Validate.notNull(startPos, "startPos");
final int nObs = observedValues.getNumberOfElements();
Validate.isTrue(nObs == sigma.getNumberOfElements(), "observedValues and sigma must be same length");
ArgumentChecker.isTrue(allowedValue.evaluate(startPos), "The start position {} is not valid for this model. Please choose a valid start position", startPos);
DoubleMatrix2D alpha;
DecompositionResult decmp;
DoubleMatrix1D theta = startPos;
double lambda = 0.0; // TODO debug if the model is linear, it will be solved in 1 step
double newChiSqr, oldChiSqr;
DoubleMatrix1D error = getError(func, observedValues, sigma, theta);
DoubleMatrix1D newError;
DoubleMatrix2D jacobian = getJacobian(jac, sigma, theta);
oldChiSqr = getChiSqr(error);
double p = getANorm(penalty, theta);
oldChiSqr += p;
// If we start at the solution we are done
if (oldChiSqr == 0.0) {
return finish(oldChiSqr, jacobian, theta, sigma);
}
DoubleMatrix1D beta = getChiSqrGrad(error, jacobian);
DoubleMatrix1D temp = (DoubleMatrix1D) _algebra.multiply(penalty, theta);
beta = (DoubleMatrix1D) _algebra.subtract(beta, temp);
for (int count = 0; count < MAX_ATTEMPTS; count++) {
alpha = getModifiedCurvatureMatrix(jacobian, lambda, penalty);
DoubleMatrix1D deltaTheta;
try {
decmp = _decomposition.evaluate(alpha);
deltaTheta = decmp.solve(beta);
} catch (final Exception e) {
throw new MathException(e);
}
DoubleMatrix1D trialTheta = (DoubleMatrix1D) _algebra.add(theta, deltaTheta);
// if the new value of theta is not in the model domain keep increasing lambda until an acceptable step is found
if (!allowedValue.evaluate(trialTheta)) {
lambda = increaseLambda(lambda);
continue;
}
newError = getError(func, observedValues, sigma, trialTheta);
p = getANorm(penalty, trialTheta);
newChiSqr = getChiSqr(newError);
newChiSqr += p;
// Check for convergence when no improvement in chiSqr occurs
if (Math.abs(newChiSqr - oldChiSqr) / (1 + oldChiSqr) < _eps) {
final DoubleMatrix2D alpha0 = lambda == 0.0 ? alpha : getModifiedCurvatureMatrix(jacobian, 0.0, penalty);
if (lambda > 0.0) {
decmp = _decomposition.evaluate(alpha0);
}
return finish(alpha0, decmp, newChiSqr, jacobian, trialTheta, sigma);