// (c) 2023-2025 Fair Isaac Corporation

import static com.dashoptimization.objects.Utils.scalarProduct;
import static com.dashoptimization.objects.Utils.sum;

import com.dashoptimization.ColumnType;
import com.dashoptimization.DefaultMessageListener;
import com.dashoptimization.IntHolder;
import com.dashoptimization.XPRSconstants;
import com.dashoptimization.XPRSenumerations;
import com.dashoptimization.objects.LinExpression;
import com.dashoptimization.objects.LinTermMap;
import com.dashoptimization.objects.Variable;
import com.dashoptimization.objects.XpressProblem;
import com.dashoptimization.objects.XpressProblem.CallbackAPI.CutRoundCallback;

/** Demonstrates how to implement cutting planes as part of a MIP
 * branch-and-bound search using the cutround callback.
 *
 * Cuts are added as user cuts using XpressProblem.addManagedCut().
 *
 * Economic lot sizing problem. Solved by adding (l,S)-inequalities
 * in a branch-and-cut heuristic (using the cutround callback).
 *
 * ELS considers production planning over a horizon of T periods. In period t,
 * t=1,...,T, there is a given demand DEMAND[p,t] that must be satisfied by
 * production produce[p,t] in period t and by inventory carried over
 * from previous periods.
 * There is a set-up cost SETUPCOST[t] associated with production in
 * period t and the total production capacity per period is limited. The unit
 * production cost in period t is PRODCOST[p,t]. There is no
 * inventory or stock-holding cost.
 *
 * A well-known class of valid inequalities for ELS are the
 * (l,S)-inequalities.  Let D(p,q) denote the demand in periods p
 * to q and y(t) be a binary variable indicating whether there is any
 * production in period t.  For each period l and each subset of periods S
 * of 1 to l, the (l,S)-inequality is
 * <pre>
 *      sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(t,l) * y(t)
 *          &gt;= D(1,l)
 * </pre>
 *
 * It says that actual production x(t) in periods included S plus maximum
 * potential production D(t,l)*y(t) in the remaining periods (those not in
 * S) must at least equal total demand in periods 1 to l.  Note that in
 * period t at most D(t,l) production is required to meet demand up to
 * period l.
 *
 * Based on the observation that
 * <pre>
 *      sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(t,l) * y(t)
 *          &gt;= sum (t=1:l) min(x(t), D(t,l) * y(t))
 *          &gt;= D(1,l)
 * </pre>
 * it is easy to develop a separation algorithm and thus a cutting plane
 * algorithm based on these (l,S)-inequalities.
 */
public class ELSManagedCuts {
    /** Tolerance for satisfiability. */
    private static final double EPS = 1e-6;

    /** Number of time periods. */
    private static final int DIM_TIMES = 15;
    /** Number of products to produce. */
    private static final int DIM_PRODUCTS = 4;

    /** Demand per product (first dim) and time period (second dim). */
    private static final int[][] DEMAND = new int[][]{
        new int[]{2, 3, 5, 3, 4, 2, 5, 4, 1, 3, 4, 2, 3, 5, 2},
        new int[]{3, 1, 2, 3, 5, 3, 1, 2, 3, 3, 4, 5, 1, 4, 1},
        new int[]{3, 5, 2, 1, 2, 1, 3, 3, 5, 2, 2, 1, 3, 2, 3},
        new int[]{2, 2, 1, 3, 2, 1, 2, 2, 3, 3, 2, 2, 3, 1, 2}
    };
    /** Setup cost. */
    private static final double[] SETUPCOST = new double[]{
        17, 14, 11, 6, 9, 6, 15, 10, 8, 7, 12, 9, 10, 8, 12
    };
    /** Production cost per product (first dim) and time period (second dim). */
    private static final double[][] PRODCOST = new double[][]{
        new double[]{5, 3, 2, 1, 3, 1, 4, 3, 2, 2, 3, 1, 2, 3, 2},
        new double[]{1, 4, 2, 3, 1, 3, 1, 2, 3, 3, 3, 4, 4, 2, 2},
        new double[]{3, 3, 3, 4, 4, 3, 3, 3, 2, 2, 1, 1, 3, 3, 3},
        new double[]{2, 2, 2, 3, 3, 3, 4, 4, 4, 3, 3, 2, 2, 2, 3}
    };
    /** Capacity. */
    private static int[] CAP = new int[]{
        12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12
    };

    /** Cut round callback for separating our cutting planes. */
    static class Callback implements CutRoundCallback {
        private Variable[][] produce;
        private Variable[][] setup;
        private double[][][] sumDemand;
        public Callback(Variable[][] produce, Variable[][] setup, double[][][] sumDemand) {
            this.produce = produce;
            this.setup = setup;
            this.sumDemand = sumDemand;
        }
        public void cutRound(XpressProblem prob, int ifxpresscuts, IntHolder action) {
            // Apply only one round of cutting on each node.
            // Because the CutRound callback is fired *before* a round of
            // cutting, the CUTROUNDS attribute will start from 0 on the first
            // invocation of the callback.
            if (prob.attributes().getCutRounds() >= 1)
                return;

            // Get the solution vector.
            // Xpress will only fire the CutRound callback when a solution is
            // available, so there is no need to check whether a solution is
            // available.
            double[] sol = prob.getCallbackSolution();

            // Search for violated constraints : the minimum of the actual
            // production produce[p][t] and the maximum potential production
            // D[p][t][l]*setup[p][t] in periods 0 to l must at least equal
            // the total demand in periods 0 to l.
            //    sum(t=1:l) min(prod[p][t], D[p][t][l]*setup[p][t]) >= D[p][0][l]
            int nCutsAdded = 0;
            for (int p = 0; p < DIM_PRODUCTS; p++) {
                for (int l = 0; l < DIM_TIMES; l++) {
                    double sum = 0.0;
                    for (int t = 0; t <= l; t++) {
                        if (produce[p][t].getValue(sol) < sumDemand[p][t][l] * setup[p][t].getValue(sol) + EPS) {
                            sum += produce[p][t].getValue(sol);
                        }
                        else {
                            sum += sumDemand[p][t][l] * setup[p][t].getValue(sol);
                        }
                    }
                    if (sum < sumDemand[p][0][l] - EPS) {
                        // Create the violated cut.
                        LinExpression cut = LinExpression.create();
                        for (int t = 0; t <= l; t++) {
                            if (produce[p][t].getValue(sol) < sumDemand[p][t][l] * setup[p][t].getValue(sol)) {
                                cut.addTerm(1.0, produce[p][t]);
                            }
                            else {
                                cut.addTerm(sumDemand[p][t][l], setup[p][t]);
                            }
                        }

                        // Give the cut to Xpress to manage.
                        // It will automatically be presolved.
                        prob.addManagedCut(true, cut.geq(sumDemand[p][0][l]));
                        ++nCutsAdded;

                        // If we modified the problem in the callback, Xpress
                        // will automatically trigger another roound of cuts,
                        // so there is no need to set the action return
                        // argument.
                    }
                }
            }
            if (nCutsAdded > 0) {
                System.out.printf("Cuts added: %d%n", nCutsAdded);
            }
        }
    }

    public static void main(String[] args) {
        try (XpressProblem prob = new XpressProblem()) {
            prob.callbacks.addMessageCallback(DefaultMessageListener::console);

            // Create an economic lot sizing problem.
            // Calculate demand D(p, s, t) as the demand for product p from
            // time s to time t(inclusive).
            double[][][] sumDemand = new double[DIM_PRODUCTS][DIM_TIMES][DIM_TIMES];
            for (int p = 0; p < DIM_PRODUCTS; p++) {
                for (int s = 0; s < DIM_TIMES; s++) {
                    double thisSumDemand = 0.0;
                    for (int t = s; t < DIM_TIMES; t++) {
                        thisSumDemand += DEMAND[p][t];
                        sumDemand[p][s][t] = thisSumDemand;
                    }
                }
            }

            Variable[][] produce = prob.addVariables(DIM_PRODUCTS, DIM_TIMES)
                .withName((p,t) -> "produce(" + p + "," + t + ")")
                .toArray();
            Variable[][] setup = prob.addVariables(DIM_PRODUCTS, DIM_TIMES)
                .withType(ColumnType.Binary)
                .withName((p,t) -> "setup(" + p + "," + t + ")")
                .toArray();

            // Add the objective function :
            // MinCost:= sum(t in TIMES) (SETUPCOST(t) * sum(p in PRODUCTS) setup(p,t) +
            //           sum(p in PRODUCTS) PRODCOST(p, t) *produce(p, t) )
            prob.setObjective(sum(DIM_TIMES, DIM_PRODUCTS,
                                  (t, p) -> setup[p][t].mul(SETUPCOST[t]).plus(produce[p][t].mul(PRODCOST[p][t]))));

            // Add constraints.
            // Production in periods 0 to t must satisfy the total demand
            // during this period of time, for all t in [0,T[
            //    forall(p in PRODUCTS, t in TIMES)
            //      Dem(t) : = sum(s in 1..t) produce(p,s) >= sum(s in 1..t) DEMAND(s)
            prob.addConstraints(DIM_PRODUCTS, DIM_TIMES,
                                (p, t) -> sum(t + 1, s -> produce[p][s]).geq(sumDemand[p][0][t]));

            // If there is production during t then there is a setup in t :
            //    forall(p in PRODUCTS, t in TIMES)
            //       ProdSetup(t) : = produce(t) <= D(t,TIMES) * setup(t)
            for (int p = 0; p < DIM_PRODUCTS; ++p) {
                for (int t = DIM_TIMES - 1; t >= 0; --t) {
                    prob.addConstraint(produce[p][t].leq(setup[p][t].mul(sumDemand[p][t][DIM_TIMES - 1])));
                }
            }

            // Capacity limits :
            //    forall(t in TIMES) Capacity(t) : = sum(p in PRODUCTS) produce(p, t) <= CAP(t)
            prob.addConstraints(DIM_TIMES,
                                t -> sum(DIM_PRODUCTS, p -> produce[p][t]).leq(CAP[t]));

            // Add a CutRound callback for separating our cuts.
            prob.callbacks.addCutRoundCallback(new Callback(produce, setup, sumDemand));

            prob.writeProb("elsmanagedcuts.lp");

            prob.optimize();

            if (prob.attributes().getSolveStatus() == XPRSenumerations.SolveStatus.COMPLETED &&
                prob.attributes().getSolStatus() == XPRSenumerations.SolStatus.OPTIMAL) {
                System.out.println("Solved problem to optimality.");
            }
            else {
                System.out.printf("Failed to solve problem with solvestatus %s and solstatus %s%n",
                                  prob.attributes().getSolveStatus(),
                                  prob.attributes().getSolStatus());
            }
        }
    }
}