simgi-0.3: src/ExtraFunctions.hs
{-# LANGUAGE ForeignFunctionInterface #-}
{-----------------------------------------------------------------
(c) 2008-2009 Markus Dittrich,
National Resource for Biomedical Supercomputing &
Carnegie Mellon University
This program 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.
This program 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 Version 3 for more details.
You should have received a copy of the GNU General Public
License along with this program; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
--------------------------------------------------------------------}
-- | definition of additional math and helper functions
module ExtraFunctions ( convert_rate
, erf
, erfc
, fact
, is_equal
, is_equal_with
, maybe_to_int
, maybe_to_positive_int
, real_exp
, to_int
) where
-- imports
import Foreign()
import Foreign.C.Types
import Prelude
-- local imports
import GenericModel (MathExpr(..))
import RpnData (RpnStack(..), RpnItem(..))
--import Debug.Trace
-- | a few constants
avogadroNum :: Double
avogadroNum = 6.0221415e23
-- | use glibc DBL_EPSILON
dbl_epsilon :: Double
dbl_epsilon = 2.2204460492503131e-16
-- | comparison function for doubles via dbl_epsion
is_equal :: Double -> Double -> Bool
is_equal x y = abs(x-y) <= abs(x) * dbl_epsilon
-- | comparison function for doubles via threshold
is_equal_with :: Double -> Double -> Double -> Bool
is_equal_with x y th = abs(x-y) <= abs(x) * th
-- | function checking if a Double can be interpreted as a non
-- negative Integer. We need this since all parsing of numbers
-- is done with Doubles but some functions only work for
-- non-negative integers such as factorial.
-- To check if we are dealing with Double, we convert to an
-- Integer via floor and the compare if the numbers are identical.
-- If yes, the number seems to be an Integer and we return it,
-- otherwise Nothing
maybe_to_positive_int :: Double -> Maybe Integer
maybe_to_positive_int x =
case (is_equal (fromInteger . floor $ x) x) && (x > 0.0) of
True -> Just $ floor x
False -> Nothing
-- | function checking if a Double can be interpreted as an
-- Integer. See is_positive_int for more detail
maybe_to_int :: Double -> Maybe Integer
maybe_to_int x =
case is_equal (fromInteger . floor $ x) x of
True -> Just $ floor x
False -> Nothing
-- | helper function for defining real powers
-- NOTE: We use glibc's pow function since it is more
-- precise than implementing it ourselves via, e.g.,
-- pow a x = exp $ x * log a
foreign import ccall "math.h pow"
c_pow :: CDouble -> CDouble -> CDouble
real_exp :: Double -> Double -> Double
real_exp a x = realToFrac $ c_pow (realToFrac a) (realToFrac x)
-- | factorial function
fact :: Integer -> Integer
fact 0 = 1
fact n = n * fact (n-1)
-- | error function erf(x)
-- we use a recursive solution of the Taylor series expansion
erf :: Double -> Double
erf x
| x == 0.0 = 0.0 -- our recursive alg. loops forever
-- in this case
| abs(x) > 2.2 = 1.0 - erfc x -- use erfc for numerical accuracy
| otherwise = 2.0 / sqrt(pi) * (erf_h x x 1.0)
where
erf_h old x_old n = let x_new = x_old * (x_next n)
tot = old + x_new
in
if abs(x_new/tot) < dbl_epsilon
then tot
else erf_h tot x_new (n+1.0)
-- Note: We need (2::Int) here to silence ghc
x_next n = -(x^(2::Int)) * (2.0*n-1.0)/(n * (2.0*n+1.0))
-- | complementary error function erfc(x) = 1 - erf(x)
-- we use a recursive solution of the continued fraction
-- expression of erfc(x) for it superior convergence
-- property. Here, we calculate the ith and (i+1)th convergent, (see
-- http://archives.math.utk.edu/articles/atuyl/confrac/intro.html)
-- and terminate when the relative difference is smaller than a
-- certain threshold.
erfc :: Double -> Double
erfc x
| abs(x) < 2.2 = 1.0 - erf(x) -- use erf(x) in [-2.2,2.2]
| signum(x) < 0 = 2.0 - erfc(-x) -- continued fraction expansion
-- only valid for x > 0
| otherwise = 1/sqrt(pi) * exp(-x^(2::Int))
* (erfc_h nc1 nc2 dc1 dc2 1.0)
where
nc1 = 1.0 :: Double -- numerator of 1st convergent
nc2 = x :: Double -- numerator of 2nd convergent
dc1 = x :: Double -- denominator of 1st convergent
dc2 = x^(2::Int)+0.5 :: Double -- denominator of 2nd convergent
erfc_h n1 n2 d1 d2 i =
let num_new = n1*i + n2*x
denom_new = d1*i + d2*x
d_old = n2/d2
d_new = num_new/denom_new
in
if abs((d_old - d_new)/d_new) < dbl_epsilon
then d_new
else erfc_h n2 num_new d2 denom_new (i+0.5)
-- | convert reaction propensities into rates if requested
-- by the user. For constants we simply multiply, for
-- rate functions we push the neccessary conversion onto
-- the stack
convert_rate :: MathExpr -> Int -> Double -> MathExpr
convert_rate theConst@(Constant c) order volume =
case order of
1 -> theConst
_ -> Constant $ c/(avogadroNum * volume)^(order-1)
convert_rate theFunc@(Function stack) order volume =
case order of
1 -> theFunc
_ -> let mult = 1.0/(avogadroNum * volume)^(order-1) in
Function . RpnStack $ (toList stack)
++ [Number mult,BinFunc (*)]
-- | convert a double to int
-- NOTE: presently, converting double -> int is done
-- via floor. Is this a good policy (once documented
-- properly)?
to_int :: (RealFrac a, Integral b) => a -> b
to_int = floor