module Data.Number.ER.Real.Arithmetic.LinearSolver
(
linearSolver
)
where
import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes
import Data.List
import Data.Maybe
import qualified Data.Map as Map
linearSolver ::
(RA.ERIntApprox ira) =>
[(Map.Map VarID ira, ira)]
->
Box ira ->
ira ->
Maybe (Box ira)
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
| width (b Map.! (widestVar b)) < tolerance =
Just b
| otherwise =
linearSolver' eqns (splitBox b ++ bs) tolerance
evalEqns box eqns =
and $ map (evalEqn box) eqns
evalEqn box (expr,cons) =
cons `RA.refines` (evalExpr expr box)
where
evalExpr expr box = Map.fold (+) 0 $ Map.unionWith (*) expr box
splitBox box =
[Map.insert k (iLg RA.\/ iMg) box, Map.insert k (iMg RA.\/ iRg) box]
where
iMg = (iLg+iRg)/2
iLg = incrementGranularity iL
iRg = incrementGranularity iR
(iL,iR) = RA.bounds i
i = box Map.! k
k = widestVar box
incrementGranularity x =
RA.setMinGranularity (RA.getGranularity x + 1) x
widestVar box =
fst $ maxm (head widthMap) $ tail widthMap
where
widthMap =
Map.assocs $ Map.map width box
maxm m [] = m
maxm m (x:xs) =
if mW < xW then
maxm x xs
else
maxm m xs
where
mW = snd m
xW = snd x
width i =
snd $ RA.bounds (iRiL)
where
(iL,iR) = RA.bounds i