Search in sources :

Example 6 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class FiniteVolumeFileWriter method convertAnalyticGeometryToRvachevFunction.

public static Expression[] convertAnalyticGeometryToRvachevFunction(GeometrySpec geoSpec) throws ExpressionException {
    SubVolume[] subVolumes = geoSpec.getSubVolumes();
    int numSubVolumes = subVolumes.length;
    Expression[] newExps = new Expression[numSubVolumes];
    if (numSubVolumes == 1) {
        newExps[0] = new Expression(-1.0);
    } else {
        for (int i = 0; i < numSubVolumes - 1; i++) {
            if (!(subVolumes[i] instanceof AnalyticSubVolume)) {
                throw new RuntimeException("Subdomain " + i + " is not a analytic subdomain.");
            }
            AnalyticSubVolume subvolume = (AnalyticSubVolume) subVolumes[i];
            try {
                newExps[i] = RvachevFunctionUtils.convertToRvachevFunction(subvolume.getExpression());
                if (newExps[numSubVolumes - 1] == null) {
                    newExps[numSubVolumes - 1] = Expression.negate(newExps[i]);
                } else {
                    newExps[numSubVolumes - 1] = Expression.max(newExps[numSubVolumes - 1], Expression.negate(newExps[i]));
                }
            } catch (ExpressionException ex) {
                throw new ExpressionException("Domain " + subvolume.getName() + ": " + ex.getMessage());
            }
        }
        for (int i = 1; i < numSubVolumes - 1; i++) {
            for (int j = 0; j < i; j++) {
                newExps[i] = Expression.max(newExps[i], Expression.negate(newExps[j]));
            }
        }
    }
    return newExps;
}
Also used : Expression(cbit.vcell.parser.Expression) SubVolume(cbit.vcell.geometry.SubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) ExpressionException(cbit.vcell.parser.ExpressionException)

Example 7 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeChomboSpec.

private void writeChomboSpec() throws ExpressionException, SolverException, PropertyVetoException, ClassNotFoundException, IOException, GeometryException, ImageException {
    if (!bChomboSolver) {
        return;
    }
    GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
    int dimension = geometrySpec.getDimension();
    if (dimension == 1) {
        throw new SolverException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " is only supported for simulations with 2D or 3D geometry.");
    }
    Simulation simulation = getSimulationTask().getSimulation();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    ChomboSolverSpec chomboSolverSpec = solverTaskDescription.getChomboSolverSpec();
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_BEGIN);
    printWriter.println(FVInputFileKeyword.DIMENSION + " " + geometrySpec.getDimension());
    Extent extent = geometrySpec.getExtent();
    Origin origin = geometrySpec.getOrigin();
    ISize isize = simulation.getMeshSpecification().getSamplingSize();
    switch(geometrySpec.getDimension()) {
        case 2:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY());
            break;
        case 3:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY() + " " + isize.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY() + " " + extent.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY() + " " + origin.getZ());
            break;
    }
    List<CompartmentSubDomain> featureList = new ArrayList<CompartmentSubDomain>();
    Enumeration<SubDomain> enum1 = simulation.getMathDescription().getSubDomains();
    while (enum1.hasMoreElements()) {
        SubDomain sd = enum1.nextElement();
        if (sd instanceof CompartmentSubDomain) {
            featureList.add((CompartmentSubDomain) sd);
        }
    }
    int numFeatures = featureList.size();
    CompartmentSubDomain[] features = featureList.toArray(new CompartmentSubDomain[0]);
    int[] phases = new int[numFeatures];
    Arrays.fill(phases, -1);
    phases[numFeatures - 1] = 0;
    int[] numAssigned = new int[] { 1 };
    assignPhases(features, numFeatures - 1, phases, numAssigned);
    Map<String, Integer> subDomainPhaseMap = new HashMap<String, Integer>();
    for (int i = 0; i < phases.length; ++i) {
        if (phases[i] == -1) {
            throw new SolverException("Failed to assign a phase to CompartmentSubdomain '" + features[i].getName() + "'. It might be caused by too coarsh a mesh.");
        }
        subDomainPhaseMap.put(features[i].getName(), phases[i]);
    }
    SubVolume[] subVolumes = geometrySpec.getSubVolumes();
    if (geometrySpec.hasImage()) {
        Geometry geometry = (Geometry) BeanUtils.cloneSerializable(simulation.getMathDescription().getGeometry());
        Geometry simGeometry = geometry;
        VCImage img = geometry.getGeometrySpec().getImage();
        int factor = Math.max(Math.max(img.getNumX(), img.getNumY()), img.getNumZ()) < 512 ? 2 : 1;
        ISize distanceMapMeshSize = new ISize(img.getNumX() * factor, img.getNumY() * factor, img.getNumZ() * factor);
        Vect3d deltaX = null;
        boolean bCellCentered = false;
        double dx = 0.5;
        double dy = 0.5;
        double dz = 0.5;
        int Nx = distanceMapMeshSize.getX();
        int Ny = distanceMapMeshSize.getY();
        int Nz = distanceMapMeshSize.getZ();
        if (dimension == 2) {
            // pad the 2D image with itself in order to obtain a 3D image used to compute the distance map
            // because the distance map algorithm is 3D only (using distance to triangles)
            byte[] oldPixels = img.getPixels();
            byte[] newPixels = new byte[oldPixels.length * 3];
            System.arraycopy(oldPixels, 0, newPixels, 0, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length * 2, oldPixels.length);
            double distX = geometry.getExtent().getX() / img.getNumX();
            double distY = geometry.getExtent().getY() / img.getNumY();
            // we set the distance on the z axis to something that makes sense
            double distZ = Math.max(distX, distY);
            Extent newExtent = new Extent(geometry.getExtent().getX(), geometry.getExtent().getY(), distZ * 3);
            VCImage newImage = new VCImageUncompressed(null, newPixels, newExtent, img.getNumX(), img.getNumY(), 3);
            // copy the pixel classes too
            ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
            for (VCPixelClass origPixelClass : geometry.getGeometrySpec().getImage().getPixelClasses()) {
                SubVolume origSubvolume = geometry.getGeometrySpec().getImageSubVolumeFromPixelValue(origPixelClass.getPixel());
                newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origPixelClass.getPixel()));
            }
            newImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
            simGeometry = new Geometry(geometry, newImage);
            Nz = 3;
        }
        GeometrySpec simGeometrySpec = simGeometry.getGeometrySpec();
        Extent simExtent = simGeometrySpec.getExtent();
        dx = simExtent.getX() / (Nx - 1);
        dy = simExtent.getY() / (Ny - 1);
        dz = simExtent.getZ() / (Nz - 1);
        if (Math.abs(dx - dy) > 0.1 * Math.max(dx, dy)) {
            dx = Math.min(dx, dy);
            dy = dx;
            Nx = (int) (simExtent.getX() / dx + 1);
            Ny = (int) (simExtent.getY() / dx + 1);
            if (dimension == 3) {
                dz = dx;
                Nz = (int) (simExtent.getZ() / dx + 1);
            }
        }
        deltaX = new Vect3d(dx, dy, dz);
        // one more point in each direction
        distanceMapMeshSize = new ISize(Nx + 1, Ny + 1, Nz + 1);
        Extent distanceMapExtent = new Extent(simExtent.getX() + dx, simExtent.getY() + dy, simExtent.getZ() + dz);
        simGeometrySpec.setExtent(distanceMapExtent);
        GeometrySurfaceDescription geoSurfaceDesc = simGeometry.getGeometrySurfaceDescription();
        geoSurfaceDesc.setVolumeSampleSize(distanceMapMeshSize);
        geoSurfaceDesc.updateAll();
        VCImage vcImage = RayCaster.sampleGeometry(simGeometry, distanceMapMeshSize, bCellCentered);
        SubvolumeSignedDistanceMap[] distanceMaps = DistanceMapGenerator.computeDistanceMaps(simGeometry, vcImage, bCellCentered);
        if (dimension == 2) {
            distanceMaps = DistanceMapGenerator.extractMiddleSlice(distanceMaps);
        }
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + simGeometrySpec.getNumSubVolumes() + " " + FVInputFileKeyword.DISTANCE_MAP);
        for (int i = 0; i < subVolumes.length; i++) {
            File distanceMapFile = new File(workingDirectory, getSimulationTask().getSimulationJobID() + "_" + subVolumes[i].getName() + DISTANCE_MAP_FILE_EXTENSION);
            writeDistanceMapFile(deltaX, distanceMaps[i], distanceMapFile);
            int phase = subDomainPhaseMap.get(subVolumes[i].getName());
            printWriter.println(subVolumes[i].getName() + " " + phase + " " + distanceMapFile.getAbsolutePath());
        }
    } else {
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + geometrySpec.getNumSubVolumes());
        Expression[] rvachevExps = convertAnalyticGeometryToRvachevFunction(geometrySpec);
        for (int i = 0; i < subVolumes.length; i++) {
            if (subVolumes[i] instanceof AnalyticSubVolume) {
                String name = subVolumes[i].getName();
                int phase = subDomainPhaseMap.get(name);
                printWriter.println(name + " " + phase + " ");
                printWriter.println(FVInputFileKeyword.IF + " " + rvachevExps[i].infix() + ";");
                printWriter.println(FVInputFileKeyword.USER + " " + ((AnalyticSubVolume) subVolumes[i]).getExpression().infix() + ";");
            }
        }
    }
    printWriter.println(FVInputFileKeyword.MAX_BOX_SIZE + " " + chomboSolverSpec.getMaxBoxSize());
    printWriter.println(FVInputFileKeyword.FILL_RATIO + " " + chomboSolverSpec.getFillRatio());
    printWriter.println(FVInputFileKeyword.RELATIVE_TOLERANCE + " " + simulation.getSolverTaskDescription().getErrorTolerance().getRelativeErrorTolerance());
    printWriter.println(FVInputFileKeyword.SAVE_VCELL_OUTPUT + " " + chomboSolverSpec.isSaveVCellOutput());
    printWriter.println(FVInputFileKeyword.SAVE_CHOMBO_OUTPUT + " " + chomboSolverSpec.isSaveChomboOutput());
    printWriter.println(FVInputFileKeyword.ACTIVATE_FEATURE_UNDER_DEVELOPMENT + " " + chomboSolverSpec.isActivateFeatureUnderDevelopment());
    printWriter.println(FVInputFileKeyword.SMALL_VOLFRAC_THRESHOLD + " " + chomboSolverSpec.getSmallVolfracThreshold());
    printWriter.println(FVInputFileKeyword.BLOCK_FACTOR + " " + chomboSolverSpec.getBlockFactor());
    printWriter.println(FVInputFileKeyword.TAGS_GROW + " " + chomboSolverSpec.getTagsGrow());
    // Refinement
    int numLevels = chomboSolverSpec.getNumRefinementLevels();
    // Refinements #Levels ratio 1, ratio 2, etc
    printWriter.print(FVInputFileKeyword.REFINEMENTS + " " + (numLevels + 1));
    List<Integer> ratios = chomboSolverSpec.getRefineRatioList();
    for (int i : ratios) {
        printWriter.print(" " + i);
    }
    // write last refinement ratio, fake
    printWriter.println(" 2");
    // membrane rois
    List<RefinementRoi> memRios = chomboSolverSpec.getMembraneRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Membrane + " " + memRios.size());
    for (RefinementRoi roi : memRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        // level tagsGrow ROIexpression
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    List<RefinementRoi> volRios = chomboSolverSpec.getVolumeRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Volume + " " + volRios.size());
    for (RefinementRoi roi : volRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    printWriter.println(FVInputFileKeyword.VIEW_LEVEL + " " + chomboSolverSpec.getViewLevel());
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_END);
    printWriter.println();
}
Also used : Origin(org.vcell.util.Origin) VCPixelClass(cbit.image.VCPixelClass) GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription) Extent(org.vcell.util.Extent) HashMap(java.util.HashMap) ISize(org.vcell.util.ISize) ArrayList(java.util.ArrayList) VCImage(cbit.image.VCImage) SubvolumeSignedDistanceMap(cbit.vcell.geometry.surface.SubvolumeSignedDistanceMap) GeometrySpec(cbit.vcell.geometry.GeometrySpec) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) RefinementRoi(org.vcell.chombo.RefinementRoi) SubVolume(cbit.vcell.geometry.SubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) SolverTaskDescription(cbit.vcell.solver.SolverTaskDescription) VCImageUncompressed(cbit.image.VCImageUncompressed) ChomboSolverSpec(org.vcell.chombo.ChomboSolverSpec) Vect3d(cbit.vcell.render.Vect3d) Geometry(cbit.vcell.geometry.Geometry) Simulation(cbit.vcell.solver.Simulation) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SolverException(cbit.vcell.solver.SolverException) File(java.io.File) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)

Example 8 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class GeometryGuiTest method getExampleGeometryCSG.

public static Geometry getExampleGeometryCSG() throws PropertyVetoException, ExpressionException, GeometryException, ImageException {
    // translated rotated cube
    CSGPrimitive cube = new CSGPrimitive("cube", CSGPrimitive.PrimitiveType.CUBE);
    CSGRotation rotatedCube = new CSGRotation("Rotation", new Vect3d(1, 2, 3), Math.PI / 4.0);
    rotatedCube.setChild(cube);
    // translated sphere
    CSGTranslation translatedSphere = new CSGTranslation("translation", new Vect3d(0.5, 0.5, 0.5));
    CSGPrimitive sphere = new CSGPrimitive("sphere", CSGPrimitive.PrimitiveType.SPHERE);
    translatedSphere.setChild(sphere);
    // union
    CSGSetOperator csgSetOperator = new CSGSetOperator("difference", OperatorType.DIFFERENCE);
    csgSetOperator.addChild(rotatedCube);
    csgSetOperator.addChild(translatedSphere);
    // scaled union
    CSGScale csgScale = new CSGScale("scale", new Vect3d(3, 3, 3));
    csgScale.setChild(csgSetOperator);
    CSGTranslation csgTranslatedUnion = new CSGTranslation("translationUnion", new Vect3d(5, 5, 5));
    csgTranslatedUnion.setChild(csgScale);
    Geometry geometry = new Geometry("csg", 3);
    CSGObject csgObject = new CSGObject(null, "obj1", 1);
    csgObject.setRoot(csgTranslatedUnion);
    geometry.getGeometrySpec().addSubVolume(new AnalyticSubVolume("background", new Expression(1.0)));
    geometry.getGeometrySpec().addSubVolume(csgObject, true);
    geometry.refreshDependencies();
    geometry.precomputeAll(new GeometryThumbnailImageFactoryAWT());
    return geometry;
}
Also used : Geometry(cbit.vcell.geometry.Geometry) GeometryThumbnailImageFactoryAWT(cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT) CSGTranslation(cbit.vcell.geometry.CSGTranslation) Expression(cbit.vcell.parser.Expression) CSGPrimitive(cbit.vcell.geometry.CSGPrimitive) CSGObject(cbit.vcell.geometry.CSGObject) CSGRotation(cbit.vcell.geometry.CSGRotation) CSGSetOperator(cbit.vcell.geometry.CSGSetOperator) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) Vect3d(cbit.vcell.render.Vect3d) CSGScale(cbit.vcell.geometry.CSGScale)

Example 9 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class GeometrySubVolumeTableModel method setValueAt.

public void setValueAt(Object aValue, int rowIndex, int columnIndex) {
    if (rowIndex < 0 || rowIndex >= getRowCount()) {
        throw new RuntimeException("GeometrySubVolumeTableModel.setValueAt(), row = " + rowIndex + " out of range [" + 0 + "," + (getRowCount() - 1) + "]");
    }
    if (columnIndex < 0 || columnIndex >= getColumnCount()) {
        throw new RuntimeException("GeometrySubVolumeTableModel.setValueAt(), column = " + columnIndex + " out of range [" + 0 + "," + (getColumnCount() - 1) + "]");
    }
    final SubVolume subVolume = getValueAt(rowIndex);
    try {
        switch(columnIndex) {
            case COLUMN_NAME:
                {
                    final String newName = (String) aValue;
                    final String oldName = subVolume.getName();
                    subVolume.setName(newName);
                    AsynchClientTask task1 = new AsynchClientTask("changing the name", AsynchClientTask.TASKTYPE_NONSWING_BLOCKING) {

                        @Override
                        public void run(Hashtable<String, Object> hashTable) throws Exception {
                            getGeometry().precomputeAll(new GeometryThumbnailImageFactoryAWT());
                            GeometrySpec gs = getGeometry().getGeometrySpec();
                            if (gs != null) {
                                gs.geometryNameChanged(oldName, newName);
                            } else {
                                if (lg.isEnabledFor(Level.WARN)) {
                                    lg.warn(getGeometry().getDescription() + " has no GeometrySpec?");
                                }
                            }
                        }
                    };
                    ClientTaskDispatcher.dispatch(ownerTable, new Hashtable<String, Object>(), new AsynchClientTask[] { task1 }, false);
                    break;
                }
            case COLUMN_VALUE:
                {
                    if (subVolume instanceof AnalyticSubVolume) {
                        final AnalyticSubVolume analyticSubVolume = (AnalyticSubVolume) subVolume;
                        if (aValue instanceof ScopedExpression) {
                            throw new RuntimeException("unexpected value type ScopedExpression");
                        } else if (aValue instanceof String) {
                            final String newExpressionString = (String) aValue;
                            analyticSubVolume.setExpression(new Expression(newExpressionString));
                            AsynchClientTask task1 = new AsynchClientTask("changing the expression", AsynchClientTask.TASKTYPE_NONSWING_BLOCKING) {

                                @Override
                                public void run(Hashtable<String, Object> hashTable) throws Exception {
                                    getGeometry().precomputeAll(new GeometryThumbnailImageFactoryAWT());
                                }
                            };
                            ClientTaskDispatcher.dispatch(ownerTable, new Hashtable<String, Object>(), new AsynchClientTask[] { task1 }, false);
                        }
                    }
                    break;
                }
        }
    } catch (Exception e) {
        e.printStackTrace(System.out);
        DialogUtils.showErrorDialog(ownerTable, e.getMessage(), e);
    }
}
Also used : AsynchClientTask(cbit.vcell.client.task.AsynchClientTask) Hashtable(java.util.Hashtable) GeometrySpec(cbit.vcell.geometry.GeometrySpec) GeometryThumbnailImageFactoryAWT(cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT) ScopedExpression(cbit.gui.ScopedExpression) ScopedExpression(cbit.gui.ScopedExpression) Expression(cbit.vcell.parser.Expression) SubVolume(cbit.vcell.geometry.SubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) CSGObject(cbit.vcell.geometry.CSGObject) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)

Example 10 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class Generate2DExpModelOpAbstract method generateModel.

public final GeneratedModelResults generateModel(double deltaX, double bleachRadius, double cellRadius, double bleachDuration, double bleachRate, double postbleachDelay, double postbleachDuration, double psfSigma, double outputTimeStep, double primaryDiffusionRate, double primaryFraction, double bleachMonitorRate, double secondaryDiffusionRate, double secondaryFraction, String extracellularName, String cytosolName, Context context) throws PropertyVetoException, ExpressionException, GeometryException, ImageException, ModelException, MappingException, MathException, MatrixException {
    double domainSize = 2.2 * cellRadius;
    Extent extent = new Extent(domainSize, domainSize, 1.0);
    Origin origin = new Origin(-extent.getX() / 2.0, -extent.getY() / 2.0, -extent.getZ() / 2.0);
    String EXTRACELLULAR_NAME = extracellularName;
    String CYTOSOL_NAME = cytosolName;
    AnalyticSubVolume cytosolSubVolume = new AnalyticSubVolume(CYTOSOL_NAME, new Expression("pow(x,2)+pow(y,2)<pow(" + cellRadius + ",2)"));
    AnalyticSubVolume extracellularSubVolume = new AnalyticSubVolume(EXTRACELLULAR_NAME, new Expression(1.0));
    Geometry geometry = new Geometry("geometry", 2);
    geometry.getGeometrySpec().setExtent(extent);
    geometry.getGeometrySpec().setOrigin(origin);
    geometry.getGeometrySpec().addSubVolume(extracellularSubVolume);
    geometry.getGeometrySpec().addSubVolume(cytosolSubVolume, true);
    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);
    SpeciesContext immobileSC = model.createSpeciesContext(cytosol);
    SpeciesContext primarySC = model.createSpeciesContext(cytosol);
    SpeciesContext secondarySC = model.createSpeciesContext(cytosol);
    // 
    // common bleaching rate for all species
    // 
    double bleachStart = 10 * outputTimeStep - bleachDuration - postbleachDelay;
    double bleachEnd = bleachStart + bleachDuration;
    Expression bleachRateExp = createBleachExpression(bleachRadius, bleachRate, bleachMonitorRate, bleachStart, bleachEnd);
    {
        SimpleReaction immobileBWM = model.createSimpleReaction(cytosol);
        GeneralKinetics immobileBWMKinetics = new GeneralKinetics(immobileBWM);
        immobileBWM.setKinetics(immobileBWMKinetics);
        immobileBWM.addReactant(immobileSC, 1);
        immobileBWMKinetics.getReactionRateParameter().setExpression(Expression.mult(bleachRateExp, new Expression(immobileSC.getName())));
    }
    {
        SimpleReaction primaryBWM = model.createSimpleReaction(cytosol);
        GeneralKinetics primaryBWMKinetics = new GeneralKinetics(primaryBWM);
        primaryBWM.setKinetics(primaryBWMKinetics);
        primaryBWM.addReactant(primarySC, 1);
        primaryBWMKinetics.getReactionRateParameter().setExpression(Expression.mult(bleachRateExp, new Expression(primarySC.getName())));
    }
    {
        SimpleReaction secondaryBWM = model.createSimpleReaction(cytosol);
        GeneralKinetics secondaryBWMKinetics = new GeneralKinetics(secondaryBWM);
        secondaryBWM.setKinetics(secondaryBWMKinetics);
        secondaryBWM.addReactant(secondarySC, 1);
        secondaryBWMKinetics.getReactionRateParameter().setExpression(Expression.mult(bleachRateExp, new Expression(secondarySC.getName())));
    }
    // create simulation context
    SimulationContext simContext = bioModel.addNewSimulationContext("simContext", SimulationContext.Application.NETWORK_DETERMINISTIC);
    simContext.setGeometry(geometry);
    FeatureMapping cytosolFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(cytosol);
    FeatureMapping extracellularFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(extracellular);
    SubVolume cytSubVolume = geometry.getGeometrySpec().getSubVolume(CYTOSOL_NAME);
    SubVolume exSubVolume = geometry.getGeometrySpec().getSubVolume(EXTRACELLULAR_NAME);
    // unused? SurfaceClass pmSurfaceClass = geometry.getGeometrySurfaceDescription().getSurfaceClass(exSubVolume, cytSubVolume);
    cytosolFeatureMapping.setGeometryClass(cytSubVolume);
    extracellularFeatureMapping.setGeometryClass(exSubVolume);
    cytosolFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    extracellularFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
    double fixedFraction = 1.0 - primaryFraction - secondaryFraction;
    SpeciesContextSpec immobileSCS = simContext.getReactionContext().getSpeciesContextSpec(immobileSC);
    immobileSCS.getInitialConditionParameter().setExpression(new Expression(fixedFraction));
    immobileSCS.getDiffusionParameter().setExpression(new Expression(0.0));
    SpeciesContextSpec primarySCS = simContext.getReactionContext().getSpeciesContextSpec(primarySC);
    primarySCS.getInitialConditionParameter().setExpression(new Expression(primaryFraction));
    primarySCS.getDiffusionParameter().setExpression(new Expression(primaryDiffusionRate));
    SpeciesContextSpec secondarySCS = simContext.getReactionContext().getSpeciesContextSpec(secondarySC);
    secondarySCS.getInitialConditionParameter().setExpression(new Expression(secondaryFraction));
    secondarySCS.getDiffusionParameter().setExpression(new Expression(secondaryDiffusionRate));
    simContext.getMicroscopeMeasurement().addFluorescentSpecies(immobileSC);
    simContext.getMicroscopeMeasurement().addFluorescentSpecies(primarySC);
    simContext.getMicroscopeMeasurement().addFluorescentSpecies(secondarySC);
    simContext.getMicroscopeMeasurement().setConvolutionKernel(new GaussianConvolutionKernel(new Expression(psfSigma), new Expression(psfSigma)));
    MathMapping mathMapping = simContext.createNewMathMapping();
    MathDescription mathDesc = mathMapping.getMathDescription();
    simContext.setMathDescription(mathDesc);
    User owner = context.getDefaultOwner();
    int meshSize = (int) (domainSize / deltaX);
    if (meshSize % 2 == 0) {
        // want an odd-sized mesh in x and y ... so centered at the origin.
        meshSize = meshSize + 1;
    }
    TimeBounds timeBounds = new TimeBounds(0.0, postbleachDuration);
    // 
    // simulation to use for data generation (output time steps as recorded by the microscope)
    // 
    double bleachBlackoutBegin = bleachStart - postbleachDelay;
    double bleachBlackoutEnd = bleachEnd + postbleachDelay;
    // ArrayList<Double> times = new ArrayList<Double>();
    // double time = 0;
    // while (time<=timeBounds.getEndingTime()){
    // if (time<=bleachBlackoutBegin || time>bleachBlackoutEnd){
    // // postbleachDelay is the time it takes to switch the filters.
    // times.add(time);
    // }
    // time += outputTimeStep.getData();
    // }
    // double[] timeArray = new double[times.size()];
    // for (int i=0;i<timeArray.length;i++){
    // timeArray[i] = times.get(i);
    // }
    // OutputTimeSpec fakeDataSimOutputTimeSpec = new ExplicitOutputTimeSpec(timeArray);
    OutputTimeSpec fakeDataSimOutputTimeSpec = new UniformOutputTimeSpec(outputTimeStep);
    KeyValue fakeDataSimKey = context.createNewKeyValue();
    SimulationVersion fakeDataSimVersion = new SimulationVersion(fakeDataSimKey, "fakeDataSim", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
    Simulation fakeDataSim = new Simulation(fakeDataSimVersion, mathDesc);
    simContext.addSimulation(fakeDataSim);
    fakeDataSim.getSolverTaskDescription().setTimeBounds(timeBounds);
    fakeDataSim.getMeshSpecification().setSamplingSize(new ISize(meshSize, meshSize, 1));
    fakeDataSim.getSolverTaskDescription().setSolverDescription(SolverDescription.SundialsPDE);
    fakeDataSim.getSolverTaskDescription().setOutputTimeSpec(fakeDataSimOutputTimeSpec);
    // 
    // simulation to use for viewing the protocol (output time steps to understand the physics)
    // 
    KeyValue fullExperimentSimKey = context.createNewKeyValue();
    SimulationVersion fullExperimentSimVersion = new SimulationVersion(fullExperimentSimKey, "fullExperiment", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
    Simulation fullExperimentSim = new Simulation(fullExperimentSimVersion, mathDesc);
    simContext.addSimulation(fullExperimentSim);
    OutputTimeSpec fullExperimentOutputTimeSpec = new UniformOutputTimeSpec(outputTimeStep / 10.0);
    fullExperimentSim.getSolverTaskDescription().setTimeBounds(timeBounds);
    fullExperimentSim.getMeshSpecification().setSamplingSize(new ISize(meshSize, meshSize, 1));
    fullExperimentSim.getSolverTaskDescription().setSolverDescription(SolverDescription.SundialsPDE);
    fullExperimentSim.getSolverTaskDescription().setOutputTimeSpec(fullExperimentOutputTimeSpec);
    GeneratedModelResults results = new GeneratedModelResults();
    results.bioModel_2D = bioModel;
    results.simulation_2D = fakeDataSim;
    results.bleachBlackoutBeginTime = bleachBlackoutBegin;
    results.bleachBlackoutEndTime = bleachBlackoutEnd;
    return results;
}
Also used : Origin(org.vcell.util.Origin) User(org.vcell.util.document.User) KeyValue(org.vcell.util.document.KeyValue) Extent(org.vcell.util.Extent) MathDescription(cbit.vcell.math.MathDescription) ISize(org.vcell.util.ISize) SpeciesContext(cbit.vcell.model.SpeciesContext) GeneralKinetics(cbit.vcell.model.GeneralKinetics) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) Feature(cbit.vcell.model.Feature) TimeBounds(cbit.vcell.solver.TimeBounds) GroupAccessNone(org.vcell.util.document.GroupAccessNone) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) OutputTimeSpec(cbit.vcell.solver.OutputTimeSpec) SimulationVersion(org.vcell.util.document.SimulationVersion) FeatureMapping(cbit.vcell.mapping.FeatureMapping) SubVolume(cbit.vcell.geometry.SubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) SimpleReaction(cbit.vcell.model.SimpleReaction) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) SimulationContext(cbit.vcell.mapping.SimulationContext) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) BigDecimal(java.math.BigDecimal) Date(java.util.Date) Geometry(cbit.vcell.geometry.Geometry) Simulation(cbit.vcell.solver.Simulation) Expression(cbit.vcell.parser.Expression) BioModel(cbit.vcell.biomodel.BioModel) BioModel(cbit.vcell.biomodel.BioModel) Model(cbit.vcell.model.Model) MathMapping(cbit.vcell.mapping.MathMapping) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)

Aggregations

AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)20 Expression (cbit.vcell.parser.Expression)14 CSGObject (cbit.vcell.geometry.CSGObject)13 SubVolume (cbit.vcell.geometry.SubVolume)11 GeometryThumbnailImageFactoryAWT (cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT)10 Geometry (cbit.vcell.geometry.Geometry)8 ImageSubVolume (cbit.vcell.geometry.ImageSubVolume)7 Extent (org.vcell.util.Extent)7 AsynchClientTask (cbit.vcell.client.task.AsynchClientTask)6 Hashtable (java.util.Hashtable)6 Origin (org.vcell.util.Origin)6 GeometrySpec (cbit.vcell.geometry.GeometrySpec)5 ExpressionException (cbit.vcell.parser.ExpressionException)5 ISize (org.vcell.util.ISize)5 ImageException (cbit.image.ImageException)4 VCImage (cbit.image.VCImage)4 BioModel (cbit.vcell.biomodel.BioModel)4 GeometrySurfaceDescription (cbit.vcell.geometry.surface.GeometrySurfaceDescription)4 SimulationContext (cbit.vcell.mapping.SimulationContext)4 GeometryClass (cbit.vcell.geometry.GeometryClass)3