limp 0.3.2.0 → 0.3.2.1
raw patch · 16 files changed
+1507/−82 lines, 16 files
Files
- limp.cabal +17/−2
- src/Numeric/Limp/Rep.hs +6/−80
- src/Numeric/Limp/Rep/Arbitrary.hs +22/−0
- src/Numeric/Limp/Rep/IntDouble.hs +29/−0
- src/Numeric/Limp/Rep/Rep.hs +64/−0
- src/Numeric/Limp/Solve/Branch/Simple.hs +84/−0
- src/Numeric/Limp/Solve/Simplex/Maps.hs +316/−0
- src/Numeric/Limp/Solve/Simplex/StandardForm.hs +227/−0
- tests/Arbitrary/Assignment.hs +33/−0
- tests/Arbitrary/Program.hs +82/−0
- tests/Arbitrary/Var.hs +38/−0
- tests/BranchExample.hs +94/−0
- tests/Convert.hs +19/−0
- tests/SimplexExample.hs +55/−0
- tests/Simplexs.hs +311/−0
- tests/Simplify.hs +110/−0
limp.cabal view
@@ -1,5 +1,5 @@ name: limp-version: 0.3.2.0+version: 0.3.2.1 synopsis: representation of Integer Linear Programs description: so far, this package just provides two representations for linear programs: "Numeric.Limp.Program", which is what I expect end-users to use, and "Numeric.Limp.Canon", which is simpler, but would be less nice for writing linear programs.@@ -23,6 +23,9 @@ hs-source-dirs: src exposed-modules: Numeric.Limp.Rep+ Numeric.Limp.Rep.Rep+ Numeric.Limp.Rep.Arbitrary+ Numeric.Limp.Rep.IntDouble Numeric.Limp.Error Numeric.Limp.Program.Bounds@@ -45,7 +48,10 @@ Numeric.Limp.Canon.Simplify.Subst Numeric.Limp.Canon.Simplify - -- other-modules: + Numeric.Limp.Solve.Simplex.StandardForm+ Numeric.Limp.Solve.Simplex.Maps+ Numeric.Limp.Solve.Branch.Simple+ build-depends: base < 5, containers == 0.5.*@@ -59,6 +65,15 @@ type: exitcode-stdio-1.0 main-is: Main.hs hs-source-dirs: tests+ other-modules:+ Arbitrary.Assignment+ Arbitrary.Program+ Arbitrary.Var+ BranchExample+ Convert+ SimplexExample+ Simplexs+ Simplify build-depends: base < 5, containers == 0.5.*,
src/Numeric/Limp/Rep.hs view
@@ -5,85 +5,11 @@ -- We bundle Z and R up into a single representation instead of abstracting over both, -- because we must be able to convert from Z to R without loss. ---module Numeric.Limp.Rep where--import Data.Map (Map)-import qualified Data.Map as M--import Data.Monoid---- | The Representation class. Requires its members @Z c@ and @R c@ to be @Num@, @Ord@ and @Eq@.------ For some reason, for type inference to work, the members must be @data@ instead of @type@.--- This gives some minor annoyances when unpacking them. See 'unwrapR' below.----class ( Num (Z c), Ord (Z c), Eq (Z c), Integral (Z c)- , Num (R c), Ord (R c), Eq (R c), RealFrac (R c)) => Rep c where-- -- | Integers- data Z c- -- | Real numbers- data R c-- -- | Convert an integer to a real. This should not lose any precision.- -- (whereas @fromIntegral 1000 :: Word8@ would lose precision)- fromZ :: Z c -> R c- fromZ = fromIntegral----- | An assignment from variables to values.--- Maps integer variables to integers, and real variables to reals.-data Assignment z r c- = Assignment (Map z (Z c)) (Map r (R c))--deriving instance (Show (Z c), Show (R c), Show z, Show r) => Show (Assignment z r c)--instance (Ord z, Ord r) => Monoid (Assignment z r c) where- mempty = Assignment M.empty M.empty- mappend (Assignment z1 r1) (Assignment z2 r2)- = Assignment (M.union z1 z2) (M.union r1 r2)----- | Retrieve value of integer variable - or 0, if there is no value.-zOf :: (Rep c, Ord z) => Assignment z r c -> z -> Z c-zOf (Assignment zs _) z- = maybe 0 id $ M.lookup z zs---- | Retrieve value of real variable - or 0, if there is no value.-rOf :: (Rep c, Ord r) => Assignment z r c -> r -> R c-rOf (Assignment _ rs) r- = maybe 0 id $ M.lookup r rs---- | Retrieve value of an integer or real variable, with result cast to a real regardless.-zrOf :: (Rep c, Ord z, Ord r) => Assignment z r c -> Either z r -> R c-zrOf a = either (fromZ . zOf a) (rOf a)--assSize :: Assignment z r c -> Int-assSize (Assignment mz mr)- = M.size mz + M.size mr----- | A representation that uses native 64-bit ints and 64-bit doubles.--- Really, this should be 32-bit ints.-data IntDouble--instance Rep IntDouble where- -- | Automatically defer numeric operations to the native int.- newtype Z IntDouble = Z Int- deriving (Ord,Eq,Integral,Real,Num,Enum)- newtype R IntDouble = R Double- deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)---- | Define show manually, so we can strip out the "Z" and "R" prefixes.-instance Show (Z IntDouble) where- show (Z i) = show i--instance Show (R IntDouble) where- show (R i) = show i--+module Numeric.Limp.Rep+ ( module Numeric.Limp.Rep.Rep+ , module Numeric.Limp.Rep.IntDouble )+ where --- | Convert a wrapped (R IntDouble) to an actual Double.-unwrapR :: R IntDouble -> Double-unwrapR (R d) = d+import Numeric.Limp.Rep.Rep+import Numeric.Limp.Rep.IntDouble
+ src/Numeric/Limp/Rep/Arbitrary.hs view
@@ -0,0 +1,22 @@+-- | Arbitrary precision number representation+module Numeric.Limp.Rep.Arbitrary where+import Numeric.Limp.Rep.Rep++-- | A representation that uses arbitrary-sized Integers and Rationals+data Arbitrary++instance Rep Arbitrary where+ -- | Automatically defer numeric operations to the native int.+ newtype Z Arbitrary = Z Integer+ deriving (Ord,Eq,Integral,Real,Num,Enum)+ newtype R Arbitrary = R Rational+ deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)++-- | Define show manually, so we can strip out the "Z" and "R" prefixes.+instance Show (Z Arbitrary) where+ show (Z i) = show i++instance Show (R Arbitrary) where+ show (R i) = show i++
+ src/Numeric/Limp/Rep/IntDouble.hs view
@@ -0,0 +1,29 @@+-- | Fixed/floating precision number representation+module Numeric.Limp.Rep.IntDouble where+import Numeric.Limp.Rep.Rep++-- | A representation that uses native 64-bit ints and 64-bit doubles.+-- Really, this should be 32-bit ints.+data IntDouble++instance Rep IntDouble where+ -- | Automatically defer numeric operations to the native int.+ newtype Z IntDouble = Z Int+ deriving (Ord,Eq,Integral,Real,Num,Enum)+ newtype R IntDouble = R Double+ deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)++-- | Define show manually, so we can strip out the "Z" and "R" prefixes.+instance Show (Z IntDouble) where+ show (Z i) = show i++instance Show (R IntDouble) where+ show (R i) = show i++++-- | Convert a wrapped (R IntDouble) to an actual Double.+unwrapR :: R IntDouble -> Double+unwrapR (R d) = d++
+ src/Numeric/Limp/Rep/Rep.hs view
@@ -0,0 +1,64 @@+-- | Representation of integers (Z) and reals (R) of similar precision.+-- Programs are abstracted over this, so that ideally in the future we could have a+-- solver that produces Integers and Rationals, instead of just Ints and Doubles.+--+-- We bundle Z and R up into a single representation instead of abstracting over both,+-- because we must be able to convert from Z to R without loss.+--+module Numeric.Limp.Rep.Rep where++import Data.Map (Map)+import qualified Data.Map as M++import Data.Monoid++-- | The Representation class. Requires its members @Z c@ and @R c@ to be @Num@, @Ord@ and @Eq@.+--+-- For some reason, for type inference to work, the members must be @data@ instead of @type@.+-- This gives some minor annoyances when unpacking them. See 'unwrapR' below.+--+class ( Num (Z c), Ord (Z c), Eq (Z c), Integral (Z c)+ , Num (R c), Ord (R c), Eq (R c), RealFrac (R c)) => Rep c where++ -- | Integers+ data Z c+ -- | Real numbers+ data R c++ -- | Convert an integer to a real. This should not lose any precision.+ -- (whereas @fromIntegral 1000 :: Word8@ would lose precision)+ fromZ :: Z c -> R c+ fromZ = fromIntegral+++-- | An assignment from variables to values.+-- Maps integer variables to integers, and real variables to reals.+data Assignment z r c+ = Assignment (Map z (Z c)) (Map r (R c))++deriving instance (Show (Z c), Show (R c), Show z, Show r) => Show (Assignment z r c)++instance (Ord z, Ord r) => Monoid (Assignment z r c) where+ mempty = Assignment M.empty M.empty+ mappend (Assignment z1 r1) (Assignment z2 r2)+ = Assignment (M.union z1 z2) (M.union r1 r2)+++-- | Retrieve value of integer variable - or 0, if there is no value.+zOf :: (Rep c, Ord z) => Assignment z r c -> z -> Z c+zOf (Assignment zs _) z+ = maybe 0 id $ M.lookup z zs++-- | Retrieve value of real variable - or 0, if there is no value.+rOf :: (Rep c, Ord r) => Assignment z r c -> r -> R c+rOf (Assignment _ rs) r+ = maybe 0 id $ M.lookup r rs++-- | Retrieve value of an integer or real variable, with result cast to a real regardless.+zrOf :: (Rep c, Ord z, Ord r) => Assignment z r c -> Either z r -> R c+zrOf a = either (fromZ . zOf a) (rOf a)++assSize :: Assignment z r c -> Int+assSize (Assignment mz mr)+ = M.size mz + M.size mr+
+ src/Numeric/Limp/Solve/Branch/Simple.hs view
@@ -0,0 +1,84 @@+-- | The simplest, stupidest possible branch and bound algorithm.+--+--+module Numeric.Limp.Solve.Branch.Simple+ (branch, makeIntegral)+ where+import Numeric.Limp.Canon.Program+import Numeric.Limp.Canon.Simplify+import Numeric.Limp.Rep++import Control.Applicative+import Control.Monad+import qualified Data.Map as M+import Data.Monoid++branch+ :: (Ord z, Ord r, Rep c)+ => (Program z r c -> Maybe (Assignment () (Either z r) c, R c))+ -> Program z r c+ -> Maybe (Assignment z r c, R c)+branch solver start_prog+ = go mempty start_prog+ where+ go ass p+ -- TODO:+ -- simp can actually change the objective function+ -- because Canon doesn't store a constant summand on the objective.+ -- we really need to return the modified summand and take that into account when+ -- choosing between two integer assignments.+ | Right (ass', p') <- simplify' ass p+ = do (assRelax,co) <- solver p'+ case makeIntegral assRelax of+ Left (var, val)+ -> branchon p' ass' (Left var) val+ Right r+ -> Just (ass' <> r, co)+ | otherwise+ = Nothing++ branchon p ass var val+ = let lo = addBound p var (Just (fromZ $ truncate val + 1), Nothing)+ up = addBound p var (Nothing, Just (fromZ $ truncate val))+ loB = go ass lo+ upB = go ass up+ in case (loB, upB) of+ (Just (a1, o1), Just (a2, o2))+ | o1 > o2+ -> Just (a1, o1)+ | otherwise+ -> Just (a2, o2)+ (Just r, Nothing)+ -> Just r+ (Nothing, Just r)+ -> Just r+ (Nothing, Nothing)+ -> Nothing+ ++ addBound p v b+ = let bs = _bounds p+ b' = maybe (Nothing,Nothing) id+ $ M.lookup v bs+ in p { _bounds = M.insert v (mergeBounds b b') bs }++makeIntegral+ :: (Ord z, Ord r, Rep c)+ => Assignment () (Either z r) c+ -> Either (z, R c)+ (Assignment z r c)+makeIntegral (Assignment _ vs)+ = uncurry Assignment+ <$> foldM go (M.empty, M.empty) (M.toList vs)+ where+ go (zs,rs) (var, val)+ = case var of+ Right r+ -> return (zs, M.insert r val rs)+ Left z+ | val' <- truncate val+ , val == fromZ val'+ -> return (M.insert z val' zs, rs)+ | otherwise+ -> Left (z, val)+
+ src/Numeric/Limp/Solve/Simplex/Maps.hs view
@@ -0,0 +1,316 @@+-- | The simplest, stupidest possible simplex algorithm.+-- The idea here is to be slow, but "obviously correct" so other algorithms+-- can be verified against it.+--+-- That's the plan, at least. For now this is just a first cut of trying to implement simplex.+--+module Numeric.Limp.Solve.Simplex.Maps+ where+import Numeric.Limp.Rep++import Numeric.Limp.Solve.Simplex.StandardForm++import Control.Arrow+import qualified Data.Map as M+import Data.Function (on)+import Data.List (minimumBy, sortBy)+++-- | Result of a single pivot attempt+data IterateResult z r c+ -- | Maximum reached!+ = Done+ -- | Pivot was made+ | Progress (Standard z r c)+ -- | No progress can be made: unbounded along the objective+ | Stuck++deriving instance (Show z, Show r, Show (R c)) => Show (IterateResult z r c)+++-- | Try to find a pivot and then perform it.+-- We're assuming, at this stage, that the existing solution is feasible.+simplex1 :: (Ord z, Ord r, Rep c)+ => Standard z r c -> IterateResult z r c+simplex1 s+ -- Check if there are any positive columns in the objective:+ = case pivotCols of+ -- if there are none, we are already at the maximum+ []+ -> Done+ -- there are some; try to find the first pivot row that works+ _+ -> go pivotCols+ where++ -- Check if there's any row worth pivoting on for this column.+ -- We're trying to see if we can increase the value of this+ -- column's variable from zero.+ go ((pc,_):pcs)+ = case pivotRowForCol s pc of+ Nothing -> go pcs+ Just pr+ -> Progress+ -- Perform the pivot.+ -- This moves the variable pr out of the basis, and pc into the basis.+ $ pivot s (pr,pc)++ -- We've tried all the pivot columns and failed.+ -- This means there's no edge we can take to increase our objective,+ -- so it must be unbounded.+ go []+ = Stuck+++ -- We want to find some positive column from the objective.+ -- In fact, find all of them and order descending.+ pivotCols+ = let ls = M.toList $ fst $ _objective s+ kvs = sortBy (compare `on` (negate . snd)) ls+ in filter ((>0) . snd) kvs+++-- | Find pivot row for given column.+-- We're trying to find a way to increase the value of+-- column from zero, and the returned row will be decreased to zero.+-- Since all variables are >= 0, we cannot return a row that would set the column to negative.+pivotRowForCol :: (Ord z, Ord r, Rep c)+ => Standard z r c+ -> StandardVar z r+ -> Maybe (StandardVar z r)+pivotRowForCol s col+ = fmap fst+ $ minBy' (compare `on` snd)+ $ concatMap (\(n,r)+ -> let rv = lookupRow r col+ o = objOfRow r+ in if rv > 0+ then [(n, o / rv)]+ else [])+ $ M.toList+ $ _constraints s++-- | Find minimum, or nothing if empty+minBy' :: (a -> a -> Ordering) -> [a] -> Maybe a+minBy' _ []+ = Nothing+minBy' f ls+ = Just $ minimumBy f ls+++-- | Perform pivot for given row and column.+-- We normalise row so that row.column = 1+--+-- > norm = row / row[column]+--+-- Then, for all other rows including the objective,+-- we want to make sure its column entry is zero:+--+-- > row' = row - row[column]*norm+--+-- In the end, this means "column" will be an identity column, or a basis column.+--+pivot :: (Ord z, Ord r, Rep c)+ => Standard z r c+ -> (StandardVar z r, StandardVar z r)+ -> Standard z r c+pivot s (pr,pc)+ = let norm = normaliseRow+ -- All other rows+ rest = filter ((/=pr) . fst) $ M.toList $ _constraints s+ in Standard+ { _constraints = M.fromList ((pc, norm) : map (id *** fixup norm) rest)+ , _objective = fixup norm $ _objective s+ , _substs = _substs s }+ where+ -- norm = row / row[column]+ normaliseRow+ | Just row@(rm, ro) <- M.lookup pr $ _constraints s+ = let c' = lookupRow row pc+ in (M.map (/c') rm, ro / c')++ -- Pivot would not be chosen if row doesn't exist..+ | otherwise+ = (M.empty, 0)++ -- row' = row - row[column]*norm+ fixup (nm,no) row@(rm,ro)+ = let co = lookupRow row pc+ in {- row' = row - co*norm -}+ ( M.unionWith (+) rm (M.map ((-co)*) nm)+ , ro - co * no )+++-- | Single phase of simplex.+-- Keep repeating until no progress can be made.+single_simplex :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Maybe (Standard z r c)+single_simplex s+ = case simplex1 s of+ Done -> Just s+ Progress s' -> single_simplex s'+ Stuck -> Nothing+++-- | Two phase:+-- first, find a satisfying solution.+-- then, solve simplex as normal.+simplex+ :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Maybe (Standard z r c)+simplex s+ = find_initial_sat s+ >>= single_simplex++-- | Find a satisfying solution.+-- if there are any rows with negative values, this means their basic values are negative+-- (which is not satisfying the x >= 0 constraint)+-- these negative-valued rows must be pivoted around using modified pivot criteria+find_initial_sat+ :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Maybe (Standard z r c)+find_initial_sat s+ = case negative_val_rows of+ [] -> Just s+ rs -> go rs+ where+ -- Find all rows with negative values+ -- because their current value is not feasible+ negative_val_rows+ = filter ((<0) . objOfRow . snd)+ $ M.toList+ $ _constraints s++ -- Find largest negative (closest to zero) to pivot on:+ -- pivoting on a negative will negate the value, setting it to positive+ min_of_row (_,(rm,_))+ = minBy' (compare `on` (negate . snd))+ $ filter ((<0) . snd)+ $ M.toList rm+++ -- There is no feasible solution+ go []+ = Nothing++ -- Try pivoting on the rows + go (r:rs)+ | Just (pc,_) <- min_of_row r+ , Just pr <- pivotRowForNegatives pc+ = simplex+ $ pivot s (pr, pc)++ | otherwise+ = go rs++ -- opposite of pivotRowForCol...+ pivotRowForNegatives col+ = fmap fst+ $ minBy' (compare `on` (negate . snd))+ $ concatMap (\(n,r)+ -> let rv = lookupRow r col+ o = objOfRow r+ in if rv < 0+ then [(n, o / rv)]+ else [])+ $ M.toList+ $ _constraints s+++ ++-- Get map of each constraint's value+assignmentAll :: (Rep c)+ => Standard z r c+ -> (M.Map (StandardVar z r) (R c), R c)+assignmentAll s+ = ( M.map val (_constraints s)+ , objOfRow (_objective s))+ where+ val (_, v)+ = v++-- Perform reverse substitution on constraint values+-- to get original values (see StandardForm)+assignment+ :: (Ord z, Ord r, Rep c)+ => Standard z r c+ -> (Assignment () (Either z r) c, R c)+assignment s+ = ( Assignment M.empty $ M.union vs' rs'+ , o )+ where+ (vs, o) = assignmentAll s++ vs' = M.fromList+ $ concatMap only_svs+ $ M.toList vs++ rs' = M.map eval $ _substs s++ eval (lin,co)+ = M.fold (+) co+ $ M.mapWithKey (\k r -> r * (maybe 0 id $ M.lookup k vs))+ $ lin++ only_svs (SV v, val)+ = [(v, val)]+ only_svs _+ = []++++-- Junk ---------------++-- | Minimise whatever variables are 'basic' in given standard+-- input must not already have an objective row "SVO",+-- because the existing objective is added as a new row with that name+minimise_basics+ :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Standard z r c+minimise_basics s+ = s+ { _objective = (M.map (const (1)) $ _constraints s, 0)+ , _constraints = M.insert SVO (_objective s) (_constraints s)+ }++-- | Find the basic variables and "price them out" of the objective function,+-- by subtracting multiples of the basic row from objective+pricing_out + :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Standard z r c+pricing_out s+ = s+ { _objective = M.foldWithKey go+ (_objective s)+ (_constraints s)+ }+ where+ go v row@(rm,ro) obj@(om,oo)+ | coeff <- lookupRow obj v+ , coeff /= 0+ , rowv <- lookupRow row v+ , mul <- -(coeff / rowv)+ = -- rowv = 1+ -- obj' = obj - (coeff/rowv)*row+ ( M.unionWith (+) om (M.map (mul*) rm)+ , oo + mul*ro )+ | otherwise+ = obj++-- | Pull the previously-hidden objective out of constraints, and use it+drop_fake_objective+ :: (Ord z, Ord r, Rep c)+ => Standard z r c -> Standard z r c+drop_fake_objective s+ | cs <- _constraints s+ , Just o <- M.lookup SVO cs+ = s+ { _objective = o+ , _constraints = M.delete SVO cs }++ | otherwise+ = s+++
+ src/Numeric/Limp/Solve/Simplex/StandardForm.hs view
@@ -0,0 +1,227 @@+-- | Standard form for programs: only equalities and all variables >= 0+-- To convert an arbitrary program to this form, we need to:+--+-- Convert unconstrained (-inf <= x <= +inf) variable into two separate parts, x+ and x-+-- wherever x occurs, it will be replaced with "x+" - "x-".+--+-- Convert variables with non-zero lower bounds (c <= x) to a new variable x', so that+-- x = x' + c+--+-- The opposite of these conversions must be performed when extracting a variable assignment+-- from the solved program.+--+-- All constraints are converted into a less-than with a constant on the right, and then+-- these less-than constraints (f <= c) have a slack variable s added such that+-- f + s == c && s >= 0+--+module Numeric.Limp.Solve.Simplex.StandardForm+ where+import Numeric.Limp.Rep+import Numeric.Limp.Canon.Constraint+import Numeric.Limp.Canon.Linear+import qualified Numeric.Limp.Canon.Program as C++import qualified Data.Map as M+import qualified Data.Set as S+++-- | A single linear function with a constant summand+type StandardRow z r c+ = (StandardLinear z r c, R c)++-- | Entire program in standard form, as well as substitutions required to extract an assignment+data Standard z r c+ = Standard+ { _objective :: StandardRow z r c+ , _constraints :: M.Map (StandardVar z r) (StandardRow z r c)+ , _substs :: StandardSubst z r c+ }+deriving instance (Show z, Show r, Show (R c)) => Show (Standard z r c)++type StandardSubst z r c+ = M.Map (Either z r) (StandardRow z r c)++type StandardLinear z r c+ = M.Map (StandardVar z r) (R c)++data StandardVar z r+ -- | A normal variable+ = SV (Either z r)++ -- | A slack variable, introduced to make less-eq constraints into equalities+ | SVS Int+ -- | Magic objective, used when hiding an existing objective as a constraint+ -- and creating a new objective+ | SVO ++ -- | When a variable has a lower bound other than 0, we replace all occurences with+ -- with a new version minus the lower bound.+ -- x >= 5+ -- ==>+ -- Lx - 5 >= 5+ -- ==>+ -- Lx >= 0+ | SVLower (Either z r)++ -- | When unconstrained variables are encountered, they are replaced with+ -- x = SVPos x - SVNeg x+ -- so both parts can be constrained to >= 0.+ | SVPos (Either z r)+ | SVNeg (Either z r)+ deriving (Eq, Ord, Show)+++-- | Sum a list of linear functions together+addLinears+ :: (Ord z, Ord r, Rep c)+ => [(StandardLinear z r c, R c)] -> (StandardLinear z r c, R c)+addLinears []+ = (M.empty, 0)+addLinears ((lin,co):rs)+ = let (lin',co') = addLinears rs+ in (M.unionWith (+) lin lin', co + co')+++-- | Perform substitution over a linear function/row+substLinear+ :: (Ord z, Ord r, Rep c)+ => StandardSubst z r c -> (StandardLinear z r c, R c) -> (StandardLinear z r c, R c)+substLinear sub (lin, co)+ = let (lin', co') = addLinears + $ map subby + $ M.toList lin+ in (lin', co + co')+ where+ subby (var, coeff)+ = case var of+ SV s+ | Just (vs,cnst) <- M.lookup s sub+ -> (M.map (*coeff) vs, -cnst * coeff)+ _+ -> (M.fromList [(var, coeff)], 0)+++-- | Convert canon program into standard form+standard :: (Ord z, Ord r, Rep c)+ => C.Program z r c+ -> Standard z r c+standard p+ = Standard+ { _objective = objective+ , _constraints = constraints+ , _substs = substs }+ where+ fv = C.varsOfProgram p+ bs = C._bounds p++ -- Objective is just negated+ objective+ = substLinear substs+ ( M.map negate+ $ standardOfLinear $ C._objective p+ , 0)++ -- Constraints are created for original program's bounds and constraints+ -- and substitution is performed.+ -- Each constraint/row receives its own slack variable.+ constraints+ = M.fromList+ $ zipWith (\c s -> (s, substLinear substs $ c s))+ ( constrs ++ bounds )+ ( map SVS [1..] )++ -- Union of all substitutions+ substs+ = M.fromList+ $ concatMap substOf+ $ S.toList fv++ -- Substitution for "x" ==> "x+" - "x-"+ negPos v+ = [(v, (M.fromList [(SVPos v, 1), (SVNeg v, -1)], 0))]++ -- Look at bounds of variables and decide+ substOf v+ = case M.lookup v bs of+ -- Unconstrained, so it can be negative+ Nothing+ -> negPos v+ Just (Nothing, Nothing)+ -> negPos v+ Just (Just 0, _)+ -> []+ -- Nonzero lower bound, so replace: v = v' + n+ Just (Just n, _)+ -> [(v, (M.fromList [(SVLower v, 1)], n)) ]+ _+ -> []++ bounds+ = concatMap linearOfBound+ $ M.toList+ $ C._bounds p++ linearOfBound (v,binds)+ = case binds of+ (_, Just n)+ -> [\s -> (M.fromList [(SV v, 1), (s, 1)], n)]+ _+ -> []++ Constraint cs = C._constraints p+ constrs+ = concatMap linearOfConstraint cs+ linearOfConstraint (C1 lo lin up)+ = let lin' = standardOfLinear lin+ in case (lo,up) of+ (Nothing,Nothing)+ -> []+ (Just lo', Nothing)+ -> [ lt lo' lin' ]+ (Nothing, Just up')+ -> [ gt up' lin' ]+ (Just lo', Just up')+ -> [ lt lo' lin'+ , gt up' lin' ]+++ lt lo lin s+ = ( M.union (M.map negate lin) (M.fromList [(s,1)])+ , negate lo )+ gt up lin s+ = ( M.union lin (M.fromList [(s, 1)])+ , up )++ standardOfLinear (Linear lin)+ = M.mapKeysMonotonic SV lin+++--- 5 <= x1 <= 40+-- ==>+-- x1 subst Lx1+5+-- Lx1 + 5 <= 40+-- ==>+-- Lx1 <= 35++-- assignmentOfMap :: Standard z r c -> M.Map (StandardVar z r) (R c) -> Assignment z r c++++-- Simple helpers ----------++-- | Get the coefficient of a variable in given row+lookupRow :: (Ord z, Ord r, Rep c)+ => StandardRow z r c+ -> StandardVar z r+ -> R c+lookupRow (r,_) v+ = case M.lookup v r of+ Nothing -> 0+ Just vv -> vv++-- | Get objective or basis value of a row+objOfRow+ :: StandardRow z r c+ -> R c+objOfRow = snd+
+ tests/Arbitrary/Assignment.hs view
@@ -0,0 +1,33 @@+module Arbitrary.Assignment where++import Numeric.Limp.Rep++import Arbitrary.Var++import Test.QuickCheck+import Control.Applicative+import Data.Map (fromList)++type Assignment' = Assignment ZVar RVar IntDouble++instance Arbitrary (Z IntDouble) where+ arbitrary = Z <$> arbitrary++instance Arbitrary (R IntDouble) where+ arbitrary = R <$> (fromIntegral <$> (arbitrary :: Gen Int))+++instance Arbitrary (Assignment ZVar RVar IntDouble) where+ arbitrary = arbitrary >>= assignment+++assignment :: Vars -> Gen Assignment'+assignment (Vars zs rs)+ = do zs' <- listOf (elements zs)+ zvs <- infiniteListOf arbitrary++ rs' <- listOf (elements rs)+ rvs <- infiniteListOf arbitrary++ return $ Assignment (fromList $ zs' `zip` zvs) (fromList $ rs' `zip` rvs)+
+ tests/Arbitrary/Program.hs view
@@ -0,0 +1,82 @@+module Arbitrary.Program where++import qualified Numeric.Limp.Program as P+import Numeric.Limp.Rep++import Arbitrary.Var+import Arbitrary.Assignment++import Test.QuickCheck+import Control.Applicative++type Program' = P.Program ZVar RVar IntDouble++data ProgramAss = ProgramAss Program' Assignment'+ deriving Show++instance Arbitrary ProgramAss where+ arbitrary+ = do a <- arbitrary+ ProgramAss <$> program a <*> assignment a++instance Arbitrary Program' where+ arbitrary = arbitrary >>= program+++program :: Vars -> Gen Program'+program vs+ = do dir <- elements [P.Minimise, P.Maximise]+ + obj <- linearR vs+ cons <- constraints vs+ bnds <- listOf (bounds vs)++ return $ P.program dir obj cons bnds+++linearR :: Vars -> Gen (P.Linear ZVar RVar IntDouble P.KR)+linearR (Vars zs rs)+ = do let vs = map Left zs ++ map Right rs+ vs' <- listOf1 (elements vs)+ cs' <- infiniteListOf arbitrary+ summand <- arbitrary+ return $ P.LR (vs' `zip` cs') summand++linearZ :: Vars -> Gen (P.Linear ZVar RVar IntDouble P.KZ)+linearZ (Vars zs _rs)+ = do vs' <- listOf1 (elements zs)+ cs' <- infiniteListOf arbitrary+ summand <- arbitrary+ return $ P.LZ (vs' `zip` cs') summand+++constraints :: Vars -> Gen (P.Constraint ZVar RVar IntDouble)+constraints vs+ = oneof+ [ (P.:==) <$> lR <*> lR+ , (P.:<=) <$> lR <*> lR+ , (P.:<) <$> lZ <*> lZ+ , (P.:>=) <$> lR <*> lR+ , (P.:>) <$> lZ <*> lZ+ , P.Between <$> lR <*> lR <*> lR+ , (P.:&&) <$> constraints vs <*> constraints vs+ , return P.CTrue ]+ where+ lR = linearR vs+ lZ = linearZ vs+++bounds :: Vars -> Gen (P.Bounds ZVar RVar IntDouble)+bounds (Vars zs rs)+ = oneof [bZ, bR]+ where+ bZ = do v <- elements zs+ a <- arbitrary+ b <- arbitrary+ return $ P.BoundZ (a,v,b)++ bR = do v <- elements rs+ a <- arbitrary+ b <- arbitrary+ return $ P.BoundR (a,v,b)+
+ tests/Arbitrary/Var.hs view
@@ -0,0 +1,38 @@+module Arbitrary.Var where+import Test.QuickCheck++data ZVar = ZVar String+ deriving (Eq,Ord)++instance Show ZVar where+ show (ZVar z) = "z$" ++ z++instance Arbitrary ZVar where+ arbitrary+ -- 26 variables should be enough for anyone!+ = do c <- elements ['a'..'z']+ return $ ZVar [c]+++data RVar = RVar String+ deriving (Eq,Ord)++instance Show RVar where+ show (RVar r) = "r$" ++ r++instance Arbitrary RVar where+ arbitrary+ = do c <- elements ['a'..'z']+ return $ RVar [c]+++data Vars = Vars [ZVar] [RVar]+ deriving Show++instance Arbitrary Vars where+ arbitrary+ = do NonEmpty zs <- arbitrary :: Gen (NonEmptyList ZVar)+ NonEmpty rs <- arbitrary :: Gen (NonEmptyList RVar)+ return $ Vars zs rs++
+ tests/BranchExample.hs view
@@ -0,0 +1,94 @@+module BranchExample where++import Numeric.Limp.Rep.Rep as R+import Numeric.Limp.Rep.Arbitrary as R+import Numeric.Limp.Program as P+import Numeric.Limp.Canon as C+import Numeric.Limp.Solve.Simplex.Maps as SM+import Numeric.Limp.Solve.Simplex.StandardForm as ST+import Numeric.Limp.Solve.Branch.Simple as B++import Numeric.Limp.Canon.Pretty+import Debug.Trace++import Control.Applicative++-- Dead simple ones -------------------------+-- x = 2+prog1 :: P.Program String String R.Arbitrary+prog1+ = P.maximise+ -- objective+ (z "x" 1)+ -- subject to+ ( z "x" 2 :<= con 5+ :&& z "x" 4 :>= con 7)+ []++-- x = 1, y = 2+prog2 :: P.Program String String R.Arbitrary+prog2+ = P.minimise+ -- objective+ (z "x" 1 .+. z "y" 1)+ -- subject to+ ( z "x" 2 :<= con 5 -- z "y" 1 .+. con 1+ :&& z "x" 1 :>= con 1 + :&& z "y" 1 :<= con 4+ :&& z "y" 1 :>= con 1)+ [ lowerZ 0 "x" + , lowerZ 0 "y" ]+++xkcd :: Direction -> P.Program String String R.Arbitrary+xkcd dir = P.program dir+ ( z1 mf .+.+ z1 ff .+.+ z1 ss .+.+ z1 hw .+.+ z1 ms .+.+ z1 sp )+ ( z mf mfp .+.+ z ff ffp .+.+ z ss ssp .+.+ z hw hwp .+.+ z ms msp .+.+ z sp spp :== con 1505 )+ [ lowerZ 0 mf+ , lowerZ 0 ff+ , lowerZ 0 ss+ , lowerZ 0 hw+ , lowerZ 0 ms+ , lowerZ 0 sp+ ]+ where+ (mf, mfp) = ("mixed-fruit", 215)+ (ff, ffp) = ("french-fries", 275)+ (ss, ssp) = ("side-salad", 335)+ (hw, hwp) = ("hot-wings", 355)+ (ms, msp) = ("mozzarella-sticks", 420)+ (sp, spp) = ("sampler-plate", 580)++test :: (Show z, Show r, Ord z, Ord r)+ => P.Program z r R.Arbitrary -> IO ()+test prog+ = let prog' = C.program prog+ + simpl p = SM.simplex $ ST.standard p++ solver p+ | st <- ST.standard p+ -- , trace (ppr show show p) True+ , Just s' <- SM.simplex st+ -- , trace ("SAT") True+ , ass <- SM.assignment s'+ = Just ass+ | otherwise+ -- , trace ("unsat") True+ = Nothing+ bb = B.branch solver+ in do + -- putStrLn (show (simpl prog'))+ -- putStrLn (show (solver prog'))+ putStrLn (show (bb prog'))+
+ tests/Convert.hs view
@@ -0,0 +1,19 @@+module Convert where++import Numeric.Limp.Program as P+import Numeric.Limp.Canon as C++import Arbitrary.Program+import Data.Monoid++import Test.Tasty.QuickCheck+import Test.Tasty.TH+++tests = $(testGroupGenerator)++prop_constraints_converted :: ProgramAss -> Bool+prop_constraints_converted (ProgramAss p a)+ = P.checkProgram a p+ == C.checkProgram a (C.program p)+
+ tests/SimplexExample.hs view
@@ -0,0 +1,55 @@+module SimplexExample where++import Numeric.Limp.Rep as R+import Numeric.Limp.Program as P+import Numeric.Limp.Canon as C+import Numeric.Limp.Solve.Simplex.Maps as SM+import Numeric.Limp.Solve.Simplex.StandardForm as ST++import Control.Monad+import qualified Data.Map as M+++data Xs = X1 | X2 | X3+ deriving (Eq, Ord, Show)++prog :: P.Program () Xs R.IntDouble+prog+ = P.maximise+ -- objective+ (r X1 60 .+. r X2 30 .+. r X3 20)+ -- subject to+ ( r X1 8 .+. r X2 6 .+. r X3 1 :<= con 48+ :&& r X1 2 .+. r X2 1.5 .+. r X3 0.5 :<= con 8+ :&& r X1 4 .+. r X2 2 .+. r X3 1.5 :<= con 20+ :&& r X2 1 :<= con 5)+ -- bounds ommitted for now+ [ lowerR 0 X1 , lowerR 0 X2 , lowerR 0 X3 ]+ -- []++test :: IO Bool+test+ = case SM.simplex $ ST.standard $ C.program prog of+ Nothing+ -> do putStrLn "Error: simplex returned Nothing"+ putStrLn (show $ ST.standard $ C.program prog)+ putStrLn (show $ SM.simplex1 $ ST.standard $ C.program prog)+ return False++ Just s+ -> do let (Assignment _ vars,obj) = SM.assignment s+ let vars' = M.toList vars+ let e_vars = [(Right X1, 2.0), (Right X3, 8.0)] :: [(Either () Xs, R IntDouble)]+ let e_obj = -280+ putStrLn "Vars:"+ putStrLn (show vars')+ putStrLn "Obj:"+ putStrLn (show obj)++ when (obj /= e_obj) $+ putStrLn ("Bad objective: should be " ++ show e_obj)+ when (vars' /= e_vars) $+ putStrLn ("Bad vars: should be " ++ show e_vars)++ return (obj == e_obj && vars' == e_vars)+
+ tests/Simplexs.hs view
@@ -0,0 +1,311 @@+module Simplexs where++import Numeric.Limp.Rep.Rep as R+import Numeric.Limp.Rep.Arbitrary as R+import Numeric.Limp.Program as P+import Numeric.Limp.Canon as C+import Numeric.Limp.Solve.Simplex.Maps as SM+import Numeric.Limp.Solve.Simplex.StandardForm as ST++import qualified Data.Map as M+++data Xs = X1 | X2 | X3+ deriving (Eq, Ord, Show)++-- Dead simple ones -------------------------+-- x1 = 10+prog1 :: P.Program () Xs R.Arbitrary+prog1+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10)+ -- bounds omitted for now+ []++-- x1 = 10+prog2 :: P.Program () Xs R.Arbitrary+prog2+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10)+ [ lowerR 0 X1 ]++-- x1 = 0+prog3 :: P.Program () Xs R.Arbitrary+prog3+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10)+ [ lowerR 0 X1 ]++-- Unbounded!+prog4 :: P.Program () Xs R.Arbitrary+prog4+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10)+ []+++-- Two constraints! --------------++-- x = 10+prog5 :: P.Program () Xs R.Arbitrary+prog5+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10+ :&& r X1 1 :>= con (-10))+ []++-- x = -10+prog6 :: P.Program () Xs R.Arbitrary+prog6+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 10+ :&& r X1 1 :>= con (-10))+ []+++-- Now two variables -------------+-- x1 = 20, x2 = 10+prog7 :: P.Program () Xs R.Arbitrary+prog7+ = P.maximise+ -- objective+ (r X1 1 .+. r X2 1)+ -- subject to+ ( r X1 1 :<= r X2 2+ :&& r X2 1 :<= con 10)+ [lowerR 0 X1, lowerR 0 X2]++-- x1 = 20, x2 = 10+prog8 :: P.Program () Xs R.Arbitrary+prog8+ = P.maximise+ -- objective+ (r X1 1 .+. r X2 1)+ -- subject to+ ( r X1 1 :<= r X2 2+ :&& r X2 1 :<= con 10)+ [] -- [lowerR 0 X1, lowerR 0 X2]++-- Something where vars=0 isn't sat ------+-- x1 = 8+prog9 :: P.Program () Xs R.Arbitrary+prog9+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :>= con 8 + :&& r X1 1 :<= con 10)+ [lowerR 0 X1]++-- x1 = 10+prog10 :: P.Program () Xs R.Arbitrary+prog10+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :>= con 8 + :&& r X1 1 :<= con 10)+ [lowerR 0 X1]++++-- An equality constraint ------------+-- x1 = 10+prog11 :: P.Program () Xs R.Arbitrary+prog11+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :== con 10 )+ [lowerR 0 X1]++-- x1 = 10+prog12 :: P.Program () Xs R.Arbitrary+prog12+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :== con 10 )+ [lowerR 0 X1]+++-- From wikipedia ----------------+-- x1 = 2.142..., x3 = 3.571...+prog13 :: P.Program () Xs R.Arbitrary+prog13+ = P.minimise+ -- objective+ (r X1 (-2) .+. r X2 (-3) .+. r X3 (-4))+ -- subject to+ ( r X1 3 .+. r X2 2 .+. r X3 1 :== con 10+ :&& r X1 2 .+. r X2 5 .+. r X3 3 :== con 15)+ [lowerR 0 X1+ ,lowerR 0 X2+ ,lowerR 0 X3]++-- x1 = 1.818..., x2 = 2.272...+prog14 :: P.Program () Xs R.Arbitrary+prog14+ = P.maximise+ -- objective+ (r X1 (-2) .+. r X2 (-3) .+. r X3 (-4))+ -- subject to+ ( r X1 3 .+. r X2 2 .+. r X3 1 :== con 10+ :&& r X1 2 .+. r X2 5 .+. r X3 3 :== con 15)+ [lowerR 0 X1+ ,lowerR 0 X2+ ,lowerR 0 X3]++-- An equality constraint on unconstrained (+-) ------------+-- x1 = 10+prog15 :: P.Program () Xs R.Arbitrary+prog15+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :== con 10 )+ []++-- A lower bound greater than zero ------------+-- x1 = 5+prog16 :: P.Program () Xs R.Arbitrary+prog16+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 30 )+ [lowerR 5 X1]++-- Lower and upper bounds -------+-- x1 = 5+prog17 :: P.Program () Xs R.Arbitrary+prog17+ = P.minimise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 30 )+ [lowerUpperR 5 X1 10]+-- x1 = 10+prog18 :: P.Program () Xs R.Arbitrary+prog18+ = P.maximise+ -- objective+ (r X1 1)+ -- subject to+ ( r X1 1 :<= con 30 )+ [lowerUpperR 5 X1 10]++-- x1 = 1, x2 = 2+prog19 :: P.Program () Xs R.Arbitrary+prog19+ = P.minimise+ (r X1 1 .+. r X2 1)+ ( r X1 2 :<= r X2 1+ :&& r X1 1 :>= con 1)+ [ lowerR 0 X1+ , lowerR 0 X2]+++-- error uncovered by branch -------+-- x1 = 1+-- x2 = 1.870...+prog20 :: P.Program () Xs R.Arbitrary+prog20+ = P.minimise+ -- x1 = mozzarella+ -- x2 = sampler plate+ (r1 X1 .+. r1 X2)+ (r X1 420 .+. r X2 580 :== con 1505)+ [ lowerR 1 X1+ , lowerUpperR 0 X2 2 ]++{-+Minimize+ 1.0 "french-fries" + 1.0 "hot-wings" + 1.0 "mixed-fruit" + 1.0 "mozzarella-sticks" + 1.0 "sampler-plate" + 1.0 "side-salad"+Subject to+ -275.0 "french-fries" - 355.0 "hot-wings" - 215.0 "mixed-fruit" - 420.0 "mozzarella-sticks" - 580.0 "sampler-plate" - 335.0 "side-salad" >= -1505.0+ -275.0 "french-fries" - 355.0 "hot-wings" - 215.0 "mixed-fruit" - 420.0 "mozzarella-sticks" - 580.0 "sampler-plate" - 335.0 "side-salad" <= -1505.0++Bounds+ 0.0 <= "french-fries"+ 0.0 <= "hot-wings"+ 0.0 <= "mixed-fruit"+ 1.0 <= "mozzarella-sticks"+ 0.0 <= "sampler-plate" <= 2.0+ 0.0 <= "side-salad"+-}++-- nonzero lower bound with non-1 coeff+-- x1 = 2.5+prog21 :: P.Program () Xs R.Arbitrary+prog21+ = P.minimise+ (r1 X1)+ (r X1 2 :>= con 5)+ [ lowerR 1 X1 ]++-- eq bound with non-1 coeff+-- x1 = 1, x2 = 3+prog22 :: P.Program () Xs R.Arbitrary+prog22+ = P.minimise+ (r1 X1 .+. r1 X2)+ (r X1 2 .+. r X2 1 :>= con 5)+ [ lowerUpperR 1 X1 1+ , lowerR 0 X2]+++std :: (Ord z, Ord r, Rep c) => P.Program z r c -> Standard z r c+std = ST.standard . C.program+++++test :: P.Program () Xs R.Arbitrary -> IO Bool+test p+ = case SM.simplex $ ST.standard $ C.program p of+ Nothing+ -> do putStrLn "Error: simplex returned Nothing"+ putStrLn (show $ ST.standard $ C.program p)+ putStrLn (show $ SM.simplex1 $ ST.standard $ C.program p)+ return False++ Just s+ -> do let (Assignment _ vars,obj) = SM.assignment s+ let vars' = M.toList vars++ putStrLn (show $ ST.standard $ C.program p)+ putStrLn (show $ SM.simplex1 $ ST.standard $ C.program p)++ putStrLn "Vars:"+ putStrLn (show vars')+ putStrLn "Obj:"+ putStrLn (show obj)++ return True+
+ tests/Simplify.hs view
@@ -0,0 +1,110 @@+module Simplify where++import Numeric.Limp.Program as P+import Numeric.Limp.Canon as C+import Numeric.Limp.Canon.Simplify as CS+import Numeric.Limp.Canon.Simplify.Subst as CS+import Numeric.Limp.Canon.Simplify.Bounder as CS+import Numeric.Limp.Canon.Simplify.Crunch as CS++import Numeric.Limp.Canon.Pretty++import Arbitrary.Assignment as Arb+import Arbitrary.Var as Arb+import Arbitrary.Program as Arb+import Data.Monoid++import Test.Tasty.QuickCheck+import Test.Tasty.TH++import Debug.Trace++tests = $(testGroupGenerator)++prop_bounder :: ProgramAss -> Property+prop_bounder (ProgramAss p a)+ = let cp = C.program p+ cp' = CS.bounderProgram cp+ valcp = C.checkProgram a cp+ valcp'+ | Right p' <- cp'+ = C.checkProgram a p'+ -- Infeasible, so assignment is false+ | otherwise+ = False+ in counterexample+ (unlines+ [ "CP: " ++ show cp+ , "CP':" ++ show cp'+ , "Val: " ++ show (valcp, valcp')])+ $ valcp == valcp'+++prop_crunch :: ProgramAss -> Property+prop_crunch (ProgramAss p a)+ = let cp = C.program p+ cp' = CS.crunchProgram cp+ valcp = C.checkProgram a cp+ valcp' = C.checkProgram a cp'+ in counterexample+ (unlines+ [ "CP: " ++ show cp+ , "CP':" ++ show cp'+ , "Val: " ++ show (valcp, valcp')])+ $ valcp == valcp'+++-- | I don't think this property is very interesting.+-- The real property should be something like:+--+-- > solve cp == solve (simplify cp)+--+prop_simplify :: Program' -> Property+prop_simplify p+ = let cp = C.program p+ simp = CS.simplify cp+ in case simp of+ Left _+ -> property True+ Right (a', cp')+ -> let valcp = C.checkProgram a' cp+ valcp' = C.checkProgram a' cp'+ in counterexample+ (unlines+ [ "CP: " ++ show cp+ , "CP':" ++ show cp'+ , "Ass:" ++ show a'+ , "Val: " ++ show (valcp, valcp')])+ $ valcp == valcp'+++prop_subst_linear :: Vars -> Property+prop_subst_linear vs+ = forAll (Arb.linearR vs) $ \f ->+ forAll (Arb.assignment vs) $ \a ->+ forAll (Arb.assignment vs) $ \b ->+ let (fc, _) = C.linear f+ (fc', c') = substLinear a fc+ in C.evalR (a <> b) fc == C.evalR b fc' + c'+++-- subst can actually make a failing program pass.+-- so this test needs to be implication, not equivalence.+prop_subst_program :: Vars -> Property+prop_subst_program vs+ = forAll (Arb.program vs) $ \f ->+ forAll (Arb.assignment vs) $ \a ->+ forAll (Arb.assignment vs) $ \b ->+ let fc = C.program f+ fc' = substProgram a fc+ both = a <> b+ valcp = C.checkProgram both fc + valcp'= C.checkProgram b fc'+ in counterexample + (unlines+ [ "CP: " ++ show fc+ , "CP':" ++ show fc'+ , "Ass:" ++ show both+ , "Val: " ++ show (valcp, valcp')])+ $ if valcp then valcp' else True+