packages feed

coinor-clp (empty) → 0.0

raw patch · 16 files changed

+1336/−0 lines, 16 filesdep +QuickCheckdep +basedep +coinor-clpsetup-changed

Dependencies added: QuickCheck, base, coinor-clp, comfort-array, deepseq, doctest-exitcode-stdio, doctest-lib, linear-programming, 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,12 @@+run-test:	update-test+	runhaskell Setup configure --user --enable-tests+	runhaskell Setup build+	runhaskell Setup haddock+	runhaskell Setup test coinor-clp-test --show-details=streaming++	runhaskell Setup configure --user -fdebug+	runhaskell Setup build++update-test:+	doctest-extract-0.1 -i src/ -o test/ --executable-main=Main.hs \+	  Numeric.COINOR.CLP Numeric.COINOR.CLP.Monad
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ cbits/support.cpp view
@@ -0,0 +1,99 @@+#define CLP_EXTERN_C+#include "coin/ClpSimplex.hpp"+#include "support.h"+++ClpPlusMinusOneMatrix *+  Clp_newPlusMinusOneMatrix (+    int numberRows,+    int numberColumns,+    bool columnOrdered,+    const int *indices,+    const CoinBigIndex *startPositive,+    const CoinBigIndex *startNegative+  ) {++  return+    new ClpPlusMinusOneMatrix (+      numberRows, numberColumns,+      columnOrdered, indices,+      startPositive, startNegative+    );+}++void+  Clp_deletePlusMinusOneMatrix (+    ClpPlusMinusOneMatrix *matrix+  ) {+  delete matrix;+}++void+  Clp_loadProblemFromMatrix (+    Clp_Simplex *model,+    const ClpMatrixBase *matrix,+    const double *collb,+    const double *colub,+    const double *obj,+    const double *rowlb,+    const double *rowub,+    const double *rowObjective+  ) {++  model->model_->loadProblem (+      *matrix,+      collb, colub, obj,+      rowlb, rowub, rowObjective+    );+}+++CoinPackedMatrix *+  Clp_newCoinPackedMatrix (+    const bool colordered,+    const int minor,+    const int major,+    const CoinBigIndex numels,+    const double *elem,+    const int *ind,+    const CoinBigIndex *start,+    const int *len+  ) {+  return+    new CoinPackedMatrix (+      colordered,+      minor, major, numels,+      elem, ind, start, len+    );+}++void+  Clp_deleteCoinPackedMatrix (+    CoinPackedMatrix *matrix+  ) {+  delete matrix;+}++void+  Clp_loadProblemFromCoinMatrix (+    Clp_Simplex *model,+    const CoinPackedMatrix *matrix,+    const double *collb,+    const double *colub,+    const double *obj,+    const double *rowlb,+    const double *rowub,+    const double *rowObjective+  ) {++  model->model_->loadProblem (+      *matrix,+      collb, colub, obj,+      rowlb, rowub, rowObjective+    );+}+++void Clp_dumpMatrix (Clp_Simplex *model) {+  model->model_->matrix()->dumpMatrix ();+}
+ coinor-clp.cabal view
@@ -0,0 +1,107 @@+Cabal-Version:    2.2+Name:             coinor-clp+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 using COIN-OR/CLP and comfort-array+Description:+  Simple interface to linear programming functions provided by COIN-OR+  using the flexible Array shape framework from @comfort-array@.+  .+  E.g. you can use @Shape.Tuple@ to convert safely+  between nested tuples and arrays with the same number of elements.+  .+  > type X = Shape.Element+  > type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+  >+  > case Shape.indexTupleFromShape (Shape.static :: PairShape) of+  >   (posIx,negIx) ->+  >     case mapSnd (mapSnd Array.toTuple) <$>+  >          LP.simplex [] [[1.*posIx, (-1).*negIx] ==. 314]+  >            (LP.Minimize,+  >             Array.fromTuple (23,42) :: Array PairShape Double)+  >       of+  >         (Right (LP.Optimal, (absol, (pos, neg)))) ->+  >           printf "absol %f,  pos %f, neg %f\n" absol pos neg+  >         _ -> fail "COINOR solver failed"+  .+  Alternatives: @comfort-glpk@, @hmatrix-glpk@, @glpk-hs@++Extra-Source-Files:+  Makefile+  include/support.h++Flag debug+  Description: Enable debug output+  Default:     False+  Manual:      True++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/coinor-clp/++Source-Repository head+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/coinor-clp/++Library+  Build-Depends:+    linear-programming >=0.0 && <0.1,+    comfort-array >=0.4 && <0.6,+    QuickCheck >=2 && <3,+    deepseq >=1.3 && <1.5,+    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+  Default-Language: Haskell98+  Hs-Source-Dirs:   src+  If flag(debug)+    Hs-Source-Dirs: src/debug-on+  Else+    Hs-Source-Dirs: src/debug-off+  Exposed-Modules:+    Numeric.COINOR.CLP+    Numeric.COINOR.CLP.Monad+  Other-Modules:+    Numeric.COINOR.CLP.FFI+    Numeric.COINOR.CLP.Private+    Numeric.COINOR.CLP.Debug++  PkgConfig-Depends: clp+  Include-Dirs: include+  Extra-Libraries: stdc+++  Cxx-Sources:+    cbits/support.cpp++Test-Suite coinor-clp-test+  Type: exitcode-stdio-1.0+  Build-Depends:+    coinor-clp,+    linear-programming,+    comfort-array >=0.5.2,+    transformers,+    non-empty,+    utility-ht,+    doctest-exitcode-stdio >=0.0 && <0.1,+    doctest-lib >=0.1 && <0.2,+    QuickCheck >=2.1 && <3,+    random >=1.0 && <1.3,+    base >=4.5 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   test+  Default-Language: Haskell98+  Main-Is: Main.hs+  Other-Modules:+    Test.Numeric.COINOR.CLP+    Test.Numeric.COINOR.CLP.Monad+    Test.Numeric.COINOR.CLP.Utility
+ include/support.h view
@@ -0,0 +1,78 @@+#ifndef COINOR_HS_SUPPORT_H+#define COINOR_HS_SUPPORT_H++#include "coin/Clp_C_Interface.h"+#include "coin/ClpPlusMinusOneMatrix.hpp"+++#ifdef __cplusplus+extern "C" {+#endif+++ClpPlusMinusOneMatrix *+  Clp_newPlusMinusOneMatrix (+    int numberRows,+    int numberColumns,+    bool columnOrdered,+    const int *indices,+    const CoinBigIndex *startPositive,+    const CoinBigIndex *startNegative+  );++void+  Clp_deletePlusMinusOneMatrix (+    ClpPlusMinusOneMatrix *matrix+  );++void+  Clp_loadProblemFromMatrix (+    Clp_Simplex *model,+    const ClpMatrixBase *matrix,+    const double *collb,+    const double *colub,+    const double *obj,+    const double *rowlb,+    const double *rowub,+    const double *rowObjective+  );+++CoinPackedMatrix *+  Clp_newCoinPackedMatrix (+    const bool colordered,+    const int minor,+    const int major,+    const CoinBigIndex numels,+    const double *elem,+    const int *ind,+    const CoinBigIndex *start,+    const int *len+  );++void+  Clp_deleteCoinPackedMatrix (+    CoinPackedMatrix *matrix+  );++void+  Clp_loadProblemFromCoinMatrix (+    Clp_Simplex *model,+    const CoinPackedMatrix *matrix,+    const double *collb,+    const double *colub,+    const double *obj,+    const double *rowlb,+    const double *rowub,+    const double *rowObjective+  );++void Clp_dumpMatrix (Clp_Simplex *model);++++#ifdef __cplusplus+} /* extern "C" */+#endif++#endif /* COINOR_HS_SUPPORT_H */
+ src/Numeric/COINOR/CLP.hs view
@@ -0,0 +1,201 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.COINOR.CLP (+   simplex,+   LP.Direction(..),+   PlusMinusOne(..),+   Term(..), (LP..*),+   Constraints,+   LP.free, (LP.<=.), (LP.>=.), (LP.==.), (LP.>=<.),+   Method, Priv.dual, Priv.primal,+   Priv.initialSolve, Priv.initialDualSolve, Priv.initialPrimalSolve,+   Priv.initialBarrierSolve, Priv.initialBarrierNoCrossSolve,+   FailureType(..),+   Result,+   ) where++import qualified Numeric.COINOR.CLP.FFI as FFI+import qualified Numeric.COINOR.CLP.Debug as Debug+import qualified Numeric.COINOR.CLP.Private as Priv+import Numeric.COINOR.CLP.Private+         (Method(runMethod), Result, FailureType(..),+          runContT, withBuffer, false,+          storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,+          storeConstraints, setOptimizationDirection, examineStatus)++import qualified Numeric.LinearProgramming.Common as LP+import Numeric.LinearProgramming.Common+         (Inequality(Inequality), Bounds,+          Term(Term), Constraints, Direction(..), Objective)++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List.HT as ListHT++import qualified Control.Monad.Trans.Cont as MC+import Control.Monad.IO.Class (liftIO)+import Control.Exception (bracket)++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Ptr (Ptr, nullPtr)+import Foreign.C.Types (CDouble)+++{- $setup+>>> import qualified Numeric.COINOR.CLP as LP+>>> import qualified Numeric.LinearProgramming.Test as TestLP+>>> import Numeric.COINOR.CLP (PlusMinusOne(..), (.*), (==.), (<=.), (>=<.))+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>>+>>> import Data.Tuple.HT (mapSnd)+>>>+>>> import qualified Test.QuickCheck as QC+>>> import Test.QuickCheck ((===), (.&&.), (.||.))+>>>+>>> type X = Shape.Element+>>> type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+>>> type TripletShape = Shape.NestedTuple Shape.TupleIndex (X,X,X)+>>>+>>> pairShape :: PairShape+>>> pairShape = Shape.static+>>>+>>> tripletShape :: TripletShape+>>> tripletShape = Shape.static+>>>+>>> approxReal :: (Ord a, Num a) => a -> a -> a -> Bool+>>> approxReal tol x y = abs (x-y) <= tol+>>>+>>> genMethod :: QC.Gen (String, LP.Method)+>>> genMethod = QC.elements $+>>>    ("dual", LP.dual) :+>>>    ("primal", LP.primal) :+>>>    ("initialSolve", LP.initialSolve) :+>>>    ("initialDualSolve", LP.initialDualSolve) :+>>>    ("initialPrimalSolve", LP.initialPrimalSolve) :+>>>    ("initialBarrierSolve", LP.initialBarrierSolve) :+>>>    -- let tests fail+>>>    -- ("initialBarrierNoCrossSolve", LP.initialBarrierNoCrossSolve) :+>>>    []+>>>+>>> forAllMethod ::+>>>    (QC.Testable prop) => (LP.Method -> prop) -> QC.Property+>>> forAllMethod prop = QC.forAllShow genMethod fst (prop . snd)+-}++++data PlusMinusOne = MinusOne | PlusOne deriving (Eq, Show)+++class Coefficient a where+   loadProblem ::+      (Shape.Indexed sh, Shape.Index sh ~ ix) =>+      sh ->+      Constraints a ix ->+      Ptr FFI.Simplex ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      Ptr CDouble -> Ptr CDouble ->+      MC.ContT () IO ()++instance Coefficient Double where+   loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr = do+      (coefficientsShape, coefficientsPtr, indexPtr, startPtr)+         <- storeConstraints shape constrs+      let createMatrix =+            FFI.newCoinPackedMatrix+               false+               (fromIntegral $ Shape.size shape)+               (fromIntegral $ length constrs)+               (fromIntegral $ Shape.size coefficientsShape)+               coefficientsPtr+               indexPtr+               startPtr+               nullPtr+      matrix <- MC.ContT $ bracket createMatrix FFI.deleteCoinPackedMatrix+      liftIO $+         FFI.loadProblemFromCoinMatrix lp matrix+            collbPtr colubPtr+            objPtr+            rowlbPtr rowubPtr+            nullPtr++instance Coefficient PlusMinusOne where+   loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr = do+      let shapeOffset = Shape.offset shape+      let coefficients =+            map+               (\(Inequality terms _bnd) ->+                  ListHT.partition (\(Term c _) -> c == PlusOne) terms)+               constrs+      indexPtr <-+         withBuffer $ Array.vectorFromList $+         concatMap+            (\(positive,negative) ->+               map fromIntegral $+                  map (\(Term _ ix) -> shapeOffset ix) positive+                  +++                  map (\(Term _ ix) -> shapeOffset ix) negative)+            coefficients+      let rowStarts =+            scanl (+) 0 $+            map (\(Inequality terms _bnd) -> length terms) constrs+      startPositivePtr <-+         withBuffer $ Array.vectorFromList $ map fromIntegral rowStarts+      startNegativePtr <-+         withBuffer $ Array.vectorFromList $+         zipWith (\k (pos,_neg) -> fromIntegral $ k + length pos)+            rowStarts coefficients+      let createMatrix =+            FFI.newPlusMinusOneMatrix+               (fromIntegral $ length constrs)+               (fromIntegral $ Shape.size shape)+               (toEnum $ fromEnum False)+               indexPtr startPositivePtr startNegativePtr+      matrix <- MC.ContT $ bracket createMatrix FFI.deletePlusMinusOneMatrix+      liftIO $+         FFI.loadProblemFromMatrix lp matrix+            collbPtr colubPtr+            objPtr+            rowlbPtr rowubPtr+            nullPtr+++{- |+>>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.dual [] [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)+Right (28.0,(5.0,0.0,4.0))+>>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.primal [y >=<. (-12,12)] [[1.*x, (-1).*y] <=. 10, [(-1).*y, (1::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)+Right (116.0,(22.0,12.0,32.0))+>>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.primal [y >=<. (-12,12)] [[PlusOne .* x, MinusOne .* y] <=. 10, [MinusOne .* y, PlusOne .* z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)+Right (116.0,(22.0,12.0,32.0))++prop> forAllMethod $ \method (QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd Array.toTuple <$> LP.simplex method [] [[1.*pos, (-1::Double).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (absol,(posResult,negResult))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False+prop> forAllMethod $ \method target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd Array.toTuple <$> LP.simplex method [] [[1.*pos, (-1::Double).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (absol,(posResult,negResult))) -> QC.counterexample (show(absol,(posResult,negResult))) $ QC.property (approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False++prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right _ -> True; _ -> False+prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right (_,sol) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False+prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right (_,sol) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.1 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False+prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Left _ -> QC.property False; Right (opt,sol) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val+-}+simplex ::+   (Coefficient a, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Method -> Bounds ix -> Constraints a ix ->+   (Direction, Objective sh) -> Result sh+simplex method bounds constrs (dir,obj) =+   unsafePerformIO $+   bracket FFI.newModel FFI.deleteModel $ \lp -> do++   Debug.initLog lp+   let shape = Array.shape obj+   runContT $ do+      objPtr <- withBuffer $ Array.map realToFrac obj+      (collbPtr,colubPtr) <-+         storeBounds $ prepareColumnBoundsArrays shape bounds+      (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs+      loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr+   setOptimizationDirection lp dir+   runMethod method lp+   examineStatus shape lp
+ src/Numeric/COINOR/CLP/FFI.hsc view
@@ -0,0 +1,159 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Numeric.COINOR.CLP.FFI where++import qualified Foreign.C.Types as C+import Foreign.Ptr (Ptr)++import Data.Int (Int32)+++type CDouble = C.CDouble+type CInt = C.CInt+type CBool = C.CBool+++#include "coin/Clp_C_Interface.h"++type BigIndex = #{type CoinBigIndex}++data CoinPackedMatrix = CoinPackedMatrix+data PlusMinusOneMatrix = PlusMinusOneMatrix+data PackedMatrix = PackedMatrix+data MatrixBase = MatrixBase+data Simplex = Simplex+++foreign import ccall "Clp_newPlusMinusOneMatrix"+   newPlusMinusOneMatrix ::+      CInt -> CInt -> CBool -> Ptr CInt ->+      Ptr BigIndex -> Ptr BigIndex ->+      IO (Ptr PlusMinusOneMatrix)++foreign import ccall "Clp_deletePlusMinusOneMatrix"+   deletePlusMinusOneMatrix ::+      Ptr PlusMinusOneMatrix -> IO ()++foreign import ccall "Clp_loadProblemFromMatrix" +   loadProblemFromMatrix ::+      Ptr Simplex ->+      Ptr matrix ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      IO ()+++foreign import ccall "Clp_newCoinPackedMatrix"+   newCoinPackedMatrix ::+      CBool -> CInt -> CInt -> BigIndex ->+      Ptr CDouble -> Ptr CInt -> Ptr BigIndex -> Ptr CInt ->+      IO (Ptr CoinPackedMatrix)++foreign import ccall "Clp_deleteCoinPackedMatrix"+   deleteCoinPackedMatrix ::+      Ptr CoinPackedMatrix -> IO ()++foreign import ccall "Clp_loadProblemFromCoinMatrix" +   loadProblemFromCoinMatrix ::+      Ptr Simplex ->+      Ptr CoinPackedMatrix ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      IO ()+++foreign import ccall "Clp_newModel"+   newModel :: IO (Ptr Simplex)++foreign import ccall "Clp_deleteModel"+   deleteModel :: Ptr Simplex -> IO ()++foreign import ccall "Clp_loadProblem"+   loadProblem ::+      Ptr Simplex ->+      CInt -> CInt ->+      Ptr BigIndex -> Ptr CInt ->+      Ptr CDouble ->+      Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      Ptr CDouble -> Ptr CDouble ->+      IO ()++foreign import ccall "Clp_setOptimizationDirection"+   setOptimizationDirection :: Ptr Simplex -> CDouble -> IO ()++foreign import ccall "Clp_chgObjCoefficients"+   chgObjCoefficients :: Ptr Simplex -> Ptr CDouble -> IO ()++foreign import ccall "Clp_addRows"+   addRows ::+      Ptr Simplex -> CInt -> Ptr CDouble -> Ptr CDouble ->+      Ptr BigIndex -> Ptr CInt -> Ptr CDouble -> IO ()++foreign import ccall "Clp_addColumns"+   addColumns ::+      Ptr Simplex -> CInt -> Ptr CDouble -> Ptr CDouble ->+      Ptr CDouble ->+      Ptr BigIndex -> Ptr CInt -> Ptr CDouble -> IO ()+++foreign import ccall "Clp_initialSolve"+   initialSolve :: Ptr Simplex -> IO CInt+foreign import ccall "Clp_initialDualSolve"+   initialDualSolve :: Ptr Simplex -> IO CInt+foreign import ccall "Clp_initialPrimalSolve"+   initialPrimalSolve :: Ptr Simplex -> IO CInt+foreign import ccall "Clp_initialBarrierSolve"+   initialBarrierSolve :: Ptr Simplex -> IO CInt+foreign import ccall "Clp_initialBarrierNoCrossSolve"+   initialBarrierNoCrossSolve :: Ptr Simplex -> IO CInt++foreign import ccall "Clp_dual"+   dual :: Ptr Simplex -> CInt -> IO CInt++foreign import ccall "Clp_primal"+   primal :: Ptr Simplex -> CInt -> IO CInt+++foreign import ccall "Clp_objectiveValue"+   objectiveValue :: Ptr Simplex -> IO CDouble++foreign import ccall "Clp_getColSolution"+   getColSolution :: Ptr Simplex -> IO (Ptr CDouble)+++foreign import ccall "Clp_status"+   status :: Ptr Simplex -> IO CInt+foreign import ccall "Clp_secondaryStatus"+   secondaryStatus :: Ptr Simplex -> IO CInt+{- Are there a numerical difficulties? -}+foreign import ccall "Clp_isAbandoned"+   isAbandoned :: Ptr Simplex -> IO CInt+{- Is optimality proven? -}+foreign import ccall "Clp_isProvenOptimal"+   isProvenOptimal :: Ptr Simplex -> IO CInt+{- Is primal infeasiblity proven? -}+foreign import ccall "Clp_isProvenPrimalInfeasible"+   isProvenPrimalInfeasible :: Ptr Simplex -> IO CInt+{- Is dual infeasiblity proven? -}+foreign import ccall "Clp_isProvenDualInfeasible"+   isProvenDualInfeasible :: Ptr Simplex -> IO CInt+{- Is the given primal objective limit reached? -}+foreign import ccall "Clp_isPrimalObjectiveLimitReached"+   isPrimalObjectiveLimitReached :: Ptr Simplex -> IO CInt+{- Is the given dual objective limit reached? -}+foreign import ccall "Clp_isDualObjectiveLimitReached"+   isDualObjectiveLimitReached :: Ptr Simplex -> IO CInt+{- Iteration limit reached? -}+foreign import ccall "Clp_isIterationLimitReached"+   isIterationLimitReached :: Ptr Simplex -> IO CInt+++foreign import ccall "Clp_setLogLevel"+   setLogLevel :: Ptr Simplex -> CInt -> IO ()++foreign import ccall "Clp_dumpMatrix"+   dumpMatrix :: Ptr Simplex -> IO ()
+ src/Numeric/COINOR/CLP/Monad.hs view
@@ -0,0 +1,208 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{- |+The monadic interface to CLP allows to optimize+with respect to multiple objectives, successively.+-}+module Numeric.COINOR.CLP.Monad (+   T,+   run,+   simplex,+   Direction(..),+   Priv.dual, Priv.primal,+   ) where++import qualified Numeric.COINOR.CLP.FFI as FFI+import qualified Numeric.COINOR.CLP.Debug as Debug+import qualified Numeric.COINOR.CLP.Private as Priv+import Numeric.COINOR.CLP.Private+         (Method(runMethod), Result,+          runContT, withBuffer,+          storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,+          storeConstraints, setOptimizationDirection, examineStatus)++import Numeric.LinearProgramming.Common+         (Bounds, Constraints, Direction(..), Objective)++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape++import qualified Control.Monad.Trans.Cont as MC+import qualified Control.Monad.Trans.Reader as MR+import Control.Monad.IO.Class (liftIO)+import Control.Monad (when)+import Control.Exception (bracket)++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Ptr (Ptr, nullPtr)+++{- $setup+>>> :set -XTypeFamilies+>>> :set -XTypeOperators+>>> import qualified Numeric.COINOR.CLP.Monad as LP+>>> import qualified Numeric.COINOR.CLP as CLP+>>> import Test.Numeric.COINOR.CLP.Utility (traverse_Lag, traverseLag)+>>> import Test.Numeric.COINOR.CLP (TripletShape, tripletShape, forAllMethod)+>>> import Numeric.COINOR.CLP (Direction, (.*), (<=.))+>>>+>>> import qualified Numeric.LinearProgramming.Monad as LPMonad+>>> import qualified Numeric.LinearProgramming.Test as TestLP+>>> import Numeric.LinearProgramming.Common (Bounds, Objective)+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>> import qualified Data.NonEmpty as NonEmpty+>>> import Data.Array.Comfort.Storable (Array)+>>> import Data.Traversable (Traversable)+>>> import Data.Foldable (Foldable)+>>>+>>> import qualified Control.Monad.Trans.Except as ME+>>>+>>> import qualified Data.List.HT as ListHT+>>> import Data.Tuple.HT (mapSnd)+>>>+>>> import Foreign.Storable (Storable)+>>>+>>> import qualified Test.QuickCheck as QC+>>>+>>>+>>> type Constraints ix = CLP.Constraints Double ix+>>>+>>>+>>> approxSuccession ::+>>>    (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>+>>>    a ->+>>>    Either CLP.FailureType (NonEmpty.T [] (a, Array sh a)) ->+>>>    Either CLP.FailureType (NonEmpty.T [] (a, Array sh a)) ->+>>>    QC.Property+>>> approxSuccession tol x y =+>>>    QC.counterexample (show x) $+>>>    QC.counterexample (show y) $+>>>    case (x,y) of+>>>       (Left sx, Left sy) -> sx==sy+>>>       (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->+>>>          let equalSol (optX, _) (optY, _) = TestLP.approxReal tol optX optY+>>>          in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys+>>>       _ -> False+>>>+>>>+>>> runSuccessive ::+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Foldable t) =>+>>>    CLP.Method ->+>>>    sh ->+>>>    Bounds ix ->+>>>    (Constraints ix, (Direction, Objective sh)) ->+>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->+>>>    Either CLP.FailureType ()+>>> runSuccessive method shape bounds (constrs,dirObj) objs =+>>>    LP.run shape bounds $ ME.runExceptT $ do+>>>       (opt, _xs) <- ME.ExceptT $ LP.simplex method constrs dirObj+>>>       traverse_Lag opt+>>>          (\prevResult (newConstr, dirObjI) -> do+>>>              (optI, _xs) <-+>>>                 ME.ExceptT $+>>>                    LP.simplex method (newConstr prevResult) dirObjI+>>>              return optI)+>>>          objs+>>>+>>> solveSuccessiveWarm ::+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>+>>>    CLP.Method ->+>>>    sh ->+>>>    Bounds ix ->+>>>    (Constraints ix, (Direction, Objective sh)) ->+>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->+>>>    Either CLP.FailureType (NonEmpty.T t (Double, Array sh Double))+>>> solveSuccessiveWarm method shape bounds (constrs,dirObj) objs =+>>>    LP.run shape bounds $ ME.runExceptT $ do+>>>       result <- ME.ExceptT $ LP.simplex method constrs dirObj+>>>       NonEmpty.Cons result <$>+>>>          traverseLag result+>>>             (\(prevOpt, _xs) (newConstr, dirObjI) ->+>>>                 ME.ExceptT $ LP.simplex method (newConstr prevOpt) dirObjI)+>>>             objs+>>>+>>> solveSuccessiveGen ::+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>+>>>    CLP.Method ->+>>>    sh ->+>>>    Bounds ix ->+>>>    (Constraints ix, (Direction, Objective sh)) ->+>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->+>>>    Either CLP.FailureType (NonEmpty.T t (Double, Array sh Double))+>>> solveSuccessiveGen method shape bounds (constrs,dirObj) objs =+>>>    LPMonad.run shape bounds $ ME.runExceptT $ do+>>>       result <-+>>>          ME.ExceptT $ LPMonad.lift (CLP.simplex method) constrs dirObj+>>>       NonEmpty.Cons result <$>+>>>          traverseLag result+>>>             (\(prevOpt, _xs) (newConstr, dirObjI) ->+>>>                 ME.ExceptT $+>>>                    LPMonad.lift (CLP.simplex method)+>>>                       (newConstr prevOpt) dirObjI)+>>>             objs+-}+++newtype T sh a = Cons (MR.ReaderT (sh, Ptr FFI.Simplex) IO 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) =+   unsafePerformIO $ runContT $ do+      lp <- MC.ContT $ bracket FFI.newModel FFI.deleteModel+      liftIO $ Debug.initLog lp+      startPtr <- withBuffer $ Array.vectorFromList [0]+      (collbPtr,colubPtr) <-+         storeBounds $ prepareColumnBoundsArrays shape bounds+      {-+      We would like to force row-major matrix layout,+      but even if we start with a CoinPackedMatrix in row-major layout,+      addRows switches back to column-major layout.+      -}+      liftIO $+         FFI.addColumns lp (fromIntegral $ Shape.size shape)+            collbPtr colubPtr nullPtr+            startPtr nullPtr nullPtr+      liftIO $ MR.runReaderT act (shape, lp)++{- |+Add new constraints to an existing problem+and simplex with a new direction and objective.++>>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.run tripletShape [] (LP.simplex LP.dual [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))+Right (28.0,(5.0,0.0,4.0))++prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (CLP.simplex method bounds constrs (dir,obj), LP.run (Array.shape origin) bounds $ LP.simplex method constrs (dir,obj)) of (Right (optA,_), Right (optB,_)) -> TestLP.approxReal 0.1 optA optB; _ -> False++prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> either (\msg -> QC.counterexample (show msg) False) (const $ QC.property True) $ runSuccessive method (Array.shape origin) bounds (constrs,dirObj) objs++prop> forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> approxSuccession 0.01 (solveSuccessiveWarm method (Array.shape origin) bounds (constrs,dirObj) objs) (solveSuccessiveGen method (Array.shape origin) bounds (constrs,dirObj) objs)+-}+simplex ::+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Method -> Constraints Double ix ->+   (Direction, Objective sh) -> T sh (Result sh)+simplex method constrs (dir,obj) = Cons $ do+   (shape, lp) <- MR.ask+   when (shape /= Array.shape obj) $+      error "COINOR.CLP.Monad.solve: objective shape mismatch"++   liftIO $ runContT $ do+      (_, coefficientsPtr, indexPtr, startPtr) <- storeConstraints shape constrs+      (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs+      objPtr <- withBuffer $ Array.map realToFrac obj+      liftIO $ do+         FFI.addRows lp (fromIntegral $ length constrs)+            rowlbPtr rowubPtr startPtr indexPtr coefficientsPtr+         FFI.chgObjCoefficients lp objPtr++   liftIO $ do+      setOptimizationDirection lp dir+      runMethod method lp+      examineStatus shape lp
+ src/Numeric/COINOR/CLP/Private.hs view
@@ -0,0 +1,154 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.COINOR.CLP.Private where++import qualified Numeric.COINOR.CLP.FFI as FFI+import Numeric.LinearProgramming.Common+         (Term(..), Bound(..), Inequality(Inequality),+          Bounds, Constraints, Direction(..))++import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayMonadic+import qualified Data.Array.Comfort.Storable.Unchecked as ArrayUnchecked+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable (Array)+import Data.Foldable (for_)+import Data.Tuple.HT (mapPair)++import qualified Control.Monad.Trans.Cont as MC+import qualified Control.Applicative.HT as AppHT+import qualified Control.Functor.HT as FuncHT+import Control.Functor.HT (void)++import Foreign.Storable (pokeElemOff, peekElemOff)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.C.Types (CDouble, CInt, CBool)++++withBuffer :: Array sh a -> MC.ContT r IO (Ptr a)+withBuffer arr =+   MC.ContT $ withForeignPtr (ArrayUnchecked.buffer arr)++runContT :: MC.ContT a IO a -> IO a+runContT act = MC.runContT act return++++false, true :: CBool+false = toEnum $ fromEnum False+true  = toEnum $ fromEnum True++positiveInfinity, negativeInfinity :: CDouble+positiveInfinity =  1/0+negativeInfinity = -1/0++prepareBounds :: Inequality a -> (a, (CDouble, CDouble))+prepareBounds (Inequality x bnd) =+   (,) x $+   case bnd of+      LessEqual up    -> (negativeInfinity, realToFrac up)+      GreaterEqual lo -> (realToFrac lo,    positiveInfinity)+      Between lo up   -> (realToFrac lo,    realToFrac up)+      Equal y         -> (realToFrac y,     realToFrac y)+      Free            -> (negativeInfinity, positiveInfinity)++prepareColumnBoundsArrays ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   sh -> Bounds ix -> (Array sh CDouble, Array sh CDouble)+prepareColumnBoundsArrays shape =+   mapPair (Array.fromBoxed, Array.fromBoxed) .+   FuncHT.unzip .+   BoxedArray.fromAssociations (0, positiveInfinity) shape .+   map prepareBounds+++type ShapeInt = Shape.ZeroBased Int++prepareRowBoundsArrays ::+   Bounds ix -> (Array ShapeInt CDouble, Array ShapeInt CDouble)+prepareRowBoundsArrays constrs =+   let shape = Shape.ZeroBased $ length constrs in+   mapPair (Array.fromList shape, Array.fromList shape) $+   unzip $ map (snd . prepareBounds) constrs++storeBounds ::+   (Array sh CDouble, Array sh CDouble) ->+   MC.ContT r IO (Ptr CDouble, Ptr CDouble)+storeBounds = AppHT.mapPair (withBuffer, withBuffer)+++storeConstraints ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   sh -> Constraints Double ix ->+   MC.ContT r IO (ShapeInt, Ptr CDouble, Ptr CInt, Ptr FFI.BigIndex)+storeConstraints shape constrs = do+   let rowStarts =+         Array.vectorFromList $ scanl (+) 0 $+         map (\(Inequality terms _bnd) -> fromIntegral $ length terms)+            constrs+   startPtr <- withBuffer rowStarts+   let shapeOffset = Shape.offset shape+   let coefficients =+         concatMap (\(Inequality terms _bnd) -> terms) constrs+   let coefficientArr =+         Array.vectorFromList $+         map (\(Term _ ix) -> fromIntegral $ shapeOffset ix) coefficients+   indexPtr <- withBuffer coefficientArr+   coefficientsPtr <-+      withBuffer $ Array.vectorFromList $+      map (\(Term c _) -> realToFrac c) coefficients+   return (Array.shape coefficientArr, coefficientsPtr, indexPtr, startPtr)+++setOptimizationDirection :: Ptr FFI.Simplex -> Direction -> IO ()+setOptimizationDirection lp dir =+   FFI.setOptimizationDirection lp $+      case dir of Minimize -> 1; Maximize -> -1+++newtype Method = Method {runMethod :: Ptr FFI.Simplex -> IO ()}++dual, primal :: Method+dual = Method $ \lp -> void $ FFI.dual lp 0+primal = Method $ \lp -> void $ FFI.primal lp 0++initialSolve, initialDualSolve, initialPrimalSolve,+   initialBarrierSolve, initialBarrierNoCrossSolve :: Method+initialSolve = Method $ void . FFI.initialSolve+initialDualSolve = Method $ void . FFI.initialDualSolve+initialPrimalSolve = Method $ void . FFI.initialPrimalSolve+initialBarrierSolve = Method $ void . FFI.initialBarrierSolve+initialBarrierNoCrossSolve = Method $ void . FFI.initialBarrierNoCrossSolve+++data FailureType =+     PrimalInfeasible+   | DualInfeasible+   | StoppedOnIterations+   | StoppedDueToErrors+   deriving (Eq, Show)++type Result sh = Either FailureType (Double, Array sh Double)++examineStatus ::+   (Shape.C sh) =>+   sh -> Ptr FFI.Simplex -> IO (Either FailureType (Double, Array sh Double))+examineStatus shape lp = do+   status <- FFI.status lp+   case status of+      0 -> do+         objVal <- FFI.objectiveValue lp+         optVec <-+            ArrayMonadic.unsafeCreateWithSize shape $ \size arrPtr -> do+               optVecPtr <- FFI.getColSolution lp+               for_ (take size [0..]) $ \k ->+                  pokeElemOff arrPtr k . realToFrac+                     =<< peekElemOff optVecPtr k+         return $ Right (realToFrac objVal, optVec)+      1 -> return $ Left PrimalInfeasible+      2 -> return $ Left DualInfeasible+      3 -> return $ Left StoppedOnIterations+      _ -> return $ Left StoppedDueToErrors
+ src/debug-off/Numeric/COINOR/CLP/Debug.hs view
@@ -0,0 +1,8 @@+module Numeric.COINOR.CLP.Debug where++import qualified Numeric.COINOR.CLP.FFI as FFI+import Foreign.Ptr (Ptr)+++initLog :: Ptr FFI.Simplex -> IO ()+initLog lp = FFI.setLogLevel lp 0
+ src/debug-on/Numeric/COINOR/CLP/Debug.hs view
@@ -0,0 +1,8 @@+module Numeric.COINOR.CLP.Debug where++import qualified Numeric.COINOR.CLP.FFI as FFI+import Foreign.Ptr (Ptr)+++initLog :: Ptr FFI.Simplex -> IO ()+initLog lp = FFI.setLogLevel lp 1
+ test/Main.hs view
@@ -0,0 +1,12 @@+-- Do not edit! Automatically created with doctest-extract.+module Main where++import qualified Test.Numeric.COINOR.CLP+import qualified Test.Numeric.COINOR.CLP.Monad++import qualified Test.DocTest.Driver as DocTest++main :: IO ()+main = DocTest.run $ do+    Test.Numeric.COINOR.CLP.test+    Test.Numeric.COINOR.CLP.Monad.test
+ test/Test/Numeric/COINOR/CLP.hs view
@@ -0,0 +1,100 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/COINOR/CLP.hs+{-# LINE 45 "src/Numeric/COINOR/CLP.hs" #-}++module Test.Numeric.COINOR.CLP where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 46 "src/Numeric/COINOR/CLP.hs" #-}+import     qualified Numeric.COINOR.CLP as LP+import     qualified Numeric.LinearProgramming.Test as TestLP+import     Numeric.COINOR.CLP (PlusMinusOne(..), (.*), (==.), (<=.), (>=<.))++import     qualified Data.Array.Comfort.Storable as Array+import     qualified Data.Array.Comfort.Shape as Shape++import     Data.Tuple.HT (mapSnd)++import     qualified Test.QuickCheck as QC+import     Test.QuickCheck ((===), (.&&.), (.||.))++type     X = Shape.Element+type     PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+type     TripletShape = Shape.NestedTuple Shape.TupleIndex (X,X,X)++pairShape     :: PairShape+pairShape     = Shape.static++tripletShape     :: TripletShape+tripletShape     = Shape.static++approxReal     :: (Ord a, Num a) => a -> a -> a -> Bool+approxReal     tol x y = abs (x-y) <= tol++genMethod     :: QC.Gen (String, LP.Method)+genMethod     = QC.elements $+       ("dual", LP.dual) :+       ("primal", LP.primal) :+       ("initialSolve", LP.initialSolve) :+       ("initialDualSolve", LP.initialDualSolve) :+       ("initialPrimalSolve", LP.initialPrimalSolve) :+       ("initialBarrierSolve", LP.initialBarrierSolve) :+       -- let tests fail+       -- ("initialBarrierNoCrossSolve", LP.initialBarrierNoCrossSolve) :+       []++forAllMethod     ::+       (QC.Testable prop) => (LP.Method -> prop) -> QC.Property+forAllMethod     prop = QC.forAllShow genMethod fst (prop . snd)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.COINOR.CLP:175: "+{-# LINE 175 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 175 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method (QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd Array.toTuple <$> LP.simplex method [] [[1.*pos, (-1::Double).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (absol,(posResult,negResult))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)+ DocTest.printPrefix "Numeric.COINOR.CLP:176: "+{-# LINE 176 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 176 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd Array.toTuple <$> LP.simplex method [] [[1.*pos, (-1::Double).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (absol,(posResult,negResult))) -> QC.counterexample (show(absol,(posResult,negResult))) $ QC.property (approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)+ DocTest.printPrefix "Numeric.COINOR.CLP:178: "+{-# LINE 178 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 178 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right _ -> True; _ -> False)+ DocTest.printPrefix "Numeric.COINOR.CLP:179: "+{-# LINE 179 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 179 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right (_,sol) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False)+ DocTest.printPrefix "Numeric.COINOR.CLP:180: "+{-# LINE 180 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 180 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Right (_,sol) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.1 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False)+ DocTest.printPrefix "Numeric.COINOR.CLP:181: "+{-# LINE 181 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.property+{-# LINE 181 "src/Numeric/COINOR/CLP.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex method bounds constrs (dir,obj) of Left _ -> QC.property False; Right (opt,sol) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val)+ DocTest.printPrefix "Numeric.COINOR.CLP:168: "+{-# LINE 168 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.example+{-# LINE 168 "src/Numeric/COINOR/CLP.hs" #-}+   (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.dual [] [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))+  [ExpectedLine [LineChunk "Right (28.0,(5.0,0.0,4.0))"]]+ DocTest.printPrefix "Numeric.COINOR.CLP:170: "+{-# LINE 170 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.example+{-# LINE 170 "src/Numeric/COINOR/CLP.hs" #-}+   (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.primal [y >=<. (-12,12)] [[1.*x, (-1).*y] <=. 10, [(-1).*y, (1::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))+  [ExpectedLine [LineChunk "Right (116.0,(22.0,12.0,32.0))"]]+ DocTest.printPrefix "Numeric.COINOR.CLP:172: "+{-# LINE 172 "src/Numeric/COINOR/CLP.hs" #-}+ DocTest.example+{-# LINE 172 "src/Numeric/COINOR/CLP.hs" #-}+   (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.simplex LP.primal [y >=<. (-12,12)] [[PlusOne .* x, MinusOne .* y] <=. 10, [MinusOne .* y, PlusOne .* z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))+  [ExpectedLine [LineChunk "Right (116.0,(22.0,12.0,32.0))"]]
+ test/Test/Numeric/COINOR/CLP/Monad.hs view
@@ -0,0 +1,137 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/COINOR/CLP/Monad.hs+{-# LINE 42 "src/Numeric/COINOR/CLP/Monad.hs" #-}++{-# OPTIONS_GHC -XTypeFamilies #-}+{-# OPTIONS_GHC -XTypeOperators #-}+module Test.Numeric.COINOR.CLP.Monad where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 45 "src/Numeric/COINOR/CLP/Monad.hs" #-}+import     qualified Numeric.COINOR.CLP.Monad as LP+import     qualified Numeric.COINOR.CLP as CLP+import     Test.Numeric.COINOR.CLP.Utility (traverse_Lag, traverseLag)+import     Test.Numeric.COINOR.CLP (TripletShape, tripletShape, forAllMethod)+import     Numeric.COINOR.CLP (Direction, (.*), (<=.))++import     qualified Numeric.LinearProgramming.Monad as LPMonad+import     qualified Numeric.LinearProgramming.Test as TestLP+import     Numeric.LinearProgramming.Common (Bounds, Objective)++import     qualified Data.Array.Comfort.Storable as Array+import     qualified Data.Array.Comfort.Shape as Shape+import     qualified Data.NonEmpty as NonEmpty+import     Data.Array.Comfort.Storable (Array)+import     Data.Traversable (Traversable)+import     Data.Foldable (Foldable)++import     qualified Control.Monad.Trans.Except as ME++import     qualified Data.List.HT as ListHT+import     Data.Tuple.HT (mapSnd)++import     Foreign.Storable (Storable)++import     qualified Test.QuickCheck as QC+++type     Constraints ix = CLP.Constraints Double ix+++approxSuccession     ::+       (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>+       a ->+       Either CLP.FailureType (NonEmpty.T [] (a, Array sh a)) ->+       Either CLP.FailureType (NonEmpty.T [] (a, Array sh a)) ->+       QC.Property+approxSuccession     tol x y =+       QC.counterexample (show x) $+       QC.counterexample (show y) $+       case (x,y) of+          (Left sx, Left sy) -> sx==sy+          (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->+             let equalSol (optX, _) (optY, _) = TestLP.approxReal tol optX optY+             in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys+          _ -> False+++runSuccessive     ::+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Foldable t) =>+       CLP.Method ->+       sh ->+       Bounds ix ->+       (Constraints ix, (Direction, Objective sh)) ->+       t (Double -> Constraints ix, (Direction, Objective sh)) ->+       Either CLP.FailureType ()+runSuccessive     method shape bounds (constrs,dirObj) objs =+       LP.run shape bounds $ ME.runExceptT $ do+          (opt, _xs) <- ME.ExceptT $ LP.simplex method constrs dirObj+          traverse_Lag opt+             (\prevResult (newConstr, dirObjI) -> do+                 (optI, _xs) <-+                    ME.ExceptT $+                       LP.simplex method (newConstr prevResult) dirObjI+                 return optI)+             objs++solveSuccessiveWarm     ::+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>+       CLP.Method ->+       sh ->+       Bounds ix ->+       (Constraints ix, (Direction, Objective sh)) ->+       t (Double -> Constraints ix, (Direction, Objective sh)) ->+       Either CLP.FailureType (NonEmpty.T t (Double, Array sh Double))+solveSuccessiveWarm     method shape bounds (constrs,dirObj) objs =+       LP.run shape bounds $ ME.runExceptT $ do+          result <- ME.ExceptT $ LP.simplex method constrs dirObj+          NonEmpty.Cons result <$>+             traverseLag result+                (\(prevOpt, _xs) (newConstr, dirObjI) ->+                    ME.ExceptT $ LP.simplex method (newConstr prevOpt) dirObjI)+                objs++solveSuccessiveGen     ::+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>+       CLP.Method ->+       sh ->+       Bounds ix ->+       (Constraints ix, (Direction, Objective sh)) ->+       t (Double -> Constraints ix, (Direction, Objective sh)) ->+       Either CLP.FailureType (NonEmpty.T t (Double, Array sh Double))+solveSuccessiveGen     method shape bounds (constrs,dirObj) objs =+       LPMonad.run shape bounds $ ME.runExceptT $ do+          result <-+             ME.ExceptT $ LPMonad.lift (CLP.simplex method) constrs dirObj+          NonEmpty.Cons result <$>+             traverseLag result+                (\(prevOpt, _xs) (newConstr, dirObjI) ->+                    ME.ExceptT $+                       LPMonad.lift (CLP.simplex method)+                          (newConstr prevOpt) dirObjI)+                objs++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:181: "+{-# LINE 181 "src/Numeric/COINOR/CLP/Monad.hs" #-}+ DocTest.property+{-# LINE 181 "src/Numeric/COINOR/CLP/Monad.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (CLP.simplex method bounds constrs (dir,obj), LP.run (Array.shape origin) bounds $ LP.simplex method constrs (dir,obj)) of (Right (optA,_), Right (optB,_)) -> TestLP.approxReal 0.1 optA optB; _ -> False)+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:183: "+{-# LINE 183 "src/Numeric/COINOR/CLP/Monad.hs" #-}+ DocTest.property+{-# LINE 183 "src/Numeric/COINOR/CLP/Monad.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> either (\msg -> QC.counterexample (show msg) False) (const $ QC.property True) $ runSuccessive method (Array.shape origin) bounds (constrs,dirObj) objs)+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:185: "+{-# LINE 185 "src/Numeric/COINOR/CLP/Monad.hs" #-}+ DocTest.property+{-# LINE 185 "src/Numeric/COINOR/CLP/Monad.hs" #-}+     (forAllMethod $ \method -> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> approxSuccession 0.01 (solveSuccessiveWarm method (Array.shape origin) bounds (constrs,dirObj) objs) (solveSuccessiveGen method (Array.shape origin) bounds (constrs,dirObj) objs))+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:178: "+{-# LINE 178 "src/Numeric/COINOR/CLP/Monad.hs" #-}+ DocTest.example+{-# LINE 178 "src/Numeric/COINOR/CLP/Monad.hs" #-}+   (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd Array.toTuple <$> LP.run tripletShape [] (LP.simplex LP.dual [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)))+  [ExpectedLine [LineChunk "Right (28.0,(5.0,0.0,4.0))"]]
+ test/Test/Numeric/COINOR/CLP/Utility.hs view
@@ -0,0 +1,23 @@+module Test.Numeric.COINOR.CLP.Utility where++import Data.Traversable (Traversable, traverse)+import Data.Foldable (Foldable, traverse_)++import qualified Control.Monad.Trans.State as MS++import Data.Tuple.HT (double)+++traverse_Lag ::+   (Foldable t, Monad m) =>+   b -> (b -> a -> m b) -> t a -> m ()+traverse_Lag b0 f =+   flip MS.evalStateT b0 .+   traverse_ (\a -> MS.StateT $ \b -> fmap double $ f b a)++traverseLag ::+   (Traversable t, Monad m) =>+   b -> (b -> a -> m b) -> t a -> m (t b)+traverseLag b0 f =+   flip MS.evalStateT b0 .+   traverse (\a -> MS.StateT $ \b -> fmap double $ f b a)