add L1 and LInf solvers

3 files changed, 407 insertions, 0 deletions

diff --git a/packages/glpk/src/C/glpk.c b/packages/glpk/src/C/glpk.c
new file mode 100644
index 0000000..bfbb435
--- /dev/null
+++ b/packages/glpk/src/C/glpk.c
@@ -0,0 +1,76 @@
+#define DVEC(A) int A##n, double*A##p
+#define DMAT(A) int A##r, int A##c, double*A##p
+
+#define AT(M,r,co) (M##p[(r)*M##c+(co)])
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <glpk.h>
+#include <math.h>
+
+/*-----------------------------------------------------*/
+
+int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
+    glp_prob *lp;
+    lp = glp_create_prob();
+    glp_set_obj_dir(lp, GLP_MAX);
+    int i,j,k;
+    int tot = cr - n;
+    glp_add_rows(lp, m);
+    glp_add_cols(lp, n);
+
+    //printf("%d %d\n",m,n);
+
+    // the first n values
+    for (k=1;k<=n;k++) {
+        glp_set_obj_coef(lp, k, AT(c, k-1, 2));
+        //printf("%d %f\n",k,AT(c, k-1, 2));
+    }
+
+    int * ia = malloc((1+tot)*sizeof(int));
+    int * ja = malloc((1+tot)*sizeof(int));
+    double * ar = malloc((1+tot)*sizeof(double));
+
+    for (k=1; k<= tot; k++) {
+        ia[k] = rint(AT(c,k-1+n,0));
+        ja[k] = rint(AT(c,k-1+n,1));
+        ar[k] =      AT(c,k-1+n,2);
+        //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
+    }
+    glp_load_matrix(lp, tot, ia, ja, ar);
+
+    int t;
+    for (i=1;i<=m;i++) {
+    switch((int)rint(AT(b,i-1,0))) {
+        case 0: { t = GLP_FR; break; }
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }
+    }
+    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
+    }
+    for (j=1;j<=n;j++) {
+    switch((int)rint(AT(b,m+j-1,0))) {
+        case 0: { t = GLP_FR; break; }
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }
+    }
+    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
+    }
+    glp_term_out(0);
+    glp_simplex(lp, NULL);
+    sp[0] = glp_get_status(lp);
+    sp[1] = glp_get_obj_val(lp);
+    for (k=1; k<=n; k++) {
+        sp[k+1] = glp_get_col_prim(lp, k);
+    }
+    glp_delete_prob(lp);
+    free(ia);
+    free(ja);
+    free(ar);
+
+    return 0;
+}
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
new file mode 100644
index 0000000..25cdab2
--- /dev/null
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -0,0 +1,264 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+{- |
+Module      :  Numeric.LinearProgramming
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+This module provides an interface to the standard simplex algorithm.
+
+For example, the following LP problem
+
+@maximize 4 x_1 - 3 x_2 + 2 x_3
+subject to
+
+2 x_1 +   x_2 <= 10
+  x_3 + 5 x_4 <= 20
+
+and
+
+x_i >= 0@
+
+can be solved as follows:
+
+@import Numeric.LinearProgramming
+
+prob = Maximize [4, -3, 2]
+
+constr1 = Sparse [ [2\#1, 1\#2] :<=: 10
+                 , [1\#2, 5\#3] :<=: 20
+                 ]
+
+\> simplex prob constr1 []
+Optimal (28.0,[5.0,0.0,4.0])@
+
+The coefficients of the constraint matrix can also be given in dense format:
+
+@constr2 = Dense [ [2,1,0] :<=: 10
+                , [0,1,5] :<=: 20
+                ]@
+
+By default all variables are bounded as @x_i >= 0@, but this can be
+changed:
+
+@\> simplex prob constr2 [ 2 :=>: 1, 3 :&: (2,7)]
+Optimal (22.6,[4.5,1.0,3.8])
+
+\> simplex prob constr2 [Free 2]
+Unbounded@
+
+The given bound for a variable completely replaces the default,
+so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.
+Multiple bounds for a variable are not allowed, instead of
+@[i :=>: a, i:<=: b]@ use @i :&: (a,b)@.
+
+-}
+
+module Numeric.LinearProgramming(
+    simplex,
+    Optimization(..),
+    Constraints(..),
+    Bounds,
+    Bound(..),
+    (#),
+    Solution(..)
+) where
+
+import Data.Packed
+import Data.Packed.Development
+import Foreign(Ptr)
+import System.IO.Unsafe(unsafePerformIO)
+import Foreign.C.Types
+import Data.List((\\),sortBy,nub)
+import Data.Function(on)
+
+--import Debug.Trace
+--debug x = trace (show x) x
+
+-----------------------------------------------------
+
+-- | Coefficient of a variable for a sparse representation of constraints.
+(#) :: Double -> Int -> (Double,Int)
+infixl 5 #
+(#) = (,)
+
+data Bound x =  x :<=: Double
+             |  x :=>: Double
+             |  x :&: (Double,Double)
+             |  x :==: Double
+             |  Free x
+             deriving Show
+
+data Solution = Undefined
+              | Feasible (Double, [Double])
+              | Infeasible (Double, [Double])
+              | NoFeasible
+              | Optimal (Double, [Double])
+              | Unbounded
+              deriving Show
+
+data Constraints = Dense  [ Bound [Double] ]
+                 | Sparse [ Bound [(Double,Int)] ]
+
+data Optimization = Maximize [Double]
+                  | Minimize [Double]
+
+type Bounds = [Bound Int]
+
+simplex :: Optimization -> Constraints -> Bounds -> Solution
+
+simplex opt (Dense  []) bnds = simplex opt (Sparse []) bnds
+simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+
+simplex opt (Dense constr) bnds = extract sg sol where
+    sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
+    n = length objfun
+    m = length constr
+    (sz, sg, objfun) = adapt opt
+
+simplex opt (Sparse constr) bnds = extract sg sol where
+    sol = simplexSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)
+    n = length objfun
+    m = length constr
+    (sz, sg, objfun) = adapt opt
+
+adapt :: Optimization -> (Int, Double, [Double])
+adapt opt = case opt of
+    Maximize x -> (size x, 1 ,x)
+    Minimize x -> (size x, -1, (map negate x))
+ where size x | null x = error "simplex: objective function with zero variables"
+              | otherwise = length x
+
+extract :: Double -> Vector Double -> Solution
+extract sg sol = r where
+    z = sg * (sol@>1)
+    v = toList $ subVector 2 (dim sol -2) sol
+    r = case round(sol@>0)::Int of
+          1 -> Undefined
+          2 -> Feasible (z,v)
+          3 -> Infeasible (z,v)
+          4 -> NoFeasible
+          5 -> Optimal (z,v)
+          6 -> Unbounded
+          _ -> error "simplex: solution type unknown"
+
+-----------------------------------------------------
+
+obj :: Bound t -> t
+obj (x :<=: _)  = x
+obj (x :=>: _)  = x
+obj (x :&: _)  = x
+obj (x :==: _) = x
+obj (Free x)   = x
+
+tb :: Bound t -> Double
+tb (_ :<=: _)  = glpUP
+tb (_ :=>: _)  = glpLO
+tb (_ :&: _)  = glpDB
+tb (_ :==: _) = glpFX
+tb (Free _)   = glpFR
+
+lb :: Bound t -> Double
+lb (_ :<=: _)     = 0
+lb (_ :=>: a)     = a
+lb (_ :&: (a,_)) = a
+lb (_ :==: a)    = a
+lb (Free _)      = 0
+
+ub :: Bound t -> Double
+ub (_ :<=: a)     = a
+ub (_ :=>: _)     = 0
+ub (_ :&: (_,a)) = a
+ub (_ :==: a)    = a
+ub (Free _)      = 0
+
+mkBound1 :: Bound t -> [Double]
+mkBound1 b = [tb b, lb b, ub b]
+
+mkBound2 :: Bound t -> (t, [Double])
+mkBound2 b = (obj b, mkBound1 b)
+
+mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double
+mkBounds n b1 b2 = fromLists (cb++vb) where
+    gv' = map obj b2
+    gv | nub gv' == gv' = gv'
+       | otherwise = error $ "simplex: duplicate bounds for vars " ++ show (gv'\\nub gv')
+    rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv
+       | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n
+    vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:=>: 0)) rv ++ map mkBound2 b2
+    cb = map mkBound1 b1
+
+mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double
+mkConstrD n f b1 | ok = fromLists (ob ++ co)
+                 | otherwise = error $ "simplex: dense constraints require "++show n
+                                     ++" variables, given " ++ show ls
+    where
+       cs = map obj b1
+       ls = map length cs
+       ok = all (==n) ls
+       den = fromLists cs
+       ob = map (([0,0]++).return) f
+       co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
+
+mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
+mkConstrS n objfun b1 = fromLists (ob ++ co) where
+    ob = map (([0,0]++).return) objfun
+    co = concat $ zipWith f [1::Int ..] cs
+    cs = map obj b1
+    f k = map (g k)
+    g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]
+              | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n
+
+-----------------------------------------------------
+
+foreign import ccall unsafe "c_simplex_sparse" c_simplex_sparse
+    :: CInt -> CInt                  -- rows and cols
+    -> CInt -> CInt -> Ptr Double    -- coeffs
+    -> CInt -> CInt -> Ptr Double    -- bounds
+    -> CInt -> Ptr Double            -- result
+    -> IO CInt                       -- exit code
+
+simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
+simplexSparse m n c b = unsafePerformIO $ do
+    s <- createVector (2+n)
+    app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
+    return s
+
+glpFR, glpLO, glpUP, glpDB, glpFX :: Double
+glpFR = 0
+glpLO = 1
+glpUP = 2
+glpDB = 3
+glpFX = 4
+
+{- Raw format of coeffs
+
+simplexSparse
+
+(12><3)
+ [ 0.0, 0.0, 10.0
+ , 0.0, 0.0,  6.0
+ , 0.0, 0.0,  4.0
+ , 1.0, 1.0,  1.0
+ , 1.0, 2.0,  1.0
+ , 1.0, 3.0,  1.0
+ , 2.0, 1.0, 10.0
+ , 2.0, 2.0,  4.0
+ , 2.0, 3.0,  5.0
+ , 3.0, 1.0,  2.0
+ , 3.0, 2.0,  2.0
+ , 3.0, 3.0,  6.0 ]
+
+bounds = (6><3)
+  [ glpUP,0,100
+  , glpUP,0,600
+  , glpUP,0,300
+  , glpLO,0,0
+  , glpLO,0,0
+  , glpLO,0,0 ]
+
+-}
+
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)])
+