use of cbit.vcell.math.FastInvariant in project vcell by virtualcell.
the class MathMapping_4_8 method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// All sizes must be set for new ODE models and ratios must be set for old ones.
simContext.checkValidity();
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
Model model = simContext.getModel();
StructureTopology structTopology = model.getStructureTopology();
//
// verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
//
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
if (sm == null || (sm instanceof FeatureMapping && getSubVolume((FeatureMapping) sm) == null)) {
throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subdomain");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
if (volFractExp != null) {
try {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + structTopology.getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0");
}
} catch (ExpressionException e) {
}
}
}
}
SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
throw new MappingException("geometry subdomain '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
// deals with model parameters
Hashtable<VolVariable, EventAssignmentInitParameter> eventVolVarHash = new Hashtable<VolVariable, EventAssignmentInitParameter>();
ModelParameter[] modelParameters = model.getModelParameters();
if (simContext.getGeometry().getDimension() == 0) {
//
// global parameters from model (that presently are constants)
//
BioEvent[] bioEvents = simContext.getBioEvents();
ArrayList<SymbolTableEntry> eventAssignTargets = new ArrayList<SymbolTableEntry>();
if (bioEvents != null && bioEvents.length > 0) {
for (BioEvent be : bioEvents) {
for (EventAssignment ea : be.getEventAssignments()) {
if (!eventAssignTargets.contains(ea.getTarget())) {
eventAssignTargets.add(ea.getTarget());
}
}
}
}
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
if (eventAssignTargets.contains(modelParameters[j])) {
EventAssignmentInitParameter eap = null;
try {
eap = addEventAssignmentInitParameter(modelParameters[j].getName(), modelParameters[j].getExpression(), PARAMETER_ROLE_EVENTASSIGN_INITCONDN, modelParameters[j].getUnitDefinition());
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
// varHash.addVariable(newFunctionOrConstant(getMathSymbol(eap, null), modelParamExpr));
VolVariable volVar = new VolVariable(modelParameters[j].getName(), nullDomain);
varHash.addVariable(volVar);
eventVolVarHash.put(volVar, eap);
} else {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), modelParamExpr));
}
}
} else {
//
for (int pass = 0; pass < 2; pass++) {
for (int j = 0; j < modelParameters.length; j++) {
Hashtable<String, Expression> structMappingVariantsHash = new Hashtable<String, Expression>();
for (int k = 0; k < structureMappings.length; k++) {
String paramVariantName = null;
Expression paramVariantExpr = null;
if (modelParameters[j].getExpression().getSymbols() == null) {
paramVariantName = modelParameters[j].getName();
paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
} else {
paramVariantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
// if the expression has symbols that do not belong in that structureMapping, do not create the variant.
Expression exp1 = modelParameters[j].getExpression();
Expression flattenedModelParamExpr = substituteGlobalParameters(exp1);
String[] symbols = flattenedModelParamExpr.getSymbols();
boolean bValid = true;
Structure sm_struct = structureMappings[k].getStructure();
if (symbols != null) {
for (int ii = 0; ii < symbols.length; ii++) {
SpeciesContext sc = model.getSpeciesContext(symbols[ii]);
if (sc != null) {
// symbol[ii] is a speciesContext, check its structure with structureMapping[k].structure. If they are the same or
// if it is the adjacent membrane(s), allow variant expression to be created. Else, continue.
Structure sp_struct = sc.getStructure();
if (sp_struct.compareEqual(sm_struct)) {
bValid = bValid && true;
} else {
// if the 2 structures are not the same, are they adjacent? then 'bValid' is true, else false.
if ((sm_struct instanceof Feature) && (sp_struct instanceof Membrane)) {
Feature sm_feature = (Feature) sm_struct;
Membrane sp_mem = (Membrane) sp_struct;
if (sp_mem.compareEqual(structTopology.getParentStructure(sm_feature)) || (structTopology.getInsideFeature(sp_mem).compareEqual(sm_feature) || structTopology.getOutsideFeature(sp_mem).compareEqual(sm_feature))) {
bValid = bValid && true;
} else {
bValid = bValid && false;
break;
}
} else if ((sm_struct instanceof Membrane) && (sp_struct instanceof Feature)) {
Feature sp_feature = (Feature) sp_struct;
Membrane sm_mem = (Membrane) sm_struct;
if (sm_mem.compareEqual(structTopology.getParentStructure(sp_feature)) || (structTopology.getInsideFeature(sm_mem).compareEqual(sp_feature) || structTopology.getOutsideFeature(sm_mem).compareEqual(sp_feature))) {
bValid = bValid && true;
} else {
bValid = bValid && false;
break;
}
} else {
bValid = bValid && false;
break;
}
}
}
}
}
if (bValid) {
if (pass == 0) {
paramVariantExpr = new Expression("VCELL_TEMPORARY_EXPRESSION_PLACEHOLDER");
} else {
paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), structureMappings[k]);
}
}
}
if (paramVariantExpr != null) {
structMappingVariantsHash.put(paramVariantName, paramVariantExpr);
}
}
globalParamVariantsHash.put(modelParameters[j], structMappingVariantsHash);
}
}
//
for (int j = 0; j < modelParameters.length; j++) {
if (modelParameters[j].getExpression().getSymbols() == null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null)));
} else {
Hashtable<String, Expression> smVariantsHash = globalParamVariantsHash.get(modelParameters[j]);
for (int k = 0; k < structureMappings.length; k++) {
String variantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
Expression variantExpr = smVariantsHash.get(variantName);
if (variantExpr != null) {
varHash.addVariable(newFunctionOrConstant(variantName, variantExpr));
}
}
}
}
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = simContext.getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(simContext.getName() + "_generated");
}
//
// volume variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof VolVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof VolVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
//
// membrane variables
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof MemVariable) {
varHash.addVariable(scm.getVariable());
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
// only calculate potential if at least one MembraneMapping has CalculateVoltage == true
//
boolean bCalculatePotential = false;
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
bCalculatePotential = true;
}
}
}
// (simContext.getGeometry().getDimension() == 0);
potentialMapping = new PotentialMapping(simContext, this);
potentialMapping.computeMath();
if (bCalculatePotential) {
//
// copy functions for currents and constants for capacitances
//
ElectricalDevice[] devices = potentialMapping.getElectricalDevices();
for (int j = 0; j < devices.length; j++) {
if (devices[j] instanceof MembraneElectricalDevice) {
MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) devices[j];
MembraneMapping memMapping = membraneElectricalDevice.getMembraneMapping();
Parameter specificCapacitanceParm = memMapping.getParameterFromRole(MembraneMapping.ROLE_SpecificCapacitance);
varHash.addVariable(new Constant(getMathSymbol(specificCapacitanceParm, memMapping), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping)));
ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
memMapping.getCalculateVoltage()) {
Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
if (capacitanceParm != null && capacitanceParm.getExpression() != null && memMapping.getCalculateVoltage()) {
StructureMappingParameter sizeParameter = membraneElectricalDevice.getMembraneMapping().getSizeParameter();
if (simContext.getGeometry().getDimension() == 0 && (sizeParameter.getExpression() == null || sizeParameter.getExpression().isZero())) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
} else {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
}
//
if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
// is Voltage Independent?
StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping)));
} else //
// membrane forced potential
//
{
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping)));
}
} else if (devices[j] instanceof CurrentClampElectricalDevice) {
CurrentClampElectricalDevice currentClampDevice = (CurrentClampElectricalDevice) devices[j];
// total current = current source (no capacitance)
Parameter totalCurrentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TotalCurrent);
Parameter currentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TransmembraneCurrent);
// Parameter dependentVoltage = currentClampDevice.getCurrentClampStimulus().getVoltageParameter();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, null), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), null)));
// varHash.addVariable(newFunctionOrConstant(getMathSymbol(dependentVoltage,null),getIdentifierSubstitutions(currentClampDevice.getDependentVoltageExpression(),dependentVoltage.getUnitDefinition(),null)));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = currentClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getExpression() != null) {
// guards against voltage parameters that are "variable".
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), null)));
}
}
} else if (devices[j] instanceof VoltageClampElectricalDevice) {
VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
// total current = current source (no capacitance)
Parameter totalCurrent = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter totalCurrentParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter voltageParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_Voltage);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrent, null), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, null), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), null)));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = voltageClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getRole() == ElectricalDevice.ROLE_UserDefined) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), null)));
}
}
}
}
} else {
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
}
}
}
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
Membrane.MembraneVoltage membraneVoltage = membraneMapping.getMembrane().getMembraneVoltage();
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membraneMapping.getMembrane());
// ElectricalDevice membraneDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i].hasCapacitance() && membraneDevices[i].getDependentVoltageExpression() == null) {
if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
if (getResolved(membraneMapping)) {
//
if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
// varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
}
} else {
//
// spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
//
varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
}
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
varHash.addVariable(initVoltageFunction);
} else {
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
varHash.addVariable(voltageFunction);
}
}
}
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), sm)));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
if (initParm != null) {
Expression initExpr = new Expression(initParm.getExpression());
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextMapping scm = getSpeciesContextMapping(speciesContextSpecs[i].getSpeciesContext());
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null && (scm.isPDERequired())) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, sm), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, sm), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, sm), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), sm)));
}
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getN_PMOLE().getName(), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getKMILLIVOLTS().getName(), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getK_GHK().getName(), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
// geometric functions
//
ModelUnitSystem modelUnitSystem = simContext.getModel().getUnitSystem();
VCUnitDefinition lengthInverseUnit = modelUnitSystem.getLengthUnit().getInverse();
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumeFraction);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_SurfaceToVolumeRatio);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
if (sm instanceof MembraneMapping && !getResolved(sm)) {
MembraneMapping mm = (MembraneMapping) sm;
parm = ((MembraneMapping) sm).getVolumeFractionParameter();
if (parm.getExpression() == null) {
throw new MappingException("volume fraction not specified for feature '" + structTopology.getInsideFeature(mm.getMembrane()).getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), sm)));
parm = mm.getSurfaceToVolumeParameter();
if (parm.getExpression() == null) {
throw new MappingException("surface to volume ratio not specified for membrane '" + mm.getMembrane().getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), lengthInverseUnit, sm)));
}
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null) {
if (simContext.getGeometry().getDimension() == 0) {
if (sizeParm.getExpression() != null) {
try {
double value = sizeParm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(sizeParm, sm), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
e.printStackTrace(System.out);
throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
}
}
} else {
String compartmentName = null;
VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
String sizeFunctionName = null;
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
if (getResolved(mm)) {
FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
compartmentName = getSubVolume(fm_inside).getName() + "_" + getSubVolume(fm_outside).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
} else {
FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
if (getSubVolume(fm_inside) == getSubVolume(fm_outside)) {
compartmentName = getSubVolume(fm_inside).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
} else {
throw new RuntimeException("unexpected structure mapping for membrane '" + mm.getMembrane().getName() + "'");
}
}
} else if (sm instanceof FeatureMapping) {
FeatureMapping fm = (FeatureMapping) sm;
compartmentName = getSubVolume(fm).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
} else {
throw new RuntimeException("structure mapping " + sm.getClass().getName() + " not yet supported");
}
Expression totalVolumeCorrection = sm.getStructureSizeCorrection(simContext, this);
Expression sizeFunctionExpression = Expression.function(sizeFunctionName, new Expression[] { new Expression("'" + compartmentName + "'") });
sizeFunctionExpression.bindExpression(mathDesc);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm)));
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
}
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], null), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), null)));
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm)));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (simContext.getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// volume subdomains
//
subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
for (int j = 0; j < subVolumes.length; j++) {
SubVolume subVolume = (SubVolume) subVolumes[j];
//
// get priority of subDomain
//
int priority;
Feature spatialFeature = getResolvedFeature(subVolume);
if (spatialFeature == null) {
if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
} else {
priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
}
} else {
// now does not have to match spatial feature, *BUT* needs to be unique
priority = j;
}
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
if (spatialFeature != null) {
FeatureMapping fm = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(spatialFeature);
subDomain.setBoundaryConditionXm(fm.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(fm.getBoundaryConditionTypeXp());
if (simContext.getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(fm.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(fm.getBoundaryConditionTypeYp());
}
if (simContext.getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(fm.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(fm.getBoundaryConditionTypeZp());
}
}
//
// create equations
//
VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
//
if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
VolVariable variable = (VolVariable) scm.getVariable();
Equation equation = null;
if ((scm.isPDERequired()) && sm instanceof FeatureMapping) {
//
if (getSubVolume((FeatureMapping) sm) == subVolume) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm)));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm)));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm)));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm)));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm)));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm)));
((PdeEquation) equation).setVelocityX((scs.getVelocityXParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityXParameter(), sm)));
((PdeEquation) equation).setVelocityY((scs.getVelocityYParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityYParameter(), sm)));
((PdeEquation) equation).setVelocityZ((scs.getVelocityZParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityZParameter(), sm)));
subDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
equation = new PdeEquation(variable, initial, rate, diffusion);
if (subDomain.getEquation(variable) == null) {
subDomain.addEquation(equation);
}
}
} else {
//
// ODE
//
SubVolume mappedSubVolume = null;
if (sm instanceof FeatureMapping) {
mappedSubVolume = getSubVolume((FeatureMapping) sm);
} else if (sm instanceof MembraneMapping) {
// membrane is mapped to that of the inside feature
FeatureMapping featureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature((Membrane) sm.getStructure()));
mappedSubVolume = getSubVolume(featureMapping);
}
if (mappedSubVolume == subVolume) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
equation = new OdeEquation(variable, initial, rate);
subDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
equation = new OdeEquation(variable, initial, rate);
if (subDomain.getEquation(variable) == null) {
subDomain.addEquation(equation);
}
}
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
VCUnitDefinition subDomainUnit = modelUnitSystem.getVolumeConcentrationUnit();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
Expression rate = getIdentifierSubstitutions(scm.getFastRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(getResolvedFeature(subVolume)));
FastRate fastRate = new FastRate(rate);
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), subDomainUnit, simContext.getGeometryContext().getStructureMapping(getResolvedFeature(subVolume)));
FastInvariant fastInvariant = new FastInvariant(rate);
fastSystem.addFastInvariant(fastInvariant);
}
}
subDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create ode's for voltages to be calculated on unresolved membranes mapped to this subVolume
//
Structure[] localStructures = getStructures(subVolume);
for (int sIndex = 0; sIndex < localStructures.length; sIndex++) {
if (localStructures[sIndex] instanceof Membrane) {
Membrane membrane = (Membrane) localStructures[sIndex];
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if (!getResolved(membraneMapping) && membraneMapping.getCalculateVoltage()) {
MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
if (capacitiveDevice.getDependentVoltageExpression() == null) {
VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping));
subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping)));
} else {
//
//
//
}
}
}
}
}
//
for (int k = 0; k < subVolumes.length; k++) {
SubVolume subVolume = (SubVolume) subVolumes[k];
//
// if there is a spatially resolved membrane surrounding this subVolume, then create a membraneSubDomain
//
structures = getStructures(subVolume);
Membrane membrane = null;
if (structures != null) {
for (int j = 0; j < structures.length; j++) {
if (structures[j] instanceof Membrane && getResolved(simContext.getGeometryContext().getStructureMapping(structures[j]))) {
membrane = (Membrane) structures[j];
}
}
}
if (membrane == null) {
continue;
}
SubVolume outerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(membrane))));
SubVolume innerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(membrane))));
if (innerSubVolume != subVolume) {
throw new MappingException("membrane " + membrane.getName() + " improperly mapped to inner subVolume " + innerSubVolume.getName());
}
//
// get priority of subDomain
//
// Feature spatialFeature = simContext.getGeometryContext().getResolvedFeature(subVolume);
// int priority = spatialFeature.getPriority();
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
SurfaceClass surfaceClass = simContext.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(innerSubVolume, outerSubVolume);
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
//
// create equations for membrane-bound molecular species
//
MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(membrane);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
//
if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
//
// independant variable, create an equation object
//
Equation equation = null;
MemVariable variable = (MemVariable) scm.getVariable();
MembraneMapping mm = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(sc.getStructure());
if (scm.isPDERequired()) {
//
if (mm.getMembrane() == membrane) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), mm));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), mm)));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), mm)));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), mm)));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), mm)));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), mm)));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), mm)));
memSubDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
equation = new PdeEquation(variable, initial, rate, diffusion);
if (memSubDomain.getEquation(variable) == null) {
memSubDomain.addEquation(equation);
}
}
} else {
//
if (mm.getMembrane() == membrane) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
equation = new OdeEquation(variable, initial, rate);
memSubDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
equation = new OdeEquation(variable, initial, rate);
if (memSubDomain.getEquation(variable) == null) {
memSubDomain.addEquation(equation);
}
}
}
}
}
//
// create dummy jump conditions for all volume variables that diffuse and/or advect
//
Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
while (enum_scm.hasMoreElements()) {
SpeciesContextMapping scm = enum_scm.nextElement();
if (scm.isPDERequired()) {
// Species species = scm.getSpeciesContext().getSpecies();
Variable var = scm.getVariable();
if (var instanceof VolVariable && (scm.isPDERequired())) {
JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) var);
if (jc == null) {
// System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
jc = new JumpCondition((VolVariable) var);
memSubDomain.addJumpCondition(jc);
}
}
}
}
//
// create jump conditions for any volume variables that bind to membrane or have explicitly defined fluxes
//
ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
if (resolvedFluxes != null) {
for (int i = 0; i < resolvedFluxes.length; i++) {
Species species = resolvedFluxes[i].getSpecies();
SpeciesContext sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getInsideFeature(membraneStructureAnalyzer.getMembrane()));
if (sc == null) {
sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getOutsideFeature(membraneStructureAnalyzer.getMembrane()));
}
SpeciesContextMapping scm = getSpeciesContextMapping(sc);
// if (scm.getVariable() instanceof VolVariable && scm.isDiffusing()){
if (scm.getVariable() instanceof VolVariable && ((MembraneStructureAnalyzer.bNoFluxIfFixed || (scm.isPDERequired())))) {
if (MembraneStructureAnalyzer.bNoFluxIfFixed && !scm.isPDERequired()) {
MembraneStructureAnalyzer.bNoFluxIfFixedExercised = true;
}
JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) scm.getVariable());
if (jc == null) {
jc = new JumpCondition((VolVariable) scm.getVariable());
memSubDomain.addJumpCondition(jc);
}
Expression inFlux = getIdentifierSubstitutions(resolvedFluxes[i].inFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
jc.setInFlux(inFlux);
Expression outFlux = getIdentifierSubstitutions(resolvedFluxes[i].outFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
jc.setOutFlux(outFlux);
} else {
throw new MappingException("APPLICATION " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + membrane.getName() + ", but doesn't diffuse in compartment " + scm.getSpeciesContext().getStructure().getName());
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = membraneStructureAnalyzer.getFastSpeciesContextMappings();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
VCUnitDefinition rateUnit = scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit);
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, membraneMapping));
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, membraneMapping));
fastSystem.addFastInvariant(fastInvariant);
}
}
memSubDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create Membrane-region equations for potential of this resolved membrane
//
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if (membraneMapping.getCalculateVoltage()) {
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membrane);
int numCapacitiveDevices = 0;
MembraneElectricalDevice capacitiveDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i] instanceof MembraneElectricalDevice) {
numCapacitiveDevices++;
capacitiveDevice = (MembraneElectricalDevice) membraneDevices[i];
}
}
if (numCapacitiveDevices != 1) {
throw new MappingException("expecting 1 capacitive electrical device on graph edge for membrane " + membrane.getName() + ", found '" + numCapacitiveDevices + "'");
}
if (mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping)) instanceof MembraneRegionVariable) {
MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping());
MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping()));
memSubDomain.addEquation(vEquation);
}
}
}
// create equations for event assign targets that are model params/strutureSize, etc.
Set<VolVariable> hashKeySet = eventVolVarHash.keySet();
Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
// working under teh assumption that we are dealing with non-spatial math, hence only one compartment domain!
SubDomain subDomain = mathDesc.getSubDomains().nextElement();
while (volVarsIter.hasNext()) {
VolVariable volVar = volVarsIter.next();
EventAssignmentInitParameter eap = eventVolVarHash.get(volVar);
Expression rateExpr = new Expression(0.0);
Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(eap, null)), rateExpr);
subDomain.addEquation(equation);
}
// events - add events to math desc and odes for event assignments that have parameters as target variables
BioEvent[] bioevents = simContext.getBioEvents();
if (bioevents != null && bioevents.length > 0) {
for (BioEvent be : bioevents) {
// transform the bioEvent trigger/delay to math Event
Expression mathTriggerExpr = getIdentifierSubstitutions(be.generateTriggerExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
Delay mathDelay = null;
if (be.getParameter(BioEventParameterType.TriggerDelay) != null) {
boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
Expression mathDelayExpr = getIdentifierSubstitutions(be.getParameter(BioEventParameterType.TriggerDelay).getExpression(), timeUnit, null);
mathDelay = new Delay(bUseValsFromTriggerTime, mathDelayExpr);
}
// now deal with (bio)event Assignment translation to math EventAssignment
ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
ArrayList<Event.EventAssignment> mathEventAssignmentsList = new ArrayList<Event.EventAssignment>();
for (EventAssignment ea : eventAssignments) {
SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
VCUnitDefinition eventAssignVarUnit = ste.getUnitDefinition();
Variable variable = varHash.getVariable(ste.getName());
Event.EventAssignment mathEA = new Event.EventAssignment(variable, getIdentifierSubstitutions(ea.getAssignmentExpression(), eventAssignVarUnit, null));
mathEventAssignmentsList.add(mathEA);
}
// use the translated trigger, delay and event assignments to create (math) event
Event mathEvent = new Event(be.getName(), mathTriggerExpr, mathDelay, mathEventAssignmentsList);
mathDesc.addEvent(mathEvent);
}
}
if (!mathDesc.isValid()) {
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
// System.out.println(mathDesc.getVCML());
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
use of cbit.vcell.math.FastInvariant in project vcell by virtualcell.
the class FastSystemAnalyzer method checkLinearity.
/**
* @exception java.lang.Exception The exception description.
*/
private void checkLinearity() throws MathException, ExpressionException {
Enumeration<Variable> enum1 = fastVarList.elements();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
//
// d invariant
// for each variable, make sure ------------- = constant;
// d Var
//
Enumeration<FastInvariant> enum_fi = fastSystem.getFastInvariants();
while (enum_fi.hasMoreElements()) {
FastInvariant fi = enum_fi.nextElement();
Expression exp = fi.getFunction().differentiate(var.getName());
exp = MathUtilities.substituteFunctions(exp, this).flatten();
if (!exp.isNumeric()) {
// If expression is in terms of 'x','y','z' - then its ok - relax the constant requirement.
String[] symbols = exp.getSymbols();
for (int i = 0; i < symbols.length; i++) {
if (!symbols[i].equals(ReservedVariable.X.getName()) && !symbols[i].equals(ReservedVariable.Y.getName()) && !symbols[i].equals(ReservedVariable.Z.getName()) && !symbols[i].equals(ReservedVariable.TIME.getName())) {
throw new MathException("FastInvariant " + fi.getFunction().toString() + " isn't linear, d/d(" + var.getName() + ") = " + exp.toString());
}
}
}
}
}
}
use of cbit.vcell.math.FastInvariant in project vcell by virtualcell.
the class MathTestingUtilities method constructOdesForSensitivity.
//
// This method is used to solve for sensitivity of variables to a given parameter.
// The mathDescription and the sensitivity parameter are passed as arguments.
// New variables and ODEs are constructed according to the rule listed below and are added to the mathDescription.
// The method returns the modified mathDescription.
//
public static MathDescription constructOdesForSensitivity(MathDescription mathDesc, Constant sensParam) throws ExpressionException, MathException, MappingException {
//
// For each ODE :
//
// dX/dt = F(X, P)
//
// (where P is the sensitivity parameter)
//
// we create two other ODEs :
//
// dX_1/dt = F(X_1, P_1) and
//
// dX_2/dt = F(X_2, P_2)
//
// where P_1 = P + epsilon, and
// P_2 = P - epsilon.
//
// We keep the initial conditions for both the new ODEs to be the same,
// i.e., X_1_init = X_2_init.
//
// Then, solving for X_1 & X_2, sensitivity of X wrt P can be computed as :
//
// dX = (X_1 - X_2)
// -- ----------- .
// dP (P_1 - P_2)
//
//
// REMOVE PRINTS AFTER CHECKING !!!
System.out.println(" \n\n------------ Old Math Description -----------------");
System.out.println(mathDesc.getVCML_database());
if (mathDesc.getGeometry().getDimension() > 0) {
throw new RuntimeException("Suppport for Spatial systems not yet implemented.");
}
VariableHash varHash = new VariableHash();
Enumeration<Variable> enumVar = mathDesc.getVariables();
while (enumVar.hasMoreElements()) {
varHash.addVariable(enumVar.nextElement());
}
//
// Get 2 values of senstivity parameter (P + epsilon) & (P - epsilon)
//
Constant epsilon = new Constant("epsilon", new Expression(sensParam.getConstantValue() * 1e-3));
Constant sensParam1 = new Constant(sensParam.getName() + "_1", new Expression(sensParam.getConstantValue() + epsilon.getConstantValue()));
Constant sensParam2 = new Constant(sensParam.getName() + "_2", new Expression(sensParam.getConstantValue() - epsilon.getConstantValue()));
//
// Iterate through each subdomain (only 1 in compartmental case), and each equation in the subdomain
//
Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
//
// Create a vector of equations to store the 2 equations for each ODE variable in the subdomain.
// Later, add it to the equations list in the subdomain.
//
Vector<Equation> equnsVector = new Vector<Equation>();
Vector<Variable> varsVector = new Vector<Variable>();
Vector<Variable> var1s = new Vector<Variable>();
Vector<Variable> var2s = new Vector<Variable>();
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
Enumeration<Equation> equationEnum = subDomain.getEquations();
Domain domain = new Domain(subDomain);
while (equationEnum.hasMoreElements()) {
Equation equation = equationEnum.nextElement();
if (equation instanceof OdeEquation) {
OdeEquation odeEquation = (OdeEquation) equation;
// Similar to substituteWithExactSolutions, to bind and substitute functions in the ODE
Expression substitutedRateExp = substituteFunctions(odeEquation.getRateExpression(), mathDesc);
String varName = odeEquation.getVariable().getName();
VolVariable var = new VolVariable(varName, domain);
varsVector.addElement(var);
//
// Create the variable var1, and get the initExpr and rateExpr from the original ODE.
// Substitute the new vars (var1 and param1) in the old initExpr and rateExpr and create a new ODE
//
String varName1 = new String("__" + varName + "_1");
Expression initExpr1 = odeEquation.getInitialExpression();
Expression rateExpr1 = new Expression(substitutedRateExp);
rateExpr1.substituteInPlace(new Expression(varName), new Expression(varName1));
rateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
VolVariable var1 = new VolVariable(varName1, domain);
var1s.addElement(var1);
OdeEquation odeEqun1 = new OdeEquation(var1, initExpr1, rateExpr1);
equnsVector.addElement(odeEqun1);
//
// Create the variable var2, and get the initExpr and rateExpr from the original ODE.
// Substitute the new vars (var2 and param2) in the old initExpr and rateExpr and create a new ODE
//
String varName2 = new String("__" + varName + "_2");
Expression initExpr2 = odeEquation.getInitialExpression();
Expression rateExpr2 = new Expression(substitutedRateExp);
rateExpr2.substituteInPlace(new Expression(varName), new Expression(varName2));
rateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
VolVariable var2 = new VolVariable(varName2, domain);
var2s.addElement(var2);
OdeEquation odeEqun2 = new OdeEquation(var2, initExpr2, rateExpr2);
equnsVector.addElement(odeEqun2);
//
// Create a function for the sensitivity function expression (X1-X2)/(P1-P2), and save in varHash
//
Expression diffVar = Expression.add(new Expression(var1.getName()), Expression.negate(new Expression(var2.getName())));
Expression diffParam = Expression.add(new Expression(sensParam1.getName()), Expression.negate(new Expression(sensParam2.getName())));
Expression sensitivityExpr = Expression.mult(diffVar, Expression.invert(diffParam));
Function sens_Func = new Function("__sens" + varName + "_wrt_" + sensParam.getName(), sensitivityExpr, domain);
varHash.addVariable(epsilon);
varHash.addVariable(sensParam1);
varHash.addVariable(sensParam2);
varHash.addVariable(var1);
varHash.addVariable(var2);
varHash.addVariable(sens_Func);
} else {
// sensitivity not implemented for PDEs or other equation types.
throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
}
}
//
// Need to substitute the new variables in the new ODEs.
// i.e., if Rate Expr for ODE_1 for variable X_1 contains variable Y, variable Z, etc.
// Rate Expr is already substituted with X_1, but it also needs substitute Y with Y_1, Z with Z_1, etc.
// So get the volume variables, from the vectors for vars, var_1s and var_2s
// Substitute the rate expressions for the newly added ODEs in equnsVector.
//
Variable[] vars = (Variable[]) BeanUtils.getArray(varsVector, Variable.class);
Variable[] var_1s = (Variable[]) BeanUtils.getArray(var1s, Variable.class);
Variable[] var_2s = (Variable[]) BeanUtils.getArray(var2s, Variable.class);
Vector<Equation> newEqunsVector = new Vector<Equation>();
for (int i = 0; i < equnsVector.size(); i++) {
Equation equn = equnsVector.elementAt(i);
Expression initEx = equn.getInitialExpression();
Expression rateEx = equn.getRateExpression();
for (int j = 0; j < vars.length; j++) {
if (equn.getVariable().getName().endsWith("_1")) {
rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_1s[j].getName()));
} else if (equn.getVariable().getName().endsWith("_2")) {
rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_2s[j].getName()));
}
}
OdeEquation odeEqun = new OdeEquation(equn.getVariable(), initEx, rateEx);
newEqunsVector.addElement(odeEqun);
}
//
for (int i = 0; i < newEqunsVector.size(); i++) {
mathDesc.getSubDomain(subDomain.getName()).addEquation((Equation) newEqunsVector.elementAt(i));
}
//
// FAST SYSTEM
// If the subdomain has a fast system, create a new fast system by substituting the high-low variables/parameters
// in the expressions for the fastInvariants and fastRates and adding them to the fast system.
//
Vector<FastInvariant> invarsVector = new Vector<FastInvariant>();
Vector<FastRate> ratesVector = new Vector<FastRate>();
Enumeration<FastInvariant> fastInvarsEnum = null;
Enumeration<FastRate> fastRatesEnum = null;
// Get the fast invariants and fast rates in the system.
FastSystem fastSystem = subDomain.getFastSystem();
if (fastSystem != null) {
fastInvarsEnum = fastSystem.getFastInvariants();
fastRatesEnum = fastSystem.getFastRates();
//
while (fastInvarsEnum.hasMoreElements()) {
FastInvariant fastInvar = fastInvarsEnum.nextElement();
Expression fastInvarExpr = fastInvar.getFunction();
fastInvarExpr = MathUtilities.substituteFunctions(fastInvarExpr, mathDesc);
Expression fastInvarExpr1 = new Expression(fastInvarExpr);
Expression fastInvarExpr2 = new Expression(fastInvarExpr);
for (int i = 0; i < vars.length; i++) {
fastInvarExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
fastInvarExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
}
fastInvarExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
FastInvariant fastInvar1 = new FastInvariant(fastInvarExpr1);
invarsVector.addElement(fastInvar1);
fastInvarExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
FastInvariant fastInvar2 = new FastInvariant(fastInvarExpr2);
invarsVector.addElement(fastInvar2);
}
// Add the newly created fast invariants to the existing list of fast invariants in the fast system.
for (int i = 0; i < invarsVector.size(); i++) {
FastInvariant inVar = (FastInvariant) invarsVector.elementAt(i);
fastSystem.addFastInvariant(inVar);
}
//
while (fastRatesEnum.hasMoreElements()) {
FastRate fastRate = fastRatesEnum.nextElement();
Expression fastRateExpr = fastRate.getFunction();
fastRateExpr = MathUtilities.substituteFunctions(fastRateExpr, mathDesc);
Expression fastRateExpr1 = new Expression(fastRateExpr);
Expression fastRateExpr2 = new Expression(fastRateExpr);
for (int i = 0; i < vars.length; i++) {
fastRateExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
fastRateExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
}
fastRateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
FastRate fastRate1 = new FastRate(fastRateExpr1);
ratesVector.addElement(fastRate1);
fastRateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
FastRate fastRate2 = new FastRate(fastRateExpr2);
ratesVector.addElement(fastRate2);
}
// Add the newly created fast rates to the existing list of fast rates in the fast system.
for (int i = 0; i < ratesVector.size(); i++) {
FastRate rate = (FastRate) ratesVector.elementAt(i);
fastSystem.addFastRate(rate);
}
}
}
// Reset all variables in mathDesc.
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
// REMOVE PRINTS AFTER CHECKING
System.out.println(" \n\n------------ New Math Description -----------------");
System.out.println(mathDesc.getVCML_database());
return mathDesc;
}
use of cbit.vcell.math.FastInvariant 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.math.FastInvariant in project vcell by virtualcell.
the class XmlReader method getFastSystem.
/**
* This method returns a FastSystemImplicit from a XML Element.
* Creation date: (5/18/2001 2:38:56 PM)
* @return cbit.vcell.math.FastSystemImplicit
* @param param org.jdom.Element
* @exception cbit.vcell.xml.XmlParseException The exception description.
*/
private FastSystem getFastSystem(Element param, MathDescription mathDesc) throws XmlParseException {
// Create a new FastSystem
FastSystem fastSystem = new FastSystem(mathDesc);
// Process the FastInvariants
Iterator<Element> iterator = param.getChildren(XMLTags.FastInvariantTag, vcNamespace).iterator();
FastInvariant fastInvariant = null;
while (iterator.hasNext()) {
Element tempElement = (Element) iterator.next();
String temp = tempElement.getText();
try {
Expression newExp = unMangleExpression(temp);
fastInvariant = new FastInvariant(newExp);
fastSystem.addFastInvariant(fastInvariant);
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding the FastInvariant " + fastInvariant + ", to a FastSystem!" + " : ", e);
}
}
// Process the FastRate
iterator = param.getChildren(XMLTags.FastRateTag, vcNamespace).iterator();
FastRate fastRate = null;
while (iterator.hasNext()) {
Element tempElement = (Element) iterator.next();
String temp = tempElement.getText();
try {
Expression newExp = unMangleExpression(temp);
fastRate = new FastRate(newExp);
fastSystem.addFastRate(fastRate);
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding the FastRate " + fastRate + ", to a FastSystem!", e);
}
}
return fastSystem;
}
Aggregations