add a TUNING.txt guide for postFringeVx.cc. Add validation steps for...

add a TUNING.txt guide for postFringeVx.cc. Add validation steps for TransformFitData and a progressive determination of the best polynomial fit of the phases

AnaFringes/TUNING.txt

+Some GuideLines to tune postFringe
+---------------------------------
+JEC: 7 Juin 18
+
+o Pour postFringe v3
+Running fitphi.pic selon
+(cd CasA; spiapp -term -exec fitphi.pic cor_CasA30mai18b ; stty sane)
+Le plot des residus des fits permets d'investiguer des problemes
+dans la plus part des cas venant de la procedure qui met bout a bout
+les segments de phase defilant entre -Pi et Pi.
+
+La routine incriminee est TransformFitData qui drecoit en entree 
+une collection (GeneralFitData) de points (freq_i, val_i) qui sont
+des medianes sur un certain nombres de points d'origine 
+(voir estimePhase les vecteurs x et y issus de medianReduct)
+
+
+Donc dans TransformFitData: 
+1) on cherche les points de transition -Pi->Pi (cas actuel) 
+ou +Pi->-Pi par la methode des differences adjacentes, puis en 
+la seuillant par une valeur Min et Max (TUNING POSSIBLE) pour 
+reconnaitre des points pathologiques que l'on suprime
+
+2) avec le nouveau set de points cleaned on reprocede de nouveau
+a des differences adjacentes. Un nouveau seuillage (TUNING POSSIBLE)
+ permet de determiner les points de transitions (vecteur idy)
+
+3) on controle si il y a des pts de transition aberrants par exemple
+en controlant que la phase avant et apres ce point doit faire un saut
+suffisant (TUNING POSSIBLE)
+
+4) une fois les points de transition valides on procede 
+progressivement au calcul des nouvelles phases en "additionnant" -2pi 
+selon le signe de la pente de variation. 
+
+Une fois les nouveaux points (freq_i, phase_i) on procede au fit
+polynomial de la relation phase(freq). Cependant, on procede
+progressivement en prenant un polynome de degre {1,...., polDegreMax}
+ avec polDegreMax controle en ligne de commande par (-degre <d>).
+Le controle se fait en calculant l'inter-quartile normalise des
+valeurs absolues des residus. Le seuil qui accepte un polynome
+est iqrn_thr (TUNING POSSIBLE). 
+
+Si aucun polynome est satisfaisant alors l'objet PolyPhi<phase> n'est
+cree dans le ficher de controle 
+
+
+fOutName =  outputDir + "/" + "PhaseParam-" + source + "-" + Num2String(FreqMin) + "-" + Num2String(FreqMax) + "-" + Type +".txt";
+  
+
+

AnaFringes/poly.cc

+#include "machdefs.h"
+
+//---- System et stdc++ include files
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <string.h>
+#include <iostream>
+#include <fstream>
+#include <complex>
+
+#include <typeinfo>
+#include <string>
+#include <vector>
+#include <map>
+#include <functional>
+#include <list>
+
+//---- Sophya include files
+#include "sopnamsp.h"
+#include "basetools.h"
+#include "systools.h"
+#include "sutils.h"
+#include "ntools.h"
+#include "array.h"
+#include "histats.h"
+
+//---- Spiapp include files
+#include "nobjmgr.h"
+#include "servnobjm.h"
+
+//---- Include files from additionnal modules
+
+//---- function to compare bits on double 
+int_8 BitCmp64(double v,int_8 flg) 
+{return ((int_8)((v<0.) ? v-0.1 : v+0.1))&flg;} 
+
+//---- function for Adding and displaying Objects  
+void Keep_Object(AnyDataObj & obj, string const & nom) 
+{ 
+  string name = nom; 
+  NamedObjMgr om; 
+  if (om.GetObj(name)) 
+    cerr << "KeepObj()/Warning Already kept object " << endl; 
+  else om.AddObj(obj, name); 
+} 
+
+void Display_Object(AnyDataObj & obj, string const & opt, string const & nom) 
+ { 
+  string name = nom; 
+  NamedObjMgr om; 
+  if (!om.GetObj(name)) 
+    om.AddObj(obj, name); 
+  om.DisplayObj(name, opt); 
+} 
+
+//---- function for getting and setting ObjectManager variables  
+void Set_ObjMgrVar(MuTyV v, string const & nom) 
+{ 
+  NamedObjMgr om; 
+  om.SetVar(nom, (string)v); 
+} 
+
+MuTyV Get_ObjMgrVar(const char * nom) 
+{ 
+  string name = nom; NamedObjMgr om; 
+  MuTyV v = om.GetVar(name); 
+  return v; 
+} 
+
+//---- Macro for Objects and variables saving
+#define KeepObj(obj) Keep_Object(obj, #obj)
+#define GetOMVar(var) Get_ObjMgrVar( #var )
+#define SetOMVar(var) Set_ObjMgrVar(var, #var )
+
+//---- Macro Displaying objects and command execution
+#define DispObj(obj, att) Display_Object(obj, att, #obj); 
+
+#define ExecCmd(cmd) srvo.ExecuteCommand(cmd); 
+
+
+
+//-------------------------------------------------//
+//----------------- User Functions ----------------//
+//-------------------------------------------------//
+
+extern "C" {
+  int Values( vector<string>& args );
+}
+
+int Values( vector<string>& args )
+{
+// Some definitions to help using spiapp;
+NamedObjMgr omg;
+Services2NObjMgr& srvo = *omg.GetServiceObj();
+
+//-------------- Object List --------------
+//Number of objects = 16
+string ___nomobj;
+
+___nomobj = "PolyPhi2H";
+Poly * ___PolyPhi2H = dynamic_cast< Poly  * >(omg.GetObj(___nomobj));
+if(___PolyPhi2H==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object PolyPhi2H... please update file");
+Poly & PolyPhi2H = (*___PolyPhi2H);
+
+___nomobj = "PolyPhi3H";
+Poly * ___PolyPhi3H = dynamic_cast< Poly  * >(omg.GetObj(___nomobj));
+if(___PolyPhi3H==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object PolyPhi3H... please update file");
+Poly & PolyPhi3H = (*___PolyPhi3H);
+
+___nomobj = "PolyPhi4H";
+Poly * ___PolyPhi4H = dynamic_cast< Poly  * >(omg.GetObj(___nomobj));
+if(___PolyPhi4H==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object PolyPhi4H... please update file");
+Poly & PolyPhi4H = (*___PolyPhi4H);
+
+
+
+//--------------------------------------------//
+//----------------- User Code ----------------//
+//--------------------------------------------//
+
+ r_8 freqMin= 1300.;
+ r_8 freqMax = 1450.;
+ int Nfreq = 1200;
+ r_8 dFreq = (freqMax-freqMin)/((r_8)Nfreq);
+
+ r_8 freqRef = 1375.091552734375;
+ 
+ TVector<r_8> freq(Nfreq);
+ freq = RegularSequence(freqMin,dFreq);
+ GeneralFitData p2H(1,Nfreq);
+ GeneralFitData p3H(1,Nfreq);
+ GeneralFitData p4H(1,Nfreq); 
+ 
+ for(int i=0;i<Nfreq;i++){
+   r_8 f=freq(i);
+   p2H.AddData1(f,PolyPhi2H(f-freqRef));
+   p3H.AddData1(f,PolyPhi3H(f-freqRef));
+   p4H.AddData1(f,PolyPhi4H(f-freqRef));
+ }
+
+ KeepObj(p2H);
+ KeepObj(p3H);
+ KeepObj(p4H);
+
+return 0;
+}

AnaFringes/postFringeV3.cc

@@ -138,9 +138,16 @@ std::string Num2String(T val){
 //---- MEDIAN : warning the nth_element reorder the input vector => input NO MORE VALID
 //------------------------
 
+//Warning; original order is not preserved
+class median_of_empty_list_exception:public std::exception{
+    virtual const char* what() const throw() {
+    return "Attempt to take the median of an empty list of numbers.  "
+      "The median of an empty list is undefined.";
+  }
+};
 template<class RandAccessIter>
 double median(RandAccessIter begin, RandAccessIter end) {
-  
+  if(begin == end){ throw median_of_empty_list_exception(); }  
   std::size_t size = end - begin;
   std::size_t middleIdx = size/2;
   RandAccessIter target = begin + middleIdx;
@@ -169,6 +176,39 @@ T median(TVector<T>& vec){
 }
 
 
+//interquartile range normalized to Gauss(0,1) => sigma estimator
+//The sigma on  iQRNorm is 1.16639*iQRNorm/sqrt(n) in case of gaussian random variable.
+//warning median does not preserve original order
+template<class RandAccessIter>
+double iQRNorm(RandAccessIter begin, RandAccessIter end) {
+  if(begin == end){ throw median_of_empty_list_exception(); }
+
+  vector<r_4> tmp1(begin,end);
+  sort(tmp1.begin(),tmp1.end());
+  vector<r_4> tmp3(tmp1);
+  std::size_t size = tmp1.end() - tmp1.begin();
+  std::size_t middleIdx = size/2;
+  double q1 = median(tmp1.begin(), tmp1.begin()+middleIdx);
+  double q3 = median(tmp3.begin()+middleIdx,tmp3.end());
+  return (q3-q1)/1.34898; //normalisation by Q3-Q1 for a Gauss(0,1) value
+}
+
+
+template <class T>
+T iQRNorm(TVector<T>& vec){
+  vector<T>tmp(vec.ConvertTostdvec());
+  return iQRNorm(tmp.begin(),tmp.end());
+}
+
+template <class T>
+T iQRNorm(const vector<T>& vec){
+  vector<T>tmp(vec);
+  return iQRNorm(tmp.begin(),tmp.end());
+}
+
+
+
+
 //input : vecin
 //input : nbre de bins sur lesquels est calcule la mediane des vecin[i]
 //output: vecout
@@ -209,7 +249,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
   std::adjacent_difference(y.begin(),y.end(),dy.begin());
   
   
-
   //  r_8 minval = *min_element(dy.begin(),dy.end(),abscomp);
   r_8 minval = median(dy);
   r_8 maxval = *max_element(dy.begin(),dy.end(),abscomp);
@@ -224,13 +263,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
   vector<bool>nonvalidPts(x.size(),false);
   for(int i=1;i<dy.size()-1;i++){
     if(fabs(dy[i])>minval*1.5 && fabs(dy[i])<0.99*maxval){
-//       cout << "Patho: i,x,y,dy[i-1],dy[i],dy[i+1]: " 
-// 	   << i << ", " 
-// 	   << x[i] << ", "  
-// 	   << y[i] << ", "  
-// 	   << dy[i-1] << ", "  
-// 	   << dy[i] << ", "  
-// 	   << dy[i+1] << "\n";  
       nonvalidPts[i]=true;
       i++;
     }
@@ -239,7 +271,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
   //clean x,y,ey vectors: proceed backward to preserve iterator
   for(int i=nonvalidPts.size()-1;i>=0;i--){
     if(nonvalidPts[i]){
-      cout << "erasing pt["<<i<<"]"<<endl;
       x.erase(x.begin()+i);
       y.erase(y.begin()+i);
       ey.erase(ey.begin()+i);
@@ -247,17 +278,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
   }
 
 
-  {
-    string outputDir = param.outputDir;
-    string fNaneDbg =  outputDir + "/redpts-" + tag + ".txt";
-    ofstream ofsDbg(fNaneDbg.c_str());
-    for(int ib=0;ib<x.size();ib++){
-      ofsDbg << x[ib] << " " << y[ib] << endl;
-    }
-    ofsDbg.close();
-  }
-
-
 
   //repeat adjacent difference with the cleaned data set
   TVector<r_8> ty(y); //for convenience use TVector
@@ -282,29 +302,59 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
 
   vector<sa_size_t> idy;
   for(sa_size_t i=0;i<tdy.Size(); i++){
-    if(fabs(tdy(i))>0.5*maxval) idy.push_back(i);
+    if(fabs(tdy(i))>0.25*maxval) idy.push_back(i);
   }
   if(idy.empty())
     throw SzMismatchError("TransformFitData: FATAL: idy vide");
 
-  cout<< "idy <";
+
+  //Confirmation of transition points JEC 7/6/18
+  //
+  cout<< "BEFORE Confirmation idy <";
   for(int i=0;i<idy.size();i++)cout<<idy[i]<< " ";
   cout << ">" << endl;
+
   
+  for(int i=idy.size()-1;i>=0;i--){ //go backward to keep iterator valid
+    int j = idy[i];
+    if(j == 0 || j == y.size()-1) { 
+      //the transition at the edge: not valid
+      cout << "Transition at edge: not validated" << endl;
+      idy.erase(idy.begin() + i);
+    } else { 
+      //Try to see if points before and after a transition point have a difference of phase close to 2pi 
+      r_8 diffPhase = fabs(y[j+1]-y[j-1]);
+      if(diffPhase < 0.5*(2*M_PI)){
+	cout << "Transition not validated: diff phase: " << diffPhase  << endl;
+	idy.erase(idy.begin() + i);
+
+      }
+    }
+  }
+
+
+  cout<< "After Confirmation idy <";
+  for(int i=0;i<idy.size();i++)cout<<idy[i]<< " ";
+  cout << ">" << endl;
+
+
+  {
+    string outputDir = param.outputDir;
+    string fNaneDbg =  outputDir + "/redpts-" + tag + ".txt";
+    ofstream ofsDbg(fNaneDbg.c_str());
+    for(int ib=0;ib<x.size();ib++){
+      ofsDbg << x[ib] << " " << y[ib] << endl;
+    }
+    ofsDbg.close();
+  }
+
   //Progressively build the one-piece curve
   for(int i=0;i<idy.size();i++){    
     for(sa_size_t j=0;j<idy[i];j++){
       ty(j) -= sigYslope*2*M_PI;
     }
   }
-  
-//   //Find a number of "2pi" to translate the curve and get a rather symmetrical curve
-//   r_8 meany = ty.Sum()/ty.Size();
-//   r_8 reste =  remainder(meany,2*M_PI);
-//   int n2pi = int((meany-reste)/(2*M_PI));
-//   ty -= n2pi*2*M_PI;
-  
-  
+    
   //Modify original data
   //Reset the structure to fill again
   data.SetDataPtr();
@@ -685,15 +735,15 @@ void estimePhase(){
     
 
 
-    vector<r_8> freqToTest;
-    vector<r_8> idxsigValid;
-    for(size_t i=0; i<dsigny.size(); i++){
-      //      if(copysign(1.0,medslope)*dsigny[i]>1.5){
-      freqToTest.push_back(freq(i)); //JEC 31/5/18 use the reduced vector
-      //	freqToTest.push_back(x[i]);
-      idxsigValid.push_back(i);
-      //      }
-    }
+//     vector<r_8> freqToTest;
+//     vector<r_8> idxsigValid;
+//     for(size_t i=0; i<dsigny.size(); i++){
+//       //      if(copysign(1.0,medslope)*dsigny[i]>1.5){
+//       freqToTest.push_back(freq(i)); //JEC 31/5/18 use the reduced vector
+//       //	freqToTest.push_back(x[i]);
+//       idxsigValid.push_back(i);
+//       //      }
+//     }
 
 
     //----------
@@ -760,46 +810,78 @@ void estimePhase(){
       }
 
 
-      int polDegre = param.poldeg;
-      Poly pol(polDegre);
+      int polDegreMax = param.poldeg;
+      
+      r_8 iqrn = M_PI;
+      r_8 iqrn_thr=0.05;
+      
+      int polDegre=0;
       r_8 freqMid = median(x);// (FreqMax+FreqMin)*0.5;
-      r_8 rcFit = dataFitReduced.PolFit(0,pol,polDegre,true,freqMid);
-      if(rcFit<0){
-	cout << "PolFit Failed: rc= " << rcFit;
-	throw MathExc("PolFit Failed:");
-      }
-      cout << "PolFit Success: dump: \n";
-      pol.Print(cout);
-      cout << endl;
+      vector<r_8>abserr;
+      
+      while(polDegre<=polDegreMax) {
+	polDegre += 1;
+
+	Poly pol(polDegre);
+	r_8 rcFit = dataFitReduced.PolFit(0,pol,polDegre,true,freqMid);
+	if(rcFit<0){
+	  cout << "PolFit Failed: rc= " << rcFit;
+	  throw MathExc("PolFit Failed:");
+	}
+	cout << "PolFit Success degre= "<< polDegre << " dump: \n";
+	pol.Print(cout);
+	cout << endl;
+      
+	//compute the residuals
+	abserr.clear();
+	for(int i=0; i<dataFit.NData(); i++){
+	  r_8 xp[1]; xp[0]=dataFit.X1(i); 
+	  abserr.push_back(fabs(dataFit.Val(i) - remainder(pol(xp[0]-freqMid),2*M_PI)));
+	}
+	r_8 iqrn = iQRNorm(abserr);
+	cout << "IQRn of AbsError:= " << iqrn << endl;
+      
+	
+	if(iqrn<=iqrn_thr) { 
+	  //save polynomial class
+	  string nameTag = "Poly"+vPhaseName[ip];
+	  fOut << PPFNameTag(nameTag) << pol;
+      
 
-      ofs << setprecision(20)
- 	  <<  vPhaseName[ip] 
- 	  << " Polynomial "
-	  << polDegre << " "
-	  << freqMid << " " ;
-      for(int i=0;i<=polDegre;i++){
-	ofs <<pol.Element(i) << " ";
-      }
-      ofs << endl;
 
-      //compute the residuals
+	  ofs << setprecision(20)
+	      <<  vPhaseName[ip] 
+	      << " Polynomial "
+	      << polDegre << " "
+	      << freqMid << " " ;
+	  for(int i=0;i<=polDegre;i++){
+	    ofs <<pol.Element(i) << " ";
+	  }
+	  ofs << endl;
 
-      for(int i=0; i<dataFit.NData(); i++){
-	r_8 xp[1]; xp[0]=dataFit.X1(i); 
-	r_8 residual = dataFit.Val(i) - remainder(pol(xp[0]-freqMid),2*M_PI);
-	dataRes.AddData1(xp[0],residual);
-      }
 
+      
+	  //save the residuals
 
+	  for(int i=0; i<dataFit.NData(); i++){
+	    r_8 xp[1]; xp[0]=dataFit.X1(i); 
+	    r_8 residual = dataFit.Val(i) - remainder(pol(xp[0]-freqMid),2*M_PI);
+	    dataRes.AddData1(xp[0],residual);
+	  }
 
-//     }//good model
-    
-    //save residuals      
-      string nameTag = "Res"+vPhaseName[ip];
-      fOut << PPFNameTag(nameTag) << dataRes;
+	  //save residuals      
+	  nameTag = "Res"+vPhaseName[ip];
+	  fOut << PPFNameTag(nameTag) << dataRes;
 
-  }//loop on Phase
+	  //jump after the while
+	  break;
+
+	}//ok polynomial
+      }//while polynomial degree
 
+      //     }//good model      
+  }//loop on Phase
+  
   ofs.close();
 
 
@@ -1164,7 +1246,7 @@ int main(int narg, char* arg[]) {
   string Type="RE"; //RE or IM
   
   //JEC 5/6/18 add polynomial degree
-  int poldeg = 9; //degree of the polynomial description of Phase in the saw-tooth case
+  int poldeg = 15; //degree of the polynomial description of Phase in the saw-tooth case
   
 
   int ka=1;
@@ -1183,7 +1265,7 @@ int main(int narg, char* arg[]) {
       cout << "-f: freq. range MHz freqMin,freqMax  (ex. 1300,1450)\n";
       cout << "-src: source name (ex. CasA18oct17)\n";
       cout << "-type: type de cross REal or IMag [RE]\n";
-      cout << "-deg: degree of the polynomial description of Phase in the saw-tooth case [9]\n";
+      cout << "-deg: max degree of the polynomial description of Phase [15]\n";
       exit(0);
     }
     else if (strcmp(arg[ka],"-debug")==0) {