packages feed

linear-programming (empty) → 0.0

raw patch · 8 files changed

+583/−0 lines, 8 filesdep +QuickCheckdep +basedep +comfort-arraysetup-changed

Dependencies added: QuickCheck, base, comfort-array, non-empty, random, transformers, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2023++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. 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.+3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS 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 AUTHORS 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.
+ Makefile view
@@ -0,0 +1,4 @@+run-test:+	runhaskell Setup configure --user --enable-tests+	runhaskell Setup build+	runhaskell Setup haddock
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ linear-programming.cabal view
@@ -0,0 +1,44 @@+Cabal-Version:    2.2+Name:             linear-programming+Version:          0.0+License:          BSD-3-Clause+License-File:     LICENSE+Author:           Henning Thielemann <haskell@henning-thielemann.de>+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>+Category:         Math+Tested-With:      GHC ==8.6.5+Build-Type:       Simple+Synopsis:         Linear Programming basic definitions+Description:+  Basic types and generic functions for use in the packages+  @coinor-clp@ and @comfort-glpk@.+Extra-Source-Files:+  Makefile++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/linear-programming/++Source-Repository head+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/linear-programming/++Library+  Build-Depends:+    comfort-array >=0.4 && <0.6,+    QuickCheck >=2.1 && <3,+    random >=1.0 && <1.3,+    transformers >=0.3 && <0.7,+    non-empty >=0.3.2 && <0.4,+    utility-ht >=0.0.16 && <0.1,+    base >=4.5 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   src+  Default-Language: Haskell98+  Exposed-Modules:+    Numeric.LinearProgramming.Common+    Numeric.LinearProgramming.Format+    Numeric.LinearProgramming.Monad+    Numeric.LinearProgramming.Test
+ src/Numeric/LinearProgramming/Common.hs view
@@ -0,0 +1,75 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.LinearProgramming.Common (+   Term(..), (.*),+   Inequality(..),+   Bound(..),+   Bounds,+   Constraints,+   Direction(..),+   Objective,+   free, (<=.), (>=.), (==.), (>=<.),+   objectiveFromTerms,+   ) where++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable (Array)++++data Term a ix = Term a ix+   deriving (Show)+++infix 7 .*++(.*) :: a -> ix -> Term a ix+(.*) = Term+++data Inequality x = Inequality x Bound+   deriving Show++data Bound =+     LessEqual Double+   | GreaterEqual Double+   | Between Double Double+   | Equal Double+   | Free+   deriving Show++instance Functor Inequality where+   fmap f (Inequality x bnd)  =  Inequality (f x) bnd++type Bounds ix = [Inequality ix]++type Constraints a ix = [Inequality [Term a ix]]++data Direction = Minimize | Maximize+   deriving (Eq, Enum, Bounded, Show)++type Objective sh = Array sh Double++++infix 4 <=., >=., >=<., ==.++(<=.), (>=.), (==.) :: x -> Double -> Inequality x+x <=. bnd = Inequality x $ LessEqual bnd+x >=. bnd = Inequality x $ GreaterEqual bnd+x ==. bnd = Inequality x $ Equal bnd++(>=<.) :: x -> (Double,Double) -> Inequality x+x >=<. bnd = Inequality x $ uncurry Between bnd++free :: x -> Inequality x+free x = Inequality x Free++++objectiveFromTerms ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   sh -> [Term Double ix] -> Objective sh+objectiveFromTerms sh =+   Array.fromAssociations 0 sh . map (\(Term x ix) -> (ix,x))
+ src/Numeric/LinearProgramming/Format.hs view
@@ -0,0 +1,97 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.LinearProgramming.Format (+   Identifier,+   mathProg,+   ) where++import qualified Numeric.LinearProgramming.Common as LP+import Numeric.LinearProgramming.Common+         (Bound(..), Inequality(Inequality),+          Bounds, Direction(..), Objective, (.*))++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List as List++import Text.Printf (printf)++import Prelude hiding (sum)++++type Term = LP.Term Double++type Constraints ix = LP.Constraints Double ix+++class Identifier ix where+   identifier :: ix -> String++instance Identifier Char where+   identifier x = [x]++instance Identifier c => Identifier [c] where+   identifier = concatMap identifier++instance Identifier Int where+   identifier = printf "x%d"++instance Identifier Integer where+   identifier = printf "x%d"+++bound :: (Identifier ix) => Inequality ix -> String+bound (Inequality ix bnd) =+   printf "var %s%s;" (identifier ix) $+   case bnd of+      LessEqual up -> printf ", <=%f" up+      GreaterEqual lo -> printf ", >=%f" lo+      Between lo up -> printf ", >=%f, <=%f" lo up+      Equal x -> printf ", =%f" x+      Free -> ""+++sum :: (Identifier ix) => [Term ix] -> String+sum [] = "0"+sum xs =+   let formatTerm (LP.Term c ix) = printf "%f*%s" c (identifier ix) in+   List.intercalate "+" $ map formatTerm xs++constraint :: (Identifier ix) => Inequality [Term ix] -> String+constraint (Inequality terms bnd) =+   let sumStr = sum terms in+   case bnd of+      LessEqual up -> printf "%s <= %f" sumStr up+      GreaterEqual lo -> printf "%f <= %s" lo sumStr+      Between lo up -> printf "%f <= %s <= %f" lo sumStr up+      Equal x -> printf "%s = %f" sumStr x+      Free -> sumStr++direction :: Direction -> String+direction Minimize = "minimize"+direction Maximize = "maximize"++objective ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Identifier ix) =>+   Objective sh -> String+objective =+   sum . map (\(ix,c) -> c .* ix) . Array.toAssociations++mathProg ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Identifier ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> [String]+mathProg bounds constrs (dir,obj) =+   map bound bounds +++   "" :+   direction dir :+   printf "value: %s;" (objective obj) :+   "" :+   "subject to" :+   zipWith+      (\k constr -> printf "constr%d: %s;" k $ constraint constr)+      [(0::Int)..] constrs +++   "" :+   "end;" :+   []
+ src/Numeric/LinearProgramming/Monad.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{- |+Generic implementation of a monad that collects constraints+over multiple stages.+It can be used to test solvers that allow for warm start+or for solvers that do not allow for warm start at all+(like GLPK's interior point solver).+-}+module Numeric.LinearProgramming.Monad (+   T,+   run,+   lift,+   ) where++import Numeric.LinearProgramming.Common++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape++import qualified Control.Monad.Trans.RWS as MRWS+import Control.Monad (when)+++newtype T sh a =+   Cons (MRWS.RWS+            (sh, Bounds (Shape.Index sh))+            ()+            (Constraints Double (Shape.Index sh))+            a)+      deriving (Functor, Applicative, Monad)+++run ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   sh -> Bounds ix -> T sh a -> a+run shape bounds (Cons act) =+   fst $ MRWS.evalRWS act (shape, bounds) []++lift ::+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Bounds ix -> Constraints Double ix -> (Direction, Objective sh) -> a) ->+   Constraints Double ix -> (Direction, Objective sh) -> T sh a+lift solver constrs dirObj@(_dir,obj) = Cons $ do+   (shape,bounds) <- MRWS.ask+   when (shape /= Array.shape obj) $+      error "LinearProgramming.Monad.solve: objective shape mismatch"+   MRWS.modify (constrs++)+   allConstrs <- MRWS.get+   return $ solver bounds allConstrs dirObj
+ src/Numeric/LinearProgramming/Test.hs view
@@ -0,0 +1,282 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.LinearProgramming.Test (+   Element,+   forAllOrigin,+   forAllProblem,+   genObjective,+   forAllObjectives,+   successiveObjectives,+   approxReal,+   approx,+   checkFeasibility,+   affineCombination,+   scalarProduct,+   ) where++import qualified Numeric.LinearProgramming.Common as LP+import Numeric.LinearProgramming.Common ((<=.), (>=.), (.*))++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((.&&.))+import System.Random (Random)++import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Ix as Ix+import Data.Array.Comfort.Storable (Array, (!))+import Data.Traversable (sequenceA, for)+import Data.Tuple.HT (mapSnd)+import Data.Maybe (fromMaybe)+import Data.Int (Int64)++import Control.Applicative (liftA2)++import Text.Printf (PrintfArg, printf)++import Foreign.Storable (Storable)++++type Term = LP.Term Double+type Constraints ix = LP.Constraints Double ix+++{- |+Generate constraints in the form of a polyhedron+which contains warrantedly the zero vector.+That is, there is an admissible solution.+In order to assert that the polyhedron is closed,+we bound all variables by a hypercube.+-}+genProblem ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (LP.Bounds ix, Constraints ix)+genProblem origin =+   liftA2 (,)+      (for (Array.toAssociations origin) $ \(ix,x) ->+         LP.Inequality ix <$>+         liftA2 LP.Between+            (doubleFromElement . (x+) <$> QC.choose (-100,-50))+            (doubleFromElement . (x+) <$> QC.choose (50,100)))+      (do+         numConstraints <- QC.choose (1,20)+         QC.vectorOf numConstraints $ do+            ixs <- QC.sublistOf $ Shape.indices $ Array.shape origin+            terms <- for ixs $ \ix -> do+               coeff <- QC.choose (-10,10)+               return (coeff, ix)+            let offset = scalarProductTerms terms origin+            let deviation = 25+            LP.Inequality+               (map (uncurry ((.*) . doubleFromElement)) terms)+               <$>+               QC.oneof (+                  (do bound <- QC.choose (offset-deviation, offset+deviation)+                      return $+                         if bound > offset+                            then LP.LessEqual    $ doubleFromElement bound+                            else LP.GreaterEqual $ doubleFromElement bound) :+                  (liftA2 LP.Between+                     (doubleFromElement <$>+                        QC.choose (offset-deviation, offset))+                     (doubleFromElement <$>+                        QC.choose (offset, offset+deviation))) :+                  []))++scalarProductTerms ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Num a) =>+   [(a,ix)] -> Array sh a -> a+scalarProductTerms terms origin =+   sum $ map (\(coeff, ix) -> coeff * origin!ix) terms++genVarShape :: QC.Gen (Shape.Range Char)+genVarShape = Shape.Range 'a' <$> QC.choose ('a','j')++genOrigin :: QC.Gen (Array (Shape.Range Char) Int64)+genOrigin = genVector =<< genVarShape++_genOrigin :: QC.Gen (Array (Shape.Range Char) Double)+_genOrigin = genVector =<< genVarShape+++_shrinkVarShape :: Shape.Range Char -> [Shape.Range Char]+_shrinkVarShape (Shape.Range from to) =+   if from<to then [Shape.Range from (pred to)] else []++shrinkOrigin ::+   (Storable a) => Array (Shape.Range Char) a -> [Array (Shape.Range Char) a]+shrinkOrigin vec =+   case Array.shape vec of+      Shape.Range from to ->+         if from<to+            then [Array.sample (Shape.Range from (pred to)) (vec!)]+            else []+++forAllOrigin ::+   (QC.Testable prop) =>+   (Array (Shape.Range Char) Int64 -> prop) -> QC.Property+forAllOrigin = QC.forAllShrink genOrigin shrinkOrigin+++class (Storable a, Random a, Num a, Ord a) => Element a where+   doubleFromElement :: a -> Double++instance Element Double where+   doubleFromElement = id++instance Element Int64 where+   doubleFromElement = fromIntegral++genObjective ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (LP.Direction, LP.Objective sh)+genObjective origin =+   liftA2 (,) QC.arbitraryBoundedEnum+      (fmap (Array.map doubleFromElement . flip asTypeOf origin) $+       genVector $ Array.shape origin)++genVector :: (Shape.Indexed sh, Element a) => sh -> QC.Gen (Array sh a)+genVector shape =+   fmap Array.fromBoxed $ sequenceA $+   BoxedArray.fromAssociations (QC.choose (-10,10)) shape []+--    BoxedArray.constant shape (QC.choose (-10,10))++shrinkProblem ::+   (LP.Bounds ix, Constraints ix) ->+   [(LP.Bounds ix, Constraints ix)]+shrinkProblem (bounds, constraints) =+   map (\shrinked -> (bounds, shrinked)) $+   filter (not . null) $ QC.shrinkList (const []) constraints++forAllProblem ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>+   (QC.Testable prop, Element a) =>+   Array sh a -> (LP.Bounds ix -> Constraints ix -> prop) -> QC.Property+forAllProblem origin =+   QC.forAllShrink (genProblem origin) shrinkProblem . uncurry+++genObjectives ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (NonEmpty.T [] (LP.Direction, [Term ix]))+genObjectives origin = do+   let shape = Array.shape origin+   let stageRange :: (Int,Int)+       stageRange = (0,3)+   stages <- for (Shape.indices shape) $ \ix -> (,) ix <$> QC.choose stageRange+   let varSets =+         fromMaybe (error "there should be at least one stage") $+         NonEmpty.fetch $+         filter (not . null) $+         map (\k -> map fst $ filter ((k==) . snd) stages) $+         Ix.range stageRange+   let asTypeOfElement :: a -> f a -> a+       asTypeOfElement = const+   for varSets $ \varSet ->+      liftA2 (,)+         QC.arbitraryBoundedEnum+         (for varSet $ \ix ->+            (.*ix) . doubleFromElement+               <$> QC.choose (-10, 10 `asTypeOfElement` origin))++shrinkObjectives ::+   NonEmpty.T [] (LP.Direction, [Term ix]) ->+   [NonEmpty.T [] (LP.Direction, [Term ix])]+shrinkObjectives (NonEmpty.Cons obj objs) =+   map (NonEmpty.Cons obj) $+   QC.shrinkList+      (\(dir,terms) ->+         map ((,) dir) $ filter (not . null) $+         QC.shrinkList (const []) terms)+      objs++forAllObjectives ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>+   (QC.Testable prop, Element a) =>+   Array sh a ->+   (NonEmpty.T [] (LP.Direction, [Term (Shape.Index sh)]) -> prop) ->+   QC.Property+forAllObjectives origin =+   QC.forAllShrink (genObjectives origin) shrinkObjectives++constraintsFromSolution ::+   Double -> (LP.Direction, x) -> Double -> [LP.Inequality x]+constraintsFromSolution tol (dir,obj) opt =+   case dir of+      LP.Minimize -> [obj <=. opt + tol]+      LP.Maximize -> [obj >=. opt - tol]++successiveObjectives ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Array sh a -> Double ->+   NonEmpty.T [] (LP.Direction, [Term ix]) ->+   ((LP.Direction, LP.Objective sh),+    [(Double -> Constraints ix, (LP.Direction, LP.Objective sh))])+successiveObjectives origin tol xs =+   let shape = Array.shape origin in+   (mapSnd (LP.objectiveFromTerms shape) $ NonEmpty.head xs,+    NonEmpty.mapAdjacent+      (\(dir,obj) y1 ->+         (constraintsFromSolution tol (dir,obj),+          mapSnd (LP.objectiveFromTerms shape) y1))+      xs)+++approxReal :: (Ord a, Num a) => a -> a -> a -> Bool+approxReal tol x y = abs (x-y) <= tol++approx :: (PrintfArg a, Ord a, Num a) => String -> a -> a -> a -> QC.Property+approx name tol x y =+   QC.counterexample (printf "%s: %f - %f" name x y) (approxReal tol x y)++++checkBound :: Double -> LP.Bound -> Double -> QC.Property+checkBound tol bound x =+   QC.counterexample (show (x, bound)) $+   case bound of+      LP.LessEqual up -> x<=up+tol+      LP.GreaterEqual lo -> x>=lo-tol+      LP.Between lo up -> lo-tol<=x && x<=up+tol+      LP.Equal y -> approxReal tol x y+      LP.Free -> True++checkBounds ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Double -> LP.Bounds ix -> Array sh Double -> QC.Property+checkBounds tol bounds sol =+   QC.conjoin $ map (\(ix,bnd) -> checkBound tol bnd (sol!ix)) $+   BoxedArray.toAssociations $+   BoxedArray.fromAssociations (LP.GreaterEqual 0) (Array.shape sol) $+   map (\(LP.Inequality ix bnd) -> (ix,bnd)) bounds++checkContraint ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Double -> LP.Inequality [LP.Term Double ix] -> Array sh Double -> QC.Property+checkContraint tol (LP.Inequality terms bnd) sol =+   checkBound tol bnd $+   scalarProductTerms (map (\(LP.Term c ix) -> (c,ix)) terms) sol++checkFeasibility ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Double -> LP.Bounds ix -> Constraints ix -> Array sh Double -> QC.Property+checkFeasibility tol bounds constrs sol =+   checkBounds tol bounds sol+   .&&.+   QC.conjoin (map (flip (checkContraint tol) sol) constrs)+++affineCombination ::+   (Shape.C sh, Eq sh, Storable a, Num a) =>+   a -> Array sh a -> Array sh a -> Array sh a+affineCombination c x y =+   Array.zipWith (+) (Array.map ((1-c)*) x) (Array.map (c*) y)++scalarProduct ::+   (Shape.C sh, Eq sh, Storable a, Num a) =>+   Array sh a -> Array sh a -> a+scalarProduct x y = Array.sum $ Array.zipWith (*) x y