use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class StructureSizeSolver method updateRelativeStructureSizes.
/**
* Insert the method's description here.
* Creation date: (5/17/2006 10:33:38 AM)
* @return double[]
* @param structName java.lang.String
* @param structSize double
*/
public static void updateRelativeStructureSizes(SimulationContext simContext) throws Exception {
if (simContext.getGeometry().getDimension() > 0) {
throw new RuntimeException("not yet supported for spatial applications");
}
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
try {
// This is rewritten in Feb 2008. Siblings and children are correctly taken into account when calculating the volume fractions.
StructureTopology structTopology = simContext.getModel().getStructureTopology();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
// calculate the sum of features' sizes inside this membrane, this is used for calculating both surface volume ratio and volume fraction.
double sumOfSubFeatures = 0;
Membrane membrane = ((MembraneMapping) structureMappings[i]).getMembrane();
Enumeration<Feature> subFeatures = structTopology.getSubFeatures(structTopology.getInsideFeature(membrane));
while (subFeatures.hasMoreElements()) {
Feature feature = subFeatures.nextElement();
sumOfSubFeatures = sumOfSubFeatures + simContext.getGeometryContext().getStructureMapping(feature).getSizeParameter().getExpression().evaluateConstant();
}
// calculate the sum of features's sizes inside the membrance's parent feature, this is used for calculating the volume fraction.
double sumOfParentMemSubFeatures = 0;
Feature parentFeature = structTopology.getOutsideFeature(membrane);
if (parentFeature != null) {
Enumeration<Feature> parentSubFeatures = structTopology.getSubFeatures(parentFeature);
while (parentSubFeatures.hasMoreElements()) {
Feature feature = parentSubFeatures.nextElement();
sumOfParentMemSubFeatures = sumOfParentMemSubFeatures + simContext.getGeometryContext().getStructureMapping(feature).getSizeParameter().getExpression().evaluateConstant();
}
// set surface volume ratio
((MembraneMapping) structureMappings[i]).getSurfaceToVolumeParameter().setExpression(new Expression(((MembraneMapping) structureMappings[i]).getSizeParameter().getExpression().evaluateConstant() / sumOfSubFeatures));
// set volume fraction
((MembraneMapping) structureMappings[i]).getVolumeFractionParameter().setExpression(new Expression(sumOfSubFeatures / sumOfParentMemSubFeatures));
}
}
}
} catch (NullPointerException e) {
e.printStackTrace(System.out);
// DialogUtils.showErrorDialog("structure sizes must all be specified");
throw new Exception("structure sizes must all be specified");
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new Exception(e.getMessage());
}
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class StructureSizeSolver method updateAbsoluteStructureSizes.
public static void updateAbsoluteStructureSizes(SimulationContext simContext, Structure struct, double structSize, VCUnitDefinition structSizeUnit) throws Exception {
StructureMapping[] structMappings = simContext.getGeometryContext().getStructureMappings();
try {
StructureTopology structTopology = simContext.getModel().getStructureTopology();
SolverParameterCollection solverParameters = new SolverParameterCollection();
ArrayList<String> unknownVars = new ArrayList<String>();
for (StructureMapping sm : structMappings) {
if (sm.getStructure() instanceof Membrane) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter svRatioParam = mm.getSurfaceToVolumeParameter();
StructureMappingParameter volFractParam = mm.getVolumeFractionParameter();
StructureMappingParameter sizeParam = mm.getSizeParameter();
solverParameters.add(new SolverParameter(svRatioParam, TokenMangler.mangleToSName("sv_" + mm.getMembrane().getName()), svRatioParam.getExpression().evaluateConstant()));
solverParameters.add(new SolverParameter(volFractParam, TokenMangler.mangleToSName("vf_" + mm.getMembrane().getName()), volFractParam.getExpression().evaluateConstant()));
}
StructureMappingParameter sizeParam = sm.getSizeParameter();
Double priorKnownValue = null;
String varName = TokenMangler.mangleToSName("size_" + sm.getStructure().getName());
if (sizeParam.getExpression() != null) {
priorKnownValue = sizeParam.getExpression().evaluateConstant();
} else if (sm.getStructure() == struct) {
priorKnownValue = structSize;
} else {
unknownVars.add(varName);
}
solverParameters.add(new SolverParameter(sizeParam, varName, priorKnownValue));
}
ArrayList<Expression> expressions = new ArrayList<Expression>();
for (int i = 0; i < structMappings.length; i++) {
if (structMappings[i] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structMappings[i];
Feature insideFeature = structTopology.getInsideFeature(membraneMapping.getMembrane());
Feature outsideFeature = structTopology.getOutsideFeature(membraneMapping.getMembrane());
StructureMappingParameter sizeParameter = membraneMapping.getSizeParameter();
StructureMappingParameter volFractParameter = membraneMapping.getVolumeFractionParameter();
StructureMappingParameter surfToVolParameter = membraneMapping.getSurfaceToVolumeParameter();
//
// EC eclosing cyt, which contains er and golgi
// "(cyt_size+ er_size + golgi_size) * cyt_svRatio - PM_size" ... implicit equation, exp == 0 is implied
//
Expression sumOfInsideVolumeExp = new Expression(0.0);
for (int j = 0; j < structMappings.length; j++) {
if (structMappings[j] instanceof FeatureMapping && structTopology.enclosedBy(structMappings[j].getStructure(), insideFeature)) {
FeatureMapping childFeatureMappingOfInside = ((FeatureMapping) structMappings[j]);
sumOfInsideVolumeExp = Expression.add(sumOfInsideVolumeExp, new Expression(solverParameters.getName(childFeatureMappingOfInside.getSizeParameter())));
}
}
Expression tempExpr = Expression.mult(sumOfInsideVolumeExp, new Expression(solverParameters.getName(surfToVolParameter)));
tempExpr = Expression.add(tempExpr, new Expression("-" + solverParameters.getName(sizeParameter)));
expressions.add(tempExpr);
//
// EC eclosing cyt, which contains er and golgi
// (EC_size + cyt_size + er_size + golgi_size) * cyt_vfRatio - (cyt_size + er_size + golgi_size) ... implicit equation, exp == 0 is implied
//
Expression sumOfParentVolumeExp = new Expression(0.0);
for (int j = 0; j < structMappings.length; j++) {
if (structMappings[j] instanceof FeatureMapping && structTopology.enclosedBy(structMappings[j].getStructure(), outsideFeature)) {
FeatureMapping childFeatureMappingOfParent = ((FeatureMapping) structMappings[j]);
sumOfParentVolumeExp = Expression.add(sumOfParentVolumeExp, new Expression(TokenMangler.mangleToSName(solverParameters.getName(childFeatureMappingOfParent.getSizeParameter()))));
}
}
Expression exp = Expression.mult(sumOfParentVolumeExp, new Expression(solverParameters.getName(volFractParameter)));
exp = Expression.add(exp, Expression.negate(sumOfInsideVolumeExp));
expressions.add(exp);
}
}
if (expressions.size() != unknownVars.size()) {
throw new RuntimeException("number of unknowns is " + unknownVars.size() + ", number of equations is " + expressions.size());
}
if (unknownVars.size() == 0 && expressions.size() == 0) {
StructureMappingParameter sizeParam = simContext.getGeometryContext().getStructureMapping(struct).getSizeParameter();
sizeParam.setExpression(new Expression(structSize));
return;
}
RationalExp[][] rowColData = new RationalExp[unknownVars.size()][unknownVars.size() + 1];
for (int row = 0; row < unknownVars.size(); row++) {
//
// verify that there is no "constant" term (without an unknown)
//
// System.out.println("equation("+row+"): "+expressions.get(row).infix());
Expression constantTerm = new Expression(expressions.get(row));
for (String var : unknownVars) {
constantTerm.substituteInPlace(new Expression(var), new Expression(0.0));
}
constantTerm = constantTerm.flatten();
//
for (int col = 0; col < unknownVars.size(); col++) {
Expression equation = new Expression(expressions.get(row));
String colVariable = unknownVars.get(col);
Expression deriv = equation.differentiate(colVariable).flatten();
String[] symbols = deriv.getSymbols();
if (symbols != null) {
for (String symbol : symbols) {
if (unknownVars.contains(symbol)) {
throw new RuntimeException("equation is not linear in the unknowns");
}
}
}
rowColData[row][col] = RationalExpUtils.getRationalExp(deriv);
}
rowColData[row][unknownVars.size()] = RationalExpUtils.getRationalExp(constantTerm).minus();
}
RationalExpMatrix rationalExpMatrix = new RationalExpMatrix(rowColData);
// rationalExpMatrix.show();
RationalExp[] solutions = rationalExpMatrix.solveLinearExpressions();
double[] solutionValues = new double[solutions.length];
for (int i = 0; i < unknownVars.size(); i++) {
Expression vcSolution = new Expression(solutions[i].infixString());
String[] symbols = vcSolution.getSymbols();
if (symbols != null) {
for (String symbol : symbols) {
SolverParameter p = solverParameters.get(symbol);
if (p.knownValue == null) {
throw new RuntimeException("solution for var " + unknownVars.get(i) + " is a function of unknown var " + p.name);
}
vcSolution.substituteInPlace(new Expression(symbol), new Expression(p.knownValue.doubleValue()));
}
}
double value = vcSolution.flatten().evaluateConstant();
// System.out.println(unknownVars.get(i)+" = "+value+" = "+solutions[i].infixString());
solutionValues[i] = value;
}
for (int i = 0; i < unknownVars.size(); i++) {
SolverParameter p = solverParameters.get(unknownVars.get(i));
p.parameter.setExpression(new Expression(solutionValues[i]));
}
//
// set the one known value (if not set already by the gui).
//
StructureMappingParameter sizeParam = simContext.getGeometryContext().getStructureMapping(struct).getSizeParameter();
sizeParam.setExpression(new Expression(structSize));
System.out.println("done");
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new Exception(e.getMessage());
}
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class SBMLExporter method addCompartments.
/**
* addCompartments comment.
* @throws XMLStreamException
* @throws SbmlException
*/
protected void addCompartments() throws XMLStreamException, SbmlException {
Model vcModel = vcBioModel.getModel();
cbit.vcell.model.Structure[] vcStructures = vcModel.getStructures();
for (int i = 0; i < vcStructures.length; i++) {
Compartment sbmlCompartment = sbmlModel.createCompartment();
sbmlCompartment.setId(TokenMangler.mangleToSName(vcStructures[i].getName()));
sbmlCompartment.setName(vcStructures[i].getName());
VCUnitDefinition sbmlSizeUnit = null;
StructureTopology structTopology = getSelectedSimContext().getModel().getStructureTopology();
Structure parentStructure = structTopology.getParentStructure(vcStructures[i]);
if (vcStructures[i] instanceof Feature) {
sbmlCompartment.setSpatialDimensions(3);
String outside = null;
if (parentStructure != null) {
outside = TokenMangler.mangleToSName(parentStructure.getName());
}
if (outside != null) {
if (outside.length() > 0) {
sbmlCompartment.setOutside(outside);
}
}
sbmlSizeUnit = sbmlExportSpec.getVolumeUnits();
UnitDefinition unitDefn = getOrCreateSBMLUnit(sbmlSizeUnit);
sbmlCompartment.setUnits(unitDefn);
} else if (vcStructures[i] instanceof Membrane) {
Membrane vcMembrane = (Membrane) vcStructures[i];
sbmlCompartment.setSpatialDimensions(2);
Feature outsideFeature = structTopology.getOutsideFeature(vcMembrane);
if (outsideFeature != null) {
sbmlCompartment.setOutside(TokenMangler.mangleToSName(outsideFeature.getName()));
sbmlSizeUnit = sbmlExportSpec.getAreaUnits();
UnitDefinition unitDefn = getOrCreateSBMLUnit(sbmlSizeUnit);
sbmlCompartment.setUnits(unitDefn);
} else if (lg.isEnabledFor(Level.WARN)) {
lg.warn(this.sbmlModel.getName() + " membrame " + vcMembrane.getName() + " has not outside feature");
}
}
sbmlCompartment.setConstant(true);
StructureMapping vcStructMapping = getSelectedSimContext().getGeometryContext().getStructureMapping(vcStructures[i]);
try {
if (vcStructMapping.getSizeParameter().getExpression() != null) {
sbmlCompartment.setSize(vcStructMapping.getSizeParameter().getExpression().evaluateConstant());
} else {
// really no need to set sizes of compartments in spatial ..... ????
// throw new RuntimeException("Compartment size not set for compartment \"" + vcStructures[i].getName() + "\" ; Please set size and try exporting again.");
}
} catch (cbit.vcell.parser.ExpressionException e) {
// If it is in the catch block, it means that the compartment size was probably not a double, but an assignment.
// Check if the expression for the compartment size is not null and add it as an assignment rule.
Expression sizeExpr = vcStructMapping.getSizeParameter().getExpression();
if (sizeExpr != null) {
ASTNode ruleFormulaNode = getFormulaFromExpression(sizeExpr);
AssignmentRule assignRule = sbmlModel.createAssignmentRule();
assignRule.setVariable(vcStructures[i].getName());
assignRule.setMath(ruleFormulaNode);
// If compartmentSize is specified by an assignment rule, the 'constant' field should be set to 'false' (default - true).
sbmlCompartment.setConstant(false);
sbmlModel.addRule(assignRule);
}
}
// Add the outside compartment of given compartment as annotation to the compartment.
// This is required later while trying to read in compartments ...
Element sbmlImportRelatedElement = null;
// if (parentStructure != null) {
// sbmlImportRelatedElement = new Element(XMLTags.VCellRelatedInfoTag, sbml_vcml_ns);
// Element compartmentElement = new Element(XMLTags.OutsideCompartmentTag, sbml_vcml_ns);
// compartmentElement.setAttribute(XMLTags.NameAttrTag, TokenMangler.mangleToSName(parentStructure.getName()));
// sbmlImportRelatedElement.addContent(compartmentElement);
// }
// Get annotation (RDF and non-RDF) for reactionStep from SBMLAnnotationUtils
sbmlAnnotationUtil.writeAnnotation(vcStructures[i], sbmlCompartment, sbmlImportRelatedElement);
// Now set notes,
sbmlAnnotationUtil.writeNotes(vcStructures[i], sbmlCompartment);
}
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class SBMLImporter method addCompartments.
protected void addCompartments(VCMetaData metaData) {
if (sbmlModel == null) {
throw new SBMLImportException("SBML model is NULL");
}
ListOf listofCompartments = sbmlModel.getListOfCompartments();
if (listofCompartments == null) {
throw new SBMLImportException("Cannot have 0 compartments in model");
}
// Using a vector here - since there can be SBML models with only
// features and no membranes.
// Hence keeping the datastructure flexible.
List<Structure> structList = new ArrayList<Structure>();
java.util.HashMap<String, Structure> structureNameMap = new java.util.HashMap<String, Structure>();
try {
int structIndx = 0;
// First pass - create the structures
for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
org.sbml.jsbml.Compartment compartment = (org.sbml.jsbml.Compartment) listofCompartments.get(i);
String compartmentName = compartment.getId();
if (!compartment.isSetSpatialDimensions() || compartment.getSpatialDimensions() == 3) {
Feature feature = new Feature(compartmentName);
structList.add(structIndx, feature);
structureNameMap.put(compartmentName, feature);
} else if (compartment.getSpatialDimensions() == 2) {
// spatial dimensions is set (see clause above)
Membrane membrane = new Membrane(compartmentName);
structList.add(structIndx, membrane);
structureNameMap.put(compartmentName, membrane);
} else {
logger.sendMessage(VCLogger.Priority.HighPriority, VCLogger.ErrorType.CompartmentError, "Cannot deal with spatial dimension : " + compartment.getSpatialDimensions() + " for compartments at this time.");
throw new SBMLImportException("Cannot deal with spatial dimension : " + compartment.getSpatialDimensions() + " for compartments at this time");
}
structIndx++;
sbmlAnnotationUtil.readAnnotation(structList.get(i), compartment);
sbmlAnnotationUtil.readNotes(structList.get(i), compartment);
}
// Second pass - connect the structures
Model model = vcBioModel.getSimulationContext(0).getModel();
for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
org.sbml.jsbml.Compartment sbmlCompartment = (org.sbml.jsbml.Compartment) listofCompartments.get(i);
String outsideCompartmentId = null;
if (sbmlCompartment.getOutside() != null && sbmlCompartment.getOutside().length() > 0) {
// compartment.getOutside returns the Sid of the 'outside'
// compartment, so get the compartment from model.
outsideCompartmentId = sbmlCompartment.getOutside();
} else {
Element sbmlImportRelatedElement = sbmlAnnotationUtil.readVCellSpecificAnnotation(sbmlCompartment);
if (sbmlImportRelatedElement != null) {
Element embeddedVCellElement = sbmlImportRelatedElement.getChild(OUTSIDE_COMP_NAME, Namespace.getNamespace(SBMLUtils.SBML_VCELL_NS));
if (embeddedVCellElement != null) {
outsideCompartmentId = embeddedVCellElement.getAttributeValue(XMLTags.NameTag);
}
}
}
if (outsideCompartmentId != null) {
Compartment outsideCompartment = sbmlModel.getCompartment(outsideCompartmentId);
Structure outsideStructure = (Structure) structureNameMap.get(outsideCompartment.getId());
Structure struct = (Structure) structureNameMap.get(sbmlCompartment.getId());
struct.setSbmlParentStructure(outsideStructure);
}
}
// set the structures in vc vcBioModel.getSimulationContext(0)
Structure[] structures = structList.toArray(new Structure[structList.size()]);
model.setStructures(structures);
// Third pass thro' the list of compartments : set the sizes on the
// structureMappings - can be done only after setting
// the structures on vcBioModel.getSimulationContext(0).
boolean allSizesSet = true;
for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
org.sbml.jsbml.Compartment compartment = (org.sbml.jsbml.Compartment) listofCompartments.get(i);
String compartmentName = compartment.getId();
if (!compartment.isSetSize()) {
// logger.sendMessage(VCLogger.Priority.MediumPriority,
// TranslationMessage.COMPARTMENT_ERROR,
// "compartment "+compartmentName+" size is not set in SBML document.");
allSizesSet = false;
} else {
double size = compartment.getSize();
// Check if size is specified by a rule
Expression sizeExpr = getValueFromAssignmentRule(compartmentName);
if (sizeExpr != null && !sizeExpr.isNumeric()) {
// We are NOT handling compartment sizes with assignment
// rules/initial Assignments that are NON-numeric at
// this time ...
logger.sendMessage(VCLogger.Priority.HighPriority, VCLogger.ErrorType.CompartmentError, "compartment " + compartmentName + " size has an assignment rule which is not a numeric value, cannot handle it at this time.");
}
// check if it is specified by initial assignment
if (sizeExpr == null) {
InitialAssignment compInitAssgnment = sbmlModel.getInitialAssignment(compartmentName);
if (compInitAssgnment != null) {
sizeExpr = getExpressionFromFormula(compInitAssgnment.getMath());
}
}
if (sizeExpr != null && !sizeExpr.isNumeric()) {
// We are NOT handling compartment sizes with assignment
// rules/initial Assignments that are NON-numeric at
// this time ...
logger.sendMessage(VCLogger.Priority.HighPriority, VCLogger.ErrorType.CompartmentError, "compartment " + compartmentName + " size has an initial assignment which is not a numeric value, cannot handle it at this time.");
}
// from 'size' attribute,
if (sizeExpr == null) {
sizeExpr = new Expression(size);
}
// Now set the size of the compartment.
Structure struct = model.getStructure(compartmentName);
StructureMapping.StructureMappingParameter mappingParam = vcBioModel.getSimulationContext(0).getGeometryContext().getStructureMapping(struct).getSizeParameter();
mappingParam.setExpression(sizeExpr);
}
}
// size is set
if (allSizesSet) {
StructureSizeSolver.updateRelativeStructureSizes(vcBioModel.getSimulationContext(0));
}
} catch (Exception e) {
e.printStackTrace(System.out);
throw new SBMLImportException("Error adding Feature to vcModel " + e.getMessage(), e);
}
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class SBMLImporter method addGeometry.
protected void addGeometry() {
// get a Geometry object via SpatialModelPlugin object.
org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = getSbmlGeometry();
if (sbmlGeometry == null) {
return;
}
int dimension = 0;
Origin vcOrigin = null;
Extent vcExtent = null;
{
// local code block
// get a CoordComponent object via the Geometry object.
ListOf<CoordinateComponent> listOfCoordComps = sbmlGeometry.getListOfCoordinateComponents();
if (listOfCoordComps == null) {
throw new RuntimeException("Cannot have 0 coordinate compartments in geometry");
}
// coord component
double ox = 0.0;
double oy = 0.0;
double oz = 0.0;
double ex = 1.0;
double ey = 1.0;
double ez = 1.0;
for (CoordinateComponent coordComponent : listOfCoordComps) {
double minValue = coordComponent.getBoundaryMinimum().getValue();
double maxValue = coordComponent.getBoundaryMaximum().getValue();
switch(coordComponent.getType()) {
case cartesianX:
{
ox = minValue;
ex = maxValue - minValue;
break;
}
case cartesianY:
{
oy = minValue;
ey = maxValue - minValue;
break;
}
case cartesianZ:
{
oz = minValue;
ez = maxValue - minValue;
break;
}
}
dimension++;
}
vcOrigin = new Origin(ox, oy, oz);
vcExtent = new Extent(ex, ey, ez);
}
// from geometry definition, find out which type of geometry : image or
// analytic or CSG
AnalyticGeometry analyticGeometryDefinition = null;
CSGeometry csGeometry = null;
SampledFieldGeometry segmentedSampledFieldGeometry = null;
SampledFieldGeometry distanceMapSampledFieldGeometry = null;
ParametricGeometry parametricGeometry = null;
for (int i = 0; i < sbmlGeometry.getListOfGeometryDefinitions().size(); i++) {
GeometryDefinition gd_temp = sbmlGeometry.getListOfGeometryDefinitions().get(i);
if (!gd_temp.isSetIsActive()) {
continue;
}
if (gd_temp instanceof AnalyticGeometry) {
analyticGeometryDefinition = (AnalyticGeometry) gd_temp;
} else if (gd_temp instanceof SampledFieldGeometry) {
SampledFieldGeometry sfg = (SampledFieldGeometry) gd_temp;
String sfn = sfg.getSampledField();
ListOf<SampledField> sampledFields = sbmlGeometry.getListOfSampledFields();
if (sampledFields.size() > 1) {
throw new RuntimeException("only one sampled field supported");
}
InterpolationKind ik = sampledFields.get(0).getInterpolationType();
switch(ik) {
case linear:
distanceMapSampledFieldGeometry = sfg;
break;
case nearestneighbor:
segmentedSampledFieldGeometry = sfg;
break;
default:
lg.warn("Unsupported " + sampledFields.get(0).getName() + " interpolation type " + ik);
}
} else if (gd_temp instanceof CSGeometry) {
csGeometry = (CSGeometry) gd_temp;
} else if (gd_temp instanceof ParametricGeometry) {
parametricGeometry = (ParametricGeometry) gd_temp;
} else {
throw new RuntimeException("unsupported geometry definition type " + gd_temp.getClass().getSimpleName());
}
}
if (analyticGeometryDefinition == null && segmentedSampledFieldGeometry == null && distanceMapSampledFieldGeometry == null && csGeometry == null) {
throw new SBMLImportException("VCell supports only Analytic, Image based (segmentd or distance map) or Constructed Solid Geometry at this time.");
}
GeometryDefinition selectedGeometryDefinition = null;
if (csGeometry != null) {
selectedGeometryDefinition = csGeometry;
} else if (analyticGeometryDefinition != null) {
selectedGeometryDefinition = analyticGeometryDefinition;
} else if (segmentedSampledFieldGeometry != null) {
selectedGeometryDefinition = segmentedSampledFieldGeometry;
} else if (distanceMapSampledFieldGeometry != null) {
selectedGeometryDefinition = distanceMapSampledFieldGeometry;
} else if (parametricGeometry != null) {
selectedGeometryDefinition = parametricGeometry;
} else {
throw new SBMLImportException("no geometry definition found");
}
Geometry vcGeometry = null;
if (selectedGeometryDefinition == analyticGeometryDefinition || selectedGeometryDefinition == csGeometry) {
vcGeometry = new Geometry("spatialGeom", dimension);
} else if (selectedGeometryDefinition == distanceMapSampledFieldGeometry || selectedGeometryDefinition == segmentedSampledFieldGeometry) {
SampledFieldGeometry sfg = (SampledFieldGeometry) selectedGeometryDefinition;
// get image from sampledFieldGeometry
// get a sampledVol object via the listOfSampledVol (from
// SampledGeometry) object.
// gcw gcw gcw
String sfn = sfg.getSampledField();
SampledField sf = null;
for (SampledField sampledField : sbmlGeometry.getListOfSampledFields()) {
if (sampledField.getSpatialId().equals(sfn)) {
sf = sampledField;
}
}
int numX = sf.getNumSamples1();
int numY = sf.getNumSamples2();
int numZ = sf.getNumSamples3();
int[] samples = new int[sf.getSamplesLength()];
StringTokenizer tokens = new StringTokenizer(sf.getSamples(), " ");
int count = 0;
while (tokens.hasMoreTokens()) {
int sample = Integer.parseInt(tokens.nextToken());
samples[count++] = sample;
}
byte[] imageInBytes = new byte[samples.length];
if (selectedGeometryDefinition == distanceMapSampledFieldGeometry) {
//
for (int i = 0; i < imageInBytes.length; i++) {
// if (interpolation(samples[i])<0){
if (samples[i] < 0) {
imageInBytes[i] = -1;
} else {
imageInBytes[i] = 1;
}
}
} else {
for (int i = 0; i < imageInBytes.length; i++) {
imageInBytes[i] = (byte) samples[i];
}
}
try {
// System.out.println("ident " + sf.getId() + " " + sf.getName());
VCImage vcImage = null;
CompressionKind ck = sf.getCompression();
DataKind dk = sf.getDataType();
if (ck == CompressionKind.deflated) {
vcImage = new VCImageCompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
} else {
switch(dk) {
case UINT8:
case UINT16:
case UINT32:
vcImage = new VCImageUncompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
default:
}
}
if (vcImage == null) {
throw new SbmlException("Unsupported type combination " + ck + ", " + dk + " for sampled field " + sf.getName());
}
vcImage.setName(sf.getId());
ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
final int numSampledVols = sampledVolumes.size();
if (numSampledVols == 0) {
throw new RuntimeException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
}
// check to see if values are uniquely integer , add set up scaling if necessary
double scaleFactor = checkPixelScaling(sampledVolumes, 1);
if (scaleFactor != 1) {
double checkScaleFactor = checkPixelScaling(sampledVolumes, scaleFactor);
VCAssert.assertTrue(checkScaleFactor != scaleFactor, "Scale factor check failed");
}
VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
// get pixel classes for geometry
for (int i = 0; i < numSampledVols; i++) {
SampledVolume sVol = sampledVolumes.get(i);
// from subVolume, get pixelClass?
final int scaled = (int) (scaleFactor * sVol.getSampledValue());
vcpixelClasses[i] = new VCPixelClass(null, sVol.getDomainType(), scaled);
}
vcImage.setPixelClasses(vcpixelClasses);
// now create image geometry
vcGeometry = new Geometry("spatialGeom", vcImage);
} catch (Exception e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to create image from SampledFieldGeometry : " + e.getMessage());
}
}
GeometrySpec vcGeometrySpec = vcGeometry.getGeometrySpec();
vcGeometrySpec.setOrigin(vcOrigin);
try {
vcGeometrySpec.setExtent(vcExtent);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new SBMLImportException("Unable to set extent on VC geometry : " + e.getMessage(), e);
}
// get listOfDomainTypes via the Geometry object.
ListOf<DomainType> listOfDomainTypes = sbmlGeometry.getListOfDomainTypes();
if (listOfDomainTypes == null || listOfDomainTypes.size() < 1) {
throw new SBMLImportException("Cannot have 0 domainTypes in geometry");
}
// get a listOfDomains via the Geometry object.
ListOf<Domain> listOfDomains = sbmlGeometry.getListOfDomains();
if (listOfDomains == null || listOfDomains.size() < 1) {
throw new SBMLImportException("Cannot have 0 domains in geometry");
}
// ListOfGeometryDefinitions listOfGeomDefns =
// sbmlGeometry.getListOfGeometryDefinitions();
// if ((listOfGeomDefns == null) ||
// (sbmlGeometry.getNumGeometryDefinitions() > 1)) {
// throw new
// RuntimeException("Can have only 1 geometry definition in geometry");
// }
// use the boolean bAnalytic to create the right kind of subvolume.
// First match the somVol=domainTypes for spDim=3. Deal witl spDim=2
// afterwards.
GeometrySurfaceDescription vcGsd = vcGeometry.getGeometrySurfaceDescription();
Vector<DomainType> surfaceClassDomainTypesVector = new Vector<DomainType>();
try {
for (DomainType dt : listOfDomainTypes) {
if (dt.getSpatialDimensions() == 3) {
// subvolume
if (selectedGeometryDefinition == analyticGeometryDefinition) {
// will set expression later - when reading in Analytic
// Volumes in GeometryDefinition
vcGeometrySpec.addSubVolume(new AnalyticSubVolume(dt.getId(), new Expression(1.0)));
} else {
// add SubVolumes later for CSG and Image-based
}
} else if (dt.getSpatialDimensions() == 2) {
surfaceClassDomainTypesVector.add(dt);
}
}
// analytic vol is needed to get the expression for subVols
if (selectedGeometryDefinition == analyticGeometryDefinition) {
// get an analyticVol object via the listOfAnalyticVol (from
// AnalyticGeometry) object.
ListOf<AnalyticVolume> aVolumes = analyticGeometryDefinition.getListOfAnalyticVolumes();
if (aVolumes.size() < 1) {
throw new SBMLImportException("Cannot have 0 Analytic volumes in analytic geometry");
}
for (AnalyticVolume analyticVol : aVolumes) {
// get subVol from VC geometry using analyticVol spatialId;
// set its expr using analyticVol's math.
SubVolume vcSubvolume = vcGeometrySpec.getSubVolume(analyticVol.getDomainType());
CastInfo<AnalyticSubVolume> ci = BeanUtils.attemptCast(AnalyticSubVolume.class, vcSubvolume);
if (!ci.isGood()) {
throw new RuntimeException("analytic volume '" + analyticVol.getId() + "' does not map to any VC subvolume.");
}
AnalyticSubVolume asv = ci.get();
try {
Expression subVolExpr = getExpressionFromFormula(analyticVol.getMath());
asv.setExpression(subVolExpr);
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new SBMLImportException("Unable to set expression on subVolume '" + asv.getName() + "'. " + e.getMessage(), e);
}
}
}
SampledFieldGeometry sfg = BeanUtils.downcast(SampledFieldGeometry.class, selectedGeometryDefinition);
if (sfg != null) {
ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
int numSampledVols = sampledVolumes.size();
if (numSampledVols == 0) {
throw new SBMLImportException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
}
VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
ImageSubVolume[] vcImageSubVols = new ImageSubVolume[numSampledVols];
// get pixel classes for geometry
int idx = 0;
for (SampledVolume sVol : sampledVolumes) {
// from subVolume, get pixelClass?
final String name = sVol.getDomainType();
final int pixelValue = SBMLUtils.ignoreZeroFraction(sVol.getSampledValue());
VCPixelClass pc = new VCPixelClass(null, name, pixelValue);
vcpixelClasses[idx] = pc;
// Create the new Image SubVolume - use index of this for
// loop as 'handle' for ImageSubVol?
ImageSubVolume isv = new ImageSubVolume(null, pc, idx);
isv.setName(name);
vcImageSubVols[idx++] = isv;
}
vcGeometry.getGeometrySpec().setSubVolumes(vcImageSubVols);
}
if (selectedGeometryDefinition == csGeometry) {
ListOf<org.sbml.jsbml.ext.spatial.CSGObject> listOfcsgObjs = csGeometry.getListOfCSGObjects();
ArrayList<org.sbml.jsbml.ext.spatial.CSGObject> sbmlCSGs = new ArrayList<org.sbml.jsbml.ext.spatial.CSGObject>(listOfcsgObjs);
// we want the CSGObj with highest ordinal to be the first
// element in the CSG subvols array.
Collections.sort(sbmlCSGs, new Comparator<org.sbml.jsbml.ext.spatial.CSGObject>() {
@Override
public int compare(org.sbml.jsbml.ext.spatial.CSGObject lhs, org.sbml.jsbml.ext.spatial.CSGObject rhs) {
// minus one to reverse sort
return -1 * Integer.compare(lhs.getOrdinal(), rhs.getOrdinal());
}
});
int n = sbmlCSGs.size();
CSGObject[] vcCSGSubVolumes = new CSGObject[n];
for (int i = 0; i < n; i++) {
org.sbml.jsbml.ext.spatial.CSGObject sbmlCSGObject = sbmlCSGs.get(i);
CSGObject vcellCSGObject = new CSGObject(null, sbmlCSGObject.getDomainType(), i);
vcellCSGObject.setRoot(getVCellCSGNode(sbmlCSGObject.getCSGNode()));
}
vcGeometry.getGeometrySpec().setSubVolumes(vcCSGSubVolumes);
}
// Call geom.geomSurfDesc.updateAll() to automatically generate
// surface classes.
// vcGsd.updateAll();
vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, true);
} catch (Exception e) {
e.printStackTrace(System.out);
throw new SBMLImportException("Unable to create VC subVolumes from SBML domainTypes : " + e.getMessage(), e);
}
// should now map each SBML domain to right VC geometric region.
GeometricRegion[] vcGeomRegions = vcGsd.getGeometricRegions();
ISize sampleSize = vcGsd.getVolumeSampleSize();
RegionInfo[] regionInfos = vcGsd.getRegionImage().getRegionInfos();
int numX = sampleSize.getX();
int numY = sampleSize.getY();
int numZ = sampleSize.getZ();
double ox = vcOrigin.getX();
double oy = vcOrigin.getY();
double oz = vcOrigin.getZ();
for (Domain domain : listOfDomains) {
String domainType = domain.getDomainType();
InteriorPoint interiorPt = domain.getListOfInteriorPoints().get(0);
if (interiorPt == null) {
DomainType currDomainType = null;
for (DomainType dt : sbmlGeometry.getListOfDomainTypes()) {
if (dt.getSpatialId().equals(domainType)) {
currDomainType = dt;
}
}
if (currDomainType.getSpatialDimensions() == 2) {
continue;
}
}
Coordinate sbmlInteriorPtCoord = new Coordinate(interiorPt.getCoord1(), interiorPt.getCoord2(), interiorPt.getCoord3());
for (int j = 0; j < vcGeomRegions.length; j++) {
if (vcGeomRegions[j] instanceof VolumeGeometricRegion) {
int regionID = ((VolumeGeometricRegion) vcGeomRegions[j]).getRegionID();
for (int k = 0; k < regionInfos.length; k++) {
// (using gemoRegion regionID).
if (regionInfos[k].getRegionIndex() == regionID) {
int volIndx = 0;
Coordinate nearestPtCoord = null;
double minDistance = Double.MAX_VALUE;
// represented by SBML 'domain[i]'.
for (int z = 0; z < numZ; z++) {
for (int y = 0; y < numY; y++) {
for (int x = 0; x < numX; x++) {
if (regionInfos[k].isIndexInRegion(volIndx)) {
double unit_z = (numZ > 1) ? ((double) z) / (numZ - 1) : 0.5;
double coordZ = oz + vcExtent.getZ() * unit_z;
double unit_y = (numY > 1) ? ((double) y) / (numY - 1) : 0.5;
double coordY = oy + vcExtent.getY() * unit_y;
double unit_x = (numX > 1) ? ((double) x) / (numX - 1) : 0.5;
double coordX = ox + vcExtent.getX() * unit_x;
// for now, find the shortest dist
// coord. Can refine algo later.
Coordinate vcCoord = new Coordinate(coordX, coordY, coordZ);
double distance = sbmlInteriorPtCoord.distanceTo(vcCoord);
if (distance < minDistance) {
minDistance = distance;
nearestPtCoord = vcCoord;
}
}
volIndx++;
}
// end - for x
}
// end - for y
}
// with domain name
if (nearestPtCoord != null) {
GeometryClass geomClassSBML = vcGeometry.getGeometryClass(domainType);
// we know vcGeometryReg[j] is a VolGeomRegion
GeometryClass geomClassVC = ((VolumeGeometricRegion) vcGeomRegions[j]).getSubVolume();
if (geomClassSBML.compareEqual(geomClassVC)) {
vcGeomRegions[j].setName(domain.getId());
}
}
}
// end if (regInfoIndx = regId)
}
// end - for regInfo
}
}
// end for - vcGeomRegions
}
// deal with surfaceClass:spDim2-domainTypes
for (int i = 0; i < surfaceClassDomainTypesVector.size(); i++) {
DomainType surfaceClassDomainType = surfaceClassDomainTypesVector.elementAt(i);
// 'surfaceClassDomainType'
for (Domain d : listOfDomains) {
if (d.getDomainType().equals(surfaceClassDomainType.getId())) {
// get the adjacent domains of this 'surface' domain
// (surface domain + its 2 adj vol domains)
Set<Domain> adjacentDomainsSet = getAssociatedAdjacentDomains(sbmlGeometry, d);
// get the domain types of the adjacent domains in SBML and
// store the corresponding subVol counterparts from VC for
// adj vol domains
Vector<SubVolume> adjacentSubVolumesVector = new Vector<SubVolume>();
Vector<VolumeGeometricRegion> adjVolGeomRegionsVector = new Vector<VolumeGeometricRegion>();
Iterator<Domain> iterator = adjacentDomainsSet.iterator();
while (iterator.hasNext()) {
Domain dom = iterator.next();
DomainType dt = getBySpatialID(sbmlGeometry.getListOfDomainTypes(), dom.getDomainType());
if (dt.getSpatialDimensions() == 3) {
// for domain type with sp. dim = 3, get
// correspoinding subVol from VC geometry.
GeometryClass gc = vcGeometry.getGeometryClass(dt.getId());
adjacentSubVolumesVector.add((SubVolume) gc);
// store volGeomRegions corresponding to this (vol)
// geomClass in adjVolGeomRegionsVector : this
// should return ONLY 1 region for subVol.
GeometricRegion[] geomRegion = vcGsd.getGeometricRegions(gc);
adjVolGeomRegionsVector.add((VolumeGeometricRegion) geomRegion[0]);
}
}
// there should be only 2 subVols in this vector
if (adjacentSubVolumesVector.size() != 2) {
throw new RuntimeException("Cannot have more or less than 2 subvolumes that are adjacent to surface (membrane) '" + d.getId() + "'");
}
// get the surface class with these 2 adj subVols. Set its
// name to that of 'surfaceClassDomainType'
SurfaceClass surfacClass = vcGsd.getSurfaceClass(adjacentSubVolumesVector.get(0), adjacentSubVolumesVector.get(1));
surfacClass.setName(surfaceClassDomainType.getSpatialId());
// get surfaceGeometricRegion that has adjVolGeomRegions as
// its adjacent vol geom regions and set its name from
// domain 'd'
SurfaceGeometricRegion surfaceGeomRegion = getAssociatedSurfaceGeometricRegion(vcGsd, adjVolGeomRegionsVector);
if (surfaceGeomRegion != null) {
surfaceGeomRegion.setName(d.getId());
}
}
// end if - domain.domainType == surfaceClassDomainType
}
// end for - numDomains
}
// structureMappings in VC from compartmentMappings in SBML
try {
// set geometry first and then set structureMappings?
vcBioModel.getSimulationContext(0).setGeometry(vcGeometry);
// update simContextName ...
vcBioModel.getSimulationContext(0).setName(vcBioModel.getSimulationContext(0).getName() + "_" + vcGeometry.getName());
Model vcModel = vcBioModel.getSimulationContext(0).getModel();
ModelUnitSystem vcModelUnitSystem = vcModel.getUnitSystem();
Vector<StructureMapping> structMappingsVector = new Vector<StructureMapping>();
SpatialCompartmentPlugin cplugin = null;
for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
Compartment c = sbmlModel.getCompartment(i);
String cname = c.getName();
cplugin = (SpatialCompartmentPlugin) c.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
CompartmentMapping compMapping = cplugin.getCompartmentMapping();
if (compMapping != null) {
// final String id = compMapping.getId();
// final String name = compMapping.getName();
CastInfo<Structure> ci = SBMLHelper.getTypedStructure(Structure.class, vcModel, cname);
if (ci.isGood()) {
Structure struct = ci.get();
String domainType = compMapping.getDomainType();
GeometryClass geometryClass = vcGeometry.getGeometryClass(domainType);
double unitSize = compMapping.getUnitSize();
Feature feat = BeanUtils.downcast(Feature.class, struct);
if (feat != null) {
FeatureMapping featureMapping = new FeatureMapping(feat, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
featureMapping.setGeometryClass(geometryClass);
if (geometryClass instanceof SubVolume) {
featureMapping.getVolumePerUnitVolumeParameter().setExpression(new Expression(unitSize));
} else if (geometryClass instanceof SurfaceClass) {
featureMapping.getVolumePerUnitAreaParameter().setExpression(new Expression(unitSize));
}
structMappingsVector.add(featureMapping);
} else if (struct instanceof Membrane) {
MembraneMapping membraneMapping = new MembraneMapping((Membrane) struct, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
membraneMapping.setGeometryClass(geometryClass);
if (geometryClass instanceof SubVolume) {
membraneMapping.getAreaPerUnitVolumeParameter().setExpression(new Expression(unitSize));
} else if (geometryClass instanceof SurfaceClass) {
membraneMapping.getAreaPerUnitAreaParameter().setExpression(new Expression(unitSize));
}
structMappingsVector.add(membraneMapping);
}
}
}
}
StructureMapping[] structMappings = structMappingsVector.toArray(new StructureMapping[0]);
vcBioModel.getSimulationContext(0).getGeometryContext().setStructureMappings(structMappings);
// if type from SBML parameter Boundary Condn is not the same as the
// boundary type of the
// structureMapping of structure of paramSpContext, set the boundary
// condn type of the structureMapping
// to the value of 'type' from SBML parameter Boundary Condn.
ListOf<Parameter> listOfGlobalParams = sbmlModel.getListOfParameters();
for (Parameter sbmlGlobalParam : sbmlModel.getListOfParameters()) {
SpatialParameterPlugin spplugin = (SpatialParameterPlugin) sbmlGlobalParam.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
ParameterType paramType = spplugin.getParamType();
if (!(paramType instanceof BoundaryCondition)) {
continue;
}
BoundaryCondition bCondn = (BoundaryCondition) paramType;
if (bCondn.isSetVariable()) {
// get the var of boundaryCondn; find appropriate spContext
// in vcell;
SpeciesContext paramSpContext = vcBioModel.getSimulationContext(0).getModel().getSpeciesContext(bCondn.getVariable());
if (paramSpContext != null) {
Structure s = paramSpContext.getStructure();
StructureMapping sm = vcBioModel.getSimulationContext(0).getGeometryContext().getStructureMapping(s);
if (sm != null) {
BoundaryConditionType bct = null;
switch(bCondn.getType()) {
case Dirichlet:
{
bct = BoundaryConditionType.DIRICHLET;
break;
}
case Neumann:
{
bct = BoundaryConditionType.NEUMANN;
break;
}
case Robin_inwardNormalGradientCoefficient:
case Robin_sum:
case Robin_valueCoefficient:
default:
throw new RuntimeException("boundary condition type " + bCondn.getType().name() + " not supported");
}
for (CoordinateComponent coordComp : getSbmlGeometry().getListOfCoordinateComponents()) {
if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMinimum().getSpatialId())) {
switch(coordComp.getType()) {
case cartesianX:
{
sm.setBoundaryConditionTypeXm(bct);
}
case cartesianY:
{
sm.setBoundaryConditionTypeYm(bct);
}
case cartesianZ:
{
sm.setBoundaryConditionTypeZm(bct);
}
}
}
if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMaximum().getSpatialId())) {
switch(coordComp.getType()) {
case cartesianX:
{
sm.setBoundaryConditionTypeXm(bct);
}
case cartesianY:
{
sm.setBoundaryConditionTypeYm(bct);
}
case cartesianZ:
{
sm.setBoundaryConditionTypeZm(bct);
}
}
}
}
} else // sm != null
{
logger.sendMessage(VCLogger.Priority.MediumPriority, VCLogger.ErrorType.OverallWarning, "No structure " + s.getName() + " requested by species context " + paramSpContext.getName());
}
}
// end if (paramSpContext != null)
}
// end if (bCondn.isSetVar())
}
// end for (sbmlModel.numParams)
vcBioModel.getSimulationContext(0).getGeometryContext().refreshStructureMappings();
vcBioModel.getSimulationContext(0).refreshSpatialObjects();
} catch (Exception e) {
e.printStackTrace(System.out);
throw new SBMLImportException("Unable to create VC structureMappings from SBML compartment mappings : " + e.getMessage(), e);
}
}
Aggregations