Examples with UniformOutputTimeSpec cbit.vcell.solver.UniformOutputTimeSpec used on opensource projects

Example 26 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec in project vcell by virtualcell.

the class SimulationListTableModel method setValueAt.

/**
 * Insert the method's description here.
 * Creation date: (7/12/2004 2:01:23 PM)
 * @param aValue java.lang.Object
 * @param rowIndex int
 * @param columnIndex int
 */
public void setValueAt(Object aValue, int rowIndex, int columnIndex) {
    Simulation simulation = getValueAt(rowIndex);
    try {
        switch(columnIndex) {
            case COLUMN_NAME:
                if (aValue instanceof String) {
                    String newName = (String) aValue;
                    if (!simulation.getName().equals(newName)) {
                        simulation.setName(newName);
                    }
                }
                break;
            case COLUMN_ENDTIME:
                if (aValue instanceof Double) {
                    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
                    double newEndTime = (Double) aValue;
                    if (newEndTime != solverTaskDescription.getTimeBounds().getEndingTime()) {
                        ClientTaskManager.changeEndTime(ownerTable, solverTaskDescription, newEndTime);
                        simulation.setIsDirty(true);
                    }
                }
                break;
            case COLUMN_OUTPUT:
                if (aValue instanceof String) {
                    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
                    OutputTimeSpec ots = solverTaskDescription.getOutputTimeSpec();
                    OutputTimeSpec newOts = null;
                    if (ots instanceof DefaultOutputTimeSpec) {
                        int newValue = Integer.parseInt((String) aValue);
                        newOts = new DefaultOutputTimeSpec(newValue, ((DefaultOutputTimeSpec) ots).getKeepAtMost());
                    } else if (ots instanceof UniformOutputTimeSpec) {
                        try {
                            boolean bValid = true;
                            double outputTime = Double.parseDouble((String) aValue);
                            if (solverTaskDescription.getOutputTimeSpec().isUniform()) {
                                if (solverTaskDescription.getSolverDescription().hasVariableTimestep()) {
                                    newOts = new UniformOutputTimeSpec(outputTime);
                                } else {
                                    double timeStep = solverTaskDescription.getTimeStep().getDefaultTimeStep();
                                    double suggestedInterval = outputTime;
                                    if (outputTime < timeStep) {
                                        suggestedInterval = timeStep;
                                        bValid = false;
                                    } else {
                                        if (!BeanUtils.isIntegerMultiple(outputTime, timeStep)) {
                                            double n = outputTime / timeStep;
                                            int intn = (int) Math.round(n);
                                            suggestedInterval = (intn * timeStep);
                                            if (suggestedInterval != outputTime) {
                                                bValid = false;
                                            }
                                        }
                                    }
                                    if (bValid) {
                                        newOts = new UniformOutputTimeSpec(outputTime);
                                    } else {
                                        String ret = PopupGenerator.showWarningDialog(ownerTable, "Output Interval", "Output Interval must " + "be integer multiple of time step.\n\nChange Output Interval to " + suggestedInterval + "?", new String[] { UserMessage.OPTION_YES, UserMessage.OPTION_NO }, UserMessage.OPTION_YES);
                                        if (ret.equals(UserMessage.OPTION_YES)) {
                                            newOts = new UniformOutputTimeSpec(suggestedInterval);
                                        }
                                    }
                                }
                            }
                        } catch (NumberFormatException ex) {
                            DialogUtils.showErrorDialog(ownerTable, "Wrong number format " + ex.getMessage().toLowerCase());
                        }
                    } else if (ots instanceof ExplicitOutputTimeSpec) {
                        newOts = ExplicitOutputTimeSpec.fromString((String) aValue);
                    }
                    if (newOts != null && !newOts.compareEqual(ots)) {
                        solverTaskDescription.setOutputTimeSpec(newOts);
                        simulation.setIsDirty(true);
                    }
                }
                break;
        }
    } catch (Exception ex) {
        ex.printStackTrace(System.out);
        PopupGenerator.showErrorDialog(ownerTable, ex.getMessage());
    }
}

Example 27 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec in project vcell by virtualcell.

the class VCellSBMLSolver method solveVCell.

public File solveVCell(String filePrefix, File outDir, String sbmlFileName, SimSpec testSpec) throws IOException, SolverException, SbmlException {
    try {
        cbit.util.xml.VCLogger logger = new LocalLogger();
        // 
        // Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
        // and import SBML  (sbml->bioModel)
        org.vcell.sbml.vcell.SBMLImporter sbmlImporter = new org.vcell.sbml.vcell.SBMLImporter(sbmlFileName, logger, false);
        BioModel bioModel = sbmlImporter.getBioModel();
        if (bRoundTrip) {
            // Round trip the bioModel (bioModel->sbml->bioModel).
            // export bioModel as sbml and save
            String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, 0, false, bioModel.getSimulationContext(0), null);
            // re-import bioModel from exported sbml
            XMLSource vcml_sbml_Src = new XMLSource(vcml_sbml);
            BioModel newBioModel = (BioModel) XmlHelper.importSBML(logger, vcml_sbml_Src, false);
            // have rest of code use the round-tripped biomodel
            bioModel = newBioModel;
        }
        // 
        // select only Application, generate math, and create a single Simulation.
        // 
        SimulationContext simContext = bioModel.getSimulationContext(0);
        MathMapping mathMapping = simContext.createNewMathMapping();
        MathDescription mathDesc = mathMapping.getMathDescription();
        simContext.setMathDescription(mathDesc);
        SimulationVersion simVersion = new SimulationVersion(new KeyValue("100"), "unnamed", null, null, null, null, null, null, null, null);
        Simulation sim = new Simulation(simVersion, mathDesc);
        sim.setName("unnamed");
        // if time factor from SBML is not 1 (i.e., it is not in secs but in minutes or hours), convert endTime to min/hr as : endTime*timeFactor
        // double endTime = testSpec.getEndTime()*timeFactor;
        double endTime = testSpec.getEndTime();
        sim.getSolverTaskDescription().setTimeBounds(new TimeBounds(0, endTime));
        TimeStep timeStep = new TimeStep();
        sim.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep.getMinimumTimeStep(), timeStep.getDefaultTimeStep(), endTime / 10000));
        sim.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((endTime - 0) / testSpec.getNumTimeSteps()));
        sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(testSpec.getAbsTolerance(), testSpec.getRelTolerance()));
        // sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
        // Generate .idaInput string
        File idaInputFile = new File(outDir, filePrefix + SimDataConstants.IDAINPUT_DATA_EXTENSION);
        PrintWriter idaPW = new java.io.PrintWriter(idaInputFile);
        SimulationJob simJob = new SimulationJob(sim, 0, null);
        SimulationTask simTask = new SimulationTask(simJob, 0);
        IDAFileWriter idaFileWriter = new IDAFileWriter(idaPW, simTask);
        idaFileWriter.write();
        idaPW.close();
        // use the idastandalone solver
        File idaOutputFile = new File(outDir, filePrefix + SimDataConstants.IDA_DATA_EXTENSION);
        // String sundialsSolverExecutable = "C:\\Developer\\Eclipse\\workspace\\VCell 4.8\\SundialsSolverStandalone_NoMessaging.exe";
        String executableName = null;
        try {
            executableName = SolverUtilities.getExes(SolverDescription.IDA)[0].getAbsolutePath();
        } catch (IOException e) {
            throw new RuntimeException("failed to get executable for solver " + SolverDescription.IDA.getDisplayLabel() + ": " + e.getMessage(), e);
        }
        Executable executable = new Executable(new String[] { executableName, idaInputFile.getAbsolutePath(), idaOutputFile.getAbsolutePath() });
        executable.start();
        /*			// Generate .cvodeInput string
			File cvodeFile = new File(outDir,filePrefix+SimDataConstants.CVODEINPUT_DATA_EXTENSION);
			PrintWriter cvodePW = new java.io.PrintWriter(cvodeFile);
			SimulationJob simJob = new SimulationJob(sim, 0, null);
		    CVodeFileWriter cvodeFileWriter = new CVodeFileWriter(cvodePW, simJob);
			cvodeFileWriter.write();
			cvodePW.close();

			// use the cvodeStandalone solver
			File cvodeOutputFile = new File(outDir,filePrefix+SimDataConstants.IDA_DATA_EXTENSION);
			String sundialsSolverExecutable = PropertyLoader.getRequiredProperty(PropertyLoader.sundialsSolverExecutableProperty);
			Executable executable = new Executable(new String[]{sundialsSolverExecutable, cvodeFile.getAbsolutePath(), cvodeOutputFile.getAbsolutePath()});
			executable.start();
*/
        // get the result
        ODESolverResultSet odeSolverResultSet = getODESolverResultSet(simJob, idaOutputFile.getPath());
        // remove CVOde input and output files ??
        idaInputFile.delete();
        idaOutputFile.delete();
        // 
        // print header
        // 
        File outputFile = new File(outDir, "results" + filePrefix + ".csv");
        java.io.PrintStream outputStream = new java.io.PrintStream(new java.io.BufferedOutputStream(new java.io.FileOutputStream(outputFile)));
        outputStream.print("time");
        for (int i = 0; i < testSpec.getVarsList().length; i++) {
            outputStream.print("," + testSpec.getVarsList()[i]);
        }
        outputStream.println();
        // 
        // extract data for time and species
        // 
        double[][] data = new double[testSpec.getVarsList().length + 1][];
        int column = odeSolverResultSet.findColumn("t");
        data[0] = odeSolverResultSet.extractColumn(column);
        int origDataLength = data[0].length;
        for (int i = 0; i < testSpec.getVarsList().length; i++) {
            column = odeSolverResultSet.findColumn(testSpec.getVarsList()[i]);
            if (column == -1) {
                Variable var = simJob.getSimulationSymbolTable().getVariable(testSpec.getVarsList()[i]);
                data[i + 1] = new double[data[0].length];
                if (var instanceof cbit.vcell.math.Constant) {
                    double value = ((cbit.vcell.math.Constant) var).getExpression().evaluateConstant();
                    for (int j = 0; j < data[i + 1].length; j++) {
                        data[i + 1][j] = value;
                    }
                } else {
                    throw new RuntimeException("Did not find " + testSpec.getVarsList()[i] + " in simulation");
                }
            } else {
                data[i + 1] = odeSolverResultSet.extractColumn(column);
            }
        }
        // 
        // for each time, print row
        // 
        int index = 0;
        double[] sampleTimes = new double[testSpec.getNumTimeSteps() + 1];
        for (int i = 0; i <= testSpec.getNumTimeSteps(); i++) {
            sampleTimes[i] = endTime * i / testSpec.getNumTimeSteps();
        }
        Model vcModel = bioModel.getModel();
        ReservedSymbol kMole = vcModel.getKMOLE();
        for (int i = 0; i < sampleTimes.length; i++) {
            // 
            while (true) {
                // 
                if (index == odeSolverResultSet.getRowCount() - 1) {
                    if (data[0][index] == sampleTimes[i]) {
                        break;
                    } else {
                        throw new RuntimeException("sampleTime does not match at last time point");
                    }
                }
                // 
                if (data[0][index + 1] > sampleTimes[i]) {
                    break;
                }
                // 
                // sampleTime must be later in our data list.
                // 
                index++;
            }
            // if data[0][index] == sampleTime no need to interpolate
            if (data[0][index] == sampleTimes[i]) {
                // if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
                // if (timeFactor != 1.0) {
                // outputStream.print(data[0][index]/timeFactor);
                // } else {
                outputStream.print(data[0][index]);
                // }
                for (int j = 0; j < testSpec.getVarsList().length; j++) {
                    // SBMLImporter.SBVCConcentrationUnits spConcUnits = speciesUnitsHash.get(testSpec.getVarsList()[j]);
                    // if (spConcUnits != null) {
                    // VCUnitDefinition sbunits = spConcUnits.getSBConcentrationUnits();
                    // VCUnitDefinition vcunits = spConcUnits.getVCConcentrationUnits();
                    // SBMLUnitParameter unitFactor = SBMLUtils.getConcUnitFactor("spConcParam", vcunits, sbunits, kMole);
                    // outputStream.print("," + data[j + 1][index] * unitFactor.getExpression().evaluateConstant()); 		//earlier, hack unitfactor = 0.000001
                    // earlier, hack unitfactor = 0.000001
                    outputStream.print("," + data[j + 1][index]);
                // }
                }
                // System.out.println("No interpolation needed!");
                outputStream.println();
            } else {
                // if data[0][index] < sampleTime, must interpolate
                double fraction = (sampleTimes[i] - data[0][index]) / (data[0][index + 1] - data[0][index]);
                // if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
                // if (timeFactor != 1.0) {
                // outputStream.print(sampleTimes[i]/timeFactor);
                // } else {
                outputStream.print(sampleTimes[i]);
                // }
                for (int j = 0; j < testSpec.getVarsList().length; j++) {
                    double interpolatedValue = 0.0;
                    double[] speciesVals = null;
                    double[] times = null;
                    // Currently using 2nd order interpolation
                    if (index == 0) {
                        // can only do 1st order interpolation
                        times = new double[] { data[0][index], data[0][index + 1] };
                        speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1] };
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    } else if (index >= 1 && index <= origDataLength - 3) {
                        double val_1 = Math.abs(sampleTimes[i] - data[0][index - 1]);
                        double val_2 = Math.abs(sampleTimes[i] - data[0][index + 2]);
                        if (val_1 < val_2) {
                            times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
                            speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
                        } else {
                            times = new double[] { data[0][index], data[0][index + 1], data[0][index + 2] };
                            speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1], data[j + 1][index + 2] };
                        }
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    } else {
                        times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
                        speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    }
                    // // Currently using 1st order interpolation
                    // times = new double[] { data[0][index], data[0][index+1] };
                    // speciesVals = new double[] { data[j+1][index], data[j+1][index+1] };
                    // interpolatedValue = taylorInterpolation(sampleTimes[i], times, speciesVals);
                    // interpolatedValue = interpolatedValue * unitFactor.getExpression().evaluateConstant(); 		//earlier, hack unitfactor = 0.000001
                    // System.out.println("Sample time: " + sampleTimes[i] + ", between time[" + index + "]=" + data[0][index]+" and time["+(index+1)+"]="+(data[0][index+1])+", interpolated = "+interpolatedValue);
                    outputStream.print("," + interpolatedValue);
                }
                outputStream.println();
            }
        }
        outputStream.close();
        return outputFile;
    } catch (Exception e) {
        e.printStackTrace(System.out);
        // File outputFile = new File(outDir,"results" + filePrefix + ".csv");
        throw new SolverException(e.getMessage());
    }
}

Example 28 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec 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 29 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec in project vcell by virtualcell.

the class OdeFileWriter method write.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
public void write(String[] parameterNames) throws Exception {
    createStateVariables();
    createSymbolTable();
    Simulation simulation = simTask.getSimulation();
    if (simulation.getSolverTaskDescription().getUseSymbolicJacobian()) {
        throw new RuntimeException("symbolic jacobian option not yet supported in interpreted Stiff solver");
    }
    writeJMSParamters();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
    ErrorTolerance errorTolerance = solverTaskDescription.getErrorTolerance();
    printWriter.println("SOLVER " + getSolverName());
    printWriter.println("STARTING_TIME " + timeBounds.getStartingTime());
    printWriter.println("ENDING_TIME " + timeBounds.getEndingTime());
    printWriter.println("RELATIVE_TOLERANCE " + errorTolerance.getRelativeErrorTolerance());
    printWriter.println("ABSOLUTE_TOLERANCE " + errorTolerance.getAbsoluteErrorTolerance());
    printWriter.println("MAX_TIME_STEP " + simulation.getSolverTaskDescription().getTimeStep().getMaximumTimeStep());
    OutputTimeSpec ots = simulation.getSolverTaskDescription().getOutputTimeSpec();
    if (ots.isDefault()) {
        printWriter.println("KEEP_EVERY " + ((DefaultOutputTimeSpec) ots).getKeepEvery());
    } else if (ots.isUniform()) {
        printWriter.println("OUTPUT_TIME_STEP " + ((UniformOutputTimeSpec) ots).getOutputTimeStep());
    } else if (ots.isExplicit()) {
        printWriter.println("OUTPUT_TIMES " + ((ExplicitOutputTimeSpec) ots).getNumTimePoints());
        printWriter.println(((ExplicitOutputTimeSpec) ots).toSpaceSeperatedMultiLinesOfString());
    }
    if (parameterNames != null && parameterNames.length != 0) {
        printWriter.println("NUM_PARAMETERS " + parameterNames.length);
        for (int i = 0; i < parameterNames.length; i++) {
            printWriter.println(parameterNames[i]);
        }
    }
    HashMap<Discontinuity, String> discontinuityNameMap = new HashMap<Discontinuity, String>();
    String eventString = null;
    if (simulation.getMathDescription().hasEvents()) {
        eventString = writeEvents(discontinuityNameMap);
    }
    String equationString = writeEquations(discontinuityNameMap);
    if (discontinuityNameMap.size() > 0) {
        printWriter.println("DISCONTINUITIES " + discontinuityNameMap.size());
        for (Discontinuity od : discontinuityNameMap.keySet()) {
            printWriter.println(discontinuityNameMap.get(od) + " " + od.getDiscontinuityExp().flatten().infix() + "; " + od.getRootFindingExp().flatten().infix() + ";");
        }
    }
    if (eventString != null) {
        printWriter.print(eventString);
    }
    printWriter.println("NUM_EQUATIONS " + getStateVariableCount());
    printWriter.println(equationString);
}

Example 30 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec in project vcell by virtualcell.

the class NFSimSolver method getMathExecutableCommand.

@Override
protected String[] getMathExecutableCommand() {
    String executableName = null;
    try {
        executableName = SolverUtilities.getExes(SolverDescription.NFSim)[0].getAbsolutePath();
    } catch (IOException e) {
        throw new RuntimeException("failed to get executable for solver " + SolverDescription.NFSim.getDisplayLabel() + ": " + e.getMessage(), e);
    }
    String inputFilename = getInputFilename();
    String outputFilename = getOutputFilename();
    String speciesOutputFilename = getSpeciesOutputFilename();
    NFsimSimulationOptions nfsso = simTask.getSimulation().getSolverTaskDescription().getNFSimSimulationOptions();
    ArrayList<String> adv = new ArrayList<String>();
    boolean observableComputationOff = nfsso.getObservableComputationOff();
    if (observableComputationOff == true) {
        // false is by default, no need to specify
        adv.add("-notf");
    }
    Integer moleculeDistance = nfsso.getMoleculeDistance();
    if (moleculeDistance != null) {
        adv.add("-utl");
        adv.add(moleculeDistance + "");
    }
    boolean aggregateBookkeeping = nfsso.getAggregateBookkeeping();
    if (aggregateBookkeeping == true || simTask.getSimulation().getMathDescription().hasSpeciesObservable()) {
        // false is by default, no need to specify
        adv.add("-cb");
    }
    Integer maxMoleculesPerType = nfsso.getMaxMoleculesPerType();
    if (maxMoleculesPerType != null) {
        adv.add("-gml");
        adv.add(maxMoleculesPerType + "");
    }
    Integer equilibrateTime = nfsso.getEquilibrateTime();
    if (equilibrateTime != null) {
        adv.add("-eq");
        adv.add(equilibrateTime + "");
    }
    boolean preventIntraBonds = nfsso.getPreventIntraBonds();
    if (preventIntraBonds == true) {
        // false is by default, no need to specify
        adv.add("-bscb");
    }
    TimeBounds tb = getSimulationJob().getSimulation().getSolverTaskDescription().getTimeBounds();
    double dtime = tb.getEndingTime() - tb.getStartingTime();
    String timeSpecOption1 = "-oSteps";
    String timeSpecOption2 = "10";
    OutputTimeSpec outputTimeSpec = getSimulationJob().getSimulation().getSolverTaskDescription().getOutputTimeSpec();
    if (outputTimeSpec instanceof DefaultOutputTimeSpec) {
        DefaultOutputTimeSpec dots = (DefaultOutputTimeSpec) outputTimeSpec;
        int steps = dots.getKeepAtMost();
        timeSpecOption1 = "-oSteps";
        timeSpecOption2 = Integer.toString(steps);
    } else if (outputTimeSpec instanceof UniformOutputTimeSpec) {
        UniformOutputTimeSpec dots = (UniformOutputTimeSpec) outputTimeSpec;
        double steps = dtime / dots.getOutputTimeStep();
        timeSpecOption1 = "-oSteps";
        int stepsi = (int) Math.round(steps);
        timeSpecOption2 = Integer.toString(stepsi);
    } else {
        throw new RuntimeException("Unsupported output time spec class");
    }
    String[] baseCommands = { "-xml", inputFilename, "-o", outputFilename, "-sim", Double.toString(dtime), "-ss", speciesOutputFilename };
    ArrayList<String> cmds = new ArrayList<String>();
    cmds.add(executableName);
    Integer seed = nfsso.getRandomSeed();
    if (seed != null) {
        cmds.add("-seed");
        cmds.add(seed.toString());
    } else {
        long randomSeed = System.currentTimeMillis();
        randomSeed = randomSeed + simTask.getSimulationJob().getJobIndex();
        // multiply with a large prime number to spread numbers that are too close and in sequence
        randomSeed = randomSeed * 89611;
        Integer rs = (int) randomSeed;
        String str = rs.toString();
        if (str.startsWith("-")) {
            // NFSim wants a positive integer, for anything else is initializing with 0
            str = str.substring(1);
        }
        cmds.add("-seed");
        cmds.add(str);
    // PrintWriter writer;
    // try {
    // writer = new PrintWriter("c:\\TEMP\\aaa\\" + randomSeed + ".txt", "UTF-8");
    // writer.println(str);
    // writer.close();
    // } catch (FileNotFoundException | UnsupportedEncodingException e) {
    // // Auto-generated catch block
    // e.printStackTrace();
    // }
    }
    cmds.add("-vcell");
    cmds.addAll(new ArrayList<String>(Arrays.asList(baseCommands)));
    cmds.add(timeSpecOption1);
    cmds.add(timeSpecOption2);
    cmds.addAll(adv);
    if (bMessaging) {
        cmds.add("-v");
    }
    return cmds.toArray(new String[cmds.size()]);
}