module Domain.Math.Polynomial.BalanceUtils
( eqView, minusView, negView
, matchLin, matchPlusCon
, cleaner, cleanerExpr
, linbal, checkForChange
, termRef, factorRef, factor1Ref, factor2Ref
, writeTerm, writeWith
, buggyBalanceRule, buggyBalanceRuleArg
, buggyBalanceExprRule
, buggyBalanceRecognizer
, collectLocal, collectGlobal
, distributeDiv, distributeTimes
, isPlusT, diffPlus
, isTimesT, diffTimes
) where
import Control.Monad
import Data.List
import Data.Maybe
import Domain.Math.Data.Polynomial
import Domain.Math.Data.Relation
import Domain.Math.Data.WithBool
import Domain.Math.Expr
import Domain.Math.Numeric.Views
import Domain.Math.Polynomial.Views
import Domain.Math.Safe
import Domain.Math.Simplification (mergeAlikeSum)
import Ideas.Common.Library
import Ideas.Utils.Prelude (fixpoint)
import Ideas.Utils.Uniplate
eqView :: View (WithBool (Equation Expr)) (WithBool (String, Rational))
eqView = makeView (either (Just . fromBool) f . fromWithBool) (fmap g)
where
f (lhs :==: rhs) = do
(s, p) <- match (polyViewWith rationalApproxView) (lhs-rhs)
case degree p of
0 -> Just $ fromBool $ coefficient 0 p == 0
1 -> Just $ singleton (s, - coefficient 0 p / coefficient 1 p)
_ -> Nothing
g (s, r) = Var s :==: fromRational r
minusView :: View Expr (Expr, Expr)
minusView = makeView isMinus (uncurry (:-:))
negView :: View Expr Expr
negView = makeView isNegate Negate
matchLin :: MonadPlus m => Expr -> m (Expr, Rational, Rational)
matchLin expr = do
(s, p) <- matchM (polyNormalForm rationalView) expr
guard (degree p == 1)
return (Var s, coefficient 1 p, coefficient 0 p)
matchPlusCon :: MonadPlus m => Expr -> m (Expr, Rational)
matchPlusCon expr =
matchM (plusView >>> second rationalView) expr
`mplus`
matchM (plusView >>> toView swapView >>> second rationalView) expr
cleaner :: WithBool (Equation Expr) -> WithBool (Equation Expr)
cleaner = fmap (fmap cleanerExpr)
cleanerExpr :: Expr -> Expr
cleanerExpr = transform f
where
f (a :/: Nat 1) = f a
f (a :/: Negate (Nat 1)) = f $ Negate a
f (Negate a :/: Negate b) = f (a/b)
f (a :/: Negate b) = f $ Negate (a/b)
f (Negate a :/: b) = f $ Negate (a/b)
f (Negate (Negate a)) = f a
f e = cleanSum (cleanProduct (simplify rationalView e))
cleanSum =
let g x y = canonical rationalView (x :+: y)
in simplifyWith (adjacent g) simpleSumView
cleanProduct =
let g x y = canonical rationalView (x :*: y)
reorder = uncurry (++) . partition (`belongsTo` rationalView)
in simplifyWith (mapSecond (adjacent g . reorder)) simpleProductView
adjacent :: (a -> a -> Maybe a) -> [a] -> [a]
adjacent f = rec
where
rec (x:y:rest) =
case f x y of
Just xy -> rec (xy:rest)
Nothing -> x:rec (y:rest)
rec xs = xs
termRef, factorRef, factor1Ref, factor2Ref :: Ref Expr
termRef = makeRef "term"
factorRef = makeRef "factor"
factor1Ref = makeRef "factor1"
factor2Ref = makeRef "factor2"
writeWith :: Trans c () -> (a -> Maybe (b, c)) -> Trans a b
writeWith wf f = transMaybe f >>> second wf >>> arr fst
writeTerm :: (a -> Maybe (b, Expr)) -> Trans a b
writeTerm = writeWith (writeRef_ termRef)
linbal :: Id
linbal = newId "algebra.equations.linear.balance"
bugbal :: IsId n => n -> Id
bugbal n = newId (linbal, "buggy", n)
checkForChange :: (MonadPlus m, Eq a) => (a -> m a) -> a -> m a
checkForChange f a = f a >>= \b -> guard (a /= b) >> return b
buggyBalanceRule :: IsId n => n -> (Equation Expr -> Maybe (Equation Expr)) -> Rule (Equation Expr)
buggyBalanceRule n = addTransRecognizer eq . buggyRule (bugbal n)
where
eq = viewEquivalent (traverseView (polyViewWith rationalView))
buggyBalanceRuleArg :: IsId n => n -> Trans (Equation Expr) (Equation Expr) -> Rule (Equation Expr)
buggyBalanceRuleArg n = addTransRecognizer eq . buggy . ruleTrans (bugbal n)
where
eq = viewEquivalent (traverseView (polyViewWith rationalView))
buggyBalanceExprRule :: IsId n => n -> (Expr -> Maybe Expr) -> Rule Expr
buggyBalanceExprRule = buggyRule . bugbal
buggyBalanceRecognizer :: IsId n => n -> Trans (a, a) () -> Rule a
buggyBalanceRecognizer n t =
addRecognizer (makeRecognizerTrans t) $ buggy $ emptyRule (bugbal n)
collectLocal :: Expr -> Expr
collectLocal = simplifyWith (mapSecond f) simpleProductView
. simplifyWith mergeAlikeSum simpleSumView
where
f xs | length ys > 1 = ys++zs
| otherwise = xs
where
(ys, zs) = partition hasNoVar xs
collectGlobal :: Expr -> Expr
collectGlobal = fixpoint (transform collectLocal)
distributeDiv :: Expr -> Expr
distributeDiv expr = fromMaybe expr $ do
(a, r) <- match (divView >>> second rationalView) expr
return $ simplifyWith (fmap (`divide` r)) simpleSumView a
where
divide x r = fromMaybe (x/fromRational r) $ do
(y, z) <- match (timesView >>> first rationalView) x
new <- y `safeDiv` r
return (fromRational new * z)
`mplus` do
(y, z) <- match (timesView >>> second rationalView) x
new <- z `safeDiv` r
return (y * fromRational new)
distributeTimes :: Expr -> Expr
distributeTimes expr = fromMaybe expr $ do
(r, a) <- match (timesView >>> first rationalView) expr
`mplus`
match (timesView >>> second rationalView >>> toView swapView) expr
return $ simplifyWith (fmap (times r)) simpleSumView a
where
times r x = fromMaybe (fromRational r*x) $ do
(a, b) <- match (divView >>> second rationalView) x
guard (b /= 0)
return (fromRational (r/b) * a)
isPlusT :: Equation Expr -> Equation Expr -> Bool
isPlusT old new = isJust (diffPlusEq old new)
diffPlusEq :: Equation Expr -> Equation Expr -> Maybe Expr
diffPlusEq (a1 :==: a2) (b1 :==: b2) = do
d1 <- diffPlus a1 b1
d2 <- diffPlus a2 b2
guard (d1 == d2)
return d1
diffPlus :: Expr -> Expr -> Maybe Expr
diffPlus a b = do
let myView = polyViewWith rationalView
(x, pa) <- matchM myView a
(y, pb) <- matchM myView b
guard (x==y)
let d = pb - pa
return $ build myView (x, d)
isTimesT :: Equation Expr -> Equation Expr -> Bool
isTimesT old new = isJust (diffTimesEq old new)
diffTimesEq :: Equation Expr -> Equation Expr -> Maybe Expr
diffTimesEq (a1 :==: a2) (b1 :==: b2) = do
d1 <- diffTimes a1 b1
d2 <- diffTimes a2 b2
guard (d1 == d2)
return d1
diffTimes :: MonadPlus m => Expr -> Expr -> m Expr
diffTimes a b = do
let myView = polyViewWith rationalView
(x, pa) <- matchM myView a
(y, pb) <- matchM myView b
guard (x==y)
if pa==0 && pb==0 then return 1 else do
d <- maybe (fail "diffTimes") return (pb `safeDiv` pa)
return $ build myView (x, d)