module Data.Number.ER.Real.Arithmetic.LinearSolver
(
linearSolver
)
where
import qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.List
import Data.Maybe
linearSolver ::
(RA.ERIntApprox ira,
DomainIntBox box varid ira,
DomainBoxMappable box box varid ira ira) =>
[(box, ira)]
->
box ->
ira ->
Maybe box
linearSolver eqns domBox tolerance =
linearSolver' eqns [domBox] tolerance
linearSolver' eqns [] tolerance =
Nothing
linearSolver' eqns (b:bs) tolerance
| not $ evalEqns b eqns =
linearSolver' eqns bs tolerance
| belowTolerance =
Just b
| otherwise =
linearSolver' eqns (splitBox b ++ bs) tolerance
where
belowTolerance =
and $ map (\d -> width d `RA.ltSingletons` tolerance) $ DBox.elems b
evalEqns box eqns =
and $ map (evalEqn box) eqns
evalEqn box (expr,cons) =
cons `RA.refines` (evalExpr expr box)
where
evalExpr expr box = sum $ DBox.elems $ DBox.intersectionWith (*) expr box
splitBox box =
[DBox.insert k (iLg RA.\/ iMg) box,
DBox.insert k (iMg RA.\/ iRg) box]
where
iMg = (iLg+iRg)/2
iLg = incrementGranularity iL
iRg = incrementGranularity iR
(iL,iR) = RA.bounds i
i = DBox.lookup "ER: LinearSolver: splitBox: " k box
k = widestVar box
incrementGranularity x =
RA.setMinGranularity (RA.getGranularity x + 1) x
widestVar box =
fst $ DBox.bestSplit box
width i =
snd $ RA.bounds (iRiL)
where
(iL,iR) = RA.bounds i