org.apache.commons.math3.random.MersenneTwister Java Examples

/**
 * Test if two clusters are significantly different in the metrics we look at for balancing.
 *
 * @param orig the utilization matrix from the original cluster
 * @param optimized the utilization matrix from the optimized cluster
 * @return The P value that the various derived resources come from the same probability distribution. The probability
 * that the null hypothesis is correct.
 */
public static double[] testDifference(double[][] orig, double[][] optimized) {
  int nResources = RawAndDerivedResource.values().length;
  if (orig.length != nResources) {
    throw new IllegalArgumentException("orig must have number of rows equal to RawAndDerivedResource.");
  }
  if (optimized.length != nResources) {
    throw new IllegalArgumentException("optimized must have number of rows equal to RawAndDerivedResource.");
  }
  if (orig[0].length != optimized[0].length) {
    throw new IllegalArgumentException("The number of brokers must be the same.");
  }

  double[] pValues = new double[orig.length];

  //TODO:  For small N we want to do statistical bootstrapping (not the same as bootstrapping data).
  for (int resourceIndex = 0; resourceIndex < nResources; resourceIndex++) {
    RandomGenerator rng = new MersenneTwister(0x5d11121018463324L);
    KolmogorovSmirnovTest kolmogorovSmirnovTest = new KolmogorovSmirnovTest(rng);
    pValues[resourceIndex] =
        kolmogorovSmirnovTest.kolmogorovSmirnovTest(orig[resourceIndex], optimized[resourceIndex]);
  }

  return pValues;
}
 

public static void drawMnist(DataSet mnist,INDArray reconstruct) throws InterruptedException {
	for(int j = 0; j < mnist.numExamples(); j++) {
		INDArray draw1 = mnist.get(j).getFeatureMatrix().mul(255);
		INDArray reconstructed2 = reconstruct.getRow(j);
		INDArray draw2 = Sampling.binomial(reconstructed2, 1, new MersenneTwister(123)).mul(255);

		DrawReconstruction d = new DrawReconstruction(draw1);
		d.title = "REAL";
		d.draw();
		DrawReconstruction d2 = new DrawReconstruction(draw2,1000,1000);
		d2.title = "TEST";
		
		d2.draw();
		Thread.sleep(1000);
		d.frame.dispose();
		d2.frame.dispose();

	}
}
 

/**
 * Initialises behaviour (determine whether the household will be a BTL investor). Households start off in social
 * housing and with their "desired bank balance" in the bank
 */
public Household(MersenneTwister prng) {
    this.prng = prng; // Passes the Model's random number generator to a private field of each instance
    home = null;
    isFirstTimeBuyer = true;
    isBankrupt = false;
    id = ++id_pool;
    age = data.Demographics.pdfHouseholdAgeAtBirth.nextDouble(this.prng);
    incomePercentile = this.prng.nextDouble();
    behaviour = new HouseholdBehaviour(this.prng, incomePercentile);
    // Find initial values for the annual and monthly gross employment income
    annualGrossEmploymentIncome = data.EmploymentIncome.getAnnualGrossEmploymentIncome(age, incomePercentile);
    monthlyGrossEmploymentIncome = annualGrossEmploymentIncome/config.constants.MONTHS_IN_YEAR;
    bankBalance = behaviour.getDesiredBankBalance(getAnnualGrossTotalIncome()); // Desired bank balance is used as initial value for actual bank balance
    monthlyGrossRentalIncome = 0.0;
}
 

/**
 * @param configFileName String with the address of the configuration file
 * @param outputFolder String with the address of the folder for storing results
 */
public Model(String configFileName, String outputFolder) {
    config = new Config(configFileName);
    prng = new MersenneTwister(config.SEED);

    government = new Government();
    demographics = new Demographics(prng);
    construction = new Construction(prng);
    centralBank = new CentralBank();
    bank = new Bank();
    households = new ArrayList<>(config.TARGET_POPULATION*2);
    houseSaleMarket = new HouseSaleMarket(prng);
    houseRentalMarket = new HouseRentalMarket(prng);

    recorder = new collectors.Recorder(outputFolder);
    transactionRecorder = new collectors.MicroDataRecorder(outputFolder);
    creditSupply = new collectors.CreditSupply(outputFolder);
    coreIndicators = new collectors.CoreIndicators();
    householdStats = new collectors.HouseholdStats();
    housingMarketStats = new collectors.HousingMarketStats(houseSaleMarket);
    rentalMarketStats = new collectors.RentalMarketStats(housingMarketStats, houseRentalMarket);

    nSimulation = 0;
}
 

/**
 * Initialise behavioural variables for a new household: propensity to save, whether the household will have the BTL
    * investor "gene" (provided its income percentile is above a certain minimum), and whether the household will be a
    * fundamentalist or a trend follower investor (provided it has received the BTL investor gene)
 *
 * @param incomePercentile Fixed income percentile for the household (assumed constant over a lifetime)
    */
HouseholdBehaviour(MersenneTwister prng, double incomePercentile) {
	this.prng = prng;  // initialize the random number generator

       // Set downpayment distributions for both first-time-buyers and owner-occupiers
       downpaymentDistFTB = new LogNormalDistribution(this.prng, config.DOWNPAYMENT_FTB_SCALE, config.DOWNPAYMENT_FTB_SHAPE);
       downpaymentDistOO = new LogNormalDistribution(this.prng, config.DOWNPAYMENT_OO_SCALE, config.DOWNPAYMENT_OO_SHAPE);
    // Compute propensity to save, so that it is constant for a given household
       propensityToSave = config.DESIRED_BANK_BALANCE_EPSILON * prng.nextGaussian();
       // Decide if household is a BTL investor and, if so, its tendency to seek capital gains or rental yields
	BTLCapGainCoefficient = 0.0;
       if(incomePercentile > config.MIN_INVESTOR_PERCENTILE &&
               prng.nextDouble() < config.getPInvestor()/config.MIN_INVESTOR_PERCENTILE) {
           BTLInvestor = true;
           if(prng.nextDouble() < config.P_FUNDAMENTALIST) {
               BTLCapGainCoefficient = config.FUNDAMENTALIST_CAP_GAIN_COEFF;
           } else {
               BTLCapGainCoefficient = config.TREND_CAP_GAIN_COEFF;
           }
       } else {
           BTLInvestor = false;
       }
}
 

/**
 * @param func Function to optimize.
 * @param startPoint Starting point.
 * @param inSigma Individual input sigma.
 * @param boundaries Upper / lower point limit.
 * @param goal Minimization or maximization.
 * @param lambda Population size used for offspring.
 * @param isActive Covariance update mechanism.
 * @param diagonalOnly Simplified covariance update.
 * @param stopValue Termination criteria for optimization.
 * @param fTol Tolerance relative error on the objective function.
 * @param pointTol Tolerance for checking that the optimum is correct.
 * @param maxEvaluations Maximum number of evaluations.
 * @param expected Expected point / value.
 */
private void doTest(MultivariateFunction func,
        double[] startPoint,
        double[] inSigma,
        double[][] boundaries,
        GoalType goal,
        int lambda,
        boolean isActive,
        int diagonalOnly, 
        double stopValue,
        double fTol,
        double pointTol,
        int maxEvaluations,
        PointValuePair expected) {
    int dim = startPoint.length;
    // test diagonalOnly = 0 - slow but normally fewer feval#
    CMAESOptimizer optim = new CMAESOptimizer( lambda, inSigma, 30000,
                                               stopValue, isActive, diagonalOnly,
                                               0, new MersenneTwister(), false);
    final double[] lB = boundaries == null ? null : boundaries[0];
    final double[] uB = boundaries == null ? null : boundaries[1];
    PointValuePair result = optim.optimize(maxEvaluations, func, goal, startPoint, lB, uB);
    Assert.assertEquals(expected.getValue(),
            result.getValue(), fTol);
    for (int i = 0; i < dim; i++) {
        Assert.assertEquals(expected.getPoint()[i],
                result.getPoint()[i], pointTol);
    }
}
 

public static List<Vector2D> createRandomPoints(int size) {
    RandomGenerator random = new MersenneTwister();

    // create the cloud container
    List<Vector2D> points = new ArrayList<Vector2D>(size);
    // fill the cloud with a random distribution of points
    for (int i = 0; i < size; i++) {
        points.add(new Vector2D(FastMath.round(random.nextDouble() * 400 + 100),
                FastMath.round(random.nextDouble() * 400 + 100)));
    }
    
    return points;
}
 

@Test
public void testMath864() {
    final CMAESOptimizer optimizer
        = new CMAESOptimizer(30000, 0, true, 10,
                             0, new MersenneTwister(), false, null);
    final MultivariateFunction fitnessFunction = new MultivariateFunction() {
            public double value(double[] parameters) {
                final double target = 1;
                final double error = target - parameters[0];
                return error * error;
            }
        };

    final double[] start = { 0 };
    final double[] lower = { -1e6 };
    final double[] upper = { 1.5 };
    final double[] sigma = { 1e-1 };
    final double[] result = optimizer.optimize(new MaxEval(10000),
                                               new ObjectiveFunction(fitnessFunction),
                                               GoalType.MINIMIZE,
                                               new CMAESOptimizer.PopulationSize(5),
                                               new CMAESOptimizer.Sigma(sigma),
                                               new InitialGuess(start),
                                               new SimpleBounds(lower, upper)).getPoint();
    Assert.assertTrue("Out of bounds (" + result[0] + " > " + upper[0] + ")",
                      result[0] <= upper[0]);
}
 

/**
 * @param func Function to optimize.
 * @param startPoint Starting point.
 * @param inSigma Individual input sigma.
 * @param boundaries Upper / lower point limit.
 * @param goal Minimization or maximization.
 * @param lambda Population size used for offspring.
 * @param isActive Covariance update mechanism.
 * @param diagonalOnly Simplified covariance update.
 * @param stopValue Termination criteria for optimization.
 * @param fTol Tolerance relative error on the objective function.
 * @param pointTol Tolerance for checking that the optimum is correct.
 * @param maxEvaluations Maximum number of evaluations.
 * @param expected Expected point / value.
 */
private void doTest(MultivariateFunction func,
        double[] startPoint,
        double[] inSigma,
        double[][] boundaries,
        GoalType goal,
        int lambda,
        boolean isActive,
        int diagonalOnly, 
        double stopValue,
        double fTol,
        double pointTol,
        int maxEvaluations,
        PointValuePair expected) {
    int dim = startPoint.length;
    // test diagonalOnly = 0 - slow but normally fewer feval#
    CMAESOptimizer optim = new CMAESOptimizer(30000, stopValue, isActive, diagonalOnly,
                                              0, new MersenneTwister(), false, null);
    final double[] lB = boundaries == null ? null : boundaries[0];
    final double[] uB = boundaries == null ? null : boundaries[1];
    PointValuePair result = boundaries == null ?
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda)) :
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new SimpleBounds(lB, uB),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda));
    // System.out.println("sol=" + Arrays.toString(result.getPoint()));
    Assert.assertEquals(expected.getValue(), result.getValue(), fTol);
    for (int i = 0; i < dim; i++) {
        Assert.assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol);
    }
}
 

@Test
public void testMath864() {
    final CMAESOptimizer optimizer
        = new CMAESOptimizer(30000, 0, true, 10,
                             0, new MersenneTwister(), false, null);
    final MultivariateFunction fitnessFunction = new MultivariateFunction() {
            public double value(double[] parameters) {
                final double target = 1;
                final double error = target - parameters[0];
                return error * error;
            }
        };

    final double[] start = { 0 };
    final double[] lower = { -1e6 };
    final double[] upper = { 1.5 };
    final double[] sigma = { 1e-1 };
    final double[] result = optimizer.optimize(new MaxEval(10000),
                                               new ObjectiveFunction(fitnessFunction),
                                               GoalType.MINIMIZE,
                                               new CMAESOptimizer.PopulationSize(5),
                                               new CMAESOptimizer.Sigma(sigma),
                                               new InitialGuess(start),
                                               new SimpleBounds(lower, upper)).getPoint();
    Assert.assertTrue("Out of bounds (" + result[0] + " > " + upper[0] + ")",
                      result[0] <= upper[0]);
}
 

/**
 * @param func Function to optimize.
 * @param startPoint Starting point.
 * @param inSigma Individual input sigma.
 * @param boundaries Upper / lower point limit.
 * @param goal Minimization or maximization.
 * @param lambda Population size used for offspring.
 * @param isActive Covariance update mechanism.
 * @param diagonalOnly Simplified covariance update.
 * @param stopValue Termination criteria for optimization.
 * @param fTol Tolerance relative error on the objective function.
 * @param pointTol Tolerance for checking that the optimum is correct.
 * @param maxEvaluations Maximum number of evaluations.
 * @param expected Expected point / value.
 */
private void doTest(MultivariateFunction func,
        double[] startPoint,
        double[] inSigma,
        double[][] boundaries,
        GoalType goal,
        int lambda,
        boolean isActive,
        int diagonalOnly, 
        double stopValue,
        double fTol,
        double pointTol,
        int maxEvaluations,
        PointValuePair expected) {
    int dim = startPoint.length;
    // test diagonalOnly = 0 - slow but normally fewer feval#
    CMAESOptimizer optim = new CMAESOptimizer(30000, stopValue, isActive, diagonalOnly,
                                              0, new MersenneTwister(), false, null);
    final double[] lB = boundaries == null ? null : boundaries[0];
    final double[] uB = boundaries == null ? null : boundaries[1];
    PointValuePair result = boundaries == null ?
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda)) :
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new SimpleBounds(lB, uB),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda));
    // System.out.println("sol=" + Arrays.toString(result.getPoint()));
    Assert.assertEquals(expected.getValue(), result.getValue(), fTol);
    for (int i = 0; i < dim; i++) {
        Assert.assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol);
    }
}
 

@Override
public MersenneTwister generator() {
	MersenneTwister random = new MersenneTwister();

	random.setSeed(seed);

	return random;
}
 

@Test
public void testMath864() {
    final CMAESOptimizer optimizer
        = new CMAESOptimizer(30000, 0, true, 10,
                             0, new MersenneTwister(), false, null);
    final MultivariateFunction fitnessFunction = new MultivariateFunction() {
            public double value(double[] parameters) {
                final double target = 1;
                final double error = target - parameters[0];
                return error * error;
            }
        };

    final double[] start = { 0 };
    final double[] lower = { -1e6 };
    final double[] upper = { 1.5 };
    final double[] sigma = { 1e-1 };
    final double[] result = optimizer.optimize(new MaxEval(10000),
                                               new ObjectiveFunction(fitnessFunction),
                                               GoalType.MINIMIZE,
                                               new CMAESOptimizer.PopulationSize(5),
                                               new CMAESOptimizer.Sigma(sigma),
                                               new InitialGuess(start),
                                               new SimpleBounds(lower, upper)).getPoint();
    Assert.assertTrue("Out of bounds (" + result[0] + " > " + upper[0] + ")",
                      result[0] <= upper[0]);
}
 

/**
 * @param func Function to optimize.
 * @param startPoint Starting point.
 * @param inSigma Individual input sigma.
 * @param boundaries Upper / lower point limit.
 * @param goal Minimization or maximization.
 * @param lambda Population size used for offspring.
 * @param isActive Covariance update mechanism.
 * @param diagonalOnly Simplified covariance update.
 * @param stopValue Termination criteria for optimization.
 * @param fTol Tolerance relative error on the objective function.
 * @param pointTol Tolerance for checking that the optimum is correct.
 * @param maxEvaluations Maximum number of evaluations.
 * @param expected Expected point / value.
 */
private void doTest(MultivariateFunction func,
        double[] startPoint,
        double[] inSigma,
        double[][] boundaries,
        GoalType goal,
        int lambda,
        boolean isActive,
        int diagonalOnly, 
        double stopValue,
        double fTol,
        double pointTol,
        int maxEvaluations,
        PointValuePair expected) {
    int dim = startPoint.length;
    // test diagonalOnly = 0 - slow but normally fewer feval#
    CMAESOptimizer optim = new CMAESOptimizer(30000, stopValue, isActive, diagonalOnly,
                                              0, new MersenneTwister(), false, null);
    final double[] lB = boundaries == null ? null : boundaries[0];
    final double[] uB = boundaries == null ? null : boundaries[1];
    PointValuePair result = boundaries == null ?
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda)) :
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new SimpleBounds(lB, uB),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda));
    // System.out.println("sol=" + Arrays.toString(result.getPoint()));
    Assert.assertEquals(expected.getValue(), result.getValue(), fTol);
    for (int i = 0; i < dim; i++) {
        Assert.assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol);
    }
}
 

@Test
public void testNormallyDistributedRandomDataShifted() {
  List<Double> values = new ArrayList<>();
  GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(0L));
  for(int i = 0;i < 1000000;++i) {
    double d = gaussian.nextNormalizedDouble() + 10;
    values.add(d);
  }
  validateEquality(values);
}
 

public static List<Vector2D> createRandomPoints(int size) {
    RandomGenerator random = new MersenneTwister();

    // create the cloud container
    List<Vector2D> points = new ArrayList<Vector2D>(size);
    // fill the cloud with a random distribution of points
    for (int i = 0; i < size; i++) {
        points.add(new Vector2D(FastMath.round(random.nextDouble() * 400 + 100),
                FastMath.round(random.nextDouble() * 400 + 100)));
    }
    
    return points;
}
 

@Test
public void testNormallyDistributedRandomDataShiftedBackwards() {
  List<Double> values = new ArrayList<>();
  GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(0L));
  for(int i = 0;i < 1000000;++i) {
    double d = gaussian.nextNormalizedDouble() - 10;
    values.add(d);
  }
  validateEquality(values);
}
 

@Test
public void testNormallyDistributedRandomDataSkewed() {
  List<Double> values = new ArrayList<>();
  GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(0L));
  for(int i = 0;i < 1000000;++i) {
    double d = (gaussian.nextNormalizedDouble()+ 10000) /1000;
    values.add(d);
  }
  validateEquality(values);
}
 

/**
 * @param func Function to optimize.
 * @param startPoint Starting point.
 * @param inSigma Individual input sigma.
 * @param boundaries Upper / lower point limit.
 * @param goal Minimization or maximization.
 * @param lambda Population size used for offspring.
 * @param isActive Covariance update mechanism.
 * @param diagonalOnly Simplified covariance update.
 * @param stopValue Termination criteria for optimization.
 * @param fTol Tolerance relative error on the objective function.
 * @param pointTol Tolerance for checking that the optimum is correct.
 * @param maxEvaluations Maximum number of evaluations.
 * @param expected Expected point / value.
 */
private void doTest(MultivariateFunction func,
        double[] startPoint,
        double[] inSigma,
        double[][] boundaries,
        GoalType goal,
        int lambda,
        boolean isActive,
        int diagonalOnly, 
        double stopValue,
        double fTol,
        double pointTol,
        int maxEvaluations,
        PointValuePair expected) {
    int dim = startPoint.length;
    // test diagonalOnly = 0 - slow but normally fewer feval#
    CMAESOptimizer optim = new CMAESOptimizer(30000, stopValue, isActive, diagonalOnly,
                                              0, new MersenneTwister(), false, null);
    final double[] lB = boundaries == null ? null : boundaries[0];
    final double[] uB = boundaries == null ? null : boundaries[1];
    PointValuePair result = boundaries == null ?
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda)) :
        optim.optimize(maxEvaluations, func, goal,
                       new InitialGuess(startPoint),
                       new SimpleBounds(lB, uB),
                       new CMAESOptimizer.Sigma(inSigma),
                       new CMAESOptimizer.PopulationSize(lambda));
    // System.out.println("sol=" + Arrays.toString(result.getPoint()));
    Assert.assertEquals(expected.getValue(), result.getValue(), fTol);
    for (int i = 0; i < dim; i++) {
        Assert.assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol);
    }
}
 

@Test
public void testMath864() {
    final CMAESOptimizer optimizer
        = new CMAESOptimizer(30000, 0, true, 10,
                             0, new MersenneTwister(), false, null);
    final MultivariateFunction fitnessFunction = new MultivariateFunction() {
            public double value(double[] parameters) {
                final double target = 1;
                final double error = target - parameters[0];
                return error * error;
            }
        };

    final double[] start = { 0 };
    final double[] lower = { -1e6 };
    final double[] upper = { 1.5 };
    final double[] sigma = { 1e-1 };
    final double[] result = optimizer.optimize(new MaxEval(10000),
                                               new ObjectiveFunction(fitnessFunction),
                                               GoalType.MINIMIZE,
                                               new CMAESOptimizer.PopulationSize(5),
                                               new CMAESOptimizer.Sigma(sigma),
                                               new InitialGuess(start),
                                               new SimpleBounds(lower, upper)).getPoint();
    Assert.assertTrue("Out of bounds (" + result[0] + " > " + upper[0] + ")",
                      result[0] <= upper[0]);
}
 

@Override
public MersenneTwister generator() {
	MersenneTwister random = new MersenneTwister();

	random.setSeed(seed);

	return random;
}
 

@Override
public Iterator<List<HStoreFile>> iterator() {
  return new Iterator<List<HStoreFile>>() {
    private GaussianRandomGenerator gen =
        new GaussianRandomGenerator(new MersenneTwister(random.nextInt()));
    private int count = 0;

    @Override
    public boolean hasNext() {
      return count < MAX_FILE_GEN_ITERS;
    }

    @Override
    public List<HStoreFile> next() {
      count += 1;
      ArrayList<HStoreFile> files = new ArrayList<>(NUM_FILES_GEN);
      for (int i = 0; i < NUM_FILES_GEN; i++) {
        files.add(createMockStoreFile(
            (int) Math.ceil(Math.max(0, gen.nextNormalizedDouble() * 32 + 32)))
        );
      }

      return files;
    }

    @Override
    public void remove() {
    }
  };
}
 

HousingMarket(MersenneTwister prng) {
    offersPQ = new PriorityQueue2D<>(new HousingMarketRecord.PQComparator()); //Priority Queue of (Price, Quality)
    // The integer passed to the ArrayList constructor is an initially declared capacity (for initial memory
    // allocation purposes), it will actually have size zero and only grow by adding elements
    // TODO: Check if this integer is too large or small, check speed penalty for using ArrayList as opposed to
    // TODO: normal arrays
    bids = new ArrayList<>(config.TARGET_POPULATION/16);
    this.prng = prng;
}
 

/***
	 * Sample from the PDF
	 * @return A random sample from the PDF
	 */
	public double nextDouble(MersenneTwister rand) {
		return(inverseCumulativeProbability(rand.nextDouble()));
//		double uniform = rand.nextDouble(); // uniform random sample on [0:1)
//		int i = (int)(uniform*(nSamples-1));
//		double remainder = uniform*(nSamples-1.0) - i;
//		return((1.0-remainder)*inverseCDF[i] + remainder*inverseCDF[i+1]);
	}
 

public StableSketch(int sksize, double alpha, long seed) {
    this.alpha = alpha;
    this.sksize = sksize;
    this.seed = seed;
    
    sk = new double[sksize];
    
    randomGen = new MersenneTwister(seed);
    stableGen = new StableRandomGenerator(randomGen, alpha, 0);
}
 

/**
 * Add element b to the stream inc times.
 * @param b The element to be added
 * @param inc How many times the element occurred
 */
public void push(byte[] b, int inc) {
    long it = hasher.hash64(b);
    
    count += inc;

    RandomGenerator r = new MersenneTwister(it);
    
    for (int i=0; i<data.length; i++) {
        double val = maxSkew(r);
        data[i] += val * (double)inc;
    }
}
 

/**
 * @param args
 */
public static void main(String[] args) throws Exception {
	MnistDataSetIterator iter = new MnistDataSetIterator(60,60000);
	@SuppressWarnings("unchecked")
	ObjectInputStream ois = new ObjectInputStream(new FileInputStream(args[0]));
	
	BasePretrainNetwork network = (BasePretrainNetwork) ois.readObject();
	
	
	DataSet test = null;
	while(iter.hasNext()) {
		INDArray reconstructed = network.transform(test.getFeatureMatrix());
		for(int i = 0; i < test.numExamples(); i++) {
			INDArray draw1 = test.get(i).getFeatureMatrix().mul(255);
			INDArray reconstructed2 = reconstructed.getRow(i);
			INDArray draw2 = Sampling.binomial(reconstructed2, 1, new MersenneTwister(123)).mul(255);

			DrawReconstruction d = new DrawReconstruction(draw1);
			d.title = "REAL";
			d.draw();
			DrawReconstruction d2 = new DrawReconstruction(draw2,100,100);
			d2.title = "TEST";
			d2.draw();
			Thread.sleep(10000);
			d.frame.dispose();
			d2.frame.dispose();
		}
	}
	
	
}
 

/**
 * Testing quantiles of normal distribution.
 * Based on feedback provided by Alessandro Gnoatto and a student of him.
 */
@Test
public void testGetQuantile() {

	final int seed = 3141;
	final int numberOfSamplePoints = 10000000;

	final MersenneTwister mersenneTwister = new MersenneTwister(seed);
	final double[] samples = new double[numberOfSamplePoints];
	for(int i = 0; i< numberOfSamplePoints; i++) {
		final double randomNumber = mersenneTwister.nextDouble();
		samples[i] = net.finmath.functions.NormalDistribution.inverseCumulativeDistribution(randomNumber);
	}

	final RandomVariable normalDistributedRandomVariable = abstractRandomVariableFactory.createRandomVariable(0.0,samples);

	final double q00 = normalDistributedRandomVariable.getQuantile(0.0);
	Assert.assertEquals(normalDistributedRandomVariable.getMin(), q00, 1E-12);

	final double q05 = normalDistributedRandomVariable.getQuantile(0.05);
	Assert.assertEquals(-1.645, q05, 1E-3);

	final double q50 = normalDistributedRandomVariable.getQuantile(0.5);
	Assert.assertEquals(0, q50, 2E-4);

	final double q95 = normalDistributedRandomVariable.getQuantile(0.95);
	Assert.assertEquals(1.645, q95, 1E-3);

	final double q99 = normalDistributedRandomVariable.getQuantile(0.99);
	Assert.assertEquals(2.33, q99, 1E-2);
}
 

/**
 * Testing quantiles of normal distribution.
 * Based on feedback provided by Alessandro Gnoatto and a student of him.
 */
@Test
public void testGetQuantile() {

	final int seed = 3141;
	final int numberOfSamplePoints = 10000000;

	final MersenneTwister mersenneTwister = new MersenneTwister(seed);
	final double[] samples = new double[numberOfSamplePoints];
	for(int i = 0; i< numberOfSamplePoints; i++) {
		final double randomNumber = mersenneTwister.nextDouble();
		samples[i] = net.finmath.functions.NormalDistribution.inverseCumulativeDistribution(randomNumber);
	}

	final RandomVariable normalDistributedRandomVariable = abstractRandomVariableFactory.createRandomVariable(0.0,samples);

	final double q00 = normalDistributedRandomVariable.getQuantile(0.0);
	Assert.assertEquals(normalDistributedRandomVariable.getMin(), q00, 1E-12);

	final double q05 = normalDistributedRandomVariable.getQuantile(0.05);
	Assert.assertEquals(-1.645, q05, 1E-3);

	final double q50 = normalDistributedRandomVariable.getQuantile(0.5);
	Assert.assertEquals(0, q50, 2E-4);

	final double q95 = normalDistributedRandomVariable.getQuantile(0.95);
	Assert.assertEquals(1.645, q95, 1E-3);

	final double q99 = normalDistributedRandomVariable.getQuantile(0.99);
	Assert.assertEquals(2.33, q99, 1E-2);
}
 

/**
 * @return a new, seeded {@link RandomGenerator}
 */
public static RandomGenerator getRandom() {
  if (useTestSeed) {
    // No need to track instances anymore
    return new MersenneTwister(TEST_SEED);
  }
  RandomGenerator random = new MersenneTwister();
  synchronized (INSTANCES) {
    INSTANCES.put(random, Boolean.TRUE);
  }
  return random;
}