diff --git a/.github/workflows/haskell-nix.yml b/.github/workflows/haskell-nix.yml new file mode 100644 index 0000000..741d389 --- /dev/null +++ b/.github/workflows/haskell-nix.yml @@ -0,0 +1,53 @@ +name: Nix CI + +on: + push: + workflow_dispatch: + +permissions: + contents: read + +jobs: + cabal: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + + - name: Check Nix flake inputs + uses: DeterminateSystems/flake-checker-action@3164002371bc90729c68af0e24d5aacf20d7c9f6 # v12 + + - name: Install Nix + uses: DeterminateSystems/nix-installer-action@c5a866b6ab867e88becbed4467b93592bce69f8a # v21 + + - name: Enable Nix cache + uses: DeterminateSystems/magic-nix-cache-action@565684385bcd71bad329742eefe8d12f2e765b39 # v13 + + - name: Build and test (cabal) + run: | + nix develop --command bash -c ' + # hsc2hs-generated test binaries crash in nix shells on Linux + # ("stack smashing detected") due to hardening flags. + # Replace hsc2hs with a wrapper that adds --cross-compile to + # avoid running those binaries entirely. + real=$(which hsc2hs) + mkdir -p /tmp/bin + printf "#!/bin/sh\nexec %s --cross-compile \"\$@\"\n" "$real" > /tmp/bin/hsc2hs + chmod +x /tmp/bin/hsc2hs + CABAL_CONFIGURE_FLAGS="--with-hsc2hs=/tmp/bin/hsc2hs" make ci + ' + + stack: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + + - name: Install Nix + uses: DeterminateSystems/nix-installer-action@c5a866b6ab867e88becbed4467b93592bce69f8a # v21 + + - name: Enable Nix cache + uses: DeterminateSystems/magic-nix-cache-action@565684385bcd71bad329742eefe8d12f2e765b39 # v13 + + - name: Build and test (stack) + run: nix develop --command make stack-ci diff --git a/.github/workflows/haskell-stack.yml b/.github/workflows/haskell-stack.yml new file mode 100644 index 0000000..a3f6fc1 --- /dev/null +++ b/.github/workflows/haskell-stack.yml @@ -0,0 +1,58 @@ +name: Haskell CI (Stack) + +on: + push: + workflow_dispatch: + +permissions: + contents: read + +jobs: + stack: + name: Stack LTS 22.44 on ${{ matrix.os }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + + steps: + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + + - uses: haskell-actions/setup@v2 + id: setup + with: + ghc-version: '9.6' + enable-stack: true + + - name: Show toolchain versions + shell: bash + run: | + ghc --numeric-version + stack --numeric-version + + - name: Cache Stack global package DB + uses: actions/cache@v4 + with: + path: ~/.stack + key: ${{ runner.os }}-stack-global-${{ hashFiles('stack.yaml') }}-${{ hashFiles('simplex-method.cabal') }} + restore-keys: | + ${{ runner.os }}-stack-global- + + - name: Cache Stack work directory + uses: actions/cache@v4 + with: + path: .stack-work + key: ${{ runner.os }}-stack-work-${{ hashFiles('stack.yaml') }}-${{ hashFiles('simplex-method.cabal') }}-${{ hashFiles('**/*.hs') }} + restore-keys: | + ${{ runner.os }}-stack-work-${{ hashFiles('stack.yaml') }}-${{ hashFiles('simplex-method.cabal') }}- + ${{ runner.os }}-stack-work- + + - name: Build dependencies + run: stack build --system-ghc --only-dependencies --test + + - name: Build the package + run: stack build --system-ghc + + - name: Run tests + run: stack test --system-ghc diff --git a/.github/workflows/haskell.yml b/.github/workflows/haskell.yml index 92d3748..e10bb7f 100644 --- a/.github/workflows/haskell.yml +++ b/.github/workflows/haskell.yml @@ -2,129 +2,82 @@ name: Haskell CI on: push: - branches: - - '*' - pull_request: - branches: [ "master" ] + workflow_dispatch: permissions: contents: read jobs: fourmolu: - runs-on: ubuntu-latest - + steps: - - uses: actions/checkout@v3 - - uses: haskell-actions/run-fourmolu@v9 + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + + - uses: haskell-actions/run-fourmolu@015e6ed9db2b3f344482792d42b7317c764eb30f # v12 with: - version: "0.14.0.0" + version: "0.15.0.0" + build: name: GHC ${{ matrix.ghc-version }} on ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: - os: [windows-latest, macos-latest, ubuntu-latest] - ghc-version: ['9.6', '9.4', '9.2', '9.0'] + os: [ubuntu-latest, macos-latest, windows-latest] + ghc-version: ['9.12', '9.10', '9.8', '9.6', '9.4', '9.2'] steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 - name: Set up GHC ${{ matrix.ghc-version }} - uses: haskell/actions/setup@v2 + uses: haskell-actions/setup@v2 id: setup with: ghc-version: ${{ matrix.ghc-version }} - enable-stack: true + cabal-version: "3.16.1.0" - - name: Installed minor versions of GHC, Cabal, and Stack + - name: Show toolchain versions 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}" + ghc --numeric-version + cabal --numeric-version + + - name: Check cabal file + run: cabal check + + - name: Update cabal package index + run: cabal update - name: Configure the build 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 .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- + cabal configure --enable-tests --enable-benchmarks --disable-documentation + cabal build --dry-run - - name: Restore %APPDATA%\stack, %LOCALAPPDATA%\Programs\stack cache (Windows) - uses: actions/cache/restore@v3 - id: cache-restore-stack-global-windows - if: runner.os == 'Windows' + - name: Restore cabal store cache + uses: actions/cache/restore@v4 + id: cache-restore with: - 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') }} + path: ${{ steps.setup.outputs.cabal-store }} + key: ${{ runner.os }}-ghc-${{ matrix.ghc-version }}-cabal-${{ hashFiles('**/plan.json') }} restore-keys: | - ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global- + ${{ runner.os }}-ghc-${{ matrix.ghc-version }}-cabal- - name: Build dependencies - run: stack build --only-dependencies - - - name: Build the package - run: stack build + run: cabal build --only-dependencies - - name: Save .stack-work cache - uses: actions/cache/save@v3 - id: cache-save-stack-work - if: steps.cache-restore-stack-work.outputs.cache-hit != 'true' - with: - 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' + - name: Save cabal store cache + uses: actions/cache/save@v4 + if: steps.cache-restore.outputs.cache-hit != 'true' with: - path: ~/.stack - key: ${{ steps.cache-restore-stack-global-unix.outputs.cache-primary-key }} + path: ${{ steps.setup.outputs.cabal-store }} + key: ${{ steps.cache-restore.outputs.cache-primary-key }} - - name: Run tests - run: stack test + - name: Build the package + run: cabal build all - - name: Check cabal file - run: cabal check + - name: Run tests + run: cabal test all - name: Build documentation - run: stack haddock \ No newline at end of file + run: cabal haddock all --disable-documentation diff --git a/.github/workflows/update-flake.yml b/.github/workflows/update-flake.yml new file mode 100644 index 0000000..010b446 --- /dev/null +++ b/.github/workflows/update-flake.yml @@ -0,0 +1,41 @@ +name: Update Nix Flake + +on: + schedule: + # Every 2 weeks on Monday at 06:00 UTC + - cron: "0 6 1,15 * *" + workflow_dispatch: + +permissions: + contents: write + pull-requests: write + +jobs: + update-flake: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + + - name: Install Nix + uses: DeterminateSystems/nix-installer-action@c5a866b6ab867e88becbed4467b93592bce69f8a # v21 + + - name: Enable Nix cache + uses: DeterminateSystems/magic-nix-cache-action@565684385bcd71bad329742eefe8d12f2e765b39 # v13 + + - name: Update flake lockfile + run: nix flake update + + - name: Create Pull Request + uses: peter-evans/create-pull-request@271a8d0340265f705b14b6d32b9829c1cb33d45e # v7.0.8 + with: + commit-message: "chore(nix): update flake lockfile" + title: "chore(nix): update flake lockfile" + body: | + Automated flake lockfile update via `nix flake update`. + + This PR was created by the scheduled update-flake workflow. + Please review CI results before merging. + branch: automation/update-flake-lock + delete-branch: true + labels: dependencies, automated diff --git a/.gitignore b/.gitignore index c8e2cad..14b3fe5 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,6 @@ *~ dist-*/ .vscode/* +.direnv/* +.envrc +cabal.project.local diff --git a/ChangeLog.md b/ChangeLog.md index 5e6fb45..a9ccb1e 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,24 @@ ## Unreleased changes +- **BREAKING CHANGE**: Restructured `VarDomain` type to support upper bounds + - Replaced `NonNegative`, `LowerBound SimplexNum`, and `Unbounded` constructors with + a single `Bounded { lowerBound :: Maybe SimplexNum, upperBound :: Maybe SimplexNum }` record + - Added smart constructors for convenience: `unbounded`, `nonNegative`, `lowerBoundOnly`, + `upperBoundOnly`, and `boundedRange` + - `Bounded Nothing Nothing` is equivalent to `Unbounded` + - `Bounded (Just 0) Nothing` is equivalent to `NonNegative` + - Upper bounds are now supported and automatically added as LEQ constraints +- Added `AddUpperBound` constructor to `VarTransform` for upper bound constraint generation +- Updated `getTransform` to return a list of transforms (can now generate both lower and upper bound transforms) +- Use Hspec for tests +- Add nix flake +- twoPhaseSimplex now takes a VarDomainMap (as the first param) + - You can specify each Var's domain using smart constructors: `nonNegative`, `unbounded`, + `lowerBoundOnly`, `upperBoundOnly`, or `boundedRange` + - If a VarDomain for a Var is undefined, it's assumed to be `unbounded` + - If you want to keep the same behaviour as before (all vars non-negative), use `nonNegative` for all Vars + ## [v0.2.0.0](https://github.com/rasheedja/LPPaver/tree/v0.2.0.0) - Setup CI diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..80ae3f6 --- /dev/null +++ b/Makefile @@ -0,0 +1,95 @@ +HS_FILES := $(shell git ls-files '*.hs') + +# BUILD_TOOL selects the Haskell build tool: cabal (default) or stack. +BUILD_TOOL ?= cabal + +# --- Shared --- + +.PHONY: format +format: + @test -n "$(HS_FILES)" || { echo "No tracked .hs files found"; exit 0; } + fourmolu -i $(HS_FILES) + +.PHONY: format-check +format-check: + @test -n "$(HS_FILES)" || { echo "No tracked .hs files found"; exit 0; } + fourmolu -m check $(HS_FILES) + +# --- Cabal --- + +.PHONY: cabal-check +cabal-check: + cabal check + +.PHONY: cabal-update +cabal-update: + cabal update + +.PHONY: cabal-configure +cabal-configure: + cabal configure --enable-tests --enable-benchmarks --disable-documentation $(CABAL_CONFIGURE_FLAGS) + cabal build --dry-run + +.PHONY: cabal-deps +cabal-deps: + cabal build --only-dependencies + +.PHONY: cabal-build +cabal-build: + cabal build all + +.PHONY: cabal-test +cabal-test: + cabal test all + +.PHONY: cabal-docs +cabal-docs: + cabal haddock all --disable-documentation + +.PHONY: cabal-ci +cabal-ci: format-check cabal-check cabal-update cabal-configure cabal-deps cabal-build cabal-test cabal-docs + +# --- Stack --- + +.PHONY: stack-update +stack-update: + stack update + +.PHONY: stack-deps +stack-deps: + stack build --only-dependencies --test --no-run-tests + +.PHONY: stack-build +stack-build: + stack build + +.PHONY: stack-test +stack-test: + stack test + +.PHONY: stack-docs +stack-docs: + stack haddock --no-haddock-deps + +.PHONY: stack-ci +stack-ci: format-check stack-update stack-deps stack-build stack-test stack-docs + +# ── Generic (delegates to BUILD_TOOL) ──────────────────────────────── + +.PHONY: update +update: $(BUILD_TOOL)-update + +.PHONY: deps +deps: $(BUILD_TOOL)-deps + +.PHONY: build +build: $(BUILD_TOOL)-build + +.PHONY: test +test: $(BUILD_TOOL)-test + +.PHONY: docs +docs: $(BUILD_TOOL)-docs + +.PHONY: ci +ci: $(BUILD_TOOL)-ci diff --git a/flake.lock b/flake.lock new file mode 100644 index 0000000..bd8d29a --- /dev/null +++ b/flake.lock @@ -0,0 +1,61 @@ +{ + "nodes": { + "flake-utils": { + "inputs": { + "systems": "systems" + }, + "locked": { + "lastModified": 1731533236, + "narHash": "sha256-l0KFg5HjrsfsO/JpG+r7fRrqm12kzFHyUHqHCVpMMbI=", + "owner": "numtide", + "repo": "flake-utils", + "rev": "11707dc2f618dd54ca8739b309ec4fc024de578b", + "type": "github" + }, + "original": { + "owner": "numtide", + "repo": "flake-utils", + "type": "github" + } + }, + "nixpkgs": { + "locked": { + "lastModified": 1775002709, + "narHash": "sha256-d3Yx83vSrN+2z/loBh4mJpyRqr9aAJqlke4TkpFmRJA=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "bcd464ccd2a1a7cd09aa2f8d4ffba83b761b1d0e", + "type": "github" + }, + "original": { + "owner": "NixOS", + "ref": "nixos-25.11", + "repo": "nixpkgs", + "type": "github" + } + }, + "root": { + "inputs": { + "flake-utils": "flake-utils", + "nixpkgs": "nixpkgs" + } + }, + "systems": { + "locked": { + "lastModified": 1681028828, + "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", + "owner": "nix-systems", + "repo": "default", + "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", + "type": "github" + }, + "original": { + "owner": "nix-systems", + "repo": "default", + "type": "github" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..a352de6 --- /dev/null +++ b/flake.nix @@ -0,0 +1,38 @@ +{ + inputs = { + nixpkgs.url = "github:NixOS/nixpkgs/nixos-25.11"; + flake-utils.url = "github:numtide/flake-utils"; + }; + + outputs = { self, nixpkgs, flake-utils }: + flake-utils.lib.eachDefaultSystem (system: + let + pkgs = import nixpkgs { inherit system; }; + hsPkgs = pkgs.haskell.packages.ghc96; + in { + devShells.default = pkgs.mkShell { + buildInputs = [ + # Haskell Tools + hsPkgs.ghc + pkgs.cabal-install + pkgs.stack + hsPkgs.fourmolu + hsPkgs.hlint + hsPkgs.haskell-language-server + # C libraries needed by Haskell deps and GHC runtime + pkgs.zlib + pkgs.zstd + pkgs.xz + pkgs.bzip2 + pkgs.libffi + pkgs.gmp + pkgs.ncurses + # Other tools + pkgs.bash # Workaround for copilot breaking in nix+direnv shells + pkgs.git + pkgs.gh + pkgs.jujutsu + ]; + }; + }); +} diff --git a/package.yaml b/package.yaml deleted file mode 100644 index f4ba4bc..0000000 --- a/package.yaml +++ /dev/null @@ -1,56 +0,0 @@ -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 >= 1.12.2 && < 1.13 - -default-extensions: - DataKinds - DeriveFunctor - DeriveGeneric - DisambiguateRecordFields - DuplicateRecordFields - FlexibleContexts - LambdaCase - OverloadedLabels - OverloadedRecordDot - OverloadedStrings - RecordWildCards - TemplateHaskell - TupleSections - TypeApplications - NamedFieldPuns - -library: - source-dirs: src - -tests: - simplex-haskell-test: - main: Spec.hs - source-dirs: test - dependencies: - - simplex-method diff --git a/simplex-method.cabal b/simplex-method.cabal index 3078198..5e4e4c8 100644 --- a/simplex-method.cabal +++ b/simplex-method.cabal @@ -1,9 +1,3 @@ -cabal-version: 1.12 - --- This file has been generated from package.yaml by hpack version 0.36.0. --- --- see: https://github.com/sol/hpack - name: simplex-method version: 0.2.0.0 synopsis: Implementation of the two-phase simplex method in exact rational arithmetic @@ -16,6 +10,7 @@ maintainer: jrasheed178@gmail.com copyright: BSD-3 license: BSD3 license-file: LICENSE +cabal-version: 1.12 build-type: Simple extra-source-files: README.md @@ -36,34 +31,78 @@ library hs-source-dirs: src default-extensions: - DataKinds DeriveFunctor DeriveGeneric DisambiguateRecordFields DuplicateRecordFields FlexibleContexts LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications NamedFieldPuns + DataKinds + DeriveFunctor + DeriveGeneric + DerivingStrategies + DisambiguateRecordFields + DuplicateRecordFields + ExtendedDefaultRules + FlexibleContexts + GeneralizedNewtypeDeriving + LambdaCase + NamedFieldPuns + OverloadedLabels + OverloadedRecordDot + OverloadedStrings + RecordWildCards + TemplateHaskell + TupleSections + TypeApplications + QuasiQuotes build-depends: 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 >=1.12.2 && <1.13 + , containers >= 0.6.5.1 && < 0.8 + , generic-lens >= 2.2 && < 2.3 + , lens >= 5.2.2 && < 5.4 + , text >= 2.0.2 && < 2.2 + , time >= 1.12.2 && < 1.15 + , monad-logger >= 0.3.40 && < 0.4 + , QuickCheck >= 2.14 && < 2.17 default-language: Haskell2010 test-suite simplex-haskell-test type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: - TestFunctions + Linear.Simplex.Solver.TwoPhaseSpec + Linear.Simplex.UtilSpec Paths_simplex_method hs-source-dirs: test default-extensions: - DataKinds DeriveFunctor DeriveGeneric DisambiguateRecordFields DuplicateRecordFields FlexibleContexts LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications NamedFieldPuns + DataKinds + DeriveFunctor + DeriveGeneric + DerivingStrategies + DisambiguateRecordFields + DuplicateRecordFields + ExtendedDefaultRules + FlexibleContexts + GeneralizedNewtypeDeriving + LambdaCase + NamedFieldPuns + OverloadedLabels + OverloadedRecordDot + OverloadedStrings + RecordWildCards + TemplateHaskell + TupleSections + TypeApplications + QuasiQuotes build-depends: 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 , simplex-method - , text >=2.0.2 && <2.1 - , time >=1.12.2 && <1.13 + , containers >= 0.6.5.1 && < 0.8 + , generic-lens >= 2.2 && < 2.3 + , lens >= 5.2.2 && < 5.4 + , text >= 2.0.2 && < 2.2 + , time >= 1.12.2 && < 1.15 + , monad-logger >= 0.3.40 && < 0.4 + , QuickCheck >= 2.14 && < 2.17 + , hspec >= 2.11.12 && < 2.12 + , hspec-expectations >= 0.8.3 && < 0.9 + , interpolatedstring-perl6 >= 1.0.2 && < 1.1 + build-tool-depends: + hspec-discover:hspec-discover >= 2.11.12 && < 2.12 default-language: Haskell2010 diff --git a/src/Linear/Simplex/Prettify.hs b/src/Linear/Simplex/Prettify.hs index b19cc44..adb62ec 100644 --- a/src/Linear/Simplex/Prettify.hs +++ b/src/Linear/Simplex/Prettify.hs @@ -12,15 +12,14 @@ -- Converts "Linear.Simplex.Types" types into human-readable 'String's module Linear.Simplex.Prettify where -import Control.Lens import Data.Generics.Labels () import Data.Map qualified as M -import Data.Ratio -import Linear.Simplex.Types +import Data.Ratio (denominator, numerator) +import Linear.Simplex.Types (ObjectiveFunction (..), PolyConstraint (..), VarLitMapSum) --- | Convert a 'VarConstMap' into a human-readable 'String' -prettyShowVarConstMap :: VarLitMapSum -> String -prettyShowVarConstMap = aux . M.toList +-- | Convert a 'VarLitMapSum' into a human-readable 'String' +prettyShowVarLitMapSum :: VarLitMapSum -> String +prettyShowVarLitMapSum = aux . M.toList where aux [] = "" aux ((vName, vCoeff) : vs) = prettyShowRational vCoeff ++ " * " ++ show vName ++ " + " ++ aux vs @@ -34,11 +33,11 @@ prettyShowVarConstMap = aux . M.toList -- | 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 (Linear.Simplex.Types.EQ vcm r) = prettyShowVarConstMap vcm ++ " == " ++ show r +prettyShowPolyConstraint (LEQ vcm r) = prettyShowVarLitMapSum vcm ++ " <= " ++ show r +prettyShowPolyConstraint (GEQ vcm r) = prettyShowVarLitMapSum vcm ++ " >= " ++ show r +prettyShowPolyConstraint (Linear.Simplex.Types.EQ vcm r) = prettyShowVarLitMapSum vcm ++ " == " ++ show r -- | Convert an 'ObjectiveFunction' into a human-readable 'String' prettyShowObjectiveFunction :: ObjectiveFunction -> String -prettyShowObjectiveFunction (Min vcm) = "min: " ++ prettyShowVarConstMap vcm -prettyShowObjectiveFunction (Max vcm) = "max: " ++ prettyShowVarConstMap vcm +prettyShowObjectiveFunction (Min vcm) = "min: " ++ prettyShowVarLitMapSum vcm +prettyShowObjectiveFunction (Max vcm) = "max: " ++ prettyShowVarLitMapSum vcm diff --git a/src/Linear/Simplex/Solver/TwoPhase.hs b/src/Linear/Simplex/Solver/TwoPhase.hs index c7dfe83..713fba0 100644 --- a/src/Linear/Simplex/Solver/TwoPhase.hs +++ b/src/Linear/Simplex/Solver/TwoPhase.hs @@ -6,27 +6,79 @@ -- Maintainer : jrasheed178@gmail.com -- Stability : experimental -- --- Module implementing the two-phase simplex method. +-- | 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 +-- 'twoPhaseSimplex' supports variable domains via its 'VarDomainMap' argument. +module Linear.Simplex.Solver.TwoPhase + ( findFeasibleSolution + , optimizeFeasibleSystem + , twoPhaseSimplex + -- Internal functions exported for testing + , preprocess + , postprocess + , computeObjective + , collectAllVars + , generateTransform + , getTransform + , applyTransforms + , applyTransform + , applyShiftToObjective + , applyShiftToConstraint + , applySplitToObjective + , applySplitToConstraint + , unapplyTransformsToVarMap + , unapplyTransformToVarMap + ) where import Prelude hiding (EQ) -import Control.Lens +import Control.Lens ((%~), (&), (.~), (<&>)) import Control.Monad (unless) import Control.Monad.IO.Class (MonadIO) -import Control.Monad.Logger -import Data.Bifunctor -import Data.List +import Control.Monad.Logger (LogLevel (LevelError, LevelInfo, LevelWarn), MonadLogger) +import Data.Bifunctor (second) import qualified Data.Map as M import Data.Maybe (fromJust, fromMaybe, mapMaybe) import Data.Ratio (denominator, numerator, (%)) +import Data.Set (Set) +import qualified Data.Set as Set import qualified Data.Text as Text import GHC.Real (Ratio) import Linear.Simplex.Types + ( Dict + , DictValue (..) + , FeasibleSystem (..) + , ObjectiveFunction (..) + , ObjectiveResult (..) + , OptimisationOutcome (..) + , PivotObjective (..) + , PolyConstraint (..) + , SimplexNum + , SimplexResult (..) + , Tableau + , TableauRow (..) + , Var + , VarDomain (..) + , VarDomainMap (..) + , VarLitMap + , VarLitMapSum + , VarTransform (..) + , nonNegative + , unbounded + ) import Linear.Simplex.Util + ( combineVarLitMapSums + , dictionaryFormToTableau + , foldVarLitMap + , insertPivotObjectiveToDict + , isMax + , logMsg + , showT + , simplifySystem + , tableauInDictionaryForm + ) -- | Find a feasible solution for the given system of 'PolyConstraint's by performing the first phase of the two-phase simplex method -- All variables in the 'PolyConstraint' must be positive. @@ -124,14 +176,17 @@ findFeasibleSolution unsimplifiedSystem = do system = simplifySystem unsimplifiedSystem maxVar = - 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 + if null system + then 0 + else + 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 [] @@ -203,7 +258,6 @@ findFeasibleSolution unsimplifiedSystem = do , 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 @@ -250,34 +304,41 @@ findFeasibleSolution unsimplifiedSystem = do -- | 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) +-- Then, the feasible system in 'Dict' form as well as a list of slack variables, a list artificial variables, and the objective variable. +-- Returns 'Optimal' with variable values if an optimal solution is found, or 'Unbounded' if the objective is unbounded. +optimizeFeasibleSystem :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> FeasibleSystem -> m OptimisationOutcome 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 + mResult <- + 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 + simplexPivot normalObjective phase1Dict + else do + logMsg LevelInfo $ + "optimizeFeasibleSystem: Artificial vars present. Pivoting system (in dict form) " + <> showT phase1Dict + <> " with objective " + <> showT adjustedObjective + simplexPivot adjustedObjective phase1Dict + case mResult of + Nothing -> do + logMsg LevelInfo "optimizeFeasibleSystem: Objective is unbounded (ratio test failed)" + pure Unbounded + Just resultDict -> do + let result = displayResults (dictionaryFormToTableau resultDict) + logMsg LevelInfo $ "optimizeFeasibleSystem: Found optimal solution: " <> showT result + pure result where - -- \| displayResults takes a 'Tableau' and returns a 'Result'. The 'Tableau' + -- \| displayResults takes a 'Tableau' and returns an 'OptimisationOutcome'. 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. + -- algorithm has been applied. The 'OptimisationOutcome' contains the values of all + -- variables appearing in the system. -- -- The function first filters out the rows of the tableau that correspond -- to the slack and artificial variables. It then extracts the values of @@ -287,12 +348,9 @@ optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..}) -- 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 -> OptimisationOutcome displayResults tableau = - Result - { objectiveVar = objectiveVar - , varValMap = extractVarVals - } + Optimal extractVarVals where extractVarVals = let tableauWithOriginalVars = @@ -376,32 +434,292 @@ optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..}) ) (M.toList objFunction.objective) --- | 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 :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> [PolyConstraint] -> m (Maybe Result) -twoPhaseSimplex objFunction unsimplifiedSystem = do +-- | Perform the two phase simplex method with variable domain information. +-- Variables not in the VarDomainMap are assumed to be Unbounded (no lower bound). +-- This function applies necessary transformations before solving and unapplies them after. +-- The returned SimplexResult contains: +-- * The feasible system (Nothing if infeasible) +-- * Results for each objective function (empty if infeasible) +twoPhaseSimplex :: + (MonadIO m, MonadLogger m) => VarDomainMap -> [ObjectiveFunction] -> [PolyConstraint] -> m SimplexResult +twoPhaseSimplex domainMap objFunctions constraints = do + logMsg LevelInfo $ + "twoPhaseSimplex: Solving system with domain map " <> showT domainMap + -- Collect original variables before any transformations + let originalVars = collectAllVars objFunctions constraints + let (transformedObjs, transformedConstraints, transforms) = preprocess objFunctions domainMap constraints logMsg LevelInfo $ - "twoPhaseSimplex: Solving system " <> showT unsimplifiedSystem <> " with objective " <> showT objFunction - phase1Result <- findFeasibleSolution unsimplifiedSystem - case phase1Result of + "twoPhaseSimplex: Applied transforms " + <> showT transforms + <> "; Transformed objectives: " + <> showT transformedObjs + <> "; Transformed constraints: " + <> showT transformedConstraints + mFeasibleSystem <- findFeasibleSolution transformedConstraints + case mFeasibleSystem of + Nothing -> do + logMsg LevelInfo "twoPhaseSimplex: No feasible solution found in phase 1" + pure $ SimplexResult Nothing [] Just feasibleSystem -> do logMsg LevelInfo $ - "twoPhaseSimplex: Feasible system found for " - <> showT unsimplifiedSystem - <> "; Feasible system: " + "twoPhaseSimplex: Feasible system found for transformed system; 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 + objResults <- optimizeAllObjectives originalVars transforms feasibleSystem (zip objFunctions transformedObjs) + logMsg LevelInfo $ "twoPhaseSimplex: All objective results: " <> showT objResults + pure $ SimplexResult (Just feasibleSystem) objResults + +-- | Optimize all objective functions over the given feasible system. +-- Returns a list of ObjectiveResult, one for each objective function. +-- The originalVars set is used to filter the result to only include original decision variables. +optimizeAllObjectives :: + (MonadIO m, MonadLogger m) => + -- | Original decision variables + Set.Set Var -> + [VarTransform] -> + FeasibleSystem -> + -- | (original, transformed) objective pairs + [(ObjectiveFunction, ObjectiveFunction)] -> + m [ObjectiveResult] +optimizeAllObjectives originalVars transforms feasibleSystem objPairs = + mapM optimizeOne objPairs + where + optimizeOne (origObj, transformedObj) = do + outcome <- optimizeFeasibleSystem transformedObj feasibleSystem + let postprocessedOutcome = postprocess originalVars transforms outcome + pure $ ObjectiveResult origObj postprocessedOutcome + +-- | Postprocess the optimisation outcome by unapplying variable transformations +-- and filtering to only include the original decision variables. +-- For Optimal outcomes, unapplies transforms to get variable values in original space. +-- For Unbounded outcomes, passes through unchanged. +postprocess :: Set.Set Var -> [VarTransform] -> OptimisationOutcome -> OptimisationOutcome +postprocess _originalVars _transforms Unbounded = Unbounded +postprocess originalVars transforms (Optimal varVals) = + let -- Unapply transforms to get variable values in original space + unappliedVarVals = unapplyTransformsToVarMap transforms varVals + -- Filter to only include original decision variables + filteredVarVals = M.filterWithKey (\k _ -> Set.member k originalVars) unappliedVarVals + in Optimal filteredVarVals + +-- | Compute the value of an objective function given variable values. +computeObjective :: ObjectiveFunction -> M.Map Var SimplexNum -> SimplexNum +computeObjective objFunction varVals = + let coeffs = case objFunction of + Max m -> m + Min m -> m + in sum $ map (\(var, coeff) -> coeff * M.findWithDefault 0 var varVals) (M.toList coeffs) + +-- | Preprocess the system by applying variable transformations based on domain information. +-- Returns the transformed objectives, constraints, and the list of transforms applied. +preprocess :: + [ObjectiveFunction] -> + VarDomainMap -> + [PolyConstraint] -> + ([ObjectiveFunction], [PolyConstraint], [VarTransform]) +preprocess objFunctions (VarDomainMap domainMap) constraints = + let -- Collect all variables in the system (from all objectives and constraints) + allVars = collectAllVars objFunctions constraints + -- Find the maximum variable to generate fresh variables + maxVar = if Set.null allVars then 0 else Set.findMax allVars + -- Generate transforms for each variable based on its domain + -- Variables not in domainMap are treated as Unbounded + (transforms, _) = foldr (generateTransform domainMap) ([], maxVar) (Set.toList allVars) + -- Apply transforms to get the transformed system + transformedObjs = map (\obj -> fst $ applyTransforms transforms obj constraints) objFunctions + (_, transformedConstraints) = case objFunctions of + [] -> (Max M.empty, applyTransformsToConstraints transforms constraints) + (obj : _) -> applyTransforms transforms obj constraints + in (transformedObjs, transformedConstraints, transforms) + +-- | Apply transforms to constraints only (used when there are no objectives) +applyTransformsToConstraints :: [VarTransform] -> [PolyConstraint] -> [PolyConstraint] +applyTransformsToConstraints transforms constraints = + snd $ applyTransforms transforms (Max M.empty) constraints + +-- | Collect all variables appearing in the objective functions and constraints +collectAllVars :: [ObjectiveFunction] -> [PolyConstraint] -> Set Var +collectAllVars objFunctions constraints = + let objVars = Set.unions $ map getObjVars objFunctions + constraintVars = Set.unions $ map getConstraintVars constraints + in Set.union objVars constraintVars + where + getObjVars :: ObjectiveFunction -> Set Var + getObjVars (Max m) = M.keysSet m + getObjVars (Min m) = M.keysSet m + + getConstraintVars :: PolyConstraint -> Set Var + getConstraintVars (LEQ m _) = M.keysSet m + getConstraintVars (GEQ m _) = M.keysSet m + getConstraintVars (EQ m _) = M.keysSet m + +-- | Generate a transform for a variable based on its domain. +-- Takes the domain map, the variable, and the current (transforms, nextFreshVar). +-- Returns updated (transforms, nextFreshVar). +generateTransform :: M.Map Var VarDomain -> Var -> ([VarTransform], Var) -> ([VarTransform], Var) +generateTransform domainMap var (transforms, nextFreshVar) = + let domain = M.findWithDefault unbounded var domainMap + (newTransforms, varOffset) = getTransform nextFreshVar var domain + in (newTransforms ++ transforms, nextFreshVar + varOffset) + +-- | Determine what transforms are needed for a variable given its domain. +-- Returns a list of transforms and the number of fresh variables consumed. +getTransform :: Var -> Var -> VarDomain -> ([VarTransform], Var) +getTransform nextFreshVar var (Bounded mLower mUpper) = + let -- Handle lower bound + (lowerTransforms, varOffset) = case mLower of + Nothing -> ([], 0) -- No lower bound: will need Split + Just l + | l == 0 -> ([], 0) -- NonNegative: no transform needed + | l > 0 -> ([AddLowerBound var l], 0) -- Positive lower bound: add constraint + | otherwise -> ([Shift var nextFreshVar l], 1) -- Negative lower bound: shift + + -- Handle upper bound (if present) + upperTransforms = case mUpper of + Nothing -> [] + Just u -> [AddUpperBound var u] + + -- If no lower bound (Nothing), we need Split transformation + -- Split replaces the variable, so upper bound would apply to the original var + -- which gets expressed as posVar - negVar + (finalTransforms, finalOffset) = case mLower of + Nothing -> + -- Unbounded: split the variable + -- Note: upperTransforms will still be added and will apply to the original variable + -- expression (posVar - negVar) via the constraint system + (Split var nextFreshVar (nextFreshVar + 1) : upperTransforms, 2) + Just _ -> + (lowerTransforms ++ upperTransforms, varOffset) + in (finalTransforms, finalOffset) + +-- | Apply all transforms to the objective function and constraints. +applyTransforms :: [VarTransform] -> ObjectiveFunction -> [PolyConstraint] -> (ObjectiveFunction, [PolyConstraint]) +applyTransforms transforms objFunction constraints = + foldr applyTransform (objFunction, constraints) transforms + +-- | Apply a single transform to the objective function and constraints. +applyTransform :: VarTransform -> (ObjectiveFunction, [PolyConstraint]) -> (ObjectiveFunction, [PolyConstraint]) +applyTransform transform (objFunction, constraints) = + case transform of + -- AddLowerBound: Add a GEQ constraint for the variable + AddLowerBound v bound -> + (objFunction, GEQ (M.singleton v 1) bound : constraints) + -- AddUpperBound: Add a LEQ constraint for the variable + AddUpperBound v bound -> + (objFunction, LEQ (M.singleton v 1) bound : constraints) + -- Shift: originalVar = shiftedVar + shiftBy (where shiftBy < 0) + -- Substitute: wherever we see originalVar, replace with shiftedVar + -- and adjust the RHS by -coeff * shiftBy + Shift origVar shiftedVar shiftBy -> + ( applyShiftToObjective origVar shiftedVar shiftBy objFunction + , map (applyShiftToConstraint origVar shiftedVar shiftBy) constraints + ) + -- Split: originalVar = posVar - negVar + -- Substitute: wherever we see originalVar with coeff c, + -- replace with posVar with coeff c and negVar with coeff -c + Split origVar posVar negVar -> + ( applySplitToObjective origVar posVar negVar objFunction + , map (applySplitToConstraint origVar posVar negVar) constraints + ) + +-- | Apply shift transformation to objective function. +-- originalVar = shiftedVar + shiftBy +-- So coefficient of originalVar becomes coefficient of shiftedVar. +-- The constant term changes but objectives don't have constants that affect optimization. +applyShiftToObjective :: Var -> Var -> SimplexNum -> ObjectiveFunction -> ObjectiveFunction +applyShiftToObjective origVar shiftedVar _shiftBy objFunction = + case objFunction of + Max m -> Max (substituteVar origVar shiftedVar m) + Min m -> Min (substituteVar origVar shiftedVar m) + where + substituteVar :: Var -> Var -> VarLitMapSum -> VarLitMapSum + substituteVar oldVar newVar m = + case M.lookup oldVar m of + Nothing -> m + Just coeff -> M.insert newVar coeff (M.delete oldVar m) + +-- | Apply shift transformation to a constraint. +-- originalVar = shiftedVar + shiftBy +-- For constraint: sum(c_i * x_i) REL rhs +-- If x_j = originalVar with coeff c_j: +-- c_j * originalVar = c_j * (shiftedVar + shiftBy) = c_j * shiftedVar + c_j * shiftBy +-- So new constraint: (replace originalVar with shiftedVar) REL (rhs - c_j * shiftBy) +applyShiftToConstraint :: Var -> Var -> SimplexNum -> PolyConstraint -> PolyConstraint +applyShiftToConstraint origVar shiftedVar shiftBy constraint = + case constraint of + LEQ m rhs -> + let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m + in LEQ newMap (rhs - rhsAdjust) + GEQ m rhs -> + let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m + in GEQ newMap (rhs - rhsAdjust) + EQ m rhs -> + let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m + in EQ newMap (rhs - rhsAdjust) + where + substituteVarInMap :: Var -> Var -> SimplexNum -> VarLitMapSum -> (VarLitMapSum, SimplexNum) + substituteVarInMap oldVar newVar shift m = + case M.lookup oldVar m of + Nothing -> (m, 0) + Just coeff -> (M.insert newVar coeff (M.delete oldVar m), coeff * shift) + +-- | Apply split transformation to objective function. +-- originalVar = posVar - negVar +-- coefficient c of originalVar becomes c for posVar and -c for negVar +applySplitToObjective :: Var -> Var -> Var -> ObjectiveFunction -> ObjectiveFunction +applySplitToObjective origVar posVar negVar objFunction = + case objFunction of + Max m -> Max (splitVar origVar posVar negVar m) + Min m -> Min (splitVar origVar posVar negVar m) + where + splitVar :: Var -> Var -> Var -> VarLitMapSum -> VarLitMapSum + splitVar oldVar pVar nVar m = + case M.lookup oldVar m of + Nothing -> m + Just coeff -> M.insert pVar coeff (M.insert nVar (-coeff) (M.delete oldVar m)) + +-- | Apply split transformation to a constraint. +-- originalVar = posVar - negVar +-- coefficient c of originalVar becomes c for posVar and -c for negVar +applySplitToConstraint :: Var -> Var -> Var -> PolyConstraint -> PolyConstraint +applySplitToConstraint origVar posVar negVar constraint = + case constraint of + LEQ m rhs -> LEQ (splitVarInMap origVar posVar negVar m) rhs + GEQ m rhs -> GEQ (splitVarInMap origVar posVar negVar m) rhs + EQ m rhs -> EQ (splitVarInMap origVar posVar negVar m) rhs + where + splitVarInMap :: Var -> Var -> Var -> VarLitMapSum -> VarLitMapSum + splitVarInMap oldVar pVar nVar m = + case M.lookup oldVar m of + Nothing -> m + Just coeff -> M.insert pVar coeff (M.insert nVar (-coeff) (M.delete oldVar m)) + +-- | Unapply transforms to convert a variable value map back to original variables. +unapplyTransformsToVarMap :: [VarTransform] -> VarLitMap -> VarLitMap +unapplyTransformsToVarMap transforms valMap = + -- Apply transforms in reverse order (since we applied them with foldr) + foldl (flip unapplyTransformToVarMap) valMap transforms + +-- | Unapply a single transform to convert a variable value map back to original variables. +unapplyTransformToVarMap :: VarTransform -> VarLitMap -> VarLitMap +unapplyTransformToVarMap transform valMap = + case transform of + -- AddLowerBound: No variable substitution was done, nothing to unapply + AddLowerBound {} -> valMap + -- AddUpperBound: No variable substitution was done, nothing to unapply + AddUpperBound {} -> valMap + -- Shift: originalVar = shiftedVar + shiftBy + -- So originalVar's value = shiftedVar's value + shiftBy + Shift origVar shiftedVar shiftBy -> + let shiftedVal = M.findWithDefault 0 shiftedVar valMap + origVal = shiftedVal + shiftBy + in M.insert origVar origVal (M.delete shiftedVar valMap) + -- Split: originalVar = posVar - negVar + -- So originalVar's value = posVar's value - negVar's value + Split origVar posVar negVar -> + let posVal = M.findWithDefault 0 posVar valMap + negVal = M.findWithDefault 0 negVar valMap + origVal = posVal - negVal + in M.insert origVar origVal (M.delete posVar (M.delete negVar valMap)) -- | 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) diff --git a/src/Linear/Simplex/Types.hs b/src/Linear/Simplex/Types.hs index 15e5d1f..83b0a9a 100644 --- a/src/Linear/Simplex/Types.hs +++ b/src/Linear/Simplex/Types.hs @@ -7,30 +7,19 @@ -- Stability : experimental module Linear.Simplex.Types where -import Control.Lens import Data.Generics.Labels () import Data.List (sort) import qualified Data.Map as M import GHC.Generics (Generic) +-- | Variable identifier type Var = Int +-- | Numeric type used in this library type SimplexNum = Rational -type SystemRow = PolyConstraint - -type System = [SystemRow] - --- 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] - +-- | A feasible system, typically produced by phase one of +-- the two-phase simplex method. data FeasibleSystem = FeasibleSystem { dict :: Dict , slackVars :: [Var] @@ -39,22 +28,34 @@ data FeasibleSystem = FeasibleSystem } deriving (Show, Read, Eq, Generic) -data Result = Result - { objectiveVar :: Var - , varValMap :: VarLitMap - -- TODO: - -- Maybe VarLitMap - -- , feasible :: Bool - -- , optimisable :: Bool +-- | The outcome of optimizing a single objective function. +data OptimisationOutcome + = -- | An optimal solution was found + Optimal {varValMap :: VarLitMap} + | -- | The objective is unbounded + Unbounded + deriving (Show, Read, Eq, Generic) + +-- | Result for a single objective function optimization. +data ObjectiveResult = ObjectiveResult + { objectiveFunction :: ObjectiveFunction + -- ^ The objective that was optimized + , outcome :: OptimisationOutcome + -- ^ The optimization outcome } deriving (Show, Read, Eq, Generic) -data SimplexMeta = SimplexMeta - { objective :: ObjectiveFunction - , feasibleSystem :: Maybe FeasibleSystem - , optimisedResult :: Maybe Result +-- | Complete result of the two-phase simplex method. +-- Contains feasibility information and results for all requested objectives. +data SimplexResult = SimplexResult + { feasibleSystem :: Maybe FeasibleSystem + -- ^ The feasible system (Nothing if infeasible) + , objectiveResults :: [ObjectiveResult] + -- ^ Results for each objective (empty if infeasible) } + deriving (Show, Read, Eq, Generic) +-- | Mapping from variable id to its numeric value/coefficient. type VarLitMap = M.Map Var SimplexNum -- | List of variables with their 'SimplexNum' coefficients. @@ -64,7 +65,7 @@ type VarLitMap = M.Map Var SimplexNum type VarLitMapSum = VarLitMap -- | For specifying constraints in a system. --- The LHS is a 'Vars', and the RHS, is a 'SimplexNum' number. +-- 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. -- @@ -76,17 +77,10 @@ data PolyConstraint deriving (Show, Read, Eq, Generic) -- | Create an objective function. --- We can either 'Max'imize or 'Min'imize a 'VarTermSum'. +-- We can either 'Max'imize or 'Min'imize a 'VarLitMapSum'. 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 @@ -94,7 +88,7 @@ data TableauRow = TableauRow } deriving (Show, Read, Eq, Generic) --- | A simplex 'Tableu' of equations. +-- | A simplex 'Tableau' of equations. -- Each entry in the map is a row. type Tableau = M.Map Var TableauRow @@ -115,9 +109,84 @@ data DictValue = DictValue -- deriving (Show, Read, Eq, Generic) type Dict = M.Map Var DictValue +-- | Objective row representation used during pivoting. +-- 'variable' is the objective basic variable and 'function'/'constant' encode +-- the objective in dictionary form. data PivotObjective = PivotObjective { variable :: Var , function :: VarLitMapSum , constant :: SimplexNum } deriving (Show, Read, Eq, Generic) + +-- | Domain specification for a variable's bounds. +-- Variables not in the VarDomainMap are assumed to be Unbounded (both bounds Nothing). +-- +-- Bounds semantics: +-- * @lowerBound = Just L@ means var >= L +-- * @lowerBound = Nothing@ means no lower bound (var can be arbitrarily negative) +-- * @upperBound = Just U@ means var <= U +-- * @upperBound = Nothing@ means no upper bound (var can be arbitrarily positive) +-- +-- Note: @Bounded Nothing Nothing@ is equivalent to unbounded. Use the smart constructors +-- ('unbounded', 'nonNegative', etc.) for clarity. +data VarDomain = Bounded + { lowerBound :: Maybe SimplexNum + -- ^ Lower bound (Nothing = -∞) + , upperBound :: Maybe SimplexNum + -- ^ Upper bound (Nothing = +∞) + } + deriving stock (Show, Read, Eq, Generic) + +-- | Smart constructor for an unbounded variable (no lower or upper bound). +-- The variable can take any real value. +unbounded :: VarDomain +unbounded = Bounded Nothing Nothing + +-- | Smart constructor for a non-negative variable (var >= 0). +-- This is the standard simplex assumption. +nonNegative :: VarDomain +nonNegative = Bounded (Just 0) Nothing + +-- | Smart constructor for a variable with only a lower bound (var >= L). +lowerBoundOnly :: SimplexNum -> VarDomain +lowerBoundOnly l = Bounded (Just l) Nothing + +-- | Smart constructor for a variable with only an upper bound (var <= U). +upperBoundOnly :: SimplexNum -> VarDomain +upperBoundOnly u = Bounded Nothing (Just u) + +-- | Smart constructor for a variable with both lower and upper bounds (L <= var <= U). +boundedRange :: SimplexNum -> SimplexNum -> VarDomain +boundedRange l u = Bounded (Just l) (Just u) + +-- | Map from variables to their domain specifications. +-- Variables not in this map are assumed to be Unbounded. +newtype VarDomainMap = VarDomainMap {unVarDomainMap :: M.Map Var VarDomain} + deriving stock (Show, Read, Eq, Generic) + +-- | Transformations applied to variables to ensure they satisfy the non-negativity requirement. +data VarTransform + = -- | var >= bound where bound > 0. Adds GEQ constraint to system. + AddLowerBound + { var :: !Var + , bound :: !SimplexNum + } + | -- | var <= bound. Adds LEQ constraint to system. + AddUpperBound + { var :: !Var + , bound :: !SimplexNum + } + | -- | originalVar = shiftedVar + shiftBy, where shiftBy < 0. After solving: originalVar = shiftedVar + shiftBy + Shift + { originalVar :: !Var + , shiftedVar :: !Var + , shiftBy :: !SimplexNum + } + | -- | originalVar = posVar - negVar, both posVar and negVar >= 0 + Split + { originalVar :: !Var + , posVar :: !Var + , negVar :: !Var + } + deriving stock (Show, Read, Eq, Generic) diff --git a/src/Linear/Simplex/Util.hs b/src/Linear/Simplex/Util.hs index 99b1495..38b415c 100644 --- a/src/Linear/Simplex/Util.hs +++ b/src/Linear/Simplex/Util.hs @@ -9,13 +9,10 @@ -- Helper functions for performing the two-phase simplex method. module Linear.Simplex.Util where -import Control.Lens import Control.Monad.IO.Class (MonadIO (..)) -import Control.Monad.Logger (LogLevel (..), LogLine, MonadLogger, logDebug, logError, logInfo, logWarn) -import Data.Bifunctor +import Control.Monad.Logger (LogLevel (..), MonadLogger, logDebug, logError, logInfo, logWarn) import Data.Generics.Labels () -import Data.Generics.Product (field) -import Data.List +import Data.List (nub, (\\)) import qualified Data.Map as Map import qualified Data.Map.Merge.Lazy as MapMerge import Data.Maybe (fromMaybe) @@ -23,6 +20,18 @@ import qualified Data.Text as T import Data.Time (getCurrentTime) import Data.Time.Format.ISO8601 (iso8601Show) import Linear.Simplex.Types + ( Dict + , DictValue (..) + , ObjectiveFunction (..) + , PivotObjective (..) + , PolyConstraint (..) + , SimplexNum + , Tableau + , TableauRow (..) + , Var + , VarLitMap + , VarLitMapSum + ) import Prelude hiding (EQ) -- | Is the given 'ObjectiveFunction' to be 'Max'imized? @@ -30,7 +39,7 @@ isMax :: ObjectiveFunction -> Bool isMax (Max _) = True isMax (Min _) = False --- | Simplifies a system of 'PolyConstraint's by first calling 'simplifyPolyConstraint', +-- | Simplifies a system of 'PolyConstraint's, -- 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] @@ -78,7 +87,7 @@ simplifySystem = nub . reduceSystem then EQ lhs rhs : reduceSystem pcs else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) --- | Converts a 'Dict' to a 'Tableau' using 'dictEntryToTableauEntry'. +-- | Converts a 'Dict' to a 'Tableau'. -- FIXME: maybe remove this line. The basic variables will have a coefficient of 1 in the 'Tableau'. dictionaryFormToTableau :: Dict -> Tableau dictionaryFormToTableau = @@ -106,15 +115,6 @@ tableauInDictionaryForm = } ) --- | 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" - Just r -> r - -- | Combines two 'VarLitMapSums together by summing values with matching keys combineVarLitMapSums :: VarLitMapSum -> VarLitMapSum -> VarLitMapSum combineVarLitMapSums = @@ -126,17 +126,6 @@ combineVarLitMapSums = 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 @@ -154,7 +143,7 @@ foldVarLitMap (vm1 : vm2 : vms) = (Just vm1VarVal, Just vm2VarVal) -> vm1VarVal + vm2VarVal (Just vm1VarVal, Nothing) -> vm1VarVal (Nothing, Just vm2VarVal) -> vm2VarVal - (Nothing, Nothing) -> error "Reached unreachable branch in foldDictValue" + (Nothing, Nothing) -> error "Reached unreachable branch in foldVarLitMap" ) ) combinedVars @@ -176,9 +165,3 @@ logMsg lvl msg = do 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) diff --git a/stack.yaml b/stack.yaml index eab5650..bfd8f93 100644 --- a/stack.yaml +++ b/stack.yaml @@ -1,9 +1,3 @@ -# 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: @@ -17,7 +11,7 @@ # # resolver: ./custom-snapshot.yaml # resolver: https://example.com/snapshots/2018-01-01.yaml -resolver: lts-21.22 +resolver: lts-22.44 # User packages to be built. # Various formats can be used as shown in the example below. @@ -30,6 +24,7 @@ resolver: lts-21.22 # - 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: @@ -39,7 +34,7 @@ packages: # - git: https://github.com/commercialhaskell/stack.git # commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a # -# extra-deps: {} +extra-deps: # Override default flag values for local packages and extra-deps # flags: {} @@ -64,5 +59,3 @@ packages: # # 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 index e8d3cc7..8d134eb 100644 --- a/stack.yaml.lock +++ b/stack.yaml.lock @@ -1,12 +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 +# https://docs.haskellstack.org/en/stable/topics/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 + sha256: 238fa745b64f91184f9aa518fe04bdde6552533d169b0da5256670df83a0f1a9 + size: 721141 + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/22/44.yaml + original: lts-22.44 diff --git a/test/Linear/Simplex/Solver/TwoPhaseSpec.hs b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs new file mode 100644 index 0000000..9cfdb82 --- /dev/null +++ b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs @@ -0,0 +1,3123 @@ +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Linear.Simplex.Solver.TwoPhaseSpec where + +import Prelude hiding (EQ) + +import Control.Monad.Logger (LogLevel (LevelInfo), filterLogger, runStdoutLoggingT) +import qualified Data.Map as M +import Data.Maybe (isJust) +import Data.Ratio ((%)) +import qualified Data.Set as Set + +import Test.Hspec (Spec, describe, expectationFailure, it, shouldBe, shouldSatisfy) +import Test.QuickCheck (NonEmptyList (..), Positive (..), property, (==>)) + +import Linear.Simplex.Prettify (prettyShowObjectiveFunction, prettyShowPolyConstraint, prettyShowVarLitMapSum) +import Linear.Simplex.Solver.TwoPhase + ( applyShiftToConstraint + , applyShiftToObjective + , applySplitToConstraint + , applySplitToObjective + , applyTransform + , applyTransforms + , collectAllVars + , computeObjective + , generateTransform + , getTransform + , postprocess + , preprocess + , twoPhaseSimplex + , unapplyTransformToVarMap + , unapplyTransformsToVarMap + ) +import Linear.Simplex.Types + ( ObjectiveFunction (..) + , ObjectiveResult (..) + , OptimisationOutcome (..) + , PolyConstraint (..) + , SimplexResult (..) + , VarDomainMap (..) + , VarTransform (..) + , boundedRange + , lowerBoundOnly + , nonNegative + , unbounded + , upperBoundOnly + ) + +spec :: Spec +spec = do + describe "twoPhaseSimplex" $ do + -- From page 50 of 'Linear and Integer Programming Made Easy' + describe "From 'Linear and Integer Programming Made Easy' (page 50)" $ do + it "Max 3x₁ + 5x₂ with LEQ constraints: obj=29, x₁=3, x₂=4" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + constraints = + [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 + , LEQ (M.fromList [(1, 1), (2, 1)]) 7 + , LEQ (M.fromList [(2, 1)]) 4 + , LEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 4)] + computeObjective obj varMap `shouldBe` 29 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min 3x₁ + 5x₂ with LEQ constraints: obj=0" $ do + let obj = Min (M.fromList [(1, 3), (2, 5)]) + constraints = + [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 + , LEQ (M.fromList [(1, 1), (2, 1)]) 7 + , LEQ (M.fromList [(2, 1)]) 4 + , LEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.empty + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max 3x₁ + 5x₂ with GEQ constraints: unbounded" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + constraints = + [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 + , GEQ (M.fromList [(1, 1), (2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) 4 + , GEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + it "Min 3x₁ + 5x₂ with GEQ constraints: obj=237/7, x₁=24/7, x₂=33/7" $ do + let obj = Min (M.fromList [(1, 3), (2, 5)]) + constraints = + [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 + , GEQ (M.fromList [(1, 1), (2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) 4 + , GEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 24 % 7), (2, 33 % 7)] + computeObjective obj varMap `shouldBe` (237 % 7) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- From https://www.eng.uwaterloo.ca/~syde05/phase1.pdf (requires two phases) + describe "From eng.uwaterloo.ca phase1.pdf (requires two phases)" $ do + it "Max x₁ - x₂ + x₃ with LEQ constraints: obj=3/5, x₂=14/5, x₃=17/5" $ do + let obj = Max (M.fromList [(1, 1), (2, -1), (3, 1)]) + constraints = + [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 14 % 5), (3, 17 % 5)] + computeObjective obj varMap `shouldBe` (3 % 5) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ - x₂ + x₃ with LEQ constraints: unbounded" $ do + let obj = Min (M.fromList [(1, 1), (2, -1), (3, 1)]) + constraints = + [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + it "Max x₁ - x₂ + x₃ with GEQ constraints: obj=1, x₁=3, x₂=2" $ do + let obj = Max (M.fromList [(1, 1), (2, -1), (3, 1)]) + constraints = + [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , GEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 2)] + computeObjective obj varMap `shouldBe` 1 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ - x₂ + x₃ with GEQ constraints: obj=-1/4, x₁=17/4, x₂=9/2" $ do + let obj = Min (M.fromList [(1, 1), (2, -1), (3, 1)]) + constraints = + [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , GEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 17 % 4), (2, 9 % 2)] + computeObjective obj varMap `shouldBe` ((-1) % 4) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- From page 49 of 'Linear and Integer Programming Made Easy' (requires two phases) + describe "From 'Linear and Integer Programming Made Easy' (page 49, requires two phases)" $ do + it "Min x₁ + x₂ + 2x₃ + x₄ with EQ constraints: obj=5, x₃=2, x₄=1" $ do + let obj = Min (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 + , EQ (M.fromList [(2, 1), (3, 1), (4, 4)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(3, 2), (4, 1)] + computeObjective obj varMap `shouldBe` 5 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₁ + x₂ + 2x₃ + x₄ with EQ constraints: obj=8, x₁=2, x₂=6" $ do + let obj = Max (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 + , EQ (M.fromList [(2, 1), (3, 1), (4, 4)]) 6 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 2), (2, 6)] + computeObjective obj varMap `shouldBe` 8 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- From page 52 of 'Linear and Integer Programming Made Easy' + describe "From 'Linear and Integer Programming Made Easy' (page 52)" $ do + it "Max -2x₃ + 2x₄ + x₅ with EQ constraints: obj=20, x₃=6, x₄=16" $ do + let obj = Max (M.fromList [(3, -2), (4, 2), (5, 1)]) + constraints = + [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 + , EQ (M.fromList [(3, 3), (4, -1), (5, 2)]) 2 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(3, 6), (4, 16)] + computeObjective obj varMap `shouldBe` 20 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min -2x₃ + 2x₄ + x₅ with EQ constraints: obj=6, x₄=2, x₅=2" $ do + let obj = Min (M.fromList [(3, -2), (4, 2), (5, 1)]) + constraints = + [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 + , EQ (M.fromList [(3, 3), (4, -1), (5, 2)]) 2 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(4, 2), (5, 2)] + computeObjective obj varMap `shouldBe` 6 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- From page 59 of 'Linear and Integer Programming Made Easy' (requires two phases) + describe "From 'Linear and Integer Programming Made Easy' (page 59, requires two phases)" $ do + it "Max 2x₁ + x₂: obj=150, x₂=150" $ do + let obj = Max (M.fromList [(1, 2), (2, 1)]) + constraints = + [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 150)] + computeObjective obj varMap `shouldBe` 150 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min 2x₁ + x₂: obj=40/3, x₂=40/3" $ do + let obj = Min (M.fromList [(1, 2), (2, 1)]) + constraints = + [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 40 % 3)] + computeObjective obj varMap `shouldBe` (40 % 3) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max 2x₁ + x₂ with GEQ constraints: unbounded" $ do + let obj = Max (M.fromList [(1, 2), (2, 1)]) + constraints = + [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + it "Min 2x₁ + x₂ with GEQ constraints: obj=75, x₁=75/2" $ do + let obj = Min (M.fromList [(1, 2), (2, 1)]) + constraints = + [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 75 % 2)] + computeObjective obj varMap `shouldBe` 75 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- From page 59 of 'Linear and Integer Programming Made Easy' + describe "From 'Linear and Integer Programming Made Easy' (page 59)" $ do + it "Min -6x₁ - 4x₂ + 2x₃: obj=-120, x₁=20" $ do + let obj = Min (M.fromList [(1, -6), (2, -4), (3, 2)]) + constraints = + [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , LEQ (M.fromList [(2, -5), (3, 5)]) 100 + , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 20)] + computeObjective obj varMap `shouldBe` (-120) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max -6x₁ - 4x₂ + 2x₃: obj=10, x₃=5" $ do + let obj = Max (M.fromList [(1, -6), (2, -4), (3, 2)]) + constraints = + [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , LEQ (M.fromList [(2, -5), (3, 5)]) 100 + , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(3, 5)] + computeObjective obj varMap `shouldBe` 10 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min -6x₁ - 4x₂ + 2x₃ with GEQ constraints: unbounded" $ do + let obj = Min (M.fromList [(1, -6), (2, -4), (3, 2)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , GEQ (M.fromList [(2, -5), (3, 5)]) 100 + , GEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + it "Max -6x₁ - 4x₂ + 2x₃ with GEQ constraints: unbounded" $ do + let obj = Max (M.fromList [(1, -6), (2, -4), (3, 2)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , GEQ (M.fromList [(2, -5), (3, 5)]) 100 + , GEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + -- From page 59 of 'Linear and Integer Programming Made Easy' + describe "From 'Linear and Integer Programming Made Easy' (page 59)" $ do + it "Max 3x₁ + 5x₂ + 2x₃: obj=250, x₂=50" $ do + let obj = Max (M.fromList [(1, 3), (2, 5), (3, 2)]) + constraints = + [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 50)] + computeObjective obj varMap `shouldBe` 250 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min 3x₁ + 5x₂ + 2x₃: obj=0" $ do + let obj = Min (M.fromList [(1, 3), (2, 5), (3, 2)]) + constraints = + [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.empty + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max 3x₁ + 5x₂ + 2x₃ with GEQ constraints: unbounded" $ do + let obj = Max (M.fromList [(1, 3), (2, 5), (3, 2)]) + constraints = + [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + SimplexResult Nothing _ -> expectationFailure "Expected unbounded but got infeasible" + _ -> expectationFailure "Expected unbounded" + + it "Min 3x₁ + 5x₂ + 2x₃ with GEQ constraints: obj=300, x₃=150" $ do + let obj = Min (M.fromList [(1, 3), (2, 5), (3, 2)]) + constraints = + [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(3, 150)] + computeObjective obj varMap `shouldBe` 300 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "Simple single/two variable tests" $ do + it "Max x₁ with x₁ <= 15: obj=15, x₁=15" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 15 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 15)] + computeObjective obj varMap `shouldBe` 15 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max 2x₁ with mixed constraints: obj=20, x₁=10, x₂=10" $ do + let obj = Max (M.fromList [(1, 2)]) + constraints = + [ LEQ (M.fromList [(1, 2)]) 20 + , GEQ (M.fromList [(2, 1)]) 10 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 10), (2, 10)] + computeObjective obj varMap `shouldBe` 20 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ with x₁ <= 15: obj=0" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 15 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.empty + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min 2x₁ with mixed constraints: obj=0, x₂=10" $ do + let obj = Min (M.fromList [(1, 2)]) + constraints = + [ LEQ (M.fromList [(1, 2)]) 20 + , GEQ (M.fromList [(2, 1)]) 10 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 10)] + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "Infeasibility tests" $ do + it "Conflicting bounds x₁ <= 15 and x₁ >= 15.01: infeasible" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 15 + , GEQ (M.fromList [(1, 1)]) 15.01 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Conflicting bounds with additional constraint: infeasible" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 15 + , GEQ (M.fromList [(1, 1)]) 15.01 + , GEQ (M.fromList [(2, 1)]) 10 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Min x₁ with duplicate GEQ constraints: obj=0, x₂=1" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1)]) 1 + , GEQ (M.fromList [(1, 1), (2, 1)]) 1 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 1 % 1)] + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Conflicting x₁+x₂ >= 2 and x₁+x₂ <= 1: infeasible" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1)]) 2 + , LEQ (M.fromList [(1, 1), (2, 1)]) 1 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + describe "LEQ/GEQ reduction bug tests" $ do + it "testLeqGeqBugMin1: obj=3, x₁=3, x₂=3" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 3)] + computeObjective obj varMap `shouldBe` 3 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "testLeqGeqBugMax1: obj=3, x₁=3, x₂=3" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 3)] + computeObjective obj varMap `shouldBe` 3 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "testLeqGeqBugMin2: obj=3, x₁=3, x₂=3" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 3)] + computeObjective obj varMap `shouldBe` 3 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "testLeqGeqBugMax2: obj=3, x₁=3, x₂=3" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 3), (2, 3)] + computeObjective obj varMap `shouldBe` 3 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + -- PolyPaver-style tests with shared parameters + describe "PolyPaver-style tests (feasible region [0,2.5]²)" $ do + let 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 + mkConstraints obj = + ( obj + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + + it "Min x₁: x₁=7/4, x₂=5/2" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 7 % 4), (2, 5 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (7 % 4) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₁: x₁=5/2, x₂=5/3" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 5 % 2), (2, 5 % 3), (3, 0)] + computeObjective obj varMap `shouldBe` (5 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₂: x₂=5/3" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 5 % 3), (1, 5 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (5 % 3) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₂: x₂=5/2" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 5 % 2), (1, 5 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (5 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "PolyPaver-style tests (infeasible region [0,1.5]²)" $ do + let 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 + mkConstraints obj = + ( obj + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + + it "Max x₁: infeasible" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Min x₁: infeasible" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Max x₂: infeasible" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Min x₂: infeasible" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + describe "PolyPaver-style tests (feasible region [0,3.5]²)" $ do + let 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 + mkConstraints obj = + ( obj + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + + it "Max x₁: x₁=7/2" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 5 % 9), (1, 7 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (7 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁: x₁=17/20" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 17 % 20), (2, 7 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (17 % 20) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₂: x₂=7/2" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 7 % 2), (1, 22 % 9)] + computeObjective obj varMap `shouldBe` (7 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₂: x₂=5/9" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 5 % 9), (1, 7 % 2), (3, 0)] + computeObjective obj varMap `shouldBe` (5 % 9) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "PolyPaver two-function tests (infeasible)" $ do + let 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 + mkConstraints obj = + ( obj + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + + it "Max x₁: infeasible" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Min x₁: infeasible" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Max x₂: infeasible" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + it "Min x₂: infeasible" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible" + + describe "PolyPaver two-function tests (feasible)" $ do + let 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 + mkConstraints obj = + ( obj + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + + it "Max x₁: x₁=5/2" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 5 % 2), (2, 45 % 22), (4, 0)] + computeObjective obj varMap `shouldBe` (5 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁: x₁=45/22" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(1, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 45 % 22), (2, 5 % 2), (4, 0)] + computeObjective obj varMap `shouldBe` (45 % 22) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₂: x₂=5/2" $ do + let (obj, constraints) = mkConstraints (Max (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 5 % 2), (1, 5 % 2), (4, 0)] + computeObjective obj varMap `shouldBe` (5 % 2) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₂: x₂=45/22" $ do + let (obj, constraints) = mkConstraints (Min (M.fromList [(2, 1)])) + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 45 % 22), (1, 5 % 2), (4, 0)] + computeObjective obj varMap `shouldBe` (45 % 22) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "QuickCheck-generated regression tests" $ do + it "testQuickCheck1: obj=-370, x₁=5/3, x₂=26" $ do + let obj = Max (M.fromList [(1, 12), (2, -15)]) + constraints = + [ EQ (M.fromList [(1, 24), (2, -2)]) (-12) + , GEQ (M.fromList [(1, -20), (2, 11)]) (-7) + , GEQ (M.fromList [(1, -28), (2, 5)]) (-8) + , GEQ (M.fromList [(1, 3), (2, 0)]) 5 + , LEQ (M.fromList [(1, -48)]) (-1) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 26), (1, 5 % 3)] + computeObjective obj varMap `shouldBe` (-370) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "testQuickCheck2: obj=-2/9, x₁=14/9, x₂=8/9" $ do + let obj = Max (M.fromList [(1, -3), (2, 5)]) + constraints = + [ LEQ (M.fromList [(1, -6), (2, 6)]) 4 + , LEQ (M.fromList [(1, 1), (2, -4), (3, 3)]) (-2) + , LEQ (M.fromList [(2, 7), (1, -4)]) 0 + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(1, 14 % 9), (2, 8 % 9)] + computeObjective obj varMap `shouldBe` ((-2) % 9) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "testQuickCheck3 (tests objective simplification): obj=-8, x₂=2" $ do + let obj = Min (M.fromList [(2, 0), (2, -4)]) + constraints = + [ GEQ (M.fromList [(1, 5), (2, 4)]) (-4) + , LEQ (M.fromList [(1, -1), (2, -1)]) 2 + , LEQ (M.fromList [(2, 1)]) 2 + , GEQ (M.fromList [(1, -5), (2, -1), (2, 1)]) (-5) + ] + allVars = collectAllVars [obj] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + varMap `shouldBe` M.fromList [(2, 2)] + computeObjective obj varMap `shouldBe` (-8) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "twoPhaseSimplex with empty constraint system" $ do + describe "Single variable with boundedRange" $ do + it "Max x₁ with 0 ≤ x₁ ≤ 10: optimal at x₁=10" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` 10 + computeObjective obj varMap `shouldBe` 10 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ with 0 ≤ x₁ ≤ 10: optimal at x₁=0" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` 0 + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₁ with 5 ≤ x₁ ≤ 15: optimal at x₁=15" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 5 15)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` 15 + computeObjective obj varMap `shouldBe` 15 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ with 5 ≤ x₁ ≤ 15: optimal at x₁=5" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 5 15)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` 5 + computeObjective obj varMap `shouldBe` 5 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max x₁ with -5 ≤ x₁ ≤ 10: optimal at x₁=10" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-5) 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` 10 + computeObjective obj varMap `shouldBe` 10 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ with -5 ≤ x₁ ≤ 10: optimal at x₁=-5" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-5) 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.findWithDefault 0 1 varMap `shouldBe` (-5) + computeObjective obj varMap `shouldBe` (-5) + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "Two variables with boundedRange" $ do + it "Max x₁ + x₂ with 0 ≤ x₁ ≤ 10, 0 ≤ x₂ ≤ 10: optimal at 20" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10), (2, boundedRange 0 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 10 + M.lookup 2 varMap `shouldBe` Just 10 + computeObjective obj varMap `shouldBe` 20 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Min x₁ + x₂ with 0 ≤ x₁ ≤ 10, 0 ≤ x₂ ≤ 10: optimal at 0" $ do + let obj = Min (M.fromList [(1, 1), (2, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10), (2, boundedRange 0 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + computeObjective obj varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + it "Max 2x₁ - x₂ with 0 ≤ x₁ ≤ 10, 0 ≤ x₂ ≤ 5: optimal at 20" $ do + let obj = Max (M.fromList [(1, 2), (2, -1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10), (2, boundedRange 0 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 10 + M.findWithDefault 0 2 varMap `shouldBe` 0 + computeObjective obj varMap `shouldBe` 20 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "NonNegative only (no upper bound), empty constraints" $ do + it "Max x₁ with x₁ ≥ 0 and no constraints: unbounded" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + _ -> expectationFailure "Expected unbounded" + + it "Min x₁ with x₁ ≥ 0 and no constraints: optimal at 0" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_s l -> l > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + M.findWithDefault 0 1 varMap `shouldBe` 0 + SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible" + _ -> expectationFailure "Unexpected result" + + describe "twoPhaseSimplex (with VarDomainMap)" $ do + it "Shift transformation with negative lower bound" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "Shift transformation finds minimum at negative bound" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 0] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-5) + _ -> expectationFailure "Unexpected result format" + + it "Split transformation for unbounded variable (max)" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 10 + , GEQ (M.fromList [(1, 1)]) (-10) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "Split transformation for unbounded variable (min)" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 10 + , GEQ (M.fromList [(1, 1)]) (-10) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-10) + _ -> expectationFailure "Unexpected result format" + + it "AddLowerBound with positive lower bound" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "AddLowerBound finds minimum at positive bound" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 5 + _ -> expectationFailure "Unexpected result format" + + it "Mixed domain types" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1), (2, 1)]) 5 + , GEQ (M.fromList [(2, 1)]) (-3) + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-2)), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let xVal = M.findWithDefault 0 1 varMap + yVal = M.findWithDefault 0 2 varMap + oVal = computeObjective obj varMap + (xVal + yVal) `shouldBe` 5 + oVal `shouldBe` 5 + _ -> expectationFailure "Unexpected result format" + + it "lowerBoundOnly 0 is equivalent to NonNegative" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + constraints = + [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 + , LEQ (M.fromList [(1, 1), (2, 1)]) 7 + , LEQ (M.fromList [(2, 1)]) 4 + , LEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + domainMap1 = VarDomainMap $ M.fromList [(1, lowerBoundOnly 0), (2, lowerBoundOnly 0)] + domainMap2 = VarDomainMap $ M.fromList [(1, nonNegative), (2, nonNegative)] + actualResult1 <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap1 [obj] constraints + actualResult2 <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap2 [obj] constraints + -- Both should produce the same optimal solution with x₁=3, x₂=4 + case (actualResult1, actualResult2) of + ( SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap1)] + , SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap2)] + ) -> do + varMap1 `shouldBe` M.fromList [(1, 3), (2, 4)] + varMap1 `shouldBe` varMap2 + _ -> expectationFailure "Expected optimal results" + + it "Infeasible system with domain constraint" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 5] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible result" + + describe "twoPhaseSimplex with upper bounds (AddUpperBound transformation)" $ do + describe "Simple single variable systems" $ do + it "Max x₁ with x₁ ≥ 0, x₁ ≤ 5 (using boundedRange): optimal at x₁=5" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 5 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ ≥ 0, x₁ ≤ 10 (using boundedRange): optimal at x₁=0" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + -- Note: non-basic variables with value 0 may not appear in varValMap + M.findWithDefault 0 1 varMap `shouldBe` 0 + _ -> expectationFailure "Unexpected result format" + + it "Max x₁ with -5 ≤ x₁ ≤ 10 (bounded range with negative lower): optimal at x₁=10" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-5) 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with -5 ≤ x₁ ≤ 10 (bounded range with negative lower): optimal at x₁=-5" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-5) 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-5) + _ -> expectationFailure "Unexpected result format" + + it "Infeasible: lower bound > upper bound" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 10 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible system" + + describe "Two variable systems with upper bounds" $ do + it "Max x₁ + x₂ with 0 ≤ x₁ ≤ 3, 0 ≤ x₂ ≤ 4: optimal at x₁=3, x₂=4" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange 0 3), (2, boundedRange 0 4)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 3 + M.lookup 2 varMap `shouldBe` Just 4 + computeObjective obj varMap `shouldBe` 7 + _ -> expectationFailure "Unexpected result format" + + it "Max 2x₁ - x₂ with -2 ≤ x₁ ≤ 5, -3 ≤ x₂ ≤ 4" $ do + -- Maximize 2x₁ - x₂: want x₁ = 5 (max), x₂ = -3 (min) + -- Optimal: 2*5 - (-3) = 13 + let obj = Max (M.fromList [(1, 2), (2, -1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-2) 5), (2, boundedRange (-3) 4)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 5 + M.lookup 2 varMap `shouldBe` Just (-3) + computeObjective obj varMap `shouldBe` 13 + _ -> expectationFailure "Unexpected result format" + + it "Mixed bounds: x₁ nonNegative, x₂ with upper bound only (unbounded below)" $ do + -- x₁ ≥ 0, x₂ ≤ 10 (no lower bound) + -- Max x₁ + x₂ with x₁ + x₂ ≤ 20 + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 20] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative), (2, upperBoundOnly 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + x1 `shouldSatisfy` (>= 0) + x2 `shouldSatisfy` (<= 10) + (x1 + x2) `shouldBe` 20 + _ -> expectationFailure "Unexpected result format" + + describe "twoPhaseSimplex with negative lower bounds (Shift transformation)" $ do + describe "Simple single variable systems" $ do + it "Max x₁ with x₁ ≤ 5, x₁ ≥ -3: optimal at upper bound x₁=5" $ do + -- Simple case: maximize x with upper bound 5 and lower bound -3 + -- Optimal should be at x₁ = 5 + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 5] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 5 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ ≤ 5, x₁ ≥ -3: optimal at lower bound x₁=-3" $ do + -- Minimize x with upper bound 5 and lower bound -3 + -- Optimal should be at x₁ = -3 + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 5] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-3) + _ -> expectationFailure "Unexpected result format" + + it "Max x₁ with x₁ ≥ -10, x₁ ≤ -2: optimal at x₁=-2" $ do + -- Both bounds are negative, maximize + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) (-2)] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-10))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-2) + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ ≥ -10, x₁ ≤ -2: optimal at x₁=-10" $ do + -- Both bounds are negative, minimize + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) (-2)] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-10))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-10) + _ -> expectationFailure "Unexpected result format" + + describe "Two variable systems with negative bounds" $ do + it "Max x₁ + x₂ with x₁ ≥ -2, x₂ ≥ -3, x₁ + x₂ ≤ 10" $ do + -- Maximize sum, both can go up to contribute to sum ≤ 10 + -- With shifts: x₁' = x₁ + 2, x₂' = x₂ + 3 + -- Constraint becomes: x₁' + x₂' ≤ 15 + -- Optimal in transformed space: x₁' + x₂' = 15 + -- After unapply: x₁ + x₂ = 15 - 5 = 10 + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-2)), (2, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + objVal = computeObjective obj varMap + -- Verify the actual objective value + objVal `shouldBe` 10 + -- Verify lower bounds are respected + x1 `shouldSatisfy` (>= (-2)) + x2 `shouldSatisfy` (>= (-3)) + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ + x₂ with x₁ ≥ -2, x₂ ≥ -3, x₁ + x₂ ≤ 10" $ do + -- Minimize sum with lower bounds -2 and -3 + -- Optimal: x₁ = -2, x₂ = -3, sum = -5 + let obj = Min (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-2)), (2, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let objVal = computeObjective obj varMap + -- Verify the actual objective value + objVal `shouldBe` (-5) + M.lookup 1 varMap `shouldBe` Just (-2) + M.lookup 2 varMap `shouldBe` Just (-3) + _ -> expectationFailure "Unexpected result format" + + it "Max 2x₁ - x₂ with x₁ ≥ -5, x₂ ≥ -4, x₁ ≤ 3, x₂ ≤ 6" $ do + -- Maximize 2x₁ - x₂: want x₁ large (up to 3) and x₂ small (down to -4) + -- Optimal: x₁ = 3, x₂ = -4, obj = 2*3 - (-4) = 10 + let obj = Max (M.fromList [(1, 2), (2, -1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 6 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5)), (2, lowerBoundOnly (-4))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + M.lookup 1 varMap `shouldBe` Just 3 + M.lookup 2 varMap `shouldBe` Just (-4) + -- Verify objective value computed from variables + (2 * x1 - x2) `shouldBe` 10 + _ -> expectationFailure "Unexpected result format" + + it "Min 2x₁ - x₂ with x₁ ≥ -5, x₂ ≥ -4, x₁ ≤ 3, x₂ ≤ 6" $ do + -- Minimize 2x₁ - x₂: want x₁ small (down to -5) and x₂ large (up to 6) + -- Optimal: x₁ = -5, x₂ = 6, obj = 2*(-5) - 6 = -16 + let obj = Min (M.fromList [(1, 2), (2, -1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 6 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5)), (2, lowerBoundOnly (-4))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + M.lookup 1 varMap `shouldBe` Just (-5) + M.lookup 2 varMap `shouldBe` Just 6 + -- Verify objective value computed from variables + (2 * x1 - x2) `shouldBe` (-16) + _ -> expectationFailure "Unexpected result format" + + describe "Systems with GEQ constraints and negative bounds" $ do + it "Max x₁ with x₁ ≥ -5, x₁ ≥ 2 (GEQ tightens bound)" $ do + -- Lower bound is -5 but GEQ constraint says x₁ ≥ 2 + -- Without upper bound, this is unbounded for Max + -- Add an upper bound via another constraint + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 2 + , LEQ (M.fromList [(1, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ ≥ -5, x₁ ≥ 2 (GEQ tightens bound)" $ do + -- Minimize with GEQ 2, so minimum is at x₁ = 2 + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 2 + , LEQ (M.fromList [(1, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 2 + _ -> expectationFailure "Unexpected result format" + + describe "Systems with EQ constraints and negative bounds" $ do + it "Max x₁ + x₂ with x₁ - x₂ = 0, x₁ ≥ -5, x₂ ≥ -5, x₁ ≤ 10" $ do + -- x₁ = x₂, maximize x₁ + x₂ = 2x₁ + -- With x₁ ≤ 10, optimal is x₁ = x₂ = 10, obj = 20 + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (2, -1)]) 0 + , LEQ (M.fromList [(1, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5)), (2, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let objVal = computeObjective obj varMap + M.lookup 1 varMap `shouldBe` Just 10 + M.lookup 2 varMap `shouldBe` Just 10 + -- Verify objective value + objVal `shouldBe` 20 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ + x₂ with x₁ - x₂ = 0, x₁ ≥ -5, x₂ ≥ -5, x₁ ≤ 10" $ do + -- x₁ = x₂, minimize x₁ + x₂ = 2x₁ + -- Lower bound is -5, so optimal is x₁ = x₂ = -5, obj = -10 + let obj = Min (M.fromList [(1, 1), (2, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (2, -1)]) 0 + , LEQ (M.fromList [(1, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5)), (2, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let objVal = computeObjective obj varMap + M.lookup 1 varMap `shouldBe` Just (-5) + M.lookup 2 varMap `shouldBe` Just (-5) + -- Verify objective value + objVal `shouldBe` (-10) + _ -> expectationFailure "Unexpected result format" + + describe "Fractional negative bounds" $ do + it "Max x₁ with x₁ ≥ -7/2, x₁ ≤ 5/2" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) (5 % 2)] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly ((-7) % 2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (5 % 2) + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ ≥ -7/2, x₁ ≤ 5/2" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) (5 % 2)] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly ((-7) % 2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just ((-7) % 2) + _ -> expectationFailure "Unexpected result format" + + describe "twoPhaseSimplex with unbounded variables (Split transformation)" $ do + describe "Simple single variable systems" $ do + it "Max x₁ with -10 ≤ x₁ ≤ 10 (unbounded var with box constraints)" $ do + -- x₁ is unbounded but constrained by -10 ≤ x₁ ≤ 10 + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 10 + , GEQ (M.fromList [(1, 1)]) (-10) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 10 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with -10 ≤ x₁ ≤ 10 (unbounded var with box constraints)" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 10 + , GEQ (M.fromList [(1, 1)]) (-10) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just (-10) + _ -> expectationFailure "Unexpected result format" + + it "unbounded variable with only upper bound: Min finds negative value" $ do + -- x₁ unbounded, only x₁ ≤ 5, minimize x₁ + -- This should be unbounded (no finite solution) since x₁ can go to -∞ + let obj = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 5] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + -- This should be unbounded (no finite optimum exists) + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + _ -> expectationFailure "Expected unbounded system" + + describe "Two variable systems with unbounded variables" $ do + it "Max x₁ + x₂ with unbounded vars, -5 ≤ x₁ ≤ 5, -3 ≤ x₂ ≤ 7" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 5 + , GEQ (M.fromList [(1, 1)]) (-5) + , LEQ (M.fromList [(2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) (-3) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 5 + M.lookup 2 varMap `shouldBe` Just 7 + let objVal = computeObjective obj varMap + objVal `shouldBe` 12 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ + x₂ with unbounded vars, -5 ≤ x₁ ≤ 5, -3 ≤ x₂ ≤ 7" $ do + let obj = Min (M.fromList [(1, 1), (2, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 5 + , GEQ (M.fromList [(1, 1)]) (-5) + , LEQ (M.fromList [(2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) (-3) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just (-5) + M.lookup 2 varMap `shouldBe` Just (-3) + let objVal = computeObjective obj varMap + objVal `shouldBe` (-8) + _ -> expectationFailure "Unexpected result format" + + it "Max x₁ - x₂ with unbounded vars: x₁ up, x₂ down" $ do + -- Maximize x₁ - x₂: want x₁ large (5) and x₂ small (-3) + let obj = Max (M.fromList [(1, 1), (2, -1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 5 + , GEQ (M.fromList [(1, 1)]) (-5) + , LEQ (M.fromList [(2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) (-3) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 5 + M.lookup 2 varMap `shouldBe` Just (-3) + let objVal = computeObjective obj varMap + objVal `shouldBe` 8 + _ -> expectationFailure "Unexpected result format" + + describe "Systems with EQ constraints and unbounded variables" $ do + it "Max x₁ with x₁ + x₂ = 10, unbounded vars, x₂ ≥ -5" $ do + -- x₁ + x₂ = 10, x₂ ≥ -5, unbounded x₁ + -- Maximize x₁: make x₂ as small as possible (-5), so x₁ = 15 + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (2, 1)]) 10 + , GEQ (M.fromList [(2, 1)]) (-5) + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just 15 + M.lookup 2 varMap `shouldBe` Just (-5) + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ with x₁ + x₂ = 10, unbounded vars, x₂ ≤ 20" $ do + -- x₁ + x₂ = 10, x₂ ≤ 20, unbounded x₁ + -- Minimize x₁: make x₂ as large as possible (20), so x₁ = -10 + let obj = Min (M.fromList [(1, 1)]) + constraints = + [ EQ (M.fromList [(1, 1), (2, 1)]) 10 + , LEQ (M.fromList [(2, 1)]) 20 + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded), (2, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just (-10) + M.lookup 2 varMap `shouldBe` Just 20 + _ -> expectationFailure "Unexpected result format" + + describe "twoPhaseSimplex with mixed domain types" $ do + describe "NonNegative, negative lower bound, and unbounded in same system" $ do + it "Max x₁ + x₂ + x₃ with x₁ ≥ 0, x₂ ≥ -5, x₃ unbounded, sum ≤ 20" $ do + -- x₁ non-negative, x₂ has lower bound -5, x₃ unbounded + -- All constrained by sum ≤ 20 and individual bounds + let obj = Max (M.fromList [(1, 1), (2, 1), (3, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1), (2, 1), (3, 1)]) 20 + , LEQ (M.fromList [(1, 1)]) 10 + , LEQ (M.fromList [(2, 1)]) 8 + , LEQ (M.fromList [(3, 1)]) 15 + , GEQ (M.fromList [(3, 1)]) (-10) + ] + domainMap = + VarDomainMap $ + M.fromList + [(1, nonNegative), (2, lowerBoundOnly (-5)), (3, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let objVal = computeObjective obj varMap + -- Verify objective value + objVal `shouldBe` 20 + _ -> expectationFailure "Unexpected result format" + + it "Min x₁ + x₂ + x₃ with x₁ ≥ 0, x₂ ≥ -5, x₃ unbounded, sum ≥ -10" $ do + -- Minimize sum with lower bound constraint + let obj = Min (M.fromList [(1, 1), (2, 1), (3, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1), (3, 1)]) (-10) + , LEQ (M.fromList [(1, 1)]) 10 + , LEQ (M.fromList [(2, 1)]) 8 + , LEQ (M.fromList [(3, 1)]) 15 + , GEQ (M.fromList [(3, 1)]) (-20) + ] + domainMap = + VarDomainMap $ + M.fromList + [(1, nonNegative), (2, lowerBoundOnly (-5)), (3, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + x3 = M.findWithDefault 0 3 varMap + objVal = computeObjective obj varMap + -- Verify constraints + x1 `shouldSatisfy` (>= 0) + x2 `shouldSatisfy` (>= (-5)) + x3 `shouldSatisfy` (>= (-20)) + -- Verify objective value + objVal `shouldBe` (-10) + _ -> expectationFailure "Unexpected result format" + + describe "Positive lower bound with other domain types" $ do + it "Max 2x₁ + 3x₂ with x₁ ≥ 2 (positive bound), x₂ ≥ -3, 2x₁ + x₂ ≤ 20" $ do + -- x₁ has positive lower bound (uses AddLowerBound) + -- x₂ has negative lower bound (uses Shift) + let obj = Max (M.fromList [(1, 2), (2, 3)]) + constraints = + [ LEQ (M.fromList [(1, 2), (2, 1)]) 20 + , LEQ (M.fromList [(2, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 2), (2, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + -- Verify constraints + x1 `shouldSatisfy` (>= 2) + x2 `shouldSatisfy` (>= (-3)) + (2 * x1 + x2) `shouldSatisfy` (<= 20) + _ -> expectationFailure "Unexpected result format" + + it "Min 2x₁ + 3x₂ with x₁ ≥ 2, x₂ ≥ -3, x₁ + x₂ ≥ 0" $ do + -- Minimize with lower bounds + -- x₁ = 2 (minimum), x₂ = -2 (to satisfy x₁ + x₂ ≥ 0) + let obj = Min (M.fromList [(1, 2), (2, 3)]) + constraints = + [ GEQ (M.fromList [(1, 1), (2, 1)]) 0 + , LEQ (M.fromList [(1, 1)]) 10 + , LEQ (M.fromList [(2, 1)]) 10 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 2), (2, lowerBoundOnly (-3))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + x1 `shouldSatisfy` (>= 2) + x2 `shouldSatisfy` (>= (-3)) + (x1 + x2) `shouldSatisfy` (>= 0) + _ -> expectationFailure "Unexpected result format" + + describe "twoPhaseSimplex edge cases and infeasibility" $ do + it "Infeasible: negative lower bound conflicts with GEQ constraint" $ do + -- x₁ ≥ -5 (domain), but x₁ ≥ 10 and x₁ ≤ 5 (constraints conflict) + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 10 + , LEQ (M.fromList [(1, 1)]) 5 + ] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible system" + + it "Infeasible: unbounded variable with conflicting constraints" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ GEQ (M.fromList [(1, 1)]) 10 + , LEQ (M.fromList [(1, 1)]) 5 + ] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> pure () + _ -> expectationFailure "Expected infeasible system" + + it "Variable at exactly zero with negative lower bound" $ do + -- x₁ ≥ -5, constraint x₁ = 0 + let obj = Max (M.fromList [(1, 1)]) + constraints = [EQ (M.fromList [(1, 1)]) 0] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 0 + _ -> expectationFailure "Unexpected result format" + + it "unbounded variable constrained to zero" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [EQ (M.fromList [(1, 1)]) 0] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> M.lookup 1 varMap `shouldBe` Just 0 + _ -> expectationFailure "Unexpected result format" + + it "Multiple variables, only some with negative bounds" $ do + -- x₁ ≥ 0 (non-negative), x₂ ≥ -10, x₃ ≥ 0 + -- Max x₁ + x₂ + x₃ with x₁ + x₂ + x₃ ≤ 15 + let obj = Max (M.fromList [(1, 1), (2, 1), (3, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1), (3, 1)]) 15] + domainMap = + VarDomainMap $ + M.fromList + [(1, nonNegative), (2, lowerBoundOnly (-10)), (3, nonNegative)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult Nothing _ -> expectationFailure "Expected a solution but got Nothing" + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let objVal = computeObjective obj varMap + -- Verify objective value + objVal `shouldBe` 15 + _ -> expectationFailure "Unexpected result format" + + -- =========================================================================== + -- Tests for internal preprocessing functions + -- =========================================================================== + + describe "collectAllVars" $ do + describe "Unit tests" $ do + it "collects variables from Max objective" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + constraints = [] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 2] + + it "collects variables from Min objective" $ do + let obj = Min (M.fromList [(3, 1), (4, -2)]) + constraints = [] + collectAllVars [obj] constraints `shouldBe` Set.fromList [3, 4] + + it "collects variables from LEQ constraint" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(2, 1), (3, 2)]) 10] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 2, 3] + + it "collects variables from GEQ constraint" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [GEQ (M.fromList [(4, 1)]) 5] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 4] + + it "collects variables from EQ constraint" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [EQ (M.fromList [(5, 2), (6, 3)]) 15] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 5, 6] + + it "collects variables from mixed constraints" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(2, 1)]) 10 + , GEQ (M.fromList [(3, 1)]) 5 + , EQ (M.fromList [(4, 1)]) 7 + ] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 2, 3, 4] + + it "handles empty objective coefficients" $ do + let obj = Max M.empty + constraints = [LEQ (M.fromList [(1, 1)]) 10] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1] + + it "handles empty constraints" $ do + let obj = Max (M.fromList [(1, 1), (2, 2)]) + constraints = [] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 2] + + it "deduplicates variables appearing in multiple places" $ do + let obj = Max (M.fromList [(1, 1), (2, 2)]) + constraints = + [ LEQ (M.fromList [(1, 3), (3, 4)]) 10 + , GEQ (M.fromList [(2, 5), (3, 6)]) 5 + ] + collectAllVars [obj] constraints `shouldBe` Set.fromList [1, 2, 3] + + describe "getTransform" $ do + describe "Unit tests" $ do + it "returns empty list for NonNegative domain" $ do + getTransform 10 1 nonNegative `shouldBe` ([], 0) + + it "returns empty list for lowerBoundOnly 0" $ do + getTransform 10 1 (lowerBoundOnly 0) `shouldBe` ([], 0) + + it "returns AddLowerBound for positive lower bound" $ do + getTransform 10 1 (lowerBoundOnly 5) `shouldBe` ([AddLowerBound 1 5], 0) + + it "returns AddLowerBound for fractional positive lower bound" $ do + getTransform 10 1 (lowerBoundOnly (3 % 2)) `shouldBe` ([AddLowerBound 1 (3 % 2)], 0) + + it "returns Shift for negative lower bound" $ do + getTransform 10 1 (lowerBoundOnly (-5)) `shouldBe` ([Shift 1 10 (-5)], 1) + + it "returns Shift for fractional negative lower bound" $ do + getTransform 10 1 (lowerBoundOnly ((-7) % 3)) `shouldBe` ([Shift 1 10 ((-7) % 3)], 1) + + it "returns Split for unbounded domain" $ do + getTransform 10 1 unbounded `shouldBe` ([Split 1 10 11], 2) + + it "returns AddUpperBound for upper bound only" $ do + getTransform 10 1 (upperBoundOnly 5) `shouldBe` ([Split 1 10 11, AddUpperBound 1 5], 2) + + it "returns AddLowerBound and AddUpperBound for bounded range" $ do + getTransform 10 1 (boundedRange 2 10) `shouldBe` ([AddLowerBound 1 2, AddUpperBound 1 10], 0) + + it "returns Shift and AddUpperBound for negative lower bound with upper bound" $ do + getTransform 10 1 (boundedRange (-5) 10) `shouldBe` ([Shift 1 10 (-5), AddUpperBound 1 10], 1) + + describe "generateTransform" $ do + describe "Unit tests" $ do + it "generates no transform for NonNegative in domain map" $ do + let domainMap = M.fromList [(1, nonNegative)] + generateTransform domainMap 1 ([], 10) `shouldBe` ([], 10) + + it "generates AddLowerBound for positive bound in domain map" $ do + let domainMap = M.fromList [(1, lowerBoundOnly 5)] + generateTransform domainMap 1 ([], 10) `shouldBe` ([AddLowerBound 1 5], 10) + + it "generates Shift for negative bound and increments fresh var" $ do + let domainMap = M.fromList [(1, lowerBoundOnly (-5))] + generateTransform domainMap 1 ([], 10) `shouldBe` ([Shift 1 10 (-5)], 11) + + it "generates Split for unbounded and increments fresh var by 2" $ do + let domainMap = M.fromList [(1, unbounded)] + generateTransform domainMap 1 ([], 10) `shouldBe` ([Split 1 10 11], 12) + + it "treats variable not in domain map as unbounded" $ do + let domainMap = M.empty + generateTransform domainMap 1 ([], 10) `shouldBe` ([Split 1 10 11], 12) + + it "accumulates transforms" $ do + let domainMap = M.fromList [(1, lowerBoundOnly 5)] + existing = [AddLowerBound 2 3] + generateTransform domainMap 1 (existing, 10) `shouldBe` ([AddLowerBound 1 5] ++ existing, 10) + + it "generates AddUpperBound for upper bound" $ do + let domainMap = M.fromList [(1, boundedRange 0 10)] + generateTransform domainMap 1 ([], 10) `shouldBe` ([AddUpperBound 1 10], 10) + + describe "applyShiftToObjective" $ do + describe "Unit tests" $ do + it "substitutes variable in Max objective" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + applyShiftToObjective 1 10 (-5) obj `shouldBe` Max (M.fromList [(10, 3), (2, 5)]) + + it "substitutes variable in Min objective" $ do + let obj = Min (M.fromList [(1, -2), (2, 4)]) + applyShiftToObjective 1 10 (-3) obj `shouldBe` Min (M.fromList [(10, -2), (2, 4)]) + + it "leaves objective unchanged if variable not present" $ do + let obj = Max (M.fromList [(2, 5), (3, 7)]) + applyShiftToObjective 1 10 (-5) obj `shouldBe` Max (M.fromList [(2, 5), (3, 7)]) + + it "preserves coefficient during substitution" $ do + let obj = Max (M.fromList [(1, 100)]) + applyShiftToObjective 1 10 (-5) obj `shouldBe` Max (M.fromList [(10, 100)]) + + describe "applyShiftToConstraint" $ do + describe "Unit tests" $ do + it "shifts LEQ constraint correctly" $ do + -- x1 = x10 + (-5), so x1 has shift -5 + -- constraint: 2*x1 <= 10 becomes 2*x10 <= 10 - 2*(-5) = 20 + let constraint = LEQ (M.fromList [(1, 2)]) 10 + applyShiftToConstraint 1 10 (-5) constraint `shouldBe` LEQ (M.fromList [(10, 2)]) 20 + + it "shifts GEQ constraint correctly" $ do + let constraint = GEQ (M.fromList [(1, 3)]) 6 + applyShiftToConstraint 1 10 (-2) constraint `shouldBe` GEQ (M.fromList [(10, 3)]) 12 + + it "shifts EQ constraint correctly" $ do + let constraint = EQ (M.fromList [(1, 4)]) 8 + applyShiftToConstraint 1 10 (-1) constraint `shouldBe` EQ (M.fromList [(10, 4)]) 12 + + it "leaves constraint unchanged if variable not present" $ do + let constraint = LEQ (M.fromList [(2, 5)]) 10 + applyShiftToConstraint 1 10 (-5) constraint `shouldBe` LEQ (M.fromList [(2, 5)]) 10 + + it "handles negative coefficients" $ do + -- x1 = x10 + (-5), constraint: -3*x1 <= 10 + -- becomes -3*x10 <= 10 - (-3)*(-5) = 10 - 15 = -5 + let constraint = LEQ (M.fromList [(1, -3)]) 10 + applyShiftToConstraint 1 10 (-5) constraint `shouldBe` LEQ (M.fromList [(10, -3)]) (-5) + + it "handles multiple variables in constraint" $ do + let constraint = LEQ (M.fromList [(1, 2), (2, 3)]) 10 + applyShiftToConstraint 1 10 (-5) constraint `shouldBe` LEQ (M.fromList [(10, 2), (2, 3)]) 20 + + describe "applySplitToObjective" $ do + describe "Unit tests" $ do + it "splits variable in Max objective" $ do + let obj = Max (M.fromList [(1, 3)]) + -- x1 = x10 - x11, so coeff 3 -> x10 gets 3, x11 gets -3 + applySplitToObjective 1 10 11 obj `shouldBe` Max (M.fromList [(10, 3), (11, -3)]) + + it "splits variable in Min objective" $ do + let obj = Min (M.fromList [(1, -2)]) + applySplitToObjective 1 10 11 obj `shouldBe` Min (M.fromList [(10, -2), (11, 2)]) + + it "leaves objective unchanged if variable not present" $ do + let obj = Max (M.fromList [(2, 5)]) + applySplitToObjective 1 10 11 obj `shouldBe` Max (M.fromList [(2, 5)]) + + it "handles multiple variables" $ do + let obj = Max (M.fromList [(1, 3), (2, 5)]) + applySplitToObjective 1 10 11 obj `shouldBe` Max (M.fromList [(10, 3), (11, -3), (2, 5)]) + + describe "applySplitToConstraint" $ do + describe "Unit tests" $ do + it "splits variable in LEQ constraint" $ do + let constraint = LEQ (M.fromList [(1, 2)]) 10 + applySplitToConstraint 1 10 11 constraint `shouldBe` LEQ (M.fromList [(10, 2), (11, -2)]) 10 + + it "splits variable in GEQ constraint" $ do + let constraint = GEQ (M.fromList [(1, 3)]) 5 + applySplitToConstraint 1 10 11 constraint `shouldBe` GEQ (M.fromList [(10, 3), (11, -3)]) 5 + + it "splits variable in EQ constraint" $ do + let constraint = EQ (M.fromList [(1, 4)]) 8 + applySplitToConstraint 1 10 11 constraint `shouldBe` EQ (M.fromList [(10, 4), (11, -4)]) 8 + + it "leaves constraint unchanged if variable not present" $ do + let constraint = LEQ (M.fromList [(2, 5)]) 10 + applySplitToConstraint 1 10 11 constraint `shouldBe` LEQ (M.fromList [(2, 5)]) 10 + + it "handles negative coefficients" $ do + let constraint = LEQ (M.fromList [(1, -3)]) 10 + applySplitToConstraint 1 10 11 constraint `shouldBe` LEQ (M.fromList [(10, -3), (11, 3)]) 10 + + it "handles multiple variables" $ do + let constraint = LEQ (M.fromList [(1, 2), (2, 3)]) 10 + applySplitToConstraint 1 10 11 constraint `shouldBe` LEQ (M.fromList [(10, 2), (11, -2), (2, 3)]) 10 + + describe "applyTransform and applyTransforms" $ do + describe "Unit tests" $ do + it "applyTransform AddLowerBound adds GEQ constraint" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + transform = AddLowerBound 1 5 + applyTransform transform (obj, constraints) + `shouldBe` (obj, [GEQ (M.singleton 1 1) 5, LEQ (M.fromList [(1, 1)]) 10]) + + it "applyTransform Shift transforms objective and constraints" $ do + let obj = Max (M.fromList [(1, 2)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + transform = Shift 1 10 (-5) + let (newObj, newConstraints) = applyTransform transform (obj, constraints) + newObj `shouldBe` Max (M.fromList [(10, 2)]) + newConstraints `shouldBe` [LEQ (M.fromList [(10, 1)]) 15] + + it "applyTransform Split transforms objective and constraints" $ do + let obj = Max (M.fromList [(1, 3)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + transform = Split 1 10 11 + let (newObj, newConstraints) = applyTransform transform (obj, constraints) + newObj `shouldBe` Max (M.fromList [(10, 3), (11, -3)]) + newConstraints `shouldBe` [LEQ (M.fromList [(10, 1), (11, -1)]) 10] + + it "applyTransforms applies multiple transforms in order" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 10] + transforms = [AddLowerBound 1 5, AddLowerBound 2 3] + let (newObj, newConstraints) = applyTransforms transforms obj constraints + newObj `shouldBe` obj + -- Two GEQ constraints should be added + length newConstraints `shouldBe` 3 + + describe "unapplyTransformToVarMap and unapplyTransformsToVarMap" $ do + describe "Unit tests" $ do + it "unapplyTransformToVarMap AddLowerBound leaves result unchanged" $ do + let varVals = M.fromList [(1, 7)] + transform = AddLowerBound 1 5 + unapplyTransformToVarMap transform varVals `shouldBe` varVals + + it "unapplyTransformToVarMap Shift recovers original variable value" $ do + -- originalVar = shiftedVar + shiftBy + -- If shiftedVar = 15 and shiftBy = -5, then originalVar = 10 + let varVals = M.fromList [(10, 15)] + transform = Shift 1 10 (-5) + let newVarVals = unapplyTransformToVarMap transform varVals + M.lookup 1 newVarVals `shouldBe` Just 10 + M.lookup 10 newVarVals `shouldBe` Nothing + + it "unapplyTransformToVarMap Split recovers original variable value" $ do + -- originalVar = posVar - negVar + -- If posVar = 8 and negVar = 3, then originalVar = 5 + let varVals = M.fromList [(10, 8), (11, 3)] + transform = Split 1 10 11 + let newVarVals = unapplyTransformToVarMap transform varVals + M.lookup 1 newVarVals `shouldBe` Just 5 + M.lookup 10 newVarVals `shouldBe` Nothing + M.lookup 11 newVarVals `shouldBe` Nothing + + it "unapplyTransformToVarMap Split handles negative original value" $ do + -- originalVar = posVar - negVar + -- If posVar = 2 and negVar = 7, then originalVar = -5 + let varVals = M.fromList [(10, 2), (11, 7)] + transform = Split 1 10 11 + let newVarVals = unapplyTransformToVarMap transform varVals + M.lookup 1 newVarVals `shouldBe` Just (-5) + + it "unapplyTransformsToVarMap applies in correct order (reverse of apply)" $ do + -- Two shifts: var 1 shifted by -5 to var 10, var 2 shifted by -3 to var 11 + let varVals = M.fromList [(10, 15), (11, 8)] + transforms = [Shift 1 10 (-5), Shift 2 11 (-3)] + let newVarVals = unapplyTransformsToVarMap transforms varVals + M.lookup 1 newVarVals `shouldBe` Just 10 + M.lookup 2 newVarVals `shouldBe` Just 5 + + describe "preprocess" $ do + describe "Unit tests" $ do + it "returns empty transforms for all NonNegative domains" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative), (2, nonNegative)] + let ([newObj], newConstraints, transforms) = preprocess [obj] domainMap constraints + transforms `shouldBe` [] + newObj `shouldBe` obj + newConstraints `shouldBe` constraints + + it "generates AddLowerBound for positive lower bounds" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly 5)] + let (_, newConstraints, transforms) = preprocess [obj] domainMap constraints + transforms `shouldBe` [AddLowerBound 1 5] + length newConstraints `shouldBe` 2 -- original + GEQ + it "generates Shift for negative lower bounds" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))] + let ([newObj], newConstraints, transforms) = preprocess [obj] domainMap constraints + length transforms `shouldBe` 1 + case head transforms of + Shift {..} -> do + originalVar `shouldBe` 1 + shiftBy `shouldBe` (-5) + _ -> expectationFailure "Expected Shift transform" + + it "generates Split for unbounded domains" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, unbounded)] + let (_, _, transforms) = preprocess [obj] domainMap constraints + length transforms `shouldBe` 1 + case head transforms of + Split {..} -> originalVar `shouldBe` 1 + _ -> expectationFailure "Expected Split transform" + + it "handles mixed domain types" $ do + let obj = Max (M.fromList [(1, 1), (2, 1), (3, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1), (3, 1)]) 10] + domainMap = + VarDomainMap $ + M.fromList + [(1, nonNegative), (2, lowerBoundOnly 5), (3, lowerBoundOnly (-3))] + let (_, _, transforms) = preprocess [obj] domainMap constraints + -- Should have AddLowerBound for var 2, Shift for var 3 + length transforms `shouldBe` 2 + + -- =========================================================================== + -- Property-based tests + -- =========================================================================== + + describe "Property-based tests" $ do + describe "collectAllVars properties" $ do + it "result is non-empty when objective is non-empty" $ + property $ + \(NonEmpty coeffs :: NonEmptyList (Int, Rational)) -> + let obj = Max (M.fromList [(abs k `mod` 100 + 1, v) | (k, v) <- coeffs]) + in not (Set.null (collectAllVars [obj] [])) + + it "result contains all objective variables" $ + property $ + \(vars :: [Int]) -> + let posVars = filter (> 0) (map abs vars) + obj = Max (M.fromList [(v, 1) | v <- take 5 posVars]) + in all (`Set.member` collectAllVars [obj] []) (M.keys $ case obj of Max m -> m; Min m -> m) + + describe "getTransform properties" $ do + it "NonNegative always produces empty list" $ + property $ + \(nextVar :: Int) (v :: Int) -> + getTransform (abs nextVar + 1) (abs v + 1) nonNegative == ([], 0) + + it "lowerBoundOnly 0 produces empty list" $ + property $ + \(nextVar :: Int) (v :: Int) -> + getTransform (abs nextVar + 1) (abs v + 1) (lowerBoundOnly 0) == ([], 0) + + it "positive lowerBoundOnly produces AddLowerBound" $ + property $ + \(Positive bound :: Positive Rational) (nextVar :: Int) (v :: Int) -> + case getTransform (abs nextVar + 1) (abs v + 1) (lowerBoundOnly bound) of + ([AddLowerBound var b], 0) -> var == abs v + 1 && b == bound + _ -> False + + it "negative lowerBoundOnly produces Shift" $ + property $ + \(Positive bound :: Positive Rational) (nextVar :: Int) (v :: Int) -> + let negBound = negate bound + in case getTransform (abs nextVar + 1) (abs v + 1) (lowerBoundOnly negBound) of + ([Shift origVar _ shiftBy], 1) -> origVar == abs v + 1 && shiftBy == negBound + _ -> False + + it "unbounded produces Split" $ + property $ + \(nextVar :: Int) (v :: Int) -> + case getTransform (abs nextVar + 1) (abs v + 1) unbounded of + ([Split origVar _ _], 2) -> origVar == abs v + 1 + _ -> False + + it "boundedRange produces both lower and upper bound transforms" $ + property $ + \(lower :: Rational) (Positive delta :: Positive Rational) (nextVar :: Int) (v :: Int) -> + let upper = lower + delta -- ensure upper > lower + in case getTransform (abs nextVar + 1) (abs v + 1) (boundedRange lower upper) of + (transforms, _) -> + any (\case AddUpperBound var u -> var == abs v + 1 && u == upper; _ -> False) transforms + + describe "applyShiftToConstraint properties" $ do + it "RHS adjustment follows formula: newRHS = oldRHS - coeff * shiftBy" $ + property $ + \(coeff :: Rational) (oldRHS :: Rational) (shiftBy :: Rational) -> + coeff + /= 0 + ==> let constraint = LEQ (M.fromList [(1, coeff)]) oldRHS + LEQ _ newRHS = applyShiftToConstraint 1 10 shiftBy constraint + in newRHS == oldRHS - coeff * shiftBy + + it "preserves constraint type (LEQ stays LEQ)" $ + property $ + \(coeff :: Rational) (rhs :: Rational) (shiftBy :: Rational) -> + coeff + /= 0 + ==> let constraint = LEQ (M.fromList [(1, coeff)]) rhs + in case applyShiftToConstraint 1 10 shiftBy constraint of + LEQ {} -> True + _ -> False + + it "preserves constraint type (GEQ stays GEQ)" $ + property $ + \(coeff :: Rational) (rhs :: Rational) (shiftBy :: Rational) -> + coeff + /= 0 + ==> let constraint = GEQ (M.fromList [(1, coeff)]) rhs + in case applyShiftToConstraint 1 10 shiftBy constraint of + GEQ {} -> True + _ -> False + + describe "applySplitToConstraint properties" $ do + it "preserves RHS value" $ + property $ + \(coeff :: Rational) (rhs :: Rational) -> + coeff + /= 0 + ==> let constraint = LEQ (M.fromList [(1, coeff)]) rhs + LEQ _ newRHS = applySplitToConstraint 1 10 11 constraint + in newRHS == rhs + + it "negVar coefficient is negation of posVar coefficient" $ + property $ + \(coeff :: Rational) (rhs :: Rational) -> + coeff + /= 0 + ==> let constraint = LEQ (M.fromList [(1, coeff)]) rhs + LEQ m _ = applySplitToConstraint 1 10 11 constraint + posCoeff = M.findWithDefault 0 10 m + negCoeff = M.findWithDefault 0 11 m + in negCoeff == negate posCoeff + + describe "unapplyTransformToVarMap Shift properties" $ do + it "recovers originalVar = shiftedVar + shiftBy" $ + property $ + \(shiftedVal :: Rational) (shiftBy :: Rational) -> + let varMap = M.fromList [(5, 100), (10, shiftedVal)] + transform = Shift 1 10 shiftBy + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 1 newVarMap == Just (shiftedVal + shiftBy) + + it "removes shifted variable from result" $ + property $ + \(shiftedVal :: Rational) (shiftBy :: Rational) -> + let varMap = M.fromList [(5, 100), (10, shiftedVal)] + transform = Shift 1 10 shiftBy + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 10 newVarMap == Nothing + + describe "unapplyTransformToVarMap Split properties" $ do + it "recovers originalVar = posVar - negVar" $ + property $ + \(posVal :: Rational) (negVal :: Rational) -> + let varMap = M.fromList [(5, 100), (10, posVal), (11, negVal)] + transform = Split 1 10 11 + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 1 newVarMap == Just (posVal - negVal) + + it "removes pos and neg variables from result" $ + property $ + \(posVal :: Rational) (negVal :: Rational) -> + let varMap = M.fromList [(5, 100), (10, posVal), (11, negVal)] + transform = Split 1 10 11 + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 10 newVarMap == Nothing + && M.lookup 11 newVarMap == Nothing + + describe "Round-trip properties" $ do + it "Shift transform and unapply is identity for variable value" $ + property $ + \(origVal :: Rational) (shiftBy :: Rational) -> + shiftBy + < 0 + ==> let shiftedVal -- Only negative shifts are valid + = + origVal - shiftBy -- shiftedVar = originalVar - shiftBy + varMap = M.fromList [(5, 100), (10, shiftedVal)] + transform = Shift 1 10 shiftBy + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 1 newVarMap == Just origVal + + it "Split with posVal=origVal and negVal=0 gives correct value for positive origVal" $ + property $ + \(Positive origVal :: Positive Rational) -> + let varMap = M.fromList [(5, 100), (10, origVal), (11, 0)] + transform = Split 1 10 11 + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 1 newVarMap == Just origVal + + it "Split with posVal=0 and negVal=-origVal gives correct value for negative origVal" $ + property $ + \(Positive origVal :: Positive Rational) -> + let negOrigVal = negate origVal + varMap = M.fromList [(5, 100), (10, 0), (11, origVal)] + transform = Split 1 10 11 + newVarMap = unapplyTransformToVarMap transform varMap + in M.lookup 1 newVarMap == Just negOrigVal + + describe "twoPhaseSimplex with upperBoundOnly domain" $ do + it "Max x\x2081 with x\x2081 \x2264 10 (upper bound only): x\x2081 unconstrained below, unbounded" $ do + -- upperBoundOnly means no lower bound (split) + upper bound + -- Max x\x2081 with only x\x2081 \x2264 10: unbounded below, but maximizing so optimal at 10 + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, upperBoundOnly 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + M.findWithDefault 0 1 varMap `shouldBe` 10 + _ -> expectationFailure "Expected optimal at upper bound" + + it "Min x\x2081 with x\x2081 \x2264 10 (upper bound only): unbounded (no lower bound)" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, upperBoundOnly 10)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ Unbounded] -> pure () + _ -> expectationFailure "Expected unbounded (no lower bound)" + + it "Max x\x2081 with x\x2081 \x2264 5 and x\x2081 + x\x2082 \x2264 8, x\x2082 upperBoundOnly 3" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 8] + domainMap = VarDomainMap $ M.fromList [(1, upperBoundOnly 5), (2, upperBoundOnly 3)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + x1 `shouldSatisfy` (<= 5) + x2 `shouldSatisfy` (<= 3) + computeObjective obj varMap `shouldBe` 8 + _ -> expectationFailure "Expected optimal result" + + it "Min x\x2081 with x\x2081 \x2264 5, x\x2081 \x2265 -3: optimal at lower bound" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [GEQ (M.fromList [(1, 1)]) (-3)] + domainMap = VarDomainMap $ M.fromList [(1, upperBoundOnly 5)] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + M.lookup 1 varMap `shouldBe` Just (-3) + _ -> expectationFailure "Expected optimal at -3" + + describe "twoPhaseSimplex with fully-negative boundedRange" $ do + it "Max x\x2081 with -10 \x2264 x\x2081 \x2264 -2: optimal at x\x2081=-2" $ do + let obj = Max (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-10) (-2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + M.lookup 1 varMap `shouldBe` Just (-2) + _ -> expectationFailure "Expected optimal at -2" + + it "Min x\x2081 with -10 \x2264 x\x2081 \x2264 -2: optimal at x\x2081=-10" $ do + let obj = Min (M.fromList [(1, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-10) (-2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> + M.lookup 1 varMap `shouldBe` Just (-10) + _ -> expectationFailure "Expected optimal at -10" + + it "Max x\x2081 + x\x2082 with -8 \x2264 x\x2081 \x2264 -1, -6 \x2264 x\x2082 \x2264 -2" $ do + let obj = Max (M.fromList [(1, 1), (2, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-8) (-1)), (2, boundedRange (-6) (-2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just (-1) + M.lookup 2 varMap `shouldBe` Just (-2) + computeObjective obj varMap `shouldBe` (-3) + _ -> expectationFailure "Expected optimal at (-1, -2)" + + it "Min x\x2081 + x\x2082 with -8 \x2264 x\x2081 \x2264 -1, -6 \x2264 x\x2082 \x2264 -2" $ do + let obj = Min (M.fromList [(1, 1), (2, 1)]) + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-8) (-1)), (2, boundedRange (-6) (-2))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + M.lookup 1 varMap `shouldBe` Just (-8) + M.lookup 2 varMap `shouldBe` Just (-6) + computeObjective obj varMap `shouldBe` (-14) + _ -> expectationFailure "Expected optimal at (-8, -6)" + + it "Max x\x2081 with -5 \x2264 x\x2081 \x2264 -1 and x\x2081 + x\x2082 \x2264 0, x\x2082 in [-5, -1]" $ do + -- With constraint x\x2081 + x\x2082 \x2264 0, and both negative ranges + let obj = Max (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1), (2, 1)]) 0] + domainMap = VarDomainMap $ M.fromList [(1, boundedRange (-5) (-1)), (2, boundedRange (-5) (-1))] + actualResult <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj] constraints + case actualResult of + SimplexResult (Just _) [ObjectiveResult _ (Optimal varMap)] -> do + let x1 = M.findWithDefault 0 1 varMap + x2 = M.findWithDefault 0 2 varMap + x1 `shouldBe` (-1) + (x1 + x2) `shouldSatisfy` (<= 0) + _ -> expectationFailure "Expected optimal result" + + describe "postprocess" $ do + it "passes through Unbounded unchanged" $ do + let originalVars = Set.fromList [1, 2] + transforms = [Shift 1 3 (-5)] + postprocess originalVars transforms Unbounded `shouldBe` Unbounded + + it "filters to original variables only" $ do + let originalVars = Set.fromList [1, 2] + transforms = [] + -- varMap has original vars and some extra (slack/artificial) + varMap = M.fromList [(1, 5), (2, 10), (3, 99), (4, 100)] + postprocess originalVars transforms (Optimal varMap) `shouldBe` Optimal (M.fromList [(1, 5), (2, 10)]) + + it "unapplies Shift transform and filters" $ do + -- Shift: original var 1 was shifted to var 3 with offset -5 + -- In transformed space: var 3 = 8, meaning original var 1 = 8 + (-5) = 3 + let originalVars = Set.fromList [1] + transforms = [Shift 1 3 (-5)] + varMap = M.fromList [(3, 8)] + postprocess originalVars transforms (Optimal varMap) `shouldBe` Optimal (M.fromList [(1, 3)]) + + it "unapplies Split transform and filters" $ do + -- Split: original var 1 was split into pos var 3 and neg var 4 + -- original var 1 = var 3 - var 4 + let originalVars = Set.fromList [1] + transforms = [Split 1 3 4] + varMap = M.fromList [(3, 7), (4, 2)] + postprocess originalVars transforms (Optimal varMap) `shouldBe` Optimal (M.fromList [(1, 5)]) + + it "unapplies multiple transforms and filters" $ do + -- var 1: Shift by -3 to var 3 + -- var 2: Split into vars 4 and 5 + let originalVars = Set.fromList [1, 2] + transforms = [Shift 1 3 (-3), Split 2 4 5] + varMap = M.fromList [(3, 10), (4, 6), (5, 1)] + postprocess originalVars transforms (Optimal varMap) + `shouldBe` Optimal (M.fromList [(1, 7), (2, 5)]) + + it "returns empty map when no original vars have values" $ do + let originalVars = Set.fromList [1] + transforms = [] + varMap = M.fromList [(2, 5)] + postprocess originalVars transforms (Optimal varMap) `shouldBe` Optimal M.empty + + describe "prettyShowVarLitMapSum" $ do + it "shows empty map as empty string" $ do + prettyShowVarLitMapSum M.empty `shouldBe` "" + + it "shows single positive coefficient" $ do + prettyShowVarLitMapSum (M.fromList [(1, 3)]) `shouldBe` "3 * 1 + " + + it "shows single negative coefficient with parentheses" $ do + prettyShowVarLitMapSum (M.fromList [(1, -2)]) `shouldBe` "(-2) * 1 + " + + it "shows multiple coefficients" $ do + let result = prettyShowVarLitMapSum (M.fromList [(1, 2), (2, 3)]) + result `shouldBe` "2 * 1 + 3 * 2 + " + + describe "prettyShowPolyConstraint" $ do + it "shows LEQ constraint" $ do + prettyShowPolyConstraint (LEQ (M.fromList [(1, 2)]) 10) `shouldBe` "2 * 1 + " ++ " <= " ++ show (10 :: Rational) + + it "shows GEQ constraint" $ do + prettyShowPolyConstraint (GEQ (M.fromList [(1, 1)]) 5) `shouldBe` "1 * 1 + " ++ " >= " ++ show (5 :: Rational) + + it "shows EQ constraint" $ do + prettyShowPolyConstraint (EQ (M.fromList [(1, 1)]) 3) `shouldBe` "1 * 1 + " ++ " == " ++ show (3 :: Rational) + + describe "prettyShowObjectiveFunction" $ do + it "shows Max objective" $ do + prettyShowObjectiveFunction (Max (M.fromList [(1, 2)])) `shouldBe` "max: 2 * 1 + " + + it "shows Min objective" $ do + prettyShowObjectiveFunction (Min (M.fromList [(1, 5)])) `shouldBe` "min: 5 * 1 + " + + describe "twoPhaseSimplex with multiple objectives" $ do + it "optimizes two objectives over the same feasible region" $ do + -- Feasible region: x₁ + x₂ ≤ 10, x₁ ≤ 6, x₂ ≤ 8, x₁,x₂ ≥ 0 + -- Max x₁: optimal at x₁=6, x₂=0 (or x₁=6, x₂=4) with obj=6 + -- Max x₂: optimal at x₁=0, x₂=8 (or x₁=2, x₂=8) with obj=8 + let obj1 = Max (M.fromList [(1, 1)]) -- Max x₁ + obj2 = Max (M.fromList [(2, 1)]) -- Max x₂ + constraints = + [ LEQ (M.fromList [(1, 1), (2, 1)]) 10 + , LEQ (M.fromList [(1, 1)]) 6 + , LEQ (M.fromList [(2, 1)]) 8 + ] + allVars = collectAllVars [obj1, obj2] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + SimplexResult mFeasibleSystem objResults <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj1, obj2] constraints + -- Should have a feasible system + mFeasibleSystem `shouldSatisfy` isJust + -- Should have two results + length objResults `shouldBe` 2 + -- First result (Max x₁) should have x₁=6 + case objResults !! 0 of + ObjectiveResult _ (Optimal varVals) -> + M.lookup 1 varVals `shouldBe` Just 6 + _ -> expectationFailure "Expected optimal result for obj1" + -- Second result (Max x₂) should have x₂=8 + case objResults !! 1 of + ObjectiveResult _ (Optimal varVals) -> + M.lookup 2 varVals `shouldBe` Just 8 + _ -> expectationFailure "Expected optimal result for obj2" + + it "handles empty objective list returning feasible system only" $ do + let constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromSet (const nonNegative) (Set.singleton 1) + SimplexResult mFeasibleSystem objResults <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [] constraints + mFeasibleSystem `shouldSatisfy` isJust + length objResults `shouldBe` 0 + + it "handles infeasible system with multiple objectives" $ do + -- x₁ ≤ 5 and x₁ ≥ 10 is infeasible + let obj1 = Max (M.fromList [(1, 1)]) + obj2 = Min (M.fromList [(1, 1)]) + constraints = + [ LEQ (M.fromList [(1, 1)]) 5 + , GEQ (M.fromList [(1, 1)]) 10 + ] + allVars = collectAllVars [obj1, obj2] constraints + domainMap = VarDomainMap $ M.fromSet (const nonNegative) allVars + SimplexResult mFeasibleSystem objResults <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj1, obj2] constraints + -- Should be infeasible + mFeasibleSystem `shouldBe` Nothing + -- No objective results when infeasible + length objResults `shouldBe` 0 + + it "optimizes Max and Min of same function over feasible region" $ do + -- Feasible region: 0 ≤ x₁ ≤ 10 + -- Max x₁: optimal at x₁=10 + -- Min x₁: optimal at x₁=0 + let obj1 = Max (M.fromList [(1, 1)]) + obj2 = Min (M.fromList [(1, 1)]) + constraints = [LEQ (M.fromList [(1, 1)]) 10] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative)] + SimplexResult mFeasibleSystem objResults <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj1, obj2] constraints + mFeasibleSystem `shouldSatisfy` isJust + length objResults `shouldBe` 2 + -- Max x₁ should be 10 + case objResults !! 0 of + ObjectiveResult _ (Optimal varVals) -> + M.lookup 1 varVals `shouldBe` Just 10 + _ -> expectationFailure "Expected optimal result for Max x₁" + -- Min x₁ should be 0 (or not present in map if zero) + case objResults !! 1 of + ObjectiveResult _ (Optimal varVals) -> + M.findWithDefault 0 1 varVals `shouldBe` 0 + _ -> expectationFailure "Expected optimal result for Min x₁" + + it "handles one unbounded objective among multiple objectives" $ do + -- x₁ with only a lower bound (non-negative) + -- Max x₁: unbounded (no upper constraint) + -- Min x₁ with x₁ ≥ 0: optimal at x₁=0 + let obj1 = Max (M.fromList [(1, 1)]) -- This will be unbounded + obj2 = Min (M.fromList [(1, 1)]) -- This will have optimal at 0 + constraints = [] + domainMap = VarDomainMap $ M.fromList [(1, nonNegative)] + SimplexResult mFeasibleSystem objResults <- + runStdoutLoggingT $ + filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ + twoPhaseSimplex domainMap [obj1, obj2] constraints + mFeasibleSystem `shouldSatisfy` isJust + length objResults `shouldBe` 2 + -- Max x₁ should be unbounded + case objResults !! 0 of + ObjectiveResult _ Unbounded -> pure () + _ -> expectationFailure "Expected unbounded result for Max x₁" + -- Min x₁ should be 0 + case objResults !! 1 of + ObjectiveResult _ (Optimal varVals) -> + M.findWithDefault 0 1 varVals `shouldBe` 0 + _ -> expectationFailure "Expected optimal result for Min x₁" diff --git a/test/Linear/Simplex/UtilSpec.hs b/test/Linear/Simplex/UtilSpec.hs new file mode 100644 index 0000000..6b74215 --- /dev/null +++ b/test/Linear/Simplex/UtilSpec.hs @@ -0,0 +1,291 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Linear.Simplex.UtilSpec where + +import Prelude hiding (EQ) + +import Control.Exception (evaluate) +import qualified Data.Map as M +import Test.Hspec (Spec, anyErrorCall, describe, expectationFailure, it, shouldBe, shouldThrow) +import Test.QuickCheck (Positive (..), property) + +import Linear.Simplex.Types + ( DictValue (..) + , ObjectiveFunction (..) + , PivotObjective (..) + , PolyConstraint (..) + , TableauRow (..) + , VarLitMapSum + ) +import Linear.Simplex.Util + ( combineVarLitMapSums + , dictionaryFormToTableau + , foldVarLitMap + , insertPivotObjectiveToDict + , isMax + , simplifySystem + , tableauInDictionaryForm + ) + +spec :: Spec +spec = do + describe "isMax" $ do + it "returns True for Max" $ do + isMax (Max (M.fromList [(1, 1)])) `shouldBe` True + + it "returns False for Min" $ do + isMax (Min (M.fromList [(1, 1)])) `shouldBe` False + + it "returns True for Max with empty coefficients" $ do + isMax (Max M.empty) `shouldBe` True + + it "returns False for Min with empty coefficients" $ do + isMax (Min M.empty) `shouldBe` False + + describe "simplifySystem" $ do + describe "Unit tests" $ do + it "returns empty list for empty input" $ do + simplifySystem [] `shouldBe` [] + + it "preserves a single LEQ constraint" $ do + let c = LEQ (M.fromList [(1, 2)]) 10 + simplifySystem [c] `shouldBe` [c] + + it "preserves a single GEQ constraint" $ do + let c = GEQ (M.fromList [(1, 1)]) 5 + simplifySystem [c] `shouldBe` [c] + + it "preserves a single EQ constraint" $ do + let c = EQ (M.fromList [(1, 3)]) 15 + simplifySystem [c] `shouldBe` [c] + + it "reduces matching LEQ and GEQ into EQ" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + simplifySystem [LEQ lhs rhs, GEQ lhs rhs] `shouldBe` [EQ lhs rhs] + + it "reduces matching GEQ and LEQ into EQ" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + simplifySystem [GEQ lhs rhs, LEQ lhs rhs] `shouldBe` [EQ lhs rhs] + + it "keeps non-matching LEQ and GEQ separate" $ do + let c1 = LEQ (M.fromList [(1, 1)]) 10 + c2 = GEQ (M.fromList [(1, 1)]) 5 + simplifySystem [c1, c2] `shouldBe` [c1, c2] + + it "removes duplicate constraints" $ do + let c = LEQ (M.fromList [(1, 2)]) 10 + simplifySystem [c, c] `shouldBe` [c] + + it "reduces LEQ matching an existing EQ" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + simplifySystem [LEQ lhs rhs, EQ lhs rhs] `shouldBe` [EQ lhs rhs] + + it "reduces GEQ matching an existing EQ" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + simplifySystem [GEQ lhs rhs, EQ lhs rhs] `shouldBe` [EQ lhs rhs] + + it "keeps EQ and removes matching LEQ from rest" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + simplifySystem [EQ lhs rhs, LEQ lhs rhs] `shouldBe` [EQ lhs rhs] + + it "preserves unrelated constraints alongside a reduction" $ do + let lhs = M.fromList [(1, 1)] + rhs = 5 + unrelated = GEQ (M.fromList [(2, 3)]) 7 + simplifySystem [LEQ lhs rhs, GEQ lhs rhs, unrelated] + `shouldBe` [EQ lhs rhs, unrelated] + + describe "Properties" $ do + it "idempotent: simplifying twice equals simplifying once" $ + property $ + \(constraints :: [(Int, Rational, Rational)]) -> + let mkConstraint (tag, coeff, rhs) = + case tag `mod` 3 of + 0 -> LEQ (M.fromList [(1, coeff)]) rhs + 1 -> GEQ (M.fromList [(1, coeff)]) rhs + _ -> EQ (M.fromList [(1, coeff)]) rhs + system = map mkConstraint constraints + in simplifySystem (simplifySystem system) == simplifySystem system + + it "never increases the number of constraints" $ + property $ + \(constraints :: [(Int, Rational, Rational)]) -> + let mkConstraint (tag, coeff, rhs) = + case tag `mod` 3 of + 0 -> LEQ (M.fromList [(1, coeff)]) rhs + 1 -> GEQ (M.fromList [(1, coeff)]) rhs + _ -> EQ (M.fromList [(1, coeff)]) rhs + system = map mkConstraint constraints + in length (simplifySystem system) <= length system + + it "matching LEQ and GEQ always produce EQ" $ + property $ + \(coeff :: Rational, rhs :: Rational) -> + let lhs = M.fromList [(1, coeff)] + in simplifySystem [LEQ lhs rhs, GEQ lhs rhs] == [EQ lhs rhs] + + describe "dictionaryFormToTableau and tableauInDictionaryForm" $ do + describe "Unit tests" $ do + it "converts a simple dictionary entry to tableau form" $ do + -- Dict: x₁ = 3 + 2*x₂ means DictValue {varMapSum = {2: 2}, constant = 3} + -- Tableau: x₁ - 2*x₂ = 3 means TableauRow {lhs = {1: 1, 2: -2}, rhs = 3} + let dict = M.fromList [(1, DictValue {varMapSum = M.fromList [(2, 2)], constant = 3})] + tableau = dictionaryFormToTableau dict + case M.lookup 1 tableau of + Just row -> do + row.lhs `shouldBe` M.fromList [(1, 1), (2, -2)] + row.rhs `shouldBe` 3 + Nothing -> expectationFailure "Expected row for var 1" + + it "converts an empty dictionary to empty tableau" $ do + dictionaryFormToTableau M.empty `shouldBe` M.empty + + it "converts a simple tableau entry to dictionary form" $ do + -- Tableau: 1*x₁ + (-2)*x₂ = 3 means x₁ = 3 + 2*x₂ + let tableau = M.fromList [(1, TableauRow {lhs = M.fromList [(1, 1), (2, -2)], rhs = 3})] + dict = tableauInDictionaryForm tableau + case M.lookup 1 dict of + Just entry -> do + entry.varMapSum `shouldBe` M.fromList [(2, 2)] + entry.constant `shouldBe` 3 + Nothing -> expectationFailure "Expected entry for var 1" + + it "converts an empty tableau to empty dictionary" $ do + tableauInDictionaryForm M.empty `shouldBe` M.empty + + describe "Round-trip properties" $ do + it "tableauInDictionaryForm . dictionaryFormToTableau is identity" $ + property $ + \(pairs :: [(Int, [(Int, Rational)], Rational)]) -> + let mkEntry (var, coeffs, c) = + let basicVar = abs var + 1 + -- Ensure no coefficient var clashes with the basic var + safeCoeffs = filter (\(v, _) -> abs v + 1000 /= basicVar) coeffs + in (basicVar, DictValue {varMapSum = M.fromList (map (\(v, coeff) -> (abs v + 1000, coeff)) safeCoeffs), constant = c}) + dict = M.fromList (map mkEntry pairs) + roundTripped = tableauInDictionaryForm (dictionaryFormToTableau dict) + in roundTripped == dict + + it "dictionaryFormToTableau . tableauInDictionaryForm is identity for well-formed tableaux" $ + property $ + \(var :: Positive Int, coeffPairs :: [(Positive Int, Rational)], c :: Rational) -> + let basicVar = getPositive var + -- Ensure basic var has coefficient 1 (well-formed tableau) + otherCoeffs = M.fromList [(getPositive v + basicVar, coeff) | (v, coeff) <- coeffPairs] + lhs = M.insert basicVar 1 otherCoeffs + tableau = M.singleton basicVar (TableauRow {lhs = lhs, rhs = c}) + roundTripped = dictionaryFormToTableau (tableauInDictionaryForm tableau) + in roundTripped == tableau + + describe "combineVarLitMapSums" $ do + describe "Unit tests" $ do + it "combines two disjoint maps" $ do + let m1 = M.fromList [(1, 3)] + m2 = M.fromList [(2, 5)] + combineVarLitMapSums m1 m2 `shouldBe` M.fromList [(1, 3), (2, 5)] + + it "sums values for overlapping keys" $ do + let m1 = M.fromList [(1, 3), (2, 4)] + m2 = M.fromList [(1, 7), (3, 1)] + combineVarLitMapSums m1 m2 `shouldBe` M.fromList [(1, 10), (2, 4), (3, 1)] + + it "returns other map when one is empty" $ do + let m = M.fromList [(1, 5)] + combineVarLitMapSums M.empty m `shouldBe` m + combineVarLitMapSums m M.empty `shouldBe` m + + it "returns empty for two empty maps" $ do + combineVarLitMapSums M.empty M.empty `shouldBe` (M.empty :: VarLitMapSum) + + describe "Properties" $ do + it "is commutative" $ + property $ + \(pairs1 :: [(Int, Rational)], pairs2 :: [(Int, Rational)]) -> + let m1 = M.fromList pairs1 + m2 = M.fromList pairs2 + in combineVarLitMapSums m1 m2 == combineVarLitMapSums m2 m1 + + it "has empty map as identity" $ + property $ + \(pairs :: [(Int, Rational)]) -> + let m = M.fromList pairs + in combineVarLitMapSums m M.empty == m + && combineVarLitMapSums M.empty m == m + + it "is associative" $ + property $ + \(p1 :: [(Int, Rational)], p2 :: [(Int, Rational)], p3 :: [(Int, Rational)]) -> + let m1 = M.fromList p1 + m2 = M.fromList p2 + m3 = M.fromList p3 + in combineVarLitMapSums (combineVarLitMapSums m1 m2) m3 + == combineVarLitMapSums m1 (combineVarLitMapSums m2 m3) + + describe "foldVarLitMap" $ do + describe "Unit tests" $ do + it "returns the single map for a singleton list" $ do + let m = M.fromList [(1, 5), (2, 3)] + foldVarLitMap [m] `shouldBe` m + + it "combines two maps" $ do + let m1 = M.fromList [(1, 3)] + m2 = M.fromList [(1, 7), (2, 4)] + foldVarLitMap [m1, m2] `shouldBe` M.fromList [(1, 10), (2, 4)] + + it "combines three maps" $ do + let m1 = M.fromList [(1, 1)] + m2 = M.fromList [(1, 2), (2, 3)] + m3 = M.fromList [(2, 7), (3, 5)] + foldVarLitMap [m1, m2, m3] `shouldBe` M.fromList [(1, 3), (2, 10), (3, 5)] + + it "throws error on empty list" $ do + evaluate (foldVarLitMap []) `shouldThrow` anyErrorCall + + describe "Properties" $ do + it "folding a singleton list is identity" $ + property $ + \(pairs :: [(Int, Rational)]) -> + let m = M.fromList pairs + in foldVarLitMap [m] == m + + it "folding two maps equals combineVarLitMapSums" $ + property $ + \(p1 :: [(Int, Rational)], p2 :: [(Int, Rational)]) -> + let m1 = M.fromList p1 + m2 = M.fromList p2 + in foldVarLitMap [m1, m2] == combineVarLitMapSums m1 m2 + + describe "insertPivotObjectiveToDict" $ do + it "inserts objective into empty dictionary" $ do + let obj = PivotObjective {variable = 1, function = M.fromList [(2, 3)], constant = 5} + result = insertPivotObjectiveToDict obj M.empty + case M.lookup 1 result of + Just entry -> do + entry.varMapSum `shouldBe` M.fromList [(2, 3)] + entry.constant `shouldBe` 5 + Nothing -> expectationFailure "Expected entry for objective variable" + + it "inserts objective into non-empty dictionary without overwriting others" $ do + let existing = M.fromList [(2, DictValue {varMapSum = M.fromList [(3, 1)], constant = 10})] + obj = PivotObjective {variable = 1, function = M.fromList [(3, 4)], constant = 7} + result = insertPivotObjectiveToDict obj existing + M.size result `shouldBe` 2 + case M.lookup 2 result of + Just entry -> entry.constant `shouldBe` 10 + Nothing -> expectationFailure "Existing entry should be preserved" + + it "overwrites existing entry with same variable" $ do + let existing = M.fromList [(1, DictValue {varMapSum = M.fromList [(2, 1)], constant = 3})] + obj = PivotObjective {variable = 1, function = M.fromList [(2, 99)], constant = 42} + result = insertPivotObjectiveToDict obj existing + M.size result `shouldBe` 1 + case M.lookup 1 result of + Just entry -> do + entry.varMapSum `shouldBe` M.fromList [(2, 99)] + entry.constant `shouldBe` 42 + Nothing -> expectationFailure "Expected updated entry" diff --git a/test/Spec.hs b/test/Spec.hs index 4a8ad55..a824f8c 100644 --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,42 +1 @@ -module Main where - -import Control.Monad -import Control.Monad.IO.Class -import Control.Monad.Logger - -import Linear.Simplex.Prettify -import Linear.Simplex.Solver.TwoPhase -import Linear.Simplex.Types -import Linear.Simplex.Util - -import TestFunctions - -main :: IO () -main = runStdoutLoggingT $ filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ runTests testsList - -runTests :: (MonadLogger m, MonadFail m, MonadIO m) => [((ObjectiveFunction, [PolyConstraint]), Maybe Result)] -> m () -runTests [] = do - liftIO $ putStrLn "All tests passed" - pure () -runTests (((testObjective, testConstraints), expectedResult) : tests) = - do - testResult <- twoPhaseSimplex testObjective testConstraints - if testResult == expectedResult - then runTests tests - else do - let msg = - "\nThe following test failed: " - <> ("\nObjective Function (Non-prettified): " ++ show testObjective) - <> ("\nConstraints (Non-prettified): " ++ show testConstraints) - <> "\n====================================" - <> ("\nObjective Function (Prettified): " ++ prettyShowObjectiveFunction testObjective) - <> "\nConstraints (Prettified): " - <> "\n" - <> concatMap (\c -> "\t" ++ prettyShowPolyConstraint c ++ "\n") testConstraints - <> "\n====================================" - <> ("\nExpected Solution (Full): " ++ show expectedResult) - <> ("\nActual Solution (Full): " ++ show testResult) - <> ("\nExpected Solution (Objective): " ++ show (extractObjectiveValue expectedResult)) - <> ("\nActual Solution (Objective): " ++ show (extractObjectiveValue testResult)) - <> "\n" - fail msg +{-# OPTIONS_GHC -F -pgmF hspec-discover #-} diff --git a/test/TestFunctions.hs b/test/TestFunctions.hs deleted file mode 100644 index b2af317..0000000 --- a/test/TestFunctions.hs +++ /dev/null @@ -1,1048 +0,0 @@ -module TestFunctions where - -import qualified Data.Map as M -import Data.Ratio -import Linear.Simplex.Types -import Prelude hiding (EQ) - -testsList :: [((ObjectiveFunction, [PolyConstraint]), Maybe Result)] -testsList = - [ (test1, Just (Result 7 (M.fromList [(7, 29), (1, 3), (2, 4)]))) - , (test2, Just (Result 7 (M.fromList [(7, 0)]))) - , (test3, Nothing) - , (test4, Just (Result 11 (M.fromList [(11, 237 % 7), (1, 24 % 7), (2, 33 % 7)]))) - , (test5, Just (Result 9 (M.fromList [(9, 3 % 5), (2, 14 % 5), (3, 17 % 5)]))) - , (test6, Nothing) - , (test7, Just (Result 8 (M.fromList [(8, 1), (2, 2), (1, 3)]))) - , (test8, Just (Result 8 (M.fromList [(8, (-1) % 4), (2, 9 % 2), (1, 17 % 4)]))) - , (test9, Just (Result 7 (M.fromList [(7, 5), (3, 2), (4, 1)]))) - , (test10, Just (Result 7 (M.fromList [(7, 8), (1, 2), (2, 6)]))) - , (test11, Just (Result 8 (M.fromList [(8, 20), (4, 16), (3, 6)]))) - , (test12, Just (Result 8 (M.fromList [(8, 6), (4, 2), (5, 2)]))) - , (test13, Just (Result 6 (M.fromList [(6, 150), (2, 150)]))) - , (test14, Just (Result 6 (M.fromList [(6, 40 % 3), (2, 40 % 3)]))) - , (test15, Nothing) - , (test16, Just (Result 6 (M.fromList [(6, 75), (1, 75 % 2)]))) - , (test17, Just (Result 7 (M.fromList [(7, (-120)), (1, 20)]))) - , (test18, Just (Result 7 (M.fromList [(7, 10), (3, 5)]))) - , (test19, Nothing) - , (test20, Nothing) - , (test21, Just (Result 7 (M.fromList [(7, 250), (2, 50)]))) - , (test22, Just (Result 7 (M.fromList [(7, 0)]))) - , (test23, Nothing) - , (test24, Just (Result 10 (M.fromList [(10, 300), (3, 150)]))) - , (test25, Just (Result 3 (M.fromList [(3, 15), (1, 15)]))) - , (test26, Just (Result 6 (M.fromList [(6, 20), (1, 10), (2, 10)]))) - , (test27, Just (Result 3 (M.fromList [(3, 0)]))) - , (test28, Just (Result 6 (M.fromList [(6, 0), (2, 10)]))) - , (test29, Nothing) - , (test30, Nothing) - , (test31, Just (Result 5 (M.fromList [(2, 1 % 1), (5, 0 % 1)]))) - , (test32, Nothing) - , (testPolyPaver1, Just (Result 12 (M.fromList [(12, 7 % 4), (2, 5 % 2), (1, 7 % 4), (3, 0)]))) - , (testPolyPaver2, Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 3), (1, 5 % 2), (3, 0)]))) - , (testPolyPaver3, Just (Result 12 (M.fromList [(12, 5 % 3), (2, 5 % 3), (1, 5 % 2), (3, 0)]))) - , (testPolyPaver4, Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 2), (1, 5 % 2), (3, 0)]))) - , (testPolyPaver5, Nothing) - , (testPolyPaver6, Nothing) - , (testPolyPaver7, Nothing) - , (testPolyPaver8, Nothing) - , (testPolyPaver9, Just (Result 12 (M.fromList [(12, 7 % 2), (2, 5 % 9), (1, 7 % 2), (3, 0)]))) - , (testPolyPaver10, Just (Result 12 (M.fromList [(12, 17 % 20), (2, 7 % 2), (1, 17 % 20), (3, 0)]))) - , (testPolyPaver11, Just (Result 12 (M.fromList [(12, 7 % 2), (2, 7 % 2), (1, 22 % 9)]))) - , (testPolyPaver12, Just (Result 12 (M.fromList [(12, 5 % 9), (2, 5 % 9), (1, 7 % 2), (3, 0)]))) - , (testPolyPaverTwoFs1, Nothing) - , (testPolyPaverTwoFs2, Nothing) - , (testPolyPaverTwoFs3, Nothing) - , (testPolyPaverTwoFs4, Nothing) - , (testPolyPaverTwoFs5, Just (Result 17 (M.fromList [(17, 5 % 2), (2, 45 % 22), (1, 5 % 2), (4, 0)]))) - , (testPolyPaverTwoFs6, Just (Result 17 (M.fromList [(17, 45 % 22), (2, 5 % 2), (1, 45 % 22), (4, 0)]))) - , (testPolyPaverTwoFs7, Just (Result 17 (M.fromList [(17, 5 % 2), (2, 5 % 2), (1, 5 % 2), (4, 0)]))) - , (testPolyPaverTwoFs8, Just (Result 17 (M.fromList [(17, 45 % 22), (2, 45 % 22), (1, 5 % 2), (4, 0)]))) - , (testLeqGeqBugMin1, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) - , (testLeqGeqBugMax1, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) - , (testLeqGeqBugMin2, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) - , (testLeqGeqBugMax2, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) - , (testQuickCheck1, Just (Result 10 (M.fromList [(10, (-370)), (2, 26), (1, 5 % 3)]))) - , (testQuickCheck2, Just (Result 8 (M.fromList [(8, (-2) % 9), (1, 14 % 9), (2, 8 % 9)]))) - , (testQuickCheck3, Just (Result 7 (M.fromList [(7, (-8)), (2, 2)]))) - ] - -testLeqGeqBugMin1 :: (ObjectiveFunction, [PolyConstraint]) -testLeqGeqBugMin1 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1)]) 3 - , LEQ (M.fromList [(1, 1)]) 3 - , GEQ (M.fromList [(2, 1)]) 3 - , LEQ (M.fromList [(2, 1)]) 3 - ] - ) - -testLeqGeqBugMax1 :: (ObjectiveFunction, [PolyConstraint]) -testLeqGeqBugMax1 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1)]) 3 - , LEQ (M.fromList [(1, 1)]) 3 - , GEQ (M.fromList [(2, 1)]) 3 - , LEQ (M.fromList [(2, 1)]) 3 - ] - ) - -testLeqGeqBugMin2 :: (ObjectiveFunction, [PolyConstraint]) -testLeqGeqBugMin2 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1)]) 3 - , LEQ (M.fromList [(1, 1)]) 3 - , GEQ (M.fromList [(2, 1)]) 3 - , LEQ (M.fromList [(2, 1)]) 3 - ] - ) - -testLeqGeqBugMax2 :: (ObjectiveFunction, [PolyConstraint]) -testLeqGeqBugMax2 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1)]) 3 - , LEQ (M.fromList [(1, 1)]) 3 - , GEQ (M.fromList [(2, 1)]) 3 - , LEQ (M.fromList [(2, 1)]) 3 - ] - ) - --- From page 50 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 29, 1 = 3, 2 = 4, -test1 :: (ObjectiveFunction, [PolyConstraint]) -test1 = - ( Max (M.fromList [(1, 3), (2, 5)]) - , - [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 - , LEQ (M.fromList [(1, 1), (2, 1)]) 7 - , LEQ (M.fromList [(2, 1)]) 4 - , LEQ (M.fromList [(1, -1), (2, 2)]) 6 - ] - ) - -test2 :: (ObjectiveFunction, [PolyConstraint]) -test2 = - ( Min (M.fromList [(1, 3), (2, 5)]) - , - [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 - , LEQ (M.fromList [(1, 1), (2, 1)]) 7 - , LEQ (M.fromList [(2, 1)]) 4 - , LEQ (M.fromList [(1, -1), (2, 2)]) 6 - ] - ) - -test3 :: (ObjectiveFunction, [PolyConstraint]) -test3 = - ( Max (M.fromList [(1, 3), (2, 5)]) - , - [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 - , GEQ (M.fromList [(1, 1), (2, 1)]) 7 - , GEQ (M.fromList [(2, 1)]) 4 - , GEQ (M.fromList [(1, -1), (2, 2)]) 6 - ] - ) - -test4 :: (ObjectiveFunction, [PolyConstraint]) -test4 = - ( Min (M.fromList [(1, 3), (2, 5)]) - , - [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 - , GEQ (M.fromList [(1, 1), (2, 1)]) 7 - , GEQ (M.fromList [(2, 1)]) 4 - , GEQ (M.fromList [(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 (M.fromList [(1, 1), (2, -1), (3, 1)]) - , - [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 - , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) - , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) - ] - ) - -test6 :: (ObjectiveFunction, [PolyConstraint]) -test6 = - ( Min (M.fromList [(1, 1), (2, -1), (3, 1)]) - , - [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 - , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) - , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) - ] - ) - -test7 :: (ObjectiveFunction, [PolyConstraint]) -test7 = - ( Max (M.fromList [(1, 1), (2, -1), (3, 1)]) - , - [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 - , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) - , GEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) - ] - ) - -test8 :: (ObjectiveFunction, [PolyConstraint]) -test8 = - ( Min (M.fromList [(1, 1), (2, -1), (3, 1)]) - , - [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 - , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) - , GEQ (M.fromList [(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 (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) - , - [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 - , EQ (M.fromList [(2, 1), (3, 1), (4, 4)]) 6 - ] - ) - -test10 :: (ObjectiveFunction, [PolyConstraint]) -test10 = - ( Max (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) - , - [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 - , EQ (M.fromList [(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 (M.fromList [(3, -2), (4, 2), (5, 1)]) - , - [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 - , EQ (M.fromList [(3, 3), (4, -1), (5, 2)]) 2 - ] - ) - -test12 :: (ObjectiveFunction, [PolyConstraint]) -test12 = - ( Min (M.fromList [(3, -2), (4, 2), (5, 1)]) - , - [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 - , EQ (M.fromList [(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 (M.fromList [(1, 2), (2, 1)]) - , - [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 - , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) - ] - ) - -test14 :: (ObjectiveFunction, [PolyConstraint]) -test14 = - ( Min (M.fromList [(1, 2), (2, 1)]) - , - [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 - , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) - ] - ) - -test15 :: (ObjectiveFunction, [PolyConstraint]) -test15 = - ( Max (M.fromList [(1, 2), (2, 1)]) - , - [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 - , GEQ (M.fromList [(1, 2), (2, -3)]) (-40) - ] - ) - -test16 :: (ObjectiveFunction, [PolyConstraint]) -test16 = - ( Min (M.fromList [(1, 2), (2, 1)]) - , - [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 - , GEQ (M.fromList [(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 (M.fromList [(1, -6), (2, -4), (3, 2)]) - , - [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 - , LEQ (M.fromList [(2, -5), (3, 5)]) 100 - , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 - ] - ) - -test18 :: (ObjectiveFunction, [PolyConstraint]) -test18 = - ( Max (M.fromList [(1, -6), (2, -4), (3, 2)]) - , - [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 - , LEQ (M.fromList [(2, -5), (3, 5)]) 100 - , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 - ] - ) - -test19 :: (ObjectiveFunction, [PolyConstraint]) -test19 = - ( Min (M.fromList [(1, -6), (2, -4), (3, 2)]) - , - [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 - , GEQ (M.fromList [(2, -5), (3, 5)]) 100 - , GEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 - ] - ) - -test20 :: (ObjectiveFunction, [PolyConstraint]) -test20 = - ( Max (M.fromList [(1, -6), (2, -4), (3, 2)]) - , - [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 - , GEQ (M.fromList [(2, -5), (3, 5)]) 100 - , GEQ (M.fromList [(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 (M.fromList [(1, 3), (2, 5), (3, 2)]) - , - [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 - , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 - , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 - ] - ) - -test22 :: (ObjectiveFunction, [PolyConstraint]) -test22 = - ( Min (M.fromList [(1, 3), (2, 5), (3, 2)]) - , - [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 - , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 - , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 - ] - ) - -test23 :: (ObjectiveFunction, [PolyConstraint]) -test23 = - ( Max (M.fromList [(1, 3), (2, 5), (3, 2)]) - , - [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 - , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 - , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 - ] - ) - -test24 :: (ObjectiveFunction, [PolyConstraint]) -test24 = - ( Min (M.fromList [(1, 3), (2, 5), (3, 2)]) - , - [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 - , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 - , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 - ] - ) - -test25 :: (ObjectiveFunction, [PolyConstraint]) -test25 = - ( Max (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, 1)]) 15 - ] - ) - -test26 :: (ObjectiveFunction, [PolyConstraint]) -test26 = - ( Max (M.fromList [(1, 2)]) - , - [ LEQ (M.fromList [(1, 2)]) 20 - , GEQ (M.fromList [(2, 1)]) 10 - ] - ) - -test27 :: (ObjectiveFunction, [PolyConstraint]) -test27 = - ( Min (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, 1)]) 15 - ] - ) - -test28 :: (ObjectiveFunction, [PolyConstraint]) -test28 = - ( Min (M.fromList [(1, 2)]) - , - [ LEQ (M.fromList [(1, 2)]) 20 - , GEQ (M.fromList [(2, 1)]) 10 - ] - ) - -test29 :: (ObjectiveFunction, [PolyConstraint]) -test29 = - ( Max (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, 1)]) 15 - , GEQ (M.fromList [(1, 1)]) 15.01 - ] - ) - -test30 :: (ObjectiveFunction, [PolyConstraint]) -test30 = - ( Max (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, 1)]) 15 - , GEQ (M.fromList [(1, 1)]) 15.01 - , GEQ (M.fromList [(2, 1)]) 10 - ] - ) - -test31 :: (ObjectiveFunction, [PolyConstraint]) -test31 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1), (2, 1)]) 1 - , GEQ (M.fromList [(1, 1), (2, 1)]) 1 - ] - ) - -test32 :: (ObjectiveFunction, [PolyConstraint]) -test32 = - ( Min (M.fromList [(1, 1)]) - , - [ GEQ (M.fromList [(1, 1), (2, 1)]) 2 - , LEQ (M.fromList [(1, 1), (2, 1)]) 1 - ] - ) - --- Tests for systems similar to those from PolyPaver2 -testPolyPaver1 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver1 = - ( Min (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) - , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(1, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 (M.fromList [(2, 1)]) - , - [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) - , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) - , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) - , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) - , GEQ (M.fromList [(1, 1)]) x1l - , LEQ (M.fromList [(1, 1)]) x1r - , GEQ (M.fromList [(2, 1)]) x2l - , LEQ (M.fromList [(2, 1)]) x2r - , LEQ (M.fromList [(3, 1)]) 0 - , LEQ (M.fromList [(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 - -testQuickCheck1 :: (ObjectiveFunction, [PolyConstraint]) -testQuickCheck1 = - ( Max (M.fromList [(1, 12), (2, -15)]) - , - [ EQ (M.fromList [(1, 24), (2, -2)]) (-12) - , GEQ (M.fromList [(1, -20), (2, 11)]) (-7) - , GEQ (M.fromList [(1, -28), (2, 5)]) (-8) - , GEQ (M.fromList [(1, 3), (2, 0)]) 5 - , LEQ (M.fromList [(1, -48)]) (-1) - ] - ) - --- Correct solution is -2/9 -testQuickCheck2 :: (ObjectiveFunction, [PolyConstraint]) -testQuickCheck2 = - ( Max (M.fromList [(1, -3), (2, 5)]) - , - [ LEQ (M.fromList [(1, -6), (2, 6)]) 4 - , LEQ (M.fromList [(1, 1), (2, -4), (3, 3)]) (-2) - , LEQ (M.fromList [(2, 7), (1, -4)]) 0 - ] - ) - --- This test will fail if the objective function is not simplified -testQuickCheck3 :: (ObjectiveFunction, [PolyConstraint]) -testQuickCheck3 = - ( Min (M.fromList [(2, 0), (2, -4)]) - , - [ GEQ (M.fromList [(1, 5), (2, 4)]) (-4) - , LEQ (M.fromList [(1, -1), (2, -1)]) 2 - , LEQ (M.fromList [(2, 1)]) 2 - , GEQ (M.fromList [(1, -5), (2, -1), (2, 1)]) (-5) - ] - )