LineOptSolver.m

classdef LineOptSolver < handle
    % LineOptSolver  Main line-opt solver. Minimizes a penalized scalar
    % objective (constraints as penalties) over the decision variables,
    % aggregating across scenarios, using the self-contained
    % opt.de.DifferentialEvolution engine (numpy-exact RNG) or an analytic/FD
    % projected-gradient path. Mirrors native-Python LineOptSolver.
 
    properties
        problem
        opt
        fixedValueMap
        evaluators        % cell of opt.LineEvaluator
        scenarioWeights
        caches            % cell of containers.Map
        iterations = 0;
        bestValue = inf;
        bestX = [];
        convergenceHistory = [];
        startTime
        deadline = inf;
    end
 
    methods
        function obj = LineOptSolver(problem, options)
            obj.problem = problem;
            obj.opt = options;
            obj.fixedValueMap = containers.Map('KeyType', 'char', 'ValueType', 'any');
            fixed = problem.getFixedVariables();
            for k = 1:numel(fixed)
                obj.fixedValueMap(fixed{k}{1}.getName()) = fixed{k}{2};
            end
            obj.evaluators = {opt.LineEvaluator(problem.getModel(), problem.getVariables(), fixed)};
            scen = problem.getScenarios();
            w = zeros(1, 1 + numel(scen));
            w(1) = 1.0;
            for i = 1:numel(scen)
                obj.evaluators{end+1} = opt.LineEvaluator(scen{i}{1}, problem.getVariables(), fixed);
                w(i+1) = scen{i}{2};
            end
            obj.scenarioWeights = w;
        end
 
        function result = solve(obj)
            obj.startTime = tic;
            obj.iterations = 0;
            obj.bestValue = inf;
            obj.bestX = [];
            obj.convergenceHistory = [];
            obj.caches = cell(1, numel(obj.evaluators));
            for i = 1:numel(obj.evaluators)
                obj.caches{i} = containers.Map('KeyType', 'char', 'ValueType', 'any');
            end
            obj.deadline = obj.opt.timeLimit;
 
            bounds = obj.evaluators{1}.getBounds();
            if isempty(bounds)
                result = obj.buildEmptyResult();
                return;
            end
            if obj.shouldUseGradient()
                result = obj.solveGradient(bounds);
            else
                result = obj.solveEvolution(bounds);
            end
        end
 
        function result = solveEvolution(obj, bounds)
            low = bounds(:, 1).';
            high = bounds(:, 2).';
            if isempty(obj.opt.seed)
                seed = mod(int64(feature('timing','cpucount')), 2^31);
            else
                seed = obj.opt.seed;
            end
            de = opt.de.DifferentialEvolution(@(x) obj.objectiveFunction(x), low, high, ...
                obj.opt.strategy, obj.opt.popsize, obj.opt.maxIterations, ...
                obj.opt.mutationLow, obj.opt.mutationHigh, obj.opt.recombination, obj.opt.tol, seed);
            de.callback = @(bx, nit) obj.deCallback(nit);
 
            timedOut = false;
            try
                r = de.solve();
                resX = obj.engineBestVector(r);
                resFun = r.fun;
            catch ME
                if strcmp(ME.identifier, 'LineOpt:TimeLimit')
                    timedOut = true;
                    if isempty(obj.bestX)
                        result = obj.buildEmptyResult(); return;
                    end
                    resX = obj.bestX; resFun = obj.bestValue;
                else
                    rethrow(ME);
                end
            end
            solveTime = toc(obj.startTime);
            result = obj.buildResult(resX, resFun, solveTime);
            if timedOut, result.terminatedBy = 'time_limit'; end
        end
 
        function stop = deCallback(obj, nit)
            obj.iterations = nit;
            obj.convergenceHistory(end+1) = obj.bestValue;
            stop = toc(obj.startTime) >= obj.deadline;
        end
 
        function x = engineBestVector(obj, r)
            if ~isempty(obj.bestX) && obj.bestValue <= r.fun
                x = obj.bestX;
            else
                x = r.x;
            end
        end
 
        function total = objectiveFunction(obj, x)
            if toc(obj.startTime) >= obj.deadline
                error('LineOpt:TimeLimit', 'time limit reached');
            end
            values = obj.evaluators{1}.decodeVariables(x);
            allValues = obj.mergeValues(values);
            objective = obj.problem.getObjective();
            pw = obj.opt.penaltyWeight;
            cons = obj.problem.getConstraints();
            scenarioValues = zeros(1, numel(obj.evaluators));
            for i = 1:numel(obj.evaluators)
                res = obj.evaluators{i}.evaluateValuesWithCache(values, obj.caches{i});
                if ~res.feasible
                    total = inf; return;
                end
                value = objective.evaluateWithPenalty(res, allValues, pw);
                for c = 1:numel(cons)
                    value = value + cons{c}.evaluate(res, allValues) * pw;
                end
                scenarioValues(i) = value;
            end
            total = obj.aggregateScenarios(scenarioValues);
            if total < obj.bestValue
                obj.bestValue = total;
                obj.bestX = x;
            end
        end
 
        function av = mergeValues(obj, values)
            % merge: start from fixed, overlay values
            av = containers.Map('KeyType', 'char', 'ValueType', 'any');
            fk = keys(obj.fixedValueMap);
            for i = 1:numel(fk), av(fk{i}) = obj.fixedValueMap(fk{i}); end
            vk = keys(values);
            for i = 1:numel(vk), av(vk{i}) = values(vk{i}); end
        end
 
        function total = aggregateScenarios(obj, values)
            if numel(values) == 1
                total = values(1); return;
            end
            if strcmp(obj.opt.scenarioAggregation, 'mean')
                total = sum(obj.scenarioWeights .* values) / sum(obj.scenarioWeights);
            else
                total = max(values);
            end
        end
 
        function tf = shouldUseGradient(obj)
            if strcmp(obj.opt.optimizer, 'gradient'), tf = true; return; end
            if strcmp(obj.opt.optimizer, 'evolution'), tf = false; return; end
            tf = obj.allContinuous();
        end
 
        function tf = allContinuous(obj)
            vars = obj.problem.getVariables();
            if isempty(vars), tf = false; return; end
            tf = true;
            for i = 1:numel(vars)
                vt = vars{i}.getVariableType();
                if ~any(strcmp(vt, {'service_rate', 'routing'}))
                    tf = false; return;
                end
            end
        end
 
        function result = solveGradient(obj, bounds)
            dim = size(bounds, 1);
            nStart = max(1, obj.opt.gradientRestarts);
            if isempty(obj.opt.seed), sd = 0; else, sd = obj.opt.seed; end
            rs = RandStream('mt19937ar', 'Seed', sd);
            starts = {0.5 * ones(1, dim)};
            for s = 2:nStart, starts{end+1} = rand(rs, 1, dim); end %#ok<AGROW>
            best = []; bestFun = inf; timedOut = false;
            for si = 1:numel(starts)
                if toc(obj.startTime) >= obj.deadline, timedOut = true; break; end
                try
                    xx = obj.projectedGradientDescent(starts{si});
                    fx = obj.objectiveFunction(xx);
                    if isfinite(fx) && fx < bestFun, bestFun = fx; best = xx; end
                catch ME
                    if strcmp(ME.identifier, 'LineOpt:TimeLimit'), timedOut = true; break; else, rethrow(ME); end
                end
            end
            if isempty(best) || (~isempty(obj.bestX) && obj.bestValue < bestFun)
                if ~isempty(obj.bestX), best = obj.bestX; bestFun = obj.bestValue;
                elseif isempty(best), result = obj.buildEmptyResult(); return; end
            end
            solveTime = toc(obj.startTime);
            result = obj.buildResult(best, bestFun, solveTime);
            if timedOut, result.terminatedBy = 'time_limit'; end
        end
 
        function x = projectedGradientDescent(obj, x0)
            x = min(max(x0, 0), 1);
            f = obj.objectiveFunction(x);
            for iter = 1:obj.opt.maxIterations
                g = obj.objectiveGradient(x);
                if norm(g) < 1e-9, break; end
                step = 1.0; improved = false;
                for ls = 1:30
                    xn = min(max(x - step * g, 0), 1);
                    fn = obj.objectiveFunction(xn);
                    if isfinite(fn) && fn < f - 1e-12
                        x = xn; f = fn; improved = true; break;
                    end
                    step = step * 0.5;
                end
                if ~improved, break; end
            end
        end
 
        function g = objectiveGradient(obj, x)
            h = obj.opt.fdStep;
            dim = numel(x);
            g = zeros(1, dim);
            f0 = [];
            for i = 1:dim
                xp = x; xm = x;
                xp(i) = min(1.0, x(i) + h);
                xm(i) = max(0.0, x(i) - h);
                fp = obj.objectiveFunction(xp);
                fm = obj.objectiveFunction(xm);
                if isfinite(fp) && isfinite(fm) && xp(i) > xm(i)
                    g(i) = (fp - fm) / (xp(i) - xm(i)); continue;
                end
                if isempty(f0), f0 = obj.objectiveFunction(x); end
                if isfinite(fp) && isfinite(f0) && xp(i) > x(i)
                    g(i) = (fp - f0) / (xp(i) - x(i));
                elseif isfinite(fm) && isfinite(f0) && x(i) > xm(i)
                    g(i) = (f0 - fm) / (x(i) - xm(i));
                else
                    g(i) = 0.0;
                end
            end
        end
 
        function result = buildResult(obj, x, objectiveValue, solveTime)
            result = opt.OptimizationResult();
            result.objectiveValue = objectiveValue;
            vals = obj.evaluators{1}.decodeVariables(x);
            vk = keys(vals);
            for i = 1:numel(vk), result.variableValues(vk{i}) = vals(vk{i}); end
            result.iterations = obj.iterations;
            result.solveTime = solveTime;
            evals = 0;
            for i = 1:numel(obj.evaluators), evals = evals + obj.evaluators{i}.getEvaluationCount(); end
            result.modelEvaluations = evals;
            result.convergenceHistory = obj.convergenceHistory;
 
            allValues = obj.mergeValues(result.variableValues);
            objective = obj.problem.getObjective();
            allCons = [objective.getConstraints(), obj.problem.getConstraints()];
            result.feasible = true;
            for i = 1:numel(obj.evaluators)
                er = obj.evaluators{i}.evaluateValuesWithCache(result.variableValues, obj.caches{i});
                if ~er.feasible, result.feasible = false; continue; end
                for c = 1:numel(allCons)
                    viol = allCons{c}.evaluate(er, allValues);
                    if viol > 0
                        nm = allCons{c}.getName();
                        prev = 0.0; if isKey(result.constraintViolations, nm), prev = result.constraintViolations(nm); end
                        result.constraintViolations(nm) = max(viol, prev);
                        result.feasible = false;
                    end
                end
            end
            elapsed = toc(obj.startTime);
            if elapsed >= obj.opt.timeLimit
                result.terminatedBy = 'time_limit';
            elseif obj.iterations >= obj.opt.maxIterations
                result.terminatedBy = 'iterations';
            else
                result.terminatedBy = 'convergence';
            end
        end
 
        function result = buildEmptyResult(obj) %#ok<MANU>
            result = opt.OptimizationResult();
            result.objectiveValue = 0.0;
            result.feasible = true;
            result.terminatedBy = 'empty';
        end
    end
end