```{-|
Module      :  Data.Number.ER.Real.LinearSolver
Description :  arbitrary precision piece-wise something function enclosures
Copyright   :  (c) Jan Duracz, Michal Konecny

Maintainer  :  mik@konecny.aow.cz
Stability   :  experimental
Portability :  portable

A simple validated solver for systems of linear equations with
interval coefficients.  It uses a naive splitting approach and is
therefore very slow.
-}
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

{-
import Data.Number.ER.Real.DefaultRepr

eq1 :: (Box IRA, IRA)
eq1 = (mkBox [1,2,1], 2)
eq2 = (mkBox [2,4,2], 4)
eq3 = (mkBox [2,4,4], 5)
eqs = [eq1,eq2,eq3]

box = mkBox \$ replicate 3 \$ (-1)RA.\/1
x1 = (-13/16)RA.\/(-3/4) :: IRA
x2 = (5/16)RA.\/(3/8) :: IRA
tol = 2^^(-20) :: IRA

mkBox :: [IRA] -> Box IRA
mkBox iras = Map.fromList \$ zip [1..] iras
-}

linearSolver ::
(RA.ERIntApprox ira) =>
[(Map.Map VarID ira, ira)]
{-^ the equations;
each equation has coefficients of linear terms
+ constant term -} ->
Box ira {-^ the domain of the variables -} ->
ira {-^ an upper bound on the size of an acceptable solution box -} ->
Maybe (Box ira)
{-^
A box containing at least one solution within the domain;
Nothing if there is no solution.
-}
linearSolver eqns domBox tolerance =
linearSolver' eqns [domBox] tolerance

linearSolver' eqns [] tolerance =
Nothing
linearSolver' eqns (b:bs) tolerance
| not \$ evalEqns b eqns = -- no solutions in the box
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

{-|
returns true iff there exists a solution to the equation in the box
-}
evalEqn box (expr,cons) =
cons `RA.refines` (evalExpr expr box)
where
evalExpr expr box = Map.fold (+) 0 \$ Map.unionWith (*) expr box

{-|
returns the list of (two) boxes resulting from splitting the widest edge
of the box in half
-}
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 (iR-iL)
where
(iL,iR) = RA.bounds i

```