//******************************************************************************
   //
   // Title:       Force Field X.
   // Description: Force Field X - Software for Molecular Biophysics.
   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
   //
   // This file is part of Force Field X.
   //
   // Force Field X is free software; you can redistribute it and/or modify it
  // under the terms of the GNU General Public License version 3 as published by
  // the Free Software Foundation.
  //
  // Force Field X 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
  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
  // Place, Suite 330, Boston, MA 02111-1307 USA
  //
  // Linking this library statically or dynamically with other modules is making a
  // combined work based on this library. Thus, the terms and conditions of the
  // GNU General Public License cover the whole combination.
  //
  // As a special exception, the copyright holders of this library give you
  // permission to link this library with independent modules to produce an
  // executable, regardless of the license terms of these independent modules, and
  // to copy and distribute the resulting executable under terms of your choice,
  // provided that you also meet, for each linked independent module, the terms
  // and conditions of the license of that module. An independent module is a
  // module which is not derived from or based on this library. If you modify this
  // library, you may extend this exception to your version of the library, but
  // you are not obligated to do so. If you do not wish to do so, delete this
  // exception statement from your version.
  //
  //******************************************************************************
  package ffx.algorithms.commands;
  
  import ffx.algorithms.cli.AlgorithmsCommand;
  import ffx.crystal.Crystal;
  import ffx.numerics.Potential;
  import ffx.numerics.math.RunningStatistics;
  import ffx.potential.ForceFieldEnergy;
  import ffx.potential.MolecularAssembly;
  import ffx.potential.parsers.DistanceMatrixFilter;
  import ffx.potential.parsers.PDBFilter;
  import ffx.potential.parsers.SystemFilter;
  import ffx.potential.parsers.XYZFilter;
  import ffx.potential.utils.Clustering.Conformation;
  import ffx.utilities.FFXBinding;
  import org.apache.commons.io.FilenameUtils;
  import org.apache.commons.math3.ml.clustering.CentroidCluster;
  import picocli.CommandLine.Command;
  import picocli.CommandLine.Option;
  import picocli.CommandLine.Parameters;
  
  import java.io.File;
  import java.io.FileWriter;
  import java.io.IOException;
  import java.util.ArrayList;
  import java.util.Collections;
  import java.util.List;
  
  import static ffx.potential.utils.Clustering.analyzeClusters;
  import static ffx.potential.utils.Clustering.hierarchicalClustering;
  import static ffx.potential.utils.Clustering.iterativeClustering;
  import static ffx.potential.utils.Clustering.kMeansClustering;
  import static java.lang.String.format;
  import static org.apache.commons.math3.util.FastMath.floorDiv;
  
  /**
   * The Cluster script clusters structures utilizing RMSD.
   * <p>
   * The Iterative Hierarchical Clustering method was contributed by Yuya Kinoshita, Koki Nishimura and
   * Masatoshi Karashima from Takeda Pharmaceuticals.
   *
   * @author Aaron J. Nessler
   * @author Mallory R. Tollefson
   * @author Michael J. Schnieders
   * <br>
   * Usage:
   * <br>
   * ffxc Cluster [options] &lt;filename&gt;
   */
  @Command(description = " Cluster structures using an RMSD matrix.", name = "Cluster")
  public class Cluster extends AlgorithmsCommand {
  
    /**
     * -a or --algorithm Clustering algorithm to use.
     * Algorithm: Multi-K-Means++ (0), Hierarchical (1), or Iterative (2)
     */
    @Option(names = {"-a", "--algorithm"}, paramLabel = "0", defaultValue = "0",
        description = "Algorithm: Multi-K-Means++ (0), Hierarchical (1), or Iterative (2)")
    private int algorithm;
  
    /**
     * -t or --trials Number of trials for Multi-K-Means or Iterative.
     */
   @Option(names = {"-t", "--trials"}, paramLabel = "1000", defaultValue = "1000",
       description = "Number of trials for Multi-K-Means or iterative method.")
   private int trials;
 
   /**
    * --tol or --tolerance Tolerance to determine equivalent points (iterative method).
    */
   @Option(names = {"--tol", "--tolerance"}, paramLabel = "0.5", defaultValue = "0.5",
       description = "Tolerance to determine if points are equivalent (iterative method).")
   private double tolerance;
 
   /**
    * -k or --clusters Maximum number of kmeans clusters.
    */
   @Option(names = {"-k", "--clusters"}, paramLabel = "0", defaultValue = "0",
       description = "Maximum number of k-means clusters.")
   private int maxClusters;
 
   /**
    * --rs or --randomSeed Set the random seed for deterministic clustering.
    */
   @Option(names = {"--rs", "--randomSeed"}, paramLabel = "-1", defaultValue = "-1",
       description = "Set the random seed for deterministic clustering (-1 uses the current time).")
   private long randomSeed;
 
   /**
    * --td or --threshold
    * RMSD value for dividing clusters from a hierarchical tree.
    */
   @Option(names = {"--td", "--threshold"}, paramLabel = "2.0", defaultValue = "2.0",
       description = "The distance value where a hierarchical tree flattened into clusters.")
   private double threshold;
 
   /**
    * -w or --write Write out an archive of a representative structure from each cluster.
    */
   @Option(names = {"-w", "--write"}, paramLabel = "false", defaultValue = "false",
       description = "Write an archive containing a representative from each cluster.")
   private boolean write;
 
   /**
    * The final argument(s) should be one or two filenames.
    */
   @Parameters(arity = "1..2", paramLabel = "files",
       description = "The RMSD distance matrix and optionally the corresponding structure archive (required for --write).")
   private List<String> filenames = null;
 
   /**
    * Generated list of clusters.
    */
   private List<CentroidCluster<Conformation>> clusterList;
 
   /**
    * Cluster constructor.
    */
   public Cluster() {
     super();
   }
 
   /**
    * Cluster constructor.
    *
    * @param binding The Binding to use.
    */
   public Cluster(FFXBinding binding) {
     super(binding);
   }
 
   /**
    * Cluster constructor that sets the command line arguments.
    *
    * @param args Command line arguments.
    */
   public Cluster(String[] args) {
     super(args);
   }
 
   /**
    * Return the Clusters.
    *
    * @return Returns the generated clusters.
    */
   public List<CentroidCluster<Conformation>> getClusterList() {
     return clusterList;
   }
 
   /**
    * Execute the script.
    */
   @Override
   public Cluster run() {
 
     // Init the context and bind variables.
     if (!init()) {
       return this;
     }
 
     if (filenames == null || filenames.isEmpty()) {
       logger.info(helpString());
       return this;
     }
 
     List<double[]> distMatrix = new ArrayList<>();
 
     String filename = filenames.get(0);
 
     DistanceMatrixFilter distanceMatrixFilter = new DistanceMatrixFilter();
     RunningStatistics runningStatistics = distanceMatrixFilter.readDistanceMatrix(filename, distMatrix);
 
     logger.info(" RMSD Distance Matrix Statistics");
     logger.info(format(" RMSD Minimum:  %8.6f", runningStatistics.getMin()));
     logger.info(format(" RMSD Maximum:  %8.6f", runningStatistics.getMax()));
     logger.info(format(" RMSD Mean:     %8.6f", runningStatistics.getMean()));
     double variance = runningStatistics.getVariance();
     if (!Double.isNaN(variance)) {
       logger.info(format(" RMSD Variance: %8.6f\n", variance));
     } else {
       logger.info("");
     }
 
     if (runningStatistics == null) {
       logger.info(format(" The distance matrix %s could not be read in.", filename));
       return this;
     }
 
     // Either use Multi-K-Means++, iterative, or Hierarchical agglomerative clustering.
     switch (algorithm) {
       case 1:
         // Hierarchical clustering.
         logger.info(" Performing Hierarchical Clustering");
         logger.info(format("  Cluster separation threshold: %6.4f A.\n", threshold));
         clusterList = hierarchicalClustering(distMatrix, threshold);
         break;
       case 2:
         // Iterative Clustering (produced by researchers at Takeda)
         logger.info(" Performing Iterative Clustering");
         clusterList = iterativeClustering(distMatrix, trials, tolerance);
         break;
       default:
         // K-Means++ and Multi K-Means++.
         logger.info(" Performing K-Means++ Clustering");
 
         // Ensure cluster are within bounds
         int dim = distMatrix.size();
         if (maxClusters <= 0 || maxClusters >= dim) {
           int newSize = floorDiv(dim, 2);
           maxClusters = newSize;
         }
         logger.info(format("  Number of clusters: %d", maxClusters));
 
         // Array representing indices of structures that represent their cluster.
         if (randomSeed == -1) {
           randomSeed = System.nanoTime();
         } else {
           logger.info(format("  Random seed:        %d", randomSeed));
         }
         logger.info(format("  Number of trials:   %d\n", trials));
         clusterList = kMeansClustering(distMatrix, maxClusters, trials, randomSeed);
     }
 
     boolean verbose = true;
     List<Integer> representativeConformer = new ArrayList<>();
     analyzeClusters(clusterList, representativeConformer, verbose);
 
     if (write) {
       writeStructures(representativeConformer);
     }
 
     return this;
   }
 
   /**
    * Write out structures corresponding to the representative of each cluster.
    *
    * @param repStructs Array list for index of representative structures.
    */
   private void writeStructures(List<Integer> repStructs) {
     String saveName = null;
     File structureFile = null;
     if (filenames.size() < 2) {
       saveName = FilenameUtils.removeExtension(filenames.get(0)) + ".arc";
       structureFile = new File(saveName);
       if (structureFile.exists()) {
         structureFile = algorithmFunctions.versionFile(structureFile);
       } else {
         saveName = FilenameUtils.removeExtension(filenames.get(0)) + ".pdb";
         structureFile = new File(saveName);
         if (structureFile.exists()) {
           structureFile = algorithmFunctions.versionFile(structureFile);
         } else {
           logger.info(" Please supply an ARC or PDB file corresponding to the distance matrix.");
           return;
         }
       }
     }
 
     // Turn off nonbonded terms.
     System.setProperty("vdwterm", "false");
 
     if (structureFile == null) {
       saveName = filenames.get(1);
       structureFile = new File(algorithmFunctions.versionFile(saveName));
     }
     File saveFile = new File(algorithmFunctions.versionFile(saveName));
     logger.info(format("\n Saving representative cluster members to %s.", structureFile.toString()));
 
     MolecularAssembly[] molecularAssemblies = algorithmFunctions.openAll(saveName);
     activeAssembly = molecularAssemblies[0];
 
     SystemFilter systemFilter = algorithmFunctions.getFilter();
     PDBFilter pdbFilter = null;
     XYZFilter xyzFilter = null;
     if (systemFilter instanceof PDBFilter) {
       pdbFilter = new PDBFilter(structureFile, activeAssembly, activeAssembly.getForceField(),
           activeAssembly.getProperties());
     } else if (systemFilter instanceof XYZFilter) {
       xyzFilter = new XYZFilter(structureFile, activeAssembly, activeAssembly.getForceField(),
           activeAssembly.getProperties());
     } else {
       logger.info(" Unknown file type.");
       return;
     }
 
     logger.info("\n Conformations: ");
     int structNum = 0;
     do {
       if (repStructs.contains(structNum)) {
         int clusterNum = repStructs.indexOf(structNum);
         // Make sure the crystal is is updated.
         Crystal crystal = activeAssembly.getCrystal();
         ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
         forceFieldEnergy.setCrystal(crystal);
         if (crystal.aperiodic()) {
           logger.info(format(" Conformation %d, Cluster %d)", structNum + 1, clusterNum + 1));
         } else {
           logger.info(format(" Conformation %d, Cluster %d (%s)", structNum + 1, clusterNum + 1,
               crystal.getUnitCell().toShortString()));
         }
         if (systemFilter instanceof PDBFilter) {
           pdbFilter.writeFile(structureFile, true, false, false);
         } else if (systemFilter instanceof XYZFilter) {
           String[] message = new String[1];
           message[0] = format(" Conformation %d, Cluster %d ", structNum + 1, clusterNum + 1);
           xyzFilter.writeFile(structureFile, true, message);
         }
       }
       structNum++;
     } while (systemFilter.readNext(false, false));
 
     if (systemFilter instanceof PDBFilter) {
       try (FileWriter fw = new FileWriter(structureFile, true)) {
         fw.append("END\n");
       } catch (IOException e) {
         logger.warning("Error appending END to PDB file: " + e.getMessage());
       }
     }
 
     systemFilter.closeReader();
   }
 
   /**
    * {@inheritDoc}
    */
   @Override
   public List<Potential> getPotentials() {
     return Collections.emptyList();
   }
 }