package iterativeclustering;

/*
 * Main.java
 *
 * Created on 12. November 2007, 21:45
 *
 * To change this template, choose Tools | Template Manager
 * and open the template in the editor.
 */

import biocomp.moltools.util.DoubleArrays;
import biocomp.moltools.util.IntArrays;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.impl.DenseDoubleMatrix1D;
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.util.StringTokenizer;
import cern.colt.matrix.doublealgo.Sorting;        



/**
 *
 * @author jan-hendrikprinz
 */
public class Main {
    
    /** Creates a new instance of Main */
    public Main() {
    }
    
    /**
     * @param args the command line arguments
     */
    public static void main(String[] args) throws FileNotFoundException, IOException {
        // TODO code application logic here
        
        
        
        int dummy = 0, state = 0;
        int fileNr = 0, numOfFiles=40, tokCount = 0;
        String nextLine;
        StringTokenizer tok;
        
        int size = 100200;    
        int lagtime = 60;
        int lengthOfSingleTrajectory = 200;
        int subTrajectories = lengthOfSingleTrajectory - lagtime - 1;

        int numOfAllTrajectories = (size / lengthOfSingleTrajectory) * (subTrajectories);
        
        int[][][] discreteTrajTransition = new int[numOfFiles][numOfAllTrajectories][2];
        double[][] discreteTrajEnergy = new double[numOfFiles][numOfAllTrajectories];
        double[][] discreteTrajTemperature = new double[numOfFiles][numOfAllTrajectories];
        
        double[] factorF = new double[numOfFiles];
        double[] factorFnew = new double[numOfFiles];
        double[][] factorFold = new double[10][numOfFiles];
        
        String outputFolder = "results/6statemodel/0K-1000K/lagtime" + Integer.toString(lagtime) + "/";
        
        // Read all data and partition into single trajectorie with one transition 
        
        // Read all Temps
        
        FileReader inputFileTrajTemp = new FileReader( "files/temps" );
        BufferedReader inputStreamTemp = new BufferedReader(inputFileTrajTemp);
        
        nextLine = inputStreamTemp.readLine();
        
        System.out.println("READ TEMPERATURE TABLE");

        tok = new StringTokenizer(nextLine, " ");
        double[] fileTemperature = new double[tok.countTokens()];
        tokCount = 0;
        while (tok.hasMoreTokens()) {
            fileTemperature[tokCount++] = 1.0 / (Double.valueOf(tok.nextToken()) * 0.001989);
        }
        
        //DoubleArrays.print(fileTemperature);
        
        System.out.println("");
        System.out.println("FINISHED TEMPERATURE TABLE");
        
        inputStreamTemp.close();
        inputFileTrajTemp.close();
        
        FileReader inputFileFFactor = new FileReader( "files/ffactor" );
        BufferedReader inputStreamFFactor = new BufferedReader(inputFileFFactor);
        
        nextLine = inputStreamFFactor.readLine();
        
        System.out.println("READ F FACTORS TABLE");

        
        tok = new StringTokenizer(nextLine, ", ");
        tokCount = 0;
        while (tok.hasMoreTokens()) {
            factorF[tokCount++] = Double.valueOf(tok.nextToken());
        }
        
        //DoubleArrays.print(fileTemperature);
        
        System.out.println("");
        System.out.println("FINISHED F FACTORS TABLE");
        
        inputStreamTemp.close();
        inputFileTrajTemp.close();
        
        
        // readAllEnergies
        
        FileReader inputFileTrajEnergies = new FileReader( "files/etot.out" );
        BufferedReader inputStreamEnergies = new BufferedReader(inputFileTrajEnergies);
        
        double[][] allEnergies = new double[numOfFiles][size];
        
        System.out.println("START READING ENERGIES");
        
        for (int ii=0; ii < size; ii++) {
            nextLine = inputStreamEnergies.readLine();

            tok = new StringTokenizer(nextLine, " ");
            double[] arr = new double[tok.countTokens()];
          
            tokCount = 0;
            while (tok.hasMoreTokens()) {
                allEnergies[tokCount++][ii] = Double.valueOf(tok.nextToken());
            }
        }
        
        for (fileNr=0; fileNr < numOfFiles; fileNr++) {
            
            for (int traj=0; traj < size / lengthOfSingleTrajectory; traj++) {
                // Iterate over all trajectories in one file
                
                double meanEnergy = 0.0;
                for (int kk=0; kk <lengthOfSingleTrajectory; kk++) {
                    meanEnergy += allEnergies[fileNr][traj * lengthOfSingleTrajectory + kk];
                }
                
                meanEnergy = meanEnergy / ((double) lengthOfSingleTrajectory);
                
                for (int ii=0; ii < subTrajectories; ii++) {
                    discreteTrajEnergy[fileNr][traj * subTrajectories + ii] = meanEnergy;
                    //System.out.println("Temp : " + fileNr + " / Traj : " + (traj * subTrajectories + ii) + " /  Energy : " +meanEnergy);
                }
            }
        }
        
        System.out.println("FINISHED READING ENERGIES");
        
        // Read All Trajectories
        
        System.out.println("START READING TRAJECTORIES");
        
        double[][] dihedrals = new double[size][];

        int phiBins = 36;
        int psiBins = 36;
        
        double phi = 0.0, psi = 0.0;

        
        int[] discreteTrajectory = new int[size];
        
        for (fileNr=0; fileNr < numOfFiles; fileNr++) {
            
            System.out.println("READ TRAJECTORIES WITH TEMPERATURE " + fileNr + " of " + numOfFiles);
            
            int currentTraj = 0;
            // Iterate over all files
            FileReader inputFileTrajectory = new FileReader( "files/" + Integer.toString(fileNr) + ".tors" );
            BufferedReader inputStreamTrajectory = new BufferedReader(inputFileTrajectory);

            nextLine = inputStreamTrajectory.readLine();

            for (int ii=0; ii < size; ii++) {
                nextLine = inputStreamTrajectory.readLine();

                tok = new StringTokenizer(nextLine, " ");
                double[] arr = new double[tok.countTokens()];
                tokCount = 0;
                tok.nextToken();
                while (tok.hasMoreTokens()) {
                    arr[tokCount++] = Double.valueOf(tok.nextToken());
                }
                dihedrals[ii] = arr;
                
                /*

                state = ((int)Math.floor(((dihedrals[ii][2] + 180.0) / 360.0) * (double)phiBins) % phiBins) * psiBins + 
                        ((int)Math.floor(((dihedrals[ii][3] + 180.0) / 360.0) * (double)psiBins) % psiBins); 

                 */
                
                // Johns Manual State Decomposition
                
                phi = dihedrals[ii][2];
                psi = dihedrals[ii][3];
                
                if((phi > 117.0)||(phi < -105.0)) {
                    if((psi > 28.0)||(psi < -124.0)) {
                        state = 0;
                    }
                    else {
                        state = 2;
                    }
                }
                
                if((phi > -105.0)&&(phi < 0.0)) {
                    if((psi > 28.0)||(psi < -124.0)) {
                        state = 1;
                    }
                    else {
                        state = 3;
                    }
                }
                
                if((phi > 0.0)&&(phi < 117.0)) {
                    if((psi > 111.0)||(psi < -5.0)) {
                        state = 4;
                    }
                    else {
                        state = 5;
                    }
                        
                }
                 
                discreteTrajectory[ii] = state;
            }
            
            // chunk each trajectory into bits of the right length and process further
            
            for (int traj=0; traj < size / lengthOfSingleTrajectory; traj++) {
                // Iterate over all trajectories in one file
                /*
                double meanEnergy = 0.0;
                for (int kk=0; kk <lengthOfSingleTrajectory; kk++) {
                    meanEnergy += allEnergies[fileNr][traj * lengthOfSingleTrajectory + kk];
                }
                
                meanEnergy = meanEnergy / ((double) lengthOfSingleTrajectory);
                */
                for (int ii=0; ii < subTrajectories; ii++) {
                    discreteTrajTransition[fileNr][currentTraj][0] = discreteTrajectory[traj * lengthOfSingleTrajectory + ii];
                    discreteTrajTransition[fileNr][currentTraj][1] = discreteTrajectory[traj * lengthOfSingleTrajectory + ii + lagtime];
                    
                    // discreteTrajEnergy[fileNr][currentTraj] = meanEnergy;
                    // discreteTrajTemperature[fileNr][currentTraj] = fileTemperature[fileNr];
                    currentTraj++;
                    /*
                    if(ii < lengthOfSingleTrajectory - lagtime - 2) {
                        meanEnergy += allEnergies[fileNr][traj * lengthOfSingleTrajectory + lagtime + ii + 1]
                                - allEnergies[fileNr][traj * lengthOfSingleTrajectory + ii];
                    }
                     **/
                    
                }
            }
            
            inputFileTrajectory.close();
            inputStreamTrajectory.close();
            
            // next file ...
        }
        
        System.out.println("FINISHED READING TRAJECTORIES");

        // use Johns Idea to compute the transition matrix
        
        // Iterate for self-consistancy
        
        //int numOfStates = phiBins*psiBins;
        
        int numOfStates = 6;
        
        double beta = 1.0;
        
        double probability = 0.0;
       
        double[][] cMatrix = new double[numOfStates][numOfStates];
        double[][] cMatrixLog = new double[numOfStates][numOfStates];
        double[][][] cMatrixLogFull = new double[numOfStates][numOfStates][numOfFiles];
        
        double cMatrixNorm = 0.0;
        double cMatrixNormLog = 0.0;
        
        double[][] vknLogMatrix = new double[numOfFiles][numOfAllTrajectories];
        
        boolean[] stateExists = new boolean[numOfStates];
        
        for(int kk = 0; kk < numOfFiles; kk++) {
            //System.out.println("Run : " + kk);
            for(int nn = 0; nn < numOfAllTrajectories; nn++) {
                
                double[] arguments2 = new double[numOfFiles];
                for(int ll = 0; ll < numOfFiles; ll++) {
                    arguments2[ll] = Math.log(numOfAllTrajectories) + factorF[ll] - fileTemperature[ll] * discreteTrajEnergy[kk][nn];
                }
                vknLogMatrix[kk][nn] = logSumExp(arguments2);
                //System.out.println(vknLogMatrix[kk][nn]);
                //vknLogMatrix[kk][nn] = 0.0;
            }
        }
        
        cMatrixNormLog = logSumExp(vknLogMatrix);

        // Scale f Factors so that the norm equals 1.0, that means the log is zero
        
        for(int kk = 0; kk < numOfFiles; kk++) {
            for(int nn = 0; nn < numOfAllTrajectories; nn++) {
                vknLogMatrix[kk][nn] = vknLogMatrix[kk][nn] - cMatrixNormLog;
            }
        }
        
        cMatrixNormLog = 0.0;

        System.out.println(cMatrixNormLog);
        
        // Compute Referenz maxlikelihood transition matrices
        
        for(int kk=0; kk<numOfFiles; kk++) {
            int[][] transitionCountMatrix = new int[numOfStates][numOfStates];
            double[][] tMatrixRef = new double[numOfStates][numOfStates];
            
            for(int nn = 0; nn < numOfAllTrajectories; nn++) {
                int st1 = discreteTrajTransition[kk][nn][0];
                int st2 = discreteTrajTransition[kk][nn][1];
                   
                transitionCountMatrix[st1][st2]++;
            }
                        
            
            for(int st1 = 0; st1 < numOfStates; st1++) {
                int rowSum = 0;
                for(int st2 = 0; st2 < numOfStates; st2++) {
                        rowSum += transitionCountMatrix[st1][st2];
                }
                if (rowSum > 0) {
                    for(int st2 = 0; st2 < numOfStates; st2++) {
                         tMatrixRef[st1][st2] = ((double)(transitionCountMatrix[st1][st2])) / ((double)rowSum);
                    }
                }
                else {    
                    for(int st2 = 0; st2 < numOfStates; st2++) {
                         tMatrixRef[st1][st2] = (1.0) / ((double)numOfStates);
                    }
                         
                } 
            }
            
            writeArrayToFile(outputFolder + "tMatrixRef" + Integer.toString(kk) + ".tMatrix", tMatrixRef);
                
        }
        
        double minTemp = 1.0;
        double maxTemp = 100.0;
        int stepTemp = 11;
        
        double[] tempTable = new double[stepTemp + 1];
        
        
        for(int temp=0; temp<stepTemp; temp+=1) {
            //beta = fileTemperature[temp];
            
            double actTemp = (((maxTemp - minTemp) * temp) / ((double)(stepTemp))) + minTemp;
            
            beta = 1.0 / (actTemp * 0.001989);
            
            
        
            // writeArrayToFile( "files/vknLogMatrix", vknLogMatrix);

            System.out.println("GET STATES & TRANSITIONS - temp : " + actTemp + "K  / beta = " + beta);

            boolean[][] doTransition = new boolean[numOfStates][numOfStates]; 
            int[][] argListCount = new int[numOfStates][numOfStates];
            
            for(int kk = 0; kk < numOfFiles; kk++) {
//                System.out.println(kk + " >>>> ");
                
                for(int nn = 0; nn < numOfAllTrajectories; nn++) {
                    int st1 = discreteTrajTransition[kk][nn][0];
                    int st2 = discreteTrajTransition[kk][nn][1];
                    argListCount[st1][st2]++;
                    argListCount[st2][st1]++;                
                    stateExists[st1] = true;
                    stateExists[st2] = true;
                    
                }
            }
            

//                System.out.println("GET STATES & TRANSITIONS");
                
                double[][][] argList = new double[numOfStates][numOfStates][];

                for(int st1 = 0; st1 < numOfStates; st1++) {
                    // System.out.println("Run : " + st1);
                    for(int st2 = 0; st2 < numOfStates; st2++) {
                        argList[st1][st2] = new double[argListCount[st1][st2]];
                        argListCount[st1][st2] = 0;
                    }
                }
            
//                System.out.println("GET ARGUMENTS");


                // System.out.println("Run : " + kk);
            for(int kk = 0; kk < numOfFiles; kk++) {
                for(int nn = 0; nn < numOfAllTrajectories; nn++) {
                    int st1 = discreteTrajTransition[kk][nn][0];
                    int st2 = discreteTrajTransition[kk][nn][1];
                    argList[st1][st2][argListCount[st1][st2]++] = - beta * discreteTrajEnergy[kk][nn] - vknLogMatrix[kk][nn] + Math.log(0.5);
                    argList[st2][st1][argListCount[st2][st1]++] = - beta * discreteTrajEnergy[kk][nn] - vknLogMatrix[kk][nn] + Math.log(0.5);
                }
            }
                
//                System.out.println("GET CMATRIX");
                /*
                for(int st1 = 0; st1 < numOfStates; st1++) {
                    for(int st2 = 0; st2 < numOfStates; st2++) {
                        if (argList[st1][st2].length > 0 ) {
                            doTransition[st1][st2] = true;
                            cMatrixLogFull[st1][st2][kk] = logSumExp(argList[st1][st2]);
                        }
                        else
                        {
                            cMatrixLogFull[st1][st2][kk] = - 1.0e100;
                        }
                    }
                }
                 */
            
                for(int st1 = 0; st1 < numOfStates; st1++) {
                    for(int st2 = 0; st2 < numOfStates; st2++) {
                        if (argList[st1][st2].length > 0 ) {
                            cMatrixLog[st1][st2] = logSumExp(argList[st1][st2]);
                        }
                        else
                        {
                            cMatrixLog[st1][st2] = - 1.0e100;
                        }
                    }
                }
            
/*
            for(int st1 = 0; st1 < numOfStates; st1++) {
                for(int st2 = 0; st2 < numOfStates; st2++) {
                    if (doTransition[st1][st2]) {
                        cMatrixLog[st1][st2] = logSumExp(cMatrixLogFull[st1][st2]);
                        for(int kk = 0; kk < numOfFiles; kk++) {
                            cMatrixLogFull[st1][st2][kk] = cMatrixLogFull[st1][st2][kk] - cMatrixLog[st1][st2];
                        }
                    }
                }
            }
 */
/*
            for(int st1 = 0; st1 < numOfStates; st1++) {
                for(int st2 = 0; st2 < numOfStates; st2++) {
                    if (argList[st1][st2].length > 0 ) {
                        //System.out.println("List for State : " + st1 + " / " + st2);
                        //DoubleArrays.print(argList[st1][st2]);
                        //System.out.println();
                        cMatrixLog[st1][st2] = logSumExp(argList[st1][st2]);
                    }
                }
            }
*/
//            writeArrayToFile( "files/CMatrixLogFull" + Integer.toString(temp) + ".cMatrixLogFull", cMatrixLogFull);

//            cMatrixLogFull = null;
            
            for(int st1 = 0; st1 < numOfStates; st1++) {
                for(int st2 = 0; st2 < numOfStates; st2++) {
                    cMatrix[st1][st2] = Math.exp(cMatrixLog[st1][st2]);
                    if (cMatrix[st1][st2] > 0.0) {
                        stateExists[st1] = true;
                    }
                }
            }

            writeArrayToFile( outputFolder + "CMatrix" + Integer.toString(temp) + ".cMatrix", cMatrix);

            System.out.println("REMOVE EMPTY STATES");
            
            // Remove Not Visited States

            int numOfVisitedStates = 0;
            int[] visitedStatesList = new int[numOfStates];

            int currentState = 0;

            for(int st1 = 0; st1 < numOfStates; st1++) {
                if(stateExists[st1]) {
                    numOfVisitedStates++;
                    visitedStatesList[st1] = currentState++;
                }
                else {
                    visitedStatesList[st1] = -1;
                }
            }

            double[] cMatrixVisitedRowSum = new double[numOfVisitedStates];
            double[][] cMatrixVisited = new double[numOfVisitedStates][numOfVisitedStates];

            double[] cMatrixLogVisitedRowSum = new double[numOfVisitedStates];
            double[][] cMatrixLogVisited = new double[numOfVisitedStates][numOfVisitedStates];


            double[][] tMatrixVisited = new double[numOfVisitedStates][numOfVisitedStates];

            for(int st1 = 0; st1 < numOfStates; st1++) {
                for(int st2 = 0; st2 < numOfStates; st2++) {
                    if (stateExists[st1]&&stateExists[st2]) {
                        cMatrixVisited[visitedStatesList[st1]][visitedStatesList[st2]] = cMatrix[st1][st2];
                        cMatrixLogVisited[visitedStatesList[st1]][visitedStatesList[st2]] = cMatrixLog[st1][st2];
                    }
                }
            }
            
            // BUILD TRANSITION MATRIX

            for(int st1 = 0; st1 < numOfVisitedStates; st1++) {
                cMatrixVisitedRowSum[st1] = 0.0;
                for(int st2 = 0; st2 < numOfVisitedStates; st2++) {
                    cMatrixVisitedRowSum[st1] += cMatrixVisited[st1][st2];
                }
                cMatrixLogVisitedRowSum[st1] = logSumExp(cMatrixLogVisited[st1]);
                
                // System.out.println("Log Matrix Sum : " + st1 + " is " + cMatrixLogVisitedRowSum[st1]);
                //if (cMatrixLogVisitedRowSum[st1] > -1000.0) {
                    for(int st2 = 0; st2 < numOfVisitedStates; st2++) {
                        tMatrixVisited[st1][st2] = Math.exp(cMatrixLogVisited[st1][st2] - cMatrixLogVisitedRowSum[st1]);
                    }
                //}
                //else
                //{
                //    stateExists[st1] = false;
                //}
            }
            
            /*
            
            System.out.println("REMOVE STATES WITH TOO LOW STATISTICS");
            
            int numOfVisitedStatesSmall = 0;
            
            int[] visitedStatesListSmall = new int[numOfStates];

            for(int st1 = 0; st1 < numOfStates; st1++) {
                if(stateExists[st1]) {
                    visitedStatesListSmall[st1] = numOfVisitedStatesSmall++;
                }
                else {
                    visitedStatesListSmall[st1] = -1;
                }
            }
            
            double[][] tMatrixVisitedSmall = new double[numOfVisitedStatesSmall][numOfVisitedStatesSmall];

            for(int st1 = 0; st1 < numOfStates; st1++) {
                for(int st2 = 0; st2 < numOfStates; st2++) {
                    if (stateExists[st1]&&stateExists[st2]) {
                        tMatrixVisitedSmall[visitedStatesListSmall[st1]][visitedStatesListSmall[st2]] = 
                                tMatrixVisited[visitedStatesList[st1]][visitedStatesList[st2]];
                    }
                }
            }
            */
            // Normalize again
            
            double avgStability = 0.0;
            
            for(int st1 = 0; st1 < numOfVisitedStates; st1++) {
                double rowSum = 0.0;
                for(int st2 = 0; st2 < numOfVisitedStates; st2++) {
                        rowSum += tMatrixVisited[st1][st2];
                }
                for(int st2 = 0; st2 < numOfVisitedStates; st2++) {
                         tMatrixVisited[st1][st2] /= rowSum;
                }
                avgStability += tMatrixVisited[st1][st1];
                
            }
            
            System.out.println("AVG Stability : " + avgStability / ((double)numOfVisitedStates));


            System.out.println(numOfVisitedStates);

            System.out.println(IntArrays.toString(visitedStatesList));

            // writeArrayToFile( "files/CMatrixVisited" + Integer.toString(temp) + ".cMatRed", cMatrixVisited);

            writeArrayToFile( outputFolder + "TMatrixVisited" + Integer.toString(temp) + ".tMatRed", tMatrixVisited);

            // DiscreteTrajectory myTraj = new DiscreteTrajectory(discreteTrajectory);

            // IntArrayDiscreteTrajectory myDiscreteTrajectory = new IntArrayDiscreteTrajectory(discreteTrajectory);        

            // IntArrayDiscreteTrajectory myDiscreteTrajectory = new IntArrayDiscreteTrajectory(new int[100]);        

            // Remove all unvisited States, might cause trouble with the eigenvalue decomposition otherwise !

            // myDiscreteTrajectory.removeEmptyStates();

            // TransitionMatrix myTransitionMatrix = new TransitionMatrix(myDiscreteTrajectory);

            /*
            TransitionMatrix myTransitionMatrix = new TransitionMatrix(tMatrixVisited, null);
            
            myTransitionMatrix.useLanczos = false;

            int numOfClusters = 1;

            double[] spec = myTransitionMatrix.spectrum(12);

            System.out.println(DoubleArrays.toString(spec));

            while ((spec[numOfClusters] > 0.95)||((spec[numOfClusters - 1] - spec[numOfClusters]) < 0.0175)) {
                numOfClusters++;
            }

            System.out.println(numOfClusters);

            int[] clustering = myTransitionMatrix.pcca(numOfClusters);

            FileWriter outputFileClustering = new FileWriter( outputFolder + "clusteringTemp" + Integer.toString(temp) + ".clustering" );
            BufferedWriter outputStreamClustering = new BufferedWriter(outputFileClustering);

            outputStreamClustering.write(IntArrays.toString(visitedStatesList)+ "\n");
            outputStreamClustering.write(IntArrays.toString(clustering) + "\n");

            outputStreamClustering.close();
            outputFileClustering.close();

            double[][] eigenVecs = myTransitionMatrix.eigenVectors(numOfClusters);

            writeArrayToFile(outputFolder + "eigenvecsTemp" + Integer.toString(temp) + ".ev", eigenVecs);
            */
        }
    }

    public static double logSumExp(double[] arguments) {
        DoubleMatrix1D myArgs = new DenseDoubleMatrix1D(arguments);
        DoubleMatrix1D myArgs2 = Sorting.quickSort.sort(myArgs);
        
        if (myArgs2.size() > 0) {
            double maxArg = myArgs2.get(myArgs2.size() - 1);
            
            //System.out.println(maxArg);
            
            //System.out.println(myArgs2.get(0));
            
            //System.out.println("logRowSum");
            //System.out.println(myArgs.toString());
        
            
            
            double expSum = 0.0;
        
            for (int ii = 0; ii < myArgs2.size() - 1; ii++) {
                expSum += Math.exp(myArgs2.get(ii) - maxArg);
            }
            
            //System.out.println(Math.log(expSum + 1.0) + maxArg);
        
            return Math.log(expSum + 1.0) + maxArg;
        }
        else {
            // minus infinity ???
            return -1.0e100;
        }
    }
    
    public static double logSumExp(double[][] arguments) {
        DoubleMatrix1D myArgs = new DenseDoubleMatrix1D(DoubleArrays.flatten(arguments));
        DoubleMatrix1D myArgs2 = Sorting.quickSort.sort(myArgs);
        
        double maxArg = myArgs2.get(myArgs2.size() - 1);
        
        // System.out.println(maxArg);
        
        double expSum = 0.0;
        
        for (int ii = 0; ii < myArgs2.size() - 1; ii++) {
            expSum += Math.exp(myArgs2.get(ii) - maxArg);
        }
        
        return Math.log(expSum + 1.0) + maxArg;
    }
    
    public static void writeArrayToFile(String fileName, int[][] array) throws IOException {
        System.out.print("WRITE " + fileName + "...");
        FileWriter outputFile = new FileWriter( fileName );
        BufferedWriter outputStream = new BufferedWriter(outputFile);

        outputStream.write(IntArrays.toString(array)+ "\n");

        outputStream.close();
        outputFile.close();
        System.out.println(" finished !");
    }
    
    public static void writeArrayToFile(String fileName, double[][] array) throws IOException {
        
        System.out.print("WRITE " + fileName + "...");
        FileWriter outputFile = new FileWriter( fileName );
        BufferedWriter outputStream = new BufferedWriter(outputFile);

        outputStream.write(DoubleArrays.toString(array)+ "\n");

        outputStream.close();
        outputFile.close();
        
        System.out.println(" finished !");
    }
    
    public static void writeArrayToFile(String fileName, double[][][] array) throws IOException {
        
        System.out.print("WRITE " + fileName + "...");
        FileWriter outputFile = new FileWriter( fileName );
        BufferedWriter outputStream = new BufferedWriter(outputFile);
        
        for(int ii=0;ii<array.length;ii++) {
            outputStream.write(DoubleArrays.toString(array[ii])+ "\n");
            outputStream.flush();
        }


        outputStream.close();
        outputFile.close();
        
        System.out.println(" finished !");
    }
    
}