{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TupleSections #-}
module Language.Drasil.Chunk.DifferentialModel (
  -- * Export Data Type
  DifferentialModel(..), ODESolverFormat(..), InitialValueProblem(..),
  -- * Input Language
  ($^^), ($*), ($+),
  -- * Constructors
  makeAODESolverFormat, makeAIVP, makeASystemDE, makeASingleDE,
  formEquations
) where

import Control.Lens (makeLenses, (^.), view)
import Language.Drasil.Chunk.Concept (ConceptChunk, dccWDS)
import Language.Drasil.UID (HasUID(uid))
import Language.Drasil.Classes (Express(..),
  ConceptDomain(..), Definition(..), Idea(..), NamedIdea(..))
import Language.Drasil.ModelExpr.Lang (ModelExpr)
import Language.Drasil.NounPhrase.Core (NP)
import Language.Drasil.Sentence (Sentence)
import Language.Drasil.Expr.Lang (Expr(..))
import Language.Drasil.Chunk.Unital (UnitalChunk)
import Language.Drasil.ModelExpr.Class (ModelExprC(nthderiv, equiv))
import Language.Drasil.Expr.Class (ExprC(..), columnVec)
import Language.Drasil.Chunk.Constrained (ConstrConcept)
import Language.Drasil.Chunk.Quantity (qw)
import Language.Drasil.Literal.Class (LiteralC(exactDbl, int))
import Data.List (find)
import Language.Drasil.WellTyped (RequiresChecking (requiredChecks))
import Language.Drasil.Space (Space, HasSpace (..))

-- | Unknown is nth order of the dependent variable 
type Unknown = Integer

-- | Term consist of a coefficient and an unknown (order)
data Term = T{
  -- | the coefficient
  Term -> Expr
_coeff :: Expr,
  -- | the order
  Term -> Unknown
_unk :: Unknown
}
makeLenses ''Term

-- | LHS is a collection of Terms
type LHS = [Term]

-- | Operation connect the dependent variable and the order
{-
  e.g. depVar $^^ d
  note: depVar is a dummy variable. It keeps the shape of the syntax.
-}
($^^) :: ConstrConcept -> Integer -> Unknown
$^^ :: ConstrConcept -> Unknown -> Unknown
($^^) ConstrConcept
_ Unknown
unk' = Unknown
unk'

-- | Operation represent multiple
{-
  e.g. exactDbl 1 $* (opProcessVariable $^^ 2), 
  exactDbl 1 is the the coefficient, 
  (opProcessVariable $^^ 2) is the 2rd order of opProcessVariable
-}
($*) :: Expr -> Unknown -> Term
$* :: Expr -> Unknown -> Term
($*) = Expr -> Unknown -> Term
T

-- | Operation represent plus (collection Terms)
{-
  e.g. [exactDbl 1 $* (opProcessVariable $^^ 2)]
       $+ (exactDbl 1 `addRe` sy qdDerivGain $* (opProcessVariable $^^ 1))
  [exactDbl 1 $* (opProcessVariable $^^ 2)] is a collection with a single Term, 
  (exactDbl 1 `addRe` sy qdDerivGain $* (opProcessVariable $^^ 1)) is the appended element
-}
($+) :: [Term] -> Term -> LHS
$+ :: [Term] -> Term -> [Term]
($+) [Term]
xs Term
x  = [Term]
xs forall a. [a] -> [a] -> [a]
++ [Term
x]

-- | Describe the structural content of a system of linear ODEs with six necessary fields
data DifferentialModel = SystemOfLinearODEs {
  -- | independent variable, often time
  DifferentialModel -> UnitalChunk
_indepVar :: UnitalChunk,
  -- | dependent variable
  DifferentialModel -> ConstrConcept
_depVar :: ConstrConcept,
  -- | coefficients matrix
  DifferentialModel -> [[Expr]]
_coefficients :: [[Expr]],
  -- | unknowns column vector (orders)
  DifferentialModel -> [Unknown]
_unknowns :: [Unknown],
  -- | constant column vector 
  DifferentialModel -> [Expr]
_dmConstants :: [Expr],
  -- | meta data
  DifferentialModel -> ConceptChunk
_dmconc :: ConceptChunk
}
makeLenses ''DifferentialModel

-- | Information for solving an initial value problem
data InitialValueProblem = IVP{
  -- | initial time
  InitialValueProblem -> Expr
initTime :: Expr,
  -- | end time
  InitialValueProblem -> Expr
finalTime :: Expr,
  -- | initial values
  InitialValueProblem -> [Expr]
initValues :: [Expr]
}

-- | Acceptable format for ODE solvers, represent the structure of X' = AX + B
-- X' is a column vector of first-order unknowns
data ODESolverFormat = X'{
  -- | represent A, the coefficient matrix with identity matrix
  ODESolverFormat -> [[Expr]]
coeffVects :: [[Expr]],
  -- | combing with the dependent variable. it represents X, the unknown column vector after reduce the highest order.
  ODESolverFormat -> [Unknown]
unknownVect :: [Integer],
  -- | represent B, the constant column vector with identity matrix
  ODESolverFormat -> [Expr]
constantVect :: [Expr]
}

-- | Finds the 'UID' of the 'ConceptChunk' used to make the 'DifferentialModel'.
instance HasUID        DifferentialModel where uid :: Lens' DifferentialModel UID
uid = Lens' DifferentialModel ConceptChunk
dmconc forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall c. HasUID c => Lens' c UID
uid
-- | Equal if 'UID's are equal.
instance Eq            DifferentialModel where DifferentialModel
a == :: DifferentialModel -> DifferentialModel -> Bool
== DifferentialModel
b = (DifferentialModel
a forall s a. s -> Getting a s a -> a
^. forall c. HasUID c => Lens' c UID
uid) forall a. Eq a => a -> a -> Bool
== (DifferentialModel
b forall s a. s -> Getting a s a -> a
^. forall c. HasUID c => Lens' c UID
uid)
-- | Finds the term ('NP') of the 'ConceptChunk' used to make the 'DifferentialModel'.
instance NamedIdea     DifferentialModel where term :: Lens' DifferentialModel NP
term = Lens' DifferentialModel ConceptChunk
dmconc forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall c. NamedIdea c => Lens' c NP
term
-- | Finds the idea contained in the 'ConceptChunk' used to make the 'DifferentialModel'.
instance Idea          DifferentialModel where getA :: DifferentialModel -> Maybe String
getA = forall c. Idea c => c -> Maybe String
getA forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Lens' DifferentialModel ConceptChunk
dmconc
-- | Finds the definition contained in the 'ConceptChunk' used to make the 'DifferentialModel'.
instance Definition    DifferentialModel where defn :: Lens' DifferentialModel Sentence
defn = Lens' DifferentialModel ConceptChunk
dmconc forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall c. Definition c => Lens' c Sentence
defn
-- | Finds the domain of the 'ConceptChunk' used to make the 'DifferentialModel'.
instance ConceptDomain DifferentialModel where cdom :: DifferentialModel -> [UID]
cdom = forall c. ConceptDomain c => c -> [UID]
cdom forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Lens' DifferentialModel ConceptChunk
dmconc
-- | Convert the 'DifferentialModel' into the model expression language.
instance Express       DifferentialModel where express :: DifferentialModel -> ModelExpr
express = DifferentialModel -> ModelExpr
formStdODE

instance RequiresChecking DifferentialModel Expr Space where
  requiredChecks :: DifferentialModel -> [(Expr, Space)]
requiredChecks DifferentialModel
dmo = forall a b. (a -> b) -> [a] -> [b]
map (, DifferentialModel
dmo forall s a. s -> Getting a s a -> a
^. (Lens' DifferentialModel ConstrConcept
depVar forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall c. HasSpace c => Getter c Space
typ)) forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Unknown] -> [Expr] -> ConstrConcept -> [Expr]
formEquations (ODESolverFormat -> [[Expr]]
coeffVects ODESolverFormat
dm) (ODESolverFormat -> [Unknown]
unknownVect ODESolverFormat
dm) (ODESolverFormat -> [Expr]
constantVect ODESolverFormat
dm) (DifferentialModel -> ConstrConcept
_depVar DifferentialModel
dmo)
    where dm :: ODESolverFormat
dm = DifferentialModel -> ODESolverFormat
makeAODESolverFormat DifferentialModel
dmo

-- | Set the expression be a system of linear ODE to Ax = b
formStdODE :: DifferentialModel -> ModelExpr
formStdODE :: DifferentialModel -> ModelExpr
formStdODE DifferentialModel
d
  | Int
size forall a. Eq a => a -> a -> Bool
== Int
1 = [Expr] -> [ModelExpr] -> [Expr] -> ModelExpr
formASingleODE (forall a. [a] -> a
head (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [[Expr]]
coefficients)) [ModelExpr]
unknownVec (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [Expr]
dmConstants)
  | Bool
otherwise = forall r. ModelExprC r => [r] -> r
equiv (ModelExpr
coeffsMatix forall r. ExprC r => r -> r -> r
$. forall r. ExprC r => [r] -> r
columnVec [ModelExpr]
unknownVec forall a. a -> [a] -> [a]
: [ModelExpr]
constantVec)
  where size :: Int
size = forall (t :: * -> *) a. Foldable t => t a -> Int
length (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [[Expr]]
coefficients)
        coeffsMatix :: ModelExpr
coeffsMatix = forall c. Express c => c -> ModelExpr
express(forall r. ExprC r => [[r]] -> r
matrix (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [[Expr]]
coefficients))
        unknownVec :: [ModelExpr]
unknownVec = [Unknown] -> ConstrConcept -> UnitalChunk -> [ModelExpr]
formAllUnknown (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [Unknown]
unknowns) (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel ConstrConcept
depVar) (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel UnitalChunk
indepVar)
        constantVec :: [ModelExpr]
constantVec = [forall c. Express c => c -> ModelExpr
express (forall r. ExprC r => [r] -> r
columnVec (DifferentialModel
d forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [Expr]
dmConstants))]

-- | Set the single ODE to a flat equation form, "left hand side" = "right hand side"
formASingleODE :: [Expr] -> [ModelExpr] -> [Expr] -> ModelExpr
formASingleODE :: [Expr] -> [ModelExpr] -> [Expr] -> ModelExpr
formASingleODE [Expr]
coeffs [ModelExpr]
unks [Expr]
consts = forall r. ModelExprC r => [r] -> r
equiv (ModelExpr
lhs forall a. a -> [a] -> [a]
: [ModelExpr]
rhs)
  where lhs :: ModelExpr
lhs = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall r. ExprC r => r -> r -> r
addRe (forall a b. (a -> b) -> [a] -> [b]
map (\(Expr, ModelExpr)
x-> forall c. Express c => c -> ModelExpr
express (forall a b. (a, b) -> a
fst (Expr, ModelExpr)
x) forall r. ExprC r => r -> r -> r
`mulRe` forall a b. (a, b) -> b
snd (Expr, ModelExpr)
x) forall a b. (a -> b) -> a -> b
$ [Expr] -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff [Expr]
coeffs [ModelExpr]
unks)
        rhs :: [ModelExpr]
rhs = forall a b. (a -> b) -> [a] -> [b]
map forall c. Express c => c -> ModelExpr
express [Expr]
consts

-- | Remove zero coefficients for the displaying purpose
filterZeroCoeff :: [Expr] -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff :: [Expr] -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff [Expr]
es [ModelExpr]
mes = forall a. (a -> Bool) -> [a] -> [a]
filter (\(Expr, ModelExpr)
x -> forall a b. (a, b) -> a
fst (Expr, ModelExpr)
x forall a. Eq a => a -> a -> Bool
/= forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [Expr]
es [ModelExpr]
mes

-- | Form all derivatives for the displaying purpose
formAllUnknown :: [Unknown] -> ConstrConcept -> UnitalChunk -> [ModelExpr]
formAllUnknown :: [Unknown] -> ConstrConcept -> UnitalChunk -> [ModelExpr]
formAllUnknown [Unknown]
unks ConstrConcept
dep UnitalChunk
ind = forall a b. (a -> b) -> [a] -> [b]
map (\Unknown
x -> Unknown -> ConstrConcept -> UnitalChunk -> ModelExpr
formAUnknown Unknown
x ConstrConcept
dep UnitalChunk
ind) [Unknown]
unks

-- | Form a derivative for the displaying purpose
formAUnknown :: Unknown -> ConstrConcept-> UnitalChunk -> ModelExpr
formAUnknown :: Unknown -> ConstrConcept -> UnitalChunk -> ModelExpr
formAUnknown Unknown
unk'' ConstrConcept
dep = forall r c.
(ModelExprC r, HasUID c, HasSymbol c) =>
Unknown -> r -> c -> r
nthderiv (forall a. Integral a => a -> Unknown
toInteger Unknown
unk'') (forall r c. (ExprC r, HasUID c, HasSymbol c) => c -> r
sy (forall q. (Quantity q, MayHaveUnit q) => q -> QuantityDict
qw ConstrConcept
dep))

-- |   Create a 'DifferentialModel' by giving a independent variable, a dependent variable a canonical matrix form, and conceptChuck.
{-
  canonical matrix form: Ax = b
    A is a known m*n matrix that contains coefficients, 
    x is an n-vector that contain derivatives of dependent variables
    b is an m-vector that contain constants
  conceptChuck: 
    uid ('String'), term ('NP'), definition ('Sentence').
-}
makeASystemDE :: UnitalChunk -> ConstrConcept -> [[Expr]] -> [Unknown] -> [Expr]-> String -> NP -> Sentence -> DifferentialModel
makeASystemDE :: UnitalChunk
-> ConstrConcept
-> [[Expr]]
-> [Unknown]
-> [Expr]
-> String
-> NP
-> Sentence
-> DifferentialModel
makeASystemDE UnitalChunk
indepVar' ConstrConcept
depVar' [[Expr]]
coeffs [Unknown]
unks [Expr]
const' String
id' NP
term' Sentence
defn'
 | forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
coeffs forall a. Eq a => a -> a -> Bool
/= forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
const' =
  forall a. HasCallStack => String -> a
error String
"Length of coefficients matrix should equal to the length of the constant vector"
 | Bool -> Bool
not forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks =
  forall a. HasCallStack => String -> a
error String
"The length of each row vector in coefficients need to equal to the length of unknowns vector"
 | Bool -> Bool
not forall a b. (a -> b) -> a -> b
$ [Unknown] -> Bool
isUnknownDescending [Unknown]
unks =
  forall a. HasCallStack => String -> a
error String
"The order of giving unknowns need to be descending"
 | Bool
otherwise = UnitalChunk
-> ConstrConcept
-> [[Expr]]
-> [Unknown]
-> [Expr]
-> ConceptChunk
-> DifferentialModel
SystemOfLinearODEs UnitalChunk
indepVar' ConstrConcept
depVar' [[Expr]]
coeffs [Unknown]
unks [Expr]
const'(String -> NP -> Sentence -> ConceptChunk
dccWDS String
id' NP
term' Sentence
defn')

-- | Create a 'DifferentialModel' by the input language
makeASingleDE :: UnitalChunk -> ConstrConcept -> LHS -> Expr-> String -> NP -> Sentence -> DifferentialModel
makeASingleDE :: UnitalChunk
-> ConstrConcept
-> [Term]
-> Expr
-> String
-> NP
-> Sentence
-> DifferentialModel
makeASingleDE UnitalChunk
indepVar'' ConstrConcept
depVar'' [Term]
lhs Expr
const'' String
id'' NP
term'' Sentence
defn''
 | forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
coeffs forall a. Eq a => a -> a -> Bool
/= forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr
const''] =
  forall a. HasCallStack => String -> a
error String
"Length of coefficients matrix should equal to the length of the constant vector"
 | Bool -> Bool
not forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks =
  forall a. HasCallStack => String -> a
error String
"The length of each row vector in coefficients need to equal to the length of unknowns vector"
 | Bool
otherwise = UnitalChunk
-> ConstrConcept
-> [[Expr]]
-> [Unknown]
-> [Expr]
-> ConceptChunk
-> DifferentialModel
SystemOfLinearODEs UnitalChunk
indepVar'' ConstrConcept
depVar'' [[Expr]]
coeffs [Unknown]
unks [Expr
const''](String -> NP -> Sentence -> ConceptChunk
dccWDS String
id'' NP
term'' Sentence
defn'')
  where unks :: [Unknown]
unks = Unknown -> ConstrConcept -> [Unknown]
createAllUnknowns([Term] -> Term
findHighestOrder [Term]
lhs forall s a. s -> Getting a s a -> a
^. Lens' Term Unknown
unk) ConstrConcept
depVar''
        coeffs :: [[Expr]]
coeffs = [[Term] -> [Unknown] -> [Expr]
createCoefficients [Term]
lhs [Unknown]
unks]

-- | Function to check whether dimension of coefficient is match with the unknown vector
isCoeffsMatchUnknowns :: [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns :: [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns [] [Unknown]
_ = forall a. HasCallStack => String -> a
error String
"Coefficients matrix can not be empty"
isCoeffsMatchUnknowns [[Expr]]
_ [] = forall a. HasCallStack => String -> a
error String
"Unknowns column vector can not be empty"
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\ [Expr]
x -> Bool -> Bool -> Bool
(&&) (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
x forall a. Eq a => a -> a -> Bool
== forall (t :: * -> *) a. Foldable t => t a -> Int
length [Unknown]
unks)) Bool
True [[Expr]]
coeffs

-- | Function to check whether the of giving the unknown vector is descending
isUnknownDescending :: [Unknown] -> Bool
isUnknownDescending :: [Unknown] -> Bool
isUnknownDescending [] = Bool
True
isUnknownDescending [Unknown
_] = Bool
True
isUnknownDescending (Unknown
x:Unknown
y:[Unknown]
xs) = Unknown
x forall a. Ord a => a -> a -> Bool
> Unknown
y Bool -> Bool -> Bool
&& [Unknown] -> Bool
isUnknownDescending [Unknown]
xs

-- | Find the order in left hand side
findHighestOrder :: LHS -> Term
findHighestOrder :: [Term] -> Term
findHighestOrder = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldr1 (\Term
x Term
y -> if Term
x forall s a. s -> Getting a s a -> a
^. Lens' Term Unknown
unk forall a. Ord a => a -> a -> Bool
>= Term
y forall s a. s -> Getting a s a -> a
^. Lens' Term Unknown
unk then Term
x else Term
y)

-- | Create all possible unknowns based on the highest order.
-- The order of the result list is from the order to zero.
createAllUnknowns :: Unknown -> ConstrConcept -> [Unknown]
createAllUnknowns :: Unknown -> ConstrConcept -> [Unknown]
createAllUnknowns Unknown
highestUnk ConstrConcept
depv
  | Unknown
highestUnk  forall a. Eq a => a -> a -> Bool
== Unknown
0  = [Unknown
highestUnk]
  | Bool
otherwise = Unknown
highestUnk forall a. a -> [a] -> [a]
: Unknown -> ConstrConcept -> [Unknown]
createAllUnknowns (Unknown
highestUnk forall a. Num a => a -> a -> a
- Unknown
1) ConstrConcept
depv

-- | Create Coefficients base on all possible unknowns
createCoefficients :: LHS -> [Unknown] -> [Expr]
createCoefficients :: [Term] -> [Unknown] -> [Expr]
createCoefficients [] [Unknown]
_ = forall a. HasCallStack => String -> a
error String
"Left hand side is an empty list"
createCoefficients [Term]
_ [] = []
createCoefficients [Term]
lhs (Unknown
x:[Unknown]
xs) = Maybe Term -> Expr
genCoefficient (Unknown -> [Term] -> Maybe Term
findCoefficient Unknown
x [Term]
lhs) forall a. a -> [a] -> [a]
: [Term] -> [Unknown] -> [Expr]
createCoefficients [Term]
lhs [Unknown]
xs

-- | Get the coefficient, if it is Nothing, return zero
genCoefficient :: Maybe Term -> Expr
genCoefficient :: Maybe Term -> Expr
genCoefficient Maybe Term
Nothing = forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0
genCoefficient (Just Term
x) = Term
x forall s a. s -> Getting a s a -> a
^. Lens' Term Expr
coeff

-- | Find the term that match with the unknown (order)
findCoefficient :: Unknown -> LHS -> Maybe Term
findCoefficient :: Unknown -> [Term] -> Maybe Term
findCoefficient Unknown
u = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find(\Term
x -> Term
x forall s a. s -> Getting a s a -> a
^. Lens' Term Unknown
unk forall a. Eq a => a -> a -> Bool
== Unknown
u)

-- | Reduce the order
transUnknowns :: [Unknown] -> [Unknown]
transUnknowns :: [Unknown] -> [Unknown]
transUnknowns = forall a. [a] -> [a]
tail

-- | Reduce the target term, move the target to the left hand side and the rest of term to the right hand side. Then, reduce its coefficient-- 
transCoefficients :: [Expr] -> [Expr]
transCoefficients :: [Expr] -> [Expr]
transCoefficients [Expr]
es
  | forall a. [a] -> a
head [Expr]
es forall a. Eq a => a -> a -> Bool
== forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1 = [Expr] -> [Expr]
mapNeg forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
tail [Expr]
es
  | Bool
otherwise = [Expr] -> [Expr]
mapNeg forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall r. ExprC r => r -> r -> r
$/ forall a. [a] -> a
head [Expr]
es) [Expr]
es
    where mapNeg :: [Expr] -> [Expr]
mapNeg = forall a b. (a -> b) -> [a] -> [b]
map (\Expr
x -> if Expr
x forall a. Eq a => a -> a -> Bool
== forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 then forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 else forall r. ExprC r => r -> r
neg Expr
x)

-- | Add the "Identity Matrix" to Coefficients
-- len is the length of the identity row,
-- index is the location of identity value (start with 0)
addIdentityCoeffs :: [[Expr]] -> Int -> Int -> [[Expr]]
addIdentityCoeffs :: [[Expr]] -> Int -> Int -> [[Expr]]
addIdentityCoeffs [[Expr]]
es Int
len Int
index
  | Int
len forall a. Eq a => a -> a -> Bool
== Int
index forall a. Num a => a -> a -> a
+ Int
1 = [[Expr]]
es
  | Bool
otherwise = [[Expr]] -> Int -> Int -> [[Expr]]
addIdentityCoeffs (Int -> Int -> [Expr]
constIdentityRowVect Int
len Int
index forall a. a -> [a] -> [a]
: [[Expr]]
es) Int
len (Int
index forall a. Num a => a -> a -> a
+ Int
1)

-- | Construct an identity row vector.
constIdentityRowVect :: Int -> Int -> [Expr]
constIdentityRowVect :: Int -> Int -> [Expr]
constIdentityRowVect Int
len Int
index = Int -> [Expr] -> [Expr]
addIdentityValue Int
index forall a b. (a -> b) -> a -> b
$ forall a. Int -> a -> [a]
replicate Int
len forall a b. (a -> b) -> a -> b
$ forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0

-- | Recreate the identity row vector with identity value 
addIdentityValue :: Int -> [Expr] -> [Expr]
addIdentityValue :: Int -> [Expr] -> [Expr]
addIdentityValue Int
n [Expr]
es = forall a b. (a, b) -> a
fst ([Expr], [Expr])
splits forall a. [a] -> [a] -> [a]
++ [forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1] forall a. [a] -> [a] -> [a]
++ forall a. [a] -> [a]
tail (forall a b. (a, b) -> b
snd ([Expr], [Expr])
splits)
  where splits :: ([Expr], [Expr])
splits = forall a. Int -> [a] -> ([a], [a])
splitAt Int
n [Expr]
es

-- | Add zeroes to Constants
-- len is the size of new constant vector
addIdentityConsts :: [Expr] -> Int -> [Expr]
addIdentityConsts :: [Expr] -> Int -> [Expr]
addIdentityConsts [Expr]
expr Int
len = forall a. Int -> a -> [a]
replicate (Int
len forall a. Num a => a -> a -> a
- Int
1) (forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) forall a. [a] -> [a] -> [a]
++ [Expr]
expr

-- | divide the leading coefficient in the constant term
divideConstant :: Expr -> Expr -> Expr
divideConstant :: Expr -> Expr -> Expr
divideConstant Expr
a Expr
b
  | Expr
b forall a. Eq a => a -> a -> Bool
== forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 = forall a. HasCallStack => String -> a
error String
"Divisor can't be zero"
  | Expr
b forall a. Eq a => a -> a -> Bool
== forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1 = Expr
a
  | Bool
otherwise       = Expr
a forall r. ExprC r => r -> r -> r
$/ Expr
b

-- | Construct an ODESolverFormat for solving the ODE.
makeAODESolverFormat :: DifferentialModel -> ODESolverFormat
makeAODESolverFormat :: DifferentialModel -> ODESolverFormat
makeAODESolverFormat DifferentialModel
dm = [[Expr]] -> [Unknown] -> [Expr] -> ODESolverFormat
X' [[Expr]]
transEs [Unknown]
transUnks [Expr]
transConsts
  where transUnks :: [Unknown]
transUnks = [Unknown] -> [Unknown]
transUnknowns forall a b. (a -> b) -> a -> b
$ DifferentialModel
dm forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [Unknown]
unknowns
        transEs :: [[Expr]]
transEs = [[Expr]] -> Int -> Int -> [[Expr]]
addIdentityCoeffs [[Expr] -> [Expr]
transCoefficients forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head (DifferentialModel
dm forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [[Expr]]
coefficients)] (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Unknown]
transUnks) Int
0
        transConsts :: [Expr]
transConsts = [Expr] -> Int -> [Expr]
addIdentityConsts [forall a. [a] -> a
head (DifferentialModel
dm forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [Expr]
dmConstants) Expr -> Expr -> Expr
`divideConstant` forall a. [a] -> a
head (forall a. [a] -> a
head (DifferentialModel
dm forall s a. s -> Getting a s a -> a
^. Lens' DifferentialModel [[Expr]]
coefficients))] (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Unknown]
transUnks)

-- | Form well-formatted ODE equations which the ODE solvers can solve.
{-
  For example:
  the original fourth-order ODE: 
  y'''' + 3y′′ − sin(t)y′ + 8y = t2

  can be re-written to

  A                 *  X      + B     = X'

  0  0      1   0      x₄       0       equation 1
  0  1      0   0      x₃       0       equation 2
  1  0      0   0      x₂       0       equation 3 
  0 -3  sin(t) -8      x₁       t^2     equation 4 

  X = x₄,x₃,x₂,x₁ and x₁ = y, x₂ = y', x₃ = y'', x₄ = y'''

  A: [[Expr]], X: [Unknown], B: [Expr]
  return [equation 1, equation 2, equation 3, equation 4]
-}
formEquations :: [[Expr]] -> [Unknown] -> [Expr] -> ConstrConcept-> [Expr]
formEquations :: [[Expr]] -> [Unknown] -> [Expr] -> ConstrConcept -> [Expr]
formEquations [] [Unknown]
_ [Expr]
_ ConstrConcept
_ = []
formEquations [[Expr]]
_ [] [Expr]
_ ConstrConcept
_ = []
formEquations [[Expr]]
_ [Unknown]
_ [] ConstrConcept
_ = []
formEquations ([Expr]
ex:[[Expr]]
exs) [Unknown]
unks (Expr
y:[Expr]
ys) ConstrConcept
depVa =
  (if Expr
y forall a. Eq a => a -> a -> Bool
== forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 then Expr
finalExpr else Expr
finalExpr forall r. ExprC r => r -> r -> r
`addRe` Expr
y) forall a. a -> [a] -> [a]
: [[Expr]] -> [Unknown] -> [Expr] -> ConstrConcept -> [Expr]
formEquations [[Expr]]
exs [Unknown]
unks [Expr]
ys ConstrConcept
depVa
  where indexUnks :: [Expr]
indexUnks = forall a b. (a -> b) -> [a] -> [b]
map (forall r. ExprC r => r -> r -> r
idx (forall r c. (ExprC r, HasUID c, HasSymbol c) => c -> r
sy ConstrConcept
depVa) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall r. LiteralC r => Unknown -> r
int) [Unknown]
unks -- create X
        filteredExprs :: [(Expr, Expr)]
filteredExprs = forall a. (a -> Bool) -> [a] -> [a]
filter (\(Expr, Expr)
x -> forall a b. (a, b) -> a
fst (Expr, Expr)
x forall a. Eq a => a -> a -> Bool
/= forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) (forall a b. [a] -> [b] -> [(a, b)]
zip [Expr]
ex [Expr]
indexUnks) -- remove zero coefficients
        termExprs :: [Expr]
termExprs = forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall r. ExprC r => r -> r -> r
mulRe) [(Expr, Expr)]
filteredExprs -- multiple coefficient with depend variables
        finalExpr :: Expr
finalExpr = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall r. ExprC r => r -> r -> r
addRe [Expr]
termExprs -- add terms together

-- Construct an InitialValueProblem.
{-
  the first Expr: start time
  the second Expr: final time
  [Expr] : initial values
-}
makeAIVP :: Expr -> Expr -> [Expr] -> InitialValueProblem
makeAIVP :: Expr -> Expr -> [Expr] -> InitialValueProblem
makeAIVP = Expr -> Expr -> [Expr] -> InitialValueProblem
IVP