src/SVMhierarchy.cpp

#include <fstream>
#include <iostream>
#include <iterator>
#include <vector>
#include <queue>

/*****************************************************************************
* This file is provided under the Creative Commons Attribution 3.0 license.
*
* You are free to share, copy, distribute, transmit, or adapt this work
* PROVIDED THAT you attribute the work to the authors listed below.
* For more information, please see the following web page:
* http://creativecommons.org/licenses/by/3.0/
*
* This file is a component of the Sleipnir library for functional genomics,
* authored by:
* Curtis Huttenhower (chuttenh@princeton.edu)
* Mark Schroeder
* Maria D. Chikina
* Olga G. Troyanskaya (ogt@princeton.edu, primary contact)
*
* If you use this library, the included executable tools, or any related
* code in your work, please cite the following publication:
* Curtis Huttenhower, Mark Schroeder, Maria D. Chikina, and
* Olga G. Troyanskaya.
* "The Sleipnir library for computational functional genomics"
*****************************************************************************/
#include "stdafx.h"
#include "cmdline.h"
#include "statistics.h"

using namespace SVMArc;
//#include "../../extlib/svm_light/svm_light/kernel.h"







int main(int iArgs, char** aszArgs) {
    gengetopt_args_info sArgs;

    CPCL PCL;
    SVMArc::CSVMSTRUCTTREE SVM;

    size_t i, j, k , iGene, jGene;
    double bestscore;
;
    ifstream ifsm;
    if (cmdline_parser(iArgs, aszArgs, &sArgs)) {
        cmdline_parser_print_help();
        return 1;
    }

    //Set Parameters
    SVM.SetLearningAlgorithm(sArgs.learning_algorithm_arg);
    SVM.SetVerbosity(sArgs.verbosity_arg);
    SVM.SetLossFunction(sArgs.loss_function_arg);

        cerr << "SetLossFunction" <<sArgs.loss_function_arg<< endl;

    if (sArgs.cross_validation_arg < 1){
        cerr << "cross_valid is <1. Must be set at least 1" << endl;
        return 1;
    }
    else if(sArgs.cross_validation_arg < 2){
        cerr << "cross_valid is set to 1. No cross validation holdouts will be run." << endl;
    }

    SVM.SetTradeoff(sArgs.tradeoff_arg);

    if (sArgs.slack_flag)
        SVM.UseSlackRescaling();
    else
        SVM.UseMarginRescaling();

    cerr << "SetRescaling" << endl;
    SVM.ReadOntology(sArgs.ontoparam_arg); // Read Ontology File
    if (!SVM.parms_check()) {
        cerr << "Parameter check not passed, see above errors" << endl;
        return 1;
    }

            cerr << "Parameter check" << endl;

    //  cout << "there are " << vecLabels.size() << " labels processed" << endl;
    size_t iFile;
    vector<string> PCLs;
    if (sArgs.input_given) {
        if (!PCL.Open(sArgs.input_arg, sArgs.skip_arg, sArgs.mmap_flag)) {
            cerr << "Could not open input PCL" << endl;
            return 1;
        }
    }

    //Read labels from file
    vector<SVMArc::SVMLabel> vecLabels;
    set<string> setLabeledGenes;
    if (sArgs.labels_given) {
        ifsm.clear();
        ifsm.open(sArgs.labels_arg);
        if (ifsm.is_open())
            vecLabels = SVM.ReadLabels(ifsm);
        else {
            cerr << "Could not read label file" << endl;
            return 1;
        }
        for (i = 0; i < vecLabels.size(); i++)
            setLabeledGenes.insert(vecLabels[i].GeneName);
    }
                cerr << "Read labels from file" << endl;


    //Training
    SAMPLE* pTrainSample;
    vector<SVMArc::SVMLabel> pTrainVector[sArgs.cross_validation_arg];
    vector<SVMArc::SVMLabel> pTestVector[sArgs.cross_validation_arg];
    vector<SVMArc::Result> AllResults;
    vector<SVMArc::Result> tmpAllResults;

    if (sArgs.model_given && sArgs.labels_given) { //learn once and write to file
        pTrainSample = SVM.CreateSample(PCL, vecLabels);
        SVM.Learn(*pTrainSample);
        SVM.WriteModel(sArgs.model_arg);
    } else if (sArgs.model_given && sArgs.output_given) { //read model and classify all
        vector<SVMLabel> vecAllLabels;

        for (size_t i = 0; i < PCL.GetGenes(); i++)
            vecAllLabels.push_back(SVMLabel(PCL.GetGene(i), 0));

        SVM.ReadModel(sArgs.model_arg);
        AllResults = SVM.Classify(PCL, vecAllLabels);
        ofstream ofsm;
        ofsm.open(sArgs.output_arg);
        if (ofsm.is_open())
            SVM.PrintResults(AllResults, ofsm);
        else {
            cerr << "Could not open output file" << endl;
        }
    } else if (sArgs.output_given && sArgs.labels_given) {
        //do learning and classifying with cross validation
        //set up training data
        if( sArgs.cross_validation_arg > 1){        
            for (i = 0; i < sArgs.cross_validation_arg; i++) {
                pTestVector[i].reserve((size_t) vecLabels.size()
                    / sArgs.cross_validation_arg + sArgs.cross_validation_arg);
                pTrainVector[i].reserve((size_t) vecLabels.size()
                    / (sArgs.cross_validation_arg)
                    * (sArgs.cross_validation_arg - 1)
                    + sArgs.cross_validation_arg);
                for (j = 0; j < vecLabels.size(); j++) {
                    if (j % sArgs.cross_validation_arg == i) {
                        pTestVector[i].push_back(vecLabels[j]);
                    } else {
                        pTrainVector[i].push_back((vecLabels[j]));
                    }
                }
            }
        }
        else{ // if you have less than 2 fold cross, no cross validation is done, all train genes are used and predicted
            // no holdout so train is the same as test gene set
            pTestVector[0].reserve((size_t) vecLabels.size() + sArgs.cross_validation_arg);
            pTrainVector[0].reserve((size_t) vecLabels.size() + sArgs.cross_validation_arg);

            for (j = 0; j < vecLabels.size(); j++) {
                pTestVector[0].push_back(vecLabels[j]);           
                pTrainVector[0].push_back(vecLabels[j]);            
            }

        }
        //set up training data done

        //set up validation data
        vector<SVMLabel> vec_allUnlabeledLabels;
        vector<Result> vec_allUnlabeledResults;
        vector<Result> vec_tmpUnlabeledResults;
        if (sArgs.all_flag) {
            vec_allUnlabeledLabels.reserve(PCL.GetGenes());
            vec_allUnlabeledResults.reserve(PCL.GetGenes());
            for (i = 0; i < PCL.GetGenes(); i++) {
                if (setLabeledGenes.find(PCL.GetGene(i))
                    == setLabeledGenes.end()) {
                        vec_allUnlabeledLabels.push_back(
                            SVMLabel(PCL.GetGene(i), 0));
                        vec_allUnlabeledResults.push_back(Result(PCL.GetGene(i)));
                }
            }
        }
        //run once
        for (i = 0; i < sArgs.cross_validation_arg; i++) {
            pTrainSample = SVM.CreateSample(PCL,
                pTrainVector[i]);

            cerr << "Cross Validation Trial " << i << endl;
            SVM.Learn(*pTrainSample);
            cerr << "Learned" << endl;
            tmpAllResults = SVM.Classify(PCL,   pTestVector[i]);
            cerr << "Classified " << tmpAllResults.size() << " examples"<< endl;
            AllResults.insert(AllResults.end(), tmpAllResults.begin(), tmpAllResults.end());
            tmpAllResults.resize(0);
            if (sArgs.all_flag) {
                vec_tmpUnlabeledResults = SVM.Classify(
                    PCL, vec_allUnlabeledLabels);
                
                if(i == 0){
                    for (j = 0; j < vec_tmpUnlabeledResults.size(); j++){
                        vec_allUnlabeledResults[j].num_class = vec_tmpUnlabeledResults[j].num_class;
                        for( k = 1; k <= vec_tmpUnlabeledResults[j].num_class; k++)
                            vec_allUnlabeledResults[j].Scores.push_back(vec_tmpUnlabeledResults[j].Scores[k]);
                    }
                }
                else{
                    for (j = 0; j < vec_tmpUnlabeledResults.size(); j++)
                        for( k = 1; k <= vec_tmpUnlabeledResults[j].num_class; k++)
                            vec_allUnlabeledResults[j].Scores[k] += vec_tmpUnlabeledResults[j].Scores[k];
                }

            }
            if (i > 0) {
                SVMArc::CSVMSTRUCTTREE::FreeSample(*pTrainSample);
            }
        }
                    cerr << "5" << endl;

        if (sArgs.all_flag) { //add the unlabeled results
            for (j = 0; j < vec_allUnlabeledResults.size(); j++)
                for( k = 1; k <= vec_allUnlabeledResults[j].num_class; k++){
                    if(k==1){
                        vec_allUnlabeledResults[j].Scores[k]/= sArgs.cross_validation_arg;
                        bestscore=vec_allUnlabeledResults[j].Scores[k];
                        vec_allUnlabeledResults[j].Value=k;
                    }else{
                        vec_allUnlabeledResults[j].Scores[k]/= sArgs.cross_validation_arg;
                        if(vec_allUnlabeledResults[j].Scores[k] < bestscore){
                            bestscore = vec_allUnlabeledResults[j].Scores[k];
                            vec_allUnlabeledResults[j].Value=k;
                        }
                    }
                }

            AllResults.insert(AllResults.end(),
                vec_allUnlabeledResults.begin(),
                vec_allUnlabeledResults.end());
                        cerr << "6" << endl;

        }
                    cerr << "7" << endl;

        ofstream ofsm;
        ofsm.clear();
        ofsm.open(sArgs.output_arg);
        SVM.PrintResults(AllResults, ofsm);
        return 0;

    } else {
        cerr << "More options are needed" << endl;
    }

}