// GPH 424.1   Planck HFI-L2 Simple Map Making
// Eric Aubourg     CEA/DAPNIA/SPP
// This version is able to produce local maps as well
// $Id: quickmap.cc,v 1.4 2003-02-24 09:18:01 aubourg Exp $

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

#include "sambainit.h"
#include "toi2map.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_1g xml_options_file" << endl;
  exit(-1);
}

// Prefs structure

struct Gph424_1_Prefs {
  struct {
    struct {
      string fname;
      string flags;
      string coord1;  // TOI name
      string coord2;  // TOI name
      unsigned long typcoord;
    } pointing;
    struct {
      string fname;
      string flags;
      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 << "  | | flags       " << input.pointing.flagfs << "\n";
    cout << "  | | coord1      " << input.pointing.coord1 << "\n";
    cout << "  | | coord2      " << input.pointing.coord2 << "\n";
    cout << "  | | typcoord    " 
	 <<  ((input.pointing.typcoord & TypCoordGal) ? "G" : "E") 
	 <<  ((input.pointing.typcoord & TypCoord1H) ? "H" : 
	      (input.pointing.typcoord & TypCoord1D) ? "D" : "R") 
	 <<  ((input.pointing.typcoord & TypCoord1C) ? "C" : "L")
	 <<  ((input.pointing.typcoord & TypCoord2H) ? "H" : 
	      (input.pointing.typcoord & TypCoord2D) ? "D" : "R") 
	 <<  ((input.pointing.typcoord & TypCoord2C) ? "C" : "L")
	 << "\n";
    cout << "  | signal\n";
    cout << "  | | fname       " << input.signal.fname << "\n";
    cout << "  | | flags       " << input.signal.flags << "\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") 
	 <<  "\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");
  if (DOMHasOption(doc, "/input/pointing/flags")) {
    prefs.input.pointing.flags=
      DOMGetString(doc, "/input/pointing/flags");
  }
  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");
  bool iscolat = DOMHasOption(doc, "/input/pointing/colat");
  if (inhd) {
    prefs.input.pointing.typcoord |= TypCoord1H|TypCoord2D;
  } else if (indd) {
    prefs.input.pointing.typcoord |= TypCoord1D|TypCoord2D;
  } else if (inrr) {
    prefs.input.pointing.typcoord |= TypCoord1R|TypCoord2R;
  } else if (prefs.input.pointing.typcoord & TypCoordEq) {
    prefs.input.pointing.typcoord = TypCoordEqStd;
  } else {
    prefs.input.pointing.typcoord = TypCoordGalStd;
  }
  if (iscolat) {
    prefs.input.pointing.typcoord |= TypCoord2C;
  }

  prefs.input.signal.fname    = DOMGetString(doc, "/input/signal/fname");
  prefs.input.signal.toiname  = DOMGetString(doc, "/input/signal/toiname");
  if (DOMHasOption(doc, "/input/signal/flags")) {
    prefs.input.signal.flags =
      DOMGetString(doc, "/input/signal/flags");
  }
  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 |= TypCoord1H | TypCoord2D;
  } else if (outdd) {
    prefs.output.typcoord |= TypCoord1D | TypCoord2D;
  } else if (outrr) {
    prefs.output.typcoord |= TypCoord1R | TypCoord2R;
  } else if (prefs.output.typcoord & TypCoordEq) {
    prefs.output.typcoord = TypCoordEqStd;
  } else {
    prefs.output.typcoord |= TypCoordGalStd;
  } 

  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();

  // Assemble pipeline

  TOIManager* mgr = TOIManager::getManager();

  // INPUT TOIs

  PIOTOIReader inPoint1(prefs.input.pointing.fname, prefs.input.pointing.coord1,
			prefs.input.pointing.flags);
  PIOTOIReader inPoint2(prefs.input.pointing.fname, prefs.input.pointing.coord2,
			prefs.input.pointing.flags);
  PIOTOIReader inSig   (prefs.input.signal.fname, prefs.input.signal.toiname,
			prefs.input.signal.flags);


  // Prepare maps for output
  
  PixelMap<r_8>* map;
  PixelMap<r_8>* wmap;
  
  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 = lmap;
    wmap = new LocalMap<r_8>(*lmap);
  } else if (prefs.output.typmap == "healpix") {
    map = new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat);
    wmap = new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat);
  } else {
    cout << "typmap " << prefs.output.typmap << " unknown" << endl;
    exit(-1);
  }
  
  TOI2Map toi2m(map,wmap);

  toi2m.SetEquinox(2000.);
  toi2m.SetCoorIn((TypAstroCoord) prefs.input.pointing.typcoord);
  toi2m.SetCoorMap((TypAstroCoord) prefs.output.typcoord);

  TOISegmented * toicoord1in = new TOISegmented(prefs.input.pointing.coord1);
  inPoint1.addOutput(prefs.input.pointing.coord1,toicoord1in);
  toi2m.addInput("Coord1In",toicoord1in);

  TOISegmented * toicoord2in = new TOISegmented(prefs.input.pointing.coord2);
  inPoint2.addOutput(prefs.input.pointing.coord2,toicoord2in);
  toi2m.addInput("Coord2In",toicoord2in);

  TOISegmented * toibolin = new TOISegmented("toi_bolo_in");
  inSig.addOutput(prefs.input.signal.toiname,toibolin);
  toi2m.addInput("BoloIn",toibolin);

  // Start processing

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

  mgr->startAll();
  mgr->waitForAll();

  // Save maps

  // Huh ? Should we create a Map(group) ? Should we create a MapObject ?
  // Should we open ? Where do we define the size ? How does this work ??

  PIOCreateMAP(prefs.output.fname, 
	       prefs.output.typcoord == TypCoordGal ? "GALACTIC" : "ECLIPTIC", 
	       "RING", 
	       prefs.output.healpix_opt.nlat);

  PIOGroup* mapGroup = PIOOpenMAP(prefs.output.fname,
				  "w");

  PIOCreateMAPObject("MAP", "PIODOUBLE", mapGroup);

  PIOWriteMAP(map->GetData(),
	      "MAP",
	      "PIODOUBLE",
	      "",
	      mapGroup);

  PIOCreateMAPObject("WMAP", "PIODOUBLE", mapGroup);

  PIOWriteMAP(wmap->GetData(),
	      "WMAP",
	      "PIODOUBLE",
	      "",
	      mapGroup);

  PIOCloseMAP(mapGroup);

  {
    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 (map->TypeOfMap() == "RING") {
      sfits << * (SphereHEALPix<r_8>*)map;
    } else if (map->TypeOfMap() == "LOCAL") {
      cout << "error, local maps not yet in PIOLib" << endl;
    } else {
      cout << "error, unknown type " << map->TypeOfMap() << endl;
    }
  }

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

  delete map;
  delete wmap;

  return 0;
}