Search in sources :

Example 1 with ElementWeights

use of cbit.vcell.opt.ElementWeights 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;
}
Also used : OptimizationException(cbit.vcell.opt.OptimizationException) ElementWeights(cbit.vcell.opt.ElementWeights) Expression(cbit.vcell.parser.Expression)

Example 2 with ElementWeights

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

the class OptXmlWriter method getDataXML.

public static Element getDataXML(ReferenceData refData) {
    Element refDataElement = null;
    if (refData instanceof SimpleReferenceData) {
        refDataElement = new Element(OptXmlTags.SimpleReferenceData_Tag);
    } else if (refData instanceof SpatialReferenceData) {
        refDataElement = new Element(OptXmlTags.SpatialReferenceData_Tag);
    }
    if (refDataElement != null) {
        // write variable declarations
        // independent variable is t and dimension is 1, these are fixed.
        Element timeVarElement = new Element(OptXmlTags.Variable_Tag);
        timeVarElement.setAttribute(OptXmlTags.VariableType_Attr, OptXmlTags.VariableType_Attr_Independent);
        timeVarElement.setAttribute(OptXmlTags.VariableName_Attr, ReservedVariable.TIME.getName());
        timeVarElement.setAttribute(OptXmlTags.VariableDimension_Attr, "1");
        refDataElement.addContent(timeVarElement);
        // check if t is at the first column
        int timeIndex = refData.findColumn(ReservedVariable.TIME.getName());
        if (timeIndex != 0) {
            throw new RuntimeException("t must be the first column");
        }
        // add all other dependent variables, recall that the dependent variables start from 2nd column onward in reference data
        for (int i = 1; i < refData.getNumDataColumns(); i++) {
            Element variableElement = new Element(OptXmlTags.Variable_Tag);
            variableElement.setAttribute(OptXmlTags.VariableType_Attr, OptXmlTags.VariableType_Attr_Dependent);
            variableElement.setAttribute(OptXmlTags.VariableName_Attr, refData.getColumnNames()[i]);
            variableElement.setAttribute(OptXmlTags.VariableDimension_Attr, refData.getDataSize() + "");
            refDataElement.addContent(variableElement);
        }
        // write data
        Element dataRowElement = new Element(OptXmlTags.Datarow_Tag);
        for (int i = 0; i < refData.getNumDataRows(); i++) {
            double[] data = refData.getDataByRow(i);
            Element rowElement = new Element(OptXmlTags.Row_Tag);
            StringBuffer rowText = new StringBuffer();
            for (int j = 0; j < data.length; j++) {
                rowText.append(data[j] + " ");
            }
            rowElement.addContent(rowText.toString());
            dataRowElement.addContent(rowElement);
        }
        refDataElement.addContent(dataRowElement);
        // write weights
        Element weightRowElement = new Element(OptXmlTags.WeightDatarow_Tag);
        // add weight type
        if (// print weights in multiple rows, each row has one-to-more elements
        refData.getWeights() instanceof ElementWeights) {
            weightRowElement.setAttribute(OptXmlTags.WeightType_Attr, OptXmlTags.WeightType_Attr_Element);
            double[][] weightData = ((ElementWeights) refData.getWeights()).getWeightData();
            // elementWeights: first dimensin is number of rows, second dimension is number of vars
            for (int i = 0; i < weightData.length; i++) {
                double[] rowData = weightData[i];
                Element rowElement = new Element(OptXmlTags.Row_Tag);
                StringBuffer rowText = new StringBuffer();
                for (int j = 0; j < rowData.length; j++) {
                    rowText.append(rowData[j] + " ");
                }
                rowElement.addContent(rowText.toString());
                weightRowElement.addContent(rowElement);
            }
        } else if (// print weights in one row, the row has one-to-more elements
        refData.getWeights() instanceof VariableWeights) {
            weightRowElement.setAttribute(OptXmlTags.WeightType_Attr, OptXmlTags.WeightType_Attr_Variable);
            double[] weightData = ((VariableWeights) refData.getWeights()).getWeightData();
            Element rowElement = new Element(OptXmlTags.Row_Tag);
            StringBuffer rowText = new StringBuffer();
            for (int j = 0; j < weightData.length; j++) {
                rowText.append(weightData[j] + " ");
            }
            rowElement.addContent(rowText.toString());
            weightRowElement.addContent(rowElement);
        } else // time weights. Print weights in multiple rows, each row has one element
        {
            weightRowElement.setAttribute(OptXmlTags.WeightType_Attr, OptXmlTags.WeightType_Attr_Time);
            double[] weightData = ((TimeWeights) refData.getWeights()).getWeightData();
            for (int j = 0; j < weightData.length; j++) {
                Element rowElement = new Element(OptXmlTags.Row_Tag);
                StringBuffer rowText = new StringBuffer();
                rowText.append(weightData[j] + " ");
                rowElement.addContent(rowText.toString());
                weightRowElement.addContent(rowElement);
            }
        }
        refDataElement.addContent(weightRowElement);
    }
    return refDataElement;
}
Also used : VariableWeights(cbit.vcell.opt.VariableWeights) ElementWeights(cbit.vcell.opt.ElementWeights) Element(org.jdom.Element) SimpleReferenceData(cbit.vcell.opt.SimpleReferenceData) SpatialReferenceData(cbit.vcell.opt.SpatialReferenceData) Constraint(cbit.vcell.opt.Constraint)

Aggregations

ElementWeights (cbit.vcell.opt.ElementWeights)2 Constraint (cbit.vcell.opt.Constraint)1 OptimizationException (cbit.vcell.opt.OptimizationException)1 SimpleReferenceData (cbit.vcell.opt.SimpleReferenceData)1 SpatialReferenceData (cbit.vcell.opt.SpatialReferenceData)1 VariableWeights (cbit.vcell.opt.VariableWeights)1 Expression (cbit.vcell.parser.Expression)1 Element (org.jdom.Element)1