diff --git a/.github/workflows/haskell.yml b/.github/workflows/haskell.yml index 7de7f26..92d3748 100644 --- a/.github/workflows/haskell.yml +++ b/.github/workflows/haskell.yml @@ -17,21 +17,17 @@ jobs: steps: - uses: actions/checkout@v3 - - uses: haskell-actions/run-fourmolu@v7 + - uses: haskell-actions/run-fourmolu@v9 + with: + version: "0.14.0.0" build: name: GHC ${{ matrix.ghc-version }} on ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: - os: [ubuntu-latest] - ghc-version: ['9.6', '9.4', '9.2', '9.0', '8.10'] - - include: - - os: windows-latest - ghc-version: '9.6' - - os: macos-latest - ghc-version: '9.6' + os: [windows-latest, macos-latest, ubuntu-latest] + ghc-version: ['9.6', '9.4', '9.2', '9.0'] steps: - uses: actions/checkout@v3 @@ -41,54 +37,94 @@ jobs: id: setup with: ghc-version: ${{ matrix.ghc-version }} - # Defaults, added for clarity: - cabal-version: 'latest' - cabal-update: true + enable-stack: true - - name: Installed minor versions of GHC and Cabal + - name: Installed minor versions of GHC, Cabal, and Stack shell: bash run: | GHC_VERSION=$(ghc --numeric-version) CABAL_VERSION=$(cabal --numeric-version) + STACK_VERSION=$(stack --numeric-version) echo "GHC_VERSION=${GHC_VERSION}" >> "${GITHUB_ENV}" echo "CABAL_VERSION=${CABAL_VERSION}" >> "${GITHUB_ENV}" + echo "STACK_VERSION=${STACK_VERSION}" >> "${GITHUB_ENV}" - name: Configure the build run: | - cabal configure --enable-tests --enable-benchmarks --disable-documentation - cabal build --dry-run + # cabal configure --enable-tests --enable-benchmarks --disable-documentation + # cabal build --dry-run + stack build --test --bench --no-haddock --dry-run # The last step generates dist-newstyle/cache/plan.json for the cache key. - - name: Restore cached dependencies + - name: Restore .stack-work cache + uses: actions/cache/restore@v3 + id: cache-restore-stack-work + with: + path: .stack-work + key: ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-work-${{ hashFiles('stack.yaml') }}-${{ hashFiles('package.yaml') }}-${{ hashFiles('**/*.hs') }} + restore-keys: | + ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-work- + + - name: Restore ~/.stack cache (Unix) + uses: actions/cache/restore@v3 + id: cache-restore-stack-global-unix + if: runner.os == 'Linux' || runner.os == 'macOS' + with: + path: ~/.stack + key: ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global-${{ hashFiles('stack.yaml') }}-${{ hashFiles('package.yaml') }} + restore-keys: | + ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global- + + - name: Restore %APPDATA%\stack, %LOCALAPPDATA%\Programs\stack cache (Windows) uses: actions/cache/restore@v3 - id: cache + id: cache-restore-stack-global-windows + if: runner.os == 'Windows' with: - path: ${{ steps.setup.outputs.cabal-store }} - key: ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-cabal-${{ env.CABAL_VERSION }}-plan-${{ hashFiles('**/plan.json') }} + path: | + ~\AppData\Roaming\stack + ~\AppData\Local\Programs\stack + key: ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global-${{ hashFiles('stack.yaml') }}-${{ hashFiles('package.yaml') }} restore-keys: | - ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-cabal-${{ env.CABAL_VERSION }}- + ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global- + + - name: Build dependencies + run: stack build --only-dependencies - - name: Install dependencies - run: cabal build all --only-dependencies + - name: Build the package + run: stack build - # Cache dependencies already here, so that we do not have to rebuild them should the subsequent steps fail. - - name: Save cached dependencies + - name: Save .stack-work cache uses: actions/cache/save@v3 - # Caches are immutable, trying to save with the same key would error. - if: ${{ !steps.cache.outputs.cache-hit - || steps.cache.outputs.cache-primary-key != steps.cache.outputs.cache-matched-key }} + id: cache-save-stack-work + if: steps.cache-restore-stack-work.outputs.cache-hit != 'true' with: - path: ${{ steps.setup.outputs.cabal-store }} - key: ${{ steps.cache.outputs.cache-primary-key }} - - - name: Build - run: cabal build all + path: .stack-work + key: ${{ steps.cache-restore-stack-work.outputs.cache-primary-key }} + + - name: Save %APPDATA%\stack, %LOCALAPPDATA%\Programs\stack cache (Windows) + uses: actions/cache/save@v3 + if: runner.os == 'Windows' + && steps.cache-restore-stack-global-windows.outputs.cache-hit != 'true' + with: + path: | + ~\AppData\Roaming\stack + ~\AppData\Local\Programs\stack + key: ${{ steps.cache-restore-stack-global-windows.outputs.cache-primary-key }} + + - name: Save ~/.stack cache (Unix) + uses: actions/cache/save@v3 + id: cache-save-stack-global + if: (runner.os == 'Linux' || runner.os == 'macOS') + && steps.cache-restore-stack-global-unix.outputs.cache-hit != 'true' + with: + path: ~/.stack + key: ${{ steps.cache-restore-stack-global-unix.outputs.cache-primary-key }} - name: Run tests - run: cabal test all + run: stack test - name: Check cabal file run: cabal check - name: Build documentation - run: cabal haddock all \ No newline at end of file + run: stack haddock \ No newline at end of file diff --git a/ChangeLog.md b/ChangeLog.md index 341ac1f..5e6fb45 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,9 +2,18 @@ ## Unreleased changes +## [v0.2.0.0](https://github.com/rasheedja/LPPaver/tree/v0.2.0.0) + - Setup CI - Use fourmolu formatter -- Switch to Cabal +- Add better types +- Use lens +- Use RecordDot syntax +- Add logging +- Improve Docs +- More Tests +- Bump Stackage LTS +- Rename Linear.Simplex.Simplex -> Linear.Simplex.TwoPhase.Simplex ## [v0.1.0.0](https://github.com/rasheedja/LPPaver/tree/v0.1.0.0) diff --git a/LICENSE b/LICENSE index ca254af..f0aec98 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright Junaid Rasheed (c) 2020-2022 +Copyright Junaid Rasheed (c) 2020-2023 All rights reserved. diff --git a/README.md b/README.md index 1bf5244..705c1be 100644 --- a/README.md +++ b/README.md @@ -4,62 +4,62 @@ ## Quick Overview -The `Linear.Simplex.Simplex` module contain both phases of the simplex method. +The `Linear.Simplex.Solver.TwoPhase` module contain both phases of the two-phase simplex method. ### Phase One Phase one is implemented by `findFeasibleSolution`: ```haskell -findFeasibleSolution :: [PolyConstraint] -> Maybe (DictionaryForm, [Integer], [Integer], Integer) +findFeasibleSolution :: [StandardConstraint] -> Maybe (DictionaryForm, [Integer], [Integer], Integer) ``` -`findFeasibleSolution` takes a list of `PolyConstraint`s. -The `PolyConstraint` type, as well as other custom types required by this library, are defined in the `Linear.Simplex.Types` module. -`PolyConstraint` is defined as: +`findFeasibleSolution` takes a list of `StandardConstraint`s. +The `StandardConstraint` type, as well as other custom types required by this library, are defined in the `Linear.Simplex.Types` module. +`StandardConstraint` is defined as: ```haskell -data PolyConstraint = - LEQ VarConstMap Rational | - GEQ VarConstMap Rational | - EQ VarConstMap Rational deriving (Show, Eq); +data StandardConstraint = + LEQ Vars Rational | + GEQ Vars Rational | + EQ Vars Rational deriving (Show, Eq); ``` -And `VarConstMap` is defined as: +And `Vars` is defined as: ```haskell -type VarConstMap = [(Integer, Rational)] +type Vars = [(Integer, Rational)] ``` -A `VarConstMap` is treated as a list of `Integer` variables mapped to their `Rational` coefficients, with an implicit `+` between each element in the list. +A `Vars` is treated as a list of `Integer` variables mapped to their `Rational` coefficients, with an implicit `+` between each element in the list. For example: `[(1, 2), (2, (-3)), (1, 3)]` is equivalent to `(2x1 + (-3x2) + 3x1)`. -And a `PolyConstraint` is an inequality/equality where the LHS is a `VarConstMap` and the RHS is a `Rational`. +And a `StandardConstraint` is an inequality/equality where the LHS is a `Vars` and the RHS is a `Rational`. For example: `LEQ [(1, 2), (2, (-3)), (1, 3)] 60` is equivalent to `(2x1 + (-3x2) + 3x1) <= 60`. -Passing a `[PolyConstraint]` to `findFeasibleSolution` will return a feasible solution if it exists as well as a list of slack variables, artificial variables, and a variable that can be safely used to represent the objective for phase two. -`Nothing` is returned if the given `[PolyConstraint]` is infeasible. +Passing a `[StandardConstraint]` to `findFeasibleSolution` will return a feasible solution if it exists as well as a list of slack variables, artificial variables, and a variable that can be safely used to represent the objective for phase two. +`Nothing` is returned if the given `[StandardConstraint]` is infeasible. The feasible system is returned as the type `DictionaryForm`: ```haskell -type DictionaryForm = [(Integer, VarConstMap)] +type DictionaryForm = [(Integer, Vars)] ``` -`DictionaryForm` can be thought of as a list of equations, where the `Integer` represents a basic variable on the LHS that is equal to the RHS represented as a `VarConstMap`. In this `VarConstMap`, the `Integer` -1 is used internally to represent a `Rational` number. +`DictionaryForm` can be thought of as a list of equations, where the `Integer` represents a basic variable on the LHS that is equal to the RHS represented as a `Vars`. In this `Vars`, the `Integer` -1 is used internally to represent a `Rational` number. ### Phase Two `optimizeFeasibleSystem` performs phase two of the simplex method, and has the type: ```haskell -data ObjectiveFunction = Max VarConstMap | Min VarConstMap deriving (Show, Eq) +data ObjectiveFunction = Max Vars | Min Vars deriving (Show, Eq) optimizeFeasibleSystem :: ObjectiveFunction -> DictionaryForm -> [Integer] -> [Integer] -> Integer -> Maybe (Integer, [(Integer, Rational)]) ``` We first pass an `ObjectiveFunction`. Then we give a feasible system in `DictionaryForm`, a list of slack variables, a list of artificial variables, and a variable to represent the objective. -`optimizeFeasibleSystem` Maximizes/Minimizes the linear equation represented as a `VarConstMap` in the given `ObjectiveFunction`. +`optimizeFeasibleSystem` Maximizes/Minimizes the linear equation represented as a `Vars` in the given `ObjectiveFunction`. The first item of the returned pair is the `Integer` variable representing the objective. The second item is a list of `Integer` variables mapped to their optimized values. If a variable is not in this list, the variable is equal to 0. @@ -70,7 +70,7 @@ If a variable is not in this list, the variable is equal to 0. It has the type: ```haskell -twoPhaseSimplex :: ObjectiveFunction -> [PolyConstraint] -> Maybe (Integer, [(Integer, Rational)]) +twoPhaseSimplex :: ObjectiveFunction -> [StandardConstraint] -> Maybe (Integer, [(Integer, Rational)]) ``` The return type is the same as that of `optimizeFeasibleSystem` @@ -87,13 +87,13 @@ There are similar functions for `DictionaryForm` as well as other custom types i ## Usage notes -You must only use positive `Integer` variables in a `VarConstMap`. +You must only use positive `Integer` variables in a `Vars`. This implementation assumes that the user only provides positive `Integer` variables; the `Integer` -1, for example, is sometimes used to represent a `Rational` number. ## Example ```haskell -exampleFunction :: (ObjectiveFunction, [PolyConstraint]) +exampleFunction :: (ObjectiveFunction, [StandardConstraint]) exampleFunction = ( Max [(1, 3), (2, 5)], -- 3x1 + 5x2 @@ -122,7 +122,7 @@ Just ``` There are many more examples in test/TestFunctions.hs. -You may use `prettyShowVarConstMap`, `prettyShowPolyConstraint`, and `prettyShowObjectiveFunction` to convert these tests into a more human-readable format. +You may use `prettyShowVarConstMap`, `prettyShowStandardConstraint`, and `prettyShowObjectiveFunction` to convert these tests into a more human-readable format. ## Issues diff --git a/fourmolu.yaml b/fourmolu.yaml index cb7e946..7d6d1f0 100644 --- a/fourmolu.yaml +++ b/fourmolu.yaml @@ -1,16 +1,16 @@ indentation: 2 -column-limit: none +column-limit: 80 function-arrows: trailing comma-style: leading -import-export-style: diff-friendly +import-export-style: leading indent-wheres: true record-brace-space: true newlines-between-decls: 1 -haddock-style: multi-line -haddock-style-module: -let-style: auto -in-style: right-align +haddock-style: single-line +haddock-style-module: single-line +let-style: inline +in-style: left-align single-constraint-parens: always unicode: never -respectful: true +respectful: false fixities: [] diff --git a/package.yaml b/package.yaml new file mode 100644 index 0000000..80225a5 --- /dev/null +++ b/package.yaml @@ -0,0 +1,64 @@ +name: simplex-method +version: 0.2.0.0 +github: "rasheedja/simplex-method" +license: BSD3 +author: "Junaid Rasheed" +maintainer: "jrasheed178@gmail.com" +copyright: "BSD-3" + +extra-source-files: +- README.md +- ChangeLog.md + +# Metadata used when publishing your package +synopsis: Implementation of the two-phase simplex method in exact rational arithmetic +category: Math, Maths, Mathematics, Optimisation, Optimization, Linear Programming + +# To avoid duplicated efforts in documentation and dealing with the +# complications of embedding Haddock markup inside cabal files, it is +# common to point users to the README.md file. +description: Please see the README on GitHub at + +dependencies: +- base >= 4.14 && < 5 +- containers >= 0.6.5.1 && < 0.7 +- generic-lens >= 2.2.0 && < 2.3 +- lens >= 5.2.2 && < 5.3 +- monad-logger >= 0.3.40 && < 0.4 +- text >= 2.0.2 && < 2.1 +- time +- hspec +- QuickCheck + +default-extensions: + DataKinds + DeriveFunctor + DeriveGeneric + DerivingStrategies + DisambiguateRecordFields + DuplicateRecordFields + FlexibleContexts + GeneralizedNewtypeDeriving + LambdaCase + OverloadedLabels + OverloadedRecordDot + OverloadedStrings + RecordWildCards + TemplateHaskell + TupleSections + TypeApplications + TypeFamilies + NamedFieldPuns + +library: + source-dirs: src + +tests: + simplex-method-test: + defaults: hspec/hspec@main + main: Spec.hs + source-dirs: test + dependencies: + - hspec + - QuickCheck + - simplex-method diff --git a/simplex-method.cabal b/simplex-method.cabal index 80ceb5b..48b78e1 100644 --- a/simplex-method.cabal +++ b/simplex-method.cabal @@ -1,11 +1,11 @@ -cabal-version: 3.6 +cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.34.4. +-- This file has been generated from package.yaml by hpack version 0.37.0. -- -- see: https://github.com/sol/hpack name: simplex-method -version: 0.1.0.0 +version: 0.2.0.0 synopsis: Implementation of the two-phase simplex method in exact rational arithmetic description: Please see the README on GitHub at category: Math, Maths, Mathematics, Optimisation, Optimization, Linear Programming @@ -14,7 +14,7 @@ bug-reports: https://github.com/rasheedja/simplex-method/issues author: Junaid Rasheed maintainer: jrasheed178@gmail.com copyright: BSD-3 -license: BSD-3-Clause +license: BSD3 license-file: LICENSE build-type: Simple extra-source-files: @@ -25,33 +25,81 @@ source-repository head type: git location: https://github.com/rasheedja/simplex-method -common common-extensions - default-extensions: - LambdaCase - TupleSections - library - import: common-extensions exposed-modules: + Comparison.Types + Linear.CanonicalForm.Solver + Linear.CanonicalForm.Types + Linear.CanonicalForm.Util + Linear.Constraint.Linear.Types + Linear.Constraint.Linear.Util + Linear.Constraint.Simple.Types + Linear.Constraint.Simple.Util + Linear.Constraint.Types + Linear.Constraint.Util + Linear.Expr.Types + Linear.Expr.Util Linear.Simplex.Prettify - Linear.Simplex.Simplex + Linear.Simplex.Solver.TwoPhase + Linear.Simplex.Solver.Types + Linear.Simplex.Standardize Linear.Simplex.Types Linear.Simplex.Util + Linear.System.Linear.Types + Linear.System.Linear.Util + Linear.System.Simple.Types + Linear.System.Simple.Util + Linear.System.Types + Linear.Tableau.Types + Linear.Term.Types + Linear.Term.Util + Linear.Var.Types + Linear.Var.Util + other-modules: + Paths_simplex_method hs-source-dirs: src + default-extensions: + DataKinds DeriveFunctor DeriveGeneric DerivingStrategies DisambiguateRecordFields DuplicateRecordFields FlexibleContexts GeneralizedNewtypeDeriving LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications TypeFamilies NamedFieldPuns build-depends: - base >=4.7 && <5 + QuickCheck + , base >=4.14 && <5 + , containers >=0.6.5.1 && <0.7 + , generic-lens >=2.2.0 && <2.3 + , hspec + , lens >=5.2.2 && <5.3 + , monad-logger >=0.3.40 && <0.4 + , text >=2.0.2 && <2.1 + , time default-language: Haskell2010 -test-suite simplex-haskell-test - import: common-extensions +test-suite simplex-method-test type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: - TestFunctions + Linear.CanonicalForm.UtilSpec + Linear.Constraint.Simple.UtilSpec + Linear.Expr.UtilSpec + Linear.System.Simple.UtilSpec + Linear.Tableau.TypesSpec + Linear.Term.UtilSpec + Linear.Var.UtilSpec + TestUtil + Paths_simplex_method hs-source-dirs: test + default-extensions: + DataKinds DeriveFunctor DeriveGeneric DerivingStrategies DisambiguateRecordFields DuplicateRecordFields FlexibleContexts GeneralizedNewtypeDeriving LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications TypeFamilies NamedFieldPuns build-depends: - base >=4.7 && <5 + QuickCheck + , base >=4.14 && <5 + , containers >=0.6.5.1 && <0.7 + , generic-lens >=2.2.0 && <2.3 + , hspec + , lens >=5.2.2 && <5.3 + , monad-logger >=0.3.40 && <0.4 , simplex-method + , text >=2.0.2 && <2.1 + , time default-language: Haskell2010 + build-tool-depends: hspec-discover:hspec-discover == 2.* diff --git a/src/Comparison/Types.hs b/src/Comparison/Types.hs new file mode 100644 index 0000000..1d860c5 --- /dev/null +++ b/src/Comparison/Types.hs @@ -0,0 +1,72 @@ +-- | +-- Module : Comparison.Types +-- Description : Types for constraints in linear programming problems +-- Copyright : (c) Junaid Rasheed, 2020-2024 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +module Comparison.Types + ( MixedComparison (..) + , getLHS + , getRHS + ) +where + +import Control.Applicative (liftA2) +import Foreign.C.Types (CBool) +import GHC.Generics (Generic) +import Test.QuickCheck (Arbitrary, arbitrary, genericShrink, oneof) + +data MixedComparison a b = a :<= b | a :>= b | a :== b + deriving (Show, Read, Eq, Generic) + +instance (Arbitrary a, Arbitrary b) => Arbitrary (MixedComparison a b) where + arbitrary = + oneof + [ liftA2 (:<=) arbitrary arbitrary + , liftA2 (:>=) arbitrary arbitrary + , liftA2 (:==) arbitrary arbitrary + ] + +getLHS :: MixedComparison a b -> a +getLHS (a :<= _) = a +getLHS (a :>= _) = a +getLHS (a :== _) = a + +getRHS :: MixedComparison a b -> b +getRHS (_ :<= b) = b +getRHS (_ :>= b) = b +getRHS (_ :== b) = b + +{- Using a class here and staying 'generic' (as in, be permissive on allowed +types) is awkward. I think it's simpler to just stick with the data type. +If we want a class, how do we best define the comparison ops? We'd need a way +for LhsType and RhsType to be compared. Maybe we can use something like +MixedTypesNum, but that's going to take some work. +-} +class MixedComparison2 c where + type LhsType c :: * + type RhsType c :: * + + lhs :: c -> LhsType c + rhs :: c -> RhsType c + + (.<=) :: c -> Bool + (.>=) :: c -> Bool + + (.==) :: c -> Bool + (.==) c = (.>=) c && (.<=) c + +data IntComparison = IntComparison Int Int + +instance MixedComparison2 IntComparison where + type LhsType IntComparison = Int + type RhsType IntComparison = Int + + lhs (IntComparison l _) = l + rhs (IntComparison _ r) = r + + (.<=) (IntComparison l r) = l <= r + (.>=) (IntComparison l r) = l >= r + + (.==) (IntComparison l r) = l == r diff --git a/src/Linear/CanonicalForm/Solver.hs b/src/Linear/CanonicalForm/Solver.hs new file mode 100644 index 0000000..17bb703 --- /dev/null +++ b/src/Linear/CanonicalForm/Solver.hs @@ -0,0 +1,9 @@ +module Linear.CanonicalForm.Solver where + +import Prelude + +import Linear.CanonicalForm.Types (CanonicalForm(..)) + +-- TwoPhase + +-- Revised diff --git a/src/Linear/CanonicalForm/Types.hs b/src/Linear/CanonicalForm/Types.hs new file mode 100644 index 0000000..3ac47df --- /dev/null +++ b/src/Linear/CanonicalForm/Types.hs @@ -0,0 +1,28 @@ +-- | +-- Module: Linear.Simplex.CanonicalForm.Types +-- Description: Types for augmented (slack) form of linear programming problems +-- Copyright: (c) Junaid Rasheed, 2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.CanonicalForm.Types where + +import qualified Data.Map as Map +import qualified Data.Set as Set +import GHC.Generics (Generic) +import Linear.Constraint.Linear.Types (LinearEquation (..)) +import Linear.Expr.Types (Expr, ExprVarsOnly) +import Linear.Expr.Util (exprVarsOnlyVars) +import Linear.System.Linear.Types (LinearSystem (..)) +import Linear.System.Simple.Types +import Linear.Var.Types (SimplexNum, Var) + +-- https://en.wikipedia.org/wiki/Linear_programming#Augmented_form_(slack_form) +data CanonicalForm = CanonicalForm + { constraints :: !LinearSystem + , originalVars :: !(Set.Set Var) + , systemVars :: !(Set.Set Var) + , systemSlackVars :: !(Set.Set Var) -- all vars are non-negative + , eliminatedVarsMap :: !(Map.Map Var Expr) + } + deriving (Show, Eq, Read, Generic) diff --git a/src/Linear/CanonicalForm/Util.hs b/src/Linear/CanonicalForm/Util.hs new file mode 100644 index 0000000..f74fd96 --- /dev/null +++ b/src/Linear/CanonicalForm/Util.hs @@ -0,0 +1,160 @@ +-- | +-- Module: Linear.Simplex.StandardForm +-- Description: Standard form of the linear program +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.CanonicalForm.Util where + +import Comparison.Types + ( MixedComparison ((:<=), (:==), (:>=)) + ) +import qualified Data.Bifunctor as Bifunctor +import qualified Data.Map as Map +import qualified Data.Maybe as Maybe +import qualified Data.Set as Set +import Linear.Constraint.Linear.Types (LinearEquation (..)) +import qualified Linear.Constraint.Linear.Util as CLU +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Constraint.Simple.Util + ( substVarSimpleConstraintExpr + ) +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.Expr.Util (exprVarsOnlyToExpr) +import Linear.CanonicalForm.Types (CanonicalForm (..)) +import Linear.System.Linear.Types (LinearSystem (..)) +import qualified Linear.System.Linear.Util as SLU +import Linear.System.Simple.Types + ( SimpleSystem (..) + , simplifySimpleSystem + ) +import qualified Linear.System.Simple.Types as SST +import Linear.System.Simple.Util (deriveBounds) +import Linear.Term.Types + ( Term (..) + , TermVarsOnly (..) + ) +import Linear.Var.Types (Bounds (..), Var(..), VarBounds) +import qualified Linear.Var.Util as LVU + +-- | Eliminate non-zero lower bounds via substitution +-- Return the system with the eliminated variables and a map of the eliminated variables to their equivalent expressions +-- First step here https://en.wikipedia.org/wiki/Simplex_algorithm#Standard_form +eliminateNonZeroLowerBounds :: + SimpleSystem -> Map.Map Var Expr -> (Map.Map Var Expr, SimpleSystem) +eliminateNonZeroLowerBounds constraints eliminatedVarsMap = aux [] constraints.unSimpleSystem + where + -- Eliminate non-zero lower bounds + aux :: + [SimpleConstraint] -> [SimpleConstraint] -> (Map.Map Var Expr, SimpleSystem) + aux _ [] = (eliminatedVarsMap, constraints) + aux checked (c : cs) = case c of + -- x >= 5 + (SimpleConstraint (ExprVarsOnly [VarTermVO var] :>= lowerBound)) -> + if lowerBound == 0 + then aux (checked ++ [c]) cs + else + let newVar = SST.nextAvailableVar constraints + -- y >= 0 + newVarLowerBound = SimpleConstraint $ ExprVarsOnly [VarTermVO newVar] :>= 0 + + -- x = y + 5 + substOldVarWith = Expr (VarTerm newVar : [ConstTerm lowerBound]) + substFn = substVarSimpleConstraintExpr var substOldVarWith + + newConstraints = + simplifySimpleSystem . SimpleSystem $ + map substFn checked ++ newVarLowerBound : map substFn cs + updatedEliminatedVarsMap = Map.insert var substOldVarWith eliminatedVarsMap + in eliminateNonZeroLowerBounds newConstraints updatedEliminatedVarsMap -- TODO: Make more efficient if needed + -- TODO: (do) Deal with == ? + -- (dont) Or remove == from the type + -- and convert to <= and >=? + _ -> aux (checked ++ [c]) cs + +-- Add slack variables... +-- Second step here https://en.wikipedia.org/wiki/Simplex_algorithm#Standard_form +-- Return system of equalities and the slack variables +-- TODO: [Var] should be a set +addSlackVars :: SimpleSystem -> ([Var], LinearSystem) +addSlackVars constraints = + let nextAvailableVar = SST.nextAvailableVar constraints + in aux constraints.unSimpleSystem nextAvailableVar [] + where + aux [] _ slackVars = (slackVars, LinearSystem []) + aux (c : cs) nextVar slackVars = case c of + (SimpleConstraint (ExprVarsOnly exprTs :<= num)) -> + let slackVar = nextVar + newNextVar = LVU.nextVar nextVar + newExpr = ExprVarsOnly $ exprTs ++ [VarTermVO slackVar] + -- slackVarLowerBound = Expr (VarTerm slackVar : []) :>= 0 + (newSlackVars, newConstraints) = aux cs newNextVar slackVars + in ( nextVar : newSlackVars + , SLU.prependLinearEquation (LinearEquation newExpr num) newConstraints + ) + (SimpleConstraint (ExprVarsOnly exprTs :>= num)) -> + let slackVar = nextVar + newNextVar = LVU.nextVar nextVar + newExpr = ExprVarsOnly $ exprTs ++ [CoeffTermVO (-1) slackVar] + -- slackVarLowerBound = Expr (VarTerm slackVar : []) :>= 0 + (newSlackVars, newConstraints) = aux cs newNextVar slackVars + in ( nextVar : newSlackVars + , SLU.prependLinearEquation (LinearEquation newExpr num) newConstraints + ) + (SimpleConstraint (expr :== num)) -> + let (newSlackVars, newConstraints) = aux cs nextVar slackVars + in ( newSlackVars + , SLU.prependLinearEquation (LinearEquation expr num) newConstraints + ) + +-- Eliminate unrestricted variables (lower bound unknown) given some bounds +-- Third step here https://en.wikipedia.org/wiki/Simplex_algorithm#Standard_form +eliminateUnrestrictedLowerBounds :: + LinearSystem -> + VarBounds -> + Map.Map Var Expr -> + (Map.Map Var Expr, LinearSystem) +eliminateUnrestrictedLowerBounds constraints varBoundMap eliminatedVarsMap = aux constraints (Map.toList varBoundMap) + where + aux :: + LinearSystem -> [(Var, Bounds)] -> (Map.Map Var Expr, LinearSystem) + aux _ [] = (eliminatedVarsMap, constraints) + aux cs ((var, Bounds Nothing _) : bounds) = + let highestVar = Maybe.fromMaybe (Var (-1)) $ SLU.findHighestVar constraints + newVarPlus = LVU.nextVar highestVar + newVarMinus = LVU.nextVar newVarPlus + -- newVarPlusLowerBound = Expr (VarTerm newVarPlus : []) :>= 0 + -- newVarMinusLowerBound = Expr (VarTerm newVarMinus : []) :>= 0 + + -- oldVar = newVarPlus - newVarMinus + substOldVarWith = ExprVarsOnly (VarTermVO newVarPlus : [CoeffTermVO (-1) newVarMinus]) + + newConstraints = + LinearSystem $ + map (CLU.substVarWith var substOldVarWith) (unLinearSystem constraints) -- TODO: simplify? + updatedEliminatedVarsMap = Map.insert var (exprVarsOnlyToExpr substOldVarWith) eliminatedVarsMap + in eliminateUnrestrictedLowerBounds + newConstraints + (Map.fromList bounds) + updatedEliminatedVarsMap + aux cs (_ : bounds) = aux cs bounds + +simpleSystemToCanonicalForm :: SimpleSystem -> CanonicalForm +simpleSystemToCanonicalForm system = + CanonicalForm + { constraints = finalSystem + , originalVars = SST.simpleSystemVars system + , systemVars = SLU.linearSystemVars finalSystem + , systemSlackVars = Set.fromList slackVars + , eliminatedVarsMap = eliminatedVarsMap + } + where + (eliminatedNonZeroLowerBoundVarsMap, system1) = eliminateNonZeroLowerBounds system Map.empty + system1Bounds = deriveBounds system1 + (slackVars, linearSystem) = addSlackVars system1 + (eliminatedVarsMap, finalSystem) = + eliminateUnrestrictedLowerBounds + linearSystem + system1Bounds + eliminatedNonZeroLowerBoundVarsMap diff --git a/src/Linear/Constraint/Linear/Types.hs b/src/Linear/Constraint/Linear/Types.hs new file mode 100644 index 0000000..e42e2a4 --- /dev/null +++ b/src/Linear/Constraint/Linear/Types.hs @@ -0,0 +1,20 @@ +-- | +-- Description: Types for linear constraints. +-- Copyright: (c) Junaid Rasheed, 2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.Constraint.Linear.Types where + +import GHC.Generics (Generic) +import Linear.Expr.Types (ExprVarsOnly) +import Linear.Var.Types (SimplexNum) + +-- lhs == rhs +-- TODO: Should I move this? Typicially, this would for a 'LinearConstraint', but I'm renaming 'Constraint' to 'LinearConstraint'. +-- TODO: Maybe I should move 'Constraint' here? +data LinearEquation = LinearEquation + { lhs :: ExprVarsOnly + , rhs :: SimplexNum + } + deriving (Show, Eq, Read, Generic) diff --git a/src/Linear/Constraint/Linear/Util.hs b/src/Linear/Constraint/Linear/Util.hs new file mode 100644 index 0000000..ed86f2d --- /dev/null +++ b/src/Linear/Constraint/Linear/Util.hs @@ -0,0 +1,29 @@ +-- | +-- Module: Linear.Constraint.Linear.Util +-- Description: Utility functions for linear constraints +-- Copyright: (c) Junaid Rasheed, 2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.Constraint.Linear.Util where + +import qualified Data.Set as Set +import Linear.Constraint.Linear.Types (LinearEquation (..)) +import Linear.Expr.Types (ExprVarsOnly) +import Linear.Expr.Util + ( exprVarsOnlyMaxVar + , exprVarsOnlyVars + , substVarExprVarsOnly + ) +import Linear.Var.Types (Var) + +-- | Get the variables in a linear equation +linearEquationVars :: LinearEquation -> Set.Set Var +linearEquationVars (LinearEquation lhs _) = exprVarsOnlyVars lhs + +findHighestVar :: LinearEquation -> Var +findHighestVar (LinearEquation lhs _) = exprVarsOnlyMaxVar lhs + +substVarWith :: + Var -> ExprVarsOnly -> LinearEquation -> LinearEquation +substVarWith var expr (LinearEquation lhs rhs) = LinearEquation (substVarExprVarsOnly var expr lhs) rhs diff --git a/src/Linear/Constraint/Simple/Types.hs b/src/Linear/Constraint/Simple/Types.hs new file mode 100644 index 0000000..7d975fb --- /dev/null +++ b/src/Linear/Constraint/Simple/Types.hs @@ -0,0 +1,24 @@ +-- | +-- Module: Linear.Constraint.Simple.Types +-- Description: Types for simple linear constraints +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Constraint.Simple.Types where + +import Comparison.Types (MixedComparison) +import GHC.Generics (Generic) +import Linear.Expr.Types (ExprVarsOnly) +import Linear.Var.Types (SimplexNum) +import Test.QuickCheck (Arbitrary (..)) + +newtype SimpleConstraint = SimpleConstraint + {unSimpleConstraint :: MixedComparison ExprVarsOnly SimplexNum} + deriving (Show, Eq, Read, Generic) + +instance Arbitrary SimpleConstraint where + arbitrary = SimpleConstraint <$> arbitrary + +class CanBeSimpleConstraint a where + toSimpleConstraint :: a -> SimpleConstraint diff --git a/src/Linear/Constraint/Simple/Util.hs b/src/Linear/Constraint/Simple/Util.hs new file mode 100644 index 0000000..223b462 --- /dev/null +++ b/src/Linear/Constraint/Simple/Util.hs @@ -0,0 +1,103 @@ +-- | +-- Module: Linear.Constraint.Simple.Util +-- Description: Utility functions for simple constraints +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Constraint.Simple.Util where + +import Comparison.Types + ( MixedComparison (..) + ) +import qualified Data.List as L +import qualified Data.Set as Set +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Constraint.Types (Constraint (..)) +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.Expr.Util + ( exprToExprVarsOnly + , exprToList + , exprVars + , exprVarsOnlyToExpr + , listToExpr + , simplifyExpr + , simplifyExprVarsOnly + , substVarExpr + , substVarExprVarsOnly + , subtractExpr + , sumExprConstTerms + , zeroConstExpr + ) +import Linear.Term.Types (Term (..), TermVarsOnly (..)) +import Linear.Var.Types (Var) + +substVarSimpleConstraintExpr :: + Var -> Expr -> SimpleConstraint -> SimpleConstraint +substVarSimpleConstraintExpr var varReplacement (SimpleConstraint (a :<= b)) = + let newExpr = substVarExpr var varReplacement (exprVarsOnlyToExpr a) + newConstraint = newExpr :<= Expr [ConstTerm b] + in constraintToSimpleConstraint $ Constraint newConstraint +substVarSimpleConstraintExpr var varReplacement (SimpleConstraint (a :>= b)) = + let newExpr = substVarExpr var varReplacement (exprVarsOnlyToExpr a) + newConstraint = newExpr :>= Expr [ConstTerm b] + in constraintToSimpleConstraint $ Constraint newConstraint +substVarSimpleConstraintExpr var varReplacement (SimpleConstraint (a :== b)) = + let newExpr = substVarExpr var varReplacement (exprVarsOnlyToExpr a) + newConstraint = newExpr :== Expr [ConstTerm b] + in constraintToSimpleConstraint $ Constraint newConstraint + +substVarSimpleConstraint :: + Var -> ExprVarsOnly -> SimpleConstraint -> SimpleConstraint +substVarSimpleConstraint var varReplacement (SimpleConstraint (a :<= b)) = SimpleConstraint $ substVarExprVarsOnly var varReplacement a :<= b +substVarSimpleConstraint var varReplacement (SimpleConstraint (a :>= b)) = SimpleConstraint $ substVarExprVarsOnly var varReplacement a :>= b +substVarSimpleConstraint var varReplacement (SimpleConstraint (a :== b)) = SimpleConstraint $ substVarExprVarsOnly var varReplacement a :== b + +constraintToSimpleConstraint :: Constraint -> SimpleConstraint +constraintToSimpleConstraint constraint = + case constraint of + Constraint (a :<= b) -> SimpleConstraint $ uncurry (:<=) (calcLhsRhs a b) + Constraint (a :>= b) -> SimpleConstraint $ uncurry (:>=) (calcLhsRhs a b) + Constraint (a :== b) -> SimpleConstraint $ uncurry (:==) (calcLhsRhs a b) + where + calcLhsRhs a b = (lhs, rhs) + where + aConsts = sumExprConstTerms a + bConsts = sumExprConstTerms b + rhs = bConsts - aConsts + + aWithoutConst = simplifyExpr . zeroConstExpr $ a + bWithoutConst = simplifyExpr . zeroConstExpr $ b + + lhs' = subtractExpr aWithoutConst bWithoutConst + lhs = case exprToExprVarsOnly lhs' of + Right exprVarsOnly -> exprVarsOnly + Left err -> + error $ + "constraintToSimpleConstraint: lhs is not ExprVarsOnly. Details: " <> err + +-- | Simplify coeff constraints by dividing the coefficient from both sides +simplifyCoeff :: SimpleConstraint -> SimpleConstraint +simplifyCoeff simpleConstraint@(SimpleConstraint (ExprVarsOnly [CoeffTermVO coeff var] :<= num)) + | coeff == 0 = simpleConstraint + | coeff > 0 = SimpleConstraint $ ExprVarsOnly [VarTermVO var] :<= (num / coeff) + | coeff < 0 = SimpleConstraint $ ExprVarsOnly [VarTermVO var] :>= (num / coeff) +simplifyCoeff simpleConstraint@(SimpleConstraint (ExprVarsOnly [CoeffTermVO coeff var] :>= num)) + | coeff == 0 = simpleConstraint + | coeff > 0 = SimpleConstraint $ ExprVarsOnly [VarTermVO var] :>= (num / coeff) + | coeff < 0 = SimpleConstraint $ ExprVarsOnly [VarTermVO var] :<= (num / coeff) +simplifyCoeff simpleConstraint@(SimpleConstraint (ExprVarsOnly [CoeffTermVO coeff var] :== num)) = + if coeff == 0 + then simpleConstraint + else SimpleConstraint $ ExprVarsOnly [VarTermVO var] :== (num / coeff) +simplifyCoeff simpleConstraint = simpleConstraint + +simplifySimpleConstraint :: SimpleConstraint -> SimpleConstraint +simplifySimpleConstraint (SimpleConstraint (expr :<= num)) = simplifyCoeff . SimpleConstraint $ simplifyExprVarsOnly expr :<= num +simplifySimpleConstraint (SimpleConstraint (expr :>= num)) = simplifyCoeff . SimpleConstraint $ simplifyExprVarsOnly expr :>= num +simplifySimpleConstraint (SimpleConstraint (expr :== num)) = simplifyCoeff . SimpleConstraint $ simplifyExprVarsOnly expr :== num + +simpleConstraintVars :: SimpleConstraint -> Set.Set Var +simpleConstraintVars (SimpleConstraint (expr :<= _)) = exprVars . exprVarsOnlyToExpr $ expr +simpleConstraintVars (SimpleConstraint (expr :>= _)) = exprVars . exprVarsOnlyToExpr $ expr +simpleConstraintVars (SimpleConstraint (expr :== _)) = exprVars . exprVarsOnlyToExpr $ expr diff --git a/src/Linear/Constraint/Types.hs b/src/Linear/Constraint/Types.hs new file mode 100644 index 0000000..e41e14d --- /dev/null +++ b/src/Linear/Constraint/Types.hs @@ -0,0 +1,26 @@ +-- | +-- Module: Linear.Constraint.Types +-- Description: Types for linear constraints +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Constraint.Types where + +import Comparison.Types + ( MixedComparison + , getLHS + , getRHS + ) +import qualified Data.Set as Set +import GHC.Generics (Generic) +import Linear.Expr.Types (Expr) +import Test.QuickCheck (Arbitrary (..)) + +-- Input +-- TODO: Consider LinearConstraint +newtype Constraint = Constraint {unConstraint :: MixedComparison Expr Expr} + deriving (Show, Eq, Read, Generic) + +instance Arbitrary Constraint where + arbitrary = Constraint <$> arbitrary diff --git a/src/Linear/Constraint/Util.hs b/src/Linear/Constraint/Util.hs new file mode 100644 index 0000000..1d50e1e --- /dev/null +++ b/src/Linear/Constraint/Util.hs @@ -0,0 +1,21 @@ +-- | +-- Module: Linear.Constraint.Util +-- Description: Utility functions for constraints +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Constraint.Util where + +import Comparison.Types + ( MixedComparison ((:<=), (:==), (:>=)) + ) +import qualified Data.Set as Set +import Linear.Constraint.Types (Constraint (..)) +import Linear.Expr.Util (exprVars) +import Linear.Var.Types (Var) + +constraintVars :: Constraint -> Set.Set Var +constraintVars (Constraint (lhs :<= rhs)) = exprVars lhs <> exprVars rhs +constraintVars (Constraint (lhs :>= rhs)) = exprVars lhs <> exprVars rhs +constraintVars (Constraint (lhs :== rhs)) = exprVars lhs <> exprVars rhs diff --git a/src/Linear/Expr/Types.hs b/src/Linear/Expr/Types.hs new file mode 100644 index 0000000..ac25fa5 --- /dev/null +++ b/src/Linear/Expr/Types.hs @@ -0,0 +1,36 @@ +-- | +-- Module: Linear.Expr.Types +-- Description: Types for linear expressions +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Expr.Types where + +import GHC.Generics (Generic) +import Linear.Term.Types (Term, TermVarsOnly) +import Test.QuickCheck (Arbitrary (..)) + +-- treat empty expr as 0 +-- Consider a version with a num instance, use + and * operators for the input +newtype Expr = Expr {unExpr :: [Term]} + deriving + ( Show + , Read + , Eq + , Generic + ) + +newtype ExprVarsOnly = ExprVarsOnly {unExprVarsOnly :: [TermVarsOnly]} + deriving + ( Show + , Read + , Eq + , Generic + ) + +instance Arbitrary Expr where + arbitrary = Expr <$> arbitrary + +instance Arbitrary ExprVarsOnly where + arbitrary = ExprVarsOnly <$> arbitrary diff --git a/src/Linear/Expr/Util.hs b/src/Linear/Expr/Util.hs new file mode 100644 index 0000000..76df7cb --- /dev/null +++ b/src/Linear/Expr/Util.hs @@ -0,0 +1,132 @@ +-- | +-- Module: Linear.Expr.Util +-- Description: Utility functions for linear expressions +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.Expr.Util where + +import Data.List.NonEmpty (NonEmpty (..)) +import qualified Data.Maybe as Maybe +import qualified Data.Set as Set +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.Term.Types (Term (..), TermVarsOnly (..)) +import Linear.Term.Util + ( negateTerm + , normalizeTerms + , normalizeTermsVarsOnly + , simplifyTerm + , termVarsOnlyToTerm + , unsafeTermToTermVarsOnly + , zeroConstTerm + ) +import Linear.Var.Types (SimplexNum, Var) + +-- | Convert an 'Expr' to a list of 'Term's. +exprToList :: Expr -> [Term] +exprToList = unExpr + +exprVarsOnlyToList :: ExprVarsOnly -> [TermVarsOnly] +exprVarsOnlyToList = unExprVarsOnly + +-- | Convert a list of 'Term's to an 'Expr'. +listToExpr :: [Term] -> Expr +listToExpr = Expr + +listToExprVarsOnly :: [TermVarsOnly] -> ExprVarsOnly +listToExprVarsOnly = ExprVarsOnly + +exprVars :: Expr -> Set.Set Var +exprVars = Set.fromList . Maybe.mapMaybe termVars . exprToList + where + termVars :: Term -> Maybe Var + termVars (ConstTerm _) = Nothing + termVars (CoeffTerm _ v) = Just v + termVars (VarTerm v) = Just v + +exprVarsOnlyVars :: ExprVarsOnly -> Set.Set Var +exprVarsOnlyVars = exprVars . exprVarsOnlyToExpr + +exprVarsOnlyMaxVar :: ExprVarsOnly -> Var +exprVarsOnlyMaxVar = maximum . exprVarsOnlyVars + +simplifyExpr :: Expr -> Expr +simplifyExpr = listToExpr . normalizeTerms . exprToList + +simplifyExprVarsOnly :: ExprVarsOnly -> ExprVarsOnly +simplifyExprVarsOnly = listToExprVarsOnly . normalizeTermsVarsOnly . exprVarsOnlyToList + +sumExprConstTerms :: Expr -> SimplexNum +sumExprConstTerms (Expr ts) = sumExprConstTerms ts + where + sumExprConstTerms = sum . Maybe.mapMaybe termConst + + termConst :: Term -> Maybe SimplexNum + termConst (ConstTerm c) = Just c + termConst _ = Nothing + +zeroConstExpr :: Expr -> Expr +zeroConstExpr (Expr ts) = Expr $ map zeroConstTerm ts + +negateExpr :: Expr -> Expr +negateExpr (Expr ts) = Expr $ map negateTerm ts + +addExpr :: Expr -> Expr -> Expr +addExpr e1 e2 = + -- Safe as Expr :+ Term is the only constructor + simplifyExpr . listToExpr $ (exprToList e1 <> exprToList e2) + +subtractExpr :: Expr -> Expr -> Expr +subtractExpr e1 e2 = addExpr e1 (negateExpr e2) + +substVarExpr :: Var -> Expr -> Expr -> Expr +substVarExpr var varReplacement = simplifyExpr . listToExpr . aux . exprToList + where + replacementTerms = exprToList varReplacement + + aux :: [Term] -> [Term] + aux [] = [] + aux (t : ts) = case t of + (VarTerm tV) -> if tV == var then aux ts ++ replacementTerms else t : aux ts + (CoeffTerm tC tV) -> + if tV == var + then + let newReplacementTerms = + map + ( simplifyTerm + . \case + (CoeffTerm rC rV) -> CoeffTerm (tC * rC) rV + (VarTerm rV) -> CoeffTerm tC rV + (ConstTerm rC) -> ConstTerm (tC * rC) + ) + replacementTerms + in aux ts ++ newReplacementTerms + else t : aux ts + (ConstTerm _) -> t : aux ts + +substVarExprVarsOnly :: Var -> ExprVarsOnly -> ExprVarsOnly -> ExprVarsOnly +substVarExprVarsOnly var varReplacement expr = + let varReplacement' = exprVarsOnlyToExpr varReplacement + expr' = exprVarsOnlyToExpr expr + result' = substVarExpr var varReplacement' expr' + in unsafeExprToExprVarsOnly result' + +unsafeExprToExprVarsOnly :: Expr -> ExprVarsOnly +unsafeExprToExprVarsOnly (Expr ts) = ExprVarsOnly (map unsafeTermToTermVarsOnly ts) + +exprToExprVarsOnly :: Expr -> Either String ExprVarsOnly +exprToExprVarsOnly expr@(Expr ts) = do + if any isConstTerm ts + then + if sumExprConstTerms expr == 0 + then Right $ ExprVarsOnly [] + else Left $ "safeExprToExprVarsOnly: Expr contains ConstTerm. Expr: " <> show expr + else Right $ unsafeExprToExprVarsOnly expr + where + isConstTerm :: Term -> Bool + isConstTerm (ConstTerm _) = True + isConstTerm _ = False + +exprVarsOnlyToExpr :: ExprVarsOnly -> Expr +exprVarsOnlyToExpr (ExprVarsOnly ts) = Expr $ map termVarsOnlyToTerm ts diff --git a/src/Linear/Simplex/Prettify.hs b/src/Linear/Simplex/Prettify.hs index c422590..647738a 100644 --- a/src/Linear/Simplex/Prettify.hs +++ b/src/Linear/Simplex/Prettify.hs @@ -1,36 +1,43 @@ -{- | -Module : Linear.Simplex.Prettify -Description : Prettifier for "Linear.Simplex.Types" types -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental - -Converts "Linear.Simplex.Types" types into human-readable 'String's --} +-- | +-- Module : Linear.Simplex.Prettify +-- Description : Prettifier for "Linear.Simplex.Types" types +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +-- +-- Converts "Linear.Simplex.Types" types into human-readable 'String's module Linear.Simplex.Prettify where -import Data.Ratio -import Linear.Simplex.Types as T +import qualified Data.Map as M +import Data.Ratio (denominator, numerator) +import Linear.Simplex.Types + ( ObjectiveFunction (Max, Min) + , VarLitMapSum + ) -- | Convert a 'VarConstMap' into a human-readable 'String' -prettyShowVarConstMap :: VarConstMap -> String -prettyShowVarConstMap [] = "" -prettyShowVarConstMap [(v, c)] = prettyShowRational c ++ " * x" ++ show v ++ "" +prettyShowVarConstMap :: VarLitMapSum -> String +prettyShowVarConstMap = aux . M.toList where - prettyShowRational r = - if r < 0 - then "(" ++ r' ++ ")" - else r' + aux [] = "" + aux ((vName, vCoeff) : vs) = prettyShowRational vCoeff ++ " * " ++ show vName ++ " + " ++ aux vs where - r' = if denominator r == 1 then show (numerator r) else show (numerator r) ++ " / " ++ show (numerator r) -prettyShowVarConstMap ((v, c) : vcs) = prettyShowVarConstMap [(v, c)] ++ " + " ++ prettyShowVarConstMap vcs + prettyShowRational r = + if r < 0 + then "(" ++ r' ++ ")" + else r' + where + r' = + if denominator r == 1 + then show (numerator r) + else show (numerator r) ++ " / " ++ show (numerator r) --- | Convert a 'PolyConstraint' into a human-readable 'String' -prettyShowPolyConstraint :: PolyConstraint -> String -prettyShowPolyConstraint (LEQ vcm r) = prettyShowVarConstMap vcm ++ " <= " ++ show r -prettyShowPolyConstraint (GEQ vcm r) = prettyShowVarConstMap vcm ++ " >= " ++ show r -prettyShowPolyConstraint (T.EQ vcm r) = prettyShowVarConstMap vcm ++ " == " ++ show r +-- | Convert a 'StandardConstraint' into a human-readable 'String' +-- prettyShowStandardConstraint :: StandardConstraint -> String +-- prettyShowStandardConstraint (LEQ vcm r) = prettyShowVarConstMap vcm ++ " <= " ++ show r +-- prettyShowStandardConstraint (GEQ vcm r) = prettyShowVarConstMap vcm ++ " >= " ++ show r +-- prettyShowStandardConstraint (Linear.Simplex.Types.EQ vcm r) = prettyShowVarConstMap vcm ++ " == " ++ show r -- | Convert an 'ObjectiveFunction' into a human-readable 'String' prettyShowObjectiveFunction :: ObjectiveFunction -> String diff --git a/src/Linear/Simplex/Simplex.hs b/src/Linear/Simplex/Simplex.hs deleted file mode 100644 index 71b2a3d..0000000 --- a/src/Linear/Simplex/Simplex.hs +++ /dev/null @@ -1,283 +0,0 @@ -{- | -Module : Linear.Simplex.Simplex -Description : Implements the twoPhaseSimplex method -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental - -Module implementing the two-phase simplex method. -'findFeasibleSolution' performs phase one of the two-phase simplex method. -'optimizeFeasibleSystem' performs phase two of the two-phase simplex method. -'twoPhaseSimplex' performs both phases of the two-phase simplex method. --} -module Linear.Simplex.Simplex (findFeasibleSolution, optimizeFeasibleSystem, twoPhaseSimplex) where - -import Data.Bifunctor -import Data.List -import Data.Maybe (fromMaybe, mapMaybe) -import Data.Ratio (denominator, numerator, (%)) -import Linear.Simplex.Types -import Linear.Simplex.Util -import Prelude hiding (EQ) - --- import Debug.Trace (trace) - -trace s a = a - -{- | Find a feasible solution for the given system of 'PolyConstraint's by performing the first phase of the two-phase simplex method - All 'Integer' variables in the 'PolyConstraint' must be positive. - If the system is infeasible, return 'Nothing' - Otherwise, return the feasible system in 'DictionaryForm' as well as a list of slack variables, a list artificial variables, and the objective variable. --} -findFeasibleSolution :: [PolyConstraint] -> Maybe (DictionaryForm, [Integer], [Integer], Integer) -findFeasibleSolution unsimplifiedSystem = - if null artificialVars -- No artificial vars, we have a feasible system - then Just (systemWithBasicVarsAsDictionary, slackVars, artificialVars, objectiveVar) - else case simplexPivot (createObjectiveDict artificialObjective objectiveVar : systemWithBasicVarsAsDictionary) of - Just phase1Dict -> - let eliminateArtificialVarsFromPhase1Tableau = map (second (filter (\(v, _) -> v `notElem` artificialVars))) phase1Dict - in case lookup objectiveVar eliminateArtificialVarsFromPhase1Tableau of - Nothing -> trace "objective row not found in phase 1 tableau" Nothing -- Should this be an error? - Just row -> - if fromMaybe 0 (lookup (-1) row) == 0 - then Just (eliminateArtificialVarsFromPhase1Tableau, slackVars, artificialVars, objectiveVar) - else trace "rhs not zero after phase 1, thus original tableau is infeasible" Nothing - Nothing -> Nothing - where - system = simplifySystem unsimplifiedSystem - - maxVar = - maximum $ - map - ( \case - LEQ vcm _ -> maximum (map fst vcm) - GEQ vcm _ -> maximum (map fst vcm) - EQ vcm _ -> maximum (map fst vcm) - ) - system - - (systemWithSlackVars, slackVars) = systemInStandardForm system maxVar [] - - maxVarWithSlackVars = if null slackVars then maxVar else maximum slackVars - - (systemWithBasicVars, artificialVars) = systemWithArtificialVars systemWithSlackVars maxVarWithSlackVars - - finalMaxVar = if null artificialVars then maxVarWithSlackVars else maximum artificialVars - - systemWithBasicVarsAsDictionary = tableauInDictionaryForm systemWithBasicVars - - artificialObjective = createArtificialObjective systemWithBasicVarsAsDictionary artificialVars - - objectiveVar = finalMaxVar + 1 - - -- Convert a system of 'PolyConstraint's to standard form; a system of only equations ('EQ'). - -- Add slack vars where necessary. - -- This may give you an infeasible system if slack vars are negative when original variables are zero. - -- If a constraint is already EQ, set the basic var to Nothing. - -- Final system is a list of equalities for the given system. - -- To be feasible, all vars must be >= 0. - systemInStandardForm :: [PolyConstraint] -> Integer -> [Integer] -> ([(Maybe Integer, PolyConstraint)], [Integer]) - systemInStandardForm [] _ sVars = ([], sVars) - systemInStandardForm (EQ v r : xs) maxVar sVars = ((Nothing, EQ v r) : newSystem, newSlackVars) - where - (newSystem, newSlackVars) = systemInStandardForm xs maxVar sVars - systemInStandardForm (LEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (v ++ [(newSlackVar, 1)]) r) : newSystem, newSlackVars) - where - newSlackVar = maxVar + 1 - (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) - systemInStandardForm (GEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (v ++ [(newSlackVar, -1)]) r) : newSystem, newSlackVars) - where - newSlackVar = maxVar + 1 - (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) - - -- Add artificial vars to a system of 'PolyConstraint's. - -- Artificial vars are added when: - -- Basic var is Nothing (When the original constraint was already an EQ). - -- Slack var is equal to a negative value (this is infeasible, all vars need to be >= 0). - -- Final system will be a feasible artificial system. - -- We keep track of artificial vars in the second item of the returned pair so they can be eliminated once phase 1 is complete. - -- If an artificial var would normally be negative, we negate the row so we can keep artificial variables equal to 1 - systemWithArtificialVars :: [(Maybe Integer, PolyConstraint)] -> Integer -> (Tableau, [Integer]) - systemWithArtificialVars [] _ = ([], []) - systemWithArtificialVars ((mVar, EQ v r) : pcs) maxVar = - case mVar of - Nothing -> - if r >= 0 - then ((newArtificialVar, (v ++ [(newArtificialVar, 1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) - else ((newArtificialVar, (v ++ [(newArtificialVar, -1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) - Just basicVar -> - case lookup basicVar v of - Just basicVarCoeff -> - if r == 0 - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else - if r > 0 - then - if basicVarCoeff >= 0 -- Should only be 1 in the standard call path - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else ((newArtificialVar, (v ++ [(newArtificialVar, 1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) -- Slack var is negative, r is positive (when original constraint was GEQ) - else -- r < 0 - - if basicVarCoeff <= 0 -- Should only be -1 in the standard call path - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else ((newArtificialVar, (v ++ [(newArtificialVar, -1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) -- Slack var is negative, r is negative (when original constraint was LEQ) - where - newArtificialVar = maxVar + 1 - - (newSystemWithNewMaxVar, artificialVarsWithNewMaxVar) = systemWithArtificialVars pcs newArtificialVar - - (newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) = systemWithArtificialVars pcs maxVar - - -- Create an artificial objective using the given 'Integer' list of artificialVars and the given 'DictionaryForm'. - -- The artificial 'ObjectiveFunction' is the negated sum of all artificial vars. - createArtificialObjective :: DictionaryForm -> [Integer] -> ObjectiveFunction - createArtificialObjective rows artificialVars = Max negatedSumWithoutArtificialVars - where - rowsToAdd = filter (\(i, _) -> i `elem` artificialVars) rows - negatedRows = map (\(_, vcm) -> map (second negate) vcm) rowsToAdd - negatedSum = foldSumVarConstMap ((sort . concat) negatedRows) - negatedSumWithoutArtificialVars = filter (\(v, _) -> v `notElem` artificialVars) negatedSum - -{- | Optimize a feasible system by performing the second phase of the two-phase simplex method. - We first pass an 'ObjectiveFunction'. - Then, the feasible system in 'DictionaryForm' as well as a list of slack variables, a list artificial variables, and the objective variable. - Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction' - and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'. --} -optimizeFeasibleSystem :: ObjectiveFunction -> DictionaryForm -> [Integer] -> [Integer] -> Integer -> Maybe (Integer, [(Integer, Rational)]) -optimizeFeasibleSystem unsimplifiedObjFunction phase1Dict slackVars artificialVars objectiveVar = - if null artificialVars - then displayResults . dictionaryFormToTableau <$> simplexPivot (createObjectiveDict objFunction objectiveVar : phase1Dict) - else displayResults . dictionaryFormToTableau <$> simplexPivot (createObjectiveDict phase2ObjFunction objectiveVar : tail phase1Dict) - where - objFunction = simplifyObjectiveFunction unsimplifiedObjFunction - - displayResults :: Tableau -> (Integer, [(Integer, Rational)]) - displayResults tableau = - ( objectiveVar - , case objFunction of - Max _ -> - map - (second snd) - $ filter (\(basicVar, _) -> basicVar `notElem` slackVars ++ artificialVars) tableau - Min _ -> - map -- We maximized -objVar, so we negate the objVar to get the final value - (\(basicVar, row) -> if basicVar == objectiveVar then (basicVar, negate (snd row)) else (basicVar, snd row)) - $ filter (\(basicVar, _) -> basicVar `notElem` slackVars ++ artificialVars) tableau - ) - - phase2Objective = - (foldSumVarConstMap . sort) $ - concatMap - ( \(var, coeff) -> - case lookup var phase1Dict of - Nothing -> [(var, coeff)] - Just row -> map (second (* coeff)) row - ) - (getObjective objFunction) - - phase2ObjFunction = if isMax objFunction then Max phase2Objective else Min phase2Objective - -{- | Perform the two phase simplex method with a given 'ObjectiveFunction' a system of 'PolyConstraint's. - Assumes the 'ObjectiveFunction' and 'PolyConstraint' is not empty. - Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction' - and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'. --} -twoPhaseSimplex :: ObjectiveFunction -> [PolyConstraint] -> Maybe (Integer, [(Integer, Rational)]) -twoPhaseSimplex objFunction unsimplifiedSystem = - case findFeasibleSolution unsimplifiedSystem of - Just r@(phase1Dict, slackVars, artificialVars, objectiveVar) -> optimizeFeasibleSystem objFunction phase1Dict slackVars artificialVars objectiveVar - Nothing -> Nothing - --- | Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'. -simplexPivot :: DictionaryForm -> Maybe DictionaryForm -simplexPivot dictionary = - trace (show dictionary) $ - case mostPositive (head dictionary) of - Nothing -> - trace - "all neg \n" - trace - (show dictionary) - Just - dictionary - Just pivotNonBasicVar -> - let mPivotBasicVar = ratioTest (tail dictionary) pivotNonBasicVar Nothing Nothing - in case mPivotBasicVar of - Nothing -> trace ("Ratio test failed on non-basic var: " ++ show pivotNonBasicVar ++ "\n" ++ show dictionary) Nothing - Just pivotBasicVar -> - trace - "one pos \n" - trace - (show dictionary) - simplexPivot - (pivot pivotBasicVar pivotNonBasicVar dictionary) - where - ratioTest :: DictionaryForm -> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer - ratioTest [] _ mCurrentMinBasicVar _ = mCurrentMinBasicVar - ratioTest ((basicVar, lp) : xs) mostNegativeVar mCurrentMinBasicVar mCurrentMin = - case lookup mostNegativeVar lp of - Nothing -> ratioTest xs mostNegativeVar mCurrentMinBasicVar mCurrentMin - Just currentCoeff -> - let rhs = fromMaybe 0 (lookup (-1) lp) - in if currentCoeff >= 0 || rhs < 0 - then -- trace (show currentCoeff) - ratioTest xs mostNegativeVar mCurrentMinBasicVar mCurrentMin -- rhs was already in right side in original tableau, so should be above zero - -- Coeff needs to be negative since it has been moved to the RHS - else case mCurrentMin of - Nothing -> ratioTest xs mostNegativeVar (Just basicVar) (Just (rhs / currentCoeff)) - Just currentMin -> - if (rhs / currentCoeff) >= currentMin - then ratioTest xs mostNegativeVar (Just basicVar) (Just (rhs / currentCoeff)) - else ratioTest xs mostNegativeVar mCurrentMinBasicVar mCurrentMin - - mostPositive :: (Integer, VarConstMap) -> Maybe Integer - mostPositive (_, lp) = - case findLargestCoeff lp Nothing of - Just (largestVar, largestCoeff) -> - if largestCoeff <= 0 - then Nothing - else Just largestVar - Nothing -> trace "No variables in first row when looking for most positive" Nothing - where - findLargestCoeff :: VarConstMap -> Maybe (Integer, Rational) -> Maybe (Integer, Rational) - findLargestCoeff [] mCurrentMax = mCurrentMax - findLargestCoeff ((var, coeff) : xs) mCurrentMax = - if var == (-1) - then findLargestCoeff xs mCurrentMax - else case mCurrentMax of - Nothing -> findLargestCoeff xs (Just (var, coeff)) - Just currentMax -> - if snd currentMax >= coeff - then findLargestCoeff xs mCurrentMax - else findLargestCoeff xs (Just (var, coeff)) - - -- Pivot a dictionary using the two given variables. - -- The first variable is the leaving (non-basic) variable. - -- The second variable is the entering (basic) variable. - -- Expects the entering variable to be present in the row containing the leaving variable. - -- Expects each row to have a unique basic variable. - -- Expects each basic variable to not appear on the RHS of any equation. - pivot :: Integer -> Integer -> DictionaryForm -> DictionaryForm - pivot leavingVariable enteringVariable rows = - case lookup enteringVariable basicRow of - Just nonBasicCoeff -> - updatedRows - where - -- Move entering variable to basis, update other variables in row appropriately - pivotEquation = (enteringVariable, map (second (/ negate nonBasicCoeff)) ((leavingVariable, -1) : filter ((enteringVariable /=) . fst) basicRow)) - -- Substitute pivot equation into other rows - updatedRows = - map - ( \(basicVar, vMap) -> - if leavingVariable == basicVar - then pivotEquation - else case lookup enteringVariable vMap of - Just subsCoeff -> (basicVar, (foldSumVarConstMap . sort) (map (second (subsCoeff *)) (snd pivotEquation) ++ filter ((enteringVariable /=) . fst) vMap)) - Nothing -> (basicVar, vMap) - ) - rows - Nothing -> trace "non basic variable not found in basic row" undefined - where - (_, basicRow) = head $ filter ((leavingVariable ==) . fst) rows diff --git a/src/Linear/Simplex/Solver/TwoPhase.hs b/src/Linear/Simplex/Solver/TwoPhase.hs new file mode 100644 index 0000000..14b1ca9 --- /dev/null +++ b/src/Linear/Simplex/Solver/TwoPhase.hs @@ -0,0 +1,680 @@ +-- | +-- Module : Linear.Simplex.Simplex.TwoPhase +-- Description : Implements the twoPhaseSimplex method +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +-- +-- Module implementing the two-phase simplex method. +-- 'findFeasibleSolution' performs phase one of the two-phase simplex method. +-- 'optimizeFeasibleSystem' performs phase two of the two-phase simplex method. +-- 'twoPhaseSimplex' performs both phases of the two-phase simplex method. +module Linear.Simplex.Solver.TwoPhase (findFeasibleSolution, optimizeFeasibleSystem, twoPhaseSimplex) where + +import Control.Lens + ( Traversable (traverse) + , (%~) + , (&) + , (.~) + , (<&>) + ) +import Control.Monad (unless) +import Control.Monad.IO.Class (MonadIO) +import Control.Monad.Logger + ( LogLevel (LevelError, LevelInfo, LevelWarn) + , MonadLogger + ) +import Data.Bifunctor (Bifunctor (second)) +import Data.List (elem, map, maximum, notElem, null, sum, (++)) +import qualified Data.Map as M +import Data.Maybe (fromJust, fromMaybe, mapMaybe) +import Data.Ratio (denominator, numerator, (%)) +import qualified Data.Text as Text +import GHC.Real (Ratio) +import Linear.Simplex.Types + ( Dict + , DictValue (..) + , FeasibleSystem (..) + , ObjectiveFunction (..) + , PivotObjective (..) + , Result (..) + , StandardConstraint + , Tableau + , TableauRow (rhs) + , VarLitMapSum + , Var + ) +import Linear.Simplex.Util + ( combineVarLitMapSums + , dictionaryFormToTableau + , foldVarLitMap + , insertPivotObjectiveToDict + , isMax + , logMsg + , showT + , tableauInDictionaryForm + ) +import Linear.Var.Types (SimplexNum) +import Prelude hiding (EQ) + +-- | Find a feasible solution for the given system of 'StandardConstraint's by performing the first phase of the two-phase simplex method +-- All variables in the 'StandardConstraint' must be positive. +-- If the system is infeasible, return 'Nothing' +-- Otherwise, return the feasible system in 'Dict' as well as a list of slack variables, a list artificial variables, and the objective variable. +findFeasibleSolution :: + (MonadIO m, MonadLogger m) => [StandardConstraint] -> m (Maybe FeasibleSystem) +findFeasibleSolution unsimplifiedSystem = do + logMsg LevelInfo $ + "findFeasibleSolution: Looking for solution for " <> showT unsimplifiedSystem + if null artificialVars -- No artificial vars, we have a feasible system + then do + logMsg + LevelInfo + "findFeasibleSolution: Feasible solution found with no artificial vars" + pure . Just $ + FeasibleSystem + systemWithBasicVarsAsDictionary + slackVars + artificialVars + objectiveVar + else do + logMsg LevelInfo $ + "findFeasibleSolution: Needed to create artificial vars. System with artificial vars (in Tableau form) " + <> showT systemWithBasicVars + mPhase1Dict <- + simplexPivot artificialPivotObjective systemWithBasicVarsAsDictionary + case mPhase1Dict of + Just phase1Dict -> do + logMsg LevelInfo $ + "findFeasibleSolution: System after pivoting with objective" + <> showT artificialPivotObjective + <> ": " + <> showT phase1Dict + let eliminateArtificialVarsFromPhase1Tableau = + M.map + ( \DictValue {..} -> + DictValue + { varMapSum = M.filterWithKey (\k _ -> k `notElem` artificialVars) varMapSum + , .. + } + ) + phase1Dict + case M.lookup objectiveVar eliminateArtificialVarsFromPhase1Tableau of + Nothing -> do + logMsg LevelWarn $ + "findFeasibleSolution: Objective row not found after eliminatiing artificial vars. This is unexpected. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + -- If the objecitve row is not found, the system is feasible iff + -- the artificial vars sum to zero. The value of an artificial + -- variable is 0 if non-basic, and the RHS of the row if basic + let artificialVarsVals = + map + ( \v -> maybe 0 (.constant) (M.lookup v eliminateArtificialVarsFromPhase1Tableau) + ) + artificialVars + let artificialVarsValsSum = sum artificialVarsVals + if artificialVarsValsSum == 0 + then do + logMsg LevelInfo $ + "findFeasibleSolution: Artifical variables sum up to 0, thus original tableau is feasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure . Just $ + FeasibleSystem + { dict = eliminateArtificialVarsFromPhase1Tableau + , slackVars = slackVars + , artificialVars = artificialVars + , objectiveVar = objectiveVar + } + else do + logMsg LevelInfo $ + "findFeasibleSolution: Artifical variables sum up to " + <> showT artificialVarsValsSum + <> ", thus original tableau is infeasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure Nothing + Just row -> + if row.constant == 0 + then do + logMsg LevelInfo $ + "findFeasibleSolution: Objective RHS is zero after pivoting, thus original tableau is feasible. feasible system (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure . Just $ + FeasibleSystem + { dict = eliminateArtificialVarsFromPhase1Tableau + , slackVars = slackVars + , artificialVars = artificialVars + , objectiveVar = objectiveVar + } + else do + unless (row.constant < 0) $ do + let errMsg = + "findFeasibleSolution: Objective RHS is negative after pivoting. This should be impossible. System without artificial vars (in Dict form) " + <> show eliminateArtificialVarsFromPhase1Tableau + logMsg LevelError $ Text.pack errMsg + error errMsg + logMsg LevelInfo $ + "findFeasibleSolution: Objective RHS not zero after phase 1, thus original tableau is infeasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure Nothing + Nothing -> do + logMsg LevelInfo $ + "findFeasibleSolution: Infeasible solution found, could not pivot with objective " + <> showT artificialPivotObjective + <> " over system (in Dict form) " + <> showT systemWithBasicVarsAsDictionary + pure Nothing + where + simplifySystem = undefined + + system = simplifySystem unsimplifiedSystem + + maxVar = undefined + -- maximum $ + -- map + -- ( \case + -- LEQ vcm _ -> maximum (map fst $ M.toList vcm) + -- GEQ vcm _ -> maximum (map fst $ M.toList vcm) + -- EQ vcm _ -> maximum (map fst $ M.toList vcm) + -- ) + -- system + + (systemWithSlackVars, slackVars) = systemInStandardForm system maxVar [] + + maxVarWithSlackVars = if null slackVars then maxVar else maximum slackVars + + (systemWithBasicVars, artificialVars) = systemWithArtificialVars systemWithSlackVars maxVarWithSlackVars + + finalMaxVar = if null artificialVars then maxVarWithSlackVars else maximum artificialVars + + systemWithBasicVarsAsDictionary = tableauInDictionaryForm systemWithBasicVars + + artificialPivotObjective = createArtificialPivotObjective systemWithBasicVarsAsDictionary artificialVars + + objectiveVar = finalMaxVar + 1 + + -- Convert a system of 'StandardConstraint's to standard form; a system of only equations ('EQ'). + -- Add slack vars where necessary. + -- This may give you an infeasible system if slack vars are negative when original variables are zero. + -- If a constraint is already EQ, set the basic var to Nothing. + -- Final system is a list of equalities for the given system. + -- To be feasible, all vars must be >= 0. + systemInStandardForm :: + [StandardConstraint] -> + Var -> + [Var] -> + ([(Maybe Var, StandardConstraint)], [Var]) + systemInStandardForm [] _ sVars = ([], sVars) + -- systemInStandardForm (EQ v r : xs) maxVar sVars = ((Nothing, EQ v r) : newSystem, newSlackVars) + -- where + -- (newSystem, newSlackVars) = systemInStandardForm xs maxVar sVars + -- systemInStandardForm (LEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (M.insert newSlackVar 1 v) r) : newSystem, newSlackVars) + -- where + -- newSlackVar = maxVar + 1 + -- (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) + -- systemInStandardForm (GEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (M.insert newSlackVar (-1) v) r) : newSystem, newSlackVars) + -- where + -- newSlackVar = maxVar + 1 + -- (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) + + -- Add artificial vars to a system of 'StandardConstraint's. + -- Artificial vars are added when: + -- Basic var is Nothing (When the original constraint was already an EQ). + -- Slack var is equal to a negative value (this is infeasible, all vars need to be >= 0). + -- Final system will be a feasible artificial system. + -- We keep track of artificial vars in the second item of the returned pair so they can be eliminated once phase 1 is complete. + -- If an artificial var would normally be negative, we negate the row so we can keep artificial variables equal to 1 + systemWithArtificialVars :: + [(Maybe Var, StandardConstraint)] -> Var -> (Tableau, [Var]) + systemWithArtificialVars [] _ = (M.empty, []) + -- systemWithArtificialVars ((mVar, EQ v r) : pcs) maxVar = + -- case mVar of + -- Nothing -> + -- if r >= 0 + -- then + -- ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar 1 v, rhs = r}) newSystemWithNewMaxVar + -- , newArtificialVar : artificialVarsWithNewMaxVar + -- ) + -- else + -- ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar (-1) v, rhs = r}) newSystemWithNewMaxVar + -- , newArtificialVar : artificialVarsWithNewMaxVar + -- ) + -- Just basicVar -> + -- case M.lookup basicVar v of + -- Just basicVarCoeff -> + -- if r == 0 + -- then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + -- else + -- if r > 0 + -- then + -- if basicVarCoeff >= 0 -- Should only be 1 in the standard call path + -- then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + -- else + -- ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar 1 v, rhs = r}) newSystemWithNewMaxVar + -- , newArtificialVar : artificialVarsWithNewMaxVar -- Slack var is negative, r is positive (when original constraint was GEQ) + -- ) + -- else -- r < 0 + + -- if basicVarCoeff <= 0 -- Should only be -1 in the standard call path + -- then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + -- else + -- ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar (-1) v, rhs = r}) newSystemWithNewMaxVar + -- , newArtificialVar : artificialVarsWithNewMaxVar -- Slack var is negative, r is negative (when original constraint was LEQ) + -- ) + -- Nothing -> error "1" -- undefined + -- where + -- newArtificialVar = maxVar + 1 + + -- (newSystemWithNewMaxVar, artificialVarsWithNewMaxVar) = systemWithArtificialVars pcs newArtificialVar + + -- (newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) = systemWithArtificialVars pcs maxVar + systemWithArtificialVars _ _ = error "systemWithArtificialVars: given system includes non-EQ constraints" + + -- \| Takes a 'Dict' and a '[Var]' as input and returns a 'PivotObjective'. + -- The 'Dict' represents the tableau of a linear program with artificial + -- variables, and '[Var]' represents the artificial variables. + + -- The function first filters out the rows of the tableau that correspond + -- to the artificial variables, and negates them. It then computes the sum + -- of the negated rows, which represents the 'PivotObjective'. + createArtificialPivotObjective :: Dict -> [Var] -> PivotObjective + createArtificialPivotObjective rows artificialVars = + PivotObjective + { variable = objectiveVar + , function = foldVarLitMap $ map (.varMapSum) negatedRowsWithoutArtificialVars + , constant = sum $ map (.constant) negatedRowsWithoutArtificialVars + } + where + -- Filter out non-artificial entries + rowsToAdd = M.filterWithKey (\k _ -> k `elem` artificialVars) rows + negatedRows = + M.map + ( \(DictValue rowVarMapSum rowConstant) -> DictValue (M.map negate rowVarMapSum) (negate rowConstant) + ) + rowsToAdd + -- Negate rows, discard keys and artificial vars since the pivot objective does not care about them + negatedRowsWithoutArtificialVars = + map + ( \(_, DictValue {..}) -> + DictValue + { varMapSum = + M.map negate $ M.filterWithKey (\k _ -> k `notElem` artificialVars) varMapSum + , constant = negate constant + } + ) + $ M.toList rowsToAdd + +-- | Optimize a feasible system by performing the second phase of the two-phase simplex method. +-- We first pass an 'ObjectiveFunction'. +-- Then, the feasible system in 'DictionaryForm' as well as a list of slack variables, a list artificial variables, and the objective variable. +-- Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction' +-- and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'. +optimizeFeasibleSystem :: + (MonadIO m, MonadLogger m) => + ObjectiveFunction -> + FeasibleSystem -> + m (Maybe Result) +optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..}) = do + logMsg LevelInfo $ + "optimizeFeasibleSystem: Optimizing feasible system " + <> showT fsys + <> " with objective " + <> showT objFunction + if null artificialVars + then do + logMsg LevelInfo $ + "optimizeFeasibleSystem: No artificial vars, system is feasible. Pivoting system (in dict form) " + <> showT phase1Dict + <> " with objective " + <> showT normalObjective + fmap (displayResults . dictionaryFormToTableau) + <$> simplexPivot normalObjective phase1Dict + else do + logMsg LevelInfo $ + "optimizeFeasibleSystem: Artificial vars present. Pivoting system (in dict form) " + <> showT phase1Dict + <> " with objective " + <> showT adjustedObjective + fmap (displayResults . dictionaryFormToTableau) + <$> simplexPivot adjustedObjective phase1Dict + where + -- \| displayResults takes a 'Tableau' and returns a 'Result'. The 'Tableau' + -- represents the final tableau of a linear program after the simplex + -- algorithm has been applied. The 'Result' contains the value of the + -- objective variable and a map of the values of all variables appearing + -- in the system, including the objective variable. + -- + -- The function first filters out the rows of the tableau that correspond + -- to the slack and artificial variables. It then extracts the values of + -- the remaining variables and stores them in a map. If the objective + -- function is a maximization problem, the map contains the values of the + -- variables as they appear in the final tableau. If the objective function + -- is a minimization problem, the map contains the values of the variables + -- as they appear in the final tableau, except for the objective variable, + -- which is negated. + displayResults :: Tableau -> Result + displayResults tableau = + Result + { objectiveVar = objectiveVar + , varValMap = extractVarVals + } + where + extractVarVals = + let tableauWithOriginalVars = + M.filterWithKey + ( \basicVarName _ -> + basicVarName `notElem` slackVars ++ artificialVars + ) + tableau + in case objFunction of + Max _ -> + M.map + ( \tableauRow -> + tableauRow.rhs + ) + tableauWithOriginalVars + Min _ -> + M.mapWithKey -- We maximized -objVar, so we negate the objVar to get the final value + ( \basicVarName tableauRow -> + if basicVarName == objectiveVar + then negate $ tableauRow.rhs + else tableauRow.rhs + ) + tableauWithOriginalVars + + -- \| Objective to use when optimising the linear program if no artificial + -- variables were necessary in the first phase. It is essentially the original + -- objective function, with a potential change of sign based on the type of + -- problem (Maximization or Minimization). + normalObjective :: PivotObjective + normalObjective = + PivotObjective + { variable = objectiveVar + , function = + if isMax objFunction + then objFunction.objective + else M.map negate objFunction.objective + , constant = 0 + } + + -- \| Objective to use when optimising the linear program if artificial + -- variables were necessary in the first phase. It is an adjustment to the + -- original objective function, where the linear coefficients are modified + -- by back-substitution of the values of the artificial variables. + adjustedObjective :: PivotObjective + adjustedObjective = + PivotObjective + { variable = objectiveVar + , function = calcVarMap + , constant = calcConstants + } + where + -- \| Compute the adjustment to the constant term of the objective + -- function. It adds up the products of the original coefficients and + -- the corresponding constant term (rhs) of each artificial variable + -- in the phase 1 'Dict'. + calcConstants :: SimplexNum + calcConstants = + sum + $ map + ( \(var, coeff) -> + let multiplyWith = if isMax objFunction then coeff else -coeff + in case M.lookup var phase1Dict of + Nothing -> 0 + Just row -> row.constant * multiplyWith + ) + $ M.toList objFunction.objective + + -- \| Compute the adjustment to the coefficients of the original + -- variables in the objective function. It performs back-substitution + -- of the variables in the original objective function using the + -- current value of each artificial variable in the phase 1 'Dict'. + calcVarMap :: VarLitMapSum + calcVarMap = + foldVarLitMap $ + map + ( M.fromList + . ( \(var, coeff) -> + let multiplyWith = if isMax objFunction then coeff else -coeff + in case M.lookup var phase1Dict of + Nothing -> + [(var, multiplyWith)] + Just row -> map (second (* multiplyWith)) (M.toList $ row.varMapSum) + ) + ) + (M.toList objFunction.objective) + +-- | Perform the two phase simplex method with a given 'ObjectiveFunction' a system of 'StandardConstraint's. +-- Assumes the 'ObjectiveFunction' and 'StandardConstraint' is not empty. +-- Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction' +-- and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'. +twoPhaseSimplex :: + (MonadIO m, MonadLogger m) => + ObjectiveFunction -> + [StandardConstraint] -> + m (Maybe Result) +twoPhaseSimplex objFunction unsimplifiedSystem = do + logMsg LevelInfo $ + "twoPhaseSimplex: Solving system " + <> showT unsimplifiedSystem + <> " with objective " + <> showT objFunction + phase1Result <- findFeasibleSolution unsimplifiedSystem + case phase1Result of + Just feasibleSystem -> do + logMsg LevelInfo $ + "twoPhaseSimplex: Feasible system found for " + <> showT unsimplifiedSystem + <> "; Feasible system: " + <> showT feasibleSystem + optimizedSystem <- optimizeFeasibleSystem objFunction feasibleSystem + logMsg LevelInfo $ + "twoPhaseSimplex: Optimized system found for " + <> showT unsimplifiedSystem + <> "; Optimized system: " + <> showT optimizedSystem + pure optimizedSystem + Nothing -> do + logMsg LevelInfo $ + "twoPhaseSimplex: Phase 1 gives infeasible result for " + <> showT unsimplifiedSystem + pure Nothing + +-- | Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'. +simplexPivot :: + (MonadIO m, MonadLogger m) => PivotObjective -> Dict -> m (Maybe Dict) +simplexPivot + objective@( PivotObjective + { variable = objectiveVar + , function = objectiveFunc + , constant = objectiveConstant + } + ) + dictionary = do + logMsg LevelInfo $ + "simplexPivot: Pivoting with objective " + <> showT objective + <> " over system (in Dict form) " + <> showT dictionary + case mostPositive objectiveFunc of + Nothing -> do + logMsg LevelInfo $ + "simplexPivot: Pivoting complete as no positive variables found in objective " + <> showT objective + <> " over system (in Dict form) " + <> showT dictionary + pure $ Just (insertPivotObjectiveToDict objective dictionary) + Just pivotNonBasicVar -> do + logMsg LevelInfo $ + "simplexPivot: Non-basic pivoting variable in objective, determined by largest coefficient = " + <> showT pivotNonBasicVar + let mPivotBasicVar = ratioTest dictionary pivotNonBasicVar Nothing Nothing + case mPivotBasicVar of + Nothing -> do + logMsg LevelInfo $ + "simplexPivot: Ratio test failed with non-basic variable " + <> showT pivotNonBasicVar + <> " over system (in Dict form) " + <> showT dictionary + pure Nothing + Just pivotBasicVar -> do + logMsg LevelInfo $ + "simplexPivot: Basic pivoting variable determined by ratio test " + <> showT pivotBasicVar + logMsg LevelInfo $ + "simplexPivot: Pivoting with basic var " + <> showT pivotBasicVar + <> ", non-basic var " + <> showT pivotNonBasicVar + <> ", objective " + <> showT objective + <> " over system (in Dict form) " + <> showT dictionary + let pivotResult = + pivot + pivotBasicVar + pivotNonBasicVar + (insertPivotObjectiveToDict objective dictionary) + pivotedObj = + let pivotedObjEntry = + fromMaybe (error "simplexPivot: Can't find objective after pivoting") $ + M.lookup objectiveVar pivotResult + in objective + & #function .~ pivotedObjEntry.varMapSum + & #constant .~ pivotedObjEntry.constant + pivotedDict = M.delete objectiveVar pivotResult + logMsg LevelInfo $ + "simplexPivot: Pivoted, Recursing with new pivoting objective " + <> showT pivotedObj + <> " for new pivoted system (in Dict form) " + <> showT pivotedDict + simplexPivot + pivotedObj + pivotedDict + where + ratioTest :: Dict -> Var -> Maybe Var -> Maybe Rational -> Maybe Var + ratioTest dict = aux (M.toList dict) + where + aux :: [(Var, DictValue)] -> Var -> Maybe Var -> Maybe Rational -> Maybe Var + aux [] _ mCurrentMinBasicVar _ = mCurrentMinBasicVar + aux (x@(basicVar, dictEquation) : xs) mostNegativeVar mCurrentMinBasicVar mCurrentMin = + case M.lookup mostNegativeVar dictEquation.varMapSum of + Nothing -> aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + Just currentCoeff -> + let dictEquationConstant = dictEquation.constant + in if currentCoeff >= 0 || dictEquationConstant < 0 + then aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + else case mCurrentMin of + Nothing -> + aux + xs + mostNegativeVar + (Just basicVar) + (Just (dictEquationConstant / currentCoeff)) + Just currentMin -> + if (dictEquationConstant / currentCoeff) >= currentMin + then + aux + xs + mostNegativeVar + (Just basicVar) + (Just (dictEquationConstant / currentCoeff)) + else aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + + mostPositive :: VarLitMapSum -> Maybe Var + mostPositive varLitMap = + case findLargestCoeff (M.toList varLitMap) Nothing of + Just (largestVarName, largestVarCoeff) -> + if largestVarCoeff <= 0 + then Nothing + else Just largestVarName + Nothing -> Nothing + where + findLargestCoeff :: + [(Var, SimplexNum)] -> Maybe (Var, SimplexNum) -> Maybe (Var, SimplexNum) + findLargestCoeff [] mCurrentMax = mCurrentMax + findLargestCoeff (v@(vName, vCoeff) : vs) mCurrentMax = + case mCurrentMax of + Nothing -> findLargestCoeff vs (Just v) + Just (_, currentMaxCoeff) -> + if currentMaxCoeff >= vCoeff + then findLargestCoeff vs mCurrentMax + else findLargestCoeff vs (Just v) + + -- Pivot a dictionary using the two given variables. + -- The first variable is the leaving (non-basic) variable. + -- The second variable is the entering (basic) variable. + -- Expects the entering variable to be present in the row containing the leaving variable. + -- Expects each row to have a unique basic variable. + -- Expects each basic variable to not appear on the RHS of any equation. + pivot :: Var -> Var -> Dict -> Dict + pivot leavingVariable enteringVariable dict = + case M.lookup enteringVariable (dictEntertingRow.varMapSum) of + Just enteringVariableCoeff -> + updatedRows + where + -- Move entering variable to basis, update other variables in row appropriately + pivotEnteringRow :: DictValue + pivotEnteringRow = + dictEntertingRow + & #varMapSum + %~ ( \basicEquation -> + -- uncurry + M.insert + leavingVariable + (-1) + (filterOutEnteringVarTerm basicEquation) + & traverse + %~ divideByNegatedEnteringVariableCoeff + ) + & #constant + %~ divideByNegatedEnteringVariableCoeff + where + newEnteringVarTerm = (leavingVariable, -1) + divideByNegatedEnteringVariableCoeff = (/ negate enteringVariableCoeff) + + -- Substitute pivot equation into other rows + updatedRows :: Dict + updatedRows = + M.fromList $ map (uncurry f2) $ M.toList dict + where + f entryVar entryVal = + if leavingVariable == entryVar + then pivotEnteringRow + else case M.lookup enteringVariable (entryVal.varMapSum) of + Just subsCoeff -> + entryVal + & #varMapSum + .~ combineVarLitMapSums + (pivotEnteringRow.varMapSum <&> (subsCoeff *)) + (filterOutEnteringVarTerm (entryVal.varMapSum)) + & #constant + .~ ((subsCoeff * (pivotEnteringRow.constant)) + entryVal.constant) + Nothing -> entryVal + + f2 :: Var -> DictValue -> (Var, DictValue) + f2 entryVar entryVal = + if leavingVariable == entryVar + then (enteringVariable, pivotEnteringRow) + else case M.lookup enteringVariable (entryVal.varMapSum) of + Just subsCoeff -> + ( entryVar + , entryVal + & #varMapSum + .~ combineVarLitMapSums + (pivotEnteringRow.varMapSum <&> (subsCoeff *)) + (filterOutEnteringVarTerm (entryVal.varMapSum)) + & #constant + .~ ((subsCoeff * (pivotEnteringRow.constant)) + entryVal.constant) + ) + Nothing -> (entryVar, entryVal) + Nothing -> error "pivot: non basic variable not found in basic row" + where + -- \| The entering row, i.e., the row in the dict which is the value of + -- leavingVariable. + dictEntertingRow = + fromMaybe + (error "pivot: Basic variable not found in Dict") + $ M.lookup leavingVariable dict + + filterOutEnteringVarTerm = M.filterWithKey (\vName _ -> vName /= enteringVariable) diff --git a/src/Linear/Simplex/Solver/Types.hs b/src/Linear/Simplex/Solver/Types.hs new file mode 100644 index 0000000..80d77c6 --- /dev/null +++ b/src/Linear/Simplex/Solver/Types.hs @@ -0,0 +1,81 @@ +module Linear.Simplex.Solver.Types where + +import qualified Data.Map as Map +import GHC.Generics (Generic) +import Linear.Expr.Types (ExprVarsOnly) +import Linear.CanonicalForm.Types (CanonicalForm) +import Linear.System.Linear.Types (LinearSystem) +import Linear.Var.Types (SimplexNum, Var) +import System.Posix.Types (CMode) + +data OptimisationDirection = Minimize | Maximize + deriving (Show, Eq, GHC.Generics.Generic) + +data Objective = Objective + { expr :: Linear.Expr.Types.ExprVarsOnly + , direction :: OptimisationDirection + } + deriving (Show, Eq, GHC.Generics.Generic) + +-- TODO: Is it useful to include the system in the result? +data Result = Result + +-- TODO: Include the canonical form? +data OptimisationResult = OptimisationResult + { varMap :: Map.Map Var SimplexNum + , objVal :: SimplexNum + } + deriving (Show, Read, Eq, GHC.Generics.Generic) + +-- class (CanBeLinearSystem s) => Solver s where +-- solve :: s -> Objective -> Result +class TwoPhaseSolver inputSystem where + firstPhase :: inputSystem -> Maybe CanonicalForm + + twoPhaseSolve :: inputSystem -> Objective -> Maybe OptimisationResult + twoPhaseSolve inputSystem obj = + let mSf = firstPhase inputSystem + in case mSf of + Nothing -> Nothing + Just sf -> Just $ systemResult $ secondPhase obj sf + where + secondPhase :: Objective -> CanonicalForm -> CanonicalForm + secondPhase = undefined + + -- This will probably be a proper function + systemResult :: CanonicalForm -> OptimisationResult + systemResult = undefined + +class CanBeStandardForm problem where + findSolution :: problem -> Maybe CanonicalForm + +-- solveStandardForm :: StandardForm -> Objective -> Maybe Result + +class LinearSystemProcessor s where + type System s :: * + +data FeasibleSystem = FeasibleSystem + { varVals :: Map.Map Var SimplexNum + , system :: LinearSystem + } + +data Model = Model {model :: Map.Map Var SimplexNum} + +data SatResult model = Unsat | Sat model + +-- s is a system +class (Monad (SatSolverMonad s)) => SatSolver s where + type SatSolverOptions s :: * + type SatSolverMonad s :: * -> * + + solve :: SatSolverOptions s -> s -> (SatSolverMonad s) (SatResult Model) + +class (Monad (OptSolverMonad s)) => OptSolver s where + type OptSolverOptions s :: * + type OptSolverMonad s :: * -> * + + optimise :: + OptSolverOptions s -> s -> Objective -> (OptSolverMonad s) (SatResult Model) + +-- class (CanBeLinearSystem s) => Solver2 s where +-- solve2 :: s -> Objective -> Result diff --git a/src/Linear/Simplex/Standardize.hs b/src/Linear/Simplex/Standardize.hs new file mode 100644 index 0000000..d8ca09c --- /dev/null +++ b/src/Linear/Simplex/Standardize.hs @@ -0,0 +1,9 @@ +module Linear.Simplex.Standardize where + +import Data.List (sort) +import qualified Data.Map as M +import GHC.Generics (Generic) + +-- Add slack vars, need type of system with only equalities + +-- Add artificial vars, can we type check this somehow? Maybe with a phantom type? Is Tableau enough? diff --git a/src/Linear/Simplex/Types.hs b/src/Linear/Simplex/Types.hs index 19cb78e..2a3c952 100644 --- a/src/Linear/Simplex/Types.hs +++ b/src/Linear/Simplex/Types.hs @@ -1,52 +1,152 @@ -{- | -Module : Linear.Simplex.Types -Description : Custom types -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental --} +-- | +-- Module : Linear.Simplex.Types +-- Description : Custom types +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental module Linear.Simplex.Types where -{- | List of 'Integer' variables with their 'Rational' coefficients. - There is an implicit addition between elements in this list. - Users must only provide positive integer variables. - - Example: [(2, 3), (6, (-1), (2, 1))] is equivalent to 3x2 + (-x6) + x2. --} -type VarConstMap = [(Integer, Rational)] - -{- | For specifying constraints in a system. - The LHS is a 'VarConstMap', and the RHS, is a 'Rational' number. - LEQ [(1, 2), (2, 1)] 3.5 is equivalent to 2x1 + x2 <= 3.5. - Users must only provide positive integer variables. - - Example: LEQ [(2, 3), (6, (-1), (2, 1))] 12.3 is equivalent to 3x2 + (-x6) + x2 <= 12.3. --} -data PolyConstraint - = LEQ VarConstMap Rational - | GEQ VarConstMap Rational - | EQ VarConstMap Rational - deriving (Show, Eq) - -{- | Create an objective function. - We can either 'Max'imize or 'Min'imize a 'VarConstMap'. --} -data ObjectiveFunction = Max VarConstMap | Min VarConstMap deriving (Show, Eq) - -{- | A 'Tableau' of equations. - Each pair in the list is a row. - The first item in the pair specifies which 'Integer' variable is basic in the equation. - The second item in the pair is an equation. - The 'VarConstMap' in the second equation is a list of variables with their coefficients. - The RHS of the equation is a 'Rational' constant. --} -type Tableau = [(Integer, (VarConstMap, Rational))] - -{- | Type representing equations. - Each pair in the list is one equation. - The first item of the pair is the basic variable, and is on the LHS of the equation with a coefficient of one. - The RHS is represented using a `VarConstMap`. - The integer variable -1 is used to represent a 'Rational' on the RHS --} -type DictionaryForm = [(Integer, VarConstMap)] +import Control.Applicative ((<|>)) +import Data.Generics.Labels () +import qualified Data.List as L +import qualified Data.List.NonEmpty as NE +import qualified Data.Map as M +import qualified Data.Maybe as Maybe +import qualified Data.Set as Set +import qualified Debug.Trace as T +import GHC.Base (liftA2) +import GHC.Generics (Generic) +import Linear.Expr.Types (Expr) +import Linear.Var.Types (SimplexNum) +import Test.QuickCheck (Arbitrary (..), genericShrink, oneof) + +type Var = Int + +data StandardFormRow = StandardFormRow + { lhs :: Expr + , rhs :: SimplexNum + } + deriving (Show, Read, Eq, Generic) + +-- | the system with slack variables +type StandardForm = [StandardFormRow] + +-- data Term = ConstTerm SimplexNum | CoeffTerm SimplexNum Var -- Consider VarTermnewVar +data StandardTerm = StandardTerm SimplexNum Var + deriving (Show, Read, Eq, Generic) + +-- Equality +data StandardConstraint = StandardConstraint [StandardTerm] SimplexNum + deriving (Show, Read, Eq, Generic) + +type StandardSystem = [StandardConstraint] + +-- Negative/no lower bounds +-- say x +-- x = x1 - x2 +-- x1, x2 >= 0 + +-- A 'Tableau' where the basic variable may be empty. +-- All non-empty basic vars are slack vars +data SystemWithSlackVarRow = SystemInStandardFormRow + { mSlackVar :: Maybe Var + -- ^ This is Nothing iff the row does not have a slack variable + , row :: TableauRow + } + +type SystemWithSlackVars = [SystemWithSlackVarRow] + +data FeasibleSystem = FeasibleSystem + { dict :: Dict + , slackVars :: [Var] + , artificialVars :: [Var] + , objectiveVar :: Var + } + deriving (Show, Read, Eq, Generic) + +data Result = Result + { objectiveVar :: Var + , varValMap :: VarLitMap + -- TODO: + -- Maybe VarLitMap + -- , feasible :: Bool + -- , optimisable :: Bool + } + deriving (Show, Read, Eq, Generic) + +data SimplexMeta = SimplexMeta + { objective :: ObjectiveFunction + , feasibleSystem :: Maybe FeasibleSystem + , optimisedResult :: Maybe Result + } + +type VarLitMap = M.Map Var SimplexNum + +-- | List of variables with their 'SimplexNum' coefficients. +-- There is an implicit addition between elements in this list. +-- +-- Example: [Var "x" 3, Var "y" -1, Var "z" 1] is equivalent to 3x + (-y) + z. +type VarLitMapSum = VarLitMap + +-- | For specifying constraints in a system. +-- The LHS is a 'VarLitMapSum', and the RHS, is a 'SimplexNum' number. +-- LEQ [(1, 2), (2, 1)] 3.5 is equivalent to 2x1 + x2 <= 3.5. +-- Users must only provide positive integer variables. +-- +-- Example: LEQ [Var "x" 3, Var "y" -1, Var "x" 1] 12.3 is equivalent to 3x + (-y) + x <= 12.3. + +-- data StandardConstraint +-- = LEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} +-- | GEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} +-- | EQ {lhs :: VarLitMapSum, rhs :: SimplexNum} +-- deriving (Show, Read, Eq, Generic) + +-- | Create an objective function. +-- We can either 'Max'imize or 'Min'imize a 'VarTermSum'. +data ObjectiveFunction + = Max {objective :: VarLitMapSum} + | Min {objective :: VarLitMapSum} + deriving (Show, Read, Eq, Generic) + +-- | TODO: Maybe we want this type +-- TODO: A better/alternative name +data Equation = Equation + { lhs :: VarLitMapSum + , rhs :: SimplexNum + } + +-- | Value for 'Tableau'. lhs = rhs. +data TableauRow = TableauRow + { lhs :: VarLitMapSum + , rhs :: SimplexNum + } + deriving (Show, Read, Eq, Generic) + +-- | A simplex 'Tableu' of equations. +-- Each entry in the map is a row. +type Tableau = M.Map Var TableauRow + +-- | Values for a 'DictEntry'. +data DictValue = DictValue + { varMapSum :: VarLitMapSum + , constant :: SimplexNum + } + deriving (Show, Read, Eq, Generic) + +-- | A simplex 'Dict' +-- One quation represents the objective function. +-- Each pair in the list is one equation in the system we're working with. +-- data Dict = Dict +-- { objective :: DictObjective +-- , entries :: DictEntries +-- } +-- deriving (Show, Read, Eq, Generic) +type Dict = M.Map Var DictValue + +data PivotObjective = PivotObjective + { variable :: Var + , function :: VarLitMapSum + , constant :: SimplexNum + } + deriving (Show, Read, Eq, Generic) diff --git a/src/Linear/Simplex/Util.hs b/src/Linear/Simplex/Util.hs index f661a76..e5a45b8 100644 --- a/src/Linear/Simplex/Util.hs +++ b/src/Linear/Simplex/Util.hs @@ -1,18 +1,45 @@ -{- | -Module : Linear.Simplex.Util -Description : Helper functions -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental - -Helper functions for performing the two-phase simplex method. --} +-- | +-- Module : Linear.Simplex.Util +-- Description : Helper functions +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +-- +-- Helper functions for performing the two-phase simplex method. module Linear.Simplex.Util where -import Data.Bifunctor -import Data.List +import Control.Monad.IO.Class (MonadIO (..)) +import Control.Monad.Logger + ( LogLevel (..) + , LogLine + , MonadLogger + , logDebug + , logError + , logInfo + , logWarn + ) +import Data.Generics.Product (field) +import Data.List (map, nub) +import qualified Data.Map as Map +import qualified Data.Map.Merge.Lazy as MapMerge +import Data.Maybe (fromMaybe) +import qualified Data.Text as T +import Data.Time (getCurrentTime) +import Data.Time.Format.ISO8601 (iso8601Show) import Linear.Simplex.Types + ( Dict + , DictValue (..) + , ObjectiveFunction (Max, Min) + , PivotObjective (constant, function, variable) + , Result (objectiveVar, varValMap) + , Tableau + , TableauRow (TableauRow, lhs, rhs) + , VarLitMap + , VarLitMapSum + , Var + ) +import Linear.Var.Types (SimplexNum) import Prelude hiding (EQ) -- | Is the given 'ObjectiveFunction' to be 'Max'imized? @@ -20,129 +47,176 @@ isMax :: ObjectiveFunction -> Bool isMax (Max _) = True isMax (Min _) = False --- | Extract the objective ('VarConstMap') from an 'ObjectiveFunction' -getObjective :: ObjectiveFunction -> VarConstMap -getObjective (Max o) = o -getObjective (Min o) = o - -{- | Simplifies a system of 'PolyConstraint's by first calling 'simplifyPolyConstraint', - then reducing 'LEQ' and 'GEQ' with same LHS and RHS (and other similar situations) into 'EQ', - and finally removing duplicate elements using 'nub'. --} -simplifySystem :: [PolyConstraint] -> [PolyConstraint] -simplifySystem = nub . reduceSystem . map simplifyPolyConstraint - where - reduceSystem :: [PolyConstraint] -> [PolyConstraint] - reduceSystem [] = [] - -- Reduce LEQ with matching GEQ and EQ into EQ - reduceSystem ((LEQ lhs rhs) : pcs) = - let matchingConstraints = - filter - ( \case - GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in if null matchingConstraints - then LEQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) - -- Reduce GEQ with matching LEQ and EQ into EQ - reduceSystem ((GEQ lhs rhs) : pcs) = - let matchingConstraints = - filter - ( \case - LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in if null matchingConstraints - then GEQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) - -- Reduce EQ with matching LEQ and GEQ into EQ - reduceSystem ((EQ lhs rhs) : pcs) = - let matchingConstraints = - filter - ( \case - LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in if null matchingConstraints - then EQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) - --- | Simplify an 'ObjectiveFunction' by first 'sort'ing and then calling 'foldSumVarConstMap' on the 'VarConstMap'. -simplifyObjectiveFunction :: ObjectiveFunction -> ObjectiveFunction -simplifyObjectiveFunction (Max varConstMap) = Max (foldSumVarConstMap (sort varConstMap)) -simplifyObjectiveFunction (Min varConstMap) = Min (foldSumVarConstMap (sort varConstMap)) - --- | Simplify a 'PolyConstraint' by first 'sort'ing and then calling 'foldSumVarConstMap' on the 'VarConstMap'. -simplifyPolyConstraint :: PolyConstraint -> PolyConstraint -simplifyPolyConstraint (LEQ varConstMap rhs) = LEQ (foldSumVarConstMap (sort varConstMap)) rhs -simplifyPolyConstraint (GEQ varConstMap rhs) = GEQ (foldSumVarConstMap (sort varConstMap)) rhs -simplifyPolyConstraint (EQ varConstMap rhs) = EQ (foldSumVarConstMap (sort varConstMap)) rhs - --- | Add a sorted list of 'VarConstMap's, folding where the variables are equal -foldSumVarConstMap :: [(Integer, Rational)] -> [(Integer, Rational)] -foldSumVarConstMap [] = [] -foldSumVarConstMap [(v, c)] = [(v, c)] -foldSumVarConstMap ((v1, c1) : (v2, c2) : vcm) = - if v1 == v2 - then - let newC = c1 + c2 - in if newC == 0 - then foldSumVarConstMap vcm - else foldSumVarConstMap $ (v1, c1 + c2) : vcm - else (v1, c1) : foldSumVarConstMap ((v2, c2) : vcm) - --- | Get a map of the value of every 'Integer' variable in a 'Tableau' -displayTableauResults :: Tableau -> [(Integer, Rational)] -displayTableauResults = map (\(basicVar, (_, rhs)) -> (basicVar, rhs)) - --- | Get a map of the value of every 'Integer' variable in a 'DictionaryForm' -displayDictionaryResults :: DictionaryForm -> [(Integer, Rational)] -displayDictionaryResults dict = displayTableauResults $ dictionaryFormToTableau dict - --- | Map the given 'Integer' variable to the given 'ObjectiveFunction', for entering into 'DictionaryForm'. -createObjectiveDict :: ObjectiveFunction -> Integer -> (Integer, VarConstMap) -createObjectiveDict (Max obj) objectiveVar = (objectiveVar, obj) -createObjectiveDict (Min obj) objectiveVar = (objectiveVar, map (second negate) obj) - -{- | Converts a 'Tableau' to 'DictionaryForm'. - We do this by isolating the basic variable on the LHS, ending up with all non basic variables and a 'Rational' constant on the RHS. - (-1) is used to represent the rational constant. --} -tableauInDictionaryForm :: Tableau -> DictionaryForm -tableauInDictionaryForm [] = [] -tableauInDictionaryForm ((basicVar, (vcm, r)) : rows) = - (basicVar, (-1, r / basicCoeff) : map (\(v, c) -> (v, negate c / basicCoeff)) nonBasicVars) : tableauInDictionaryForm rows - where - basicCoeff = if null basicVars then 1 else snd $ head basicVars - (basicVars, nonBasicVars) = partition (\(v, _) -> v == basicVar) vcm - -{- | Converts a 'DictionaryForm' to a 'Tableau'. - This is done by moving all non-basic variables from the right to the left. - The rational constant (represented by the 'Integer' variable -1) stays on the right. - The basic variables will have a coefficient of 1 in the 'Tableau'. --} -dictionaryFormToTableau :: DictionaryForm -> Tableau -dictionaryFormToTableau [] = [] -dictionaryFormToTableau ((basicVar, row) : rows) = - (basicVar, ((basicVar, 1) : map (second negate) nonBasicVars, r)) : dictionaryFormToTableau rows - where - (rationalConstant, nonBasicVars) = partition (\(v, _) -> v == (-1)) row - r = if null rationalConstant then 0 else (snd . head) rationalConstant -- If there is no rational constant found in the right side, the rational constant is 0. - -{- | If this function is given 'Nothing', return 'Nothing'. - Otherwise, we 'lookup' the 'Integer' given in the first item of the pair in the map given in the second item of the pair. - This is typically used to extract the value of the 'ObjectiveFunction' after calling 'Linear.Simplex.Simplex.twoPhaseSimplex'. --} -extractObjectiveValue :: Maybe (Integer, [(Integer, Rational)]) -> Maybe Rational -extractObjectiveValue Nothing = Nothing -extractObjectiveValue (Just (objVar, results)) = - case lookup objVar results of +-- | Simplifies a system of 'StandardConstraint's by first calling 'simplifyStandardConstraint', +-- then reducing 'LEQ' and 'GEQ' with same LHS and RHS (and other similar situations) into 'EQ', +-- and finally removing duplicate elements using 'nub'. +-- simplifySystem :: [StandardConstraint] -> [StandardConstraint] +-- simplifySystem = nub . reduceSystem +-- where +-- reduceSystem :: [StandardConstraint] -> [StandardConstraint] +-- reduceSystem [] = [] +-- -- Reduce LEQ with matching GEQ and EQ into EQ +-- reduceSystem ((LEQ lhs rhs) : pcs) = +-- let matchingConstraints = +-- filter +-- ( \case +-- GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- _ -> False +-- ) +-- pcs +-- in if null matchingConstraints +-- then LEQ lhs rhs : reduceSystem pcs +-- else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) +-- -- Reduce GEQ with matching LEQ and EQ into EQ +-- reduceSystem ((GEQ lhs rhs) : pcs) = +-- let matchingConstraints = +-- filter +-- ( \case +-- LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- _ -> False +-- ) +-- pcs +-- in if null matchingConstraints +-- then GEQ lhs rhs : reduceSystem pcs +-- else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) +-- -- Reduce EQ with matching LEQ and GEQ into EQ +-- reduceSystem ((EQ lhs rhs) : pcs) = +-- let matchingConstraints = +-- filter +-- ( \case +-- LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' +-- _ -> False +-- ) +-- pcs +-- in if null matchingConstraints +-- then EQ lhs rhs : reduceSystem pcs +-- else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) + +-- | Converts a 'Dict' to a 'Tableau' using 'dictEntryToTableauEntry'. +-- FIXME: maybe remove this line. The basic variables will have a coefficient of 1 in the 'Tableau'. +dictionaryFormToTableau :: Dict -> Tableau +dictionaryFormToTableau = + Map.mapWithKey + ( \basicVar (DictValue {..}) -> + TableauRow + { lhs = Map.insert basicVar 1 $ negate <$> varMapSum + , rhs = constant + } + ) + +-- | Converts a 'Tableau' to a 'Dict'. +-- We do this by isolating the basic variable on the LHS, ending up with all non basic variables and a 'SimplexNum' constant on the RHS. +tableauInDictionaryForm :: Tableau -> Dict +tableauInDictionaryForm = + Map.mapWithKey + ( \basicVar (TableauRow {..}) -> + let basicVarCoeff = fromMaybe 1 $ Map.lookup basicVar lhs + in DictValue + { varMapSum = + Map.map + (\c -> negate c / basicVarCoeff) + $ Map.delete basicVar lhs + , constant = rhs / basicVarCoeff + } + ) + +-- | If this function is given 'Nothing', return 'Nothing'. +-- Otherwise, we 'lookup' the 'Integer' given in the first item of the pair in the map given in the second item of the pair. +-- This is typically used to extract the value of the 'ObjectiveFunction' after calling 'Linear.Simplex.Solver.TwoPhase.twoPhaseSimplex'. +extractObjectiveValue :: Maybe Result -> Maybe SimplexNum +extractObjectiveValue = fmap $ \result -> + case Map.lookup result.objectiveVar result.varValMap of Nothing -> error "Objective not found in results when extracting objective value" - r -> r + Just r -> r + +-- | Combines two 'VarLitMapSums together by summing values with matching keys +combineVarLitMapSums :: VarLitMapSum -> VarLitMapSum -> VarLitMapSum +combineVarLitMapSums = + MapMerge.merge + (MapMerge.mapMaybeMissing keepVal) + (MapMerge.mapMaybeMissing keepVal) + (MapMerge.zipWithMaybeMatched sumVals) + where + keepVal = const pure + sumVals k v1 v2 = Just $ v1 + v2 + +foldDictValue :: [DictValue] -> DictValue +foldDictValue [] = error "Empty list of DictValues given to foldDictValue" +foldDictValue [x] = x +foldDictValue + ( DictValue {varMapSum = vm1, constant = c1} + : DictValue {varMapSum = vm2, constant = c2} + : dvs + ) = + let combinedDictValue = + DictValue + { varMapSum = foldVarLitMap [vm1, vm2] + , constant = c1 + c2 + } + in foldDictValue $ combinedDictValue : dvs + +foldVarLitMap :: [VarLitMap] -> VarLitMap +foldVarLitMap [] = error "Empty list of VarLitMaps given to foldVarLitMap" +foldVarLitMap [x] = x +foldVarLitMap (vm1 : vm2 : vms) = + let combinedVars = nub $ Map.keys vm1 <> Map.keys vm2 + + combinedVarMap = + Map.fromList $ + map + ( \var -> + let mVm1VarVal = Map.lookup var vm1 + mVm2VarVal = Map.lookup var vm2 + in ( var + , case (mVm1VarVal, mVm2VarVal) of + (Just vm1VarVal, Just vm2VarVal) -> vm1VarVal + vm2VarVal + (Just vm1VarVal, Nothing) -> vm1VarVal + (Nothing, Just vm2VarVal) -> vm2VarVal + (Nothing, Nothing) -> error "Reached unreachable branch in foldDictValue" + ) + ) + combinedVars + in foldVarLitMap $ combinedVarMap : vms + +insertPivotObjectiveToDict :: PivotObjective -> Dict -> Dict +insertPivotObjectiveToDict objective = + Map.insert + objective.variable + (DictValue {varMapSum = objective.function, constant = objective.constant}) + +showT :: (Show a) => a -> T.Text +showT = T.pack . show + +logMsg :: (MonadIO m, MonadLogger m) => LogLevel -> T.Text -> m () +logMsg lvl msg = do + currTime <- T.pack . iso8601Show <$> liftIO getCurrentTime + let msgToLog = currTime <> ": " <> msg + case lvl of + LevelDebug -> $logDebug msgToLog + LevelInfo -> $logInfo msgToLog + LevelWarn -> $logWarn msgToLog + LevelError -> $logError msgToLog + LevelOther otherLvl -> error "logMsg: LevelOther is not implemented" + +extractTableauValues :: Tableau -> Map.Map Var SimplexNum +extractTableauValues = Map.map (.rhs) + +extractDictValues :: Dict -> Map.Map Var SimplexNum +extractDictValues = Map.map (.constant) + +-- scanConstraintVars :: Constraint -> [Var] +-- scanConstraintVars (Var x <= c) = x : scanConstraintVars c +-- scanConstraintVars (c <= Var x) = x : scanConstraintVars c +-- scanConstraintVars (Var x >= c) = x : scanConstraintVars c +-- scanConstraintVars (c >= Var x) = x : scanConstraintVars c +-- scanConstraintVars (Var x == c) = x : scanConstraintVars c +-- scanConstraintVars (c == Var x) = x : scanConstraintVars c +-- scanConstraintVars (c1 <= c2) = scanConstraintVars c1 <> scanConstraintVars c2 +-- scanConstraintVars (c1 >= c2) = scanConstraintVars c1 <> scanConstraintVars c2 +-- scanConstraintVars (c1 == c2) = scanConstraintVars c1 <> scanConstraintVars c2 + +-- scanConstraintVarBounds :: Constraint -> VarBounds +-- scanConstraintVarBounds (Var x <= c) = Map.singleton x (Just c, Nothing) diff --git a/src/Linear/System/Linear/Types.hs b/src/Linear/System/Linear/Types.hs new file mode 100644 index 0000000..30d06a7 --- /dev/null +++ b/src/Linear/System/Linear/Types.hs @@ -0,0 +1,31 @@ +-- | +-- Module: Linear.System.Linear.Types +-- Description: Types for linear programming problems +-- Copyright: (c) Junaid Rasheed, 2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.System.Linear.Types where + +import GHC.Generics (Generic) +import Linear.Constraint.Linear.Types (LinearEquation) +import Linear.Expr.Types (Expr) + +-- TODO: name this system of equations or something +-- TODO: OR, should I just get rid of this? +newtype LinearSystem = LinearSystem {unLinearSystem :: [LinearEquation]} + deriving (Show, Eq, Read, Generic) + +{- When would I ever want this? Do I need things to be able to be +able to turn into linear systems? Yes. But do I want other people +to be able to do that? It would be nice, but I think we do that in +the future if people actually want it. +-} +-- class CanBeLinearSystem a where +-- toLinearSystem :: a -> LinearSystem + +-- instance CanBeLinearSystem LinearSystem where +-- toLinearSystem = id + +-- instance CanBeLinearSystem LinearEquation where +-- toLinearSystem id = LinearSystem [id] diff --git a/src/Linear/System/Linear/Util.hs b/src/Linear/System/Linear/Util.hs new file mode 100644 index 0000000..4f3cc8d --- /dev/null +++ b/src/Linear/System/Linear/Util.hs @@ -0,0 +1,29 @@ +-- | +-- Module: Linear.System.Linear.Util +-- Description: Utility functions for linear programming systems +-- Copyright: (c) Junaid Rasheed, 2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.System.Linear.Util where + +import qualified Data.Set as Set +import Linear.Constraint.Linear.Types (LinearEquation (..)) +import qualified Linear.Constraint.Linear.Util as CLU +import Linear.System.Linear.Types (LinearSystem (..)) +import Linear.Var.Types (Var, VarBounds) + +-- | Prepend a linear equation to a linear system +prependLinearEquation :: LinearEquation -> LinearSystem -> LinearSystem +prependLinearEquation eq (LinearSystem eqs) = LinearSystem (eq : eqs) + +-- | Append a linear equation to a linear system +appendLinearEquation :: LinearEquation -> LinearSystem -> LinearSystem +appendLinearEquation eq (LinearSystem eqs) = LinearSystem (eqs ++ [eq]) + +findHighestVar :: LinearSystem -> Maybe Var +findHighestVar (LinearSystem []) = Nothing +findHighestVar (LinearSystem eqs) = Just $ maximum $ map CLU.findHighestVar eqs + +linearSystemVars :: LinearSystem -> Set.Set Var +linearSystemVars = Set.unions . map CLU.linearEquationVars . unLinearSystem diff --git a/src/Linear/System/Simple/Types.hs b/src/Linear/System/Simple/Types.hs new file mode 100644 index 0000000..f2e2540 --- /dev/null +++ b/src/Linear/System/Simple/Types.hs @@ -0,0 +1,57 @@ +-- | +-- Module: Linear.System.Simple.Types +-- Description: Types for the Simplex system +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.System.Simple.Types where + +import Comparison.Types (getLHS) +import qualified Data.Set as Set +import GHC.Generics (Generic) +import Linear.Constraint.Simple.Types (SimpleConstraint) +import Linear.Constraint.Simple.Util + ( simpleConstraintVars + , simplifySimpleConstraint + ) +import Linear.Expr.Util (exprVarsOnlyToList) +import Linear.System.Types (System) +import Linear.Term.Types (TermVarsOnly (..)) +import Linear.Var.Types (Var(..)) +import Test.QuickCheck (Arbitrary (..)) + +-- TODO: Use a more descriptive name +newtype SimpleSystem = SimpleSystem {unSimpleSystem :: [SimpleConstraint]} + deriving (Show, Eq, Read, Generic) + +instance Arbitrary SimpleSystem where + arbitrary = SimpleSystem <$> arbitrary + +simplifySimpleSystem :: SimpleSystem -> SimpleSystem +simplifySimpleSystem = SimpleSystem . map simplifySimpleConstraint . unSimpleSystem + +simpleSystemVars :: SimpleSystem -> Set.Set Var +simpleSystemVars = Set.unions . map simpleConstraintVars . unSimpleSystem + +findHighestVar :: SimpleSystem -> Maybe Var +findHighestVar simpleSystem = + let vars = simpleSystemVars simpleSystem + in if Set.null vars + then Nothing + else Just $ Set.findMax vars + +nextAvailableVar :: SimpleSystem -> Var +nextAvailableVar simpleSystem = + case findHighestVar simpleSystem of + Just v -> Var $ v.unVar + 1 + Nothing -> Var 0 + +class CanBeSimpleSystem a where + toSimpleSystem :: a -> SimpleSystem + +instance CanBeSimpleSystem SimpleSystem where + toSimpleSystem = id + +instance CanBeSimpleSystem System where + toSimpleSystem = undefined diff --git a/src/Linear/System/Simple/Util.hs b/src/Linear/System/Simple/Util.hs new file mode 100644 index 0000000..cfe5b2c --- /dev/null +++ b/src/Linear/System/Simple/Util.hs @@ -0,0 +1,64 @@ +-- | +-- Module: Linear.System.Simple.Util +-- Description: Utility functions for the Simplex method +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD3 +-- Maintainer: jrasheed178@gmail.com +-- Stability: experimental +module Linear.System.Simple.Util where + +import Comparison.Types + ( MixedComparison ((:<=), (:==), (:>=)) + ) +import qualified Data.Map as M +import qualified Data.Set as Set +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.System.Simple.Types + ( SimpleSystem (..) + , simpleSystemVars + ) +import Linear.Term.Types (Term (..), TermVarsOnly (..)) +import Linear.Var.Types (Bounds (..), VarBounds) + +-- | Derive bounds for all variables in a system +deriveBounds :: SimpleSystem -> VarBounds +deriveBounds simpleSystem = foldr updateBounds initialVarBounds simpleSystem.unSimpleSystem + where + systemVars = simpleSystemVars simpleSystem + initialVarBounds = M.fromList [(v, Bounds Nothing Nothing) | v <- Set.toList systemVars] + + updateBounds :: SimpleConstraint -> VarBounds -> VarBounds + updateBounds (SimpleConstraint (ExprVarsOnly [VarTermVO var] :<= num)) = M.insertWith mergeBounds var (Bounds Nothing (Just num)) + updateBounds (SimpleConstraint (ExprVarsOnly [VarTermVO var] :>= num)) = M.insertWith mergeBounds var (Bounds (Just num) Nothing) + updateBounds (SimpleConstraint (ExprVarsOnly [VarTermVO var] :== num)) = M.insertWith mergeBounds var (Bounds (Just num) (Just num)) + updateBounds _ = id + + -- \| Merge two bounds, very simple + mergeBounds :: Bounds -> Bounds -> Bounds + mergeBounds (Bounds l1 u1) (Bounds l2 u2) = Bounds (mergeLower l1 l2) (mergeUpper u1 u2) + where + mergeLower Nothing b = b + mergeLower a Nothing = a + mergeLower (Just a) (Just b) = Just (max a b) + + mergeUpper Nothing b = b + mergeUpper a Nothing = a + mergeUpper (Just a) (Just b) = Just (min a b) + +-- Eliminate inequalities which are outside the bounds +-- precondition: no zero coefficients +removeObviousInequalities :: SimpleSystem -> VarBounds -> SimpleSystem +removeObviousInequalities constraints bounds = + SimpleSystem $ + filter + ( \case + (SimpleConstraint (ExprVarsOnly [VarTermVO var] :<= num)) -> case M.lookup var bounds of + Just (Bounds _ (Just upper)) -> num <= upper + _ -> True + (SimpleConstraint (ExprVarsOnly [VarTermVO var] :>= num)) -> case M.lookup var bounds of + Just (Bounds (Just lower) _) -> num >= lower + _ -> True + _ -> True + ) + constraints.unSimpleSystem diff --git a/src/Linear/System/Types.hs b/src/Linear/System/Types.hs new file mode 100644 index 0000000..c279f49 --- /dev/null +++ b/src/Linear/System/Types.hs @@ -0,0 +1,7 @@ +module Linear.System.Types where + +import Linear.Constraint.Types (Constraint) + +-- TODO: create Sytem type, list of Constraints +newtype System = System {unSystem :: [Constraint]} + deriving (Show, Eq, Read) diff --git a/src/Linear/Tableau/Types.hs b/src/Linear/Tableau/Types.hs new file mode 100644 index 0000000..b2d5b8f --- /dev/null +++ b/src/Linear/Tableau/Types.hs @@ -0,0 +1,89 @@ +module Linear.Tableau.Types where + +import qualified Data.Set as Set +import Test.QuickCheck (Arbitrary(..)) +import qualified Test.QuickCheck as QC + +import Linear.Var.Types (SimplexNum) +import Linear.Simplex.Solver.Types (Objective) +import Linear.CanonicalForm.Types (CanonicalForm) + +-- TODO: Consider a type where the size of the list is specified +-- | A row for a @Tableau@ +data TableauRow = TableauRow { coeffs :: ![SimplexNum] + -- ^ Variable coefficients + , rhs :: !SimplexNum + -- ^ Right-hand side constants + } + deriving stock (Show, Read, Eq) + +instance Arbitrary TableauRow where + arbitrary = + TableauRow <$> QC.listOf1 arbitrary <*> arbitrary + +-- | Type representing a Simplex Tableau +data Tableau = Tableau { rows :: ![TableauRow] + , basicVars :: !(Set.Set Int) -- TODO: document the sizes are the same + -- ^ Column index of the basic vars for each row + } + deriving stock (Show, Read, Eq) + +instance Arbitrary Tableau where + arbitrary = do + n <- QC.choose (1, 5) + rows <- QC.vectorOf n arbitrary + basicVars <- Set.fromList <$> QC.vectorOf n (QC.choose (0, 10)) + pure $ Tableau rows basicVars + +rowCoeffsLength :: TableauRow -> Int +rowCoeffsLength row = length row.coeffs + +rowLength :: TableauRow -> Int +rowLength row = 1 + rowCoeffsLength row + +-- | Perform a pivot operation on a @Tableau@ +pivot :: Int -- ^ Entering column index (0-based) + -> Int -- ^ Leaving row index (0-based) + -> Tableau + -> Tableau +pivot enteringCol leavingRowIdx tableau = + let + -- Split tableau into leaving row and others + (beforeRows, pivotRow:afterRows) = splitAt leavingRowIdx tableau.rows + pivotElemCoeff = getCoeff enteringCol pivotRow + + -- Normalize pivot row + normalizedRow = TableauRow + { coeffs = map (/ pivotElemCoeff) pivotRow.coeffs + , rhs = pivotRow.rhs / pivotElemCoeff + } + + -- Update all other rows + adjustedRows = map (adjustRow normalizedRow enteringCol) (beforeRows ++ afterRows) + (adjustedBefore, adjustedAfter) = splitAt leavingRowIdx adjustedRows + newRows = adjustedBefore ++ [normalizedRow] ++ adjustedAfter + + -- Update basic variables + newBasicVars = updateBasicVar leavingRowIdx enteringCol tableau.basicVars + in + Tableau + { rows = newRows + , basicVars = newBasicVars + } + where + getCoeff col (TableauRow cs _) = cs !! col + + -- Eliminate entering col using gaussian elimination + adjustRow normalizedRow enterCol' row = + let multiplier = getCoeff enterCol' row + adjustedCoeffs = zipWith (-) (coeffs row) (map (* multiplier) (coeffs normalizedRow)) + adjustedRHS = rhs row - multiplier * rhs normalizedRow + in TableauRow adjustedCoeffs adjustedRHS + + updateBasicVar idx newVar vars = + let (before, after) = (Set.take idx vars, Set.drop (idx + 1) vars) + in before <> Set.fromList [newVar] <> after + +mkTableau :: Objective -> CanonicalForm -> Tableau +mkTableau = undefined + diff --git a/src/Linear/Term/Types.hs b/src/Linear/Term/Types.hs new file mode 100644 index 0000000..a8f2236 --- /dev/null +++ b/src/Linear/Term/Types.hs @@ -0,0 +1,38 @@ +-- | +-- Module : Linear.Term.Types +-- Description : Types for linear terms +-- Copright : (c) Junaid Rasheed, 2020-2024 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +module Linear.Term.Types where + +import GHC.Generics (Generic) +import Linear.Var.Types (SimplexNum, Var) +import Test.QuickCheck (Arbitrary (..), genericShrink, oneof) + +data Term + = ConstTerm {constant :: SimplexNum} + | CoeffTerm {coeff :: SimplexNum, var :: Var} + | VarTerm {var :: Var} + deriving (Show, Read, Eq, Ord, Generic) + +data TermVarsOnly + = VarTermVO {var :: Var} + | CoeffTermVO {coeff :: SimplexNum, var :: Var} + deriving (Show, Read, Eq, Generic) + +instance Arbitrary Term where + arbitrary = + oneof + [ ConstTerm <$> arbitrary + , CoeffTerm <$> arbitrary <*> arbitrary + , VarTerm <$> arbitrary + ] + +instance Arbitrary TermVarsOnly where + arbitrary = + oneof + [ VarTermVO <$> arbitrary + , CoeffTermVO <$> arbitrary <*> arbitrary + ] diff --git a/src/Linear/Term/Util.hs b/src/Linear/Term/Util.hs new file mode 100644 index 0000000..406293b --- /dev/null +++ b/src/Linear/Term/Util.hs @@ -0,0 +1,103 @@ +-- | +-- Module : Linear.Expr.Types +-- Description : Util functions for terms +-- Copyright : (c) Junaid Rasheed, 2020-2024 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +module Linear.Term.Util where + +import qualified Data.List as L +import Linear.Term.Types + ( Term (..) + , TermVarsOnly (..) + ) +import Linear.Var.Types (Var) + +simplifyTerm :: Term -> Term +simplifyTerm (CoeffTerm 0 _) = ConstTerm 0 +simplifyTerm (CoeffTerm 1 v) = VarTerm v +simplifyTerm t = t + +negateTerm :: Term -> Term +negateTerm (ConstTerm c) = ConstTerm (-c) +negateTerm (CoeffTerm (-1) v) = VarTerm v +negateTerm (CoeffTerm c v) = CoeffTerm (-c) v +negateTerm (VarTerm v) = CoeffTerm (-1) v + +-- Function to get all vars from a Term +termVar :: Term -> Maybe Var +termVar (VarTerm v) = Just v +termVar (CoeffTerm _ v) = Just v +termVar _ = Nothing + +zeroConstTerm :: Term -> Term +zeroConstTerm (ConstTerm _) = ConstTerm 0 +zeroConstTerm t = t + +isConstTerm :: Term -> Bool +isConstTerm (ConstTerm _) = True +isConstTerm _ = False + +-- | Normalize a list of 'Term's where each term is added together. +normalizeTerms :: [Term] -> [Term] +normalizeTerms = + L.sort + . map simplifyTerm + . combineTerms + . L.sortBy orderForCombineTerms + . map (varToCoeff . simplifyTerm) + where + orderForCombineTerms :: Term -> Term -> Ordering + orderForCombineTerms _ (VarTerm _) = error "Unexpected VarTerm in orderForCombineTerms" + orderForCombineTerms (VarTerm _) _ = error "Unexpected VarTerm in orderForCombineTerms" + orderForCombineTerms (ConstTerm c1) (ConstTerm c2) = compare c1 c2 + orderForCombineTerms (CoeffTerm c1 v1) (CoeffTerm c2 v2) = + case compare v1 v2 of + EQ -> compare c1 c2 + x -> x + orderForCombineTerms (ConstTerm _) (CoeffTerm _ _) = LT + orderForCombineTerms (CoeffTerm _ _) (ConstTerm _) = GT + + varToCoeff :: Term -> Term + varToCoeff (VarTerm v) = CoeffTerm 1 v + varToCoeff t = t + + combineTerms :: [Term] -> [Term] + combineTerms [] = [] + combineTerms [ConstTerm 0] = [] + combineTerms [CoeffTerm 0 _] = [] + combineTerms [x] = [x] + combineTerms allXs@(x1 : x2 : xs) = + case (x1, x2) of + (ConstTerm 0, _) -> combineTerms (x2 : xs) + (_, ConstTerm 0) -> combineTerms (x1 : xs) + (CoeffTerm 0 _, _) -> combineTerms (x2 : xs) + (_, CoeffTerm 0 _) -> combineTerms (x1 : xs) + (ConstTerm c1, ConstTerm c2) -> + if c1 + c2 == 0 + then combineTerms xs + else combineTerms (ConstTerm (c1 + c2) : xs) + (CoeffTerm c1 v1, CoeffTerm c2 v2) -> + if v1 == v2 + then combineTerms (CoeffTerm (c1 + c2) v1 : xs) + else x1 : combineTerms (x2 : xs) + _otherwise -> x1 : combineTerms (x2 : xs) + +normalizeTermsVarsOnly :: [TermVarsOnly] -> [TermVarsOnly] +normalizeTermsVarsOnly = map unsafeTermToTermVarsOnly . normalizeTerms . map termVarsOnlyToTerm + +termToTermVarsOnly :: Term -> Either String TermVarsOnly +termToTermVarsOnly (VarTerm v) = Right $ VarTermVO v +termToTermVarsOnly (CoeffTerm c v) = Right $ CoeffTermVO c v +termToTermVarsOnly (ConstTerm _) = Left "termToTermVarsOnly: ConstTerm not allowed" + +unsafeTermToTermVarsOnly :: Term -> TermVarsOnly +unsafeTermToTermVarsOnly t = + case termToTermVarsOnly t of + Right x -> x + Left e -> error e + +termVarsOnlyToTerm :: TermVarsOnly -> Term +termVarsOnlyToTerm (VarTermVO v) = VarTerm v +termVarsOnlyToTerm (CoeffTermVO c v) = CoeffTerm c v diff --git a/src/Linear/Var/Types.hs b/src/Linear/Var/Types.hs new file mode 100644 index 0000000..37d9df2 --- /dev/null +++ b/src/Linear/Var/Types.hs @@ -0,0 +1,19 @@ +module Linear.Var.Types where + +import qualified Data.Map as M +import GHC.Generics (Generic) +import Test.QuickCheck (Arbitrary) + +type SimplexNum = Rational + +newtype Var = Var { unVar :: Int } + deriving (Show, Read, Eq, Ord, Generic) + deriving newtype (Arbitrary) + +data Bounds = Bounds + { lowerBound :: !(Maybe SimplexNum) + , upperBound :: !(Maybe SimplexNum) + } + deriving (Show, Read, Eq, Generic) + +type VarBounds = M.Map Var Bounds diff --git a/src/Linear/Var/Util.hs b/src/Linear/Var/Util.hs new file mode 100644 index 0000000..ff95de9 --- /dev/null +++ b/src/Linear/Var/Util.hs @@ -0,0 +1,18 @@ +module Linear.Var.Util where + +import qualified Data.Map as M +import Linear.Var.Types (Var(..), Bounds (..), VarBounds) + +validateBounds :: VarBounds -> Bool +validateBounds boundsMap = all soundBounds $ M.toList boundsMap + where + soundBounds (_, Bounds lowerBound upperBound) = + case (lowerBound, upperBound) of + (Just l, Just u) -> l <= u + (_, _) -> True + +nextVar :: Var -> Var +nextVar = Var . succ . unVar + +prevVar :: Var -> Var +prevVar = Var . pred . unVar diff --git a/stack.yaml b/stack.yaml new file mode 100644 index 0000000..eab5650 --- /dev/null +++ b/stack.yaml @@ -0,0 +1,68 @@ +# This file was automatically generated by 'stack init' +# +# Some commonly used options have been documented as comments in this file. +# For advanced use and comprehensive documentation of the format, please see: +# https://docs.haskellstack.org/en/stable/yaml_configuration/ + +# Resolver to choose a 'specific' stackage snapshot or a compiler version. +# A snapshot resolver dictates the compiler version and the set of packages +# to be used for project dependencies. For example: +# +# resolver: lts-3.5 +# resolver: nightly-2015-09-21 +# resolver: ghc-7.10.2 +# +# The location of a snapshot can be provided as a file or url. Stack assumes +# a snapshot provided as a file might change, whereas a url resource does not. +# +# resolver: ./custom-snapshot.yaml +# resolver: https://example.com/snapshots/2018-01-01.yaml +resolver: lts-21.22 + +# User packages to be built. +# Various formats can be used as shown in the example below. +# +# packages: +# - some-directory +# - https://example.com/foo/bar/baz-0.0.2.tar.gz +# subdirs: +# - auto-update +# - wai +packages: +- . +# Dependency packages to be pulled from upstream that are not in the resolver. +# These entries can reference officially published versions as well as +# forks / in-progress versions pinned to a git hash. For example: +# +# extra-deps: +# - acme-missiles-0.3 +# - git: https://github.com/commercialhaskell/stack.git +# commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a +# +# extra-deps: {} + +# Override default flag values for local packages and extra-deps +# flags: {} + +# Extra package databases containing global packages +# extra-package-dbs: [] + +# Control whether we use the GHC we find on the path +# system-ghc: true +# +# Require a specific version of stack, using version ranges +# require-stack-version: -any # Default +# require-stack-version: ">=2.5" +# +# Override the architecture used by stack, especially useful on Windows +# arch: i386 +# arch: x86_64 +# +# Extra directories used by stack for building +# extra-include-dirs: [/path/to/dir] +# extra-lib-dirs: [/path/to/dir] +# +# Allow a newer minor version of GHC than the snapshot specifies +# compiler-check: newer-minor + +system-ghc: true diff --git a/stack.yaml.lock b/stack.yaml.lock new file mode 100644 index 0000000..e8d3cc7 --- /dev/null +++ b/stack.yaml.lock @@ -0,0 +1,12 @@ +# This file was autogenerated by Stack. +# You should not edit this file by hand. +# For more information, please see the documentation at: +# https://docs.haskellstack.org/en/stable/lock_files + +packages: [] +snapshots: +- completed: + sha256: afd5ba64ab602cabc2d3942d3d7e7dd6311bc626dcb415b901eaf576cb62f0ea + size: 640060 + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/21/22.yaml + original: lts-21.22 diff --git a/test/Linear/CanonicalForm/UtilSpec.hs b/test/Linear/CanonicalForm/UtilSpec.hs new file mode 100644 index 0000000..13d20fb --- /dev/null +++ b/test/Linear/CanonicalForm/UtilSpec.hs @@ -0,0 +1,415 @@ +module Linear.CanonicalForm.UtilSpec where + +import Comparison.Types + ( MixedComparison ((:<=), (:==), (:>=)) + , getLHS + ) +import Control.Monad (forM) +import Data.Functor ((<&>)) +import qualified Data.List as List +import qualified Data.Map as Map +import qualified Data.Maybe as Maybe +import qualified Data.Set as Set +import qualified Debug.Trace as T +import Linear.Constraint.Linear.Types (LinearEquation (..)) +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.CanonicalForm.Util + ( addSlackVars + , eliminateNonZeroLowerBounds + , eliminateUnrestrictedLowerBounds + ) +import Linear.System.Linear.Types (LinearSystem (..)) +import Linear.System.Simple.Types (SimpleSystem (..)) +import Linear.System.Simple.Util (deriveBounds) +import Linear.Term.Types + ( Term (..) + , TermVarsOnly (..) + ) +import Linear.Var.Types + ( Var (..) + ) +import Test.Hspec (Spec, describe, it, shouldBe) +import Test.QuickCheck (Testable (property), withMaxSuccess) + +-- data Term = ConstTerm SimplexNum | CoeffTerm SimplexNum Var | VarTerm Var -- Consider VarTerm Var - note, we must consider normalizing this: Considered. It makes going to standard form easier due to type safety +-- deriving (Show, Read, Eq, Ord, Generic) + +-- TODO: consider type NumberConstraint = MixedComparison SimplexNum SimplexNum +spec :: Spec +spec = describe "Slack Form Transformations" $ do + it + "eliminateNonZeroLowerBounds does not do anything when all lower bounds are zero" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + expectedEliminatedVarExprMap = Map.empty + updatedSystem `shouldBe` expectedSimpleSystem + updatedBounds `shouldBe` expectedEliminatedVarExprMap + it "eliminateNonZeroLowerBounds correctly eliminates positive lower bounds" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 2)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + expectedEliminatedVarExprMap = Map.fromList [(Var 0, Expr (VarTerm (Var 2) : [ConstTerm 1]))] + updatedSystem `shouldBe` expectedSimpleSystem + updatedBounds `shouldBe` expectedEliminatedVarExprMap + it "eliminateNonZeroLowerBounds correctly eliminates negative lower bounds" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= (-1) + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 2)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + ] + expectedEliminatedVarExprMap = Map.fromList [(Var 0, Expr (VarTerm (Var 2) : [ConstTerm (-1)]))] + updatedSystem `shouldBe` expectedSimpleSystem + updatedBounds `shouldBe` expectedEliminatedVarExprMap + it + "eliminateNonZeroLowerBounds correctly eliminates positive and negative lower bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= (-1) + ] + (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 2)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 3)] :>= 0 + ] + expectedEliminatedVarExprMap = + Map.fromList + [ (Var 0, Expr (VarTerm (Var 2) : [ConstTerm 1])) + , (Var 1, Expr (VarTerm (Var 3) : [ConstTerm (-1)])) + ] + updatedSystem `shouldBe` expectedSimpleSystem + updatedBounds `shouldBe` expectedEliminatedVarExprMap + it + "eliminateNonZeroLowerBounds correctly substitutes vars with non-zero lower bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 0) : [VarTermVO (Var 1)]) :>= 1 + ] + (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 2)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [VarTermVO (Var 2)]) :>= 0 + ] + expectedEliminatedVarExprMap = Map.fromList [(Var 0, Expr (VarTerm (Var 2) : [ConstTerm 1]))] + updatedSystem `shouldBe` expectedSimpleSystem + updatedBounds `shouldBe` expectedEliminatedVarExprMap + it "eliminateNonZeroLowerBounds property based test lower bounds" $ withMaxSuccess 5 $ property $ \simpleSystem -> do + let (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + all + ( \case + SimpleConstraint (ExprVarsOnly [VarTermVO (Var _)] :>= num) -> num == 0 + _ -> True + ) + updatedSystem.unSimpleSystem + it "eliminateNonZeroLowerBounds property based test map" $ withMaxSuccess 5 $ property $ \simpleSystem -> do + let (updatedBounds, updatedSystem) = eliminateNonZeroLowerBounds simpleSystem Map.empty + all + ( \(var, _) -> + any + ( \(SimpleConstraint constraint) -> + let getVars _a = [] + lhs = getLHS constraint + allVars = getVars lhs + in var `notElem` allVars + ) + updatedSystem.unSimpleSystem + ) + (Map.toList updatedBounds) + it + "addSlackVars correctly transforms inequalities to equalities (wikipedia case)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 2) : [CoeffTermVO 2 (Var 3)]) :<= 3 -- x_2 + 2x_3 <= 3 + , SimpleConstraint $ ExprVarsOnly (CoeffTermVO (-1) (Var 4) : [CoeffTermVO 3 (Var 5)]) :>= 2 -- -x_4 + 3x_5 >= 2 + ] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 2) : [CoeffTermVO 2 (Var 3), VarTermVO (Var 6)])) 3 -- x_2 + 2x_3 + x_6 = 3 + , LinearEquation + (ExprVarsOnly (CoeffTermVO (-1) (Var 4) : [CoeffTermVO 3 (Var 5), CoeffTermVO (-1) (Var 7)])) + 2 -- -x_4 + 3x_5 + x_7 = 2 + ] + expectedSlackVars = [Var 6, Var 7] + (slackVars, updatedSystem) = addSlackVars simpleSystem + updatedSystem `shouldBe` expectedSystem + slackVars `shouldBe` expectedSlackVars + it + "addSlackVars correctly transforms inequalities to equalities (test case 1)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2)]) :<= 4 -- x_1 + 2x_2 <= 4 + , SimpleConstraint $ ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)]) :>= 3 -- -x_3 + 2x_4 >= 3 + ] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2), VarTermVO (Var 5)])) 4 -- x_1 + 2x_2 + x_5 = 4 + , LinearEquation + (ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4), CoeffTermVO (-1) (Var 6)])) + 3 -- -x_3 + 2x_4 - x_6 = 3 + ] + expectedSlackVars = [Var 5, Var 6] + (slackVars, updatedSystem) = addSlackVars simpleSystem + updatedSystem `shouldBe` expectedSystem + slackVars `shouldBe` expectedSlackVars + it + "addSlackVars correctly transforms inequalities to equalities (test case 2)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2)]) :<= 5 -- x_1 + 2x_2 <= 5 + , SimpleConstraint $ ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)]) :>= 4 -- -x_3 + 2x_4 >= 4 + ] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2), VarTermVO (Var 5)])) 5 -- x_1 + 2x_2 + x_5 = 5 + , LinearEquation + (ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4), CoeffTermVO (-1) (Var 6)])) + 4 -- -x_3 + 2x_4 - x_6 = 4 + ] + expectedSlackVars = [Var 5, Var 6] + (slackVars, updatedSystem) = addSlackVars simpleSystem + updatedSystem `shouldBe` expectedSystem + slackVars `shouldBe` expectedSlackVars + it + "addSlackVars correctly transforms inequalities to equalities (test case 3)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2)]) :<= 5 -- x_1 + 2x_2 <= 5 + , SimpleConstraint $ ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)]) :== 4 -- -x_3 + 2x_4 = 4 + ] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2), VarTermVO (Var 5)])) 5 -- x_1 + 2x_2 + x_5 = 5 + , LinearEquation (ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)])) 4 -- -x_3 + 2x_4 = 4 + ] + expectedSlackVars = [Var 5] + (slackVars, updatedSystem) = addSlackVars simpleSystem + updatedSystem `shouldBe` expectedSystem + slackVars `shouldBe` expectedSlackVars + it + "addSlackVars correctly transforms inequalities to equalities (test case 4)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2)]) :== 5 -- x_1 + 2x_2 = 5 + , SimpleConstraint $ ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)]) :== 4 -- -x_3 + 2x_4 = 4 + ] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 2)])) 5 -- x_1 + 2x_2 = 5 + , LinearEquation (ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [CoeffTermVO 2 (Var 4)])) 4 -- -x_3 + 2x_4 = 4 + ] + expectedSlackVars = [] + (slackVars, updatedSystem) = addSlackVars simpleSystem + updatedSystem `shouldBe` expectedSystem + slackVars `shouldBe` expectedSlackVars + it + "eliminateUnrestrictedLowerBounds correctly eliminates unrestricted lower bounds (wikipedia case)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 -- x_1 >= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [VarTermVO (Var 2)]) :>= 0 -- x_1 + x_2 >= 0 + ] + systemBounds = deriveBounds simpleSystem + (eliminatedNonZeroLowerBounds, systemWithoutNonZeroLowerBounds) = eliminateNonZeroLowerBounds simpleSystem Map.empty + (slackVars, systemWithSlackVars) = addSlackVars systemWithoutNonZeroLowerBounds + (updatedEliminatedVarsMap, updatedSystem) = + eliminateUnrestrictedLowerBounds + systemWithSlackVars + systemBounds + eliminatedNonZeroLowerBounds + expectedSlackVars = [Var 3, Var 4] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (CoeffTermVO (-1) (Var 3) : [VarTermVO (Var 1)])) 0 -- -x_3 + x_1 = 0 + , LinearEquation + ( ExprVarsOnly + (CoeffTermVO (-1) (Var 4) : [CoeffTermVO (-1) (Var 6), VarTermVO (Var 1), VarTermVO (Var 5)]) + ) + 0 -- -x_4 - x_6 + x_1 + x_5 = 0 + ] + expectedEliminatedVarExprMap = Map.fromList [(Var 2, Expr (VarTerm (Var 5) : [CoeffTerm (-1) (Var 6)]))] + + slackVars `shouldBe` expectedSlackVars + updatedSystem `shouldBe` expectedSystem + updatedEliminatedVarsMap `shouldBe` expectedEliminatedVarExprMap + it + "eliminateUnrestrictedLowerBounds correctly eliminates unrestricted lower bounds (test case 2)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 -- x_1 >= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [VarTermVO (Var 2), VarTermVO (Var 3)]) :>= 0 -- x_1 + x_2 + x_3 >= 0 + ] + systemBounds = deriveBounds simpleSystem + (eliminatedNonZeroLowerBounds, systemWithoutNonZeroLowerBounds) = eliminateNonZeroLowerBounds simpleSystem Map.empty + (slackVars, systemWithSlackVars) = addSlackVars systemWithoutNonZeroLowerBounds + (updatedEliminatedVarsMap, updatedSystem) = + eliminateUnrestrictedLowerBounds + systemWithSlackVars + systemBounds + eliminatedNonZeroLowerBounds + expectedSlackVars = [Var 4, Var 5] + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (CoeffTermVO (-1) (Var 4) : [VarTermVO (Var 1)])) 0 -- -x_4 + x_1 = 0 + , LinearEquation + ( ExprVarsOnly + ( CoeffTermVO (-1) (Var 5) + : [CoeffTermVO (-1) (Var 7), CoeffTermVO (-1) (Var 9), VarTermVO (Var 1), VarTermVO (Var 6), VarTermVO (Var 8)] + ) + ) + 0 -- -x_5 - x_7 - x_9 + x_1 + x_6 + x_8 = 0 + ] + expectedEliminatedVarExprMap = + Map.fromList + [ (Var 2, Expr (VarTerm (Var 6) : [CoeffTerm (-1) (Var 7)])) + , (Var 3, Expr (VarTerm (Var 8) : [CoeffTerm (-1) (Var 9)])) + ] + slackVars `shouldBe` expectedSlackVars + updatedSystem `shouldBe` expectedSystem + updatedEliminatedVarsMap `shouldBe` expectedEliminatedVarExprMap + + it + "eliminateUnrestrictedLowerBounds correctly eliminates unrestricted lower bounds (test case 3)" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 -- x_1 >= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [VarTermVO (Var 2)]) :>= 0 -- x_1 + x_2 >= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 2) : [VarTermVO (Var 3)]) :>= 0 -- x_2 + x_3 >= 0 + ] + systemBounds = deriveBounds simpleSystem + (eliminatedNonZeroLowerBounds, systemWithoutNonZeroLowerBounds) = eliminateNonZeroLowerBounds simpleSystem Map.empty + (slackVars, systemWithSlackVars) = addSlackVars systemWithoutNonZeroLowerBounds + expectedSlackVars = [Var 4, Var 5, Var 6] + (updatedEliminatedVarsMap, updatedSystem) = + eliminateUnrestrictedLowerBounds + systemWithSlackVars + systemBounds + eliminatedNonZeroLowerBounds + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (CoeffTermVO (-1) (Var 4) : [VarTermVO (Var 1)])) 0 -- -x_4 + x_1 = 0 + , LinearEquation + ( ExprVarsOnly + (CoeffTermVO (-1) (Var 5) : [CoeffTermVO (-1) (Var 8), VarTermVO (Var 1), VarTermVO (Var 7)]) + ) + 0 -- -x_5 - x_8 + x_1 + x_7 = 0 + , LinearEquation + ( ExprVarsOnly + ( CoeffTermVO (-1) (Var 6) + : [CoeffTermVO (-1) (Var 8), CoeffTermVO (-1) (Var 10), VarTermVO (Var 7), VarTermVO (Var 9)] + ) + ) + 0 -- -x_6 - x_8 - x_10 + x_7 + x_9 = 0 + ] + expectedEliminatedVarExprMap = + Map.fromList + [ (Var 2, Expr (VarTerm (Var 7) : [CoeffTerm (-1) (Var 8)])) + , (Var 3, Expr (VarTerm (Var 9) : [CoeffTerm (-1) (Var 10)])) + ] + slackVars `shouldBe` expectedSlackVars + updatedSystem `shouldBe` expectedSystem + updatedEliminatedVarsMap `shouldBe` expectedEliminatedVarExprMap + it + "eliminateUnrestrictedLowerBounds correctly eliminates non-zero lower bounds for all variables" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO 2 (Var 1)]) :>= 0 -- x_1 + 2x_1 >= 0 + , SimpleConstraint $ ExprVarsOnly (VarTermVO (Var 2) : [CoeffTermVO 3 (Var 1)]) :>= 0 -- x_2 + 3x_1 >= 0 + ] + systemBounds = deriveBounds simpleSystem + (eliminatedNonZeroLowerBounds, systemWithoutNonZeroLowerBounds) = eliminateNonZeroLowerBounds simpleSystem Map.empty + (slackVars, systemWithSlackVars) = addSlackVars systemWithoutNonZeroLowerBounds + expectedSlackVars = [Var 3, Var 4] + (updatedEliminatedVarsMap, updatedSystem) = + eliminateUnrestrictedLowerBounds + systemWithSlackVars + systemBounds + eliminatedNonZeroLowerBounds + expectedSystem = + LinearSystem + [ LinearEquation + (ExprVarsOnly (CoeffTermVO (-3) (Var 6) : [CoeffTermVO (-1) (Var 3), CoeffTermVO 3 (Var 5)])) + 0 -- -3x_6 - x_3 + 3x_5 = 0 + , LinearEquation + ( ExprVarsOnly + ( CoeffTermVO (-3) (Var 6) + : [CoeffTermVO (-1) (Var 4), CoeffTermVO (-1) (Var 8), CoeffTermVO 3 (Var 5), VarTermVO (Var 7)] + ) + ) + 0 -- -3x_6 - x_4 - x_8 + 3x_5 + x_7 = 0 + ] + expectedEliminatedVarExprMap = + Map.fromList + [ (Var 1, Expr (VarTerm (Var 5) : [CoeffTerm (-1) (Var 6)])) + , (Var 2, Expr (VarTerm (Var 7) : [CoeffTerm (-1) (Var 8)])) + ] + slackVars `shouldBe` expectedSlackVars + updatedSystem `shouldBe` expectedSystem + updatedEliminatedVarsMap `shouldBe` expectedEliminatedVarExprMap + + it + "eliminateUnrestrictedLowerBounds correctly handles all variables with zero lower bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 -- x_1 >= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 2)] :>= 0 -- x_2 >= 0 + ] + systemBounds = deriveBounds simpleSystem + (eliminatedNonZeroLowerBounds, systemWithoutNonZeroLowerBounds) = eliminateNonZeroLowerBounds simpleSystem Map.empty + (slackVars, systemWithSlackVars) = addSlackVars systemWithoutNonZeroLowerBounds + expectedSlackVars = [Var 3, Var 4] + (updatedEliminatedVarsMap, updatedSystem) = + eliminateUnrestrictedLowerBounds + systemWithSlackVars + systemBounds + eliminatedNonZeroLowerBounds + expectedSystem = + LinearSystem + [ LinearEquation (ExprVarsOnly (VarTermVO (Var 1) : [CoeffTermVO (-1) (Var 3)])) 0 -- x_1 - x_3 = 0 + , LinearEquation (ExprVarsOnly (VarTermVO (Var 2) : [CoeffTermVO (-1) (Var 4)])) 0 -- x_2 - x_4 = 0 + ] + expectedEliminatedVarExprMap = Map.empty + slackVars `shouldBe` expectedSlackVars + updatedSystem `shouldBe` expectedSystem + updatedEliminatedVarsMap `shouldBe` expectedEliminatedVarExprMap diff --git a/test/Linear/Constraint/Simple/UtilSpec.hs b/test/Linear/Constraint/Simple/UtilSpec.hs new file mode 100644 index 0000000..bb1d78c --- /dev/null +++ b/test/Linear/Constraint/Simple/UtilSpec.hs @@ -0,0 +1,163 @@ +-- | +-- Module: Linear.Constraint.Simple.TypesSpec +-- Description: Tests for Linear.Constraint.Simple.Types +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.Constraint.Simple.UtilSpec where + +import qualified Data.Map as Map +import qualified Data.Set as Set +import Linear.Constraint.Simple.Util + ( constraintToSimpleConstraint + , simpleConstraintVars + , simplifyCoeff + , simplifySimpleConstraint + , substVarSimpleConstraint + , substVarSimpleConstraintExpr + ) +import Linear.Constraint.Util (constraintVars) +import Linear.Expr.Types (Expr (..), ExprVarsOnly (..)) +import Linear.Expr.Util (exprVars, exprVarsOnlyVars) +import Linear.Term.Types (Term (..), TermVarsOnly (..)) +import Test.Hspec (Spec, describe) +import Test.Hspec.QuickCheck (prop) +import Test.QuickCheck (counterexample, elements) +import TestUtil + ( evalConstraint + , evalExpr + , evalExprVarsOnly + , evalSimpleConstraint + , genVarMap + ) +import Prelude + +spec :: Spec +spec = do + describe "SimpleConstraint" $ do + prop + "substVarSimpleConstraintExpr with a constant is the same as evaluating with the variable mapped to the constant" + $ \simpleConstraint c -> do + let vars = Set.toList $ simpleConstraintVars simpleConstraint + var <- elements vars + let varReplacement = Expr (ConstTerm c : []) + initialVarMap <- genVarMap vars + let varMap = Map.insert var c initialVarMap + substitutedSimpleConstraint = substVarSimpleConstraintExpr var varReplacement simpleConstraint + substitutedSimpleConstraintEval = evalSimpleConstraint varMap substitutedSimpleConstraint + simpleConstraintEval = evalSimpleConstraint varMap simpleConstraint + pure + $ counterexample + ( "simpleConstraint: " + <> show simpleConstraint + <> "\nvar: " + <> show var + <> "\nconst: " + <> show c + <> "\nvarMap: " + <> show varMap + <> "\nvarReplacement: " + <> show varReplacement + <> "\nsubstitutedSimpleConstraint: " + <> show substitutedSimpleConstraint + <> "\nsubstitutedSimpleConstraintEval: " + <> show substitutedSimpleConstraintEval + <> "\nsimpleConstraintEval: " + <> show simpleConstraintEval + ) + $ substitutedSimpleConstraintEval == simpleConstraintEval + prop + "substVarSimpleConstraintExpr with an expr is the same as evaluating with the variable mapped to the expr" + $ \simpleConstraint exprReplacement -> do + let vars = + Set.toList $ simpleConstraintVars simpleConstraint <> exprVars exprReplacement + var <- elements vars + initialVarMap <- genVarMap vars + let exprReplacementEval = evalExpr initialVarMap exprReplacement + varMap = Map.insert var exprReplacementEval initialVarMap + substitutedSimpleConstraint = substVarSimpleConstraintExpr var exprReplacement simpleConstraint + simpleConstraintEval = evalSimpleConstraint varMap simpleConstraint + substitutedSimpleConstraintEval = evalSimpleConstraint initialVarMap substitutedSimpleConstraint + pure + $ counterexample + ( "simpleConstraint: " + <> show simpleConstraint + <> "\nvar: " + <> show var + <> "\nexprReplacement: " + <> show exprReplacement + <> "\ninitialVarMap: " + <> show initialVarMap + <> "\nexprReplacementEval: " + <> show exprReplacementEval + <> "\nvarMap: " + <> show varMap + <> "\nsubstitutedSimpleConstraint: " + <> show substitutedSimpleConstraint + <> "\nsimpleConstraintEval: " + <> show simpleConstraintEval + <> "\nsubstitutedSimpleConstraintEval: " + <> show substitutedSimpleConstraintEval + ) + $ substitutedSimpleConstraintEval == simpleConstraintEval + prop "constraintToSimpleConstraint leads to the same evaluation" $ \constraint -> do + let vars = Set.toList $ constraintVars constraint + varMap <- genVarMap vars + let simpleConstraint = constraintToSimpleConstraint constraint + constraintEval = evalConstraint varMap constraint + simpleConstraintEval = evalSimpleConstraint varMap simpleConstraint + pure + $ counterexample + ( "constraint: " + <> show constraint + <> "\nsimpleConstraint: " + <> show simpleConstraint + <> "\ninitialVarMap: " + <> show varMap + <> "\nconstraintEval: " + <> show constraintEval + <> "\nsimpleConstraintEval" + <> show simpleConstraintEval + ) + $ constraintEval == simpleConstraintEval + prop "simplifyCoeff leads to the same evaluation" $ \simpleConstraint -> do + let vars = Set.toList $ simpleConstraintVars simpleConstraint + varMap <- genVarMap vars + let simplifiedSimpleConstraint = simplifyCoeff simpleConstraint + simpleConstraintEval = evalSimpleConstraint varMap simpleConstraint + simplifiedSimpleConstraintEval = evalSimpleConstraint varMap simplifiedSimpleConstraint + pure + $ counterexample + ( "simpleConstraint: " + <> show simpleConstraint + <> "\nsimplifiedSimpleConstraint: " + <> show simplifiedSimpleConstraint + <> "\ninitialVarMap: " + <> show varMap + <> "\nsimpleConstraintEval: " + <> show simpleConstraintEval + <> "\nsimplifiedSimpleConstraintEval: " + <> show simplifiedSimpleConstraintEval + ) + $ simpleConstraintEval == simplifiedSimpleConstraintEval + prop "simplifySimpleConstraint leads to the same evaluation" $ \simpleConstraint -> do + let vars = Set.toList $ simpleConstraintVars simpleConstraint + varMap <- genVarMap vars + let simplifiedSimpleConstraint = simplifySimpleConstraint simpleConstraint + simpleConstraintEval = evalSimpleConstraint varMap simpleConstraint + simplifiedSimpleConstraintEval = evalSimpleConstraint varMap simplifiedSimpleConstraint + pure + $ counterexample + ( "simpleConstraint: " + <> show simpleConstraint + <> "\nsimplifiedSimpleConstraint: " + <> show simplifiedSimpleConstraint + <> "\ninitialVarMap: " + <> show varMap + <> "\nsimpleConstraintEval: " + <> show simpleConstraintEval + <> "\nsimplifiedSimpleConstraintEval: " + <> show simplifiedSimpleConstraintEval + ) + $ simpleConstraintEval == simplifiedSimpleConstraintEval diff --git a/test/Linear/Expr/UtilSpec.hs b/test/Linear/Expr/UtilSpec.hs new file mode 100644 index 0000000..5cc925a --- /dev/null +++ b/test/Linear/Expr/UtilSpec.hs @@ -0,0 +1,268 @@ +-- | +-- Module: Linear.Expr.TypesSpec +-- Description: Tests for Linear.Expr.Types +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.Expr.UtilSpec where + +import qualified Data.List.NonEmpty as NE +import qualified Data.Map as Map +import qualified Data.Set as Set +import Linear.Expr.Types (Expr (..)) +import Linear.Expr.Util + ( addExpr + , exprToList + , exprVars + , listToExpr + , negateExpr + , simplifyExpr + , substVarExpr + , subtractExpr + , sumExprConstTerms + , zeroConstExpr + ) +import Linear.Term.Types (Term (..)) +import Test.Hspec (Spec, describe) +import Test.Hspec.QuickCheck (prop) +import Test.QuickCheck (counterexample, elements) +import TestUtil (evalExpr, genVarMap) + +spec :: Spec +spec = do + describe "Expr" $ do + prop "simplifying leads to same evaluation" $ \expr -> do + let vars = Set.toList $ exprVars expr + varMap <- genVarMap vars + pure $ evalExpr varMap (simplifyExpr expr) == evalExpr varMap expr + prop "simplifying twice is the same as simplifying once" $ \expr -> do + -- expr: (Expr (CoeffTerm (1 % 1) 0) :+ CoeffTerm (1 % 2) (-1)) :+ CoeffTerm (2 % 1) (-1) + -- 1 = x + -- -1 = y + -- 1x + 0.5y + 2y + -- simplifiedExpr: (Expr (CoeffTerm (1 % 2) (-1)) :+ CoeffTerm (2 % 1) (-1)) :+ VarTerm 0 + -- 0.5y + 2y + x + -- simplifiedTwiceExpr: Expr (CoeffTerm (5 % 2) (-1)) :+ VarTerm 0 + -- 2.5y + x + let simplifiedExpr = simplifyExpr expr + let simplifiedTwiceExpr = simplifyExpr simplifiedExpr + counterexample + ( "expr: " + <> show expr + <> "\nsimplifiedExpr: " + <> show simplifiedExpr + <> "\nsimplifiedTwiceExpr: " + <> show simplifiedTwiceExpr + ) + $ simplifiedTwiceExpr == simplifiedExpr + prop "composing listToExpr and exprToList is the same as id" $ \expr -> do + let exprConvertedTwice = listToExpr (exprToList expr) + counterexample + ( "expr: " + <> show expr + <> "\nexprConvertedTwice: " + <> show exprConvertedTwice + ) + $ exprConvertedTwice == expr + prop "summing c terms is the same as evaluating the expr with zero coefficients" $ \expr -> do + let vars = Set.toList $ exprVars expr + varMap = Map.fromList $ map (,0) vars + constSum = sumExprConstTerms expr + exprEval = evalExpr varMap expr + counterexample + ( "expr: " + <> show expr + <> "\nvarMap: " + <> show varMap + <> "\nconstSum: " + <> show constSum + <> "\nexprEval: " + <> show exprEval + ) + $ constSum == exprEval + prop "negating and evaluating is the same as negating the evaluation" $ \expr -> do + let vars = Set.toList $ exprVars expr + varMap <- genVarMap vars + let negatedExpr = negateExpr expr + exprEval = evalExpr varMap expr + exprEvalNegated = negate exprEval + negatedExprEval = evalExpr varMap negatedExpr + pure + $ counterexample + ( "expr: " + <> show expr + <> "\nnegatedExpr: " + <> show negatedExpr + <> "\nvarMap: " + <> show varMap + <> "\nexprEval: " + <> show exprEval + <> "\nexprEvalNegated: " + <> show exprEvalNegated + <> "\nnegatedExprEval: " + <> show negatedExprEval + ) + $ exprEvalNegated == negatedExprEval + prop "negating a simplified expression twice is the same as not negating" $ \expr -> do + let simplifiedExpr = simplifyExpr expr + negatedTwiceSimpleExpr = negateExpr (negateExpr simplifiedExpr) + counterexample + ( "expr: " + <> show expr + <> "\nsimplifiedExpr: " + <> show simplifiedExpr + <> "\nnegatedTwiceSimpleExpr: " + <> show negatedTwiceSimpleExpr + ) + $ negatedTwiceSimpleExpr == simplifiedExpr + prop "addExpr is the same as evaluating the sum of the exprs" $ \expr1 expr2 -> do + let vars = Set.toList $ exprVars expr1 <> exprVars expr2 + varMap <- genVarMap vars + let addExpr1Expr2 = addExpr expr1 expr2 + addExpr1Expr2Eval = evalExpr varMap addExpr1Expr2 + expr1Eval = evalExpr varMap expr1 + expr2Eval = evalExpr varMap expr2 + sumExpr1EvalExpr2Eval = expr1Eval + expr2Eval + pure + $ counterexample + ( "expr1: " + <> show expr1 + <> "\nexpr2: " + <> show expr2 + <> "\nvarMap:" + <> show varMap + <> "\naddExpr1Expr2: " + <> show addExpr1Expr2 + <> "\naddExpr1Expr2Eval: " + <> show addExpr1Expr2Eval + <> "\nexpr1Eval: " + <> show expr1Eval + <> "\nexpr2Eval: " + <> show expr2Eval + <> "\nexpr1EvalPlusExpr2Eval: " + <> show sumExpr1EvalExpr2Eval + ) + $ addExpr1Expr2Eval == sumExpr1EvalExpr2Eval + prop "subtractExpr is the same as evaluating the difference of the exprs" $ \expr1 expr2 -> do + let vars = Set.toList $ exprVars expr1 <> exprVars expr2 + varMap <- genVarMap vars + let subtractExpr1Expr2 = subtractExpr expr1 expr2 + subtractExpr1Expr2Eval = evalExpr varMap subtractExpr1Expr2 + expr1Eval = evalExpr varMap expr1 + expr2Eval = evalExpr varMap expr2 + diffExpr1EvalExpr2Eval = expr1Eval - expr2Eval + pure + $ counterexample + ( "expr1: " + <> show expr1 + <> "\nexpr2: " + <> show expr2 + <> "\nvarMap:" + <> show varMap + <> "\nsubtractExpr1Expr2: " + <> show subtractExpr1Expr2 + <> "\nsubtractExpr1Expr2Eval: " + <> show subtractExpr1Expr2Eval + <> "\nexpr1Eval: " + <> show expr1Eval + <> "\nexpr2Eval: " + <> show expr2Eval + <> "\nexpr1EvalMinusExpr2Eval: " + <> show diffExpr1EvalExpr2Eval + ) + $ subtractExpr1Expr2Eval == diffExpr1EvalExpr2Eval + prop "substVarExpr with the same variable is the same as simplfying" $ \expr -> do + let vars = Set.toList $ exprVars expr + var <- elements vars + varMap <- genVarMap vars + let varReplacement = Expr (VarTerm var : []) + exprSubst = substVarExpr var varReplacement expr + exprSimplified = simplifyExpr expr + exprSubstEval = evalExpr varMap exprSubst + exprSimplifiedEval = evalExpr varMap exprSimplified + pure + $ counterexample + ( "expr: " + <> show expr + <> "\nvar: " + <> show var + <> "\nvarMap: " + <> show varMap + <> "\nvarReplacement: " + <> show varReplacement + <> "\nsubstVarExpr: " + <> show exprSubst + <> "\nsimplifyExpr: " + <> show exprSimplified + <> "\nexprSubstEval: " + <> show exprSubstEval + <> "\nexprSimplifiedEval: " + <> show exprSimplifiedEval + ) + $ exprSubstEval == exprSimplifiedEval + prop + "substVarExpr with a constant is the same as evaluating with the variable mapped to the constant" + $ \expr c -> do + let varReplacement = Expr (ConstTerm c : []) + let vars = Set.toList $ exprVars expr + var <- elements vars + initialVarMap <- genVarMap vars + let varMap = Map.insert var c initialVarMap + substitutedExpr = substVarExpr var varReplacement expr + substitutedExprEval = evalExpr varMap substitutedExpr + exprEval = evalExpr varMap expr + pure + $ counterexample + ( "expr: " + <> show expr + <> "\nvar: " + <> show var + <> "\nconst: " + <> show c + <> "\nvarMap: " + <> show varMap + <> "\nvarReplacement: " + <> show varReplacement + <> "\nsubstitutedExpr: " + <> show substitutedExpr + <> "\nsubstitutedExprEval: " + <> show substitutedExprEval + <> "\nexprEval: " + <> show exprEval + ) + $ evalExpr varMap (substVarExpr var varReplacement expr) == evalExpr varMap expr + prop + "substVarExpr with an expr is the same as evaluating with the variable mapped to the expr" + $ \expr exprReplacement -> do + let vars = Set.toList $ exprVars expr <> exprVars exprReplacement + var <- elements vars + initialVarMap <- genVarMap vars + let exprReplacementEval = evalExpr initialVarMap exprReplacement + varMap = Map.insert var exprReplacementEval initialVarMap + substitutedExpr = substVarExpr var exprReplacement expr + exprEval = evalExpr varMap expr + substitutedExprEval = evalExpr initialVarMap substitutedExpr + pure + $ counterexample + ( "expr: " + <> show expr + <> "\nvar: " + <> show var + <> "\nexprReplacement: " + <> show exprReplacement + <> "\ninitialVarMap: " + <> show initialVarMap + <> "\nexprReplacementEval: " + <> show exprReplacementEval + <> "\nvarMap: " + <> show varMap + <> "\nsubstExpr: " + <> show substitutedExpr + <> "\nexprEval: " + <> show exprEval + <> "\nsubstExprEval: " + <> show substitutedExprEval + ) + $ substitutedExprEval == exprEval + prop "zeroConstExpr correctly zeroes constant terms in expressions" $ \expr -> sumExprConstTerms (zeroConstExpr expr) == 0 diff --git a/test/Linear/System/Simple/UtilSpec.hs b/test/Linear/System/Simple/UtilSpec.hs new file mode 100644 index 0000000..bfde917 --- /dev/null +++ b/test/Linear/System/Simple/UtilSpec.hs @@ -0,0 +1,287 @@ +-- | +-- Module: Linear.System.Simple.TypesSpec +-- Description: Tests for Linear.System.Simple.Types +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed =)) + ) +import Data.List.NonEmpty (NonEmpty (..)) +import qualified Data.Map as Map +import qualified Data.Set as Set +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Expr.Types (ExprVarsOnly (..)) +import Linear.System.Simple.Types + ( SimpleSystem (SimpleSystem) + , findHighestVar + , simpleSystemVars + , simplifySimpleSystem + ) +import Linear.System.Simple.Util + ( deriveBounds + , removeObviousInequalities + ) +import Linear.Term.Types (TermVarsOnly (..)) +import Linear.Var.Types (Bounds (..), Var(..)) +import Linear.Var.Util (validateBounds) +import Test.Hspec (Spec, describe, it, shouldBe) +import Test.Hspec.QuickCheck (prop) +import Test.QuickCheck (counterexample) +import TestUtil (evalSimpleSystem, genVarMap) + +spec :: Spec +spec = do + describe "SimpleSystem" $ do + prop "simplifySimpleSystem leads to the same evaluation" $ \simpleSystem -> do + let vars = Set.toList $ simpleSystemVars simpleSystem + varMap <- genVarMap vars + let simplifiedSimpleSystem = simplifySimpleSystem simpleSystem + simpleSystemEval = evalSimpleSystem varMap simpleSystem + simplifiedSimpleSystemEval = evalSimpleSystem varMap simplifiedSimpleSystem + pure + $ counterexample + ( "simpleSystem: " + <> show simpleSystem + <> "\nsimplifiedSimpleSystem: " + <> show simplifiedSimpleSystem + <> "\ninitialVarMap: " + <> show varMap + <> "\nsimpleSystemEval: " + <> show simpleSystemEval + <> "\nsimplifiedSimpleSystemEval: " + <> show simplifiedSimpleSystemEval + ) + $ simpleSystemEval == simplifiedSimpleSystemEval + it "findHighestVar finds the highest variable in a simple system" $ do + let simpleSystem1 = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :<= 1 + ] + simpleSystem100 = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 50)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 100)] :<= 1 + ] + simpleSystem10 = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var (-10))] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 10)] :<= 1 + ] + simpleSystemMinus10 = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var (-10))] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var (-20))] :<= 1 + ] + + findHighestVar (SimpleSystem []) `shouldBe` Nothing + findHighestVar simpleSystem1 `shouldBe` Just (Var 1) + findHighestVar simpleSystem100 `shouldBe` Just (Var 100) + findHighestVar simpleSystem10 `shouldBe` Just (Var 10) + findHighestVar simpleSystemMinus10 `shouldBe` Just (Var (-10)) + describe "Bounds" + $ it + "validateBounds finds that deriving bounds for a system where -1 <= x <= 1 has valid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= (-1) + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just (-1)) (Just 1))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` True + it + "validateBounds finds that deriving bounds for a system where 0 <= x <= 1 has valid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 1))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` True + it + "validateBounds finds that deriving bounds for a system where 1 <= x <= 1 has valid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 1) (Just 1))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` True + it + "validateBounds finds that deriving bounds for a system where 1 <= x <= 0 has invalid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 0 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 1) (Just 0))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` False + it + "validateBounds finds that deriving bounds for a system where 0 <= x <= 1 and 1 <= y <= 3 has valid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :<= 3 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 1)), (Var 1, Bounds (Just 1) (Just 3))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` True + it + "validateBounds finds that deriving bounds for a system where 1 <= x <= 0 and 3 <= y <= 1 has invalid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 3 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :<= 1 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 1) (Just 0)), (Var 1, Bounds (Just 3) (Just 1))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` False + it + "validateBounds finds that deriving bounds for a system where 1 <= x <= 0 and 1 <= y <= 3 has invalid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :<= 3 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 1) (Just 0)), (Var 1, Bounds (Just 1) (Just 3))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` False + it + "validateBounds finds that deriving bounds for a system where 0 <= x <= 1 and 3 <= y <= 1 has invalid bounds" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 1 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :>= 3 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 1)] :<= 1 + ] + derivedBounds = deriveBounds simpleSystem + expectedBounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 1)), (Var 1, Bounds (Just 3) (Just 1))] + derivedBounds `shouldBe` expectedBounds + validateBounds derivedBounds `shouldBe` False + it "removeObviousInequalities removes x <= 3 when bounds say x <= 2" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 3 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 2))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it "removeObviousInequalities does not remove x <= 2 when bounds say x <= 2" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 2))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it "removeObviousInequalities removes x >= 3 when bounds say x >= 4" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 4 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 3 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 4) (Just 5))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 4] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it "removeObviousInequalities does not remove x >= 4 when bounds say x >= 4" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 4 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 4) (Just 5))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 4] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it + "removeObviousInequalities does not remove 0 <= x <= 2 when bounds say 0 <= x <= 2" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 2))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + ] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it + "removeObviousInequalities removes upper bound of 0 <= x <= 2 when bounds say 0 <= x <= 1" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 1))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it + "removeObviousInequalities removes lower bound of 0 <= x <= 2 when bounds say 1 <= x <= 2" + $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :>= 0 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 1) (Just 2))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = SimpleSystem [SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem + it "removeUselssSystemBounds only removes constraints of the form x <= c" $ do + let simpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + , SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 3 + , SimpleConstraint $ ExprVarsOnly [CoeffTermVO 2 (Var 0)] :<= 6 + ] + bounds = Map.fromList [(Var 0, Bounds (Just 0) (Just 2))] + simplifiedSimpleSystem = removeObviousInequalities simpleSystem bounds + expectedSimpleSystem = + SimpleSystem + [ SimpleConstraint $ ExprVarsOnly [VarTermVO (Var 0)] :<= 2 + , SimpleConstraint $ ExprVarsOnly [CoeffTermVO 2 (Var 0)] :<= 6 + ] + simplifiedSimpleSystem `shouldBe` expectedSimpleSystem diff --git a/test/Linear/Tableau/TypesSpec.hs b/test/Linear/Tableau/TypesSpec.hs new file mode 100644 index 0000000..2600bbc --- /dev/null +++ b/test/Linear/Tableau/TypesSpec.hs @@ -0,0 +1,78 @@ +module Linear.Tableau.TypesSpec where + +import Data.List (nub) +import qualified Data.Set as Set +import Linear.Tableau.Types + +import Test.Hspec +import Test.Hspec.QuickCheck (prop) +import Test.QuickCheck + +spec :: Spec +spec = + describe "Tableau Pivoting" $ do + let -- Sample initial tableau + initialTableau = Tableau + { rows = + [ TableauRow { coeffs = [1, 2], rhs = 3 } -- x₁ + 2x₂ = 3 + , TableauRow { coeffs = [4, 5], rhs = 6 } -- 4x₁ +5x₂ = 6 + ] + , basicVars = Set.fromList [2, 3] -- x₃ basic in row 0, x₄ in row 1 + } + + it "handles basic pivot operation" $ do + let result = pivot 0 0 initialTableau + result.rows `shouldBe` + [ TableauRow { coeffs = [1, 2], rhs = 3 } -- x₁ remains same + , TableauRow { coeffs = [0, -3], rhs = -6 } -- 0x₁ -3x₂ = -6 + ] + result.basicVars `shouldBe` Set.fromList [0, 3] -- x₁ now basic in row 0 + + it "handles pivot with negative coefficients" $ do + let tab = Tableau + { rows = [ TableauRow { coeffs = [-2, 4], rhs = 6 } + , TableauRow { coeffs = [3, -1], rhs = 2 } + ] + , basicVars = Set.fromList [1, 2] + } + result = pivot 0 0 tab + + result.rows `shouldBe` + [ TableauRow { coeffs = [1, -2], rhs = -3 } -- Normalized row + , TableauRow { coeffs = [0, 5], rhs = 11 } -- Eliminated row: 3 - (3 * -3) = 0, -1 - (3 * -2) = 5 + ] + result.basicVars `shouldBe` Set.fromList [0, 2] + + it "maintains row count" $ -- TODO: Why am I testing this? + property $ \(Positive n) -> + let tab = Tableau { rows = (replicate n (TableauRow { coeffs = [], rhs = 0 })), basicVars = Set.fromList (replicate n 0) } + in length (rows (pivot 0 0 tab)) `shouldBe` n + + it "handles fractional results" $ do + let tab = Tableau + { rows = [ TableauRow { coeffs = [2, 4], rhs = 6 } ] + , basicVars = Set.fromList [0] + } + result = pivot 0 0 tab + + coeffs (head (rows result)) `shouldBe` [1, 2] + rhs (head (rows result)) `shouldBe` 3 + + it "pivots the wikipeda worked example correctly" $ do + let tab = Tableau + { rows = [ TableauRow { coeffs = [1, 2, 3, 4, 0, 0], rhs = 0 } + , TableauRow { coeffs = [0, 3, 2, 1, 1, 0], rhs = 10 } + , TableauRow { coeffs = [0, 2, 5, 3, 0, 1], rhs = 15 } + ] + , basicVars = Set.fromList [0, 4, 5] + } + let step1 = pivot 3 2 tab + let step1Expected = Tableau + { rows = [ TableauRow { coeffs = [1, -(2 / 3), -(11 / 3), 0, 0, -(4 / 3)], rhs = -20 } + , TableauRow { coeffs = [0, 7/3, 1/3, 0, 1,- (1 / 3)], rhs = 5 } + , TableauRow { coeffs = [0, 2/3, 5/3, 1, 0, 1/3], rhs = 5 } + ] + , basicVars = Set.fromList [0, 4, 3] + } + step1 `shouldBe` step1Expected + diff --git a/test/Linear/Term/UtilSpec.hs b/test/Linear/Term/UtilSpec.hs new file mode 100644 index 0000000..b276c29 --- /dev/null +++ b/test/Linear/Term/UtilSpec.hs @@ -0,0 +1,133 @@ +-- | +-- Module: Linear.Term.TypesSpec +-- Description: Tests for Linear.Term.Types +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module Linear.Term.UtilSpec where + +import qualified Data.Either as Either +import qualified Data.List as List +import qualified Data.Map as Map +import Linear.Term.Types + ( Term (..) + , TermVarsOnly (..) + ) +import Linear.Term.Util + ( isConstTerm + , negateTerm + , normalizeTerms + , simplifyTerm + , termToTermVarsOnly + , termVar + , unsafeTermToTermVarsOnly + , zeroConstTerm + ) +import Test.Hspec (Spec, describe, it, shouldBe, shouldSatisfy) +import Test.Hspec.QuickCheck (prop) +import Test.QuickCheck (counterexample) +import TestUtil (evalTerm, genVarMap) +import Prelude +import Linear.Var.Types (Var(..)) + +spec :: Spec +spec = describe "Term" $ do + prop "simplifying leads to same evaluation" $ \term -> do + varMap <- maybe (pure Map.empty) (genVarMap . List.singleton) $ termVar term + let simplifiedTerm = simplifyTerm term + termEval = evalTerm varMap term + simplifiedTermEval = evalTerm varMap simplifiedTerm + pure + $ counterexample + ( "term: " + <> show term + <> "simplifiedTerm: " + <> show simplifiedTerm + <> "\nvarMap: " + <> show varMap + <> "\ntermEval: " + <> show termEval + <> "\nsimplifiedTermEval: " + <> show simplifiedTermEval + ) + $ evalTerm varMap (simplifyTerm term) == evalTerm varMap term + prop "simplifyTerm is idempotent" $ \term -> do + let simplifiedTerm = simplifyTerm term + simplifiedTwiceTerm = simplifyTerm simplifiedTerm + counterexample + ( "term: " + <> show term + <> "\nsimplifiedTerm: " + <> show simplifiedTerm + <> "\nsimplifiedTwiceTerm: " + <> show simplifiedTwiceTerm + ) + $ simplifiedTwiceTerm == simplifiedTerm + prop "negating and evaluating is the same as negating the evaluation" $ \term -> do + varMap <- maybe (pure $ Map.empty) (genVarMap . List.singleton) $ termVar term + let negatedTerm = negateTerm term + termEval = evalTerm varMap term + negatedTermEval = evalTerm varMap negatedTerm + pure + $ counterexample + ( "term: " + <> show term + <> "\nnegatedTerm: " + <> show negatedTerm + <> "\nvarMap: " + <> show varMap + <> "\ntermEval: " + <> show termEval + <> "\nnegatedTermEval: " + <> show negatedTermEval + ) + $ negate termEval == negatedTermEval + prop "negating twice is the same as not negating" $ \term -> do + let simplifiedTerm = simplifyTerm term + negatedTwiceSimpleTerm = negateTerm (negateTerm simplifiedTerm) + counterexample + ( "term: " + <> show term + <> "\nsimplifiedTerm: " + <> show simplifiedTerm + <> "\nnegatedTwiceSimpleTerm: " + <> show negatedTwiceSimpleTerm + ) + $ negatedTwiceSimpleTerm == simplifiedTerm + prop "zeroConstTerm correctly zeroes constant terms" $ \term -> do + let termZeroedConsts = zeroConstTerm term + counterexample + ( "term: " + <> show term + <> "\ntermZeroedConsts: " + <> show termZeroedConsts + ) + $ case term of + ConstTerm _ -> termZeroedConsts == ConstTerm 0 + _ -> termZeroedConsts == term + it "isConstTerm correctly identifies constant terms" $ do + isConstTerm (ConstTerm 0) `shouldBe` True + isConstTerm (ConstTerm 1) `shouldBe` True + isConstTerm (CoeffTerm 1 (Var 1)) `shouldBe` False + isConstTerm (VarTerm (Var 1)) `shouldBe` False + it "termToTermVarsOnly correctly converts terms to vars only" $ do + termToTermVarsOnly (ConstTerm 0) `shouldSatisfy` Either.isLeft + termToTermVarsOnly (ConstTerm 1) `shouldSatisfy` Either.isLeft + termToTermVarsOnly (CoeffTerm 1 (Var 1)) `shouldBe` Right (CoeffTermVO 1 (Var 1)) + termToTermVarsOnly (VarTerm (Var 1)) `shouldBe` Right (VarTermVO (Var 1)) + it "unsafeTermToTermVarsOnly correctly converts terms without vars" $ do + unsafeTermToTermVarsOnly (CoeffTerm 1 (Var 1)) `shouldBe` (CoeffTermVO 1 (Var 1)) + unsafeTermToTermVarsOnly (VarTerm (Var 1)) `shouldBe` (VarTermVO (Var 1)) + prop "normalizeTerms is idempotent" $ \terms -> do + let normalizedTerms = normalizeTerms terms + normalizedTwiceTerms = normalizeTerms normalizedTerms + counterexample + ( "terms: " + <> show terms + <> "\nnormalizedTerms: " + <> show normalizedTerms + <> "\nnormalizedTwiceTerms: " + <> show normalizedTwiceTerms + ) + $ normalizedTwiceTerms == normalizedTerms diff --git a/test/Linear/Var/UtilSpec.hs b/test/Linear/Var/UtilSpec.hs new file mode 100644 index 0000000..58c76cf --- /dev/null +++ b/test/Linear/Var/UtilSpec.hs @@ -0,0 +1,20 @@ +module Linear.Var.UtilSpec where + +import qualified Data.Map as Map +import Linear.Var.Types ( Bounds(..), Var(..) ) +import Linear.Var.Util (validateBounds) +import Test.Hspec (Spec, describe, it, shouldBe) + +spec :: Spec +spec = describe "Bounds" $ do + it "validateBounds returns true for valid bounds" $ do + validateBounds (Map.fromList [(Var 1, Bounds (Just 1) (Just 2))]) `shouldBe` True + validateBounds (Map.fromList [(Var 1, Bounds (Just 1.1) (Just 1.2))]) + `shouldBe` True + validateBounds (Map.fromList [(Var 1, Bounds (Just 1) Nothing)]) `shouldBe` True + validateBounds (Map.fromList [(Var 1, Bounds Nothing (Just 2))]) `shouldBe` True + validateBounds (Map.fromList [(Var 1, Bounds Nothing Nothing)]) `shouldBe` True + it "validateBounds returns false for invalid bounds" $ do + validateBounds (Map.fromList [(Var 1, Bounds (Just 2) (Just 1))]) `shouldBe` False + validateBounds (Map.fromList [(Var 1, Bounds (Just 1.2) (Just 1.1))]) + `shouldBe` False diff --git a/test/Spec.hs b/test/Spec.hs index f1c9919..a824f8c 100644 --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,28 +1 @@ -module Main where - -import Linear.Simplex.Prettify -import Linear.Simplex.Simplex -import Linear.Simplex.Util -import TestFunctions - -main :: IO () -main = runTests testsList - -runTests [] = putStrLn "All tests passed" -runTests (((testObjective, testConstraints), expectedResult) : tests) = - let testResult = twoPhaseSimplex testObjective testConstraints - in if testResult == expectedResult - then runTests tests - else do - putStrLn "The following test failed: \n" - putStrLn ("Objective Function (Non-prettified): " ++ show testObjective) - putStrLn ("Constraints (Non-prettified): " ++ show testConstraints) - putStrLn "====================================\n" - putStrLn ("Objective Function (Prettified): " ++ prettyShowObjectiveFunction testObjective) - putStrLn "Constraints (Prettified): " - putStrLn (concatMap ((\c -> "\t" ++ prettyShowPolyConstraint c ++ "\n")) testConstraints) - putStrLn "====================================\n" - putStrLn ("Expected Solution (Full): " ++ show expectedResult) - putStrLn ("Actual Solution (Full): " ++ show testResult) - putStrLn ("Expected Solution (Objective): " ++ show (extractObjectiveValue expectedResult)) - putStrLn ("Actual Solution (Objective): " ++ show (extractObjectiveValue testResult)) +{-# OPTIONS_GHC -F -pgmF hspec-discover #-} diff --git a/test/TestFunctions.hs b/test/TestFunctions.hs deleted file mode 100644 index cc0852b..0000000 --- a/test/TestFunctions.hs +++ /dev/null @@ -1,1024 +0,0 @@ -module TestFunctions where - -import Data.Ratio -import Linear.Simplex.Types -import Prelude hiding (EQ) - -testsList :: [((ObjectiveFunction, [PolyConstraint]), Maybe (Integer, [(Integer, Rational)]))] -testsList = - [ (test1, Just (7, [(7, 29 % 1), (1, 3 % 1), (2, 4 % 1)])) - , (test2, Just (7, [(7, 0 % 1)])) - , (test3, Nothing) - , (test4, Just (11, [(11, 237 % 7), (1, 24 % 7), (2, 33 % 7)])) - , (test5, Just (9, [(9, 3 % 5), (2, 14 % 5), (3, 17 % 5)])) - , (test6, Nothing) - , (test7, Just (8, [(8, 1 % 1), (2, 2 % 1), (1, 3 % 1)])) - , (test8, Just (8, [(8, (-1) % 4), (2, 9 % 2), (1, 17 % 4)])) - , (test9, Just (7, [(7, 5 % 1), (3, 2 % 1), (4, 1 % 1)])) - , (test10, Just (7, [(7, 8 % 1), (1, 2 % 1), (2, 6 % 1)])) - , (test11, Just (8, [(8, 20 % 1), (4, 16 % 1), (3, 6 % 1)])) - , (test12, Just (8, [(8, 6 % 1), (4, 2 % 1), (5, 2 % 1)])) - , (test13, Just (6, [(6, 150 % 1), (2, 150 % 1)])) - , (test14, Just (6, [(6, 40 % 3), (2, 40 % 3)])) - , (test15, Nothing) - , (test16, Just (6, [(6, 75 % 1), (1, 75 % 2)])) - , (test17, Just (7, [(7, (-120) % 1), (1, 20 % 1)])) - , (test18, Just (7, [(7, 10 % 1), (3, 5 % 1)])) - , (test19, Nothing) - , (test20, Nothing) - , (test21, Just (7, [(7, 250 % 1), (2, 50 % 1)])) - , (test22, Just (7, [(7, 0 % 1)])) - , (test23, Nothing) - , (test24, Just (10, [(10, 300 % 1), (3, 150 % 1)])) - , (test25, Just (3, [(3, 15 % 1), (1, 15 % 1)])) - , (test26, Just (6, [(6, 20 % 1), (1, 10 % 1), (2, 10 % 1)])) - , (test27, Just (3, [(3, 0 % 1)])) - , (test28, Just (6, [(6, 0 % 1), (2, 10 % 1)])) - , (test29, Nothing) - , (test30, Nothing) - , (testPolyPaver1, Just (12, [(12, 7 % 4), (2, 5 % 2), (1, 7 % 4), (3, 0 % 1)])) - , (testPolyPaver2, Just (12, [(12, 5 % 2), (2, 5 % 3), (1, 5 % 2), (3, 0 % 1)])) - , (testPolyPaver3, Just (12, [(12, 5 % 3), (2, 5 % 3), (1, 5 % 2), (3, 0 % 1)])) - , (testPolyPaver4, Just (12, [(12, 5 % 2), (2, 5 % 2), (1, 5 % 2), (3, 0 % 1)])) - , (testPolyPaver5, Nothing) - , (testPolyPaver6, Nothing) - , (testPolyPaver7, Nothing) - , (testPolyPaver8, Nothing) - , (testPolyPaver9, Just (12, [(12, 7 % 2), (2, 5 % 9), (1, 7 % 2), (3, 0 % 1)])) - , (testPolyPaver10, Just (12, [(12, 17 % 20), (2, 7 % 2), (1, 17 % 20), (3, 0 % 1)])) - , (testPolyPaver11, Just (12, [(12, 7 % 2), (2, 7 % 2), (1, 22 % 9)])) - , (testPolyPaver12, Just (12, [(12, 5 % 9), (2, 5 % 9), (1, 7 % 2), (3, 0 % 1)])) - , (testPolyPaverTwoFs1, Nothing) - , (testPolyPaverTwoFs2, Nothing) - , (testPolyPaverTwoFs3, Nothing) - , (testPolyPaverTwoFs4, Nothing) - , (testPolyPaverTwoFs5, Just (17, [(17, 5 % 2), (2, 45 % 22), (1, 5 % 2), (4, 0 % 1)])) - , (testPolyPaverTwoFs6, Just (17, [(17, 45 % 22), (2, 5 % 2), (1, 45 % 22), (4, 0 % 1)])) - , (testPolyPaverTwoFs7, Just (17, [(17, 5 % 2), (2, 5 % 2), (1, 5 % 2), (4, 0 % 1)])) - , (testPolyPaverTwoFs8, Just (17, [(17, 45 % 22), (2, 45 % 22), (1, 5 % 2), (4, 0 % 1)])) - , (testLeqGeqBugMin1, Just (5, [(5, 3 % 1), (1, 3 % 1), (2, 3 % 1)])) - , (testLeqGeqBugMax1, Just (5, [(5, 3 % 1), (1, 3 % 1), (2, 3 % 1)])) - , (testLeqGeqBugMin2, Just (5, [(5, 3 % 1), (1, 3 % 1), (2, 3 % 1)])) - , (testLeqGeqBugMax2, Just (5, [(5, 3 % 1), (1, 3 % 1), (2, 3 % 1)])) - , (testQuickCheck1, Just (10, [(10, (-370) % 1), (2, 26 % 1), (1, 5 % 3)])) - , (testQuickCheck2, Just (8, [(8, (-2) % 9), (1, 14 % 9), (2, 8 % 9)])) - , (testQuickCheck3, Just (7, [(7, (-8) % 1), (2, 2 % 1)])) - ] - -testLeqGeqBugMin1 = - ( Min [(1, 1)] - , - [ GEQ [(1, 1 % 1)] (3 % 1) - , LEQ [(1, 1 % 1)] (3 % 1) - , GEQ [(2, 1 % 1)] (3 % 1) - , LEQ [(2, 1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMax1 = - ( Min [(1, 1)] - , - [ GEQ [(1, 1 % 1)] (3 % 1) - , LEQ [(1, 1 % 1)] (3 % 1) - , GEQ [(2, 1 % 1)] (3 % 1) - , LEQ [(2, 1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMin2 = - ( Min [(1, 1)] - , - [ GEQ [(1, 1 % 1)] (3 % 1) - , LEQ [(1, 1 % 1)] (3 % 1) - , GEQ [(2, 1 % 1)] (3 % 1) - , LEQ [(2, 1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMax2 = - ( Min [(1, 1)] - , - [ GEQ [(1, 1 % 1)] (3 % 1) - , LEQ [(1, 1 % 1)] (3 % 1) - , GEQ [(2, 1 % 1)] (3 % 1) - , LEQ [(2, 1 % 1)] (3 % 1) - ] - ) - --- From page 50 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 29, 1 = 3, 2 = 4, -test1 :: (ObjectiveFunction, [PolyConstraint]) -test1 = - ( Max [(1, 3), (2, 5)] - , - [ LEQ [(1, 3), (2, 1)] 15 - , LEQ [(1, 1), (2, 1)] 7 - , LEQ [(2, 1)] 4 - , LEQ [(1, -1), (2, 2)] 6 - ] - ) - -test2 :: (ObjectiveFunction, [PolyConstraint]) -test2 = - ( Min [(1, 3), (2, 5)] - , - [ LEQ [(1, 3), (2, 1)] 15 - , LEQ [(1, 1), (2, 1)] 7 - , LEQ [(2, 1)] 4 - , LEQ [(1, -1), (2, 2)] 6 - ] - ) - -test3 :: (ObjectiveFunction, [PolyConstraint]) -test3 = - ( Max [(1, 3), (2, 5)] - , - [ GEQ [(1, 3), (2, 1)] 15 - , GEQ [(1, 1), (2, 1)] 7 - , GEQ [(2, 1)] 4 - , GEQ [(1, -1), (2, 2)] 6 - ] - ) - -test4 :: (ObjectiveFunction, [PolyConstraint]) -test4 = - ( Min [(1, 3), (2, 5)] - , - [ GEQ [(1, 3), (2, 1)] 15 - , GEQ [(1, 1), (2, 1)] 7 - , GEQ [(2, 1)] 4 - , GEQ [(1, -1), (2, 2)] 6 - ] - ) - --- From https://www.eng.uwaterloo.ca/~syde05/phase1.pdf --- Solution: obj = 3/5, 2 = 14/5, 3 = 17/5 --- requires two phases -test5 :: (ObjectiveFunction, [PolyConstraint]) -test5 = - ( Max [(1, 1), (2, -1), (3, 1)] - , - [ LEQ [(1, 2), (2, -1), (3, 2)] 4 - , LEQ [(1, 2), (2, -3), (3, 1)] (-5) - , LEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - -test6 :: (ObjectiveFunction, [PolyConstraint]) -test6 = - ( Min [(1, 1), (2, -1), (3, 1)] - , - [ LEQ [(1, 2), (2, -1), (3, 2)] 4 - , LEQ [(1, 2), (2, -3), (3, 1)] (-5) - , LEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - -test7 :: (ObjectiveFunction, [PolyConstraint]) -test7 = - ( Max [(1, 1), (2, -1), (3, 1)] - , - [ GEQ [(1, 2), (2, -1), (3, 2)] 4 - , GEQ [(1, 2), (2, -3), (3, 1)] (-5) - , GEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - -test8 :: (ObjectiveFunction, [PolyConstraint]) -test8 = - ( Min [(1, 1), (2, -1), (3, 1)] - , - [ GEQ [(1, 2), (2, -1), (3, 2)] 4 - , GEQ [(1, 2), (2, -3), (3, 1)] (-5) - , GEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - --- From page 49 of 'Linear and Integer Programming Made Easy' --- Solution: obj = -5, 3 = 2, 4 = 1, objVar was negated so actual val is 5 wa --- requires two phases -test9 :: (ObjectiveFunction, [PolyConstraint]) -test9 = - ( Min [(1, 1), (2, 1), (3, 2), (4, 1)] - , - [ EQ [(1, 1), (3, 2), (4, -2)] 2 - , EQ [(2, 1), (3, 1), (4, 4)] 6 - ] - ) - -test10 :: (ObjectiveFunction, [PolyConstraint]) -test10 = - ( Max [(1, 1), (2, 1), (3, 2), (4, 1)] - , - [ EQ [(1, 1), (3, 2), (4, -2)] 2 - , EQ [(2, 1), (3, 1), (4, 4)] 6 - ] - ) - --- Adapted from page 52 of 'Linear and Integer Programming Made Easy' --- Removed variables which do not appear in the system (these should be artificial variables) --- Solution: obj = 20, 3 = 6, 4 = 16 wq -test11 :: (ObjectiveFunction, [PolyConstraint]) -test11 = - ( Max [(3, -2), (4, 2), (5, 1)] - , - [ EQ [(3, -2), (4, 1), (5, 1)] 4 - , EQ [(3, 3), (4, -1), (5, 2)] 2 - ] - ) - -test12 :: (ObjectiveFunction, [PolyConstraint]) -test12 = - ( Min [(3, -2), (4, 2), (5, 1)] - , - [ EQ [(3, -2), (4, 1), (5, 1)] 4 - , EQ [(3, 3), (4, -1), (5, 2)] 2 - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 150, 1 = 0, 2 = 150 --- requires two phases -test13 :: (ObjectiveFunction, [PolyConstraint]) -test13 = - ( Max [(1, 2), (2, 1)] - , - [ LEQ [(1, 4), (2, 1)] 150 - , LEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test14 :: (ObjectiveFunction, [PolyConstraint]) -test14 = - ( Min [(1, 2), (2, 1)] - , - [ LEQ [(1, 4), (2, 1)] 150 - , LEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test15 :: (ObjectiveFunction, [PolyConstraint]) -test15 = - ( Max [(1, 2), (2, 1)] - , - [ GEQ [(1, 4), (2, 1)] 150 - , GEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test16 :: (ObjectiveFunction, [PolyConstraint]) -test16 = - ( Min [(1, 2), (2, 1)] - , - [ GEQ [(1, 4), (2, 1)] 150 - , GEQ [(1, 2), (2, -3)] (-40) - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 120, 1 = 20, 2 = 0, 3 = 0, objVar was negated so actual val is -120 -test17 :: (ObjectiveFunction, [PolyConstraint]) -test17 = - ( Min [(1, -6), (2, -4), (3, 2)] - , - [ LEQ [(1, 1), (2, 1), (3, 4)] 20 - , LEQ [(2, -5), (3, 5)] 100 - , LEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test18 :: (ObjectiveFunction, [PolyConstraint]) -test18 = - ( Max [(1, -6), (2, -4), (3, 2)] - , - [ LEQ [(1, 1), (2, 1), (3, 4)] 20 - , LEQ [(2, -5), (3, 5)] 100 - , LEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test19 :: (ObjectiveFunction, [PolyConstraint]) -test19 = - ( Min [(1, -6), (2, -4), (3, 2)] - , - [ GEQ [(1, 1), (2, 1), (3, 4)] 20 - , GEQ [(2, -5), (3, 5)] 100 - , GEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test20 :: (ObjectiveFunction, [PolyConstraint]) -test20 = - ( Max [(1, -6), (2, -4), (3, 2)] - , - [ GEQ [(1, 1), (2, 1), (3, 4)] 20 - , GEQ [(2, -5), (3, 5)] 100 - , GEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 250, 1 = 0, 2 = 50, 3 = 0 -test21 :: (ObjectiveFunction, [PolyConstraint]) -test21 = - ( Max [(1, 3), (2, 5), (3, 2)] - , - [ LEQ [(1, 5), (2, 1), (3, 4)] 50 - , LEQ [(1, 1), (2, -1), (3, 1)] 150 - , LEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test22 :: (ObjectiveFunction, [PolyConstraint]) -test22 = - ( Min [(1, 3), (2, 5), (3, 2)] - , - [ LEQ [(1, 5), (2, 1), (3, 4)] 50 - , LEQ [(1, 1), (2, -1), (3, 1)] 150 - , LEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test23 :: (ObjectiveFunction, [PolyConstraint]) -test23 = - ( Max [(1, 3), (2, 5), (3, 2)] - , - [ GEQ [(1, 5), (2, 1), (3, 4)] 50 - , GEQ [(1, 1), (2, -1), (3, 1)] 150 - , GEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test24 :: (ObjectiveFunction, [PolyConstraint]) -test24 = - ( Min [(1, 3), (2, 5), (3, 2)] - , - [ GEQ [(1, 5), (2, 1), (3, 4)] 50 - , GEQ [(1, 1), (2, -1), (3, 1)] 150 - , GEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test25 :: (ObjectiveFunction, [PolyConstraint]) -test25 = - ( Max [(1, 1)] - , - [ LEQ [(1, 1)] 15 - ] - ) - -test26 :: (ObjectiveFunction, [PolyConstraint]) -test26 = - ( Max [(1, 2)] - , - [ LEQ [(1, 2)] 20 - , GEQ [(2, 1)] 10 - ] - ) - -test27 :: (ObjectiveFunction, [PolyConstraint]) -test27 = - ( Min [(1, 1)] - , - [ LEQ [(1, 1)] 15 - ] - ) - -test28 :: (ObjectiveFunction, [PolyConstraint]) -test28 = - ( Min [(1, 2)] - , - [ LEQ [(1, 2)] 20 - , GEQ [(2, 1)] 10 - ] - ) - -test29 :: (ObjectiveFunction, [PolyConstraint]) -test29 = - ( Max [(1, 1)] - , - [ LEQ [(1, 1)] 15 - , GEQ [(1, 1)] 15.01 - ] - ) - -test30 :: (ObjectiveFunction, [PolyConstraint]) -test30 = - ( Max [(1, 1)] - , - [ LEQ [(1, 1)] 15 - , GEQ [(1, 1)] 15.01 - , GEQ [(2, 1)] 10 - ] - ) - --- Tests for systems similar to those from PolyPaver2 -testPolyPaver1 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver1 = - ( Min [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver2 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver2 = - ( Max [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver3 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver3 = - ( Min [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver4 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver4 = - ( Max [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver5 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver5 = - ( Max [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver6 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver6 = - ( Min [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver7 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver7 = - ( Max [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver8 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver8 = - ( Min [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver9 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver9 = - ( Max [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver10 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver10 = - ( Min [(1, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver11 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver11 = - ( Max [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver12 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver12 = - ( Min [(2, 1)] - , - [ LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)) -- -5 - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaverTwoFs1 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs1 = - ( Max [(1, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs2 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs2 = - ( Min [(1, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs3 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs3 = - ( Max [(2, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs4 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs4 = - ( Min [(2, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs5 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs5 = - ( Max [(1, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l -- don't need variable >= 0, already assumed - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs6 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs6 = - ( Min [(1, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l -- don't need variable >= 0, already assumed - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs7 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs7 = - ( Max [(2, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l -- don't need variable >= 0, already assumed - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs8 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs8 = - ( Min [(2, 1)] - , - [ LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)) -- -4, This will need an artificial variable - , GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)) - , LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)) - , GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)) - , GEQ [(1, 1)] x1l -- don't need variable >= 0, already assumed - , LEQ [(1, 1)] x1r - , GEQ [(2, 1)] x2l - , LEQ [(2, 1)] x2r - , LEQ [(3, 1)] 0 - , LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - --- Test cases produced by old simplex-haskell/SoPlex QuickCheck prop - --- SoPlex gives -400 for the following system but -370 is the optimized solution --- simplex-haskell gives -370 --- SoPlex gives -370 if we simplify the system before sending it to SoPlex -testQuickCheck1 = - ( Max [(1, -6), (1, -8), (1, 9), (1, 10), (1, 8), (2, -15), (1, 13), (1, -14), (2, 0)] - , - [ EQ [(1, 5), (1, 6), (2, -2), (1, 7), (1, 6), (2, 0)] (-12) - , GEQ [(1, 11), (1, 0), (1, -5), (1, -12), (1, -14), (2, 11)] (-7) - , GEQ [(1, -12), (1, -7), (1, -2), (2, -9), (1, 3), (1, 5), (1, -15), (2, 14)] (-8) - , GEQ [(1, 13), (1, 1), (1, -11), (2, 0)] 5 - , LEQ [(1, -10), (1, -14), (1, 4), (1, -2), (1, -10), (1, -5), (1, -11)] (-1) - ] - ) - --- If we do not call simplifyPolyConstraints before we start the simplex algorithm, the following return a wrong solution --- Correct solution is -2/9 -testQuickCheck2 = - ( Max [(1, -3), (2, 5)] - , - [ LEQ [(2, -1), (1, -6), (2, 7)] 4 - , LEQ [(1, 1), (2, -4), (3, 3)] (-2) - , LEQ [(2, 6), (1, -4), (2, 1)] 0 - ] - ) - --- This test will fail if the objective function is not simplified -testQuickCheck3 = - ( Min [(2, 0), (2, -4)] - , - [ GEQ [(1, 5), (2, 4)] (-4) - , LEQ [(1, -1), (2, -1)] 2 - , LEQ [(2, 1)] 2 - , GEQ [(1, -5), (2, -1), (2, 1)] (-5) - ] - ) diff --git a/test/TestUtil.hs b/test/TestUtil.hs new file mode 100644 index 0000000..bf753d0 --- /dev/null +++ b/test/TestUtil.hs @@ -0,0 +1,80 @@ +-- | +-- Module: TestUtil +-- Description: Utility functions for testing +-- Copyright: (c) Junaid Rasheed, 2020-2024 +-- License: BSD-3 +-- Maintainer: Junaid Rasheed +-- Stability: experimental +module TestUtil where + +import Comparison.Types + ( MixedComparison ((:<=), (:==), (:>=)) + ) +import Control.Monad (forM) +import qualified Data.Map as Map +import Linear.Constraint.Simple.Types (SimpleConstraint (..)) +import Linear.Constraint.Types + ( Constraint (..) + ) +import Linear.Expr.Types (Expr, ExprVarsOnly) +import Linear.Expr.Util (exprToList, exprVarsOnlyToExpr) +import Linear.System.Simple.Types (SimpleSystem (..)) +import Linear.Term.Types + ( Term (..) + , TermVarsOnly + ) +import Linear.Term.Util (termVarsOnlyToTerm) +import Linear.Var.Types (SimplexNum, Var) +import Test.QuickCheck (Arbitrary (..), Gen) +import Prelude + +evalTerm :: Map.Map Var SimplexNum -> Linear.Term.Types.Term -> SimplexNum +evalTerm _ (Linear.Term.Types.ConstTerm c) = c +evalTerm varMap (Linear.Term.Types.CoeffTerm c v) = + c + * Map.findWithDefault + (error $ "evalTerm: " <> show v <> " not found in varMap " <> show varMap) + v + varMap +evalTerm varMap (Linear.Term.Types.VarTerm v) = + Map.findWithDefault + (error $ "evalTerm: " <> show v <> " not found in varMap " <> show varMap) + v + varMap + +evalTermVarsOnly :: Map.Map Var SimplexNum -> TermVarsOnly -> SimplexNum +evalTermVarsOnly varMap terms = evalTerm varMap $ termVarsOnlyToTerm terms + +evalExpr :: Map.Map Var SimplexNum -> Expr -> SimplexNum +evalExpr varMap expr = sum $ map (evalTerm varMap) $ exprToList expr + +evalExprVarsOnly :: Map.Map Var SimplexNum -> ExprVarsOnly -> SimplexNum +evalExprVarsOnly varMap = evalExpr varMap . exprVarsOnlyToExpr + +evalConstraint :: Map.Map Var SimplexNum -> Constraint -> Bool +evalConstraint varMap (Constraint (lhs :<= rhs)) = evalExpr varMap lhs <= evalExpr varMap rhs +evalConstraint varMap (Constraint (lhs :>= rhs)) = evalExpr varMap lhs >= evalExpr varMap rhs +evalConstraint varMap (Constraint (lhs :== rhs)) = evalExpr varMap lhs == evalExpr varMap rhs + +evalSimpleConstraint :: Map.Map Var SimplexNum -> SimpleConstraint -> Bool +evalSimpleConstraint varMap (SimpleConstraint (lhs :<= rhs)) = evalExprVarsOnly varMap lhs <= rhs +evalSimpleConstraint varMap (SimpleConstraint (lhs :>= rhs)) = evalExprVarsOnly varMap lhs >= rhs +evalSimpleConstraint varMap (SimpleConstraint (lhs :== rhs)) = evalExprVarsOnly varMap lhs == rhs + +evalSimpleSystem :: Map.Map Var SimplexNum -> SimpleSystem -> Bool +evalSimpleSystem varMap = all (evalSimpleConstraint varMap) . unSimpleSystem + +genVarMap :: [Var] -> Gen (Map.Map Var SimplexNum) +genVarMap vars = do + varVals <- forM vars $ const arbitrary + pure $ Map.fromList $ zip vars varVals + +isConstExpr :: Expr -> Bool +isConstExpr expr = + let listExpr = exprToList expr + in all + ( \case + ConstTerm _ -> True + _ -> False + ) + listExpr