use of cbit.vcell.matrix.MatrixException in project vcell by virtualcell.
the class MassActionSolver method solveMassAction.
public static MassActionFunction solveMassAction(Parameter optionalForwardRateParameter, Parameter optionalReverseRateParameter, Expression orgExp, ReactionStep rs) throws ExpressionException, ModelException, DivideByZeroException {
MassActionFunction maFunc = new MassActionFunction();
// get reactants, products, overlaps, non-overlap reactants and non-overlap products
ArrayList<ReactionParticipant> reactants = new ArrayList<ReactionParticipant>();
ArrayList<ReactionParticipant> products = new ArrayList<ReactionParticipant>();
ReactionParticipant[] rp = rs.getReactionParticipants();
// should use this one to compare functional equavalent since this duplicatedExp has all params substituted.
Expression duplicatedExp = substituteParameters(orgExp, false);
// separate the reactants and products, fluxes, catalysts
String rxnName = rs.getName();
// reaction with membrane current can not be transformed to mass action
if (rs.getPhysicsOptions() == ReactionStep.PHYSICS_MOLECULAR_AND_ELECTRICAL || rs.getPhysicsOptions() == ReactionStep.PHYSICS_ELECTRICAL_ONLY) {
throw new ModelException("Kinetics of reaction " + rxnName + " has membrane current. It can not be automatically transformed to Mass Action kinetics.");
}
for (int i = 0; i < rp.length; i++) {
if (rp[i] instanceof Reactant) {
reactants.add(rp[i]);
} else if (rp[i] instanceof Product) {
products.add(rp[i]);
} else if (rp[i] instanceof Catalyst) {
String catalystName = rp[i].getSpeciesContext().getName();
// only if duplictedExp is not a non-linear function of catalystName.
if (duplicatedExp.hasSymbol(catalystName)) {
// Only when catalyst appears in ReactionRate, we add catalyst to both reactant and product
// the stoichiometry should be set to 1.
ReactionParticipant catalystRP = new ReactionParticipant(null, rs, rp[i].getSpeciesContext(), 1) {
public boolean compareEqual(Matchable obj) {
ReactionParticipant rp = (ReactionParticipant) obj;
if (rp == null) {
return false;
}
if (!Compare.isEqual(getSpecies(), rp.getSpecies())) {
return false;
}
if (!Compare.isEqual(getStructure(), rp.getStructure())) {
return false;
}
if (getStoichiometry() != rp.getStoichiometry()) {
return false;
}
return true;
}
@Override
public void writeTokens(PrintWriter pw) {
}
@Override
public void fromTokens(CommentStringTokenizer tokens, Model model) throws Exception {
}
};
products.add(catalystRP);
reactants.add(catalystRP);
}
}
}
/**
* The code below is going to solve reaction with kinetics that are NOT Massaction. Or Massaction with catalysts involved.
*/
//
// 2x2 rational matrix
//
// lets define
// J() is the substituted total rate expression
// P() as the theoretical "product" term with Kf = 1
// R() as the theoretical "reactant" term with Kr = 1
//
// Kf * R([1 1 1]) - Kr * P([1 1 1]) = J([1 1 1])
// Kf * R([2 3 4]) - Kr * P([2 3 4]) = J([2 3 4])
//
// in matrix form,
//
// | |
// | R([1 1 1]) -P([1 1 1]) J([1 1 1]) |
// | R([2 3 4]) -P([2 3 4]) J([2 3 4]) |
// | |
//
//
// example: Kf*A*B^2*C - Kr*C*A
// J() = forwardRate*43/Kabc*A*B^2*C - (myC*5-2)*C*A
// P() = C*A
// R() = A*B^2*C
//
// | |
// | R([1 1 1]) -P([1 1 1]) J([1 1 1]) |
// | R([2 3 4]) -P([2 3 4]) J([2 3 4]) |
// | |
//
// | |
// | R([1 1 1])*R([2 3 4]) -P([1 1 1])*R([2 3 4]) J([1 1 1])*R([2 3 4]) |
// | R([2 3 4])*R([1 1 1]) -P([2 3 4])*R([1 1 1]) J([2 3 4])*R([1 1 1]) |
// | |
//
//
// | |
// | R([1 1 1])*R([2 3 4]) -P([1 1 1])*R([2 3 4]) J([1 1 1])*R([2 3 4]) |
// | 0 -P([2 3 4])*R([1 1 1])+P([1 1 1])*R([2 3 4]) J([2 3 4])*R([1 1 1])-J([1 1 1])*R([2 3 4]) |
// | |
//
//
//
//
Expression forwardExp = null;
Expression reverseExp = null;
Expression R_exp = new Expression(1);
if (reactants.size() > 0) {
R_exp = new Expression(1.0);
for (ReactionParticipant reactant : reactants) {
R_exp = Expression.mult(R_exp, Expression.power(new Expression(reactant.getName()), new Expression(reactant.getStoichiometry())));
}
}
Expression P_exp = new Expression(1);
if (products.size() > 0) {
P_exp = new Expression(1.0);
for (ReactionParticipant product : products) {
P_exp = Expression.mult(P_exp, Expression.power(new Expression(product.getName()), new Expression(product.getStoichiometry())));
}
}
HashSet<String> reactionParticipantNames = new HashSet<String>();
for (ReactionParticipant reactionParticipant : rs.getReactionParticipants()) {
reactionParticipantNames.add(reactionParticipant.getName());
}
Expression R_1 = new Expression(R_exp);
Expression R_2 = new Expression(R_exp);
Expression P_1 = new Expression(P_exp);
Expression P_2 = new Expression(P_exp);
Expression J_1 = new Expression(duplicatedExp);
Expression J_2 = new Expression(duplicatedExp);
int index = 0;
for (String rpName : reactionParticipantNames) {
R_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
R_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
P_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
P_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
J_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
J_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
index++;
}
R_1 = R_1.flatten();
R_2 = R_2.flatten();
P_1 = P_1.flatten();
P_2 = P_2.flatten();
J_1 = J_1.flatten();
J_2 = J_2.flatten();
// e.g. A+B <-> A+B, A+2B <-> A+2B
if (ExpressionUtils.functionallyEquivalent(R_exp, P_exp, false, 1e-8, 1e-8)) {
throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Identical reactants and products not supported for stochastic models."));
}
if (R_exp.compareEqual(new Expression(1)) && P_exp.compareEqual(new Expression(1))) {
// no reactants, no products ... nothing to do
forwardExp = null;
reverseExp = null;
} else {
// both reactants and products
RationalExpMatrix matrix = new RationalExpMatrix(2, 3);
matrix.set_elem(0, 0, RationalExpUtils.getRationalExp(R_1, true));
matrix.set_elem(0, 1, RationalExpUtils.getRationalExp(Expression.negate(P_1), true));
matrix.set_elem(0, 2, RationalExpUtils.getRationalExp(J_1, true));
matrix.set_elem(1, 0, RationalExpUtils.getRationalExp(R_2, true));
matrix.set_elem(1, 1, RationalExpUtils.getRationalExp(Expression.negate(P_2), true));
matrix.set_elem(1, 2, RationalExpUtils.getRationalExp(J_2, true));
RationalExp[] solution;
try {
// matrix.show();
solution = matrix.solveLinearExpressions();
// solution[0] is forward rate.
solution[0] = solution[0].simplify();
// solution[1] is reverse rate.
solution[1] = solution[1].simplify();
} catch (MatrixException e) {
e.printStackTrace();
throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "MatrixException: " + e.getMessage()));
}
forwardExp = new Expression(solution[0].infixString());
reverseExp = new Expression(solution[1].infixString());
// for massAction, if there is no reactant in the forward rate, the forwardExp should be set to null to avoid writing jump process for the forward reaction.
if (R_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardExp = null;
}
// for massAction, if there is no product in the reverse rate, the reverseExp should be set to null to avoid writing jump process for the reverse reaction.
if (P_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
reverseExp = null;
}
}
if (forwardExp != null) {
forwardExp.bindExpression(rs);
}
if (reverseExp != null) {
reverseExp.bindExpression(rs);
}
// Reconstruct the rate based on the extracted forward rate and reverse rate. If the reconstructed rate is not equivalent to the original rate,
// it means the original rate is not in the form of Kf*r1^n1*r2^n2-Kr*p1^m1*p2^m2.
Expression constructedExp = reconstructedRate(forwardExp, reverseExp, reactants, products, rs.getNameScope());
Expression orgExp_withoutCatalyst = removeCatalystFromExp(duplicatedExp, rs);
Expression constructedExp_withoutCatalyst = removeCatalystFromExp(constructedExp, rs);
if (!ExpressionUtils.functionallyEquivalent(orgExp_withoutCatalyst, constructedExp_withoutCatalyst, false, 1e-8, 1e-8)) {
throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Mathmatical form incompatible with mass action."));
}
// check if forward rate constant and reverse rate constant both can be evaluated to constants(numbers) after substituting all parameters.
if (forwardExp != null) {
Expression forwardExpCopy = new Expression(forwardExp);
try {
substituteParameters(forwardExpCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in forward rate '" + forwardExp.infix() + "', " + e.getMessage()));
}
//
if (optionalForwardRateParameter != null) {
Expression forwardRateParameterCopy = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
try {
substituteParameters(forwardRateParameterCopy, true).evaluateConstant();
if (forwardExpCopy.compareEqual(forwardRateParameterCopy)) {
forwardExp = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
}
} catch (ExpressionException e) {
// not expecting a problem because forwardExpCopy didn't have a problem, but in any case it is ok to swallow this exception
e.printStackTrace(System.out);
}
}
}
if (reverseExp != null) {
Expression reverseExpCopy = new Expression(reverseExp);
try {
substituteParameters(reverseExpCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in reverse rate '" + reverseExp.infix() + "', " + e.getMessage()));
}
//
if (optionalReverseRateParameter != null) {
Expression reverseRateParameterCopy = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
try {
substituteParameters(reverseRateParameterCopy, true).evaluateConstant();
if (reverseExpCopy.compareEqual(reverseRateParameterCopy)) {
reverseExp = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
}
} catch (ExpressionException e) {
// not expecting a problem because reverseExpCopy didn't have a problem, but in any case it is ok to swallow this exception
e.printStackTrace(System.out);
}
}
}
maFunc.setForwardRate(forwardExp);
maFunc.setReverseRate(reverseExp);
maFunc.setReactants(reactants);
maFunc.setProducts(products);
return maFunc;
}
use of cbit.vcell.matrix.MatrixException in project vcell by virtualcell.
the class FastSystemAnalyzer method refreshInvarianceMatrix.
/**
*/
private void refreshInvarianceMatrix() throws MathException, ExpressionException {
//
// the invariance's are expressed in matrix form
//
// |a a a a a -1 0 0| |x1| |0|
// |a a a a a 0 -1 0| * |x2| = |0|
// |a a a a a 0 0 -1| |x3| |0|
// |x4|
// |x5|
// |c1|
// |c2|
// |c3|
//
Variable[] vars = new Variable[fastVarList.size()];
fastVarList.copyInto(vars);
int numVars = fastVarList.size();
int rows = fastSystem.getNumFastInvariants();
int cols = numVars + fastSystem.getNumFastInvariants();
RationalExpMatrix matrix = new RationalExpMatrix(rows, cols);
Enumeration<FastInvariant> fastInvariantsEnum = fastSystem.getFastInvariants();
for (int i = 0; i < rows && fastInvariantsEnum.hasMoreElements(); i++) {
FastInvariant fi = (FastInvariant) fastInvariantsEnum.nextElement();
Expression function = fi.getFunction();
for (int j = 0; j < numVars; j++) {
Variable var = (Variable) fastVarList.elementAt(j);
Expression exp = function.differentiate(var.getName());
exp.bindExpression(null);
exp = exp.flatten();
RationalExp coeffRationalExp = RationalExpUtils.getRationalExp(exp);
matrix.set_elem(i, j, coeffRationalExp);
}
matrix.set_elem(i, numVars + i, -1);
}
// Print
System.out.println("origMatrix");
matrix.show();
//
// gaussian elimination on the matrix give the following representation
// note that some column pivoting (variable re-ordering) is sometimes required to
// determine N-r dependent vars
//
// |10i0iccc|
// |01i0iccc| where (c)'s are the coefficients for constants of invariances
// |00i1iccc| (i)'s are the coefficients for dependent vars in terms of independent vars
//
// Print
System.out.println("reducedMatrix");
if (rows > 0) {
try {
matrix.gaussianElimination(new RationalExpMatrix(rows, rows));
} catch (MatrixException e) {
e.printStackTrace(System.out);
throw new MathException(e.getMessage());
}
}
matrix.show();
for (int i = 0; i < vars.length; i++) {
System.out.print(vars[i].getName() + " ");
}
System.out.println("");
//
for (int i = 0; i < rows; i++) {
for (int j = 0; j < rows; j++) {
RationalExp rexp = matrix.get(i, j).simplify();
matrix.set_elem(i, j, rexp);
}
}
for (int i = 0; i < rows; i++) {
//
if (!matrix.get(i, i).isConstant() || matrix.get(i, i).getConstant().doubleValue() != 1) {
for (int j = i + 1; j < numVars; j++) {
if (matrix.get(i, j).isConstant() && matrix.get(i, j).getConstant().doubleValue() == 1.0) {
for (int ii = 0; ii < rows; ii++) {
RationalExp temp = matrix.get(ii, i);
matrix.set_elem(ii, i, matrix.get(ii, j));
matrix.set_elem(ii, j, temp);
}
Variable tempVar = vars[i];
vars[i] = vars[j];
vars[j] = tempVar;
break;
}
}
}
}
// Print
for (int i = 0; i < vars.length; i++) {
System.out.print(vars[i].getName() + " ");
}
System.out.println("");
matrix.show();
//
// separate into dependent and indepent variables, and chop off identity matrix (left N-r columns)
//
// T |iiccc| T
// [x1 x2 x4] = -1 * |iiccc| * [x3 x5 c1 c2 c3]
// |iiccc|
//
//
int numInvariants = fastSystem.getNumFastInvariants();
dependentVarList.removeAllElements();
for (int i = 0; i < numInvariants; i++) {
dependentVarList.addElement(vars[i]);
}
independentVarList.removeAllElements();
for (int i = numInvariants; i < vars.length; i++) {
independentVarList.addElement(vars[i]);
}
int new_cols = independentVarList.size() + numInvariants;
dependencyMatrix = new RationalExpMatrix(rows, new_cols);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < new_cols; j++) {
RationalExp rexp = matrix.get(i, j + dependentVarList.size()).simplify().minus();
dependencyMatrix.set_elem(i, j, rexp);
}
}
// Print
System.out.println("\n\nDEPENDENCY MATRIX");
dependencyMatrix.show();
System.out.print("dependent vars: ");
for (int i = 0; i < dependentVarList.size(); i++) {
System.out.print(((Variable) dependentVarList.elementAt(i)).getName() + " ");
}
System.out.println("");
System.out.print("independent vars: ");
for (int i = 0; i < independentVarList.size(); i++) {
System.out.print(((Variable) independentVarList.elementAt(i)).getName() + " ");
}
System.out.println("");
}
use of cbit.vcell.matrix.MatrixException 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.matrix.MatrixException in project vcell by virtualcell.
the class ModelOptimizationMapping method getRemappedReferenceData.
/**
* Gets the constraintData property (cbit.vcell.opt.ConstraintData) value.
* @return The constraintData property value.
* @see #setConstraintData
*/
private ReferenceData getRemappedReferenceData(MathMapping mathMapping) throws MappingException {
if (modelOptimizationSpec.getReferenceData() == null) {
return null;
}
//
// make sure time is mapped
//
ReferenceData refData = modelOptimizationSpec.getReferenceData();
ReferenceDataMappingSpec[] refDataMappingSpecs = modelOptimizationSpec.getReferenceDataMappingSpecs();
RowColumnResultSet rowColResultSet = new RowColumnResultSet();
Vector<SymbolTableEntry> modelObjectList = new Vector<SymbolTableEntry>();
Vector<double[]> dataList = new Vector<double[]>();
//
// find bound columns, (time is always mapped to the first column)
//
int mappedColumnCount = 0;
for (int i = 0; i < refDataMappingSpecs.length; i++) {
SymbolTableEntry modelObject = refDataMappingSpecs[i].getModelObject();
if (modelObject != null) {
int mappedColumnIndex = mappedColumnCount;
if (modelObject instanceof Model.ReservedSymbol && ((ReservedSymbol) modelObject).isTime()) {
mappedColumnIndex = 0;
}
String origRefDataColumnName = refDataMappingSpecs[i].getReferenceDataColumnName();
int origRefDataColumnIndex = refData.findColumn(origRefDataColumnName);
if (origRefDataColumnIndex < 0) {
throw new RuntimeException("reference data column named '" + origRefDataColumnName + "' not found");
}
double[] columnData = refData.getDataByColumn(origRefDataColumnIndex);
if (modelObjectList.contains(modelObject)) {
throw new RuntimeException("multiple reference data columns mapped to same model object '" + modelObject.getName() + "'");
}
modelObjectList.insertElementAt(modelObject, mappedColumnIndex);
dataList.insertElementAt(columnData, mappedColumnIndex);
mappedColumnCount++;
}
}
//
if (modelObjectList.size() == 0) {
throw new RuntimeException("reference data was not associated with model");
}
if (modelObjectList.size() == 1) {
throw new RuntimeException("reference data was not associated with model, must map time and at least one other column");
}
boolean bFoundTimeVar = false;
for (SymbolTableEntry ste : modelObjectList) {
if (ste instanceof Model.ReservedSymbol && ((ReservedSymbol) ste).isTime()) {
bFoundTimeVar = true;
break;
}
}
if (!bFoundTimeVar) {
throw new RuntimeException("must map time column of reference data to model");
}
//
for (int i = 0; i < modelObjectList.size(); i++) {
SymbolTableEntry modelObject = (SymbolTableEntry) modelObjectList.elementAt(i);
try {
// Find by name because MathSybolMapping has different 'objects' than refDataMapping 'objects'
Variable variable = mathMapping.getMathSymbolMapping().findVariableByName(modelObject.getName());
if (variable != null) {
String symbol = variable.getName();
rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
} else if (modelObject instanceof Model.ReservedSymbol && ((Model.ReservedSymbol) modelObject).isTime()) {
Model.ReservedSymbol time = (Model.ReservedSymbol) modelObject;
String symbol = time.getName();
rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
}
} catch (MathException | MatrixException | ExpressionException | ModelException e) {
e.printStackTrace();
throw new MappingException(e.getMessage(), e);
}
}
//
// populate data columns (time and rest)
//
double[] weights = new double[rowColResultSet.getColumnDescriptionsCount()];
weights[0] = 1.0;
int numRows = ((double[]) dataList.elementAt(0)).length;
int numColumns = modelObjectList.size();
for (int j = 0; j < numRows; j++) {
double[] row = new double[numColumns];
for (int i = 0; i < numColumns; i++) {
row[i] = ((double[]) dataList.elementAt(i))[j];
if (i > 0) {
weights[i] += row[i] * row[i];
}
}
rowColResultSet.addRow(row);
}
for (int i = 0; i < numColumns; i++) {
if (weights[i] == 0) {
weights[i] = 1;
} else {
weights[i] = 1 / weights[i];
}
}
SimpleReferenceData remappedRefData = new SimpleReferenceData(rowColResultSet, weights);
return remappedRefData;
}
use of cbit.vcell.matrix.MatrixException in project vcell by virtualcell.
the class MathVerifier method testMathGeneration.
public MathGenerationResults testMathGeneration(KeyValue simContextKey) throws SQLException, ObjectNotFoundException, DataAccessException, XmlParseException, MappingException, MathException, MatrixException, ExpressionException, ModelException, PropertyVetoException {
User adminUser = new User(PropertyLoader.ADMINISTRATOR_ACCOUNT, new org.vcell.util.document.KeyValue(PropertyLoader.ADMINISTRATOR_ID));
if (lg.isTraceEnabled())
lg.trace("Testing SimContext with key '" + simContextKey + "'");
// get biomodel refs
java.sql.Connection con = null;
java.sql.Statement stmt = null;
con = conFactory.getConnection(new Object());
cbit.vcell.modeldb.BioModelSimContextLinkTable bmscTable = cbit.vcell.modeldb.BioModelSimContextLinkTable.table;
cbit.vcell.modeldb.BioModelTable bmTable = cbit.vcell.modeldb.BioModelTable.table;
cbit.vcell.modeldb.UserTable userTable = cbit.vcell.modeldb.UserTable.table;
String sql = "SELECT " + bmscTable.bioModelRef.getQualifiedColName() + "," + bmTable.ownerRef.getQualifiedColName() + "," + userTable.userid.getQualifiedColName() + " FROM " + bmscTable.getTableName() + "," + bmTable.getTableName() + "," + userTable.getTableName() + " WHERE " + bmscTable.simContextRef.getQualifiedColName() + " = " + simContextKey + " AND " + bmTable.id.getQualifiedColName() + " = " + bmscTable.bioModelRef.getQualifiedColName() + " AND " + bmTable.ownerRef.getQualifiedColName() + " = " + userTable.id.getQualifiedColName();
ArrayList<KeyValue> bioModelKeys = new ArrayList<KeyValue>();
stmt = con.createStatement();
User owner = null;
try {
ResultSet rset = stmt.executeQuery(sql);
while (rset.next()) {
KeyValue key = new KeyValue(rset.getBigDecimal(bmscTable.bioModelRef.getUnqualifiedColName()));
bioModelKeys.add(key);
KeyValue ownerRef = new KeyValue(rset.getBigDecimal(bmTable.ownerRef.getUnqualifiedColName()));
String userid = rset.getString(userTable.userid.getUnqualifiedColName());
owner = new User(userid, ownerRef);
}
} finally {
if (stmt != null) {
stmt.close();
}
con.close();
}
// use the first biomodel...
if (bioModelKeys.size() == 0) {
throw new RuntimeException("zombie simContext ... no biomodels");
}
BioModelInfo bioModelInfo = dbServerImpl.getBioModelInfo(owner, bioModelKeys.get(0));
//
// read in the BioModel from the database
//
BigString bioModelXML = dbServerImpl.getBioModelXML(owner, bioModelInfo.getVersion().getVersionKey());
BioModel bioModelFromDB = XmlHelper.XMLToBioModel(new XMLSource(bioModelXML.toString()));
BioModel bioModelNewMath = XmlHelper.XMLToBioModel(new XMLSource(bioModelXML.toString()));
bioModelFromDB.refreshDependencies();
bioModelNewMath.refreshDependencies();
//
// get all Simulations for this model
//
Simulation[] modelSimsFromDB = bioModelFromDB.getSimulations();
//
// ---> only for the SimContext we started with...
// recompute mathDescription, and verify it is equivalent
// then check each associated simulation to ensure math overrides are applied in an equivalent manner also.
//
SimulationContext[] simContextsFromDB = bioModelFromDB.getSimulationContexts();
SimulationContext[] simContextsNewMath = bioModelNewMath.getSimulationContexts();
SimulationContext simContextFromDB = null;
SimulationContext simContextNewMath = null;
for (int k = 0; k < simContextsFromDB.length; k++) {
// find it...
if (simContextsFromDB[k].getKey().equals(simContextKey)) {
simContextFromDB = simContextsFromDB[k];
simContextNewMath = simContextsNewMath[k];
break;
}
}
if (simContextFromDB == null) {
throw new RuntimeException("BioModel referred to by this SimContext does not contain this SimContext");
} else {
MathDescription origMathDesc = simContextFromDB.getMathDescription();
//
try {
if (simContextNewMath.getGeometry().getDimension() > 0 && simContextNewMath.getGeometry().getGeometrySurfaceDescription().getGeometricRegions() == null) {
simContextNewMath.getGeometry().getGeometrySurfaceDescription().updateAll();
}
} catch (Exception e) {
e.printStackTrace(System.out);
}
//
// updated mathDescription loaded into copy of bioModel, then test for equivalence.
//
cbit.vcell.mapping.MathMapping mathMapping = simContextNewMath.createNewMathMapping();
MathDescription mathDesc_latest = mathMapping.getMathDescription();
MathMapping_4_8 mathMapping_4_8 = new MathMapping_4_8(simContextNewMath);
MathDescription mathDesc_4_8 = mathMapping_4_8.getMathDescription();
String issueString = null;
org.vcell.util.Issue[] issues = mathMapping.getIssues();
if (issues != null && issues.length > 0) {
StringBuffer buffer = new StringBuffer("Issues(" + issues.length + "):\n");
for (int l = 0; l < issues.length; l++) {
buffer.append(" <<" + issues[l].toString() + ">>\n");
}
issueString = buffer.toString();
}
simContextNewMath.setMathDescription(mathDesc_latest);
MathCompareResults mathCompareResults_latest = MathDescription.testEquivalency(SimulationSymbolTable.createMathSymbolTableFactory(), origMathDesc, mathDesc_latest);
MathCompareResults mathCompareResults_4_8 = null;
try {
mathCompareResults_4_8 = MathDescription.testEquivalency(SimulationSymbolTable.createMathSymbolTableFactory(), origMathDesc, mathDesc_4_8);
} catch (Exception e) {
e.printStackTrace(System.out);
mathCompareResults_4_8 = new MathCompareResults(Decision.MathDifferent_FAILURE_UNKNOWN, e.getMessage());
}
return new MathGenerationResults(bioModelFromDB, simContextFromDB, origMathDesc, mathDesc_latest, mathCompareResults_latest, mathDesc_4_8, mathCompareResults_4_8);
}
}
Aggregations