/****************************************************************************
 *
 * Copyright (c) 2009 Pengyi Yang & Liang Hsu 
 *     All rights reserved.  
 * 
 * Filename: BPSO.java
 * Abstract: A particle swarm optimization example which is used to find the maximum value of the function:
 *           f(x,y) = (1-x)*(1-x)*Math.exp(-x*x-(y+1)*(y+1))-(x-x*x*x-y*y*y)*Math.exp(-x*x-y*y)
 *           The maximum value is about 1.6705
 *           
 * Version:  1.0
 * 
 * Author:   Pengyi Yang            
 *          (mikesilver7@msn.com)  
 *          
 * Date:     2009/1/15
 *
 * Notes:  This program is free software; you can redistribute it and/or modify
 *         it under the terms of the GNU General Public License as published by
 *         the Free Software Foundation; either version 2 of the License, or
 *         (at your option) any later version.
 *
 *         This program is distributed in the hope that it will be useful,
 *         but WITHOUT ANY WARRANTY; without even the implied warranty of
 *         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *         GNU General Public License for more details.
 *
 *         You should have received a copy of the GNU General Public License
 *         along with this program; if not, write to the Free Software
 *         Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 * 
 * History: 
     <Date>     <Author>     <Modification> 
    2009/2/26   Liang Hsu     fix some bugs(|v|<Vmax)   
    2009/3/05   Pengyi Yang   modify the algorithm for imbalance data sampling
    2009/3/07   Liang Hsu     Adjust parameters
    2009/3/10   Liang Hsu     Add Mutation
    2009/3/10   Liang Hsu     Modified BPSO
 ****************************************************************************/

// output to test
import java.util.*;
import java.text.DecimalFormat;
import java.io.*;
import java.io.File;   
import java.io.FileOutputStream;   
import java.io.PrintStream; 

import weka.classifiers.Classifier;
import weka.classifiers.Evaluation;
import weka.classifiers.bayes.NaiveBayes;
import weka.classifiers.functions.Logistic;
import weka.classifiers.functions.RBFNetwork;
import weka.classifiers.functions.SMO;
import weka.classifiers.lazy.IBk;
import weka.classifiers.trees.J48;
import weka.classifiers.trees.RandomForest;
import weka.core.*;

public class BPSO
{
	// hold the input file name
	private String fileName;
	// hold input data set
	private Instances dataset;
	// hold the internal training set
	private Instances internalTrain;
	// hold the internal test set
	private Instances internalTest;
	// whether display detail
	private boolean detail;
	// hold control samples (samples with label "0")
	private String[] control;
	
	private int ITERATIONS = 150;
	private int POPULATIONSIZE = 100;
	// control size, will vary with different data sets
	private int CONTROLSIZE;
	
	
	// w is the inertia weight
	private final double w = 0.689343;

	
	// c1 and c2 are the cognitive and social acceleration constants respectively.
	private final double c1 = 1.42694; 
	private final double c2 = c1;


	// velocity bound  
	//0.9820 = 1/(1+exp(4))
	//0.0180 =  1/(1+exp(-4))
	private final double Vmax = 0.9820; 
	private final double Vmin = 0.0180; 
	
	//Particles
	private int[][] particles ;
	
	// The so far best position of each particle combined into 1-D fitness
	private double[] localBest;
	
	// The so far best position of the swarm combined into 1-D fitness 
	private double globalBest;
		
	// The current velocities associate with each particle
	private double[][] velocity;
	
	// The so far the best position of each particle (local optimal)
	private int[][] localBestOfParticles;

	//Mutation probability
	private final double mutation = 0.3; 
	
	// The so far the best position of the swarm (global optimal)
	private int[] globalBestOfParticles;;
	
	//Genotype 
	private double [][] genotype; 
	
	private Random random ;
	
	private double sum;
	
	// format output
	DecimalFormat dec; 
	
	public BPSO(String fileName, boolean detail)
	{
		this.fileName = fileName;
		this.detail = detail;
		this.random = new Random(System.currentTimeMillis());
		//this.ITERATIONS = this.ITERATIONS ;
		
		this.sum = 0.0;
		/* load in the training data set */
		try 
		{
			dataset = new Instances(new BufferedReader(new FileReader(this.fileName)));
			dataset.setClassIndex(dataset.numAttributes() - 1);
		}
		catch(IOException ioe)
		{
			ioe.printStackTrace();
		}
			
		// create a copy of original data set for folding
		Instances randData = new Instances(dataset);
		
		// dividing the data set with 3-fold stratify measure
		randData.stratify(3);
		
		
		for (int fold = 0; fold < 3; fold++) {
		
			// using the first 2 fold as internal training set. the last fold be the internal test set.
			internalTrain = randData.trainCV(3, fold);
			internalTest = randData.testCV(3, fold);
		
			//System.out.println("fold " + fold +":");
			//System.out.println(internalTest.toString());
			
			CONTROLSIZE = 0;
			
			// figure out how many control sample the internal training set has
			for (int i = 0; i < internalTrain.numInstances(); i++) 
			{
				if (internalTrain.instance(i).classValue() == 0)
				{
					CONTROLSIZE++;
				}
				
			}
		
			this.dec = new DecimalFormat("##.######");
			this.localBest = new double[POPULATIONSIZE];
			this.localBestOfParticles = new int[POPULATIONSIZE][CONTROLSIZE];
			this.globalBest = Double.MIN_VALUE;
			this.globalBestOfParticles = new int[CONTROLSIZE];
			this.velocity = new double[POPULATIONSIZE][CONTROLSIZE];
			this.particles = new int[POPULATIONSIZE][CONTROLSIZE];
			this.genotype = new double[POPULATIONSIZE][CONTROLSIZE];
		
			System.out.println("fold = " + fold);
			System.out.println("w = " + this.w);
			System.out.println("c1,c2 = " + this.c1);
			System.out.println("ITERATIONS = "  + this.ITERATIONS);
			System.out.println("POPULATIONSIZE = "  + this.POPULATIONSIZE);
		
			initialization();
			findMaxValue();
		
			saveResults(fold);
		}
	}
	
	// Initiate the particle positions and their velocities
	public void initialization()
	{	
		// print the initialization results
		System.out.println("initialization: ");
		
		// initiate particle position
		for (int x = 0; x < POPULATIONSIZE; x++)
		{
			for(int y = 0; y < CONTROLSIZE; y++)
			{
				if (random.nextDouble() > 0.5)
				{
					this.particles[x][y] = 1;
				}
				else
				{
					this.particles[x][y] = 0;
				}
				
				this.genotype[x][y] = 0.0;
			}
		}

						
		// initiate globalBest and localBest
		for(int count = 0; count < POPULATIONSIZE; count++)
		{
			this.localBest[count] = Double.MIN_VALUE;
		}
		
		this.globalBest = Double.MIN_VALUE;
		
		if (detail == true)
			printMatrix(particles);
	}
	
	// finding the best results with the given iterations
	public void findMaxValue()
	{	
		for (int iterator = 0; iterator < ITERATIONS; iterator++) 
		{
			
			System.out.println("iteration: " + iterator);
			
			/** 
			 * update the local bests and global best positions and their fitness 
			 */
			for (int x = 0; x < POPULATIONSIZE; x++) 
			{
				// evaluate the fitness of a given particle
				double fitness = evaluate(x);
				
				// update local values
				if (this.localBest[x] < fitness) 
				{
					for (int y = 0; y < this.CONTROLSIZE; y++)
					{
						this.localBestOfParticles[x][y] = this.particles[x][y];
					}				
					this.localBest[x] = fitness;
				}
				
				// update global values
				if (this.globalBest < fitness) 
				{
					for (int y = 0; y < this.CONTROLSIZE; y++)
					{
						this.globalBestOfParticles[y] = this.particles[x][y];
					}									
					this.globalBest = fitness;
				}
			}
			
			/**
			 * update each particle's velocity and position
			 */
			for (int x = 0; x < this.POPULATIONSIZE; x++) 
			{
				int postion = random.nextInt(CONTROLSIZE);
				for (int y = 0; y < this.CONTROLSIZE; y++) 
				{			
					// r1 and r2 are the uniform random numbers generated in the range of [0, 1]
					double r1 = random.nextDouble();
					double r2 = random.nextDouble();
					double r3 = random.nextDouble();
					
					// updating equations that govern the BPSO
					this.velocity[x][y] = w * this.velocity[x][y] + c1 * r1 * (this.localBestOfParticles[x][y] - this.particles[x][y]) + c2 * r2 * (this.globalBestOfParticles[y] - this.particles[x][y]);
					
					this.genotype[x][y] += this.velocity[x][y];
					/*if (this.velocity[x][y] > this.Vmax)
					{
						this.velocity[x][y] = this.Vmax;
					}
		
					if (this.velocity[x][y] < this.Vmin)
					{
						this.velocity[x][y] = this.Vmin;
					}*/
					
					if (this.genotype[x][y] > this.Vmax)
					{
						this.genotype[x][y] = this.Vmax;
					}
		
					if (this.genotype[x][y] < this.Vmin)
					{
						this.genotype[x][y] = this.Vmin;
					}
					
					//Mutation
					if (r3 < mutation && postion == y)
					{
						this.genotype[x][y] = -this.genotype[x][y];
					}
					this.particles[x][y] = simgiod(this.genotype[x][y]);
				}
			}
			
			if (detail == true)
				printMatrix(particles);
		}
	}


	public int simgiod(double velocity)
	{
		double  randomNumber = random.nextDouble();
		if (randomNumber >= 1/( 1+ Math.exp(-velocity)))