diff --git a/Data/SVM.hs b/Data/SVM.hs --- a/Data/SVM.hs +++ b/Data/SVM.hs @@ -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 - diff --git a/Data/SVM/Raw.hsc b/Data/SVM/Raw.hsc --- a/Data/SVM/Raw.hsc +++ b/Data/SVM/Raw.hsc @@ -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 diff --git a/HSvm.cabal b/HSvm.cabal --- a/HSvm.cabal +++ b/HSvm.cabal @@ -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 +Maintainer: Pavel Ryzhov License: BSD3 License-File: LICENSE Author: Paolo Losi 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 diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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 diff --git a/cbits/svm.cpp b/cbits/svm.cpp --- a/cbits/svm.cpp +++ b/cbits/svm.cpp @@ -5,23 +5,25 @@ #include #include #include +#include +#include #include "svm.h" int libsvm_version = LIBSVM_VERSION; typedef float Qfloat; typedef signed char schar; #ifndef min -template inline T min(T x,T y) { return (x static inline T min(T x,T y) { return (x inline T max(T x,T y) { return (x>y)?x:y; } +template static inline T max(T x,T y) { return (x>y)?x:y; } #endif -template inline void swap(T& x, T& y) { T t=x; x=y; y=t; } -template inline void clone(T*& dst, S* src, int n) +template static inline void swap(T& x, T& y) { T t=x; x=y; y=t; } +template 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= 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= 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= 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*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*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*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 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;tl-(end-begin); subprob.x = Malloc(struct svm_node*,subprob.l); subprob.y = Malloc(double,subprob.l); - + k=0; for(j=0;jx[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;il;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;isvm_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;il;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;iC; for(i=0;inr_weight;i++) - { + { int j; for(j=0;jweight_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;inr_class = nr_class; - + model->label = Malloc(int,nr_class); for(i=0;ilabel[i] = label[i]; - + model->rho = Malloc(double,nr_class*(nr_class-1)/2); for(i=0;irho[i] = f[i].rho; @@ -2247,21 +2281,26 @@ int nSV = 0; for(int j=0;jnSV[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;iSV[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;ksv_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;iprobability && + 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;jx[perm[j]],prob_estimates); - free(prob_estimates); + free(prob_estimates); } else for(j=begin;jx[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;il;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;il;i++) + for(i=0;il;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;iSV[i],model->param); @@ -2489,6 +2553,10 @@ for(i=1;inSV[i-1]; + int *vote = Malloc(int,nr_class); + for(i=0;inSV[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 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 0) - ++vote[i]; - else - ++vote[j]; - } - - int vote_max_idx = 0; - for(i=1;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;ilabel[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;irho[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;irho[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;ilabel[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;iprobA[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;iprobB[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;inSV[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;inr_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;inr_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;iprobA!=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; } diff --git a/cbits/svm.h b/cbits/svm.h --- a/cbits/svm.h +++ b/cbits/svm.h @@ -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