# include "relphot.h"

# define USE_IMAGE_ID 1

static off_t       **bin;     // link from catalog,measure to image	     
static off_t       **clist;   // catalog which supplied measurement on image 
static off_t       **mlist;   // measure reference for measurement on image  
static off_t        *Nlist;   // number of measurements on image	     
static off_t        *NLIST;   // allocated number of measurements on image   

static Image        *image;   // list of available images  
static off_t        Nimage;   // number of available images

// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
// index

# if USE_IMAGE_ID
static off_t        *imageIDs; // list of all image IDs
static off_t        *imageIdx; // list of index for image IDs 
# else
static unsigned int *start;
static unsigned int *stop;
# endif

Image *getimages (off_t *N) {

  *N = Nimage;
  return (image);
}

Image *getimage (off_t N) {
  return (&image[N]);
}

void initImages (Image *input, off_t N) {

  off_t i;

  image = input;
  Nimage = N;

# if USE_IMAGE_ID
  ALLOCATE (imageIDs, off_t, Nimage);
  ALLOCATE (imageIdx, off_t, Nimage);

  for (i = 0; i < Nimage; i++) {
    imageIdx[i] = i;
    imageIDs[i] = image[i].imageID;
  }
  llsortpair (imageIDs, imageIdx, Nimage);
# else
  ALLOCATE (start,   unsigned, Nimage);
  ALLOCATE (stop,    unsigned, Nimage);

  for (i = 0; i < Nimage; i++) {
    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
  }
# endif
}

off_t getImageByID (off_t ID) {

  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
  // use bisection to find the specified image ID

# if USE_IMAGE_ID
  off_t Nlo, Nhi, N;

  Nlo = 0; Nhi = Nimage;
  while (Nhi - Nlo > 10) {
    N = 0.5*(Nlo + Nhi);
    if (imageIDs[N] < ID) {
      Nlo = MAX(N, 0);
    } else {
      Nhi = MIN(N + 1, Nimage);
    }
  }

  for (N = Nlo; N < Nhi; N++) {
    if (imageIDs[N] == ID)
      return (imageIdx[N]);
  }
# endif

  return (-1);
}

void initImageBins (Catalog *catalog, int Ncatalog) {

  off_t i, j;

  ALLOCATE (bin, off_t *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
  }

  ALLOCATE (Nlist, off_t, Nimage);
  ALLOCATE (NLIST, off_t, Nimage);
  ALLOCATE (clist, off_t *, Nimage);
  ALLOCATE (mlist, off_t *, Nimage);

  for (i = 0; i < Nimage; i++) {
    Nlist[i] = 0;
    NLIST[i] = 100;
    ALLOCATE (clist[i], off_t, NLIST[i]);
    ALLOCATE (mlist[i], off_t, NLIST[i]);
  }
}

void freeImageBins (int Ncatalog) {

  off_t i;

  for (i = 0; i < Ncatalog; i++) {
    free (bin[i]);
  }
  free (bin);
  for (i = 0; i < Nimage; i++) {
    free (clist[i]);
    free (mlist[i]);
  }
  free (clist);
  free (mlist);
}

/* select all image equivalent to the current photcode */
void findImages (Catalog *catalog, int Ncatalog) {

  off_t j;
  int i, ecode;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
      if (photcode[0].code != ecode) continue;
      matchImage (catalog, j, i);
    }
  }
}

int findCCD (off_t idx, off_t meas, int cat, Measure *measure) {

  int ccdnum;
  double X, Y;
  char *pname, *filter, *p, base[256];

  /* identify the ccd on the basis of the photcode name */
  pname = GetPhotcodeNamebyCode (image[idx].photcode);
  filter = photcode[0].name;
  sprintf (base, "%s.%s.", MOSAICNAME, filter);
  if (strncmp (pname, base, strlen (base))) return (FALSE);
  p = pname + strlen(base);

  /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR 
     for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
     and nn is an integer -> ccdnum */
      
  /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode

  while (!isdigit(*p) && *p) p++;
  if (*p == 0) return (FALSE);
  ccdnum = atoi (p);

  /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
     constructed, there will be null values for undefined ccdnums */

  // old code to add this measurement to the grid cell for this chip 
  // ave = measure[0].averef;
  // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
  // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
  // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);

  // XXX we can now use these values (but need to be careful about old formats)
  X = measure[0].Xccd;
  Y = measure[0].Yccd;
  setGridMeasure (meas, cat, X, Y, ccdnum);

  return (TRUE);
}

# if USE_IMAGE_ID
void matchImage (Catalog *catalog, off_t meas, int cat) {

  off_t idx, ID;
  int status;
  Measure *measure;
  
  measure = &catalog[cat].measure[meas];

  ID = measure[0].imageID;
  idx = getImageByID (ID);
  if (idx == -1) {
    fprintf (stderr, "can't match detection to image?\n");
    abort();
  }

  if (USE_GRID) {
    status = findCCD (idx, meas, cat, measure);
    if (!status) {
      fprintf (stderr, "failed to determine CCD for %lld, %d\n", (long long) meas, cat);
      abort();
    }
  }

  // index for (catalog, measure) -> image
  bin[cat][meas] = idx;

  // index for image, Nentry -> catalog
  clist[idx][Nlist[idx]] = cat;

  // index for image, Nentry -> measure
  mlist[idx][Nlist[idx]] = meas;
  Nlist[idx] ++;

  if (Nlist[idx] == NLIST[idx]) {
    NLIST[idx] += 100;
    REALLOCATE (clist[idx], off_t, NLIST[idx]);
    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
  }	

  return;
}

# else
// this is the time-based match
void matchImage (Catalog *catalog, off_t meas, int cat) {

  off_t i;
  int status;
  Measure *measure;
  
  measure = &catalog[cat].measure[meas];
  for (i = 0; i < Nimage; i++) {
    if (image[0].photcode == -1) continue;
    if (measure[0].photcode != image[i].photcode) continue;
    if (measure[0].t < start[i]) continue;
    if (measure[0].t > stop[i]) continue;
    
    if (USE_GRID) {
      status = findCCD(i, meas, cat);
      if (!status) continue;
    }

    // index for (catalog, measure) -> image
    bin[cat][meas] = i;

    // index for image, Nentry -> catalog
    clist[i][Nlist[i]] = cat;

    // index for image, Nentry -> measure
    mlist[i][Nlist[i]] = meas;
    Nlist[i] ++;

    if (Nlist[i] == NLIST[i]) {
      NLIST[i] += 100;
      REALLOCATE (clist[i], off_t, NLIST[i]);
      REALLOCATE (mlist[i], off_t, NLIST[i]);
    }	
    return;
  }
  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
}
# endif

off_t getImageEntry (off_t meas, int cat) {

  off_t i;

  i = bin[cat][meas];
  return (i);
}

float getMcal (off_t meas, int cat) {

  off_t i;
  float value;

  i = bin[cat][meas];
  if (i == -1) return (NAN);

  if (image[i].flags & IMAGE_BAD) return (NAN);  
  value = image[i].Mcal;
  return (value);
}

Coords *getCoords (off_t meas, int cat) {

  off_t i;

  i = bin[cat][meas];
  if (i == -1) return (NULL);
  return (&image[i].coords);
}

/* determine Mcal values for all images */
void setMcal (Catalog *catalog, int PoorImages) {

  off_t i, j, m, c, n, N, Nmax, mark, bad;
  int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys;
  float Msys, Mrel, Mmos, Mgrid;
  double *list, *dlist;
  StatType stats;

  if (FREEZE_IMAGES) return;

  if (PoorImages) {
    IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
  }

  Nmax = 0;
  for (i = 0; i < Nimage; i++) {
    Nmax = MAX (Nmax, Nlist[i]);
  }
  ALLOCATE (list, double, Nmax);
  ALLOCATE (dlist, double, Nmax);

  Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;

  for (i = 0; i < Nimage; i++) {
    
    /* on PoorImages run, skip good images */
    if (PoorImages) {
      bad = image[i].flags & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
      if (!bad) continue;
    }      

    N = 0;
    for (j = 0; j < Nlist[i]; j++) {
      
      m = mlist[i][j];
      c = clist[i][j];
      
      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
	  Nbad++;
	  continue;
      }
      Mmos  = getMmos  (m, c);
      if (isnan(Mmos)) {
	  Nmos ++;
	  continue;
      }
      Mgrid = getMgrid (m, c);
      if (isnan(Mgrid)) {
	  Ngrid++;
	  continue;
      }
      Mrel  = getMrel  (catalog, m, c);
      if (isnan(Mrel)) {
	  Nrel ++;
	  continue;
      }
      
      n = catalog[c].measure[m].averef;
      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
      if (isnan(Msys)) {
	Nsys++;
	continue;
      }
      list[N] = Msys - Mrel - Mmos - Mgrid;
      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
      N++;
    }
    /* Nlist[i] is all measurements, N is good measurements */

    /* too few good measurements or too many bad measurements */
    if (!PoorImages) {
      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
      if (mark) {
	image[i].flags |= ID_IMAGE_FEW;
	Nfew ++;
      } else {
	image[i].flags &= ~ID_IMAGE_FEW;
      }      
    }
    
    liststats (list, dlist, N, &stats);
    image[i].Mcal  = stats.mean;
    image[i].dMcal = stats.sigma;
    image[i].Xm    = 100.0*log10(stats.chisq);
  }
  free (list);
  free (dlist);

  fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys);

  if (PoorImages) {
    IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }
  return;
}

/* mark image if: abs(Mcal) too large, dMcal too large */
void clean_images () {

  int mark, Nmark;
  off_t i, N;
  double *mlist, *slist, *dlist;
  double MaxOffset, MaxScatter, MedOffset;
  StatType stats;

  if (FREEZE_IMAGES) return;

  if (VERBOSE) fprintf (stderr, "marking poor images\n");

  ALLOCATE (mlist, double, Nimage);
  ALLOCATE (slist, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  for (i = N = 0; i < Nimage; i++) {
    if (image[i].flags & IMAGE_BAD) continue;
    mlist[N] = fabs (image[i].Mcal);
    slist[N] = image[i].dMcal;
    dlist[N] = 1;
    N++;
  }
  initstats ("MEAN");
  liststats (mlist, dlist, N, &stats);
  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
  MedOffset = stats.median;
  liststats (slist, dlist, N, &stats);
  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
  
  Nmark = 0;
  for (i = 0; i < Nimage; i++) {
    mark = FALSE;
    image[i].flags &= ~ID_IMAGE_POOR;
    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
    if (mark) { 
      Nmark ++;
      image[i].flags |= ID_IMAGE_POOR;
    } else {
      image[i].flags &= ~ID_IMAGE_POOR;
    }
  }

  fprintf (stderr, "%d images marked poor\n", Nmark);
  initstats (STATMODE);
  free (mlist);
  free (slist);
  free (dlist);
}

void plot_images () {

  off_t i, bin;
  double *xlist, *Mlist, *dlist;
  Graphdata graphdata;

  if (FREEZE_IMAGES) return;

  ALLOCATE (xlist, double, Nimage);
  ALLOCATE (dlist, double, Nimage);
  ALLOCATE (Mlist, double, Nimage);

  /**** dMcal vs airmass ****/
  for (i = 0; i < Nimage; i++) {
    Mlist[i] = image[i].Mcal;
    dlist[i] = image[i].dMcal;
    xlist[i] = image[i].secz;
  }

  plot_defaults (&graphdata);
  graphdata.ymin = PlotdMmin; 
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
  plot_defaults (&graphdata);
  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", "%s.Mcal.dMcal.png", OUTROOT);

# define NBIN 200
  REALLOCATE (xlist, double, NBIN);
  REALLOCATE (Mlist, double, NBIN);

  /**** dMcal histgram ****/
  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nimage; i++) {
    bin = image[i].dMcal / 0.00025;
    bin = MAX (0, MIN (NBIN - 1, bin));
    Mlist[bin] += 1.0;
  }

  plot_defaults (&graphdata);
  graphdata.style = 1;
  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT);

  free (dlist);
  free (xlist);
  free (Mlist);
}

StatType statsImageN (Catalog *catalog) {

  off_t i, j, m, c, n, N;
  double *list, *dlist;
  float Mcal, Mmos, Mgrid;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {
    if (image[i].flags & IMAGE_BAD)  continue;

    N = 0;
    for (j = 0; j < Nlist[i]; j++) {

      m = mlist[i][j];
      c = clist[i][j];

      Mcal  = getMcal  (m, c);
      if (isnan(Mcal)) continue;
      Mmos  = getMmos  (m, c);
      if (isnan(Mmos)) continue;
      Mgrid = getMgrid (m, c);
      if (isnan(Mgrid)) continue;
      N++;
    }
    list[n] = N;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImageX (Catalog *catalog) {

  off_t i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].flags & IMAGE_BAD)  continue;

    list[n] = pow (10.0, 0.01*image[i].Xm);
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImageM (Catalog *catalog) {

  off_t i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].flags & IMAGE_BAD)  continue;

    list[n] = image[i].Mcal;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImagedM (Catalog *catalog) {

  off_t i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].flags & IMAGE_BAD)  continue;

    list[n] = image[i].dMcal;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

