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;
}
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;
}
Aggregations