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

    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