// GPH 424.1   Planck HFI-L2 Simple Map Making
// Eric Aubourg     CEA/DAPNIA/SPP
// $Id: quickmap_p.cc,v 1.1 2001-11-29 18:19:32 aubourg Exp $

#include <stdlib.h>
#include <stdio.h>
#include "toi.h"
#include "toiprocessor.h"
#include "fitstoirdr.h"
#include "fitstoiwtr.h"
#include "toimanager.h"
#include "toisegment.h"

#include "sambainit.h"
#include "toi2gmap.h"
#include "fitsspherehealpix.h"
#include "localmap.h"
#include "fitslocalmap.h"

#include <util/PlatformUtils.hpp>
#include <sax/SAXException.hpp>
#include <sax/SAXParseException.hpp>
#include <parsers/DOMParser.hpp>
#include <dom/DOM_DOMException.hpp>


void usage(void);

void usage(void) {
  cerr << "usage: gph424_1gp xml_options_file" << endl;
  exit(-1);
}

// Prefs structure

struct Gph424_1_Prefs {
  struct {
    struct {
      string fname;
      string coord1;  // TOI name
      string coord2;  // TOI name
      unsigned long typcoord;
    } pointing;
    struct {
      string fname;
      string toiname;
    } signal;
  } input;

  struct {
    string fname;
    string wfname; // weight
    unsigned long typcoord;
    string typmap;  // local, healpix...
    struct {
      int nlat;
    } healpix_opt;
    struct {
      int npx, npy;
      double org1,org2;
      double anglex, angley;
    } local_opt;
  } output;

  int nproc;
  long snb, sne;

  void dump() {
    cout << "*** GPH424-1 settings\n";
    cout << "  input\n";
    cout << "  | pointing\n";
    cout << "  | | fname       " << input.pointing.fname << "\n";
    cout << "  | | coord1      " << input.pointing.coord1 << "\n";
    cout << "  | | coord2      " << input.pointing.coord2 << "\n";
    cout << "  | | typcoord    " 
	 <<  ((input.pointing.typcoord & TypCoordGal) ? "G" : "E") 
	 <<  ((input.pointing.typcoord & TypCoordHD) ? "H" : 
	      (input.pointing.typcoord & TypCoordDD) ? "D" : "R") << "\n";
    cout << "  | signal\n";
    cout << "  | | fname       " << input.signal.fname << "\n";
    cout << "  | | toiname     " << input.signal.toiname << "\n";
    cout << "  output\n";
    cout << "  | fname         " << output.fname << "\n";
    cout << "  | wfname        " << output.wfname << "\n";
    cout << "  | typcoord      " 
	 <<  ((output.typcoord & TypCoordGal) ? "G" : "E") 
	 <<  ((output.typcoord & TypCoordHD) ? "H" : 
	      (input.pointing.typcoord & TypCoordDD) ? "D" : "R")  << "\n";
    cout << "  | typmap        " << output.typmap << "\n";
    if (output.typmap == "healpix") {
      cout << "  | | nlat           " << output.healpix_opt.nlat << "\n";
    } else if (output.typmap == "local") {
      cout << "  | | npix           " << output.local_opt.npx    << ", " << output.local_opt.npy << "\n";
      cout << "  | | origin         " << output.local_opt.org1   << ", " << output.local_opt.org2 << "\n";
      cout << "  | | extension      " << output.local_opt.anglex << ", " << output.local_opt.angley << "\n";
    }
    cout << "  nproc           " << nproc << "\n";
    if (sne >= 0 && snb >= 0) {
      cout << "  samples         " << snb << " " << sne << "\n";
    }
    cout << "***" << endl;
  }
};

DOM_Node DOMGetNode(DOM_Document& doc, string path) {
  DOM_Element el = doc.getDocumentElement();

  while (path != "") {
    if (path[0] == '/') {
      path = path.substr(1);
    }
    if (path == "") break;
    int x = path.find('/');
    string pathElem;
    if (x == -1) {
      pathElem = path;
      path = "";
    } else {
      pathElem = path.substr(0, x);
      path = path.substr(x);
    }
    DOM_Node n = el.getElementsByTagName(pathElem.c_str()).item(0);
    if (n==0) return n;
    el = (DOM_Element&) n;
  }
  return el;
}

string DOMGetString(DOM_Document& doc, string path, string def="") {
  DOM_Node n = DOMGetNode(doc, path);
  if (n==0) return def;
  char* s = n.getFirstChild().getNodeValue().transcode();
  char* p=s; while ((*p)==' ') p++;
  char* q=p+strlen(p)-1;
  while (q>p && *q==' ') {
    *q = '\0';
    q--;
  }
  string ss = p;
  free(s);
  return ss;
}

bool DOMHasOption(DOM_Document& doc, string path) {
  DOM_Node n = DOMGetNode(doc, path);
  return n != 0;
}

void parsePrefs(string prefFile, Gph424_1_Prefs& prefs) {
  try {
    XMLPlatformUtils::Initialize();
  } catch (const XMLException& toCatch) {
    cerr << "Error during XML initialization! :\n"
	 << toCatch.getMessage() << endl;
    exit(-1);
  }

  DOMParser* parser = new DOMParser;
  
  try {
    parser->parse(prefFile.c_str());
  } catch (const XMLException& toCatch) {
    cerr << "\nError during parsing: '" << prefFile << "'\n"
	 << "Exception message is:  \n"
	 << toCatch.getMessage() << "\n" << endl;
    exit(-1);
  } catch (const DOM_DOMException& toCatch) {
    cerr << "\nDOM Error during parsing: '" << prefFile << "'\n"
	 << "DOMException code is:  \n"
	 << toCatch.code << "\n" << endl;
    exit(-1);
  } catch (...) {
    cerr << "\nUnexpected exception during parsing: '" << prefFile << "'\n";
    exit(-1);  
  }

  DOM_Document doc = parser->getDocument();
  
  prefs.input.pointing.fname  = DOMGetString(doc, "/input/pointing/fname");
  prefs.input.pointing.coord1 = DOMGetString(doc, "/input/pointing/coord1","phi");
  prefs.input.pointing.coord2 = DOMGetString(doc, "/input/pointing/coord2","theta");
  bool inequ = DOMHasOption(doc, "/input/pointing/equatorial");
  bool ingal = DOMHasOption(doc, "/input/pointing/galactic");
  prefs.input.pointing.typcoord = ingal ? TypCoordGal : TypCoordEq;
  bool inhd  = DOMHasOption(doc, "/input/pointing/hourdeg");
  bool indd  = DOMHasOption(doc, "/input/pointing/degdeg");
  bool inrr  = DOMHasOption(doc, "/input/pointing/radian");
  if (inhd) {
    prefs.input.pointing.typcoord |= TypCoordHD;
  } else if (indd) {
    prefs.input.pointing.typcoord |= TypCoordDD;
  } else if (inrr) {
    prefs.input.pointing.typcoord |= TypCoordRR;
  } else if (prefs.input.pointing.typcoord & TypCoordEq) {
    prefs.input.pointing.typcoord |= TypCoordHD;
  } else {
    prefs.input.pointing.typcoord |= TypCoordDD;
  }

  prefs.input.signal.fname    = DOMGetString(doc, "/input/signal/fname");
  prefs.input.signal.toiname  = DOMGetString(doc, "/input/signal/toiname");

  prefs.output.fname          = DOMGetString(doc, "/output/fname");
  prefs.output.wfname         = DOMGetString(doc, "/output/wfname");
  bool outequ = DOMHasOption(doc, "/output/equatorial");
  bool outgal = DOMHasOption(doc, "/output/galactic");
  prefs.output.typcoord = outgal ? TypCoordGal : TypCoordEq;
  bool outhd  = DOMHasOption(doc, "/output/hourdeg");
  bool outdd  = DOMHasOption(doc, "/output/degdeg");
  bool outrr  = DOMHasOption(doc, "/output/radian");
  if (outhd) {
    prefs.output.typcoord |= TypCoordHD;
  } else if (outdd) {
    prefs.output.typcoord |= TypCoordDD;
  } else if (outrr) {
    prefs.output.typcoord |= TypCoordRR;
  } else if (prefs.output.typcoord & TypCoordEq) {
    prefs.output.typcoord |= TypCoordHD;
  } else {
    prefs.output.typcoord |= TypCoordDD;
  } 

  prefs.output.typmap = DOMHasOption(doc, "/output/local") ? "local" :
    (DOMHasOption(doc, "/output/healpix") ? "healpix" : "healpix");
  prefs.output.healpix_opt.nlat      = atoi(DOMGetString(doc, "/output/healpix/nlat", "128").c_str());

  string s                    = DOMGetString(doc, "/output/local/npx");
  prefs.output.local_opt.npx = prefs.output.local_opt.npy = -1;
  int nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.npx, &prefs.output.local_opt.npy);
  if (nc == 1) prefs.output.local_opt.npy = prefs.output.local_opt.npx;

  s                           = DOMGetString(doc, "/output/local/center");
  prefs.output.local_opt.org1 = prefs.output.local_opt.org2 = 0;
  sscanf(s.c_str(), "%lg %lg", &prefs.output.local_opt.org1, &prefs.output.local_opt.org2);  

  s                           = DOMGetString(doc, "/output/local/angle");
  prefs.output.local_opt.anglex = prefs.output.local_opt.angley = 0;
  nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.anglex, &prefs.output.local_opt.angley);
  if (nc == 1) prefs.output.local_opt.angley = prefs.output.local_opt.anglex;

  prefs.nproc                 = atoi(DOMGetString(doc, "/nproc", "1").c_str());

  s                           = DOMGetString(doc, "/samples", "-1 -1");
  sscanf(s.c_str(), "%ld %ld", &prefs.snb, &prefs.sne);
}

int main(int argc, char** argv) {

  // Read prefs
  Gph424_1_Prefs prefs;
  parsePrefs(argv[1], prefs);
  prefs.dump();

  int np = prefs.nproc;

  // Assemble pipeline

  TOIManager* mgr = TOIManager::getManager();

  // Two input files for each data chunk

  vector<FITSTOIReader*> inPoint;
  vector<FITSTOIReader*> inSig;
  for (int i=0; i<np; i++) {
    inPoint.push_back(new FITSTOIReader(prefs.input.pointing.fname));
    inSig.push_back(new FITSTOIReader(prefs.input.signal.fname));
    inPoint[i]->setImplicitSN();
    inSig[i]->setImplicitSN();
  }

  // Prepare maps for output

  vector<PixelMap<r_8>*> map;
  vector<PixelMap<r_8>*> wmap;

  for (int i=0; i<np; i++) {
    if (prefs.output.typmap == "local") {
      LocalMap<r_8>* lmap = 
	new LocalMap<r_8>(prefs.output.local_opt.npx, prefs.output.local_opt.npy);
      lmap->SetOrigin(prefs.output.local_opt.org1, prefs.output.local_opt.org2);
      lmap->SetSize(prefs.output.local_opt.anglex, prefs.output.local_opt.angley);
      map.push_back(lmap);
      wmap.push_back(new LocalMap<r_8>(*lmap));
    } else if (prefs.output.typmap == "healpix") {
      map.push_back(new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat));
      wmap.push_back(new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat));
    } else {
      cout << "typmap " << prefs.output.typmap << " unknown" << endl;
      exit(-1);
    }
  }

  // One map projector per data chunk

  vector<TOI2GMap*> toi2m;
  {for (int i=0; i<np; i++) {
    toi2m.push_back(new TOI2GMap(map[i],wmap[i]));
    toi2m[i]->SetEquinox(2000.);
    toi2m[i]->SetCoorIn((TypAstroCoord) prefs.input.pointing.typcoord);
    toi2m[i]->SetCoorOut((TypAstroCoord) prefs.output.typcoord);
  }}


  if (prefs.snb >= 0 && prefs.sne >= 0) {
    mgr->setRequestedSample(prefs.snb, prefs.sne);
  }

  int sn0 = inPoint[0]->getMinOut();
  int sn1 = inPoint[0]->getMaxOut();

  cout << "Range to be processed : " << sn0 << " - " << sn1 << endl;

  int chunksize = (sn1-sn0)/np;
  {for (int i=0; i<np; i++) {
    int isn0 = sn0 + chunksize*i;
    int isn1 = sn0 + chunksize*(i+1) - 1;
    if (i == np-1) isn1 = sn1;
    cout << "Processor " << i << " range " << isn0 << " - " << isn1 << endl;
    inPoint[i]->setMinSn(isn0);
    inSig[i]->setMinSn(isn0);
    inPoint[i]->setMaxSn(isn1);
    inSig[i]->setMaxSn(isn1);
  }}

  // Connect pipes

  {for (int i=0; i<np; i++) {
    TOISegmented * toicoord1in = new TOISegmented(prefs.input.pointing.coord1);
    inPoint[i]->addOutput(prefs.input.pointing.coord1,toicoord1in);
    toi2m[i]->addInput("Coord1In",toicoord1in);

    TOISegmented * toicoord2in = new TOISegmented(prefs.input.pointing.coord2);
    inPoint[i]->addOutput(prefs.input.pointing.coord2,toicoord2in);
    toi2m[i]->addInput("Coord2In",toicoord2in);

    TOISegmented * toibolin = new TOISegmented("toi_bolo_in");
    inSig[i]->addOutput(prefs.input.signal.toiname,toibolin);
    toi2m[i]->addInput("BoloIn",toibolin);
  }}

  // Start processing

  {for (int i=0; i<np; i++) {
    inPoint[i]->start();
    inSig[i]->start();
    toi2m[i]->start();
  }}
  mgr->joinAll();

  // Merge maps

  PixelMap<r_8>* smap = map[0];
  PixelMap<r_8>* swmap = wmap[0];

  {for (int i=1; i<np; i++) {
    // There should be a += operator in all maps, not only in local maps
    // *smap += *map[i];
    // *swmap += *wmap[i];
    int n = smap->NbPixels();
    for (int j=0; j<n; j++) {
      (*smap)(j) += (*map[i])(j);
      (*swmap)(j) += (*wmap[i])(j);
    }
  }}

  // Save maps

  {
    FitsOutFile sfits(prefs.output.fname,FitsFile::clear);
    cout<<"Creating map fits file "<<prefs.output.fname<<endl;
    // Would need some generic mechanism here $CHECK$ $$TBD$$
    if (smap->TypeOfMap() == "RING") {
      sfits << * (SphereHEALPix<r_8>*)smap;
    } else if (smap->TypeOfMap() == "LOCAL") {
      sfits << * (LocalMap<r_8>*)smap;
    } else {
      cout << "error, unknown type " << smap->TypeOfMap() << endl;
    }
  }

  {
    FitsOutFile swfits(prefs.output.wfname,FitsFile::clear);
    cout<<"Creating map weight fits file "<<prefs.output.wfname<<endl;
    //    swfits << *wmap;
    if (swmap->TypeOfMap() == "RING") {
      swfits << * (SphereHEALPix<r_8>*)swmap;
    } else if (swmap->TypeOfMap() == "LOCAL") {
      swfits << * (LocalMap<r_8>*)swmap;
    } else {
      cout << "error, unknown type " << swmap->TypeOfMap() << endl;
    }
  }

  return 0;
}