svmpredict.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "../svm.h"

#include "mex.h"
#include "svm_model_matlab.h"

#ifdef MX_API_VER
#if MX_API_VER < 0x07030000
typedef int mwIndex;
#endif
#endif

#define CMD_LEN 2048

int print_null(const char *s,...) {}
int (*info)(const char *fmt,...) = &mexPrintf;

void read_sparse_instance(const mxArray *prhs, int index, struct svm_node *x)
{
        int i, j, low, high;
        mwIndex *ir, *jc;
        double *samples;

        ir = mxGetIr(prhs);
        jc = mxGetJc(prhs);
        samples = mxGetPr(prhs);

        // each column is one instance
        j = 0;
        low = (int)jc[index], high = (int)jc[index+1];
        for(i=low;i<high;i++)
        {
                x[j].index = (int)ir[i] + 1;
                x[j].value = samples[i];
                j++;
        }
        x[j].index = -1;
}

static void fake_answer(mxArray *plhs[])
{
        plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
        plhs[1] = mxCreateDoubleMatrix(0, 0, mxREAL);
        plhs[2] = mxCreateDoubleMatrix(0, 0, mxREAL);
}

void predict(mxArray *plhs[], const mxArray *prhs[], struct svm_model *model, const int predict_probability)
{
        int label_vector_row_num, label_vector_col_num;
        int feature_number, testing_instance_number;
        int instance_index;
        double *ptr_instance, *ptr_label, *ptr_predict_label; 
        double *ptr_prob_estimates, *ptr_dec_values, *ptr;
        struct svm_node *x;
        mxArray *pplhs[1]; // transposed instance sparse matrix

        int correct = 0;
        int total = 0;
        double error = 0;
        double sump = 0, sumt = 0, sumpp = 0, sumtt = 0, sumpt = 0;

        int svm_type=svm_get_svm_type(model);
        int nr_class=svm_get_nr_class(model);
        double *prob_estimates=NULL;

        // prhs[1] = testing instance matrix
        feature_number = (int)mxGetN(prhs[1]);
        testing_instance_number = (int)mxGetM(prhs[1]);
        label_vector_row_num = (int)mxGetM(prhs[0]);
        label_vector_col_num = (int)mxGetN(prhs[0]);

        if(label_vector_row_num!=testing_instance_number)
        {
                mexPrintf("Length of label vector does not match # of instances.\n");
                fake_answer(plhs);
                return;
        }
        if(label_vector_col_num!=1)
        {
                mexPrintf("label (1st argument) should be a vector (# of column is 1).\n");
                fake_answer(plhs);
                return;
        }

        ptr_instance = mxGetPr(prhs[1]);
        ptr_label    = mxGetPr(prhs[0]);

        // transpose instance matrix
        if(mxIsSparse(prhs[1]))
        {
                if(model->param.kernel_type == PRECOMPUTED)
                {
                        // precomputed kernel requires dense matrix, so we make one
                        mxArray *rhs[1], *lhs[1];
                        rhs[0] = mxDuplicateArray(prhs[1]);
                        if(mexCallMATLAB(1, lhs, 1, rhs, "full"))
                        {
                                mexPrintf("Error: cannot full testing instance matrix\n");
                                fake_answer(plhs);
                                return;
                        }
                        ptr_instance = mxGetPr(lhs[0]);
                        mxDestroyArray(rhs[0]);
                }
                else
                {
                        mxArray *pprhs[1];
                        pprhs[0] = mxDuplicateArray(prhs[1]);
                        if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose"))
                        {
                                mexPrintf("Error: cannot transpose testing instance matrix\n");
                                fake_answer(plhs);
                                return;
                        }
                }
        }

        if(predict_probability)
        {
                if(svm_type==NU_SVR || svm_type==EPSILON_SVR)
                        info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=%g\n",svm_get_svr_probability(model));
                else
                        prob_estimates = (double *) malloc(nr_class*sizeof(double));
        }

        plhs[0] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL);
        if(predict_probability)
        {
                // prob estimates are in plhs[2]
                if(svm_type==C_SVC || svm_type==NU_SVC)
                        plhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class, mxREAL);
                else
                        plhs[2] = mxCreateDoubleMatrix(0, 0, mxREAL);
        }
        else
        {
                // decision values are in plhs[2]
                if(svm_type == ONE_CLASS ||
                   svm_type == EPSILON_SVR ||
                   svm_type == NU_SVR ||
                   nr_class == 1) // if only one class in training data, decision values are still returned.
                        plhs[2] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL);
                else
                        plhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class*(nr_class-1)/2, mxREAL);
        }

        ptr_predict_label = mxGetPr(plhs[0]);
        ptr_prob_estimates = mxGetPr(plhs[2]);
        ptr_dec_values = mxGetPr(plhs[2]);
        x = (struct svm_node*)malloc((feature_number+1)*sizeof(struct svm_node) );
        for(instance_index=0;instance_index<testing_instance_number;instance_index++)
        {
                int i;
                double target_label, predict_label;

                target_label = ptr_label[instance_index];

                if(mxIsSparse(prhs[1]) && model->param.kernel_type != PRECOMPUTED) // prhs[1]^T is still sparse
                        read_sparse_instance(pplhs[0], instance_index, x);
                else
                {
                        for(i=0;i<feature_number;i++)
                        {
                                x[i].index = i+1;
                                x[i].value = ptr_instance[testing_instance_number*i+instance_index];
                        }
                        x[feature_number].index = -1;
                }

                if(predict_probability)
                {
                        if(svm_type==C_SVC || svm_type==NU_SVC)
                        {
                                predict_label = svm_predict_probability(model, x, prob_estimates);
                                ptr_predict_label[instance_index] = predict_label;
                                for(i=0;i<nr_class;i++)
                                        ptr_prob_estimates[instance_index + i * testing_instance_number] = prob_estimates[i];
                        } else {
                                predict_label = svm_predict(model,x);
                                ptr_predict_label[instance_index] = predict_label;
                        }
                }
                else
                {
                        if(svm_type == ONE_CLASS ||
                           svm_type == EPSILON_SVR ||
                           svm_type == NU_SVR)
                        {
                                double res;
                                predict_label = svm_predict_values(model, x, &res);
                                ptr_dec_values[instance_index] = res;
                        }
                        else
                        {
                                double *dec_values = (double *) malloc(sizeof(double) * nr_class*(nr_class-1)/2);
                                predict_label = svm_predict_values(model, x, dec_values);
                                if(nr_class == 1) 
                                        ptr_dec_values[instance_index] = 1;
                                else
                                        for(i=0;i<(nr_class*(nr_class-1))/2;i++)
                                                ptr_dec_values[instance_index + i * testing_instance_number] = dec_values[i];
                                free(dec_values);
                        }
                        ptr_predict_label[instance_index] = predict_label;
                }

                if(predict_label == target_label)
                        ++correct;
                error += (predict_label-target_label)*(predict_label-target_label);
                sump += predict_label;
                sumt += target_label;
                sumpp += predict_label*predict_label;
                sumtt += target_label*target_label;
                sumpt += predict_label*target_label;
                ++total;
        }
        if(svm_type==NU_SVR || svm_type==EPSILON_SVR)
        {
                info("Mean squared error = %g (regression)\n",error/total);
                info("Squared correlation coefficient = %g (regression)\n",
                        ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/
                        ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt))
                        );
        }
        else
                info("Accuracy = %g%% (%d/%d) (classification)\n",
                        (double)correct/total*100,correct,total);

        // return accuracy, mean squared error, squared correlation coefficient
        plhs[1] = mxCreateDoubleMatrix(3, 1, mxREAL);
        ptr = mxGetPr(plhs[1]);
        ptr[0] = (double)correct/total*100;
        ptr[1] = error/total;
        ptr[2] = ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/
                                ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt));

        free(x);
        if(prob_estimates != NULL)
                free(prob_estimates);
}

void exit_with_help()
{
        mexPrintf(
                "Usage: [predicted_label, accuracy, decision_values/prob_estimates] = svmpredict(testing_label_vector, testing_instance_matrix, model, 'libsvm_options')\n"
                "Parameters:\n"
                "  model: SVM model structure from svmtrain.\n"
                "  libsvm_options:\n"
                "    -b probability_estimates: whether to predict probability estimates, 0 or 1 (default 0); one-class SVM not supported yet\n"
                "    -q : quiet mode (no outputs)\n"
                "Returns:\n"
                "  predicted_label: SVM prediction output vector.\n"
                "  accuracy: a vector with accuracy, mean squared error, squared correlation coefficient.\n"
                "  prob_estimates: If selected, probability estimate vector.\n"
        );
}

void mexFunction( int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[] )
{
        int prob_estimate_flag = 0;
        struct svm_model *model;
        info = &mexPrintf;

        if(nrhs > 4 || nrhs < 3)
        {
                exit_with_help();
                fake_answer(plhs);
                return;
        }

        if(!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) {
                mexPrintf("Error: label vector and instance matrix must be double\n");
                fake_answer(plhs);
                return;
        }

        if(mxIsStruct(prhs[2]))
        {
                const char *error_msg;

                // parse options
                if(nrhs==4)
                {
                        int i, argc = 1;
                        char cmd[CMD_LEN], *argv[CMD_LEN/2];

                        // put options in argv[]
                        mxGetString(prhs[3], cmd,  mxGetN(prhs[3]) + 1);
                        if((argv[argc] = strtok(cmd, " ")) != NULL)
                                while((argv[++argc] = strtok(NULL, " ")) != NULL)
                                        ;

                        for(i=1;i<argc;i++)
                        {
                                if(argv[i][0] != '-') break;
                                if((++i>=argc) && argv[i-1][1] != 'q')
                                {
                                        exit_with_help();
                                        fake_answer(plhs);
                                        return;
                                }
                                switch(argv[i-1][1])
                                {
                                        case 'b':
                                                prob_estimate_flag = atoi(argv[i]);
                                                break;
                                        case 'q':
                                                i--;
                                                info = &print_null;
                                                break;
                                        default:
                                                mexPrintf("Unknown option: -%c\n", argv[i-1][1]);
                                                exit_with_help();
                                                fake_answer(plhs);
                                                return;
                                }
                        }
                }

                model = matlab_matrix_to_model(prhs[2], &error_msg);
                if (model == NULL)
                {
                        mexPrintf("Error: can't read model: %s\n", error_msg);
                        fake_answer(plhs);
                        return;
                }

                if(prob_estimate_flag)
                {
                        if(svm_check_probability_model(model)==0)
                        {
                                mexPrintf("Model does not support probabiliy estimates\n");
                                fake_answer(plhs);
                                svm_free_and_destroy_model(&model);
                                return;
                        }
                }
                else
                {
                        if(svm_check_probability_model(model)!=0)
                                info("Model supports probability estimates, but disabled in predicton.\n");
                }

                predict(plhs, prhs, model, prob_estimate_flag);
                // destroy model
                svm_free_and_destroy_model(&model);
        }
        else
        {
                mexPrintf("model file should be a struct array\n");
                fake_answer(plhs);
        }

        return;
}