# include "relphot.h"

static int Nmax;
static double *list;
static double *dlist;

void initMrel (Catalog *catalog, int Ncatalog) {

  int i, j;
  
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].average[j].Nm);
    }
  }

  ALLOCATE (list, double, MAX (1, Nmax));
  ALLOCATE (dlist, double, MAX (1, Nmax));
}  

float getMrel (Catalog *catalog, int meas, int cat) {

  int ave;
  float value;

  ave = catalog[cat].measure[meas].averef;
  if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG);  
 
  value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M_PS;
  return (value);
}

int setMrel (Catalog *catalog, int Ncatalog) {

  int i, j, k, m, N;
  float Msys, Mcal, Mmos, Mgrid;
  StatType stats;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].average[j].code & STAR_BAD) continue;  
      m = catalog[i].average[j].offset;

      // XXX allow REF stars to be included in the calculation
      // this should be optionally set, and should allow for 
      // REF stars to be downweighted by more than their reported
      // errors.  how such info is carried is unclear...
      N = 0;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
	Mcal  = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) {
	  Mcal = Mmos = Mgrid = 0;
	} else {
	  Mmos  = getMmos  (m, i);
	  if (Mmos == NO_MAG) continue;
	  if (Mmos == NO_IMAGE) continue;
	  if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
	}

	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	list[N] = Msys - Mcal - Mmos - Mgrid;
	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
	N++;
      }
      if (N < STAR_TOOFEW) { /* too few measurements */
	catalog[i].average[j].code |= ID_STAR_FEW;
      } else {
	catalog[i].average[j].code &= ~ID_STAR_FEW;
      }	

      liststats (list, dlist, N, &stats);

      if (PhotPrimary) {
	catalog[i].average[j].M = stats.mean;
	catalog[i].average[j].dM = stats.sigma;
	catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
      } else {
	catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
	catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = stats.sigma;
	catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
      }      
    }
  }
  return (TRUE);
}

int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) {

  int i, j, k, m, N, Nmax;
  float Msys, Mcal, Mmos, Mgrid;
  double *list, *dlist;
  StatType stats;

  /* Nmeasure is now different, need to reallocate */
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].average[j].Nm);
    }
  }
  ALLOCATE (list, double, MAX (1, Nmax));
  ALLOCATE (dlist, double, MAX (1, Nmax));

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* skip stars already calibrated */
      if (catalog[i].found[j]) continue;  

      N = 0;
      m = catalog[i].average[j].offset;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
	Mcal  = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) {
	  Mcal = Mmos = Mgrid = 0;
	} else {
	  Mmos  = getMmos  (m, i);
	  if (Mmos == NO_MAG) continue;
	  if (Mmos == NO_IMAGE) continue;
	  if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
	}

	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	list[N] = Msys - Mcal - Mmos - Mgrid;
	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
	N++;
      }
      if (N < 1) continue;

      liststats (list, dlist, N, &stats);
      if (mark) catalog[i].found[j] = TRUE;

      /* use sigma or error in dM for output? */
      if (PhotPrimary) {
	catalog[i].average[j].M = stats.mean;
	catalog[i].average[j].dM = MAX (stats.sigma, stats.error);
	catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
      } else {
	catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
	catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = MAX (stats.error, stats.sigma);
	catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
      }      
    }
  }

  free (list);
  free (dlist);
  return (TRUE);
}

/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_NOCAL */
int setMcalOutput (Catalog *catalog, int Ncatalog) {

  int i, j, k, m;
  float Mcal, Mmos, Mgrid;

  MEAS_BAD = ID_MEAS_NOCAL;
  IMAGE_BAD = ID_IMAGE_NOCAL;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      m = catalog[i].average[j].offset;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
	Mcal  = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) continue;
	Mmos  = getMmos  (m, i);
	if (Mmos == NO_MAG) continue;
	if (Mmos == NO_IMAGE) continue;
	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
	catalog[i].measure[m].Mcal_PS = Mcal + Mmos + Mgrid;
      }
    }
  }
  return (TRUE);
}

void clean_stars (Catalog *catalog, int Ncatalog) {

  int i, j, Ndel, Nave, Ntot, mark;
  float dM, Xm;
  double Chisq, MaxScatter, MaxChisq;
  double *xlist, *slist, *dlist;
  StatType stats;

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

  /* find Xm median -> ChiSq lim must be > median */
  for (i = Ntot = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }
  ALLOCATE (xlist, double, Ntot);
  ALLOCATE (slist, double, Ntot);
  ALLOCATE (dlist, double, Ntot);
  for (i = Ntot = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      if (catalog[i].average[j].code & STAR_BAD) continue;
      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
      if (Xm == -1) continue;
      Chisq = pow (10.0, 0.01*Xm);
      xlist[Ntot] = Chisq;
      slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
      dlist[Ntot] = 1;
      Ntot ++;
    }
  }
  
  initstats ("MEAN");
  liststats (xlist, dlist, Ntot, &stats);
  MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
  liststats (slist, dlist, Ntot, &stats);
  MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
  fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);

  Ndel = Nave = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
      Xm = PhotPrimary ? catalog[i].average[j].Xm    : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
      Chisq = pow (10.0, 0.01*Xm);
      mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq);
      if (mark) {
	catalog[i].average[j].code |= ID_STAR_POOR;
	Ndel ++;
      } else {
	catalog[i].average[j].code &= ~ID_STAR_POOR;
      }
      Nave ++;
    }
  }
  fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
  initstats (STATMODE);
  free (xlist);
  free (slist);
  free (dlist);
}

void clean_measures (Catalog *catalog, int Ncatalog, int final) {

  int i, j, k, m, N, Ndel, Nave, Nmax, image_bad, TOOFEW;
  int *ilist;
  double *tlist, *list, *dlist, Ns;
  float Msys, Mcal, Mmos, Mgrid;
  StatType stats;

  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
  /* Nmeasure is now different, need to reallocate */
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].average[j].Nm);
    }
  }
  ALLOCATE (list, double, MAX (1, Nmax));
  ALLOCATE (dlist, double, MAX (1, Nmax));
  ALLOCATE (ilist, int, Nmax);
  ALLOCATE (tlist, double, Nmax);
  
  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
  TOOFEW = MAX (5, STAR_TOOFEW);

  Ns = 3;
  Ndel = Nave = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* skip bad stars to prevent them from becoming good (on inner sample) */
      if (catalog[i].average[j].code & STAR_BAD) continue;  

      /* on final processing, skip stars already measured */
      if (final && catalog[i].found[j]) continue;  

      /* accumulate list of valid measurements */
      m = catalog[i].average[j].offset;
      N = 0;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	/* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
	Mcal  = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) continue;
	Mmos  = getMmos  (m, i);
	if (Mmos == NO_MAG) continue;
	if (Mmos == NO_IMAGE) continue;
	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;

	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	list[N] = Msys - Mcal - Mmos - Mgrid;
	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
	N++;
      }
      if (N < TOOFEW) continue;

      /* 3-sigma clip based on stats of inner 50% */
      initstats ("INNER_MEAN");
      liststats (list, dlist, N, &stats);
      stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
      for (k = m = 0; k < N; k++) {
	if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
	  list[m] = list[k];
	  m++;
	}
      }
      initstats ("MEAN");
      liststats (list, dlist, m, &stats);
      stats.sigma = MAX (MIN_ERROR, stats.sigma);

      /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
      image_bad = IMAGE_BAD;
      IMAGE_BAD = ID_IMAGE_NOCAL;
      m = catalog[i].average[j].offset;
      N = 0;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	/* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
	Mcal  = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) continue;
	Mmos  = getMmos  (m, i);
	if (Mmos == NO_MAG) continue;
	if (Mmos == NO_IMAGE) continue;
	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;

	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	list[N] = Msys - Mcal - Mmos - Mgrid;
	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
	ilist[N] = m;
	N++;
	Nave ++;
      }
      if (N < TOOFEW) continue;

      /* mark bad measures */
      for (k = 0; k < N; k++) {
	if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
	  catalog[i].measure[ilist[k]].flags |= ID_MEAS_POOR;
	  Ndel ++;
	}
      }
      IMAGE_BAD = image_bad;
    }
  }
  initstats (STATMODE);
  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
  free (ilist);
  free (tlist);
}

StatType statsStarN (Catalog *catalog, int Ncatalog) {

  int i, j, k, m, n, N, Ntot;
  double *list, *dlist;
  float Mcal, Mmos;
  StatType stats;

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

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

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].average[j].code & STAR_BAD) continue;  
      m = catalog[i].average[j].offset;

      N = 0;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	Mcal = getMcal  (m, i);
	if (Mcal == NO_MAG) continue;
	if (Mcal == NO_IMAGE) continue;
	Mmos = getMmos  (m, i);
	if (Mmos == NO_MAG) continue;
	if (Mmos == NO_IMAGE) continue;
	if (getMgrid (m, i) == NO_MAG) continue;
	N++;
      }
      
      list[n] = N;
      dlist[n] = 1;
      n++;
    }
  }

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

StatType statsStarX (Catalog *catalog, int Ncatalog) {

  int i, j, n, Ntot, Xm;
  double *list, *dlist;
  StatType stats;

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

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

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].average[j].code & STAR_BAD) continue;  

      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
      if (Xm == NO_MAG) continue;
      list[n] = pow (10.0, 0.01*Xm);
      dlist[n] = 1;
      n++;
    }
  }

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

StatType statsStarS (Catalog *catalog, int Ncatalog) {

  int i, j, n, Ntot;
  double *list, *dlist;
  float dM;
  StatType stats;

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

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

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].average[j].code & STAR_BAD) continue;  

      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
      list[n] = dM;
      dlist[n] = 1;
      n++;
    }
  }

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

void plot_stars (Catalog *catalog, int Ncatalog) {

  int i, j, bin;
  float dMrel;
  double *xlist, *Mlist;
  Graphdata graphdata;

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

  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      if (catalog[i].average[j].code & STAR_BAD) continue;  
      dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
      bin = dMrel / 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, "dMrel hist", "dMhist.png");

  free (xlist);
  free (Mlist);
}

void plot_chisq (Catalog *catalog, int Ncatalog) {

  int i, j, N, Ntotal, value;
  double *xlist, *ylist;
  Graphdata graphdata;

  Ntotal = 0;
  for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage;

  ALLOCATE (xlist, double, Ntotal);
  ALLOCATE (ylist, double, Ntotal);

  N = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      if (catalog[i].average[j].code & STAR_BAD) continue;
      xlist[N] = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
      value    = PhotPrimary ? catalog[i].average[j].Xm   : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
      if (value == NO_MAG) continue;
      ylist[N] = 0.01*value;
      N++;
    }
  }

  plot_defaults (&graphdata);
  graphdata.ymin = -3.0;
  plot_list (&graphdata, xlist, ylist, N, "chisq", NULL);
  free (xlist);
  free (ylist);
}

void plot_star_coords (Catalog *catalog, int Ncatalog) {

  int i, j, N;
  double *xlist, *ylist;
  double Xmin, Ymin, Xmax, Ymax;
  Graphdata graphdata;

  N = 0; 
  for (i = 0; i < Ncatalog; i++) {
    N += catalog[i].Naverage;
  }
  ALLOCATE (xlist, double, N);
  ALLOCATE (ylist, double, N);

  N = 0;
  Xmin = Ymin = +360.0;
  Xmax = Ymax = -360.0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      xlist[N] = catalog[i].average[j].R;
      ylist[N] = catalog[i].average[j].D;
      N++;
    }
  }
  plot_defaults (&graphdata);
  plot_list (&graphdata, xlist, ylist, N, "coords", NULL);

  free (xlist);
  free (ylist);
}

