StringVector theta,
Environment rho,
DoubleVector dir) {
if(dir.length() != theta.length()) {
throw new EvalException("'dir' is not a numeric vector of the correct length");
}
SEXP responseExpr = context.evaluate(modelExpr, rho);
if(!(responseExpr instanceof Vector)) {
throw new EvalException("Expected numeric response from model");
}
Vector response = (Vector)responseExpr;
for(int i=0;i!=response.length();++i) {
if(!DoubleVector.isFinite(response.getElementAsDouble(i))) {
throw new EvalException("Missing value or an infinity produced when evaluating the model");
}
}
double parameterValues[][] = new double[theta.length()][];
int totalParameterValueCount = 0;
for (int i = 0; i < theta.length(); i++) {
String parameterName = theta.getElementAsString(i);
SEXP parameterValue = rho.findVariable(Symbol.get(parameterName));
if (!(parameterValue instanceof AtomicVector)) {
throw new EvalException("variable '%s' is not numeric", parameterName);
}
parameterValues[i] = ((AtomicVector) parameterValue).toDoubleArray();
totalParameterValueCount += parameterValues[i].length;
}
double eps = Math.sqrt(DoubleVector.EPSILON);
DoubleMatrixBuilder gradient = new DoubleMatrixBuilder(response.length(), totalParameterValueCount);
int gradientColumn = 0;
for (int parameterIndex = 0; parameterIndex < theta.length(); parameterIndex++) {
double direction = dir.getElementAsDouble(parameterIndex);
// Parameters can be multivariate, not just single values, so loop over each
// element in this parameter
for (int parameterValueIndex = 0; parameterValueIndex < parameterValues[parameterIndex].length;
parameterValueIndex++) {
// update this individual parameter value
double startingParameterValue = parameterValues[parameterIndex][parameterValueIndex];
double absoluteParameterValue = Math.abs(startingParameterValue);
double delta = (absoluteParameterValue == 0) ? eps : absoluteParameterValue * eps;
parameterValues[parameterIndex][parameterValueIndex] +=
direction * delta;
// update this parameter in the model's environment
rho.setVariable(theta.getElementAsString(parameterIndex), new DoubleArrayVector(parameterValues[parameterIndex]));
// compute the new response given this updated parameter
DoubleVector responseDelta = (DoubleVector) context.evaluate(modelExpr, rho);
for (int k = 0; k < response.length(); k++) {
if (!DoubleVector.isFinite(responseDelta.getElementAsDouble(k))) {
throw new EvalException("Missing value or an infinity produced when evaluating the model");
}
double difference = responseDelta.getElementAsDouble(k) - response.getElementAsDouble(k);
double relativeDifference = difference / delta;
gradient.set(k, gradientColumn, direction * relativeDifference);
}