tsrif.cpp

1 gps基本算法

      srif.doSequential = false;

      // starting value
      sol(0)   = sol(1)   = 0.0;
      doIt(data,sol,cov);
      cout << endl << " ---------------------------------------------------------\n"
         << " do it again, starting at (4,4) to get the false minimum.\n";
      sol(0)   = sol(1)   = 4.0;
      data = D;
      srif.zeroAll();
      srif.iterationsLimit = 40;
      doIt(data,sol,cov);
   }
   else if(test == 11) {
      //{
      //   ifstream ifs("globaltemp.data");
      //   if(!ifs) {
      //      cout << "could not open globaltemp.data .. abort\n";
      //      return 0;
      //   }
      //   const int BUFF_SIZE=100;
      //   char buffer[BUFF_SIZE];
      //   int year=0;
      //   double temp;
      //   vector<int> years;
      //   vector<double> temps;
      //   while(ifs.getline(buffer,BUFF_SIZE)) {
      //      if(buffer[0]=='#') continue;
      //      string line=buffer;
      //      StringUtils::stripTrailing(buffer,'\r');
      //      string::size_type p;
      //      vector<string> fs;
      //      while(line.size() > 0) {
      //         p = line.find(" ");
      //         if(p==string::npos) p=line.size();
      //         if(p != 0) fs.push_back(line.substr(0,p));
      //         if(p >= line.size()) break;
      //         line.erase(0,p+1);
      //      }
      //      if(StringUtils::asInt(fs[0]) == year) continue;    // percentages line
      //      year = StringUtils::asInt(fs[0]);
      //      temp = StringUtils::asDouble(fs[13]);
      //      years.push_back(year);
      //      temps.push_back(temp);
      //      cout << year << " " << fixed << setprecision(3) << setw(6)
      //         << temp << endl;
      //   }
      //   ifs.close();
      //   for(i=0; i<years.size(); i++) cout << "," << years[i]; cout << endl;
      //   for(i=0; i<temps.size(); i++) cout << "," << temps[i]; cout << endl;
      //}
      cout << " a simple but noisy polynomial fit" << endl;
      cout << " to global temperature anomaly data." << endl;
      cout << " do it in 8 batches (7 of 20, then 1 of 9) points each." << endl;
      cout << " cf. www.cru.uea.ac.uk/cru/data/temperature" << endl;
      NL += "c0"; NL += "c1"; NL += "c2"; NL += "c3"; NL += "c4"; NL += "c5";
      N = NL.size();
      M = batchlen = batchsize = 20;
      Vector<double> f(M);
      Matrix<double> Coef(8,N);  // save the results from each batch
      data.resize(M);
      truth.resize(N);
      sol.resize(N);
      partials.resize(M,N);
      truth = 0.0;         // don't know it

      srif = SRIFilter(NL);
      srif.doLinearize = false;
      srif.doSequential = true;

      sol = 0.0;
      Coef = 0.0;
      for(batch=0; batch<8; batch++) {
         if(batch == 7) {
            batchlen = 9;
            data.resize(batchlen);         // actually not necessary...
         }
         for(i=0; i<batchlen; i++) data(i) = temps11[batch*batchsize+i];
         cout << "\n -------- Batch " << batch+1 << " --------" << endl;
         doIt(data,sol,cov);
         for(i=0; i<N; i++) Coef(batch,i)=sol(i);
      }

      // print out the results, and generate a plot
      cout << "\n Coefficients (batch vs coefficients) :\n" << fixed
         << setw(10) << setprecision(3) << Coef << endl;

      ofstream ofs("tsrif11.dat",ios::out);
      if(!ofs) {
         cout << " Could not open tsrif11.dat .. abort plot\n";
      }
      else {
         // print all the data and all the solutions -- do in one big batch
         batch = 0;         // this necessary for LSFunc to operate correctly
         M = batchlen = batchsize = 149;
         Matrix<double> F(M,8);
         f.resize(M);
         partials.resize(M,N);
         for(j=0; j<8; j++) {
            sol = Coef.rowCopy(j);              // pick out the jth solution
            LSFunc(sol,f,partials);             // evaluate at all 149 times
            //for(i=0; i<M; i++) F(i,j)=f(i);   // use slice instead
            //                    (matrix,col index,slice=(start,length,stride))
            MatrixColSlice<double> Fcol(F,j,std::slice(0,M,1));
            Fcol = f;                           // copy into (the slice of) F
            //MatrixRowSlice<double> Frow(F,17,std::slice(0,7,1));
            //cout << "Col slice is\n" << Fcol << endl;
            //cout << "Row slice is\n" << Frow << endl;
         }
         // print all the evaluated polynomials at all the times
         for(i=0; i<M; i++) {
            // times and data
            ofs << setw(4) << years11[i] << " " << f63 << temps11[i];
            // solutions
            for(j=0; j<8; j++) ofs << " " << f63 << F(i,j);
            ofs << endl;
         }
         ofs.close();
         cout << " Created tsrif11.dat" << endl;
         ofs.open("tsrif11.gp",ios::out);
         if(!ofs) {
            cout << " Could not open tsrif11.gp .. abort plot\n";
         }
         else {
            ofs << "set title \"Global temperature anomaly - "
               << "sequential fits of order " << N-1 << " in 7 batches of 20 points "
               << "and 1 of 9 ... tsrif(11)\"\n";
            ofs << "set xlabel \"Year\"\n";
            ofs << "set ylabel \"Temperature anomaly\"\n";
            ofs << "unset mouse\n";
#ifndef _WIN32
            ofs << "set term x11 enhanced font \"luxi sans,17\"\n";  // linux only
#endif
            ofs << "set key left\n";
            ofs << "#set xrange [40:60]\n";
            ofs << "set yrange [-0.6:0.6]\n";
            ofs << "set style line 1 lt 8 lw 2\n";
            ofs << "plot \"tsrif11.dat\" using 1:2 t \"dT\" with points\n";
            ofs << "replot \"tsrif11.dat\" using 1:3 t \"fit1\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:4 t \"fit2\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:5 t \"fit3\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:6 t \"fit4\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:7 t \"fit5\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:8 t \"fit6\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:9 t \"fit7\" with lines\n";
            ofs << "replot \"tsrif11.dat\" using 1:10 t \"fit8\" with lines ls 1\n";
            ofs.close();
            cout << " Created file tsrif11.gp .. try gnuplot tsrif11.gp\n";
         }
      }
   }
   else if(test == 12) {
      cout << " a simple but noisy polynomial fit";
      cout << " to global temperature anomaly data." << endl;
      cout << " do it in one big batch of 149 points,";
      cout << " and plot the data." << endl;
      cout << " cf. www.cru.uea.ac.uk/cru/data/temperature" << endl;
      NL += "c0"; NL += "c1"; NL += "c2"; NL += "c3"; NL += "c4"; NL += "c5";
      //NL += "c6"; NL += "c7"; NL += "c8";
      //NL += "c9"; NL += "c10"; NL += "c11";
      batchlen=batchsize=M=149; N=NL.size();
      Vector<double> f(M);
      data.resize(M);
      truth.resize(N);
      sol.resize(N);
      partials.resize(M,N);
      truth = 0.0;         // don't know it

      srif = SRIFilter(NL);
      srif.doLinearize = false;
      srif.doSequential = false;

      sol = 0.0;
      batch = 0;     // for LSFunc
      for(i=0; i<M; i++) data(i) = temps11[i];
      doIt(data,sol,cov);

      ofstream ofs("tsrif12.dat",ios::out);
      if(!ofs) {
         cout << " Could not open tsrif12.dat .. abort plot\n";
      }
      else {
         LSFunc(sol,f,partials);
         for(i=0; i<M; i++) {
            ofs << setw(4) << years11[i]
               << " " << f63 << temps11[i]
               << " " << f63 << f[i] << endl;
         }
         ofs.close();
         cout << " Created tsrif12.dat" << endl;
         ofs.open("tsrif12.gp",ios::out);
         if(!ofs) {
            cout << " Could not open tsrif12.gp .. abort plot\n";
         }
         else {
            ofs << "set title \"Global temperature anomaly - fit of order "
               << N-1 << " tsrif(12)\"\n";
            ofs << "set xlabel \"Year\"\n";
            ofs << "set ylabel \"Temperature anomaly\"\n";
            ofs << "unset mouse\n";
#ifndef _WIN32
            ofs << "set term x11 enhanced font \"luxi sans,17\"\n";  // linux only
#endif
            ofs << "set key left\n";
            ofs << "#set xrange [40:60]\n";
            ofs << "#set yrange [250:290]\n";
            ofs << "plot \"tsrif12.dat\" using 1:2 t \"dT\" with points\n";
            ofs << "replot \"tsrif12.dat\" using 1:3 t \"fit\" with lines\n";
            ofs.close();
            cout << "\n Created file tsrif12.gp .. try gnuplot tsrif12.gp\n";
         }
      }
   }
   else if(test == 13 || test == 14) {
         // see discussion  below under test == 14
      double tau=50.0,sig=0.1,B=0.0;       // data residuals are about 0.1
         // input may come from command line
      if(inputsig > 0) sig=inputsig;
      if(inputtau > 0) tau=inputtau;

      cout << "Try fitting the global temperature data of options 11 and 12\n";
      cout << "another way. Divide the timeline into N equal parts. In each part,\n";
      cout << "fit the data to a constant (13). Then do it again (14), but with\n";
      cout << "a priori information which correlates the solution elements\n";
      cout << "using a first order Markov (random walk) process.\n";

      NL += "c0"; NL += "c1"; NL += "c2"; NL += "c3"; NL += "c4"; NL += "c5";
      NL += "c6"; NL += "c7"; NL += "c8"; NL += "c9";
      N = NL.size();

      batchlen = batchsize = M = 150;
      data.resize(M);
      truth.resize(N);
      sol.resize(N);
      partials.resize(M,N);
      for(i=0; i<N; i++) {
         truth(i) = 0.0;
         for(int j=i*(150/N); j<(i+1)*(150/N); j++)
            truth(i) += temps11[j];
         truth(i) /= (150/N);
      }

      srif = SRIFilter(NL);
      srif.doLinearize = false;
      srif.doSequential = false;

      sol = 0.0;
      batch = 0;     // for LSFunc
      for(i=0; i<M; i++) data(i) = temps11[i];

      if(test == 13) cout << " Don't correlate the state elements" << endl;
      if(test == 14) {
         cout << " Add correlation to the state elements" << endl;
         // State = 0
         // Cov = (sig2,B2=const, e=exp(-dt/tau) where dt is
         //                                the time spacing of the state elements....
         // [   1    e    e^2   ]
         // [   e    1    e     ] * sig2 + B2
         // [   e^2  e    1     ]
         // sig2 = variance = (sigma)^2, sigma = uncertainty on state element
         // tau = time constant => stiffer for larger tau
         double d,ex;
         ex = exp(-(M/N)/tau);             // dt = M/N years
         Matrix<double> apCov(N,N);
         Vector<double> apSt(N,0.0);       // apriori state is the 'expected' values
         ident(apCov);
         for(i=0; i<N-1; i++) {            // row i
            d = ex;
            for(j=i+1; j<N; j++) {         // col j
               apCov(i,j) = apCov(j,i) = d;
               d *= ex;
            }
         }
         apCov *= sig*sig;
         //apCov += B*B;

         LabelledMatrix LC(NL,apCov);
         LC.setw(10).setprecision(6);
         LC.message(" apCov");
         cout << LC << endl;

         srif.addAPriori(apCov,apSt);
      }

      doIt(data,sol,cov);

      ofstream ofs("tsrif13.dat",ios::out);
      if(!ofs) {
         cout << " Could not open tsrif13.dat .. abort plot\n";
      }
      else {
         for(i=0; i<M; i++) {
            k = i/(M/N);
            ofs << setw(4) << years11[i]
               << " " << f63 << temps11[i]
               << " " << f63 << sol(k) << endl;
         }
         ofs.close();
         cout << " Created tsrif13.dat" << endl;
         ofs.open("tsrif13.gp",ios::out);
         if(!ofs) {
            cout << " Could not open tsrif13.gp .. abort plot\n";
         }
         else {
            ofs << "set title \"Global temperature anomaly : piecewise";
            if(test==13) ofs << " fit (13)\"\n";
            if(test==14) ofs << " constrained fit, tau="
               << setprecision(1) << tau << "yrs, sig="
               << setprecision(3) << sig << "C"
               << " (14)\"\n";
            ofs << "set xlabel \"Year\"\n";
            ofs << "set ylabel \"Temperature anomaly (deg C)\"\n";
            ofs << "unset mouse\n";
#ifndef _WIN32
            ofs << "set term x11 enhanced font \"luxi sans,17\"\n";  // linux only
#endif
            ofs << "set key left\n";
            ofs << "#set xrange [40:60]\n";
            ofs << "#set yrange [250:290]\n";
            ofs << "plot \"tsrif13.dat\" using 1:2 t \"dT\" with points\n";
            ofs << "replot \"tsrif13.dat\" using 1:3 t \"fit\" with linespoints\n";
            ofs.close();
            cout << "\n Created file tsrif13.gp .. try gnuplot tsrif13.gp\n";
         }
      }
   }
   else if(test == 15 || test == 16) {
      cout << "Dataset is " << dataset << endl;
      if(dataset==1) { msg=msg1; n15 = &n151[0]; y15 = &y151[0]; M15=M151; }
      if(dataset==2) { msg=msg2; n15 = &n152[0]; y15 = &y152[0]; M15=M152; }
      if(dataset==3) { msg=msg3; n15 = &n153[0]; y15 = &y153[0]; M15=M153; }
      if(dataset==4) { msg=msg4; n15 = &n154[0]; y15 = &y154[0]; M15=M154; }
      NL += "Bias"; NL += "Linear"; NL += "Quad";
      N=NL.size();
      batchlen=batchsize=M=M15;
      if(test == 15) cout << " a weighted polynomial fit";
      if(test == 16) cout << " a robust polynomial fit";
      cout << " to " << M  << " double difference phase data points." << endl;
      cout << msg << endl;

      Vector<double> f(M);
      Vector<double> res,wt(M);
      Stats<double> stat;
      PolyFit<double> PF;

      data.resize(M);
      truth.resize(N);
      sol.resize(N);
      partials.resize(M,N);
      truth = 0.0;         // don't know it

      // configure srif
      srif = SRIFilter(NL);
      srif.iterationsLimit = 20;
      srif.doLinearize = false;
      srif.doSequential = false;

      // ---------- do it robust