HSvm 0.1.0.2.90 → 0.1.0.3.22
raw patch · 6 files changed
+632/−427 lines, 6 filesPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
- Data.SVM: Model :: (ForeignPtr CSvmModel) -> Model
- Data.SVM: checkParam :: Ptr CSvmProblem -> Ptr CSvmParameter -> IO ()
- Data.SVM: defaultExtra :: ExtraParam
- Data.SVM: freeCSVmProblem :: Ptr CSvmProblem -> IO ()
- Data.SVM: mergeAlgo :: Algorithm -> CSvmParameter -> CSvmParameter
- Data.SVM: mergeExtra :: ExtraParam -> CSvmParameter -> CSvmParameter
- Data.SVM: mergeKernel :: KernelType -> CSvmParameter -> CSvmParameter
- Data.SVM: newCSvmNodeArray :: Vector -> IO (Ptr CSvmNode)
- Data.SVM: newCSvmProblem :: Problem -> IO (Ptr CSvmProblem)
- Data.SVM: newtype Model
- Data.SVM: withParam :: ExtraParam -> Algorithm -> KernelType -> (Ptr CSvmParameter -> IO a) -> IO a
- Data.SVM: withProblem :: Problem -> (Ptr CSvmProblem -> IO a) -> IO a
+ Data.SVM: data Model
- Data.SVM: loadModel :: FilePath -> IO (Model)
+ Data.SVM: loadModel :: FilePath -> IO Model
- Data.SVM: predict :: Model -> Vector -> Double
+ Data.SVM: predict :: Model -> Vector -> IO Double
- Data.SVM: train :: Algorithm -> KernelType -> Problem -> IO (Model)
+ Data.SVM: train :: Algorithm -> KernelType -> Problem -> IO Model
- Data.SVM: train' :: ExtraParam -> Algorithm -> KernelType -> Problem -> IO (Model)
+ Data.SVM: train' :: ExtraParam -> Algorithm -> KernelType -> Problem -> IO Model
- Data.SVM.Raw: c_clone_model_support_vectors :: Ptr CSvmModel -> IO ()
+ Data.SVM.Raw: c_clone_model_support_vectors :: Ptr CSvmModel -> IO CInt
Files
- Data/SVM.hs +101/−75
- Data/SVM/Raw.hsc +59/−58
- HSvm.cabal +6/−5
- LICENSE +1/−1
- cbits/svm.cpp +431/−284
- cbits/svm.h +34/−4
Data/SVM.hs view
@@ -1,32 +1,48 @@-module Data.SVM where+module Data.SVM+ ( Vector+ , Problem+ , KernelType (..)+ , Algorithm (..)+ , ExtraParam (..)+ , Model+ , train+ , train'+ , crossValidate+ , crossValidate'+ , loadModel+ , saveModel+ , predict+ ) where -- TODO limitare l'export -- TODO verificare l'import -import Control.Arrow ((***))-import Control.Monad (when, liftM)-import Control.Exception -import Data.IntMap (IntMap, toList)-import qualified Data.IntMap as M-import Foreign.Storable (poke, peek)-import Foreign.Marshal.Alloc (malloc, alloca, free)-import Foreign.Marshal.Array-import Foreign.ForeignPtr -import Foreign.Ptr (Ptr, nullPtr)-import Foreign.C.String -import System.IO.Unsafe-import qualified Data.SVM.Raw as R-import Data.SVM.Raw (CSvmModel, CSvmProblem(..), CSvmNode(..), CSvmParameter,- c_svm_train, c_svm_cross_validation,- c_svm_destroy_model, c_svm_check_parameter,- c_svm_load_model, c_svm_save_model, c_svm_predict,- c_clone_model_support_vectors, defaultCParam)+-- import Control.Arrow ((***))+import Control.Exception+import Control.Monad (liftM, when)+import Data.IntMap (IntMap, toList)+import qualified Data.IntMap as M+import Data.SVM.Raw (CSvmModel, CSvmNode (..), CSvmParameter,+ CSvmProblem (..),+ c_clone_model_support_vectors,+ c_svm_check_parameter,+ c_svm_cross_validation,+ c_svm_destroy_model, c_svm_load_model,+ c_svm_predict, c_svm_save_model,+ c_svm_train, defaultCParam)+import qualified Data.SVM.Raw as R+import Foreign.C.String+import Foreign.ForeignPtr+import Foreign.Marshal.Alloc (alloca, free, malloc)+import Foreign.Marshal.Array+import Foreign.Ptr (Ptr, nullPtr)+import Foreign.Storable (peek, poke) type Vector = IntMap Double type Problem = [(Double, Vector)] newtype Model = Model (ForeignPtr CSvmModel) -data KernelType = Linear +data KernelType = Linear | RBF { gamma :: Double } | Sigmoid { gamma :: Double, coef0 :: Double } | Poly { gamma :: Double, coef0 :: Double, degree :: Int}@@ -37,77 +53,89 @@ | EpsilonSvr { epsilon :: Double, c :: Double } | OneClassSvm { nu :: Double } -data ExtraParam = ExtraParam {cacheSize :: Double, - shrinking :: Int, +data ExtraParam = ExtraParam {cacheSize :: Double,+ shrinking :: Int, probability :: Int} -defaultExtra = ExtraParam {cacheSize = 100, shrinking = 1, probability = 0}+defaultExtra :: ExtraParam+defaultExtra = ExtraParam {cacheSize = 1000, shrinking = 1, probability = 0} mergeKernel :: KernelType -> CSvmParameter -> CSvmParameter mergeKernel Linear p = p { R.kernel_type = R.linear } mergeKernel (RBF g) p = p { R.kernel_type = R.rbf, R.gamma = realToFrac g }-mergeKernel (Sigmoid g c) p = p { R.kernel_type = R.sigmoid,- R.gamma = realToFrac g, - R.coef0 = realToFrac c }-mergeKernel (Poly g c d) p = p { R.kernel_type = R.poly, - R.gamma = realToFrac g, - R.coef0 = realToFrac c, +mergeKernel (Sigmoid g cf) p = p { R.kernel_type = R.sigmoid,+ R.gamma = realToFrac g,+ R.coef0 = realToFrac cf }+mergeKernel (Poly g cf d) p = p { R.kernel_type = R.poly,+ R.gamma = realToFrac g,+ R.coef0 = realToFrac cf, R.degree = fromIntegral d} mergeAlgo :: Algorithm -> CSvmParameter -> CSvmParameter-mergeAlgo (CSvc c) p = p { R.svm_type = R.cSvc, - R.c = realToFrac c }-mergeAlgo (NuSvc nu) p = p { R.svm_type = R.nuSvc, - R.nu = realToFrac nu }-mergeAlgo (NuSvr nu c) p = p { R.svm_type = R.nuSvr, - R.nu = realToFrac nu, - R.c = realToFrac c }-mergeAlgo (EpsilonSvr e c) p = p { R.svm_type = R.epsilonSvr, - R.eps = realToFrac e, - R.c = realToFrac c }+mergeAlgo (CSvc cf) p = p { R.svm_type = R.cSvc,+ R.c = realToFrac cf }+mergeAlgo (NuSvc n) p = p { R.svm_type = R.nuSvc,+ R.nu = realToFrac n }+mergeAlgo (NuSvr n cf) p = p { R.svm_type = R.nuSvr,+ R.nu = realToFrac n,+ R.c = realToFrac cf }+mergeAlgo (EpsilonSvr e cf) p = p { R.svm_type = R.epsilonSvr,+ R.eps = realToFrac e,+ R.c = realToFrac cf }+mergeAlgo (OneClassSvm n) p = p { R.svm_type = R.oneClass,+ R.nu = realToFrac n } -mergeExtra (ExtraParam c s pr) p = p { R.cache_size = realToFrac c,+mergeExtra :: ExtraParam -> CSvmParameter -> CSvmParameter+mergeExtra (ExtraParam cf s pr) p = p { R.cache_size = realToFrac cf, R.shrinking = fromIntegral s, R.probability = fromIntegral pr } ------------------------------------------------------------------------------- +convertToNodeArray :: Vector -> [CSvmNode]+convertToNodeArray = map convertNode . toList . M.filter (/= 0)+ where+ convertNode (key, val) = CSvmNode (fromIntegral key) (realToFrac val)++endMarker :: CSvmNode+endMarker = CSvmNode (-1) 0.0+ newCSvmNodeArray :: Vector -> IO (Ptr CSvmNode)-newCSvmNodeArray v = newArray (convertVector v ++ [CSvmNode (-1) 0])- where convertVector :: Vector -> [CSvmNode]- convertVector = map convertNode . toList . M.filter (/= 0)- convertNode = uncurry CSvmNode . (fromIntegral *** realToFrac)+newCSvmNodeArray v = newArray0 endMarker (convertToNodeArray v) +withCSvmNodeArray :: Vector -> (Ptr CSvmNode -> IO a) -> IO a+withCSvmNodeArray v = withArray0 endMarker (convertToNodeArray v)+ newCSvmProblem :: Problem -> IO (Ptr CSvmProblem)-newCSvmProblem lvs = do nodePtrList <- mapM newCSvmNodeArray $ map snd lvs+newCSvmProblem lvs = do nodePtrList <- mapM (newCSvmNodeArray . snd) lvs nodePtrPtr <- newArray nodePtrList labelPtr <- newArray . map realToFrac $ map fst lvs- let l = fromIntegral . length $ lvs+ let z = fromIntegral . length $ lvs ptr <- malloc- poke ptr $ CSvmProblem l labelPtr nodePtrPtr+ poke ptr $ CSvmProblem z labelPtr nodePtrPtr return ptr freeCSVmProblem :: Ptr CSvmProblem -> IO () freeCSVmProblem ptr = do prob <- peek ptr free $ y prob vecList <- peekArray (fromIntegral $ l prob) (x prob)- mapM_ free vecList + mapM_ free vecList free $ x prob free ptr withProblem :: Problem -> (Ptr CSvmProblem -> IO a) -> IO a-withProblem prob = bracket (newCSvmProblem prob) freeCSVmProblem +withProblem prob = bracket (newCSvmProblem prob) freeCSVmProblem --- -withParam :: ExtraParam - -> Algorithm - -> KernelType - -> (Ptr CSvmParameter -> IO a) +withParam :: ExtraParam+ -> Algorithm+ -> KernelType+ -> (Ptr CSvmParameter -> IO a) -> IO a-withParam extra algo kern f = - let merge = mergeAlgo algo . mergeKernel kern . mergeExtra extra +withParam extra algo kern f =+ let merge = mergeAlgo algo . mergeKernel kern . mergeExtra extra param = merge defaultCParam in alloca $ \paramPtr -> poke paramPtr param >> f paramPtr @@ -118,29 +146,27 @@ -- -train' :: ExtraParam -> Algorithm -> KernelType -> Problem -> IO (Model)-train' extra algo kern prob = +train' :: ExtraParam -> Algorithm -> KernelType -> Problem -> IO Model+train' extra algo kern prob = withProblem prob $ \probPtr -> withParam extra algo kern $ \paramPtr -> do checkParam probPtr paramPtr- -- rereadParam <- peek paramPtr- -- print rereadParam modelPtr <- c_svm_train probPtr paramPtr- c_clone_model_support_vectors modelPtr+ _ <- c_clone_model_support_vectors modelPtr modelForeignPtr <- newForeignPtr c_svm_destroy_model modelPtr return $ Model modelForeignPtr -- | The 'train' function allows training a 'Model' starting from a 'Problem' -- by specifying an 'Algorithm' and a 'KernelType'-train :: Algorithm -> KernelType -> Problem -> IO (Model)+train :: Algorithm -> KernelType -> Problem -> IO Model train = train' defaultExtra -crossValidate' :: ExtraParam - -> Algorithm - -> KernelType - -> Problem - -> Int +crossValidate' :: ExtraParam+ -> Algorithm+ -> KernelType+ -> Problem+ -> Int -> IO [Double] crossValidate' extra algo kern prob nFold = withProblem prob $ \probPtr ->@@ -152,27 +178,27 @@ c_svm_cross_validation probPtr paramPtr c_nFold targetPtr map realToFrac `liftM` peekArray probLen targetPtr +crossValidate :: Algorithm -> KernelType -> Problem -> Int -> IO [Double] crossValidate = crossValidate' defaultExtra ----------------------------------------------------------------------- saveModel :: Model -> FilePath -> IO ()-saveModel (Model modelForeignPtr) path = +saveModel (Model modelForeignPtr) path = withForeignPtr modelForeignPtr $ \modelPtr -> do pathString <- newCString path ret <- c_svm_save_model pathString modelPtr- when (ret /= 0) $ error "svm: error saving the model"+ when (ret /= 0) $ error $ "svm: error saving the model:" ++ show ret -loadModel :: FilePath -> IO (Model)+loadModel :: FilePath -> IO Model loadModel path = do- modelPtr <- c_svm_load_model =<< newCString path + modelPtr <- c_svm_load_model =<< newCString path Model `liftM` newForeignPtr c_svm_destroy_model modelPtr ----predict :: Model -> Vector -> Double-predict (Model modelForeignPtr) vector = unsafePerformIO action+predict :: Model -> Vector -> IO Double+predict (Model modelForeignPtr) vector = action where action :: IO Double- action = withForeignPtr modelForeignPtr $ \modelPtr -> - bracket (newCSvmNodeArray vector) free $ \vectorPtr ->+ action = withForeignPtr modelForeignPtr $ \modelPtr ->+ withCSvmNodeArray vector $ \vectorPtr -> return . realToFrac . c_svm_predict modelPtr $ vectorPtr-
Data/SVM/Raw.hsc view
@@ -1,4 +1,4 @@-{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving,+{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving, EmptyDataDecls #-} #include "svm.h"@@ -14,22 +14,22 @@ -- TODO verificare l'import import Foreign.Storable (Storable(..), peekByteOff, pokeByteOff)-import Foreign.C.Types (CDouble (..), CInt (..) )+import Foreign.C.Types (CDouble (..), CInt (..)) import Foreign.C.String (CString) import Foreign.Ptr(nullPtr, Ptr) import Foreign.ForeignPtr (FinalizerPtr) -data CSvmNode = CSvmNode {+data CSvmNode = CSvmNode { index:: CInt,- value:: CDouble+ value:: CDouble } instance Storable CSvmNode where sizeOf _ = #size struct svm_node alignment _ = #alignment struct svm_node- peek ptr = do index <- (#peek struct svm_node, index) ptr- value <- (#peek struct svm_node, value) ptr- return $ CSvmNode index value+ peek ptr = do idx <- (#peek struct svm_node, index) ptr+ val <- (#peek struct svm_node, value) ptr+ return $ CSvmNode idx val poke ptr (CSvmNode i v) = do (#poke struct svm_node, index) ptr i (#poke struct svm_node, value) ptr v @@ -37,18 +37,18 @@ l:: CInt, y:: Ptr CDouble, x:: Ptr (Ptr CSvmNode)-}+} instance Storable CSvmProblem where sizeOf _ = #size struct svm_problem alignment _ = #alignment struct svm_problem- peek ptr = do l <- (#peek struct svm_problem, l) ptr- y <- (#peek struct svm_problem, y) ptr- x <- (#peek struct svm_problem, x) ptr- return $ CSvmProblem l y x- poke ptr (CSvmProblem l y x) = do (#poke struct svm_problem, l) ptr l- (#poke struct svm_problem, y) ptr y- (#poke struct svm_problem, x) ptr x+ peek ptr = do lp <- (#peek struct svm_problem, l) ptr+ yp <- (#peek struct svm_problem, y) ptr+ xp <- (#peek struct svm_problem, x) ptr+ return $ CSvmProblem lp yp xp+ poke ptr (CSvmProblem lp yp xp) = do (#poke struct svm_problem, l) ptr lp+ (#poke struct svm_problem, y) ptr yp+ (#poke struct svm_problem, x) ptr xp -- TODO esportare solo il tipo e non il costruttore?@@ -56,7 +56,7 @@ deriving (Storable, Show) #enum CSvmType, CSvmType, C_SVC, NU_SVC, ONE_CLASS, EPSILON_SVR, NU_SVR -newtype CKernelType = CKernelType {unCKernelType :: CInt}+newtype CKernelType = CKernelType {unCKernelType :: CInt} deriving (Storable, Show) #enum CKernelType, CKernelType, LINEAR, POLY, RBF, SIGMOID, PRECOMPUTED @@ -78,64 +78,65 @@ probability :: CInt } deriving Show -defaultCParam = CSvmParameter cSvc rbf 3 0 0 100 1e-3 1+defaultCParam :: CSvmParameter+defaultCParam = CSvmParameter cSvc rbf 3 0 0 100 1e-3 1 0 nullPtr nullPtr 0.5 0.1 1 0 instance Storable CSvmParameter where sizeOf _ = #size struct svm_parameter alignment _ = #alignment struct svm_parameter- peek ptr = do svm_type <- (#peek struct svm_parameter, svm_type) ptr- kernel_type <- (#peek struct svm_parameter, kernel_type) ptr- degree <- (#peek struct svm_parameter, degree) ptr- gamma <- (#peek struct svm_parameter, gamma) ptr- coef0 <- (#peek struct svm_parameter, coef0) ptr- cache_size <- (#peek struct svm_parameter, cache_size) ptr- eps <- (#peek struct svm_parameter, eps) ptr- c <- (#peek struct svm_parameter, C) ptr- nr_weight <- (#peek struct svm_parameter, nr_weight) ptr- weight_label <- (#peek struct svm_parameter, weight_label) ptr- weight <- (#peek struct svm_parameter, weight) ptr- nu <- (#peek struct svm_parameter, nu) ptr- p <- (#peek struct svm_parameter, p) ptr- shrinking <- (#peek struct svm_parameter, degree) ptr- probability <- (#peek struct svm_parameter, probability) ptr- return $ CSvmParameter svm_type kernel_type degree- gamma coef0 cache_size eps c nr_weight- weight_label weight nu p shrinking probability- poke ptr (CSvmParameter svm_type kernel_type degree- gamma coef0 cache_size eps c nr_weight- weight_label weight nu p shrinking probability) =- do (#poke struct svm_parameter, svm_type) ptr svm_type- (#poke struct svm_parameter, kernel_type) ptr kernel_type- (#poke struct svm_parameter, degree) ptr degree- (#poke struct svm_parameter, gamma) ptr gamma- (#poke struct svm_parameter, coef0) ptr coef0- (#poke struct svm_parameter, cache_size) ptr cache_size- (#poke struct svm_parameter, eps) ptr eps- (#poke struct svm_parameter, C) ptr c- (#poke struct svm_parameter, nr_weight) ptr nr_weight- (#poke struct svm_parameter, weight_label) ptr weight_label- (#poke struct svm_parameter, weight) ptr weight- (#poke struct svm_parameter, nu) ptr nu- (#poke struct svm_parameter, p) ptr p- (#poke struct svm_parameter, shrinking) ptr shrinking- (#poke struct svm_parameter, probability) ptr probability+ peek ptr = do svm_type_p <- (#peek struct svm_parameter, svm_type) ptr+ kernel_type_p <- (#peek struct svm_parameter, kernel_type) ptr+ degree_p <- (#peek struct svm_parameter, degree) ptr+ gamma_p <- (#peek struct svm_parameter, gamma) ptr+ coef0_p <- (#peek struct svm_parameter, coef0) ptr+ cache_size_p <- (#peek struct svm_parameter, cache_size) ptr+ eps_p <- (#peek struct svm_parameter, eps) ptr+ c_p <- (#peek struct svm_parameter, C) ptr+ nr_weight_p <- (#peek struct svm_parameter, nr_weight) ptr+ weight_label_p <- (#peek struct svm_parameter, weight_label) ptr+ weight_p <- (#peek struct svm_parameter, weight) ptr+ nu_p <- (#peek struct svm_parameter, nu) ptr+ p_p <- (#peek struct svm_parameter, p) ptr+ shrinking_p <- (#peek struct svm_parameter, degree) ptr+ probability_p <- (#peek struct svm_parameter, probability) ptr+ return $ CSvmParameter svm_type_p kernel_type_p degree_p + gamma_p coef0_p cache_size_p eps_p c_p nr_weight_p+ weight_label_p weight_p nu_p p_p shrinking_p probability_p+ poke ptr (CSvmParameter svm_type_p kernel_type_p degree_p+ gamma_p coef0_p cache_size_p eps_p c_p nr_weight_p+ weight_label_p weight_p nu_p p_p shrinking_p probability_p) =+ do (#poke struct svm_parameter, svm_type) ptr svm_type_p+ (#poke struct svm_parameter, kernel_type) ptr kernel_type_p+ (#poke struct svm_parameter, degree) ptr degree_p+ (#poke struct svm_parameter, gamma) ptr gamma_p+ (#poke struct svm_parameter, coef0) ptr coef0_p+ (#poke struct svm_parameter, cache_size) ptr cache_size_p+ (#poke struct svm_parameter, eps) ptr eps_p+ (#poke struct svm_parameter, C) ptr c_p+ (#poke struct svm_parameter, nr_weight) ptr nr_weight_p+ (#poke struct svm_parameter, weight_label) ptr weight_label_p+ (#poke struct svm_parameter, weight) ptr weight_p+ (#poke struct svm_parameter, nu) ptr nu_p+ (#poke struct svm_parameter, p) ptr p_p+ (#poke struct svm_parameter, shrinking) ptr shrinking_p+ (#poke struct svm_parameter, probability) ptr probability_p data CSvmModel --- TODO cambiare il return type da+-- TODO cambiare il return type da foreign import ccall unsafe "svm.h svm_train" c_svm_train :: Ptr CSvmProblem -> Ptr CSvmParameter -> IO (Ptr CSvmModel)--foreign import ccall unsafe "svm.h svm_cross_validation" c_svm_cross_validation:: Ptr CSvmProblem -> Ptr CSvmParameter -> CInt -> Ptr CDouble -> IO ()+ +foreign import ccall unsafe "svm.h svm_cross_validation" c_svm_cross_validation:: Ptr CSvmProblem -> Ptr CSvmParameter -> CInt -> Ptr CDouble -> IO () foreign import ccall unsafe "svm.h svm_predict" c_svm_predict :: Ptr CSvmModel -> Ptr CSvmNode -> CDouble foreign import ccall unsafe "svm.h svm_save_model" c_svm_save_model :: CString -> Ptr CSvmModel -> IO CInt foreign import ccall unsafe "svm.h svm_load_model" c_svm_load_model :: CString -> IO (Ptr CSvmModel)-+ foreign import ccall unsafe "svm.h svm_check_parameter" c_svm_check_parameter :: Ptr CSvmProblem -> Ptr CSvmParameter -> CString foreign import ccall unsafe "svm.h &svm_destroy_model" c_svm_destroy_model :: FinalizerPtr CSvmModel -foreign import ccall unsafe "svm.h clone_model_support_vectors" c_clone_model_support_vectors :: Ptr CSvmModel -> IO ()+foreign import ccall unsafe "svm.h clone_model_support_vectors" c_clone_model_support_vectors :: Ptr CSvmModel -> IO CInt
HSvm.cabal view
@@ -1,13 +1,13 @@ Name: HSvm-Version: 0.1.0.2.90+Version: 0.1.0.3.22 Copyright: (c) 2009 Paolo Losi, 2017 Pavel Ryzhov-Maintainer: Paolo Losi <paolo.losi@gmail.com>+Maintainer: Pavel Ryzhov <paul@paulrz.cz> License: BSD3 License-File: LICENSE Author: Paolo Losi <paolo.losi@gmail.com> Category: Datamining, Classification Synopsis: Haskell Bindings for libsvm-Stability: alpha+Stability: provisional Build-Type: Simple Cabal-Version: >= 1.6 Extra-Source-Files: cbits/svm.cpp cbits/svm.h@@ -17,10 +17,11 @@ location: https://github.com/paulrzcz/HSvm.git Library- Build-Depends: base >= 4 && < 5, containers+ Build-Depends: base >= 4 && < 5+ , containers Exposed-modules: Data.SVM, Data.SVM.Raw Includes: svm.h Include-Dirs: cbits C-Sources: cbits/svm.cpp Extra-Libraries: stdc++- Ghc-Options: -Wall+ Ghc-Options: -Wall
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2009, Paolo Losi+Copyright (c) 2009, Paolo Losi, Pavel Ryzhov All rights reserved. Redistribution and use in source and binary forms, with or without
cbits/svm.cpp view
@@ -5,23 +5,25 @@ #include <float.h> #include <string.h> #include <stdarg.h>+#include <limits.h>+#include <locale.h> #include "svm.h" int libsvm_version = LIBSVM_VERSION; typedef float Qfloat; typedef signed char schar; #ifndef min-template <class T> inline T min(T x,T y) { return (x<y)?x:y; }+template <class T> static inline T min(T x,T y) { return (x<y)?x:y; } #endif #ifndef max-template <class T> inline T max(T x,T y) { return (x>y)?x:y; }+template <class T> static inline T max(T x,T y) { return (x>y)?x:y; } #endif-template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }-template <class S, class T> inline void clone(T*& dst, S* src, int n)+template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }+template <class S, class T> static inline void clone(T*& dst, S* src, int n) { dst = new T[n]; memcpy((void *)dst,(void *)src,sizeof(T)*n); }-inline double powi(double base, int times)+static inline double powi(double base, int times) { double tmp = base, ret = 1.0; @@ -41,7 +43,7 @@ fputs(s,stdout); fflush(stdout); }-void (*svm_print_string) (const char *) = &print_string_stdout;+static void (*svm_print_string) (const char *) = &print_string_stdout; #if 1 static void info(const char *fmt,...) {@@ -72,7 +74,7 @@ // return some position p where [p,len) need to be filled // (p >= len if nothing needs to be filled) int get_data(const int index, Qfloat **data, int len);- void swap_index(int i, int j); + void swap_index(int i, int j); private: int l; long int size;@@ -192,7 +194,7 @@ class QMatrix { public: virtual Qfloat *get_Q(int column, int len) const = 0;- virtual Qfloat *get_QD() const = 0;+ virtual double *get_QD() const = 0; virtual void swap_index(int i, int j) const = 0; virtual ~QMatrix() {} };@@ -205,7 +207,7 @@ static double k_function(const svm_node *x, const svm_node *y, const svm_parameter& param); virtual Qfloat *get_Q(int column, int len) const = 0;- virtual Qfloat *get_QD() const = 0;+ virtual double *get_QD() const = 0; virtual void swap_index(int i, int j) const // no so const... { swap(x[i],x[j]);@@ -306,11 +308,19 @@ ++py; else ++px;- } + } } return sum; } +void print_node(const svm_node *p) {+ while(p->index != -1) {+ printf("%d:%.8g ",p->index,p->value);+ p++;+ }+ printf("\n");+}+ double Kernel::k_function(const svm_node *x, const svm_node *y, const svm_parameter& param) {@@ -322,6 +332,12 @@ return powi(param.gamma*dot(x,y)+param.coef0,param.degree); case RBF: {+ // printf("x: ");+ // print_node(x);+ // printf("y: ");+ // print_node(y);+ //+ // printf("k_function - rbf: "); double sum = 0; while(x->index != -1 && y->index !=-1) {@@ -335,7 +351,7 @@ else { if(x->index > y->index)- { + { sum += y->value * y->value; ++y; }@@ -358,7 +374,9 @@ sum += y->value * y->value; ++y; }- ++ // printf(" sum=%g\n", sum);+ return exp(-param.gamma*sum); } case SIGMOID:@@ -366,7 +384,7 @@ case PRECOMPUTED: //x: test (validation), y: SV return x[(int)(y->value)].value; default:- return 0; // Unreachable + return 0; // Unreachable } } @@ -412,7 +430,7 @@ char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE double *alpha; const QMatrix *Q;- const Qfloat *QD;+ const double *QD; double eps; double Cp,Cn; double *p;@@ -442,7 +460,7 @@ virtual double calculate_rho(); virtual void do_shrinking(); private:- bool be_shrunk(int i, double Gmax1, double Gmax2); + bool be_shrunk(int i, double Gmax1, double Gmax2); }; void Solver::swap_index(int i, int j)@@ -474,7 +492,7 @@ nr_free++; if(2*nr_free < active_size)- info("\nWarning: using -h 0 may be faster\n");+ info("\nWARNING: using -h 0 may be faster\n"); if (nr_free*l > 2*active_size*(l-active_size)) {@@ -556,9 +574,10 @@ // optimization step int iter = 0;+ int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l); int counter = min(l,1000)+1; - while(1)+ while(iter < max_iter) { // show progress and do shrinking @@ -582,11 +601,11 @@ else counter = 1; // do shrinking next iteration }- + ++iter; // update alpha[i] and alpha[j], handle bounds carefully- + const Qfloat *Q_i = Q.get_Q(i,active_size); const Qfloat *Q_j = Q.get_Q(j,active_size); @@ -598,14 +617,14 @@ if(y[i]!=y[j]) {- double quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];+ double quad_coef = QD[i]+QD[j]+2*Q_i[j]; if (quad_coef <= 0) quad_coef = TAU; double delta = (-G[i]-G[j])/quad_coef; double diff = alpha[i] - alpha[j]; alpha[i] += delta; alpha[j] += delta;- + if(diff > 0) { if(alpha[j] < 0)@@ -641,7 +660,7 @@ } else {- double quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];+ double quad_coef = QD[i]+QD[j]-2*Q_i[j]; if (quad_coef <= 0) quad_coef = TAU; double delta = (G[i]-G[j])/quad_coef;@@ -687,7 +706,7 @@ double delta_alpha_i = alpha[i] - old_alpha_i; double delta_alpha_j = alpha[j] - old_alpha_j;- + for(int k=0;k<active_size;k++) { G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;@@ -725,6 +744,18 @@ } } + if(iter >= max_iter)+ {+ if(active_size < l)+ {+ // reconstruct the whole gradient to calculate objective value+ reconstruct_gradient();+ active_size = l;+ info("*");+ }+ fprintf(stderr,"\nWARNING: reaching max number of iterations\n");+ }+ // calculate rho si->rho = calculate_rho();@@ -775,7 +806,7 @@ // j: minimizes the decrease of obj value // (if quadratic coefficeint <= 0, replace it with tau) // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)- + double Gmax = -INF; double Gmax2 = -INF; int Gmax_idx = -1;@@ -783,7 +814,7 @@ double obj_diff_min = INF; for(int t=0;t<active_size;t++)- if(y[t]==+1) + if(y[t]==+1) { if(!is_upper_bound(t)) if(-G[t] >= Gmax)@@ -818,8 +849,8 @@ Gmax2 = G[j]; if (grad_diff > 0) {- double obj_diff; - double quad_coef=Q_i[i]+QD[j]-2.0*y[i]*Q_i[j];+ double obj_diff;+ double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else@@ -842,8 +873,8 @@ Gmax2 = -G[j]; if (grad_diff > 0) {- double obj_diff; - double quad_coef=Q_i[i]+QD[j]+2.0*y[i]*Q_i[j];+ double obj_diff;+ double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else@@ -859,7 +890,7 @@ } } - if(Gmax+Gmax2 < eps)+ if(Gmax+Gmax2 < eps || Gmin_idx == -1) return 1; out_i = Gmax_idx;@@ -880,7 +911,7 @@ { if(y[i]==+1) return(G[i] > Gmax2);- else + else return(G[i] > Gmax1); } else@@ -896,27 +927,27 @@ // find maximal violating pair first for(i=0;i<active_size;i++) {- if(y[i]==+1) + if(y[i]==+1) {- if(!is_upper_bound(i)) + if(!is_upper_bound(i)) { if(-G[i] >= Gmax1) Gmax1 = -G[i]; }- if(!is_lower_bound(i)) + if(!is_lower_bound(i)) { if(G[i] >= Gmax2) Gmax2 = G[i]; } }- else + else {- if(!is_upper_bound(i)) + if(!is_upper_bound(i)) { if(-G[i] >= Gmax2) Gmax2 = -G[i]; }- if(!is_lower_bound(i)) + if(!is_lower_bound(i)) { if(G[i] >= Gmax1) Gmax1 = G[i];@@ -924,7 +955,7 @@ } } - if(unshrink == false && Gmax1 + Gmax2 <= eps*10) + if(unshrink == false && Gmax1 + Gmax2 <= eps*10) { unshrink = true; reconstruct_gradient();@@ -991,7 +1022,7 @@ // // additional constraint: e^T \alpha = constant //-class Solver_NU : public Solver+class Solver_NU: public Solver { public: Solver_NU() {}@@ -1063,15 +1094,15 @@ { if(y[j]==+1) {- if (!is_lower_bound(j)) + if (!is_lower_bound(j)) { double grad_diff=Gmaxp+G[j]; if (G[j] >= Gmaxp2) Gmaxp2 = G[j]; if (grad_diff > 0) {- double obj_diff; - double quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];+ double obj_diff;+ double quad_coef = QD[ip]+QD[j]-2*Q_ip[j]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else@@ -1094,8 +1125,8 @@ Gmaxn2 = -G[j]; if (grad_diff > 0) {- double obj_diff; - double quad_coef = Q_in[in]+QD[j]-2*Q_in[j];+ double obj_diff;+ double quad_coef = QD[in]+QD[j]-2*Q_in[j]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else@@ -1111,7 +1142,7 @@ } } - if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)+ if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps || Gmin_idx == -1) return 1; if (y[Gmin_idx] == +1)@@ -1129,14 +1160,14 @@ { if(y[i]==+1) return(-G[i] > Gmax1);- else + else return(-G[i] > Gmax4); } else if(is_lower_bound(i)) { if(y[i]==+1) return(G[i] > Gmax2);- else + else return(G[i] > Gmax3); } else@@ -1165,14 +1196,14 @@ if(!is_lower_bound(i)) { if(y[i]==+1)- { + { if(G[i] > Gmax2) Gmax2 = G[i]; } else if(G[i] > Gmax3) Gmax3 = G[i]; } } - if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) + if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) { unshrink = true; reconstruct_gradient();@@ -1235,12 +1266,12 @@ r1 = sum_free1/nr_free1; else r1 = (ub1+lb1)/2;- + if(nr_free2 > 0) r2 = sum_free2/nr_free2; else r2 = (ub2+lb2)/2;- + si->r = (r1+r2)/2; return (r1-r2)/2; }@@ -1249,18 +1280,18 @@ // Q matrices for various formulations // class SVC_Q: public Kernel-{ +{ public: SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_) :Kernel(prob.l, prob.x, param) { clone(y,y_,prob.l); cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));- QD = new Qfloat[prob.l];+ QD = new double[prob.l]; for(int i=0;i<prob.l;i++)- QD[i]= (Qfloat)(this->*kernel_function)(i,i);+ QD[i] = (this->*kernel_function)(i,i); }- + Qfloat *get_Q(int i, int len) const { Qfloat *data;@@ -1273,7 +1304,7 @@ return data; } - Qfloat *get_QD() const+ double *get_QD() const { return QD; }@@ -1295,7 +1326,7 @@ private: schar *y; Cache *cache;- Qfloat *QD;+ double *QD; }; class ONE_CLASS_Q: public Kernel@@ -1305,11 +1336,11 @@ :Kernel(prob.l, prob.x, param) { cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));- QD = new Qfloat[prob.l];+ QD = new double[prob.l]; for(int i=0;i<prob.l;i++)- QD[i]= (Qfloat)(this->*kernel_function)(i,i);+ QD[i] = (this->*kernel_function)(i,i); }- + Qfloat *get_Q(int i, int len) const { Qfloat *data;@@ -1322,7 +1353,7 @@ return data; } - Qfloat *get_QD() const+ double *get_QD() const { return QD; }@@ -1341,18 +1372,18 @@ } private: Cache *cache;- Qfloat *QD;+ double *QD; }; class SVR_Q: public Kernel-{ +{ public: SVR_Q(const svm_problem& prob, const svm_parameter& param) :Kernel(prob.l, prob.x, param) { l = prob.l; cache = new Cache(l,(long int)(param.cache_size*(1<<20)));- QD = new Qfloat[2*l];+ QD = new double[2*l]; sign = new schar[2*l]; index = new int[2*l]; for(int k=0;k<l;k++)@@ -1361,8 +1392,8 @@ sign[k+l] = -1; index[k] = k; index[k+l] = k;- QD[k]= (Qfloat)(this->*kernel_function)(k,k);- QD[k+l]=QD[k];+ QD[k] = (this->*kernel_function)(k,k);+ QD[k+l] = QD[k]; } buffer[0] = new Qfloat[2*l]; buffer[1] = new Qfloat[2*l];@@ -1375,7 +1406,7 @@ swap(index[i],index[j]); swap(QD[i],QD[j]); }- + Qfloat *get_Q(int i, int len) const { Qfloat *data;@@ -1395,7 +1426,7 @@ return buf; } - Qfloat *get_QD() const+ double *get_QD() const { return QD; }@@ -1416,7 +1447,7 @@ int *index; mutable int next_buffer; Qfloat *buffer[2];- Qfloat *QD;+ double *QD; }; //@@ -1436,7 +1467,7 @@ { alpha[i] = 0; minus_ones[i] = -1;- if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;+ if(prob->y[i] > 0) y[i] = +1; else y[i] = -1; } Solver s;@@ -1626,10 +1657,10 @@ struct decision_function { double *alpha;- double rho; + double rho; }; -decision_function svm_train_one(+static decision_function svm_train_one( const svm_problem *prob, const svm_parameter *param, double Cp, double Cn) {@@ -1686,33 +1717,9 @@ return f; } -//-// svm_model-//-struct svm_model-{- svm_parameter param; // parameter- int nr_class; // number of classes, = 2 in regression/one class svm- int l; // total #SV- svm_node **SV; // SVs (SV[l])- double **sv_coef; // coefficients for SVs in decision functions (sv_coef[k-1][l])- double *rho; // constants in decision functions (rho[k*(k-1)/2])- double *probA; // pariwise probability information- double *probB;-- // for classification only-- int *label; // label of each class (label[k])- int *nSV; // number of SVs for each class (nSV[k])- // nSV[0] + nSV[1] + ... + nSV[k-1] = l- // XXX- int free_sv; // 1 if svm_model is created by svm_load_model- // 0 if svm_model is created by svm_train-};- // Platt's binary SVM Probablistic Output: an improvement from Lin et al.-void sigmoid_train(- int l, const double *dec_values, const double *labels, +static void sigmoid_train(+ int l, const double *dec_values, const double *labels, double& A, double& B) { double prior1=0, prior0 = 0;@@ -1721,7 +1728,7 @@ for (i=0;i<l;i++) if (labels[i] > 0) prior1+=1; else prior0+=1;- + int max_iter=100; // Maximal number of iterations double min_step=1e-10; // Minimal step taken in line search double sigma=1e-12; // For numerically strict PD of Hessian@@ -1731,8 +1738,8 @@ double *t=Malloc(double,l); double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize; double newA,newB,newf,d1,d2;- int iter; - + int iter;+ // Initial Point and Initial Fun Value A=0.0; B=log((prior0+1.0)/(prior1+1.0)); double fval = 0.0;@@ -1824,9 +1831,10 @@ free(t); } -double sigmoid_predict(double decision_value, double A, double B)+static double sigmoid_predict(double decision_value, double A, double B) { double fApB = decision_value*A+B;+ // 1-p used later; avoid catastrophic cancellation if (fApB >= 0) return exp(-fApB)/(1.0+exp(-fApB)); else@@ -1834,14 +1842,14 @@ } // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng-void multiclass_probability(int k, double **r, double *p)+static void multiclass_probability(int k, double **r, double *p) { int t,j; int iter = 0, max_iter=max(100,k); double **Q=Malloc(double *,k); double *Qp=Malloc(double,k); double pQp, eps=0.005/k;- + for (t=0;t<k;t++) { p[t]=1.0/k; // Valid if k = 1@@ -1877,7 +1885,7 @@ max_error=error; } if (max_error<eps) break;- + for (t=0;t<k;t++) { double diff=(-Qp[t]+pQp)/Q[t][t];@@ -1898,7 +1906,7 @@ } // Cross-validation decision values for probability estimates-void svm_binary_svc_probability(+static void svm_binary_svc_probability( const svm_problem *prob, const svm_parameter *param, double Cp, double Cn, double& probA, double& probB) {@@ -1924,7 +1932,7 @@ subprob.l = prob->l-(end-begin); subprob.x = Malloc(struct svm_node*,subprob.l); subprob.y = Malloc(double,subprob.l);- + k=0; for(j=0;j<begin;j++) {@@ -1969,23 +1977,23 @@ struct svm_model *submodel = svm_train(&subprob,&subparam); for(j=begin;j<end;j++) {- svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); + svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); // ensure +1 -1 order; reason not using CV subroutine dec_values[perm[j]] *= submodel->label[0];- } - svm_destroy_model(submodel);+ }+ svm_free_and_destroy_model(&submodel); svm_destroy_param(&subparam); } free(subprob.x); free(subprob.y);- } + } sigmoid_train(prob->l,dec_values,prob->y,probA,probB); free(dec_values); free(perm); } -// Return parameter of a Laplace distribution -double svm_svr_probability(+// Return parameter of a Laplace distribution+static double svm_svr_probability( const svm_problem *prob, const svm_parameter *param) { int i;@@ -2000,15 +2008,15 @@ { ymv[i]=prob->y[i]-ymv[i]; mae += fabs(ymv[i]);- } + } mae /= prob->l; double std=sqrt(2*mae*mae); int count=0; mae=0; for(i=0;i<prob->l;i++)- if (fabs(ymv[i]) > 5*std) + if (fabs(ymv[i]) > 5*std) count=count+1;- else + else mae+=fabs(ymv[i]); mae /= (prob->l-count); info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);@@ -2019,14 +2027,14 @@ // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data // perm, length l, must be allocated before calling this subroutine-void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)+static void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm) { int l = prob->l; int max_nr_class = 16; int nr_class = 0; int *label = Malloc(int,max_nr_class); int *count = Malloc(int,max_nr_class);- int *data_label = Malloc(int,l); + int *data_label = Malloc(int,l); int i; for(i=0;i<l;i++)@@ -2056,6 +2064,24 @@ } } + //+ // Labels are ordered by their first occurrence in the training set.+ // However, for two-class sets with -1/+1 labels and -1 appears first,+ // we swap labels to ensure that internally the binary SVM has positive data corresponding to the +1 instances.+ //+ if (nr_class == 2 && label[0] == -1 && label[1] == 1)+ {+ swap(label[0],label[1]);+ swap(count[0],count[1]);+ for(i=0;i<l;i++)+ {+ if(data_label[i] == 0)+ data_label[i] = 1;+ else+ data_label[i] = 0;+ }+ }+ int *start = Malloc(int,nr_class); start[0] = 0; for(i=1;i<nr_class;i++)@@ -2089,6 +2115,9 @@ param->svm_type == EPSILON_SVR || param->svm_type == NU_SVR) {+ printf("SVM type %d\n", param->svm_type);+ printf("Kernel type %d\n", param->kernel_type);+ printf("Gamma %f\n", param->gamma); // regression or one-class-svm model->nr_class = 2; model->label = NULL;@@ -2096,7 +2125,7 @@ model->probA = NULL; model->probB = NULL; model->sv_coef = Malloc(double *,1); - if(param->probability && + if(param->probability && (param->svm_type == EPSILON_SVR || param->svm_type == NU_SVR)) {@@ -2115,14 +2144,16 @@ model->l = nSV; model->SV = Malloc(svm_node *,nSV); model->sv_coef[0] = Malloc(double,nSV);+ model->sv_indices = Malloc(int,nSV); int j = 0; for(i=0;i<prob->l;i++) if(fabs(f.alpha[i]) > 0) { model->SV[j] = prob->x[i]; model->sv_coef[0][j] = f.alpha[i];+ model->sv_indices[j] = i+1; ++j;- } + } free(f.alpha); }@@ -2137,7 +2168,10 @@ int *perm = Malloc(int,l); // group training data of the same class- svm_group_classes(prob,&nr_class,&label,&start,&count,perm); + svm_group_classes(prob,&nr_class,&label,&start,&count,perm);+ if(nr_class == 1)+ info("WARNING: training data in only one class. See README for details.\n");+ svm_node **x = Malloc(svm_node *,l); int i; for(i=0;i<l;i++)@@ -2149,19 +2183,19 @@ for(i=0;i<nr_class;i++) weighted_C[i] = param->C; for(i=0;i<param->nr_weight;i++)- { + { int j; for(j=0;j<nr_class;j++) if(param->weight_label[i] == label[j]) break; if(j == nr_class)- fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);+ fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]); else weighted_C[j] *= param->weight[i]; } // train k*(k-1)/2 models- + bool *nonzero = Malloc(bool,l); for(i=0;i<l;i++) nonzero[i] = false;@@ -2214,11 +2248,11 @@ // build output model->nr_class = nr_class;- + model->label = Malloc(int,nr_class); for(i=0;i<nr_class;i++) model->label[i] = label[i];- + model->rho = Malloc(double,nr_class*(nr_class-1)/2); for(i=0;i<nr_class*(nr_class-1)/2;i++) model->rho[i] = f[i].rho;@@ -2247,21 +2281,26 @@ int nSV = 0; for(int j=0;j<count[i];j++) if(nonzero[start[i]+j])- { + { ++nSV; ++total_sv; } model->nSV[i] = nSV; nz_count[i] = nSV; }- + info("Total nSV = %d\n",total_sv); model->l = total_sv; model->SV = Malloc(svm_node *,total_sv);+ model->sv_indices = Malloc(int,total_sv); p = 0; for(i=0;i<l;i++)- if(nonzero[i]) model->SV[p++] = x[i];+ if(nonzero[i])+ {+ model->SV[p] = x[i];+ model->sv_indices[p++] = perm[i] + 1;+ } int *nz_start = Malloc(int,nr_class); nz_start[0] = 0;@@ -2284,7 +2323,7 @@ int sj = start[j]; int ci = count[i]; int cj = count[j];- + int q = nz_start[i]; int k; for(k=0;k<ci;k++)@@ -2296,7 +2335,7 @@ model->sv_coef[i][q++] = f[p].alpha[ci+k]; ++p; }- + free(label); free(probA); free(probB);@@ -2319,11 +2358,16 @@ void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target) { int i;- int *fold_start = Malloc(int,nr_fold+1);+ int *fold_start; int l = prob->l; int *perm = Malloc(int,l); int nr_class;-+ if (nr_fold > l)+ {+ nr_fold = l;+ fprintf(stderr,"WARNING: # folds > # data. Will use # folds = # data instead (i.e., leave-one-out cross validation)\n");+ }+ fold_start = Malloc(int,nr_fold+1); // stratified cv may not give leave-one-out rate // Each class to l folds -> some folds may have zero elements if((param->svm_type == C_SVC ||@@ -2340,7 +2384,7 @@ int *index = Malloc(int,l); for(i=0;i<l;i++) index[i]=perm[i];- for (c=0; c<nr_class; c++) + for (c=0; c<nr_class; c++) for(i=0;i<count[c];i++) { int j = i+rand()%(count[c]-i);@@ -2369,9 +2413,9 @@ fold_start[0]=0; for (i=1;i<=nr_fold;i++) fold_start[i] = fold_start[i-1]+fold_count[i-1];- free(start); + free(start); free(label);- free(count); + free(count); free(index); free(fold_count); }@@ -2397,7 +2441,7 @@ subprob.l = l-(end-begin); subprob.x = Malloc(struct svm_node*,subprob.l); subprob.y = Malloc(double,subprob.l);- + k=0; for(j=0;j<begin;j++) {@@ -2412,23 +2456,23 @@ ++k; } struct svm_model *submodel = svm_train(&subprob,param);- if(param->probability && + if(param->probability && (param->svm_type == C_SVC || param->svm_type == NU_SVC)) { double *prob_estimates=Malloc(double,svm_get_nr_class(submodel)); for(j=begin;j<end;j++) target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);- free(prob_estimates); + free(prob_estimates); } else for(j=begin;j<end;j++) target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);- svm_destroy_model(submodel);+ svm_free_and_destroy_model(&submodel); free(subprob.x); free(subprob.y);- } + } free(fold_start);- free(perm); + free(perm); } @@ -2449,6 +2493,18 @@ label[i] = model->label[i]; } +void svm_get_sv_indices(const svm_model *model, int* indices)+{+ if (model->sv_indices != NULL)+ for(int i=0;i<model->l;i++)+ indices[i] = model->sv_indices[i];+}++int svm_get_nr_sv(const svm_model *model)+{+ return model->l;+}+ double svm_get_svr_probability(const svm_model *model) { if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&@@ -2461,25 +2517,33 @@ } } -void svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)+double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values) {+ int i; if(model->param.svm_type == ONE_CLASS || model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) { double *sv_coef = model->sv_coef[0]; double sum = 0;- for(int i=0;i<model->l;i++)+ for(i=0;i<model->l;i++) { sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);+ } sum -= model->rho[0]; *dec_values = sum;++ printf("%g\n", sum);++ if(model->param.svm_type == ONE_CLASS)+ return (sum>0)?1:-1;+ else+ return sum; } else {- int i; int nr_class = model->nr_class; int l = model->l;- + double *kvalue = Malloc(double,l); for(i=0;i<l;i++) kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);@@ -2489,6 +2553,10 @@ for(i=1;i<nr_class;i++) start[i] = start[i-1]+model->nSV[i-1]; + int *vote = Malloc(int,nr_class);+ for(i=0;i<nr_class;i++)+ vote[i] = 0;+ int p=0; for(i=0;i<nr_class;i++) for(int j=i+1;j<nr_class;j++)@@ -2498,7 +2566,7 @@ int sj = start[j]; int ci = model->nSV[i]; int cj = model->nSV[j];- + int k; double *coef1 = model->sv_coef[j-1]; double *coef2 = model->sv_coef[i];@@ -2508,56 +2576,40 @@ sum += coef2[sj+k] * kvalue[sj+k]; sum -= model->rho[p]; dec_values[p] = sum;++ if(dec_values[p] > 0)+ ++vote[i];+ else+ ++vote[j]; p++; } + int vote_max_idx = 0;+ for(i=1;i<nr_class;i++)+ if(vote[i] > vote[vote_max_idx])+ vote_max_idx = i;+ free(kvalue); free(start);+ free(vote);+ return model->label[vote_max_idx]; } } double svm_predict(const svm_model *model, const svm_node *x) {+ print_node(x);+ int nr_class = model->nr_class;+ double *dec_values; if(model->param.svm_type == ONE_CLASS || model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR)- {- double res;- svm_predict_values(model, x, &res);- - if(model->param.svm_type == ONE_CLASS)- return (res>0)?1:-1;- else- return res;- }+ dec_values = Malloc(double, 1); else- {- int i;- int nr_class = model->nr_class;- double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);- svm_predict_values(model, x, dec_values);-- int *vote = Malloc(int,nr_class);- for(i=0;i<nr_class;i++)- vote[i] = 0;- int pos=0;- for(i=0;i<nr_class;i++)- for(int j=i+1;j<nr_class;j++)- {- if(dec_values[pos++] > 0)- ++vote[i];- else- ++vote[j];- }-- int vote_max_idx = 0;- for(i=1;i<nr_class;i++)- if(vote[i] > vote[vote_max_idx])- vote_max_idx = i;- free(vote);- free(dec_values);- return model->label[vote_max_idx];- }+ dec_values = Malloc(double, nr_class*(nr_class-1)/2);+ double pred_result = svm_predict_values(model, x, dec_values);+ free(dec_values);+ return pred_result; } double svm_predict_probability(@@ -2583,7 +2635,13 @@ pairwise_prob[j][i]=1-pairwise_prob[i][j]; k++; }- multiclass_probability(nr_class,pairwise_prob,prob_estimates);+ if (nr_class == 2)+ {+ prob_estimates[0] = pairwise_prob[0][1];+ prob_estimates[1] = pairwise_prob[1][0];+ }+ else+ multiclass_probability(nr_class,pairwise_prob,prob_estimates); int prob_max_idx = 0; for(i=1;i<nr_class;i++)@@ -2592,19 +2650,19 @@ for(i=0;i<nr_class;i++) free(pairwise_prob[i]); free(dec_values);- free(pairwise_prob); + free(pairwise_prob); return model->label[prob_max_idx]; }- else + else return svm_predict(model, x); } -const char *svm_type_table[] =+static const char *svm_type_table[] = { "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL }; -const char *kernel_type_table[]=+static const char *kernel_type_table[]= { "linear","polynomial","rbf","sigmoid","precomputed",NULL };@@ -2614,6 +2672,12 @@ FILE *fp = fopen(model_file_name,"w"); if(fp==NULL) return -1; + char *old_locale = setlocale(LC_ALL, NULL);+ if (old_locale) {+ old_locale = strdup(old_locale);+ }+ setlocale(LC_ALL, "C");+ const svm_parameter& param = model->param; fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);@@ -2632,14 +2696,14 @@ int l = model->l; fprintf(fp, "nr_class %d\n", nr_class); fprintf(fp, "total_sv %d\n",l);- + { fprintf(fp, "rho"); for(int i=0;i<nr_class*(nr_class-1)/2;i++) fprintf(fp," %g",model->rho[i]); fprintf(fp, "\n"); }- + if(model->label) { fprintf(fp, "label");@@ -2692,8 +2756,12 @@ } fprintf(fp, "\n"); }++ setlocale(LC_ALL, old_locale);+ free(old_locale);+ if (ferror(fp) != 0 || fclose(fp) != 0) return -1;- else return 0;+ else return 0; } static char *line = NULL;@@ -2717,29 +2785,30 @@ return line; } -svm_model *svm_load_model(const char *model_file_name)+//+// FSCANF helps to handle fscanf failures.+// Its do-while block avoids the ambiguity when+// if (...)+// FSCANF();+// is used+//+#define FSCANF(_stream, _format, _var) do{ if (fscanf(_stream, _format, _var) != 1) return false; }while(0)+bool read_model_header(FILE *fp, svm_model* model) {- FILE *fp = fopen(model_file_name,"rb");- if(fp==NULL) return NULL;- - // read parameters-- svm_model *model = Malloc(svm_model,1); svm_parameter& param = model->param;- model->rho = NULL;- model->probA = NULL;- model->probB = NULL;- model->label = NULL;- model->nSV = NULL;+ // parameters for training only won't be assigned, but arrays are assigned as NULL for safety+ param.nr_weight = 0;+ param.weight_label = NULL;+ param.weight = NULL; char cmd[81]; while(1) {- fscanf(fp,"%80s",cmd);+ FSCANF(fp,"%80s",cmd); if(strcmp(cmd,"svm_type")==0) {- fscanf(fp,"%80s",cmd);+ FSCANF(fp,"%80s",cmd); int i; for(i=0;svm_type_table[i];i++) {@@ -2752,16 +2821,12 @@ if(svm_type_table[i] == NULL) { fprintf(stderr,"unknown svm type.\n");- free(model->rho);- free(model->label);- free(model->nSV);- free(model);- return NULL;+ return false; } } else if(strcmp(cmd,"kernel_type")==0)- { - fscanf(fp,"%80s",cmd);+ {+ FSCANF(fp,"%80s",cmd); int i; for(i=0;kernel_type_table[i];i++) {@@ -2774,78 +2839,108 @@ if(kernel_type_table[i] == NULL) { fprintf(stderr,"unknown kernel function.\n");- free(model->rho);- free(model->label);- free(model->nSV);- free(model);- return NULL;+ return false; } } else if(strcmp(cmd,"degree")==0)- fscanf(fp,"%d",¶m.degree);+ FSCANF(fp,"%d",¶m.degree); else if(strcmp(cmd,"gamma")==0)- fscanf(fp,"%lf",¶m.gamma);+ FSCANF(fp,"%lf",¶m.gamma); else if(strcmp(cmd,"coef0")==0)- fscanf(fp,"%lf",¶m.coef0);+ FSCANF(fp,"%lf",¶m.coef0); else if(strcmp(cmd,"nr_class")==0)- fscanf(fp,"%d",&model->nr_class);+ FSCANF(fp,"%d",&model->nr_class); else if(strcmp(cmd,"total_sv")==0)- fscanf(fp,"%d",&model->l);+ FSCANF(fp,"%d",&model->l); else if(strcmp(cmd,"rho")==0) { int n = model->nr_class * (model->nr_class-1)/2; model->rho = Malloc(double,n); for(int i=0;i<n;i++)- fscanf(fp,"%lf",&model->rho[i]);+ FSCANF(fp,"%lf",&model->rho[i]); } else if(strcmp(cmd,"label")==0) { int n = model->nr_class; model->label = Malloc(int,n); for(int i=0;i<n;i++)- fscanf(fp,"%d",&model->label[i]);+ FSCANF(fp,"%d",&model->label[i]); } else if(strcmp(cmd,"probA")==0) { int n = model->nr_class * (model->nr_class-1)/2; model->probA = Malloc(double,n); for(int i=0;i<n;i++)- fscanf(fp,"%lf",&model->probA[i]);+ FSCANF(fp,"%lf",&model->probA[i]); } else if(strcmp(cmd,"probB")==0) { int n = model->nr_class * (model->nr_class-1)/2; model->probB = Malloc(double,n); for(int i=0;i<n;i++)- fscanf(fp,"%lf",&model->probB[i]);+ FSCANF(fp,"%lf",&model->probB[i]); } else if(strcmp(cmd,"nr_sv")==0) { int n = model->nr_class; model->nSV = Malloc(int,n); for(int i=0;i<n;i++)- fscanf(fp,"%d",&model->nSV[i]);+ FSCANF(fp,"%d",&model->nSV[i]); } else if(strcmp(cmd,"SV")==0) { while(1) { int c = getc(fp);- if(c==EOF || c=='\n') break; + if(c==EOF || c=='\n') break; } break; } else { fprintf(stderr,"unknown text in model file: [%s]\n",cmd);- free(model->rho);- free(model->label);- free(model->nSV);- free(model);- return NULL;+ return false; } } + return true;++}++svm_model *svm_load_model(const char *model_file_name)+{+ FILE *fp = fopen(model_file_name,"rb");+ if(fp==NULL) return NULL;++ char *old_locale = setlocale(LC_ALL, NULL);+ if (old_locale) {+ old_locale = strdup(old_locale);+ }+ setlocale(LC_ALL, "C");++ // read parameters++ svm_model *model = Malloc(svm_model,1);+ model->rho = NULL;+ model->probA = NULL;+ model->probB = NULL;+ model->sv_indices = NULL;+ model->label = NULL;+ model->nSV = NULL;++ // read header+ if (!read_model_header(fp, model))+ {+ fprintf(stderr, "ERROR: fscanf failed to read model\n");+ setlocale(LC_ALL, old_locale);+ free(old_locale);+ free(model->rho);+ free(model->label);+ free(model->nSV);+ free(model);+ return NULL;+ }+ // read sv_coef and SV int elements = 0;@@ -2910,6 +3005,9 @@ } free(line); + setlocale(LC_ALL, old_locale);+ free(old_locale);+ if (ferror(fp) != 0 || fclose(fp) != 0) return NULL; @@ -2917,58 +3015,58 @@ return model; } -void svm_destroy_model(svm_model* model)+void svm_free_model_content(svm_model* model_ptr) {- if(model->free_sv && model->l > 0)- free((void *)(model->SV[0]));- for(int i=0;i<model->nr_class-1;i++)- free(model->sv_coef[i]);- free(model->SV);- free(model->sv_coef);- free(model->rho);- free(model->label);- free(model->probA);- free(model->probB);- free(model->nSV);- free(model);+ if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)+ free((void *)(model_ptr->SV[0]));+ if(model_ptr->sv_coef)+ {+ for(int i=0;i<model_ptr->nr_class-1;i++)+ free(model_ptr->sv_coef[i]);+ }++ free(model_ptr->SV);+ model_ptr->SV = NULL;++ free(model_ptr->sv_coef);+ model_ptr->sv_coef = NULL;++ free(model_ptr->rho);+ model_ptr->rho = NULL;++ free(model_ptr->label);+ model_ptr->label= NULL;++ free(model_ptr->probA);+ model_ptr->probA = NULL;++ free(model_ptr->probB);+ model_ptr->probB= NULL;++ free(model_ptr->sv_indices);+ model_ptr->sv_indices = NULL;++ free(model_ptr->nSV);+ model_ptr->nSV = NULL; } -int clone_model_support_vectors(svm_model* model)+void svm_destroy_model(struct svm_model *model_ptr) {- svm_node **new_sv_array;- svm_node *new_node_array, *curr_node, *new_node;- int i, node_n = 0; - if(model->free_sv) return 0;- model->free_sv = 1;-- for(i = 0; i < model->l; i++) {- curr_node = model->SV[i];- node_n++;- while (curr_node->index != -1) {- node_n++;- curr_node++;- }- }- new_sv_array = Malloc(svm_node *, model->l);- new_node_array = Malloc(svm_node, node_n);- if(new_sv_array == NULL || new_node_array == NULL) return 1;+ if(model_ptr != NULL)+ {+ svm_free_model_content(model_ptr);+ free(model_ptr);+ }+} - new_node = new_node_array;- for(i = 0; i < model->l; i++) {- curr_node = model->SV[i];- new_sv_array[i] = new_node;- while(curr_node->index != -1) {- new_node->index = curr_node->index;- new_node->value = curr_node->value;- new_node++; curr_node++;- }- new_node->index = -1;- new_node->value = 0;- new_node++;- }- free(model->SV);- model->SV = new_sv_array;- return 0;+void svm_free_and_destroy_model(svm_model** model_ptr_ptr)+{+ if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)+ {+ svm_free_model_content(*model_ptr_ptr);+ free(*model_ptr_ptr);+ *model_ptr_ptr = NULL;+ } } void svm_destroy_param(svm_parameter* param)@@ -2988,9 +3086,9 @@ svm_type != EPSILON_SVR && svm_type != NU_SVR) return "unknown svm type";- + // kernel_type, degree- + int kernel_type = param->kernel_type; if(kernel_type != LINEAR && kernel_type != POLY &&@@ -2999,6 +3097,9 @@ kernel_type != PRECOMPUTED) return "unknown kernel type"; + if(param->gamma < 0)+ return "gamma < 0";+ if(param->degree < 0) return "degree of polynomial kernel < 0"; @@ -3040,7 +3141,7 @@ // check whether nu-svc is feasible- + if(svm_type == NU_SVC) { int l = prob->l;@@ -3073,7 +3174,7 @@ ++nr_class; } }- + for(i=0;i<nr_class;i++) { int n1 = count[i];@@ -3101,4 +3202,50 @@ model->probA!=NULL && model->probB!=NULL) || ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) && model->probA!=NULL);+}++void svm_set_print_string_function(void (*print_func)(const char *))+{+ if(print_func == NULL)+ svm_print_string = &print_string_stdout;+ else+ svm_print_string = print_func;+}++int clone_model_support_vectors(svm_model* model)+{+ svm_node **new_sv_array;+ svm_node *new_node_array, *curr_node, *new_node;+ int i, node_n = 0;+ if(model->free_sv) return 0;+ model->free_sv = 1;++ for(i = 0; i < model->l; i++) {+ curr_node = model->SV[i];+ node_n++;+ while (curr_node->index != -1) {+ node_n++;+ curr_node++;+ }+ }+ new_sv_array = Malloc(svm_node *, model->l);+ new_node_array = Malloc(svm_node, node_n);+ if(new_sv_array == NULL || new_node_array == NULL) return 1;++ new_node = new_node_array;+ for(i = 0; i < model->l; i++) {+ curr_node = model->SV[i];+ new_sv_array[i] = new_node;+ while(curr_node->index != -1) {+ new_node->index = curr_node->index;+ new_node->value = curr_node->value;+ new_node++; curr_node++;+ }+ new_node->index = -1;+ new_node->value = 0;+ new_node++;+ }+ free(model->SV);+ model->SV = new_sv_array;+ return 0; }
cbits/svm.h view
@@ -1,7 +1,7 @@ #ifndef _LIBSVM_H #define _LIBSVM_H -#define LIBSVM_VERSION 289+#define LIBSVM_VERSION 322 #ifdef __cplusplus extern "C" {@@ -46,6 +46,31 @@ int probability; /* do probability estimates */ }; +//+// svm_model+//+struct svm_model+{+ struct svm_parameter param; /* parameter */+ int nr_class; /* number of classes, = 2 in regression/one class svm */+ int l; /* total #SV */+ struct svm_node **SV; /* SVs (SV[l]) */+ double **sv_coef; /* coefficients for SVs in decision functions (sv_coef[k-1][l]) */+ double *rho; /* constants in decision functions (rho[k*(k-1)/2]) */+ double *probA; /* pariwise probability information */+ double *probB;+ int *sv_indices; /* sv_indices[0,...,nSV-1] are values in [1,...,num_traning_data] to indicate SVs in the training set */++ /* for classification only */++ int *label; /* label of each class (label[k]) */+ int *nSV; /* number of SVs for each class (nSV[k]) */+ /* nSV[0] + nSV[1] + ... + nSV[k-1] = l */+ /* XXX */+ int free_sv; /* 1 if svm_model is created by svm_load_model*/+ /* 0 if svm_model is created by svm_train */+};+ struct svm_model *svm_train(const struct svm_problem *prob, const struct svm_parameter *param); void svm_cross_validation(const struct svm_problem *prob, const struct svm_parameter *param, int nr_fold, double *target); @@ -55,20 +80,25 @@ int svm_get_svm_type(const struct svm_model *model); int svm_get_nr_class(const struct svm_model *model); void svm_get_labels(const struct svm_model *model, int *label);+void svm_get_sv_indices(const struct svm_model *model, int *sv_indices);+int svm_get_nr_sv(const struct svm_model *model); double svm_get_svr_probability(const struct svm_model *model); -void svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values);+double svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values); double svm_predict(const struct svm_model *model, const struct svm_node *x); double svm_predict_probability(const struct svm_model *model, const struct svm_node *x, double* prob_estimates); -void svm_destroy_model(struct svm_model *model);+void svm_destroy_model(struct svm_model *model_ptr);+void svm_free_model_content(struct svm_model *model_ptr);+void svm_free_and_destroy_model(struct svm_model **model_ptr_ptr); void svm_destroy_param(struct svm_parameter *param); const char *svm_check_parameter(const struct svm_problem *prob, const struct svm_parameter *param); int svm_check_probability_model(const struct svm_model *model); -extern void (*svm_print_string) (const char *);+void svm_set_print_string_function(void (*print_func)(const char *)); +// added by HSvm int clone_model_support_vectors(struct svm_model *model); #ifdef __cplusplus