use of cbit.vcell.model.Product in project vcell by virtualcell.
the class NetworkTransformer method transform.
private void transform(SimulationContext simContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback, NetworkGenerationRequirements networkGenerationRequirements) {
String msg = "Generating network: flattening...";
mathMappingCallback.setMessage(msg);
TaskCallbackMessage tcm = new TaskCallbackMessage(TaskCallbackStatus.Clean, "");
simContext.appendToConsole(tcm);
tcm = new TaskCallbackMessage(TaskCallbackStatus.TaskStart, msg);
simContext.appendToConsole(tcm);
long startTime = System.currentTimeMillis();
System.out.println("Convert to bngl, execute BNG, retrieve the results.");
try {
BNGOutputSpec outputSpec = generateNetwork(simContext, mathMappingCallback, networkGenerationRequirements);
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
long endTime = System.currentTimeMillis();
long elapsedTime = endTime - startTime;
System.out.println(" " + elapsedTime + " milliseconds");
Model model = transformedSimulationContext.getModel();
ReactionContext reactionContext = transformedSimulationContext.getReactionContext();
// ---- Parameters -----------------------------------------------------------------------------------------------
startTime = System.currentTimeMillis();
for (int i = 0; i < outputSpec.getBNGParams().length; i++) {
BNGParameter p = outputSpec.getBNGParams()[i];
// System.out.println(i+1 + ":\t\t"+ p.toString());
if (model.getRbmModelContainer().getParameter(p.getName()) != null) {
// if it's already there we don't try to add it again; this should be true for all of them!
continue;
}
String s = p.getName();
FakeSeedSpeciesInitialConditionsParameter fakeICParam = FakeSeedSpeciesInitialConditionsParameter.fromString(s);
if (speciesEquivalenceMap.containsKey(fakeICParam)) {
// we get rid of the fake parameters we use as keys
continue;
}
FakeReactionRuleRateParameter fakeKineticParam = FakeReactionRuleRateParameter.fromString(s);
if (fakeKineticParam != null) {
System.out.println("found fakeKineticParam " + fakeKineticParam.fakeParameterName);
// we get rid of the fake parameters we use as keys
continue;
}
throw new RuntimeException("unexpected parameter " + p.getName() + " in internal BNG processing");
// Expression exp = new Expression(p.getValue());
// exp.bindExpression(model.getRbmModelContainer().getSymbolTable());
// model.getRbmModelContainer().addParameter(p.getName(), exp, model.getUnitSystem().getInstance_TBD());
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGParams().length + " parameters to model, " + elapsedTime + " ms";
System.out.println(msg);
// ---- Species ------------------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding species...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f);
startTime = System.currentTimeMillis();
System.out.println("\nSpecies :");
// the reactions will need this map to recover the names of species knowing only the networkFileIndex
HashMap<Integer, String> speciesMap = new HashMap<Integer, String>();
LinkedHashMap<String, Species> sMap = new LinkedHashMap<String, Species>();
LinkedHashMap<String, SpeciesContext> scMap = new LinkedHashMap<String, SpeciesContext>();
LinkedHashMap<String, BNGSpecies> crossMap = new LinkedHashMap<String, BNGSpecies>();
List<SpeciesContext> noMapForThese = new ArrayList<SpeciesContext>();
// final int decimalTickCount = Math.max(outputSpec.getBNGSpecies().length/10, 1);
for (int i = 0; i < outputSpec.getBNGSpecies().length; i++) {
BNGSpecies s = outputSpec.getBNGSpecies()[i];
// System.out.println(i+1 + ":\t\t"+ s.toString());
String key = s.getConcentration().infix();
FakeSeedSpeciesInitialConditionsParameter fakeParam = FakeSeedSpeciesInitialConditionsParameter.fromString(key);
if (fakeParam != null) {
Pair<SpeciesContext, Expression> value = speciesEquivalenceMap.get(fakeParam);
// the species context of the original model
SpeciesContext originalsc = value.one;
Expression initial = value.two;
// replace the fake initial condition with the real one
s.setConcentration(initial);
// we'll have to find the species context from the cloned model which correspond to the original species
SpeciesContext sc = model.getSpeciesContext(originalsc.getName());
// System.out.println(sc.getName() + ", " + sc.getSpecies().getCommonName() + " ...is one of the original seed species.");
// existing name
speciesMap.put(s.getNetworkFileIndex(), sc.getName());
sMap.put(sc.getName(), sc.getSpecies());
scMap.put(sc.getName(), sc);
crossMap.put(sc.getName(), s);
noMapForThese.add(sc);
continue;
}
// all these species are new!
// generate unique name for the species
int count = 0;
String speciesName = null;
String nameRoot = "s";
String speciesPatternNameString = s.extractName();
while (true) {
speciesName = nameRoot + count;
if (Model.isNameUnused(speciesName, model) && !sMap.containsKey(speciesName) && !scMap.containsKey(speciesName)) {
break;
}
count++;
}
// newly created name
speciesMap.put(s.getNetworkFileIndex(), speciesName);
SpeciesContext speciesContext;
if (s.hasCompartment()) {
String speciesPatternCompartmentString = s.extractCompartment();
speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(speciesPatternCompartmentString), null);
} else {
speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(0), null);
}
speciesContext.setName(speciesName);
try {
if (speciesPatternNameString != null) {
SpeciesPattern sp = RbmUtils.parseSpeciesPattern(speciesPatternNameString, model);
speciesContext.setSpeciesPattern(sp);
}
} catch (ParseException e) {
e.printStackTrace();
throw new RuntimeException("Bad format for species pattern string: " + e.getMessage());
}
// speciesContext.setSpeciesPatternString(speciesPatternString);
// model.addSpecies(speciesContext.getSpecies());
// model.addSpeciesContext(speciesContext);
sMap.put(speciesName, speciesContext.getSpecies());
scMap.put(speciesName, speciesContext);
crossMap.put(speciesName, s);
// }
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
// if(i%50 == 0) {
// System.out.println(i+"");
// }
// if(i%decimalTickCount == 0) {
// int multiplier = i/decimalTickCount;
// float progress = progressFractionQuota/4.0f + progressFractionQuotaSpecies*multiplier;
// mathMappingCallback.setProgressFraction(progress);
// }
}
for (SpeciesContext sc1 : model.getSpeciesContexts()) {
boolean found = false;
for (Map.Entry<String, SpeciesContext> entry : scMap.entrySet()) {
SpeciesContext sc2 = entry.getValue();
if (sc1.getName().equals(sc2.getName())) {
found = true;
// System.out.println("found species context " + sc1.getName() + " of species " + sc1.getSpecies().getCommonName() + " // " + sc2.getSpecies().getCommonName());
break;
}
}
if (found == false) {
// we add to the map the species context and the species which exist in the model but which are not in the map yet
// the only ones in this situation should be plain species which were not given to bngl for flattening (they are flat already)
// System.out.println("species context " + sc1.getName() + " not found in the map. Adding it.");
scMap.put(sc1.getName(), sc1);
sMap.put(sc1.getName(), sc1.getSpecies());
noMapForThese.add(sc1);
}
}
for (Species s1 : model.getSpecies()) {
boolean found = false;
for (Map.Entry<String, Species> entry : sMap.entrySet()) {
Species s2 = entry.getValue();
if (s1.getCommonName().equals(s2.getCommonName())) {
found = true;
// System.out.println("found species " + s1.getCommonName());
break;
}
}
if (found == false) {
System.err.println("species " + s1.getCommonName() + " not found in the map!");
}
}
SpeciesContext[] sca = new SpeciesContext[scMap.size()];
scMap.values().toArray(sca);
Species[] sa = new HashSet<Species>(sMap.values()).toArray(new Species[0]);
model.setSpecies(sa);
model.setSpeciesContexts(sca);
boolean isSpatial = transformedSimulationContext.getGeometry().getDimension() > 0;
for (SpeciesContext sc : sca) {
if (noMapForThese.contains(sc)) {
continue;
}
SpeciesContextSpec scs = reactionContext.getSpeciesContextSpec(sc);
Parameter param = scs.getParameter(SpeciesContextSpec.ROLE_InitialConcentration);
BNGSpecies s = crossMap.get(sc.getName());
param.setExpression(s.getConcentration());
SpeciesContext origSpeciesContext = simContext.getModel().getSpeciesContext(s.getName());
if (origSpeciesContext != null) {
ModelEntityMapping em = new ModelEntityMapping(origSpeciesContext, sc);
entityMappings.add(em);
} else {
ModelEntityMapping em = new ModelEntityMapping(new GeneratedSpeciesSymbolTableEntry(sc), sc);
if (isSpatial) {
scs.initializeForSpatial();
}
entityMappings.add(em);
}
}
// for(SpeciesContext sc : sca) { // clean all the species patterns from the flattened species, we have no sp now
// sc.setSpeciesPattern(null);
// }
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGSpecies().length + " species to model, " + elapsedTime + " ms";
System.out.println(msg);
// ---- Reactions -----------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding reactions...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f * 3.0f);
startTime = System.currentTimeMillis();
System.out.println("\nReactions :");
Map<String, HashSet<String>> ruleKeyMap = new HashMap<String, HashSet<String>>();
Map<String, BNGReaction> directBNGReactionsMap = new HashMap<String, BNGReaction>();
Map<String, BNGReaction> reverseBNGReactionsMap = new HashMap<String, BNGReaction>();
for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
BNGReaction r = outputSpec.getBNGReactions()[i];
if (!r.isRuleReversed()) {
// direct
directBNGReactionsMap.put(r.getKey(), r);
} else {
reverseBNGReactionsMap.put(r.getKey(), r);
}
//
// for each rule name, store set of keySets (number of unique keysets are number of generated reactions from this ruleName).
//
HashSet<String> keySet = ruleKeyMap.get(r.getRuleName());
if (keySet == null) {
keySet = new HashSet<String>();
ruleKeyMap.put(r.getRuleName(), keySet);
}
keySet.add(r.getKey());
}
Map<String, ReactionStep> reactionStepMap = new HashMap<String, ReactionStep>();
for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
BNGReaction bngReaction = outputSpec.getBNGReactions()[i];
// System.out.println(i+1 + ":\t\t"+ r.writeReaction());
String baseName = bngReaction.getRuleName();
String reactionName = null;
HashSet<String> keySetsForThisRule = ruleKeyMap.get(bngReaction.getRuleName());
if (keySetsForThisRule.size() == 1 && model.getReactionStep(bngReaction.getRuleName()) == null && !reactionStepMap.containsKey(bngReaction.getRuleName())) {
// we can reuse the reaction rule labels
reactionName = bngReaction.getRuleName();
} else {
reactionName = bngReaction.getRuleName() + "_0";
while (true) {
if (model.getReactionStep(reactionName) == null && !reactionStepMap.containsKey(reactionName)) {
// we can reuse the reaction rule labels
break;
}
reactionName = TokenMangler.getNextEnumeratedToken(reactionName);
}
}
//
if (directBNGReactionsMap.containsValue(bngReaction)) {
BNGReaction forwardBNGReaction = bngReaction;
BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
String name = forwardBNGReaction.getRuleName();
if (name.endsWith(ReactionRule.DirectHalf)) {
name = name.substring(0, name.indexOf(ReactionRule.DirectHalf));
}
if (name.endsWith(ReactionRule.InverseHalf)) {
name = name.substring(0, name.indexOf(ReactionRule.InverseHalf));
}
ReactionRule rr = model.getRbmModelContainer().getReactionRule(name);
Structure structure = rr.getStructure();
boolean bReversible = reverseBNGReaction != null;
SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
for (int j = 0; j < forwardBNGReaction.getReactants().length; j++) {
BNGSpecies s = forwardBNGReaction.getReactants()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Reactant reactant = sr.getReactant(scName);
if (reactant == null) {
int stoichiometry = 1;
sr.addReactant(sc, stoichiometry);
} else {
int stoichiometry = reactant.getStoichiometry();
stoichiometry += 1;
reactant.setStoichiometry(stoichiometry);
}
}
for (int j = 0; j < forwardBNGReaction.getProducts().length; j++) {
BNGSpecies s = forwardBNGReaction.getProducts()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Product product = sr.getProduct(scName);
if (product == null) {
int stoichiometry = 1;
sr.addProduct(sc, stoichiometry);
} else {
int stoichiometry = product.getStoichiometry();
stoichiometry += 1;
product.setStoichiometry(stoichiometry);
}
}
MassActionKinetics targetKinetics = new MassActionKinetics(sr);
sr.setKinetics(targetKinetics);
KineticsParameter kforward = targetKinetics.getForwardRateParameter();
KineticsParameter kreverse = targetKinetics.getReverseRateParameter();
String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
if (!kforward.getName().equals(kforwardNewName)) {
targetKinetics.renameParameter(kforward.getName(), kforwardNewName);
kforward = targetKinetics.getForwardRateParameter();
}
final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
if (!kreverse.getName().equals(kreverseNewName)) {
targetKinetics.renameParameter(kreverse.getName(), kreverseNewName);
kreverse = targetKinetics.getReverseRateParameter();
}
applyKineticsExpressions(forwardBNGReaction, kforward, targetKinetics);
if (reverseBNGReaction != null) {
applyKineticsExpressions(reverseBNGReaction, kreverse, targetKinetics);
}
// String fieldParameterName = kforward.getName();
// fieldParameterName += "_" + r.getRuleName();
// kforward.setName(fieldParameterName);
reactionStepMap.put(reactionName, sr);
} else if (reverseBNGReactionsMap.containsValue(bngReaction) && !directBNGReactionsMap.containsKey(bngReaction.getKey())) {
// reverse only (must be irreversible)
BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
ReactionRule rr = model.getRbmModelContainer().getReactionRule(reverseBNGReaction.extractRuleName());
Structure structure = rr.getStructure();
boolean bReversible = false;
SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
for (int j = 0; j < reverseBNGReaction.getReactants().length; j++) {
BNGSpecies s = reverseBNGReaction.getReactants()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Reactant reactant = sr.getReactant(scName);
if (reactant == null) {
int stoichiometry = 1;
sr.addReactant(sc, stoichiometry);
} else {
int stoichiometry = reactant.getStoichiometry();
stoichiometry += 1;
reactant.setStoichiometry(stoichiometry);
}
}
for (int j = 0; j < reverseBNGReaction.getProducts().length; j++) {
BNGSpecies s = reverseBNGReaction.getProducts()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Product product = sr.getProduct(scName);
if (product == null) {
int stoichiometry = 1;
sr.addProduct(sc, stoichiometry);
} else {
int stoichiometry = product.getStoichiometry();
stoichiometry += 1;
product.setStoichiometry(stoichiometry);
}
}
MassActionKinetics k = new MassActionKinetics(sr);
sr.setKinetics(k);
KineticsParameter kforward = k.getForwardRateParameter();
KineticsParameter kreverse = k.getReverseRateParameter();
String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
if (!kforward.getName().equals(kforwardNewName)) {
k.renameParameter(kforward.getName(), kforwardNewName);
kforward = k.getForwardRateParameter();
}
final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
if (!kreverse.getName().equals(kreverseNewName)) {
k.renameParameter(kreverse.getName(), kreverseNewName);
kreverse = k.getReverseRateParameter();
}
applyKineticsExpressions(reverseBNGReaction, kforward, k);
// String fieldParameterName = kforward.getName();
// fieldParameterName += "_" + r.getRuleName();
// kforward.setName(fieldParameterName);
reactionStepMap.put(reactionName, sr);
}
}
for (ReactionStep rs : model.getReactionSteps()) {
reactionStepMap.put(rs.getName(), rs);
}
ReactionStep[] reactionSteps = new ReactionStep[reactionStepMap.size()];
reactionStepMap.values().toArray(reactionSteps);
model.setReactionSteps(reactionSteps);
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGReactions().length + " reactions to model, " + elapsedTime + " ms";
System.out.println(msg);
// clean all the reaction rules
model.getRbmModelContainer().getReactionRuleList().clear();
// ---- Observables -------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding observables...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 8.0f * 7.0f);
startTime = System.currentTimeMillis();
System.out.println("\nObservables :");
RbmModelContainer rbmmc = model.getRbmModelContainer();
for (int i = 0; i < outputSpec.getObservableGroups().length; i++) {
ObservableGroup o = outputSpec.getObservableGroups()[i];
if (rbmmc.getParameter(o.getObservableGroupName()) != null) {
System.out.println(" ...already exists.");
// if it's already there we don't try to add it again; this should be true for all of them!
continue;
}
ArrayList<Expression> terms = new ArrayList<Expression>();
for (int j = 0; j < o.getListofSpecies().length; j++) {
Expression term = Expression.mult(new Expression(o.getSpeciesMultiplicity()[j]), new Expression(speciesMap.get(o.getListofSpecies()[j].getNetworkFileIndex())));
terms.add(term);
}
Expression exp = Expression.add(terms.toArray(new Expression[terms.size()])).flatten();
exp.bindExpression(rbmmc.getSymbolTable());
RbmObservable originalObservable = rbmmc.getObservable(o.getObservableGroupName());
VCUnitDefinition observableUnitDefinition = originalObservable.getUnitDefinition();
rbmmc.removeObservable(originalObservable);
Parameter newParameter = rbmmc.addParameter(o.getObservableGroupName(), exp, observableUnitDefinition);
RbmObservable origObservable = simContext.getModel().getRbmModelContainer().getObservable(o.getObservableGroupName());
ModelEntityMapping em = new ModelEntityMapping(origObservable, newParameter);
entityMappings.add(em);
}
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getObservableGroups().length + " observables to model, " + elapsedTime + " ms";
System.out.println(msg);
} catch (PropertyVetoException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ExpressionBindingException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ModelException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ExpressionException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ClassNotFoundException ex) {
throw new RuntimeException(ex.getMessage());
} catch (IOException ex) {
throw new RuntimeException(ex.getMessage());
}
System.out.println("Done transforming");
msg = "Generating math...";
System.out.println(msg);
mathMappingCallback.setMessage(msg);
mathMappingCallback.setProgressFraction(progressFractionQuota);
}
use of cbit.vcell.model.Product 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 ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
use of cbit.vcell.model.Product in project vcell by virtualcell.
the class RulebasedTransformer method transform.
private void transform(SimulationContext originalSimContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback) throws PropertyVetoException {
Model newModel = transformedSimulationContext.getModel();
Model originalModel = originalSimContext.getModel();
ModelEntityMapping em = null;
// list of rules created from the reactions; we apply the symmetry factor computed by bionetgen only to these
Set<ReactionRule> fromReactions = new HashSet<>();
for (SpeciesContext newSpeciesContext : newModel.getSpeciesContexts()) {
final SpeciesContext originalSpeciesContext = originalModel.getSpeciesContext(newSpeciesContext.getName());
// map new and old species contexts
em = new ModelEntityMapping(originalSpeciesContext, newSpeciesContext);
entityMappings.add(em);
if (newSpeciesContext.hasSpeciesPattern()) {
// it's perfect already and can't be improved
continue;
}
try {
MolecularType newmt = newModel.getRbmModelContainer().createMolecularType();
newModel.getRbmModelContainer().addMolecularType(newmt, false);
MolecularTypePattern newmtp_sc = new MolecularTypePattern(newmt);
SpeciesPattern newsp_sc = new SpeciesPattern();
newsp_sc.addMolecularTypePattern(newmtp_sc);
newSpeciesContext.setSpeciesPattern(newsp_sc);
RbmObservable newo = new RbmObservable(newModel, "O0_" + newmt.getName() + "_tot", newSpeciesContext.getStructure(), RbmObservable.ObservableType.Molecules);
MolecularTypePattern newmtp_ob = new MolecularTypePattern(newmt);
SpeciesPattern newsp_ob = new SpeciesPattern();
newsp_ob.addMolecularTypePattern(newmtp_ob);
newo.addSpeciesPattern(newsp_ob);
newModel.getRbmModelContainer().addObservable(newo);
// map new observable to old species context
em = new ModelEntityMapping(originalSpeciesContext, newo);
entityMappings.add(em);
} catch (ModelException e) {
e.printStackTrace();
throw new RuntimeException("unable to transform species context: " + e.getMessage());
} catch (PropertyVetoException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
ReactionSpec[] reactionSpecs = transformedSimulationContext.getReactionContext().getReactionSpecs();
for (ReactionSpec reactionSpec : reactionSpecs) {
if (reactionSpec.isExcluded()) {
// we create rules only from those reactions which are not excluded
continue;
}
ReactionStep rs = reactionSpec.getReactionStep();
String name = rs.getName();
String mangled = TokenMangler.fixTokenStrict(name);
mangled = newModel.getReactionName(mangled);
Kinetics k = rs.getKinetics();
if (!(k instanceof MassActionKinetics)) {
throw new RuntimeException("Only Mass Action Kinetics supported at this time, reaction \"" + rs.getName() + "\" uses kinetic law type \"" + rs.getKinetics().getName() + "\"");
}
boolean bReversible = rs.isReversible();
ReactionRule rr = new ReactionRule(newModel, mangled, rs.getStructure(), bReversible);
fromReactions.add(rr);
MassActionKinetics massActionKinetics = (MassActionKinetics) k;
List<Reactant> rList = rs.getReactants();
List<Product> pList = rs.getProducts();
// counting the stoichiometry - 2A+B means 3 reactants
int numReactants = 0;
for (Reactant r : rList) {
numReactants += r.getStoichiometry();
if (numReactants > 2) {
String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
throw new RuntimeException(message);
}
}
int numProducts = 0;
for (Product p : pList) {
numProducts += p.getStoichiometry();
if (bReversible && numProducts > 2) {
String message = "NFSim doesn't support more than 2 products within a reversible reaction: " + name;
throw new RuntimeException(message);
}
}
RateLawType rateLawType = RateLawType.MassAction;
RbmKineticLaw kineticLaw = new RbmKineticLaw(rr, rateLawType);
try {
String forwardRateName = massActionKinetics.getForwardRateParameter().getName();
Expression forwardRateExp = massActionKinetics.getForwardRateParameter().getExpression();
String reverseRateName = massActionKinetics.getReverseRateParameter().getName();
Expression reverseRateExp = massActionKinetics.getReverseRateParameter().getExpression();
LocalParameter fR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
fR.setName(forwardRateName);
LocalParameter rR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
rR.setName(reverseRateName);
if (rs.hasReactant()) {
kineticLaw.setParameterValue(fR, forwardRateExp, true);
}
if (rs.hasProduct()) {
kineticLaw.setParameterValue(rR, reverseRateExp, true);
}
//
for (KineticsParameter reaction_p : massActionKinetics.getKineticsParameters()) {
if (reaction_p.getRole() == Kinetics.ROLE_UserDefined) {
LocalParameter rule_p = kineticLaw.getLocalParameter(reaction_p.getName());
if (rule_p == null) {
//
// after lazy parameter creation we didn't find a user-defined rule parameter with this same name.
//
// there must be a global symbol with the same name, that the local reaction parameter has overridden.
//
ParameterContext.LocalProxyParameter rule_proxy_parameter = null;
for (ProxyParameter proxyParameter : kineticLaw.getProxyParameters()) {
if (proxyParameter.getName().equals(reaction_p.getName())) {
rule_proxy_parameter = (LocalProxyParameter) proxyParameter;
}
}
if (rule_proxy_parameter != null) {
// we want to convert to local
boolean bConvertToGlobal = false;
kineticLaw.convertParameterType(rule_proxy_parameter, bConvertToGlobal);
} else {
// could find neither local parameter nor proxy parameter
throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " didn't map to a reactionRule parameter");
}
} else if (rule_p.getRole() == RbmKineticLawParameterType.UserDefined) {
kineticLaw.setParameterValue(rule_p, reaction_p.getExpression(), true);
rule_p.setUnitDefinition(reaction_p.getUnitDefinition());
} else {
throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " mapped to a reactionRule parameter with unexpected role " + rule_p.getRole().getDescription());
}
}
}
} catch (ExpressionException e) {
e.printStackTrace();
throw new RuntimeException("Problem attempting to set RbmKineticLaw expression: " + e.getMessage());
}
rr.setKineticLaw(kineticLaw);
KineticsParameter[] kpList = k.getKineticsParameters();
ModelParameter[] mpList = rs.getModel().getModelParameters();
ModelParameter mp = rs.getModel().getModelParameter(kpList[0].getName());
ReactionParticipant[] rpList = rs.getReactionParticipants();
for (ReactionParticipant p : rpList) {
if (p instanceof Reactant) {
int stoichiometry = p.getStoichiometry();
for (int i = 0; i < stoichiometry; i++) {
SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
ReactantPattern reactantPattern = new ReactantPattern(speciesPattern, p.getStructure());
rr.addReactant(reactantPattern);
}
} else if (p instanceof Product) {
int stoichiometry = p.getStoichiometry();
for (int i = 0; i < stoichiometry; i++) {
SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
ProductPattern productPattern = new ProductPattern(speciesPattern, p.getStructure());
rr.addProduct(productPattern);
}
}
}
// commented code below is probably obsolete, we verify (above) in the reaction the number of participants,
// no need to do it again in the corresponding rule
// if(rr.getReactantPatterns().size() > 2) {
// String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
// throw new RuntimeException(message);
// }
// if(rr.getProductPatterns().size() > 2) {
// String message = "NFSim doesn't support more than 2 products within a reaction: " + name;
// throw new RuntimeException(message);
// }
newModel.removeReactionStep(rs);
newModel.getRbmModelContainer().addReactionRule(rr);
}
for (ReactionRuleSpec rrs : transformedSimulationContext.getReactionContext().getReactionRuleSpecs()) {
if (rrs == null) {
continue;
}
ReactionRule rr = rrs.getReactionRule();
if (rrs.isExcluded()) {
// delete those rules which are disabled (excluded) in the Specifications / Reaction table
newModel.getRbmModelContainer().removeReactionRule(rr);
continue;
}
}
// now that we generated the rules we can delete the reaction steps they're coming from
for (ReactionStep rs : newModel.getReactionSteps()) {
newModel.removeReactionStep(rs);
}
try {
// we invoke bngl just for the purpose of generating the xml file, which we'll then use to extract the symmetry factor
generateNetwork(transformedSimulationContext, fromReactions, mathMappingCallback);
} catch (ClassNotFoundException | IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
System.out.println("Finished RuleBased Transformer.");
}
use of cbit.vcell.model.Product in project vcell by virtualcell.
the class ModelProcessEquation method parseReaction.
public static ReactionParticipant[] parseReaction(ReactionStep reactionStep, Model model, String equationString) throws ExpressionException, PropertyVetoException {
int gotoIndex = equationString.indexOf(REACTION_GOESTO);
if (gotoIndex < 1 && equationString.length() == 0) {
throw new ExpressionException("Syntax error! " + REACTION_GOESTO + " not found. (e.g. a+b->c)");
}
if (reactionStep == null) {
return null;
}
String leftHand = equationString.substring(0, gotoIndex);
String rightHand = equationString.substring(gotoIndex + REACTION_GOESTO.length());
StringTokenizer st = new StringTokenizer(leftHand, "+");
ArrayList<ReactionParticipant> rplist = new ArrayList<ReactionParticipant>();
HashMap<String, SpeciesContext> speciesContextMap = new HashMap<String, SpeciesContext>();
Structure rxnStructure = reactionStep.getStructure();
while (st.hasMoreElements()) {
String nextToken = st.nextToken().trim();
if (nextToken.length() == 0) {
continue;
}
int stoichiIndex = 0;
while (true) {
if (Character.isDigit(nextToken.charAt(stoichiIndex))) {
stoichiIndex++;
} else {
break;
}
}
int stoichi = 1;
String tmp = nextToken.substring(0, stoichiIndex);
if (tmp.length() > 0) {
stoichi = Integer.parseInt(tmp);
}
String var = nextToken.substring(stoichiIndex).trim();
SpeciesContext sc = model.getSpeciesContext(var);
if (sc == null) {
sc = speciesContextMap.get(var);
if (sc == null) {
Species species = model.getSpecies(var);
if (species == null) {
species = new Species(var, null);
}
sc = new SpeciesContext(species, rxnStructure);
sc.setName(var);
speciesContextMap.put(var, sc);
}
}
// if (reactionStep instanceof SimpleReaction) {
rplist.add(new Reactant(null, (SimpleReaction) reactionStep, sc, stoichi));
// } else if (reactionStep instanceof FluxReaction) {
// rplist.add(new Flux(null, (FluxReaction) reactionStep, sc));
// }
}
st = new StringTokenizer(rightHand, "+");
while (st.hasMoreElements()) {
String nextToken = st.nextToken().trim();
if (nextToken.length() == 0) {
continue;
}
int stoichiIndex = 0;
while (true) {
if (Character.isDigit(nextToken.charAt(stoichiIndex))) {
stoichiIndex++;
} else {
break;
}
}
int stoichi = 1;
String tmp = nextToken.substring(0, stoichiIndex);
if (tmp.length() > 0) {
stoichi = Integer.parseInt(tmp);
}
String var = nextToken.substring(stoichiIndex);
SpeciesContext sc = model.getSpeciesContext(var);
if (sc == null) {
sc = speciesContextMap.get(var);
if (sc == null) {
Species species = model.getSpecies(var);
if (species == null) {
species = new Species(var, null);
}
sc = new SpeciesContext(species, rxnStructure);
sc.setName(var);
speciesContextMap.put(var, sc);
}
}
// if (reactionStep instanceof SimpleReaction) {
rplist.add(new Product(null, (SimpleReaction) reactionStep, sc, stoichi));
// } else if (reactionStep instanceof FluxReaction) {
// rplist.add(new Flux(null, (FluxReaction) reactionStep, sc));
// }
}
return rplist.toArray(new ReactionParticipant[0]);
}
use of cbit.vcell.model.Product in project vcell by virtualcell.
the class BioModelEditor method setRightBottomPanelOnSelection.
@Override
protected void setRightBottomPanelOnSelection(Object[] selections) {
if (selections == null) {
return;
}
JComponent bottomComponent = rightBottomEmptyPanel;
int destComponentIndex = DocumentEditorTabID.object_properties.ordinal();
boolean bShowInDatabaseProperties = false;
boolean bShowPathway = false;
if (selections.length == 1) {
Object singleSelection = selections[0];
if (singleSelection instanceof ReactionStep) {
bottomComponent = getReactionPropertiesPanel();
} else if (singleSelection instanceof ReactionRule) {
bottomComponent = getReactionRulePropertiesPanel();
} else if (singleSelection instanceof SpeciesContext) {
bottomComponent = getSpeciesPropertiesPanel();
} else if (singleSelection instanceof MolecularType) {
bottomComponent = getMolecularTypePropertiesPanel();
} else if (singleSelection instanceof RbmObservable) {
bottomComponent = getObservablePropertiesPanel();
} else if (singleSelection instanceof Structure) {
bottomComponent = getStructurePropertiesPanel();
getStructurePropertiesPanel().setModel(bioModel.getModel());
} else if (singleSelection instanceof Parameter) {
bottomComponent = getParameterPropertiesPanel();
} else if (singleSelection instanceof SimulationContext) {
bottomComponent = getApplicationPropertiesPanel();
} else if (singleSelection instanceof ParameterEstimationTask) {
bottomComponent = parameterEstimationTaskPropertiesPanel;
} else if (singleSelection instanceof Product || singleSelection instanceof Reactant) {
bottomComponent = getReactionParticipantPropertiesPanel();
} else if (singleSelection instanceof BioModelInfo) {
bShowInDatabaseProperties = true;
bottomComponent = bioModelMetaDataPanel;
} else if (singleSelection instanceof MathModelInfo) {
bShowInDatabaseProperties = true;
bottomComponent = mathModelMetaDataPanel;
} else if (singleSelection instanceof GeometryInfo) {
bShowInDatabaseProperties = true;
bottomComponent = geometryMetaDataPanel;
} else if (singleSelection instanceof SpeciesContextSpec) {
bottomComponent = getSpeciesContextSpecPanel();
} else if (singleSelection instanceof ReactionSpec) {
bottomComponent = getKineticsTypeTemplatePanel();
} else if (singleSelection instanceof ReactionRuleSpec) {
//
bottomComponent = getReactionRuleSpecPropertiesPanel();
} else if (singleSelection instanceof BioModelsNetModelInfo) {
bShowInDatabaseProperties = true;
bottomComponent = getBioModelsNetPropertiesPanel();
} else if (singleSelection instanceof Simulation) {
bottomComponent = getSimulationSummaryPanel();
} else if (singleSelection instanceof DataSymbol) {
bottomComponent = getDataSymbolsSpecPanel();
} else if (singleSelection instanceof BioEvent) {
bottomComponent = getEventPanel();
} else if (singleSelection instanceof SpatialObject) {
bottomComponent = getSpatialObjectPropertyPanel();
} else if (singleSelection instanceof SpatialProcess) {
bottomComponent = getSpatialProcessPropertyPanel();
} else if (singleSelection instanceof BioPaxObject) {
bottomComponent = bioPaxObjectPropertiesPanel;
} else if (singleSelection instanceof BioModel || singleSelection instanceof VCMetaData) {
bottomComponent = bioModelEditorAnnotationPanel;
} else if (singleSelection instanceof PathwayData) {
bShowPathway = true;
bottomComponent = getBioModelEditorPathwayPanel();
} else if (singleSelection instanceof Model) {
} else if (singleSelection instanceof RuleParticipantSignature) {
bottomComponent = getReactionRuleParticipantSignaturePropertiesPanel();
} else if (singleSelection instanceof CSGObject) {
bottomComponent = csgObjectPropertiesPanel;
csgObjectPropertiesPanel.setSimulationContext(getSelectedSimulationContext());
} else if (singleSelection instanceof DocumentEditorTreeFolderNode) {
DocumentEditorTreeFolderClass folderClass = ((DocumentEditorTreeFolderNode) singleSelection).getFolderClass();
if ((folderClass == DocumentEditorTreeFolderClass.REACTIONS_NODE) && !(singleSelection instanceof ReactionRule)) {
bottomComponent = getReactionPropertiesPanel();
} else if ((folderClass == DocumentEditorTreeFolderClass.REACTIONS_NODE) && (singleSelection instanceof ReactionRule)) {
bottomComponent = getReactionRulePropertiesPanel();
} else if (folderClass == DocumentEditorTreeFolderClass.STRUCTURES_NODE) {
bottomComponent = getStructurePropertiesPanel();
} else if (folderClass == DocumentEditorTreeFolderClass.SPECIES_NODE) {
bottomComponent = getSpeciesPropertiesPanel();
} else if (folderClass == DocumentEditorTreeFolderClass.MOLECULAR_TYPES_NODE) {
bottomComponent = getMolecularTypePropertiesPanel();
} else if (folderClass == DocumentEditorTreeFolderClass.OBSERVABLES_NODE) {
bottomComponent = getObservablePropertiesPanel();
} else if (folderClass == DocumentEditorTreeFolderClass.APPLICATIONS_NODE) {
bottomComponent = getApplicationsPropertiesPanel();
getApplicationsPropertiesPanel().setBioModel(bioModel);
} else if (folderClass == DocumentEditorTreeFolderClass.PARAMETER_ESTIMATION_NODE) {
bottomComponent = parameterEstimationTaskPropertiesPanel;
}
}
}
if (bShowPathway) {
for (destComponentIndex = 0; destComponentIndex < rightBottomTabbedPane.getTabCount(); destComponentIndex++) {
if (rightBottomTabbedPane.getComponentAt(destComponentIndex) == bottomComponent) {
break;
}
}
String tabTitle = "Pathway Preview";
if (rightBottomTabbedPane.getTabCount() == destComponentIndex) {
rightBottomTabbedPane.addTab(tabTitle, new TabCloseIcon(), bottomComponent);
}
} else if (bShowInDatabaseProperties) {
for (destComponentIndex = 0; destComponentIndex < rightBottomTabbedPane.getTabCount(); destComponentIndex++) {
Component c = rightBottomTabbedPane.getComponentAt(destComponentIndex);
if (c == bioModelMetaDataPanel || c == mathModelMetaDataPanel || c == geometryMetaDataPanel || c == getBioModelsNetPropertiesPanel()) {
break;
}
}
if (rightBottomTabbedPane.getTabCount() == destComponentIndex) {
rightBottomTabbedPane.addTab(DATABASE_PROPERTIES_TAB_TITLE, new TabCloseIcon(), bottomComponent);
}
}
if (rightBottomTabbedPane.getComponentAt(destComponentIndex) != bottomComponent) {
bottomComponent.setBorder(GuiConstants.TAB_PANEL_BORDER);
rightBottomTabbedPane.setComponentAt(destComponentIndex, bottomComponent);
rightSplitPane.repaint();
}
if (rightBottomTabbedPane.getSelectedComponent() != bottomComponent) {
rightBottomTabbedPane.setSelectedComponent(bottomComponent);
}
}
Aggregations