org.apache.commons.math.linear.Array2DRowRealMatrix 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 Array2DRowRealMatrix(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

@Test
public void testLeastSquares1()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 3.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 4.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * <p>Compute the "hat" matrix.
 * </p>
 * <p>The hat matrix is defined in terms of the design matrix X
 *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
 * </p>
 * <p>The implementation here uses the QR decomposition to compute the
 * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
 * p-dimensional identity matrix augmented by 0's.  This computational
 * formula is from "The Hat Matrix in Regression and ANOVA",
 * David C. Hoaglin and Roy E. Welsch,
 * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
 *
 * @return the hat matrix
 */
public RealMatrix calculateHat() {
    // Create augmented identity matrix
    RealMatrix Q = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = Q.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j =0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }

    // Compute and return Hat matrix
    return Q.multiply(augI).multiply(Q.transpose());
}
 

@Test
public void testLeastSquares1()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 3.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 4.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * 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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * 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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * 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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * <p>Compute the "hat" matrix.
 * </p>
 * <p>The hat matrix is defined in terms of the design matrix X
 *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
 * </p>
 * <p>The implementation here uses the QR decomposition to compute the
 * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
 * p-dimensional identity matrix augmented by 0's.  This computational
 * formula is from "The Hat Matrix in Regression and ANOVA",
 * David C. Hoaglin and Roy E. Welsch,
 * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
 *
 * @return the hat matrix
 */
public RealMatrix calculateHat() {
    // Create augmented identity matrix
    RealMatrix Q = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = Q.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j =0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }

    // Compute and return Hat matrix
    return Q.multiply(augI).multiply(Q.transpose());
}
 

/**
 * 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 Array2DRowRealMatrix(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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

@Test
public void testLeastSquares1()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 3.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 4.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

@Test
public void testLeastSquares3() {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
                { 1, 0 },
                { 0, 1 }
            }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
            public double[] value(double[] variables) {
                return factors.operate(variables);
            }
        }, new double[] { 2, -3 }, new Array2DRowRealMatrix(new double [][] {
                { 1, 1.2 }, { 1.2, 2 }
            }));
    SimplexOptimizer optimizer = new SimplexOptimizer(-1, 1e-6);
    optimizer.setSimplex(new NelderMeadSimplex(2));
    RealPointValuePair optimum =
        optimizer.optimize(200, ls, GoalType.MINIMIZE, new double[] { 10, 10 });
    Assert.assertEquals( 2, optimum.getPointRef()[0], 2e-3);
    Assert.assertEquals(-3, optimum.getPointRef()[1], 8e-4);
    Assert.assertTrue(optimizer.getEvaluations() > 60);
    Assert.assertTrue(optimizer.getEvaluations() < 80);
    Assert.assertTrue(optimum.getValue() < 1e-6);
}
 

/**
 * 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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

@Test
public void testLeastSquares2()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 }, new double[] { 10.0, 0.1 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 5.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 8.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * 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 Array2DRowRealMatrix(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 Array2DRowRealMatrix(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

@Test
public void testLeastSquares1()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 3.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 4.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * 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 Array2DRowRealMatrix(createTableau(goalType == GoalType.MAXIMIZE));
    initialize();
}
 

@Test
public void testLeastSquares1()
    throws FunctionEvaluationException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
                { 1.0, 0.0 },
                { 0.0, 1.0 }
            }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
            public double[] value(double[] variables) {
                return factors.operate(variables);
            }
        }, new double[] { 2.0, -3.0 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxEvaluations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 3.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 4.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * 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 Array2DRowRealMatrix(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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

/**
 * 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 Array2DRowRealMatrix(matrix);
    this.numArtificialVariables = 0;
}
 

@Test
public void testLeastSquares2()
throws FunctionEvaluationException, ConvergenceException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
            { 1.0, 0.0 },
            { 0.0, 1.0 }
        }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
        public double[] value(double[] variables) {
            return factors.operate(variables);
        }
    }, new double[] { 2.0, -3.0 }, new double[] { 10.0, 0.1 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxIterations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 5.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 8.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/** Reinitialize the instance.
 * <p>Beware that all arrays <em>must</em> be references to integrator
 * arrays, in order to ensure proper update without copy.</p>
 * @param time time at which all arrays are defined
 * @param stepSize step size used in the scaled and nordsieck arrays
 * @param scaledDerivative reference to the integrator array holding the first
 * scaled derivative
 * @param nordsieckVector reference to the integrator matrix holding the
 * nordsieck vector
 */
public void reinitialize(final double time, final double stepSize,
                         final double[] scaledDerivative,
                         final Array2DRowRealMatrix nordsieckVector) {
    this.referenceTime = time;
    this.scalingH      = stepSize;
    this.scaled        = scaledDerivative;
    this.nordsieck     = nordsieckVector;

    // make sure the state and derivatives will depend on the new arrays
    setInterpolatedTime(getInterpolatedTime());

}
 

/**
 * Loads model x and y sample data from a flat array of data, overriding any previous sample.
 * Assumes that rows are concatenated with y values first in each row.
 *
 * @param data input data array
 * @param nobs number of observations (rows)
 * @param nvars number of independent variables (columns, not counting y)
 */
public void newSampleData(double[] data, int nobs, int nvars) {
    double[] y = new double[nobs];
    double[][] x = new double[nobs][nvars + 1];
    int pointer = 0;
    for (int i = 0; i < nobs; i++) {
        y[i] = data[pointer++];
        x[i][0] = 1.0d;
        for (int j = 1; j < nvars + 1; j++) {
            x[i][j] = data[pointer++];
        }
    }
    this.X = new Array2DRowRealMatrix(x);
    this.Y = new ArrayRealVector(y);
}
 

/** Initialize the high order scaled derivatives at step start.
 * @param first first scaled derivative at step start
 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
 * will be modified
 * @return high order derivatives at step start
 */
public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
                                                 final double[][] multistep) {
    for (int i = 0; i < multistep.length; ++i) {
        final double[] msI = multistep[i];
        for (int j = 0; j < first.length; ++j) {
            msI[j] -= first[j];
        }
    }
    return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
}
 

/**
 * @param m
 *            Input matrix.
 * @return Row matrix representing the sums of the rows.
 */
private static RealMatrix sumRows(final RealMatrix m) {
    double[][] d = new double[1][m.getColumnDimension()];
    for (int c = 0; c < m.getColumnDimension(); c++) {
        double sum = 0;
        for (int r = 0; r < m.getRowDimension(); r++) {
            sum += m.getEntry(r, c);
        }
        d[0][c] = sum;
    }
    return new Array2DRowRealMatrix(d, false);
}
 

/** {@inheritDoc} */
@Override
public void readExternal(final ObjectInput in)
    throws IOException, ClassNotFoundException {

    // read the base class 
    final double t = readBaseExternal(in);

    // read the local attributes
    scalingH      = in.readDouble();
    referenceTime = in.readDouble();

    final int n = (currentState == null) ? -1 : currentState.length;
    final boolean hasScaled = in.readBoolean();
    if (hasScaled) {
        scaled = new double[n];
        for (int j = 0; j < n; ++j) {
            scaled[j] = in.readDouble();
        }
    } else {
        scaled = null;
    }

    final boolean hasNordsieck = in.readBoolean();
    if (hasNordsieck) {
        nordsieck = (Array2DRowRealMatrix) in.readObject();
    } else {
        nordsieck = null;
    }

    if (hasScaled && hasNordsieck) {
        // we can now set the interpolated time and state
        stateVariation = new double[n];
        setInterpolatedTime(t);
    } else {
        stateVariation = null;
    }

}
 

/**
 * @param n
 *            Number of rows.
 * @param m
 *            Number of columns.
 * @return n X m matrix of 1.0-values.
 */
private static RealMatrix ones(int n, int m) {
    double[][] d = new double[n][m];
    for (int r = 0; r < n; r++) {
        Arrays.fill(d[r], 1.0);
    }
    return new Array2DRowRealMatrix(d, false);
}
 

/**
 * @param n
 *            Number of rows.
 * @param m
 *            Number of columns.
 * @return n X m matrix of 0.0-values, diagonal has values 1.0.
 */
private static RealMatrix eye(int n, int m) {
    double[][] d = new double[n][m];
    for (int r = 0; r < n; r++)
        if (r < m)
            d[r][r] = 1;
    return new Array2DRowRealMatrix(d, false);
}