// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using Optimizer.Objects;
using ColumnType = Optimizer.ColumnType;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.LinExpression;

namespace XpressExamples
{
    /// <summary>Load balancing of train wagons.</summary>
    /// <remarks>
    ///   second version, using heuristic solution as start solution for MIP
    /// </remarks>
    class Wagon
    {
        /* Box weights */
        private static readonly int[] WEIGHT = { 34, 6, 8, 17, 16, 5, 13, 21, 25, 31, 14, 13, 33, 9, 25, 25 };
        private static readonly int NBOXES = WEIGHT.Length;       /* Number of boxes                       */
        private static readonly int NWAGONS = 3;                  /* Number of wagons                      */
        private static readonly int WMAX = 100;                   /* Weight limit of the wagons            */
        private static readonly int[] HeurSol = new int[NBOXES];  /* Heuristic solution: for each box      */

        /* Variables */
        private static Variable[,] load;
        private static Variable maxweight;

        private static void model(XpressProblem p)
        {
            // Create load[box,wagon] (binary)
            load = p.AddVariables(NBOXES, NWAGONS)
                .WithType(ColumnType.Binary)
                .WithName((b, w) => $"load_{b + 1}_{w + 1}")
                .ToArray();

            // Create maxweight (continuous with lb=ceil((sum(b in BOXES) WEIGHT(b))/NBOXES)
            maxweight = p.AddVariable(
                Math.Ceiling(WEIGHT.Sum() / (double)NBOXES),
                Double.PositiveInfinity,
                ColumnType.Continuous,
                "maxweight"
            );

            // Every box into one wagon: forall(b in BOXES) sum(w in WAGONS) load(b,w) = 1
            p.AddConstraints(NBOXES,
                    b => Sum(NWAGONS, w => load[b, w]) == 1.0);

            // Limit the weight loaded into every wagon: forall(w in WAGONS) sum(b in BOXES) WEIGHT(b)*load(b,w) <= maxweight
            p.AddConstraints(NWAGONS,
                    w => Sum(NBOXES, b => WEIGHT[b] * load[b, w])
                    .Leq(maxweight).SetName($"wmilit_{w}")
                    );

            // Minimize maxweight
            p.SetObjective(maxweight, Optimizer.ObjSense.Minimize);
            p.WriteProb("Wagon.lp", "l");
        }

        /** LPT (Longest processing time) heuristic:
         * One at a time, place the heaviest unassigned
         * box onto the wagon with the least load
         * @return A dense vector with the heuristic solution.
         */
        private static double heuristic()
        {
            int[] ORDERW = new int[NBOXES];   // Box indices sorted in decreasing weight order
            int[] CurNum = new int[NWAGONS];          // For each wagon w, this is the number of boxes currently loaded
            int[] CurWeight = new int[NWAGONS];       // For each wagon w, this is the current weight, i.e. the sum of weights of loaded boxes
            int[,] Load = new int[NWAGONS, NBOXES];  // Load[w][i] (for i=0..CurNum[w]-1) contains the box index of the i-th box loaded on wagon w

            // Copy the box indices into array ORDERW and sort them in decreasing
            // order of box weights (the sorted indices are returned in array ORDERW)
            for (int b = 0; b < NBOXES; b++) ORDERW[b] = b;
            Array.Sort(ORDERW, (i1, i2) => WEIGHT[i2].CompareTo(WEIGHT[i1]));

            // Distribute the loads to the wagons using the LPT heuristic
            for (int b = 0; b < NBOXES; b++)
            {
                int v = 0;                          /* Find wagon v with the smallest load */
                for (int w = 0; w < NWAGONS; w++) if (CurWeight[w] <= CurWeight[v]) v = w;
                Load[v, CurNum[v]] = ORDERW[b];     /* Add current box to wagon v */
                CurNum[v]++;                        /* Increase the counter of boxes on v */
                CurWeight[v] += WEIGHT[ORDERW[b]];  /* Update current weight of the wagon */
            }

            // Calculate the solution value
            double heurobj = 0;                   /* heuristic solution objective value (max wagon weight) */
            for (int w = 0; w < NWAGONS; w++) if (CurWeight[w] > heurobj) heurobj = CurWeight[w];

            // Solution printing
            Console.WriteLine("Heuristic solution:");
            Console.WriteLine("Max weight: ", heurobj);

            for (int w = 0; w < NWAGONS; w++)
            {
                Console.Write($" {w + 1}:");
                for (int i = 0; i < CurNum[w]; i++) Console.Write(" " + (Load[w, i] + 1));
                Console.WriteLine($" (total weight: {CurWeight[w]:%d})");
            }

            // Save the heuristic solution into the HeurSol array
            for (int w = 0; w < NWAGONS; w++) for (int i = 0; i < CurNum[w]; i++) HeurSol[Load[w, i]] = w;
            return heurobj;
        }

        private static void optimization(XpressProblem p)
        {
            // Get the solution from the heuristic solution we have found
            // Send the solution to the optimizer
            Variable[] heurmipvar = new Variable[NBOXES];
            double[] heurmipsol = new double[NBOXES];
            for (int b = 0; b < NBOXES; b++)
            {
                heurmipvar[b] = load[b, HeurSol[b]];
                heurmipsol[b] = 1.0;
            }
            p.AddMipSol(heurmipsol, heurmipvar, "heuristic");
            p.MipOptimize();
            if (p.SolStatus != Optimizer.SolStatus.Optimal &&
                    p.SolStatus != Optimizer.SolStatus.Feasible)
                throw new Exception("failed to optimize with status " + p.SolStatus);

            Console.WriteLine(
                    $"Problem status:\n\tSolve status: {p.SolveStatus}\n\tLP status: {p.LPStatus}\n\tMIP status: {p.MIPStatus}\n\tSol status: {p.SolStatus}");

            // An integer solution has been found
            if (p.MIPStatus == Optimizer.MIPStatus.Optimal
                    || p.MIPStatus == Optimizer.MIPStatus.Solution)
            {
                double[] sol = p.GetSolution();
                Console.WriteLine($"Optimal solution:\n Max weight: {p.ObjVal}");
                for (int w = 0; w < NWAGONS; w++)
                {
                    double tot_weight = 0;
                    Console.Write($" {w + 1}:");
                    for (int b = 0; b < NBOXES; b++) if (load[b, w].GetValue(sol) * WEIGHT[b] > .5)
                        {
                            Console.Write(" {0:f}", load[b, w].GetValue(sol) * WEIGHT[b]);
                            tot_weight += load[b, w].GetValue(sol) * WEIGHT[b];
                        }
                    Console.WriteLine(" (total weight: {0:f}", tot_weight);
                }
            }
        }

        public static void Main(string[] args)
        {
            if (heuristic() <= WMAX)
            {
                Console.WriteLine("Heuristic solution fits capacity limits");
            }
            else
            {
                using (XpressProblem prob = new XpressProblem())
                {
                    /* Suppress outputs */
                    model(prob);
                    optimization(prob);
                }
            }
        }
    }
}