-- |
-- Module      : Main
-- Description : Command Line Interface to analyse CREST results.
-- Copyright   : Phillip Seeber, 2021
-- License     : AGPL-3
-- Maintainer  : phillip.seeber@uni-jena.de
-- Stability   : experimental
-- Portability : POSIX, Windows
module Main (main) where

import ConClusion.Chemistry.Topology hiding (xyz)
import qualified ConClusion.Numeric.Statistics as Statistics
import Data.Aeson
import Data.Attoparsec.Text hiding (D)
import Data.Foldable
import qualified Data.IntSet as IntSet
import Data.Massiv.Array as Massiv hiding (B, D)
import qualified Data.Massiv.Array as Massiv
import Data.Text.Lazy (toStrict)
import Data.Text.Lazy.Builder as TB
import Data.Version (showVersion)
import Formatting hiding (char, (%))
import qualified Formatting as F
import Optics hiding (view)
import Paths_ConClusion (version)
import RIO hiding (Lens', Vector, lens, view, (^.))
import qualified RIO.Text as Text
import System.Console.CmdArgs
import System.IO.Unsafe (unsafePerformIO)

main :: IO ()
main = do
  app <- cmdArgs conClusionCmd >>= cmdToEnv
  runRIO app $ do
    logInfo $ "ConClusion, version " <> (displayShow . showVersion $ version) <> "\n"

    -- Processing.
    (featureMat, _trj) <- processTrajectory
    pcaMat <- dimReduction featureMat
    _clusters <- clusterFeatures pcaMat

    logInfo "Done. Good Luck with the results!"

----------------------------------------------------------------------------------------------------
-- Processing Logic

-- | Read the trajectory and calculate the specified features.
{-# SCC processTrajectory #-}
processTrajectory :: (HasXYZ env, HasDim env, HasLogFunc env) => RIO env (Matrix DL Double, Trajectory)
processTrajectory = do
  -- Obtain environment
  file <- view xyzL
  dims <- view dimL
  let outputFile = "features.dat"

  -- Logging
  {- ORMOLU_DISABLE -}
  logInfo $
    "+-----------------------+\n\
    \| Trajectory Processing |\n\
    \+-----------------------+\n\n\
    \  Trajectory File     : "  <> displayShow file <> "\n\
    \  Feature Matrix File : "  <> displayShow outputFile <> "\n\
    \  Features            :\n" <> displayFeatures dims
  {- ORMOLU_ENABLE -}

  -- Processing.
  raw <- readFileUtf8 file
  trj <- handleFailure $ parse' trajectory raw
  features <- handleFailure $ getFeatures dims trj

  -- Output.
  toGnuPlot Nothing (Just . Massiv.fromList Seq $ dims) (compute features)
    >>= writeFileUtf8 outputFile

  logInfo "\n\n\n"
  return (features, trj)

-- | Possibly perform a PCA.
{-# SCC dimReduction #-}
dimReduction :: (HasPrincipalComponentAnalysis env, HasLogFunc env) => Matrix DL Double -> RIO env (Matrix U Double)
dimReduction mat = do
  doPCA <- view pcaL
  case doPCA of
    Nothing -> return . compute $ featureMat
    Just PrincipalComponentAnalysis {keep} -> do
      let outputFile = "pca.dat"

      -- Logging.
      {- ORMOLU_DISABLE -}
      logInfo
        "+------------------------------+\n\
        \| Principal Component Analysis |\n\
        \+------------------------------+\n"
      {- ORMOLU_ENABLE -}

      -- Processing
      pcaData <- Statistics.pca keep featureMat
      let behaviourPercent = Statistics.remaining pcaData

      -- Logging.
      {- ORMOLU_DISABLE -}
      logInfo $
        "  Dimensions         : " <> display keep <> "\n\
        \  Eigenvalues        : " <> (displayEigenvalues . Statistics.allEigenValues $ pcaData) <> "\n\
        \  Mean Squared Error : " <> (displayMSE . Statistics.mse $ pcaData) <> "\n\
        \  Behaviour Captured : " <> displayPercent behaviourPercent <> "\n\
        \  PCA Matrix File    : " <> displayShow outputFile
      {- ORMOLU_ENABLE -}

      -- Output
      toGnuPlot Nothing Nothing (Statistics.y pcaData) >>= writeFileUtf8 outputFile

      logInfo "\n\n\n"
      return . Statistics.y $ pcaData
  where
    featureMat = compute @U mat

-- | Cluster analysis of the (dimensionalty reduced) feature matrix.
{-# SCC clusterFeatures #-}
clusterFeatures ::
  (HasClustering env, HasLogFunc env) => Matrix U Double -> RIO env Statistics.Clusters
clusterFeatures mat = do
  doClustering <- view clusteringL
  case doClustering of
    Nothing -> let Sz (_ :. n) = size mat in return $ makeArray @Massiv.B @Ix1 @IntSet Par (Sz n) $ \i -> IntSet.singleton i
    Just cl -> do
      let outputFile = "cluster.dat"
      -- Logging
      {- ORMOLU_DISABLE -}
      logInfo
        "+------------+\n\
        \| Clustering |\n\
        \+------------+\n"
      {- ORMOLU_ENABLE -}

      clusters <- case cl of
        DBScan {distance, measure, minSize} -> clusterDB distance measure minSize
        HCA {distance, measure, joinstrat, forcemin} -> clusterHCA distance measure joinstrat forcemin

      {- ORMOLU_DISABLE -}
      logInfo $
        "  Number of Clusters    : " <> (display . Massiv.elemsCount $ clusters) <> "\n\
        \  Cluster Data File     : " <> displayShow outputFile <> "\n\
        \  Cluster               : " <> displayClusters clusters
      {- ORMOLU_ENABLE -}

      -- Output
      toGnuPlot (Just clusters) Nothing mat >>= writeFileUtf8 outputFile

      return clusters
  where
    -- Make the distance function for clustering.
    distFn m = case m of
      Lr i -> Statistics.lpNorm i
      Manhattan -> Statistics.manhattan
      Euclidean -> Statistics.euclidean
      Mahalanobis -> Statistics.mahalanobis

    -- Cluster data with DBScan algorithm.
    clusterDB d m s = do
      clusters <- Statistics.dbscan (distFn m) s d mat
      {- ORMOLU_DISABLE -}
      logInfo $
        "  Algorithm            : DBScan\n\
        \  Distance Measure     : " <> displayShow m <> "\n\
        \  Search Distance ε    : " <> display d <> "\n\
        \  Minimal Cluster Size : " <> display s
      {- ORMOLU_ENABLE -}
      return clusters

    -- Cluster data with Hierarchical Cluster Analysis.
    clusterHCA d m j fm = do
      dendrogram <- Statistics.hca (distFn m) j mat
      let allClusters = Statistics.cutDendroAt dendrogram d
          dendroFile = "dendrogram.json"
          clusters = case fm of
            Nothing -> allClusters
            Just minSize -> compute . sfilter (\c -> IntSet.size c >= minSize) $ allClusters

      -- Write dendrogram file to disk.
      liftIO $ encodeFile dendroFile dendrogram

      {- ORMOLU_DISABLE -}
      logInfo $
        "  Algorithm             : Hierarchical Cluster Analysis\n\
        \  Distance Measure      : " <> displayShow m <> "\n\
        \  Search Distance       : " <> display d <> "\n\
        \  Cluster Join Strategy : " <> displayShow j <> "\n" <>
        (case fm of
          Nothing -> mempty
          Just minSize -> "  Minimal Cluster Size  : " <> display minSize <> "\n") <>
        "  Dendrogram File       : " <> displayShow dendroFile
      {- ORMOLU_ENABLE -}
      return clusters

----------------------------------------------------------------------------------------------------
-- Command line interfaces

-- | Command line arguments for ConClusion.
data ConClusionCmd = ConClusionCmd
  { -- | Filepath to the CREST xyz trajectory file.
    xyz :: FilePath,
    -- | Dimensions to process.
    dim :: String,
    -- | If PCA should be used to reduce dimensionalty of the problem. If so, the number of
    -- dimensions to keep.
    pca :: Maybe Int,
    -- | Cluster algorithm to be used for clustering, if any.
    cluster :: Maybe String,
    -- | Distance measure of points.
    measure :: String,
    -- | Cluster joining strategy for hierarchical clustering.
    joinstrat :: String,
    -- | Search distance in case of DBScan or distance for the Dendogram cut in case of hierarchical
    -- clustering.
    distance :: Double,
    -- | Mininmal size of a cluster in case of DBScan.
    minsize :: Int,
    -- | If to force a minimum size of clusters in HCA.
    forcemin :: Bool,
    -- | Do verbose logging.
    verbose :: Bool
  }
  deriving (Show, Data, Typeable)

conClusionCmd :: ConClusionCmd
conClusionCmd =
  ConClusionCmd
    { xyz =
        "crest_conformers.xyz"
          &= typFile
          &= help "Path to CREST trajectory.",
      dim =
        "e"
          &= help
            "Expression of dimensions to analyse:\n\
            \  - \"e\"         Energy\n\
            \  - \"b m n\"     Bond length between the atoms with\n\
            \                indices m and n\n\
            \  - \"a m n o\"   Angle between the atoms with indices\n\
            \                m, n and o\n\
            \  - \"d m n o p\" Sinus of the dihedral between the atoms\n\
            \                m, n, o and p\n\
            \Multiple features may be combined by giving a comma separated list of them.\n\
            \Default: e",
      pca =
        Nothing
          &= typ "INT"
          &= help
            "Perform a principal component analysis on the features.\
            \Optionally give the number of dimensions you want to keep.\n\
            \Default: 2"
          &= opt (2 :: Int),
      cluster =
        Nothing
          &= help
            "Cluster algorithm to cluster conformeres.\n\
            \  - \"dbscan\" DBScan algorithm. Depends on \n\
            \             \"measure\", \"distance\" and \"minSize\"\n\
            \  - \"hca\"    Hierarchical Cluster Analysis. Depends on\n\
            \             \"measure\", \"joinstrat\" and \"distance\"\n\
            \Default: dbscan"
          &= opt ("dbscan" :: String),
      measure =
        "mahalanobis"
          &= help
            "Distance measure for distance between conformeres.\n\
            \  - \"lrX\"         General L_r norm, where X is a value\n\
            \                  between 3 and 9\n\
            \  - \"manhattan\"   Manhattan distances\n\
            \  - \"euclidean\"   Euclidean distances\n\
            \  - \"mahalanobis\" Mahalanobis distances\n\
            \Default: mahalanobis",
      joinstrat =
        "ward"
          &= help
            "Strategy for inter-cluster distance in hierarchical clustering.\n\
            \  - \"single\"     Single-Linkage (minimal distance of\n\
            \                 points between clusters)\n\
            \  - \"complete\"   Complete-Linkage (maximal distance of\n\
            \                 points between clusters)\n\
            \  - \"median\"      Median distance linkage between\n\
            \                 all points between clusters)\n\
            \  - \"upgma\"      Average-Group-Linkage (average distance\n\
            \                 of all points of joined clusters)\n\
            \  - \"wpgma\"      Average-Linkage (average distance\n\
            \                 of all points between clusters)\n\
            \  - \"centroid\" Distance of cluster centroids\n\
            \  - \"ward\"     Ward distance of minimal residues\n\
            \Default: ward",
      distance =
        0.1
          &= typ "FLOAT"
          &= help
            "Search distance in DBScan or Dendogram cut distance in hierarchical clustering.\n\
            \Default: 0.1",
      minsize =
        5
          &= typ "INT"
          &= help
            "Minimal size of a cluster in DBScan.\n\
            \Default: 5",
      forcemin =
        False
          &= help
            "Force a minimal cluster size in hierarchical clustering.\n\
            \Default: False",
      verbose =
        False
          &= help "Verbose logging."
    }
    &= summary
      "Analysis of conformere ensembles.\n\
      \Processes an annotated XYZ trajectory (as generated by CREST, for example), where the energy\
      \is provided in the comment line. The analysis is a three-step procedure:\n\
      \  1. Calculate selected features from the ensemble, that are meant to\n\
      \     characterise the conformeres.\n\
      \     Relevant settings: \"--dim\"\n\
      \  2. (Optional) Perform a PCA analysis on the feature matrix to reduce the\n\
      \     number of dimensions.\n\
      \     Relevant settings: \"--pca\"\n\
      \  3. (Optional) Cluster the feature matrix into groups of conformeres. If a PCA\n\
      \     has been performed, use the feature matrix with reduced dimensionalty.\n\
      \     Relevant settins: \"--cluster\", \"--measure\", \"joinstrat\", \"--distance\",\n\
      \     \"--minsize\""

----------------------------------------------------------------------------------------------------
-- Runtime Environment.

-- Reader classes
class HasXYZ a where
  xyzL :: Lens' a FilePath

class HasDim a where
  dimL :: Lens' a [Feature]

class HasPrincipalComponentAnalysis a where
  pcaL :: Lens' a (Maybe PrincipalComponentAnalysis)

class HasClustering a where
  clusteringL :: Lens' a (Maybe Clustering)

-- | Full runtime environment of ConClusion with all settings combined.
data App = App
  { -- | Filepath to the CREST xyz trajectory file.
    xyz :: FilePath,
    -- | Features that are to be analysed.
    dim :: [Feature],
    -- | Settings for principal component analysis.
    pca :: Maybe PrincipalComponentAnalysis,
    -- | Settings for Clustering.
    clustering :: Maybe Clustering,
    -- | RIO's 'logFunc'.
    logFunc :: LogFunc
  }

-- Lenses
instance (k ~ A_Lens, a ~ FilePath, b ~ a) => LabelOptic "xyz" k App App a b where
  labelOptic = lens (xyz :: App -> FilePath) $ \s b -> (s {xyz = b} :: App)

instance (k ~ A_Lens, a ~ [Feature], b ~ a) => LabelOptic "dim" k App App a b where
  labelOptic = lens (dim :: App -> [Feature]) $ \s b -> (s {dim = b} :: App)

instance (k ~ A_Lens, a ~ Maybe PrincipalComponentAnalysis, b ~ a) => LabelOptic "pca" k App App a b where
  labelOptic = lens (pca :: App -> Maybe PrincipalComponentAnalysis) $ \s b -> (s {pca = b} :: App)

instance (k ~ A_Lens, a ~ Maybe Clustering, b ~ a) => LabelOptic "clustering" k App App a b where
  labelOptic = lens clustering $ \s b -> s {clustering = b}

instance (k ~ A_Lens, a ~ LogFunc, b ~ a) => LabelOptic "logFunc" k App App a b where
  labelOptic = lens logFunc $ \s b -> s {logFunc = b}

-- Reader Classes
instance HasXYZ App where
  xyzL = castOptic #xyz

instance HasDim App where
  dimL = castOptic #dim

instance HasPrincipalComponentAnalysis App where
  pcaL = castOptic #pca

instance HasClustering App where
  clusteringL = castOptic #clustering

instance HasLogFunc App where
  logFuncL = toLensVL #logFunc

-- | Settings for principal component analysis.
newtype PrincipalComponentAnalysis = PrincipalComponentAnalysis
  { -- | Number of dimensions to keep from PCA.
    keep :: Int
  }

-- Lenses
instance (k ~ A_Lens, a ~ Int, b ~ a) => LabelOptic "keep" k PrincipalComponentAnalysis PrincipalComponentAnalysis a b where
  labelOptic = lens keep $ \s b -> s {keep = b}

-- | Settings for clustering.
data Clustering
  = -- | DBScan.
    DBScan
      { -- | Search distance around a point.
        distance :: Double,
        -- | Distance measure in clustering.
        measure :: Measure,
        -- | Minimal size of a cluster.
        minSize :: Int
      }
  | -- | Hierarchical clustering.
    HCA
      { -- | Cut distance in the dendogram to select the number of clusters.
        distance :: Double,
        -- | Distance measure in clustering.
        measure :: Measure,
        -- | Join strategy for clusters.
        joinstrat :: Statistics.JoinStrat Double,
        -- | A minimum cluster size might be forced in post processing for HCA.
        forcemin :: Maybe Int
      }

-- | Distance measure between clusters.
data Measure
  = Lr Int
  | Manhattan
  | Euclidean
  | Mahalanobis
  deriving (Show)

----------------------------------------------------------------------------------------------------
-- Utilities

-- | Conversion of command line arguments into runtime environment.
cmdToEnv :: ConClusionCmd -> IO App
cmdToEnv ConClusionCmd {xyz, dim, pca, cluster, measure, joinstrat, distance, minsize, forcemin, verbose} = do
  -- Dimensionalty/Feature construction.
  dimApp <- case parseOnly dimParser (Text.pack dim) of
    Left err -> throwM $ ArgumentException "dim" err "a comma separated list of dimension specifications"
    Right res -> return res

  -- PCA construction.
  pcaSettings <- case pca of
    Nothing -> return Nothing
    Just x ->
      if x >= 1
        then return . Just $ PrincipalComponentAnalysis {keep = x}
        else throwM $ ArgumentException "pca" (show x) "a positive integer smaller than the number of clusters."

  -- Clustering construction.
  clDistance <-
    if distance <= 0
      then throwM $ ArgumentException "distance" (show distance) "a positive number"
      else return distance
  clMeasure <- case measure of
    "lr3" -> return $ Lr 3
    "lr4" -> return $ Lr 4
    "lr5" -> return $ Lr 5
    "lr6" -> return $ Lr 6
    "lr7" -> return $ Lr 7
    "lr8" -> return $ Lr 8
    "lr9" -> return $ Lr 9
    "manhattan" -> return Manhattan
    "euclidean" -> return Euclidean
    "mahalanobis" -> return Mahalanobis
    _ -> throwM $ ArgumentException "measure" measure "one of \"manhattan | euclidean | mahalanobis | lrX\" with X = [ 3 .. 9]"
  clMinsize <-
    if minsize <= 0
      then throwM $ ArgumentException "minsize" (show minsize) ">= 1"
      else return minsize
  clJoinstrat <- case joinstrat of
    "single" -> return Statistics.SingleLinkage
    "complete" -> return Statistics.CompleteLinkage
    "median" -> return Statistics.Median
    "upgma" -> return Statistics.UPGMA
    "wpgma" -> return Statistics.WPGMA
    "centroid" -> return Statistics.Centroid
    "ward" -> return Statistics.Ward
    _ -> throwM $ ArgumentException "joinstrat" joinstrat "one of \"single | complete | upgma | wpgma | centroid | ward\""
  clustering <- case cluster of
    Just "dbscan" ->
      return . Just $
        DBScan
          { distance = clDistance,
            measure = clMeasure,
            minSize = clMinsize
          }
    Just "hca" ->
      return . Just $
        HCA
          { distance = clDistance,
            measure = clMeasure,
            joinstrat = clJoinstrat,
            forcemin = if forcemin then Just clMinsize else Nothing
          }
    Nothing -> return Nothing
    Just c -> throwM $ ArgumentException "cluster" c "one of \"dbscan | hca\""

  -- Log function setup.
  logOptions <-
    id
      . setLogUseTime False
      . setLogUseColor True
      . setLogVerboseFormat False
      . setLogUseLoc False
      <$> logOptionsHandle stdout verbose
  (lf, _ :: IO ()) <- newLogFunc logOptions
  return
    App
      { xyz = xyz,
        dim = dimApp,
        pca = pcaSettings,
        clustering = clustering,
        logFunc = lf
      }

-- | Parser for dimensionalty expressions.
dimParser :: Parser [Feature]
dimParser = many1 $ do
  skipMany (char ' ')
  fType <- eParser <|> bParser <|> aParser <|> dParser
  skipMany (char ' ')
  _ <- option ',' (char ',')
  return fType
  where
    eParser = char 'e' *> return Energy
    bParser = do
      _ <- char 'b'
      skipMany (char ' ')
      a <- decimal
      skipMany (char ' ')
      b <- decimal
      return . Bond $ B a b
    aParser = do
      _ <- char 'a'
      skipMany (char ' ')
      a <- decimal
      skipMany (char ' ')
      b <- decimal
      skipMany (char ' ')
      c <- decimal
      return . Angle $ A a b c
    dParser = do
      _ <- char 'd'
      skipMany (char ' ')
      a <- decimal
      skipMany (char ' ')
      b <- decimal
      skipMany (char ' ')
      c <- decimal
      skipMany (char ' ')
      d <- decimal
      return . Dihedral $ D a b c d

-- | Substitution for RIO's @view@ for optics.
view :: (Is k A_Getter, MonadReader s f) => Optic' k is s b -> f b
view l = (^. l) <$> ask

-- | Parsers lifted to 'MonadThrow'.
parse' :: MonadThrow m => Parser x -> Text -> m x
parse' p t = case parseOnly p t of
  Left err -> throwM $ ParseException err
  Right res -> return res

-- | Logging a problem as an error and then rethrow the underlying exception.
handleFailure :: (Exception e, HasLogFunc env) => Either e r -> RIO env r
handleFailure (Right r) = return r
handleFailure (Left e) = do
  logError . displayShow $ e
  throwM e

-- | Display features for a trajectory in a nice way.
displayFeatures :: Foldable f => f Feature -> Utf8Builder
displayFeatures features =
  foldl
    ( \acc f ->
        acc
          <> ( case f of
                 Energy -> "    - E"
                 Bond (B a b) -> "    - B " <> display a <> " " <> display b
                 Angle (A a b c) -> "    - A " <> display a <> " " <> display b <> " " <> display c
                 Dihedral (D a b c d) -> "    - D " <> display a <> " " <> display b <> " " <> display c <> " " <> display d
             )
          <> "\n"
    )
    mempty
    features

-- | Display a percent value.
displayPercent :: Double -> Utf8Builder
displayPercent v = display $ sformat (fixed 1 F.% " %") v

-- | Display a mean squared error.
displayMSE :: Double -> Utf8Builder
displayMSE v = display $ sformat (fixed 2) v

-- | Display eigenvalues
displayEigenvalues :: Vector U Double -> Utf8Builder
displayEigenvalues vec = display . Massiv.fold . Massiv.map fmtD $ vec
  where
    fmtD d = format (right 6 ' ' F.%. fixed 2) d

-- | Display best cluster members.
displayBestMember :: Vector U Int -> Utf8Builder
displayBestMember vec =
  display . toLazyText . Massiv.fold . Massiv.map (\i -> bformat ("  - " F.% int F.% "\n") i) $ vec

-- | Converts a matrix to a Gnuplot matrix. Prints a header line as a comment with the features, if
-- available. Adds the assigned cluster number in the first column if available.
{-# SCC toGnuPlot #-}
toGnuPlot ::
  ( MonadThrow m
  ) =>
  Maybe Statistics.Clusters ->
  Maybe (Vector Massiv.B Feature) ->
  Matrix U Double ->
  m Text
toGnuPlot clusters features mat = do
  -- Build the lines representing the points. If cluster information is available, the first column
  -- is the cluster, that this point is assigned to.
  valueLines <- case clusters of
    Nothing -> return . Massiv.fold . Massiv.map (writeLine . Left) . innerSlices $ mat
    Just cl -> do
      let points = innerSlices mat :: Vector Massiv.D (Vector M Double)
          labeledCl = compute . Massiv.imap (,) $ cl :: Vector Massiv.B (Int, IntSet)

      -- Anotate all points with their assignment to a cluster number. If a point is from DBScan, it
      -- is not necessarily assigned to a cluster.
      annoPoints <- Massiv.iforM @Massiv.B points $ \pointIx vec -> do
        let filteredCl = Massiv.dropWhile (\(_, clPoints) -> not $ pointIx `IntSet.member` clPoints) labeledCl
        let assignedCl = filteredCl !? 0 :: Maybe (Int, IntSet)
            clNumber = fst <$> assignedCl :: Maybe Int
        return (clNumber, vec)

      -- return annoClusters
      return . Massiv.fold . Massiv.map (writeLine . Right) $ annoPoints

  featureLine <- case features of
    Nothing -> return mempty
    Just f -> do
      let Sz nFeatureLabels = actualFeatureSize f
          featureLine = writeFeatures f
      if nFeatures == nFeatureLabels
        then return featureLine
        else throwM $ SizeMismatchException (Sz nFeatureLabels) (Sz nFeatures)
  return . toStrict . TB.toLazyText $ featureLine <> valueLines
  where
    Sz (nFeatures :. _) = size mat

    -- The column width in the Gnuplot file.
    cw = 14

    -- Formatter for doubles.
    dForm a = bformat (left cw ' ' F.%. fixed 4) a

    -- Formatter for Integers.
    iForm = fitLeft 3 F.%. left 3 ' ' F.%. int

    -- Formatter for Features
    fForm f = case f of
      Energy -> bformat (left cw ' ' F.%. builder) "E"
      Bond (B a b) -> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm)) "B" a b
      Angle (A a b c) -> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A" a b c
      Dihedral (D a b c d) ->
        bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A sin" a b c d
          <> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A cos" a b c d

    -- Print a vector of numeric values.
    doubleVecB :: Source r Ix1 Double => Vector r Double -> TB.Builder
    doubleVecB v = Massiv.fold . Massiv.map dForm $ v

    -- Writer for a value line.
    writeLine ::
      (Source r Ix1 Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder
    writeLine (Left vec) = doubleVecB vec <> "\n"
    writeLine (Right (cl, vec)) = bformat iForm (fromMaybe (-1) cl) <> doubleVecB vec <> "\n"

    -- Writer for the Feature line.
    writeFeatures :: Vector Massiv.B Feature -> TB.Builder
    writeFeatures vec = "#" <> (Massiv.fold . Massiv.map fForm $ vec) <> "\n"

    -- Actual feature size.
    actualFeatureSize vec =
      unsafePerformIO $
        Massiv.foldlP
          ( \acc f -> case f of
              Dihedral _ -> 2 + acc
              _ -> 1 + acc
          )
          0
          (+)
          0
          vec

-- | Pretty printer for clusters.
displayClusters :: Statistics.Clusters -> Utf8Builder
displayClusters clusters = display . toLazyText . Massiv.fold . Massiv.map printCl $ clusters
  where
    -- Print clusters as list of indices.
    iForm = fitLeft 5 F.%. left 5 ' ' F.%. int

    -- Printer of a single cluster.
    printCl :: IntSet -> TB.Builder
    printCl cl =
      fst $
        IntSet.foldl
          ( \(bAcc, counter) i ->
              let this = bformat iForm i <> if counter `mod` (16 :: Int) == 0 then "\n      " else mempty
               in (bAcc <> this, counter + 1)
          )
          ("\n\n    - ", 1)
          cl

----------------------------------------------------------------------------------------------------
-- Exceptions

-- | Exception thrown on unexcpected command line arguments.
data ArgumentException = ArgumentException
  { argName :: String,
    is :: String,
    should :: String
  }

instance Show ArgumentException where
  show ArgumentException {argName, is, should} =
    "ArgumentException: Field \""
      <> argName
      <> "\" got \""
      <> is
      <> "\" but must be "
      <> should
      <> "."

instance Exception ArgumentException

-- | Exceptions thrown on failed parsing.
newtype ParseException = ParseException String

instance Show ParseException where
  show (ParseException err) = "ParseException: \"" <> err <> "\""

instance Exception ParseException