// Copyright (C) Stichting Deltares 2024. All rights reserved.
// 
// This file is part of the Dam Engine.
// 
// The Dam Engine is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// 
// All names, logos, and references to "Deltares" are registered trademarks of
// Stichting Deltares and remain full property of Stichting Deltares at all times.
// All rights reserved.

using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.DamEngine.Data.Geometry;
using Deltares.DamEngine.Data.Geotechnics;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.Calculators.KernelWrappers.Common;

public static class SoilProfile2DHelper
{
    public enum LayerType
    {
        BottomAquiferCluster,
        InBetweenAquiferCluster,
        HighestAquiferCluster,
        LowestLayer
    }

    private enum BoundaryType
    {
        Top,
        Bottom
    }

    private enum PolyLineConnectionType
    {
        SamePoint,
        SameX
    }

    private const double deviationX = 1e-06;
    private const double toleranceAlmostEqual = 1e-09;

    /// <summary>
    /// Find the intersection point between the aquifer (bottom aquifer or in-between aquifer) and the surface line if relevant.
    /// Start the search at <see cref="xStart"/> in the inward direction.
    /// Return the first intersection point found as <see cref="intersectionPoint"/>. 
    /// </summary>
    /// <param name="aquiferType">The type of aquifer: bottom aquifer or in-between aquifer. </param>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <param name="xStart">The X coordinate of the starting point for the search.</param>
    /// <param name="intersectionPoint">The intersection point of the surface line (often the ditch) with the aquifer.</param>
    /// <returns><c>true</c> if the aquifer intersects the surface line; otherwise, <c>false</c>.</returns>
    public static bool IsSurfaceLineIntersectedByAquifer(LayerType aquiferType, SoilProfile2D soilProfile, double xStart, out GeometryPoint intersectionPoint)
    {
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile);
        foreach (double xCoordinate in xCoordinates.Where(xCoordinate => xCoordinate.IsGreaterThanOrEqualTo(xStart)))
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate);
            switch (aquiferType)
            {
                case LayerType.BottomAquiferCluster when crossSection.Layers.All(layer => layer.IsAquifer):
                case LayerType.InBetweenAquiferCluster when crossSection.Layers[0].IsAquifer && crossSection.Layers.Any(layer => !layer.IsAquifer):
                    intersectionPoint = new GeometryPoint
                    {
                        X = xCoordinate,
                        Z = crossSection.TopLevel
                    };
                    return true;
            }
        }
        intersectionPoint = null;
        return false;
    }

    /// <summary>
    /// Check if at least one isolated in-between aquifer is present.  
    /// </summary>
    /// <param name="soilProfile"></param>
    /// <returns>True if an isolated in between aquifer is present, otherwise false.</returns>
    public static bool IsAtLeastOneIsolatedInBetweenAquiferPresent(SoilProfile2D soilProfile)
    {
        if (IsAquitardPresentWithAquiferInnerLoop(soilProfile))
        {
            return true;
        }

        // The inner loops are removed otherwise in case of "Aquitard inner loop in aquifer" the method
        // IsAtLeastOneIsolatedInBetweenAquiferPresent will return True because the number of aquifers in the first crossSection
        // cutting the inner loop is higher than in the previousCrossSection outside the inner loop.
        SoilProfile2D soilProfileWithoutInnerLoops = CreateSoilProfileWithoutInnerLoops(soilProfile);
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfileWithoutInnerLoops);
        var previousCount = 0;
        SoilProfile1D previousCrossSection = null;
        for (var i = 0; i < xCoordinates.Length; i++)
        {
            SoilProfile1D crossSection = soilProfileWithoutInnerLoops.GetSoilProfile1D(xCoordinates[i]);
            int currentCount = crossSection.GetInBetweenAquiferClusters?.Count ?? 0;

            if (i > 0 && (currentCount != previousCount) && IsAtLeastOneIsolatedInBetweenAquiferPresent(currentCount, previousCount, crossSection, previousCrossSection))
            {
                return true;
            }

            previousCount = currentCount;
            previousCrossSection = crossSection;
        }

        return false;
    }

    private static SoilProfile2D CreateSoilProfileWithoutInnerLoops(SoilProfile2D soilProfile)
    {
        SoilProfile2D soilProfileWithoutInnerLoops = soilProfile.Clone();

        foreach (GeometrySurface geometrySurface in soilProfileWithoutInnerLoops.Geometry.Surfaces.Where(surface => surface.InnerLoops.Count > 0))
        {
            geometrySurface.InnerLoops.Clear();
        }

        foreach (SoilLayer2D surface in soilProfileWithoutInnerLoops.Surfaces.Where(surface => surface.GeometrySurface.InnerLoops.Count > 0))
        {
            surface.GeometrySurface.InnerLoops.Clear();
        }

        return soilProfileWithoutInnerLoops;
    }

    private static bool IsAtLeastOneIsolatedInBetweenAquiferPresent(int currentCount, int previousCount, SoilProfile1D currentCrossSection, SoilProfile1D previousCrossSection)
    {
        // If more than one new or one less in-between aquifers are found, at least one is an isolated aquifer.
        if (Math.Abs(currentCount - previousCount) > 1)
        {
            return true;
        }

        if (currentCount <= 0 || currentCrossSection.GetInBetweenAquiferClusters == null)
        {
            return false;
        }

        // If one new in-between aquifer is found, the top layer in the previous cross-section must be an aquifer.
        // If not, it is an isolated aquifer.
        if (currentCount - previousCount == 1 && !previousCrossSection.Layers[0].IsAquifer)
        {
            return true;
        }

        // If one less in-between aquifer is found, the top layer in the current section must be an aquifer.
        // If not, it is an isolated aquifer.
        return previousCount - currentCount == 1 && !currentCrossSection.Layers[0].IsAquifer;
    }

    private static bool IsAquitardPresentWithAquiferInnerLoop(SoilProfile2D soilProfile2D)
    {
        foreach (GeometrySurface surface in soilProfile2D.Geometry.Surfaces.Where(surface => surface.InnerLoops.Count > 0))
        {
            SoilLayer2D soilLayerOfOuterLoop = soilProfile2D.Surfaces.First(s => s.GeometrySurface == surface);
            if (!soilLayerOfOuterLoop.IsAquifer)
            {
                foreach (GeometryLoop innerLoop in surface.InnerLoops)
                {
                    GeometrySurface geometrySurfaceOfInnerLoop = soilProfile2D.Geometry.Surfaces.First(s => s.OuterLoop.HasSameCurves(innerLoop));
                    SoilLayer2D soilLayerOfInnerLoop = soilProfile2D.Surfaces.First(s => s.GeometrySurface == geometrySurfaceOfInnerLoop);
                    if (soilLayerOfInnerLoop.IsAquifer)
                    {
                        return true;
                    }
                }
            }
        }

        return false;
    }

    /// <summary>
    /// Determines all the points of the layer boundary for the specified layer type:
    /// - For LayerType = BottomAquiferCluster, layerBoundaryTop is the top boundary of the bottom aquifer (= waternet line for PL3)
    /// - For LayerType = HighestAquiferCluster, layerBoundaryTop is the top of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet)
    /// - For LayerType = LowestLayer, layerBoundaryBottom is the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet)
    /// </summary>
    /// <param name="layerType">The type of layer boundary.</param>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <param name="layerBoundaryTop">All the points of the layer top boundary.</param>
    /// <param name="layerBoundaryBottom">All the points of the layer bottom boundary.</param>
    /// <returns>All the points of the layer boundary.</returns>
    internal static void DetermineAquiferLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, out Point2D[] layerBoundaryTop, out Point2D[] layerBoundaryBottom)
    {
        double[] xCoordinates = DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(layerType, soilProfile);
        var layerBoundaryTopPoints = new List<Point2D>();
        var layerBoundaryBottomPoints = new List<Point2D>();
        foreach (double xCoordinate in xCoordinates)
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate);
            double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection, -1);
            double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection, -1);

            if (!double.IsNaN(currentAquiferTop) && !double.IsNaN(currentAquiferBottom))
            {
                layerBoundaryTopPoints.Add(new Point2D(xCoordinate, currentAquiferTop));
                layerBoundaryBottomPoints.Add(new Point2D(xCoordinate, currentAquiferBottom));
            }
        }
        
        // Perform a short validation that the coordinates are fully defined from the beginning to the end 
        // of the profile. If not, the layer is not continuous.
        if (!IsLayerBoundaryContinuous(layerBoundaryTopPoints, xCoordinates.First(), xCoordinates.Last())
            || !IsLayerBoundaryContinuous(layerBoundaryBottomPoints, xCoordinates.First(), xCoordinates.Last()))
        {
            layerBoundaryTop = null;
            layerBoundaryBottom = null;
            return;
        }

        layerBoundaryTop = layerBoundaryTopPoints.ToArray();
        layerBoundaryBottom = layerBoundaryBottomPoints.ToArray();
    }

    private static void AddAquiferContainingPointIfNotYetAdded(SoilProfile2D soilProfile, double xCoord, double zCoord, ref List<GeometrySurface> relevantAquifers)
    {
        foreach (GeometrySurface surface in soilProfile.Geometry.Surfaces)
        {
            var polygon = new GeometryLoop();
            foreach (Point2D point in surface.OuterLoop.CalcPoints)
            {
                polygon.CalcPoints.Add(point);
            }

            if (Routines2D.CheckIfPointIsInPolygon(polygon, xCoord, zCoord) == PointInPolygon.OutsidePolygon)
            {
                continue;
            }

            if (!relevantAquifers.Contains(surface))
            {
                relevantAquifers.Add(surface);
            }

            break;
        }
    }

    /// <summary>
    /// Determines all the points of the layer boundaries (top and bottom) of all the in-between aquifers.
    /// </summary>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <param name="aquifersBoundaryTop">All the points of the layer top boundary of all the in-between aquifers.</param>
    /// <param name="aquifersBoundaryBottom">All the points of the layer bottom boundary of all the in-between aquifers.</param>
    /// <returns>All the points of the layer boundary of all the in-between aquifers.</returns>
    internal static void DetermineInBetweenAquifersLayerBoundaryPoints(SoilProfile2D soilProfile, out List<GeometryPointString> aquifersBoundaryTop, out List<GeometryPointString> aquifersBoundaryBottom)
    {
        // The inner loops are removed because they can't be in-between aquifers
        SoilProfile2D soilProfileWithoutInnerLoops = CreateSoilProfileWithoutInnerLoops(soilProfile);
        double[] xCoordinates = DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(LayerType.InBetweenAquiferCluster, soilProfileWithoutInnerLoops);

        var relevantTopAquifers = new List<GeometrySurface>();
        var relevantBottomAquifers = new List<GeometrySurface>();
        foreach (double xCoordinate in xCoordinates)
        {
            SoilProfile1D crossSection = soilProfileWithoutInnerLoops.GetSoilProfile1D(xCoordinate);
            if (crossSection.GetInBetweenAquiferClusters.Count > 0)
            {
                for (var i = 0; i < crossSection.GetInBetweenAquiferClusters.Count; i++)
                {
                    double zLevelTop = crossSection.GetInBetweenAquiferClusters[i].Item1.TopLevel - GeometryConstants.Accuracy;
                    double zLevelBottom = crossSection.GetInBetweenAquiferClusters[i].Item2.BottomLevel + GeometryConstants.Accuracy;
                    AddAquiferContainingPointIfNotYetAdded(soilProfileWithoutInnerLoops, xCoordinate, zLevelTop, ref relevantTopAquifers);
                    AddAquiferContainingPointIfNotYetAdded(soilProfileWithoutInnerLoops, xCoordinate, zLevelBottom, ref relevantBottomAquifers);
                }
            }
        }

        List<GeometryPointString> isolatedLayerBoundaryTop = [];
        List<GeometryPointString> isolatedLayerBoundaryBottom = [];
        aquifersBoundaryTop = [];
        aquifersBoundaryBottom = [];
        if (relevantTopAquifers != null && relevantBottomAquifers != null)
        {
            isolatedLayerBoundaryTop.AddRange(relevantTopAquifers.Select(surface => surface.DetermineTopGeometrySurface()));
            isolatedLayerBoundaryBottom.AddRange(relevantBottomAquifers.Select(surface => surface.DetermineBottomGeometrySurface()));
            aquifersBoundaryTop = ConnectPolyLines(isolatedLayerBoundaryTop);
            aquifersBoundaryBottom = ConnectPolyLines(isolatedLayerBoundaryBottom);
        }
    }

    private static double[] DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(LayerType layerType, SoilProfile2D soilProfile)
    {
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile);
        var xCoordinatesAll = new List<double>();
        foreach (double xCoordinate in xCoordinates)
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate);
            int count = layerType == LayerType.InBetweenAquiferCluster ? crossSection.GetInBetweenAquiferClusters.Count : 1;
            if (count == 0)
            {
                xCoordinatesAll.Add(xCoordinate);
            }
            else
            {
                for (var i = 0; i < count; i++)
                {
                    if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, i))
                    {
                        xCoordinatesAll.Add(xCoordinate - deviationX);
                        xCoordinatesAll.Add(xCoordinate + deviationX);
                    }
                    else
                    {
                        xCoordinatesAll.Add(xCoordinate);
                    }
                }
            }
        }
        return xCoordinatesAll.ToArray();
    }

    internal static List<GeometryPointString> ConnectPolyLines(List<GeometryPointString> polyLines)
    {
        List<GeometryPointString> connectedPolyLinesSamePoint = ConnectPolyLines(PolyLineConnectionType.SamePoint, polyLines);
        List<GeometryPointString> connectedPolyLinesSameX = ConnectPolyLines(PolyLineConnectionType.SameX, connectedPolyLinesSamePoint);
        return connectedPolyLinesSameX;
    }

    private static bool CanPolyLinesBeConnected(PolyLineConnectionType connectionType, GeometryPointString currentPolyLine, GeometryPointString nextPolyLine, List<GeometryPointString> polyLines)
    {
        if (connectionType == PolyLineConnectionType.SamePoint)
        {
            return currentPolyLine.CalcPoints.Last().LocationEquals(nextPolyLine.CalcPoints.First());
        }
        bool isPolyLineUnique = polyLines.FindAll(p => p.CalcPoints.First().X.IsNearEqual(nextPolyLine.CalcPoints.First().X, toleranceAlmostEqual)).Count == 1;
        return currentPolyLine.CalcPoints.Last().X.IsNearEqual(nextPolyLine.CalcPoints.First().X, toleranceAlmostEqual) && isPolyLineUnique;

    }

    /// <summary>
    /// Connect the poly-lines in the list when possible. Two types of connection are available:
    /// - SamePoint = the end point of a poly-line coincides with the start point of another poly-line
    /// - SameX = the X_end point of a poly-line coincides with the X_start point of another poly-line (connection is done only is
    /// this other poly-line is unique)
    /// </summary>
    /// <param name="connectionType">The type of connection for the poly-lines (same point or same X).</param>
    /// <param name="polyLines">The list of poly-lines that must be connected when possible.</param>
    /// <returns>A list of connected poly-lines.</returns>
    private static List<GeometryPointString> ConnectPolyLines(PolyLineConnectionType connectionType, List<GeometryPointString> polyLines)
    {
        var connectedPolyLines = new List<GeometryPointString>();
        var usedPolyLines = new HashSet<GeometryPointString>();

        foreach (GeometryPointString polyLine in polyLines)
        {
            if (usedPolyLines.Contains(polyLine))
                continue;

            var currentPolyLine = new GeometryPointString();
            currentPolyLine.CalcPoints.AddRange(polyLine.CalcPoints);
            usedPolyLines.Add(polyLine);

            bool isConnectionFound;
            do
            {
                isConnectionFound = IsConnectionFound(connectionType, polyLines, usedPolyLines, currentPolyLine);
            } while (isConnectionFound);

            currentPolyLine.SyncPoints();
            connectedPolyLines.Add(currentPolyLine);
        }
        
        return connectedPolyLines;
    }

    private static bool IsConnectionFound(PolyLineConnectionType connectionType, List<GeometryPointString> polyLines, HashSet<GeometryPointString> usedPolyLines, GeometryPointString currentPolyLine)
    {
        var isConnectionFound = false;
        foreach (GeometryPointString nextPolyLine in polyLines)
        {
            if (usedPolyLines.Contains(nextPolyLine))
                continue;
                    
            if (CanPolyLinesBeConnected(connectionType, currentPolyLine, nextPolyLine, polyLines))
            {
                currentPolyLine.CalcPoints.AddRange(connectionType == PolyLineConnectionType.SamePoint ? nextPolyLine.CalcPoints.Skip(1) : nextPolyLine.CalcPoints);
                usedPolyLines.Add(nextPolyLine);
                isConnectionFound = true;
                break;
            }
        }

        return isConnectionFound;
    }

    /// <summary>
    /// Determine all the xCoordinates to make cross-sections. 
    /// </summary>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <returns>All the xCoordinates of the soil profile 2D.</returns>
    private static double[] DetermineAllXCoordinatesOfSoilProfile(SoilProfile2D soilProfile)
    {
        IEnumerable<Point2D> points = soilProfile.Surfaces.SelectMany(surf => surf.GeometrySurface.OuterLoop.CalcPoints);
        double[] xCoordinates = points.Select(point => point.X).OrderBy(x => x).Distinct().ToArray();
        return xCoordinates;
    }
    
    private static bool IsLayerBoundaryContinuous(IEnumerable<Point2D> boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary)
    {
        IEnumerable<Point2D> point2Ds = boundaryPoints.ToList();
        return point2Ds.Any() && (Math.Abs(point2Ds.First().X - leftGeometryBoundary) < toleranceAlmostEqual)
                              && Math.Abs(point2Ds.Last().X - rightGeometryBoundary) < toleranceAlmostEqual;
    }
    
    private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer)
    {
        SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D, indexInBetweenAquifer);
        if (layer != null)
        {
            return boundaryType == BoundaryType.Top ? layer.TopLevel : layer.BottomLevel;
        }

        return double.NaN;
    }

    private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer)
    {
        switch (layerType)
        {
            case LayerType.BottomAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.BottomAquiferLayer : soilProfile1D.DeepestAquiferLayer;
            case LayerType.InBetweenAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item1 : soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item2;
            case LayerType.HighestAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.GetHighestAquifer() : soilProfile1D.GetLowestLayerOfHighestAquiferCluster();
            case LayerType.LowestLayer:
                return soilProfile1D.Layers.Last();
        }

        return null;
    }

    private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer)
    {
        SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX);
        SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX);
        if (crossSectionLeft == null || crossSectionRight == null)
        {
            return false;
        }

        double aquiferLeftTop = GetLevel(layerType, BoundaryType.Top, crossSectionLeft, indexInBetweenAquifer);
        double aquiferLefBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionLeft, indexInBetweenAquifer);
        double aquiferRightTop = GetLevel(layerType, BoundaryType.Top, crossSectionRight, indexInBetweenAquifer);
        double aquiferRightBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionRight, indexInBetweenAquifer);
        if (!double.IsNaN(aquiferLeftTop) && !double.IsNaN(aquiferLefBottom) && !double.IsNaN(aquiferRightTop) && !double.IsNaN(aquiferRightBottom))
        {
            return Math.Abs(aquiferLeftTop - aquiferRightTop) > GeometryConstants.Accuracy || Math.Abs(aquiferLefBottom - aquiferRightBottom) > GeometryConstants.Accuracy;
        }

        return false;
    }
}