using System;
using System.Collections.Generic;
using System.Linq;

namespace HydraulicSimulator.Models
{
    /// <summary>
    /// Resultado de la simulación de la red hidráulica
    /// </summary>
    public class SolutionResult
    {
        public bool Converged { get; set; }
        public Dictionary<string, double> Flows { get; set; } = new Dictionary<string, double>();
        public Dictionary<string, double> Pressures { get; set; } = new Dictionary<string, double>();
        public int Iterations { get; set; }
        public double Residual { get; set; }
        public string ErrorMessage { get; set; } = "";

        public SolutionResult(bool converged = false)
        {
            Converged = converged;
        }
    }

    /// <summary>
    /// Red hidráulica y solver
    /// </summary>
    public class HydraulicNetwork
    {
        public Fluid Fluid { get; set; }
        public Dictionary<string, Node> Nodes { get; set; } = new Dictionary<string, Node>();
        public List<Branch> Branches { get; set; } = new List<Branch>();

        public HydraulicNetwork(Fluid fluid = null)
        {
            Fluid = fluid ?? Fluid.Water20C;
        }

        public void AddNode(string name, double? fixedP = null)
        {
            if (Nodes.ContainsKey(name))
                return;

            if (fixedP == null)
            {
                Nodes[name] = new Node(name, fixedP: false, p: 0.0);
            }
            else
            {
                Nodes[name] = new Node(name, fixedP: true, p: fixedP.Value);
            }
        }

        public void AddElement(Element element, string n1, string n2, string name = "")
        {
            if (!Nodes.ContainsKey(n1) || !Nodes.ContainsKey(n2))
                throw new ArgumentException("Nodos inexistentes");

            var branchName = string.IsNullOrEmpty(name) ? $"{n1}->{n2}" : name;
            Branches.Add(new Branch(n1, n2, new List<Element> { element }, branchName));
        }

        public void AddBranch(string n1, string n2, List<Element> elements, string name = "")
        {
            if (!Nodes.ContainsKey(n1) || !Nodes.ContainsKey(n2))
                throw new ArgumentException("Nodos inexistentes");

            Branches.Add(new Branch(n1, n2, elements, name));
        }

        public SolutionResult Solve(int maxIterations = 100, double tolerance = 1e-3, 
                                   double relaxationFactor = 0.1, bool verbose = false)
        {
            try
            {
                var names = Nodes.Keys.ToList();
                var free = names.Where(n => !Nodes[n].FixedP).ToList();
                var fixedNodes = names.Where(n => Nodes[n].FixedP).ToList();
                var idxFree = free.Select((n, i) => new { n, i }).ToDictionary(x => x.n, x => x.i);

                // Inicial: presiones ~ promedio de fijos
                double pRef = 0.0;
                if (fixedNodes.Any())
                {
                    pRef = fixedNodes.Average(n => Nodes[n].P);
                }

                foreach (var n in free)
                {
                    Nodes[n].P = pRef;
                }

                // Inicializar caudales a valores pequeños
                foreach (var b in Branches)
                {
                    b.Q = 0.001; // m³/s inicial pequeño
                }

                double normR = 0.0;
                int it = 0;

                // Iteración global sobre presiones nodales
                for (it = 0; it < maxIterations; it++)
                {
                    // 1) con presiones actuales, resolvés q de cada rama
                    foreach (var b in Branches)
                    {
                        var dpTarget = Nodes[b.N1].P - Nodes[b.N2].P;
                        b.SolveFlowGivenDp(dpTarget, Fluid);
                    }

                    // 2) residuos de continuidad en nodos libres
                    var R = new double[free.Count];
                    var G = new double[free.Count, free.Count]; // Jacobiano nodal

                    foreach (var b in Branches)
                    {
                        var i = b.N1;
                        var j = b.N2;
                        // sensibilidad dq/d(dp) = 1 / d(ΔP)/dq
                        var ddpDq = b.DdpDq(b.Q, Fluid);

                        // Evitar división por cero y valores demasiado grandes
                        ddpDq = Math.Max(1e-10, Math.Min(1e10, Math.Abs(ddpDq)));
                        var k = 1.0 / ddpDq;

                        // aporte a residuos (signo: positiva saliendo de nodo n1)
                        if (idxFree.ContainsKey(i))
                        {
                            var idx_i = idxFree[i];
                            R[idx_i] += b.Q;
                            G[idx_i, idx_i] += k;
                            if (idxFree.ContainsKey(j))
                            {
                                var idx_j = idxFree[j];
                                G[idx_i, idx_j] -= k;
                            }
                        }
                        if (idxFree.ContainsKey(j))
                        {
                            var idx_j = idxFree[j];
                            R[idx_j] -= b.Q;
                            G[idx_j, idx_j] += k;
                            if (idxFree.ContainsKey(i))
                            {
                                var idx_i = idxFree[i];
                                G[idx_j, idx_i] -= k;
                            }
                        }
                    }

                    normR = R.Length > 0 ? R.Max(Math.Abs) : 0.0;
                    if (verbose)
                        Console.WriteLine($"it {it}: |R|_inf={normR:E3}");

                    if (normR < tolerance)
                        break;

                    if (free.Count > 0)
                    {
                        // 3) resolver actualización de presiones
                        var dpUpdate = SolveLinearSystem(G, R, free.Count);

                        // Limitar la actualización de presión
                        var maxDp = 50000.0; // 50 kPa máximo por iteración
                        for (int k = 0; k < dpUpdate.Length; k++)
                        {
                            dpUpdate[k] = Math.Max(-maxDp, Math.Min(maxDp, dpUpdate[k]));
                        }

                        // amortiguación más conservadora
                        var relax = it < 10 ? relaxationFactor : relaxationFactor * 0.5;
                        for (int k = 0; k < dpUpdate.Length; k++)
                        {
                            dpUpdate[k] *= relax;
                        }

                        foreach (var kvp in idxFree)
                        {
                            Nodes[kvp.Key].P += dpUpdate[kvp.Value];
                        }
                    }
                }

                // Preparar resultado
                var result = new SolutionResult(normR < tolerance)
                {
                    Iterations = it,
                    Residual = normR
                };

                // Llenar flows y pressures
                foreach (var b in Branches)
                {
                    result.Flows[b.Name] = b.Q;
                }

                foreach (var kvp in Nodes)
                {
                    result.Pressures[kvp.Key] = kvp.Value.P;
                }

                return result;
            }
            catch (Exception ex)
            {
                return new SolutionResult(false)
                {
                    ErrorMessage = ex.Message,
                    Iterations = 0,
                    Residual = double.MaxValue
                };
            }
        }

        private double[] SolveLinearSystem(double[,] G, double[] R, int n)
        {
            // Implementación simple de eliminación gaussiana con pivoteo parcial
            var A = new double[n, n + 1];
            
            // Copiar G y -R al sistema aumentado
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    A[i, j] = G[i, j];
                }
                A[i, n] = -R[i];
            }

            // Regularización
            for (int i = 0; i < n; i++)
            {
                A[i, i] += 1e-6;
            }

            // Eliminación gaussiana
            for (int i = 0; i < n; i++)
            {
                // Encontrar el pivote
                int maxRow = i;
                for (int k = i + 1; k < n; k++)
                {
                    if (Math.Abs(A[k, i]) > Math.Abs(A[maxRow, i]))
                        maxRow = k;
                }

                // Intercambiar filas
                for (int k = i; k <= n; k++)
                {
                    var temp = A[maxRow, k];
                    A[maxRow, k] = A[i, k];
                    A[i, k] = temp;
                }

                // Hacer todos los elementos debajo del pivote igual a 0
                for (int k = i + 1; k < n; k++)
                {
                    if (Math.Abs(A[i, i]) < 1e-12)
                        A[i, i] = 1e-12;
                    
                    var c = A[k, i] / A[i, i];
                    for (int j = i; j <= n; j++)
                    {
                        if (i == j)
                            A[k, j] = 0;
                        else
                            A[k, j] -= c * A[i, j];
                    }
                }
            }

            // Sustitución hacia atrás
            var solution = new double[n];
            for (int i = n - 1; i >= 0; i--)
            {
                solution[i] = A[i, n];
                for (int j = i + 1; j < n; j++)
                {
                    solution[i] -= A[i, j] * solution[j];
                }
                
                if (Math.Abs(A[i, i]) < 1e-12)
                    A[i, i] = 1e-12;
                    
                solution[i] = solution[i] / A[i, i];
            }

            return solution;
        }

        /// <summary>
        /// Genera un reporte de la solución
        /// </summary>
        public void Report()
        {
            Console.WriteLine("== Nodos (Pa) ==");
            foreach (var kvp in Nodes)
            {
                var node = kvp.Value;
                var kind = node.FixedP ? "FIX" : "FREE";
                Console.WriteLine($"{node.Name,10}: {node.P,12:F1}  [{kind}]");
            }

            Console.WriteLine("\n== Ramas (q m³/s) ==");
            foreach (var b in Branches)
            {
                Console.WriteLine($"{b.Name,15}: {b.Q,10:E6}");
            }
        }
    }
}