{-# 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 $+ sy qdDerivGain $* (opProcessVariable $^^ 1))
  [exactDbl 1 $* (opProcessVariable $^^ 2)] is a collection with a single Term, 
  (exactDbl 1 $+ sy qdDerivGain $* (opProcessVariable $^^ 1)) is the appended element
-}
($++) :: [Term] -> Term -> LHS
$++ :: [Term] -> Term -> [Term]
($++) [Term]
xs Term
x  = [Term]
xs [Term] -> [Term] -> [Term]
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 :: Getter DifferentialModel UID
uid = (ConceptChunk -> f ConceptChunk)
-> DifferentialModel -> f DifferentialModel
Lens' DifferentialModel ConceptChunk
dmconc ((ConceptChunk -> f ConceptChunk)
 -> DifferentialModel -> f DifferentialModel)
-> ((UID -> f UID) -> ConceptChunk -> f ConceptChunk)
-> (UID -> f UID)
-> DifferentialModel
-> f DifferentialModel
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (UID -> f UID) -> ConceptChunk -> f ConceptChunk
forall c. HasUID c => Getter c UID
Getter ConceptChunk UID
uid
-- | Equal if 'UID's are equal.
instance Eq            DifferentialModel where DifferentialModel
a == :: DifferentialModel -> DifferentialModel -> Bool
== DifferentialModel
b = (DifferentialModel
a DifferentialModel -> Getting UID DifferentialModel UID -> UID
forall s a. s -> Getting a s a -> a
^. Getting UID DifferentialModel UID
forall c. HasUID c => Getter c UID
Getter DifferentialModel UID
uid) UID -> UID -> Bool
forall a. Eq a => a -> a -> Bool
== (DifferentialModel
b DifferentialModel -> Getting UID DifferentialModel UID -> UID
forall s a. s -> Getting a s a -> a
^. Getting UID DifferentialModel UID
forall c. HasUID c => Getter c UID
Getter DifferentialModel UID
uid)
-- | Finds the term ('NP') of the 'ConceptChunk' used to make the 'DifferentialModel'.
instance NamedIdea     DifferentialModel where term :: Lens' DifferentialModel NP
term = (ConceptChunk -> f ConceptChunk)
-> DifferentialModel -> f DifferentialModel
Lens' DifferentialModel ConceptChunk
dmconc ((ConceptChunk -> f ConceptChunk)
 -> DifferentialModel -> f DifferentialModel)
-> ((NP -> f NP) -> ConceptChunk -> f ConceptChunk)
-> (NP -> f NP)
-> DifferentialModel
-> f DifferentialModel
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (NP -> f NP) -> ConceptChunk -> f ConceptChunk
forall c. NamedIdea c => Lens' c NP
Lens' ConceptChunk NP
term
-- | Finds the idea contained in the 'ConceptChunk' used to make the 'DifferentialModel'.
instance Idea          DifferentialModel where getA :: DifferentialModel -> Maybe String
getA = ConceptChunk -> Maybe String
forall c. Idea c => c -> Maybe String
getA (ConceptChunk -> Maybe String)
-> (DifferentialModel -> ConceptChunk)
-> DifferentialModel
-> Maybe String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting ConceptChunk DifferentialModel ConceptChunk
-> DifferentialModel -> ConceptChunk
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting ConceptChunk DifferentialModel ConceptChunk
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 = (ConceptChunk -> f ConceptChunk)
-> DifferentialModel -> f DifferentialModel
Lens' DifferentialModel ConceptChunk
dmconc ((ConceptChunk -> f ConceptChunk)
 -> DifferentialModel -> f DifferentialModel)
-> ((Sentence -> f Sentence) -> ConceptChunk -> f ConceptChunk)
-> (Sentence -> f Sentence)
-> DifferentialModel
-> f DifferentialModel
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Sentence -> f Sentence) -> ConceptChunk -> f ConceptChunk
forall c. Definition c => Lens' c Sentence
Lens' ConceptChunk Sentence
defn
-- | Finds the domain of the 'ConceptChunk' used to make the 'DifferentialModel'.
instance ConceptDomain DifferentialModel where cdom :: DifferentialModel -> [UID]
cdom = ConceptChunk -> [UID]
forall c. ConceptDomain c => c -> [UID]
cdom (ConceptChunk -> [UID])
-> (DifferentialModel -> ConceptChunk)
-> DifferentialModel
-> [UID]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting ConceptChunk DifferentialModel ConceptChunk
-> DifferentialModel -> ConceptChunk
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting ConceptChunk DifferentialModel ConceptChunk
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 = (Expr -> (Expr, Space)) -> [Expr] -> [(Expr, Space)]
forall a b. (a -> b) -> [a] -> [b]
map (, DifferentialModel
dmo DifferentialModel -> Getting Space DifferentialModel Space -> Space
forall s a. s -> Getting a s a -> a
^. ((ConstrConcept -> Const Space ConstrConcept)
-> DifferentialModel -> Const Space DifferentialModel
Lens' DifferentialModel ConstrConcept
depVar ((ConstrConcept -> Const Space ConstrConcept)
 -> DifferentialModel -> Const Space DifferentialModel)
-> ((Space -> Const Space Space)
    -> ConstrConcept -> Const Space ConstrConcept)
-> Getting Space DifferentialModel Space
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Space -> Const Space Space)
-> ConstrConcept -> Const Space ConstrConcept
forall c. HasSpace c => Getter c Space
Getter ConstrConcept Space
typ)) ([Expr] -> [(Expr, Space)]) -> [Expr] -> [(Expr, Space)]
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = [Expr] -> [ModelExpr] -> [Expr] -> ModelExpr
formASingleODE ([[Expr]] -> [Expr]
forall a. HasCallStack => [a] -> a
head (DifferentialModel
d DifferentialModel
-> Getting [[Expr]] DifferentialModel [[Expr]] -> [[Expr]]
forall s a. s -> Getting a s a -> a
^. Getting [[Expr]] DifferentialModel [[Expr]]
Lens' DifferentialModel [[Expr]]
coefficients)) [ModelExpr]
unknownVec (DifferentialModel
d DifferentialModel
-> Getting [Expr] DifferentialModel [Expr] -> [Expr]
forall s a. s -> Getting a s a -> a
^. Getting [Expr] DifferentialModel [Expr]
Lens' DifferentialModel [Expr]
dmConstants)
  | Bool
otherwise = [ModelExpr] -> ModelExpr
forall r. ModelExprC r => [r] -> r
equiv (ModelExpr
coeffsMatix ModelExpr -> ModelExpr -> ModelExpr
forall r. ExprC r => r -> r -> r
$. [ModelExpr] -> ModelExpr
forall r. ExprC r => [r] -> r
columnVec [ModelExpr]
unknownVec ModelExpr -> [ModelExpr] -> [ModelExpr]
forall a. a -> [a] -> [a]
: [ModelExpr]
constantVec)
  where size :: Int
size = [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (DifferentialModel
d DifferentialModel
-> Getting [[Expr]] DifferentialModel [[Expr]] -> [[Expr]]
forall s a. s -> Getting a s a -> a
^. Getting [[Expr]] DifferentialModel [[Expr]]
Lens' DifferentialModel [[Expr]]
coefficients)
        coeffsMatix :: ModelExpr
coeffsMatix = Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express([[Expr]] -> Expr
forall r. ExprC r => [[r]] -> r
matrix (DifferentialModel
d DifferentialModel
-> Getting [[Expr]] DifferentialModel [[Expr]] -> [[Expr]]
forall s a. s -> Getting a s a -> a
^. Getting [[Expr]] DifferentialModel [[Expr]]
Lens' DifferentialModel [[Expr]]
coefficients))
        unknownVec :: [ModelExpr]
unknownVec = [Unknown] -> ConstrConcept -> UnitalChunk -> [ModelExpr]
formAllUnknown (DifferentialModel
d DifferentialModel
-> Getting [Unknown] DifferentialModel [Unknown] -> [Unknown]
forall s a. s -> Getting a s a -> a
^. Getting [Unknown] DifferentialModel [Unknown]
Lens' DifferentialModel [Unknown]
unknowns) (DifferentialModel
d DifferentialModel
-> Getting ConstrConcept DifferentialModel ConstrConcept
-> ConstrConcept
forall s a. s -> Getting a s a -> a
^. Getting ConstrConcept DifferentialModel ConstrConcept
Lens' DifferentialModel ConstrConcept
depVar) (DifferentialModel
d DifferentialModel
-> Getting UnitalChunk DifferentialModel UnitalChunk -> UnitalChunk
forall s a. s -> Getting a s a -> a
^. Getting UnitalChunk DifferentialModel UnitalChunk
Lens' DifferentialModel UnitalChunk
indepVar)
        constantVec :: [ModelExpr]
constantVec = [Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express ([Expr] -> Expr
forall r. ExprC r => [r] -> r
columnVec (DifferentialModel
d DifferentialModel
-> Getting [Expr] DifferentialModel [Expr] -> [Expr]
forall s a. s -> Getting a s a -> a
^. Getting [Expr] DifferentialModel [Expr]
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 = [ModelExpr] -> ModelExpr
forall r. ModelExprC r => [r] -> r
equiv (ModelExpr
lhs ModelExpr -> [ModelExpr] -> [ModelExpr]
forall a. a -> [a] -> [a]
: [ModelExpr]
rhs)
  where lhs :: ModelExpr
lhs = (ModelExpr -> ModelExpr -> ModelExpr) -> [ModelExpr] -> ModelExpr
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 ModelExpr -> ModelExpr -> ModelExpr
forall r. ExprC r => r -> r -> r
($+) (((Expr, ModelExpr) -> ModelExpr)
-> [(Expr, ModelExpr)] -> [ModelExpr]
forall a b. (a -> b) -> [a] -> [b]
map (\(Expr, ModelExpr)
x-> Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express ((Expr, ModelExpr) -> Expr
forall a b. (a, b) -> a
fst (Expr, ModelExpr)
x) ModelExpr -> ModelExpr -> ModelExpr
forall r. ExprC r => r -> r -> r
$* (Expr, ModelExpr) -> ModelExpr
forall a b. (a, b) -> b
snd (Expr, ModelExpr)
x) ([(Expr, ModelExpr)] -> [ModelExpr])
-> [(Expr, ModelExpr)] -> [ModelExpr]
forall a b. (a -> b) -> a -> b
$ [Expr] -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff [Expr]
coeffs [ModelExpr]
unks)
        rhs :: [ModelExpr]
rhs = (Expr -> ModelExpr) -> [Expr] -> [ModelExpr]
forall a b. (a -> b) -> [a] -> [b]
map Expr -> ModelExpr
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 = ((Expr, ModelExpr) -> Bool)
-> [(Expr, ModelExpr)] -> [(Expr, ModelExpr)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Expr, ModelExpr)
x -> (Expr, ModelExpr) -> Expr
forall a b. (a, b) -> a
fst (Expr, ModelExpr)
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
/= Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) ([(Expr, ModelExpr)] -> [(Expr, ModelExpr)])
-> [(Expr, ModelExpr)] -> [(Expr, ModelExpr)]
forall a b. (a -> b) -> a -> b
$ [Expr] -> [ModelExpr] -> [(Expr, ModelExpr)]
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 = (Unknown -> ModelExpr) -> [Unknown] -> [ModelExpr]
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 = Unknown -> ModelExpr -> UnitalChunk -> ModelExpr
forall c.
(HasUID c, HasSymbol c) =>
Unknown -> ModelExpr -> c -> ModelExpr
forall r c.
(ModelExprC r, HasUID c, HasSymbol c) =>
Unknown -> r -> c -> r
nthderiv (Unknown -> Unknown
forall a. Integral a => a -> Unknown
toInteger Unknown
unk'') (QuantityDict -> ModelExpr
forall c. (HasUID c, HasSymbol c) => c -> ModelExpr
forall r c. (ExprC r, HasUID c, HasSymbol c) => c -> r
sy (ConstrConcept -> QuantityDict
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'
 | [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
coeffs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
const' =
  String -> DifferentialModel
forall a. HasCallStack => String -> a
error String
"Length of coefficients matrix should equal to the length of the constant vector"
 | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks =
  String -> DifferentialModel
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 (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ [Unknown] -> Bool
isUnknownDescending [Unknown]
unks =
  String -> DifferentialModel
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''
 | [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
coeffs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr
const''] =
  String -> DifferentialModel
forall a. HasCallStack => String -> a
error String
"Length of coefficients matrix should equal to the length of the constant vector"
 | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Unknown] -> Bool
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks =
  String -> DifferentialModel
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 Term -> Getting Unknown Term Unknown -> Unknown
forall s a. s -> Getting a s a -> a
^. Getting Unknown Term Unknown
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]
_ = String -> Bool
forall a. HasCallStack => String -> a
error String
"Coefficients matrix can not be empty"
isCoeffsMatchUnknowns [[Expr]]
_ [] = String -> Bool
forall a. HasCallStack => String -> a
error String
"Unknowns column vector can not be empty"
isCoeffsMatchUnknowns [[Expr]]
coeffs [Unknown]
unks = ([Expr] -> Bool -> Bool) -> Bool -> [[Expr]] -> Bool
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\ [Expr]
x -> Bool -> Bool -> Bool
(&&) ([Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Unknown] -> Int
forall a. [a] -> Int
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 Unknown -> Unknown -> Bool
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 = (Term -> Term -> Term) -> [Term] -> Term
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldr1 (\Term
x Term
y -> if Term
x Term -> Getting Unknown Term Unknown -> Unknown
forall s a. s -> Getting a s a -> a
^. Getting Unknown Term Unknown
Lens' Term Unknown
unk Unknown -> Unknown -> Bool
forall a. Ord a => a -> a -> Bool
>= Term
y Term -> Getting Unknown Term Unknown -> Unknown
forall s a. s -> Getting a s a -> a
^. Getting Unknown Term Unknown
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  Unknown -> Unknown -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown
0  = [Unknown
highestUnk]
  | Bool
otherwise = Unknown
highestUnk Unknown -> [Unknown] -> [Unknown]
forall a. a -> [a] -> [a]
: Unknown -> ConstrConcept -> [Unknown]
createAllUnknowns (Unknown
highestUnk Unknown -> Unknown -> Unknown
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]
_ = String -> [Expr]
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) Expr -> [Expr] -> [Expr]
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 = Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0
genCoefficient (Just Term
x) = Term
x Term -> Getting Expr Term Expr -> Expr
forall s a. s -> Getting a s a -> a
^. Getting Expr Term Expr
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 = (Term -> Bool) -> [Term] -> Maybe Term
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find(\Term
x -> Term
x Term -> Getting Unknown Term Unknown -> Unknown
forall s a. s -> Getting a s a -> a
^. Getting Unknown Term Unknown
Lens' Term Unknown
unk Unknown -> Unknown -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown
u)

-- | Reduce the order
transUnknowns :: [Unknown] -> [Unknown]
transUnknowns :: [Unknown] -> [Unknown]
transUnknowns = [Unknown] -> [Unknown]
forall a. HasCallStack => [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
  | [Expr] -> Expr
forall a. HasCallStack => [a] -> a
head [Expr]
es Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1 = [Expr] -> [Expr]
mapNeg ([Expr] -> [Expr]) -> [Expr] -> [Expr]
forall a b. (a -> b) -> a -> b
$ [Expr] -> [Expr]
forall a. HasCallStack => [a] -> [a]
tail [Expr]
es
  | Bool
otherwise = [Expr] -> [Expr]
mapNeg ([Expr] -> [Expr]) -> [Expr] -> [Expr]
forall a b. (a -> b) -> a -> b
$ [Expr] -> [Expr]
forall a. HasCallStack => [a] -> [a]
tail ([Expr] -> [Expr]) -> [Expr] -> [Expr]
forall a b. (a -> b) -> a -> b
$ (Expr -> Expr) -> [Expr] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map (Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
$/ [Expr] -> Expr
forall a. HasCallStack => [a] -> a
head [Expr]
es) [Expr]
es
    where mapNeg :: [Expr] -> [Expr]
mapNeg = (Expr -> Expr) -> [Expr] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map (\Expr
x -> if Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 then Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 else Expr -> Expr
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 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
index Int -> Int -> Int
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 [Expr] -> [[Expr]] -> [[Expr]]
forall a. a -> [a] -> [a]
: [[Expr]]
es) Int
len (Int
index Int -> Int -> Int
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 ([Expr] -> [Expr]) -> [Expr] -> [Expr]
forall a b. (a -> b) -> a -> b
$ Int -> Expr -> [Expr]
forall a. Int -> a -> [a]
replicate Int
len (Expr -> [Expr]) -> Expr -> [Expr]
forall a b. (a -> b) -> a -> b
$ Unknown -> Expr
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 = ([Expr], [Expr]) -> [Expr]
forall a b. (a, b) -> a
fst ([Expr], [Expr])
splits [Expr] -> [Expr] -> [Expr]
forall a. [a] -> [a] -> [a]
++ [Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1] [Expr] -> [Expr] -> [Expr]
forall a. [a] -> [a] -> [a]
++ [Expr] -> [Expr]
forall a. HasCallStack => [a] -> [a]
tail (([Expr], [Expr]) -> [Expr]
forall a b. (a, b) -> b
snd ([Expr], [Expr])
splits)
  where splits :: ([Expr], [Expr])
splits = Int -> [Expr] -> ([Expr], [Expr])
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 = Int -> Expr -> [Expr]
forall a. Int -> a -> [a]
replicate (Int
len Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) [Expr] -> [Expr] -> [Expr]
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 Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 = String -> Expr
forall a. HasCallStack => String -> a
error String
"Divisor can't be zero"
  | Expr
b Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1 = Expr
a
  | Bool
otherwise       = Expr
a Expr -> Expr -> Expr
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 ([Unknown] -> [Unknown]) -> [Unknown] -> [Unknown]
forall a b. (a -> b) -> a -> b
$ DifferentialModel
dm DifferentialModel
-> Getting [Unknown] DifferentialModel [Unknown] -> [Unknown]
forall s a. s -> Getting a s a -> a
^. Getting [Unknown] DifferentialModel [Unknown]
Lens' DifferentialModel [Unknown]
unknowns
        transEs :: [[Expr]]
transEs = [[Expr]] -> Int -> Int -> [[Expr]]
addIdentityCoeffs [[Expr] -> [Expr]
transCoefficients ([Expr] -> [Expr]) -> [Expr] -> [Expr]
forall a b. (a -> b) -> a -> b
$ [[Expr]] -> [Expr]
forall a. HasCallStack => [a] -> a
head (DifferentialModel
dm DifferentialModel
-> Getting [[Expr]] DifferentialModel [[Expr]] -> [[Expr]]
forall s a. s -> Getting a s a -> a
^. Getting [[Expr]] DifferentialModel [[Expr]]
Lens' DifferentialModel [[Expr]]
coefficients)] ([Unknown] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Unknown]
transUnks) Int
0
        transConsts :: [Expr]
transConsts = [Expr] -> Int -> [Expr]
addIdentityConsts [[Expr] -> Expr
forall a. HasCallStack => [a] -> a
head (DifferentialModel
dm DifferentialModel
-> Getting [Expr] DifferentialModel [Expr] -> [Expr]
forall s a. s -> Getting a s a -> a
^. Getting [Expr] DifferentialModel [Expr]
Lens' DifferentialModel [Expr]
dmConstants) Expr -> Expr -> Expr
`divideConstant` [Expr] -> Expr
forall a. HasCallStack => [a] -> a
head ([[Expr]] -> [Expr]
forall a. HasCallStack => [a] -> a
head (DifferentialModel
dm DifferentialModel
-> Getting [[Expr]] DifferentialModel [[Expr]] -> [[Expr]]
forall s a. s -> Getting a s a -> a
^. Getting [[Expr]] DifferentialModel [[Expr]]
Lens' DifferentialModel [[Expr]]
coefficients))] ([Unknown] -> Int
forall a. [a] -> Int
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 Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0 then Expr
finalExpr else Expr
finalExpr Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
$+ Expr
y) Expr -> [Expr] -> [Expr]
forall a. a -> [a] -> [a]
: [[Expr]] -> [Unknown] -> [Expr] -> ConstrConcept -> [Expr]
formEquations [[Expr]]
exs [Unknown]
unks [Expr]
ys ConstrConcept
depVa
  where indexUnks :: [Expr]
indexUnks = (Unknown -> Expr) -> [Unknown] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map (Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
idx (ConstrConcept -> Expr
forall c. (HasUID c, HasSymbol c) => c -> Expr
forall r c. (ExprC r, HasUID c, HasSymbol c) => c -> r
sy ConstrConcept
depVa) (Expr -> Expr) -> (Unknown -> Expr) -> Unknown -> Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unknown -> Expr
forall r. LiteralC r => Unknown -> r
int) [Unknown]
unks -- create X
        filteredExprs :: [(Expr, Expr)]
filteredExprs = ((Expr, Expr) -> Bool) -> [(Expr, Expr)] -> [(Expr, Expr)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Expr, Expr)
x -> (Expr, Expr) -> Expr
forall a b. (a, b) -> a
fst (Expr, Expr)
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
/= Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0) ([Expr] -> [Expr] -> [(Expr, Expr)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Expr]
ex [Expr]
indexUnks) -- remove zero coefficients
        termExprs :: [Expr]
termExprs = ((Expr, Expr) -> Expr) -> [(Expr, Expr)] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map ((Expr -> Expr -> Expr) -> (Expr, Expr) -> Expr
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
($*)) [(Expr, Expr)]
filteredExprs -- multiple coefficient with depend variables
        finalExpr :: Expr
finalExpr = (Expr -> Expr -> Expr) -> [Expr] -> Expr
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
($+) [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