package networkinference;

/**
 *
 * @author Ilias Foudalis
 */


import java.io.*;
import com.jmatio.io.*;
import com.jmatio.types.MLArray;
import com.jmatio.types.MLDouble;
import java.util.*;

public class NetworkInference {
    /*
     * Latitude arrays for data sets of different resolutions.
     * When adding a new one also upddate the code on setCurrentLatitudes() function
     */

    double latitudes45[] = {88.5,84.5,80.5,76.5,72.5,68.5,64.5,60.5,56.5,52.5,48.5,44.5,40.5,36.5,32.5,28.5,24.5,20.5,16.5,12.5,8.5,4.5,0.5,-4.5,-8.5,-12.5,-16.5,-20.5,-24.5,-28.5,-32.5,-36.5,-40.5,-44.5,-48.5,-52.5,-56.5,-60.5,-64.5,-68.5,-72.5,-76.5,-80.5,-84.5,-88.5};
    double latitudes72[] = {88.75000,86.25000,83.75000,81.25000,78.75000,76.25000,73.75000,71.25000,68.75000,66.25000,63.75000,61.25000,58.75000,56.25000,53.75000,51.25000,48.75000,46.25000,43.75000,41.25000,38.75000,36.25000,33.75000,31.25000,28.75000,26.25000,23.75000,21.25000,18.75000,16.25000,13.75000,11.25000,8.750000,6.250000,3.750000,1.250000,-1.250000,-3.750000,-6.250000,-8.750000,-11.25000,-13.75000,-16.25000,-18.75000,-21.25000,-23.75000,-26.25000,-28.75000,-31.25000,-33.75000,-36.25000,-38.75000,-41.25000,-43.75000,-46.25000,-48.75000,-51.25000,-53.75000,-56.25000,-58.75000,-61.25000,-63.75000,-66.25000,-68.75000,-71.25000,-73.75000,-76.25000,-78.75000,-81.25000,-83.75000,-86.25000,-88.75000};
    double latitudes90[] = {-89,-87,-85,-83,-81,-79,-77,-75,-73,-71,-69,-67,-65,-63,-61,-59,-57,-55,-53,-51,-49,-47,-45,-43,-41,-39,-37,-35,-33,-31,-29,-27,-25,-23,-21,-19,-17,-15,-13,-11,-9,-7,-5,-3,-1,1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79,81,83,85,87,89};    
    double latitudes180[] = {89.50000,88.50000,87.50000,86.50000,85.50000,84.50000,83.50000,82.50000,81.50000,80.50000,79.50000,78.50000,77.50000,76.50000,75.50000,74.50000,73.50000,72.50000,71.50000,70.50000,69.50000,68.50000,67.50000,66.50000,65.50000,64.50000,63.50000,62.50000,61.50000,60.50000,59.50000,58.50000,57.50000,56.50000,55.50000,54.50000,53.50000,52.50000,51.50000,50.50000,49.50000,48.50000,47.50000,46.50000,45.50000,44.50000,43.50000,42.50000,41.50000,40.50000,39.50000,38.50000,37.50000,36.50000,35.50000,34.50000,33.50000,32.50000,31.50000,30.50000,29.50000,28.50000,27.50000,26.50000,25.50000,24.50000,23.50000,22.50000,21.50000,20.50000,19.50000,18.50000,17.50000,16.50000,15.50000,14.50000,13.50000,12.50000,11.50000,10.50000,9.500000,8.500000,7.500000,6.500000,5.500000,4.500000,3.500000,2.500000,1.500000,0.5000000,-0.5000000,-1.500000,-2.500000,-3.500000,-4.500000,-5.500000,-6.500000,-7.500000,-8.500000,-9.500000,-10.50000,-11.50000,-12.50000,-13.50000,-14.50000,-15.50000,-16.50000,-17.50000,-18.50000,-19.50000,-20.50000,-21.50000,-22.50000,-23.50000,-24.50000,-25.50000,-26.50000,-27.50000,-28.50000,-29.50000,-30.50000,-31.50000,-32.50000,-33.50000,-34.50000,-35.50000,-36.50000,-37.50000,-38.50000,-39.50000,-40.50000,-41.50000,-42.50000,-43.50000,-44.50000,-45.50000,-46.50000,-47.50000,-48.50000,-49.50000,-50.50000,-51.50000,-52.50000,-53.50000,-54.50000,-55.50000,-56.50000,-57.50000,-58.50000,-59.50000,-60.50000,-61.50000,-62.50000,-63.50000,-64.50000,-65.50000,-66.50000,-67.50000,-68.50000,-69.50000,-70.50000,-71.50000,-72.50000,-73.50000,-74.50000,-75.50000,-76.50000,-77.50000,-78.50000,-79.50000,-80.50000,-81.50000,-82.50000,-83.50000,-84.50000,-85.50000,-86.50000,-87.50000,-88.50000,-89.50000};
        
    double []currentLatitudes;

    

        
    double [][][] dataArray;
    Vector<Area> areas;
    int dimX, dimY, dimZ; //(lat x lon x time)
    
    PublicFunctions pf = new PublicFunctions();
    
    
    public static void main(String[] args) {        
        
        System.out.println("Usage: ");
        System.out.println("java -Xmx#g -jar NetworkInference.jar grided_timeseries.m area_file.txt ");
        String matFile = args[0];//"hadsst_1979_2005_win_noTre.mat";
        String areaFile = args[1];//"allAreas_merged_had_new.txt";
                
        NetworkInference ni = new NetworkInference(matFile, areaFile);
        
    }
    
  
    
    public NetworkInference(String matFile, String areaFile)
    {
                
        loadFileMatlab(matFile);
        loadAreaFile(areaFile);
                
        System.out.println("Loaded "+areas.size()+" areas.");
            
        setCurrentLatitudes();
        setTrueAreaSize();
        constructAreaCummAnomalyVectors();        
        constructLinks();
        
        
                        
        //and now the network is constructed!
        exportGraphCLM(areas,"",false);
    
        scoreDecomposition();                        
        printAllAreas();
        //exportGraphCLM(areas,"forLink",true);
        //and var metrics output here
        exportCorrelations();        
        areaSizeDistributions();
        exportLinkWeights();
        degreeVsStrength();
                                
        //links from the X strongest areas
        constructColorMapLinksStrength(areas,"",6);                                 
    }
    
    
    
    
    public void constructLinks()
    {
        //just making sure that no areas have 0 std (if a grid needs to be reversed such situations can occur)
       Vector<Area> toRemove = new Vector<Area>();
        for(int i = 0; i < areas.size(); i++)
        {
            Area alpha = areas.get(i);
            double powe = getVariance(alpha);
            if(powe == 0)
                toRemove.add(areas.get(i));
        }
        areas.removeAll(toRemove);
             
        //compute all edges;
        Vector<Edge> allPossibleEdges = new Vector<Edge>();
        for(int i = 0; i < areas.size(); i++)
        {
            Area fromArea = areas.get(i);
            for(int j = 0; j < areas.size(); j++)
            {
                Area toArea = areas.get(j);
                if(fromArea.areaID > toArea.areaID)
                {                    
                    double correl = (pearsonCorrel(toArea.cummulativeAnnomalies, fromArea.cummulativeAnnomalies));
                    Edge e = new Edge(fromArea.areaID,toArea.areaID,Math.abs(correl));
                    e.realCorrel = correl;
                    allPossibleEdges.add(e);                                        
                }
            }
        }
                
         //add the edges to their respective areas
         for(int i = 0; i < allPossibleEdges.size(); i++)
         {            
            Edge e = allPossibleEdges.get(i);          
            Area from = new Area(e.from);
            Area to = new Area(e.to);
            from = areas.get(areas.indexOf(from));
            to = areas.get(areas.indexOf(to));
            Edge e1 = new Edge(from.areaID,to.areaID,e.absCorrel);
            e1.realCorrel = e.realCorrel;        
            Edge e2 = new Edge(to.areaID,from.areaID,e.absCorrel);
            e2.realCorrel = e.realCorrel;            
            from.edges.add(e1);
            to.edges.add(e2);                                    
        }
        //now we need to compute the realweights
         
         //first of all set area power, power = log(sqrt(signal.power))
         //also compute totalPower
         //double totalPower = 0.0;
         
         for(int i = 0; i < areas.size(); i++)
         {
             double power = getVariance(areas.get(i));                                  
             power = Math.sqrt(power); 
             
             try
             {
               PrintWriter out = new PrintWriter(new FileWriter("stdAndSize.txt",true));
               out.println(power+" "+areas.get(i).nodes.size());
               out.close();
             }
             catch(IOException ioe)
             {}
             //totalPower+=power;             
             areas.get(i).power = power;            
         }         
         
         
         //now set link weights
         for(int i = 0; i < areas.size(); i++)
         {
             Area alpha = areas.get(i);
             for(int j = 0; j < alpha.edges.size(); j++)
             {
                 Edge e = alpha.edges.get(j);
                 Area beta = new Area(e.to);
                 beta = areas.get(areas.indexOf(beta));
                 //double realWeight = (alpha.power*beta.power*e.realCorrel)/totalPower;
                 //double realWeight = ((alpha.power+beta.power)/totalPower)*e.realCorrel;
                 double realWeight = (alpha.power*beta.power)*e.realCorrel;
                 areas.get(i).edges.get(j).realWeight = realWeight;
                 areas.get(i).edges.get(j).absWeight = Math.abs(realWeight);                         
             }
         }
         //finally set area strength
         for(int i = 0; i < areas.size(); i++)
         {
             Area alpha = areas.get(i);
             double strength = 0.0;
             for(int j = 0; j < alpha.edges.size(); j++)
                 strength+=alpha.edges.get(j).absWeight;
             areas.get(i).strength = strength;
         }
         //and print a color map
         try
        {
          PrintWriter out = new PrintWriter(new FileWriter("map_strength.txt"));
          double map[][] = new double[dimX][dimY];
          for(int i = 0; i < areas.size(); i++)
          {
              Area alpha = areas.get(i);
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);
                  map[pos[0]][pos[1]] = alpha.strength;
              }
          }
          for(int i = 0; i < dimX; i++)
          {
              for(int j = 0; j < dimY; j++)
                  out.print(map[i][j]+" ");
              out.println();
          }
          out.close();
        }
        catch(IOException ioe)
        {}
    }
    
    
    
    
    
    public double getVariance(Area alpha)
    {
        double signPower = 0.0;
        for(int i = 0; i < alpha.cummulativeAnnomalies.length; i++)
            signPower += Math.pow(alpha.cummulativeAnnomalies[i],2);
        signPower=signPower/(double)(dimZ-1);
        return signPower;
    }
    
    
    public void constructAreaCummAnomalyVectors()
    {
                //init vectors.
        for(int i = 0; i < areas.size(); i++)        
            areas.get(i).cummulativeAnnomalies = new double[dimZ];                        
        
        //now start constructing the vectors.
        for(int i = 0; i < areas.size(); i++)
        {
            Area alpha = areas.get(i);            
            for(int j = 0; j < alpha.nodes.size(); j++)
            {                
                //find the position of the cell
                int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);                                                               
                //add the data to my cummulative vector
                double[] timeSeries = dataArray[pos[0]][pos[1]];
                if(timeSeries[0] != 100 && timeSeries[0]!=-1000000)//mask values (e.g. for land)
                {
                    /*
                     * If one wants to compare data sets of different resolution
                     *the cumulative anomaly formula below needs to be modified
                     * to account for the "real" cell size
                     * Example:
                     * Convert cumulative anomalies of a 180x360 to a 90x144
                     * A cell in an 180x360 resolution has a size 1x1 (degrees)
                     * A cell in an 90x144 resolution has a suze 2x2.5 (degrees)
                     * Thus the cumulative anomaly of a180x360 cell needs to be 5 times smaller.
                     * To account of the latter the cumulative anomaly formula needs to be modified as:
                     * alpha.cummulativeAnnomalies[z] += (timeSeries[z]/5.0)*Math.abs(Math.cos(Math.toRadians(currentLatitudes[pos[0]])));
                     * 
                     */
                    
                    for(int z = 0; z < dimZ; z++)
                        alpha.cummulativeAnnomalies[z] += (timeSeries[z])*Math.abs(Math.cos(Math.toRadians(currentLatitudes[pos[0]])));
                }
            }            
        }    
    }
        
    
    public void loadAreaFile(String filename)
    {
        areas = new Vector<Area>();                  
        try
        {
            BufferedReader in = new BufferedReader(new FileReader(filename));
            String line = in.readLine();
            try
            {
                while(line!=null)
                {                                        
                    int areaID = Integer.parseInt(line.trim());                    
                    //get the cells
                    String cellIDs[] = in.readLine().split(" ");
                    //construct the area
                    Area alpha = new Area(areaID);
                    for(int i = 0; i < cellIDs.length; i++)
                        alpha.nodes.add(new Integer(Integer.parseInt(cellIDs[i])));
                    //add the area
                    areas.add(alpha);
                    //read new area id
                    line = in.readLine();
                }
            }
            catch(NullPointerException npe)
            {}
        }
        catch(IOException ioe)
        {
            System.out.println("Error while parsing the areas file.");
        }        
    }
    
    public void loadFileMatlab(String filename)
    {
        try
        {

            MatFileReader mat = new MatFileReader(filename);
            HashMap<String,MLArray> contents = new HashMap<String,MLArray>();
            contents = (HashMap)mat.getContent();//get all tha arrays of the mat file.
            //print their names.
            Set keys = contents.keySet();
            Object[] keynames = keys.toArray();
            for(int i = 0; i < keynames.length; i++)
                System.out.println(keynames[i].toString());
            //Now I just want the first array.
            boolean isdouble = (contents.get(keynames[0].toString())).isDouble();
            MLDouble mldouble = (MLDouble)(contents.get(keynames[0].toString()));
            double [][] temp = mldouble.getArray();
            
            int[] dimIds = mldouble.getDimensions();
            dimX = dimIds[0]; dimY = dimIds[1]; dimZ = dimIds[2];
            //reshape temp
            dataArray = new double[dimX][dimY][dimZ];
            
            for(int i = 0; i < dimX; i++)
            {
                int timePos = 0, lonPos = 0;;
                for(int j = 0; j < temp[i].length; j++)
                {
                    dataArray[i][lonPos][timePos] = temp[i][j];
                    lonPos++;
                    if(lonPos%(dimY)==0) //-1 cause have counts from 0.
                    {
                        lonPos = 0;
                        timePos = timePos+1;
                    }
                }
            }
            temp = null;
            System.gc();
        }
        catch(IOException io)
        {
            System.out.println("Error parsing .mat file.");
        }
    }
    
    
    
    public double pearsonCorrel(double[] scores1, double[] scores2)
    {
        double sum_sq_x = 0;
        double sum_sq_y = 0;
        double sum_coproduct = 0;
        double mean_x = scores1[0];
        double mean_y = scores2[0];
        for(int i=2;i<scores1.length+1;i+=1){
            double sweep =Double.valueOf(i-1)/i;
            double delta_x = scores1[i-1]-mean_x;
            double delta_y = scores2[i-1]-mean_y;
            sum_sq_x += delta_x * delta_x * sweep;
            sum_sq_y += delta_y * delta_y * sweep;
            sum_coproduct += delta_x * delta_y * sweep;
            mean_x += delta_x / i;
            mean_y += delta_y / i;
        }
        double pop_sd_x = (double) Math.sqrt(sum_sq_x/scores1.length);
        double pop_sd_y = (double) Math.sqrt(sum_sq_y/scores1.length);
        double cov_x_y = sum_coproduct / scores1.length;
        double result = cov_x_y / (pop_sd_x*pop_sd_y);
        return result;
    }
    
    public void markAreaCentersGCdist()
    {
        
        //the area center is defined as the node in the area
        //which minimized the great circle distance from all other nodes                        
        double longStep = 360.0/(double)dimY;
        for(int i = 0; i < areas.size(); i++)
        {
            Area alpha = areas.get(i);
            //quick check
            if(alpha.nodes.size() <= 2)
                alpha.centerNode = alpha.nodes.get(0);
            else
            {
                int bestCandidate = -1;
                double minDistanceSoFar = Double.MAX_VALUE;
                for(int x = 0; x < alpha.nodes.size(); x++)
                {
                    double currentDistance = 0.0;
                    int nodeFrom = alpha.nodes.get(x);
                    int []posFrom = pf.getPosition(nodeFrom+1, dimX, dimY);
                    double latFromDeg = currentLatitudes[posFrom[0]];
                    double lonFromDeg = posFrom[1]*longStep;
                    for(int y = 0; y < alpha.nodes.size(); y++)
                    {
                        int nodeTo = alpha.nodes.get(y);
                        if(nodeTo != nodeFrom)
                        {
                            int []posTo = pf.getPosition(nodeTo+1, dimX, dimY);
                            double latToDeg = currentLatitudes[posTo[0]];
                            double lonToDeg = posTo[1]*longStep;
                            currentDistance+=pf.distance3(Math.toRadians(latFromDeg), Math.toRadians(lonFromDeg), Math.toRadians(latToDeg), Math.toRadians(lonToDeg));                            
                        }                        
                    }
                    if(minDistanceSoFar > currentDistance)
                    {
                        minDistanceSoFar = currentDistance;
                        bestCandidate = alpha.nodes.get(x);
                    }
                }//
                areas.get(i).centerNode = bestCandidate;
                
            }
        }             
    }
    
    public void gcdLinkAndSize()
    {
        //a) sort the areas by size
        AreaSizeComparatorTrueSizeRev mycomp = new AreaSizeComparatorTrueSizeRev();
        EdgeStrengthComparator mycompEdge = new EdgeStrengthComparator();        
        java.util.Collections.sort(areas,mycomp);
        //gut, now for each area...
        try
        {
            PrintWriter out = new PrintWriter(new FileWriter("gcdSizeAndWeight.txt"));
            Vector<Edge> edges = new Vector<Edge>();
            for(int i = 0; i < areas.size(); i++)
            {
                Area alpha = areas.get(i);
                //get the links and put them to a vector
                edges.removeAllElements();
                edges.addAll(alpha.edges);
                java.util.Collections.sort(edges,mycompEdge); //sort the edges
                //for each edge print 
                // area size (from) (gcd) (lw)
                for(int j = 0; j < edges.size(); j++)
                {
                    Edge e = edges.get(j);
                    Area beta = new Area(e.to);
                    beta = areas.get(areas.indexOf(beta));
                    double gcd = 0;
                    int posFrom[] = pf.getPosition(alpha.centerNode+1, dimX, dimY);
                    int posTo[] = pf.getPosition(beta.centerNode+1, dimX, dimY);
                    double latFrom = Math.toRadians(currentLatitudes[posFrom[0]]);
                    double latTo = Math.toRadians(currentLatitudes[posTo[0]]);
                    double stepLon = 360.0/(double)dimY;
                    double lonTo = Math.toRadians(stepLon*posTo[1]);
                    double lonFrom  = Math.toRadians(stepLon*posFrom[1]);
                    gcd = pf.distance3(latFrom, lonFrom, latTo, lonTo);
                    out.println(i+" "+e.realWeight+" "+gcd);
                }
            
            }
            out.close();
        }
        catch(IOException ioe)
        {}
    }
    
    
    public void setTrueAreaSize()
    {
        Area alpha;
        double sum = 0.0;
        
        double scalingFactor = 0.0;
        for(int i = 0; i < dimX; i++)
            for(int j = 0; j < dimY; j++)
                scalingFactor += Math.cos(Math.toRadians(currentLatitudes[i]));
        scalingFactor = scalingFactor/10000;
          System.out.println(scalingFactor);
          
         
        
        for(int i = 0; i < areas.size(); i++)
        {
            alpha = areas.get(i);
            double areaSize = 0.0;
            for(int j = 0; j < alpha.nodes.size(); j++)
            {
                int pos[] = pf.getPosition(alpha.nodes.get(j).intValue()+1, dimX, dimY);
                double degreeLat = currentLatitudes[pos[0]];                
                areaSize+= Math.cos(Math.toRadians(degreeLat))*scalingFactor;
            }
            areas.get(areas.indexOf(alpha)).trueSize = areaSize;
            sum+= areaSize;
        }    
    }
    
    public void setCurrentLatitudes()
    {
        if(dimX == 72)
            currentLatitudes = latitudes72;
        else if(dimX == 180)
            currentLatitudes = latitudes180;
        else if(dimX == 90)
            currentLatitudes = latitudes90;
        else if(dimX == 45)
            currentLatitudes = latitudes45;
    }
    
    

    public void exportGraphCLM(Vector<Area> areas, String fileDelim, boolean normalize)
    {
        try
        {
            double maxLinkWeight = 1.0;//Double.MIN_VALUE;
            if(normalize)
            {
                maxLinkWeight = Double.MIN_VALUE;
                for(int i = 0; i < areas.size(); i++)
                {
                    Area alpha = areas.get(i);
                    for(int j = 0; j < alpha.edges.size(); j++)
                    {
                        double lw = alpha.edges.get(j).absWeight;
                        if(lw > maxLinkWeight)
                            maxLinkWeight = lw;
                    }
                }
            }
             
             
            
            
            PrintWriter out = new PrintWriter(new FileWriter("finalGraph"+fileDelim+".clm"));
            out.println("Vertices\t"+areas.size());
            int numOfEdges = 0;
            Area alpha;
            for(int i = 0; i < areas.size(); i++)
            {
                alpha = areas.get(i);
                numOfEdges+=alpha.edges.size();
                out.println(alpha.areaID+"\t");
                for(int j = 0; j < alpha.nodes.size(); j++)
                {
                    out.print(alpha.nodes.get(j)+" ");
                }
                out.print("\r\n");
            }//Finished with nodes.
            out.println("Edges\t"+numOfEdges);
            for(int i = 0; i < areas.size(); i++)
            {
                alpha = areas.get(i);
                for(int j = 0; j < alpha.edges.size(); j++)
                {
                    out.println(alpha.areaID+"\t"+alpha.edges.get(j).to+"\t"+(alpha.edges.get(j).realWeight/maxLinkWeight));
                }
            }
            out.close();
        }
        catch(IOException ioe)
        {
            ioe.printStackTrace();
        }
        //and now on .net
        try
        {
            for(int i = 0; i < areas.size(); i++)
            {
                Area alpha = areas.get(i);
                double power = getVariance(alpha);
                power = Math.sqrt(power);
                areas.get(i).power = power;
            }
            double totalPower = 0.0;
            for(int i = 0; i < areas.size(); i++)            
                totalPower+=areas.get(i).power;
            
            
            PrintWriter out = new PrintWriter(new FileWriter("finalNet"+fileDelim+".net"));
            java.util.HashMap<Integer,Integer> map = new java.util.HashMap();
            out.println("*Vertices "+areas.size());
            for(int i = 0; i < areas.size(); i++)
            {
                map.put(areas.get(i).areaID, (i+1));                
                out.println((i+1)+" "+areas.get(i).areaID);
            }
            out.println("*Arcs");
            //count how many edges
            int nedges = 0;
            for(int i = 0; i < areas.size(); i++)
                nedges+= areas.get(i).edges.size();
            nedges = nedges/2;//don't want to count twice
            out.println("*Edges "+nedges);
            Area from;
            for(int i = 0; i < areas.size(); i++)
            {
                from = areas.get(i);
                for(int j = 0; j <  from.edges.size(); j++)
                    if(map.get(from.areaID) < map.get(from.edges.get(j).to))
                        out.println(map.get(from.areaID)+" "+map.get(from.edges.get(j).to)+" "+Math.abs(from.edges.get(j).absWeight));                        
            }
            out.close();                        
        }
        catch(IOException ioe)
        {
            ioe.printStackTrace();
        }
    }
    
    /////////////////////////////////////
    ////METRICS AND STATS WILL GO HERE
    ////////////////////////////////////
    public void areaSizeDistributions()
    {
        try
        {
          PrintWriter outCellSize = new PrintWriter(new FileWriter("areaSize.txt"));
          PrintWriter outTrueSize = new PrintWriter(new FileWriter("trueAreaSize.txt"));
          Area alpha;
          for(int i = 0; i < areas.size(); i++)
          {
              alpha = areas.get(i);
              outCellSize.println(alpha.nodes.size());
              outTrueSize.println(alpha.trueSize);
          }
          outCellSize.close();
          outTrueSize.close();

        }
        catch(IOException ioe)
        {
            ioe.printStackTrace();
        }
    }
    
        /*
     * Output 
     1) all degrees
     2) all strengths
     3) degree of an area vs strength
     4) degree vs true area size
     4) strength vs true area size
     *
     * REQUIRED -> compute true area sizes first
     */
    public void degreeVsStrength()
    {
        try
        {                    
          PrintWriter outStr = new PrintWriter(new FileWriter("AreaStrengths.txt"));          
          PrintWriter outStrSize = new PrintWriter(new FileWriter("strengthVsTrueSize.txt"));
          double maxWeight = Double.MIN_VALUE,minWeight = Double.MAX_VALUE;
          
          Area alpha;
          for(int i = 0; i < areas.size(); i++)
          {
              alpha = areas.get(i);
              double areaDegree = alpha.edges.size();
              double areaStrength = 0.0;
              double size = alpha.trueSize;
              for(int j = 0; j < alpha.edges.size(); j++)
                  areaStrength+= Math.abs(alpha.edges.get(j).absWeight);
                            
              outStr.println(areaStrength);              
              outStrSize.println(areaStrength+" "+size);
          }          
          outStr.close();          
          outStrSize.close();
        }
        catch(IOException ioe)
        {
            ioe.printStackTrace();
        }
    }
    
    public void printAllAreas()
    {
        try
        {
          double[][] map = new double[dimX][dimY];
          PrintWriter out = new PrintWriter(new FileWriter("mapwithAllAreas.txt"));
          for(int i = 0; i < areas.size(); i++)
          {
              Area alpha = areas.get(i);
              double id = alpha.color;
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1,dimX,dimY);
                  map[pos[0]][pos[1]] = id;//alpha.areaID;
              }
          }
          for(int x = 0; x < dimX; x++)
          {
              for(int y = 0 ; y < dimY; y++)
                  out.print(map[x][y]+" ");
              out.println();
          }
          out.close();
        }
        catch(IOException ioe)
        {}
    }
    

    
    
    public void exportCorrelations()
    {
        try
        {
          PrintWriter out = new PrintWriter(new FileWriter("allCorrelations.txt"));
          for(int i = 0; i < areas.size(); i++)
          {
              Area alpha = areas.get(i);
              for(int j = 0; j < alpha.edges.size(); j++)
              {
                  Edge e = alpha.edges.get(j);
                  if(e.to > alpha.areaID)
                      out.println(e.realCorrel);
              }
          }
          out.close();
        }
        catch(IOException ioe)
        {}
    }

    
    public void exportLinkWeights()
    {
        try
        {
          PrintWriter out = new PrintWriter(new FileWriter("linkWeightsAreas.txt"));  
          double maxWeight = Double.MIN_VALUE,minWeight = Double.MAX_VALUE;
          
          Area alpha;
          
          for(int i = 0; i < areas.size(); i++)
          {
              alpha = areas.get(i);
              for(int j = 0; j < alpha.edges.size(); j++)
              {                  
                  Edge edge = alpha.edges.get(j);
                  if(edge.to < alpha.areaID)//do not count twice
                  {
                        double weight = Math.abs(edge.realWeight);
                        if(weight > maxWeight)
                            maxWeight = weight;
                        if(weight < minWeight)
                            minWeight = weight;
                        out.println(weight);              
                  }
              }
          }
          out.close();
          System.out.println("Max weight between areas: "+maxWeight);
          System.out.println("Min weight between areas: "+minWeight);
         
        }
        catch(IOException ioe)
        {
            ioe.printStackTrace();
        }
    }
    
    
    //What is this?
    //A color map depicting for a gaiven area A the areas who are connected to A.
    //The color of each area will be eqyal to the link weight from A to that area.
    //Which areas? the largest X
    public void constructColorMapLinks(int howmany)
    {
        AreaSizeComparatorTrueSize mycomp = new AreaSizeComparatorTrueSize();
        java.util.Collections.sort(areas,mycomp);        
        for(int i = 0; i < howmany; i++)
        {
                        
            Area alpha = areas.get(i);
            try
            {
              String map[][] = new String[dimX][dimY];  
              for(int x = 0; x < dimX; x++)
                  for(int y = 0; y < dimY; y++)
                      map[x][y]="NaN";
              
              PrintWriter out = new  PrintWriter(new FileWriter("linkColormap_area_"+alpha.areaID+".txt"));

              for(int j = 0; j < alpha.edges.size(); j++)
              {
                  Area temp = new Area(alpha.edges.get(j).to);
                  temp = areas.get(areas.indexOf(temp));
                  double weight = Math.abs(alpha.edges.get(j).realWeight);

                  boolean possign = true;
                  if(alpha.edges.get(j).realWeight < 0)
                      possign = false;
                  for(int x = 0; x < temp.nodes.size(); x++)
                  {
                      int pos[] = pf.getPosition(temp.nodes.get(x)+1, dimX, dimY);
                      if(possign)
                      {
                        map[pos[0]][pos[1]] = Double.toString(weight);
                      }
                      else
                      {
                        map[pos[0]][pos[1]] = Double.toString(-weight);                    
                      }
                  }
              }
              //now what -> give alpha a unique color and "print it too"
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);
                  map[pos[0]][pos[1]] = "NaN";                                    
              }
              //finished print the maps
              for(int x  = 0; x < dimX; x++)
              {
                  for(int y = 0; y < dimY; y++)                  
                      out.print(map[x][y]+" ");                                      
                  out.println();                  
              }
              out.close();              
            }
            catch(IOException ioe)
            {}
        }        
    }
    
    //What is this?
    //A color map depicting for a gaiven area A the areas who are connected to A.
    //The color of each area will be eqyal to the link weight from A to that area.
    //Which areas? the largest X
    public void constructColorMapLinksStrength(Vector<Area> areas, String delimiter, int howmany)
    {
        
        
        //need to compute the strength of each area (in abs. value)
        for(int i = 0; i < areas.size(); i++)
        {
            double strength = 0.0;
            Area alpha = areas.get(i);
            for(int j = 0; j < alpha.edges.size(); j++)
                strength += Math.abs(alpha.edges.get(j).absWeight);
            areas.get(i).strength = strength;
        }
        
        AreaStrengthComparator mycomp = new AreaStrengthComparator();
        java.util.Collections.sort(areas,mycomp);
        
        for(int i = 0; i < howmany; i++)
        {
                        
            Area alpha = areas.get(i);
            try
            {
              String map[][] = new String[dimX][dimY];  
              for(int x = 0; x < dimX; x++)
                  for(int y = 0; y < dimY; y++)
                      map[x][y]="NaN";              
              PrintWriter out = new  PrintWriter(new FileWriter(delimiter+"str_linkColormap_area_"+i+".txt"));              
              for(int j = 0; j < alpha.edges.size(); j++)
              {
                  Area temp = new Area(alpha.edges.get(j).to);
                  temp = areas.get(areas.indexOf(temp));
                  double weight = Math.abs(alpha.edges.get(j).realWeight);
                
                  boolean possign = true;
                  if(alpha.edges.get(j).realWeight < 0)
                      possign = false;
                  for(int x = 0; x < temp.nodes.size(); x++)
                  {
                      int pos[] = pf.getPosition(temp.nodes.get(x)+1, dimX, dimY);
                      if(possign)                      
                        map[pos[0]][pos[1]] = Double.toString(weight);                                      
                      else                      
                        map[pos[0]][pos[1]] = Double.toString(-weight);                                                                    
                  }
              }
              //now what -> give alpha a unique color and "print it too"
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);
                  map[pos[0]][pos[1]] = "NaN";                                    
              }
              //finished print the maps
              for(int x  = 0; x < dimX; x++)
              {
                  for(int y = 0; y < dimY; y++)                  
                      out.print(map[x][y]+" ");                                      
                  out.println();                  
              }
              out.close();              
            }
            catch(IOException ioe)
            {}
        }        
    }
    
    //return the ids of the neighboring (geographically connected) areas to alpha.
    public Vector<Integer> findNeighbors(Area alpha)
    {
        Vector<Integer> neighbors = new Vector<Integer>();
        for(int i = 0; i < alpha.nodes.size(); i++)
        {
            Vector<Integer> neighCells = pf.getNeighbors(alpha.nodes.get(i), dimX, dimY);
            for(int j = 0; j < areas.size(); j++)
            {
                boolean foundit = false;
                Area tempArea = areas.get(j);
                if(tempArea.areaID != alpha.areaID );
                {
                    for(int z = 0; z < neighCells.size(); z++)
                    {
                        if(tempArea.nodes.contains(neighCells.get(z)))
                        {
                            foundit= true;
                            break;
                        }
                    }
                    if(foundit && !neighbors.contains(tempArea.areaID))
                        neighbors.add(new Integer(tempArea.areaID));
                }
            }            
        }
        return neighbors;
        
    }
           
    
    public void scoreDecomposition()
    {         
        //again construct a dummy linked list representing the graph
                //Algorithm
        //for a particular strength s
        //   1) Remove all nodes with strength < s
        //   2) Recompute remaining node strength
        //      3) If nodes are found with strength < s, repeat 1
       //       4) else
       //    5) Nodes will be put in s-core
        Hashtable<Integer,Integer> mapAreaToNodeID = new Hashtable<Integer,Integer>();
        Hashtable<Integer,Integer> mapNodeToAreaID = new Hashtable<Integer,Integer>();
        Hashtable[] linkedListGraph = new Hashtable[areas.size()];
        Vector<Node> nodes = new Vector<Node>();
        NodeStrengthComparator nodeStrComp = new NodeStrengthComparator();
        for(int i = 0; i < areas.size(); i++)
        {
            mapAreaToNodeID.put(areas.get(i).areaID, i);
            mapNodeToAreaID.put(i, areas.get(i).areaID);
        }
        //construct linked list graph
        for(int i = 0; i < linkedListGraph.length; i++)//init
            linkedListGraph[i] = new Hashtable<Integer,Double>();
        for(int i = 0; i < areas.size(); i ++)
        {
            Area alpha = areas.get(i);
            for(int j = 0; j < alpha.edges.size(); j++)
            {
                Edge e = alpha.edges.get(j);
                int from = mapAreaToNodeID.get(alpha.areaID);
                int to = mapAreaToNodeID.get(e.to);
                linkedListGraph[from].put(to, e.absWeight);
            }            
        }
        //costruct the nodes and their strength
        double totalStr = 0.0;
        for(int i = 0; i < linkedListGraph.length; i++)
        {
            double strength = 0.0;
            java.util.Enumeration<Double> weights = linkedListGraph[i].elements();
            while(weights.hasMoreElements())
                strength+=weights.nextElement();
            Node n = new Node(i);
            n.strength = strength;
            totalStr+=strength;
            nodes.add(n);                    
        }
        java.util.Collections.sort(nodes,nodeStrComp);
        
        
        
        
        //aaaand we are ready to start
        boolean updateScore = true;
        double currentScore = 0.0;
        
        
        while(true)
        {
            java.util.Collections.sort(nodes,nodeStrComp);
            //remove the node with the lowest strength            
            Node n = nodes.remove(0);
            //find the area and set the areas score            
            if(updateScore)
            {
               currentScore = n.strength;
               updateScore = false;
            }
            Area dummyArea = new Area(mapNodeToAreaID.get(n.nodeid));
            areas.get(areas.indexOf(dummyArea)).score = currentScore;
            if(nodes.isEmpty())
                break;
            //now recompute the strength of the nodes
            Integer nodeToRemove = n.nodeid;
            for(int i = 0; i < linkedListGraph.length; i++)
                linkedListGraph[i].remove(nodeToRemove);
            linkedListGraph[nodeToRemove].clear();
            for(int i = 0; i < linkedListGraph.length; i++)
            {
                if(linkedListGraph[i].size() > 0)
                {
                    java.util.Enumeration<Double> weights =  linkedListGraph[i].elements();
                    double strength = 0.0;
                    while(weights.hasMoreElements())
                        strength+=weights.nextElement();                    
                    Node dummy = new Node(i);
                    nodes.get(nodes.indexOf(dummy)).strength = strength;
                }
            }
            java.util.Collections.sort(nodes,nodeStrComp);
            
            double minStr = nodes.get(0).strength;
            if(minStr > currentScore)
                updateScore=true;

        }            
        
        
        //thats it, now print a color map
        try
        {
          PrintWriter out = new PrintWriter(new FileWriter("map_score.txt"));
          double map[][] = new double[dimX][dimY];
          for(int i = 0; i < areas.size(); i++)
          {
              Area alpha = areas.get(i);
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);
                  map[pos[0]][pos[1]] = alpha.score;
              }
          }
          for(int i = 0; i < dimX; i++)
          {
              for(int j = 0; j < dimY; j++)
                  out.print(map[i][j]+" ");
              out.println();
          }
          out.close();
          out = new PrintWriter(new FileWriter("map_score_ranked.txt"));
          map = new double[dimX][dimY];
          AreaScoreComparator areaScoreComp = new AreaScoreComparator();
          java.util.Collections.sort(areas,areaScoreComp);
          int rank = 1;
          for(int i = 0; i < areas.size(); i++)
          {
              Area alpha = areas.get(i);
              for(int j = 0; j < alpha.nodes.size(); j++)
              {
                  int pos[] = pf.getPosition(alpha.nodes.get(j)+1, dimX, dimY);
                  map[pos[0]][pos[1]] = rank;
              }
              //get the next area
              if(i < areas.size()-1)
              {
                  Area beta = areas.get(i+1);
                  if(beta.score < alpha.score)
                      rank++;
              }
          }
          for(int i = 0; i < dimX; i++)
          {
              for(int j = 0; j < dimY; j++)
                  out.print(map[i][j]+" ");
              out.println();
          }
          out.close();
        }
        catch(IOException ioe)
        {}
        
           
    }                                                        
}

class AreaStrengthComparator implements Comparator<Area>
{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.strength < beta.strength)
            return 1;
        else if(alpha.strength > beta.strength)
            return -1;
        else
            return 0;
                
    }
    
}

class EdgeStrengthComparator implements Comparator<Edge>
{
    public int compare(Edge alpha, Edge beta)
    {
        if(alpha.absWeight < beta.absWeight)
            return -1;
        else if(alpha.absWeight > beta.absWeight)
            return 1;
        else
            return 0;
                
    }
    
}
class AreaDegreeComparator implements Comparator<Area>{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.edges.size() > beta.edges.size())
            return 1;
        else if(alpha.edges.size() < beta.edges.size())
            return -1;
        else return 0;
    }
}


class AreaSizeComparatorTrueSize implements Comparator<Area>
{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.trueSize > beta.trueSize)
            return -1;
        else if(alpha.trueSize < beta.trueSize)
            return 1;
        else
            return 0;
    }
}

class AreaSizeComparatorTrueSizeRev implements Comparator<Area>
{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.trueSize > beta.trueSize)
            return 1;
        else if(alpha.trueSize < beta.trueSize)
            return -1;
        else
            return 0;
    }
}

class AreaSizeComparatorCells implements Comparator<Area>
{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.nodes.size() > beta.nodes.size())
            return -1;
        else if(alpha.nodes.size() < beta.nodes.size())
            return 1;
        else return 0;
                    
    }
}




class Node //used on sCore on correlations (it is actually an area)
{
    int nodeid;
    double strength;//of correlations
    
    public Node(int nodeid)
    {
        this.nodeid = nodeid;
    }
    
    public Node()
    {
        
    }

    @Override
    public int hashCode() {
        int hash = 7;
        hash = 43 * hash + this.nodeid;
        return hash;
    }
    
    @Override
    public boolean equals(Object o)
    {
        Node n = (Node)o;
        if(!(n instanceof Node))
            return false;
        else if(n.nodeid == this.nodeid)
            return true;
        else
            return false;
    }
}

//sort the dipoles (from max to min str)
class NodeStrengthComparator implements Comparator<Node>
{
    public int compare(Node alpha, Node beta)
    {
        if(alpha.strength < beta.strength)
            return  -1;
        else if(alpha.strength > beta.strength)
            return 1;
        else
            return 0;
    }
    
}

class AreaScoreComparator implements Comparator<Area>
{
    public int compare(Area alpha, Area beta)
    {
        if(alpha.score < beta.score)
            return 1;
        else if(alpha.score > beta.score)
            return -1;
        else return 0;
    }
}
