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

Example 76 with MathException

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

the class SmoldynFileWriter method writeInitialConcentration.

private int writeInitialConcentration(ParticleInitialConditionConcentration initialConcentration, SubDomain subDomain, Variable variable, String variableName, StringBuilder sb) throws ExpressionException, MathException {
    SimpleSymbolTable simpleSymbolTable = new SimpleSymbolTable(new String[] { ReservedVariable.X.getName(), ReservedVariable.Y.getName(), ReservedVariable.Z.getName() });
    Expression disExpression = new Expression(initialConcentration.getDistribution());
    disExpression.bindExpression(simulationSymbolTable);
    disExpression = simulationSymbolTable.substituteFunctions(disExpression).flatten();
    disExpression.bindExpression(simpleSymbolTable);
    double[] values = new double[3];
    if (dimension == 1) {
        if (disExpression.getSymbolBinding(ReservedVariable.Y.getName()) != null || disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'y' or 'z'", dimension, variable, disExpression));
        }
    } else if (dimension == 2) {
        if (disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'z'", dimension, variable, disExpression));
        }
    }
    int totalCount = 0;
    StringBuilder localsb = new StringBuilder();
    if (subDomain instanceof CompartmentSubDomain) {
        MeshSpecification meshSpecification = simulation.getMeshSpecification();
        ISize sampleSize = meshSpecification.getSamplingSize();
        int numX = sampleSize.getX();
        int numY = dimension < 2 ? 1 : sampleSize.getY();
        int numZ = dimension < 3 ? 1 : sampleSize.getZ();
        boolean bCellCentered = simulation.hasCellCenteredMesh();
        double dx = meshSpecification.getDx(bCellCentered);
        double dy = meshSpecification.getDy(bCellCentered);
        double dz = meshSpecification.getDz(bCellCentered);
        Origin origin = resampledGeometry.getGeometrySpec().getOrigin();
        double ox = origin.getX();
        double oy = origin.getY();
        double oz = origin.getZ();
        Extent extent = resampledGeometry.getExtent();
        double ex = extent.getX();
        double ey = extent.getY();
        double ez = extent.getZ();
        int offset = 0;
        for (int k = 0; k < numZ; k++) {
            double centerz = oz + k * dz;
            double loz = Math.max(oz, centerz - dz / 2);
            double hiz = Math.min(oz + ez, centerz + dz / 2);
            double lz = hiz - loz;
            values[2] = centerz;
            for (int j = 0; j < numY; j++) {
                double centery = oy + j * dy;
                double loy = Math.max(oy, centery - dy / 2);
                double hiy = Math.min(oy + ey, centery + dy / 2);
                values[1] = centery;
                double ly = hiy - loy;
                for (int i = 0; i < numX; i++) {
                    int regionIndex = resampledGeometry.getGeometrySurfaceDescription().getRegionImage().getRegionInfoFromOffset(offset).getRegionIndex();
                    offset++;
                    GeometricRegion region = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions(regionIndex);
                    if (region instanceof VolumeGeometricRegion) {
                        if (!((VolumeGeometricRegion) region).getSubVolume().getName().equals(subDomain.getName())) {
                            continue;
                        }
                    }
                    double centerx = ox + i * dx;
                    double lox = Math.max(ox, centerx - dx / 2);
                    double hix = Math.min(ox + ex, centerx + dx / 2);
                    double lx = hix - lox;
                    values[0] = centerx;
                    double volume = lx;
                    if (dimension > 1) {
                        volume *= ly;
                        if (dimension > 2) {
                            volume *= lz;
                        }
                    }
                    double expectedCount = disExpression.evaluateVector(values) * volume;
                    if (expectedCount <= 0) {
                        continue;
                    }
                    long count = dist.nextPoisson(expectedCount);
                    if (count <= 0) {
                        continue;
                    }
                    totalCount += count;
                    localsb.append(SmoldynVCellMapper.SmoldynKeyword.mol + " " + count + " " + variableName + " " + (float) lox + "-" + (float) hix);
                    if (lg.isDebugEnabled()) {
                        lg.debug("Component subdomain " + variableName + " count " + count);
                    }
                    if (dimension > 1) {
                        localsb.append(" " + loy + "-" + hiy);
                        if (dimension > 2) {
                            localsb.append(" " + loz + "-" + hiz);
                        }
                    }
                    localsb.append("\n");
                }
            }
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.compartment_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    } else if (subDomain instanceof MembraneSubDomain) {
        ArrayList<TrianglePanel> trianglePanelList = membraneSubdomainTriangleMap.get(subDomain);
        for (TrianglePanel trianglePanel : trianglePanelList) {
            Triangle triangle = trianglePanel.triangle;
            switch(dimension) {
                case 1:
                    values[0] = triangle.getNodes(0).getX();
                    break;
                case 2:
                    {
                        double centroidX = triangle.getNodes(0).getX();
                        double centroidY = triangle.getNodes(0).getY();
                        if (triangle.getNodes(0).getX() == triangle.getNodes(1).getX() && triangle.getNodes(0).getY() == triangle.getNodes(1).getY()) {
                            centroidX += triangle.getNodes(2).getX();
                            centroidY += triangle.getNodes(2).getY();
                        } else {
                            centroidX += triangle.getNodes(1).getX();
                            centroidY += triangle.getNodes(1).getY();
                        }
                        values[0] = centroidX / 2;
                        values[1] = centroidY / 2;
                        break;
                    }
                case 3:
                    {
                        double centroidX = triangle.getNodes(0).getX() + triangle.getNodes(1).getX() + triangle.getNodes(2).getX();
                        double centroidY = triangle.getNodes(0).getY() + triangle.getNodes(1).getY() + triangle.getNodes(2).getY();
                        double centroidZ = triangle.getNodes(0).getZ() + triangle.getNodes(1).getZ() + triangle.getNodes(2).getZ();
                        values[0] = centroidX / 3;
                        values[1] = centroidY / 3;
                        values[2] = centroidZ / 3;
                        break;
                    }
            }
            double expectedCount = disExpression.evaluateVector(values) * triangle.getArea();
            if (expectedCount <= 0) {
                continue;
            }
            long count = dist.nextPoisson(expectedCount);
            if (count <= 0) {
                continue;
            }
            totalCount += count;
            if (lg.isDebugEnabled()) {
                lg.debug("Membrane subdomain " + subDomain.getName() + ' ' + variableName + " count " + count);
            }
            localsb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol + " " + count + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + trianglePanel.name + "\n");
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.all + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    }
    if (lg.isDebugEnabled()) {
        lg.debug("Subdomain " + subDomain.getName() + ' ' + variableName + " total count " + totalCount);
    }
    return totalCount;
}

Example 77 with MathException

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

the class CartesianMeshFileReader method readFromFiles.

public CartesianMesh readFromFiles(VCellSimFiles vcellSimFiles) throws IOException, MathException {
    // 
    // read meshFile and parse into 'mesh' object
    // 
    BufferedReader meshReader = null;
    BufferedReader meshMetricsReader = null;
    try {
        meshReader = new BufferedReader(new FileReader(vcellSimFiles.cartesianMeshFile));
        CommentStringTokenizer meshST = new CommentStringTokenizer(meshReader);
        CommentStringTokenizer membraneMeshMetricsST = null;
        if (vcellSimFiles.meshMetricsFile != null) {
            meshMetricsReader = new BufferedReader(new FileReader(vcellSimFiles.meshMetricsFile));
            membraneMeshMetricsST = new CommentStringTokenizer(meshMetricsReader);
        }
        MembraneMeshMetrics membraneMeshMetrics = null;
        SubdomainInfo subdomainInfo = null;
        if (membraneMeshMetricsST != null) {
            membraneMeshMetrics = readMembraneMeshMetrics(membraneMeshMetricsST);
        }
        if (vcellSimFiles.subdomainFile != null) {
            subdomainInfo = SubdomainInfo.read(vcellSimFiles.subdomainFile);
        }
        CartesianMesh mesh = readCartesianMesh(meshST, membraneMeshMetrics, subdomainInfo);
        return mesh;
    } finally {
        if (meshReader != null) {
            try {
                meshReader.close();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        if (meshMetricsReader != null) {
            try {
                meshMetricsReader.close();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
    }
}

Example 78 with MathException

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

the class ChomboVtkFileWriter method writeEmptyMeshFiles.

public File[] writeEmptyMeshFiles(ChomboFiles chomboFiles, File destinationDirectory, ProgressListener progressListener) throws IOException, MathException, DataAccessException {
    ArrayList<File> meshFiles = new ArrayList<File>();
    int timeIndex = 0;
    HashMap<String, VisMesh> domainMeshMap = new HashMap<String, VisMesh>();
    ChomboDataset chomboDataset;
    try {
        chomboDataset = ChomboFileReader.readDataset(chomboFiles, chomboFiles.getTimeIndices().get(timeIndex));
    } catch (Exception e) {
        throw new DataAccessException("failed to read Chombo Dataset: " + e.getMessage(), e);
    }
    for (ChomboCombinedVolumeMembraneDomain chomboCombinedVolumeMembraneDomain : chomboDataset.getCombinedVolumeMembraneDomains()) {
        ChomboMeshData chomboMeshData = chomboCombinedVolumeMembraneDomain.getChomboMeshData();
        ChomboMeshMapping chomboMeshMapping = new ChomboMeshMapping();
        VisMesh visMesh = domainMeshMap.get(chomboCombinedVolumeMembraneDomain.getVolumeDomainName());
        if (visMesh == null) {
            visMesh = chomboMeshMapping.fromMeshData(chomboMeshData, chomboCombinedVolumeMembraneDomain);
            domainMeshMap.put(chomboCombinedVolumeMembraneDomain.getVolumeDomainName(), visMesh);
        }
        String volumeDomainName = chomboCombinedVolumeMembraneDomain.getVolumeDomainName();
        // 
        // write volume mesh file
        // 
        {
            File volumeMeshFile = getVtuMeshFileName(destinationDirectory, chomboFiles, volumeDomainName);
            File chomboIndexDataFile = getChomboIndexDataFileName(destinationDirectory, chomboFiles, volumeDomainName);
            VtkService.getInstance().writeChomboVolumeVtkGridAndIndexData(visMesh, volumeDomainName, volumeMeshFile, chomboIndexDataFile);
            meshFiles.add(volumeMeshFile);
        }
        if (chomboMeshData.getMembraneVarData().size() > 0) {
            // 
            // write membrane mesh file
            // 
            String membraneDomainName = chomboCombinedVolumeMembraneDomain.getMembraneDomainName();
            File membraneMeshFile = getVtuMeshFileName(destinationDirectory, chomboFiles, membraneDomainName);
            File membraneIndexDataFile = getChomboIndexDataFileName(destinationDirectory, chomboFiles, membraneDomainName);
            VtkService.getInstance().writeChomboMembraneVtkGridAndIndexData(visMesh, membraneDomainName, membraneMeshFile, membraneIndexDataFile);
            meshFiles.add(membraneMeshFile);
        }
    }
    return meshFiles.toArray(new File[0]);
}