Examples with Expression cbit.vcell.parser.Expression used on opensource projects

Example 91 with Expression

use of cbit.vcell.parser.Expression in project vcell by virtualcell.

the class FRAPDataAnalysis method fitRecovery_diffusionOnly.

/**
 * Method fitRecovery2.
 * @param frapData, the original image info.
 * @param arg_bleachType, the gaussian spot or half cell bleaching types.
 * @return FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults
 * @throws ExpressionException
 */
public static FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults fitRecovery_diffusionOnly(FRAPData frapData, int arg_bleachType, int startIndexForRecovery) throws ExpressionException, OptimizationException, IOException, IllegalArgumentException {
    // int startIndexForRecovery = getRecoveryIndex(frapData);
    // 
    // get unnormalized average background fluorescence at each time point
    // 
    double[] temp_background = frapData.getAvgBackGroundIntensity();
    // the prebleachAvg has backgroud subtracted.
    double[] preBleachAvgXYZ = FrapDataUtils.calculatePreBleachAverageXYZ(frapData, startIndexForRecovery);
    // temp_fluor has subtracted background and divided by prebleach average.
    // get average intensity under the bleached area according to each time point
    double[] temp_fluor = getAverageROIIntensity(frapData, frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_BLEACHED.name()), preBleachAvgXYZ, temp_background);
    double[] temp_time = frapData.getImageDataset().getImageTimeStamps();
    // get nomalized preBleachAverage under bleached area.
    double preBleachAverage_bleachedArea = 0.0;
    for (int i = 0; i < startIndexForRecovery; i++) {
        preBleachAverage_bleachedArea += temp_fluor[i];
    }
    preBleachAverage_bleachedArea /= startIndexForRecovery;
    // get number of pixels in bleached region(non ROI region pixels are saved as 0)
    ROI bleachedROI_2D = frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_BLEACHED.name());
    long numPixelsInBleachedROI = bleachedROI_2D.getRoiImages()[0].getNumXYZ() - bleachedROI_2D.getRoiImages()[0].countPixelsByValue((short) 0);
    // assume ROI is a circle, A = Pi*R^2
    // so R = sqrt(A/Pi)
    double imagePixelArea = frapData.getImageDataset().getAllImages()[0].getPixelAreaXY();
    double area = imagePixelArea * numPixelsInBleachedROI;
    // Radius of ROI(assume that ROI is a circle)
    double bleachRadius = Math.sqrt(area / Math.PI);
    // assume cell is a circle, A = Pi*R^2
    // so R = sqrt(A/Pi)
    area = frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name()).getRoiImages()[0].getPixelAreaXY() * (frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name()).getRoiImages()[0].getNumXYZ() - bleachedROI_2D.getRoiImages()[0].countPixelsByValue((short) 0));
    // Radius of ROI(assume that ROI is a circle)
    double cellRadius = Math.sqrt(area / Math.PI);
    // average intensity under bleached area. The time points start from the first post bleach
    double[] fluor = new double[temp_fluor.length - startIndexForRecovery];
    // Time points stat from the first post bleach
    double[] time = new double[temp_time.length - startIndexForRecovery];
    System.arraycopy(temp_fluor, startIndexForRecovery, fluor, 0, fluor.length);
    System.arraycopy(temp_time, startIndexForRecovery, time, 0, time.length);
    FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults diffAnalysisResults = new FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults();
    // frapDataAnalysisResults.setBleachWidth(bleachRadius);
    // frapDataAnalysisResults.setStartingIndexForRecovery(startIndexForRecovery);
    // frapDataAnalysisResults.setBleachRegionData(temp_fluor);// average intensity under bleached region accroding to time points
    diffAnalysisResults.setBleachType(arg_bleachType);
    // curve fitting according to different bleach types
    double[] inputParamValues = null;
    double[] outputParamValues = null;
    // 
    // Bleach while monitoring fit
    // 
    double[] tempCellROIAverage = getAverageROIIntensity(frapData, frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name()), preBleachAvgXYZ, temp_background);
    double[] cellROIAverage = new double[tempCellROIAverage.length - startIndexForRecovery];
    // Cell Avg. points start from the first post bleach
    System.arraycopy(tempCellROIAverage, startIndexForRecovery, cellROIAverage, 0, cellROIAverage.length);
    outputParamValues = new double[1];
    Expression bleachWhileMonitorFitExpression = CurveFitting.fitBleachWhileMonitoring(time, cellROIAverage, outputParamValues);
    double bleachWhileMonitoringRate = outputParamValues[0];
    diffAnalysisResults.setBleachWhileMonitoringTau(bleachWhileMonitoringRate);
    diffAnalysisResults.setFitBleachWhileMonitorExpression(bleachWhileMonitorFitExpression.flatten());
    // to fit If, Io and recovery tau by fitting expressions based on bleaching types
    if (arg_bleachType == FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_GaussianSpot) {
        double cellAreaBleached = getCellAreaBleachedFraction(frapData);
        // the input parameter is the bleach while monitoring rate
        inputParamValues = new double[] { bleachWhileMonitoringRate };
        // the array is used to get If, Io, and tau back.
        outputParamValues = new double[3];
        Expression fittedCurve = CurveFitting.fitRecovery_diffOnly(time, fluor, FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_GaussianSpot, inputParamValues, outputParamValues);
        // get diffusion rate
        double fittedRecoveryTau = outputParamValues[2];
        Expression diffExp = new Expression(FRAPDataAnalysis.gaussianSpot_DiffFunc);
        diffExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_w), new Expression(bleachRadius));
        diffExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_tau), new Expression(fittedRecoveryTau));
        double fittedDiffusionRate = diffExp.evaluateConstant();
        // double fittedDiffusionRate = bleachRadius*bleachRadius/(4.0*fittedRecoveryTau);
        // get mobile fraction
        double If = outputParamValues[0];
        double Io = outputParamValues[1];
        Expression mFracExp = new Expression(FRAPDataAnalysis.gaussianSpot_mobileFracFunc);
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.para_If.getName()), new Expression(If));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.para_Io.getName()), new Expression(Io));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_preBlchAvg), new Expression(preBleachAverage_bleachedArea));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_fB), new Expression(cellAreaBleached));
        double mobileFrac = mFracExp.evaluateConstant();
        // sometimes the mobile fraction goes beyond [0,1], we have to restrict the mobile fraction value.
        mobileFrac = Math.min(1, mobileFrac);
        mobileFrac = Math.max(0, mobileFrac);
        // set frap data analysis results
        diffAnalysisResults.setDiffFitExpression(fittedCurve.flatten());
        diffAnalysisResults.setRecoveryTau(fittedRecoveryTau);
        diffAnalysisResults.setRecoveryDiffusionRate(fittedDiffusionRate);
        diffAnalysisResults.setMobilefraction(mobileFrac);
    } else if (arg_bleachType == FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_HalfCell) {
        double cellAreaBleached = getCellAreaBleachedFraction(frapData);
        // the input parameter is the bleach while monitoring rate
        inputParamValues = new double[] { bleachWhileMonitoringRate };
        // the array is used get If, Io, and tau back.
        outputParamValues = new double[3];
        Expression fittedCurve = CurveFitting.fitRecovery_diffOnly(time, fluor, FrapDataAnalysisResults.DiffusionOnlyAnalysisRestults.BleachType_HalfCell, inputParamValues, outputParamValues);
        // get diffusion rate
        double fittedRecoveryTau = outputParamValues[2];
        Expression diffExp = new Expression(FRAPDataAnalysis.halfCell_DiffFunc);
        diffExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_r), new Expression(cellRadius));
        diffExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_tau), new Expression(fittedRecoveryTau));
        diffExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_Pi), new Expression(Math.PI));
        double fittedDiffusionRate = diffExp.evaluateConstant();
        // double fittedDiffusionRate = (cellRadius*cellRadius)/(fittedRecoveryTau*Math.PI*Math.PI);
        // get mobile fraction
        double If = outputParamValues[0];
        double Io = outputParamValues[1];
        Expression mFracExp = new Expression(FRAPDataAnalysis.gaussianSpot_mobileFracFunc);
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.para_If.getName()), new Expression(If));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.para_Io.getName()), new Expression(Io));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_preBlchAvg), new Expression(preBleachAverage_bleachedArea));
        mFracExp.substituteInPlace(new Expression(FRAPDataAnalysis.symbol_fB), new Expression(cellAreaBleached));
        double mobileFrac = mFracExp.evaluateConstant();
        // set diffusion only analysis results
        diffAnalysisResults.setDiffFitExpression(fittedCurve.flatten());
        diffAnalysisResults.setRecoveryTau(fittedRecoveryTau);
        diffAnalysisResults.setRecoveryDiffusionRate(fittedDiffusionRate);
        diffAnalysisResults.setMobilefraction(mobileFrac);
    } else {
        throw new IllegalArgumentException("Unknown Bleach Type " + arg_bleachType);
    }
    return diffAnalysisResults;
}

Example 92 with Expression

use of cbit.vcell.parser.Expression in project vcell by virtualcell.

the class FRAPDataAnalysis method fitRecovery_reacOffRateOnly.

/**
 * Method fitRecovery2.
 * @param frapData, the original image info.
 * @param fixedParameter: the fixed parameter from profile likelihood distribution analysis.
 * @param measurementError: the measurementError is used as weights for calculating objectiveFunction errors.
 * @return FrapDataAnalysisResults.ReactionOnlyAnalysisRestults
 * @throws ExpressionException
 */
public static FrapDataAnalysisResults.ReactionOnlyAnalysisRestults fitRecovery_reacOffRateOnly(FRAPData frapData, Parameter fixedParam, double[][] measurementError, int startIndexForRecovery) throws ExpressionException, OptimizationException, IOException {
    // int startIndexForRecovery = getRecoveryIndex(frapData);
    // 
    // get unnormalized average background fluorescence at each time point
    // 
    double[] temp_background = frapData.getAvgBackGroundIntensity();
    // the prebleachAvg has backgroud subtracted.
    double[] preBleachAvgXYZ = FrapDataUtils.calculatePreBleachAverageXYZ(frapData, startIndexForRecovery);
    // tempBeachedAverage and tempCellROIAverage have subtracted background and divided by prebleach average. the following array has data for the full time duration.
    double[] tempBeachedAverage = getAverageROIIntensity(frapData, frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_BLEACHED.name()), preBleachAvgXYZ, temp_background);
    double[] tempCellROIAverage = getAverageROIIntensity(frapData, frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name()), preBleachAvgXYZ, temp_background);
    double[] temp_time = frapData.getImageDataset().getImageTimeStamps();
    // get bleached and cell roi data starting from first post bleach
    // The time points start from the first post bleach
    double[] bleachedAverage = new double[tempBeachedAverage.length - startIndexForRecovery];
    // time points start from the first post bleach
    double[] cellROIAverage = new double[tempCellROIAverage.length - startIndexForRecovery];
    // Time points stat from the first post bleach
    double[] time = new double[temp_time.length - startIndexForRecovery];
    System.arraycopy(tempBeachedAverage, startIndexForRecovery, bleachedAverage, 0, bleachedAverage.length);
    System.arraycopy(tempCellROIAverage, startIndexForRecovery, cellROIAverage, 0, cellROIAverage.length);
    System.arraycopy(temp_time, startIndexForRecovery, time, 0, time.length);
    // initialize reaction off rate analysis results
    FrapDataAnalysisResults.ReactionOnlyAnalysisRestults offRateAnalysisResults = new FrapDataAnalysisResults.ReactionOnlyAnalysisRestults();
    /**
     *curve fitting
     */
    // index 0: cell ROI intensity average at time 0, I_cell_ini. index 1: bleached intensity average at time 0, I_bleached_ini
    double[] inputParamValues = new double[] { cellROIAverage[0], bleachedAverage[0] };
    // if fixed parameter is null, then outputs are bwm rate, koff rate, fitting parameter A. Otherwise outputs are two out of three.
    double[] outputParamValues = null;
    // create data array,first col is cell average, second col is bleached average.
    double[][] fitData = new double[2][];
    fitData[0] = cellROIAverage;
    fitData[1] = bleachedAverage;
    // create element weights array, first col is cell data weights, second col is bleached data weights
    double[][] weightData = new double[time.length][2];
    double[] cellROIMeasurementError = measurementError[FRAPData.VFRAP_ROI_ENUM.ROI_CELL.ordinal()];
    double[] bleachedROIMeasurementError = measurementError[FRAPData.VFRAP_ROI_ENUM.ROI_BLEACHED.ordinal()];
    for (// elementWeight, first dimension is number of rows(time points), second dimension is number of variables(fitDatasets)
    int i = 0; // elementWeight, first dimension is number of rows(time points), second dimension is number of variables(fitDatasets)
    i < time.length; // elementWeight, first dimension is number of rows(time points), second dimension is number of variables(fitDatasets)
    i++) {
        weightData[i][0] = 1 / (cellROIMeasurementError[i] * cellROIMeasurementError[i]);
        weightData[i][1] = 1 / (bleachedROIMeasurementError[i] * bleachedROIMeasurementError[i]);
    }
    ElementWeights eleWeights = new ElementWeights(weightData);
    if (fixedParam == null) {
        // call curvefitting
        // bwmrate, fitting parameter A & reaction off rate
        outputParamValues = new double[3];
        double error = CurveFitting.fitRecovery_reacKoffRateOnly(time, fitData, inputParamValues, outputParamValues, fixedParam, eleWeights);
        // set objective function value
        offRateAnalysisResults.setObjectiveFunctionValue(error);
        // set reaction off rate analysis results
        offRateAnalysisResults.setBleachWhileMonitoringTau(outputParamValues[0]);
        offRateAnalysisResults.setFittingParamA(outputParamValues[1]);
        offRateAnalysisResults.setOffRate(outputParamValues[2]);
        // set cell intensity expression ( only t left in expression)
        Expression cellIntensityExp = new Expression(FRAPOptFunctions.FUNC_CELL_INTENSITY);
        cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(cellROIAverage[0]));
        // subsitute bwmRate
        cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(outputParamValues[0]));
        // undo time shift
        cellIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
        offRateAnalysisResults.setFitBleachWhileMonitorExpression(cellIntensityExp);
        // set bleached region intensity expression ( only t left in expression)
        Expression bleachIntensityExp = new Expression(FRAPOptFunctions.FUNC_RECOVERY_BLEACH_REACTION_DOMINANT);
        bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(bleachedAverage[0]));
        // subsitute bwmRate
        bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(outputParamValues[0]));
        // subsitute parameter A
        bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_A), new Expression(outputParamValues[1]));
        // reaction off rate
        bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_KOFF), new Expression(outputParamValues[2]));
        // undo time shift
        bleachIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
        offRateAnalysisResults.setOffRateFitExpression(bleachIntensityExp);
    } else {
        if (fixedParam != null && fixedParam.getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_BLEACH_MONITOR_RATE])) {
            // fitting parameter A & reaction off rate
            outputParamValues = new double[2];
            double error = CurveFitting.fitRecovery_reacKoffRateOnly(time, fitData, inputParamValues, outputParamValues, fixedParam, eleWeights);
            // set objective function value
            offRateAnalysisResults.setObjectiveFunctionValue(error);
            // set reaction off rate analysis results
            offRateAnalysisResults.setBleachWhileMonitoringTau(fixedParam.getInitialGuess());
            offRateAnalysisResults.setFittingParamA(outputParamValues[0]);
            offRateAnalysisResults.setOffRate(outputParamValues[1]);
            // set cell intensity expression ( only t left in expression)
            Expression cellIntensityExp = new Expression(FRAPOptFunctions.FUNC_CELL_INTENSITY);
            cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(cellROIAverage[0]));
            // subsitute bwmRate
            cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(fixedParam.getInitialGuess()));
            // undo time shift
            cellIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
            offRateAnalysisResults.setFitBleachWhileMonitorExpression(cellIntensityExp);
            // set bleached region intensity expression ( only t left in expression)
            Expression bleachIntensityExp = new Expression(FRAPOptFunctions.FUNC_RECOVERY_BLEACH_REACTION_DOMINANT);
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(bleachedAverage[0]));
            // subsitute bwmRate
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(fixedParam.getInitialGuess()));
            // subsitute parameter A
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_A), new Expression(outputParamValues[0]));
            // reaction off rate
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_KOFF), new Expression(outputParamValues[1]));
            // undo time shift
            bleachIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
            offRateAnalysisResults.setOffRateFitExpression(bleachIntensityExp);
        } else if (fixedParam != null && fixedParam.getName().equals(FRAPModel.MODEL_PARAMETER_NAMES[FRAPModel.INDEX_OFF_RATE])) {
            // bwmRate & fitting parameter A
            outputParamValues = new double[2];
            double error = CurveFitting.fitRecovery_reacKoffRateOnly(time, fitData, inputParamValues, outputParamValues, fixedParam, eleWeights);
            // set objective function value
            offRateAnalysisResults.setObjectiveFunctionValue(error);
            // set reaction off rate analysis results
            offRateAnalysisResults.setBleachWhileMonitoringTau(outputParamValues[0]);
            offRateAnalysisResults.setFittingParamA(outputParamValues[1]);
            offRateAnalysisResults.setOffRate(fixedParam.getInitialGuess());
            // set cell intensity expression ( only t left in expression)
            Expression cellIntensityExp = new Expression(FRAPOptFunctions.FUNC_CELL_INTENSITY);
            cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inicell), new Expression(cellROIAverage[0]));
            // subsitute bwmRate
            cellIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(outputParamValues[0]));
            // undo time shift
            cellIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
            offRateAnalysisResults.setFitBleachWhileMonitorExpression(cellIntensityExp);
            // set bleached region intensity expression ( only t left in expression)
            Expression bleachIntensityExp = new Expression(FRAPOptFunctions.FUNC_RECOVERY_BLEACH_REACTION_DOMINANT);
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(bleachedAverage[0]));
            // subsitute bwmRate
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(outputParamValues[0]));
            // subsitute parameter A
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_A), new Expression(outputParamValues[1]));
            // reaction off rate
            bleachIntensityExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_KOFF), new Expression(fixedParam.getInitialGuess()));
            // undo time shift
            bleachIntensityExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + time[0]));
            offRateAnalysisResults.setOffRateFitExpression(bleachIntensityExp);
        } else {
            throw new OptimizationException("Unknown fixed parameter:" + fixedParam.getName());
        }
    }
    return offRateAnalysisResults;
}

Example 93 with Expression

use of cbit.vcell.parser.Expression in project vcell by virtualcell.

the class FRAPOptFunctions method getRecoveryExpressionWithCurrentParameters.

public Expression getRecoveryExpressionWithCurrentParameters(Parameter[] currentParams) throws ExpressionException {
    FRAPData frapData = getExpFrapStudy().getFrapData();
    double[] frapDataTimeStamps = frapData.getImageDataset().getImageTimeStamps();
    // Experiment - Cell ROI Average
    double[] temp_background = frapData.getAvgBackGroundIntensity();
    int startIndexRecovery = getExpFrapStudy().getStartingIndexForRecovery();
    double[] preBleachAvgXYZ = FrapDataUtils.calculatePreBleachAverageXYZ(frapData, startIndexRecovery);
    double[] bleachRegionData = FRAPDataAnalysis.getAverageROIIntensity(frapData, frapData.getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_BLEACHED.name()), preBleachAvgXYZ, temp_background);
    ;
    Expression bleachedAvgExp = new Expression(FRAPOptFunctions.FUNC_RECOVERY_BLEACH_REACTION_DOMINANT);
    // substitute parameter values
    bleachedAvgExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_I_inibleached), new Expression(bleachRegionData[startIndexRecovery]));
    bleachedAvgExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_BWM_RATE), new Expression(currentParams[FRAPModel.INDEX_BLEACH_MONITOR_RATE].getInitialGuess()));
    bleachedAvgExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_A), new Expression(currentParams[FRAPModel.INDEX_BINDING_SITE_CONCENTRATION].getInitialGuess()));
    bleachedAvgExp.substituteInPlace(new Expression(FRAPOptFunctions.SYMBOL_KOFF), new Expression(currentParams[FRAPModel.INDEX_OFF_RATE].getInitialGuess()));
    // time shift
    bleachedAvgExp.substituteInPlace(new Expression(ReservedVariable.TIME.getName()), new Expression(ReservedVariable.TIME.getName() + "-" + frapDataTimeStamps[startIndexRecovery]));
    return bleachedAvgExp;
}

Example 94 with Expression

use of cbit.vcell.parser.Expression in project vcell by virtualcell.

the class FRAPOptFunctions method getFitData.

// fitExp contains time parameter only, all other parameter are substituted with parameter values.
public double[][] getFitData(Parameter[] currentParams) throws DivideByZeroException, ExpressionException {
    Expression fitExp = getRecoveryExpressionWithCurrentParameters(currentParams);
    double[] frapDataTimeStamps = getExpFrapStudy().getFrapData().getImageDataset().getImageTimeStamps();
    int startIndexRecovery = getExpFrapStudy().getStartingIndexForRecovery();
    double[] truncatedTimes = new double[frapDataTimeStamps.length - startIndexRecovery];
    for (int i = startIndexRecovery; i < frapDataTimeStamps.length; i++) {
        truncatedTimes[i - startIndexRecovery] = frapDataTimeStamps[i];
    }
    // System.out.println("error" + getWeightedError(fitExp.flatten()));
    return createData(fitExp.flatten(), truncatedTimes);
}

Example 95 with Expression

use of cbit.vcell.parser.Expression in project vcell by virtualcell.

the class FRAPStudy method createNewRefBioModel.

public static BioModel createNewRefBioModel(FRAPStudy sourceFrapStudy, String baseDiffusionRate, TimeStep tStep, KeyValue simKey, User owner, FieldDataIdentifierSpec psfFDIS, int startingIndexForRecovery) throws Exception {
    if (owner == null) {
        throw new Exception("Owner is not defined");
    }
    ROI cellROI_2D = sourceFrapStudy.getFrapData().getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name());
    Extent extent = sourceFrapStudy.getFrapData().getImageDataset().getExtent();
    TimeBounds timeBounds = FRAPOptData.getEstimatedRefTimeBound(sourceFrapStudy);
    double timeStepVal = FRAPOptData.REFERENCE_DIFF_DELTAT;
    int numX = cellROI_2D.getRoiImages()[0].getNumX();
    int numY = cellROI_2D.getRoiImages()[0].getNumY();
    int numZ = cellROI_2D.getRoiImages().length;
    short[] shortPixels = cellROI_2D.getRoiImages()[0].getPixels();
    byte[] bytePixels = new byte[numX * numY * numZ];
    final byte EXTRACELLULAR_PIXVAL = 0;
    final byte CYTOSOL_PIXVAL = 1;
    for (int i = 0; i < bytePixels.length; i++) {
        if (shortPixels[i] != 0) {
            bytePixels[i] = CYTOSOL_PIXVAL;
        }
    }
    VCImage maskImage;
    try {
        maskImage = new VCImageUncompressed(null, bytePixels, extent, numX, numY, numZ);
    } catch (ImageException e) {
        e.printStackTrace();
        throw new RuntimeException("failed to create mask image for geometry");
    }
    Geometry geometry = new Geometry("geometry", maskImage);
    if (geometry.getGeometrySpec().getNumSubVolumes() != 2) {
        throw new Exception("Cell ROI has no ExtraCellular.");
    }
    int subVolume0PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(0)).getPixelValue();
    geometry.getGeometrySpec().getSubVolume(0).setName((subVolume0PixVal == EXTRACELLULAR_PIXVAL ? EXTRACELLULAR_NAME : CYTOSOL_NAME));
    int subVolume1PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(1)).getPixelValue();
    geometry.getGeometrySpec().getSubVolume(1).setName((subVolume1PixVal == CYTOSOL_PIXVAL ? CYTOSOL_NAME : EXTRACELLULAR_NAME));
    geometry.getGeometrySurfaceDescription().updateAll();
    BioModel bioModel = new BioModel(null);
    bioModel.setName("unnamed");
    Model model = new Model("model");
    bioModel.setModel(model);
    Feature extracellular = model.addFeature(EXTRACELLULAR_NAME);
    Feature cytosol = model.addFeature(CYTOSOL_NAME);
    Membrane plasmaMembrane = model.addMembrane(PLASMAMEMBRANE_NAME);
    String roiDataName = FRAPStudy.ROI_EXTDATA_NAME;
    final int ONE_DIFFUSION_SPECIES_COUNT = 1;
    final int MOBILE_SPECIES_INDEX = 0;
    Expression[] diffusionConstants = new Expression[ONE_DIFFUSION_SPECIES_COUNT];
    Species[] species = new Species[ONE_DIFFUSION_SPECIES_COUNT];
    SpeciesContext[] speciesContexts = new SpeciesContext[ONE_DIFFUSION_SPECIES_COUNT];
    Expression[] initialConditions = new Expression[ONE_DIFFUSION_SPECIES_COUNT];
    // Mobile Species
    diffusionConstants[MOBILE_SPECIES_INDEX] = new Expression(baseDiffusionRate);
    species[MOBILE_SPECIES_INDEX] = new Species(SPECIES_NAME_PREFIX_MOBILE, "Mobile bleachable species");
    speciesContexts[MOBILE_SPECIES_INDEX] = new SpeciesContext(null, species[MOBILE_SPECIES_INDEX].getCommonName(), species[MOBILE_SPECIES_INDEX], cytosol);
    FieldFunctionArguments postBleach_first = new FieldFunctionArguments(roiDataName, "postbleach_first", new Expression(0), VariableType.VOLUME);
    FieldFunctionArguments prebleach_avg = new FieldFunctionArguments(roiDataName, "prebleach_avg", new Expression(0), VariableType.VOLUME);
    Expression expPostBleach_first = new Expression(postBleach_first.infix());
    Expression expPreBleach_avg = new Expression(prebleach_avg.infix());
    initialConditions[MOBILE_SPECIES_INDEX] = Expression.div(expPostBleach_first, expPreBleach_avg);
    SimulationContext simContext = new SimulationContext(bioModel.getModel(), geometry);
    bioModel.addSimulationContext(simContext);
    FeatureMapping cytosolFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(cytosol);
    FeatureMapping extracellularFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(extracellular);
    MembraneMapping plasmaMembraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(plasmaMembrane);
    SubVolume cytSubVolume = geometry.getGeometrySpec().getSubVolume(CYTOSOL_NAME);
    SubVolume exSubVolume = geometry.getGeometrySpec().getSubVolume(EXTRACELLULAR_NAME);
    SurfaceClass pmSurfaceClass = geometry.getGeometrySurfaceDescription().getSurfaceClass(exSubVolume, cytSubVolume);
    cytosolFeatureMapping.setGeometryClass(cytSubVolume);
    extracellularFeatureMapping.setGeometryClass(exSubVolume);
    plasmaMembraneMapping.setGeometryClass(pmSurfaceClass);
    cytosolFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    extracellularFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    plasmaMembraneMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    for (int i = 0; i < initialConditions.length; i++) {
        model.addSpecies(species[i]);
        model.addSpeciesContext(speciesContexts[i]);
    }
    for (int i = 0; i < speciesContexts.length; i++) {
        SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i]);
        scs.getInitialConditionParameter().setExpression(initialConditions[i]);
        scs.getDiffusionParameter().setExpression(diffusionConstants[i]);
    }
    MathMapping mathMapping = simContext.createNewMathMapping();
    MathDescription mathDesc = mathMapping.getMathDescription();
    // Add PSF function
    mathDesc.addVariable(new Function(Simulation.PSF_FUNCTION_NAME, new Expression(psfFDIS.getFieldFuncArgs().infix()), null));
    simContext.setMathDescription(mathDesc);
    SimulationVersion simVersion = new SimulationVersion(simKey, "sim1", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
    Simulation newSimulation = new Simulation(simVersion, simContext.getMathDescription());
    newSimulation.getSolverTaskDescription().setSolverDescription(SolverDescription.FiniteVolumeStandalone);
    simContext.addSimulation(newSimulation);
    newSimulation.getSolverTaskDescription().setTimeBounds(timeBounds);
    newSimulation.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec(timeStepVal));
    newSimulation.getMeshSpecification().setSamplingSize(cellROI_2D.getISize());
    newSimulation.getSolverTaskDescription().setTimeStep(new TimeStep(timeStepVal, timeStepVal, timeStepVal));
    return bioModel;
}