use of cbit.vcell.mapping.FastSystemAnalyzer in project vcell by virtualcell.
the class DefaultODESolver method initialize.
/**
* This method was created by a SmartGuide.
* @exception SolverException The exception description.
*/
protected void initialize() throws SolverException {
SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
Simulation sim = simSymbolTable.getSimulation();
try {
// create a fast system if necessary
fieldFastAlgebraicSystem = null;
if (getSubDomain().getFastSystem() != null) {
fieldFastAlgebraicSystem = new FastAlgebraicSystem(new FastSystemAnalyzer(getSubDomain().getFastSystem(), simSymbolTable));
}
// refreshIdentifiers();
fieldIdentifiers = createIdentifiers();
fieldSensVariables = createSensitivityVariables();
// refreshStateVariables();
fieldStateVariables = createStateVariables();
//
// allocate ValueVectors object
fieldValueVectors = new ValueVectors(getValueVectorCount(), fieldIdentifiers.size());
// initialize indexes of variables
fieldVariableIndexes = new int[getStateVariableCount()];
for (int i = 0; i < getStateVariableCount(); i++) {
fieldVariableIndexes[i] = getStateVariable(i).getVariable().getIndex();
}
// initialize constants
double[] initialValues = getValueVector(0);
for (int i = 0; i < fieldIdentifiers.size(); i++) {
if (fieldIdentifiers.elementAt(i) instanceof Constant) {
Constant constant = (Constant) fieldIdentifiers.elementAt(i);
constant.bind(simSymbolTable);
if (constant.isConstant()) {
// constant.getValue();
initialValues[constant.getIndex()] = constant.getExpression().evaluateConstant();
} else {
throw new SolverException("cannot evaluate constant '" + constant.getName() + "' = " + constant.getExpression());
}
}
}
// initialize variables
for (int i = 0; i < getStateVariableCount(); i++) {
initialValues[getVariableIndex(i)] = getStateVariable(i).evaluateIC(initialValues);
}
fieldODESolverResultSet = createODESolverResultSet();
// reset - in the ** default ** solvers we don't pick up from where we left off, we can override that behaviour in integrate() if ever necessary
fieldCurrentTime = sim.getSolverTaskDescription().getTimeBounds().getStartingTime();
} catch (ExpressionException expressionException) {
expressionException.printStackTrace(System.out);
throw new SolverException(expressionException.getMessage());
} catch (MathException mathException) {
mathException.printStackTrace(System.out);
throw new SolverException(mathException.getMessage());
}
}
use of cbit.vcell.mapping.FastSystemAnalyzer in project vcell by virtualcell.
the class IDAFileWriter method writeEquations.
/**
* Insert the method's description here.
* Creation date: (3/8/00 10:31:52 PM)
*/
protected String writeEquations(HashMap<Discontinuity, String> discontinuityNameMap) throws MathException, ExpressionException {
Simulation simulation = simTask.getSimulation();
StringBuffer sb = new StringBuffer();
MathDescription mathDescription = simulation.getMathDescription();
if (mathDescription.hasFastSystems()) {
//
// define vector of original variables
//
SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDescription.getSubDomains().nextElement();
FastSystem fastSystem = subDomain.getFastSystem();
FastSystemAnalyzer fs_Analyzer = new FastSystemAnalyzer(fastSystem, simSymbolTable);
int numIndependent = fs_Analyzer.getNumIndependentVariables();
int systemDim = mathDescription.getStateVariableNames().size();
int numDependent = systemDim - numIndependent;
//
// get all variables from fast system (dependent and independent)
//
HashSet<String> fastSystemVarHash = new HashSet<String>();
Enumeration<Variable> dependentfastSystemVarEnum = fs_Analyzer.getDependentVariables();
while (dependentfastSystemVarEnum.hasMoreElements()) {
fastSystemVarHash.add(dependentfastSystemVarEnum.nextElement().getName());
}
Enumeration<Variable> independentfastSystemVarEnum = fs_Analyzer.getIndependentVariables();
while (independentfastSystemVarEnum.hasMoreElements()) {
fastSystemVarHash.add(independentfastSystemVarEnum.nextElement().getName());
}
//
// get all equations including for variables that are not in the fastSystem (ode equations for "slow system")
//
RationalExpMatrix origInitVector = new RationalExpMatrix(systemDim, 1);
RationalExpMatrix origSlowRateVector = new RationalExpMatrix(systemDim, 1);
RationalExpMatrix origVarColumnVector = new RationalExpMatrix(systemDim, 1);
Enumeration<Equation> enumEquations = subDomain.getEquations();
int varIndex = 0;
while (enumEquations.hasMoreElements()) {
Equation equation = enumEquations.nextElement();
origVarColumnVector.set_elem(varIndex, 0, new RationalExp(equation.getVariable().getName()));
Expression rateExpr = equation.getRateExpression();
rateExpr.bindExpression(varsSymbolTable);
rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable);
origSlowRateVector.set_elem(varIndex, 0, new RationalExp("(" + rateExpr.flatten().infix() + ")"));
Expression initExpr = new Expression(equation.getInitialExpression());
initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable).flatten();
origInitVector.set_elem(varIndex, 0, new RationalExp("(" + initExpr.flatten().infix() + ")"));
varIndex++;
}
//
// make symbolic matrix for fast invariants (from FastSystem's fastInvariants as well as a new fast invariant for each variable not included in the fast system.
//
RationalExpMatrix fastInvarianceMatrix = new RationalExpMatrix(numDependent, systemDim);
int row = 0;
for (int i = 0; i < origVarColumnVector.getNumRows(); i++) {
//
if (!fastSystemVarHash.contains(origVarColumnVector.get(i, 0).infixString())) {
fastInvarianceMatrix.set_elem(row, i, RationalExp.ONE);
row++;
}
}
Enumeration<FastInvariant> enumFastInvariants = fastSystem.getFastInvariants();
while (enumFastInvariants.hasMoreElements()) {
FastInvariant fastInvariant = enumFastInvariants.nextElement();
Expression fastInvariantExpression = fastInvariant.getFunction();
for (int col = 0; col < systemDim; col++) {
Expression coeff = fastInvariantExpression.differentiate(origVarColumnVector.get(col, 0).infixString()).flatten();
coeff = simSymbolTable.substituteFunctions(coeff);
fastInvarianceMatrix.set_elem(row, col, RationalExpUtils.getRationalExp(coeff));
}
row++;
}
for (int i = 0; i < systemDim; i++) {
sb.append("VAR " + origVarColumnVector.get(i, 0).infixString() + " INIT " + origInitVector.get(i, 0).infixString() + ";\n");
}
RationalExpMatrix fullMatrix = null;
RationalExpMatrix inverseFullMatrix = null;
RationalExpMatrix newSlowRateVector = null;
try {
RationalExpMatrix fastMat = ((RationalExpMatrix) fastInvarianceMatrix.transpose().findNullSpace());
fullMatrix = new RationalExpMatrix(systemDim, systemDim);
row = 0;
for (int i = 0; i < fastInvarianceMatrix.getNumRows(); i++) {
for (int col = 0; col < systemDim; col++) {
fullMatrix.set_elem(row, col, fastInvarianceMatrix.get(i, col));
}
row++;
}
for (int i = 0; i < fastMat.getNumRows(); i++) {
for (int col = 0; col < systemDim; col++) {
fullMatrix.set_elem(row, col, fastMat.get(i, col));
}
row++;
}
inverseFullMatrix = new RationalExpMatrix(systemDim, systemDim);
inverseFullMatrix.identity();
RationalExpMatrix copyOfFullMatrix = new RationalExpMatrix(fullMatrix);
copyOfFullMatrix.gaussianElimination(inverseFullMatrix);
newSlowRateVector = new RationalExpMatrix(numDependent, 1);
newSlowRateVector.matmul(fastInvarianceMatrix, origSlowRateVector);
} catch (MatrixException ex) {
ex.printStackTrace();
throw new MathException(ex.getMessage());
}
sb.append("TRANSFORM\n");
for (row = 0; row < systemDim; row++) {
for (int col = 0; col < systemDim; col++) {
sb.append(fullMatrix.get(row, col).getConstant().doubleValue() + " ");
}
sb.append("\n");
}
sb.append("INVERSETRANSFORM\n");
for (row = 0; row < systemDim; row++) {
for (int col = 0; col < systemDim; col++) {
sb.append(inverseFullMatrix.get(row, col).getConstant().doubleValue() + " ");
}
sb.append("\n");
}
int numDifferential = numDependent;
int numAlgebraic = numIndependent;
sb.append("RHS DIFFERENTIAL " + numDifferential + " ALGEBRAIC " + numAlgebraic + "\n");
int equationIndex = 0;
while (equationIndex < numDependent) {
// print row of mass matrix followed by slow rate corresponding to fast invariant
Expression slowRateExp = new Expression(newSlowRateVector.get(equationIndex, 0).infixString()).flatten();
slowRateExp.bindExpression(simSymbolTable);
slowRateExp = MathUtilities.substituteFunctions(slowRateExp, varsSymbolTable).flatten();
Vector<Discontinuity> v = slowRateExp.getDiscontinuities();
for (Discontinuity od : v) {
od = getSubsitutedAndFlattened(od, varsSymbolTable);
String dname = discontinuityNameMap.get(od);
if (dname == null) {
dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
discontinuityNameMap.put(od, dname);
}
slowRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
}
sb.append(slowRateExp.infix() + ";\n");
equationIndex++;
}
Enumeration<FastRate> enumFastRates = fastSystem.getFastRates();
while (enumFastRates.hasMoreElements()) {
// print the fastRate for this row
Expression fastRateExp = new Expression(enumFastRates.nextElement().getFunction());
fastRateExp = MathUtilities.substituteFunctions(fastRateExp, varsSymbolTable).flatten();
Vector<Discontinuity> v = fastRateExp.getDiscontinuities();
for (Discontinuity od : v) {
od = getSubsitutedAndFlattened(od, varsSymbolTable);
String dname = discontinuityNameMap.get(od);
if (dname == null) {
dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
discontinuityNameMap.put(od, dname);
}
fastRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
}
sb.append(fastRateExp.flatten().infix() + ";\n");
equationIndex++;
}
} else {
for (int i = 0; i < getStateVariableCount(); i++) {
StateVariable stateVar = getStateVariable(i);
Expression initExpr = new Expression(stateVar.getInitialRateExpression());
initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable);
initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
sb.append("VAR " + stateVar.getVariable().getName() + " INIT " + initExpr.flatten().infix() + ";\n");
}
sb.append("TRANSFORM\n");
for (int row = 0; row < getStateVariableCount(); row++) {
for (int col = 0; col < getStateVariableCount(); col++) {
sb.append((row == col ? 1 : 0) + " ");
}
sb.append("\n");
}
sb.append("INVERSETRANSFORM\n");
for (int row = 0; row < getStateVariableCount(); row++) {
for (int col = 0; col < getStateVariableCount(); col++) {
sb.append((row == col ? 1 : 0) + " ");
}
sb.append("\n");
}
sb.append("RHS DIFFERENTIAL " + getStateVariableCount() + " ALGEBRAIC 0\n");
for (int i = 0; i < getStateVariableCount(); i++) {
StateVariable stateVar = getStateVariable(i);
Expression rateExpr = new Expression(stateVar.getRateExpression());
rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable).flatten();
Vector<Discontinuity> v = rateExpr.getDiscontinuities();
for (Discontinuity od : v) {
od = getSubsitutedAndFlattened(od, varsSymbolTable);
String dname = discontinuityNameMap.get(od);
if (dname == null) {
dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
discontinuityNameMap.put(od, dname);
}
rateExpr.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
}
sb.append(rateExpr.infix() + ";\n");
}
}
return sb.toString();
}
use of cbit.vcell.mapping.FastSystemAnalyzer in project vcell by virtualcell.
the class FiniteVolumeFileWriter method writeFastSystem.
/*private void writeDataProcessor() throws DataAccessException, IOException, MathException, DivideByZeroException, ExpressionException {
Simulation simulation = simTask.getSimulation();
DataProcessingInstructions dpi = simulation.getDataProcessingInstructions();
if (dpi == null) {
printWriter.println("DATA_PROCESSOR_BEGIN " + DataProcessingInstructions.ROI_TIME_SERIES);
printWriter.println("DATA_PROCESSOR_END");
printWriter.println();
return;
}
FieldDataIdentifierSpec fdis = dpi.getSampleImageFieldData(simulation.getVersion().getOwner());
if (fdis == null) {
throw new DataAccessException("Can't find sample image in data processing instructions");
}
String secondarySimDataDir = PropertyLoader.getProperty(PropertyLoader.secondarySimDataDirProperty, null);
DataSetControllerImpl dsci = new DataSetControllerImpl(new NullSessionLog(),null,userDirectory.getParentFile(),secondarySimDataDir == null ? null : new File(secondarySimDataDir));
CartesianMesh origMesh = dsci.getMesh(fdis.getExternalDataIdentifier());
SimDataBlock simDataBlock = dsci.getSimDataBlock(null,fdis.getExternalDataIdentifier(), fdis.getFieldFuncArgs().getVariableName(), fdis.getFieldFuncArgs().getTime().evaluateConstant());
VariableType varType = fdis.getFieldFuncArgs().getVariableType();
VariableType dataVarType = simDataBlock.getVariableType();
if (!varType.equals(VariableType.UNKNOWN) && !varType.equals(dataVarType)) {
throw new IllegalArgumentException("field function variable type (" + varType.getTypeName() + ") doesn't match real variable type (" + dataVarType.getTypeName() + ")");
}
double[] origData = simDataBlock.getData();
String filename = SimulationJob.createSimulationJobID(Simulation.createSimulationID(simulation.getKey()), simulationJob.getJobIndex()) + FieldDataIdentifierSpec.getDefaultFieldDataFileNameForSimulation(fdis.getFieldFuncArgs());
File fdatFile = new File(userDirectory, filename);
DataSet.writeNew(fdatFile,
new String[] {fdis.getFieldFuncArgs().getVariableName()},
new VariableType[]{simDataBlock.getVariableType()},
new ISize(origMesh.getSizeX(),origMesh.getSizeY(),origMesh.getSizeZ()),
new double[][]{origData});
printWriter.println("DATA_PROCESSOR_BEGIN " + dpi.getScriptName());
printWriter.println(dpi.getScriptInput());
printWriter.println("SampleImageFile " + fdis.getFieldFuncArgs().getVariableName() + " " + fdis.getFieldFuncArgs().getTime().infix() + " " + fdatFile);
printWriter.println("DATA_PROCESSOR_END");
printWriter.println();
}
*/
/**
*# fast system dimension num_dependents
*FAST_SYSTEM_BEGIN 2 2
*INDEPENDENT_VARIALBES rf r
*DEPENDENT_VARIALBES rB rfB
*
*PSEUDO_CONSTANT_BEGIN
*__C0 (rfB + rf);
*__C1 (r + rB);
*PSEUDO_CONSTANT_END
*
*FAST_RATE_BEGIN
* - ((0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0) * rf) - (0.1 * ( - rf + __C0)));
*((0.02 * r * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) - (0.1 * ( - r + __C1)));
*FAST_RATE_END
*
*FAST_DEPENDENCY_BEGIN
*rB ( - r + __C1);
*rfB ( - rf + __C0);
*FAST_DEPENDENCY_END
*
*JACOBIAN_BEGIN
* - (0.1 + (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf) + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)));
* - (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf);
*(0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))));
*(0.1 + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) + (0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))))));
*JACOBIAN_END
*
*FAST_SYSTEM_END
* @throws ExpressionException
* @throws MathException
*/
private void writeFastSystem(SubDomain subDomain) throws MathException, ExpressionException {
VariableDomain variableDomain = (subDomain instanceof CompartmentSubDomain) ? VariableDomain.VARIABLEDOMAIN_VOLUME : VariableDomain.VARIABLEDOMAIN_MEMBRANE;
FastSystem fastSystem = subDomain.getFastSystem();
if (fastSystem == null) {
return;
}
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, simTask.getSimulationJob().getSimulationSymbolTable());
int numIndep = fs_analyzer.getNumIndependentVariables();
int numDep = fs_analyzer.getNumDependentVariables();
int numPseudo = fs_analyzer.getNumPseudoConstants();
printWriter.println("# fast system dimension num_dependents");
printWriter.println("FAST_SYSTEM_BEGIN " + numIndep + " " + numDep);
if (numIndep != 0) {
printWriter.print("INDEPENDENT_VARIALBES ");
Enumeration<Variable> enum1 = fs_analyzer.getIndependentVariables();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
printWriter.print(var.getName() + " ");
}
printWriter.println();
}
if (numDep != 0) {
printWriter.print("DEPENDENT_VARIALBES ");
Enumeration<Variable> enum1 = fs_analyzer.getDependentVariables();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
printWriter.print(var.getName() + " ");
}
printWriter.println();
}
printWriter.println();
if (numPseudo != 0) {
printWriter.println("PSEUDO_CONSTANT_BEGIN");
Enumeration<PseudoConstant> enum1 = fs_analyzer.getPseudoConstants();
while (enum1.hasMoreElements()) {
PseudoConstant pc = enum1.nextElement();
printWriter.println(pc.getName() + " " + subsituteExpression(pc.getPseudoExpression(), fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("PSEUDO_CONSTANT_END");
printWriter.println();
}
if (numIndep != 0) {
printWriter.println("FAST_RATE_BEGIN");
Enumeration<Expression> enum1 = fs_analyzer.getFastRateExpressions();
while (enum1.hasMoreElements()) {
Expression exp = enum1.nextElement();
printWriter.println(subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("FAST_RATE_END");
printWriter.println();
}
if (numDep != 0) {
printWriter.println("FAST_DEPENDENCY_BEGIN");
Enumeration<Expression> enum_exp = fs_analyzer.getDependencyExps();
Enumeration<Variable> enum_var = fs_analyzer.getDependentVariables();
while (enum_exp.hasMoreElements()) {
Expression exp = enum_exp.nextElement();
Variable depVar = enum_var.nextElement();
printWriter.println(depVar.getName() + " " + subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("FAST_DEPENDENCY_END");
printWriter.println();
}
if (numIndep != 0) {
printWriter.println("JACOBIAN_BEGIN");
Enumeration<Expression> enum_fre = fs_analyzer.getFastRateExpressions();
while (enum_fre.hasMoreElements()) {
Expression fre = enum_fre.nextElement();
Enumeration<Variable> enum_var = fs_analyzer.getIndependentVariables();
while (enum_var.hasMoreElements()) {
Variable var = enum_var.nextElement();
Expression exp = subsituteExpression(fre, fs_analyzer, variableDomain).flatten();
Expression differential = exp.differentiate(var.getName());
printWriter.println(subsituteExpression(differential, fs_analyzer, variableDomain).infix() + ";");
}
}
printWriter.println("JACOBIAN_END");
printWriter.println();
}
printWriter.println("FAST_SYSTEM_END");
printWriter.println();
}
Aggregations