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

import Control.Lens (makeLenses, (^.), view)
import Data.List (find)
import Data.List.NonEmpty (NonEmpty((:|)))
import qualified Data.List.NonEmpty as NE

import Drasil.Database (HasUID(uid), HasChunkRefs(..), mkUid)

import Language.Drasil
  (ConceptChunk, cncpt''', Express(..), Definition(..), Idea(..), NamedIdea(..)
  , ModelExpr, NP, Sentence, Expr, ModelExprC(nthderiv, equiv), DefinedQuantityDict
  , ExprC(..), columnVec, ConstrConcept, LiteralC(exactDbl, int), RequiresChecking (requiredChecks)
  , 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 -> DefinedQuantityDict
_indepVar :: DefinedQuantityDict,
  -- | dependent variable
  DifferentialModel -> ConstrConcept
_depVar :: ConstrConcept,
  -- | coefficients matrix
  DifferentialModel -> NonEmpty (NonEmpty Expr)
_coefficients :: NonEmpty (NonEmpty Expr),
  -- | unknowns column vector (orders)
  DifferentialModel -> [Unknown]
_unknowns :: [Unknown],
  -- | constant column vector
  DifferentialModel -> NonEmpty Expr
_dmConstants :: NonEmpty 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]
}

instance HasChunkRefs DifferentialModel where
  chunkRefs :: DifferentialModel -> Set UID
chunkRefs DifferentialModel
dm = [Set UID] -> Set UID
forall a. Monoid a => [a] -> a
mconcat
    [ DefinedQuantityDict -> Set UID
forall a. HasChunkRefs a => a -> Set UID
chunkRefs (DifferentialModel
dm DifferentialModel
-> Getting
     DefinedQuantityDict DifferentialModel DefinedQuantityDict
-> DefinedQuantityDict
forall s a. s -> Getting a s a -> a
^. Getting DefinedQuantityDict DifferentialModel DefinedQuantityDict
Lens' DifferentialModel DefinedQuantityDict
indepVar)
    , ConstrConcept -> Set UID
forall a. HasChunkRefs a => a -> Set UID
chunkRefs (DifferentialModel
dm DifferentialModel
-> Getting ConstrConcept DifferentialModel ConstrConcept
-> ConstrConcept
forall s a. s -> Getting a s a -> a
^. Getting ConstrConcept DifferentialModel ConstrConcept
Lens' DifferentialModel ConstrConcept
depVar)
    , ConceptChunk -> Set UID
forall a. HasChunkRefs a => a -> Set UID
chunkRefs (DifferentialModel
dm DifferentialModel
-> Getting ConceptChunk DifferentialModel ConceptChunk
-> ConceptChunk
forall s a. s -> Getting a s a -> a
^. Getting ConceptChunk DifferentialModel ConceptChunk
Lens' DifferentialModel ConceptChunk
dmconc)
    ]
  {-# INLINABLE chunkRefs #-}

-- | 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
-- | 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 = NonEmpty Expr -> [ModelExpr] -> NonEmpty Expr -> ModelExpr
formASingleODE (NonEmpty (NonEmpty Expr) -> NonEmpty Expr
forall a. NonEmpty a -> a
NE.head (DifferentialModel
d DifferentialModel
-> Getting
     (NonEmpty (NonEmpty Expr))
     DifferentialModel
     (NonEmpty (NonEmpty Expr))
-> NonEmpty (NonEmpty Expr)
forall s a. s -> Getting a s a -> a
^. Getting
  (NonEmpty (NonEmpty Expr))
  DifferentialModel
  (NonEmpty (NonEmpty Expr))
Lens' DifferentialModel (NonEmpty (NonEmpty Expr))
coefficients)) [ModelExpr]
unknownVec (DifferentialModel
d DifferentialModel
-> Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
-> NonEmpty Expr
forall s a. s -> Getting a s a -> a
^. Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
Lens' DifferentialModel (NonEmpty 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 = NonEmpty (NonEmpty Expr) -> Int
forall a. NonEmpty a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (DifferentialModel
d DifferentialModel
-> Getting
     (NonEmpty (NonEmpty Expr))
     DifferentialModel
     (NonEmpty (NonEmpty Expr))
-> NonEmpty (NonEmpty Expr)
forall s a. s -> Getting a s a -> a
^. Getting
  (NonEmpty (NonEmpty Expr))
  DifferentialModel
  (NonEmpty (NonEmpty Expr))
Lens' DifferentialModel (NonEmpty (NonEmpty Expr))
coefficients)
    coeffsMatix :: ModelExpr
coeffsMatix = Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express([[Expr]] -> Expr
forall r. ExprC r => [[r]] -> r
matrix ((NonEmpty Expr -> [Expr]) -> [NonEmpty Expr] -> [[Expr]]
forall a b. (a -> b) -> [a] -> [b]
map NonEmpty Expr -> [Expr]
forall a. NonEmpty a -> [a]
NE.toList (NonEmpty (NonEmpty Expr) -> [NonEmpty Expr]
forall a. NonEmpty a -> [a]
NE.toList (NonEmpty (NonEmpty Expr) -> [NonEmpty Expr])
-> NonEmpty (NonEmpty Expr) -> [NonEmpty Expr]
forall a b. (a -> b) -> a -> b
$ DifferentialModel
d DifferentialModel
-> Getting
     (NonEmpty (NonEmpty Expr))
     DifferentialModel
     (NonEmpty (NonEmpty Expr))
-> NonEmpty (NonEmpty Expr)
forall s a. s -> Getting a s a -> a
^. Getting
  (NonEmpty (NonEmpty Expr))
  DifferentialModel
  (NonEmpty (NonEmpty Expr))
Lens' DifferentialModel (NonEmpty (NonEmpty Expr))
coefficients)))
    unknownVec :: [ModelExpr]
unknownVec = [Unknown] -> ConstrConcept -> DefinedQuantityDict -> [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
     DefinedQuantityDict DifferentialModel DefinedQuantityDict
-> DefinedQuantityDict
forall s a. s -> Getting a s a -> a
^. Getting DefinedQuantityDict DifferentialModel DefinedQuantityDict
Lens' DifferentialModel DefinedQuantityDict
indepVar)
    constantVec :: [ModelExpr]
constantVec = [Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express ([Expr] -> Expr
forall r. ExprC r => [r] -> r
columnVec (NonEmpty Expr -> [Expr]
forall a. NonEmpty a -> [a]
NE.toList (NonEmpty Expr -> [Expr]) -> NonEmpty Expr -> [Expr]
forall a b. (a -> b) -> a -> b
$ DifferentialModel
d DifferentialModel
-> Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
-> NonEmpty Expr
forall s a. s -> Getting a s a -> a
^. Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
Lens' DifferentialModel (NonEmpty Expr)
dmConstants))]

-- | Set the single ODE to a flat equation form, "left hand side" = "right hand side"
formASingleODE :: NonEmpty Expr -> [ModelExpr] -> NonEmpty Expr -> ModelExpr
formASingleODE :: NonEmpty Expr -> [ModelExpr] -> NonEmpty Expr -> ModelExpr
formASingleODE NonEmpty Expr
coeffs [ModelExpr]
unks NonEmpty 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
x,ModelExpr
y) -> Expr -> ModelExpr
forall c. Express c => c -> ModelExpr
express Expr
x ModelExpr -> ModelExpr -> ModelExpr
forall r. ExprC r => r -> r -> r
$* ModelExpr
y) ([(Expr, ModelExpr)] -> [ModelExpr])
-> [(Expr, ModelExpr)] -> [ModelExpr]
forall a b. (a -> b) -> a -> b
$ NonEmpty Expr -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff NonEmpty 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] -> [ModelExpr]) -> [Expr] -> [ModelExpr]
forall a b. (a -> b) -> a -> b
$ NonEmpty Expr -> [Expr]
forall a. NonEmpty a -> [a]
NE.toList NonEmpty Expr
consts

-- | Remove zero coefficients for the displaying purpose
filterZeroCoeff :: NonEmpty Expr -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff :: NonEmpty Expr -> [ModelExpr] -> [(Expr, ModelExpr)]
filterZeroCoeff NonEmpty 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 (NonEmpty Expr -> [Expr]
forall a. NonEmpty a -> [a]
NE.toList NonEmpty Expr
es) [ModelExpr]
mes

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

-- | Form a derivative for the displaying purpose
formAUnknown :: Unknown -> ConstrConcept -> DefinedQuantityDict -> ModelExpr
formAUnknown :: Unknown -> ConstrConcept -> DefinedQuantityDict -> ModelExpr
formAUnknown Unknown
unk'' ConstrConcept
dep = Unknown -> ModelExpr -> DefinedQuantityDict -> ModelExpr
forall c.
(IsChunk c, HasSymbol c) =>
Unknown -> ModelExpr -> c -> ModelExpr
forall r c.
(ModelExprC r, IsChunk c, HasSymbol c) =>
Unknown -> r -> c -> r
nthderiv (Unknown -> Unknown
forall a. Integral a => a -> Unknown
toInteger Unknown
unk'') (ConstrConcept -> ModelExpr
forall c. (IsChunk c, HasSymbol c) => c -> ModelExpr
forall r c. (ExprC r, IsChunk c, HasSymbol c) => c -> r
sy 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 :: DefinedQuantityDict -> ConstrConcept -> NonEmpty (NonEmpty Expr) -> [Unknown] -> NonEmpty Expr -> String -> NP -> Sentence -> DifferentialModel
makeASystemDE :: DefinedQuantityDict
-> ConstrConcept
-> NonEmpty (NonEmpty Expr)
-> [Unknown]
-> NonEmpty Expr
-> String
-> NP
-> Sentence
-> DifferentialModel
makeASystemDE DefinedQuantityDict
indepVar' ConstrConcept
depVar' NonEmpty (NonEmpty Expr)
coeffs [Unknown]
unks NonEmpty Expr
const' String
id' NP
term' Sentence
defn'
 | NonEmpty (NonEmpty Expr) -> Int
forall a. NonEmpty a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length NonEmpty (NonEmpty Expr)
coeffs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= NonEmpty Expr -> Int
forall a. NonEmpty a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length NonEmpty 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
$ NonEmpty (NonEmpty Expr) -> [Unknown] -> Bool
isCoeffsMatchUnknowns NonEmpty (NonEmpty 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 = DefinedQuantityDict
-> ConstrConcept
-> NonEmpty (NonEmpty Expr)
-> [Unknown]
-> NonEmpty Expr
-> ConceptChunk
-> DifferentialModel
SystemOfLinearODEs DefinedQuantityDict
indepVar' ConstrConcept
depVar' NonEmpty (NonEmpty Expr)
coeffs [Unknown]
unks NonEmpty Expr
const' (UID -> NP -> Sentence -> ConceptChunk
cncpt''' (String -> UID
mkUid String
id') NP
term' Sentence
defn')

-- | Create a 'DifferentialModel' by the input language
makeASingleDE :: DefinedQuantityDict -> ConstrConcept -> LHS -> Expr-> String -> NP -> Sentence -> DifferentialModel
makeASingleDE :: DefinedQuantityDict
-> ConstrConcept
-> [Term]
-> Expr
-> String
-> NP
-> Sentence
-> DifferentialModel
makeASingleDE DefinedQuantityDict
indepVar'' ConstrConcept
depVar'' [Term]
lhs Expr
const'' String
id'' NP
term'' Sentence
defn''
 | NonEmpty (NonEmpty Expr) -> Int
forall a. NonEmpty a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length NonEmpty (NonEmpty 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
$ NonEmpty (NonEmpty Expr) -> [Unknown] -> Bool
isCoeffsMatchUnknowns NonEmpty (NonEmpty 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 = DefinedQuantityDict
-> ConstrConcept
-> NonEmpty (NonEmpty Expr)
-> [Unknown]
-> NonEmpty Expr
-> ConceptChunk
-> DifferentialModel
SystemOfLinearODEs DefinedQuantityDict
indepVar'' ConstrConcept
depVar'' NonEmpty (NonEmpty Expr)
coeffs [Unknown]
unks (Expr
const'' Expr -> [Expr] -> NonEmpty Expr
forall a. a -> [a] -> NonEmpty a
:| []) (UID -> NP -> Sentence -> ConceptChunk
cncpt''' (String -> UID
mkUid 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 :: NonEmpty (NonEmpty Expr)
coeffs = [Term] -> [Unknown] -> NonEmpty Expr
createCoefficients [Term]
lhs [Unknown]
unks NonEmpty Expr -> [NonEmpty Expr] -> NonEmpty (NonEmpty Expr)
forall a. a -> [a] -> NonEmpty a
:| []

-- | Function to check whether dimension of coefficient is match with the unknown vector
isCoeffsMatchUnknowns :: NonEmpty (NonEmpty Expr) -> [Unknown] -> Bool
isCoeffsMatchUnknowns :: NonEmpty (NonEmpty Expr) -> [Unknown] -> Bool
isCoeffsMatchUnknowns NonEmpty (NonEmpty Expr)
_ [] = String -> Bool
forall a. HasCallStack => String -> a
error String
"Unknowns column vector can not be empty"
isCoeffsMatchUnknowns NonEmpty (NonEmpty Expr)
coeffs [Unknown]
unks = (NonEmpty Expr -> Bool -> Bool)
-> Bool -> NonEmpty (NonEmpty Expr) -> Bool
forall a b. (a -> b -> b) -> b -> NonEmpty a -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\ NonEmpty Expr
x -> Bool -> Bool -> Bool
(&&) (NonEmpty Expr -> Int
forall a. NonEmpty a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length NonEmpty 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 NonEmpty (NonEmpty 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] -> NonEmpty Expr
createCoefficients :: [Term] -> [Unknown] -> NonEmpty Expr
createCoefficients [] [Unknown]
_ = String -> NonEmpty Expr
forall a. HasCallStack => String -> a
error String
"Left hand side is an empty list"
createCoefficients [Term]
_ [] = String -> NonEmpty Expr
forall a. HasCallStack => String -> a
error String
"No unknowns"
createCoefficients [Term]
lhs (Unknown
x:[Unknown]
xs) =
  Maybe Term -> Expr
genCoefficient (Unknown -> [Term] -> Maybe Term
findCoefficient Unknown
x [Term]
lhs) Expr -> [Expr] -> NonEmpty Expr
forall a. a -> [a] -> NonEmpty a
:|
  (Unknown -> Expr) -> [Unknown] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map (\Unknown
z ->  Maybe Term -> Expr
genCoefficient (Unknown -> [Term] -> Maybe Term
findCoefficient Unknown
z [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 [] = String -> [Unknown]
forall a. HasCallStack => String -> a
error String
"impossible, there are always unknowns"
transUnknowns (Unknown
_ : [Unknown]
us) = [Unknown]
us

-- | 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 :: NonEmpty Expr -> [Expr]
transCoefficients :: NonEmpty Expr -> [Expr]
transCoefficients (Expr
e :| [Expr]
es) =
  (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
== Expr
zero then Expr
zero else Expr -> Expr
forall r. ExprC r => r -> r
neg Expr
x Expr -> Expr -> Expr
forall r. ExprC r => r -> r -> r
$/ Expr
e) [Expr]
es
  where zero :: Expr
zero = Unknown -> Expr
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
0

-- | 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]
front [Expr] -> [Expr] -> [Expr]
forall a. [a] -> [a] -> [a]
++ [Expr] -> [Expr]
forall {a}. LiteralC a => [a] -> [a]
ident [Expr]
back
  where
    ([Expr]
front, [Expr]
back) = Int -> [Expr] -> ([Expr], [Expr])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
n [Expr]
es
    ident :: [a] -> [a]
ident [] = String -> [a]
forall a. HasCallStack => String -> a
error String
"second half should not be empty"
    ident (a
_ : [a]
xs) = Unknown -> a
forall r. LiteralC r => Unknown -> r
exactDbl Unknown
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xs

-- | 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"
  | 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 [NonEmpty Expr -> [Expr]
transCoefficients (NonEmpty Expr -> [Expr]) -> NonEmpty Expr -> [Expr]
forall a b. (a -> b) -> a -> b
$ NonEmpty (NonEmpty Expr) -> NonEmpty Expr
forall a. NonEmpty a -> a
NE.head (DifferentialModel
dm DifferentialModel
-> Getting
     (NonEmpty (NonEmpty Expr))
     DifferentialModel
     (NonEmpty (NonEmpty Expr))
-> NonEmpty (NonEmpty Expr)
forall s a. s -> Getting a s a -> a
^. Getting
  (NonEmpty (NonEmpty Expr))
  DifferentialModel
  (NonEmpty (NonEmpty Expr))
Lens' DifferentialModel (NonEmpty (NonEmpty 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 [NonEmpty Expr -> Expr
forall a. NonEmpty a -> a
NE.head (DifferentialModel
dm DifferentialModel
-> Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
-> NonEmpty Expr
forall s a. s -> Getting a s a -> a
^. Getting (NonEmpty Expr) DifferentialModel (NonEmpty Expr)
Lens' DifferentialModel (NonEmpty Expr)
dmConstants) Expr -> Expr -> Expr
`divideConstant` NonEmpty Expr -> Expr
forall a. NonEmpty a -> a
NE.head (NonEmpty (NonEmpty Expr) -> NonEmpty Expr
forall a. NonEmpty a -> a
NE.head (DifferentialModel
dm DifferentialModel
-> Getting
     (NonEmpty (NonEmpty Expr))
     DifferentialModel
     (NonEmpty (NonEmpty Expr))
-> NonEmpty (NonEmpty Expr)
forall s a. s -> Getting a s a -> a
^. Getting
  (NonEmpty (NonEmpty Expr))
  DifferentialModel
  (NonEmpty (NonEmpty Expr))
Lens' DifferentialModel (NonEmpty (NonEmpty 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. (IsChunk c, HasSymbol c) => c -> Expr
forall r c. (ExprC r, IsChunk 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