use of cbit.vcell.opt.OptimizationResultSet in project vcell by virtualcell.
the class CurveFitting method fitRecovery_diffOnly.
/*
* @para: time, time points since the first post bleach.
* @para: flour, average intensities under bleached region according to time points since the first post bleach.
* @para: bleachType, gaussian spot of half cell bleaching. model expressions vary based on bleaching type.
* @para: inputparam, bleaching while monitoring rate, will be substituted in the model expression.
* @para: outputParam, the array which will pass results back.
*/
public static Expression fitRecovery_diffOnly(double[] time, double[] normalized_fluor, int bleachType, double[] inputparam, double[] outputParam) throws ExpressionException, OptimizationException, IOException {
if (time.length != normalized_fluor.length) {
throw new RuntimeException("Fluorecence and time arrays must be the same length");
}
// normaliztion for time by subtracting the starting time: time[0]
double[] normalized_time = new double[time.length];
for (int i = 0; i < time.length; i++) {
normalized_time[i] = time[i] - time[0];
}
Expression modelExp = null;
OptimizationResultSet optResultSet = null;
OptSolverResultSet optSolverResultSet = null;
String[] paramNames = null;
double[] paramValues = null;
if (bleachType == FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_GaussianSpot || bleachType == FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_HalfCell) {
if (bleachType == FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_GaussianSpot) {
Expression muExp = new Expression(FRAPDataAnalysis.gaussianSpot_MuFunc);
modelExp = new Expression(FRAPDataAnalysis.gaussianSpot_IntensityFunc);
modelExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_u), muExp);
} else {
Expression muExp = new Expression(FRAPDataAnalysis.halfCell_MuFunc);
modelExp = new Expression(FRAPDataAnalysis.halfCell_IntensityFunc);
modelExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_u), muExp);
}
// inputparam[0] is the bleach while monitoring rate.
double bleachWhileMonitoringRate = inputparam[0];
modelExp = modelExp.getSubstitutedExpression(new Expression(FRAPDataAnalysis.symbol_bwmRate), new Expression(bleachWhileMonitoringRate));
Parameter[] parameters = new Parameter[] { FRAPDataAnalysis.para_If, FRAPDataAnalysis.para_Io, FRAPDataAnalysis.para_tau };
bindExpressionToParametersAndTime(modelExp, parameters);
// estimate parameters by minimizing the errors between function values and reference data
optResultSet = solve(modelExp.flatten(), parameters, normalized_time, normalized_fluor);
optSolverResultSet = optResultSet.getOptSolverResultSet();
paramNames = optSolverResultSet.getParameterNames();
paramValues = optSolverResultSet.getBestEstimates();
// copy into "output" buffer from parameter values.
for (int i = 0; i < paramValues.length; i++) {
outputParam[i] = paramValues[i];
}
} else {
throw new IllegalArgumentException("Unknown bleach type " + bleachType);
}
for (int i = 0; i < paramNames.length; i++) {
System.out.println("finally: " + paramNames[i] + " = " + paramValues[i]);
}
// investigate the return information
processReturnCode(OptimizationSolverSpec.SOLVERTYPE_POWELL, optSolverResultSet);
// construct recovery under bleached region by diffusion only expression
Expression fit = new Expression(modelExp);
System.out.println("fit before subsituting parameters:" + fit.infix());
// substitute parameter values
for (int i = 0; i < paramValues.length; i++) {
fit.substituteInPlace(new Expression(paramNames[i]), new Expression(paramValues[i]));
}
// undo time shift
fit.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
// undo fluorescence normalization
System.out.println("fit equation after unnorm:" + fit.infix());
return fit;
}
use of cbit.vcell.opt.OptimizationResultSet in project vcell by virtualcell.
the class CurveFitting method fitBleachWhileMonitoring.
/*
* @para: time, time points since the first post bleach
* @para: flour, average intensities under bleached region according to time points since the first post bleach
* @para: parameterValues, the array which will pass results back
*/
public static Expression fitBleachWhileMonitoring(double[] time, double[] normalized_fluor, double[] outputParam, Weights weights) throws ExpressionException, OptimizationException, IOException {
if (time.length != normalized_fluor.length) {
throw new RuntimeException("Fluorecence and time arrays must be the same length");
}
// normaliztion for time by subtracting the starting time: time[0]
double[] normalized_time = new double[time.length];
for (int i = 0; i < time.length; i++) {
normalized_time[i] = time[i] - time[0];
}
Expression modelExp = null;
OptimizationResultSet optResultSet = null;
String[] paramNames = null;
double[] paramValues = null;
double cellFirstPostBleach = normalized_fluor[0];
modelExp = new Expression(FRAPOptFunctions.FUNC_CELL_INTENSITY);
modelExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(cellFirstPostBleach));
// initialize starting guess, arguments in Parameter are name, Lower Bound, Upper Bound, Scale, Initial Guess
Parameter[] parameters = new Parameter[] { new Parameter(FRAPOptFunctions.SYMBOL_BWM_RATE, FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getLowerBound(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getUpperBound(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getScale(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getInitialGuess()) };
bindExpressionToParametersAndTime(modelExp, parameters);
// estimate blech while monitoring rate by minimizing the error between funtion values and reference data
if (weights == null) {
optResultSet = solve(modelExp, parameters, normalized_time, normalized_fluor);
} else {
optResultSet = solve(modelExp, parameters, normalized_time, normalized_fluor, weights);
}
OptSolverResultSet optSolverResultSet = optResultSet.getOptSolverResultSet();
paramNames = optSolverResultSet.getParameterNames();
paramValues = optSolverResultSet.getBestEstimates();
// copy into "output" buffer for parameter values.
for (int i = 0; i < paramValues.length; i++) {
outputParam[i] = paramValues[i];
}
System.out.println(optSolverResultSet.getOptimizationStatus().toString());
for (int i = 0; i < paramNames.length; i++) {
System.out.println("finally: " + paramNames[i] + " = " + paramValues[i]);
}
// investigate the return information
processReturnCode(OptimizationSolverSpec.SOLVERTYPE_POWELL, optSolverResultSet);
//
// construct final equation
//
Expression fit = new Expression(modelExp);
System.out.println("bleach while monitoring fit before subsituting parameters:" + fit.infix());
//
for (int i = 0; i < paramValues.length; i++) {
fit.substituteInPlace(new Expression(paramNames[i]), new Expression(paramValues[i]));
}
//
// undo time shift
//
fit.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
//
// undo fluorescence normalization
//
System.out.println("bleach while monitoring fit equation after unnorm:" + fit.infix());
return fit;
}
use of cbit.vcell.opt.OptimizationResultSet in project vcell by virtualcell.
the class CurveFitting method fitRecovery_reacKoffRateOnly.
/*
* @para: time, time points since the first post bleach.
* @para: normalized_data, first dimension index 0:average intensities under cell region, first dimension index 1: average intensities under bleached region
* @para: inputparam, index 0:cellROIAvg at time 0(first post bleach). index 1:bleachedROIAvg at time 0(first post bleach).
* @para: outputParam, results for 2 or 3 parameters, depending on if there is any fixed parameter.
*/
public static double fitRecovery_reacKoffRateOnly(double[] time, double[][] normalized_data, double[] inputparam, double[] outputParam, Parameter fixedParameter, Weights weights) throws ExpressionException, OptimizationException, IOException {
if (normalized_data != null && normalized_data.length > 0) {
for (int i = 0; i < normalized_data.length; i++) {
if (normalized_data[i] != null && time.length != normalized_data[i].length) {
throw new RuntimeException("Fluorecence and time arrays must be the same length");
}
}
}
// normaliztion for time by subtracting the starting time: time[0]
double[] normalized_time = new double[time.length];
for (int i = 0; i < time.length; i++) {
normalized_time[i] = time[i] - time[0];
}
// initiate variables
OptimizationResultSet optResultSet = null;
Expression bwmRateExp = new Expression(FRAPOptFunctions.FUNC_CELL_INTENSITY);
Expression koffRateExp = new Expression(FRAPOptFunctions.FUNC_RECOVERY_BLEACH_REACTION_DOMINANT);
// get parameters to be estimated(use all of the 3 if there is no fixed parameter. Or 2 out of the 3 if there is a fixed parameter)
Parameter bwmParam = new Parameter(FRAPOptFunctions.SYMBOL_BWM_RATE, FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getLowerBound(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getUpperBound(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getScale(), FRAPModel.REF_BLEACH_WHILE_MONITOR_PARAM.getInitialGuess());
Parameter koffParam = new Parameter(FRAPOptFunctions.SYMBOL_KOFF, FRAPModel.REF_REACTION_OFF_RATE.getLowerBound(), FRAPModel.REF_REACTION_OFF_RATE.getUpperBound(), FRAPModel.REF_REACTION_OFF_RATE.getScale(), FRAPModel.REF_REACTION_OFF_RATE.getInitialGuess());
Parameter fittingParamA = new Parameter(FRAPOptFunctions.SYMBOL_A, /*binding site concentration is reused to store fitting parameter A, but the name can not be reused*/
FRAPModel.REF_BS_CONCENTRATION_OR_A.getLowerBound(), FRAPModel.REF_BS_CONCENTRATION_OR_A.getUpperBound(), FRAPModel.REF_BS_CONCENTRATION_OR_A.getScale(), FRAPModel.REF_BS_CONCENTRATION_OR_A.getInitialGuess());
// get column names for reference data to be constructed
String[] columnNames = new String[] { ReservedVariable.TIME.getName(), "cellIntensityAvg", "bleachIntensityAvg" };
if (fixedParameter == null) {
// get fitting parameter array
Parameter[] parameters = new Parameter[] { bwmParam, fittingParamA, koffParam };
// get expression pairs
bwmRateExp = bwmRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(inputparam[0]));
bindExpressionToParametersAndTime(bwmRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair bwmRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(bwmRateExp, 1);
koffRateExp = koffRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(inputparam[1]));
bindExpressionToParametersAndTime(koffRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair koffRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(koffRateExp, 2);
ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = new ExplicitFitObjectiveFunction.ExpressionDataPair[] { bwmRateExpDataPair, koffRateExpDataPair };
// solve
optResultSet = solve(expDataPairs, parameters, normalized_time, normalized_data, columnNames, weights);
} else {
if (fixedParameter != null && fixedParameter.getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_BLEACH_MONITOR_RATE])) {
// get fitting parameter array
Parameter[] parameters = new Parameter[] { fittingParamA, koffParam };
// get expression pairs
bwmRateExp = bwmRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(inputparam[0]));
bwmRateExp = bwmRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(fixedParameter.getInitialGuess()));
bindExpressionToParametersAndTime(bwmRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair bwmRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(bwmRateExp, 1);
koffRateExp = koffRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(inputparam[1]));
koffRateExp = koffRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(fixedParameter.getInitialGuess()));
bindExpressionToParametersAndTime(koffRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair koffRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(koffRateExp, 2);
ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = new ExplicitFitObjectiveFunction.ExpressionDataPair[] { bwmRateExpDataPair, koffRateExpDataPair };
// solve
optResultSet = solve(expDataPairs, parameters, normalized_time, normalized_data, columnNames, weights);
} else if (fixedParameter != null && fixedParameter.getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_OFF_RATE])) {
// get fitting parameter array
Parameter[] parameters = new Parameter[] { bwmParam, fittingParamA };
// get expression pairs
bwmRateExp = bwmRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(inputparam[0]));
bindExpressionToParametersAndTime(bwmRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair bwmRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(bwmRateExp, 1);
koffRateExp = koffRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(inputparam[1]));
koffRateExp = koffRateExp.getSubstitutedExpression(new Expression(FRAPOptFunctions.SYMBOL_KOFF), new Expression(fixedParameter.getInitialGuess()));
bindExpressionToParametersAndTime(koffRateExp, parameters);
ExplicitFitObjectiveFunction.ExpressionDataPair koffRateExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(koffRateExp, 2);
ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = new ExplicitFitObjectiveFunction.ExpressionDataPair[] { bwmRateExpDataPair, koffRateExpDataPair };
// solve
optResultSet = solve(expDataPairs, parameters, normalized_time, normalized_data, columnNames, weights);
}
}
// get output parameter values
OptSolverResultSet optSolverResultSet = optResultSet.getOptSolverResultSet();
String[] paramNames = optSolverResultSet.getParameterNames();
double[] paramValues = optSolverResultSet.getBestEstimates();
// copy into "output" buffer from parameter values.
for (int i = 0; i < paramValues.length; i++) {
outputParam[i] = paramValues[i];
}
// for debug purpose
for (int i = 0; i < paramNames.length; i++) {
System.out.println("finally: " + paramNames[i] + " = " + paramValues[i]);
}
// investigate the return information
processReturnCode(OptimizationSolverSpec.SOLVERTYPE_POWELL, optSolverResultSet);
// return objective function value
return optSolverResultSet.getLeastObjectiveFunctionValue();
}
use of cbit.vcell.opt.OptimizationResultSet in project vcell by virtualcell.
the class CurveFitting method solveSpatial.
public static OptimizationResultSet solveSpatial(Simulation sim, Parameter[] parameters, SpatialReferenceData refData, File dataDir, FieldDataIdentifierSpec[] fieldDataIDSs) throws ExpressionException, OptimizationException, IOException {
// choose optimization solver, currently we have Powell
PowellOptimizationSolver optService = new PowellOptimizationSolver();
OptimizationSpec optSpec = new OptimizationSpec();
// send to optimization service
optSpec.setObjectiveFunction(new PdeObjectiveFunction(sim.getMathDescription(), refData, dataDir, fieldDataIDSs));
double[] parameterValues = new double[parameters.length];
for (int i = 0; i < parameters.length; i++) {
parameterValues[i] = parameters[i].getInitialGuess();
System.out.println("initial " + parameters[i].getName() + " = " + parameterValues[i] + ";\tLB = " + parameters[i].getLowerBound() + ";\tUB = " + parameters[i].getUpperBound());
}
// get the initial guess to send it to the f() function. ....
for (int i = 0; i < parameters.length; i++) {
optSpec.addParameter(parameters[i]);
}
// Parameters in OptimizationSolverSpec are solver type and objective function change tolerance.
OptimizationSolverSpec optSolverSpec = new OptimizationSolverSpec(OptimizationSolverSpec.SOLVERTYPE_POWELL, 0.000001);
OptSolverCallbacks optSolverCallbacks = new DefaultOptSolverCallbacks();
OptimizationResultSet optResultSet = null;
optResultSet = optService.solve(optSpec, optSolverSpec, optSolverCallbacks);
OptSolverResultSet optSolverResultSet = optResultSet.getOptSolverResultSet();
String[] paramNames = optSolverResultSet.getParameterNames();
double[] paramValues = optSolverResultSet.getBestEstimates();
for (int i = 0; i < paramNames.length; i++) {
System.out.println("finally: " + paramNames[i] + " = " + paramValues[i]);
}
return optResultSet;
}
use of cbit.vcell.opt.OptimizationResultSet in project vcell by virtualcell.
the class OptimizationService method optimize.
public static OptimizationResultSet optimize(ParameterEstimationTask parameterEstimationTask) throws Exception {
// if (OperatingSystemInfo.getInstance().isMac()){
// throw new RuntimeException("parameter estimation not currently available on Mac\n\n try Windows or Linux.\n\n coming soon on Mac.");
// }
copasiOptCallbacks.reset();
updateMath(parameterEstimationTask.getSimulationContext(), NetworkGenerationRequirements.ComputeFullStandardTimeout);
StringBuffer issueText = new StringBuffer();
java.util.Vector<Issue> issueList = new java.util.Vector<Issue>();
IssueContext issueContext = new IssueContext();
parameterEstimationTask.gatherIssues(issueContext, issueList);
boolean bFailed = false;
for (int i = 0; i < issueList.size(); i++) {
Issue issue = (Issue) issueList.elementAt(i);
issueText.append(issue.getMessage() + "\n");
if (issue.getSeverity() == Issue.SEVERITY_ERROR) {
bFailed = true;
break;
}
}
if (bFailed) {
throw new OptimizationException(issueText.toString());
}
parameterEstimationTask.refreshMappings();
OptimizationResultSet optResultSet = CopasiOptimizationSolver.solveLocalPython(new ParameterEstimationTaskSimulatorIDA(), parameterEstimationTask, copasiOptCallbacks, callback);
return optResultSet;
}
Aggregations