#include <stdio.h>
#include <string.h>
#include <curses.h>
#include "structs.h"
#include "nrutil.h"

#define DONE 1

struct image *mask;      /* stamp-sized region of bad pixel mask */
struct image model, cmodel, tpsf, sincwin;
struct derivs *df;   /* spaces reserved for derivative images */
struct inpars input;
float **y1a, **y2a, **y12a;
float xskycent, yskycent;
struct bicubic_grid **cg;

struct image *sigma_img (struct image *data, struct image *mask);
float sigcheck (struct image *sigma);
void count_pars (struct fitpars *fpar, int *nfree, int *nfixed, int *unused);
double LevMar (struct image *d, struct image *n, struct image *psf, 
    struct fitpars *fpar, struct cons *constr, struct convpars *cpar, 
    double **covar, double *sigma, int *ndof, int type);
struct bicubic_grid **get_cgrid (struct image *psf);
void fprintpar (int lowx, int lowy, struct fitpars *par, FILE *ftype, 
    char galfitlog[], double *sigma, float *coerr, double chisq, int ndof);
void delete_work_arrays (struct fitpars *fpar, struct derivs *);
struct derivs *create_work_arrays (struct fitpars *fpar, struct image *model,
    int convolve, struct convpars *cpar);
void point_fpar_at_df (struct fitpars *fpar, struct derivs *df);
void strong_links (struct fitpars *fpar, struct derivs *df, 
    struct cons *constr);
int geterr (double *chi2, int ndof, struct image *d, struct image *n, struct
    image *psf, struct fitpars *fpar, struct convpars *cpar, float
    *coerr, double *sigma);
void outmenu (struct fitpars *fpar, char *galfitlog, float xoffset, 
    float yoffset, float chisq, int ndof);
void outmodel (struct fitpars *fpar, struct convpars *cpar,
    struct image *d, struct image *model, struct image *sigma, 
    struct derivs *df, double *chisq, int *ndof, int nfree, int nfixed, 
    int offx, int offy);
int (*pfunc) (const char *, ...);
void free_cgrid (struct bicubic_grid **m, long nrl, long nrh, long ncl, 
    long nch);
void assign_err (struct fitpars *fpar, double *sigma, float chi2nu);
void outfile (char galfitlog[]);

void search_min (struct image *data, struct image *sigma, struct image *psf,
    struct image *badpix, struct fitpars *fpar, struct cons *constr,
    struct convpars *cpar) 
{
    void offset_center (struct fitpars *fpar, struct convpars *cpar, int xoff,
							int yoff);
    struct image *to_mini (struct image *);

    FILE *logfile;
    struct image *d, *s;  /* sub-region of data and noise image to fit */
    int i, ndof=0, totpars=0, nfree=0, nfixed=0, unused=0;
    int lowx=0, lowy=0, highx=0, highy=0;
    double **covar, *sig;
    float *coerr, sigval;
    double chisq=0.;
    long tnaxes[3];
    struct fitpars *fptr;
    extern struct sampars sample;
    extern char device[];

    /* When doing more than one object in an image, can modify this *\
    \* subroutine to loop through all the objects                   */

 /************************************************************************\
 *                                                                        *
 *    Set up all the image matrices for fitting once and for all          *
 *                                                                        *
 \************************************************************************/

    /*-----------------------------------------------------------------\
    |  First cut the section of the image that we want to fit from     |
    |  the original data/sigma images into new "mini-images."  In      |
    |  general the mini-images will be larger than the convolution     |
    |  region.  The mini-images are regions the user specified         |
    |  to fit with no convolution padding added.                       |
    \-----------------------------------------------------------------*/

    d = to_mini (data);

    free_matrix (data->z, 1, data->naxes[2], 1, data->naxes[1]);

    mask = to_mini (badpix);

    if (sigma->err == 0)
        s = to_mini (sigma);
    else
	s = sigma_img (d, mask);

    if (sigval = sigcheck(s))
        pfunc ("-- WARNING: Sigma image has 0. values.  Set to a low value of %.2e. \n", sigval);

    /*-------------------------------------------------\
    |  Offset the object centers and the convolution   |
    |  box relative to the stamp-sized image           |
    \-------------------------------------------------*/

    sscanf (d->imgsect, "[%d:%d,%d:%d]", &lowx, &highx, &lowy, &highy);
    offset_center (fpar, cpar, -(lowx-1), -(lowy-1));

    xskycent = (d->naxes[1] + 1)/ 2.;
    yskycent = (d->naxes[2] + 1)/ 2.;

    /*------------------------------------------\
    |  Prepare for bicubic interpolation of PSF |
    \------------------------------------------*/

    if (strncmp (psf->name, "none", 4) != 0) { 
        tpsf.naxes[1] = psf->naxes[1];
        tpsf.naxes[2] = psf->naxes[2];

	tpsf.z = matrix (1, tpsf.naxes[2], 1, tpsf.naxes[1]);
	cg = get_cgrid (psf);
    }

    /*------------------------------------------------------------\
    |  Create storage space for the derivative and model images.  |
    |  Store them as global.  These are real (as opposed to       |
    |  complex).  They're the same size as the postage sized      |
    |  image and will NOT be 2^N in size and will not have any    |
    |  padding.  "model" and "df" will have the size of the       |
    |  entire fitting region.                                     |
    \------------------------------------------------------------*/

    model = *d;
    df = create_work_arrays (fpar, &model, 0, cpar);

    /*---------------------------------------------------------------------\
    |  Now point fpar->dp at the corresponding component in the df         |
    |  structure.  This is quite useful when you need to know what         |
    |  derivative structure corresponds to which set of parameters on the  |
    |  fly.  This happens often when dealing with truncation parameters.   |
    \---------------------------------------------------------------------*/

    point_fpar_at_df (fpar, df);

    if (model.naxes[2] > 0 && model.naxes[1] > 0)
        model.z = matrix (1, model.naxes[2], 1, model.naxes[1]);


 /************************************************\
 *  Create covariance matrix and then minimize!   *
 \************************************************/

    count_pars (fpar, &nfree, &nfixed, &unused);
    totpars = nfree + nfixed + unused;

    covar = dmatrix (1, totpars, 1, totpars);
    coerr = vector (1, totpars);
    sig = dvector (1, totpars);

    if (!input.create) {
        chisq = LevMar (d, s, psf, fpar, constr, cpar, covar, sig, &ndof, 1);

        /*********\
        *  DONE!   \*************************************************\
        *  Now create model, output the data, model, residual,       *
        *  images into a 3-layered block.                            *
        \************************************************************/

        assign_err (fpar, sig, chisq/ndof);
        pfunc ("\n\nFit summary is now being saved into `fit.log'.\n\n");
    } else
        pfunc ("\n\n-- Creating output image....\n\n");


    if (strncmp (device, "curses", 6) == 0)
        refresh();

    outfile (d->logfile);
    strcpy (model.logfile, d->logfile);

    outmodel (fpar, cpar, d, &model, s, df, &chisq, &ndof, nfree, nfixed,
	(lowx-1), (lowy-1));

 /*********************************************************\
 *  Restore x and y offsets, then output the initial and   *
 *  final parameters into files.                           *
 \*********************************************************/

    xskycent += (lowx-1);
    yskycent += (lowy-1);
    offset_center (fpar, cpar, (lowx-1), (lowy-1));

    if (! input.create) {
        outmenu (fpar, d->logfile, 0., 0., chisq, ndof);

        if ( (logfile = fopen ("fit.log", "a")) == (FILE *) 0)
            logfile = fopen ("fit.log", "w");
        fprintpar (1, 1, fpar, logfile, d->logfile, sig, coerr, chisq, ndof);
        fclose (logfile);
    };

 /***********************************************\
 *  Now free up all the memory grabbed up.       *
 \***********************************************/

    delete_work_arrays (fpar, df);

    if (model.naxes[2] > 0 && model.naxes[1] > 0)
        free_matrix (model.z, 1, model.naxes[2], 1, model.naxes[1]);
    free_vector (coerr, 1, totpars);
    free_dmatrix (covar, 1, totpars, 1, totpars);
    free_matrix (mask->z, 1, mask->naxes[2], 1, mask->naxes[1]);
    free_matrix (d->z, 1, d->naxes[2], 1, d->naxes[1]);
    free_matrix (s->z, 1, s->naxes[2], 1, s->naxes[1]);
    free_dvector (sig, 1, totpars);
    if (strncmp (psf->name, "none", 4) != 0)  
        free_matrix (tpsf.z, 1, tpsf.naxes[2], 1, tpsf.naxes[1]);
    if (strncmp (psf->name, "none", 4) != 0)
        free_cgrid (cg, 0, psf->naxes[2]+1, 0, psf->naxes[1]+1);

    free (mask);
    free (s);
    free (d);

    return;
}


/*****************************************************************************/


struct image *to_mini (struct image *img)
{
    extern int (*pfunc)(const char *, ...);
    extern char device[];
    int i, j, xmin=0, xmax=0, ymin=0, ymax=0;
    struct image *new;

    new = (struct image *) calloc (1, (size_t)(sizeof (struct image)));
    if (new == (struct image *) 0) {
	pfunc ("Error in allocating memory for mini-image....");
        if (strncmp(device, "curses", 6)== 0)
            refresh();
	exit (1);
    };

    *new = *img;
    sscanf (img->imgsect, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
    new->naxes[1] = xmax - xmin + 1;
    new->naxes[2] = ymax - ymin + 1;

    new->z = matrix(1, new->naxes[2], 1, new->naxes[1]);

    xmax = IMIN (xmax, img->naxes[1]);   /*  In case the image size is      */
    ymax = IMIN (ymax, img->naxes[2]);   /*  smaller than the fitting size  */
					 /*  for whatever reason....        */

    for (j=ymin; j <= ymax; ++j){
	for (i=xmin; i <= xmax; ++i)
	    new->z[j-ymin+1][i-xmin+1] = img->z[j][i];
    }; 
    return (new);
}


/*****************************************************************************/

void offset_center (struct fitpars *fpar, struct convpars *cpar, int xoff,
							int yoff)
{
    while (fpar != NULL) {
        if (strncmp (fpar->objtype, "sky", 3) != 0) {
            fpar->a[1] += xoff;
            fpar->a[2] += yoff;
        };
        fpar = fpar->next;
    };
    cpar->cent[1] += xoff;
    cpar->cent[2] += yoff;
}