namespace GerbilManagerWebAPI.Genetics
{
///
/// Computes the inbreeding coefficient (Inzuchtkoeffizient) using Wright's path
/// method:
///
/// F = Σ over common ancestors, over every valid path, of
/// (1/2)^(n1 + n2 + 1) * (1 + F_A)
///
/// where, for a given common ancestor A, n1 is the number of links from the sire
/// up to A and n2 the number of links from the dam up to A, and F_A is A's own
/// inbreeding coefficient (computed recursively at full available pedigree depth).
/// A path counts only if the sire-side and dam-side chains share no individual
/// other than A — this is what prevents double counting.
///
/// This type is deliberately free of any database / EF dependency: it works purely
/// from a parent lookup, so it can be unit-tested with hand-built pedigrees.
///
public sealed class PedigreeCalculator
{
private readonly IReadOnlyDictionary _parents;
private readonly IReadOnlyDictionary _names;
// F_A memo + re-entrancy guard (the latter only matters if the data somehow
// contains a cycle, which is biologically impossible but cheap to defend).
private readonly Dictionary _fCache = new();
private readonly HashSet _fInProgress = new();
private readonly Dictionary _depthCache = new();
public PedigreeCalculator(
IReadOnlyDictionary parents,
IReadOnlyDictionary names)
{
_parents = parents;
_names = names;
}
/// True if the given individual exists in the pedigree.
public bool Contains(Guid id) => _parents.ContainsKey(id);
/// Inbreeding coefficient of an existing individual.
public InbreedingResult ForIndividual(Guid id)
{
var (father, mother) = ParentsOf(id);
return ForOffspringOf(father, mother);
}
///
/// Inbreeding coefficient of the (possibly hypothetical) offspring of the given
/// sire and dam — equivalently, the kinship between the two parents.
///
public InbreedingResult ForOffspringOf(Guid? sire, Guid? dam)
{
int generations = GenerationsAvailable(sire, dam);
if (sire is not Guid s || dam is not Guid d)
{
// A parent is unknown → no inbreeding can be established.
return new InbreedingResult(0.0, 0.0, generations, Array.Empty());
}
var contributions = Coefficient(s, d);
double f = contributions.Values.Sum();
var commonAncestors = contributions
.Select(kv => new CommonAncestorContribution(kv.Key, NameOf(kv.Key), kv.Value))
.OrderByDescending(c => c.Contribution)
.ThenBy(c => c.Name, StringComparer.OrdinalIgnoreCase)
.ToList();
return new InbreedingResult(f, f * 100.0, generations, commonAncestors);
}
///
/// Per-common-ancestor contributions to the kinship of
/// and . The sum of the values is the coefficient.
///
private Dictionary Coefficient(Guid sire, Guid dam)
{
var sirePaths = EnumerateAncestorPaths(sire);
var damPaths = EnumerateAncestorPaths(dam);
// Group the dam's paths by the ancestor they terminate at for quick lookup.
var damPathsByAncestor = damPaths
.GroupBy(p => p[^1])
.ToDictionary(g => g.Key, g => g.ToList());
var contributions = new Dictionary();
foreach (var sirePath in sirePaths)
{
var ancestor = sirePath[^1];
if (!damPathsByAncestor.TryGetValue(ancestor, out var matchingDamPaths))
{
continue;
}
var sireNodes = new HashSet(sirePath);
foreach (var damPath in matchingDamPaths)
{
// The two chains may meet only at the common ancestor itself.
if (damPath.Any(node => node != ancestor && sireNodes.Contains(node)))
{
continue;
}
int n1 = sirePath.Count - 1;
int n2 = damPath.Count - 1;
double term = Math.Pow(0.5, n1 + n2 + 1) * (1.0 + F(ancestor));
contributions[ancestor] = contributions.GetValueOrDefault(ancestor) + term;
}
}
return contributions;
}
/// Inbreeding coefficient F of a single individual (memoized).
private double F(Guid id)
{
if (_fCache.TryGetValue(id, out var cached))
{
return cached;
}
// Defend against accidental cycles in the data.
if (!_fInProgress.Add(id))
{
return 0.0;
}
double f = 0.0;
var (father, mother) = ParentsOf(id);
if (father is Guid s && mother is Guid d)
{
f = Coefficient(s, d).Values.Sum();
}
_fInProgress.Remove(id);
_fCache[id] = f;
return f;
}
///
/// Every chain of parent links starting at , including
/// the trivial chain [start] itself. Each returned list runs from the start
/// individual down to one of its ancestors (inclusive). A visited set keeps the
/// walk finite even if the data contains a cycle.
///
private List> EnumerateAncestorPaths(Guid start)
{
var all = new List>();
var path = new List();
void Walk(Guid node)
{
path.Add(node);
all.Add(new List(path));
var (father, mother) = ParentsOf(node);
foreach (var parent in new[] { father, mother })
{
if (parent is Guid p && !path.Contains(p))
{
Walk(p);
}
}
path.RemoveAt(path.Count - 1);
}
Walk(start);
return all;
}
///
/// Deepest generation of known ancestry for the offspring of the two parents
/// (parents = 1, grandparents = 2, …). 0 when neither parent is known.
///
private int GenerationsAvailable(Guid? sire, Guid? dam)
{
int FromParent(Guid? parent) => parent is Guid p ? 1 + AncestorDepth(p) : 0;
return Math.Max(FromParent(sire), FromParent(dam));
}
/// Number of ancestral generations above (0 for a founder).
private int AncestorDepth(Guid id)
{
if (_depthCache.TryGetValue(id, out var cached))
{
return cached;
}
// Temporary 0 breaks any cycle while the real value is computed.
_depthCache[id] = 0;
var (father, mother) = ParentsOf(id);
int depth = 0;
if (father is Guid f)
{
depth = Math.Max(depth, 1 + AncestorDepth(f));
}
if (mother is Guid m)
{
depth = Math.Max(depth, 1 + AncestorDepth(m));
}
_depthCache[id] = depth;
return depth;
}
private (Guid? Father, Guid? Mother) ParentsOf(Guid id) =>
_parents.TryGetValue(id, out var p) ? p : (null, null);
private string NameOf(Guid id) =>
_names.TryGetValue(id, out var name) ? name : id.ToString();
}
}