public class CSVDataJava extends ScaledData {
	/**
	 * Parse given csvfile into set of input and target values.
	 * 
	 * @param csvfile file of comma-separated values, last value in each line is
	 *        the target value
	 */
	public CSVDataJava(String csvfile) {
		super(DataSizeRetreiver.num_fitness_cases(csvfile), DataSizeRetreiver.num_terminals(csvfile));
		BufferedReader f;
		try {
                    f = new BufferedReader(new FileReader(csvfile));
                    String[] token;
                    try {
                            /**
                             * Keep track of max/min target values (Vladislavleva suggested) to
                             * perform approximate scaling.
                             */
                            int fitnessCaseIndex = 0;
                            while (f.ready() && fitnessCaseIndex < numberOfFitnessCases) {
                                    token = f.readLine().split(",");
                                    for (int i = 0; i < token.length - 1; i++) {
                                            this.fitnessCases[fitnessCaseIndex][i] = Double.valueOf(token[i]);
                                            if(fitnessCases[fitnessCaseIndex][i] < minFeatures[i]) minFeatures[i] = fitnessCases[fitnessCaseIndex][i];
                                            if(fitnessCases[fitnessCaseIndex][i] > maxFeatures[i]) maxFeatures[i] = fitnessCases[fitnessCaseIndex][i];
                                    }
                                    Double val = Double.valueOf(token[token.length - 1]);
                                    addTargetValue(val, fitnessCaseIndex);
                                    fitnessCaseIndex++;
                            }
                            this.scaleTarget();
                    } catch (NumberFormatException e) {
                    } catch (IOException e) {
                    }
                } catch (FileNotFoundException e) {
                    System.err.print("File not found Exception in CSVDataJava constructor");
                }
	}
	
        /**
         * Display length/width of data
         */
	public void printDataInfo() {
		System.out.println("We have " + this.numberOfFitnessCases + " fitness cases and " + this.numberOfFeatures+ " values");
	}
} public interface DataJava {
	/**
	 * Get the current set of input data, where each row of input corresponds to
	 * a list of values in the return.
	 * 
	 * @return Two Dimensional Array; equivalent to the matrix of
	 *         input data, where each training case corresponds to a row.
	 */
	public double[][] getInputValues();

	/**
	 * Get the target values corresponding to the training cases. The ith targt
	 * value is the output corresponding to the ith row of the input values (the
	 * ith training case).
	 * 
	 * @return One Dimensional Array; equivalent to the column of output
	 *         values, where each value is the output for the given row of
	 *         input.
	 */
	public double[] getTargetValues();

	/**
	 * Returns the precomputed mean on the target values. Used in the ModelScaler
         * @return mean target value
	 */
	public Double getTargetMean();
	
	/**
	 * Return the maximum of all the target values.
	 * @return
	 */
	public Double getTargetMax();

	/**
	 * Return the minimum of all target values.
	 * @return
	 */
	public Double getTargetMin();

	/**
	 * Get the target values, scaled to be in the range [0,1]. This is done
	 * according to the following formula.
	 * <p>
	 * <code>y_i_scaled = (y_i-getTargetMin())/(getTargetMax()-getTargetMin())</code>
	 * <p>
	 * Otherwise, this has the same properties as {@link #getTargetValues()}.
	 * 
	 * @return
	 */
	public double[] getScaledTargetValues();
        
        /**
         * @return the number Of Fitness Cases
         */
        public int getNumberOfFitnessCases();

        /**
         * @return the number Of Features
         */
        public int getNumberOfFeatures();
} public abstract class FitnessFunction {

	// used to access this fitness function's position in the
	// gp.Individual.fitnesses HashMap. Each subclass must
	// define this value
	//public static String FITNESS_KEY;

	// should this value be maximized (true) or minimized (false)?
	// used by non-dominated sort to convert all fitness functions
	// to minimizers
	// public Boolean isMaximizingFunction;

	/**
	 * Evaluate a single individual
	 * 
	 * @param pop
	 */

	/**
	 * Evaluate each individual in a population
	 * 
	 * @param pop
	 */
	public abstract void evalPop(Population pop);

        /**
         * Determine whether the fitness function needs to be minimized or maximized
         * @return flag for minimization or maximization
         */
	public abstract Boolean isMaximizingFunction();

} public class SubtreeComplexityFitness extends FitnessFunction {

        /**
         * Fitness key
         */
	public static final String FITNESS_KEY = Parameters.Operators.SUBTREE_COMPLEXITY_FITNESS;

        /**
         * whether to minimize/maximize function
         */
	public Boolean isMaximizingFunction = false;

        /**
         * this function provides discrete fitness values (integer)
         */
	public Boolean discreteFitness = true;

	/**
	 * Should this fitness function be minimized (i.e. mean squared error) or
	 * maximized?
         * @return 
	 */
        @Override
	public Boolean isMaximizingFunction() {
		return this.isMaximizingFunction;
	}

	/**
	 * Compute the subtree complexity of the individual
         * @param ind
	 */
	public void eval(Individual ind) {
		Tree t = (Tree) ind.getGenotype();
		Integer complexity = t.getSubtreeComplexity();
		ind.setFitness(SubtreeComplexityFitness.FITNESS_KEY, (double) complexity);
	}

	@Override
	public void evalPop(Population pop) {
		for (Individual individual : pop) {
			this.eval(individual);
		}
	}
} public class SRLARSJava extends FitnessFunction {

    private final DataJava data;
    private int pow;
    private final boolean USE_INT;
    
    /**
     * Fitness key
     */
    public static String FITNESS_KEY = Parameters.Operators.SR_JAVA_FITNESS;
    
    /**
     * Maximizing function
     */
    public Boolean isMaximizingFunction = true;
    private int numThreads;
    /**
     * Create a new fitness operator, using the provided data, for assessing
     * individual solutions to Symbolic Regression problems. There are two
     * parameters for this fitness evaluation:
     * <p>
     * <ul>
     * <li>The power <i>p</i> to use in computing the mean of the errors. See
     * {@link #getMeanFromP(int)} for recognized values and more information.
     * This parameter is set to the value specified by the key
     * {@value algorithm.Parameters.Names#MEAN_POW} in the properties file. The
     * default value is {@value algorithm.Parameters.Defaults#MEAN_POW}.
     * <li>A boolean flag specifying if the predictions should be forced to
     * integers because the output variable is integer-valued (ie you might want
     * to set this to true if your output is integer-valued). This parameter is
     * set to the boolean at the key
     * {@value algorithm.Parameters.Names#COERCE_TO_INT}. If there is no value
     * specified, the default is
     * {@value algorithm.Parameters.Defaults#COERCE_TO_INT}.
     * </ul>
     * 
     * @param data
     *            The dataset (training cases, output variable) to use in
     *            computing the fitness of individuals.
     */
    public SRLARSJava(DataJava data) {
        this(data, 2, false,1);
    }
    
    /**
     * Constructor that allows to specify the number of threads for parallel evaluation
     * @param aData
     * @param aPow
     * @param is_int
     * @param anumThreads 
     */
    public SRLARSJava(DataJava aData, int aPow, boolean is_int,int anumThreads) {
        this.data = aData;
        pow = aPow;
        USE_INT = is_int;
        numThreads = anumThreads;
    }
    /**
     * Should this fitness function be minimized (i.e. mean squared error) or
     * maximized?
     * @return 
     */
    @Override
    public Boolean isMaximizingFunction() {
        return this.isMaximizingFunction;
    }

    /**
     * Simple "factory"-like method for returning the proper generalized mean
     * object to use, given the parameter p (see
     * http://en.wikipedia.org/wiki/Generalized_mean). Since we are only
     * concerned with means of errors, we don't utilize all the valid values of
     * p defined for generalized means.
     * 
     * @param p
     *            power to use in computing a mean
     * @return an instance of {@link PowerMean} if p > 1; an instance of
     *         {@link ArithmeticMean} if p == 1; otherwise, an instance of
     *         {@link Maximum}.
     */
    public static Mean getMeanFromP(int p) {
        if (p == 1) {
            return new ArithmeticMean();
        } else if (p > 1) {
            return new PowerMean(p);
        } else {
            return new Maximum();
        }
    }

    /**
     * @param ind
     * @throws java.lang.Exception
     * @see Function
     */
    public void eval(Individual ind) throws Exception {
        
        Tree genotype = (Tree) ind.getGenotype();
        Mean MEAN_FUNC = getMeanFromP(pow);
        Function func = genotype.generate();
        List<Double> d;
        ArrayList<Double> interVals;
        double[][] inputValuesAux = data.getInputValues();
        Tree tAux = (Tree) ind.getGenotype();
        //double[] targetAux = data.getScaledTargetValues();
        double[] targetAux = data.getTargetValues();
        float[][] intermediateValues = new float[data.getNumberOfFitnessCases()][tAux.getSize()];
        for (int i = 0; i < data.getNumberOfFitnessCases(); i++) {
            d = new ArrayList<Double>();
            for (int j = 0; j < data.getNumberOfFeatures(); j++) {
                d.add(j, inputValuesAux[i][j]);
            }
            interVals = new ArrayList<Double>();
            Double val = func.evalIntermediate(d,interVals);
            for(int t=0;t<interVals.size();t++){
                intermediateValues[i][t] = interVals.get(t).floatValue();
            }
            d.clear();
            interVals.clear();
        }

        /*
         * LassoFitGenerator is initialized
         */
        LassoFitGenerator fitGenerator = new LassoFitGenerator();
        int numObservations = data.getNumberOfFitnessCases();
        fitGenerator.init(tAux.getSize(), numObservations);
        for (int i = 0; i < numObservations; i++) {
            fitGenerator.setObservationValues(i,intermediateValues[i] );
            fitGenerator.setTarget(i, targetAux[i]);
        }

        /*
         * Generate the Lasso fit. The -1 arguments means that
         * there would be no limit on the maximum number of 
         * features per model
         */
        LassoFit fit = fitGenerator.fit(-1);

        
        // We pick the first value of lamda that includes all the features in the model
        int indexWeights = 0;
        int usedVars=0;
        for(int i=0;i<fit.lambdas.length;i++){
            if(fit.nonZeroWeights[i]>usedVars){
                indexWeights = i;
                usedVars = fit.nonZeroWeights[i];
            }
        }
        ArrayList<String> alWeights = new ArrayList<String>();
        double lassoIntercept = 0;
        double fitness = 0;
        //if(fit.nonZeroWeights[indexWeights]==tAux.getSize()){
            //double[] lassoWeights = fit.compressedWeights[indexWeights];
            double[] lassoWeights = fit.getWeights(indexWeights);
            for(int j=0;j<lassoWeights.length;j++){
                alWeights.add(Double.toString(lassoWeights[j]));
            }
            lassoIntercept = fit.intercepts[indexWeights];
            //SRPhenotype phenotype_tmp = new SRPhenotype();
            for (int i = 0; i < data.getTargetValues().length; i++) {
                double prediction = 0;
                for(int j=0;j<lassoWeights.length;j++){
                    prediction += intermediateValues[i][j]*lassoWeights[j];
                }
                prediction += lassoIntercept;
                //phenotype_tmp.addNewDataValue(prediction);
                if (this.USE_INT) {
                    prediction = Math.round(prediction);
                }
                MEAN_FUNC.addValue(Math.abs(targetAux[i] - prediction));
            }
            /*final ScaledTransform scalePrediction = new ScaledTransform(phenotype_tmp.min, phenotype_tmp.max);
            final ScaledTransform scaleTarget = new ScaledTransform(this.data.getTargetMin(), this.data.getTargetMax());
            for (int i = 0; i < this.getTarget().length; i++) {
                Double scaled_val = scalePrediction.scaleValue(phenotype_tmp.getDataValue(i));
                if (this.USE_INT) {
                    int rounded_val = scaleTarget.unScaleValue(scaled_val).intValue();
                    scaled_val = scaleTarget.scaleValue((double) rounded_val);
                }
                MEAN_FUNC.addValue(Math.abs(targetAux[i] - scaled_val));
            }
            phenotype_tmp = null;
            */
            Double error = MEAN_FUNC.getMean();
            fitness=errorToFitness(error);
            
        /*}else{
            for(int j=0;j<tAux.getSize();j++){
                alWeights.add("0");
            }
        }*/
        ind.setWeights(alWeights);
        ind.setLassoIntercept(Double.toString(lassoIntercept));
        ind.setFitness(SRLARSJava.FITNESS_KEY, fitness);
}
    
    /**
     * Transform errors to fitness values. For errors, smaller values are
     * better, while for fitness, values closer to 1 are better. This particular
     * transformation also places a greater emphasis on small changes close to
     * the solution (fitness == 1.0 represents a complete solution). However,
     * this transformation also assumes that the error is in the range [0,1].
     * 
     * @param error
     *            Error on training set; value in range [0,1].
     * @return 0 if error is not a number (NaN) or outside of the range [0,1].
     *         Otherwise, return the fitness.
     */
    private Double errorToFitness(Double error) {
        if(error==0){
            return 1.0;
        }else if (error.isNaN() || error < 0.0 ) {
            return 0.0;
        } else {
            return (1.0) / error;
        }
    }

    @Override
    public void evalPop(Population pop) {
        
        ArrayList<SRJavaThread> alThreads = new ArrayList<SRJavaThread>();
        for(int i=0;i<numThreads;i++){
            SRJavaThread threadAux = new SRJavaThread(i, pop,numThreads);
            alThreads.add(threadAux);
        }
        
        for(int i=0;i<numThreads;i++){
            SRJavaThread threadAux = alThreads.get(i);
            threadAux.start();
        }
        
        for(int i=0;i<numThreads;i++){
            SRJavaThread threadAux = alThreads.get(i);
            try {
                threadAux.join();
            } catch (InterruptedException ex) {
                Logger.getLogger(SRLARSJava.class.getName()).log(Level.SEVERE, null, ex);
            }
        }
        

    }

    /**
     * @return the fitnessCases
     */
    public double[][] getFitnessCases() {
        return data.getInputValues();
    }

    /**
     * @return the scaled_target
     */
    public double[] getTarget() {
        return data.getScaledTargetValues();
    }

    /**
     * Utility class for computing the element-wise normalization of a vector.
     * That is, for a given vector <code>A</code>, where all values are in the
     * range <code>[min(A), max(A)]</code>, we want an easy way to compute the
     * scaled vector <code>B</code> such that every value in B is in the range
     * <code>[0, 1]</code>.
     * <p>
     * In particular, if the ith element of <code>A</code> has value
     * <code>a_i</code>, then the ith element of <code>B</code> will have the
     * value <code>(a_i-min(A))/(max(A)-min(A))</code>.
     * <p>
     * Also, provide the inverse transformation, to recover the values of
     * <code>A</code>.
     * 
     * @author Owen Derby
     * @see ScaledData
     */
    private class ScaledTransform {
        private final double min, range;

        public ScaledTransform(double min, double max) {
            this.min = min;
            this.range = max - min;
        }

        public Double scaleValue(Double val) {
            return (val - min) / range;
        }

        public Double unScaleValue(Double scaled_val) {
            return (scaled_val * range) + min;
        }
    }
    
    /**
     * Auxiliary thread class to perform parallel evaluation of the population
     * Each thread evaluates a subset of the populaiton
     */
    public class SRJavaThread extends Thread{
        private int indexThread, totalThreads;
        private Population pop;
        
        /**
         * Constructor of the thread class
         * @param anIndex
         * @param aPop
         * @param aTotalThreads 
         */
        public SRJavaThread(int anIndex, Population aPop,int aTotalThreads){
            indexThread = anIndex;
            pop = aPop;
            totalThreads = aTotalThreads;
        }

        
        /**
         * Each thread evaluates a subset of the population.
         * The subset is selected with the indices of the thread and that of individuals
         */
        @Override
        public void run(){
            int indexIndi = 0;
            for (Individual individual : pop) {
                if(indexIndi%totalThreads==indexThread){
                    try {
                        eval(individual);
                    } catch (Exception ex) {
                        Logger.getLogger(SRLARSJava.class.getName()).log(Level.SEVERE, null, ex);
                    }
                }
                indexIndi++;
            }
        }
     }
        
} public class TreeGenerator {

    private MersenneTwisterFast rng;
    private final List<String> terminals;
    private final List<String> functions;

    /**
     * Create the factory, setting the terminal and function sets to use. If no
     * function sets are provided, the defaults are used from
     * {@link algorithm.Parameters.Defaults}.
     * 
     * @param r An instance of {@link MersenneTwisterFast} to use when selecting
     *        labels for nodes in trees.
     * @param funcset The set of functions to select internal node labels from.
     * @param termset The set of terminals to select terminal node labels from.
     */
    public TreeGenerator(MersenneTwisterFast r, List<String> funcset,List<String> termset) {
        this.rng = r;
        if (funcset == null) {
            System.out.println("No function set provided - using defaults!");
            functions = Parameters.Defaults.FUNCTIONS;
        } else {
            functions = funcset;
        }
        if (termset == null) {
            System.out.println("No terminal set provided - using defaults!");
            terminals = Parameters.Defaults.TERMINALS;
        } else {
            terminals = termset;
        }
    }

    /**
     * Generate an expression tree coding for the sum of the variables of the problem
     * @param termset
     * @return tree coding the sum of all the variables of the problem
     */
    public Tree generateLinearModel(List<String> termset) {
        ArrayList<String> featureList = new ArrayList<String>();
        featureList.addAll(termset);
        TreeNode holder = new TreeNode(null, "holder");
        TreeNode root = new TreeNode(holder, "null");
        generateLinearModel(root, featureList);
        holder.addChild(root);
        Tree t = new Tree(holder);
        return t;
    }
    
    /**
     * Generate a linear model only with the variables in the list
     * @param node
     * @param featureList 
     */
    public void generateLinearModel(TreeNode node,List<String> featureList) {
        // See node as a placeholder: we overwrite its label,
        // children, but we don't change its parent or id.
        node.children.clear();
        if (featureList.size() == 1) {
            // can't go deeper, make it a terminal
            node.label = featureList.get(0);
        } else {
            node.label = "+";
            TreeNode leftNode = new TreeNode(node, "null");
            int indexFeature = rng.nextInt(featureList.size());
            leftNode.label = featureList.get(indexFeature);
            featureList.remove(indexFeature);
            node.addChild(leftNode);
            
            TreeNode rightNode = new TreeNode(node, "null");
            generateLinearModel(rightNode, featureList);
            node.addChild(rightNode);
        }
        node.resetAbove();
    }
    
    /**
     * Generate a tree by either the grow or the fill method (as described by
     * Koza). If using the fill method, the new tree will be a full depth with
     * every terminal at maxDepth. If using the grow method, then the tree is
     * "grown" by randomly choosing either a terminal or function for each node,
     * and recursing for each function (non-terminal) selected. This continues
     * until all terminals are selected or maxDepth is reached.
     * 
     * @param maxDepth the maximum depth the new tree can reach.
     * @param full Boolean value indicating to use fill or grow method. If true,
     *        the fill method will be used. If false, the grow method.
     * @return the newly generated instance of the tree.
     */
    public Tree generateTree(int maxDepth, boolean full) {
        TreeNode holder = new TreeNode(null, "holder");
        TreeNode root = new TreeNode(holder, "null");
        generate(root, maxDepth, full);
        holder.addChild(root);
        Tree t = new Tree(holder);
        return t;
    }

    /**
     * Constructor a tree from a string. The string is assumed to be a LISP-like
     * S-expression, using prefix notation for operators.
     * 
     * @param input The S-expression string.
     * @return new tree representing the inpu.
     */
    public static Tree generateTree(String input) {
        // Make sure the string is tokenizable
        // FIXME allow other delimiters?
        input = input.replace("(", " ( ");
        input = input.replace("[", " [ ");
        input = input.replace(")", " ) ");
        input = input.replace("]", " ] ");

        StringTokenizer st = new StringTokenizer(input);
        TreeNode holder = new TreeNode(null, "holder");
        parseString(holder, st);
        return new Tree(holder);
    }
    
    /**
     * Constructor a tree from a string. The string is assumed to be a LISP-like
     * S-expression, using prefix notation for operators.
     * 
     * @param input The S-expression string.
     * @return new tree representing the inpu.
     */
    public static Tree generateTreeConstants(String input) {
        // Make sure the string is tokenizable
        // FIXME allow other delimiters?
        input = input.replace("(", " ( ");
        input = input.replace("[", " [ ");
        input = input.replace(")", " ) ");
        input = input.replace("]", " ] ");

        StringTokenizer st = new StringTokenizer(input);
        TreeNode holder = new TreeNode(null, "holder");
        parseString(holder, st);
        return new Tree(holder);
    }

    /**
     * Given a node in a tree, generate an entire new subtree below, using
     * either the fill or grow method (see {@link #generateTree(int, boolean)}).
     * Note that the given node will have it's label overwritten and all it's
     * current children cleared.
     * 
     * @param node TreeNode to use as the root of the new subtree.
     * @param depth The allowed depth of the generated subtree.
     * @param full Whether to use the fill or grow method for generating the
     *        subtree.
     */
    public void generate(TreeNode node, int depth, boolean full) {
        // See node as a placeholder: we overwrite its label,
        // children, but we don't change its parent or id.
        node.children.clear();
        if (depth <= 0) {
            // can't go deeper, make it a terminal
            node.label = terminals.get(rng.nextInt(terminals.size()));
        } else {
            int label_idx;
            if (full) {
                // want full tree, so make it a function
                label_idx = rng.nextInt(functions.size());
            } else {
                // growing tree, randomly choose function or terminal
                label_idx = rng.nextInt(functions.size() + terminals.size());
            }
            if (label_idx < functions.size()) {
                node.label = functions.get(label_idx);
                // create and recurse on each argument (child) of function
                for (int i = 0; i < arity(node.label); i++) {
                    TreeNode newNode = new TreeNode(node, "null");
                    node.addChild(newNode);
                    generate(newNode, depth - 1, full);
                }
            } else {
                node.label = terminals.get(label_idx - functions.size());
            }
        }
        node.resetAbove();
    }
    
    /**
     * Given a node in a tree, generate an entire new subtree below, using
     * either the fill or grow method (see {@link #generateTree(int, boolean)}).
     * Note that the given node will have it's label overwritten and all it's
     * current children cleared.
     * 
     * @param node TreeNode to use as the root of the new subtree.
     * @param depth The allowed depth of the generated subtree.
     * @param full Whether to use the fill or grow method for generating the
     *        subtree.
     * @throws evogpj.gp.GPException
     */
    public void generateForMutationWithConstants(TreeNode node, int depth, boolean full) throws GPException {
        // See node as a placeholder: we overwrite its label,
        // children, but we don't change its parent or id.
        
        if (node.children.isEmpty()) {
            node.children.clear();
            if(rng.nextBoolean()){// swap variable
                node.label = terminals.get(rng.nextInt(terminals.size()));
                node.resetCoeff();
            }else{ // change constant
                if(rng.nextBoolean()){// increase coefficient a 10%
                    node.increaseCoeff();
                }else{// decrease coefficient a 10%
                    node.decreaseCoeff();
                }
                
            }
        } else {
            generate(node, depth, full);
        }
        node.resetAbove();
    }

	/**
	 * Compute the arity of the function represented by the string.
	 * 
	 * @param label string encoding a function (ie "sin" or "*")
	 * @return the arity of the function (number of arguments the function
	 *         takes)
	 * @see math.Function#getArityFromLabel(String)
	 */
	public static int arity(String label) {
		// Boolean problems functions:
		if (label.equals("if")) {
			return 3;
		} else if (label.equals("and")) {
			return 2;
		} else if (label.equals("or")) {
			return 2;
		} else if (label.equals("nand")) {
			return 2;
		} else if (label.equals("nor")) {
			return 2;
		} else if (label.equals("not")) {
			return 1;
		}
		// Symbolic Regression functions
		return Function.getArityFromLabel(label);
	}

    private static void parseString(TreeNode parent, StringTokenizer st) {

        while (st.hasMoreTokens()) {
            String currTok = st.nextToken().trim();

            if (currTok.equals("")) {                    
            } else if (currTok.equals("(") || currTok.equals("[")) {
                // The next token is the parent of a new subtree
                currTok = st.nextToken().trim();
                TreeNode newNode = new TreeNode(parent, currTok);
                parent.addChild(newNode);
                parseString(newNode, st);
            } else if (currTok.equals(")") || currTok.equals("]")) {
                // Finish this subtree
                return;
            } else {
                // An ordinary child node: add it to parent and continue.
                TreeNode newNode = new TreeNode(parent, currTok);
                parent.addChild(newNode);
            }
        }
    }

} public class GPException extends Exception {
        
        /**
         * Constructor for GP exceptions
         * @param string 
         */
	public GPException(String string) {
            super(string);
	}

	/**
	 * 
	 */
	private static final long serialVersionUID = 2346788405132405602L;

} public class Individual implements Serializable {
    private static final long serialVersionUID = -1448696047207647168L;
    private final Genotype genotype;
    // store fitness values as key:value pairs to support multiple objectives
    private LinkedHashMap<String, Double> fitnesses;
    // distance from the origin
    private Double euclideanDistance;
    // the "hypervolume" of a particular individual
    private Double crowdingDistance;
    private Integer dominationCount;
    private double threshold;
    private double crossValAreaROC;	
    private double crossValFitness;
    private double scaledMAE;
    private double scaledMSE;
    private double RT_Cost;
    private double minTrainOutput, maxTrainOutput;
    private ArrayList<Double> estimatedDensPos, estimatedDensNeg;
    ArrayList<String> weights;
    String lassoIntercept;
    /**
     * Create an individual with the given genotype. The new individuals
     * phenotype and fitness are left unspecified.
     * 
     * @param genotype an Instance of a {@link Genotype}
     */
    public Individual(Genotype genotype) {
        this.genotype = genotype;
        this.fitnesses = new LinkedHashMap<String, Double>();
        this.euclideanDistance = Double.MAX_VALUE;
        this.crowdingDistance = 0.0;
        this.dominationCount = 0;
        threshold = 0;
        crossValAreaROC = 0;    
        RT_Cost = 0;
    }

    @SuppressWarnings("unchecked")
    private Individual(Individual i) {
        this.crowdingDistance = i.crowdingDistance;
        this.dominationCount = i.dominationCount;
        this.genotype = i.genotype.copy();
        this.fitnesses = (LinkedHashMap<String, Double>) i.fitnesses.clone();
        this.euclideanDistance = i.euclideanDistance;
        this.threshold = i.threshold;
        this .crossValAreaROC = i.crossValAreaROC;     
        RT_Cost = 0;
        this.weights = i.weights;
    }

    /**
     * deep copy of an individual. THe new individual will be identical to the
     * old in every way, but will be completely independent. So any changes made
     * to the old do not affect the new, and vice versa.
     * 
     * @return new individual
     */
    public Individual copy() {
        return new Individual(this);
    }
    
    /**
     * Return Genotype
     * @return 
     */
    public Genotype getGenotype() {
        return this.genotype;
    }

    /**
     * return fitness names
     * @return 
     */
    public Set<String> getFitnessNames() {
        return this.fitnesses.keySet();
    }
	
    /**
     * Return the LinkedHashMap storing fitness values
     * @return
     */
    public LinkedHashMap<String, Double> getFitnesses() {
        return this.fitnesses;
    }

    /**
     * return fitness function name
     * @param key
     * @return 
     */
    public Double getFitness(String key) {
        return this.fitnesses.get(key);
    }

    /**
     * Overload getFitness to support returning simply the first fitness value if prompted
     * @return
     */
    public Double getFitness() {
        String first = getFirstFitnessKey();
        return getFitness(first);
    }

    /**
     * @return the first fitness key from this individual's stored fitness values
     */
    public String getFirstFitnessKey() {
        Set<String> keys = getFitnesses().keySet();
        return keys.iterator().next();
    }

    /**
     * Set fitness functions
     * @param d 
     */
    public void setFitnesses(LinkedHashMap<String, Double> d) {
        this.fitnesses = d;
    }

    /**
     * Update a particular fitness value
     * @param key
     * @param d the new fitness 
     */
    public void setFitness(String key, Double d) {
        fitnesses.put(key, d);
    }

    /**
     * Obtain the memorized euclidean distance of fitness
     * @return
     */
    public Double getEuclideanDistance() {
        return euclideanDistance;
    }

    /**
     * Calculate the euclidean distance of the fitnesses of this individual from
     * the origin.
     * @param fitnessFunctions
     * @param standardizedMins
     * @param standardizedRanges
     */
    public void calculateEuclideanDistance(LinkedHashMap<String, FitnessFunction> fitnessFunctions,LinkedHashMap<String, Double> standardizedMins,
                    LinkedHashMap<String, Double> standardizedRanges) {
        // reset euclidean distance to 0
        euclideanDistance = 0.0;
        for (String fitnessKey : fitnesses.keySet()) {
            // get fitness converted to minimization if necessary
            Double standardizedFitness = FitnessComparisonStandardizer.getFitnessForMinimization(this, fitnessKey,fitnessFunctions);
            // normalize 
            Double normalizedStandardizedFitness = (standardizedFitness - standardizedMins.get(fitnessKey)) / standardizedRanges.get(fitnessKey);
            // add to euclidean distance
            euclideanDistance += Math.pow(normalizedStandardizedFitness,2);
        }
    }
	
    /**
     * return crowding distance
     * @return 
     */
    public Double getCrowdingDistance() {
        return this.crowdingDistance;
    }

    /**
     * Set crowding distance
     * @param newCrowdingDistance 
     */
    public void setCrowdingDistance(Double newCrowdingDistance) {
        this.crowdingDistance = newCrowdingDistance;
    }

    /**
     * Update the crowding distance of this individual with information from a new fitness function
     * @param localCrowdingDistance
     */
    public void updateCrowdingDistance(Double localCrowdingDistance) {
        crowdingDistance *= localCrowdingDistance;
    }

    /**
     * return domination count
     * @return 
     */
    public Integer getDominationCount() {
        return dominationCount;
    }

    /**
     * increment domination count
     */
    public void incrementDominationCount() {
        dominationCount++;
    }

    /**
     * set domination count
     * @param dominationCount 
     */
    public void setDominationCount(Integer dominationCount) {
        this.dominationCount = dominationCount;
    }

    /**
     * get threshold
     * @return 
     */
    public double getThreshold(){
        return threshold;
    }
    
    /**
     * set threshold
     * @param aThreshold 
     */
    public void setThreshold(double aThreshold){
        threshold = aThreshold;
    }	

    /**
     * get area under the ROC curve on validation data
     * @return 
     */
    public double getCrossValAreaROC(){
        return crossValAreaROC;
    }

    /**
     * set auc on validation data
     * @param anArea 
     */
    public void setCrossValAreaROC(double anArea){
        crossValAreaROC = anArea;
    }

    /**
     * Clear any nonessential memorized values
     */
    public void reset() {
        euclideanDistance = Double.MAX_VALUE;
        crowdingDistance = 0.0;
        dominationCount = 0;
        RT_Cost = 0;
    }
	
    /**
     * return whether 2 individuals are equal
     * @param i
     * @return 
     */
    public Boolean equals(Individual i){
        if (!this.getGenotype().equals(i.getGenotype())){
            return false;
        } else{ // g and p are equal, check fitnesses
            return (this.getFitnesses().equals(i.getFitnesses()));
        }
    }

    /**
     * @return unscaled MATLAB infix string
     */
    @Override
    public String toString() {
        return this.genotype.toString();
    }

    /**
     * @param shouldRound
     * @return scaled MATLAB infix string, with appropriate rounding if necessary
     */
    public String toScaledString(boolean shouldRound) {
        String scaledModel = ((Tree) this.genotype).toScaledString();
        return (shouldRound) ? "round(" + scaledModel + ")" : scaledModel ;
    }
    
    /**
     * get final string representation
     * @param shouldRound
     * @return 
     */
    public String toFinalScaledString(boolean shouldRound) {
        String scaledModel = ((Tree) this.genotype).toFinalScaledString();
        return (shouldRound) ? "round(" + scaledModel + ")" : scaledModel ;
    }

    /**
     * get scaled string representation
     * @return 
     */
    public String toScaledString() {
        return toScaledString(false);
    }
    
    /**
     * get final scaled representation
     * @return 
     */
    public String toFinalScaledString() {
        return toFinalScaledString(false);
    }

    /**
     * set fitness on validatio data
     * @param aFitness 
     */
    public void setScaledCrossValFitness(double aFitness){
        this.crossValFitness = aFitness;
    }

    /**
     * get fitness on validation data
     * @return 
     */
    public double getCrossValFitness(){
        return crossValFitness;
    }

    /**
     * set MSE of the scaled model
     * @param aMSE 
     */
    public void setScaledMSE(double aMSE){
        scaledMSE = aMSE;
    }

    /**
     * get MSE of the scaled model
     * @return 
     */
    public double getScaledMSE(){
        return scaledMSE;
    }
    
    /**
     * set MAE of scaled model
     * @param aMAE 
     */
    public void setScaledMAE(double aMAE){
        scaledMAE = aMAE;
    }

    /**
     * get MAE of the scaled value
     * @return 
     */
    public double getScaledMAE(){
        return scaledMAE;
    }


    /**
     * @return the minTrainOutput
     */
    public double getMinTrainOutput() {
        return minTrainOutput;
    }

    /**
     * @param minTrainOutput the minTrainOutput to set
     */
    public void setMinTrainOutput(double minTrainOutput) {
        this.minTrainOutput = minTrainOutput;
    }

    /**
     * @return the maxTrainOutput
     */
    public double getMaxTrainOutput() {
        return maxTrainOutput;
    }

    /**
     * @param maxTrainOutput the maxTrainOutput to set
     */
    public void setMaxTrainOutput(double maxTrainOutput) {
        this.maxTrainOutput = maxTrainOutput;
    }
    

    /**
     * Set LASSO weights
     * @param aWeights 
     */
    public void setWeights(ArrayList<String> aWeights){
        weights = aWeights;
    }

    /**
     * get LASSO weights
     * @return 
     */
    public ArrayList<String> getWeights(){
        return weights;
    }
    
    /**
     * get intercept
     * @return 
     */
    public String getLassoIntercept(){
        return lassoIntercept;
    }
    
    /**
     * set intercept
     * @param aLassoIntercept 
     */
    public void setLassoIntercept(String aLassoIntercept){
        lassoIntercept = aLassoIntercept;
    }
} public strictfp class MersenneTwisterFast implements Serializable, Cloneable
    {
    // Serialization
    private static final long serialVersionUID = -8219700664442619525L;  // locked as of Version 15
    
    // Period parameters
    private static final int N = 624;
    private static final int M = 397;
    private static final int MATRIX_A = 0x9908b0df;   //    private static final * constant vector a
    private static final int UPPER_MASK = 0x80000000; // most significant w-r bits
    private static final int LOWER_MASK = 0x7fffffff; // least significant r bits


    // Tempering parameters
    private static final int TEMPERING_MASK_B = 0x9d2c5680;
    private static final int TEMPERING_MASK_C = 0xefc60000;
    
    private int mt[]; // the array for the state vector
    private int mti; // mti==N+1 means mt[N] is not initialized
    private int mag01[];
    
    // a good initial seed (of int size, though stored in a long)
    //private static final long GOOD_SEED = 4357;

    private double __nextNextGaussian;
    private boolean __haveNextNextGaussian;
    
    /* We're overriding all internal data, to my knowledge, so this should be okay */
    public Object clone()
        {
        try
            {
            MersenneTwisterFast f = (MersenneTwisterFast)(super.clone());
            f.mt = (int[])(mt.clone());
            f.mag01 = (int[])(mag01.clone());
            return f;
            }
        catch (CloneNotSupportedException e) { throw new InternalError(); } // should never happen
        }
    
    public boolean stateEquals(Object o)
        {
        if (o==this) return true;
        if (o == null || !(o instanceof MersenneTwisterFast))
            return false;
        MersenneTwisterFast other = (MersenneTwisterFast) o;
        if (mti != other.mti) return false;
        for(int x=0;x<mag01.length;x++)
            if (mag01[x] != other.mag01[x]) return false;
        for(int x=0;x<mt.length;x++)
            if (mt[x] != other.mt[x]) return false;
        return true;
        }

    /** Reads the entire state of the MersenneTwister RNG from the stream */
    public void readState(DataInputStream stream) throws IOException
        {
        int len = mt.length;
        for(int x=0;x<len;x++) mt[x] = stream.readInt();
        
        len = mag01.length;
        for(int x=0;x<len;x++) mag01[x] = stream.readInt();
        
        mti = stream.readInt();
        __nextNextGaussian = stream.readDouble();
        __haveNextNextGaussian = stream.readBoolean();
        }
        
    /** Writes the entire state of the MersenneTwister RNG to the stream */
    public void writeState(DataOutputStream stream) throws IOException
        {
        int len = mt.length;
        for(int x=0;x<len;x++) stream.writeInt(mt[x]);
        
        len = mag01.length;
        for(int x=0;x<len;x++) stream.writeInt(mag01[x]);
        
        stream.writeInt(mti);
        stream.writeDouble(__nextNextGaussian);
        stream.writeBoolean(__haveNextNextGaussian);
        }

    /**
     * Constructor using the default seed.
     */
    public MersenneTwisterFast()
        {
        this(System.currentTimeMillis());
        }
    
    /**
     * Constructor using a given seed.  Though you pass this seed in
     * as a long, it's best to make sure it's actually an integer.
     *
     */
    public MersenneTwisterFast(long seed)
        {
        setSeed(seed);
        }
    

    /**
     * Constructor using an array of integers as seed.
     * Your array must have a non-zero length.  Only the first 624 integers
     * in the array are used; if the array is shorter than this then
     * integers are repeatedly used in a wrap-around fashion.
     */
    public MersenneTwisterFast(int[] array)
        {
        setSeed(array);
        }


    /**
     * Initalize the pseudo random number generator.  Don't
     * pass in a long that's bigger than an int (Mersenne Twister
     * only uses the first 32 bits for its seed).   
     */

    synchronized public void setSeed(long seed)
        {
        // Due to a bug in java.util.Random clear up to 1.2, we're
        // doing our own Gaussian variable.
        __haveNextNextGaussian = false;

        mt = new int[N];
        
        mag01 = new int[2];
        mag01[0] = 0x0;
        mag01[1] = MATRIX_A;

        mt[0]= (int)(seed & 0xffffffff);
        for (mti=1; mti<N; mti++) 
            {
            mt[mti] = 
                (1812433253 * (mt[mti-1] ^ (mt[mti-1] >>> 30)) + mti); 
            /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
            /* In the previous versions, MSBs of the seed affect   */
            /* only MSBs of the array mt[].                        */
            /* 2002/01/09 modified by Makoto Matsumoto             */
            // mt[mti] &= 0xffffffff;
            /* for >32 bit machines */
            }
        }


    /**
     * Sets the seed of the MersenneTwister using an array of integers.
     * Your array must have a non-zero length.  Only the first 624 integers
     * in the array are used; if the array is shorter than this then
     * integers are repeatedly used in a wrap-around fashion.
     */

    synchronized public void setSeed(int[] array)
        {
        if (array.length == 0)
            throw new IllegalArgumentException("Array length must be greater than zero");
        int i, j, k;
        setSeed(19650218);
        i=1; j=0;
        k = (N>array.length ? N : array.length);
        for (; k!=0; k--) 
            {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >>> 30)) * 1664525)) + array[j] + j; /* non linear */
            // mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
            i++;
            j++;
            if (i>=N) { mt[0] = mt[N-1]; i=1; }
            if (j>=array.length) j=0;
            }
        for (k=N-1; k!=0; k--) 
            {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >>> 30)) * 1566083941)) - i; /* non linear */
            // mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
            i++;
            if (i>=N) 
                {
                mt[0] = mt[N-1]; i=1; 
                }
            }
        mt[0] = 0x80000000; /* MSB is 1; assuring non-zero initial array */ 
        }


    public int nextInt()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return y;
        }



    public short nextShort()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (short)(y >>> 16);
        }



    public char nextChar()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (char)(y >>> 16);
        }


    public boolean nextBoolean()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (boolean)((y >>> 31) != 0);
        }



    /** This generates a coin flip with a probability <tt>probability</tt>
        of returning true, else returning false.  <tt>probability</tt> must
        be between 0.0 and 1.0, inclusive.   Not as precise a random real
        event as nextBoolean(double), but twice as fast. To explicitly
        use this, remember you may need to cast to float first. */

    public boolean nextBoolean(float probability)
        {
        int y;
        
        if (probability < 0.0f || probability > 1.0f)
            throw new IllegalArgumentException ("probability must be between 0.0 and 1.0 inclusive.");
        if (probability==0.0f) return false;            // fix half-open issues
        else if (probability==1.0f) return true;        // fix half-open issues
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (y >>> 8) / ((float)(1 << 24)) < probability;
        }


    /** This generates a coin flip with a probability <tt>probability</tt>
        of returning true, else returning false.  <tt>probability</tt> must
        be between 0.0 and 1.0, inclusive. */

    public boolean nextBoolean(double probability)
        {
        int y;
        int z;

        if (probability < 0.0 || probability > 1.0)
            throw new IllegalArgumentException ("probability must be between 0.0 and 1.0 inclusive.");
        if (probability==0.0) return false;             // fix half-open issues
        else if (probability==1.0) return true; // fix half-open issues
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
            
            mti = 0;
            }
        
        z = mt[mti++];
        z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
        z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
        z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
        z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
        
        /* derived from nextDouble documentation in jdk 1.2 docs, see top */
        return ((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53) < probability;
        }


    public byte nextByte()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (byte)(y >>> 24);
        }


    public void nextBytes(byte[] bytes)
        {
        int y;
        
        for (int x=0;x<bytes.length;x++)
            {
            if (mti >= N)   // generate N words at one time
                {
                int kk;
                final int[] mt = this.mt; // locals are slightly faster 
                final int[] mag01 = this.mag01; // locals are slightly faster 
                
                for (kk = 0; kk < N - M; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                for (; kk < N-1; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
                
                mti = 0;
                }
            
            y = mt[mti++];
            y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
            y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
            y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
            y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

            bytes[x] = (byte)(y >>> 24);
            }
        }


    public long nextLong()
        {
        int y;
        int z;

        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
            
            mti = 0;
            }
        
        z = mt[mti++];
        z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
        z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
        z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
        z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
        
        return (((long)y) << 32) + (long)z;
        }



    /** Returns a long drawn uniformly from 0 to n-1.  Suffice it to say,
        n must be > 0, or an IllegalArgumentException is raised. */
    public long nextLong(long n)
        {
        if (n<=0)
            throw new IllegalArgumentException("n must be positive, got: " + n);
        
        long bits, val;
        do 
            {
            int y;
            int z;
    
            if (mti >= N)   // generate N words at one time
                {
                int kk;
                final int[] mt = this.mt; // locals are slightly faster 
                final int[] mag01 = this.mag01; // locals are slightly faster 
            
                for (kk = 0; kk < N - M; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                for (; kk < N-1; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
    
                mti = 0;
                }
    
            y = mt[mti++];
            y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
            y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
            y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
            y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)
    
            if (mti >= N)   // generate N words at one time
                {
                int kk;
                final int[] mt = this.mt; // locals are slightly faster 
                final int[] mag01 = this.mag01; // locals are slightly faster 
                
                for (kk = 0; kk < N - M; kk++)
                    {
                    z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
                    }
                for (; kk < N-1; kk++)
                    {
                    z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
                    }
                z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
                
                mti = 0;
                }
            
            z = mt[mti++];
            z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
            z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
            z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
            z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
            
            bits = (((((long)y) << 32) + (long)z) >>> 1);
            val = bits % n;
            } while (bits - val + (n-1) < 0);
        return val;
        }

    /** Returns a random double in the half-open range from [0.0,1.0).  Thus 0.0 is a valid
        result but 1.0 is not. */
    public double nextDouble()
        {
        int y;
        int z;

        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
                }
            z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
            
            mti = 0;
            }
        
        z = mt[mti++];
        z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
        z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
        z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
        z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
        
        /* derived from nextDouble documentation in jdk 1.2 docs, see top */
        return ((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53);
        }



    /** Returns a double in the range from 0.0 to 1.0, possibly inclusive of 0.0 and 1.0 themselves.  Thus:

        <p><table border=0>
        <th><td>Expression<td>Interval
        <tr><td>nextDouble(false, false)<td>(0.0, 1.0)
        <tr><td>nextDouble(true, false)<td>[0.0, 1.0)
        <tr><td>nextDouble(false, true)<td>(0.0, 1.0]
        <tr><td>nextDouble(true, true)<td>[0.0, 1.0]
        </table>
        
        <p>This version preserves all possible random values in the double range.
    */
    public double nextDouble(boolean includeZero, boolean includeOne)
        {
        double d = 0.0;
        do
            {
            d = nextDouble();                           // grab a value, initially from half-open [0.0, 1.0)
            if (includeOne && nextBoolean()) d += 1.0;  // if includeOne, with 1/2 probability, push to [1.0, 2.0)
            } 
        while ( (d > 1.0) ||                            // everything above 1.0 is always invalid
            (!includeZero && d == 0.0));            // if we're not including zero, 0.0 is invalid
        return d;
        }



    public double nextGaussian()
        {
        if (__haveNextNextGaussian)
            {
            __haveNextNextGaussian = false;
            return __nextNextGaussian;
            } 
        else 
            {
            double v1, v2, s;
            do 
                { 
                int y;
                int z;
                int a;
                int b;
                    
                if (mti >= N)   // generate N words at one time
                    {
                    int kk;
                    final int[] mt = this.mt; // locals are slightly faster 
                    final int[] mag01 = this.mag01; // locals are slightly faster 
                    
                    for (kk = 0; kk < N - M; kk++)
                        {
                        y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                        }
                    for (; kk < N-1; kk++)
                        {
                        y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                        }
                    y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                    mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
                    
                    mti = 0;
                    }
                
                y = mt[mti++];
                y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
                y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
                y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
                y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)
                
                if (mti >= N)   // generate N words at one time
                    {
                    int kk;
                    final int[] mt = this.mt; // locals are slightly faster 
                    final int[] mag01 = this.mag01; // locals are slightly faster 
                    
                    for (kk = 0; kk < N - M; kk++)
                        {
                        z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
                        }
                    for (; kk < N-1; kk++)
                        {
                        z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
                        }
                    z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                    mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
                    
                    mti = 0;
                    }
                
                z = mt[mti++];
                z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
                z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
                z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
                z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
                
                if (mti >= N)   // generate N words at one time
                    {
                    int kk;
                    final int[] mt = this.mt; // locals are slightly faster 
                    final int[] mag01 = this.mag01; // locals are slightly faster 
                    
                    for (kk = 0; kk < N - M; kk++)
                        {
                        a = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+M] ^ (a >>> 1) ^ mag01[a & 0x1];
                        }
                    for (; kk < N-1; kk++)
                        {
                        a = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+(M-N)] ^ (a >>> 1) ^ mag01[a & 0x1];
                        }
                    a = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                    mt[N-1] = mt[M-1] ^ (a >>> 1) ^ mag01[a & 0x1];
                    
                    mti = 0;
                    }
                
                a = mt[mti++];
                a ^= a >>> 11;                          // TEMPERING_SHIFT_U(a)
                a ^= (a << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(a)
                a ^= (a << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(a)
                a ^= (a >>> 18);                        // TEMPERING_SHIFT_L(a)
                
                if (mti >= N)   // generate N words at one time
                    {
                    int kk;
                    final int[] mt = this.mt; // locals are slightly faster 
                    final int[] mag01 = this.mag01; // locals are slightly faster 
                    
                    for (kk = 0; kk < N - M; kk++)
                        {
                        b = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+M] ^ (b >>> 1) ^ mag01[b & 0x1];
                        }
                    for (; kk < N-1; kk++)
                        {
                        b = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                        mt[kk] = mt[kk+(M-N)] ^ (b >>> 1) ^ mag01[b & 0x1];
                        }
                    b = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                    mt[N-1] = mt[M-1] ^ (b >>> 1) ^ mag01[b & 0x1];
                    
                    mti = 0;
                    }
                
                b = mt[mti++];
                b ^= b >>> 11;                          // TEMPERING_SHIFT_U(b)
                b ^= (b << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(b)
                b ^= (b << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(b)
                b ^= (b >>> 18);                        // TEMPERING_SHIFT_L(b)
                
                /* derived from nextDouble documentation in jdk 1.2 docs, see top */
                v1 = 2 *
                    (((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53))
                    - 1;
                v2 = 2 * (((((long)(a >>> 6)) << 27) + (b >>> 5)) / (double)(1L << 53))
                    - 1;
                s = v1 * v1 + v2 * v2;
                } while (s >= 1 || s==0);
            double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
            __nextNextGaussian = v2 * multiplier;
            __haveNextNextGaussian = true;
            return v1 * multiplier;
            }
        }
    
    
    
    

    /** Returns a random float in the half-open range from [0.0f,1.0f).  Thus 0.0f is a valid
        result but 1.0f is not. */
    public float nextFloat()
        {
        int y;
        
        if (mti >= N)   // generate N words at one time
            {
            int kk;
            final int[] mt = this.mt; // locals are slightly faster 
            final int[] mag01 = this.mag01; // locals are slightly faster 
            
            for (kk = 0; kk < N - M; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            for (; kk < N-1; kk++)
                {
                y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                }
            y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];

            mti = 0;
            }
  
        y = mt[mti++];
        y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
        y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
        y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
        y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)

        return (y >>> 8) / ((float)(1 << 24));
        }


    /** Returns a float in the range from 0.0f to 1.0f, possibly inclusive of 0.0f and 1.0f themselves.  Thus:

        <p><table border=0>
        <th><td>Expression<td>Interval
        <tr><td>nextFloat(false, false)<td>(0.0f, 1.0f)
        <tr><td>nextFloat(true, false)<td>[0.0f, 1.0f)
        <tr><td>nextFloat(false, true)<td>(0.0f, 1.0f]
        <tr><td>nextFloat(true, true)<td>[0.0f, 1.0f]
        </table>
        
        <p>This version preserves all possible random values in the float range.
    */
    public float nextFloat(boolean includeZero, boolean includeOne)
        {
        float d = 0.0f;
        do
            {
            d = nextFloat();                            // grab a value, initially from half-open [0.0f, 1.0f)
            if (includeOne && nextBoolean()) d += 1.0f; // if includeOne, with 1/2 probability, push to [1.0f, 2.0f)
            } 
        while ( (d > 1.0f) ||                           // everything above 1.0f is always invalid
            (!includeZero && d == 0.0f));           // if we're not including zero, 0.0f is invalid
        return d;
        }



    /** Returns an integer drawn uniformly from 0 to n-1.  Suffice it to say,
        n must be > 0, or an IllegalArgumentException is raised. */
    public int nextInt(int n)
        {
        if (n<=0)
            throw new IllegalArgumentException("n must be positive, got: " + n);
        
        if ((n & -n) == n)  // i.e., n is a power of 2
            {
            int y;
        
            if (mti >= N)   // generate N words at one time
                {
                int kk;
                final int[] mt = this.mt; // locals are slightly faster 
                final int[] mag01 = this.mag01; // locals are slightly faster 
                
                for (kk = 0; kk < N - M; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                for (; kk < N-1; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
                
                mti = 0;
                }
            
            y = mt[mti++];
            y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
            y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
            y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
            y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)
            
            return (int)((n * (long) (y >>> 1) ) >> 31);
            }
        
        int bits, val;
        do 
            {
            int y;
            
            if (mti >= N)   // generate N words at one time
                {
                int kk;
                final int[] mt = this.mt; // locals are slightly faster 
                final int[] mag01 = this.mag01; // locals are slightly faster 
                
                for (kk = 0; kk < N - M; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                for (; kk < N-1; kk++)
                    {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
                    mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
                    }
                y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
                
                mti = 0;
                }
            
            y = mt[mti++];
            y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
            y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
            y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
            y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)
        
            bits = (y >>> 1);
            val = bits % n;
            } while(bits - val + (n-1) < 0);
        return val;
        }
    } public class Population extends ArrayList<Individual> implements Serializable {
	private static final long serialVersionUID = 6111020814262385165L;

        /**
         * default constructor
         */
	public Population() {
		super();
	}

	/**
         * a constructor which yields the combination of populations
         * @param populations 
         */
        public Population(Population... populations) {
            super();
            for (Population population : populations) {
                this.addAll(population);
            }
	}

        /**
         * equals method
         * @param other
         * @return 
         */
	@Override
	public boolean equals(Object other) {
		if (!this.getClass().equals(other.getClass()))
			return false;
		ArrayList<Individual> otherL = (ArrayList<Individual>) other;
		if (otherL.size() != this.size())
			return false;
		for (int i = 0; i < otherL.size(); i++) {
			if (!this.get(i).equals(otherL.get(i)))
				return false;
		}
		return true;
	}

        /**
         * perform a deep copy of the population
         * @return 
         */
	public Population deepCopy() {
		ByteArrayOutputStream bos = new ByteArrayOutputStream();
		ObjectOutputStream oos;
		try {
			oos = new ObjectOutputStream(bos);
			oos.writeObject(this);
			ObjectInputStream ois = new ObjectInputStream(new ByteArrayInputStream(bos.toByteArray()));
			return (Population) ois.readObject();
		} catch (IOException e) {
			System.err.println("IO Exception in deepCopy");
		} catch (ClassNotFoundException e) {
			System.err.println("Class not found Exception in deepCopy");
		}
		return null;
	}

	/**
	 * Calculates the euclidean distance of all individuals in the first front. Others will keep a distance of Double.MAX_VALUE
	 * @param fitnessFunctions
	 */
	public void calculateEuclideanDistances(LinkedHashMap<String, FitnessFunction> fitnessFunctions) {
		// first get the mins and maxes for the first front only
		LinkedHashMap<String, Double>[] minMax = getMinMax(this, fitnessFunctions, true);
		LinkedHashMap<String, Double> mins = minMax[0];
		LinkedHashMap<String, Double> maxes = minMax[1];
		// convert the mins and maxes to standardized form
		for (String key : mins.keySet()) {
			// get new standardized min and max. swap values if this fitness function isn't minimizing already
			if (fitnessFunctions.get(key).isMaximizingFunction()) { // swap min and max since they've both been inverted
				Double standardizedMin = FitnessComparisonStandardizer.getFitnessForMinimization(maxes, key, fitnessFunctions);
				Double standardizedMax = FitnessComparisonStandardizer.getFitnessForMinimization(mins, key, fitnessFunctions);
				mins.put(key, standardizedMin);
				maxes.put(key, standardizedMax);
			} else {
				Double standardizedMin = FitnessComparisonStandardizer.getFitnessForMinimization(mins, key, fitnessFunctions);
				Double standardizedMax = FitnessComparisonStandardizer.getFitnessForMinimization(maxes, key, fitnessFunctions);
				mins.put(key, standardizedMin);
				maxes.put(key, standardizedMax);
			}
		}
		// create the ranges needed for scaling
		LinkedHashMap<String, Double> ranges = new LinkedHashMap<String, Double>();
		for (String key : mins.keySet()) {
			ranges.put(key, Math.abs(maxes.get(key) - mins.get(key)));
		}
		// compute euclidean distances for the first front only
		for (Individual individual : this) {
			if (individual.getDominationCount() > 0) continue;
			individual.calculateEuclideanDistance(fitnessFunctions, mins, ranges);
		}
	}
	
        /**
	 * Find the min and max value for each fitness function
	 * 
	 * @param pop
         * @param fitnessFunctions
	 * @param onlyFirstFront if true, mins/maxes calculated only for individuals in the first front
         * @return 
	 */
	public static LinkedHashMap<String, Double>[] getMinMax(Population pop,LinkedHashMap<String, FitnessFunction> fitnessFunctions, Boolean onlyFirstFront) {
		LinkedHashMap<String, Double>[] minMax = new LinkedHashMap[2];
		LinkedHashMap<String, Double> mins = new LinkedHashMap<String, Double>();
		LinkedHashMap<String, Double> maxes = new LinkedHashMap<String, Double>();
		
		// establish order of fitness functions in mins and maxes
		for (String id : fitnessFunctions.keySet()) {
			mins.put(id, pop.get(0).getFitness(id));
			maxes.put(id, pop.get(0).getFitness(id));
		}

		// find mins and maxes
		for (Individual i : pop) {
			if (onlyFirstFront && i.getDominationCount() > 0) continue;
			for (String funcName : i.getFitnessNames()) {
				Double iFitness = i.getFitness(funcName);
				if (iFitness < mins.get(funcName)) // lower min
					mins.put(funcName, iFitness);
				if (iFitness > maxes.get(funcName)) // higher max
					maxes.put(funcName, iFitness);
			}
		}
		minMax[0] = mins;
		minMax[1] = maxes;
		return minMax;
	}
	
        /**
         * Non-dominating sort of the population
         */
        public void sort() {
            Collections.sort(this, new DominationCrowdingSortComparator(false));
	}

    /**
     * 
     * @return 
     */
    @Override
    public int hashCode() {
        int hash = 7;
        return hash;
    }
        
        /**
         * Comparator to compare the domination relationship between individuals of the population
         */
	public class DominationCrowdingSortComparator implements Comparator<Individual> {
		// can't use crowding distance for comparison if it's not being computed a level above
		boolean crowdingDistanceEnabled;

                /**
                 * default constructor
                 */
		public DominationCrowdingSortComparator() {
			this(false);
		}
                
                /**
                 * whether to use crowding distance to break ties
                 * @param _crowdingDistanceEnabled 
                 */
		public DominationCrowdingSortComparator(boolean _crowdingDistanceEnabled) {
			crowdingDistanceEnabled = _crowdingDistanceEnabled;
		}
		
                /**
                 * compare two individuals
                 * @param a
                 * @param b
                 * @return 
                 */
		@Override
		public int compare(Individual a, Individual b) {
			Integer ad = a.getDominationCount();
			Integer bd = b.getDominationCount();
			// sort by domination count: lower is better
			if (ad > bd) return 1;
			else if (ad < bd) return -1;
			// don't proceed to crowding distance sorting if we're not computing crowding distance
			if (!crowdingDistanceEnabled) return 0;
			// now sort by crowding distance: higher is better
			Double ac = a.getCrowdingDistance();
			Double bc = b.getCrowdingDistance();
			if (ac < bc) return 1;
			else if (ac > bc) return -1;
			else return 0;
		}
	}
	
        /**
         * this class sorts individuals according to their fitness in validation data
         */
	public class CrossValSortComparator implements Comparator<Individual>{
                /**
                 * default constructor
                 */
		public CrossValSortComparator() {
		}
		
		/**
                 * compare 2 individuals
                 * @param a
                 * @param b
                 * @return 
                 */
                @Override
		public int compare(Individual a, Individual b) {
			double aCVF = a.getCrossValFitness();
			double bCVF = b.getCrossValFitness();
			// sort by domination count: lower is better
                        if (aCVF == bCVF) {
                            return 0;
                        } else if (aCVF > bCVF){
                            return 1;
                        }else {
                            return -1;
                        }
		}
	}
        
        /**
         * sort individuals according to non-domination
         * @param crowdingSortEnabled 
         */
	public void sort(boolean crowdingSortEnabled) {
            Collections.sort(this, new DominationCrowdingSortComparator(crowdingSortEnabled));
	}
        
        /**
         * sort individuals according to the fitness on validation
         */
        public void sortCrossVal() {
            Collections.sort(this, new CrossValSortComparator());
	}
                
} public interface Crossover {

	/**
	 * Perform crossover between two individuals, returning their offspring. The
	 * two individuals may be thought of as the parents. The returned population
	 * may contain one or more offspring, depending on the crossover
	 * implementation. Neither parents shall be mutated - copies of their genome
	 * should be recombined to form offspring.
	 * 
	 * @param ind1 The first parent
	 * @param ind2 The second parent
	 * @return population The offspring. We return a population because there
	 *         may be more than one offspring returned.
	 * @throws GPException if this operator is not compatible with the genotype
	 *         of one or more of the parents
	 */
	public abstract Population crossOver(Individual ind1, Individual ind2)
			throws GPException;
} public class CrowdedTournamentSelection extends TournamentSelection {

    /**
     * Constructor
     * @param rand
     * @param props 
     */
    public CrowdedTournamentSelection(MersenneTwisterFast rand,Properties props) {
        super(rand, props);
    }

    /**
     * Perform crowded tournament selection
     * 
     * Note: this depends on the NonDominationRank, which is memoized from the
     * last call to NonDominatedSort. It also depends on the crowding distance,
     * which is memoized from the call to CrowdedSort.computeCrowdingDistances
     * @param pop
     * @return 
     */
    @Override
    public Individual select(Population pop) {
        int n = pop.size();
        Individual best, challenger;
        best = pop.get(rand.nextInt(n));
        for (int j = 0; j < TOURNEY_SIZE - 1; j++) {
            challenger = pop.get(rand.nextInt(n));
            // challenger wins if it dominates best
            if (challenger.getDominationCount() < best.getDominationCount()) {
                best = challenger;
                // or if neither dominates the other (same nondom rank) and
                // challenger has higher crowding distance
            } else if ((challenger.getDominationCount().equals(best.getDominationCount()))&& (challenger.getCrowdingDistance() > best.getCrowdingDistance())) {
                best = challenger;
            }
        }
        return best;
    }
} public interface Initialize {

	/**
	 * Instantiate a new population of individuals
	 * 
	 * @param popSize Number of individuals in the new population.
	 * @return the new population.
	 */
	public abstract Population initialize(int popSize);
} public interface Mutate {

	/**
	 * Produce a new individual by randomly mutating i. Implementing class
	 * encodes particular mutation procedure. Individual i is not mutated.
	 * Instead, a new individual is cloned from i, then mutated.
	 * 
	 * @param i Individual to mutate from.
	 * @return New individual
	 * @throws GPException
	 */
	public abstract Individual mutate(Individual i) throws GPException;
} public interface Select {

	/**
	 * Given a population, select an individual and return it.
	 * 
	 * @param pop the population to select an individual from.
	 * @return the selected individual.
	 */
	public abstract Individual select(Population pop);
} public class SinglePointKozaCrossover extends SinglePointUniformCrossover {

	private final double USE_FUNCTION;

	/**
	 * A new mutation operator. The probability for selecting to crossover at an
	 * internal node is specified by the
	 * {@value algorithm.Parameters.Names#KOZA_FUNC_RATE} key in the properties
	 * file. If no value is specified, then the default value of
	 * {@value algorithm.Parameters.Defaults#KOZA_FUNC_RATE} is used.
	 * 
	 * @param rand random number generator.
	 * @param props properties of the system.
	 */
	public SinglePointKozaCrossover(MersenneTwisterFast rand, Properties props) {
		super(rand, props);
		if (props.containsKey(Parameters.Names.KOZA_FUNC_RATE)) {
			USE_FUNCTION = Double.valueOf(props
					.getProperty(Parameters.Names.KOZA_FUNC_RATE));
		} else {
			USE_FUNCTION = Parameters.Defaults.KOZA_FUNC_RATE;
		}
	}

	/**
	 * Implement our biased point selection.
         * @return 
	 */
	@Override
	protected TreeNode selectXOverPt(Tree t) {
		List<TreeNode> nodes = t.getRoot().depthFirstTraversal();
		// Separate nodes into function nodes and terminal nodes.
		List<TreeNode> function_nodes = new ArrayList<TreeNode>();
		List<TreeNode> term_nodes = new ArrayList<TreeNode>();
		for (TreeNode n : nodes) {
			if (n.children.isEmpty())
				term_nodes.add(n);
			else
				function_nodes.add(n);
		}
		// determine which set to select from, then uniformly select from that
		// set.
		if (!function_nodes.isEmpty() && rand.nextDouble() <= USE_FUNCTION)
			return function_nodes.get(rand.nextInt(function_nodes.size()));
		else
			return term_nodes.get(rand.nextInt(term_nodes.size()));
	}

} public class SinglePointUniformCrossover extends RandomOperator implements Crossover {

    private final int TREE_XOVER_MAX_DEPTH;
    private final int TREE_XOVER_TRIES;

    /**
     * Create crossover operator which swaps subtrees at uniformly selected
     * points in two parent individuals. There are two parameters for this
     * operator.
     * <ul>
     * <li>The maximum depth allowed in new trees, specified by the key
     * {@value algorithm.Parameters.Names#TREE_XOVER_MAX_DEPTH}, which defaults
     * to {@value algorithm.Parameters.Defaults#TREE_XOVER_MAX_DEPTH}.
     * <li>The maximum number of times to try the crossover before giving up. A
     * crossover attempt fails if either of the new trees is too deep. This
     * parameter is specified by the key
     * {@value algorithm.Parameters.Names#TREE_XOVER_TRIES}, which defaults to
     * {@value algorithm.Parameters.Defaults#TREE_XOVER_TRIES}.
     * 
     * @param rand
     * @param props
     */
    public SinglePointUniformCrossover(MersenneTwisterFast rand,Properties props) {
        super(rand);
        if (props.containsKey(Parameters.Names.TREE_XOVER_MAX_DEPTH))
            TREE_XOVER_MAX_DEPTH = Integer.valueOf(props.getProperty(Parameters.Names.TREE_XOVER_MAX_DEPTH));
        else
            TREE_XOVER_MAX_DEPTH = Parameters.Defaults.TREE_XOVER_MAX_DEPTH;
        if (props.containsKey(Parameters.Names.TREE_XOVER_TRIES))
            TREE_XOVER_TRIES = Integer.valueOf(props.getProperty(Parameters.Names.TREE_XOVER_TRIES));
        else
            TREE_XOVER_TRIES = Parameters.Defaults.TREE_XOVER_TRIES;
    }

    @Override
    public Population crossOver(Individual ind1, Individual ind2) throws GPException {
            if (!(ind1.getGenotype() instanceof Tree && ind2.getGenotype() instanceof Tree)) {
                throw new GPException("attempting SinglePointUniformCrossover of two genotypes not of type Tree");
            }
		Tree c1, c2;
		int tries = 0;
		do {
			c2 = (Tree) ind2.getGenotype().copy();
			c1 = (Tree) ind1.getGenotype().copy();
			// pick a xover pt in this Tree by uniform sampling of
			// depthFirstTraversal, and find that node's index among its
			// siblings
			TreeNode xoverPt1 = selectXOverPt(c1);
			int xoverPt1idxInChildren = xoverPt1.parent.children
					.indexOf(xoverPt1);

			// same for other Tree
			TreeNode xoverPt2 = selectXOverPt(c2);
			int xoverPt2idxInChildren = xoverPt2.parent.children
					.indexOf(xoverPt2);

			// other.xoverpt = this.xoverpt, and fix up parent link
			xoverPt2.parent.children.set(xoverPt2idxInChildren, xoverPt1);
			TreeNode tmpParent = xoverPt1.parent;
			xoverPt1.parent = xoverPt2.parent;

			// this.xoverpt = other.xoverpt, and fix up parent link
			tmpParent.children.set(xoverPt1idxInChildren, xoverPt2);
			xoverPt2.parent = tmpParent;

			// reset cached values
			xoverPt1.reset();
			xoverPt2.reset();
			tries++;
		} while ((c1.getDepth() > TREE_XOVER_MAX_DEPTH || c2.getDepth() > TREE_XOVER_MAX_DEPTH) && tries < TREE_XOVER_TRIES);
		Population twoPop = new Population();
		if (tries >= TREE_XOVER_TRIES) {
			// System.out.println("failed to xover properly"+ c1.getDepth() +
			// " " +c2.getDepth());

			// One of the children might be of permissible depth...
			if (c1.getDepth() > TREE_XOVER_MAX_DEPTH)
				twoPop.add(ind1.copy());
			else
				twoPop.add(new Individual(c1));

			if (c1.getDepth() > TREE_XOVER_MAX_DEPTH)
				twoPop.add(ind2.copy());
			else
				twoPop.add(new Individual(c2));
		} else {
			Individual i1 = new Individual(c1);
			Individual i2 = new Individual(c2);
			// clear any nonessential memoized values
			i1.reset();
			i2.reset();
			twoPop.add(i1);
			twoPop.add(i2);
		}
		return twoPop;
	}

	/**
	 * Select point (node) uniformly in the given tree.
	 * 
	 * @param t Tree to select from
	 * @return Chosen node
	 */
	protected TreeNode selectXOverPt(Tree t) {
		List<TreeNode> nodes = t.getRoot().depthFirstTraversal();
		int nodeCount = nodes.size();
		int xoverPt1idx = rand.nextInt(nodeCount);
		return nodes.get(xoverPt1idx);
	}

} public class SubtreeMutate extends RandomOperator implements Mutate {

	private final int TREE_MUTATE_MAX_DEPTH;
	private final TreeGenerator treeGen;

    /**
     * Construct a new mutation operator. The maximum depth of the new tree of
     * the returned individual is specified by the value at the key
     * {@value algorithm.Parameters.Names#TREE_MUTATE_MAX_DEPTH} in the
     * properties file. If no value is specified, the default depth limit is
     * {@value algorithm.Parameters.Defaults#TREE_MUTATE_MAX_DEPTH}.
     * 
     * @param rand random number generator.
     * @param props object encoding system properties.
     * @param TGen generator to use for growing new subtrees.
     */
    public SubtreeMutate(MersenneTwisterFast rand, Properties props,TreeGenerator TGen) {
        super(rand);
        if (props.containsKey(Parameters.Names.TREE_MUTATE_MAX_DEPTH))
            TREE_MUTATE_MAX_DEPTH = Integer.valueOf(props.getProperty(Parameters.Names.TREE_MUTATE_MAX_DEPTH));
        else
            TREE_MUTATE_MAX_DEPTH = Parameters.Defaults.TREE_MUTATE_MAX_DEPTH;
        treeGen = TGen;
    }

    @Override
    public Individual mutate(Individual i) throws GPException {
        if (!(i.getGenotype() instanceof Tree)) {
            throw new GPException("attempting SubtreeMutate of genotype not of type Tree");
        }
        Tree copy = (Tree) i.getGenotype().copy();
        ArrayList<TreeNode> treeNodes = copy.getRoot().depthFirstTraversal();
        int whichNode = rand.nextInt(treeNodes.size());
        TreeNode n = treeNodes.get(whichNode);
        // System.out.println("Node picked is: " + whichNode + "; " +
        // n.toString());
        int curDepth = n.getDepth();
        // grow a new subtree at the selected point
        // if selected point is variable
        treeGen.generate(n, TREE_MUTATE_MAX_DEPTH - curDepth, false);
        // System.out.println("New subtree there is: " + n.toStringAsTree());

        return new Individual(copy);
    }

} public class TournamentSelection extends RandomOperator implements Select {

    /**
     * Size of the tournament to run.
     */
    protected final int TOURNEY_SIZE;


    /**
     * Create a select operator, of the specified tournament size.
     * 
     * @param rand instance of {@link MersenneTwisterFast} to use for random
     *        number generation.
     * @param props
     */
    public TournamentSelection(MersenneTwisterFast rand, Properties props) {
        super(rand);
        if (props.containsKey(Parameters.Names.TOURNEY_SIZE))
            TOURNEY_SIZE = Integer.valueOf(props.getProperty(Parameters.Names.TOURNEY_SIZE));
        else
            TOURNEY_SIZE = Parameters.Defaults.TOURNEY_SIZE;
    }

    @Override
    public Individual select(Population pop) {
        int n = pop.size();
        // want newPop of size n, but have to account for elitism
        Individual best, challenger;
        best = pop.get(rand.nextInt(n));
        if (best.getDominationCount() == null) {
            System.err.println("Individual has no non-domination rank");
        }
        for (int j = 0; j < TOURNEY_SIZE - 1; j++) {
            challenger = pop.get(rand.nextInt(n));
            if (challenger.getDominationCount() == null) {
                System.err.println("Individual has no non-domination rank");
            }
            if (challenger.getDominationCount() < best.getDominationCount())
                best = challenger;
            // if the individuals are in the same front, randomize selection
            else if (challenger.getDominationCount() == best.getDominationCount()) {
                if (rand.nextBoolean(0.5)) {
                        best = challenger;
                }
            }
        }
        return best;
    }

} public class TreeInitialize extends RandomOperator implements Initialize {
	private final int TREE_INITIAL_MAX_DEPTH;
	private final TreeGenerator treeGen;

	/**
	 * Create new intialize operator. The generated trees are limited by a
	 * maximum allowed depth. This limit is specified by the value at the key
	 * {@value algorithm.Parameters.Names#TREE_MUTATE_MAX_DEPTH} in the
	 * properties file. If there is no value specified, the default limit of
	 * {@value algorithm.Parameters.Defaults#TREE_INIT_MAX_DEPTH} is used.
	 * 
	 * @param rand random number generator instance.
	 * @param props properties file.
	 * @param TGen generator for creating the new trees.
	 */
	public TreeInitialize(MersenneTwisterFast rand, Properties props,TreeGenerator TGen) {
		super(rand);
		if (props.containsKey(Parameters.Names.TREE_INIT_MAX_DEPTH))
                    TREE_INITIAL_MAX_DEPTH = Integer.valueOf(props.getProperty(Parameters.Names.TREE_INIT_MAX_DEPTH));
		else
                    TREE_INITIAL_MAX_DEPTH = Parameters.Defaults.TREE_INIT_MAX_DEPTH;
		this.treeGen = TGen;
	}

	@Override
	public Population initialize(int popSize) {
		Population ret = new Population();
		// number of individuals to generate per depth, per type (grow or full)
		int inds_per_d = popSize / 2 / TREE_INITIAL_MAX_DEPTH;

		for (int depth = 1; depth <= TREE_INITIAL_MAX_DEPTH; depth++) {
			for (int i = 0; i < inds_per_d; i++) {
                            Tree t = treeGen.generateTree(depth, true);
                            Individual ind = new Individual(t);
                            ret.add(ind);
                            //ret.add(new Individual(treeGen.generateTree(depth, true)));
                            ret.add(new Individual(treeGen.generateTree(depth, false)));
			}
		}
		return ret;
	}

} public class CrowdingSort {

        /**
         * arbitrary large distance assigned to boundary solutions
         */
        public static final Double BOUNDARY_DISTANCE = 1000000000.0;

	/**
	 * Calculate cuboid region which is linearly proportional to crowding distance
         * @param p
         * @param f
	 */
	public static void computeCrowdingDistances(Population p, LinkedHashMap<String, FitnessFunction> f) {
		// for each fitness function:
		// first use comparator to generate sort indices. Pass population to
		// comparator, use it to access the population when making comparisons
		// get difference between first and last individual for normalization
		// then for each individual, update the crowding distance with the distance between neighbors, divided by max-min
		for (String fitnessFuncName : f.keySet()) {
			IntegerComparator ic = new IntegerComparator(p, fitnessFuncName);
			ArrayList<Integer> indices = getRange(0, p.size());
			Collections.sort(indices, ic);
			for (int metaIndex = 0; metaIndex < indices.size(); metaIndex++) {
				Integer previousMetaIndex = metaIndex - 1;
				Integer nextMetaIndex = metaIndex + 1;
				Double localCrowdingDistance;
				if ((previousMetaIndex < 0) || (nextMetaIndex >= p.size()))
					localCrowdingDistance = CrowdingSort.BOUNDARY_DISTANCE;
				else {
                                    Individual prev = p.get(indices.get(previousMetaIndex));
                                    Individual next = p.get(indices.get(nextMetaIndex));
                                    localCrowdingDistance = Math.abs(next.getFitness(fitnessFuncName) - prev.getFitness(fitnessFuncName));
				}
				// update the crowding distance for this individual
				Individual current = p.get(indices.get(metaIndex));
				current.updateCrowdingDistance(localCrowdingDistance);
			}
		}
	}

        /**
         * 
         * @param min
         * @param max
         * @return 
         */
	public static ArrayList<Integer> getRange(Integer min, Integer max) {
            ArrayList<Integer> r = new ArrayList<Integer>();
            for (int i = min; i < max; i++)
                    r.add(i);
            return r;
	}

        /**
         * Comparator for ranking values
         */
	public static class IntegerComparator implements Comparator<Integer> {
		Population p;
		String fitnessFuncName;

                /**
                 * 
                 * @param _p
                 * @param _fitnessFuncName 
                 */
		public IntegerComparator(Population _p, String _fitnessFuncName) {
			p = _p;
			fitnessFuncName = _fitnessFuncName;
		}
		
                /**
                 * 
                 * @param index1
                 * @param index2
                 * @return 
                 */
		@Override
		public int compare(Integer index1, Integer index2) {
			// fetch corresponding individuals from Population
			Individual i1 = p.get(index1);
			Individual i2 = p.get(index2);
			Double f1 = i1.getFitness(fitnessFuncName);
			Double f2 = i2.getFitness(fitnessFuncName);
			return (f1 < f2) ? -1 : (f1 > f2) ? 1 : 0;
		}
	}
} public class DominatedCount extends Operator {

    /**
     * Performs a O(mn^2) non-dominated sort, where m is the number of objective
     * functions and n is the population size
     * 
     * @param p
     * @param f
     * @throws DominationException
     */
    public static void countDominated(Population p,LinkedHashMap<String, FitnessFunction> f) throws DominationException {
        for (int i = 0; i < p.size(); i++) {		
            Individual a = p.get(i);
            for (int j = i + 1; j < p.size(); j++) {
                Individual b = p.get(j);
                Boolean aDomb = domination(a, b, f);
                Boolean bDoma = domination(b, a, f);
                if (!aDomb && !bDoma && a.getFitnesses().equals(b.getFitnesses())) { // happens if a and b are identical
                    // in this case whoever happens to be located earlier in the population wins
                    if (i <= j) b.incrementDominationCount();
                    else a.incrementDominationCount();
                } else if (aDomb) {
                    b.incrementDominationCount();
                } else if (bDoma) {
                    a.incrementDominationCount();
                }
            }
        }
    }

    /**
     * Domination exception
     */
    public static class DominationException extends Exception {
        
        private static final long serialVersionUID = 243416542135464L;
        
        /**
         * 
         * @param s 
         */
        public DominationException(String s) {super(s);}
    }

    /**
     * The domination operator returns true if A dominates B, which occurs when
     * (1) A is at least as good as B in all objectives, and (2) A is strictly
     * better than B in at least one objective
     * 
     * @param a
     * @param b
     * @param fitnessFunctions
     * @return 
     * @throws evogpj.sort.DominatedCount.DominationException
     */
    @SuppressWarnings("static-access")
    public static boolean domination(Individual a, Individual b, LinkedHashMap<String, FitnessFunction> fitnessFunctions )throws DominationException {
        Set<String> aFuncNames = a.getFitnessNames();
        Set<String> bFuncNames = b.getFitnessNames();
        if (aFuncNames.size() != bFuncNames.size())
            throw new DominationException(("Error: individuals' fitnesses are not of same " + "length: a=\"%s\" and b=\"%s\"").format(a.toString(),b.toString()));
        Boolean strictlyBetter = false;
        for (String fitnessFunctionName : aFuncNames) {
            Double af = FitnessComparisonStandardizer.getFitnessForMinimization(a, fitnessFunctionName, fitnessFunctions);
            Double bf = FitnessComparisonStandardizer.getFitnessForMinimization(b, fitnessFunctionName, fitnessFunctions);
            if (!(af <= bf)) // domination violated
                return false;
            if (af < bf)
                strictlyBetter = true;
        }
        return strictlyBetter;
    }
} public static class DominationException extends Exception {
        
    private static final long serialVersionUID = 243416542135464L;
        
    /**
     * 
     * @param s 
     */
    public DominationException(String s) {super(s);}
}
import evogpj.evaluation.java.CSVDataJava;
import evogpj.evaluation.java.DataJava;
import evogpj.evaluation.FitnessFunction;
import evogpj.evaluation.java.SubtreeComplexityFitness;
import evogpj.evaluation.java.SRLARSJava;
import evogpj.genotype.TreeGenerator;
import evogpj.gp.GPException;
import evogpj.gp.Individual;
import evogpj.gp.MersenneTwisterFast;
import evogpj.gp.Population;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Properties;
import evogpj.operator.Crossover;
import evogpj.operator.CrowdedTournamentSelection;
import evogpj.operator.Initialize;
import evogpj.operator.Mutate;
import evogpj.operator.Select;
import evogpj.operator.SinglePointKozaCrossover;
import evogpj.operator.SinglePointUniformCrossover;
import evogpj.operator.SubtreeMutate;
import evogpj.operator.TournamentSelection;
import evogpj.operator.TreeInitialize;
import evogpj.sort.CrowdingSort;
import evogpj.sort.DominatedCount;
import evogpj.sort.DominatedCount.DominationException;
     * Create all the operators from the loaded params. Seed is the seed to use
     * for the rng. If specified, d_in is some DataJava to use. Otherwise, d_in
     * should be null and fitness will load in the appropriate data.
     * 
     * @param props
     * @param seed
     * @throws java.io.IOException
     * 
     */
    protected void create_operators(Properties props, long seed) throws IOException {
        System.out.println("Running evogpj with seed: " + seed);
        rand = new MersenneTwisterFast(seed);
        fitnessFunctions = splitFitnessOperators(FITNESS);
        for (String fitnessOperatorName : fitnessFunctions.keySet()) {
            if (fitnessOperatorName.equals(Parameters.Operators.SR_JAVA_FITNESS)) {
                DataJava data = new CSVDataJava(PROBLEM);
                minTarget = data.getTargetMin();
                maxTarget = data.getTargetMax();
                if (TERM_SET == null) {
                    TERM_SET = new ArrayList<String>();
                    for (int i = 0; i < data.getNumberOfFeatures(); i++) TERM_SET.add("X" + (i + 1));
                    System.out.println(TERM_SET);
                }
                fitnessFunctions.put(fitnessOperatorName,new SRLARSJava(data, MEAN_POW, COERCE_TO_INT,EXTERNAL_THREADS));
                //modelScalerJava = new SRModelScalerJava(data);
            } else if (fitnessOperatorName.equals(Parameters.Operators.SUBTREE_COMPLEXITY_FITNESS)) {
                fitnessFunctions.put(fitnessOperatorName,new SubtreeComplexityFitness());
            } else {
                System.err.format("Invalid fitness function %s specified for problem type %s%n",fitnessOperatorName);
            }
        }

        TreeGenerator treeGen = new TreeGenerator(rand, FUNC_SET, TERM_SET);
        if (INITIALIZE.equals(Parameters.Operators.TREE_INITIALIZE)) {
            initialize = new TreeInitialize(rand, props, treeGen);
        } else {
            System.err.format("Invalid initialize function %s specified%n",INITIALIZE);
            System.exit(-1);
        }

        // Set up operators.
        if (SELECT.equals(Parameters.Operators.TOURNEY_SELECT)) {
            select = new TournamentSelection(rand, props);
        } else if (SELECT.equals(Parameters.Operators.CROWD_SELECT)) {
            select = new CrowdedTournamentSelection(rand, props);
        } else {
            System.err.format("Invalid select function %s specified%n", SELECT);
            System.exit(-1);
        }

        mutate = new SubtreeMutate(rand, props, treeGen);
        //mutate = new SubtreeMutateConstants(rand, props, treeGen);

        if (XOVER.equals(Parameters.Operators.SPU_XOVER)) {
            xover = new SinglePointUniformCrossover(rand, props);
        } else if (XOVER.equals(Parameters.Operators.SPK_XOVER)) {
            xover = new SinglePointKozaCrossover(rand, props);
        } else {
            System.err.format("Invalid crossover function %s specified%n",XOVER);
            System.exit(-1);
        }

        //Tree t = treeGen.generateLinearModel(TERM_SET);
        //Individual linearModelInd = new Individual(t);
        // to set up equalization operator, we need to evaluate all the
        // individuals first
        pop = initialize.initialize(POP_SIZE);
        //pop.set(0,linearModelInd);
        // initialize totalPop to simply the initial population
        for (FitnessFunction f : fitnessFunctions.values())
            f.evalPop(pop);
        // calculate domination counts of initial population for tournament selection
        try {
            DominatedCount.countDominated(pop, fitnessFunctions);
        } catch (DominationException e) {
            System.exit(-1);
        }
        // save first front of initial population
        paretoFront = new Population();
        for (int index = 0; index < pop.size(); index++) {
            Individual individual = pop.get(index);
            if (individual.getDominationCount().equals(0))
                paretoFront.add(individual);
        }
        // calculate crowding distances of initial population for crowding sort
        if (SELECT.equals(Parameters.Operators.CROWD_SELECT)) {
            CrowdingSort.computeCrowdingDistances(pop, fitnessFunctions);
        }
    }

    /**
     * Accept potential migrants into the population
     * @param migrants
     */
    protected void acceptMigrants(Population migrants) {
            pop.addAll(migrants);
    }
	
    /**
     * This is the heart of the algorithm. This corresponds to running the
     * {@link #pop} forward one generation
     * <p>
     * Basically while we still need to produce offspring, we select an
     * individual (or two) as parent(s) and perform a genetic operator, chosen
     * at random according to the parameters, to apply to the parent(s) to
     * produce children. Then evaluate the fitness of the new child(ren) and if
     * they are accepted by the equalizer, add them to the next generation.
     * <p>
     * The application of operators is mutually exclusive. That is, for each
     * iteration of this algorithm, we will choose exactly one of crossover,
     * mutation and replication. However, which one we choose is determined by
     * sampling from the distribution specified by the mutation and crossover
     * rates.
     * 
     * @returns a LinkedHashMap mapping fitness function name to the best
     *          individual for that fitness function
     * @throws GPException
     *             if any of the operators receive a individual with an
     *             unexpected genotype, this is an error.
     */
    protected void step() throws GPException {