/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

#include <math.h>
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <unistd.h>
#include <stdint.h>
#include <limits.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <assert.h>

#include "starutil.h"
#include "codefile.h"
#include "mathutil.h"
#include "quadfile.h"
#include "kdtree.h"
#include "fitsioutils.h"
#include "anqfits.h"
#include "starkd.h"
#include "boilerplate.h"
#include "errors.h"
#include "log.h"
#include "quad-utils.h"
#include "quad-builder.h"
#include "allquads.h"

allquads_t* allquads_init() {
    allquads_t* aq = calloc(1, sizeof(allquads_t));
    aq->dimquads = 4;
    return aq;
}

void allquads_free(allquads_t* aq) {
    free(aq);
}

static void add_quad(quadbuilder_t* qb, unsigned int* quad, void* token) {
    allquads_t* aq = token;
    if (log_get_level() > LOG_VERB) {
        int k;
        debug("quad: ");
        for (k=0; k<qb->dimquads; k++)
            debug("%-6i ", quad[k]);
        logverb("\n");
    }
    quad_write_const(aq->codes, aq->quads, quad, aq->starkd, qb->dimquads, aq->dimcodes);
}

static anbool check_AB(quadbuilder_t* qb, pquad_t* pq, void* token) {
    allquads_t* aq = token;
    debug("check_AB: iA=%i, iB=%i, aq->starA=%i\n", pq->iA, pq->iB, aq->starA);
    return (pq->iA == aq->starA);
}

/*
 #include "codefile.h"
 static bool check_full_quad(quadbuilder_t* qb, unsigned int* quad, int nstars, void* token) {
 allquads_t* aq = token;
 double code[4];
 double xyz[12];
 double R1, R2;
 startree_get(aq->starkd, quad[0], xyz+0);
 startree_get(aq->starkd, quad[1], xyz+3);
 startree_get(aq->starkd, quad[2], xyz+6);
 startree_get(aq->starkd, quad[3], xyz+9);
 codefile_compute_star_code(xyz, code, 4);
 logmsg("Code: %g, %g, %g, %g\n", code[0], code[1], code[2], code[3]);
 R1 = hypot(code[0] - 0.5, code[1] - 0.5);
 R2 = hypot(code[2] - 0.5, code[3] - 0.5);
 logmsg("Radii %g, %g\n", R1, R2);
 return TRUE;
 }
 */

int allquads_create_quads(allquads_t* aq) {
    quadbuilder_t* qb;
    int i, N;
    double* xyz;

    qb = quadbuilder_init();

    qb->quadd2_high = aq->quad_d2_upper;
    qb->quadd2_low  = aq->quad_d2_lower;
    qb->check_scale_high = aq->use_d2_upper;
    qb->check_scale_low  = aq->use_d2_lower;

    qb->dimquads = aq->dimquads;
    qb->add_quad = add_quad;
    qb->add_quad_token = aq;

    N = startree_N(aq->starkd);

    logverb("check scale high: %i\n", qb->check_scale_high);

    if (!qb->check_scale_high) {
        int* inds;
        inds = malloc(N * sizeof(int));
        for (i=0; i<N; i++)
            inds[i] = i;
        xyz = malloc(3 * N * sizeof(double));
        kdtree_copy_data_double(aq->starkd->tree, 0, N, xyz);

        qb->starxyz = xyz;
        qb->starinds = inds;
        qb->Nstars = N;

        quadbuilder_create(qb);

        free(xyz);
        free(inds);
    } else {
        int nq;
        int lastgrass = 0;

        /*
         xyz = malloc(3 * N * sizeof(double));
         kdtree_copy_data_double(aq->starkd->tree, 0, N, xyz);
         */

        // star A = i
        nq = aq->quads->numquads;
        for (i=0; i<N; i++) {
            double xyzA[3];
            int* inds;
            int NR;

            int grass = (i*80 / N);
            if (grass != lastgrass) {
                printf(".");
                fflush(stdout);
                lastgrass = grass;
            }

            startree_get(aq->starkd, i, xyzA);

            startree_search_for(aq->starkd, xyzA, aq->quad_d2_upper,
                                &xyz, NULL, &inds, &NR);

            /*
             startree_search_for(aq->starkd, xyzA, aq->quad_d2_upper,
             NULL, NULL, &inds, &NR);
             */

            logverb("Star %i of %i: found %i stars in range\n", i+1, N, NR);
            aq->starA = i;
            qb->starxyz = xyz;
            qb->starinds = inds;
            qb->Nstars = NR;
            qb->check_AB_stars = check_AB;
            qb->check_AB_stars_token = aq;
            //qb->check_full_quad = check_full_quad;
            //qb->check_full_quad_token = aq;

            quadbuilder_create(qb);

            logverb("Star %i of %i: wrote %i quads for this star, total %i so far.\n", i+1, N, aq->quads->numquads - nq, aq->quads->numquads);
            free(inds);
            free(xyz);
        }
        //
        //free(xyz);

        printf("\n");
    }

    quadbuilder_free(qb);
    return 0;
}

int allquads_close(allquads_t* aq) {
    startree_close(aq->starkd);

    // fix output file headers.
    if (quadfile_fix_header(aq->quads) ||
        quadfile_close(aq->quads)) {
        ERROR("Couldn't write quad output file");
        return -1;
    }
    if (codefile_fix_header(aq->codes) ||
        codefile_close(aq->codes)) {
        ERROR("Couldn't write code output file");
        return -1;
    }
    return 0;
}

int allquads_open_outputs(allquads_t* aq) {
    int hp, hpnside;
    qfits_header* hdr;

    printf("Reading star kdtree %s ...\n", aq->skdtfn);
    aq->starkd = startree_open(aq->skdtfn);
    if (!aq->starkd) {
        ERROR("Failed to open star kdtree %s\n", aq->skdtfn);
        return -1;
    }
    printf("Star tree contains %i objects.\n", startree_N(aq->starkd));

    printf("Will write to quad file %s and code file %s\n", aq->quadfn, aq->codefn);
    aq->quads = quadfile_open_for_writing(aq->quadfn);
    if (!aq->quads) {
        ERROR("Couldn't open file %s to write quads.\n", aq->quadfn);
        return -1;
    }
    aq->codes = codefile_open_for_writing(aq->codefn);
    if (!aq->codes) {
        ERROR("Couldn't open file %s to write codes.\n", aq->quadfn);
        return -1;
    }

    aq->quads->dimquads = aq->dimquads;
    aq->codes->dimcodes = aq->dimcodes;

    if (aq->id) {
        aq->quads->indexid = aq->id;
        aq->codes->indexid = aq->id;
    }

    // get the "HEALPIX" header from the skdt and put it in the code and quad headers.
    hp = qfits_header_getint(startree_header(aq->starkd), "HEALPIX", -1);
    if (hp == -1) {
        logmsg("Warning: skdt does not contain \"HEALPIX\" header.  Code and quad files will not contain this header either.\n");
    }
    aq->quads->healpix = hp;
    aq->codes->healpix = hp;
    // likewise "HPNSIDE"
    hpnside = qfits_header_getint(startree_header(aq->starkd), "HPNSIDE", 1);
    aq->quads->hpnside = hpnside;
    aq->codes->hpnside = hpnside;

    hdr = quadfile_get_header(aq->quads);
    qfits_header_add(hdr, "CXDX", "T", "All codes have the property cx<=dx.", NULL);
    qfits_header_add(hdr, "CXDXLT1", "T", "All codes have the property cx+dx<=1.", NULL);
    qfits_header_add(hdr, "MIDHALF", "T", "All codes have the property cx+dx<=1.", NULL);
    qfits_header_add(hdr, "CIRCLE", "T", "Codes live in the circle, not the box.", NULL);

    hdr = codefile_get_header(aq->codes);
    qfits_header_add(hdr, "CXDX", "T", "All codes have the property cx<=dx.", NULL);
    qfits_header_add(hdr, "CXDXLT1", "T", "All codes have the property cx+dx<=1.", NULL);
    qfits_header_add(hdr, "MIDHALF", "T", "All codes have the property cx+dx<=1.", NULL);
    qfits_header_add(hdr, "CIRCLE", "T", "Codes live in the circle, not the box.", NULL);

    if (quadfile_write_header(aq->quads)) {
        ERROR("Couldn't write headers to quads file %s\n", aq->quadfn);
        return -1;
    }
    if (codefile_write_header(aq->codes)) {
        ERROR("Couldn't write headers to code file %s\n", aq->codefn);
        return -1;
    }

    if (!aq->use_d2_lower)
        aq->quad_d2_lower = 0.0;
    if (!aq->use_d2_upper)
        aq->quad_d2_upper = 10.0;

    aq->codes->numstars = startree_N(aq->starkd);
    aq->codes->index_scale_upper = distsq2rad(aq->quad_d2_upper);
    aq->codes->index_scale_lower = distsq2rad(aq->quad_d2_lower);

    aq->quads->numstars = aq->codes->numstars;
    aq->quads->index_scale_upper = aq->codes->index_scale_upper;
    aq->quads->index_scale_lower = aq->codes->index_scale_lower;
    return 0;
}