Examples with FeatureMapping cbit.vcell.mapping.FeatureMapping used on opensource projects

Example 16 with FeatureMapping

use of cbit.vcell.mapping.FeatureMapping in project vcell by virtualcell.

the class XmlReader method getFeatureMapping.

/**
 * This method retuns a FeatureMapping object from a XML representation.
 * Creation date: (5/7/2001 4:12:03 PM)
 * @return cbit.vcell.mapping.FeatureMapping
 * @param param org.jdom.Element
 */
private FeatureMapping getFeatureMapping(Element param, SimulationContext simulationContext) throws XmlParseException {
    // Retrieve attributes
    String featurename = unMangle(param.getAttributeValue(XMLTags.FeatureAttrTag));
    String geometryClassName = param.getAttributeValue(XMLTags.SubVolumeAttrTag);
    if (geometryClassName != null) {
        geometryClassName = unMangle(geometryClassName);
    } else {
        geometryClassName = param.getAttributeValue(XMLTags.GeometryClassAttrTag);
        if (geometryClassName != null) {
            geometryClassName = unMangle(geometryClassName);
        }
    }
    Feature featureref = (Feature) simulationContext.getModel().getStructure(featurename);
    if (featureref == null) {
        throw new XmlParseException("The Feature " + featurename + " could not be resolved!");
    }
    // *** Create new Feature Mapping ****
    FeatureMapping feamap = new FeatureMapping(featureref, simulationContext, simulationContext.getModel().getUnitSystem());
    // Set Size
    if (param.getAttributeValue(XMLTags.SizeTag) != null) {
        String size = unMangle(param.getAttributeValue(XMLTags.SizeTag));
        try {
            feamap.getSizeParameter().setExpression(unMangleExpression(size));
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            throw new XmlParseException("An expressionException was fired when setting the size Expression " + size + " to a featureMapping!", e);
        }
    } else {
        try {
            feamap.getSizeParameter().setExpression(null);
        } catch (Exception e) {
            e.printStackTrace();
            throw new RuntimeException("unexpected exception while setting structure size: " + e.getMessage());
        }
    }
    // Set Volume/unit_area if it exists
    if (param.getAttributeValue(XMLTags.VolumePerUnitAreaTag) != null) {
        String volPerUnitArea = unMangle(param.getAttributeValue(XMLTags.VolumePerUnitAreaTag));
        try {
            feamap.getVolumePerUnitAreaParameter().setExpression(unMangleExpression(volPerUnitArea));
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            throw new XmlParseException("An expressionException was fired when setting the VolumePerUnitArea Expression " + volPerUnitArea + " to a featureMapping!", e);
        }
    }
    // Set Volume/unitVol if it exists
    if (param.getAttributeValue(XMLTags.VolumePerUnitVolumeTag) != null) {
        String volPerUnitVol = unMangle(param.getAttributeValue(XMLTags.VolumePerUnitVolumeTag));
        try {
            feamap.getVolumePerUnitVolumeParameter().setExpression(unMangleExpression(volPerUnitVol));
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            throw new XmlParseException("An expressionException was fired when setting the size Expression " + volPerUnitVol + " to a featureMapping!", e);
        }
    }
    if (geometryClassName != null) {
        GeometryClass[] geometryClasses = simulationContext.getGeometry().getGeometryClasses();
        for (int i = 0; i < geometryClasses.length; i++) {
            if (geometryClasses[i].getName().equals(geometryClassName)) {
                try {
                    feamap.setGeometryClass(geometryClasses[i]);
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new XmlParseException("A propertyVetoException was fired when trying to set the subvolume or surface " + geometryClassName + " to a MembraneMapping!", e);
                }
            }
        }
    }
    // Set Boundary conditions
    Element tempElement = param.getChild(XMLTags.BoundariesTypesTag, vcNamespace);
    // Xm
    String temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueXm);
    BoundaryConditionType bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeXm(bct);
    // Xp
    temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueXp);
    bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeXp(bct);
    // Ym
    temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueYm);
    bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeYm(bct);
    // Yp
    temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueYp);
    bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeYp(bct);
    // Zm
    temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueZm);
    bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeZm(bct);
    // Zp
    temp = tempElement.getAttributeValue(XMLTags.BoundaryAttrValueZp);
    bct = new BoundaryConditionType(temp);
    feamap.setBoundaryConditionTypeZp(bct);
    return feamap;
}

Example 17 with FeatureMapping

use of cbit.vcell.mapping.FeatureMapping in project vcell by virtualcell.

the class FRAPStudy method createNewSimBioModel.

public static BioModel createNewSimBioModel(FRAPStudy sourceFrapStudy, Parameter[] params, TimeStep tStep, KeyValue simKey, User owner, int startingIndexForRecovery) throws Exception {
    if (owner == null) {
        throw new Exception("Owner is not defined");
    }
    ROI cellROI_2D = sourceFrapStudy.getFrapData().getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name());
    double df = params[FRAPModel.INDEX_PRIMARY_DIFF_RATE].getInitialGuess();
    double ff = params[FRAPModel.INDEX_PRIMARY_FRACTION].getInitialGuess();
    double bwmRate = params[FRAPModel.INDEX_BLEACH_MONITOR_RATE].getInitialGuess();
    double dc = 0;
    double fc = 0;
    double bs = 0;
    double onRate = 0;
    double offRate = 0;
    if (params.length == FRAPModel.NUM_MODEL_PARAMETERS_TWO_DIFF) {
        dc = params[FRAPModel.INDEX_SECONDARY_DIFF_RATE].getInitialGuess();
        fc = params[FRAPModel.INDEX_SECONDARY_FRACTION].getInitialGuess();
    } else if (params.length == FRAPModel.NUM_MODEL_PARAMETERS_BINDING) {
        dc = params[FRAPModel.INDEX_SECONDARY_DIFF_RATE].getInitialGuess();
        fc = params[FRAPModel.INDEX_SECONDARY_FRACTION].getInitialGuess();
        bs = params[FRAPModel.INDEX_BINDING_SITE_CONCENTRATION].getInitialGuess();
        onRate = params[FRAPModel.INDEX_ON_RATE].getInitialGuess();
        offRate = params[FRAPModel.INDEX_OFF_RATE].getInitialGuess();
    }
    // immobile fraction
    double fimm = 1 - ff - fc;
    if (fimm < FRAPOptimizationUtils.epsilon && fimm > (0 - FRAPOptimizationUtils.epsilon)) {
        fimm = 0;
    }
    if (fimm < (1 + FRAPOptimizationUtils.epsilon) && fimm > (1 - FRAPOptimizationUtils.epsilon)) {
        fimm = 1;
    }
    Extent extent = sourceFrapStudy.getFrapData().getImageDataset().getExtent();
    double[] timeStamps = sourceFrapStudy.getFrapData().getImageDataset().getImageTimeStamps();
    TimeBounds timeBounds = new TimeBounds(0.0, timeStamps[timeStamps.length - 1] - timeStamps[startingIndexForRecovery]);
    double timeStepVal = timeStamps[startingIndexForRecovery + 1] - timeStamps[startingIndexForRecovery];
    int numX = cellROI_2D.getRoiImages()[0].getNumX();
    int numY = cellROI_2D.getRoiImages()[0].getNumY();
    int numZ = cellROI_2D.getRoiImages().length;
    short[] shortPixels = cellROI_2D.getRoiImages()[0].getPixels();
    byte[] bytePixels = new byte[numX * numY * numZ];
    final byte EXTRACELLULAR_PIXVAL = 0;
    final byte CYTOSOL_PIXVAL = 1;
    for (int i = 0; i < bytePixels.length; i++) {
        if (shortPixels[i] != 0) {
            bytePixels[i] = CYTOSOL_PIXVAL;
        }
    }
    VCImage maskImage;
    try {
        maskImage = new VCImageUncompressed(null, bytePixels, extent, numX, numY, numZ);
    } catch (ImageException e) {
        e.printStackTrace();
        throw new RuntimeException("failed to create mask image for geometry");
    }
    Geometry geometry = new Geometry("geometry", maskImage);
    if (geometry.getGeometrySpec().getNumSubVolumes() != 2) {
        throw new Exception("Cell ROI has no ExtraCellular.");
    }
    int subVolume0PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(0)).getPixelValue();
    geometry.getGeometrySpec().getSubVolume(0).setName((subVolume0PixVal == EXTRACELLULAR_PIXVAL ? EXTRACELLULAR_NAME : CYTOSOL_NAME));
    int subVolume1PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(1)).getPixelValue();
    geometry.getGeometrySpec().getSubVolume(1).setName((subVolume1PixVal == CYTOSOL_PIXVAL ? CYTOSOL_NAME : EXTRACELLULAR_NAME));
    geometry.getGeometrySurfaceDescription().updateAll();
    BioModel bioModel = new BioModel(null);
    bioModel.setName("unnamed");
    Model model = new Model("model");
    bioModel.setModel(model);
    model.addFeature(EXTRACELLULAR_NAME);
    Feature extracellular = (Feature) model.getStructure(EXTRACELLULAR_NAME);
    model.addFeature(CYTOSOL_NAME);
    Feature cytosol = (Feature) model.getStructure(CYTOSOL_NAME);
    // Membrane mem = model.addMembrane(EXTRACELLULAR_CYTOSOL_MEM_NAME);
    // model.getStructureTopology().setInsideFeature(mem, cytosol);
    // model.getStructureTopology().setOutsideFeature(mem, extracellular);
    String roiDataName = FRAPStudy.ROI_EXTDATA_NAME;
    final int SPECIES_COUNT = 4;
    final int FREE_SPECIES_INDEX = 0;
    final int BS_SPECIES_INDEX = 1;
    final int COMPLEX_SPECIES_INDEX = 2;
    final int IMMOBILE_SPECIES_INDEX = 3;
    Expression[] diffusionConstants = null;
    Species[] species = null;
    SpeciesContext[] speciesContexts = null;
    Expression[] initialConditions = null;
    diffusionConstants = new Expression[SPECIES_COUNT];
    species = new Species[SPECIES_COUNT];
    speciesContexts = new SpeciesContext[SPECIES_COUNT];
    initialConditions = new Expression[SPECIES_COUNT];
    // total initial condition
    FieldFunctionArguments postBleach_first = new FieldFunctionArguments(roiDataName, "postbleach_first", new Expression(0), VariableType.VOLUME);
    FieldFunctionArguments prebleach_avg = new FieldFunctionArguments(roiDataName, "prebleach_avg", new Expression(0), VariableType.VOLUME);
    Expression expPostBleach_first = new Expression(postBleach_first.infix());
    Expression expPreBleach_avg = new Expression(prebleach_avg.infix());
    Expression totalIniCondition = Expression.div(expPostBleach_first, expPreBleach_avg);
    // Free Species
    diffusionConstants[FREE_SPECIES_INDEX] = new Expression(df);
    species[FREE_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_MOBILE, "Mobile bleachable species");
    speciesContexts[FREE_SPECIES_INDEX] = new SpeciesContext(null, species[FREE_SPECIES_INDEX].getCommonName(), species[FREE_SPECIES_INDEX], cytosol);
    initialConditions[FREE_SPECIES_INDEX] = Expression.mult(new Expression(ff), totalIniCondition);
    // Immobile Species (No diffusion)
    // Set very small diffusion rate on immobile to force evaluation as state variable (instead of FieldData function)
    // If left as a function errors occur because functions involving FieldData require a database connection
    final String IMMOBILE_DIFFUSION_KLUDGE = "1e-14";
    diffusionConstants[IMMOBILE_SPECIES_INDEX] = new Expression(IMMOBILE_DIFFUSION_KLUDGE);
    species[IMMOBILE_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_IMMOBILE, "Immobile bleachable species");
    speciesContexts[IMMOBILE_SPECIES_INDEX] = new SpeciesContext(null, species[IMMOBILE_SPECIES_INDEX].getCommonName(), species[IMMOBILE_SPECIES_INDEX], cytosol);
    initialConditions[IMMOBILE_SPECIES_INDEX] = Expression.mult(new Expression(fimm), totalIniCondition);
    // BS Species
    diffusionConstants[BS_SPECIES_INDEX] = new Expression(IMMOBILE_DIFFUSION_KLUDGE);
    species[BS_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_BINDING_SITE, "Binding Site species");
    speciesContexts[BS_SPECIES_INDEX] = new SpeciesContext(null, species[BS_SPECIES_INDEX].getCommonName(), species[BS_SPECIES_INDEX], cytosol);
    initialConditions[BS_SPECIES_INDEX] = Expression.mult(new Expression(bs), totalIniCondition);
    // Complex species
    diffusionConstants[COMPLEX_SPECIES_INDEX] = new Expression(dc);
    species[COMPLEX_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_SLOW_MOBILE, "Slower mobile bleachable species");
    speciesContexts[COMPLEX_SPECIES_INDEX] = new SpeciesContext(null, species[COMPLEX_SPECIES_INDEX].getCommonName(), species[COMPLEX_SPECIES_INDEX], cytosol);
    initialConditions[COMPLEX_SPECIES_INDEX] = Expression.mult(new Expression(fc), totalIniCondition);
    // add reactions to species if there is bleachWhileMonitoring rate.
    for (int i = 0; i < initialConditions.length; i++) {
        model.addSpecies(species[i]);
        model.addSpeciesContext(speciesContexts[i]);
        // reaction with BMW rate, which should not be applied to binding site
        if (!(species[i].getCommonName().equals(FRAPStudy.SPECIES_NAME_PREFIX_BINDING_SITE))) {
            SimpleReaction simpleReaction = new SimpleReaction(model, cytosol, speciesContexts[i].getName() + "_bleach", true);
            model.addReactionStep(simpleReaction);
            simpleReaction.addReactant(speciesContexts[i], 1);
            MassActionKinetics massActionKinetics = new MassActionKinetics(simpleReaction);
            simpleReaction.setKinetics(massActionKinetics);
            KineticsParameter kforward = massActionKinetics.getForwardRateParameter();
            simpleReaction.getKinetics().setParameterValue(kforward, new Expression(new Double(bwmRate)));
        }
    }
    // add the binding reaction: F + BS <-> C
    SimpleReaction simpleReaction2 = new SimpleReaction(model, cytosol, "reac_binding", true);
    model.addReactionStep(simpleReaction2);
    simpleReaction2.addReactant(speciesContexts[FREE_SPECIES_INDEX], 1);
    simpleReaction2.addReactant(speciesContexts[BS_SPECIES_INDEX], 1);
    simpleReaction2.addProduct(speciesContexts[COMPLEX_SPECIES_INDEX], 1);
    MassActionKinetics massActionKinetics = new MassActionKinetics(simpleReaction2);
    simpleReaction2.setKinetics(massActionKinetics);
    KineticsParameter kforward = massActionKinetics.getForwardRateParameter();
    KineticsParameter kreverse = massActionKinetics.getReverseRateParameter();
    simpleReaction2.getKinetics().setParameterValue(kforward, new Expression(new Double(onRate)));
    simpleReaction2.getKinetics().setParameterValue(kreverse, new Expression(new Double(offRate)));
    // create simulation context
    SimulationContext simContext = new SimulationContext(bioModel.getModel(), geometry);
    bioModel.addSimulationContext(simContext);
    FeatureMapping cytosolFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(cytosol);
    FeatureMapping extracellularFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(extracellular);
    // Membrane plasmaMembrane = model.getStructureTopology().getMembrane(cytosol, extracellular);
    // MembraneMapping plasmaMembraneMapping = (MembraneMapping)simContext.getGeometryContext().getStructureMapping(plasmaMembrane);
    SubVolume cytSubVolume = geometry.getGeometrySpec().getSubVolume(CYTOSOL_NAME);
    SubVolume exSubVolume = geometry.getGeometrySpec().getSubVolume(EXTRACELLULAR_NAME);
    SurfaceClass pmSurfaceClass = geometry.getGeometrySurfaceDescription().getSurfaceClass(exSubVolume, cytSubVolume);
    cytosolFeatureMapping.setGeometryClass(cytSubVolume);
    extracellularFeatureMapping.setGeometryClass(exSubVolume);
    // plasmaMembraneMapping.setGeometryClass(pmSurfaceClass);
    cytosolFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    extracellularFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    for (int i = 0; i < speciesContexts.length; i++) {
        SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i]);
        scs.getInitialConditionParameter().setExpression(initialConditions[i]);
        scs.getDiffusionParameter().setExpression(diffusionConstants[i]);
    }
    MathMapping mathMapping = simContext.createNewMathMapping();
    MathDescription mathDesc = mathMapping.getMathDescription();
    // Add total fluorescence as function of mobile(optional: and slower mobile) and immobile fractions
    mathDesc.addVariable(new Function(FRAPStudy.SPECIES_NAME_PREFIX_COMBINED, new Expression(species[FREE_SPECIES_INDEX].getCommonName() + "+" + species[COMPLEX_SPECIES_INDEX].getCommonName() + "+" + species[IMMOBILE_SPECIES_INDEX].getCommonName()), null));
    simContext.setMathDescription(mathDesc);
    SimulationVersion simVersion = new SimulationVersion(simKey, "sim1", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
    Simulation newSimulation = new Simulation(simVersion, mathDesc);
    simContext.addSimulation(newSimulation);
    newSimulation.getSolverTaskDescription().setTimeBounds(timeBounds);
    newSimulation.getMeshSpecification().setSamplingSize(cellROI_2D.getISize());
    // newSimulation.getSolverTaskDescription().setTimeStep(timeStep); // Sundials doesn't need time step
    newSimulation.getSolverTaskDescription().setSolverDescription(SolverDescription.SundialsPDE);
    // use exp time step as output time spec
    newSimulation.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec(timeStepVal));
    return bioModel;
}

Example 18 with FeatureMapping

use of cbit.vcell.mapping.FeatureMapping in project vcell by virtualcell.

the class ITextWriter method writeStructureMapping.

// boundary types ignored.
protected void writeStructureMapping(Section simContextSection, SimulationContext sc) throws DocumentException {
    GeometryContext gc = sc.getGeometryContext();
    if (gc == null) {
        return;
    }
    Section structMappSection = null;
    /*try {
			ByteArrayOutputStream bos = generateStructureMappingImage(sc);
			com.lowagie.text.Image structMapImage = com.lowagie.text.Image.getInstance(Toolkit.getDefaultToolkit().createImage(bos.toByteArray()), null);
			structMappSection = simContextSection.addSection("Structure Mapping For " + sc.getName(), simContextSection.numberDepth() + 1);
			structMappSection.add(structMapImage);
		} catch (Exception e) {
			System.err.println("Unable to add structure mapping image to report.");
			e.printStackTrace();
		}*/
    StructureMapping[] structMap = gc.getStructureMappings();
    Table structMapTable = null;
    ModelUnitSystem modelUnitSystem = sc.getModel().getUnitSystem();
    for (int i = 0; i < structMap.length; i++) {
        if (!(structMap[i] instanceof FeatureMapping)) {
            continue;
        }
        if (i == 0) {
            structMapTable = getTable(5, 100, 1, 3, 3);
            structMapTable.addCell(createCell("Structure Mapping", getBold(DEF_HEADER_FONT_SIZE), 5, 1, Element.ALIGN_CENTER, true));
            structMapTable.addCell(createHeaderCell("Structure", getBold(), 1));
            structMapTable.addCell(createHeaderCell("Subdomain", getBold(), 1));
            structMapTable.addCell(createHeaderCell("Resolved (T/F)", getBold(), 1));
            structMapTable.addCell(createHeaderCell("Surf/Vol", getBold(), 1));
            structMapTable.addCell(createHeaderCell("VolFract", getBold(), 1));
            structMapTable.endHeaders();
        }
        String structName = structMap[i].getStructure().getName();
        SubVolume subVol = (SubVolume) ((FeatureMapping) structMap[i]).getGeometryClass();
        String subDomain = "";
        if (subVol != null) {
            subDomain = subVol.getName();
        }
        // ((FeatureMapping)structMap[i]).getResolved();
        boolean isResolved = false;
        String surfVolStr = "", volFractStr = "";
        MembraneMapping mm = (MembraneMapping) gc.getStructureMapping(sc.getModel().getStructureTopology().getParentStructure(structMap[i].getStructure()));
        if (mm != null) {
            StructureMapping.StructureMappingParameter smp = mm.getSurfaceToVolumeParameter();
            if (smp != null) {
                Expression tempExp = smp.getExpression();
                VCUnitDefinition tempUnit = smp.getUnitDefinition();
                if (tempExp != null) {
                    surfVolStr = tempExp.infix();
                    if (tempUnit != null && !modelUnitSystem.getInstance_DIMENSIONLESS().compareEqual(tempUnit)) {
                        // no need to add '1' for dimensionless unit
                        surfVolStr += " " + tempUnit.getSymbolUnicode();
                    }
                }
            }
            smp = mm.getVolumeFractionParameter();
            if (smp != null) {
                Expression tempExp = smp.getExpression();
                VCUnitDefinition tempUnit = smp.getUnitDefinition();
                if (tempExp != null) {
                    volFractStr = tempExp.infix();
                    if (tempUnit != null && !modelUnitSystem.getInstance_DIMENSIONLESS().compareEqual(tempUnit)) {
                        volFractStr += " " + tempUnit.getSymbolUnicode();
                    }
                }
            }
        }
        structMapTable.addCell(createCell(structName, getFont()));
        structMapTable.addCell(createCell(subDomain, getFont()));
        structMapTable.addCell(createCell((isResolved ? " T " : " F "), getFont()));
        structMapTable.addCell(createCell(surfVolStr, getFont()));
        structMapTable.addCell(createCell(volFractStr, getFont()));
    }
    if (structMapTable != null) {
        if (structMappSection == null) {
            structMappSection = simContextSection.addSection("Structure Mapping For " + sc.getName(), simContextSection.numberDepth() + 1);
        }
        structMappSection.add(structMapTable);
    }
}

Example 19 with FeatureMapping

use of cbit.vcell.mapping.FeatureMapping in project vcell by virtualcell.

the class StructureSizeSolver method updateUnitStructureSizes.

/**
 * 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 updateUnitStructureSizes(SimulationContext simContext, GeometryClass geometryClass) {
    if (simContext.getGeometryContext().getGeometry().getDimension() == 0) {
        return;
    }
    StructureMapping[] myStructMappings = simContext.getGeometryContext().getStructureMappings(geometryClass);
    if (myStructMappings != null && myStructMappings.length == 1) {
        // if the unitSizeParameter is dimensionless, then features are mapped to SubVolumes or Membranes are mapped to surfaces (should sum to 1)
        boolean bDimensionless = myStructMappings[0].getUnitSizeParameter().getUnitDefinition().isEquivalent(simContext.getModel().getUnitSystem().getInstance_DIMENSIONLESS());
        if (bDimensionless) {
            try {
                myStructMappings[0].getUnitSizeParameter().setExpression(new Expression(1.0));
                return;
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new RuntimeException(e.getMessage());
            }
        }
    }
    if (myStructMappings != null && myStructMappings.length == 0) {
        // nothing to solve, there are no mappings for this geometryClass
        return;
    }
    StructureMapping[] structMappings = simContext.getGeometryContext().getStructureMappings();
    try {
        ConstraintContainerImpl ccImpl = new ConstraintContainerImpl();
        Structure struct = null;
        Expression totalVolExpr = new Expression(0.0);
        StructureTopology structureTopology = simContext.getModel().getStructureTopology();
        for (int i = 0; i < structMappings.length; i++) {
            if (structMappings[i].getGeometryClass() != geometryClass) {
                continue;
            }
            // new model with unit sizes already
            if (structMappings[i].getUnitSizeParameter() != null && structMappings[i].getUnitSizeParameter().getExpression() != null) {
                return;
            }
            if (struct == null) {
                struct = structMappings[i].getStructure();
            }
            if (structMappings[i] instanceof MembraneMapping) {
                MembraneMapping membraneMapping = (MembraneMapping) structMappings[i];
                Membrane membrane = membraneMapping.getMembrane();
                String membraneSizeName = TokenMangler.mangleToSName(membrane.getName() + "_size");
                ccImpl.addSimpleBound(new SimpleBounds(membraneSizeName, new RealInterval(0, 100000), AbstractConstraint.PHYSICAL_LIMIT, "definition"));
                Feature insideFeature = structureTopology.getInsideFeature(membrane);
                String volFractName = TokenMangler.mangleToSName(insideFeature.getName() + "_volFract");
                String svRatioName = TokenMangler.mangleToSName(insideFeature.getName() + "_svRatio");
                StructureMapping.StructureMappingParameter volFractParameter = membraneMapping.getVolumeFractionParameter();
                double volFractValue = volFractParameter.getExpression().evaluateConstant();
                ccImpl.addSimpleBound(new SimpleBounds(volFractName, new RealInterval(volFractValue, volFractValue), AbstractConstraint.MODELING_ASSUMPTION, "from model"));
                StructureMapping.StructureMappingParameter surfToVolParameter = membraneMapping.getSurfaceToVolumeParameter();
                double svRatioValue = surfToVolParameter.getExpression().evaluateConstant();
                ccImpl.addSimpleBound(new SimpleBounds(svRatioName, new RealInterval(svRatioValue, svRatioValue), AbstractConstraint.MODELING_ASSUMPTION, "from model"));
                // membrane mapped to volume
                if (geometryClass instanceof SubVolume) {
                    // 
                    // EC eclosing cyt, which contains er and golgi
                    // "(cyt_size+ er_size + golgi_size) * cyt_svRatio - PM_size == 0"
                    // 
                    Expression sumOfInsideVolumeExp = new Expression(0.0);
                    for (int j = 0; j < structMappings.length; j++) {
                        if (structMappings[j] instanceof FeatureMapping && structureTopology.enclosedBy(structMappings[j].getStructure(), insideFeature)) {
                            Feature childFeatureOfInside = ((FeatureMapping) structMappings[j]).getFeature();
                            if (simContext.getGeometryContext().getStructureMapping(childFeatureOfInside).getGeometryClass() == geometryClass) {
                                sumOfInsideVolumeExp = Expression.add(sumOfInsideVolumeExp, new Expression(TokenMangler.mangleToSName(childFeatureOfInside.getName() + "_size")));
                            }
                        }
                    }
                    Expression tempExpr = Expression.mult(sumOfInsideVolumeExp, new Expression(svRatioName));
                    tempExpr = Expression.add(tempExpr, new Expression("-" + membraneSizeName));
                    ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(tempExpr.infix() + "==0"), AbstractConstraint.MODELING_ASSUMPTION, "svRatio definition"));
                    // 
                    // EC eclosing cyt, which contains er and golgi
                    // (EC_size + cyt_size + er_size + golgi_size) * cyt_vfRatio - (cyt_size + er_size + golgi_size) == 0
                    // 
                    Feature outsideFeature = structureTopology.getOutsideFeature(membrane);
                    Expression sumOfParentVolumeExp = new Expression(0.0);
                    for (int j = 0; j < structMappings.length; j++) {
                        if (structMappings[j] instanceof FeatureMapping && structureTopology.enclosedBy(structMappings[j].getStructure(), outsideFeature)) {
                            Feature childFeatureOfParent = ((FeatureMapping) structMappings[j]).getFeature();
                            if (simContext.getGeometryContext().getStructureMapping(childFeatureOfParent).getGeometryClass() == geometryClass) {
                                sumOfParentVolumeExp = Expression.add(sumOfParentVolumeExp, new Expression(TokenMangler.mangleToSName(childFeatureOfParent.getName() + "_size")));
                            }
                        }
                    }
                    Expression exp = Expression.mult(sumOfParentVolumeExp, new Expression(volFractName));
                    exp = Expression.add(exp, Expression.negate(sumOfInsideVolumeExp));
                    ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(exp.infix() + "==0.0"), AbstractConstraint.MODELING_ASSUMPTION, "volFract definition"));
                }
            } else if (structMappings[i] instanceof FeatureMapping) {
                FeatureMapping featureMapping = (FeatureMapping) structMappings[i];
                String featureSizeName = TokenMangler.mangleToSName(featureMapping.getFeature().getName() + "_size");
                totalVolExpr = Expression.add(totalVolExpr, new Expression(featureSizeName));
                ccImpl.addSimpleBound(new SimpleBounds(featureSizeName, new RealInterval(0, 1), AbstractConstraint.PHYSICAL_LIMIT, "definition"));
            }
        }
        if (geometryClass instanceof SubVolume) {
            ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(totalVolExpr.infix() + "==1.0"), AbstractConstraint.MODELING_ASSUMPTION, "total volume"));
        }
        // ccImpl.show();
        ConstraintSolver constraintSolver = new ConstraintSolver(ccImpl);
        constraintSolver.resetIntervals();
        int numTimesNarrowed = 0;
        RealInterval[] lastSolution = null;
        boolean bChanged = true;
        while (constraintSolver.narrow() && bChanged && numTimesNarrowed < 125) {
            numTimesNarrowed++;
            bChanged = false;
            RealInterval[] thisSolution = constraintSolver.getIntervals();
            if (lastSolution != null) {
                for (int i = 0; i < thisSolution.length; i++) {
                    if (!thisSolution[i].equals(lastSolution[i])) {
                        bChanged = true;
                    }
                }
            } else {
                bChanged = true;
            }
            lastSolution = thisSolution;
        }
        System.out.println("num of times narrowed = " + numTimesNarrowed);
        if (numTimesNarrowed > 0) {
            String[] symbols = constraintSolver.getSymbols();
            net.sourceforge.interval.ia_math.RealInterval[] solution = constraintSolver.getIntervals();
            double totalArea = 0;
            double totalVolume = 0;
            for (int i = 0; i < symbols.length; i++) {
                System.out.println("solution[" + i + "] \"" + symbols[i] + "\" = " + solution[i]);
                for (int j = 0; j < structMappings.length; j++) {
                    if (symbols[i].equals(TokenMangler.mangleToSName(structMappings[j].getStructure().getName() + "_size"))) {
                        if (!Double.isInfinite(solution[i].lo()) && !Double.isInfinite(solution[i].hi())) {
                            double value = (solution[i].lo() + solution[i].hi()) / 2;
                            Expression exp = new Expression(value);
                            if (structMappings[j] instanceof FeatureMapping) {
                                FeatureMapping fm = (FeatureMapping) structMappings[j];
                                totalVolume += value;
                                if (geometryClass instanceof SubVolume) {
                                    fm.getVolumePerUnitVolumeParameter().setExpression(exp);
                                } else if (geometryClass instanceof SurfaceClass) {
                                    fm.getVolumePerUnitAreaParameter().setExpression(exp);
                                }
                            } else if (structMappings[j] instanceof MembraneMapping) {
                                MembraneMapping mm = (MembraneMapping) structMappings[j];
                                totalArea += value;
                                if (geometryClass instanceof SubVolume) {
                                    mm.getAreaPerUnitVolumeParameter().setExpression(exp);
                                } else if (geometryClass instanceof SurfaceClass) {
                                    mm.getAreaPerUnitAreaParameter().setExpression(exp);
                                }
                            }
                        }
                    }
                }
            }
            // 
            // normalize all so that total volume is 1.0 for subVolumes or
            // total area is 1.0 for surfaceClasses
            // 
            double scaleFactor = 1;
            if (geometryClass instanceof SubVolume) {
                scaleFactor = totalVolume;
            } else if (geometryClass instanceof SurfaceClass) {
                scaleFactor = totalArea;
            } else {
                throw new RuntimeException("unexpected GeometryClass");
            }
            for (int j = 0; j < structMappings.length; j++) {
                if (structMappings[j].getGeometryClass() == geometryClass) {
                    if (structMappings[j] instanceof FeatureMapping) {
                        FeatureMapping fm = (FeatureMapping) structMappings[j];
                        if (geometryClass instanceof SubVolume) {
                            fm.getVolumePerUnitVolumeParameter().setExpression(new Expression(fm.getVolumePerUnitVolumeParameter().getExpression().evaluateConstant() / scaleFactor));
                        } else if (geometryClass instanceof SurfaceClass) {
                            fm.getVolumePerUnitAreaParameter().setExpression(new Expression(fm.getVolumePerUnitAreaParameter().getExpression().evaluateConstant() / scaleFactor));
                        }
                    } else if (structMappings[j] instanceof MembraneMapping) {
                        MembraneMapping mm = (MembraneMapping) structMappings[j];
                        if (geometryClass instanceof SubVolume) {
                            mm.getAreaPerUnitVolumeParameter().setExpression(new Expression(mm.getAreaPerUnitVolumeParameter().getExpression().evaluateConstant() / scaleFactor));
                        } else if (geometryClass instanceof SurfaceClass) {
                            mm.getAreaPerUnitAreaParameter().setExpression(new Expression(mm.getAreaPerUnitAreaParameter().getExpression().evaluateConstant() / scaleFactor));
                        }
                    }
                }
            }
        } else {
            throw new RuntimeException("cannot solve for size");
        }
    } catch (ExpressionException e) {
        e.printStackTrace(System.out);
        throw new RuntimeException(e.getMessage());
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new RuntimeException(e.getMessage());
    }
}

Example 20 with FeatureMapping

use of cbit.vcell.mapping.FeatureMapping in project vcell by virtualcell.

the class SimulationContextDbDriver method assignStructureMappingsSQL.

/**
 * This method was created in VisualAge.
 * @param simContext cbit.vcell.mapping.SimulationContext
 */
private void assignStructureMappingsSQL(QueryHashtable dbc, Connection con, KeyValue simContextKey, SimulationContext simContext) throws SQLException, DataAccessException {
    String sql;
    sql = " SELECT " + "*" + " FROM " + structureMappingTable.getTableName() + " WHERE " + structureMappingTable.simContextRef + " = " + simContextKey;
    Statement stmt = con.createStatement();
    try {
        // log.print(sql);
        ResultSet rset = stmt.executeQuery(sql);
        while (rset.next()) {
            BigDecimal subvolumeRefBigDecimal = rset.getBigDecimal(structureMappingTable.subVolumeRef.toString());
            KeyValue subVolumeRef = (subvolumeRefBigDecimal == null ? null : new KeyValue(subvolumeRefBigDecimal));
            BigDecimal surfaceClassRefBigDecimal = rset.getBigDecimal(structureMappingTable.surfaceClassRef.toString());
            KeyValue surfaceClassRef = (surfaceClassRefBigDecimal == null ? null : new KeyValue(surfaceClassRefBigDecimal));
            KeyValue structureRef = new KeyValue(rset.getBigDecimal(structureMappingTable.structRef.toString()));
            // 
            // lookup structure and subVolume from SimulationContext by their keys
            // 
            // DBCache will not always give same instance consistently (usually this is
            // fixed up later in the ReferenceResolver at the Client).
            // 
            Structure theStructure = null;
            Structure[] structureArray = simContext.getModel().getStructures();
            for (int i = 0; i < structureArray.length; i++) {
                Structure structure = structureArray[i];
                if (structure.getKey().compareEqual(structureRef)) {
                    theStructure = structure;
                    break;
                }
            }
            if (theStructure == null) {
                throw new DataAccessException("Can't match structure and subvolume");
            }
            GeometryClass theGeometryClass = null;
            KeyValue geometryClassKey = (subVolumeRef == null ? surfaceClassRef : subVolumeRef);
            if (geometryClassKey != null) {
                GeometryClass[] geometryClasses = simContext.getGeometry().getGeometryClasses();
                for (int i = 0; i < geometryClasses.length; i++) {
                    if (geometryClasses[i].getKey().compareEqual(geometryClassKey)) {
                        theGeometryClass = geometryClasses[i];
                        break;
                    }
                }
                if (theGeometryClass == null) {
                    throw new DataAccessException("Can't find Geometryclass");
                }
            }
            Expression sizeExpression = null;
            String sizeExpressionS = rset.getString(StructureMappingTable.table.sizeExp.getUnqualifiedColName());
            if (!rset.wasNull() && sizeExpressionS != null && sizeExpressionS.length() > 0) {
                try {
                    sizeExpressionS = TokenMangler.getSQLRestoredString(sizeExpressionS);
                    sizeExpression = new Expression(sizeExpressionS);
                } catch (ExpressionException e) {
                    e.printStackTrace();
                    throw new DataAccessException("SimulationContextDbDriver.assignStructureMappingSQL : Couldn't parse non-null size expression for Structure " + theStructure.getName());
                }
            }
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(theStructure);
            try {
                sm.getSizeParameter().setExpression(sizeExpression);
            } catch (Exception e1) {
                throw new DataAccessException("SimulationContextDbDriver.assignStructureMappingSQL : Couldn't set size expression '" + sizeExpressionS + "'for Structure " + theStructure.getName());
            }
            try {
                sm.setGeometryClass(theGeometryClass);
            } catch (PropertyVetoException e) {
                lg.error(e.getMessage(), e);
                throw new DataAccessException(e.getMessage());
            }
            if (sm instanceof FeatureMapping) {
                FeatureMapping fm = (FeatureMapping) sm;
                String boundaryTypeXmString = rset.getString(structureMappingTable.boundaryTypeXm.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeXm(new BoundaryConditionType(boundaryTypeXmString));
                }
                String boundaryTypeXpString = rset.getString(structureMappingTable.boundaryTypeXp.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeXp(new BoundaryConditionType(boundaryTypeXpString));
                }
                String boundaryTypeYmString = rset.getString(structureMappingTable.boundaryTypeYm.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeYm(new BoundaryConditionType(boundaryTypeYmString));
                }
                String boundaryTypeYpString = rset.getString(structureMappingTable.boundaryTypeYp.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeYp(new BoundaryConditionType(boundaryTypeYpString));
                }
                String boundaryTypeZmString = rset.getString(structureMappingTable.boundaryTypeZm.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeZm(new BoundaryConditionType(boundaryTypeZmString));
                }
                String boundaryTypeZpString = rset.getString(structureMappingTable.boundaryTypeZp.toString());
                if (!rset.wasNull()) {
                    fm.setBoundaryConditionTypeZp(new BoundaryConditionType(boundaryTypeZpString));
                }
                String volPerUnitArea = rset.getString(structureMappingTable.volPerUnitAreaExp.toString());
                if (!rset.wasNull()) {
                    try {
                        fm.getVolumePerUnitAreaParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(volPerUnitArea)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in surfaceToVol expression: " + e.getMessage());
                    }
                }
                String volPerUnitVol = rset.getString(structureMappingTable.volPerUnitVolExp.toString());
                if (!rset.wasNull()) {
                    try {
                        fm.getVolumePerUnitVolumeParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(volPerUnitVol)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in surfaceToVol expression: " + e.getMessage());
                    }
                }
            } else if (sm instanceof MembraneMapping) {
                MembraneMapping mm = (MembraneMapping) sm;
                String surfToVolString = rset.getString(structureMappingTable.surfToVolExp.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getSurfaceToVolumeParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(surfToVolString)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in surfaceToVol expression: " + e.getMessage());
                    }
                }
                String volFractString = rset.getString(structureMappingTable.volFractExp.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getVolumeFractionParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(volFractString)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in volFract expression: " + e.getMessage());
                    }
                }
                boolean bCalculateVoltage = rset.getBoolean(structureMappingTable.bCalculateVoltage.toString());
                if (!rset.wasNull()) {
                    mm.setCalculateVoltage(bCalculateVoltage);
                }
                java.math.BigDecimal specificCapacitance = rset.getBigDecimal(structureMappingTable.specificCap.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getSpecificCapacitanceParameter().setExpression(new Expression(specificCapacitance.doubleValue()));
                    } catch (ExpressionBindingException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("error setting membrane specific capacitance: " + e.getMessage());
                    }
                }
                String initialVoltageString = rset.getString(structureMappingTable.initialVoltage.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getInitialVoltageParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(initialVoltageString)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("database parse error in initial membrane voltage: " + e.getMessage());
                    }
                }
                String areaPerUnitArea = rset.getString(structureMappingTable.areaPerUnitAreaExp.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getAreaPerUnitAreaParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(areaPerUnitArea)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in surfaceToVol expression: " + e.getMessage());
                    }
                }
                String areaPerUnitVol = rset.getString(structureMappingTable.areaPerUnitVolExp.toString());
                if (!rset.wasNull()) {
                    try {
                        mm.getAreaPerUnitVolumeParameter().setExpression(new Expression(TokenMangler.getSQLRestoredString(areaPerUnitVol)));
                    } catch (ExpressionException e) {
                        e.printStackTrace(System.out);
                        throw new DataAccessException("parse error in surfaceToVol expression: " + e.getMessage());
                    }
                }
            } else {
                throw new DataAccessException("unknown structureMapping type");
            }
        // System.out.println("Structure Key = " + theStructure + " - " + "SubVolume Key " + theSubVolume.getKey());
        }
    } finally {
        stmt.close();
    }
}