Examples with MathDescription cbit.vcell.math.MathDescription used on opensource projects

Example 91 with MathDescription

use of cbit.vcell.math.MathDescription 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 92 with MathDescription

use of cbit.vcell.math.MathDescription in project vcell by virtualcell.

the class MathOverrides method updateFromMathDescription.

void updateFromMathDescription() {
    MathDescription mathDescription = getSimulation().getMathDescription();
    // 
    // get list of names of constants in this math
    // 
    Enumeration<Constant> enumeration = mathDescription.getConstants();
    java.util.HashSet<String> mathDescriptionHash = new java.util.HashSet<String>();
    while (enumeration.hasMoreElements()) {
        Constant constant = enumeration.nextElement();
        mathDescriptionHash.add(constant.getName());
    }
    // 
    // for any elements in this MathOverrides but not in the new MathDescription
    // 
    // 1) try to "repair" overridden constants for automatically generated constant names (via math generation)
    // which have changed due to changes to Math generation
    // 
    // 2) if not repaired, will be reported as an issue.
    // 
    HashMap<String, String> renamedMap = new HashMap<String, String>();
    boolean bNameRepaired = true;
    while (bNameRepaired) {
        bNameRepaired = false;
        Enumeration<String> mathOverrideNamesEnum = getOverridesHash().keys();
        while (mathOverrideNamesEnum.hasMoreElements()) {
            String name = mathOverrideNamesEnum.nextElement();
            if (!mathDescriptionHash.contains(name)) {
                // 
                if (name.endsWith(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old)) {
                    String name_repaired_to_uM = name.replace(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old, DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old_uM);
                    if (mathDescriptionHash.contains(name_repaired_to_uM)) {
                        Element element = overridesHash.remove(name);
                        element.name = name_repaired_to_uM;
                        overridesHash.put(name_repaired_to_uM, element);
                        renamedMap.put(name, name_repaired_to_uM);
                        removeConstant(name);
                        bNameRepaired = true;
                        break;
                    }
                }
                // 
                if (name.endsWith(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old)) {
                    String name_repaired_to_molecule_per_um2 = name.replace(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old, DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old_molecules_um_2);
                    if (mathDescriptionHash.contains(name_repaired_to_molecule_per_um2)) {
                        Element element = overridesHash.remove(name);
                        element.name = name_repaired_to_molecule_per_um2;
                        overridesHash.put(name_repaired_to_molecule_per_um2, element);
                        renamedMap.put(name, name_repaired_to_molecule_per_um2);
                        removeConstant(name);
                        bNameRepaired = true;
                        break;
                    }
                }
                // 
                if (name.endsWith(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old_molecules_per_um2)) {
                    String name_repaired_to_molecule_per_um2 = name.replace(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old_molecules_per_um2, DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONCENTRATION_old_molecules_um_2);
                    if (mathDescriptionHash.contains(name_repaired_to_molecule_per_um2)) {
                        Element element = overridesHash.remove(name);
                        element.name = name_repaired_to_molecule_per_um2;
                        overridesHash.put(name_repaired_to_molecule_per_um2, element);
                        renamedMap.put(name, name_repaired_to_molecule_per_um2);
                        removeConstant(name);
                        bNameRepaired = true;
                        break;
                    }
                }
            }
        }
    }
    // 
    for (String origName : renamedMap.keySet()) {
        for (Element element : overridesHash.values()) {
            if (element.actualValue != null && element.actualValue.hasSymbol(origName)) {
                try {
                    element.actualValue = element.actualValue.getSubstitutedExpression(new Expression(origName), new Expression(mathDescription.getVariable(renamedMap.get(origName)), null));
                } catch (ExpressionException e) {
                    // won't happen ... here because new Expression(origName) throws a parse exception ... but it must be able to parse.
                    e.printStackTrace();
                }
            }
        }
    }
    refreshDependencies();
}

Example 93 with MathDescription

use of cbit.vcell.math.MathDescription in project vcell by virtualcell.

the class MathOverrides method removeUnusedOverrides.

/**
 * explicit user action invokes this method (currently pressing a button), this is not automatic cleanup of old/obsolete models ... we want the user to address the issues.
 */
public void removeUnusedOverrides() {
    MathDescription mathDescription = getSimulation().getMathDescription();
    // 
    // get list of names of constants in this math
    // 
    Enumeration<Constant> enumeration = mathDescription.getConstants();
    java.util.HashSet<String> mathDescriptionHash = new java.util.HashSet<String>();
    while (enumeration.hasMoreElements()) {
        Constant constant = enumeration.nextElement();
        mathDescriptionHash.add(constant.getName());
    }
    // 
    // for any elements in this MathOverrides but not in the new MathDescription, add an "error" issue
    // 
    ArrayList<String> overridesToDelete = new ArrayList<String>();
    Enumeration<String> mathOverrideNamesEnum = getOverridesHash().keys();
    while (mathOverrideNamesEnum.hasMoreElements()) {
        String name = mathOverrideNamesEnum.nextElement();
        if (!mathDescriptionHash.contains(name)) {
            // constant 'name' no longer part of mathDescription
            overridesToDelete.add(name);
        } else if (getSimulation().getSimulationOwner() != null) {
            IssueContext issueContext = new IssueContext(ContextType.Simulation, getSimulation(), null);
            if (getSimulation().getSimulationOwner().gatherIssueForMathOverride(issueContext, getSimulation(), name) != null) {
                // constant 'name' is part of math, but simulation owner doesn't want you to override it.
                overridesToDelete.add(name);
            }
        }
    }
    for (String deletedName : overridesToDelete) {
        removeConstant(deletedName);
    }
}

Example 94 with MathDescription

use of cbit.vcell.math.MathDescription in project vcell by virtualcell.

the class MathOverrides method verifyExpression.

/**
 * Insert the method's description here.
 * Creation date: (9/22/2005 3:24:36 PM)
 */
private void verifyExpression(Constant value, boolean checkScan) throws ExpressionBindingException {
    Expression exp = value.getExpression();
    String[] symbols = exp.getSymbols();
    MathDescription mathDescription = getSimulation().getMathDescription();
    exp.bindExpression(mathDescription);
    if (symbols != null) {
        for (int i = 0; i < symbols.length; i++) {
            // 
            // expression must be a function of another Simulation parameter
            // 
            Variable variable = mathDescription.getVariable(symbols[i]);
            if (!(variable != null && variable instanceof Constant)) {
                throw new ExpressionBindingException("identifier " + symbols[i] + " is not a constant. " + "Parameter overrides should only refer to constants.");
            }
            if (checkScan && isScan(symbols[i])) {
                throw new ExpressionBindingException("Parameter overrides cannot depend on another scanned parameter " + symbols[i]);
            }
            // 
            if (symbols[i].equals(value.getName())) {
                throw new ExpressionBindingException("Parameter overrides can not be recursive definition, can't use identifier " + value.getName() + " in expression for " + value.getName());
            }
        }
    }
}

Example 95 with MathDescription

use of cbit.vcell.math.MathDescription in project vcell by virtualcell.

the class MathOverrides method gatherIssues.

/**
 * This method was created in VisualAge.
 * @return boolean
 */
public void gatherIssues(IssueContext issueContext, List<Issue> issueList) {
    MathDescription mathDescription = getSimulation().getMathDescription();
    // 
    // get list of names of constants in this math
    // 
    Enumeration<Constant> enumeration = mathDescription.getConstants();
    java.util.HashSet<String> mathDescriptionHash = new java.util.HashSet<String>();
    while (enumeration.hasMoreElements()) {
        Constant constant = enumeration.nextElement();
        mathDescriptionHash.add(constant.getName());
    }
    // 
    // for any elements in this MathOverrides but not in the new MathDescription, add an "error" issue
    // 
    Enumeration<String> mathOverrideNamesEnum = getOverridesHash().keys();
    while (mathOverrideNamesEnum.hasMoreElements()) {
        String name = mathOverrideNamesEnum.nextElement();
        if (!mathDescriptionHash.contains(name)) {
            Issue issue = new Issue(getSimulation(), issueContext, IssueCategory.Simulation_Override_NotFound, VCellErrorMessages.getErrorMessage(VCellErrorMessages.SIMULATION_OVERRIDE_NOTFOUND, name, getSimulation().getName()), Issue.SEVERITY_ERROR);
            issueList.add(issue);
        }
        Variable var = mathDescription.getVariable(name);
        if (getSimulation().getSimulationOwner() != null) {
            Issue issue = getSimulation().getSimulationOwner().gatherIssueForMathOverride(issueContext, getSimulation(), name);
            if (issue != null) {
                issueList.add(issue);
            }
        }
    }
}