hmatrix-glpk 0.3.1 → 0.4.0

raw patch · 6 files changed

+473/−346 lines, 6 filesdep ~basedep ~hmatrixPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base, hmatrix

API changes (from Hackage documentation)

Files

• hmatrix-glpk.cabal +6/−5
• lib/Numeric/LinearProgramming.hs +0/−265
• lib/Numeric/LinearProgramming/glpk.c +0/−76
• src/C/glpk.c +76/−0
• src/Numeric/LinearProgramming.hs +271/−0
• src/Numeric/LinearProgramming/L1.hs +120/−0

@@ -1,5 +1,5 @@ Name:               hmatrix-glpk-Version:            0.3.1+Version:            0.4.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -11,7 +11,7 @@  Simple interface to linear programming functions provided by GLPK.  Category:           Math-tested-with:        GHC ==7.4+tested-with:        GHC ==7.8  cabal-version:      >=1.6 build-type:         Simple@@ -22,13 +22,14 @@                         examples/simplex4.hs  library-    Build-Depends:      base >= 3 && < 5, hmatrix >= 0.8.3+    Build-Depends:      base <5, hmatrix >= 0.16 -    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 

@@ -1,265 +0,0 @@-{-# LANGUAGE ForeignFunctionInterface #-}--{- |-Module      :  Numeric.LinearProgramming-Copyright   :  (c) Alberto Ruiz 2010-License     :  GPL--Maintainer  :  Alberto Ruiz (aruiz at um dot es)-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 Numeric.LinearAlgebra hiding (i)-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)-    let fi = fromIntegral-    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 ]---}-

@@ -1,76 +0,0 @@-#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;-}

@@ -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;+}

@@ -0,0 +1,271 @@+{-# 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 ]++-}+

@@ -0,0 +1,120 @@+{- |+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 (+    l1Solve, l1SolveGT,+    l1SolveO, lInfSolveO,+    l1SolveU,+) where++import Numeric.LinearAlgebra+import Numeric.LinearProgramming++-- | L_inf solution of overconstrained system Ax=b.+--+-- @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.+--+-- @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.+--+-- @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)])++--------------------------------------------------------------------------------+-- | Solution in the L_1 norm, with L_1 regularization, of a linear system @Ax=b@.+--+-- @argmin_x  λ||x||_1 + ||Ax-b||_1@+l1Solve+    :: Double        -- ^ λ+    -> Matrix Double -- ^ A+    -> Vector Double -- ^ b+    -> Vector Double -- ^ x+l1Solve λ a b = fromList (take n x)+  where+    n = cols a+    m = rows a+    as = toRows a+    bs = toList b+    c1Res = zipWith3 (mkR (1)) as bs [1..m]+    c2Res = zipWith3 (mkR (-1)) as bs [1..m]+    mkR sign a_i b_i k = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(k+2*n)]) :<=: (sign * b_i)+    c1Sol = map (\k ->  [ 1#k, -1#k+n] :<=: 0) [1..n]+    c2Sol = map (\k ->  [-1#k, -1#k+n] :<=: 0) [1..n]+    p = Sparse (c1Res++c2Res++c1Sol++c2Sol)+    cost = replicate n 0 ++ replicate n λ ++ replicate m 1+    Optimal (_j,x) = simplex (Minimize cost) p (map Free [1..(2*n+m)])++--------------------------------------------------------------------------------++-- | Solution in the L_1 norm, with L_1 regularization, of a system of linear inequalities @Ax>=b@.+--+-- @argmin_x  λ||x||_1 + ||step(b-Ax)||_1@+l1SolveGT+    :: Double        -- ^ λ+    -> Matrix Double -- ^ A+    -> Vector Double -- ^ b+    -> Vector Double -- ^ x+l1SolveGT λ a b = fromList (take n x)+  where+    n = cols a+    m = rows a+    as = toRows a+    bs = toList b+    cRes = zipWith3 mkR as bs [1..m]+    mkR a_i b_i k = (zipWith (#) (toList a_i) [1..] ++ [1#(k+2*n)]) :>=: (b_i)+    c1Sol = map (\k ->  [ 1#k, -1#k+n] :<=: 0) [1..n]+    c2Sol = map (\k ->  [-1#k, -1#k+n] :<=: 0) [1..n]+    p = Sparse (cRes++c1Sol++c2Sol)+    cost = replicate n 0 ++ replicate n λ ++ replicate m 1+    Optimal (_j,x) = simplex (Minimize cost) p (map Free [1..(2*n)])++--------------------------------------------------------------------------------++