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

/** {@inheritDoc} */
@Override
protected RealPointValuePair doOptimize() {
    final double[] lowerBound = getLowerBound();
    final double[] upperBound = getUpperBound();

    // Validity checks.
    setup(lowerBound, upperBound);

    isMinimize = (getGoalType() == GoalType.MINIMIZE);
    currentBest = new ArrayRealVector(getStartPoint());

    final double value = bobyqa(lowerBound, upperBound);

    return new RealPointValuePair(currentBest.getDataRef(),
                                  isMinimize ? value : -value);
}
 

/**
 * 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);
}
 

public static double computeSimilarityScore(double[] vec1, double[] vec2,
    SimilarityMethodType simMethodType) throws IllegalArgumentException {

  double simScore = 0.0;

  try {

    if (simMethodType == SimilarityMethodType.DOT) {
      double[] vec1_norm = new double[vec1.length];
      double[] vec2_norm = new double[vec2.length];

      double div1 = 0.0, div2 = 0.0;
      for (int i = 0; i < vec1.length; ++i) {
        div1 += vec1[i] * vec1[i];
        div2 += vec2[i] * vec2[i];
      }
      for (int i = 0; i < vec1.length; ++i) {
        vec1_norm[i] = vec1[i] / Math.sqrt(div1);
        vec2_norm[i] = vec2[i] / Math.sqrt(div2);
      }
      simScore = (new ArrayRealVector(vec1_norm)).dotProduct(vec2_norm);
    } else if (simMethodType == SimilarityMethodType.PEARSON) {
      simScore = new PearsonsCorrelation().correlation(vec1, vec2);
    }
  } catch (IllegalArgumentException e) {
    throw new IllegalArgumentException("Failed to compute similarity score for vec1.length="
        + vec1.length + " and vec2.length=" + vec2.length);
  }

  return simScore;
}
 

/** test solve */
public void testSolve() {
    DecompositionSolver solver =
        new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1, 2, 3 }, { 0, -5, 1 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { -8.0 / 25.0, -263.0 / 75.0, -29.0 / 75.0 },
            { 19.0 / 25.0,   78.0 / 25.0,  49.0 / 25.0 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using Array2DRowRealMatrix
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealMatrix with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * 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);
}
 

public RepeatingLSH(List<LSH> lshList) throws MathException
{
  super(lshList.get(0).getDim(), lshList.get(0).getRandomGenerator());
  this.lshList = lshList;
  RandomGenerator rg = lshList.get(0).getRandomGenerator();
  RandomData rd = new RandomDataImpl(rg);
  /*
   * Compute a random vector of lshList.size() with each component taken from U(0,10)
   */
  randomVec = new ArrayRealVector(lshList.size());
  for(int i = 0; i < randomVec.getDimension();++i)
  {
    randomVec.setEntry(i, rd.nextUniform(0, 10.0));
  }
}
 

/**
 * Loads new y sample data, overriding any previous data.
 *
 * @param y the array representing the y sample
 * @throws IllegalArgumentException if y is null or empty
 */
protected void newYSampleData(double[] y) {
    if (y == null) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.NULL_NOT_ALLOWED);
    }
    if (y.length == 0) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.NO_DATA);
    }
    this.Y = new ArrayRealVector(y);
}
 

/** test solve */
public void testSolve() {
    DecompositionSolver solver =
        new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1, 2, 3 }, { 0, -5, 1 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { -8.0 / 25.0, -263.0 / 75.0, -29.0 / 75.0 },
            { 19.0 / 25.0,   78.0 / 25.0,  49.0 / 25.0 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using Array2DRowRealMatrix
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealMatrix with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * 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);
}
 

/** test solve */
public void testSolve() {
    DecompositionSolver solver =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1, 0 }, { 2, -5 }, { 3, 1 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 19, -71 }, { -6, 22 }, { -2, 9 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using ArrayRealVector
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealVector with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * 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);
}
 

/**
 * 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);
}
 

/** test solve */
public void testSolve() {
    DecompositionSolver solver =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1, 0 }, { 2, -5 }, { 3, 1 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 19, -71 }, { -6, 22 }, { -2, 9 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using ArrayRealVector
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealVector with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * 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);
}
 

/**
 * Loads new y sample data, overriding any previous data.
 *
 * @param y the array representing the y sample
 * @throws IllegalArgumentException if y is null or empty
 */
protected void newYSampleData(double[] y) {
    if (y == null) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.NULL_NOT_ALLOWED);
    }
    if (y.length == 0) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.NO_DATA);
    }
    this.Y = new ArrayRealVector(y);
}
 

/**
 * 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);
}
 

/**
 * 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);
}
 

/**
 * 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);
}
 

/**
 * @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()));
    }

}
 

/**
 * 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)));

}
 

/** test solve */
public void testSolve() throws MathException {
    DecompositionSolver solver =
        new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            {   78,  -13,    1 },
            {  414,  -62,   -1 },
            { 1312, -202,  -37 },
            { 2989, -542,  145 },
            { 5510, -1465, 201 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1,  0,  1 },
            { 0,  1,  1 },
            { 2,  1, -4 },
            { 2,  2,  2 },
            { 5, -3,  0 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using ArrayRealVector
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealVector with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * @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 an array has zero-length.
 * @throws NullArgumentException if an argument is {@code null}.
 */
public MicrosphereInterpolatingFunction(double[][] xval,
                                        double[] yval,
                                        int brightnessExponent,
                                        int microsphereElements,
                                        UnitSphereRandomVectorGenerator rand) {
    if (xval == null ||
        yval == null) {
        throw new NullArgumentException();
    }
    if (xval.length == 0) {
        throw new NoDataException();
    }
    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }
    if (xval[0] == null) {
        throw new NullArgumentException();
    }

    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 == null) {
            throw new NullArgumentException();
        }
        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()));
    }
}
 

/** test solve */
public void testSolve() {
    RealMatrix m = MatrixUtils.createRealMatrix(new double[][] {
            { 91,  5, 29, 32, 40, 14 },
            {  5, 34, -1,  0,  2, -1 },
            { 29, -1, 12,  9, 21,  8 },
            { 32,  0,  9, 14,  9,  0 },
            { 40,  2, 21,  9, 51, 19 },
            { 14, -1,  8,  0, 19, 14 }
    });
    DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1561, 269, 188 },
            {   69, -21,  70 },
            {  739, 108,  63 },
            {  324,  86,  59 },
            { 1624, 194, 107 },
            {  796,  69,  36 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1,   2, 1 },
            { 2,  -1, 2 },
            { 4,   2, 3 },
            { 8,  -1, 0 },
            { 16,  2, 0 },
            { 32, -1, 0 }
    });

    // using RealMatrix
    assertEquals(0, es.solve(b).subtract(xRef).getNorm(), 2.0e-12);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(es.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

    // using Array2DRowRealMatrix
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     es.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

    // using RealMatrix with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     es.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

}
 

/**
 * 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)));

}
 

/** 
 * 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)));

}
 

/**
 * @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 an array has zero-length.
 * @throws NullArgumentException if an argument is {@code null}.
 */
public MicrosphereInterpolatingFunction(double[][] xval,
                                        double[] yval,
                                        int brightnessExponent,
                                        int microsphereElements,
                                        UnitSphereRandomVectorGenerator rand) {
    if (xval == null ||
        yval == null) {
        throw new NullArgumentException();
    }
    if (xval.length == 0) {
        throw new NoDataException();
    }
    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }
    if (xval[0] == null) {
        throw new NullArgumentException();
    }

    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 == null) {
            throw new NullArgumentException();
        }
        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()));
    }
}
 

public FortranArray(ArrayRealVector data) {
    super(data, false);
}
 

/** test solve */
public void testSolve() throws MathException {
    DecompositionSolver solver =
        new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            {   78,  -13,    1 },
            {  414,  -62,   -1 },
            { 1312, -202,  -37 },
            { 2989, -542,  145 },
            { 5510, -1465, 201 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1,  0,  1 },
            { 0,  1,  1 },
            { 2,  1, -4 },
            { 2,  2,  2 },
            { 5, -3,  0 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using ArrayRealVector
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealVector with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

@Test
public void testConstant() {
    double constantValue = 10d;
    double measurementNoise = 0.1d;
    double processNoise = 1e-5d;

    // A = [ 1 ]
    RealMatrix A = new Array2DRowRealMatrix(new double[] { 1d });
    // no control input
    RealMatrix B = null;
    // H = [ 1 ]
    RealMatrix H = new Array2DRowRealMatrix(new double[] { 1d });
    // x = [ 10 ]
    RealVector x = new ArrayRealVector(new double[] { constantValue });
    // Q = [ 1e-5 ]
    RealMatrix Q = new Array2DRowRealMatrix(new double[] { processNoise });
    // R = [ 0.1 ]
    RealMatrix R = new Array2DRowRealMatrix(new double[] { measurementNoise });

    ProcessModel pm
        = new DefaultProcessModel(A, B, Q,
                                  new ArrayRealVector(new double[] { constantValue }), null);
    MeasurementModel mm = new DefaultMeasurementModel(H, R);
    KalmanFilter filter = new KalmanFilter(pm, mm);

    Assert.assertEquals(1, filter.getMeasurementDimension());
    Assert.assertEquals(1, filter.getStateDimension());

    assertMatrixEquals(Q.getData(), filter.getErrorCovariance());

    // check the initial state
    double[] expectedInitialState = new double[] { constantValue };
    assertVectorEquals(expectedInitialState, filter.getStateEstimation());

    RealVector pNoise = new ArrayRealVector(1);
    RealVector mNoise = new ArrayRealVector(1);

    RandomGenerator rand = new JDKRandomGenerator();
    // iterate 60 steps
    for (int i = 0; i < 60; i++) {
        filter.predict();

        // Simulate the process
        pNoise.setEntry(0, processNoise * rand.nextGaussian());

        // x = A * x + p_noise
        x = A.operate(x).add(pNoise);

        // Simulate the measurement
        mNoise.setEntry(0, measurementNoise * rand.nextGaussian());

        // z = H * x + m_noise
        RealVector z = H.operate(x).add(mNoise);

        filter.correct(z);

        // state estimate should be larger than measurement noise
        double diff = Math.abs(constantValue - filter.getStateEstimation()[0]);
        // System.out.println(diff);
        Assert.assertTrue(MathUtils.compareTo(diff, measurementNoise, 1e-6) < 0);
    }

    // error covariance should be already very low (< 0.02)
    Assert.assertTrue(MathUtils.compareTo(filter.getErrorCovariance()[0][0],
                                          0.02d, 1e-6) < 0);
}
 

/** test solve */
public void testSolve() {
    RealMatrix m = MatrixUtils.createRealMatrix(new double[][] {
            { 91,  5, 29, 32, 40, 14 },
            {  5, 34, -1,  0,  2, -1 },
            { 29, -1, 12,  9, 21,  8 },
            { 32,  0,  9, 14,  9,  0 },
            { 40,  2, 21,  9, 51, 19 },
            { 14, -1,  8,  0, 19, 14 }
    });
    DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1561, 269, 188 },
            {   69, -21,  70 },
            {  739, 108,  63 },
            {  324,  86,  59 },
            { 1624, 194, 107 },
            {  796,  69,  36 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1,   2, 1 },
            { 2,  -1, 2 },
            { 4,   2, 3 },
            { 8,  -1, 0 },
            { 16,  2, 0 },
            { 32, -1, 0 }
    });

    // using RealMatrix
    assertEquals(0, es.solve(b).subtract(xRef).getNorm(), 2.0e-12);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(es.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

    // using Array2DRowRealMatrix
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     es.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

    // using RealMatrix with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     es.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     2.0e-11);
    }

}