Examples with UpdateFunction edu.cmu.tetrad.gene.tetrad.gene.history.UpdateFunction used on opensource projects

Example 1 with UpdateFunction

use of edu.cmu.tetrad.gene.tetrad.gene.history.UpdateFunction in project tetrad by cmu-phil.

the class MeasurementSimulator method simulate.

/**
 * Simulates (optionally) neither, either, or both of two three-dimensionaly
 * data sets, rawData and measuredData. For the form of the first set, see
 * <code>getRawData</code>; for the form of the second set, see
 * <code>getMeasuredData</code>.  The idea of the simulation is that cells
 * are grown in different dishes and periodically sampled. When they are
 * sampled, their DNA is ground up and then pipetted in aggregate onto a
 * number of microarray chips. Each indivual cell is simulated separately;
 * the simulated cells are grouped into dishes. Expression levels of cells
 * in each dish are aggregated and then random noise is added to this
 * aggregate expression level representing various sources of noise: sample
 * to sample variation, chip to chip variation, and pixellation measurement
 * variation. For more information, see the Genetic Simulator spec, which
 * this method implements.
 *
 * @see #getRawData
 * @see #getMeasuredData
 */
public void simulate(GeneHistory history) {
    setHistory(history);
    history.reset();
    // Get the update function.
    UpdateFunction updateFunction = history.getUpdateFunction();
    int numFactors = updateFunction.getNumFactors();
    // Set up the dish model.
    DishModel dishModel = new DishModel(getNumDishes(), getDishDishVariability());
    // Set up the history.
    history.setInitSync(isInitSync());
    history.setDishModel(dishModel);
    // We will use two data cubes to store data for display
    // purposes: a measurement cube and a raw data cube. We first
    // construct the measurement cube, with two extra columns for
    // (a) dish number and (b) chip number.
    int numChips = getNumSamplesPerDish();
    if (isMeasuredDataSaved()) {
        int numSteps = timeSteps.length;
        int numSamples = getNumDishes() * numChips;
        this.measuredData = new double[numFactors][numSteps][numSamples];
    }
    // We next construct the raw data cube, with two extra columns
    // for (a) individual number and (b) dish number. Raw data is
    // saved out only if the 'rawDataSave' parameter is true. If
    // it's not true, these objects are not constructed. (They can
    // be extremely large.) References to the cube and the initial
    // columns need to be saved in this scope for future use (but
    // may be null).
    this.rawData = null;
    if (isRawDataSaved()) {
        int numSteps = timeSteps.length;
        int numCells = getNumDishes() * getNumCellsPerDish();
        this.rawData = new double[numFactors][numSteps][numCells];
    }
    // There are three sources of error--sample to sample, chip to
    // chip, and pixel digitalization.  A certain number of
    // samples are drawn from each dish, and these are deposited
    // on the same number of chips. There is an error associated
    // with each chip and with each sample. (Pixel digitalization
    // error is drawn later.)
    double sampleSd = getSampleSampleVariability();
    double chipSd = getChipChipVariability();
    double pixelSd = getPixelDigitalization();
    Distribution sampleErrorDist = new Normal(0.0, sampleSd);
    Distribution chipErrorDist = new Normal(0.0, chipSd);
    Distribution pixelErrorDist = new Normal(0.0, pixelSd);
    double[][] chipErrors = new double[getNumDishes()][numChips];
    double[][] sampleErrors = new double[getNumDishes()][numChips];
    for (int d = 0; d < getNumDishes(); d++) {
        for (int ch = 0; ch < numChips; ch++) {
            chipErrors[d][ch] = chipErrorDist.nextRandom();
            sampleErrors[d][ch] = sampleErrorDist.nextRandom();
        }
    }
    // Make two double[] arrays, one to store data generated for a
    // particular cell (in a particular dish) and one to store
    // aggregations of these data.
    double[][] cellData = new double[timeSteps.length][numFactors];
    double[][] aggregation = new double[cellData.length][cellData[0].length];
    // 'aggregation'. First, iterate over dishes.
    for (int d = 0; d < getNumDishes(); d++) {
        dishModel.setDishNumber(d);
        // Reset the aggregation.
        for (int sIndex = 0; sIndex < timeSteps.length; sIndex++) {
            Arrays.fill(aggregation[sIndex], 0);
        }
        // Next, iterate over the cells in a dish.
        for (int c = 0; c < getNumCellsPerDish(); c++) {
            // 12/01/01
            if ((c + 1) % 50 == 0) {
                this.dishNumber = d;
                this.cellNumber = c;
            // System.out.println("Dish # " + (d + 1) + ", Cell # " +
            // (c + 1));
            }
            // Reset the cell data.
            for (int sIndex = 0; sIndex < timeSteps.length; sIndex++) {
                Arrays.fill(cellData[sIndex], 0);
            }
            // Initialize the history array.
            history.initialize();
            // Generate data for one cell, storing only those time
            // steps that are indicated in the timeSteps
            // array. Note that the timeSteps array is 1-indexed,
            // while the 's' variable is 0-indexed. (Might want to
            // fix this in a future version.)
            int s = -1;
            int sIndex = 0;
            // for the steps in timeSteps[].
            while (++s < getStepsGenerated()) {
                // being synchronized.
                if (s > 0) {
                    history.update();
                }
                if (s == timeSteps[sIndex] - 1) {
                    double[][] historyArray = history.getHistoryArray();
                    // the factors (genes).
                    for (int f = 0; f < numFactors; f++) {
                        // Copy data from the getModel time step in
                        // the history to the cellData[][] array.
                        cellData[sIndex][f] = historyArray[0][f];
                        // Antilog it if necessary.
                        if (isAntilogCalculated()) {
                            cellData[sIndex][f] = Math.exp(cellData[sIndex][f]);
                        }
                        // be null.)
                        if (isRawDataSaved()) {
                            int row = d * getNumCellsPerDish() + c;
                            /*this.rawData[f][sIndex][row] =
                                historyArray[0][f];*/
                            this.rawData[f][sIndex][row] = cellData[sIndex][f];
                        }
                    }
                    if (++sIndex >= timeSteps.length) {
                        break;
                    }
                }
            // END if (s == timeSteps[sIndex])
            }
            // Aggregate data.
            for (int i = 0; i < cellData.length; i++) {
                for (int j = 0; j < cellData[0].length; j++) {
                    aggregation[i][j] += cellData[i][j];
                }
            }
        }
        if (isMeasuredDataSaved()) {
            // generated.)
            for (int sIndex = 0; sIndex < timeSteps.length; sIndex++) {
                for (int f = 0; f < numFactors; f++) {
                    for (int ch = 0; ch < numChips; ch++) {
                        double average = aggregation[sIndex][f] / getNumCellsPerDish();
                        double pixelError = pixelErrorDist.nextRandom();
                        double measurement = average + sampleErrors[d][ch] + chipErrors[d][ch] + pixelError;
                        this.measuredData[f][sIndex][d * numChips + ch] = measurement;
                    }
                // END for (int ch = 0; ch < numChips; ch++)
                }
            // END for (int f = 0; f < history.getNumFacto...
            }
        // END for (int sIndex = 0; sIndex < timeSteps...
        }
    // END if (measuredDataSaved)
    }
    return;
}