Hooking Your Solver to ojAlgo

import static org.ojalgo.function.constant.PrimitiveMath.*;

import java.util.ArrayList;
import java.util.List;

import org.ojalgo.OjAlgoUtils;
import org.ojalgo.equation.Equation;
import org.ojalgo.matrix.store.PhysicalStore;
import org.ojalgo.matrix.store.R064Store;
import org.ojalgo.matrix.task.iterative.ConjugateGradientSolver;
import org.ojalgo.netio.BasicLogger;
import org.ojalgo.optimisation.Expression;
import org.ojalgo.optimisation.ExpressionsBasedModel;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.Variable;
import org.ojalgo.optimisation.convex.ConvexSolver;
import org.ojalgo.optimisation.integer.IntegerSolver;
import org.ojalgo.optimisation.linear.LinearSolver;
import org.ojalgo.structure.Structure1D.IntIndex;
import org.ojalgo.structure.Structure2D.IntRowColumn;

/**
 * An example of how to:
 * <ol>
 * <li>implement a simple solver, and
 * <li>hook that up to be used with {@link ExpressionsBasedModel}.
 * </ol>
 *
 * @see https://www.ojalgo.org/2025/02/hooking-your-solver-to-ojalgo/
 */
public class HookingUpSolver {

    /**
     * The {@link Optimisation.Integration} interface actually declares a few more methods than what we
     * implement here. Since we extend the abstract {@link ExpressionsBasedModel.Integration} class, and use
     * it the standard way, several things are taken care of for us.
     */
    static final class UnconstrainedIntegration extends ExpressionsBasedModel.Integration<UnconstrainedSolver> {

        /**
         * Turn parameter scaling on/off
         */
        private static final boolean ADJUSTED = false;

        /**
         * Build the solver instance.
         */
        @Override
        public HookingUpSolver.UnconstrainedSolver build(final ExpressionsBasedModel model) {

            List<Variable> variables = model.getVariables();

            Expression objective = model.objective();

            /*
             * In this case we know that the objective function is quadratic, and that there are no
             * constraints. We can use this information to build a simple equation system. Essentially we take
             * the partial derivatives of the objective function with respect to each variable, and create a
             * system of linear equations by equating each of them to 0.0.
             */

            List<Equation> equations = new ArrayList<>(variables.size());
            for (int i = 0, limit = variables.size(); i < limit; i++) {
                equations.add(Equation.sparse(i, limit));
            }

            for (IntRowColumn quadraticKey : objective.getQuadraticKeySet()) {
                double quadraticValue = objective.doubleValue(quadraticKey, ADJUSTED);
                equations.get(quadraticKey.row).add(quadraticKey.column, quadraticValue);
                equations.get(quadraticKey.column).add(quadraticKey.row, quadraticValue);
            }

            for (IntIndex linearKey : objective.getLinearKeySet()) {
                double linearValue = objective.doubleValue(linearKey, ADJUSTED);
                equations.get(linearKey.index).setRHS(-linearValue);
            }

            return new UnconstrainedSolver(equations, model.options);
        }

        /**
         * Can the solver, built by this integration, handle the model?
         */
        @Override
        public boolean isCapable(final ExpressionsBasedModel model) {

            /*
             * Verify that the model is unconstrained and that the objective function is quadratic (not just
             * linear). ExpressionsBasedModel cannot model expressions of higher order than 2, and we check
             * that at least one factor is quadratic – that's NOT enough to guarantee that the problem is
             * solvable by this solver, but it's a good start. To be thorough we should check the eigenvalues
             * of the Hessian.
             */
            return model.constraints().count() == 0 && model.objective().isAnyQuadraticFactorNonZero();
        }

    }

    /**
     * Solves unconstrained optimisation problems. A very simple implementation that only works for (some)
     * quadratic problems modelled in {@link ExpressionsBasedModel}.
     * <p>
     * Actually this solver solves an equation system, and doesn't know where that comes from. Transforming
     * the unconstrained optimisation model into an equation system that this solver can solve is done by the
     * {@link UnconstrainedIntegration} class.
     * <p>
     * If this solver should be usable without that integration class (without {@link ExpressionsBasedModel}),
     * it would need some other way to build the equation system from a multivariate function or similar.
     */
    static final class UnconstrainedSolver implements Optimisation.Solver {

        /**
         * Requires the equation system body matrix to be symmetric and positive-definite. So this wont be
         * able to solve anything we throw at it
         */
        private final ConjugateGradientSolver myDelegateSolver = new ConjugateGradientSolver();
        private final List<Equation> myEquations;

        UnconstrainedSolver(final List<Equation> equations, final Optimisation.Options options) {
            super();
            myEquations = equations;
            myDelegateSolver.configurator().debug(options.logger_appender);
        }

        @Override
        public Result solve(final Result kickStarter) {

            PhysicalStore<Double> solution = R064Store.FACTORY.column(kickStarter);

            /*
             * Solve the equation system, updating the solution vector, returning an error/accuracy measure.
             */
            double error = myDelegateSolver.resolve(myEquations, solution);

            /*
             * Is the solution good enough?
             */
            State state = error < HUNDREDTH ? Optimisation.State.OPTIMAL : Optimisation.State.APPROXIMATE;

            /*
             * The objective function value
             */
            double value = this.evaluate(solution);

            return Result.wrap(solution).withState(state).withValue(value);
        }

        /**
         * Evaluate the objective function at the given solution.
         */
        private double evaluate(final PhysicalStore<Double> solution) {

            double tmpVal;
            double retVal = ZERO;

            for (Equation equation : myEquations) {
                tmpVal = equation.dot(solution) * HALF;
                tmpVal -= equation.getRHS();
                retVal += tmpVal * solution.doubleValue(equation.index);
            }

            return retVal;
        }

    }

    /**
     * Integrations absolutely must be stateless, and you typically only ever need one instance.
     */
    static final ExpressionsBasedModel.Integration<UnconstrainedSolver> INTEGRATION = new UnconstrainedIntegration();

    /**
     * Solve a small problem you'll find towards the end of this document:
     * https://www.ece.mcmaster.ca/~xwu/part4.pdf
     * <p>
     * Using the custom solver and integration defined above.
     */
    public static void main(final String[] args) {

        BasicLogger.debug();
        BasicLogger.debug(HookingUpSolver.class);
        BasicLogger.debug(OjAlgoUtils.getTitle());
        BasicLogger.debug(OjAlgoUtils.getDate());
        BasicLogger.debug();

        /*
         * Register the solver integration with ExpressionsBasedModel. This is a one-time operation – you only
         * need to do this once per JVM.
         */
        ExpressionsBasedModel.addIntegration(INTEGRATION);
        /*
         * Provided it's capable to solve the problem, ExpressionsBasedModel will now use the
         * UnconstrainedSolver rather than the default set of solvers. The integration's isCapable method
         * determines if the solver will be used.
         */

        ExpressionsBasedModel model = new ExpressionsBasedModel();

        /*
         * max: 2xy + 2x - x^2 - 2y^2
         */

        Variable x = model.addVariable("x");
        Variable y = model.addVariable("y");

        Expression expression = model.addExpression().weight(1);
        /*
         * Setting a weight on the expression is what turns it into (part of) the objective function. In a
         * more complex case the objective function could be a weighted combination of several expressions.
         * Here we only have one expression but still need to set a non-zero weight.
         */
        expression.set(x, y, 2); // +2xy
        expression.set(x, 2); // +2x
        expression.set(x, x, -1); // -x^2
        expression.set(y, y, -2); // -2y^2

        BasicLogger.debug("Expected result: {}", Optimisation.Result.of(2.0, Optimisation.State.OPTIMAL, 2.0, 1.0));
        BasicLogger.debug();

        /*
         * Turn on debug logging of the solver. With each iteration the current solution is printed as well as
         * a measure of the current error and the iteration number.
         */
        model.options.debug(UnconstrainedSolver.class);

        Optimisation.Result result = model.maximise();

        BasicLogger.debug();
        BasicLogger.debug("Actual result: {}", result);

        /*
         * ojAlgo's buil-in solvers are integrated with ExpressionsBasedModel the same way any other solver
         * is. The main built-in solvers are:
         */
        ExpressionsBasedModel.addIntegration(LinearSolver.INTEGRATION);
        ExpressionsBasedModel.addIntegration(ConvexSolver.INTEGRATION);
        ExpressionsBasedModel.addIntegration(IntegerSolver.INTEGRATION);
        /*
         * To learn how more complex integrations are built, study the source code of those.
         */
        /*
         * The only difference between the built-in solvers and custom ones is that you never have worry about
         * adding/registering built-in solver integrations. Adding them explicitly, like we just did above,
         * works fine but is entirely redundant. You only need to add/register solvers you want to use rather
         * than the built-in ones.
         */

    }

}