```{-|
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 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
--import qualified Data.Map as Map

-- the following is code for unit testing
{-

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,
DomainIntBox box varid ira,
DomainBoxMappable box box varid ira ira) =>
[(box, ira)]
{-^ the equations;
each equation has coefficients of linear terms
+ constant term -} ->
box {-^ the domain of the variables -} ->
ira {-^ an upper bound on the size of an acceptable solution box -} ->
Maybe box
{-^
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
| 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

{-|
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 = sum \$ DBox.elems \$ DBox.intersectionWith (*) expr box

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

```