packages feed

mfsolve (empty) → 0.1.0

raw patch · 5 files changed

+1016/−0 lines, 5 filesdep +basedep +hashabledep +mfsolvesetup-changed

Dependencies added: base, hashable, mfsolve, tasty, tasty-hunit, unordered-containers

Files

+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) 2014, Kristof Bastiaensen+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++* Redistributions of source code must retain the above copyright notice, this+  list of conditions and the following disclaimer.++* Redistributions in binary form must reproduce the above copyright notice,+  this list of conditions and the following disclaimer in the documentation+  and/or other materials provided with the distribution.++* Neither the name of the {organization} nor the names of its+  contributors may be used to endorse or promote products derived from+  this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Math/MFSolve.hs view
@@ -0,0 +1,766 @@+{-# LANGUAGE DeriveGeneric, PatternGuards, ViewPatterns #-}++{-|+Module      : Math.MFSolve+Description : Equation solver and calculator à la metafont+Copyright   : (c) Kristof Bastiaensen, 2015+License     : BSD-3+Maintainer  : kristof@resonata.be+Stability   : unstable+Portability : portable++This module implements an equation solver that solves and evaluates+expressions on the fly.  It is based on Prof. D.E.Knuth's+/metafont/.  The goal of mfsolve is to make the solver useful in an+interactive program, by enhancing the bidirectionality of the solver.+Like metafont, it can solve linear equations, and evaluate nonlinear+expressions.  In addition to metafont, it also solves for angles, and+makes the solution independend of the order of the equations.++The `Expr` datatype allows for calculations with constants and unknown+variables.  The `Dependencies` datatype contains all dependencies and known equations.++=== Examples:++Let's define some variables.  The `SimpleVar` type is a simple wrapper+around `String` to provide nice output.++> let [x, y, t, a] = map (makeVariable . SimpleVar) ["x", "y", "t", "a"]++Solve linear equations:++> showVars $ solveEqs emptyDeps+> [ 2*x + y === 5,+>   x - y   === 1]++> x = 2.0+> y = 1.0++Solve for angle (pi/4):++> showVars $ solveEqs emptyDeps+> [ sin(t) === 1/sqrt(2) ]++> t = 0.7853981633974484++Solve for angle (pi/3) and amplitude:++> showVars $ solveEqs emptyDeps+> [ a*sin(x) === sqrt 3,+>   a*cos(x) === 1 ]++> x = 1.0471975511965979+> a = 2.0++Allow nonlinear expression with unknown variables:++> showVars $ solveEqs emptyDeps+> [ sin(sqrt(x)) === y,+>   x === 2]++>x = 2.0+>y = 0.9877659459927355++Find the angle and amplitude when using a rotation matrix:++> showVars $ solveEqs emptyDeps+> [ a*cos t*x - a*sin t*y === 30,+>   a*sin t*x + a*cos t*y === 40,+>   x === 10,+>   y === 10 ]++> x = 10.0+> y = 10.0+> t = 0.14189705460416402+> a = 3.5355339059327373++-}++module Math.MFSolve+       (SimpleExpr(..), Expr, LinExpr(..), UnaryOp(..), BinaryOp(..),+        Dependencies, DepError(..), SimpleVar(..),+        getKnown, knownVars, varDefined, nonlinearEqs, dependendVars,+        simpleExpr, emptyDeps, makeVariable, makeConstant,+        (===), (=&=), solveEqs, showVars)+where+import qualified Data.HashMap.Strict as M+import qualified Data.HashSet as H+import GHC.Generics+import Data.Hashable+import Data.Maybe+import Data.List+import Data.Function(on)+import Control.Monad++infixr 1 === , =&=++--  | _labeled_ black box mathematical functions+data UnaryFun n = UnaryFun UnaryOp (n -> n)++data BinaryOp = Add | Mul+              deriving Eq+data UnaryOp =+  Sin | Abs | Recip | Signum |+  Exp | Log | Cos | Cosh | Atanh |+  Tan | Sinh | Asin | Acos | Asinh | Acosh | Atan+  deriving (Eq, Generic)++-- | A simplified datatype representing an expression+data SimpleExpr v n =+  SEBin BinaryOp (SimpleExpr v n) (SimpleExpr v n) |+  SEUn UnaryOp (SimpleExpr v n) |+  Var v |+  Const n++newtype SimpleVar = SimpleVar String+                  deriving (Eq, Ord, Generic)++-- | An mathematical expression of several variables.+data Expr v n = Expr (LinExpr v n) [TrigTerm v n] [NonLinExpr v n]+                deriving Generic++-- | A linear expression of several variables.+-- For example: @2*a + 3*b + 2@ would be represented as+-- @LinExpr 2 [(a, 2), (b, 3)]@.+data LinExpr v n = LinExpr n [(v, n)]+                 deriving (Generic, Eq, Show)+type Period v n = [(v, n)]+type Phase n = n+type Amplitude v n = LinExpr v n++-- A sum of sinewaves with the same period (a linear sum of several+-- variables), but possibly different (constant) phase.  For example+-- @(2*a+b) sin (x+y) + 2*b*sin(x+y+pi)@ would be represented by:+-- @TrigTerm [(x,1),(y,1)] [(0, LinExpr 0 [(a, 2), (b, 1)]),+-- (pi, LinExpr 0 [(b, 2)])@+type TrigTerm v n = (Period v n, [(Phase n, Amplitude v n)])++-- Any other term+data NonLinExpr v n = +  UnaryApp (UnaryFun n) (Expr v n) |+  MulExp (Expr v n) (Expr v n) |+  SinExp (Expr v n)+  deriving Generic++-- | An angular function of the form @c + n*sin(theta + alpha)@+-- where @theta@, and @n@ are linear terms, @alpha@ and @c@ are constants.+type LinearMap v n = M.HashMap v (LinExpr v n)+type TrigEq v n = (Period v n, Amplitude v n, Phase n, n)+type TrigEq2 v n = M.HashMap (Period v n)+                   (M.HashMap v (Expr v n))++instance (Hashable v, Hashable n) => Hashable (LinExpr v n)+instance (Hashable v, Hashable n) => Hashable (NonLinExpr v n)+instance (Hashable n) => Hashable (UnaryFun n) where+  hashWithSalt s (UnaryFun o _) = hashWithSalt s o+instance Hashable UnaryOp+instance (Hashable v, Hashable n) => Hashable (Expr v n)+instance Hashable SimpleVar++-- | A simple String wrapper, which will print formulas more cleanly.+instance Show SimpleVar where+  show (SimpleVar s) = s++-- | An opaque datatype containing the dependencies of each variable.+-- A variable who's dependency is just a number is called /known/.  A+-- variables which depends on other variables is called /dependend/.+-- A variable which is neither known or dependend is called+-- /independend/.  A variable can only depend on other /independend/+-- variables.  It also contains nonlinear equations which it couldn't+-- reduce to a linear equation yet.+data Dependencies v n = Dependencies+                        (M.HashMap v (H.HashSet v))+                        (LinearMap v n)+                        [TrigEq v n]+                        (TrigEq2 v n)+                        [Expr v n]+                        +-- | An error type for '===', '=&=' and 'solveEq':+data DepError n =+  -- | 'InconsistentEq' @a@: The equation was reduced to the+  -- impossible equation `a == 0` for nonzero a, which means the+  -- equation is inconsistent with previous equations.+  InconsistentEq n |+  -- | 'RedundantEq': The equation was reduced to the redundant equation 0 == 0, which+  -- means it doesn't add any information.+  RedundantEq++instance (Ord n, Num n, Eq n, Show v, Show n) => Show (Expr v n) where+  show e = show (simpleExpr e)++withParens :: (Show t1, Show t, Ord t1, Num t1, Eq t1) => SimpleExpr t t1 -> [BinaryOp] -> String+withParens e@(SEBin op _ _) ops+  | op `elem` ops = "(" ++ show e ++ ")"+withParens e _ = show e++instance (Show v, Ord n, Show n, Num n, Eq n) => Show (SimpleExpr v n) where+  show (Var v) = show v+  show (Const n) = show n+  show (SEBin Add e1 (SEBin Mul (Const e2) e3))+    | e2 < 0 =+      show e1 ++ " - " ++ show (SEBin Mul (Const (negate e2)) e3)+  show (SEBin Add e1 e2) =+    show e1 ++ " + " ++ show e2+  show (SEBin Mul (Const 1) e) = show e+  show (SEBin Mul e (Const 1)) = show e+  show (SEBin Mul e1 (SEUn Recip e2)) =+    withParens e1 [Add] ++ "/" ++ withParens e2 [Add, Mul]+  show (SEBin Mul e1 e2) =+    withParens e1 [Add] ++ "*" ++ withParens e2 [Add]+  show (SEUn Exp (SEBin Mul (SEUn Log e1) e2)) =+    withParens e1 [Add, Mul] ++ "**" ++ withParens e2 [Add, Mul]+  show (SEUn op e) = show op ++ "(" ++ show e ++ ")"++instance Show BinaryOp where+  show Add = "+"+  show Mul = "*"++instance Show UnaryOp where+  show Sin = "sin"+  show Abs = "abs"+  show Recip = "1/"+  show Signum = "sign"+  show Exp = "exp"+  show Log = "log"+  show Cos = "cos"+  show Cosh = "cosh"+  show Atanh = "atanh"+  show Tan = "tan"+  show Sinh = "sinh"+  show Asin = "asin"+  show Acos = "acos"+  show Asinh = "asinh"+  show Acosh = "acosh"+  show Atan = "atan"++instance (Floating n, Ord n, Ord v) => Num (Expr v n) where+  (+) = addExpr+  (*) = mulExpr+  negate = mulExpr (makeConstant (-1))+  abs = unExpr (UnaryFun Abs abs)+  signum = unExpr (UnaryFun Signum signum)+  fromInteger = makeConstant . fromInteger++instance (Floating n, Ord n, Ord v) => Fractional (Expr v n) where+  recip = unExpr (UnaryFun Recip (1.0/))+  fromRational = makeConstant . fromRational++instance (Floating n, Ord n, Ord v) => Floating (Expr v n) where+  pi = makeConstant pi+  exp = unExpr (UnaryFun Exp exp)+  log = unExpr (UnaryFun Log log)+  sin = sinExpr+  cos a = sinExpr (a + makeConstant (pi/2))+  cosh = unExpr (UnaryFun Cosh cosh)+  atanh = unExpr (UnaryFun Atanh atanh)+  tan = unExpr (UnaryFun Tan tan)+  sinh = unExpr (UnaryFun Sinh sinh)+  asin = unExpr (UnaryFun Asin asin)+  acos = unExpr (UnaryFun Acos acos)+  asinh = unExpr (UnaryFun Asinh asinh)+  acosh = unExpr (UnaryFun Acosh acosh)+  atan = unExpr (UnaryFun Atan atan)++instance (Show n, Floating n, Ord n, Ord v, Show v) =>Show (Dependencies v n) where+  show dep@(Dependencies _ lin _ _ _) = +    unlines (map showLin (M.toList lin) +++             map showNl (nonlinearEqs dep))+    where showLin (v, e) = show v ++ " = " ++ show (linExpr e)+          showNl e = show e ++ " = 0"++instance (Show n) => Show (DepError n) where+  show (InconsistentEq a) =+    "Inconsistent equations, off by " ++ show a+  show RedundantEq =+    "Redundant Equation."+++addSimple :: (Num t1, Eq t1) => SimpleExpr t t1 -> SimpleExpr t t1 -> SimpleExpr t t1+addSimple (Const 0) e = e+addSimple e (Const 0) = e+addSimple e1 e2 = SEBin Add e1 e2++linToSimple :: (Num t1, Eq t1) => LinExpr t t1 -> SimpleExpr t t1+linToSimple (LinExpr v t) =+  Const v `addSimple`+  foldr (addSimple.mul) (Const 0) t+  where+    mul (v2, 1) = Var v2+    mul (v2, c) = SEBin Mul (Const c) (Var v2)+    +trigToSimple :: (Num n, Eq n) => TrigTerm v n -> SimpleExpr v n+trigToSimple (theta, t) =+  foldr (addSimple.makeSin) (Const 0) t+  where+    makeSin (alpha, n) =+      SEBin Mul (linToSimple n)+      (SEUn Sin angle) where+        angle | alpha == 0 =+                linToSimple (LinExpr 0 theta)+              | otherwise =+                SEBin Add (linToSimple (LinExpr 0 theta))+                (Const alpha)++nonlinToSimple :: (Num n, Eq n) => NonLinExpr v n -> SimpleExpr v n+nonlinToSimple (UnaryApp (UnaryFun o _) e) =+  SEUn o (simpleExpr e)+nonlinToSimple (MulExp e1 e2) =+  SEBin Mul (simpleExpr e1) (simpleExpr e2)+nonlinToSimple (SinExp e) =+  SEUn Sin (simpleExpr e)++-- | Convert an `Expr` to a `SimpleExpr`.+simpleExpr :: (Num n, Eq n) => Expr v n -> SimpleExpr v n+simpleExpr (Expr lin trig nonlin) =+  linToSimple lin `addSimple`+  foldr (addSimple.trigToSimple)+  (Const 0) trig `addSimple`+  foldr (addSimple.nonlinToSimple)+  (Const 0) nonlin++zeroTerm :: (Num n) => LinExpr v n+zeroTerm = LinExpr 0 []++linExpr :: LinExpr v n -> Expr v n+linExpr lt = Expr lt [] []++-- | Create an expression from a constant+makeConstant :: n -> Expr v n+makeConstant c = linExpr (LinExpr c [])++-- | Create an expression from a variable+makeVariable :: Num n => v -> Expr v n+makeVariable v = linExpr (LinExpr 0 [(v, 1)])++trigExpr :: (Num n) => [TrigTerm v n] -> Expr v n+trigExpr t = Expr zeroTerm t []++nonlinExpr :: Num n => [NonLinExpr v n] -> Expr v n+nonlinExpr = Expr zeroTerm []++getConst :: LinExpr t a -> Maybe a+getConst (LinExpr a []) = Just a+getConst _ = Nothing++getLin :: Expr t n -> Maybe (LinExpr t n)+getLin (Expr lt [] []) = Just lt+getLin _ = Nothing++getConstExpr :: Expr t b -> Maybe b+getConstExpr e = getLin e >>= getConst++addLin :: (Ord v, Num n, Eq n) => LinExpr v n -> LinExpr v n -> LinExpr v n+addLin (LinExpr c1 terms1) (LinExpr c2 terms2) =+  LinExpr (c1+c2) terms3 where+    terms3 = filter ((/= 0) . snd) $+             merge terms1 terms2 (+)++addExpr :: (Ord n, Ord v, Floating n) => Expr v n -> Expr v n -> Expr v n+addExpr (Expr lt1 trig1 nl1) (Expr lt2 trig2 nl2) =+  Expr (addLin lt1 lt2) trig3 (nl1++nl2)+  where+    trig3 = merge trig1 trig2 addTrigTerms++-- merge two association lists, by combining equal keys with+-- the given function, and keeping keys sorted.+merge :: Ord k => [(k, v)] -> [(k, v)] -> (v -> v -> v) -> [(k, v)]+merge [] l _ = l+merge l [] _ = l+merge (a@(k,v):as) (b@(k2,v2):bs) f = case compare k k2 of+  LT -> a: merge as (b:bs) f+  EQ -> (k, f v v2): merge as bs f+  GT -> b: merge (a:as) bs f++-- add trigonometric terms with the same period+addTrigTerms :: (Ord a, Ord t, Floating a) => [(a, LinExpr t a)] -> [(a, LinExpr t a)] -> [(a, LinExpr t a)]+addTrigTerms [] p = p+addTrigTerms terms terms2 =+  foldr mergeTerms terms terms2+  where+    mergeTerms (alpha, n) ((beta, m):rest) =+      case addTrigTerm alpha n beta m of+       Just (_, LinExpr 0 []) -> rest+       Just (gamma, o) ->+         mergeTerms (gamma, o) rest+       Nothing -> (beta, m) : mergeTerms (alpha, n) rest+    mergeTerms a [] = [a]++addTrigTerm :: (Ord a, Ord t, Floating a) => a -> LinExpr t a -> a -> LinExpr t a -> Maybe (a, LinExpr t a)+addTrigTerm alpha n beta m+  | alpha == beta =+    Just (alpha, addLin n m)+  | Just r <- termIsMultiple n m =+      let gamma = atan (divident/divisor) ++                  (if divisor < 0 then pi else 0)+          divident = r*sin alpha + sin beta+          divisor = r*cos alpha + cos beta+          o = sqrt(divident*divident + divisor*divisor)+      in Just (gamma, mulLinExpr o m)+  | otherwise = Nothing++-- compare if the linear term is a multiple of the other, within roundoff                +termIsMultiple :: (Ord a, Fractional a, Eq t) => LinExpr t a -> LinExpr t a -> Maybe a+termIsMultiple (LinExpr _ _) (LinExpr 0 []) = Nothing+termIsMultiple (LinExpr 0 []) (LinExpr _ _) = Nothing+termIsMultiple (LinExpr 0 r1@((_, d1):_)) (LinExpr 0 r2@((_, d2):_))+  | compareBy r1 r2 (compareTerm (d1/d2)) =+      Just (d1/d2)+termIsMultiple (LinExpr c1 r1) (LinExpr c2 r2)+  | compareBy r1 r2 (compareTerm (c1/c2)) =+      Just (c1/c2)+  | otherwise = Nothing++compareTerm :: (Ord a1, Fractional a1, Eq a) => a1 -> (a, a1) -> (a, a1) -> Bool+compareTerm ratio (v3,c3) (v4, c4) = +  v3 == v4 && (abs(c3 - (c4 * ratio)) <= abs c3*2e-50)++compareBy :: [a] -> [b] -> (a -> b -> Bool) -> Bool+compareBy [] [] _ = True+compareBy (e:l) (e2:l2) f =+  f e e2 && compareBy l l2 f+compareBy _ _ _ = False+        +-- multiply a linear term by a constant.+mulLinExpr :: Num n => n -> LinExpr v n -> LinExpr v n+mulLinExpr x (LinExpr e terms) =+  LinExpr (e*x) $ map (fmap (*x)) terms++-- multiply all sines with the constant+-- constant multiplier+mulConstTrig :: (Ord n, Num n) => n -> TrigTerm v n -> TrigTerm v n+mulConstTrig c (theta, terms) =  (theta, tt) where+  tt = map (fmap (mulLinExpr c)) terms++mulLinTrig :: (Ord n, Ord v, Floating n) => LinExpr v n -> TrigTerm v n -> Expr v n+mulLinTrig lt (theta, terms) =+  -- linear multiplier+  foldr ((+).mul1) 0 terms+  where+    -- constant amplitude+    mul1 (alpha, LinExpr c []) =+      trigExpr [(theta, [(alpha, mulLinExpr c lt)])]+    -- linear amplitude+    mul1 t =+      nonlinExpr [MulExp (trigExpr [(theta, [t])])+                      (Expr lt [] [])]++-- constant * (linear + trig)+mulExpr :: (Ord a, Ord t, Floating a) => Expr t a -> Expr t a -> Expr t a+mulExpr (getConstExpr -> Just c) (Expr lt2 trig []) =+  Expr (mulLinExpr c lt2)+  (map (mulConstTrig c) trig) []++mulExpr (Expr lt2 trig []) (getConstExpr -> Just c) =+  Expr (mulLinExpr c lt2)+  (map (mulConstTrig c) trig) []++-- linear * (constant + trig)+mulExpr (getLin -> Just lt) (Expr (getConst -> Just c) trig []) =+  linExpr (mulLinExpr c lt) ++  foldr ((+).mulLinTrig lt) 0 trig++mulExpr (Expr (getConst -> Just c) trig []) (getLin -> Just lt) =+  linExpr (mulLinExpr c lt) ++  foldr ((+).mulLinTrig lt) 0 trig++-- anything else+mulExpr e1 e2 = nonlinExpr [MulExp e1 e2]+      +sinExpr :: Floating n => Expr v n -> Expr v n+sinExpr (Expr (LinExpr c t) [] [])+  | null t = makeConstant (sin c)+  | otherwise = trigExpr [(t, [(c, LinExpr 1 [])])]+sinExpr e = nonlinExpr [SinExp e]++unExpr :: Num n => UnaryFun n -> Expr v n -> Expr v n+unExpr (UnaryFun _ f) e+  | Just c <- getConstExpr e = makeConstant (f c)+unExpr f e = nonlinExpr [UnaryApp f e]++substVarLin :: (Ord v, Num n, Eq n) => (v -> Maybe (LinExpr v n)) -> LinExpr v n -> LinExpr v n+substVarLin s (LinExpr a terms) =+  let substOne (v, c) =+        maybe (LinExpr 0 [(v, c)]) (mulLinExpr c) (s v)+  in foldr (addLin.substOne) (LinExpr a []) terms++substVarNonLin :: (Ord n, Ord v, Floating n) => (v -> Maybe (LinExpr v n)) -> NonLinExpr v n -> Expr v n+substVarNonLin s (UnaryApp f e1) =+  unExpr f (subst s e1)+substVarNonLin s (MulExp e1 e2) =+  subst s e1 * subst s e2+substVarNonLin s (SinExp e1) =+  sin (subst s e1)++substVarTrig :: (Ord v, Ord n, Floating n) => (v -> Maybe (LinExpr v n)) -> ([(v, n)], [(n, LinExpr v n)]) -> Expr v n+substVarTrig s (period, terms) =+  let period2 = linExpr $ substVarLin s (LinExpr 0 period)+      terms2 = map (fmap $ linExpr.substVarLin s) terms+  in foldr (\(p,a) -> (+ (a * sin (makeConstant p + period2))))+     0 terms2++subst :: (Ord n, Ord v, Floating n) => (v -> Maybe (LinExpr v n)) -> Expr v n -> Expr v n+subst s (Expr lt trig nl) =+  linExpr (substVarLin s lt) ++  foldr ((+).substVarTrig s) 0 trig ++  foldr ((+).substVarNonLin s) 0 nl++-- | An empty set of dependencies.+emptyDeps :: Dependencies v n+emptyDeps = Dependencies M.empty M.empty [] M.empty []++simpleSubst :: Eq a => a -> b -> a -> Maybe b+simpleSubst x y z+  | x == z = Just y+  | otherwise = Nothing++-- | Make the expressions on both sides equal, and add the result to the Set of+-- dependencies.+(===) :: (Hashable n, Hashable v, RealFrac (Phase n), Ord v,+          Floating n) => Expr v n -> Expr v n+         -> Dependencies v n+         -> Either (DepError n) (Dependencies v n)+(===) lhs rhs deps = addEq deps (lhs - rhs)++addEqs :: (Hashable v, Hashable n, RealFrac (Phase n), Ord v, Floating n) => Dependencies v n -> [Expr v n] -> Either (DepError n) (Dependencies v n)+addEqs = foldM addEq++addEq :: (Hashable n, Hashable v, RealFrac (Phase n), Ord v,+          Floating n) =>+         Dependencies v n+         -> Expr v n -> Either (DepError n) (Dependencies v n)+addEq deps@(Dependencies _ lin _ _ _) expr =+  addEq0 deps $+  -- substitute known and dependend variables+  subst (flip M.lookup lin) expr+  +-- the following alternative would continue after a redundant error.+-- However a redundant expression is supposed to be an error in metafont.+-- +-- addEqs dep [] = Right dep+-- addEqs dep (e:r) =+--   case addEq0 dep e of+--    Left (InconsistentEq c) ->+--      Left $ InconsistentEq c+--    Left RedundantEq ->+--      addEqs dep r+--    Right newdep ->+--      case addEqs newdep r of+--       Left (InconsistentEq c) ->+--         Left $ InconsistentEq c+--       Left RedundantEq -> Right newdep+--       Right newerdep -> Right newerdep++-- This one is by Cale Gibbard: ++select :: [a] -> [(a, [a])]+select [] = []+select (x:xs) =+  (x,xs) : [(y,x:ys) | (y,ys) <- select xs]+  +addEq0 :: (Hashable v, Hashable n, RealFrac (Phase n), Ord v, Floating n) => Dependencies v n -> Expr v n -> Either (DepError n) (Dependencies v n)+-- adding a constant equation+addEq0 _  (getConstExpr -> Just c) =+  if c == 0 then Left RedundantEq+  else Left (InconsistentEq c)++-- adding a linear equation+addEq0 (Dependencies vdep lin trig trig2 nonlin) (Expr lt [] []) =+  let (v, _, lt2) = splitMax lt+      -- variables that depend on v+      depVars = fromMaybe H.empty (M.lookup v vdep)+      -- substitute v in all dependend variables+      -- and add v as dependend variable+      lin' = M.insert v lt2 $+             H.foldl' (flip $ M.adjust $ substVarLin $+                       simpleSubst v lt2) lin depVars+      -- independend variables substituted+      ltVars = case lt2 of+        LinExpr _ vars -> map fst vars+      -- add dependency link from independend variables to the+      -- substituted equations and v, and remove v (since it has+      -- become dependend, so no variable can depend on it).+      depVars2 = H.insert v depVars+      vdep' = H.foldl'+              (\mp k -> M.insertWith H.union k depVars2 mp)+              (M.delete v vdep) (H.fromList ltVars)+      -- simply substitute var in nonlinear eqs+      nonlin' = map (subst (simpleSubst v lt2)) nonlin+      -- substitute and evaluate trigonometric equations+      trigSubst (p, a, ph, c) =+        subst (simpleSubst v lt2) $+        sin (linExpr $ LinExpr ph p) *+        linExpr a + makeConstant c+      newTrig = map trigSubst trig+      trigSubst2 (v2, ex) =+        subst (simpleSubst v lt2) $+        makeVariable v2 - ex+      newTrig2 =+        map trigSubst2 $+        concatMap M.toList $+        M.elems trig2+      +  in addEqs (Dependencies vdep' lin' [] M.empty []) (newTrig++newTrig2++nonlin')++-- adding a sine equation+addEq0 deps@(Dependencies vdep lin trig trig2 nl)+  (Expr (LinExpr c lt) [(theta, [(alpha, getConst -> Just n)])] []) =+  if null lt then+    -- reduce a sine to linear equation+    addEq0 deps (linExpr $ LinExpr (alpha - asin (-c/n)) theta)+  else+    -- add a variable dependency on the sine equation+    case M.lookup theta trig2 of+     -- no sine with same period+     Nothing -> addSin (LinExpr c lt) alpha n+     Just map2 ->+       case foldr ((+).doSubst)+            (makeConstant c ++             makeConstant n *+             sin (linExpr $ LinExpr alpha theta))+            lt of+        Expr lt2 [(_, [(alpha2, getConst -> Just n2)])] []+          | isNothing(getConst lt2)+          -> addSin lt2 alpha2 n2+        e2 -> addEq0 deps e2+       where+         doSubst (v,c2) = case M.lookup v map2 of+           Nothing -> makeVariable v * makeConstant c2+           Just e2 -> e2 * makeConstant c2+  where+    addSin l' a' n' =+      let (v, c', r) = splitMax l'+          trig2' = M.insertWith M.union theta+                   (M.singleton v $+                    Expr r [(theta, [(a', LinExpr (n'/negate c') [])])] [])+                   trig2+      in Right $ Dependencies  vdep lin trig trig2' nl++--  adding the first sine equation+addEq0 (Dependencies d lin [] trig2 nl) (Expr (getConst -> Just c) [(theta, [(alpha, n)])] []) =+  Right $ Dependencies d lin [(theta, n, alpha, c)] trig2 nl++-- try reducing this equation with another sine equation+addEq0 (Dependencies deps lin trig trig2 nl)+  (Expr (getConst -> Just x) [(theta, [(a, n)])] []) =+  case mapMaybe similarTrig $ select trig of+   -- no matching equation found+   [] -> Right $ Dependencies deps lin ((theta, n, a, x):trig) trig2 nl+   -- solve for angle and amplitude, and add resulting linear+   -- equations+   l -> addEqs (Dependencies deps lin rest trig2 nl) [lin1, lin2]+     where+       ((b,y), rest) = maximumBy (maxTrig `on` fst) l+       maxTrig (t1,_) (t2,_) = +         compare ((t1-a)`dmod`pi) ((t2-a)`dmod`pi)+       d      = sin(a-b)+       e      = y*cos(a-b)-x+       theta2 = atan (-y*d/e)-b ++                (if (d*e) < 0 then pi else 0)+       n2     = sqrt(y*y + e*e/(d*d))+       lin1   = linExpr $ LinExpr (-theta2) theta+       lin2   = linExpr n - makeConstant n2+  where+    similarTrig ((t,m,b,y),rest)+      | Just r <- termIsMultiple m n,+        t == theta &&+        (b-a) `dmod` pi > pi/8 =+          Just ((b,y/r),rest)+      | otherwise = Nothing++-- just add any other equation to the list of nonlinear equations+addEq0 (Dependencies d lin trig trig2 nonlin) e =+  Right $ Dependencies d lin trig trig2 (e:nonlin)++dmod :: RealFrac a => a -> a -> a+dmod a b = abs((a/b) - fromInteger (round (a/b)) * b)++-- put the variable with the maximum coefficient on the lhs of the+-- equation+splitMax :: (Ord b, Fractional b, Eq v) => LinExpr v b -> (v, b, LinExpr v b)+splitMax (LinExpr c t) =+  let (v,c2) = maximumBy (compare `on` (abs.snd)) t+  in (v, c2,+      LinExpr (-c/c2) $+      map (fmap (negate.(/c2))) $+      filter ((/= v).fst) t)+      +-- | Return True if the variable is known or dependend.+varDefined :: (Eq v, Hashable v) => Dependencies v n -> v -> Bool+varDefined (Dependencies _ dep _ _ _) v =+  case M.lookup v dep of+    Nothing -> False+    _ -> True++-- | Return all dependend variables with their dependencies.+dependendVars :: (Eq n) => Dependencies v n -> [(v, LinExpr v n)]+dependendVars (Dependencies _ lin _ _ _) =+  filter (notConst.snd) (M.toList lin)+  where+    notConst (LinExpr _ []) = False+    notConst _ = True+  ++-- | Return all known variables.+knownVars :: Dependencies v n -> [(v, n)]+knownVars (Dependencies _ lin _ _ _) =+  mapMaybe knownVar $ M.toList lin+  where+    knownVar (v, LinExpr n []) = Just (v, n)+    knownVar _ = Nothing++-- -- | Return all independend variables.+-- freeVars :: (Eq v, Hashable v) => Dependencies n v -> [v]+-- freeVars (Dependencies dep) =+--   HS.toList $ M.foldl' addVars HS.empty dep+--   where addVars s (LinExpr _ a) =+--           HS.union s $ HS.fromList $ map fst a++-- | Return the value of the variable, or a list of variables+-- it depends on.  Only linear dependencies are shown.+getKnown :: (Eq v, Hashable v) => Dependencies v n -> v -> Either [v] n+getKnown (Dependencies _ lin _ _ _) var =+  case M.lookup var lin of+    Nothing -> Left  []+    Just (LinExpr a []) ->+      Right a+    Just (LinExpr _ v) ->+      Left $ map fst v++trigToExpr :: (Ord n, Ord v, Floating n) => TrigEq v n -> Expr v n+trigToExpr (p, a, ph, c) =+  linExpr a * sin(linExpr $ LinExpr ph p) ++  makeConstant c++-- | Give all nonlinear equations as an `Expr` equal to 0.+nonlinearEqs :: (Ord n, Ord v, Floating n) => Dependencies v n -> [Expr v n]+nonlinearEqs  (Dependencies _ _ trig trig2 nl) =+  map trigToExpr trig +++  map (\(v, e) -> makeVariable v - e) +  (concatMap M.toList (M.elems trig2)) +++  nl+  +-- | Make the pairs of expressions on both sides equal, and add the+-- result to the Set of dependencies.  No error is signaled if the+-- equation for one of the sides is redundant for example in (x, 0) ==+-- (y, 0).+(=&=) :: (Hashable n, Hashable v, RealFrac (Phase n), Ord v, Floating n) => (Expr v n, Expr v n) -> (Expr v n, Expr v n) -> Dependencies v n -> Either (DepError n) (Dependencies v n)+(=&=) (a, b) (c, d) dep = +  case (a === c) dep of+    Left RedundantEq ->+      (b === d) dep+    Right res ->+      case (b === d) res of+        Left RedundantEq -> Right res+        Right res2 -> Right res2+        err -> err+    err -> err++-- | Solve a list of equations in order.  Returns either a new set of+-- dependencies, or signals an error.+solveEqs :: Dependencies v n -> [Dependencies v n -> Either (DepError n) (Dependencies v n)] -> Either (DepError n) (Dependencies v n)+solveEqs = foldM $ flip ($)++-- | Show all variables and equations.+showVars :: (Show n, Show v, Show a, Ord n, Ord v, Floating n) => Either (DepError a) (Dependencies v n) -> IO ()+showVars (Left e) = print e+showVars (Right dep) = print dep
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ mfsolve.cabal view
@@ -0,0 +1,39 @@+Name:		mfsolve+Version: 	0.1.0+Synopsis:	Equation solver and calculator à la metafont+Category: 	Math+Copyright: 	Kristof Bastiaensen (2015)+Stability:	Unstable+License:	BSD3+License-file:	LICENSE+Author:		Kristof Bastiaensen+Maintainer:	Kristof Bastiaensen+Bug-Reports: 	https://github.com/kuribas/mfsolve/issues+Build-type:	Simple+Cabal-version:	>=1.8+Description:  An equation solver and calculator in the spirit of Metafont.+  .+  Like metafont, it can solve linear equations, and evaluate nonlinear expressions.  In addition to metafont, it also solves for angles, and makes the solution independend of the order of the equations.+ +source-repository head+  type:		git+  location:	https://github.com/kuribas/mfsolve++Library+  Ghc-options: -Wall+  Build-depends: base >= 3 && < 5, unordered-containers > 0.2, hashable >= 0.1.2+  Exposed-Modules:+    Math.MFSolve++test-suite test+  type: exitcode-stdio-1.0+  hs-source-dirs:+    tests+  main-is:+    test.hs+  build-depends:+    base >= 4 && < 5,+    tasty >= 0.8,+    tasty-hunit >= 0.9,+    mfsolve+
+ tests/test.hs view
@@ -0,0 +1,182 @@+import Test.Tasty+import Test.Tasty.HUnit+import Math.MFSolve+import Data.List+import Data.Maybe++data TestExpr = TestExpr+                ((SimpleVar -> Double) -> Double)+                (Expr SimpleVar Double)++type TestEq = ((SimpleVar -> Double) -> String,+               (Expr SimpleVar Double,+                Expr SimpleVar Double))++testVar :: String -> TestExpr+testVar str = TestExpr ($ SimpleVar str) (makeVariable $ SimpleVar str)++testConst :: Double -> TestExpr+testConst n = TestExpr (const n) (makeConstant n)++testBin :: (Double -> Double -> Double)+        -> (Expr SimpleVar Double+            -> Expr SimpleVar Double -> Expr SimpleVar Double)+        -> TestExpr -> TestExpr -> TestExpr+testBin f f2(TestExpr a b) (TestExpr c d) =+  TestExpr (\s -> f (a s) (c s)) (f2 b d)++testUn :: (Double -> Double)+       -> (Expr SimpleVar Double -> Expr SimpleVar Double)+       -> TestExpr -> TestExpr+testUn f f2(TestExpr a b) =+  TestExpr (f.a) (f2 b)++testEq :: Double -> Double -> String -> String -> String+testEq a b as bs+  | abs(a-b) <= 1e-10 =+      ""+  | otherwise =+      "substituting solutions in:\n"+      ++ as ++ " = " ++ bs ++ "\n" +++      "gives: " ++ show a ++ " = " ++ show b ++ "\n"++infixr 1 ?= +(?=) :: TestExpr -> TestExpr -> TestEq++TestExpr a b ?= TestExpr c d =+  (\s -> testEq (a s) (c s) (show b) (show d),+   (b,d))++instance Num TestExpr where+  (+) = testBin (+) (+)+  (*) = testBin (*) (*)+  signum = testUn signum signum+  abs = testUn abs abs+  fromInteger = testConst . fromInteger+  negate = testUn negate negate++instance Fractional TestExpr where+  recip = testUn recip recip+  fromRational = testConst . fromRational++instance Floating TestExpr where+  pi = testConst pi+  exp = testUn exp exp+  log = testUn log log+  sin = testUn sin sin +  cos = testUn cos cos+  cosh = testUn cosh cosh+  atanh = testUn atanh atanh+  tan = testUn tan tan+  sinh = testUn sinh sinh+  asin = testUn asin asin+  acos = testUn acos acos+  asinh = testUn asinh asinh+  acosh = testUn acosh acosh+  atan = testUn atan atan++a, b, c, d, x, y :: TestExpr+[a, b, c, d, x, y] = map testVar ["a", "b", "c", "d", "x", "y"]++randomSol :: [(SimpleVar, Double)]+randomSol = [+  (SimpleVar "a", 0.897),+  (SimpleVar "b", 0.905),+  (SimpleVar "c", -0.585),+  (SimpleVar "d", 0.018),+  (SimpleVar "x", -0.3628),+  (SimpleVar "y", 0.887)]++trySolve :: [TestEq]  -> Assertion+trySolve eqs =+  assertString $ fromMaybe "" $+  find (not.null) $ map solveOne $+  permutations eqs++solveOne :: [TestEq] -> String+solveOne eqs =  +  case solveEqs emptyDeps $+       map (uncurry (===) . snd) eqs+  of+   Left RedundantEq ->+     "Found redundant equation"+   Left (InconsistentEq n) ->+     "Equation off by "++ show n+   Right d+     | not $ null $ nonlinearEqs d ->+       "Some nonlinear equations were unevaluated:\n" +++       show d+     | not $ null $ dependendVars d ->+       "Some linear dependencies were left:\n" +++       show d+     | otherwise ->+       let kv = knownVars d ++ randomSol+           sols = map (($ (fromMaybe 0 . (`lookup` kv))).fst) eqs+       in case find (not.null) sols of+           Nothing -> ""+           Just er ->+             "Solution didn't match equations:\n" +++             show d ++ er+         +tests :: TestTree+tests = testGroup "Tests" [unitTests]++unitTests = testGroup "Unit tests" [+  testCase "Linear system of equations" $+  trySolve [+    3*a + 2*b - c ?= 1,+    2*a - 2*b + 4 ?= -2,+    -a + b/2 - c ?= 0],++  testCase "Solve single sine" $+  trySolve [+    2*sin(0.3*a - 2*b + pi/8) ?= -1.42,+    3*a+2*b ?= 7],++  testCase "Adding sinewaves with same period" $+  trySolve [+    a ?= 2*b,+    2*sin(a + 0.1) + 3*cos(2*b + 0.1) ?= -0.524],+    +  testCase "Adding sinewaves with same period (using substitution)" $+  trySolve [+    2*x + 3*y ?= -0.524,+    a ?= 2*b,+    x ?= sin(a + 0.1),+    y ?= cos(2*b + 0.1)],++  testCase "Solving for angle and amplitude" $+  trySolve [+    a*sin (b+c+0.1) ?= 5,+    a*cos (b+c+0.2) ?= 6,+    b+2*c ?= 3],++  testCase "Mixing nonlinear and trigonometric functions" $+  trySolve [+    sqrt a * (b + sin(2*c)) ?= 0.243,+    a - b ?= 1,+    a + b ?= 3],+  +  testCase "Simplifying nonlinear into trigonometric." $+  trySolve [+    2*sin(2*c*sqrt b) ?= a,+    a - b ?= 1,+    a + b ?= 3],++  testCase "Simplifying trigonometric into linear" $+  trySolve [+    sin(a+b*2) ?= c,+    2*a+4*b ?= pi/8,+    b ?= 2*c],++  testCase "Rotation by rotation matrix." $+  trySolve [+    cos a*x - sin a*y ?= 10,+    sin a*x + cos a*y ?= 13,+    x ?= d*sin b,+    y ?= d*cos b,+    b ?= pi/3+    ]+  ]++main = defaultMain tests