abcranger.cpp

#include "readstatobs.hpp"
#include "readreftable.hpp"
#include "cxxopts.hpp"
#include "EstimParam.hpp"
#include "ModelChoice.hpp"

int main(const int argc, const char* argv[]) {
    size_t nref, chosenscen;
    std::string headerfile,reftablefile,statobsfile;

    try {
        size_t nref;
        cxxopts::Options options(argv[0], " - ABC Random Forest - Model choice or parameter estimation command line options");

        options
            .positional_help("[optional args]")
            .show_positional_help();

        options.add_options()
            ("h,header","Header file",cxxopts::value<std::string>()->default_value("headerRF.txt"))
            ("r,reftable","Reftable file",cxxopts::value<std::string>()->default_value("reftableRF.bin"))
            ("b,statobs","Statobs file",cxxopts::value<std::string>()->default_value("statobsRF.txt"))
            ("o,output","Prefix output (modelchoice_out or estimparam_out by default)",cxxopts::value<std::string>())
            ("n,nref","Number of samples, 0 means all",cxxopts::value<size_t>()->default_value("0"))
            ("m,minnodesize","Minimal node size. 0 means 1 for classification or 5 for regression",cxxopts::value<size_t>()->default_value("0"))
            ("t,ntree","Number of trees",cxxopts::value<size_t>()->default_value("500"))
            ("j,threads","Number of threads, 0 means all",cxxopts::value<size_t>()->default_value("0"))
            ("s,seed","Seed, generated by default",cxxopts::value<size_t>()->default_value("0"))
            ("c,noisecolumns","Number of noise columns",cxxopts::value<size_t>()->default_value("5"))
            ("nolinear","Disable LDA for model choice or PLS for parameter estimation")
            ("plsmaxvar","Percentage of maximum explained Y-variance for retaining pls axis",cxxopts::value<double>()->default_value("0.9"))
            ("chosenscen","Chosen scenario (mandatory for parameter estimation)", cxxopts::value<size_t>())
            ("noob","number of oob testing samples (mandatory for parameter estimation)",cxxopts::value<size_t>())
            ("allpost","calculate all posteriors per model not just the selected one")
            ("parameter","name of the parameter of interest (mandatory for parameter estimation)",cxxopts::value<std::string>())
            ("g,groups","Groups of models",cxxopts::value<std::string>())
            ("save","save forest in ranger format")
            ("help", "Print help")
            ;
        auto opts = options.parse(argc,argv);

        if (opts.count("help")) {
          std::cout << options.help({"", "Group"}) << std::endl;
            exit(0);
        }

        headerfile = opts["h"].as<std::string>();
        reftablefile = opts["r"].as<std::string>();
        statobsfile = opts["b"].as<std::string>();
        nref = opts["n"].as<size_t>();

        if (opts.count("chosenscen") != 0 ||
            opts.count("parameter") != 0 ||
            opts.count("noob") != 0)
        {
            std::cout << "> Parameter Estimation <" << std::endl;
            if (opts.count("chosenscen") == 0 ||
                opts.count("parameter") == 0 ||
                opts.count("noob") == 0)
            {
                std::cout << "Error : please provide noob, parameter AND chosenscen arguments for parameter estimation." << std::endl;
                exit(1);
            }

            chosenscen = static_cast<double>(opts["chosenscen"].as<size_t>());

            auto myread = readreftable_scen(headerfile, reftablefile, chosenscen, nref);
            auto origobs = readStatObs(statobsfile);
            size_t nstat = myread.stats_names.size();
            size_t num_samples = origobs.size() / nstat;
            if (((origobs.size() % nstat) != 0) || (num_samples < 1)) {
                std::cout << "wrong number of summary statistics in statobs file." << std::endl;
                exit(1);           
            }       
            MatrixXd statobs = Map<MatrixXd>(origobs.data(),nstat,num_samples).transpose();
            auto res = EstimParam_fun(myread, statobs, opts);
        } else {
            std::cout << "> Model Choice <" << std::endl;
            
            auto myread = readreftable(headerfile, reftablefile, nref, false, opts.count("g") > 0 ? opts["g"].as<std::string>() : "");
            auto origobs = readStatObs(statobsfile);
            size_t nstat = myread.stats_names.size();
            size_t num_samples = origobs.size() / nstat;
            if (((origobs.size() % nstat) != 0) || (num_samples < 1)) {
                std::cout << "wrong number of summary statistics in statobs file." << std::endl;
                exit(1);           
            }       
            MatrixXd statobs = Map<MatrixXd>(origobs.data(),nstat,num_samples).transpose();
            auto res = ModelChoice_fun(myread,statobs,opts);

        }



    } catch (const cxxopts::exceptions::exception& e)
      {
        std::cout << "error parsing options: " << e.what() << std::endl;
        exit(1);
    }
}