path: root/svd.c

/////////////////////////

#include <iostream>
using namespace std;

void usage() {
  cout << R"EoF(Program to fit to perform Singular Value Decomposition (SVD).
Typical usage:
  ./svd -m -1 -ifile whatever.csv -svd eigen.csv

Options include:
  -h                    # print this message (and immediately exit)
  -ifile   $fn          # input data file
  -svd     $fn          # write the results to $fn
  -meta    NN           # use NNth line of input header as metadata
                        # NN=-1 means use *last* line of input header

To write to stdout, use "-" for the $fn.
)EoF";
}

#include <iomanip>      /* for setw */
#include <vector>
#include "arg_parser.h"
#include <string>
#include <sstream>
#include <map>
#include "parse_csv.h"
#include "poissonier.h"
#include <armadillo>


// Format all the elements of a C++ vector,
// creating a string suitable for printing.
template<class T>
string dump(vector<T> const arg) {
  stringstream rslt;
  string sep = "";
  for (T val : arg) {
    rslt << sep << val;
    sep = ", ";
  }
  return rslt.str();
}

// Format all the elements of a counted array,
// creating a string suitable for printing.
template<class T>
string dump(int const nn, T const arg[]) {
  stringstream rslt;
  string sep = "";
  for (int ii = 0; ii < nn; ii++) {
    rslt << sep << arg[ii];
    sep = ",";
  }
  return rslt.str();
}

// numerically accurate sinh(arg)/arg
template<class T>
inline T sinhc(T arg) {
  if (abs(arg) < 1e-4) return 1 + arg*arg/6;
  return sinh(arg)/arg;
}

struct evendor {
  double t;               // time
  double n;               // number of events
};

// Intensity, usually denoted lambda,
// i.e. integral of rate(t) dt over the given interval.
double inty_multi_exp(vector<double> const& prm,
        double const t1, double const t2, vector<double> const norm) {
  double rslt(0);
  int nn = prm.size();
  for (int ii = 0; ii < nn; ii+=2){
    double amp(prm[ii] * norm[ii]);
    double rate(1+ii < nn ? prm[1+ii] * norm[1+ii] : 0.);
    double mid((t2 + t1) / 2.);
    double tau(t2 - t1);
    rslt += amp * exp(-rate*mid) * tau * sinhc(rate*tau/2.);
  }
  return rslt;
}

// Add something to a particular component of a vector.
template<class T>
vector<T> goose(vector<T> const& vvv, int const ii, T const delta) {
  vector<T> rslt(vvv);
  rslt[ii] += delta;
  return rslt;
}

typedef vector<vector<string> > VVS;

void dump_matrix(ostream& ouch, arma::mat const &mmm,
        int const pr = 3, int const texmode=0){
  int wd(pr+6);
  int dim(mmm.n_rows);
  for (int ii = 0; ii < dim; ii++) {
    for (int jj = 0; jj < dim; jj++) {
      if (texmode) {
      if ( mmm(ii, jj) > 3 )
        ouch << "\\y{"; /* } */
        else ouch << "\\x{"; /* } */
        ouch  << setprecision(pr) << setw(wd) << mmm(ii, jj)
                /* { */ << "} ";
      } else {
        ouch << setprecision(pr) << setw(wd) << mmm(ii, jj) << ", ";
      }
    }
    if (texmode) {
      ouch << " \\z" << endl;
    } else {
      ouch << endl;
    }
  }
}

string dump(arma::mat const mmm, string const prefix="") {
  int NR = mmm.n_rows;
  int NC = mmm.n_cols;
  stringstream rslt;
  for (int ii = 0; ii < NR; ii++) {
    string sep = prefix;
    for (int jj = 0; jj < NC; jj++) {
      rslt << sep << fixed << setw(11) << mmm(ii, jj);
      sep = ", ";
    }
    rslt << endl;
  }
  return rslt.str();
}

string dump(arma::vec vvv, arma::mat const mmm) {
  int NR = mmm.n_rows;
  int NC = mmm.n_cols;
  stringstream rslt;
  for (int ii = 0; ii < NR; ii++) {
    rslt << fixed << setw(11) << vvv(ii)
         << ",";              // skip a column
    for (int jj = 0; jj < NC; jj++) {
      rslt << "," << fixed << setw(11) << mmm(ii, jj);
    }
    rslt << endl;
  }
  return rslt.str();
}

string dumpx(arma::vec vvv, arma::mat const mmm) {
  int NR = mmm.n_rows;
  int NC = mmm.n_cols;
  stringstream rslt;
  for (int ii = 0; ii < NR; ii++) {
    double tmp = vvv(ii);
    rslt << fixed << setw(11) << tmp;
    rslt << ",";              // skip a column
    for (int jj = 0; jj < NC; jj++) {
      rslt << "," << fixed << setw(11) << mmm(ii, jj);
    }
    rslt << endl;
  }
  return rslt.str();
}

void do_svd(string const ifile, string const ofile,
        int const gotmeta, int whichmeta){
  if (ifile == "") {
    if (ofile != "") throw invalid_argument("do_fit needs an input file");
    return;
  }

  ifstream inx;
  if (ifile != "-") {
    inx.open(ifile);
    if (! inx.good()) {
      cerr << "Cannot open input '" << ifile << "'" << endl;
      exit(1);
    }
  }

  istream& in(ifile == "-" ? cin : inx);
  VVS alldata;
  alldata = readCSV<string>(in);
  cout << alldata.size() << endl;
  int hh = count_header(alldata);
  cout << hh << endl;
  VVS aoa(alldata.begin()+hh, alldata.end());
  int NR = aoa.size();
  if (NR == 0) throw invalid_argument("no data???");
  int NC = aoa[0].size();
  cout << "NR: " << NR
       << "  NC: " << NC << endl;
  arma::mat maha(NR, NC, arma::fill::zeros);
  for (int ii = 0; ii < NR; ii++) {
    for (int jj = 0; jj < NC; jj++) {
      maha(ii, jj) = stod(aoa[ii][jj]);
    }
  }

  if (ofile != "") {
    ofstream xxout;
    if (ofile != "-") {
      xxout.open(ofile);
    }
    ostream& xout(ofile != "-" ? xxout : cout);

    string meta;
    if (gotmeta) {
      if (whichmeta < 0) whichmeta += hh;
      meta = dump(alldata[whichmeta]);
    }

    xout << ",," << meta << endl;
    xout << dump(maha, ",,");   // already includes newlines

    arma::mat U, V;
    arma::vec eival;
    arma::svd(U, eival, V, maha);

    xout << endl;
    xout << "cost²,," << meta << endl;
    xout << dumpx(eival, V) << endl;
  }
}

int main(int argc, char * const * argv) {
  poissonier poi;
  string seed;
  string ifile;
  arg_parser arg(argc, argv);
  arg.fold_case = 1;
  string progname = arg.nextRaw();
  string svdfile;
  int gotmeta(0);
  int whichmeta(0);

  for (;;){
    string raw_arg = arg.nextRaw();         // get next keyword
    if (arg.fail()) break;
    int where = arg.cur.length();
    if (where){
      if (arg.cur[where-1] == ':') arg.cur = arg.cur.substr(0,where-1);
    }
    if (0) {}
    else if (arg.prefix("help") || arg.prefix("-h")) {
      usage();
      exit(0);
    }
    else if (arg.prefix("-meta")) {
      gotmeta++;
      arg >> whichmeta;
    }
    else if (arg.prefix("-ifile")) arg >> ifile;
    else if (arg.prefix("-seed")) arg >> seed;
    else if (arg.prefix("-svd")) arg >> svdfile;
    else {
      cerr << "Unrecognized command '" << raw_arg << "'" << endl;
      exit(1);
    }
  }
  if (seed == "--") seed = poi.randomize();
  else if (seed != "") poi.seed(seed);

  if (!(svdfile.length()
     ))
         svdfile = "-";
  do_svd(ifile, svdfile, gotmeta, whichmeta);
}