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

/**
 * Converts a test string to a {@link LinearConstraint}.
 * Ex: x0 + x1 + x2 + x3 - x12 = 0
 */
private LinearConstraint equationFromString(int numCoefficients, String s) {
    Relationship relationship;
    if (s.contains(">=")) {
        relationship = Relationship.GEQ;
    } else if (s.contains("<=")) {
        relationship = Relationship.LEQ;
    } else if (s.contains("=")) {
        relationship = Relationship.EQ;
    } else {
        throw new IllegalArgumentException();
    }

    String[] equationParts = s.split("[>|<]?=");
    double rhs = Double.parseDouble(equationParts[1].trim());

    RealVector lhs = new ArrayRealVector(numCoefficients);
    String left = equationParts[0].replaceAll(" ?x", "");
    String[] coefficients = left.split(" ");
    for (String coefficient : coefficients) {
        double value = coefficient.charAt(0) == '-' ? -1 : 1;
        int index = Integer.parseInt(coefficient.replaceFirst("[+|-]", "").trim());
        lhs.setEntry(index, value);
    }
    return new LinearConstraint(lhs, relationship, rhs);
}
 

/**
 * Predict the internal state estimation one time step ahead.
 *
 * @param u
 *            the control vector
 * @throws DimensionMismatchException
 *             if the dimension of the control vector does not fit
 */
public void predict(final RealVector u) {
    // sanity checks
    if (u != null &&
        u.getDimension() != controlMatrix.getColumnDimension()) {
        throw new DimensionMismatchException(u.getDimension(),
                                             controlMatrix.getColumnDimension());
    }

    // project the state estimation ahead (a priori state)
    // xHat(k)- = A * xHat(k-1) + B * u(k-1)
    stateEstimation = transitionMatrix.operate(stateEstimation);

    // add control input if it is available
    if (u != null) {
        stateEstimation = stateEstimation.add(controlMatrix.operate(u));
    }

    // project the error covariance ahead
    // P(k)- = A * P(k-1) * A' + Q
    errorCovariance = transitionMatrix.multiply(errorCovariance)
            .multiply(transitionMatrixT)
            .add(processModel.getProcessNoise());
}
 

/**
 * Verifies that setting X, Y and covariance separately has the same effect as newSample(X,Y,cov).
 */
@Test
public void testNewSample2() throws Exception {
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    double[][] covariance = MatrixUtils.createRealIdentityMatrix(4).scalarMultiply(2).getData();
    GLSMultipleLinearRegression regression = new GLSMultipleLinearRegression();
    regression.newSampleData(y, x, covariance);
    RealMatrix combinedX = regression.X.copy();
    RealVector combinedY = regression.Y.copy();
    RealMatrix combinedCovInv = regression.getOmegaInverse();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    assertEquals(combinedX, regression.X);
    assertEquals(combinedY, regression.Y);
    assertEquals(combinedCovInv, regression.getOmegaInverse());
}
 

/**
 * Verifies that setting X, Y and covariance separately has the same effect as newSample(X,Y,cov).
 */
@Test
public void testNewSample2() throws Exception {
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    double[][] covariance = MatrixUtils.createRealIdentityMatrix(4).scalarMultiply(2).getData();
    GLSMultipleLinearRegression regression = new GLSMultipleLinearRegression();
    regression.newSampleData(y, x, covariance);
    RealMatrix combinedX = regression.X.copy();
    RealVector combinedY = regression.Y.copy();
    RealMatrix combinedCovInv = regression.getOmegaInverse();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(combinedX, regression.X);
    Assert.assertEquals(combinedY, regression.Y);
    Assert.assertEquals(combinedCovInv, regression.getOmegaInverse());
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * @param xval the arguments for the interpolation points.
 * {@code xval[i][0]} is the first component of interpolation point
 * {@code i}, {@code xval[i][1]} is the second component, and so on
 * until {@code xval[i][d-1]}, the last component of that interpolation
 * point (where {@code dimension} is thus the dimension of the sampled
 * space).
 * @param yval the values for the interpolation points
 * @param brightnessExponent Brightness dimming factor.
 * @param microsphereElements Number of surface elements of the
 * microsphere.
 * @param rand Unit vector generator for creating the microsphere.
 * @throws DimensionMismatchException if the lengths of {@code yval} and
 * {@code xval} (equal to {@code n}, the number of interpolation points)
 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
 * have lengths different from {@code dimension}.
 * @throws NoDataException if there are no data (xval null or zero length)
 */
public MicrosphereInterpolatingFunction(double[][] xval,
                                        double[] yval,
                                        int brightnessExponent,
                                        int microsphereElements,
                                        UnitSphereRandomVectorGenerator rand)
    throws DimensionMismatchException, NoDataException {
    if (xval.length == 0 || xval[0] == null) {
        throw new NoDataException();
    }

    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }

    dimension = xval[0].length;
    this.brightnessExponent = brightnessExponent;

    // Copy data samples.
    samples = new HashMap<RealVector, Double>(yval.length);
    for (int i = 0; i < xval.length; ++i) {
        final double[] xvalI = xval[i];
        if ( xvalI.length != dimension) {
            throw new DimensionMismatchException(xvalI.length, dimension);
        }

        samples.put(new ArrayRealVector(xvalI), yval[i]);
    }

    microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
    // Generate the microsphere, assuming that a fairly large number of
    // randomly generated normals will represent a sphere.
    for (int i = 0; i < microsphereElements; i++) {
        microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
    }

}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoefficientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

@Test
public void testL1UDFSparse() throws Exception
{

  setMemorySettings();
  RandomGenerator rg = new JDKRandomGenerator();
  rg.setSeed(0);
  RandomData rd = new RandomDataImpl(rg);
  int n = 1000;
  List<RealVector> vectors = LSHTest.getVectors(rd, 1000, n);
  PigTest test = createPigTestFromString(l1SparseTest);
  writeLinesToFile("input", getSparseLines(vectors));
  List<RealVector> queries = LSHTest.getVectors(rd, 1000, 10);
  writeLinesToFile("queries", getSparseLines(queries));
  test.runScript();
  List<Tuple> neighbors = this.getLinesForAlias(test, "NEIGHBOR_CNT");
  Assert.assertEquals( queries.size(), neighbors.size() );
  for(long cnt : getCounts(neighbors))
  {
    Assert.assertTrue(cnt >= 3);
  }
  Distance d = new Distance()
  {

    @Override
    public double distance(RealVector v1, RealVector v2) {
      return L1.distance(v1, v2);
    }

  };
  verifyPoints(neighbors, d, 1000);
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Store the illumination and index of the brightest sample.
 * @param illuminationFromSample illumination received from sample
 * @param sample current sample illuminating the element
 */
void store(final double illuminationFromSample,
           final Map.Entry<RealVector, Double> sample) {
    if (illuminationFromSample > this.brightestIllumination) {
        this.brightestIllumination = illuminationFromSample;
        this.brightestSample = sample;
    }
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/**
 * Store the illumination and index of the brightest sample.
 * @param illuminationFromSample illumination received from sample
 * @param sample current sample illuminating the element
 */
void store(final double illuminationFromSample,
           final Map.Entry<RealVector, Double> sample) {
    if (illuminationFromSample > this.brightestIllumination) {
        this.brightestIllumination = illuminationFromSample;
        this.brightestSample = sample;
    }
}
 

/**
 * Store the illumination and index of the brightest sample.
 * @param illuminationFromSample illumination received from sample
 * @param sample current sample illuminating the element
 */
void store(final double illuminationFromSample,
           final Map.Entry<RealVector, Double> sample) {
    if (illuminationFromSample > this.brightestIllumination) {
        this.brightestIllumination = illuminationFromSample;
        this.brightestSample = sample;
    }
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * @param xval the arguments for the interpolation points.
 * {@code xval[i][0]} is the first component of interpolation point
 * {@code i}, {@code xval[i][1]} is the second component, and so on
 * until {@code xval[i][d-1]}, the last component of that interpolation
 * point (where {@code dimension} is thus the dimension of the sampled
 * space).
 * @param yval the values for the interpolation points
 * @param brightnessExponent Brightness dimming factor.
 * @param microsphereElements Number of surface elements of the
 * microsphere.
 * @param rand Unit vector generator for creating the microsphere.
 * @throws DimensionMismatchException if the lengths of {@code yval} and
 * {@code xval} (equal to {@code n}, the number of interpolation points)
 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
 * have lengths different from {@code dimension}.
 * @throws IllegalArgumentException if there are no data (xval null or zero length)
 */
public MicrosphereInterpolatingFunction(double[][] xval,
                                        double[] yval,
                                        int brightnessExponent,
                                        int microsphereElements,
                                        UnitSphereRandomVectorGenerator rand)
    throws DimensionMismatchException, IllegalArgumentException {
    if (xval.length == 0 || xval[0] == null) {
        throw MathRuntimeException.createIllegalArgumentException("no data");
    }

    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }

    dimension = xval[0].length;
    this.brightnessExponent = brightnessExponent;

    // Copy data samples.
    samples = new HashMap<RealVector, Double>(yval.length);
    for (int i = 0; i < xval.length; ++i) {
        final double[] xvalI = xval[i];
        if ( xvalI.length != dimension) {
            throw new DimensionMismatchException(xvalI.length, dimension);
        }

        samples.put(new ArrayRealVector(xvalI), yval[i]);
    }

    microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
    // Generate the microsphere, assuming that a fairly large number of
    // randomly generated normals will represent a sphere.
    for (int i = 0; i < microsphereElements; i++) {
        microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
    }

}
 

/**
 * Compute the LSH for a given vector.
 */
public long apply(RealVector vector)
{
  /*
   * The hash is just floor(<v, a>/w)
   */
   double ret = b;
  
   for(int i = 0;i < dim;++i)
   {
      ret += vector.getEntry(i)*a[i];
   }
   return (long)Math.floor(ret/w);
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * @param xval the arguments for the interpolation points.
 * {@code xval[i][0]} is the first component of interpolation point
 * {@code i}, {@code xval[i][1]} is the second component, and so on
 * until {@code xval[i][d-1]}, the last component of that interpolation
 * point (where {@code dimension} is thus the dimension of the sampled
 * space).
 * @param yval the values for the interpolation points
 * @param brightnessExponent Brightness dimming factor.
 * @param microsphereElements Number of surface elements of the
 * microsphere.
 * @param rand Unit vector generator for creating the microsphere.
 * @throws DimensionMismatchException if the lengths of {@code yval} and
 * {@code xval} (equal to {@code n}, the number of interpolation points)
 * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]},
 * have lengths different from {@code dimension}.
 * @throws IllegalArgumentException if there are no data (xval null or zero length)
 */
public MicrosphereInterpolatingFunction(double[][] xval,
                                        double[] yval,
                                        int brightnessExponent,
                                        int microsphereElements,
                                        UnitSphereRandomVectorGenerator rand)
    throws DimensionMismatchException, IllegalArgumentException {
    if (xval.length == 0 || xval[0] == null) {
        throw MathRuntimeException.createIllegalArgumentException("no data");
    }

    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }

    dimension = xval[0].length;
    this.brightnessExponent = brightnessExponent;

    // Copy data samples.
    samples = new HashMap<RealVector, Double>(yval.length);
    for (int i = 0; i < xval.length; ++i) {
        final double[] xvalI = xval[i];
        if ( xvalI.length != dimension) {
            throw new DimensionMismatchException(xvalI.length, dimension);
        }

        samples.put(new ArrayRealVector(xvalI), yval[i]);
    }

    microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements);
    // Generate the microsphere, assuming that a fairly large number of
    // randomly generated normals will represent a sphere.
    for (int i = 0; i < microsphereElements; i++) {
        microsphere.add(new MicrosphereSurfaceElement(rand.nextVector()));
    }

}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Verifies that newSampleData methods consistently insert unitary columns
 * in design matrix.  Confirms the fix for MATH-411.
 */
@Test
public void testNewSample() throws Exception {
    double[] design = new double[] {
      1, 19, 22, 33,
      2, 20, 30, 40,
      3, 25, 35, 45,
      4, 27, 37, 47
    };
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    AbstractMultipleLinearRegression regression = createRegression();
    regression.newSampleData(design, 4, 3);
    RealMatrix flatX = regression.X.copy();
    RealVector flatY = regression.Y.copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    assertEquals(flatX, regression.X);
    assertEquals(flatY, regression.Y);
    
    // No intercept
    regression.setNoIntercept(true);
    regression.newSampleData(design, 4, 3);
    flatX = regression.X.copy();
    flatY = regression.Y.copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    assertEquals(flatX, regression.X);
    assertEquals(flatY, regression.Y);
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.getData()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * {@inheritDoc}
 */
public double[] estimateRegressionParameters() {
    RealVector b = calculateBeta();
    return b.getData();
}
 

/**
 * 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 not 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
 *
 */
@Override
public void estimate(EstimationProblem problem)
throws EstimationException {

    initializeEstimate(problem);

    // work matrices
    double[] grad             = new double[parameters.length];
    ArrayRealVector bDecrement = new ArrayRealVector(parameters.length);
    double[] bDecrementData   = bDecrement.getDataRef();
    RealMatrix wGradGradT     = MatrixUtils.createRealMatrix(parameters.length, parameters.length);

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

        // build the linear problem
        incrementJacobianEvaluationsCounter();
        RealVector b = new ArrayRealVector(parameters.length);
        RealMatrix a = MatrixUtils.createRealMatrix(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] = weight * residual * grad[j];
                }

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

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

            }
        }

        try {

            // solve the linearized least squares problem
            RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b);

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

        } catch(InvalidMatrixException e) {
            throw new EstimationException("unable to solve: singular problem");
        }


        previous = cost;
        updateResidualsAndCost();

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

}
 

/**
 * Create the tableau by itself.
 * @param maximize if true, goal is to maximize the objective function
 * @return created tableau
 */
protected RealMatrix createTableau(final boolean maximize) {

    // create a matrix of the correct size
    int width = numDecisionVariables + numSlackVariables +
    numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
    int height = constraints.size() + getNumObjectiveFunctions();
    Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);

    // initialize the objective function rows
    if (getNumObjectiveFunctions() == 2) {
        matrix.setEntry(0, 0, -1);
    }
    int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
    matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);
    RealVector objectiveCoefficients =
        maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
    copyArray(objectiveCoefficients.getData(), matrix.getDataRef()[zIndex]);
    matrix.setEntry(zIndex, width - 1,
        maximize ? f.getConstantTerm() : -1 * f.getConstantTerm());

    if (!restrictToNonNegative) {
        matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
            getInvertedCoeffiecientSum(objectiveCoefficients));
    }

    // initialize the constraint rows
    int slackVar = 0;
    int artificialVar = 0;
    for (int i = 0; i < constraints.size(); i++) {
        LinearConstraint constraint = constraints.get(i);
        int row = getNumObjectiveFunctions() + i;

        // decision variable coefficients
        copyArray(constraint.getCoefficients().getData(), matrix.getDataRef()[row]);

        // x-
        if (!restrictToNonNegative) {
            matrix.setEntry(row, getSlackVariableOffset() - 1,
                getInvertedCoeffiecientSum(constraint.getCoefficients()));
        }

        // RHS
        matrix.setEntry(row, width - 1, constraint.getValue());

        // slack variables
        if (constraint.getRelationship() == Relationship.LEQ) {
            matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);  // slack
        } else if (constraint.getRelationship() == Relationship.GEQ) {
            matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess
        }

        // artificial variables
        if ((constraint.getRelationship() == Relationship.EQ) ||
                (constraint.getRelationship() == Relationship.GEQ)) {
            matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
            matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
            matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
        }
    }

    return matrix;
}
 

/**
 * {@inheritDoc}
 */
public double[] estimateRegressionParameters() {
    RealVector b = calculateBeta();
    return b.getData();
}
 

/**
 * Create the tableau by itself.
 * @param maximize if true, goal is to maximize the objective function
 * @return created tableau
 */
protected RealMatrix createTableau(final boolean maximize) {

    // create a matrix of the correct size
    int width = numDecisionVariables + numSlackVariables +
    numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
    int height = constraints.size() + getNumObjectiveFunctions();
    Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);

    // initialize the objective function rows
    if (getNumObjectiveFunctions() == 2) {
        matrix.setEntry(0, 0, -1);
    }
    int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
    matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);
    RealVector objectiveCoefficients =
        maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
    copyArray(objectiveCoefficients.toArray(), matrix.getDataRef()[zIndex]);
    matrix.setEntry(zIndex, width - 1,
        maximize ? f.getConstantTerm() : -1 * f.getConstantTerm());

    if (!restrictToNonNegative) {
        matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
            getInvertedCoefficientSum(objectiveCoefficients));
    }

    // initialize the constraint rows
    int slackVar = 0;
    int artificialVar = 0;
    for (int i = 0; i < constraints.size(); i++) {
        LinearConstraint constraint = constraints.get(i);
        int row = getNumObjectiveFunctions() + i;

        // decision variable coefficients
        copyArray(constraint.getCoefficients().toArray(), matrix.getDataRef()[row]);

        // x-
        if (!restrictToNonNegative) {
            matrix.setEntry(row, getSlackVariableOffset() - 1,
                getInvertedCoefficientSum(constraint.getCoefficients()));
        }

        // RHS
        matrix.setEntry(row, width - 1, constraint.getValue());

        // slack variables
        if (constraint.getRelationship() == Relationship.LEQ) {
            matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);  // slack
        } else if (constraint.getRelationship() == Relationship.GEQ) {
            matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess
        }

        // artificial variables
        if ((constraint.getRelationship() == Relationship.EQ) ||
                (constraint.getRelationship() == Relationship.GEQ)) {
            matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
            matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
            matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
        }
    }

    return matrix;
}