chemical-equation-0.0.1: src/LAPACK.hs

module Main where

import LAPACK.Common
import Common

import qualified Numeric.LAPACK.Linear.LowerUpper as LU
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Shape as ExtShape
import qualified Numeric.Netlib.Class as Class
import Numeric.LAPACK.Matrix ((|||), (##/#), (.*#))

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Shape ((::+)((::+)))

import Text.Printf (printf)

import qualified Options.Applicative as OP



nullspacePermutable ::
   (ExtShape.Permutable height, Eq height) =>
   (ExtShape.Permutable width, Eq width) =>
   (Class.Floating a) =>
   Matrix.Tall height width a ->
   Matrix.General height Matrix.ShapeInt a
nullspacePermutable mat =
   let lu = LU.fromMatrix mat

       n = Shape.size $ Matrix.height mat
       k = Shape.size $ Matrix.width mat

       lSplit =
         Matrix.mapHeight
            (const $ Matrix.width mat ::+ Matrix.shapeInt (n-k)) $
         Matrix.fromFull $ Matrix.toFull $ LU.extractL lu

       det = Triangular.determinant $ Triangular.takeUpper $ LU.extractU lu

       {- FixMe:
       Triangular.takeLower does not know,
       that the matrix has unit diagonal
       and it packs the triangle, which is unnecessary.
       -}
       sol =
         Matrix.negate (Matrix.takeBottom lSplit)
            ##/# Triangular.takeLower (toSquare (Matrix.takeTop lSplit))

   in LU.multiplyP LU.NonInverted lu $
      Matrix.mapHeight (const $ Matrix.height mat) $
      Matrix.transpose $
      det .*# (sol ||| Matrix.fromFull (Square.identityFromHeight sol))

nullspace ::
   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   Matrix.Tall height width a ->
   Matrix.General height Matrix.ShapeInt a
nullspace mat =
   Matrix.mapHeight ExtShape.deconsIntIndexed $
   nullspacePermutable $
   Matrix.mapHeight ExtShape.IntIndexed $
   Matrix.mapWidth ExtShape.IntIndexed $
   mat



main :: IO ()
main = do
   (verbosity,format,args) <- OP.execParser $ info parser
   (reactantSet,reactantsWithoutDuplicates)
      <- preprocessArguments verbosity args

   {-
   LU decomposition for wide matrix fails for C2H5OH CH4 O2 CO2 H2O.
   We must work on transposed matrix in order to benefit from row permutations.
   -}
   {-
   FixMe:
   HCl NaOH NaCl H2O
   has 4 reactants and 4 elements, though there is one solution,
   The matrix has rank 3.
   -}
   tall <-
      case Matrix.caseTallWide $ matrixFromReactantSet reactantSet of
         Right wide -> return $ Matrix.transpose wide
         Left tall ->
            fail $
            printf
               ("there must be more reactants than elements, "
                  ++ "but we have %d reactant(s) and %d element(s)")
               (Shape.size $ Matrix.width tall)
               (Shape.size $ Matrix.height tall)

   displaySolutions verbosity format reactantsWithoutDuplicates $
      map (cancelIntegerVector . fmap round . Array.toMap) $
      Matrix.toColumns $ nullspace tall