Convert libsvm format files to LP files

svm2lp.hs

import Control.Monad
import qualified Data.Foldable as F
import Data.Default.Class
import Data.List.Split
import Data.Maybe
import Data.IntMap (IntMap)
import qualified Data.IntMap as IntMap
import Data.Map (Map)
import qualified Data.Map as Map
import qualified ToySolver.Data.MIP as MIP
import System.Environment
import System.FilePath
import System.IO
import Text.Printf

type Problem = [(Int, IntMap Double)]

-- http://ntucsu.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/
loadFile :: FilePath -> IO Problem
loadFile fname = do
  s <- readFile fname
  return $ map f (lines s)
  where
    f :: String -> (Int, IntMap Double)
    f s =
      case words s of
        (y : xs) -> (read (dropWhile ('+'==) y), IntMap.fromList [(read v, read val) | x <- xs, let [v,val] = splitOn ":" x])

primal :: Maybe Double -> Problem -> MIP.Problem
primal c prob
  = def
  { MIP.objectiveFunction = def
      { MIP.objDir = MIP.OptMin
      , MIP.objExpr =
         sum [MIP.constExpr (1/2) * wj * wj | wj <- fmap MIP.varExpr $ IntMap.elems w] +
         sum [MIP.constExpr (realToFrac (fromJust c)) * xi_i | isJust c, xi_i <- fmap MIP.varExpr xi]
      }
  , MIP.constraints =
      [ def{ MIP.constrBody =
               ( MIP.Expr (IntMap.elems (IntMap.intersectionWith (\w_j xs_ij -> MIP.Term (fromIntegral y_i * realToFrac xs_ij) [w_j]) w xs_i))
                 - (MIP.constExpr (fromIntegral y_i) *  MIP.varExpr b - (if isJust c then MIP.varExpr xi_i else 0))
               , MIP.Ge
               , 1
               ) }
      | ((y_i, xs_i), xi_i) <- zip prob xi
      ]
  , MIP.varType = Map.fromList [(x, MIP.ContinuousVariable) | x <- b  : [w_j | w_j <- IntMap.elems w] ++ [xi_i | isJust c, xi_i <- xi]]
  , MIP.varBounds =
      Map.unions
      [ Map.singleton b (MIP.NegInf, MIP.PosInf)
      , Map.fromList [(w_j, (MIP.NegInf, MIP.PosInf)) | w_j <- IntMap.elems w]
      , Map.fromList [(xi_i, (0, MIP.PosInf)) | isJust c, xi_i <- xi]
      ]
  }
  where
    m = length prob
    n = fst $ IntMap.findMax $ IntMap.unions (map snd prob)
    w = IntMap.fromList [(j, MIP.toVar ("w_" ++ show j)) | j <- [1..n]]
    b = MIP.toVar "b"
    xi = [MIP.toVar ("xi_" ++ show i) | i <- [1..m]]

dual
  :: Maybe Double
  -> (IntMap Double -> IntMap Double -> Double)
  -> Problem
  -> MIP.Problem
dual c kernel prob
  = def
  { MIP.objectiveFunction = def
      { MIP.objDir = MIP.OptMax
      , MIP.objExpr = MIP.Expr $
          [MIP.Term 1 [a_i] | a_i <- a] ++
          [ MIP.Term (- (1/2) * fromIntegral (y_i * y_j) * realToFrac (kernel xs_i xs_j)) [a_i, a_j]
          | ((y_i, xs_i), a_i) <- zip prob a
          , ((y_j, xs_j), a_j) <- zip prob a
          ]
      }
  , MIP.constraints =
      [ def{ MIP.constrBody =
               ( MIP.Expr [ MIP.Term (fromIntegral y_i) [a_i]  | ((y_i, xs_i), a_i) <- zip prob a]
               , MIP.Eql
               , 0
               ) }
      ]
  , MIP.varType = Map.fromList [(a_i, MIP.ContinuousVariable) | a_i <- a]
  , MIP.varBounds = Map.fromList [(a_i, (0, if isJust c then MIP.Finite (realToFrac (fromJust c)) else MIP.PosInf)) | a_i <- a]
  }
  where
    m = length prob
    n = fst $ IntMap.findMax $ IntMap.unions (map snd prob)
    a = [MIP.toVar ("a_" ++ show i) | i <- [1..m]]

dot :: IntMap Double -> IntMap Double -> Double
dot a b = sum $ IntMap.elems $ IntMap.intersectionWith (*) a b

gaussian :: Double -> IntMap Double -> IntMap Double -> Double
gaussian sigma a b
  = exp (- F.sum (IntMap.map (**2) (IntMap.unionWith (+) a (IntMap.map negate b))) / (2 * sigma**2))

convert fname = do
  putStrLn fname
  hFlush stdout
  prob <- loadFile fname
  forM_ [(Nothing, ""), (Just 1, "C1-"), (Just 10, "C10-")] $ \(c, cstr) -> do
    let mip = primal c prob
    MIP.writeLPFile (printf "%s-%sprimal.lp" fname cstr) mip
    forM_ [(dot, ""), (gaussian 1, "-gaussian-sigma1"), (gaussian 10, "-gaussian_sigma10")] $ \(kernel, kstr) -> do
      let mip2 = dual c kernel prob
      MIP.writeLPFile (printf "%s-%sdual%s.lp" fname cstr kstr) mip2

main = do
  args <- getArgs
  forM_ args convert