diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -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 diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/cbits/support.cpp b/cbits/support.cpp new file mode 100644 --- /dev/null +++ b/cbits/support.cpp @@ -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 (); +} diff --git a/coinor-clp.cabal b/coinor-clp.cabal new file mode 100644 --- /dev/null +++ b/coinor-clp.cabal @@ -0,0 +1,107 @@ +Cabal-Version: 2.2 +Name: coinor-clp +Version: 0.0 +License: BSD-3-Clause +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +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 diff --git a/include/support.h b/include/support.h new file mode 100644 --- /dev/null +++ b/include/support.h @@ -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 */ diff --git a/src/Numeric/COINOR/CLP.hs b/src/Numeric/COINOR/CLP.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/COINOR/CLP.hs @@ -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 diff --git a/src/Numeric/COINOR/CLP/FFI.hsc b/src/Numeric/COINOR/CLP/FFI.hsc new file mode 100644 --- /dev/null +++ b/src/Numeric/COINOR/CLP/FFI.hsc @@ -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 () diff --git a/src/Numeric/COINOR/CLP/Monad.hs b/src/Numeric/COINOR/CLP/Monad.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/COINOR/CLP/Monad.hs @@ -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 diff --git a/src/Numeric/COINOR/CLP/Private.hs b/src/Numeric/COINOR/CLP/Private.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/COINOR/CLP/Private.hs @@ -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 diff --git a/src/debug-off/Numeric/COINOR/CLP/Debug.hs b/src/debug-off/Numeric/COINOR/CLP/Debug.hs new file mode 100644 --- /dev/null +++ b/src/debug-off/Numeric/COINOR/CLP/Debug.hs @@ -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 diff --git a/src/debug-on/Numeric/COINOR/CLP/Debug.hs b/src/debug-on/Numeric/COINOR/CLP/Debug.hs new file mode 100644 --- /dev/null +++ b/src/debug-on/Numeric/COINOR/CLP/Debug.hs @@ -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 diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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 diff --git a/test/Test/Numeric/COINOR/CLP.hs b/test/Test/Numeric/COINOR/CLP.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/COINOR/CLP.hs @@ -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))"]] diff --git a/test/Test/Numeric/COINOR/CLP/Monad.hs b/test/Test/Numeric/COINOR/CLP/Monad.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/COINOR/CLP/Monad.hs @@ -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))"]] diff --git a/test/Test/Numeric/COINOR/CLP/Utility.hs b/test/Test/Numeric/COINOR/CLP/Utility.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/COINOR/CLP/Utility.hs @@ -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)