# include "dvosplit.h"
# define NROWS 100000 /* ~10MB per block for measures */
# define DNOUT 1000

AveLinks *split_averages (Catalog *incatalog, SkyList *outlist, Catalog *outcatalogs) {

  double inR, inD;
  int n, block, ave, cat, averef, Nblocks, Ncat, Nout, Nsecfilt;
  int *outref, *outcat, *outmem;
  AveLinks *avelinks;

  ALLOCATE (outref, int, incatalog[0].Naves_disk);
  ALLOCATE (outcat, int, incatalog[0].Naves_disk);
  ALLOCATE (outmem, int, outlist[0].Nregions);

  Nsecfilt = GetPhotcodeNsecfilt ();

  // allocate enough space for these output buffers: use Nsecfilt + 1 incase the file
  // contains primary photcodes, which will increase Nsecfilt by one.
  for (cat = 0; cat < outlist[0].Nregions; cat++) {
    outmem[cat] = DNOUT;
    REALLOCATE (outcatalogs[cat].average, Average, outmem[cat]);
    REALLOCATE (outcatalogs[cat].secfilt, SecFilt, outmem[cat]*(Nsecfilt + 1));
  }

  // split out the average & secfilt entries:
  incatalog[0].catflags = LOAD_AVES | LOAD_SECF;
  Nblocks = incatalog[0].Naves_disk / NROWS;
  if (incatalog[0].Naves_disk % NROWS) Nblocks ++;
  for (block = 0; block < Nblocks; block++) {

    // read up to NROWS at a time
    dvo_catalog_load_segment (incatalog, VERBOSE, block*NROWS, NROWS);
    fprintf (stderr, "splitting %s (averages) .. %d of %d\n", incatalog[0].filename, block, Nblocks);

    assert (block*NROWS == incatalog[0].Naves_off);

    // distribute data to the output catalogs
    for (ave = 0; ave < incatalog[0].Naverage; ave++) {
      averef = ave + incatalog[0].Naves_off;
	
      inR = incatalog[0].average[ave].R;
      inD = incatalog[0].average[ave].D;

      // which of the outcatalogs contains this coordinate?

      Ncat = -1;
      for (cat = 0; cat < outlist[0].Nregions; cat++) {
	if (inR < outlist[0].regions[cat][0].Rmin) continue;
	if (inR > outlist[0].regions[cat][0].Rmax) continue;
	if (inD < outlist[0].regions[cat][0].Dmin) continue;
	if (inD > outlist[0].regions[cat][0].Dmax) continue;
	Ncat = cat;
	break;
      }

      if (Ncat == -1) {
	fprintf (stderr, "WARNING: missed %d (%f, %f)\n", averef, inR, inD);
	continue;
      }

      Nout = outcatalogs[Ncat].Naverage;
      outref[averef] = Nout + outcatalogs[Ncat].Naves_off;
      outcat[averef] = Ncat;

      // assign the value to the next element of the output catalog
      outcatalogs[Ncat].average[Nout] = incatalog[0].average[ave];
      outcatalogs[Ncat].Naverage ++;

      // update secfilt at the same time
      for (n = 0; n < Nsecfilt; n++) {
	outcatalogs[Ncat].secfilt[Nout*Nsecfilt + n] = incatalog[0].secfilt[ave*Nsecfilt + n];
	outcatalogs[Ncat].Nsecf_mem++;
      }

      if (outcatalogs[Ncat].Naverage >= outmem[Ncat]) {
	outmem[Ncat] += DNOUT;
	REALLOCATE (outcatalogs[Ncat].average, Average, outmem[Ncat]);
	REALLOCATE (outcatalogs[Ncat].secfilt, SecFilt, outmem[Ncat]*(Nsecfilt + 1));
      }
    }
    dvo_catalog_free_data (incatalog);

    // double check the values of Naverage, Nsecf_mem?

    // write out the new values
    for (cat = 0; cat < outlist[0].Nregions; cat++) {
      outcatalogs[cat].catflags = LOAD_AVES | LOAD_SECF;

      dvo_catalog_save (&outcatalogs[cat], VERBOSE);
      // fprintf (stderr, "secfilt: %d %d %d %d\n", outcatalogs[cat].Nsecf_mem, outcatalogs[cat].Nsecf_disk, outcatalogs[cat].Nsecf_off, outcatalogs[cat].Naverage, outcatalogs[cat].Nsecfilt);

      // advance the pointers and free the current data
      // XXX these should be done within save segment:
      outcatalogs[cat].Naves_disk += outcatalogs[cat].Naverage;
      outcatalogs[cat].Naves_off  += outcatalogs[cat].Naverage;
      outcatalogs[cat].Nsecf_disk += outcatalogs[cat].Nsecfilt * outcatalogs[cat].Naverage;
      outcatalogs[cat].Nsecf_off  += outcatalogs[cat].Nsecfilt * outcatalogs[cat].Naverage;
      outcatalogs[cat].Nsecfilt    = Nsecfilt;

      outcatalogs[cat].Naverage    = 0;
      outcatalogs[cat].Nsecf_mem   = 0;
    }
  }

  for (cat = 0; cat < outlist[0].Nregions; cat++) {
    dvo_catalog_free_data (&outcatalogs[cat]);
  }

  free (outmem);

  ALLOCATE (avelinks, AveLinks, 1);
  avelinks[0].outref = outref;
  avelinks[0].outcat = outcat;


  return (avelinks);
}
