From 6fd3f70ca2ce5460fc3c6eab138c6a25b7bcfc69 Mon Sep 17 00:00:00 2001
From: John Denker <jsd@av8n.com>
Date: Mon, 18 Oct 2021 17:01:35 -0700
Subject: smarter about metadata; get rid of unsanitary magic numbers (u5)

---
 poiss.c | 41 +++++++++++++++++++++++++----------------
 1 file changed, 25 insertions(+), 16 deletions(-)

diff --git a/poiss.c b/poiss.c
index a286184..dace05e 100644
--- a/poiss.c
+++ b/poiss.c
@@ -8,6 +8,9 @@ void usage() {
 Uses actual Poisson probability (i.e. not Gaussians, and
 therefore not least squares).
 
+Typical usage:
+:; ./poiss -if my_data.csv -fit cu.csv
+
 Options include:
   -h                    # print this message (and immediately exit)
   -ifile   $fn          # input file, date for -fit
@@ -38,7 +41,9 @@ See ./doit for a working example, including the SS --meld stanzas.
 #include "poissonier.h"
 
 // Used for equal weighting:
-vector<double> const u5(5, 1.0);
+vector<double> uvec(int const ddd) {
+  return vector<double>(ddd, 1.0);
+}
 
 // Class to describe a sitch, i.e. a situation, i.e.
 // everything the objective function needs to know:
@@ -257,12 +262,14 @@ struct guesser {
           << endl;
 
     vector<double> slow_model{slow_r, slow_a, bl_r};
+    int nprm = slow_model.size();
 
     double fast_t0 = evt[0].t;
     double fast_n0 = evt[0].n;
     double fast_t1 = evt[fast_i].t;
     double fast_n1 = evt[fast_i].n;
-    double model_inty  = inty_multi_exp(slow_model, fast_t0, fast_t1, u5);
+    double model_inty = inty_multi_exp(slow_model,
+                 fast_t0, fast_t1, uvec(nprm));
     double fast_tot = fast_n1 - fast_n0;
     double fast_dn = fast_tot - model_inty;
     double fast_r = fast_dn * fast_a;
@@ -296,21 +303,20 @@ struct guesser {
     }
     ostream& xout(ofile != "-" ? xxout : cout);
 
-    size_t ss = sitch.t1.size();
-    xout << ss << endl;
-    xout << sitch.t1[ss-1]
-          << "  " << sitch.t2[ss-1]
-          << "  " << sitch.obs[ss-1]
+    size_t nprm = sitch.t1.size();
+    xout << nprm << endl;
+    xout << sitch.t1[nprm-1]
+          << "  " << sitch.t2[nprm-1]
+          << "  " << sitch.obs[nprm-1]
           << endl;
 
-    vector<double> prm(u5);
+    vector<double> prm(uvec(nprm));
 
     xout << dump(model) << endl;
 
     xout << "bl/norm,bl,limp" << endl;
     double bl0 = prm[4];
     void* context = (void*)(&sitch);
-    int nprm(prm.size());
 
     for (double bl = bl0*.9; bl <= bl0*1.1; bl += bl/100.) {
       prm[4] = bl;
@@ -377,9 +383,9 @@ void do_fit(poissonier& poi, string const ifile,
   ggg.more();
   ggg.sitch.norm = ggg.model;
   ggg.show(gfile);
+  int nprm = ggg.model.size();
 
-  vector<double> prm = u5;
-  int nprm(prm.size());
+  vector<double> prm = uvec(nprm);
   vector<double> const lower(nprm, 0.5);
   vector<double> const upper(nprm, 1.5);
 
@@ -402,6 +408,7 @@ void do_fit(poissonier& poi, string const ifile,
   }
 
   if (OK && ofile != "") {
+    string meta = "fast.amp, fast.dk, slow.amp, slow.dk, bl";
     ofstream xxout;
     if (ofile != "-") {
       xxout.open(ofile);
@@ -421,7 +428,7 @@ void do_fit(poissonier& poi, string const ifile,
           << ",   perdof:, " << limp_end / DoF
           << endl;
     xout << endl;
-    xout << ",fast.amp, fast.dk, slow.amp, slow.dk, bl" << endl;
+    xout << "," << meta << endl;
     xout << "Normed:, " << dump(found) << endl;
     vector<double> combi(nprm);
     for (int ii = 0; ii < nprm; ii++) {
@@ -437,12 +444,13 @@ void do_fit(poissonier& poi, string const ifile,
 // In particular, calculate the Mahalanobis metric
 // i.e. the second derivative (Hessian) of the log improbability.
     vector<vector<double> > covar(nprm, vector<double>(nprm));
-    vector<double> probe(u5);
+    vector<double> probe(uvec(nprm));
     sitcher sitch(ggg.sitch);
     sitch.norm = combi;
     void* context = (void*)(&sitch);
     double delta(0.01);
     xout << endl;
+    xout << meta << endl;
     for (int ii = 0; ii < nprm; ii++) {
       string sep = "";
       for (int jj = 0; jj < nprm; jj++) {
@@ -486,12 +494,13 @@ void do_dk(poissonier& poi, string const ofile) {
 
   double step(1);
   vector<double> prm( {30, .1, 10, .01, 1} );
+  int nprm(prm.size());
 
   int dt(5);
   for (int ii = 0; ii < 500; ii+= dt) {
     double tt1(step*ii);
     double tt2(tt1 + dt);
-    double intensity(inty_multi_exp(prm, tt1, tt2, u5));
+    double intensity(inty_multi_exp(prm, tt1, tt2, uvec(nprm)));
     int obs = poi.sample(intensity);
     sitch.t1.push_back(tt1);
     sitch.t2.push_back(tt2);
@@ -509,13 +518,13 @@ void do_dk(poissonier& poi, string const ofile) {
 // and the observations (to the extent they exist)
   for (int ii = 0; ii < 500; ii++) {
     double tt1(smooth.t1[ii]);
-    double sminty(inst_multi_exp(prm, tt1, u5));
+    double sminty(inst_multi_exp(prm, tt1, uvec(nprm)));
     xout << tt1 << "," << sminty;
     if (ii < (int) sitch.t1.size()) {
       double fake(max(0.+sitch.obs[ii], 0.2));   // to facilitate log axes
       double tt1(sitch.t1[ii]);
       double tt2(sitch.t2[ii]);
-      double minty(inty_multi_exp(prm, tt1, tt2, u5));
+      double minty(inty_multi_exp(prm, tt1, tt2, uvec(nprm)));
       xout << "," << tt1
            << "," << tt2
            << ",x"
-- 
cgit v1.2.3