org.apache.commons.math.linear.RealMatrixImpl Java Examples

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Build a tableau for a linear problem.
 * @param f linear objective function
 * @param constraints linear constraints
 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 * or {@link GoalType#MINIMIZE}
 * @param restrictToNonNegative whether to restrict the variables to non-negative values
 * @param epsilon amount of error to accept in floating point comparisons
 */
SimplexTableau(final LinearObjectiveFunction f,
               final Collection<LinearConstraint> constraints,
               final GoalType goalType, final boolean restrictToNonNegative,
               final double epsilon) {
    this.f                      = f;
    this.constraints            = constraints;
    this.restrictToNonNegative  = restrictToNonNegative;
    this.epsilon                = epsilon;
    this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
    this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                  getConstraintTypeCounts(Relationship.GEQ);
    this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

/**
 * Get the covariance matrix.
 * @return covariance matrix
 */
public RealMatrix getResult() {

    int dimension = sums.length;
    RealMatrixImpl result = new RealMatrixImpl(dimension, dimension);

    if (n > 1) {
        double[][] resultData = result.getDataRef();
        double c = 1.0 / (n * (isBiasCorrected ? (n - 1) : n));
        int k = 0;
        for (int i = 0; i < dimension; ++i) {
            for (int j = 0; j <= i; ++j) {
                double e = c * (n * productsSums[k++] - sums[i] * sums[j]);
                resultData[i][j] = e;
                resultData[j][i] = e;
            }
        }
    }

    return result;

}
 

/**
 * Removes the phase 1 objective function and artificial variables from this tableau.
 */
protected void discardArtificialVariables() {
    if (numArtificialVariables == 0) {
        return;
    }
    int width = getWidth() - numArtificialVariables - 1;
    int height = getHeight() - 1;
    double[][] matrix = new double[height][width];
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width - 1; j++) {
            matrix[i][j] = getEntry(i + 1, j + 1);
        }
        matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
    }
    this.tableau = new RealMatrixImpl(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {
 
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int rows = problem.getMeasurements().length;
    final int cols = problem.getAllParameters().length;
    final int max  = cols * rows;
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += cols) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        return new RealMatrixImpl(jTj).inverse().getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem",
                                      new Object[0]);
    }

}
 

/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {
 
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int rows = problem.getMeasurements().length;
    final int cols = problem.getUnboundParameters().length;
    final int max  = cols * rows;
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += cols) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        return new RealMatrixImpl(jTj).inverse().getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem",
                                      new Object[0]);
    }

}
 

/**
 * Get the covariance matrix of estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {
 
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int rows = problem.getMeasurements().length;
    final int cols = problem.getAllParameters().length;
    final int max  = cols * rows;
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += cols) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        return new RealMatrixImpl(jTj).inverse().getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem",
                                      new Object[0]);
    }

}
 

/** 
 * Solve an estimation problem using a least squares criterion.
 *
 * <p>This method set the unbound parameters of the given problem
 * starting from their current values through several iterations. At
 * each step, the unbound parameters are changed in order to
 * minimize a weighted least square criterion based on the
 * measurements of the problem.</p>
 *
 * <p>The iterations are stopped either when the criterion goes
 * below a physical threshold under which improvement are considered
 * useless or when the algorithm is unable to improve it (even if it
 * is still high). The first condition that is met stops the
 * iterations. If the convergence it nos reached before the maximum
 * number of iterations, an {@link EstimationException} is
 * thrown.</p>
 *
 * @param problem estimation problem to solve
 * @exception EstimationException if the problem cannot be solved
 *
 * @see EstimationProblem
 *
 */
public void estimate(EstimationProblem problem)
throws EstimationException {

    initializeEstimate(problem);

    // work matrices
    double[] grad             = new double[parameters.length];
    RealMatrixImpl bDecrement = new RealMatrixImpl(parameters.length, 1);
    double[][] bDecrementData = bDecrement.getDataRef();
    RealMatrixImpl wGradGradT = new RealMatrixImpl(parameters.length, parameters.length);
    double[][] wggData        = wGradGradT.getDataRef();

    // iterate until convergence is reached
    double previous = Double.POSITIVE_INFINITY;
    do {

        // build the linear problem
        incrementJacobianEvaluationsCounter();
        RealMatrix b = new RealMatrixImpl(parameters.length, 1);
        RealMatrix a = new RealMatrixImpl(parameters.length, parameters.length);
        for (int i = 0; i < measurements.length; ++i) {
            if (! measurements [i].isIgnored()) {

                double weight   = measurements[i].getWeight();
                double residual = measurements[i].getResidual();

                // compute the normal equation
                for (int j = 0; j < parameters.length; ++j) {
                    grad[j] = measurements[i].getPartial(parameters[j]);
                    bDecrementData[j][0] = weight * residual * grad[j];
                }

                // build the contribution matrix for measurement i
                for (int k = 0; k < parameters.length; ++k) {
                    double[] wggRow = wggData[k];
                    double gk = grad[k];
                    for (int l = 0; l < parameters.length; ++l) {
                        wggRow[l] =  weight * gk * grad[l];
                    }
                }

                // update the matrices
                a = a.add(wGradGradT);
                b = b.add(bDecrement);

            }
        }

        try {

            // solve the linearized least squares problem
            RealMatrix dX = a.solve(b);

            // update the estimated parameters
            for (int i = 0; i < parameters.length; ++i) {
                parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i, 0));
            }

        } catch(InvalidMatrixException e) {
            throw new EstimationException("unable to solve: singular problem", new Object[0]);
        }


        previous = cost;
        updateResidualsAndCost();

    } while ((getCostEvaluations() < 2) ||
             (Math.abs(previous - cost) > (cost * steadyStateThreshold) &&
              (Math.abs(cost) > convergence)));

}