Example 1 with MacroscopicRateConstant
use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.
the class NFsimXMLWriter method writeNFsimXML.
public static Element writeNFsimXML(SimulationTask origSimTask, long randomSeed, NFsimSimulationOptions nfsimSimulationOptions, boolean bUseLocationMarks) throws SolverException {
try {
System.out.println("VCML ORIGINAL .... START\n" + origSimTask.getSimulation().getMathDescription().getVCML_database() + "\nVCML ORIGINAL .... END\n====================\n");
} catch (MathException e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
SimulationTask clonedSimTask = null;
try {
clonedSimTask = (SimulationTask) BeanUtils.cloneSerializable(origSimTask);
} catch (Exception eee) {
throw new SolverException("failed to clone mathDescription while preparing NFSim input: " + eee.getMessage(), eee);
}
MathDescription clonedMathDesc = clonedSimTask.getSimulation().getMathDescription();
if (bUseLocationMarks) {
try {
//
// get list of Compartment Names (stored in locations).
//
ArrayList<String> locations = new ArrayList<String>();
Enumeration<Variable> varEnum = clonedMathDesc.getVariables();
ArrayList<VolumeParticleSpeciesPattern> volumeParticleSpeciesPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
while (varEnum.hasMoreElements()) {
Variable var = varEnum.nextElement();
if (var instanceof VolumeParticleSpeciesPattern) {
VolumeParticleSpeciesPattern speciesPattern = (VolumeParticleSpeciesPattern) var;
if (!locations.contains(speciesPattern.getLocationName())) {
locations.add(speciesPattern.getLocationName());
}
volumeParticleSpeciesPatterns.add(speciesPattern);
}
}
//
for (ParticleMolecularType particleMolecularType : clonedMathDesc.getParticleMolecularTypes()) {
String pmcLocationName = RbmUtils.SiteStruct;
String pmcLocationId = particleMolecularType.getName() + "_" + RbmUtils.SiteStruct;
ParticleMolecularComponent locationComponent = new ParticleMolecularComponent(pmcLocationId, pmcLocationName);
for (String location : locations) {
locationComponent.addComponentStateDefinition(new ParticleComponentStateDefinition(location));
}
particleMolecularType.insertMolecularComponent(0, locationComponent);
String pmcMarkName = RbmUtils.SiteProduct;
String pmcMarkId = particleMolecularType.getName() + "_" + RbmUtils.SiteProduct;
ParticleMolecularComponent markComponent = new ParticleMolecularComponent(pmcMarkId, pmcMarkName);
markComponent.addComponentStateDefinition(new ParticleComponentStateDefinition("0"));
markComponent.addComponentStateDefinition(new ParticleComponentStateDefinition("1"));
particleMolecularType.insertMolecularComponent(1, markComponent);
}
//
for (VolumeParticleSpeciesPattern speciesPattern : volumeParticleSpeciesPatterns) {
for (ParticleMolecularTypePattern molTypePattern : speciesPattern.getParticleMolecularTypePatterns()) {
//
// add location component to pattern ... state=<location>
//
{
final ParticleMolecularComponent locationComponentDefinition = molTypePattern.getMolecularType().getComponentList().get(0);
ParticleMolecularComponentPattern locationPattern = new ParticleMolecularComponentPattern(locationComponentDefinition);
ParticleComponentStateDefinition locationStateDefinition = null;
for (ParticleComponentStateDefinition stateDef : locationComponentDefinition.getComponentStateDefinitions()) {
if (stateDef.getName().equals(speciesPattern.getLocationName())) {
locationStateDefinition = stateDef;
}
}
ParticleComponentStatePattern locationStatePattern = new ParticleComponentStatePattern(locationStateDefinition);
locationPattern.setComponentStatePattern(locationStatePattern);
locationPattern.setBondType(ParticleBondType.None);
locationPattern.setBondId(-1);
molTypePattern.insertMolecularComponentPattern(0, locationPattern);
}
//
// add mark component to pattern ... state="0" (for observables and reactants ... later we will clone and use "1" for products).
{
final ParticleMolecularComponent markComponentDefinition = molTypePattern.getMolecularType().getComponentList().get(1);
ParticleMolecularComponentPattern markPattern = new ParticleMolecularComponentPattern(markComponentDefinition);
final int clearStateIndex = 0;
final int setStateIndex = 1;
ParticleComponentStateDefinition markStateClearedDefinition = markComponentDefinition.getComponentStateDefinitions().get(clearStateIndex);
ParticleComponentStatePattern markStatePattern = new ParticleComponentStatePattern(markStateClearedDefinition);
markPattern.setComponentStatePattern(markStatePattern);
markPattern.setBondType(ParticleBondType.None);
markPattern.setBondId(-1);
molTypePattern.insertMolecularComponentPattern(1, markPattern);
}
}
}
//
// when processing ParticleJumpProcesses, we add a new "product" species pattern (by cloning the original speciesPattern)
// and setting the mark site to "1", change name to name+"_PRODUCT", and add to math model if it doesn't already exist.
//
// cloned the "standard" reactant/observable speciesPattern, set the mark for all molecules, and add to mathDesc.
//
CompartmentSubDomain subDomain = (CompartmentSubDomain) clonedMathDesc.getSubDomains().nextElement();
for (ParticleJumpProcess particleJumpProcess : subDomain.getParticleJumpProcesses()) {
for (Action action : particleJumpProcess.getActions()) {
if (action.getOperation().equals(Action.ACTION_CREATE)) {
VolumeParticleSpeciesPattern volumeParticleSpeciesPattern = (VolumeParticleSpeciesPattern) action.getVar();
String newSpeciesPatternName = volumeParticleSpeciesPattern.getName() + "_" + particleJumpProcess.getName();
VolumeParticleSpeciesPattern productPattern = new VolumeParticleSpeciesPattern(volumeParticleSpeciesPattern, newSpeciesPatternName);
// VolumeParticleSpeciesPattern productPattern = (VolumeParticleSpeciesPattern) BeanUtils.cloneSerializable(volumeParticleSpeciesPattern);
for (ParticleMolecularTypePattern productMolTypePattern : productPattern.getParticleMolecularTypePatterns()) {
ParticleComponentStateDefinition markSet = productMolTypePattern.getMolecularType().getComponentList().get(1).getComponentStateDefinitions().get(1);
productMolTypePattern.getMolecularComponentPatternList().get(1).setComponentStatePattern(new ParticleComponentStatePattern(markSet));
}
System.out.println(productPattern.getName());
if (clonedMathDesc.getVariable(productPattern.getName()) == null) {
clonedMathDesc.addVariable(productPattern);
}
action.setVar(productPattern);
}
}
}
try {
System.out.println("===============================\n ----------- VCML HACKED .... START\n" + clonedMathDesc.getVCML_database() + "\nVCML HACKED .... END\n====================\n");
} catch (MathException e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
} catch (Exception e) {
throw new SolverException("failed to apply location mark transformation: " + e.getMessage(), e);
}
}
Element sbmlElement = new Element("sbml");
Element modelElement = new Element("model");
modelElement.setAttribute("id", "nameless");
SimulationSymbolTable simulationSymbolTable = new SimulationSymbolTable(clonedSimTask.getSimulation(), clonedSimTask.getSimulationJob().getJobIndex());
Element listOfParametersElement = getListOfParameters(clonedMathDesc, simulationSymbolTable);
Element listOfMoleculeTypesElement = getListOfMoleculeTypes(clonedMathDesc);
Element listOfSpeciesElement = getListOfSpecies(clonedMathDesc, simulationSymbolTable);
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) clonedMathDesc.getSubDomains().nextElement();
Element listOfReactionRules = new Element("ListOfReactionRules");
for (int reactionRuleIndex = 0; reactionRuleIndex < compartmentSubDomain.getParticleJumpProcesses().size(); reactionRuleIndex++) {
ParticleJumpProcess particleJumpProcess = compartmentSubDomain.getParticleJumpProcesses().get(reactionRuleIndex);
MathRuleFactory mathRuleFactory = new MathRuleFactory();
MathRuleEntry rule = mathRuleFactory.createRuleEntry(particleJumpProcess, reactionRuleIndex);
RuleAnalysisReport report = RuleAnalysis.analyze(rule, true);
// remember, we have to add RateLaw
Element reactionRuleElement = RuleAnalysis.getNFSimXML(rule, report);
// ArrayList<MolecularTypeOfReactionParticipant> currentReactantElementsOfReaction = new ArrayList<MolecularTypeOfReactionParticipant>();
// ArrayList<ComponentOfMolecularTypeOfReactionParticipant> currentComponentOfReactantElementsOfReaction = new ArrayList<ComponentOfMolecularTypeOfReactionParticipant>();
// ArrayList<MolecularTypeOfReactionParticipant> currentProductElementsOfReaction = new ArrayList<MolecularTypeOfReactionParticipant>();
// ArrayList<ComponentOfMolecularTypeOfReactionParticipant> currentComponentOfProductElementsOfReaction = new ArrayList<ComponentOfMolecularTypeOfReactionParticipant>();
// currentMappingOfReactionParticipants.clear();
// reactionProductBondSites.clear();
// reactionReactantBondSites.clear();
//
// Element reactionRuleElement = new Element("ReactionRule");
// String reactionRuleID = "RR" + (reactionRuleIndex + 1);
// reactionRuleElement.setAttribute("id",reactionRuleID);
// reactionRuleElement.setAttribute("name",particleJumpProcess.getName());
// reactionRuleElement.setAttribute("symmetry_factor","1");
// reactionRule.resolveBonds();
//
// ArrayList<VolumeParticleSpeciesPattern> selectedPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
// for (ParticleVariable particleVariable : particleJumpProcess.getParticleVariables()){
// if (!(particleVariable instanceof VolumeParticleSpeciesPattern)){
// throw new SolverException("expecting only "+VolumeParticleSpeciesPattern.class.getSimpleName()+"s for "+ParticleJumpProcess.class.getSimpleName()+" "+particleJumpProcess.getName());
// }
// selectedPatterns.add((VolumeParticleSpeciesPattern) particleVariable);
// }
// ArrayList<VolumeParticleSpeciesPattern> createdPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
// HashSet<VolumeParticleSpeciesPattern> destroyedPatterns = new HashSet<VolumeParticleSpeciesPattern>();
// for (Action action : particleJumpProcess.getActions()){
// if (!(action.getVar() instanceof VolumeParticleSpeciesPattern)){
// throw new SolverException("expecting only "+VolumeParticleSpeciesPattern.class.getSimpleName()+"s for "+ParticleJumpProcess.class.getSimpleName()+" "+particleJumpProcess.getName());
// }
// if (action.getOperation().equals(Action.ACTION_CREATE)){
// createdPatterns.add((VolumeParticleSpeciesPattern) action.getVar());
// }else if (action.getOperation().equals(Action.ACTION_DESTROY)){
// destroyedPatterns.add((VolumeParticleSpeciesPattern) action.getVar());
// }else{
// throw new RuntimeException("unexpected action operation "+action.getOperation()+" for jump process "+particleJumpProcess.getName());
// }
// }
//
// Element listOfReactantPatternsElement = new Element("ListOfReactantPatterns");
// for(int reactantPatternIndex=0; reactantPatternIndex < selectedPatterns.size(); reactantPatternIndex++) {
// VolumeParticleSpeciesPattern reactantSpeciesPattern = selectedPatterns.get(reactantPatternIndex);
// String reactantPatternID = "RP" + (reactantPatternIndex + 1);
// patternReactantBondSites.clear();
// Element reactantPatternElement = getReactionParticipantPattern1(reactionRuleID, reactantPatternID, reactantSpeciesPattern,
// currentReactantElementsOfReaction, currentComponentOfReactantElementsOfReaction, "ReactantPattern");
// listOfReactantPatternsElement.addContent(reactantPatternElement);
// reactionReactantBondSites.addAll(patternReactantBondSites);
// }
// reactionRuleElement.addContent(listOfReactantPatternsElement);
//
// Element listOfProductPatternsElement = new Element("ListOfProductPatterns");
// ArrayList<VolumeParticleSpeciesPattern> productSpeciesPatterns = new ArrayList<VolumeParticleSpeciesPattern>(selectedPatterns);
// productSpeciesPatterns.removeAll(destroyedPatterns);
// productSpeciesPatterns.addAll(createdPatterns);
// // for products, add all "created" species from Actions and all "particles" that are selected but not destroyed
// for(int productPatternIndex=0; productPatternIndex < productSpeciesPatterns.size(); productPatternIndex++) {
// VolumeParticleSpeciesPattern productSpeciesPattern = productSpeciesPatterns.get(productPatternIndex);
// String productPatternID = "PP" + (productPatternIndex + 1);
// patternProductBondSites.clear();
// Element productPatternElement = getReactionParticipantPattern1(reactionRuleID, productPatternID, productSpeciesPattern,
// currentProductElementsOfReaction, currentComponentOfProductElementsOfReaction, "ProductPattern");
// listOfProductPatternsElement.addContent(productPatternElement);
// reactionProductBondSites.addAll(patternProductBondSites);
// }
// reactionRuleElement.addContent(listOfProductPatternsElement);
// <RateLaw id="RR1_RateLaw" type="Ele" totalrate="0">
// <ListOfRateConstants>
// <RateConstant value="kon"/>
// </ListOfRateConstants>
// </RateLaw>
Element rateLawElement = new Element("RateLaw");
rateLawElement.setAttribute("id", RuleAnalysis.getID(rule));
String rateConstantValue = null;
JumpProcessRateDefinition particleProbabilityRate = particleJumpProcess.getParticleRateDefinition();
if (particleProbabilityRate.getExpressions().length > 0) {
JumpProcessRateDefinition particleRateDefinition = particleJumpProcess.getParticleRateDefinition();
Expression expression = null;
if (particleRateDefinition instanceof MacroscopicRateConstant) {
expression = ((MacroscopicRateConstant) particleProbabilityRate).getExpression();
} else {
throw new SolverException("ParticleRateDefinition type " + particleRateDefinition.getClass().getSimpleName() + " not supported");
}
rateConstantValue = expression.infixBng();
// all rates constants are being flattened and given reserved names
Expression substitutedValExpr = null;
try {
substitutedValExpr = simulationSymbolTable.substituteFunctions(expression);
} catch (MathException | ExpressionException e) {
e.printStackTrace(System.out);
throw new SolverException("ParticleJumpProcess " + particleJumpProcess.getName() + " substitution failed : exp = \"" + expression.infix() + "\": " + e.getMessage());
}
Double value = null;
try {
value = substitutedValExpr.evaluateConstant();
Element parameterElement = new Element("Parameter");
String id = "K_reserved_" + reactionRuleIndex;
parameterElement.setAttribute("id", id);
if (value != null) {
parameterElement.setAttribute("type", "Constant");
parameterElement.setAttribute("value", value.toString());
parameterElement.addContent(new Comment(rateConstantValue));
rateConstantValue = id;
listOfParametersElement.addContent(parameterElement);
}
} catch (ExpressionException e) {
System.out.println("ParticleJumpProcess " + particleJumpProcess.getName() + " = " + substitutedValExpr.infix() + " does not have a constant value");
}
}
if (isFunction(rateConstantValue, clonedMathDesc, simulationSymbolTable)) {
rateLawElement.setAttribute("type", "Function");
rateLawElement.setAttribute("totalrate", "0");
rateLawElement.setAttribute("name", rateConstantValue);
} else {
rateLawElement.setAttribute("type", "Ele");
rateLawElement.setAttribute("totalrate", "0");
Element listOfRateConstantsElement = new Element("ListOfRateConstants");
Element rateConstantElement = new Element("RateConstant");
// System.out.println(" --- " + particleJumpProcess.getParticleRateDefinition().getExpressions());
if (particleProbabilityRate.getExpressions().length > 0) {
rateConstantElement.setAttribute("value", rateConstantValue);
}
listOfRateConstantsElement.addContent(rateConstantElement);
rateLawElement.addContent(listOfRateConstantsElement);
}
reactionRuleElement.addContent(rateLawElement);
// // <Map>
// // <MapItem sourceID="RR1_RP1_M1" targetID="RR1_PP1_M1"/>
// // <MapItem sourceID="RR1_RP1_M1_C1" targetID="RR1_PP1_M1_C1"/>
// // <MapItem sourceID="RR1_RP1_M1_C2" targetID="RR1_PP1_M1_C2"/>
// // <MapItem sourceID="RR1_RP2_M1" targetID="RR1_PP1_M2"/>
// // <MapItem sourceID="RR1_RP2_M1_C1" targetID="RR1_PP1_M2_C1"/>
// // </Map>
// Element mapElement = new Element("Map");
// System.out.println("----------------------------------------------------------------------");
// for(MolecularTypeOfReactionParticipant p : currentReactantElementsOfReaction) {
// System.out.println(p.moleculeName + ", " + p.elementID);
// }
// for(ComponentOfMolecularTypeOfReactionParticipant c : currentComponentOfReactantElementsOfReaction) {
// System.out.println(c.moleculeName + ", " + c.componentName + ", " + c.elementID);
// }
// System.out.println("----------------------------------------------------------------------");
// for(MolecularTypeOfReactionParticipant p : currentProductElementsOfReaction) {
// System.out.println(p.moleculeName + ", " + p.elementID);
// }
// for(ComponentOfMolecularTypeOfReactionParticipant c : currentComponentOfProductElementsOfReaction) {
// System.out.println(c.moleculeName + ", " + c.componentName + ", " + c.elementID);
// }
// System.out.println("----------------------------------------------------------------------");
//
// List<MolecularTypeOfReactionParticipant> cloneOfReactants = new ArrayList<MolecularTypeOfReactionParticipant>(currentReactantElementsOfReaction);
// List<MolecularTypeOfReactionParticipant> cloneOfProducts = new ArrayList<MolecularTypeOfReactionParticipant>(currentProductElementsOfReaction);
// for(Iterator<MolecularTypeOfReactionParticipant> itReactant = cloneOfReactants.iterator(); itReactant.hasNext();) { // participants
// MolecularTypeOfReactionParticipant reactant = itReactant.next();
// boolean foundProduct = false;
// for(Iterator<MolecularTypeOfReactionParticipant> itProduct = cloneOfProducts.iterator(); itProduct.hasNext();) {
// MolecularTypeOfReactionParticipant product = itProduct.next();
// if(reactant.find(product)) {
// MappingOfReactionParticipants m = new MappingOfReactionParticipants(reactant.elementID, product.elementID, "");
// currentMappingOfReactionParticipants.add(m );
// itProduct.remove();
// foundProduct = true;
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// if(foundProduct == false) {
// System.out.println("Did not found a match for reactant " + reactant.moleculeName + ", " + reactant.elementID);
// }
// itReactant.remove(); // found or not, we remove the reactant
// }
// if(!currentProductElementsOfReaction.isEmpty()) {
// for(MolecularTypeOfReactionParticipant p : currentProductElementsOfReaction) {
// System.out.println("Did not found a match for product " + p.moleculeName + ", " + p.elementID);
// }
// }
// for(Iterator<ComponentOfMolecularTypeOfReactionParticipant> itReactant = currentComponentOfReactantElementsOfReaction.iterator(); itReactant.hasNext();) { // components
// ComponentOfMolecularTypeOfReactionParticipant reactant = itReactant.next();
// boolean foundProduct = false;
// for(Iterator<ComponentOfMolecularTypeOfReactionParticipant> itProduct = currentComponentOfProductElementsOfReaction.iterator(); itProduct.hasNext();) {
// ComponentOfMolecularTypeOfReactionParticipant product = itProduct.next();
// String state = "";
// if(reactant.find(product)) {
// if(!reactant.state.equals(product.state)) {
// state = product.state;
// }
// MappingOfReactionParticipants m = new MappingOfReactionParticipants(reactant.elementID, product.elementID, state);
// currentMappingOfReactionParticipants.add(m );
// itProduct.remove();
// foundProduct = true;
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// if(foundProduct == false) {
// System.out.println("Did not found a match for reactant " + reactant.moleculeName + ", " + reactant.elementID);
// }
// itReactant.remove(); // found or not, we remove the reactant
// }
// if(!currentComponentOfProductElementsOfReaction.isEmpty()) {
// for(ComponentOfMolecularTypeOfReactionParticipant p : currentComponentOfProductElementsOfReaction) {
// System.out.println("Did not found a match for product " + p.moleculeName + ", " + p.elementID);
// }
// }
// for(Iterator<MappingOfReactionParticipants> it = currentMappingOfReactionParticipants.iterator(); it.hasNext();) {
// MappingOfReactionParticipants m = it.next();
// Element mapItemElement = new Element("MapItem");
// mapItemElement.setAttribute("sourceID", m.reactantElementID);
// mapItemElement.setAttribute("targetID", m.productElementID);
// mapElement.addContent(mapItemElement);
// }
// reactionRuleElement.addContent(mapElement);
//
// // <ListOfOperations>
// // <AddBond site1="RR1_RP1_M1_C1" site2="RR1_RP2_M1_C1"/>
// // <StateChange site="RR0_RP0_M0_C2" finalState="Y"/>
// // </ListOfOperations>
// Element listOfOperationsElement = new Element("ListOfOperations");
//
// // AddBond elements
// // add any bond in the product which is not present in the reactant
// Iterator<BondSites> it = patternProductBondSites.iterator();
// while (it.hasNext()) {
// BondSites bs = it.next();
// String reactantS1 = MappingOfReactionParticipants.findMatchingReactant(bs.component1, currentMappingOfReactionParticipants);
// String reactantS2 = MappingOfReactionParticipants.findMatchingReactant(bs.component2, currentMappingOfReactionParticipants);
// // we check if the bonds in the product existed already in the reactant, in which case they were not "added" in this reaction
// BondSites candidate = new BondSites(reactantS1, reactantS2);
// boolean preExistent = false;
// for(BondSites bsReactant : reactionReactantBondSites) {
// if(bsReactant.equals(candidate)) {
// preExistent = true;
// break;
// }
// }
// if(preExistent == true) {
// continue; // we don't add preexisting bonds
// }
// Element addBondElement = new Element("AddBond");
// addBondElement.setAttribute("site1", reactantS1);
// addBondElement.setAttribute("site2", reactantS2);
// listOfOperationsElement.addContent(addBondElement);
// }
// // StateChange elements
// for(Iterator<MappingOfReactionParticipants> it1 = currentMappingOfReactionParticipants.iterator(); it1.hasNext();) {
// MappingOfReactionParticipants m = it1.next();
// if(!m.componentFinalState.equals("")) { // state has changed if it's different from ""
// Element stateChangeElement = new Element("StateChange");
// stateChangeElement.setAttribute("site", m.reactantElementID);
// stateChangeElement.setAttribute("finalState", m.componentFinalState);
// listOfOperationsElement.addContent(stateChangeElement);
// }
// }
// // eliminate all the common entries (molecule types) in reactants and products
// // what's left in reactants was deleted, what's left in products was added
// List<MolecularTypeOfReactionParticipant> commonParticipants = new ArrayList<MolecularTypeOfReactionParticipant>();
// for(Iterator<MolecularTypeOfReactionParticipant> itReactant = currentReactantElementsOfReaction.iterator(); itReactant.hasNext();) { // participants
// MolecularTypeOfReactionParticipant reactant = itReactant.next();
// for(Iterator<MolecularTypeOfReactionParticipant> itProduct = currentProductElementsOfReaction.iterator(); itProduct.hasNext();) {
// MolecularTypeOfReactionParticipant product = itProduct.next();
// if(reactant.find(product)) {
// // commonParticipants contains the reactant molecules with a equivalent molecule in the product (meaning they are not in the "Deleted" category)
// commonParticipants.add(reactant);
// itReactant.remove();
// itProduct.remove();
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// }
// // DeleteBond element
// // there is no need to mention deletion of bond if the particleSpeciesPattern
// // or the MolecularType involved in the bond are deleted as well
// // We only keep those "Deleted" bonds which belong to the molecules (of the reactant) present in commonParticipants
// // Both components (sites) of the bond need to have their molecules in commonParticipants
// boolean foundMoleculeForComponent1 = false;
// boolean foundMoleculeForComponent2 = false;
// HashSet<BondSites> cloneOfReactantBondSites = new HashSet<BondSites>(patternReactantBondSites);
// Iterator<BondSites> itbs = cloneOfReactantBondSites.iterator();
// while (itbs.hasNext()) {
// BondSites bs = itbs.next();
// String bondComponent1MoleculeId = BondSites.extractMoleculeId(bs.component1);
// String bondComponent2MoleculeId = BondSites.extractMoleculeId(bs.component2);
// for(MolecularTypeOfReactionParticipant commonReactionMoleculeule : commonParticipants) {
// String commonReactantPatternId = commonReactionMoleculeule.elementID;
// if(bondComponent1MoleculeId.equals(commonReactantPatternId)) {
// foundMoleculeForComponent1 = true;
// }
// if(bondComponent2MoleculeId.equals(commonReactantPatternId)) {
// foundMoleculeForComponent2 = true;
// }
// }
// if(!foundMoleculeForComponent1 || !foundMoleculeForComponent2) {
// // at least one of bond's molecule is not in common, hence we don't need to report the deletion of this bond
// itbs.remove();
// }
// }
// // the clone has now all the deleted bonds whose molecules have not been deleted
// itbs = cloneOfReactantBondSites.iterator();
// while (itbs.hasNext()) {
// BondSites bs = itbs.next();
// Element addBondElement = new Element("DeleteBond");
// addBondElement.setAttribute("site1", bs.component1);
// addBondElement.setAttribute("site2", bs.component2);
// listOfOperationsElement.addContent(addBondElement);
// }
// // Add MolecularType element
// for(MolecularTypeOfReactionParticipant molecule : currentProductElementsOfReaction) {
// System.out.println("created molecule: " + molecule.elementID + "' " + molecule.moleculeName);
// Element addMolecularTypePatternElement = new Element("Add");
// addMolecularTypePatternElement.setAttribute("id", molecule.elementID);
// listOfOperationsElement.addContent(addMolecularTypePatternElement);
// }
// // Delete MolecularType element
// // if the reactant pattern of the molecule being deleted still exists as part of the common, then we only delete the molecule
// // if the reactant pattern of the molecule being deleted is not as part of the common, then it's gone completely and we delete the reactant pattern
// ArrayList<String> patternsToDelete = new ArrayList<String>();
// for(MolecularTypeOfReactionParticipant molecule : currentReactantElementsOfReaction) {
// String reactantPatternId = molecule.extractReactantPatternId();
// boolean found = false;
// for(MolecularTypeOfReactionParticipant common : commonParticipants) {
// String commonId = common.extractReactantPatternId();
// if(reactantPatternId.equals(commonId)) {
// found = true;
// break; // some other molecule of this pattern still there, we don't delete the pattern
// }
// }
// if(found == true) { // some other molecule of this pattern still there, we don't delete the pattern
// System.out.println("deleted molecule: " + molecule.elementID + "' " + molecule.moleculeName);
// Element addMolecularTypePatternElement = new Element("Delete");
// addMolecularTypePatternElement.setAttribute("id", molecule.elementID);
// addMolecularTypePatternElement.setAttribute("DeleteMolecules", "0");
// listOfOperationsElement.addContent(addMolecularTypePatternElement);
// } else { // no molecule of this pattern left, we delete the pattern
// if(patternsToDelete.contains(reactantPatternId)) {
// // nothing to do, we're already deleting this pattern
// break;
// } else {
// patternsToDelete.add(reactantPatternId);
// System.out.println("deleted pattern: " + reactantPatternId);
// Element addParticleSpeciesPatternElement = new Element("Delete");
// addParticleSpeciesPatternElement.setAttribute("id", reactantPatternId);
// addParticleSpeciesPatternElement.setAttribute("DeleteMolecules", "0");
// listOfOperationsElement.addContent(addParticleSpeciesPatternElement);
// }
// }
// }
// reactionRuleElement.addContent(listOfOperationsElement);
listOfReactionRules.addContent(reactionRuleElement);
}
Element listOfObservablesElement = getListOfObservables(clonedMathDesc);
Element listOfFunctionsElement = getListOfFunctions(clonedMathDesc, simulationSymbolTable);
modelElement.addContent(listOfParametersElement);
modelElement.addContent(listOfMoleculeTypesElement);
modelElement.addContent(listOfSpeciesElement);
modelElement.addContent(listOfReactionRules);
modelElement.addContent(listOfObservablesElement);
modelElement.addContent(listOfFunctionsElement);
sbmlElement.addContent(modelElement);
// // return e1;
return sbmlElement;
}
Also used :
Action(cbit.vcell.math.Action)
SimulationTask(cbit.vcell.messaging.server.SimulationTask)
Variable(cbit.vcell.math.Variable)
MathDescription(cbit.vcell.math.MathDescription)
Element(org.jdom.Element)
ArrayList(java.util.ArrayList)
ParticleMolecularComponent(cbit.vcell.math.ParticleMolecularComponent)
ExpressionException(cbit.vcell.parser.ExpressionException)
ParticleComponentStateDefinition(cbit.vcell.math.ParticleComponentStateDefinition)
ParticleMolecularComponentPattern(cbit.vcell.math.ParticleMolecularComponentPattern)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
MathRuleFactory(cbit.vcell.math.MathRuleFactory)
ParticleMolecularType(cbit.vcell.math.ParticleMolecularType)
RuleAnalysisReport(org.vcell.model.rbm.RuleAnalysisReport)
Comment(org.jdom.Comment)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
ParticleComponentStatePattern(cbit.vcell.math.ParticleComponentStatePattern)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
VolumeParticleSpeciesPattern(cbit.vcell.math.VolumeParticleSpeciesPattern)
SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable)
SolverException(cbit.vcell.solver.SolverException)
ExpressionException(cbit.vcell.parser.ExpressionException)
MathException(cbit.vcell.math.MathException)
ParticleMolecularTypePattern(cbit.vcell.math.ParticleMolecularTypePattern)
MathRuleEntry(cbit.vcell.math.MathRuleFactory.MathRuleEntry)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SolverException(cbit.vcell.solver.SolverException)
Example 2 with MacroscopicRateConstant
use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.
the class SmoldynFileWriter method writeReactions.
private void writeReactions() throws ExpressionException, MathException {
printWriter.println("# reactions");
Enumeration<SubDomain> subdomains = mathDesc.getSubDomains();
while (subdomains.hasMoreElements()) {
SubDomain subdomain = subdomains.nextElement();
for (ParticleJumpProcess pjp : subdomain.getParticleJumpProcesses()) {
ArrayList<Variable> reactants = new ArrayList<Variable>();
ArrayList<Variable> products = new ArrayList<Variable>();
for (Action a : pjp.getActions().toArray(new Action[pjp.getActions().size()])) {
if (a.getOperation().equals(Action.ACTION_CREATE)) {
products.add(a.getVar());
} else if (a.getOperation().equals(Action.ACTION_DESTROY)) {
reactants.add(a.getVar());
}
}
Expression rateDefinition = null;
JumpProcessRateDefinition jprd = pjp.getParticleRateDefinition();
if (jprd instanceof MacroscopicRateConstant) {
rateDefinition = subsituteFlatten(((MacroscopicRateConstant) jprd).getExpression());
} else if (jprd instanceof InteractionRadius) {
rateDefinition = subsituteFlatten(((InteractionRadius) jprd).getExpression());
} else {
new RuntimeException("The jump process rate definition is not supported");
}
if (rateDefinition.isZero()) {
continue;
}
if (mathDesc.isSpatialHybrid()) {
String[] symbols = rateDefinition.getSymbols();
if (symbols != null) {
if (subdomain instanceof MembraneSubDomain) {
rateDefinition = new Expression(FiniteVolumeFileWriter.replaceVolumeVariable(getSimulationTask(), (MembraneSubDomain) subdomain, rateDefinition));
}
}
} else {
try {
rateDefinition.evaluateConstant();
} catch (ExpressionException ex) {
throw new ExpressionException("reaction rate for jump process " + pjp.getName() + " is not a constant. Constants are required for all reaction rates.");
}
}
// Smoldyn takes maximum 2nd order reaction.
if (reactants.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 reactants.");
}
if (products.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 products.");
}
String rateDefinitionStr = simulation.getMathDescription().isSpatialHybrid() ? rateDefinition.infix() + ";" : rateDefinition.evaluateConstant() + "";
if (subdomain instanceof CompartmentSubDomain) {
// 0th order reaction, product limited to one and we'll let the reaction know where it happens
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_cmpt + " " + subdomain.getName() + " " + pjp.getName() + " ");
} else {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction + " " + /* + subdomain.getName() + " "*/
pjp.getName() + " ");
}
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (subdomain instanceof MembraneSubDomain) {
// 0th order reaction, product limited to one and it can be on mem or in vol
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // consuming of a species to nothing, limited to one reactant
if (reactants.size() == 1 && products.size() == 0) {
if (// consuming a mem species in mem reaction
getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // it equals to adsorption, species A from vol adsorbed to mem as again species A, and then we kill the speceis A on mem.
if (getVolumeVariableCount(reactants) == 1) {
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
String speciesName = reactants.get(0).getName();
String killMolCmd = "cmd " + SmoldynVCellMapper.SmoldynKeyword.E + " " + SmoldynVCellMapper.SmoldynKeyword.killmol + " " + speciesName + "(" + SmoldynKeyword.up + ")";
killMolCommands.add(killMolCmd);
}
} else // Use rate command for any membrane reactions with 1 reactant and 1 product
if ((reactants.size() == 1) && (products.size() == 1)) {
// Membrane reaction (1 react to 1 product).
if (getMembraneVariableCount(products) == 1 && getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // Other single molecular reactions
{
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else // membrane reactions which are not one to one, or 0th order, or consuming species
{
if (// membrane reaction has one membrane bound reactant
(getMembraneVariableCount(reactants) == 1)) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (// bimolecular membrane reaction
getMembraneVariableCount(reactants) == 2) {
if (jprd instanceof InteractionRadius) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionByInteractionRadius(reactants, products, subdomain, rateDefinitionStr, pjp.getName());
} else {
// throw new MathException("Error with reaction: " + pjp.getName() + ".\nVCell Spatial stochastic modeling requires macroscopic or microscopic kinetics for bimolecular membrane reactions.");
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else if (getMembraneVariableCount(reactants) == 0) {
throw new MathException("Error with reaction: " + pjp.getName() + ".\nIn VCell spatial stochastic modeling, the membrane reaction requires at least one membrane bound reactant.");
}
}
}
}
}
printWriter.println();
}
Also used :
Action(cbit.vcell.math.Action)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReservedVariable(cbit.vcell.math.ReservedVariable)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ArrayList(java.util.ArrayList)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Example 3 with MacroscopicRateConstant
use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.
the class XmlReader method getParticleJumpProcess.
private ParticleJumpProcess getParticleJumpProcess(Element param, MathDescription md) throws XmlParseException {
// name
String name = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
ProcessSymmetryFactor processSymmetryFactor = null;
Attribute symmetryFactorAttr = param.getAttribute(XMLTags.ProcessSymmetryFactorAttrTag);
if (symmetryFactorAttr != null) {
processSymmetryFactor = new ProcessSymmetryFactor(Double.parseDouble(symmetryFactorAttr.getValue()));
}
// selected particle
List<ParticleVariable> varList = new ArrayList<ParticleVariable>();
Iterator<Element> iterator = param.getChildren(XMLTags.SelectedParticleTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
String varname = unMangle(tempelement.getAttributeValue(XMLTags.NameAttrTag));
Variable var = md.getVariable(varname);
if (!(var instanceof ParticleVariable)) {
throw new XmlParseException("Not a ParticleVariable in ParticleJumpProcess.");
}
varList.add((ParticleVariable) var);
}
// probability rate
JumpProcessRateDefinition jprd = null;
// for old models
Element pb = param.getChild(XMLTags.ParticleProbabilityRateTag, vcNamespace);
if (pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new MacroscopicRateConstant(exp);
} else // for new models
{
pb = param.getChild(XMLTags.MacroscopicRateConstantTag, vcNamespace);
if (// jump process rate defined by macroscopic rate constant
pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new MacroscopicRateConstant(exp);
} else // jump process rate defined by binding radius
{
pb = param.getChild(XMLTags.InteractionRadiusTag, vcNamespace);
if (pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new InteractionRadius(exp);
}
}
}
// add actions
List<Action> actionList = new ArrayList<Action>();
iterator = param.getChildren(XMLTags.ActionTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
try {
actionList.add(getAction(tempelement, md));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
} catch (ExpressionException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
ParticleJumpProcess jump = new ParticleJumpProcess(name, varList, jprd, actionList, processSymmetryFactor);
return jump;
}
Also used :
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
Action(cbit.vcell.math.Action)
FilamentVariable(cbit.vcell.math.FilamentVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
StochVolVariable(cbit.vcell.math.StochVolVariable)
RandomVariable(cbit.vcell.math.RandomVariable)
VolumeRandomVariable(cbit.vcell.math.VolumeRandomVariable)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
VolVariable(cbit.vcell.math.VolVariable)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
PointVariable(cbit.vcell.math.PointVariable)
MembraneRandomVariable(cbit.vcell.math.MembraneRandomVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MemVariable(cbit.vcell.math.MemVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
InteractionRadius(cbit.vcell.math.InteractionRadius)
Attribute(org.jdom.Attribute)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Element(org.jdom.Element)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ArrayList(java.util.ArrayList)
ExpressionException(cbit.vcell.parser.ExpressionException)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Example 4 with MacroscopicRateConstant
use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
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);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().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(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
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(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().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 = getSimulationContext().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.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
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.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
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.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
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.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
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.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
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.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// 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(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().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");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
rate constants are important in calculating the probability rate.
for Mass Action, we use KForward and KReverse,
for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
*/
Expression forwardRate = null;
Expression reverseRate = null;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
}
if (maFunc != null) {
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ReactionStep(cbit.vcell.model.ReactionStep)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Kinetics(cbit.vcell.model.Kinetics)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
GeometryClass(cbit.vcell.geometry.GeometryClass)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Pair(org.vcell.util.Pair)
ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MathException(cbit.vcell.math.MathException)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
MassActionSolver(cbit.vcell.model.MassActionSolver)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Example 5 with MacroscopicRateConstant
use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.
the class ParticleMathMapping method combineHybrid.
private void combineHybrid() throws MappingException, ExpressionException, MatrixException, MathException, ModelException {
ArrayList<SpeciesContext> continuousSpecies = new ArrayList<SpeciesContext>();
ArrayList<ParticleVariable> continuousSpeciesParticleVars = new ArrayList<ParticleVariable>();
ArrayList<SpeciesContext> stochSpecies = new ArrayList<SpeciesContext>();
//
// categorize speciesContexts as continuous and stochastic
//
SpeciesContextSpec[] scsArray = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
continuousSpecies = new ArrayList<SpeciesContext>();
stochSpecies = new ArrayList<SpeciesContext>();
for (SpeciesContextSpec speciesContextSpec : scsArray) {
if (!getSimulationContext().isStoch() || speciesContextSpec.isForceContinuous()) {
continuousSpecies.add(speciesContextSpec.getSpeciesContext());
Variable variable = getMathSymbolMapping().getVariable(speciesContextSpec.getSpeciesContext());
if (variable instanceof ParticleVariable) {
continuousSpeciesParticleVars.add((ParticleVariable) variable);
}
} else {
stochSpecies.add(speciesContextSpec.getSpeciesContext());
}
}
if (continuousSpecies.isEmpty()) {
return;
}
//
// create continuous mathDescription ... add stochastic variables and processes to the continuous Math and use this.
//
DiffEquMathMapping mathMapping = new DiffEquMathMapping(getSimulationContext(), callback, networkGenerationRequirements);
mathMapping.refresh(null);
MathDescription contMathDesc = mathMapping.getMathDescription();
//
// get list of all continuous variables
//
HashMap<String, Variable> allContinuousVars = new HashMap<String, Variable>();
Enumeration<Variable> enumVar = contMathDesc.getVariables();
while (enumVar.hasMoreElements()) {
Variable var = enumVar.nextElement();
allContinuousVars.put(var.getName(), var);
}
//
// replace those continuous variables and equations for stochastic speciesContexts
// with the particleVariables and particleProperties
// (ParticleJumpProcesses removed later)
//
ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
for (SpeciesContext stochSpeciesContext : stochSpecies) {
Variable contVar = mathMapping.getMathSymbolMapping().getVariable(stochSpeciesContext);
Variable stochVar = getMathSymbolMapping().getVariable(stochSpeciesContext);
allContinuousVars.put(stochVar.getName(), stochVar);
//
// replace continuous "concentration" VolVariable/MemVariable for this particle with a Function for concentration
//
allContinuousVars.remove(contVar);
VCUnitDefinition sizeUnit = unitSystem.getLengthUnit().raiseTo(new RationalNumber(stochSpeciesContext.getStructure().getDimension()));
VCUnitDefinition stochasticDensityUnit = unitSystem.getStochasticSubstanceUnit().divideBy(sizeUnit);
VCUnitDefinition continuousDensityUnit = unitSystem.getConcentrationUnit(stochSpeciesContext.getStructure());
if (stochasticDensityUnit.isEquivalent(continuousDensityUnit)) {
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), new Expression(stochVar, getNameScope()), contVar.getDomain()));
} else {
Expression conversionFactorExp = getUnitFactor(continuousDensityUnit.divideBy(stochasticDensityUnit));
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), Expression.mult(new Expression(stochVar, getNameScope()), conversionFactorExp), contVar.getDomain()));
}
//
// remove continuous equation
//
Enumeration<SubDomain> contSubDomains = contMathDesc.getSubDomains();
while (contSubDomains.hasMoreElements()) {
SubDomain contSubDomain = contSubDomains.nextElement();
contSubDomain.removeEquation(contVar);
if (contSubDomain instanceof MembraneSubDomain) {
((MembraneSubDomain) contSubDomain).removeJumpCondition(contVar);
}
}
//
// remove all continuous variables for speciesContextSpec parameters (e.g. initial conditions, diffusion rates, boundary conditions, velocities)
//
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(stochSpeciesContext);
Parameter[] scsParameters = scs.getParameters();
for (Parameter parameter : scsParameters) {
Variable continuousScsParmVariable = mathMapping.getMathSymbolMapping().getVariable(parameter);
allContinuousVars.remove(continuousScsParmVariable);
}
//
// copy ParticleJumpProcess and ParticleProperties to the continuous math
//
SubDomain contSubDomain = contMathDesc.getSubDomain(contVar.getDomain().getName());
SubDomain stochSubDomain = mathDesc.getSubDomain(stochVar.getDomain().getName());
ParticleProperties particleProperties = stochSubDomain.getParticleProperties(stochVar);
contSubDomain.addParticleProperties(particleProperties);
}
//
// add all ParticleJumpProcesses to the continuous model
//
Enumeration<SubDomain> enumStochSubdomains = mathDesc.getSubDomains();
while (enumStochSubdomains.hasMoreElements()) {
SubDomain stochSubdomain = enumStochSubdomains.nextElement();
SubDomain contSubdomain = contMathDesc.getSubDomain(stochSubdomain.getName());
for (ParticleJumpProcess particleJumpProcess : stochSubdomain.getParticleJumpProcesses()) {
//
// modify "selection list" (particleVariables), probability rate, and actions if referenced particleVariable is to be "forced continuous"
//
ParticleVariable[] selectedParticles = particleJumpProcess.getParticleVariables();
for (ParticleVariable particleVariable : selectedParticles) {
if (continuousSpeciesParticleVars.contains(particleVariable)) {
particleJumpProcess.remove(particleVariable);
JumpProcessRateDefinition jumpProcessRateDefinition = particleJumpProcess.getParticleRateDefinition();
if (jumpProcessRateDefinition instanceof MacroscopicRateConstant) {
MacroscopicRateConstant macroscopicRateConstant = (MacroscopicRateConstant) jumpProcessRateDefinition;
macroscopicRateConstant.setExpression(Expression.mult(macroscopicRateConstant.getExpression(), new Expression(particleVariable, null)));
} else if (jumpProcessRateDefinition instanceof InteractionRadius) {
throw new MappingException("cannot adjust interaction radius for reaction process " + particleJumpProcess.getName() + ", particle " + particleVariable.getName() + " is continuous");
} else {
throw new MappingException("rate definition type " + jumpProcessRateDefinition.getClass().getSimpleName() + " not yet implemented for hybrid PDE/Particle math generation");
}
}
Iterator<Action> iterAction = particleJumpProcess.getActions().iterator();
while (iterAction.hasNext()) {
Action action = iterAction.next();
if (continuousSpeciesParticleVars.contains(action.getVar())) {
iterAction.remove();
}
}
}
if (!particleJumpProcess.getActions().isEmpty()) {
contSubdomain.addParticleJumpProcess(particleJumpProcess);
}
}
}
//
for (MathMappingParameter mathMappingParameter : fieldMathMappingParameters) {
if (mathMappingParameter instanceof UnitFactorParameter) {
String name = mathMappingParameter.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, newFunctionOrConstant(name, mathMappingParameter.getExpression(), null));
}
}
}
//
// add constants and functions from the particle math that aren't already defined in the continuous math
//
Enumeration<Variable> enumVars = mathDesc.getVariables();
while (enumVars.hasMoreElements()) {
Variable var = enumVars.nextElement();
if (var instanceof Constant || var instanceof Function) {
String name = var.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, var);
}
}
}
contMathDesc.setAllVariables(allContinuousVars.values().toArray(new Variable[0]));
mathDesc = contMathDesc;
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
System.out.println(mathDesc.getVCML_database());
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 ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used :
GeometryClass(cbit.vcell.geometry.GeometryClass)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
MathDescription(cbit.vcell.math.MathDescription)
HashMap(java.util.HashMap)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
ArrayList(java.util.ArrayList)
SpeciesContext(cbit.vcell.model.SpeciesContext)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Function(cbit.vcell.math.Function)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Aggregations
JumpProcessRateDefinition (cbit.vcell.math.JumpProcessRateDefinition)7
MacroscopicRateConstant (cbit.vcell.math.MacroscopicRateConstant)7
Action (cbit.vcell.math.Action)6
ParticleJumpProcess (cbit.vcell.math.ParticleJumpProcess)6
Expression (cbit.vcell.parser.Expression)6
InteractionRadius (cbit.vcell.math.InteractionRadius)5
MembraneParticleVariable (cbit.vcell.math.MembraneParticleVariable)5
ParticleVariable (cbit.vcell.math.ParticleVariable)5
Variable (cbit.vcell.math.Variable)5
VolumeParticleVariable (cbit.vcell.math.VolumeParticleVariable)5
ArrayList (java.util.ArrayList)5
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)4
MathException (cbit.vcell.math.MathException)4
ExpressionException (cbit.vcell.parser.ExpressionException)4
MathDescription (cbit.vcell.math.MathDescription)3
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)3
ProcessSymmetryFactor (cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)3
SubDomain (cbit.vcell.math.SubDomain)3
Element (org.jdom.Element)3
GeometryClass (cbit.vcell.geometry.GeometryClass)2
