Merge branch 'master' of github.com:albertoruiz/hmatrix

8 files changed, 178 insertions, 72 deletions

diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index dda47cb..599f69e 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -13,6 +13,7 @@ typedef float  complex TCF;
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <stdint.h>
 
 #define MACRO(B) do {B} while (0)
 #define ERROR(CODE) MACRO(return CODE;)
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 9adc023..257ad73 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -241,7 +241,7 @@ instance Container Vector (Complex Float)
 
 ---------------------------------------------------------------
 
-instance (Container Vector a) => Container Matrix a
+instance (Fractional a, Element a, Container Vector a) => Container Matrix a
   where
     size' = size
     scale' x = liftMatrix (scale' x)
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 037396d..3398e6a 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -52,7 +52,7 @@ module Numeric.LinearAlgebra.Static(
     linSolve, (<\>),
     -- * Factorizations
     svd, withCompactSVD, svdTall, svdFlat, Eigen(..),
-    withNullspace, qr,
+    withNullspace, qr, chol,
     -- * Misc
     mean,
     Disp(..), Domain(..),
@@ -68,7 +68,7 @@ import Numeric.LinearAlgebra.HMatrix hiding (
     row,col,vector,matrix,linspace,toRows,toColumns,
     (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
     eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
-    qr,size,app,mul,dot)
+    qr,size,app,mul,dot,chol)
 import qualified Numeric.LinearAlgebra.HMatrix as LA
 import Data.Proxy(Proxy)
 import Numeric.LinearAlgebra.Static.Internal
@@ -306,6 +306,9 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n)
       where
         (l,v) = LA.eig m
 
+chol :: KnownNat n => Sym n -> Sq n
+chol (extract . unSym -> m) = mkL $ LA.cholSH m
+
 --------------------------------------------------------------------------------
 
 withNullspace
diff --git a/packages/glpk/examples/simplex2.hs b/packages/glpk/examples/simplex2.hs
index e9e8859..0d83ca9 100644
--- a/packages/glpk/examples/simplex2.hs
+++ b/packages/glpk/examples/simplex2.hs
@@ -10,9 +10,14 @@ constr2 = Dense [ [2,1,0] :<=: 10
                 , [0,1,5] :<=: 20
                 ]
 
+constr3 = General [ [1#1, 1#1, 1#2] :<=: 10
+                  , [1#2, 5#3] :<=: 20
+                  ]
+
 main = do
     print $ simplex prob constr1 []
     print $ simplex prob constr2 []
+    print $ simplex prob constr3 []
     print $ simplex prob constr2 [ 2 :>=: 1, 3 :&: (2,7)]
     print $ simplex prob constr2 [ Free 2 ]
 
diff --git a/packages/glpk/examples/simplex5.hs b/packages/glpk/examples/simplex5.hs
new file mode 100644
index 0000000..ecbcdaa
--- /dev/null
+++ b/packages/glpk/examples/simplex5.hs
@@ -0,0 +1,27 @@
+import Numeric.LinearProgramming
+
+-- This is a linear program from the paper "Picking vs. Guessing Secrets: A Game-theoretic Analysis"
+
+gamma = 100000 :: Double
+sigma = 1 :: Double
+n = 64 :: Int
+cost_fun :: Int -> Double
+cost_fun i = (fromIntegral i) / (fromIntegral n)
+size_fun :: Int -> Double
+size_fun i = 2^(fromIntegral i)
+
+prob = Minimize $ map cost_fun [1..n]
+bnds = [i :&: (0,1) | i <- [1..n]]
+
+constr1 = [[1 # i | i <- [1..n]] :==: 1] ++
+          [[1/(size_fun i) # i,
+            -1/(size_fun (i+1)) # i+1] :>=: 0 | i <- [1..n-1]] ++
+          [(
+            [gamma#i | i <- [1..k]] ++
+            (concat [[sigma*(size_fun i) # j | j <- [1..i-1]] | i <- [1..k]]) ++
+            [((size_fun i) - 1)/2 # i | i <- [1..k]])
+           :<=: (sigma * (foldr (+) 0 (map size_fun [1..k]))) | k <- [1..n]]
+
+main = do
+  print $ simplex prob (General constr1) bnds -- NoFeasible
+  print $ exact   prob (General constr1) bnds -- solution found
diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index cd761e0..a9859f9 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -20,9 +20,10 @@ extra-source-files:     examples/simplex1.hs
                         examples/simplex2.hs
                         examples/simplex3.hs
                         examples/simplex4.hs
+                        examples/simplex5.hs
 
 library
-    Build-Depends:      base <5, hmatrix >= 0.16
+    Build-Depends:      base <5, hmatrix >= 0.16, containers >= 0.5.4.0
 
     hs-source-dirs:     src
 
diff --git a/packages/glpk/src/C/glpk.c b/packages/glpk/src/C/glpk.c
index bfbb435..86b1277 100644
--- a/packages/glpk/src/C/glpk.c
+++ b/packages/glpk/src/C/glpk.c
@@ -10,67 +10,71 @@
 
 /*-----------------------------------------------------*/
 
-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);
+#define C_X_SPARSE(X)                                           \
+  int c_##X##_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_##X(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;                                                         \
+  }                                                                     \
 
-    //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;
-}
+C_X_SPARSE(simplex);
+C_X_SPARSE(exact);
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index b0537cc..6a0c47d 100644
--- a/packages/glpk/src/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -49,6 +49,14 @@ constr2 = Dense [ [2,1,0] :<=: 10
                 ]
 @
 
+Note that when using sparse constraints, coefficients cannot appear more than once in each constraint. You can alternatively use General which will automatically sum any duplicate coefficients when necessary.
+
+@
+constr3 = General [ [1\#1, 1\#1, 1\#2] :<=: 10
+                  , [1\#2, 5\#3] :<=: 20
+                  ]
+@
+
 By default all variables are bounded as @x_i >= 0@, but this can be
 changed:
 
@@ -67,6 +75,8 @@ Multiple bounds for a variable are not allowed, instead of
 
 module Numeric.LinearProgramming(
     simplex,
+    exact,
+    sparseOfGeneral,
     Optimization(..),
     Constraints(..),
     Bounds,
@@ -82,13 +92,14 @@ import System.IO.Unsafe(unsafePerformIO)
 import Foreign.C.Types
 import Data.List((\\),sortBy,nub)
 import Data.Function(on)
+import qualified Data.Map.Strict as Map
 
 --import Debug.Trace
 --debug x = trace (show x) x
 
 -----------------------------------------------------
 
--- | Coefficient of a variable for a sparse representation of constraints.
+-- | Coefficient of a variable for a sparse and general representations of constraints.
 (#) :: Double -> Int -> (Double,Int)
 infixl 5 #
 (#) = (,)
@@ -108,18 +119,29 @@ data Solution = Undefined
               | Unbounded
               deriving Show
 
-data Constraints = Dense  [ Bound [Double] ]
-                 | Sparse [ Bound [(Double,Int)] ]
+data Constraints = Dense   [ Bound [Double] ]
+                 | Sparse  [ Bound [(Double,Int)] ]
+                 | General [ Bound [(Double,Int)] ]
 
 data Optimization = Maximize [Double]
                   | Minimize [Double]
 
 type Bounds = [Bound Int]
 
+-- | Convert a system of General constraints to one with unique coefficients. 
+sparseOfGeneral :: Constraints -> Constraints
+sparseOfGeneral (General cs) =
+    Sparse $ map (\bl -> 
+                      let cl = obj bl in
+                      let cl_unique = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in
+                      withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] cl_unique)) cs
+sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected"
+
 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   []) bnds = simplex opt (Sparse []) bnds
+simplex opt (Sparse  []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+simplex opt (General []) 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)
@@ -133,6 +155,29 @@ simplex opt (Sparse constr) bnds = extract sg sol where
     m = length constr
     (sz, sg, objfun) = adapt opt
 
+simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds
+
+-- | Simplex method with exact internal arithmetic. See glp_exact in glpk documentation for more information.
+exact :: Optimization -> Constraints -> Bounds -> Solution
+
+exact opt (Dense   []) bnds = exact opt (Sparse []) bnds
+exact opt (Sparse  []) bnds = exact opt (Sparse [Free [0#1]]) bnds
+exact opt (General []) bnds = exact opt (Sparse [Free [0#1]]) bnds
+
+exact opt (Dense constr) bnds = extract sg sol where
+    sol = exactSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
+    n = length objfun
+    m = length constr
+    (sz, sg, objfun) = adapt opt
+
+exact opt (Sparse constr) bnds = extract sg sol where
+    sol = exactSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)
+    n = length objfun
+    m = length constr
+    (sz, sg, objfun) = adapt opt
+
+exact opt constr@(General _) bnds = exact opt (sparseOfGeneral constr) bnds
+
 adapt :: Optimization -> (Int, Double, [Double])
 adapt opt = case opt of
     Maximize x -> (size x, 1 ,x)
@@ -162,6 +207,13 @@ obj (x :&: _)  = x
 obj (x :==: _) = x
 obj (Free x)   = x
 
+withObj :: Bound t -> t -> Bound t
+withObj (_ :<=: b) x = (x :<=: b)
+withObj (_ :>=: b) x = (x :>=: b)
+withObj (_ :&: b) x = (x :&: b)
+withObj (_ :==: b) x = (x :==: b)
+withObj (Free _) x = Free x
+
 tb :: Bound t -> Double
 tb (_ :<=: _)  = glpUP
 tb (_ :>=: _)  = glpLO
@@ -235,6 +287,19 @@ simplexSparse m n c b = unsafePerformIO $ do
     app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
     return s
 
+foreign import ccall unsafe "c_exact_sparse" c_exact_sparse
+    :: CInt -> CInt                  -- rows and cols
+    -> CInt -> CInt -> Ptr Double    -- coeffs
+    -> CInt -> CInt -> Ptr Double    -- bounds
+    -> CInt -> Ptr Double            -- result
+    -> IO CInt                       -- exit code
+
+exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
+exactSparse m n c b = unsafePerformIO $ do
+    s <- createVector (2+n)
+    app3 (c_exact_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_exact_sparse"
+    return s
+
 glpFR, glpLO, glpUP, glpDB, glpFX :: Double
 glpFR = 0
 glpLO = 1