add L1 and LInf solvers

4 files changed, 70 insertions, 2 deletions

diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index 5f3b34e..533ee2d 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -24,11 +24,12 @@ extra-source-files:     examples/simplex1.hs
 library
     Build-Depends:      base, hmatrix
 
-    hs-source-dirs:     lib
+    hs-source-dirs:     src
 
     Exposed-modules:    Numeric.LinearProgramming
+                        Numeric.LinearProgramming.L1
 
-    c-sources:          lib/Numeric/LinearProgramming/glpk.c
+    c-sources:          src/C/glpk.c
 
     ghc-options:  -Wall
 
diff --git a/packages/glpk/lib/Numeric/LinearProgramming/glpk.c b/packages/glpk/src/C/glpk.c
index bfbb435..bfbb435 100644
--- a/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
+++ b/packages/glpk/src/C/glpk.c
diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index 25cdab2..25cdab2 100644
--- a/packages/glpk/lib/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
diff --git a/packages/glpk/src/Numeric/LinearProgramming/L1.hs b/packages/glpk/src/Numeric/LinearProgramming/L1.hs
new file mode 100644
index 0000000..2e8f94a
--- /dev/null
+++ b/packages/glpk/src/Numeric/LinearProgramming/L1.hs
@@ -0,0 +1,67 @@
+{- |
+Module      :  Numeric.LinearProgramming.L1
+Copyright   :  (c) Alberto Ruiz 2011-14
+Stability   :  provisional
+
+Linear system solvers in the L_1 norm using linear programming.
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.LinearProgramming.L1 (
+    l1SolveO, lInfSolveO,
+    l1SolveU,
+) where
+
+import Numeric.LinearAlgebra
+import Numeric.LinearProgramming
+
+-- | L_Inf solution of overconstrained system Ax=b.
+--
+-- Find argmin x ||Ax-b||_inf
+lInfSolveO :: Matrix Double -> Vector Double -> Vector Double
+lInfSolveO a b = fromList (take n x)
+  where
+    n = cols a
+    as = toRows a
+    bs = toList b
+    c1 = zipWith (mk (1)) as bs
+    c2 = zipWith (mk (-1)) as bs
+    mk sign a_i b_i = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(n+1)]) :<=: (sign * b_i)
+    p = Sparse (c1++c2)
+    Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ [1])) p (map Free [1..(n+1)])
+
+
+-- | L_1 solution of overconstrained system Ax=b.
+--
+-- Find argmin x ||Ax-b||_1.
+l1SolveO :: Matrix Double -> Vector Double -> Vector Double
+l1SolveO a b = fromList (take n x)
+  where
+    n = cols a
+    m = rows a
+    as = toRows a
+    bs = toList b
+    ks = [1..]
+    c1 = zipWith3 (mk (1)) as bs ks
+    c2 = zipWith3 (mk (-1)) as bs ks
+    mk sign a_i b_i k = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(k+n)]) :<=: (sign * b_i)
+    p = Sparse (c1++c2)
+    Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate m 1)) p (map Free [1..(n+m)])
+
+
+
+-- | L1 solution of underconstrained linear system Ax=b.
+--
+-- Find argmin x ||x||_1 such that Ax=b.
+l1SolveU :: Matrix Double -> Vector Double -> Vector Double
+l1SolveU a y = fromList (take n x)
+  where
+    n = cols a
+    c1 = map (\k ->  [ 1#k, -1#k+n] :<=: 0) [1..n]
+    c2 = map (\k ->  [-1#k, -1#k+n] :<=: 0) [1..n]
+    c3 = zipWith (:==:) (map sp $ toRows a) (toList y)
+    sp v = zipWith (#) (toList v) [1..]
+    p = Sparse (c1 ++ c2 ++ c3)
+    Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate n 1)) p (map Free [1..(2*n)])
+