public enum DISTRIBUTION {
	Lognormal, Normal, Student, Uniform
} public interface MultivariateDistribution {
	
	public double[] getRealization();
	
	public double[] getRealization(double[] u);
	
	public int getDim();
	
	public List<double[]> getSample(int size);
	
	public void setGenerator(Random generator);
	
	public Random getGenerator();
	
	public double[] getMean();
	
	public double[] getVariance();
	
	public double[][] getCov();
	
	public double[] getSkewness();
	
	public double[] getKurtosis();
	
	public DISTRIBUTION getType();

} public class MultivariateLognormal extends MultivariateNormal {
	
	double[] _logmean;
	double[][] _logcov;
	double[] _logvar;
	double[] _logskew;
	double[] _logkurt;
	
	public MultivariateLognormal(double[] mean, double[][] covariance, Random generator) {
		super(mean,covariance,generator);
		_logmean = new double[_dim];
		_logcov = new double[_dim][_dim];
		_logskew = new double[_dim];
		_logkurt = new double[_dim];
		for (int i=0; i<mean.length; i++) {
			_logmean[i] = Math.exp(mean[i]+covariance[i][i]/2);
			for (int j=0; j<mean.length; j++)
				_logcov[i][j] = (Math.exp(covariance[i][j])-1)*Math.exp(2*mean[i]+covariance[i][j]);
			_logskew[i] = (Math.exp(covariance[i][i])+2)*Math.sqrt(Math.exp(covariance[i][i])-1);
			_logkurt[i] = Math.exp(4*covariance[i][i])+2*Math.exp(3*covariance[i][i])+3*Math.exp(2*covariance[i][i])-6;
		}
		_logvar = new double[_dim];
		for (int i=0; i<_dim; i++)
			_logvar[i] = _logcov[i][i];
	}
	
	@Override
	public double[] transform(double[] z) {
		double[] x = super.transform(z);
		for (int i=0; i<_dim; i++)
			x[i] = Math.exp(x[i]);
		return x;
	}
	
	@Override
	public double[] getMean() {
		return _logmean;
	}

	@Override
	public double[][] getCov() {
		return _logcov;
	}
	
	@Override
	public double[] getVariance() {
		return _logvar;
	}

	@Override
	public double[] getSkewness() {
		return _logskew;
	}

	@Override
	public double[] getKurtosis() {
		return _logkurt;
	}

	@Override
	public DISTRIBUTION getType() {
		return DISTRIBUTION.Lognormal;
	}
	
	public double[] getNormalMean() {
		return super._mean;
	}
	
	public double[][] getNormalCov() {
		return super._cov;
	}

} public class MultivariateNormal implements MultivariateDistribution {
	
	double[][] _chol;
	double[][] _cov;
	double[] _var;
	double[] _mean;
	int _dim;
	Random _generator;
	
	public MultivariateNormal(double[] mean, double[][] covariance, Random generator) {
		if (covariance.length != mean.length || covariance[0].length != mean.length)
	         throw new IllegalArgumentException("Covariance matrix dimensions invalid.");
		_mean = mean;
		_chol = choleskyDecomposition(covariance);
		_cov = covariance;
		_dim = mean.length;
		_var = new double[_dim];
		for (int i=0; i<_dim; i++)
			_var[i] = _cov[i][i];
		_generator = generator;
	}
	
	public static void main (String... args) {
		int dim = 1;
		double[] mean = new double[dim];
		double[][] correl = new double[dim][dim];
		for (int i=0; i<dim; i++) {
			mean[i] = 10;
			correl[i][i] = 2.5*2.5;
		}
		ReductionMethod mm = new MonteCarlo(new MultivariateNormal(mean,correl,new Random()));
		Map<double[],Double> map = mm.getScenarios(20);
		for (double[] x : map.keySet()) {
			for (int i=0; i<x.length; i++)
				System.out.print(x[i]+"\t");
			System.out.println();
		}
	}
	
	
	public void setSeed(long n) {
		_generator.setSeed(n);
	}
	
	//cholesky-banachiewicz algorithm
	static double[][] choleskyDecomposition(double[][] a){
		int m = a.length;
		double[][] l = new double[m][m]; //automatically initialzed to 0's
		for(int i = 0; i< m; i++){
			for(int k = 0; k < (i+1); k++){
				double sum = 0;
				for(int j = 0; j < k; j++){
					sum += l[i][j] * l[k][j];
				}
				l[i][k] = (i == k) ? Math.sqrt(a[i][i] - sum) :
					(1.0 / l[k][k] * (a[i][k] - sum));
			}
		}
		return l;
	}
	
	public double[] getRealization() {
		double[] z = new double[_dim];
		for (int i=0; i<_dim; i++)
			z[i] = _generator.nextGaussian();
		return transform(z);
	}
	
	public double[] transform(double[] z) {
		double[] x = new double[_chol.length];
		for (int i=0; i<_dim; i++) {
			double sum = 0.0;
			for (int j=0; j<i+1; j++)
				sum += z[j]*_chol[i][j];
			x[i] = _mean[i] + sum;
		}
		return x;
	}

	@Override
	public int getDim() {
		return _dim;
	}

	@Override
	public List<double[]> getSample(int size) {
		List<double[]> list = new LinkedList<>();
		for (int i=0; i<size; i++)
			list.add(getRealization());
		return list;
	}

	@Override
	public void setGenerator(Random generator) {
		_generator = generator;
	}

	@Override
	public Random getGenerator() {
		return _generator;
	}

	@Override
	public double[] getMean() {
		return _mean;
	}

	@Override
	public double[][] getCov() {
		return _cov;
	}
	
	public double[][] getCholesky() {
		return _chol;
	}

	@Override
	public double[] getSkewness() {
		return new double[_dim];
	}

	@Override
	public double[] getKurtosis() {
		return new double[_dim];
	}

	@Override
	public DISTRIBUTION getType() {
		return DISTRIBUTION.Normal;
	}

	@Override
	public double[] getRealization(double[] u) {
		double[] z = new double[_dim];
		for (int i=0; i<_dim; i++)
			z[i] = NormalDistQuick.inverseF01(u[i]);
		return transform(z);
	}

	@Override
	public double[] getVariance() {
		return _var;
	}
	
	public boolean isIndependent() {
		for (int i=0; i<_dim; i++) {
			for (int j=i+1; j<_dim; j++) {
				if (_cov[i][j]!=0)
					return false;
			}
		}
		return true;
	}

} public class MultivariateStudent extends MultivariateNormal {
	
	double[][] _cov;
	double[][] _correlation;
	double[] _kurt;
	double[] _scale;
	int[] _df;
	
	public MultivariateStudent(double[] mean, double[] scale, double[][] correlation, int[] df, Random generator) {
		super(mean,correlation,generator);
		_df = df;
		_cov = new double[_dim][_dim];
		_kurt = new double[_dim];
		_scale = scale;
		_correlation = correlation;
		_var = new double[_dim];
		for (int i=0; i<_dim; i++) {
			for (int j=0; j<_dim; j++)
				_cov[i][j] = scale[i]*correlation[i][j]*Math.sqrt(df[i]/(df[i]-2)*df[j]/(df[j]-2));
			_kurt[i] = 6/(df[i]-4);
			_var[i] = _cov[i][i];
		}
	}
	
	@Override
	public double[] getRealization() {
		double[] z = new double[_dim];
		for (int i=0; i<_dim; i++) {
			double u,v,w;
			do {
				u = 2*_generator.nextDouble()-1;
				v = 2*_generator.nextDouble()-1;
				w = u*u+v*v;
			} while (w>1);
			double c = u*u/w;
			double r = _df[i]*(Math.pow(w,-2./_df[i])-1);
			z[i] = Math.sqrt(r*c)*(_generator.nextBoolean()?1:-1);
		}
		return transform(z);
	}
	
	@Override
	public double[] getRealization(double[] u) {
		double[] z = new double[_dim];
		for (int i=0; i<_dim; i++)
			z[i] = StudentDistQuick.inverseF(_df[i], u[i]);
		return transform(z);
	}
	
	@Override
	public double[] transform(double[] z) {
		double[] x = new double[_chol.length];
		for (int i=0; i<_dim; i++) {
			double sum = 0.0;
			for (int j=0; j<i+1; j++)
				sum += z[j]*_chol[i][j];
			x[i] = _mean[i] + sum*_scale[i];
		}
		return x;
	}
	
	public double[] getScale() {
		return _scale;
	}
	
	public int[] getDf() {
		return _df;
	}
	
	public double[][] getCorrelation() {
		return _correlation;
	}

	@Override
	public double[][] getCov() {
		return _cov;
	}

	@Override
	public double[] getKurtosis() {
		return _kurt;
	}

	@Override
	public DISTRIBUTION getType() {
		return DISTRIBUTION.Student;
	}

} public class MultivariateUniform extends MultivariateNormal {
	
	double[] _min;
	double[] _max;
	double[] _correl2;
	
	public MultivariateUniform(double[] min, double[] max, double[][] correlation, Random generator) {
		super(new double[min.length],adjustCorrelation(correlation),generator);
		_min = min;
		_max = max;
		_mean = new double[_dim];
		_cov = new double[_dim][_dim];
		for (int i=0; i<_dim; i++) {
			_mean[i] = min[i] + (max[i]-min[i])/2;
			_cov[i][i] = 1./12*Math.pow(max[i]-min[i],2);
		}
		for (int i=0; i<_dim; i++) {
			for (int j=i+1; j<_dim; j++)  {
				_cov[i][j] = correlation[i][j]*Math.sqrt(_cov[i][i]*_cov[j][j]);
				_cov[j][i] = _cov[i][j];
			}
			_var[i] = _cov[i][i];
		}
	}
	
	static double[][] adjustCorrelation(double[][] correlation) {
		int dim = correlation.length;
		double[][] correl2 = new double[dim][dim];
		for (int i=0; i<dim; i++)  {
			correl2[i][i] = 1;
			for (int j=i+1; j<dim; j++) {
				correl2[i][j] = 2*Math.sin(Math.PI/6.0*correlation[i][j]);
				correl2[j][i] = correl2[i][j];
			}
		}
		return correl2;
	}
	
	public MultivariateUniform(double[] mean, double[][] covariance, Random generator) {
		super(new double[mean.length],covariance,generator);
		double[][] correl = new double[_dim][_dim];
		_min = new double[mean.length];
		_max = new double[mean.length];
		_mean = new double[_dim];
		for (int i=0; i<mean.length; i++) {
			_min[i] = mean[i] - Math.sqrt(3*covariance[i][i]);
			_max[i] = mean[i] + Math.sqrt(3*covariance[i][i]);
			_mean[i] = _min[i] + (_max[i]-_min[i])/2;
			for (int j=0; j<_dim; j++) 
				correl[i][j] = _cov[i][j]/Math.sqrt(_cov[i][i]*_cov[j][j]);
		}
		_chol = choleskyDecomposition(adjustCorrelation(correl));
	}
	
	@Override
	public double[] transform(double[] z) {
		double[] x = new double[_dim];
		for (int i=0; i<_dim; i++) {
			x[i] = 0.0;
			for (int j=0; j<i+1; j++)
				x[i] += z[j]*_chol[i][j];
			x[i] = NormalDistQuick.cdf01(x[i])*(_max[i]-_min[i])+_min[i];
		}
		return x;
	}
	
	@Override
	public double[] getMean() {
		return _mean;
	}

	@Override
	public double[][] getCov() {
		return _cov;
	}

	@Override
	public double[] getSkewness() {
		return new double[_dim];
	}

	@Override
	public double[] getKurtosis() {
		double[] kurt = new double[_dim];
		Arrays.fill(kurt, -6./5);
		return kurt;
	}

	@Override
	public DISTRIBUTION getType() {
		return DISTRIBUTION.Uniform;
	}

} public enum METHOD {
	MonteCarlo, QuasiMonteCarlo, MomentMatching, QuantizationLearning, QuantizationGrids, Scenred2, VoronoiCellSampling
} public class MomentMatching extends ReductionMethod {
	
	static String _sep = System.getProperty("file.separator");
	static String _mmpath = System.getProperty("user.dir")+"/bin/";
	static String _userdir = System.getProperty("user.dir")+_sep;
	long _time = System.nanoTime();
	int _dim;
	double[] _mean;
	double[] _stdev;
	double[] _skewness;
	double[] _kurtosis;
	double[][] _correl;
	boolean _debugmode = false;
	long _seed;
	
	public MomentMatching(double[] mean, double[][] covariance, long seed) {
		_dim = mean.length;
		if (covariance.length != mean.length || covariance[0].length != mean.length)
	         throw new IllegalArgumentException("Covariance matrix dimensions invalid.");
		_mean = mean;
		_skewness = new double[_dim];
		_kurtosis = new double[_dim];
		_correl = correlation(mean,covariance);
		_stdev = stdev(covariance);
		_seed = seed;
	}
	
	public MomentMatching(double[] mean, double[][] covariance, double[] skewness, double[] kurtosis, long seed) {
		_dim = mean.length;
		if (covariance.length != mean.length || covariance[0].length != mean.length)
	         throw new IllegalArgumentException("Covariance matrix dimensions invalid.");
		_mean = mean;
		_skewness = skewness;
		_kurtosis = kurtosis;
		_correl = correlation(mean,covariance);
		_stdev = stdev(covariance);
		_seed = seed;
	}
	
	static double[][] correlation (double[] mean, double[][] covariance) {
		double[][] correl = new double[covariance.length][covariance.length];
		for (int i=0; i<covariance.length; i++) {
			correl[i][i] = 1.;
			for (int j=i+1; j<covariance.length; j++) {
				correl[i][j] = covariance[i][j]/(Math.sqrt(covariance[i][i])*Math.sqrt(covariance[j][j]));
				correl[j][i] = correl[i][j];
			}
		}	
		return correl;
	}
	
	static double[] stdev (double[][] covariance) {
		double[] stdev = new double[covariance.length];
		for (int i=0; i<covariance.length; i++)
			stdev[i] = Math.sqrt(covariance[i][i]);
		return stdev;
	}
	
	@Override
	public synchronized Map<double[], Double> getScenarios(int numScen) {
		Map<double[],Double> map = null;
		StringBuilder sb = new StringBuilder();
		int seed = new Xorshift(_seed).nextInt();
		try {
			writeMoments(_dim);
			writeCorrelation(_dim);
			String filename = null;
			if(System.getProperty("os.name").startsWith("Mac")){
				filename = "scengen_HKW.mac";
			}
			else if(System.getProperty("os.name").startsWith("Linux")){
				filename = "scengen_HKW.linux";
			}
			else if(System.getProperty("os.name").startsWith("Windows")){
				filename = "scen-gen_HKW.exe";
			}
			else if(System.getProperty("os.name").startsWith("Microsoft")){
				filename = "scen-gen_HKW.exe";
			}
			else {
				throw new IllegalStateException("No compatible binary found for your OS.");
			}
			Process p = Runtime.getRuntime().exec(_mmpath+filename+" "+numScen+" -f 1 -r "+seed);
			BufferedReader reader=new BufferedReader(new InputStreamReader(p.getInputStream())); 
			String line = reader.readLine();
			while(line!=null) { 
				sb.append(line);
				line=reader.readLine(); 
			}
			map = readOutput(numScen,_dim);
		} catch (IOException e) { 
			System.out.println(sb.toString()); 
			e.printStackTrace(); }
		finally {
			new File(_userdir+"out_scen.txt").delete();
			new File(_userdir+"tg_moms.txt").delete();
			new File(_userdir+"tg_corrs.txt").delete();
		}
		return map;
	}
	
	void writeMoments(int dim) throws IOException {
		FileOutputStream stream = new FileOutputStream(_userdir+"tg_moms.txt");
		BufferedWriter br = new BufferedWriter(new OutputStreamWriter(stream));
		br.write("4\n");
		br.write(dim+"\n\n");
		StringBuilder sb = new StringBuilder();
		for (int i=0; i<dim; i++)
			sb.append(_mean[i]+"\t");
		sb.append("\n");
		for (int i=0; i<dim; i++)
			sb.append(_stdev[i]+"\t");
		sb.append("\n");
		for (int i=0; i<dim; i++)
			sb.append(_skewness[i]+"\t");
		sb.append("\n");
		for (int i=0; i<dim; i++)
			sb.append(_kurtosis[i]+"\t");
		br.write(sb.toString());
		br.close();
	}
	
	void writeCorrelation(int dim) throws IOException {
		FileOutputStream stream = new FileOutputStream(_userdir+"tg_corrs.txt");
		BufferedWriter br = new BufferedWriter(new OutputStreamWriter(stream));
		br.write(dim+"\n");
		br.write(dim+"\n");
		StringBuilder sb = new StringBuilder();
		for (int i=0; i<dim; i++) {
			sb.append("\n");
			for (int j=0; j<dim; j++)
				sb.append(_correl[i][j]+"\t");
		}
		br.write(sb.toString());
		br.close();
	}
	
	Map<double[], Double> readOutput(int numScen, int dim) throws IOException  {
		Map<double[],Double> map = new LinkedHashMap<>();
		FileInputStream stream = new FileInputStream(_userdir+"out_scen.txt");
		BufferedReader br = new BufferedReader(new InputStreamReader(stream));
		for (int i=0; i<4; i++) br.readLine();
		String s;
		double[] values = null;
		while ((s = br.readLine())!=null) {
			String[] sar = s.split("\\s+");
//			System.out.println(Arrays.toString(sar));
			values = new double[dim];
			for (int i=0; i<dim; i++)
				values[i] = Double.parseDouble(sar[i+1]);
			map.put(values,Double.parseDouble(sar[0]));
		}
		br.close();
		return map;
	}
	
	

} public class MonteCarlo extends ReductionMethod {
	
	MultivariateDistribution _mvDist;
	
	public MonteCarlo (MultivariateDistribution mvDist) {
		_mvDist = mvDist;
		
	}
	
	@Override
	public Map<double[], Double> getScenarios(int numScen) {
		Map<double[],Double> resultMap = new LinkedHashMap<>();
		for (int i=0; i<numScen; i++)
			resultMap.put(_mvDist.getRealization(),1./numScen);
		return resultMap;
	}
} public class QuantizationGrids extends ReductionMethod {

	MultivariateNormal _normal;
	int _dim;
	
	public QuantizationGrids (MultivariateDistribution mvDist) {
		if (mvDist.getDim()>10 || mvDist.getType()!=DISTRIBUTION.Normal)
			throw new IllegalArgumentException("VectorQuantization1 accepts only multivariate normal distributions with at most 10 dimensions.");
		_dim = mvDist.getDim();
		_normal = (MultivariateNormal) mvDist;
	}
	
	public Map<double[],Double> getScenarios(int numScen) {
		Map<double[],Double> map = new LinkedHashMap<>();
		try {
			String filename = String.format("%d_%d_%s",numScen,_dim,_dim>1?"nopti":"dualopti");
			InputStream stream = ClassLoader.getSystemClassLoader().getResourceAsStream(filename);
			BufferedReader br = new BufferedReader(new InputStreamReader(stream));
			String s;
			double[] values = null;
			int shift = _dim>1?1:0;
			while ((s = br.readLine())!=null) {
				String[] sar = s.split("\\s+");
				values = new double[_dim];
				for (int i=0; i<_dim; i++)
					values[i] = Double.parseDouble(sar[i+1+shift]);
				map.put(values,Double.parseDouble(sar[shift]));
			}
			//remove last line
			map.remove(values);
			br.close();
		} catch (IOException e) {e.printStackTrace();}
		return fitMultivariate(map);
	}
	
	public Map<double[],Double> fitMultivariate(Map<double[],Double> scen) {
		double[] mean = _normal.getMean();
		double[][] chol = _normal.getCholesky();
		Map<double[],Double> map = new LinkedHashMap<>();
		for (double[] d: scen.keySet()) {
			double[] x = new double[_dim];
			for (int i=0; i<_dim; i++) {
				double sum = 0.0;
				for (int j=0; j<i+1; j++)
					sum += d[j]*chol[i][j];
				x[i] = mean[i] + sum;
			}
			map.put(x,scen.get(d));
		}
		return map;
	}
	
	
	
	
} public class QuantizationLearning extends ReductionMethod {
	
	MultivariateDistribution _mvdist;
	
	public QuantizationLearning (MultivariateDistribution dist) {
		_mvdist = dist;
	}
	
	@Override
	public Map<double[], Double> getScenarios(int numScen) {
		int dim = _mvdist.getDim();
		double a = 100*numScen;
		int sampleSize = 10000*numScen; 
		double[][] scenarios = new double[numScen][];
		double[] probs = new double[numScen];
		Arrays.fill(probs, 1);
		MonteCarlo rand = new MonteCarlo(_mvdist);
		Iterator<double[]> iter = rand.getScenarios(numScen).keySet().iterator();
		for (int j=0; j<numScen; j++) 
			scenarios[j] = iter.next();
		for (int i=0; i<sampleSize; i++) {
			double[] point = _mvdist.getRealization();
			int nearest = 0;
			double minDist = Double.MAX_VALUE;
			for (int j=0; j<numScen; j++) {
				double dist = Metrics.squaredDistance(scenarios[j], point);
				if (dist<minDist) {
					minDist = dist;
					nearest = j;
				}
			}
			for (int d=0; d<dim; d++) 
				scenarios[nearest][d] += a/(a+i+1)*(point[d]-scenarios[nearest][d]);
			for (int j=0; j<numScen; j++) {
				if (j==nearest)
					probs[j] += (1.-probs[j])/(i+1.);
				else
					probs[j] -= probs[j]/(i+1.);
			}
    	}
		Map<double[],Double> map = new LinkedHashMap<>();
		for (int j=0; j<numScen; j++) 
			map.put(scenarios[j],probs[j]);
		return map;
	}

} public class QuasiMonteCarlo extends ReductionMethod {
	
	MultivariateDistribution _mvDist;
	static RandomStream rand = new MRG32k3a();
	
	public QuasiMonteCarlo (MultivariateDistribution mvDist) {
		_mvDist = mvDist;		
	}
	
	@Override
	public Map<double[], Double> getScenarios(int numScen) {
		int dim = _mvDist.getDim();
		SobolSequence sobol = new SobolSequence(numScen+1,dim);
		sobol.leftMatrixScramble(rand);
		PointSetIterator iter = sobol.iterator();
		iter.resetToNextPoint();
		Map<double[],Double> resultMap = new LinkedHashMap<>();
		for (int i=0; i<numScen; i++) {
			double[] u = new double[dim];
			iter.nextPoint(u, dim);
			resultMap.put(_mvDist.getRealization(u),1./numScen);
		}
		return resultMap;
	}
} public abstract class ReductionMethod {
	
	public abstract Map<double[],Double> getScenarios(int numScen);
	
	/**
	 * @param dim
	 * @param scenarios
	 * @return Returns the scenario along one particular dimension.
	 */
	public static double[] getScenarios(int dim, Map<double[],Double> scenarios) {
		double[] scen = new double[scenarios.size()];
		int i=0;
		for (double[] x : scenarios.keySet())
			scen[i++] = x[dim];
		return scen;
	}
	
	public static double[] getWeights(Map<double[],Double> scenarios) {
		double[] weights = new double[scenarios.size()];
		int i=0;
		for (double[] x : scenarios.keySet())
			weights[i++] = scenarios.get(x);
		return weights;
	}

} public class Scenred2 extends ReductionMethod {
	
	static String _sep = System.getProperty("file.separator");
	static String _scenredpath = System.getProperty("user.dir")+"/bin/";
	static String _userdir = System.getProperty("user.dir")+_sep;
	
	static int _timeLimit = 600; //seconds
	static int _metricType = 3; //1:transport, 2:fortet-mourier ,3:wasserstein
	static int _reductionMethod = 1; //1:forward, 2:backward
	static int _pNorm = 2; //0: max, 1:sum, 2:euclidean
	static int _order = 2; //metric order 
	static int _dim;
	static int _size;
	Map<double[],Double> _input;
	
	public Scenred2(Map<double[],Double> input, int dim) {
		_dim = dim;
		_input = input;
	}
	
	@Override
	public synchronized Map<double[], Double> getScenarios(int numScen) {
		Map<double[],Double> map = null;
		StringBuilder sb = new StringBuilder();
		try {
			writeCommand();
			writeParam(numScen);
			writeSample();
			Process p = null;
			if(System.getProperty("os.name").startsWith("Windows"))
				p = Runtime.getRuntime().exec(_scenredpath+"scenred2.exe command.txt -nogams");
			else {
				throw new IllegalStateException("No compatible binary found for your OS.");
			}
			BufferedReader reader=new BufferedReader(new InputStreamReader(p.getInputStream())); 
			String line = reader.readLine();
			while(line!=null) { 
				sb.append(line);
				line=reader.readLine(); 
			}
			map = readOutput(numScen,_dim);
		} catch (IOException e) { 
			System.out.println(sb.toString()); 
			e.printStackTrace(); 
		}
		finally {
			new File(_userdir+"scenred2_in.txt").delete();
			new File(_userdir+"scenred2_out.txt").delete();
			new File(_userdir+"scenred2_param.txt").delete();
		}
		return map;
	}
	
	void writeSample() throws IOException {
		FileOutputStream stream = new FileOutputStream(_userdir+"scenred2_in.txt");
		BufferedWriter br = new BufferedWriter(new OutputStreamWriter(stream));
		br.write(String.format("TYPE FAN\n"
				+ "TIME 2\n"
				+ "SCEN %d\n"
				+ "RANDOM %d\n"
				+ "DATA\n",_input.size(),_dim));
		StringBuilder sb = new StringBuilder();
		for (double[] d : _input.keySet()) {
			sb.append(_input.get(d)+"\n");
			for (int j=0; j<_dim; j++)
				sb.append(1.0+"\t");
			sb.append("\n");
			for (int j=0; j<_dim; j++)
				sb.append(d[j]+"\t");
			sb.append("\n\n");
		}
		sb.append("END");
		br.write(sb.toString());
		br.close();
	}
	
	void writeCommand() throws IOException {
		FileOutputStream stream = new FileOutputStream(_userdir+"command.txt");
		BufferedWriter br = new BufferedWriter(new OutputStreamWriter(stream));
		br.write(String.format("report_level 1\n"
				+ "runtime_limit %d\n"
				+ "read_scen %sscenred2_in.txt\n"
				+ "scen_red %sscenred2_param.txt\n"
				+ "out_scen %sscenred2_out.txt",_timeLimit,_userdir,_userdir,_userdir));
		br.close();
	}
	
	void writeParam(int numScen) throws IOException {
		FileOutputStream stream = new FileOutputStream(_userdir+"scenred2_param.txt");
		BufferedWriter br = new BufferedWriter(new OutputStreamWriter(stream));
		br.write(String.format("red_num_leaves %d\n"
				+ "metric_type %d\n"
				+ "p_norm %d\n"
				+ "reduction_method %d\n"
				+ "order %d"
				+ "\n",numScen,_metricType,_pNorm,_reductionMethod,_order));
		br.close();
	}
	
	Map<double[], Double> readOutput(int numScen, int dim) throws IOException  {
		Map<double[],Double> map = new LinkedHashMap<>();
		FileInputStream stream = new FileInputStream(_userdir+"scenred2_out.txt");
		BufferedReader br = new BufferedReader(new InputStreamReader(stream));
//		for (int i=0; i<9; i++) br.readLine();
		String s;
		double[] values = null;
		while (!(s = br.readLine()).matches("DATA"));
		while ((s = br.readLine())!=null) {
			double prob = Double.parseDouble(s);
			s = br.readLine();
			s = br.readLine();
			String[] sar = s.split("\\s+");
			values = new double[dim];
			for (int i=0; i<dim; i++)
				values[i] = Double.parseDouble(sar[i]);
			map.put(values,prob);
			s = br.readLine();
			if (s.matches("END"))
				break;
		}
		br.close();
		return map;
	}
	
	

} public class VoronoiCellSampling extends QuantizationLearning {
	
	public VoronoiCellSampling(MultivariateDistribution dist) {
		super(dist);

	}
	
	public Map<double[], Double> getScenarios(int numScen) {
		int dim = _mvdist.getDim();
		double a = 100*numScen;
		int sampleSize = 10000*numScen; 
		double[][] scenarios1 = new double[numScen][];
		double[][] scenarios2 = new double[numScen][];
		double[] probs = new double[numScen];
		Arrays.fill(probs, 1);
		MonteCarlo rand = new MonteCarlo(_mvdist);
		Iterator<double[]> iter = rand.getScenarios(numScen).keySet().iterator();
		for (int j=0; j<numScen; j++) {
			double[] x = iter.next();
			scenarios1[j] = x;
			scenarios2[j] = x;
		}
		for (int i=0; i<sampleSize; i++) {
			double[] point = _mvdist.getRealization();
			int nearest = 0;
			double minDist = Double.MAX_VALUE;
			for (int j=0; j<numScen; j++) {
				double dist = Metrics.squaredDistance(scenarios1[j], point);
				if (dist<minDist) {
					minDist = dist;
					nearest = j;
				}
			}
			scenarios2[nearest] = point;
			for (int d=0; d<dim; d++) 
				scenarios1[nearest][d] += a/(a+i+1)*(point[d]-scenarios1[nearest][d]);
			for (int j=0; j<numScen; j++) {
				if (j==nearest)
					probs[j] += (1.-probs[j])/(i+1.);
				else
					probs[j] -= probs[j]/(i+1.);
			}
    	}
		Map<double[],Double> map = new LinkedHashMap<>();
		for (int j=0; j<numScen; j++) {
			map.put(scenarios2[j],probs[j]);
		}
		return map;
	}
	

}
import java.util.Map;
import java.util.Random;
import scengen.distribution.DISTRIBUTION;
import scengen.distribution.MultivariateDistribution;
import scengen.distribution.MultivariateLognormal;
import scengen.distribution.MultivariateNormal;
import scengen.distribution.MultivariateStudent;
import scengen.distribution.MultivariateUniform;
import scengen.methods.METHOD;
import scengen.methods.MomentMatching;
import scengen.methods.MonteCarlo;
import scengen.methods.QuantizationGrids;
import scengen.methods.QuantizationLearning;
import scengen.methods.QuasiMonteCarlo;
import scengen.methods.ReductionMethod;
import scengen.methods.Scenred2;
import scengen.methods.VoronoiCellSampling;
package scengen.paper;




public class RunGenerator {

	public static void main(String... args) {
		int dim = 2;
		int numScen = 10;
		print(makeScen(dim,numScen,DISTRIBUTION.Normal,METHOD.MomentMatching,new Random()),dim);
	}

	public static Map<double[],Double> makeScen(int dim, int numScen, DISTRIBUTION dist, METHOD method, Random rand) {
		double[] mean = new double[dim];
		double[] scale = new double[dim];
		int[] df = new int[dim];
		double[][] cov = new double[dim][dim];
		for (int i=0; i<dim; i++) {
			mean[i] = dist==DISTRIBUTION.Lognormal?1:0;
			scale[i] = 1;
			df[i] = 5;
			cov[i][i] = 0.5;
		}
		MultivariateDistribution mvdist = null;
		ReductionMethod pb = null;
		switch (dist) {
		case Normal: mvdist = new MultivariateNormal(mean,cov,rand); 
			break;
		case Lognormal: mvdist = new MultivariateLognormal(mean,cov,rand);
			break;
		case Uniform: mvdist = new MultivariateUniform(mean,cov,rand);
			break;
		case Student: mvdist = new MultivariateStudent(mean,scale,cov,df,rand);