Examples with ArrayDouble ucar.ma2.ArrayDouble used on opensource projects

Example 1 with ArrayDouble

use of ucar.ma2.ArrayDouble in project vcell by virtualcell.

the class NetCDFEvaluator method printAllHistograms.

public void printAllHistograms(int timePointNo) throws IOException {
    if (ncreader != null) {
        try {
            // get all species' names
            String[] names = ncreader.getSpeciesNames_val();
            // get the data for all species at a specific time point
            ArrayDouble data = (ArrayDouble) getDataOverTrials(timePointNo);
            // shape[0]:num of trial, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                System.out.println(names[0] + ":");
                // get one specie's values at a specific time point over trials
                double[] val = new double[shape[0]];
                for (int j = 0; j < shape[0]; j++) {
                    val[j] = temData.get(j);
                // System.out.println("\t"+val[j]);
                }
                Point2D[] histogram = generateHistogram(val);
                for (int k = 0; k < histogram.length; k++) {
                    System.out.println(histogram[k].getX() + "\t" + histogram[k].getY());
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                for (// go through species one by one
                int i = 0; // go through species one by one
                i < shape[1]; // go through species one by one
                i++) {
                    System.out.println(names[i] + ":");
                    // get one specie's values at a specific time point over trials
                    double[] val = new double[shape[0]];
                    for (int j = 0; j < shape[0]; j++) {
                        val[j] = temData.get(j, i);
                    // System.out.println("\t"+val[j]);
                    }
                    Point2D[] histogram = generateHistogram(val);
                    for (int k = 0; k < histogram.length; k++) {
                        System.out.println(histogram[k].getX() + "\t" + histogram[k].getY());
                    }
                }
            }
        } catch (Exception e) {
            e.printStackTrace(System.err);
            throw new IOException("Can not get species' names from model.");
        }
    }
}

Example 2 with ArrayDouble

use of ucar.ma2.ArrayDouble in project vcell by virtualcell.

the class NetCDFWriter method writeHybridInputFile.

/**
 * Write the model to a NetCDF file which serves as an input for stoch hybrid simulator.
 * To write to a NetCDF file is a bit complicated. First, we have to create a NetCDF-3
 * file. And then feed in the data.
 * Creation date: (5/22/2007 5:36:03 PM)
 */
public void writeHybridInputFile(String[] parameterNames) throws Exception, cbit.vcell.parser.ExpressionException, IOException, MathException, InvalidRangeException {
    Simulation simulation = simTask.getSimulation();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (initialize()) {
        // we need to get model and control information first
        NetcdfFileWriteable ncfile = NetcdfFileWriteable.createNew(filename, false);
        // Model info. will be extracted from subDomain of mathDescription
        java.util.Enumeration<SubDomain> e = simulation.getMathDescription().getSubDomains();
        // remember we are dealing with compartmental model here. only 1 subdomain.
        SubDomain subDomain = e.nextElement();
        JumpProcess[] reactions = (JumpProcess[]) subDomain.getJumpProcesses().toArray(new JumpProcess[subDomain.getJumpProcesses().size()]);
        // get species variable names
        Variable[] variables = simSymbolTable.getVariables();
        String[] speciesNames = new String[variables.length];
        for (int i = 0; i < variables.length; i++) speciesNames[i] = variables[i].getName();
        // the probabilities for reactions
        Expression[] probs = new Expression[reactions.length];
        for (int i = 0; i < reactions.length; i++) {
            probs[i] = simSymbolTable.substituteFunctions(reactions[i].getProbabilityRate());
            probs[i] = probs[i].flatten();
        }
        VarIniCondition[] varInis = (VarIniCondition[]) subDomain.getVarIniConditions().toArray(new VarIniCondition[subDomain.getVarIniConditions().size()]);
        // the non-constant stoch variables
        Vector<Variable> vars = new Vector<Variable>();
        for (int i = 0; i < varInis.length; i++) {
            if (varInis[i].getVar() instanceof StochVolVariable) {
                vars.addElement(varInis[i].getVar());
            }
        }
        // get reaction rate law types and rate constants
        ReactionRateLaw[] reactionRateLaws = getReactionRateLaws(probs);
        SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
        TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
        UniformOutputTimeSpec timeSpec = (UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec();
        UniformOutputTimeSpec outputTimeSpec = ((UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec());
        NonspatialStochSimOptions stochOpt = solverTaskDescription.getStochOpt();
        // create an empty NetCDF-3 file
        // define dimensions
        /* these sizes must match the buffers allocated in corresponding Fortran code -- see globalvariables.f90
			in numerics Hy3S/src directory */
        Dimension numTrial = ncfile.addDimension("NumTrials", (int) stochOpt.getNumOfTrials());
        Dimension numSpecies = ncfile.addDimension("NumSpecies", vars.size());
        Dimension numReactions = ncfile.addDimension("NumReactions", subDomain.getJumpProcesses().size());
        int outPoints = ((int) ((timeBounds.getEndingTime() - timeBounds.getStartingTime()) / outputTimeSpec.getOutputTimeStep())) + 1;
        Dimension numTimePoints = ncfile.addDimension("NumTimePoints", outPoints);
        Dimension numModels = ncfile.addDimension("NumModels", 1);
        Dimension numMaxDepList = ncfile.addDimension("NumMaxDepList", 6);
        Dimension numMaxStoichList = ncfile.addDimension("NumMaxStoichList", 25);
        Dimension stringLen = ncfile.addDimension("StringLen", 72);
        // define variables
        // jms info
        ArrayList<Dimension> dims = new ArrayList<Dimension>();
        dims.add(stringLen);
        if (bMessaging) {
            ncfile.addVariable("JMS_BROKER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_USER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_PASSWORD", DataType.CHAR, dims);
            ncfile.addVariable("JMS_QUEUE", DataType.CHAR, dims);
            ncfile.addVariable("JMS_TOPIC", DataType.CHAR, dims);
            ncfile.addVariable("VCELL_USER", DataType.CHAR, dims);
            ncfile.addVariable("SIMULATION_KEY", DataType.INT, new ArrayList<Dimension>());
            ncfile.addVariable("JOB_INDEX", DataType.INT, new ArrayList<Dimension>());
        }
        // scalars
        ncfile.addVariable("TStart", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("TEnd", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("SaveTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("Volume", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTimeSD", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("ExpType", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastTrial", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastModel", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("MaxNumModels", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("NumModels", DataType.INT, new ArrayList<Dimension>());
        // variables with at least 1 dimension
        ArrayList<Dimension> dimspecies = new ArrayList<Dimension>();
        dimspecies.add(numSpecies);
        ArrayList<Dimension> dimreactions = new ArrayList<Dimension>();
        dimreactions.add(numReactions);
        ncfile.addVariable("SpeciesSplitOnDivision", DataType.INT, dimspecies);
        ncfile.addVariable("SaveSpeciesData", DataType.INT, dimspecies);
        ncfile.addVariable("Reaction_Rate_Laws", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_DListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_StoichListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_OptionalData", DataType.INT, dimreactions);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxStoichList);
        ncfile.addVariable("Reaction_StoichCoeff", DataType.INT, dims);
        ncfile.addVariable("Reaction_StoichSpecies", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_DepList", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(stringLen);
        ncfile.addVariable("Reaction_names", DataType.CHAR, dims);
        dims.clear();
        dims.add(numSpecies);
        dims.add(stringLen);
        ncfile.addVariable("Species_names", DataType.CHAR, dims);
        ncfile.addVariable("SpeciesIC", DataType.INT, dimspecies);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_Rate_Constants", DataType.DOUBLE, dims);
        // create the file
        try {
            ncfile.create();
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error creating hybrid file " + filename + ": " + ioe.getMessage());
        }
        // write data to the NetCDF file
        try {
            // write jms info
            if (bMessaging) {
                ArrayChar.D1 jmsString = new ArrayChar.D1(stringLen.getLength());
                String jmshost = PropertyLoader.getRequiredProperty(PropertyLoader.jmsHostExternal);
                // 
                // Used for new REST HTTP messaging api (USE THIS WHEN Hyrbid Solvers are compiled).
                // 
                // String jmsrestport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsRestPortExternal);
                // String jmsurl = jmshost+":"+jmsrestport;
                // 
                // connect to messaging using legacy AMQP protocol instead of new REST api.  Needed for legacy pre-compiled solvers.
                // 
                String jmsport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsPortExternal);
                String jmsurl = "failover:(tcp://" + jmshost + ":" + jmsport + ")";
                jmsString.setString(jmsurl);
                ncfile.write("JMS_BROKER", jmsString);
                jmsString.setString(PropertyLoader.getRequiredProperty(PropertyLoader.jmsUser));
                ncfile.write("JMS_USER", jmsString);
                String jmsPassword = PropertyLoader.getSecretValue(PropertyLoader.jmsPasswordValue, PropertyLoader.jmsPasswordFile);
                jmsString.setString(jmsPassword);
                ncfile.write("JMS_PASSWORD", jmsString);
                jmsString.setString(VCellQueue.WorkerEventQueue.getName());
                ncfile.write("JMS_QUEUE", jmsString);
                jmsString.setString(VCellTopic.ServiceControlTopic.getName());
                ncfile.write("JMS_TOPIC", jmsString);
                jmsString.setString(simulation.getVersion().getOwner().getName());
                ncfile.write("VCELL_USER", jmsString);
                ArrayInt.D0 scalarJMS = new ArrayInt.D0();
                scalarJMS.set(Integer.parseInt(simulation.getVersion().getVersionKey() + ""));
                ncfile.write("SIMULATION_KEY", scalarJMS);
                scalarJMS.set(simTask.getSimulationJob().getJobIndex());
                ncfile.write("JOB_INDEX", scalarJMS);
            }
            ArrayDouble.D0 scalarDouble = new ArrayDouble.D0();
            // TStart, TEnd, SaveTime
            if ((timeBounds.getEndingTime() > timeBounds.getStartingTime()) && (outputTimeSpec.getOutputTimeStep() > 0)) {
                scalarDouble.set(timeBounds.getStartingTime());
                ncfile.write("TStart", scalarDouble);
                scalarDouble.set(timeBounds.getEndingTime());
                ncfile.write("TEnd", scalarDouble);
                scalarDouble.set(outputTimeSpec.getOutputTimeStep());
                ncfile.write("SaveTime", scalarDouble);
            } else {
                System.err.println("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
                throw new RuntimeException("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
            }
            // Volume
            // we set volume to 1. This model file cannot support multi-compartmental sizes.
            // When writting the rate constants, we must take the volume into account according to the reaction type.
            scalarDouble.set(1);
            ncfile.write("Volume", scalarDouble);
            // CellGrowthTime, CellGrowthTimeSD,
            scalarDouble.set(0);
            ncfile.write("CellGrowthTime", scalarDouble);
            ncfile.write("CellGrowthTimeSD", scalarDouble);
            // ExpType, Last Trial, Last Model, MaxNumModels, NumModels
            ArrayInt.D0 scalarInt = new ArrayInt.D0();
            scalarInt.set(0);
            ncfile.write("LastTrial", scalarInt);
            ncfile.write("LastModel", scalarInt);
            scalarInt.set(1);
            ncfile.write("ExpType", scalarInt);
            ncfile.write("MaxNumModels", scalarInt);
            ncfile.write("NumModels", scalarInt);
            // SpeciesSplitOnDivision
            ArrayInt A1 = new ArrayInt.D1(numSpecies.getLength());
            Index idx = A1.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A1.setInt(idx.set(i), 0);
            }
            ncfile.write("SpeciesSplitOnDivision", new int[1], A1);
            // SaveSpeciesData
            ArrayInt A2 = new ArrayInt.D1(numSpecies.getLength());
            idx = A2.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A2.setInt(idx.set(i), 1);
            }
            ncfile.write("SaveSpeciesData", new int[1], A2);
            // Reaction_Rate_Laws
            ArrayInt A3 = new ArrayInt.D1(numReactions.getLength());
            idx = A3.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A3.setInt(idx.set(i), reactionRateLaws[i].getLawType());
            }
            ncfile.write("Reaction_Rate_Laws", new int[1], A3);
            // Reaction_DListLen
            ArrayInt A4 = new ArrayInt.D1(numReactions.getLength());
            idx = A4.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_0)
                    A4.setInt(idx.set(i), 0);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_1) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_1substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_1substrate))
                    A4.setInt(idx.set(i), 1);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_2substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_2substrate))
                    A4.setInt(idx.set(i), 2);
                else if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_3substrate)
                    A4.setInt(idx.set(i), 3);
            }
            ncfile.write("Reaction_DListLen", new int[1], A4);
            // Reaction_StoichListLen
            ArrayInt A5 = new ArrayInt.D1(numReactions.getLength());
            idx = A5.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A5.setInt(idx.set(i), reactions[i].getActions().size());
            }
            ncfile.write("Reaction_StoichListLen", new int[1], A5);
            // Reaction_OptionalData
            ArrayInt A6 = new ArrayInt.D1(numReactions.getLength());
            idx = A6.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A6.setInt(idx.set(i), 0);
            }
            ncfile.write("Reaction_OptionalData", new int[1], A6);
            // Reaction_StoichCoeff
            ArrayInt A7 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A7.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                Action[] actions = (Action[]) reactions[i].getActions().toArray(new Action[reactions[i].getActions().size()]);
                for (int j = 0; j < actions.length; j++) {
                    try {
                        actions[j].getOperand().evaluateConstant();
                        int coeff = (int) Math.round(actions[j].getOperand().evaluateConstant());
                        A7.setInt(idx.set(i, j), coeff);
                    } catch (ExpressionException ex) {
                        ex.printStackTrace(System.err);
                        throw new ExpressionException(ex.getMessage());
                    }
                }
            }
            ncfile.write("Reaction_StoichCoeff", new int[2], A7);
            // Reaction_StoichSpecies
            ArrayInt A8 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A8.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ArrayList<Action> actions = reactions[i].getActions();
                for (int j = 0; j < actions.size(); j++) {
                    A8.setInt(idx.set(i, j), getVariableIndex(((Action) actions.get(j)).getVar().getName(), vars));
                }
            }
            ncfile.write("Reaction_StoichSpecies", new int[2], A8);
            // Reaction_DepList
            ArrayInt A9 = new ArrayInt.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A9.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                Hashtable<String, Integer> tem = varInProbOrderHash[i];
                Enumeration<String> varnames = tem.keys();
                if (rl.getLawType() == ReactionRateLaw.order_0) {
                // don't do anything here.
                } else if ((rl.getLawType() == ReactionRateLaw.order_1) || (rl.getLawType() == ReactionRateLaw.order_2_1substrate) || (rl.getLawType() == ReactionRateLaw.order_3_1substrate) || (rl.getLawType() == ReactionRateLaw.order_2_2substrate) || (rl.getLawType() == ReactionRateLaw.order_3_3substrate)) {
                    int j = 0;
                    while (varnames.hasMoreElements()) {
                        String name = varnames.nextElement();
                        A9.setInt(idx.set(i, j), getVariableIndex(name, vars));
                        j++;
                    }
                } else if (rl.getLawType() == ReactionRateLaw.order_3_2substrate) {
                    int order = 0;
                    String highOrderName = "";
                    String lowOrderName = "";
                    // we must make sure to put the higher order species first.
                    while (varnames.hasMoreElements()) {
                        lowOrderName = varnames.nextElement();
                        if (tem.get(lowOrderName) > order) {
                            String s = highOrderName;
                            highOrderName = lowOrderName;
                            lowOrderName = s;
                            order = tem.get(highOrderName);
                        }
                    }
                    A9.setInt(idx.set(i, 0), getVariableIndex(highOrderName, vars));
                    A9.setInt(idx.set(i, 1), getVariableIndex(lowOrderName, vars));
                }
            }
            ncfile.write("Reaction_DepList", new int[2], A9);
            // Reaction_names
            ArrayChar A10 = new ArrayChar.D2(numReactions.getLength(), stringLen.getLength());
            for (int i = 0; i < numReactions.getLength(); i++) {
                String name = reactions[i].getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A10.setString(i, name);
                } else
                    throw new RuntimeException("Name of Reaction:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Reaction_names", A10);
            // Species_names
            ArrayChar A11 = new ArrayChar.D2(numSpecies.getLength(), stringLen.getLength());
            for (int i = 0; i < numSpecies.getLength(); i++) {
                String name = vars.elementAt(i).getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A11.setString(i, name);
                } else
                    throw new RuntimeException("Name of Species:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Species_names", A11);
            // Species Initial Condition (in number of molecules).
            // Species iniCondition are sampled from a poisson distribution(which has a mean of the current iniExp value)
            RandomDataGenerator dist = new RandomDataGenerator();
            if (stochOpt.isUseCustomSeed()) {
                Integer randomSeed = stochOpt.getCustomSeed();
                if (randomSeed != null) {
                    dist.reSeed(randomSeed);
                }
            }
            ArrayLong A12 = new ArrayLong.D1(numSpecies.getLength());
            idx = A12.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                try {
                    VarIniCondition varIniCondition = subDomain.getVarIniCondition(vars.elementAt(i));
                    Expression varIniExp = varIniCondition.getIniVal();
                    varIniExp.bindExpression(simSymbolTable);
                    varIniExp = simSymbolTable.substituteFunctions(varIniExp).flatten();
                    double expectedCount = varIniExp.evaluateConstant();
                    long varCount = 0;
                    if (varIniCondition instanceof VarIniCount) {
                        varCount = (long) expectedCount;
                    } else {
                        if (expectedCount > 0) {
                            varCount = dist.nextPoisson(expectedCount);
                        }
                    }
                    A12.setLong(idx.set(i), varCount);
                } catch (ExpressionException ex) {
                    ex.printStackTrace(System.err);
                    throw new ExpressionException(ex.getMessage());
                }
            }
            ncfile.write("SpeciesIC", new int[1], A12);
            // Reaction_Rate_Constants(NumReactions, NumMaxDepList) ;
            ArrayDouble A13 = new ArrayDouble.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A13.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                A13.setDouble(idx.set(i, 0), rl.getRateConstant());
            }
            ncfile.write("Reaction_Rate_Constants", A13);
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error writing hybrid input file " + filename + ": " + ioe.getMessage());
        } catch (InvalidRangeException ire) {
            ire.printStackTrace(System.err);
            throw new InvalidRangeException("Error writing hybrid input file " + filename + ": " + ire.getMessage());
        }
        try {
            ncfile.close();
        } catch (IOException ioe) {
            throw new IOException("Error closing file " + filename + ". " + ioe.getMessage());
        }
    }
}

Example 3 with ArrayDouble

use of ucar.ma2.ArrayDouble in project vcell by virtualcell.

the class HybridSolver method getHybridSolverResultSet.

/**
 * Get data columns and function columns to be ready for plot.
 * Creation date: (8/15/2006 11:36:43 AM)
 * @return cbit.vcell.solver.stoch.StochSolverResultSet
 */
private ODESolverResultSet getHybridSolverResultSet() {
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    ODESolverResultSet stSolverResultSet = new ODESolverResultSet();
    try {
        String filename = getBaseName() + NETCDF_DATA_EXTENSION;
        NetCDFEvaluator ncEva = new NetCDFEvaluator();
        NetCDFReader ncReader = null;
        try {
            ncEva.setNetCDFTarget(filename);
            ncReader = ncEva.getNetCDFReader();
        } catch (Exception e) {
            e.printStackTrace(System.err);
            throw new RuntimeException("Cannot open simulation result file: " + filename + "!");
        }
        // Read result according to trial number
        if (ncReader.getNumTrials() == 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("t"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no time points
            ArrayDouble data = (ArrayDouble) ncEva.getTimeSeriesData(1);
            double[] timePoints = ncReader.getTimePoints();
            System.out.println("time points length is " + timePoints.length);
            // shape[0]:num of timepoints, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                System.out.println("one species in time series data and size is " + temData.getSize());
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                System.out.println("multiple species in time series, the length of time series is :" + data.getShape()[0] + ", and the total number of speceis is: " + data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else if (ncReader.getNumTrials() > 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("TrialNo"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no trial numbers
            ArrayDouble data = (ArrayDouble) ncEva.getDataOverTrials(ncReader.getTimePoints().length - 1);
            int[] trialNum = ncEva.getNetCDFReader().getTrialNumbers();
            // System.out.println("total trials are "+trialNum.length);
            // shape[0]:number of trials, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                // System.out.println("one species over trials, size is: "+temData.getSize());
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                // System.out.println("multiple species in multiple trials, the length of trials is :"+data.getShape()[0]+", and the total number of speceis is: "+data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else {
            throw new RuntimeException("Number of trials should be a countable positive value, from 1 to N.");
        }
    } catch (Exception e) {
        e.printStackTrace(System.err);
        throw new RuntimeException("Problem encountered in parsing hybrid simulation results.\n" + e.getMessage());
    }
    /*
	 *Add appropriate Function columns to result set if the stochastic simulation is to display the trajectory.
	 *No function columns for the results of multiple stochastic trials.
	 *In stochastic simulation the functions include probability functions and clamped variable.
	 */
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (simSymbolTable.getSimulation().getSolverTaskDescription().getStochOpt().getNumOfTrials() == 1) {
        Function[] functions = simSymbolTable.getFunctions();
        for (int i = 0; i < functions.length; i++) {
            if (SimulationSymbolTable.isFunctionSaved(functions[i])) {
                Expression exp1 = new Expression(functions[i].getExpression());
                try {
                    exp1 = simSymbolTable.substituteFunctions(exp1);
                } catch (MathException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                }
                try {
                    FunctionColumnDescription cd = new FunctionColumnDescription(exp1.flatten(), functions[i].getName(), null, functions[i].getName(), false);
                    stSolverResultSet.addFunctionColumn(cd);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
    return stSolverResultSet;
}

Example 4 with ArrayDouble

use of ucar.ma2.ArrayDouble in project vcell by virtualcell.

the class ODESimData method readNCDataFile.

public static ODESimData readNCDataFile(VCDataIdentifier vcdId, File dataFile, File functionsFile) throws DataAccessException {
    // read ida file
    System.out.println("reading NetCDF file : " + dataFile);
    ODESimData odeSimData = new ODESimData();
    odeSimData.formatID = NETCDF_DATA_FORMAT_ID;
    odeSimData.mathName = vcdId.getID();
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    try {
        NetCDFEvaluator ncEva = new NetCDFEvaluator();
        NetCDFReader ncReader = null;
        try {
            ncEva.setNetCDFTarget(dataFile.getAbsolutePath());
            ncReader = ncEva.getNetCDFReader();
        } catch (Exception e) {
            e.printStackTrace(System.err);
            throw new RuntimeException("Cannot open simulation result file: " + dataFile.getAbsolutePath() + "!");
        }
        // Read result according to trial number
        if (ncReader.getNumTrials() == 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            odeSimData.addDataColumn(new ODESolverResultSetColumnDescription("t"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                odeSimData.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no time points
            ArrayDouble data = (ArrayDouble) ncEva.getTimeSeriesData(1);
            double[] timePoints = ncReader.getTimePoints();
            System.out.println("time points length is " + timePoints.length);
            // shape[0]:num of timepoints, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                System.out.println("one species in time series data and size is " + temData.getSize());
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[odeSimData.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < odeSimData.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    odeSimData.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                System.out.println("multiple species in time series, the length of time series is :" + data.getShape()[0] + ", and the total number of speceis is: " + data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[odeSimData.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < odeSimData.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    odeSimData.addRow(values);
                }
            }
        } else if (ncReader.getNumTrials() > 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            odeSimData.addDataColumn(new ODESolverResultSetColumnDescription("TrialNo"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                odeSimData.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no trial numbers
            ArrayDouble data = (ArrayDouble) ncEva.getDataOverTrials(ncReader.getTimePoints().length - 1);
            int[] trialNum = ncEva.getNetCDFReader().getTrialNumbers();
            // System.out.println("total trials are "+trialNum.length);
            // shape[0]:number of trials, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                // System.out.println("one species over trials, size is: "+temData.getSize());
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[odeSimData.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < odeSimData.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    odeSimData.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                // System.out.println("multiple species in multiple trials, the length of trials is :"+data.getShape()[0]+", and the total number of speceis is: "+data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[odeSimData.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < odeSimData.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    odeSimData.addRow(values);
                }
            }
        } else {
            throw new RuntimeException("Number of trials should be a countable positive value, from 1 to N.");
        }
    } catch (Exception e) {
        e.printStackTrace(System.err);
        throw new RuntimeException("Problem encountered in parsing hybrid simulation results.\n" + e.getMessage());
    }
    if (!odeSimData.getColumnDescriptions(0).getName().equals(SimDataConstants.HISTOGRAM_INDEX_NAME)) {
        Vector<AnnotatedFunction> funcList;
        try {
            funcList = FunctionFileGenerator.readFunctionsFile(functionsFile, vcdId.getID());
            for (AnnotatedFunction func : funcList) {
                try {
                    Expression expression = new Expression(func.getExpression());
                    odeSimData.addFunctionColumn(new FunctionColumnDescription(expression, func.getName(), null, func.getName(), false));
                } catch (ExpressionException e) {
                    throw new RuntimeException("Could not add function " + func.getName() + " to annotatedFunctionList");
                }
            }
        } catch (FileNotFoundException e1) {
            e1.printStackTrace(System.out);
            throw new DataAccessException(e1.getMessage());
        } catch (IOException e1) {
            e1.printStackTrace(System.out);
            throw new DataAccessException(e1.getMessage());
        }
    }
    return odeSimData;
}