Search in sources :

Example 6 with OptimizationResultSet

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;
}
Also used : Expression(cbit.vcell.parser.Expression) OptimizationResultSet(cbit.vcell.opt.OptimizationResultSet) Parameter(cbit.vcell.opt.Parameter) OptSolverResultSet(cbit.vcell.opt.OptSolverResultSet)

Example 7 with OptimizationResultSet

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;
}
Also used : Expression(cbit.vcell.parser.Expression) OptimizationResultSet(cbit.vcell.opt.OptimizationResultSet) Parameter(cbit.vcell.opt.Parameter) OptSolverResultSet(cbit.vcell.opt.OptSolverResultSet)

Example 8 with OptimizationResultSet

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();
}
Also used : OptimizationResultSet(cbit.vcell.opt.OptimizationResultSet) ExplicitFitObjectiveFunction(cbit.vcell.opt.ExplicitFitObjectiveFunction) OptSolverResultSet(cbit.vcell.opt.OptSolverResultSet) Expression(cbit.vcell.parser.Expression) Parameter(cbit.vcell.opt.Parameter)

Example 9 with OptimizationResultSet

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;
}
Also used : PowellOptimizationSolver(cbit.vcell.opt.solvers.PowellOptimizationSolver) DefaultOptSolverCallbacks(org.vcell.optimization.DefaultOptSolverCallbacks) OptSolverCallbacks(org.vcell.optimization.OptSolverCallbacks) OptimizationResultSet(cbit.vcell.opt.OptimizationResultSet) PdeObjectiveFunction(cbit.vcell.opt.PdeObjectiveFunction) DefaultOptSolverCallbacks(org.vcell.optimization.DefaultOptSolverCallbacks) OptimizationSolverSpec(cbit.vcell.opt.OptimizationSolverSpec) OptimizationSpec(cbit.vcell.opt.OptimizationSpec) OptSolverResultSet(cbit.vcell.opt.OptSolverResultSet)

Example 10 with OptimizationResultSet

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;
}
Also used : OptimizationException(cbit.vcell.opt.OptimizationException) Issue(org.vcell.util.Issue) OptimizationResultSet(cbit.vcell.opt.OptimizationResultSet) IssueContext(org.vcell.util.IssueContext) ParameterEstimationTaskSimulatorIDA(org.vcell.optimization.ParameterEstimationTaskSimulatorIDA)

Aggregations

OptimizationResultSet (cbit.vcell.opt.OptimizationResultSet)20 OptSolverResultSet (cbit.vcell.opt.OptSolverResultSet)12 OptimizationSolverSpec (cbit.vcell.opt.OptimizationSolverSpec)10 PowellOptimizationSolver (cbit.vcell.opt.solvers.PowellOptimizationSolver)8 OptimizationSpec (cbit.vcell.opt.OptimizationSpec)7 DefaultOptSolverCallbacks (org.vcell.optimization.DefaultOptSolverCallbacks)7 OptSolverCallbacks (org.vcell.optimization.OptSolverCallbacks)7 OptimizationException (cbit.vcell.opt.OptimizationException)6 Parameter (cbit.vcell.opt.Parameter)6 SimpleReferenceData (cbit.vcell.opt.SimpleReferenceData)5 ExplicitFitObjectiveFunction (cbit.vcell.opt.ExplicitFitObjectiveFunction)4 OptRunResultSet (cbit.vcell.opt.OptSolverResultSet.OptRunResultSet)4 Expression (cbit.vcell.parser.Expression)4 CopasiOptimizationMethod (cbit.vcell.opt.CopasiOptimizationMethod)3 OptimizationStatus (cbit.vcell.opt.OptimizationStatus)3 DefaultScalarFunction (cbit.function.DefaultScalarFunction)2 RowColumnResultSet (cbit.vcell.math.RowColumnResultSet)2 Parameter (cbit.vcell.model.Parameter)2 CopasiOptimizationParameter (cbit.vcell.opt.CopasiOptimizationParameter)2 ImplicitObjectiveFunction (cbit.vcell.opt.ImplicitObjectiveFunction)2