packages feed

ideas-0.7: src/Domain/Math/Power/Equation/NormViews.hs

-----------------------------------------------------------------------------
-- Copyright 2010, Open Universiteit Nederland. This file is distributed 
-- under the terms of the GNU General Public License. For more information, 
-- see the file "LICENSE.txt", which is included in the distribution.
-----------------------------------------------------------------------------
-- |
-- Maintainer  :  alex.gerdes@ou.nl
-- Stability   :  provisional
-- Portability :  portable (depends on ghc)
--
-----------------------------------------------------------------------------

module Domain.Math.Power.Equation.NormViews
{-   ( normPowerEqApproxView
   , normPowerEqView
   , normExpEqView
   , normLogEqView
--   , normLogView
   ) -} where

import Common.Classes
import Common.View
import Common.Rewriting hiding (rewrite)
import Control.Arrow ( (>>^) )
import Control.Monad
import Data.List
import Data.Maybe
import Data.Ratio
import Domain.Math.Approximation
import Domain.Math.Data.OrList
import Domain.Math.Data.PrimeFactors
import Domain.Math.Data.Relation
import Domain.Math.Expr
import Domain.Math.Numeric.Views
import Domain.Math.Polynomial.CleanUp
import Domain.Math.Polynomial.Views
import Domain.Math.Power.NormViews
import Domain.Math.Power.Utils
import Domain.Math.Power.Views
import Domain.Math.Simplification hiding (simplify, simplifyWith)

import Common.Uniplate

-- Change to configurable strategy!
normPowerEqApproxView :: Int -> View (Relation Expr) (Expr, Expr)
normPowerEqApproxView d = makeView f (uncurry (.~=.))
   where
     f rel = case relationType rel of 
      EqualTo       -> fmap (second (simplifyWith (precision d) doubleView)) 
                     $ match (equationView >>> normPowerEqView) rel 
      Approximately -> return (leftHandSide rel, rightHandSide rel)
      _             -> Nothing

normPowerEqView :: View (Equation Expr) (Expr, Expr) -- with x>0!
normPowerEqView = makeView f (uncurry (:==:))
  where
    f expr = do
      -- selected var to the left, the rest to the right
      (lhs :==: rhs) <- varLeft expr >>= constRight
      -- match power
      (c, ax)        <- match (timesView <&> (identity >>^ (,) 1)) $
                          simplify normPowerView lhs
      (a, x)         <- match myPowerView ax
      -- simplify, scale and take root
      let y = cleanUpExpr $ (rhs ./. c) .^. (1 ./. x)
      return (a, simplify rationalView y)

    myPowerView =  powerView 
               <&> (rootView >>> second (makeView (\a->Just (1 ./. a)) (1 ./.)))
               <&> (identity >>^ \a->(a,1))

normPowerEqView' :: View (Equation Expr) (Expr, Expr) -- with x>0!
normPowerEqView' = makeView f (uncurry (:==:))
  where
    f expr = do
      -- selected var to the left, the rest to the right
      (lhs :==: rhs) <- varLeft expr >>= constRight
      -- match power
      (c, (a, x))    <- match unitPowerView lhs
      -- simplify, scale and take root
      let y = cleanUpExpr $ (rhs ./. c) .^. (1 ./. x)
      return (a, simplify myRationalView y)

constRight :: Equation Expr -> Maybe (Equation Expr)
constRight (lhs :==: rhs) = do
  (vs, cs) <- fmap (partition hasSomeVar) (match sumView lhs)
  let rhs' = rhs .+. build sumView (map neg cs)
  return $ negateEq $ build sumView vs :==: simplifyWith mergeAlikeSum sumView rhs'

negateEq :: Equation Expr -> Equation Expr
negateEq (lhs :==: rhs) = 
  case lhs of
    Negate lhs' -> lhs' :==: neg rhs
    _           -> lhs  :==: rhs

varLeft :: Equation Expr -> Maybe (Equation Expr)
varLeft (lhs :==: rhs) = do
  (vs, cs) <- fmap (partition hasSomeVar) (match sumView rhs)
  return $ lhs .+. build sumView (map neg vs) :==: build sumView cs

scaleLeft :: Equation Expr -> Maybe (Equation Expr)
scaleLeft (lhs :==: rhs) = 
  match timesView lhs >>= \(c, x) -> return $ 
    x :==: simplifyWith (second mergeAlikeProduct) productView (rhs ./. c)

normExpEqView :: View (Equation Expr) (String, Rational)
normExpEqView = makeView f id >>> linearEquationView
  where
    try g a = fromMaybe a $ g a
    f e = do
      let (l :==: r) = try scaleLeft $ try constRight e
      return $ case match powerView l of
        Just (b, x) -> x :==: simplify normLogView (logBase b r)
        Nothing     -> l :==: r

normLogEqView :: View (OrList (Equation Expr)) (OrList (Equation Expr))
normLogEqView = makeView (liftM g . switch . fmap f) id  -- AG: needs to be replaced by higherOrderEqView
  where
    f expr@(lhs :==: rhs) = return $
      case match logView lhs of
        Just (b, x) -> x :==: simplify myRationalView (b .^. rhs)
        Nothing     -> expr
    g = fmap (fmap (simplify myRationalView) . simplify normPowerEqView) 
      . simplify quadraticEquationsView 

-- liftToOrListView :: View a b -> View (OrList a) (OrList b)
-- liftToOrListView v = makeView (switch . fmap (match v)) ()

normLogView :: View Expr Expr
normLogView = makeView g id
  where
    g expr = 
      case expr of 
        Sym s [x, y] 
          | isLogSymbol s -> do
              b <- match integerView x
              let divExp (be, n) = return $ f be y ./. fromInteger n
              maybe (Just $ f b y) divExp $ greatestPower b
          | otherwise -> Nothing
        _ -> Nothing
    f b expr= 
      case expr of
        Nat 1     -> Nat 0
        Nat n     
          | n == b    -> Nat 1
          | otherwise -> maybe (logBase (fromInteger b) (fromInteger n)) Nat 
                       $ lookup b (allPowers n)
        e1 :*: e2 -> f b e1 .+. f b e2
        e1 :/: e2 -> f b e1 .-. f b e2
        Sqrt e    -> f b (e .^. (1 ./. 2))
        Negate e  -> Negate $ f b e
        Sym s [x,y]
          | isPowerSymbol s -> y .*. f b x
          | isRootSymbol  s -> f b (x .^. (1 ./. y))
        _         -> expr

myRationalView :: View Expr Rational
myRationalView = makeView (return . rewrite simplerPower) id >>> rationalView

simplerPower :: Expr -> Maybe Expr
simplerPower expr = 
  case expr of      
    Sqrt x -> simplerPower $ x .^. (1/2)
    Sym s [x, y]
      | isRootSymbol s  -> simplerPower $ x .^. (1/y)
      | isPowerSymbol s -> f
      | otherwise -> Nothing
        where f | y == 0 || x == 1 = Just 1
                | y == 1 = Just x
                | x == 0 = Just 0
                | otherwise =
                  -- geheel getal
                  liftM fromRational (match rationalView expr) 
                  `mplus`
                  -- wortel
                  do 
                    ry <- match rationalView y
                    rx <- match rationalView x
                    guard $ numerator ry == 1 && denominator rx == 1
                    liftM fromInteger $ takeRoot (numerator rx) (denominator ry)
                  `mplus`
                  -- (a/b)^y -> a^x/b^y
                  do
                    (a, b) <- match divView x
                    return $ build divView (a .^. y, b .^. y)
    _ -> Nothing

-- myRationalView = makeView (exprToNum f) id >>> rationalView
--   where
--     f s [x, y] 
--       | isDivideSymbol s = 
--           fracDiv x y
--       | isPowerSymbol s = do
--           ry <- match rationalView y
--           rx <- match rationalView x
--           if      ry == 0 then return 1                      -- 0
--           else if ry == 1 then return rx                     -- 1
--           else if denominator ry == 1 then            -- geheel getal
--             let a = x Prelude.^ abs (numerator ry)
--             in return $ if numerator ry < 0 then 1 / a else a
--           else if numerator ry == 1 then              -- breuk / root
--             if denominator ry > 1 then 
--               if denominator rx == 1 then
--                 takeRoot (numerator rx) (denominator ry)   -- breuk/root
--               else
--                 f powerSymbol [numerator rx, ] / f powerSymbol []
--             else
--               take
--           else                                       -- no calculation
--       | isRootSymbol s = do
--           n <- match integerView y
--           b <- match integerView x
--           liftM fromInteger $ lookup b $ map swap (allPowers n)
--     f _ _ = Nothing