diff --git a/hmatrix-glpk.cabal b/hmatrix-glpk.cabal --- a/hmatrix-glpk.cabal +++ b/hmatrix-glpk.cabal @@ -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 diff --git a/lib/Numeric/LinearProgramming.hs b/lib/Numeric/LinearProgramming.hs deleted file mode 100644 --- a/lib/Numeric/LinearProgramming.hs +++ /dev/null @@ -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 ] - --} - diff --git a/lib/Numeric/LinearProgramming/glpk.c b/lib/Numeric/LinearProgramming/glpk.c deleted file mode 100644 --- a/lib/Numeric/LinearProgramming/glpk.c +++ /dev/null @@ -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 -#include -#include -#include - -/*-----------------------------------------------------*/ - -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/src/C/glpk.c b/src/C/glpk.c new file mode 100644 --- /dev/null +++ b/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 +#include +#include +#include + +/*-----------------------------------------------------*/ + +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/src/Numeric/LinearProgramming.hs b/src/Numeric/LinearProgramming.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearProgramming.hs @@ -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 ] + +-} + diff --git a/src/Numeric/LinearProgramming/L1.hs b/src/Numeric/LinearProgramming/L1.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearProgramming/L1.hs @@ -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)]) + +-------------------------------------------------------------------------------- + +