Example usage for org.apache.commons.math.optimization.linear SimplexSolver SimplexSolver

List of usage examples for org.apache.commons.math.optimization.linear SimplexSolver SimplexSolver

Introduction

In this page you can find the example usage for org.apache.commons.math.optimization.linear SimplexSolver SimplexSolver.

Prototype

public SimplexSolver(final double epsilon) 

Source Link

Document

Build a simplex solver with a specified accepted amount of error

Usage

From source file:net.sf.tweety.math.opt.solver.ApacheCommonsSimplex.java

@Override
public Map<Variable, Term> solve(ConstraintSatisfactionProblem problem) {
    if (!problem.isLinear())
        throw new IllegalArgumentException("Simplex algorithm is for linear problems only.");
    //this.log.info("Wrapping optimization problem for calling the Apache Commons Simplex algorithm.");            
    // 1.) bring all constraints in linear and normalized form
    Set<Statement> constraints = new HashSet<Statement>();
    for (Statement s : problem)
        constraints.add(s.toNormalizedForm().toLinearForm());
    // 2.) for every constraint we need an extra variable
    int numVariables = problem.getVariables().size();
    // 3.) define mappings from variables to indices
    int index = 0;
    Map<Variable, Integer> origVars2Idx = new HashMap<Variable, Integer>();
    for (Variable v : problem.getVariables())
        origVars2Idx.put(v, index++);/*from  w  w  w.jav  a 2  s  . c o  m*/
    // 4.) Check for target function (for constraint satisfaction problems
    //      its empty
    double[] coefficientsTarget = new double[numVariables];
    int i = 0;
    for (; i < numVariables; i++)
        coefficientsTarget[i] = 0;
    double constTerm = 0;
    if (problem instanceof OptimizationProblem) {
        // bring target function in linear form
        Sum t = ((OptimizationProblem) problem).getTargetFunction().toLinearForm();
        for (Term summand : t.getTerms()) {
            // as t is in linear form every summand is a product
            Product p = (Product) summand;
            if (p.getTerms().size() == 1) {
                // p consists of just a constant term
                Constant c = (Constant) p.getTerms().get(0);
                if (c instanceof FloatConstant)
                    constTerm += ((FloatConstant) c).getValue();
                else
                    constTerm += ((IntegerConstant) c).getValue();
            } else {
                // p consists of a variable and a constant
                Variable v = (Variable) ((p.getTerms().get(0) instanceof Variable) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                Constant c = (Constant) ((p.getTerms().get(0) instanceof Constant) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                double coefficient = (c instanceof FloatConstant) ? (((FloatConstant) c).getValue())
                        : (((IntegerConstant) c).getValue());
                coefficientsTarget[origVars2Idx.get(v)] += coefficient;
            }
        }
    }
    LinearObjectiveFunction target = new LinearObjectiveFunction(coefficientsTarget, constTerm);
    // 5.) Represent the constraints
    Set<LinearConstraint> finalConstraints = new HashSet<LinearConstraint>();
    for (Statement s : constraints) {
        double[] coefficientsConstraint = new double[numVariables];
        for (i = 0; i < numVariables; i++)
            coefficientsConstraint[i] = 0;
        // as s is in linear form go through the summands
        Sum leftTerm = (Sum) s.getLeftTerm();
        double rest = 0;
        for (Term summand : leftTerm.getTerms()) {
            // as s is in linear form every summand is a product
            Product p = (Product) summand;
            if (p.getTerms().size() == 1) {
                // p consists of just a constant term
                Constant c = (Constant) p.getTerms().get(0);
                if (c instanceof FloatConstant)
                    rest += ((FloatConstant) c).getValue();
                else
                    rest += ((IntegerConstant) c).getValue();
            } else {
                // p consists of a variable and a constant
                Variable v = (Variable) ((p.getTerms().get(0) instanceof Variable) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                Constant c = (Constant) ((p.getTerms().get(0) instanceof Constant) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                double coefficient = (c instanceof FloatConstant) ? (((FloatConstant) c).getValue())
                        : (((IntegerConstant) c).getValue());
                coefficientsConstraint[origVars2Idx.get(v)] += coefficient;
            }
        }
        Relationship r = Relationship.EQ;
        if (s instanceof Inequation)
            r = Relationship.GEQ;
        finalConstraints.add(new LinearConstraint(coefficientsConstraint, r, -rest));
    }
    // 6.) Optimize.
    try {
        //this.log.info("Calling the Apache Commons Simplex algorithm.");
        SimplexSolver solver = new SimplexSolver(0.01);
        solver.setMaxIterations(this.MAXITERATIONS);
        RealPointValuePair r = null;
        if (problem instanceof OptimizationProblem) {
            int type = ((OptimizationProblem) problem).getType();
            r = solver.optimize(target, finalConstraints,
                    (type == OptimizationProblem.MINIMIZE) ? (GoalType.MINIMIZE) : (GoalType.MAXIMIZE),
                    this.onlyPositive);
        } else
            r = solver.optimize(target, finalConstraints, GoalType.MINIMIZE, this.onlyPositive);
        //this.log.info("Parsing output from the Apache Commons Simplex algorithm.");
        Map<Variable, Term> result = new HashMap<Variable, Term>();
        for (Variable v : origVars2Idx.keySet())
            result.put(v, new FloatConstant(r.getPoint()[origVars2Idx.get(v)]));
        return result;
    } catch (OptimizationException e) {
        //log.error(e.getMessage());
        throw new ProblemInconsistentException();
    }
}

From source file:rb.app.RBnSCMCSystem.java

public double get_SCM_LB() {
    //return 0.01;

    double min_J_obj = 0.;
    double[] min_Jlocal_obj = new double[n_mubar];

    // Sort the indices of mu_bar based on distance from current_parameters
    List<Integer> sortedmubarIndices = getSorted_CJ_Indices(mu_bar);
    int icount = 0;
    //while ((min_J_obj<=0) && (icount < sortedmubarIndices.size())){
    while ((min_J_obj <= 0) && (icount < sortedmubarIndices.size())) {
        int imubar = sortedmubarIndices.get(icount);

        // First, declare the constraints
        Collection constraints = new ArrayList();

        // Add bounding box constraints for the get_Q_a() variables
        for (int q = 0; q < get_Q_a(); q++) {
            double[] index = new double[get_Q_a() * 2];
            index[q] = 1.;/*  w  w  w .j  a v a 2 s  .co m*/

            constraints.add(new LinearConstraint(index, Relationship.GEQ, B_min[q] / beta_bar[imubar]));
            constraints.add(new LinearConstraint(index, Relationship.LEQ, B_max[q] / beta_bar[imubar]));

            index[q] = 0.;
            index[q + get_Q_a()] = 1.;

            constraints.add(new LinearConstraint(index, Relationship.GEQ, B_min[q] / beta_bar[imubar]));
            constraints.add(new LinearConstraint(index, Relationship.LEQ, B_max[q] / beta_bar[imubar]));
        }

        // Save the current_parameters since we'll change them in the loop below
        save_current_parameters();

        // Add the constraint rows
        if (n_muhat[imubar] > 0) {
            for (int imuhat = 0; imuhat < n_muhat[imubar]; imuhat++) {
                current_parameters = mu_hat[imubar].get(imuhat);

                double[] constraint_row = new double[get_Q_a() * 2];
                for (int q = 0; q < get_Q_a(); q++) {
                    Complex theta_q_a = complex_eval_theta_q_a(q);
                    constraint_row[q] = theta_q_a.getReal() * beta_bar[imubar];
                    constraint_row[q + get_Q_a()] = theta_q_a.getImaginary() * beta_bar[imubar];
                }

                constraints
                        .add(new LinearConstraint(constraint_row, Relationship.GEQ, beta_hat[imubar][imuhat]));
            }
        }

        // Now load the original parameters back into current_parameters
        // in order to set the coefficients of the objective function
        reload_current_parameters();

        // Create objective function object
        double[] objectiveFn = new double[get_Q_a() * 2];
        for (int q = 0; q < get_Q_a(); q++) {
            Complex theta_q_a = complex_eval_theta_q_a(q);
            objectiveFn[q] = theta_q_a.getReal() * beta_bar[imubar];
            objectiveFn[q + get_Q_a()] = theta_q_a.getImaginary() * beta_bar[imubar];
        }
        LinearObjectiveFunction f = new LinearObjectiveFunction(objectiveFn, 0.);

        try {
            SimplexSolver solver = new SimplexSolver(1e-6);
            RealPointValuePair opt_pair = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
            min_Jlocal_obj[icount] = opt_pair.getValue();
        } catch (OptimizationException e) {
            Log.e("RBSCMSYSTEM_TAG", "Optimal solution not found");
            e.printStackTrace();
        } catch (Exception e) {
            Log.e("RBSCMSYSTEM_TAG", "Exception occurred during SCM_LB calculation");
            e.printStackTrace();
        }

        min_J_obj = min_J_obj > min_Jlocal_obj[icount] ? min_J_obj : min_Jlocal_obj[icount];
        icount++;
    }
    return min_J_obj;
}