GerbilManager/GerbilManagerWebAPI/Genetics/PedigreeCalculator.cs

namespace GerbilManagerWebAPI.Genetics
{
    /// <summary>
    /// 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.
    /// </summary>
    public sealed class PedigreeCalculator
    {
        private readonly IReadOnlyDictionary<Guid, (Guid? Father, Guid? Mother)> _parents;
        private readonly IReadOnlyDictionary<Guid, string> _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<Guid, double> _fCache = new();
        private readonly HashSet<Guid> _fInProgress = new();
        private readonly Dictionary<Guid, int> _depthCache = new();

        public PedigreeCalculator(
            IReadOnlyDictionary<Guid, (Guid? Father, Guid? Mother)> parents,
            IReadOnlyDictionary<Guid, string> names)
        {
            _parents = parents;
            _names = names;
        }

        /// <summary>True if the given individual exists in the pedigree.</summary>
        public bool Contains(Guid id) => _parents.ContainsKey(id);

        /// <summary>Inbreeding coefficient of an existing individual.</summary>
        public InbreedingResult ForIndividual(Guid id)
        {
            var (father, mother) = ParentsOf(id);
            return ForOffspringOf(father, mother);
        }

        /// <summary>
        /// Inbreeding coefficient of the (possibly hypothetical) offspring of the given
        /// sire and dam — equivalently, the kinship between the two parents.
        /// </summary>
        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<CommonAncestorContribution>());
            }

            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);
        }

        /// <summary>
        /// Per-common-ancestor contributions to the kinship of <paramref name="sire"/>
        /// and <paramref name="dam"/>. The sum of the values is the coefficient.
        /// </summary>
        private Dictionary<Guid, double> 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<Guid, double>();

            foreach (var sirePath in sirePaths)
            {
                var ancestor = sirePath[^1];
                if (!damPathsByAncestor.TryGetValue(ancestor, out var matchingDamPaths))
                {
                    continue;
                }

                var sireNodes = new HashSet<Guid>(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;
        }

        /// <summary>Inbreeding coefficient F of a single individual (memoized).</summary>
        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;
        }

        /// <summary>
        /// Every chain of parent links starting at <paramref name="start"/>, 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.
        /// </summary>
        private List<List<Guid>> EnumerateAncestorPaths(Guid start)
        {
            var all = new List<List<Guid>>();
            var path = new List<Guid>();

            void Walk(Guid node)
            {
                path.Add(node);
                all.Add(new List<Guid>(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;
        }

        /// <summary>
        /// Deepest generation of known ancestry for the offspring of the two parents
        /// (parents = 1, grandparents = 2, …). 0 when neither parent is known.
        /// </summary>
        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));
        }

        /// <summary>Number of ancestral generations above <paramref name="id"/> (0 for a founder).</summary>
        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();
    }
}