use of cbit.vcell.math.MathException in project vcell by virtualcell.
the class NFsimXMLWriter method getListOfSpecies.
private static Element getListOfSpecies(MathDescription mathDesc, SimulationSymbolTable simulationSymbolTable) throws SolverException {
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) mathDesc.getSubDomains().nextElement();
//
// NFsim expects a list of seed species. These are concrete species patterns that have a "ParticleProperties" element defined (for initial conditions).
//
Element listOfSpeciesElement = new Element("ListOfSpecies");
for (int speciesIndex = 0; speciesIndex < compartmentSubDomain.getParticleProperties().size(); speciesIndex++) {
// seedSpecies.getSpeciesPattern().resolveBonds();
ParticleProperties particleProperties = compartmentSubDomain.getParticleProperties().get(speciesIndex);
Element speciesElement = new Element("Species");
String speciesID = "S" + (speciesIndex + 1);
ParticleSpeciesPattern seedSpecies = (ParticleSpeciesPattern) particleProperties.getVariable();
speciesElement.setAttribute("id", speciesID);
speciesElement.setAttribute("name", seedSpecies.getName());
List<ParticleInitialCondition> particleInitialConditions = particleProperties.getParticleInitialConditions();
if (particleInitialConditions.size() != 1) {
throw new SolverException("multiple particle initial conditions not expected for " + ParticleSpeciesPattern.class.getSimpleName() + " " + seedSpecies.getName());
}
// the initial conditions must be a count in math (ParticleInitialConditionCount)
if (!(particleInitialConditions.get(0) instanceof ParticleInitialConditionCount)) {
throw new SolverException("expecting initial count for " + ParticleSpeciesPattern.class.getSimpleName() + " " + seedSpecies.getName());
}
ParticleInitialConditionCount initialCount = (ParticleInitialConditionCount) particleInitialConditions.get(0);
try {
double value = evaluateConstant(initialCount.getCount(), simulationSymbolTable);
Integer maxMoleculesPerType = simulationSymbolTable.getSimulation().getSolverTaskDescription().getNFSimSimulationOptions().getMaxMoleculesPerType();
if (maxMoleculesPerType == null) {
maxMoleculesPerType = NFsimSimulationOptions.DefaultMaxMoleculesPerType;
}
if (maxMoleculesPerType.doubleValue() < value) {
String eMessage = "The Initial count for Species '" + seedSpecies.getName() + "' is " + BigDecimal.valueOf(value).toBigInteger();
eMessage += ", which is higher than the limit of " + maxMoleculesPerType + ".\n";
eMessage += "Please do one of the following: \n- reduce the Initial Condition value for this Species or reduce the compartment size\n";
eMessage += "- increase the maximal number of Molecules per Molecular Type in the Advanced Solver Options panel.";
throw new RuntimeException(eMessage);
}
speciesElement.setAttribute("concentration", Double.toString(value));
} catch (ExpressionException | MathException e) {
e.printStackTrace();
throw new SolverException("error processing initial count of " + ParticleSpeciesPattern.class.getSimpleName() + " " + seedSpecies.getName() + ": " + e.getMessage());
}
HashMap<Bond, BondSites> bondSiteMapping = new HashMap<Bond, BondSites>();
Element listOfMoleculesElement = getListOfMolecules(speciesID, seedSpecies, bondSiteMapping);
speciesElement.addContent(listOfMoleculesElement);
if (bondSiteMapping.size() > 0) {
Element listOfBondsElement = getListOfBonds(bondSiteMapping);
speciesElement.addContent(listOfBondsElement);
}
listOfSpeciesElement.addContent(speciesElement);
}
return listOfSpeciesElement;
}
use of cbit.vcell.math.MathException in project vcell by virtualcell.
the class SmoldynFileWriter method writeMolecules.
private void writeMolecules() throws ExpressionException, MathException {
// write molecules
StringBuilder sb = new StringBuilder();
int max_mol = 0;
Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
while (subDomainEnumeration.hasMoreElements()) {
SubDomain subDomain = subDomainEnumeration.nextElement();
for (ParticleProperties particleProperties : subDomain.getParticleProperties()) {
ArrayList<ParticleInitialCondition> particleInitialConditions = particleProperties.getParticleInitialConditions();
String variableName = getVariableName(particleProperties.getVariable(), subDomain);
for (ParticleInitialCondition pic : particleInitialConditions) {
if (pic instanceof ParticleInitialConditionCount) {
max_mol += writeInitialCount((ParticleInitialConditionCount) pic, subDomain, variableName, sb);
} else if (pic instanceof ParticleInitialConditionConcentration) {
max_mol += writeInitialConcentration((ParticleInitialConditionConcentration) pic, subDomain, particleProperties.getVariable(), variableName, sb);
}
}
if (lg.isDebugEnabled()) {
lg.debug("subdomain " + subDomain.getName() + ' ' + variableName + " processed, maximum mol estimate now " + max_mol);
}
}
}
if (max_mol > MAX_MOLECULE_LIMIT) {
throw new MathException(VCellErrorMessages.getSmoldynMaxMolReachedErrorMessage((long) max_mol, MAX_MOLECULE_LIMIT));
}
int max_adjusted = max_mol * MOLECULE_MAX_COEFFICIENT;
if (max_adjusted < MIN_MOLECULE_LIMIT) {
if (lg.isInfoEnabled()) {
lg.info("adjusting computed max " + max_adjusted + " to minimum " + MIN_MOLECULE_LIMIT);
}
max_adjusted = MIN_MOLECULE_LIMIT;
}
if (max_adjusted > MAX_MOLECULE_LIMIT) {
if (lg.isInfoEnabled()) {
lg.info("adjusting computed max " + max_adjusted + " to maximum " + MAX_MOLECULE_LIMIT);
}
max_adjusted = MAX_MOLECULE_LIMIT;
}
printWriter.println("# molecules");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_mol + " " + max_adjusted);
printWriter.println(sb);
}
use of cbit.vcell.math.MathException 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();
}
use of cbit.vcell.math.MathException in project vcell by virtualcell.
the class SmoldynSolver method initialize.
/**
* This method takes the place of the old runUnsteady()...
*/
protected void initialize() throws SolverException {
if (lg.isTraceEnabled())
lg.trace("SmoldynSolver.initialize()");
fireSolverStarting(SimulationMessage.MESSAGE_SOLVEREVENT_STARTING_INIT);
writeFunctionsFile();
// write subdomains file
try {
SubdomainInfo.write(new File(getBaseName() + SimDataConstants.SUBDOMAINS_FILE_SUFFIX), simTask.getSimulation().getMathDescription());
} catch (IOException e1) {
e1.printStackTrace();
throw new SolverException(e1.getMessage());
} catch (MathException e1) {
e1.printStackTrace();
throw new SolverException(e1.getMessage());
}
String inputFilename = getInputFilename();
if (lg.isTraceEnabled())
lg.trace("SmoldynSolver.initialize() baseName = " + getBaseName());
setSolverStatus(new SolverStatus(SolverStatus.SOLVER_RUNNING, SimulationMessage.MESSAGE_SOLVER_RUNNING_INPUT_FILE));
fireSolverStarting(SimulationMessage.MESSAGE_SOLVEREVENT_STARTING_INPUT_FILE);
PrintWriter pw = null;
try {
pw = new PrintWriter(inputFilename);
if (SystemUtils.IS_OS_WINDOWS) {
//
// the windows executable is compiled under cygwin (or mingw) and expect Unix style line termination characters.
// so SmoldynInput file is written with platform default Java convention (different on windows) to a String.
// Then the string is translated to Unix style before written to the file.
// smoldyn is particularly sensitive to this issue, other compiled solvers are more tolerant.
//
StringWriter stringWriter = new StringWriter();
PrintWriter pw2 = new PrintWriter(stringWriter);
SmoldynFileWriter stFileWriter = new SmoldynFileWriter(pw2, false, getBaseName(), simTask, bMessaging);
stFileWriter.write();
String fileContents = stringWriter.getBuffer().toString();
fileContents = fileContents.replace("\r\n", "\n");
pw.write(fileContents);
} else {
//
// for linux or macos, no translation is necessary.
//
SmoldynFileWriter stFileWriter = new SmoldynFileWriter(pw, false, getBaseName(), simTask, bMessaging);
stFileWriter.write();
}
} catch (Exception e) {
setSolverStatus(new SolverStatus(SolverStatus.SOLVER_ABORTED, SimulationMessage.solverAborted("Could not generate input file: " + e.getMessage())));
e.printStackTrace(System.out);
throw new SolverException(e.getMessage());
} finally {
if (pw != null) {
pw.close();
}
}
setSolverStatus(new SolverStatus(SolverStatus.SOLVER_RUNNING, SimulationMessage.MESSAGE_SOLVER_RUNNING_START));
// get executable path+name.
setMathExecutable(new MathExecutable(getMathExecutableCommand(), getSaveDirectory()));
// setMathExecutable(new cbit.vcell.solvers.MathExecutable(executableName + " gibson " + getBaseName() + ".stochInput" + " " + getBaseName() + ".stoch"));
}
use of cbit.vcell.math.MathException in project vcell by virtualcell.
the class ModelOptimizationMapping method getRemappedReferenceData.
/**
* Gets the constraintData property (cbit.vcell.opt.ConstraintData) value.
* @return The constraintData property value.
* @see #setConstraintData
*/
private ReferenceData getRemappedReferenceData(MathMapping mathMapping) throws MappingException {
if (modelOptimizationSpec.getReferenceData() == null) {
return null;
}
//
// make sure time is mapped
//
ReferenceData refData = modelOptimizationSpec.getReferenceData();
ReferenceDataMappingSpec[] refDataMappingSpecs = modelOptimizationSpec.getReferenceDataMappingSpecs();
RowColumnResultSet rowColResultSet = new RowColumnResultSet();
Vector<SymbolTableEntry> modelObjectList = new Vector<SymbolTableEntry>();
Vector<double[]> dataList = new Vector<double[]>();
//
// find bound columns, (time is always mapped to the first column)
//
int mappedColumnCount = 0;
for (int i = 0; i < refDataMappingSpecs.length; i++) {
SymbolTableEntry modelObject = refDataMappingSpecs[i].getModelObject();
if (modelObject != null) {
int mappedColumnIndex = mappedColumnCount;
if (modelObject instanceof Model.ReservedSymbol && ((ReservedSymbol) modelObject).isTime()) {
mappedColumnIndex = 0;
}
String origRefDataColumnName = refDataMappingSpecs[i].getReferenceDataColumnName();
int origRefDataColumnIndex = refData.findColumn(origRefDataColumnName);
if (origRefDataColumnIndex < 0) {
throw new RuntimeException("reference data column named '" + origRefDataColumnName + "' not found");
}
double[] columnData = refData.getDataByColumn(origRefDataColumnIndex);
if (modelObjectList.contains(modelObject)) {
throw new RuntimeException("multiple reference data columns mapped to same model object '" + modelObject.getName() + "'");
}
modelObjectList.insertElementAt(modelObject, mappedColumnIndex);
dataList.insertElementAt(columnData, mappedColumnIndex);
mappedColumnCount++;
}
}
//
if (modelObjectList.size() == 0) {
throw new RuntimeException("reference data was not associated with model");
}
if (modelObjectList.size() == 1) {
throw new RuntimeException("reference data was not associated with model, must map time and at least one other column");
}
boolean bFoundTimeVar = false;
for (SymbolTableEntry ste : modelObjectList) {
if (ste instanceof Model.ReservedSymbol && ((ReservedSymbol) ste).isTime()) {
bFoundTimeVar = true;
break;
}
}
if (!bFoundTimeVar) {
throw new RuntimeException("must map time column of reference data to model");
}
//
for (int i = 0; i < modelObjectList.size(); i++) {
SymbolTableEntry modelObject = (SymbolTableEntry) modelObjectList.elementAt(i);
try {
// Find by name because MathSybolMapping has different 'objects' than refDataMapping 'objects'
Variable variable = mathMapping.getMathSymbolMapping().findVariableByName(modelObject.getName());
if (variable != null) {
String symbol = variable.getName();
rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
} else if (modelObject instanceof Model.ReservedSymbol && ((Model.ReservedSymbol) modelObject).isTime()) {
Model.ReservedSymbol time = (Model.ReservedSymbol) modelObject;
String symbol = time.getName();
rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
}
} catch (MathException | MatrixException | ExpressionException | ModelException e) {
e.printStackTrace();
throw new MappingException(e.getMessage(), e);
}
}
//
// populate data columns (time and rest)
//
double[] weights = new double[rowColResultSet.getColumnDescriptionsCount()];
weights[0] = 1.0;
int numRows = ((double[]) dataList.elementAt(0)).length;
int numColumns = modelObjectList.size();
for (int j = 0; j < numRows; j++) {
double[] row = new double[numColumns];
for (int i = 0; i < numColumns; i++) {
row[i] = ((double[]) dataList.elementAt(i))[j];
if (i > 0) {
weights[i] += row[i] * row[i];
}
}
rowColResultSet.addRow(row);
}
for (int i = 0; i < numColumns; i++) {
if (weights[i] == 0) {
weights[i] = 1;
} else {
weights[i] = 1 / weights[i];
}
}
SimpleReferenceData remappedRefData = new SimpleReferenceData(rowColResultSet, weights);
return remappedRefData;
}
Aggregations