diff --git a/.github/workflows/haskell.yml b/.github/workflows/haskell.yml
index 92d3748..bbf33cf 100644
--- a/.github/workflows/haskell.yml
+++ b/.github/workflows/haskell.yml
@@ -2,22 +2,18 @@ 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@3b7702b41516aa428dfe6e295dc73476ae58f69e # v11
         with:
           version: "0.14.0.0"
   build:
@@ -27,104 +23,58 @@ jobs:
       fail-fast: false
       matrix:
         os: [windows-latest, macos-latest, ubuntu-latest]
-        ghc-version: ['9.6', '9.4', '9.2', '9.0']
+        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@96f3dafd067155f32643c2a0757ab71d2910e2c2 # v2.8.0
         id: setup
         with:
           ghc-version: ${{ matrix.ghc-version }}
-          enable-stack: true
 
       - name: Installed minor versions of GHC, Cabal, and Stack
         shell: bash
         run: |
           GHC_VERSION=$(ghc --numeric-version)
           CABAL_VERSION=$(cabal --numeric-version)
-          STACK_VERSION=$(stack --numeric-version)
           echo "GHC_VERSION=${GHC_VERSION}"     >> "${GITHUB_ENV}"
           echo "CABAL_VERSION=${CABAL_VERSION}" >> "${GITHUB_ENV}"
-          echo "STACK_VERSION=${STACK_VERSION}" >> "${GITHUB_ENV}"
+
+      - name: Check cabal file
+        run: cabal check
 
       - 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-
-
-      - name: Restore %APPDATA%\stack, %LOCALAPPDATA%\Programs\stack cache (Windows)
-        uses: actions/cache/restore@v3
-        id: cache-restore-stack-global-windows
-        if: runner.os == 'Windows'
+          cabal configure --enable-tests --enable-benchmarks --disable-documentation
+          cabal build --dry-run
+
+      - name: Restore cached dependencies
+        uses: actions/cache/restore@5a3ec84eff668545956fd18022155c47e93e2684 # v4.2.3
+        id: cache
+        env:
+          key: ${{ runner.os }}-ghc-${{ steps.setup.outputs.ghc-version }}-cabal-${{ steps.setup.outputs.cabal-version }}
         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') }}
-          restore-keys: |
-            ${{ runner.os }}-ghc-${{ env.GHC_VERSION }}-stack-${{ env.STACK_VERSION }}-stack-global-
+          path: ${{ steps.setup.outputs.cabal-store }}
+          key: ${{ env.key }}-plan-${{ hashFiles('**/plan.json') }}
+          restore-keys: ${{ env.key }}-
 
       - name: Build dependencies
-        run: stack build --only-dependencies
+        run: cabal build --only-dependencies
 
-      - name: Build the package
-        run: stack build
-
-      - 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'
+      - name: Save cached dependencies
+        uses: actions/cache/save@5a3ec84eff668545956fd18022155c47e93e2684 # v4.2.3
+        if: steps.cache.outputs.cache-hit != 'true'
         with:
-          path: |
-            ~\AppData\Roaming\stack
-            ~\AppData\Local\Programs\stack    
-          key: ${{ steps.cache-restore-stack-global-windows.outputs.cache-primary-key }}
-  
-      - name: Save ~/.stack cache (Unix)
-        uses: actions/cache/save@v3
-        id: cache-save-stack-global
-        if: (runner.os == 'Linux' || runner.os == 'macOS') 
-            && steps.cache-restore-stack-global-unix.outputs.cache-hit != 'true'
-        with:
-          path: ~/.stack
-          key: ${{ steps.cache-restore-stack-global-unix.outputs.cache-primary-key }}
+          path: ${{ steps.setup.outputs.cabal-store }}
+          key: ${{ steps.cache.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/.gitignore b/.gitignore
index c8e2cad..3b65193 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,5 @@
 *~
 dist-*/
 .vscode/*
+.direnv/*
+.envrc
diff --git a/ChangeLog.md b/ChangeLog.md
index 5e6fb45..62325e2 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -2,6 +2,9 @@
 
 ## Unreleased changes
 
+- Use Hspec for tests
+- Add nix flake
+
 ## [v0.2.0.0](https://github.com/rasheedja/LPPaver/tree/v0.2.0.0)
 
 - Setup CI
diff --git a/flake.lock b/flake.lock
new file mode 100644
index 0000000..7d8c09d
--- /dev/null
+++ b/flake.lock
@@ -0,0 +1,633 @@
+{
+  "nodes": {
+    "HTTP": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1451647621,
+        "narHash": "sha256-oHIyw3x0iKBexEo49YeUDV1k74ZtyYKGR2gNJXXRxts=",
+        "owner": "phadej",
+        "repo": "HTTP",
+        "rev": "9bc0996d412fef1787449d841277ef663ad9a915",
+        "type": "github"
+      },
+      "original": {
+        "owner": "phadej",
+        "repo": "HTTP",
+        "type": "github"
+      }
+    },
+    "cabal-32": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1603716527,
+        "narHash": "sha256-X0TFfdD4KZpwl0Zr6x+PLxUt/VyKQfX7ylXHdmZIL+w=",
+        "owner": "haskell",
+        "repo": "cabal",
+        "rev": "48bf10787e27364730dd37a42b603cee8d6af7ee",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "3.2",
+        "repo": "cabal",
+        "type": "github"
+      }
+    },
+    "cabal-34": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1645834128,
+        "narHash": "sha256-wG3d+dOt14z8+ydz4SL7pwGfe7SiimxcD/LOuPCV6xM=",
+        "owner": "haskell",
+        "repo": "cabal",
+        "rev": "5ff598c67f53f7c4f48e31d722ba37172230c462",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "3.4",
+        "repo": "cabal",
+        "type": "github"
+      }
+    },
+    "cabal-36": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1669081697,
+        "narHash": "sha256-I5or+V7LZvMxfbYgZATU4awzkicBwwok4mVoje+sGmU=",
+        "owner": "haskell",
+        "repo": "cabal",
+        "rev": "8fd619e33d34924a94e691c5fea2c42f0fc7f144",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "3.6",
+        "repo": "cabal",
+        "type": "github"
+      }
+    },
+    "cardano-shell": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1608537748,
+        "narHash": "sha256-PulY1GfiMgKVnBci3ex4ptk2UNYMXqGjJOxcPy2KYT4=",
+        "owner": "input-output-hk",
+        "repo": "cardano-shell",
+        "rev": "9392c75087cb9a3d453998f4230930dea3a95725",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "repo": "cardano-shell",
+        "type": "github"
+      }
+    },
+    "flake-compat": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1672831974,
+        "narHash": "sha256-z9k3MfslLjWQfnjBtEtJZdq3H7kyi2kQtUThfTgdRk0=",
+        "owner": "input-output-hk",
+        "repo": "flake-compat",
+        "rev": "45f2638735f8cdc40fe302742b79f248d23eb368",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "ref": "hkm/gitlab-fix",
+        "repo": "flake-compat",
+        "type": "github"
+      }
+    },
+    "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"
+      }
+    },
+    "hackage": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1758846310,
+        "narHash": "sha256-kVnn9TScof8n41p7LqwvBvoLlfFhLDkjrP+aOAhmQ9k=",
+        "owner": "input-output-hk",
+        "repo": "hackage.nix",
+        "rev": "173aca690d454916a2d1ab5a7d13b593240fa0f5",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "repo": "hackage.nix",
+        "type": "github"
+      }
+    },
+    "hackage-for-stackage": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1758846300,
+        "narHash": "sha256-uS0e51ny5rGdI5HiOttTYMjGyOqBSoraXDWCY7gFc9g=",
+        "owner": "input-output-hk",
+        "repo": "hackage.nix",
+        "rev": "813f87b29c01a70bf479ff7c72b240d7d6a3fe16",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "ref": "for-stackage",
+        "repo": "hackage.nix",
+        "type": "github"
+      }
+    },
+    "hackage-internal": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1750307553,
+        "narHash": "sha256-iiafNoeLHwlSLQTyvy8nPe2t6g5AV4PPcpMeH/2/DLs=",
+        "owner": "input-output-hk",
+        "repo": "hackage.nix",
+        "rev": "f7867baa8817fab296528f4a4ec39d1c7c4da4f3",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "repo": "hackage.nix",
+        "type": "github"
+      }
+    },
+    "haskellNix": {
+      "inputs": {
+        "HTTP": "HTTP",
+        "cabal-32": "cabal-32",
+        "cabal-34": "cabal-34",
+        "cabal-36": "cabal-36",
+        "cardano-shell": "cardano-shell",
+        "flake-compat": "flake-compat",
+        "hackage": "hackage",
+        "hackage-for-stackage": "hackage-for-stackage",
+        "hackage-internal": "hackage-internal",
+        "hls": "hls",
+        "hls-1.10": "hls-1.10",
+        "hls-2.0": "hls-2.0",
+        "hls-2.10": "hls-2.10",
+        "hls-2.11": "hls-2.11",
+        "hls-2.2": "hls-2.2",
+        "hls-2.3": "hls-2.3",
+        "hls-2.4": "hls-2.4",
+        "hls-2.5": "hls-2.5",
+        "hls-2.6": "hls-2.6",
+        "hls-2.7": "hls-2.7",
+        "hls-2.8": "hls-2.8",
+        "hls-2.9": "hls-2.9",
+        "hpc-coveralls": "hpc-coveralls",
+        "iserv-proxy": "iserv-proxy",
+        "nixpkgs": [
+          "haskellNix",
+          "nixpkgs-unstable"
+        ],
+        "nixpkgs-2305": "nixpkgs-2305",
+        "nixpkgs-2311": "nixpkgs-2311",
+        "nixpkgs-2405": "nixpkgs-2405",
+        "nixpkgs-2411": "nixpkgs-2411",
+        "nixpkgs-2505": "nixpkgs-2505",
+        "nixpkgs-unstable": "nixpkgs-unstable",
+        "old-ghc-nix": "old-ghc-nix",
+        "stackage": "stackage"
+      },
+      "locked": {
+        "lastModified": 1758847890,
+        "narHash": "sha256-rGX7RF8Au5ZJJSqlQivsl4seyEslI/K3OnEC9ulLwNM=",
+        "owner": "input-output-hk",
+        "repo": "haskell.nix",
+        "rev": "46abef90b4101ff9253a574cf6fbdc74b78a5863",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "repo": "haskell.nix",
+        "type": "github"
+      }
+    },
+    "hls": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1741604408,
+        "narHash": "sha256-tuq3+Ip70yu89GswZ7DSINBpwRprnWnl6xDYnS4GOsc=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "682d6894c94087da5e566771f25311c47e145359",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-1.10": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1680000865,
+        "narHash": "sha256-rc7iiUAcrHxwRM/s0ErEsSPxOR3u8t7DvFeWlMycWgo=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "b08691db779f7a35ff322b71e72a12f6e3376fd9",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "1.10.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.0": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1687698105,
+        "narHash": "sha256-OHXlgRzs/kuJH8q7Sxh507H+0Rb8b7VOiPAjcY9sM1k=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "783905f211ac63edf982dd1889c671653327e441",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.0.0.1",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.10": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1743069404,
+        "narHash": "sha256-q4kDFyJDDeoGqfEtrZRx4iqMVEC2MOzCToWsFY+TOzY=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "2318c61db3a01e03700bd4b05665662929b7fe8b",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.10.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.11": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1747306193,
+        "narHash": "sha256-/MmtpF8+FyQlwfKHqHK05BdsxC9LHV70d/FiMM7pzBM=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "46ef4523ea4949f47f6d2752476239f1c6d806fe",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.11.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.2": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1693064058,
+        "narHash": "sha256-8DGIyz5GjuCFmohY6Fa79hHA/p1iIqubfJUTGQElbNk=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "b30f4b6cf5822f3112c35d14a0cba51f3fe23b85",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.2.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.3": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1695910642,
+        "narHash": "sha256-tR58doOs3DncFehHwCLczJgntyG/zlsSd7DgDgMPOkI=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "458ccdb55c9ea22cd5d13ec3051aaefb295321be",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.3.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.4": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1699862708,
+        "narHash": "sha256-YHXSkdz53zd0fYGIYOgLt6HrA0eaRJi9mXVqDgmvrjk=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "54507ef7e85fa8e9d0eb9a669832a3287ffccd57",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.4.0.1",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.5": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1701080174,
+        "narHash": "sha256-fyiR9TaHGJIIR0UmcCb73Xv9TJq3ht2ioxQ2mT7kVdc=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "27f8c3d3892e38edaef5bea3870161815c4d014c",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.5.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.6": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1705325287,
+        "narHash": "sha256-+P87oLdlPyMw8Mgoul7HMWdEvWP/fNlo8jyNtwME8E8=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "6e0b342fa0327e628610f2711f8c3e4eaaa08b1e",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.6.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.7": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1708965829,
+        "narHash": "sha256-LfJ+TBcBFq/XKoiNI7pc4VoHg4WmuzsFxYJ3Fu+Jf+M=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "50322b0a4aefb27adc5ec42f5055aaa8f8e38001",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.7.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.8": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1715153580,
+        "narHash": "sha256-Vi/iUt2pWyUJlo9VrYgTcbRviWE0cFO6rmGi9rmALw0=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "dd1be1beb16700de59e0d6801957290bcf956a0a",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.8.0.0",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hls-2.9": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1719993701,
+        "narHash": "sha256-wy348++MiMm/xwtI9M3vVpqj2qfGgnDcZIGXw8sF1sA=",
+        "owner": "haskell",
+        "repo": "haskell-language-server",
+        "rev": "90319a7e62ab93ab65a95f8f2bcf537e34dae76a",
+        "type": "github"
+      },
+      "original": {
+        "owner": "haskell",
+        "ref": "2.9.0.1",
+        "repo": "haskell-language-server",
+        "type": "github"
+      }
+    },
+    "hpc-coveralls": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1607498076,
+        "narHash": "sha256-8uqsEtivphgZWYeUo5RDUhp6bO9j2vaaProQxHBltQk=",
+        "owner": "sevanspowell",
+        "repo": "hpc-coveralls",
+        "rev": "14df0f7d229f4cd2e79f8eabb1a740097fdfa430",
+        "type": "github"
+      },
+      "original": {
+        "owner": "sevanspowell",
+        "repo": "hpc-coveralls",
+        "type": "github"
+      }
+    },
+    "iserv-proxy": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1755243078,
+        "narHash": "sha256-GLbl1YaohKdpzZVJFRdcI1O1oE3F3uBer4lFv3Yy0l8=",
+        "owner": "stable-haskell",
+        "repo": "iserv-proxy",
+        "rev": "150605195cb7183a6fb7bed82f23fedf37c6f52a",
+        "type": "github"
+      },
+      "original": {
+        "owner": "stable-haskell",
+        "ref": "iserv-syms",
+        "repo": "iserv-proxy",
+        "type": "github"
+      }
+    },
+    "nixpkgs-2305": {
+      "locked": {
+        "lastModified": 1705033721,
+        "narHash": "sha256-K5eJHmL1/kev6WuqyqqbS1cdNnSidIZ3jeqJ7GbrYnQ=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "a1982c92d8980a0114372973cbdfe0a307f1bdea",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-23.05-darwin",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-2311": {
+      "locked": {
+        "lastModified": 1719957072,
+        "narHash": "sha256-gvFhEf5nszouwLAkT9nWsDzocUTqLWHuL++dvNjMp9I=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "7144d6241f02d171d25fba3edeaf15e0f2592105",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-23.11-darwin",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-2405": {
+      "locked": {
+        "lastModified": 1735564410,
+        "narHash": "sha256-HB/FA0+1gpSs8+/boEavrGJH+Eq08/R2wWNph1sM1Dg=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "1e7a8f391f1a490460760065fa0630b5520f9cf8",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-24.05-darwin",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-2411": {
+      "locked": {
+        "lastModified": 1748037224,
+        "narHash": "sha256-92vihpZr6dwEMV6g98M5kHZIttrWahb9iRPBm1atcPk=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "f09dede81861f3a83f7f06641ead34f02f37597f",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-24.11-darwin",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-2505": {
+      "locked": {
+        "lastModified": 1754477006,
+        "narHash": "sha256-suIgZZHXdb4ca9nN4MIcmdjeN+ZWsTwCtYAG4HExqAo=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "4896699973299bffae27d0d9828226983544d9e9",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-25.05-darwin",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "nixpkgs-unstable": {
+      "locked": {
+        "lastModified": 1754393734,
+        "narHash": "sha256-fbnmAwTQkuXHKBlcL5Nq1sMAzd3GFqCOQgEQw6Hy0Ak=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "a683adc19ff5228af548c6539dbc3440509bfed3",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixpkgs-unstable",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "old-ghc-nix": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1631092763,
+        "narHash": "sha256-sIKgO+z7tj4lw3u6oBZxqIhDrzSkvpHtv0Kki+lh9Fg=",
+        "owner": "angerman",
+        "repo": "old-ghc-nix",
+        "rev": "af48a7a7353e418119b6dfe3cd1463a657f342b8",
+        "type": "github"
+      },
+      "original": {
+        "owner": "angerman",
+        "ref": "master",
+        "repo": "old-ghc-nix",
+        "type": "github"
+      }
+    },
+    "root": {
+      "inputs": {
+        "flake-utils": "flake-utils",
+        "haskellNix": "haskellNix",
+        "nixpkgs": [
+          "haskellNix",
+          "nixpkgs-2505"
+        ]
+      }
+    },
+    "stackage": {
+      "flake": false,
+      "locked": {
+        "lastModified": 1758845522,
+        "narHash": "sha256-SgkvlWF9a+Qrkn791ZOiUVt3wuZXRJ06YjpTZMRy+R8=",
+        "owner": "input-output-hk",
+        "repo": "stackage.nix",
+        "rev": "e2f097d435e38fb6e649efa4a95e214a506a1da5",
+        "type": "github"
+      },
+      "original": {
+        "owner": "input-output-hk",
+        "repo": "stackage.nix",
+        "type": "github"
+      }
+    },
+    "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..69be272
--- /dev/null
+++ b/flake.nix
@@ -0,0 +1,45 @@
+{
+  inputs = {
+    haskellNix.url = "github:input-output-hk/haskell.nix";
+    nixpkgs.follows = "haskellNix/nixpkgs-2505";
+    flake-utils.url = "github:numtide/flake-utils";
+  };
+
+  outputs = { self, nixpkgs, flake-utils, haskellNix }:
+    flake-utils.lib.eachDefaultSystem (system:
+      let
+        overlays = [ haskellNix.overlay ];
+        pkgs = import nixpkgs {
+          inherit system overlays;
+          inherit (haskellNix) config;
+        };
+
+        project = pkgs.haskell-nix.cabalProject' {
+          src = ./.;
+          compiler-nix-name = "ghc967";
+          shell = {
+            tools = {
+              cabal = "3.16.0.0";
+              hlint = "3.8";
+              haskell-language-server = "2.11.0.0";
+              fourmolu = "0.17.0.0";
+            };
+            buildInputs = with pkgs; [
+              # system dependencies go here
+            ];
+          };
+        };
+
+        flake = project.flake {};
+      in flake);
+  
+  # --- Flake Local Nix Configuration ---
+  nixConfig = {
+    # This sets the flake to use the IOG nix cache.
+    # Nix should ask for permission before using it,
+    # but remove it here if you do not want it to.
+    extra-substituters = ["https://cache.iog.io"];
+    extra-trusted-public-keys = ["hydra.iohk.io:f/Ea+s+dFdN+3Y/G+FDgSq+a5NEWhJGzdjvKNGv0/EQ="];
+    allow-import-from-derivation = "true";
+  };
+}
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 <https://github.com/rasheedja/simplex-method#readme>
-
-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..f3e9673 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,77 @@ 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.16.0 && < 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
       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.16.0 && < 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/Solver/TwoPhase.hs b/src/Linear/Simplex/Solver/TwoPhase.hs
index c7dfe83..cac4ef8 100644
--- a/src/Linear/Simplex/Solver/TwoPhase.hs
+++ b/src/Linear/Simplex/Solver/TwoPhase.hs
@@ -6,11 +6,32 @@
 -- 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'' performs both phases with variable domain support.
+module Linear.Simplex.Solver.TwoPhase 
+  ( findFeasibleSolution
+  , optimizeFeasibleSystem
+  , twoPhaseSimplex
+  , twoPhaseSimplex'
+  -- Internal functions exported for testing
+  , preprocess
+  , postprocess
+  , computeObjective
+  , collectAllVars
+  , generateTransform
+  , getTransform
+  , applyTransforms
+  , applyTransform
+  , applyShiftToObjective
+  , applyShiftToConstraint
+  , applySplitToObjective
+  , applySplitToConstraint
+  , unapplyTransforms
+  , unapplyTransform
+  ) where
 
 import Prelude hiding (EQ)
 
@@ -24,6 +45,8 @@ import qualified Data.Map as M
 import Data.Maybe (fromJust, fromMaybe, mapMaybe)
 import Data.Ratio (denominator, numerator, (%))
 import qualified Data.Text as Text
+import Data.Set (Set)
+import qualified Data.Set as Set
 import GHC.Real (Ratio)
 import Linear.Simplex.Types
 import Linear.Simplex.Util
@@ -403,6 +426,240 @@ twoPhaseSimplex objFunction unsimplifiedSystem = do
       logMsg LevelInfo $ "twoPhaseSimplex: Phase 1 gives infeasible result for " <> showT unsimplifiedSystem
       pure Nothing
 
+-- | 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 Result contains variable values and objective value in the original space.
+-- TODO: use this as twoPhaseSimplex, add instructions in CHANGELOG for old users
+twoPhaseSimplex' :: (MonadIO m, MonadLogger m) => VarDomainMap -> ObjectiveFunction -> [PolyConstraint] -> m (Maybe Result)
+twoPhaseSimplex' domainMap objFunction constraints = do
+  logMsg LevelInfo $
+    "twoPhaseSimplex': Solving system with domain map " <> showT domainMap
+  let (transformedObj, transformedConstraints, transforms) = preprocess objFunction domainMap constraints 
+  logMsg LevelInfo $
+    "twoPhaseSimplex': Applied transforms " <> showT transforms
+      <> "; Transformed objective: " <> showT transformedObj
+      <> "; Transformed constraints: " <> showT transformedConstraints
+  mResult <- twoPhaseSimplex transformedObj transformedConstraints
+  case mResult of
+    Nothing -> do
+      logMsg LevelInfo "twoPhaseSimplex': No solution found"
+      pure Nothing
+    Just result -> do
+      let finalResult = postprocess objFunction transforms result
+      logMsg LevelInfo $
+        "twoPhaseSimplex': Postprocessed result: " <> showT finalResult
+      pure (Just finalResult)
+
+-- | Postprocess the result by unapplying variable transformations and computing
+-- the objective value in the original space.
+postprocess :: ObjectiveFunction -> [VarTransform] -> Result -> Result
+postprocess objFunction transforms result =
+  let -- First unapply transforms to get variable values in original space
+      unappliedResult = unapplyTransforms transforms result
+      -- Then compute the objective value using the original objective function
+      objVal = computeObjective objFunction unappliedResult.varValMap
+      -- Update the objective value in the result
+      finalVarValMap = M.insert unappliedResult.objectiveVar objVal unappliedResult.varValMap
+  in unappliedResult { varValMap = finalVarValMap }
+
+-- | 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 objective, constraints, and the list of transforms applied.
+preprocess :: ObjectiveFunction 
+           -> VarDomainMap 
+           -> [PolyConstraint]
+           -> (ObjectiveFunction, [PolyConstraint], [VarTransform])
+preprocess objFunction (VarDomainMap domainMap) constraints =
+  let -- Collect all variables in the system
+      allVars = collectAllVars objFunction 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
+      (transformedObj, transformedConstraints) = applyTransforms transforms objFunction constraints
+  in (transformedObj, transformedConstraints, transforms)
+
+-- | Collect all variables appearing in the objective function and constraints
+collectAllVars :: ObjectiveFunction -> [PolyConstraint] -> Set Var
+collectAllVars objFunction constraints =
+  let objVars = case objFunction of
+        Max m -> M.keysSet m
+        Min m -> M.keysSet m
+      constraintVars = Set.unions $ map getConstraintVars constraints
+  in Set.union objVars constraintVars
+  where
+    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
+  in case getTransform nextFreshVar var domain of
+       Nothing -> (transforms, nextFreshVar)
+       Just t@(AddLowerBound {}) -> (t : transforms, nextFreshVar)
+       Just t@(Shift {}) -> (t : transforms, nextFreshVar + 1)
+       Just t@(Split {}) -> (t : transforms, nextFreshVar + 2)
+
+-- | Determine what transform (if any) is needed for a variable given its domain.
+getTransform :: Var -> Var -> VarDomain -> Maybe VarTransform
+getTransform nextFreshVar var domain =
+  case domain of
+    NonNegative -> Nothing
+    
+    LowerBound l 
+      | l == 0 -> Nothing
+      | l > 0  -> Just $ AddLowerBound var l
+      | otherwise -> Just $ Shift var nextFreshVar l  -- l < 0, need to shift
+    
+    Unbounded ->
+      Just $ Split var nextFreshVar (nextFreshVar + 1)
+
+-- | 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)
+    
+    -- 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 the result back to original variables.
+unapplyTransforms :: [VarTransform] -> Result -> Result
+unapplyTransforms transforms result =
+  -- Apply transforms in reverse order (since we applied them with foldr)
+  foldl (flip unapplyTransform) result transforms
+
+-- | Unapply a single transform to convert result back to original variable.
+unapplyTransform :: VarTransform -> Result -> Result
+unapplyTransform transform result@(Result {varValMap = valMap, ..}) =
+  case transform of
+    -- AddLowerBound: No variable substitution was done, nothing to unapply
+    AddLowerBound {} -> result
+    
+    -- 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
+          newMap = M.insert origVar origVal (M.delete shiftedVar valMap)
+      in result { varValMap = newMap }
+    
+    -- 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
+          newMap = M.insert origVar origVal (M.delete posVar (M.delete negVar valMap))
+      in result { varValMap = newMap }
+
 -- | Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'.
 simplexPivot :: (MonadIO m, MonadLogger m) => PivotObjective -> Dict -> m (Maybe Dict)
 simplexPivot objective@(PivotObjective {variable = objectiveVar, function = objectiveFunc, constant = objectiveConstant}) dictionary = do
diff --git a/src/Linear/Simplex/Types.hs b/src/Linear/Simplex/Types.hs
index 15e5d1f..3d2ea63 100644
--- a/src/Linear/Simplex/Types.hs
+++ b/src/Linear/Simplex/Types.hs
@@ -121,3 +121,39 @@ data PivotObjective = PivotObjective
   , constant :: SimplexNum
   }
   deriving (Show, Read, Eq, Generic)
+
+-- | Domain specification for a variable's lower bound.
+-- Note: This only concerns lower bounds. Upper bounds are handled via constraints.
+-- Variables not in the VarDomainMap are assumed to be Unbounded.
+data VarDomain 
+  = NonNegative           -- ^ var >= 0 (standard simplex assumption, no transformation needed)
+  | LowerBound SimplexNum -- ^ var >= L for some L (if L < 0: shift, if L > 0: add constraint)
+  | Unbounded             -- ^ No lower bound (split into difference of two non-negative vars)
+  -- TODO: Upperbound can still be useful, can negate it to get a loewr bound, can add it to the constraints
+  deriving stock (Show, Read, Eq, Generic)
+
+-- | 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 
+  = AddLowerBound 
+      { var :: !Var
+      , bound :: !SimplexNum 
+      } -- ^ var >= bound where bound > 0. Adds GEQ constraint to system.
+  | Shift 
+      { originalVar :: !Var
+      , shiftedVar :: !Var  
+      , shiftBy :: !SimplexNum 
+      } -- ^ originalVar = shiftedVar + shiftBy, where shiftBy < 0. After solving: originalVar = shiftedVar + shiftBy
+  | Split 
+      { originalVar :: !Var
+      , posVar :: !Var
+      , negVar :: !Var 
+      } -- ^ originalVar = posVar - negVar, both posVar and negVar >= 0
+  deriving stock (Show, Read, Eq, Generic)
+                  
+
+
diff --git a/stack.yaml b/stack.yaml
deleted file mode 100644
index eab5650..0000000
--- a/stack.yaml
+++ /dev/null
@@ -1,68 +0,0 @@
-# This file was automatically generated by 'stack init'
-#
-# Some commonly used options have been documented as comments in this file.
-# For advanced use and comprehensive documentation of the format, please see:
-# https://docs.haskellstack.org/en/stable/yaml_configuration/
-
-# Resolver to choose a 'specific' stackage snapshot or a compiler version.
-# A snapshot resolver dictates the compiler version and the set of packages
-# to be used for project dependencies. For example:
-#
-# resolver: lts-3.5
-# resolver: nightly-2015-09-21
-# resolver: ghc-7.10.2
-#
-# The location of a snapshot can be provided as a file or url. Stack assumes
-# a snapshot provided as a file might change, whereas a url resource does not.
-#
-# resolver: ./custom-snapshot.yaml
-# resolver: https://example.com/snapshots/2018-01-01.yaml
-resolver: lts-21.22
-
-# User packages to be built.
-# Various formats can be used as shown in the example below.
-#
-# packages:
-# - some-directory
-# - https://example.com/foo/bar/baz-0.0.2.tar.gz
-#   subdirs:
-#   - auto-update
-#   - wai
-packages:
-- .
-# Dependency packages to be pulled from upstream that are not in the resolver.
-# These entries can reference officially published versions as well as
-# forks / in-progress versions pinned to a git hash. For example:
-#
-# extra-deps:
-# - acme-missiles-0.3
-# - git: https://github.com/commercialhaskell/stack.git
-#   commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a
-#
-# extra-deps: {}
-
-# Override default flag values for local packages and extra-deps
-# flags: {}
-
-# Extra package databases containing global packages
-# extra-package-dbs: []
-
-# Control whether we use the GHC we find on the path
-# system-ghc: true
-#
-# Require a specific version of stack, using version ranges
-# require-stack-version: -any # Default
-# require-stack-version: ">=2.5"
-#
-# Override the architecture used by stack, especially useful on Windows
-# arch: i386
-# arch: x86_64
-#
-# Extra directories used by stack for building
-# extra-include-dirs: [/path/to/dir]
-# extra-lib-dirs: [/path/to/dir]
-#
-# Allow a newer minor version of GHC than the snapshot specifies
-# compiler-check: newer-minor
-
-system-ghc: true
diff --git a/stack.yaml.lock b/stack.yaml.lock
deleted file mode 100644
index e8d3cc7..0000000
--- a/stack.yaml.lock
+++ /dev/null
@@ -1,12 +0,0 @@
-# This file was autogenerated by Stack.
-# You should not edit this file by hand.
-# For more information, please see the documentation at:
-#   https://docs.haskellstack.org/en/stable/lock_files
-
-packages: []
-snapshots:
-- completed:
-    sha256: afd5ba64ab602cabc2d3942d3d7e7dd6311bc626dcb415b901eaf576cb62f0ea
-    size: 640060
-    url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/21/22.yaml
-  original: lts-21.22
diff --git a/test/Linear/Simplex/Solver/TwoPhaseSpec.hs b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs
new file mode 100644
index 0000000..8d0ca39
--- /dev/null
+++ b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs
@@ -0,0 +1,1823 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+module Linear.Simplex.Solver.TwoPhaseSpec where
+
+import Prelude hiding (EQ)
+
+import Control.Monad.IO.Class
+import Control.Monad.Logger
+import qualified Data.Map as M
+import qualified Data.Set as Set
+import Data.Ratio
+
+import Text.InterpolatedString.Perl6
+
+import Test.Hspec
+import Test.Hspec.Expectations.Contrib (annotate)
+import Test.QuickCheck hiding (Result)
+import qualified Linear.Simplex.Types as T
+
+import Linear.Simplex.Prettify
+import Linear.Simplex.Solver.TwoPhase
+import Linear.Simplex.Types hiding (NonNegative)
+import Linear.Simplex.Util
+
+-- | Helper to run a test case and check result
+runTest :: (ObjectiveFunction, [PolyConstraint]) -> Maybe Result -> IO ()
+runTest (obj, constraints) expectedResult = do
+  actualResult <-
+    runStdoutLoggingT $
+      filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+        twoPhaseSimplex obj constraints
+  let prettyObj = prettyShowObjectiveFunction obj
+      prettyConstraints = map prettyShowPolyConstraint constraints
+      expectedObjVal = extractObjectiveValue expectedResult
+      actualObjVal = extractObjectiveValue actualResult
+      -- HACK: Verify NonNegative twoPhaseSimplex' NonNegative == twoPhaseSimplex
+      allVars = collectAllVars obj constraints
+      domainMap = VarDomainMap $ M.fromSet (const T.NonNegative) allVars
+  actualResult' <-
+    runStdoutLoggingT $
+    filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+      twoPhaseSimplex' domainMap obj constraints
+  let actualObjVal' = extractObjectiveValue actualResult'
+  annotate
+    [qc|
+
+Objective Function (Non-prettified): {obj}
+Constraints        (Non-prettified): {constraints}
+====================================
+Objective Function     (Prettified): {prettyObj}
+Constraints            (Prettified): {prettyConstraints}
+====================================
+Expected Solution            (Full): {expectedResult}
+Actual Solution              (Full): {actualResult}
+Expected Solution       (Objective): {expectedObjVal}
+Actual Solution         (Objective): {actualObjVal}
+====================================
+Actual Solution'             (Full): {actualResult'}
+Actual Solution'        (Objective): {actualObjVal'}
+    |]
+    $ do
+      actualResult `shouldBe` expectedResult
+      -- TODO: worth removing twoPhaseSimplex?
+      actualResult' `shouldBe` expectedResult
+
+
+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 testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 29), (1, 3), (2, 4)])))
+
+      it "Min 3x₁ + 5x₂ with LEQ constraints: obj=0" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 0)])))
+
+      it "Max 3x₁ + 5x₂ with GEQ constraints: infeasible" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Min 3x₁ + 5x₂ with GEQ constraints: obj=237/7, x₁=24/7, x₂=33/7" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 11 (M.fromList [(11, 237 % 7), (1, 24 % 7), (2, 33 % 7)])))
+
+    -- 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 testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase (Just (Result 9 (M.fromList [(9, 3 % 5), (2, 14 % 5), (3, 17 % 5)])))
+
+      it "Min x₁ - x₂ + x₃ with LEQ constraints: infeasible" $ do
+        let testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Max x₁ - x₂ + x₃ with GEQ constraints: obj=1, x₁=3, x₂=2" $ do
+        let testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase (Just (Result 8 (M.fromList [(8, 1), (2, 2), (1, 3)])))
+
+      it "Min x₁ - x₂ + x₃ with GEQ constraints: obj=-1/4, x₁=17/4, x₂=9/2" $ do
+        let testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase (Just (Result 8 (M.fromList [(8, (-1) % 4), (2, 9 % 2), (1, 17 % 4)])))
+
+    -- 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 testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 5), (3, 2), (4, 1)])))
+
+      it "Max x₁ + x₂ + 2x₃ + x₄ with EQ constraints: obj=8, x₁=2, x₂=6" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 8), (1, 2), (2, 6)])))
+
+    -- 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 testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 8 (M.fromList [(8, 20), (4, 16), (3, 6)])))
+
+      it "Min -2x₃ + 2x₄ + x₅ with EQ constraints: obj=6, x₄=2, x₅=2" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 8 (M.fromList [(8, 6), (4, 2), (5, 2)])))
+
+    -- 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 testCase =
+              ( Max (M.fromList [(1, 2), (2, 1)])
+              , [ LEQ (M.fromList [(1, 4), (2, 1)]) 150
+                , LEQ (M.fromList [(1, 2), (2, -3)]) (-40)
+                ]
+              )
+        runTest testCase (Just (Result 6 (M.fromList [(6, 150), (2, 150)])))
+
+      it "Min 2x₁ + x₂: obj=40/3, x₂=40/3" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 2), (2, 1)])
+              , [ LEQ (M.fromList [(1, 4), (2, 1)]) 150
+                , LEQ (M.fromList [(1, 2), (2, -3)]) (-40)
+                ]
+              )
+        runTest testCase (Just (Result 6 (M.fromList [(6, 40 % 3), (2, 40 % 3)])))
+
+      it "Max 2x₁ + x₂ with GEQ constraints: infeasible" $ do
+        let testCase =
+              ( Max (M.fromList [(1, 2), (2, 1)])
+              , [ GEQ (M.fromList [(1, 4), (2, 1)]) 150
+                , GEQ (M.fromList [(1, 2), (2, -3)]) (-40)
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Min 2x₁ + x₂ with GEQ constraints: obj=75, x₁=75/2" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 2), (2, 1)])
+              , [ GEQ (M.fromList [(1, 4), (2, 1)]) 150
+                , GEQ (M.fromList [(1, 2), (2, -3)]) (-40)
+                ]
+              )
+        runTest testCase (Just (Result 6 (M.fromList [(6, 75), (1, 75 % 2)])))
+
+    -- 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 testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, (-120)), (1, 20)])))
+
+      it "Max -6x₁ - 4x₂ + 2x₃: obj=10, x₃=5" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 10), (3, 5)])))
+
+      it "Min -6x₁ - 4x₂ + 2x₃ with GEQ constraints: infeasible" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Max -6x₁ - 4x₂ + 2x₃ with GEQ constraints: infeasible" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase Nothing
+
+    -- 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 testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 250), (2, 50)])))
+
+      it "Min 3x₁ + 5x₂ + 2x₃: obj=0" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, 0)])))
+
+      it "Max 3x₁ + 5x₂ + 2x₃ with GEQ constraints: infeasible" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Min 3x₁ + 5x₂ + 2x₃ with GEQ constraints: obj=300, x₃=150" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 10 (M.fromList [(10, 300), (3, 150)])))
+
+    describe "Simple single/two variable tests" $ do
+      it "Max x₁ with x₁ <= 15: obj=15, x₁=15" $ do
+        let testCase =
+              ( Max (M.fromList [(1, 1)])
+              , [ LEQ (M.fromList [(1, 1)]) 15
+                ]
+              )
+        runTest testCase (Just (Result 3 (M.fromList [(3, 15), (1, 15)])))
+
+      it "Max 2x₁ with mixed constraints: obj=20, x₁=10, x₂=10" $ do
+        let testCase =
+              ( Max (M.fromList [(1, 2)])
+              , [ LEQ (M.fromList [(1, 2)]) 20
+                , GEQ (M.fromList [(2, 1)]) 10
+                ]
+              )
+        runTest testCase (Just (Result 6 (M.fromList [(6, 20), (1, 10), (2, 10)])))
+
+      it "Min x₁ with x₁ <= 15: obj=0" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 1)])
+              , [ LEQ (M.fromList [(1, 1)]) 15
+                ]
+              )
+        runTest testCase (Just (Result 3 (M.fromList [(3, 0)])))
+
+      it "Min 2x₁ with mixed constraints: obj=0, x₂=10" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 2)])
+              , [ LEQ (M.fromList [(1, 2)]) 20
+                , GEQ (M.fromList [(2, 1)]) 10
+                ]
+              )
+        runTest testCase (Just (Result 6 (M.fromList [(6, 0), (2, 10)])))
+
+    describe "Infeasibility tests" $ do
+      it "Conflicting bounds x₁ <= 15 and x₁ >= 15.01: infeasible" $ do
+        let testCase =
+              ( Max (M.fromList [(1, 1)])
+              , [ LEQ (M.fromList [(1, 1)]) 15
+                , GEQ (M.fromList [(1, 1)]) 15.01
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Conflicting bounds with additional constraint: infeasible" $ do
+        let testCase =
+              ( Max (M.fromList [(1, 1)])
+              , [ LEQ (M.fromList [(1, 1)]) 15
+                , GEQ (M.fromList [(1, 1)]) 15.01
+                , GEQ (M.fromList [(2, 1)]) 10
+                ]
+              )
+        runTest testCase Nothing
+
+      it "Min x₁ with duplicate GEQ constraints: obj=0, x₂=1" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 1)])
+              , [ GEQ (M.fromList [(1, 1), (2, 1)]) 1
+                , GEQ (M.fromList [(1, 1), (2, 1)]) 1
+                ]
+              )
+        runTest testCase (Just (Result 5 (M.fromList [(2, 1 % 1), (5, 0 % 1)])))
+
+      it "Conflicting x₁+x₂ >= 2 and x₁+x₂ <= 1: infeasible" $ do
+        let testCase =
+              ( Min (M.fromList [(1, 1)])
+              , [ GEQ (M.fromList [(1, 1), (2, 1)]) 2
+                , LEQ (M.fromList [(1, 1), (2, 1)]) 1
+                ]
+              )
+        runTest testCase Nothing
+
+    describe "LEQ/GEQ reduction bug tests" $ do
+      it "testLeqGeqBugMin1: obj=3, x₁=3, x₂=3" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)])))
+
+      it "testLeqGeqBugMax1: obj=3, x₁=3, x₂=3" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)])))
+
+      it "testLeqGeqBugMin2: obj=3, x₁=3, x₂=3" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)])))
+
+      it "testLeqGeqBugMax2: obj=3, x₁=3, x₂=3" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)])))
+
+    -- 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
+        runTest (mkConstraints (Min (M.fromList [(1, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 7 % 4), (2, 5 % 2), (1, 7 % 4), (3, 0)])))
+
+      it "Max x₁: x₁=5/2, x₂=5/3" $ do
+        runTest (mkConstraints (Max (M.fromList [(1, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 3), (1, 5 % 2), (3, 0)])))
+
+      it "Min x₂: x₂=5/3" $ do
+        runTest (mkConstraints (Min (M.fromList [(2, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 5 % 3), (2, 5 % 3), (1, 5 % 2), (3, 0)])))
+
+      it "Max x₂: x₂=5/2" $ do
+        runTest (mkConstraints (Max (M.fromList [(2, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 2), (1, 5 % 2), (3, 0)])))
+
+    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
+        runTest (mkConstraints (Max (M.fromList [(1, 1)]))) Nothing
+
+      it "Min x₁: infeasible" $ do
+        runTest (mkConstraints (Min (M.fromList [(1, 1)]))) Nothing
+
+      it "Max x₂: infeasible" $ do
+        runTest (mkConstraints (Max (M.fromList [(2, 1)]))) Nothing
+
+      it "Min x₂: infeasible" $ do
+        runTest (mkConstraints (Min (M.fromList [(2, 1)]))) Nothing
+
+    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
+        runTest (mkConstraints (Max (M.fromList [(1, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 7 % 2), (2, 5 % 9), (1, 7 % 2), (3, 0)])))
+
+      it "Min x₁: x₁=17/20" $ do
+        runTest (mkConstraints (Min (M.fromList [(1, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 17 % 20), (2, 7 % 2), (1, 17 % 20), (3, 0)])))
+
+      it "Max x₂: x₂=7/2" $ do
+        runTest (mkConstraints (Max (M.fromList [(2, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 7 % 2), (2, 7 % 2), (1, 22 % 9)])))
+
+      it "Min x₂: x₂=5/9" $ do
+        runTest (mkConstraints (Min (M.fromList [(2, 1)]))) 
+                (Just (Result 12 (M.fromList [(12, 5 % 9), (2, 5 % 9), (1, 7 % 2), (3, 0)])))
+
+    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
+        runTest (mkConstraints (Max (M.fromList [(1, 1)]))) Nothing
+
+      it "Min x₁: infeasible" $ do
+        runTest (mkConstraints (Min (M.fromList [(1, 1)]))) Nothing
+
+      it "Max x₂: infeasible" $ do
+        runTest (mkConstraints (Max (M.fromList [(2, 1)]))) Nothing
+
+      it "Min x₂: infeasible" $ do
+        runTest (mkConstraints (Min (M.fromList [(2, 1)]))) Nothing
+
+    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
+        runTest (mkConstraints (Max (M.fromList [(1, 1)]))) 
+                (Just (Result 17 (M.fromList [(17, 5 % 2), (2, 45 % 22), (1, 5 % 2), (4, 0)])))
+
+      it "Min x₁: x₁=45/22" $ do
+        runTest (mkConstraints (Min (M.fromList [(1, 1)]))) 
+                (Just (Result 17 (M.fromList [(17, 45 % 22), (2, 5 % 2), (1, 45 % 22), (4, 0)])))
+
+      it "Max x₂: x₂=5/2" $ do
+        runTest (mkConstraints (Max (M.fromList [(2, 1)]))) 
+                (Just (Result 17 (M.fromList [(17, 5 % 2), (2, 5 % 2), (1, 5 % 2), (4, 0)])))
+
+      it "Min x₂: x₂=45/22" $ do
+        runTest (mkConstraints (Min (M.fromList [(2, 1)]))) 
+                (Just (Result 17 (M.fromList [(17, 45 % 22), (2, 45 % 22), (1, 5 % 2), (4, 0)])))
+
+    describe "QuickCheck-generated regression tests" $ do
+      it "testQuickCheck1: obj=-370, x₁=5/3, x₂=26" $ do
+        let testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase (Just (Result 10 (M.fromList [(10, (-370)), (2, 26), (1, 5 % 3)])))
+
+      it "testQuickCheck2: obj=-2/9, x₁=14/9, x₂=8/9" $ do
+        let testCase =
+              ( 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
+                ]
+              )
+        runTest testCase (Just (Result 8 (M.fromList [(8, (-2) % 9), (1, 14 % 9), (2, 8 % 9)])))
+
+      it "testQuickCheck3 (tests objective simplification): obj=-8, x₂=2" $ do
+        let testCase =
+              ( 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)
+                ]
+              )
+        runTest testCase (Just (Result 7 (M.fromList [(7, (-8)), (2, 2)])))
+
+  describe "twoPhaseSimplex' (with VarDomainMap)" $ do
+    it "NonNegative domain gives same result as twoPhaseSimplex" $ do
+      -- TODO: redundant if we keep the runTest hack
+      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
+            ]
+          domainMap = VarDomainMap $ M.fromList [(1, T.NonNegative), (2, T.NonNegative)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      actualResult `shouldBe` Just (Result 7 (M.fromList [(7, 29), (1, 3), (2, 4)]))
+
+    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, LowerBound (-5))]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 10
+
+    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, LowerBound (-5))]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-5)
+
+    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
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 10
+
+    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
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-10)
+
+    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, LowerBound 5)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 10
+
+    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, LowerBound 5)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 5
+
+    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, LowerBound (-2)), (2, Unbounded)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> do
+          let xVal = M.findWithDefault 0 1 result.varValMap
+              yVal = M.findWithDefault 0 2 result.varValMap
+              oVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+          (xVal + yVal) `shouldBe` 5
+          oVal `shouldBe` 5
+
+    it "LowerBound 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, LowerBound 0), (2, LowerBound 0)]
+          domainMap2 = VarDomainMap $ M.fromList [(1, T.NonNegative), (2, T.NonNegative)]
+      actualResult1 <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap1 obj constraints
+      actualResult2 <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap2 obj constraints
+      actualResult1 `shouldBe` Just (Result 7 (M.fromList [(7, 29), (1, 3), (2, 4)]))
+      actualResult1 `shouldBe` actualResult2
+
+    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, LowerBound 10)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      actualResult `shouldBe` Nothing
+
+  describe "twoPhaseSimplex' with negative lower bound s (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, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just 5
+
+      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, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-3)
+
+      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, LowerBound (-10))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-2)
+
+      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, LowerBound (-10))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-10)
+
+    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, LowerBound (-2)), (2, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+                objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            -- Verify the actual objective value
+            objVal `shouldBe` 10
+            -- Verify lower bounds are respected
+            x1 `shouldSatisfy` (>= (-2))
+            x2 `shouldSatisfy` (>= (-3))
+
+      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, LowerBound (-2)), (2, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            -- Verify the actual objective value
+            objVal `shouldBe` (-5)
+            M.lookup 1 result.varValMap `shouldBe` Just (-2)
+            M.lookup 2 result.varValMap `shouldBe` Just (-3)
+
+      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, LowerBound (-5)), (2, LowerBound (-4))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+            M.lookup 1 result.varValMap `shouldBe` Just 3
+            M.lookup 2 result.varValMap `shouldBe` Just (-4)
+            -- Verify objective value computed from variables
+            (2 * x1 - x2) `shouldBe` 10
+
+      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, LowerBound (-5)), (2, LowerBound (-4))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+            M.lookup 1 result.varValMap `shouldBe` Just (-5)
+            M.lookup 2 result.varValMap `shouldBe` Just 6
+            -- Verify objective value computed from variables
+            (2 * x1 - x2) `shouldBe` (-16)
+
+    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, LowerBound (-5))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just 10
+
+      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, LowerBound (-5))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just 2
+
+    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, LowerBound (-5)), (2, LowerBound (-5))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            M.lookup 1 result.varValMap `shouldBe` Just 10
+            M.lookup 2 result.varValMap `shouldBe` Just 10
+            -- Verify objective value
+            objVal `shouldBe` 20
+
+      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, LowerBound (-5)), (2, LowerBound (-5))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            M.lookup 1 result.varValMap `shouldBe` Just (-5)
+            M.lookup 2 result.varValMap `shouldBe` Just (-5)
+            -- Verify objective value
+            objVal `shouldBe` (-10)
+
+    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, LowerBound ((-7) % 2))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just (5 % 2)
+
+      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, LowerBound ((-7) % 2))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just ((-7) % 2)
+
+  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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just 10
+
+      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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> M.lookup 1 result.varValMap `shouldBe` Just (-10)
+
+      it "Unbounded variable with only upper bound: Min finds negative value" $ do
+        -- x₁ unbounded, only x₁ ≤ 5, minimize x₁
+        -- This should be unbounded (no 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 (infeasible for optimization)
+        actualResult `shouldBe` Nothing
+
+    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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            M.lookup 1 result.varValMap `shouldBe` Just 5
+            M.lookup 2 result.varValMap `shouldBe` Just 7
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            objVal `shouldBe` 12
+
+      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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            M.lookup 1 result.varValMap `shouldBe` Just (-5)
+            M.lookup 2 result.varValMap `shouldBe` Just (-3)
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            objVal `shouldBe` (-8)
+
+      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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            M.lookup 1 result.varValMap `shouldBe` Just 5
+            M.lookup 2 result.varValMap `shouldBe` Just (-3)
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            objVal `shouldBe` 8
+
+    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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            M.lookup 1 result.varValMap `shouldBe` Just 15
+            M.lookup 2 result.varValMap `shouldBe` Just (-5)
+
+      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
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            M.lookup 1 result.varValMap `shouldBe` Just (-10)
+            M.lookup 2 result.varValMap `shouldBe` Just 20
+
+  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, T.NonNegative), (2, LowerBound (-5)), (3, Unbounded)]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            -- Verify objective value
+            objVal `shouldBe` 20
+
+      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, T.NonNegative), (2, LowerBound (-5)), (3, Unbounded)]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+                x3 = M.findWithDefault 0 3 result.varValMap
+                objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+            -- Verify constraints
+            x1 `shouldSatisfy` (>= 0)
+            x2 `shouldSatisfy` (>= (-5))
+            x3 `shouldSatisfy` (>= (-20))
+            -- Verify objective value
+            objVal `shouldBe` (-10)
+
+    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, LowerBound 2), (2, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+            -- Verify constraints
+            x1 `shouldSatisfy` (>= 2)
+            x2 `shouldSatisfy` (>= (-3))
+            (2 * x1 + x2) `shouldSatisfy` (<= 20)
+
+      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, LowerBound 2), (2, LowerBound (-3))]
+        actualResult <-
+          runStdoutLoggingT $
+            filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+              twoPhaseSimplex' domainMap obj constraints
+        case actualResult of
+          Nothing -> expectationFailure "Expected a solution but got Nothing"
+          Just result -> do
+            let x1 = M.findWithDefault 0 1 result.varValMap
+                x2 = M.findWithDefault 0 2 result.varValMap
+            x1 `shouldSatisfy` (>= 2)
+            x2 `shouldSatisfy` (>= (-3))
+            (x1 + x2) `shouldSatisfy` (>= 0)
+
+  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, LowerBound (-5))]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      actualResult `shouldBe` Nothing
+
+    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
+      actualResult `shouldBe` Nothing
+
+    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, LowerBound (-5))]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 0
+
+    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
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> M.lookup 1 result.varValMap `shouldBe` Just 0
+
+    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, T.NonNegative), (2, LowerBound (-10)), (3, T.NonNegative)]
+      actualResult <-
+        runStdoutLoggingT $
+          filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $
+            twoPhaseSimplex' domainMap obj constraints
+      case actualResult of
+        Nothing -> expectationFailure "Expected a solution but got Nothing"
+        Just result -> do
+          let objVal = M.findWithDefault 0 result.objectiveVar result.varValMap
+          -- Verify objective value
+          objVal `shouldBe` 15
+
+  -- ===========================================================================
+  -- 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 Nothing for NonNegative domain" $ do
+        getTransform 10 1 T.NonNegative `shouldBe` Nothing
+
+      it "returns Nothing for LowerBound 0" $ do
+        getTransform 10 1 (LowerBound 0) `shouldBe` Nothing
+
+      it "returns AddLowerBound for positive lower bound" $ do
+        getTransform 10 1 (LowerBound 5) `shouldBe` Just (AddLowerBound 1 5)
+
+      it "returns AddLowerBound for fractional positive lower bound" $ do
+        getTransform 10 1 (LowerBound (3 % 2)) `shouldBe` Just (AddLowerBound 1 (3 % 2))
+
+      it "returns Shift for negative lower bound" $ do
+        getTransform 10 1 (LowerBound (-5)) `shouldBe` Just (Shift 1 10 (-5))
+
+      it "returns Shift for fractional negative lower bound" $ do
+        getTransform 10 1 (LowerBound ((-7) % 3)) `shouldBe` Just (Shift 1 10 ((-7) % 3))
+
+      it "returns Split for Unbounded domain" $ do
+        getTransform 10 1 Unbounded `shouldBe` Just (Split 1 10 11)
+
+  describe "generateTransform" $ do
+    describe "Unit tests" $ do
+      it "generates no transform for NonNegative in domain map" $ do
+        let domainMap = M.fromList [(1, T.NonNegative)]
+        generateTransform domainMap 1 ([], 10) `shouldBe` ([], 10)
+
+      it "generates AddLowerBound for positive bound in domain map" $ do
+        let domainMap = M.fromList [(1, LowerBound 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, LowerBound (-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, LowerBound 5)]
+            existing = [AddLowerBound 2 3]
+        generateTransform domainMap 1 (existing, 10) `shouldBe` ([AddLowerBound 1 5, AddLowerBound 2 3], 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 "unapplyTransform and unapplyTransforms" $ do
+    describe "Unit tests" $ do
+      it "unapplyTransform AddLowerBound leaves result unchanged" $ do
+        let result = Result 5 (M.fromList [(5, 10), (1, 7)])
+            transform = AddLowerBound 1 5
+        unapplyTransform transform result `shouldBe` result
+
+      it "unapplyTransform Shift recovers original variable value" $ do
+        -- originalVar = shiftedVar + shiftBy
+        -- If shiftedVar = 15 and shiftBy = -5, then originalVar = 10
+        let result = Result 5 (M.fromList [(5, 100), (10, 15)])
+            transform = Shift 1 10 (-5)
+        let newResult = unapplyTransform transform result
+        M.lookup 1 (varValMap newResult) `shouldBe` Just 10
+        M.lookup 10 (varValMap newResult) `shouldBe` Nothing
+
+      it "unapplyTransform Split recovers original variable value" $ do
+        -- originalVar = posVar - negVar
+        -- If posVar = 8 and negVar = 3, then originalVar = 5
+        let result = Result 5 (M.fromList [(5, 100), (10, 8), (11, 3)])
+            transform = Split 1 10 11
+        let newResult = unapplyTransform transform result
+        M.lookup 1 (varValMap newResult) `shouldBe` Just 5
+        M.lookup 10 (varValMap newResult) `shouldBe` Nothing
+        M.lookup 11 (varValMap newResult) `shouldBe` Nothing
+
+      it "unapplyTransform Split handles negative original value" $ do
+        -- originalVar = posVar - negVar
+        -- If posVar = 2 and negVar = 7, then originalVar = -5
+        let result = Result 5 (M.fromList [(5, 100), (10, 2), (11, 7)])
+            transform = Split 1 10 11
+        let newResult = unapplyTransform transform result
+        M.lookup 1 (varValMap newResult) `shouldBe` Just (-5)
+
+      it "unapplyTransforms 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 result = Result 5 (M.fromList [(5, 100), (10, 15), (11, 8)])
+            transforms = [Shift 1 10 (-5), Shift 2 11 (-3)]
+        let newResult = unapplyTransforms transforms result
+        M.lookup 1 (varValMap newResult) `shouldBe` Just 10
+        M.lookup 2 (varValMap newResult) `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, T.NonNegative), (2, T.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, LowerBound 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, LowerBound (-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, T.NonNegative), (2, LowerBound 5), (3, LowerBound (-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 Nothing" $ property $
+        \(nextVar :: Int) (v :: Int) ->
+          getTransform (abs nextVar + 1) (abs v + 1) T.NonNegative == Nothing
+
+      it "LowerBound 0 produces Nothing" $ property $
+        \(nextVar :: Int) (v :: Int) ->
+          getTransform (abs nextVar + 1) (abs v + 1) (LowerBound 0) == Nothing
+
+      it "positive LowerBound produces AddLowerBound" $ property $
+        \(Positive bound :: Positive Rational) (nextVar :: Int) (v :: Int) ->
+          case getTransform (abs nextVar + 1) (abs v + 1) (LowerBound bound) of
+            Just (AddLowerBound var b) -> var == abs v + 1 && b == bound
+            _ -> False
+
+      it "negative LowerBound produces Shift" $ property $
+        \(Positive bound :: Positive Rational) (nextVar :: Int) (v :: Int) ->
+          let negBound = negate bound
+          in case getTransform (abs nextVar + 1) (abs v + 1) (LowerBound negBound) of
+            Just (Shift origVar _ shiftBy) -> 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
+            Just (Split origVar _ _) -> origVar == abs v + 1
+            _ -> False
+
+    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 "unapplyTransform Shift properties" $ do
+      it "recovers originalVar = shiftedVar + shiftBy" $ property $
+        \(shiftedVal :: Rational) (shiftBy :: Rational) ->
+          let result = Result 5 (M.fromList [(5, 100), (10, shiftedVal)])
+              transform = Shift 1 10 shiftBy
+              newResult = unapplyTransform transform result
+          in M.lookup 1 (varValMap newResult) == Just (shiftedVal + shiftBy)
+
+      it "removes shifted variable from result" $ property $
+        \(shiftedVal :: Rational) (shiftBy :: Rational) ->
+          let result = Result 5 (M.fromList [(5, 100), (10, shiftedVal)])
+              transform = Shift 1 10 shiftBy
+              newResult = unapplyTransform transform result
+          in M.lookup 10 (varValMap newResult) == Nothing
+
+    describe "unapplyTransform Split properties" $ do
+      it "recovers originalVar = posVar - negVar" $ property $
+        \(posVal :: Rational) (negVal :: Rational) ->
+          let result = Result 5 (M.fromList [(5, 100), (10, posVal), (11, negVal)])
+              transform = Split 1 10 11
+              newResult = unapplyTransform transform result
+          in M.lookup 1 (varValMap newResult) == Just (posVal - negVal)
+
+      it "removes pos and neg variables from result" $ property $
+        \(posVal :: Rational) (negVal :: Rational) ->
+          let result = Result 5 (M.fromList [(5, 100), (10, posVal), (11, negVal)])
+              transform = Split 1 10 11
+              newResult = unapplyTransform transform result
+          in M.lookup 10 (varValMap newResult) == Nothing && 
+             M.lookup 11 (varValMap newResult) == Nothing
+
+    describe "Round-trip properties" $ do
+      it "Shift transform and unapply is identity for variable value" $ property $
+        \(origVal :: Rational) (shiftBy :: Rational) ->
+          shiftBy < 0 ==>  -- Only negative shifts are valid
+            let shiftedVal = origVal - shiftBy  -- shiftedVar = originalVar - shiftBy
+                result = Result 5 (M.fromList [(5, 100), (10, shiftedVal)])
+                transform = Shift 1 10 shiftBy
+                newResult = unapplyTransform transform result
+            in M.lookup 1 (varValMap newResult) == Just origVal
+
+      it "Split with posVal=origVal and negVal=0 gives correct value for positive origVal" $ property $
+        \(Positive origVal :: Positive Rational) ->
+          let result = Result 5 (M.fromList [(5, 100), (10, origVal), (11, 0)])
+              transform = Split 1 10 11
+              newResult = unapplyTransform transform result
+          in M.lookup 1 (varValMap newResult) == 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
+              result = Result 5 (M.fromList [(5, 100), (10, 0), (11, origVal)])
+              transform = Split 1 10 11
+              newResult = unapplyTransform transform result
+          in M.lookup 1 (varValMap newResult) == Just negOrigVal
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)
-    ]
-  )