diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/chemical-equation.cabal b/chemical-equation.cabal new file mode 100644 --- /dev/null +++ b/chemical-equation.cabal @@ -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 diff --git a/src/Common.hs b/src/Common.hs new file mode 100644 --- /dev/null +++ b/src/Common.hs @@ -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 diff --git a/src/FLINT.hs b/src/FLINT.hs new file mode 100644 --- /dev/null +++ b/src/FLINT.hs @@ -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 diff --git a/src/FLINT/Example.hs b/src/FLINT/Example.hs new file mode 100644 --- /dev/null +++ b/src/FLINT/Example.hs @@ -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 "" diff --git a/src/FLINT/Integer.hs b/src/FLINT/Integer.hs new file mode 100644 --- /dev/null +++ b/src/FLINT/Integer.hs @@ -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 () diff --git a/src/FLINT/MatrixInteger.hs b/src/FLINT/MatrixInteger.hs new file mode 100644 --- /dev/null +++ b/src/FLINT/MatrixInteger.hs @@ -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 diff --git a/src/FLINT/Type.hsc b/src/FLINT/Type.hsc new file mode 100644 --- /dev/null +++ b/src/FLINT/Type.hsc @@ -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} diff --git a/src/FLINT/Utility.hs b/src/FLINT/Utility.hs new file mode 100644 --- /dev/null +++ b/src/FLINT/Utility.hs @@ -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 diff --git a/src/LAPACK.hs b/src/LAPACK.hs new file mode 100644 --- /dev/null +++ b/src/LAPACK.hs @@ -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 diff --git a/src/LAPACK/Common.hs b/src/LAPACK/Common.hs new file mode 100644 --- /dev/null +++ b/src/LAPACK/Common.hs @@ -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 diff --git a/src/Main.hs b/src/Main.hs new file mode 100644 --- /dev/null +++ b/src/Main.hs @@ -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