Examples with Parameter cbit.vcell.opt.Parameter used on opensource projects

Example 21 with Parameter

use of cbit.vcell.opt.Parameter in project vcell by virtualcell.

the class CurveFitting method bindExpressionToParametersAndTime.

private static void bindExpressionToParametersAndTime(Expression modelExp, Parameter[] parameters) throws ExpressionBindingException {
    ArrayList<String> symbols = new ArrayList<String>();
    symbols.add("t");
    for (Parameter p : parameters) {
        symbols.add(p.getName());
    }
    modelExp.bindExpression(new SimpleSymbolTable(symbols.toArray(new String[0])));
}

Example 22 with Parameter

use of cbit.vcell.opt.Parameter 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;
}

Example 23 with Parameter

use of cbit.vcell.opt.Parameter in project vcell by virtualcell.

the class ReactionDominantTest method solve.

private OptimizationResultSet solve() throws ExpressionException, IOException {
    Expression Fbleached_bleachFastExp = new Expression(Fbleached_bleachFast);
    Expression OmegaExp = new Expression(strOmega);
    Expression Omega_cExp = new Expression(strOmega_c);
    Expression WExp = new Expression(strW);
    Fbleached_bleachFastExp.substituteInPlace(OmegaExp, new Expression(Omega));
    Fbleached_bleachFastExp.substituteInPlace(Omega_cExp, new Expression(Omega_c));
    Fbleached_bleachFastExp.substituteInPlace(WExp, new Expression(W));
    Parameter[] parameters = new Parameter[] { para_alpha, para_Kon_star, para_Koff };
    // Expression Fbleached_bleachFast = Fbleached_bleachFastExp.flatten();
    // choose optimization solver, currently we have Powell and CFSQP
    PowellOptimizationSolver optService = new PowellOptimizationSolver();
    OptimizationSpec optSpec = new OptimizationSpec();
    // create simple reference data
    double[][] realData = new double[2][t.length];
    for (int i = 0; i < t.length; i++) {
        realData[0][i] = t[i];
        realData[1][i] = data_bleached[i];
    }
    String[] colNames = new String[] { "t", "intensity" };
    double[] weights = new double[] { 1.0, 1.0 };
    SimpleReferenceData refData = new SimpleReferenceData(colNames, weights, realData);
    // send to optimization service
    ExplicitFitObjectiveFunction.ExpressionDataPair oneExpDataPair = new ExplicitFitObjectiveFunction.ExpressionDataPair(Fbleached_bleachFastExp.flatten(), 1);
    ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = new ExplicitFitObjectiveFunction.ExpressionDataPair[] { oneExpDataPair };
    optSpec.setObjectiveFunction(new ExplicitFitObjectiveFunction(expDataPairs, refData));
    optSpec.setComputeProfileDistributions(true);
    // 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);
    String[] paramNames = optResultSet.getOptSolverResultSet().getParameterNames();
    ArrayList<ProfileDistribution> profileDistributionList = optResultSet.getOptSolverResultSet().getProfileDistributionList();
    for (int i = 0; i < profileDistributionList.size(); i++) {
        outputProfileLikelihood(profileDistributionList.get(i), i, paramNames[i], new File(fileDir));
    }
    return optResultSet;
}

Example 24 with Parameter

use of cbit.vcell.opt.Parameter in project vcell by virtualcell.

the class FRAPModel method createReacBindingParametersFromDiffusionParameters.

public static Parameter[] createReacBindingParametersFromDiffusionParameters(Parameter[] origParams) {
    Parameter[] params = new Parameter[NUM_MODEL_PARAMETERS_BINDING];
    Parameter primaryDiff = new Parameter(MODEL_PARAMETER_NAMES[INDEX_PRIMARY_DIFF_RATE], REF_DIFFUSION_RATE_PARAM.getLowerBound(), REF_DIFFUSION_RATE_PARAM.getUpperBound(), REF_DIFFUSION_RATE_PARAM.getScale(), origParams[INDEX_PRIMARY_DIFF_RATE].getInitialGuess());
    Parameter primaryFrac = new Parameter(MODEL_PARAMETER_NAMES[INDEX_PRIMARY_FRACTION], REF_MOBILE_FRACTION_PARAM.getLowerBound(), REF_MOBILE_FRACTION_PARAM.getUpperBound(), REF_MOBILE_FRACTION_PARAM.getScale(), origParams[INDEX_PRIMARY_FRACTION].getInitialGuess());
    Parameter bleachWhileMonitoringRate = new Parameter(MODEL_PARAMETER_NAMES[INDEX_BLEACH_MONITOR_RATE], REF_BLEACH_WHILE_MONITOR_PARAM.getLowerBound(), REF_BLEACH_WHILE_MONITOR_PARAM.getUpperBound(), REF_BLEACH_WHILE_MONITOR_PARAM.getScale(), origParams[INDEX_BLEACH_MONITOR_RATE].getInitialGuess());
    Parameter secondaryDiff = null;
    Parameter secondaryFrac = null;
    secondaryDiff = new Parameter(MODEL_PARAMETER_NAMES[INDEX_SECONDARY_DIFF_RATE], REF_DIFFUSION_RATE_PARAM.getLowerBound(), REF_DIFFUSION_RATE_PARAM.getUpperBound(), REF_DIFFUSION_RATE_PARAM.getScale(), 0);
    secondaryFrac = new Parameter(MODEL_PARAMETER_NAMES[INDEX_SECONDARY_FRACTION], REF_MOBILE_FRACTION_PARAM.getLowerBound(), REF_MOBILE_FRACTION_PARAM.getUpperBound(), REF_MOBILE_FRACTION_PARAM.getScale(), 0);
    Parameter bsConcentration = new Parameter(MODEL_PARAMETER_NAMES[INDEX_BINDING_SITE_CONCENTRATION], 0, 1, 1, 0);
    Parameter onReacRate = new Parameter(MODEL_PARAMETER_NAMES[INDEX_ON_RATE], 0, 1e6, 1, 0);
    Parameter offReacRate = new Parameter(MODEL_PARAMETER_NAMES[INDEX_OFF_RATE], 0, 1e6, 1, 0);
    params = new Parameter[FRAPModel.NUM_MODEL_PARAMETERS_BINDING];
    params[FRAPModel.INDEX_PRIMARY_DIFF_RATE] = primaryDiff;
    params[FRAPModel.INDEX_PRIMARY_FRACTION] = primaryFrac;
    params[FRAPModel.INDEX_BLEACH_MONITOR_RATE] = bleachWhileMonitoringRate;
    params[FRAPModel.INDEX_SECONDARY_DIFF_RATE] = secondaryDiff;
    params[FRAPModel.INDEX_SECONDARY_FRACTION] = secondaryFrac;
    params[FRAPModel.INDEX_BINDING_SITE_CONCENTRATION] = bsConcentration;
    params[FRAPModel.INDEX_ON_RATE] = onReacRate;
    params[FRAPModel.INDEX_OFF_RATE] = offReacRate;
    return params;
}

Example 25 with Parameter

use of cbit.vcell.opt.Parameter in project vcell by virtualcell.

the class FRAPOptFunctions method evaluateParameters.

// profileData array contains only two profile distribution, one for bleachWhileMonitoringRate and another for reaction off rate
public ProfileData[] evaluateParameters(Parameter[] currentParams, ClientTaskStatusSupport clientTaskStatusSupport) throws Exception {
    int totalParamLen = currentParams.length;
    int resultDataCounter = 0;
    ProfileData[] resultData = new ProfileData[NUM_PARAM_ESTIMATED];
    FRAPStudy frapStudy = getExpFrapStudy();
    for (int j = 0; j < totalParamLen; j++) {
        // only bleach while monitoring rate and reaction off rate need to evaluated
        if (currentParams[j] != null && (currentParams[j].getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_BLEACH_MONITOR_RATE]) || currentParams[j].getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_OFF_RATE]))) {
            ProfileData profileData = new ProfileData();
            // add the fixed parameter to profileData, output exp data and opt results
            Parameter[] newBestParameters = getBestParamters(frapStudy.getFrapData(), null, true);
            double iniError = getOffRateResults().getObjectiveFunctionValue();
            // fixed parameter(make sure parameter shall not be smaller than epsilon)
            Parameter fixedParam = newBestParameters[j];
            if (// log function cannot take 0 as parameter
            fixedParam.getInitialGuess() == 0) {
                fixedParam = new Parameter(fixedParam.getName(), fixedParam.getLowerBound(), fixedParam.getUpperBound(), fixedParam.getScale(), FRAPOptimizationUtils.epsilon);
            }
            if (clientTaskStatusSupport != null) {
                clientTaskStatusSupport.setMessage("<html>Evaluating confidence intervals of \'" + fixedParam.getName() + "\' <br> of finding reaction off rate model.</html>");
                // start evaluation of a parameter.
                clientTaskStatusSupport.setProgress(0);
            }
            ProfileDataElement pde = new ProfileDataElement(fixedParam.getName(), Math.log10(fixedParam.getInitialGuess()), iniError, newBestParameters);
            profileData.addElement(pde);
            // increase
            int iterationCount = 1;
            double paramLogVal = Math.log10(fixedParam.getInitialGuess());
            double lastError = iniError;
            boolean isBoundReached = false;
            double incrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter];
            int stepIncreaseCount = 0;
            while (true) {
                if (// if exceeds the maximum iterations, break;
                iterationCount > FRAPOptData.MAX_ITERATION) {
                    break;
                }
                if (isBoundReached) {
                    break;
                }
                paramLogVal = paramLogVal + incrementStep;
                double paramVal = Math.pow(10, paramLogVal);
                if (paramVal > (fixedParam.getUpperBound() - FRAPOptimizationUtils.epsilon)) {
                    paramVal = fixedParam.getUpperBound();
                    paramLogVal = Math.log10(fixedParam.getUpperBound());
                    isBoundReached = true;
                }
                Parameter increasedParam = new Parameter(fixedParam.getName(), fixedParam.getLowerBound(), fixedParam.getUpperBound(), fixedParam.getScale(), paramVal);
                // getBestParameters returns the whole set of parameters including the fixed parameters
                Parameter[] newParameters = getBestParamters(frapStudy.getFrapData(), increasedParam, true);
                double error = getOffRateResults().getObjectiveFunctionValue();
                pde = new ProfileDataElement(increasedParam.getName(), paramLogVal, error, newParameters);
                profileData.addElement(pde);
                // check if the we run enough to get confidence intervals(99% @6.635, we plus 10 over the min error)
                if (error > (iniError + 10)) {
                    break;
                }
                if (Math.abs((error - lastError) / lastError) < FRAPOptData.MIN_LIKELIHOOD_CHANGE) {
                    stepIncreaseCount++;
                    incrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter] * Math.pow(2, stepIncreaseCount);
                } else {
                    if (stepIncreaseCount > 1) {
                        incrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter] / Math.pow(2, stepIncreaseCount);
                        stepIncreaseCount--;
                    }
                }
                if (clientTaskStatusSupport.isInterrupted()) {
                    throw UserCancelException.CANCEL_GENERIC;
                }
                lastError = error;
                iterationCount++;
                clientTaskStatusSupport.setProgress((int) ((iterationCount * 1.0 / FRAPOptData.MAX_ITERATION) * 0.5 * 100));
            }
            // half way through evaluation of a parameter.
            clientTaskStatusSupport.setProgress(50);
            // decrease
            iterationCount = 1;
            paramLogVal = Math.log10(fixedParam.getInitialGuess());
            lastError = iniError;
            isBoundReached = false;
            double decrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter];
            stepIncreaseCount = 0;
            while (true) {
                if (// if exceeds the maximum iterations, break;
                iterationCount > FRAPOptData.MAX_ITERATION) {
                    break;
                }
                if (isBoundReached) {
                    break;
                }
                paramLogVal = paramLogVal - decrementStep;
                double paramVal = Math.pow(10, paramLogVal);
                // System.out.println("paramVal:" + paramVal);
                if (paramVal < (fixedParam.getLowerBound() + FRAPOptimizationUtils.epsilon)) {
                    paramVal = fixedParam.getLowerBound();
                    paramLogVal = Math.log10(FRAPOptimizationUtils.epsilon);
                    isBoundReached = true;
                }
                Parameter decreasedParam = new Parameter(fixedParam.getName(), fixedParam.getLowerBound(), fixedParam.getUpperBound(), fixedParam.getScale(), paramVal);
                // getBestParameters returns the whole set of parameters including the fixed parameters
                Parameter[] newParameters = getBestParamters(frapStudy.getFrapData(), decreasedParam, true);
                double error = getOffRateResults().getObjectiveFunctionValue();
                pde = new ProfileDataElement(decreasedParam.getName(), paramLogVal, error, newParameters);
                profileData.addElement(0, pde);
                if (error > (iniError + 10)) {
                    break;
                }
                if (Math.abs((error - lastError) / lastError) < FRAPOptData.MIN_LIKELIHOOD_CHANGE) {
                    stepIncreaseCount++;
                    decrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter] * Math.pow(2, stepIncreaseCount);
                } else {
                    if (stepIncreaseCount > 1) {
                        incrementStep = FRAPOptData.DEFAULT_CI_STEPS_OFF_RATE[resultDataCounter] / Math.pow(2, stepIncreaseCount);
                        stepIncreaseCount--;
                    }
                }
                if (clientTaskStatusSupport.isInterrupted()) {
                    throw UserCancelException.CANCEL_GENERIC;
                }
                lastError = error;
                iterationCount++;
                clientTaskStatusSupport.setProgress((int) (((iterationCount + FRAPOptData.MAX_ITERATION) * 1.0 / FRAPOptData.MAX_ITERATION) * 0.5 * 100));
            }
            resultData[resultDataCounter++] = profileData;
            // finish evaluation of a parameter
            clientTaskStatusSupport.setProgress(100);
        } else {
            continue;
        }
    }
    // this message is specifically set for batchrun, the message will stay in the status panel. It doesn't affect single run,which disappears quickly that user won't notice.
    clientTaskStatusSupport.setMessage("Evaluating confidence intervals ...");
    return resultData;
}