packages feed

chemical-equation (empty) → 0.0

raw patch · 13 files changed

+795/−0 lines, 13 filesdep +basedep +comfort-arraydep +containerssetup-changed

Dependencies added: base, comfort-array, containers, lapack, netlib-ffi, non-empty, numeric-quest, parsec, shell-utility, transformers, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2026, Henning Thielemann++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Henning Thielemann nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ chemical-equation.cabal view
@@ -0,0 +1,139 @@+Cabal-Version:       2.2+Name:                chemical-equation+Version:             0.0+Synopsis:            Balance chemical equations+Description:+  Balance coefficients of chemical equations.+  You only give the list of reactants,+  the program will automatically find out+  on which side of the equation the reactants are placed.+  However, the first reactant is always placed on the left side+  and the order of all reactants is preserved as good as possible.+  .+  > $ balance-chemical-equation C O2 CO2+  > C + O2 <=> CO2+  .+  > $ balance-chemical-equation CH4 O2 CO2 H2O+  > CH4 + 2 O2 <=> CO2 + 2 H2O+  .+  > $ balance-chemical-equation C8H18 O2 CO2 H2O+  > 2 C8H18 + 25 O2 <=> 16 CO2 + 18 H2O+  .+  > $ balance-chemical-equation N2 H2 NH3+  > N2 + 3 H2 <=> 2 NH3+  .+  > $ balance-chemical-equation HCl NaOH NaCl H2O+  > HCl + NaOH <=> NaCl + H2O+  .+  If you mix multiple reactions+  then the program takes them apart, again.+  However, the resulting partial reactions may not be chemically sensible.+  .+  > $ balance-chemical-equation C CH4 O2 CO2 H2O+  > 2 C + 2 H2O <=> CH4 + CO2+  > C + O2 <=> CO2+  .+  The program knows nothing about chemistry and has no dictionary of elements.+  It cannot check whether the reactants are existing chemical substances+  or whether the reactions can happen in the real world.+  It accepts any capital letter followed+  by many lower-case letters as element symbols.+  .+  Mathematically speaking, balancing a reaction equation+  means computing the nullspace of an integer matrix,+  where all coefficients in the nullspace vectors are in turn integers.+  The matrix has a row for every element and a column for every reactant+  and every cell contains+  the multiplicity of an element in the according reactant.+  Our solver transforms this matrix to Reduced Row Echelon form,+  a generalization of Gauss elimination.+Homepage:            https://hub.darcs.net/thielema/chemical-equation+License:             BSD-3-Clause+License-File:        LICENSE+Author:              Henning Thielemann+Maintainer:          haskell@henning-thielemann.de+Category:            Chemistry+Build-Type:          Simple++Flag flint+  Description: Implementation using FLINT library+  Manual:  True+  Default: False++Flag lapack+  Description: Enable experimental implementation using LAPACK+  Manual:  True+  Default: False++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/chemical-equation++Source-Repository head+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/chemical-equation++Executable balance-chemical-equation+  Build-Depends:+    numeric-quest >=0.2.1 && <0.3,+    shell-utility >=0.1 && <0.2,+    comfort-array >=0.5 && <0.6,+    parsec >=3.1 && <3.2,+    transformers >=0.3 && <0.7,+    containers >=0.4 && <0.9,+    non-empty >=0.3 && <0.4,+    utility-ht >=0.0.11 && <0.1,+    base >=4.5 && <5+  Default-Language:    Haskell2010+  GHC-Options:         -Wall+  Hs-Source-Dirs:      src+  Main-is:             Main.hs+  Other-Modules:+    Common++Executable balance-chemical-equation-flint+  If flag(flint)+    Build-Depends:+      shell-utility >=0.1 && <0.2,+      comfort-array >=0.5 && <0.6,+      parsec >=3.1 && <3.2,+      transformers >=0.3 && <0.7,+      containers >=0.4 && <0.9,+      non-empty >=0.3 && <0.4,+      utility-ht >=0.0.11 && <0.1,+      base >=4.5 && <5+  Else+    Buildable: False+  PkgConfig-Depends:   flint >=3 && <4+  Default-Language:    Haskell2010+  GHC-Options:         -Wall+  Hs-Source-Dirs:      src+  Main-is:             FLINT.hs+  Other-Modules:+    FLINT.Example+    FLINT.Utility+    FLINT.MatrixInteger+    FLINT.Integer+    FLINT.Type+    Common++Executable balance-chemical-equation-lapack+  If flag(lapack)+    Build-Depends:+      shell-utility >=0.1 && <0.2,+      lapack >=0.4 && <0.6,+      netlib-ffi >=0.1 && <0.2,+      comfort-array >=0.5 && <0.6,+      parsec >=3.1 && <3.2,+      containers >=0.4 && <0.9,+      non-empty >=0.3 && <0.4,+      utility-ht >=0.0.11 && <0.1,+      base >=4.5 && <5+  Else+    Buildable: False+  Default-Language:    Haskell2010+  GHC-Options:         -Wall+  Hs-Source-Dirs:      src+  Main-is:             LAPACK.hs+  Other-Modules:       LAPACK.Common, Common
+ src/Common.hs view
@@ -0,0 +1,138 @@+module Common where++import qualified Data.Traversable as Trav+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import qualified Data.Set as Set+import qualified Data.List.HT as ListHT+import qualified Data.List as List+import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import Data.Map (Map)+import Data.Set (Set)+import Data.Maybe (mapMaybe, isNothing)+import Data.Tuple.HT (mapPair, mapSnd, swap)++import qualified Text.Parsec as Parsec+import Text.Parsec.String (Parser)++import Text.Printf (printf)++import qualified Shell.Utility.Log as Log+import Shell.Utility.Verbosity (Verbosity)++import Control.Monad (when)+import Control.Applicative (liftA2, (<$>))++++type Reactant = Map String Integer+++parseReactant :: Parser Reactant+parseReactant =+   fmap (Map.fromListWith (+)) $+   Parsec.many1 $ do+      element <- liftA2 (:) Parsec.upper (Parsec.many Parsec.lower)+      kStr <- Parsec.many Parsec.digit+      case (kStr, read kStr) of+         ("", _) -> return (element, 1)+         {- ToDo:+         should be an unrecoverable parser error+         in order to not generate confusing additional errors,+         e.g. when parsing C0+         -}+         (_, 0) ->+            Parsec.unexpected $+            printf "element %s has multiplicity zero" element+         (_, k) -> return (element, k)++preprocessArguments ::+   Verbosity -> [String] -> IO (Set Reactant, [(Reactant, String)])+preprocessArguments verbosity args = do+   let (brokenArgs,reactants) =+         ListHT.unzipEithers $+         map (\arg -> flip (,) arg <$> Parsec.parse parseReactant arg arg)+            args+   case brokenArgs of+      _:_ -> fail $ concatMap ("\n\n"++) $ map show brokenArgs+      [] -> do+         let m = Map.fromListWith NonEmptyC.append $+                     map (mapSnd NonEmpty.singleton) reactants+         let duplicatesGroups =+               Map.filter+                  (\xs -> case xs of NonEmpty.Cons _ (_:_) -> True; _ -> False)+                  m+         Fold.for_ duplicatesGroups $ \duplicates ->+            Log.warn verbosity $+               printf "duplicate reactants: %s\n"+                  (List.intercalate ", " $ NonEmpty.flatten duplicates)++         return+            (Map.keysSet m,+               map snd $ filter fst $ snd $+               Trav.mapAccumL+                  (\seenSoFar p@(reactant,_name) ->+                     (Set.insert reactant seenSoFar,+                        (Set.notMember reactant seenSoFar, p)))+                  Set.empty reactants)++++mapFromListWithoutDuplicates :: (Ord k) => [(k,a)] -> Either [k] (Map k a)+mapFromListWithoutDuplicates xs =+   let m = Map.fromListWith (\_ _ -> Nothing) $ map (mapSnd Just) xs in+   case Trav.sequence m of+      Just table -> Right table+      Nothing -> Left $ Map.keys $ Map.filter isNothing m+++cancelIntegerVector :: (Functor f, Foldable f) => f Integer -> f Integer+cancelIntegerVector v =+   let d = Fold.foldl gcd 0 v+   in fmap (flip div d) v+++integerLinearCombination :: Map Reactant Integer -> Map String Integer+integerLinearCombination =+   Map.unionsWith (+) . Map.elems .+   Map.mapWithKey (\r k -> fmap (k*) r)++formatEquation :: [(Reactant,String)] -> Map Reactant Integer -> String+formatEquation rs reactantMap =+   let tagged =+         mapMaybe+            (\(r,name) ->+               let n = reactantMap Map.! r in+               case compare n 0 of+                  EQ -> Nothing+                  GT -> Just $ (True,(n,name))+                  LT -> Just $ (False,(-n,name)))+            rs+       (lhs,rhs) =+         (case tagged of ((True,_):_) -> id; _ -> swap) $+         mapPair (map snd, map snd) $+         List.partition fst tagged+       formatSum =+         List.intercalate " + " .+         map (\(n,str) -> if n==1 then str else printf "%d %s" n str)+   in  formatSum lhs ++ " <=> " ++ formatSum rhs++displaySolutions ::+   Verbosity -> [(Reactant, String)] -> [Map Reactant Integer] -> IO ()+displaySolutions verbosity reactants sols = do+   Fold.for_ sols $ \sol -> do+      putStrLn . formatEquation reactants $ sol+      let probe = integerLinearCombination sol+      when (Fold.any (0/=) probe) $+         Log.warn verbosity $ "probe failed\n" ++ show probe ++ "\n"++   let unused =+         if null sols+         then Set.fromList $ map fst reactants+         else Map.keysSet $+               Map.filter id $ Map.unionsWith (&&) $ map (fmap (0==)) sols+   when (not $ Set.null unused) $+      Log.warn verbosity $+      printf "Unused reactants: %s\n" $ List.intercalate ", " $ map snd $+      filter (flip Set.member unused . fst) $ reactants
+ src/FLINT.hs view
@@ -0,0 +1,70 @@+module Main where++import Common++import qualified FLINT.Utility as F+import qualified FLINT.MatrixInteger as FmpzMat+import qualified FLINT.Type as FLINT+import FLINT.Example ()++import Foreign.Ptr (Ptr)++import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Traversable as Trav+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import qualified Data.Set as Set+import Data.Set (Set)++import qualified Shell.Utility.Verbosity as Verbosity++import System.Environment (getArgs)++import Control.Monad.IO.Class (liftIO)+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)++++matrixFromReactantSet :: Set Reactant -> ContT r IO (Ptr FLINT.Matrix)+matrixFromReactantSet reactantSet = do+   let elements = foldMap Map.keysSet reactantSet+   mat <-+      ContT $+      F.withMatrix+         (fromIntegral $ Set.size elements)+         (fromIntegral $ Set.size reactantSet)+   liftIO $ do+      FmpzMat.zero mat+      Fold.for_ (zip [0..] $ Set.toAscList reactantSet) $ \(j,reactant) ->+         Fold.sequence_ $ flip Map.mapWithKey reactant $ \element n -> do+            let i = fromIntegral $ Shape.offset elements element+            ptr <- FmpzMat.entry mat i j+            F.integerFlintFromPrelude ptr n+   return mat+++{- |+The solution consists essentially of calling 'FmpzMat.nullspace'.+All the Haskell code is needed for transferring data+between FLINT and Haskell data structures,+parsing and validating arguments and presenting the results.+-}+main :: IO ()+main = do+   let verbosity = Verbosity.normal+   (reactantSet,reactantsWithoutDuplicates)+      <- preprocessArguments verbosity =<< getArgs++   sols <- evalContT $ do+      mat <- matrixFromReactantSet reactantSet+      n <- liftIO $ FmpzMat.ncols mat+      nullBasis <- ContT $ F.withMatrix n n+      nullspaceDim <- liftIO $ FmpzMat.nullspace nullBasis mat+      let reactantIndices = Map.fromList $ zip (Set.toAscList reactantSet) [0..]+      liftIO $+         Trav.for [0 .. nullspaceDim-1] $ \j ->+            Trav.for reactantIndices $ \i ->+               F.integerPreludeFromFlint =<< FmpzMat.entry nullBasis i j++   displaySolutions verbosity reactantsWithoutDuplicates $+      map cancelIntegerVector sols
+ src/FLINT/Example.hs view
@@ -0,0 +1,47 @@+module FLINT.Example where++import qualified FLINT.Utility as F+import qualified FLINT.MatrixInteger as FmpzMat++import qualified Data.Foldable as Fold+import Text.Printf (printf)+++example :: IO ()+example =+   let m = 3; n = 4 in+   F.withMatrix m n $ \mat ->+   F.withMatrix n n $ \null_basis ->++   do+   {- Set up example matrix:+   [1  2  3  4]+   [2  4  6  8]+   [3  6  9  12]+   (rank 1, nullspace dimension 3)+   -}+   FmpzMat.set_si mat 0 0 1+   FmpzMat.set_si mat 0 1 2+   FmpzMat.set_si mat 0 2 3+   FmpzMat.set_si mat 0 3 4++   FmpzMat.set_si mat 1 0 2+   FmpzMat.set_si mat 1 1 4+   FmpzMat.set_si mat 1 2 6+   FmpzMat.set_si mat 1 3 8++   FmpzMat.set_si mat 2 0 3+   FmpzMat.set_si mat 2 1 6+   FmpzMat.set_si mat 2 2 9+   FmpzMat.set_si mat 2 3 12++   nullspace_dim <- FmpzMat.nullspace null_basis mat++   printf "Matrix rank: %d\n" (toInteger (m + n - nullspace_dim))+   printf "Nullspace dimension: %d\n" (toInteger nullspace_dim)+   printf "\nNullspace basis vectors:\n"++   Fold.for_ [0 .. n-1] $ \i -> do+      Fold.for_ [0 .. nullspace_dim-1] $ \j ->+         printf "%s " =<< F.decimalFromInteger =<< FmpzMat.entry null_basis i j+      putStrLn ""
+ src/FLINT/Integer.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module FLINT.Integer where++import qualified FLINT.Type as FLINT++import qualified Foreign.C.Types as C+import Foreign.C.String (CString)+import Foreign.C.Types (CInt, CSize)+import Foreign.Ptr (Ptr)+++foreign import ccall "flint/fmpz.h fmpz_init"+   init :: Ptr FLINT.Integer -> IO ()++foreign import ccall "flint/fmpz.h fmpz_clear"+   clear :: Ptr FLINT.Integer -> IO ()++foreign import ccall "flint/fmpz.h fmpz_set"+   set :: Ptr FLINT.Integer -> Ptr FLINT.Integer -> IO ()++foreign import ccall "flint/fmpz.h fmpz_set_si"+   set_si :: Ptr FLINT.Integer -> C.CLong -> IO ()++foreign import ccall "flint/fmpz.h fmpz_print"+   print :: Ptr FLINT.Integer -> IO ()++foreign import ccall "flint/fmpz.h fmpz_sizeinbase"+   sizeinbase :: Ptr FLINT.Integer -> CInt -> IO CSize++foreign import ccall "flint/fmpz.h fmpz_get_str"+   get_str :: CString -> CInt -> Ptr FLINT.Integer -> IO CString++foreign import ccall "flint/fmpz.h fmpz_set_str"+   set_str :: Ptr FLINT.Integer -> CString -> CInt -> IO CInt++foreign import ccall "flint/fmpz.h fmpz_sgn"+   sgn :: Ptr FLINT.Integer -> IO CInt++foreign import ccall "flint/fmpz.h fmpz_abs"+   abs :: Ptr FLINT.Integer -> Ptr FLINT.Integer -> IO ()++foreign import ccall "flint/fmpz.h fmpz_add"+   add :: Ptr FLINT.Integer -> Ptr FLINT.Integer -> Ptr FLINT.Integer -> IO ()
+ src/FLINT/MatrixInteger.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module FLINT.MatrixInteger where++import qualified FLINT.Integer as Fmpz+import qualified FLINT.Type as FLINT++import qualified Foreign.C.Types as C+import Foreign.C.Types (CLong)+import Foreign.Ptr (Ptr)++++foreign import ccall "flint/fmpz_mat.h fmpz_mat_init"+   init :: Ptr FLINT.Matrix -> CLong -> CLong -> IO ()++foreign import ccall "flint/fmpz_mat.h fmpz_mat_clear"+   clear :: Ptr FLINT.Matrix -> IO ()++foreign import ccall "flint/fmpz_mat.h fmpz_mat_entry"+   entry :: Ptr FLINT.Matrix -> CLong -> CLong -> IO (Ptr FLINT.Integer)++foreign import ccall "flint/fmpz_mat.h fmpz_mat_nrows"+   nrows :: Ptr FLINT.Matrix -> IO C.CLong++foreign import ccall "flint/fmpz_mat.h fmpz_mat_ncols"+   ncols :: Ptr FLINT.Matrix -> IO C.CLong++foreign import ccall "flint/fmpz_mat.h fmpz_mat_zero"+   zero :: Ptr FLINT.Matrix -> IO ()++foreign import ccall "flint/fmpz_mat.h fmpz_mat_nullspace"+   nullspace :: Ptr FLINT.Matrix -> Ptr FLINT.Matrix -> IO CLong+++set_si :: Ptr FLINT.Matrix -> CLong -> CLong -> CLong -> IO ()+set_si mat i j x = do+   ptr <- entry mat i j+   Fmpz.set_si ptr x
+ src/FLINT/Type.hsc view
@@ -0,0 +1,16 @@+module FLINT.Type where++import Prelude hiding (Integer)+++#include "flint/fmpz.h"+#include "flint/fmpz_mat.h"++data Integer = Integer+data Matrix = Matrix++integerSize :: Int+integerSize = #{size fmpz_t}++matrixSize :: Int+matrixSize = #{size fmpz_mat_t}
+ src/FLINT/Utility.hs view
@@ -0,0 +1,72 @@+module FLINT.Utility where++import qualified FLINT.MatrixInteger as FmpzMat+import qualified FLINT.Integer as Fmpz+import qualified FLINT.Type as FLINT++import Text.Printf (printf)++import Control.Exception (bracket_)+import Control.Monad (void)+import Control.Applicative ((<$>))++import Foreign.Marshal.Alloc (allocaBytes)+import Foreign.C.String (peekCString, withCString)+import Foreign.C.Types (CLong)+import Foreign.Ptr (Ptr)+++withInteger :: (Ptr FLINT.Integer -> IO a) -> IO a+withInteger act =+   allocaBytes FLINT.integerSize $ \ptr ->+   bracket_ (Fmpz.init ptr) (Fmpz.clear ptr) (act ptr)++withMatrix :: CLong -> CLong -> (Ptr FLINT.Matrix -> IO a) -> IO a+withMatrix m n act =+   allocaBytes FLINT.matrixSize $ \mat ->+   bracket_ (FmpzMat.init mat m n) (FmpzMat.clear mat) (act mat)++decimalFromInteger :: Ptr FLINT.Integer -> IO String+decimalFromInteger ptr = do+   let base = 10+   num_digits <- Fmpz.sizeinbase ptr base+   allocaBytes (fromIntegral num_digits + 2) $ \buffer -> do+      void $ Fmpz.get_str buffer base ptr+      peekCString buffer++naturalPreludeFromFlint :: Ptr FLINT.Integer -> IO Integer+naturalPreludeFromFlint ptr = do+   let base = 16+   num_digits <- Fmpz.sizeinbase ptr base+   allocaBytes (fromIntegral num_digits + 2) $ \buffer -> do+      void $ Fmpz.get_str buffer base ptr+      read . ("0x"++) <$> peekCString buffer++integerPreludeFromFlint :: Ptr FLINT.Integer -> IO Integer+integerPreludeFromFlint ptr =+   withInteger $ \absInt -> do+   Fmpz.abs absInt ptr+   sgn <- Fmpz.sgn ptr+   let base = 16+   num_digits <- Fmpz.sizeinbase ptr base+   allocaBytes (fromIntegral num_digits + 2) $ \buffer -> do+      void $ Fmpz.get_str buffer base absInt+      (toInteger sgn *) . read . ("0x"++) <$> peekCString buffer++-- | return sign and turn integer into absolute value+splitSign :: Ptr FLINT.Integer -> IO Ordering+splitSign ptr = do+   sgn <- Fmpz.sgn ptr+   Fmpz.abs ptr ptr+   return $ compare sgn 0++-- | returns success+integerFromDecimal :: Ptr FLINT.Integer -> String -> IO Bool+integerFromDecimal ptr str =+   withCString str $ \strPtr ->+      (0==) <$> Fmpz.set_str ptr strPtr 10++integerFlintFromPrelude :: Ptr FLINT.Integer -> Integer -> IO ()+integerFlintFromPrelude ptr n =+   withCString (printf "%x" n) $ \strPtr ->+      void $ Fmpz.set_str ptr strPtr 16
+ src/LAPACK.hs view
@@ -0,0 +1,101 @@+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 Shell.Utility.Verbosity as Verbosity++import System.Environment (getArgs)++++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+   let verbosity = Verbosity.normal+   (reactantSet,reactantsWithoutDuplicates)+      <- preprocessArguments verbosity =<< getArgs++   {-+   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 reactantsWithoutDuplicates $+      map (cancelIntegerVector . fmap round . Array.toMap) $+      Matrix.toColumns $ nullspace tall
+ src/LAPACK/Common.hs view
@@ -0,0 +1,31 @@+module LAPACK.Common where++import Common (Reactant)++import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.Map as Map+import Data.Set (Set)++++toSquare :: (Shape.C sh) => Matrix.General sh sh a -> Matrix.Square sh a+toSquare a =+   ArrMatrix.liftOmni1+      (Array.reshape (MatrixShape.square (ArrMatrix.order a) (Matrix.height a)))+      a+++matrixFromReactantSet ::+   Set Reactant -> Matrix.General (Set String) (Set Reactant) Double+matrixFromReactantSet reactantSet =+   let elements = foldMap Map.keysSet reactantSet+   in Matrix.fromColumnArray elements $+      fmap (Array.fromAssociations 0 elements . Map.toList . fmap fromInteger) $+      BoxedArray.indices reactantSet
+ src/Main.hs view
@@ -0,0 +1,68 @@+module Main where++import Common++import qualified RowEchelon++import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import qualified Data.Set as Set+import qualified Data.List.HT as ListHT+import qualified Data.List as List+import qualified Data.Ratio as Ratio+import Data.Array.Comfort.Boxed (Array)+import Data.Set (Set)+import Data.Ratio ((%))++import qualified Shell.Utility.Verbosity as Verbosity++import System.Environment (getArgs)++++matrixFromReactantSet ::+   Set Reactant -> Array (Set String, Set Reactant) Rational+matrixFromReactantSet reactantSet =+   let elements = foldMap Map.keysSet reactantSet+   in Array.fromAssociations 0 (elements, reactantSet) $ do+         r <- Set.toList reactantSet+         (e,n) <- Map.toList r+         return ((e,r), fromInteger n)+++gcdRat :: Rational -> Rational -> Rational+gcdRat x y =+   gcd (Ratio.numerator x) (Ratio.numerator y)+   %+   lcm (Ratio.denominator x) (Ratio.denominator y)++integerFromRational :: Rational -> Integer+integerFromRational r =+   let n = round r in+   if r == n%1 then n else error "gcdRat failed"++normalizeRationalVector :: (Functor f, Foldable f) => f Rational -> f Integer+normalizeRationalVector v =+   let d = Fold.foldl gcdRat 0 v+   in fmap (integerFromRational . (/d)) v++++main :: IO ()+main = do+   let verbosity = Verbosity.normal+   (reactantSet,reactantsWithoutDuplicates)+      <- preprocessArguments verbosity =<< getArgs++   let mat = matrixFromReactantSet reactantSet+   let sols =+         map (Map.fromList . zip (Set.toAscList reactantSet)) $+         List.transpose $ RowEchelon.matrixRows $+         RowEchelon.nullspace $+         RowEchelon.matrixFromRows $+         ListHT.sliceVertical (Set.size reactantSet) $+         Array.toList mat++   displaySolutions verbosity reactantsWithoutDuplicates $+      map normalizeRationalVector sols