struct covfunconfig {

  /* The possible covariance kernels are:
     1 - time&space squared error,
     2 - space only,
     3 - periodic.
     4 - Matérn with nu=5/2.
     5 - windfield-informed kernel
  */
  int covftype;

  /* covariance function arguments and their inverses (squared for l)
     - not all of these are set for all types of kernels */
  float tau, llat, llon, lt, lperiodic, invllat2, invllon2, invlt, invlperiodic2, rho;

  /* Lower and higher limits for those parameters */
  float tau_low, llat_low, llon_low, lt_low, lperiodic_low, rho_low;
  float tau_high, llat_high, llon_high, lt_high, lperiodic_high, rho_high;

  /* nmaxobs is the max number of obs that we consider - we thin data
     to this number if there are more. */
  int nmaxobs;

  /* This is only for the Matérn kernels; to precompute the b-factors */
  float *bfact;

  /* This is the list of pointers to the kernels for a compound kernel */
  struct covfunconfig **kernels;
};


struct config {
  /* Restrict area to these coordinates */
  double minlat, maxlat;
  double minlon, maxlon;

  /* Resolution of the grid */
  double gridres;

  /* meanfun and compound covfunconfig */
  struct mfstruct *mfs;
  struct covfunconfig *cfc;

  /* Number of lats, lons, days, and total grid points */
  size_t nlat, nlon, maxdays, ngp;

  /* Arrays for the latitude values of the grid points */
  double *lat, *lon;

  /* Arrays for the cartesian unit ball locations for the
     longitude-latitude points. */
  double *x, *y, *z;

  /* For training data, if we do not do a gridded GP, we also save the
     time, dv, and dv_unc variables in the config. For the regular
     GP, these will be set to NULL. */
  double *t, *dv, *dv_unc, *u, *v;

  /* Minimum correlation under which observations are no longer
     included; rad is the distance on unit ball surface that this
     corresponds to. max_rad_equator_degrees is the same, but in
     degrees when you travel the latitudes, or equivalently, on the
     equator. max_tdist refers to the amount of days beyond which the
     correlation would be too weak for observations to be any more
     included. */
  double min_cov, max_rad, max_rad_equator_degrees, max_tdist;

  /* Time in seconds from Jan 1st 1970 */
  double first_day_noon_unix;

  /* The amount of memory we allocate when we need more (this linear
     scheme is ok for this code since we don't allocate more in
     practice often */
  size_t malloc_chunk;

  /* nmaxobs is the maximum size of the K-matrix to be considered for
     Cholesky - not counting the point where the GP is evaluated */
  size_t nmaxobs;

  /* Grid index increment for each sweep of data. Done so that we can
     calculate large grids and not run out of memory*/
  int gi_inc;

  /* Mean distance between observations to be added (or proportional
     to it. To include everything, make this very tiny. To have sparse
     observations, have a big one here. These are unit ball
     distances. Only a random fraction of observations is added for
     performance reasons */
  double obs_dist;

  /* mode is 1 for gridGP, 2 for find_parameters() */
  int mode;

  /* 0 for xco2, 1 for wind u component, and 2 for wind v component */
  int dependent_variable;

  /* use_wind: do we read and use wind information? */
  int use_wind;

  /* List of days for which GP will be run */
  int *daylist;
  int ndays_in_daylist;
};


struct daydata {
  /* Data array containing dv, dv_unc, and time */
  float **data;

  /* Data arrays for latitudes and longitudes */
  float *lat, *lon;

  /* Number of observations in the grid added */
  size_t ngobs;

  /* The time (in seconds since 1970) of the first observation in
     file. This is used for deducing the current day */
  double first_datatime;

  /* Set to 3 (dv, dv_unc, time) if no wind is used, 5 if winds
     are used. */
  int nvars;
};


struct cholstruct {
  /* This struct contains the covariance matrix data (nobs^2 vector),
     number of observations, and the vector with which the predicted
     observation can be recovered */
  size_t nobs, max_tot_ksize;

  /* Prior variance for predicted point */
  float prior_var;

  /* The nobs*nobs between-observations covariance matrix, as an array
     as required by lapacke */
  double *K;

  /* This is the observations vector, of size nobs */
  double *dv;

  /* This vector contains covariances between the grid point where the
     Gaussian process is evaluated and the observations */
  double *Kxx;

  /* meanvalue is used for holding the mean values of the observations
     for each grid point */
  double meanvalue;

  /* x contains a helper array of size E.nmaxobs for each gridpoint
     to calculate the concentration in a thread-safe manner. */
  double *x;

  /* obsind gives the index of the observations in S->gridx[gi], so
     that accessing the location is possible afterwards. This was
     earlier needed for interpolating observations from dynamic to
     static GP, and can also be used for debugging */
  size_t *obsind;

  /* For recording what winds were used if covftype == 5 */
  float uwind, vwind;
};


struct staticGP {
  /* This struct is close to state, but we don't need to track
     time so this should be easier */

  /* dv, dv_uncertainty arrays of length E.ngp */
  float *result_dv, *result_dv_unc;

  /* For calculating the cholesky decompositions at each point */
  struct cholstruct **C;
};


struct state {
  /* For each point in grid, contains an array of Cartesian
     coordinates of the observations within radius for each gridpoint

     Difference between state and config is that data in state (S in
     code) may change but config data does not, once experiment has
     started. */
  
  /* Arrays containing anything that has to do with the observations
     that we want to deal with. */
  float *x;
  float *y;
  float *z;
  float *t;
  float *dv;
  float *dv_unc;

  /* Wind fields */
  float *u;
  float *v;

  /* How many obs we had per day */
  size_t **griddaycounts;

  /* Same as E->maxdays */
  int total_days;

  /* Increment indexes for doing part of the data at a time */
  size_t gi0, gi1;

  /* Number of observation that we have in S, so length of S->x etc.*/
  size_t nobs;

  /* These are used for bookkeeping how many failures we got. */
  size_t allocated_for_all_obs, nancount, total_choleskys;

  /* Current unix time is used for observation correlation calclations */
  float first_day_noon;

  /* How much memory is allocated for each grid point on each day*/
  size_t **allocated;

  /* What indexes in arrays x,y,z etc are close enough for day d and
     grid index gi? FIXME This changed: Each grid index and day has
     now a set of indexes that refer to observations on that day in
     that location. */
  size_t ***gridobsind;

  /* List of indexes for each day about which grid indexes have at
     least one observation. This is used to facilitate reading, but
     does not work with sampling */
  size_t **obs_grid_mask;
  size_t *obs_grid_mask_counts;

  /* Value to override E->min_cov for sampling - see comment in
     initialize_state for more */
  float min_cov_override;

  /* These are the gridded wind data, if winds from e.g. a previous GP
     are used. */
  float *uwinds_gridded;
  float *vwinds_gridded;
};


void initialize_cholstruct(struct cholstruct *C, size_t nobs) {
  /* This function just mallocs according to the number of
     observations given */

  C->nobs = nobs;
  if (!nobs) {
    C->K = NULL;
    C->Kxx = NULL;
    C->dv = NULL;
    C->x = NULL;
  } else {
    C->K = calloc(C->nobs*C->nobs, sizeof(double));
    C->Kxx = calloc(C->nobs, sizeof(double));
    C->dv = calloc(C->nobs, sizeof(double));
    C->x = calloc(C->nobs, sizeof(double));
    C->obsind = calloc(C->nobs, sizeof(size_t));
  }
  C->uwind = 0;
  C->vwind = 0;
}

struct state* initialize_state(struct state *S, struct config *E, size_t gi0, size_t gi1, float min_cov_override) {

  size_t i, j;

  S = malloc(sizeof(struct state));

  /* Allocate memory for the arrays holding the data, and for the
     arrays holding each day's datapoints for each information about
     data points that are close to these points


     Arguments:

     - *S is the state pointer that'll contain the initialized state
     - *E is the config with respect to which the state is created
     - gi0 and gi1 give the start and end of the grid points that
       belong to this sweep

     - min_cov_override: If we draw samples from the posterior, it
       becomes very memory intensive easily since all sampled points
       are added to list of points nearby. Additionally, due to
       Gaussian process's screening property, only nearby observations
       matter in the end. That's why E->max_rad is impractical, and we
       use the float given in sampling_mode to override that.  If 0,
       then E->min_cov is honored.
  */

  /* In the new one S->gridx[gi][somedaydata...] would be
     S->x[gi][day...] but starting point needs not be found by
     counting previous days, and otoh, many days need to be iterated
     over */

  S->min_cov_override = min_cov_override;

  S->x = NULL;
  S->y = NULL;
  S->z = NULL;
  S->t = NULL;
  S->dv = NULL;
  S->dv_unc = NULL;
  S->u = NULL;
  S->v = NULL;

  S->uwinds_gridded = NULL;
  S->vwinds_gridded = NULL;

  S->gridobsind = malloc(E->gi_inc * sizeof(size_t*));
  S->griddaycounts = malloc(E->gi_inc * sizeof(size_t*));
  S->allocated = malloc(E->gi_inc * sizeof(size_t*));

  S->obs_grid_mask = malloc(E->maxdays * sizeof(size_t*));
  S->obs_grid_mask_counts = malloc(E->maxdays * sizeof(size_t*));

  for (i=0; i<E->gi_inc; i++) {
    /* Array for each grid index in current batch */
    S->gridobsind[i] = calloc(E->maxdays, sizeof(size_t*));
    S->griddaycounts[i] = calloc(E->maxdays, sizeof(size_t));
    S->allocated[i] = calloc(E->maxdays, sizeof(size_t));
    for (j=0; j<E->maxdays; j++) {
      S->gridobsind[i][j] = NULL;
      S->griddaycounts[i][j] = 0;
      S->allocated[i][j] = 0;
    }
  }

  for (j=0; j<E->maxdays; j++) {
    S->obs_grid_mask[j] = calloc(E->ngp, sizeof(size_t));
    S->obs_grid_mask_counts[j] = 0;
  }


  /* Initialize the rest of the variables */
  S->total_days = E->maxdays;
  S->gi0 = gi0;
  S->gi1 = gi1;
  S->nancount = 0;
  S->total_choleskys = 0;
  S->allocated_for_all_obs = 0;
  S->nobs = 0;
  /* The first day's noon is in unix time given below. The way to get
     these is with the date-command in terminal; e.g.
     
     date --date="2014-09-06T12:00:00-00:00" +%s
     
     To prevent overflow, we subtract 1,400,000,000 from these
     values. */
  S->first_day_noon = E->first_day_noon_unix - 1400000000.;

  printf("time at first day's noon set to %f\n", S->first_day_noon);
  fflush(stdout);

  return S;
}

void rebuild_S_obs_grid_mask(struct state *S, struct config *E) {
  int i, j;

  /* First empty the whole thing */
  for (j=0; j<E->maxdays; j++) {
    S->obs_grid_mask_counts[j] = 0;
    for (i=0; i<E->ngp; i++) {
      S->obs_grid_mask[j][i] = 0;
    }
  }

  for (i=0; i<E->ngp; i++) {
    for (j=0; j<E->maxdays; j++) {
      if (S->griddaycounts[i][j]) {
	S->obs_grid_mask[j][S->obs_grid_mask_counts[j]++] = i;
      }
    }
  }
}

void teardown_daydata(struct daydata *D) {
  size_t i;

  if (D->lat) {free(D->lat);}
  if (D->lon) {free(D->lon);}

  if (D->data) {
    for (i=0; i<D->nvars; i++) {
      if (D->data[i]) {
	free(D->data[i]);
      }
    }
    free(D->data);
  }

  free(D);
}

void teardown_state(struct state *S, struct config *E) {
  int i, j;

  printf("Failure rate calculating Cholesky decompositions: %f%% (%zu/%zu)\n",
	 100.*S->nancount/(float)S->total_choleskys, S->nancount, S->total_choleskys);

  /* Free allocated memory in the range where we calculated */
  for (i=0; i<E->gi_inc; i++) {
    if (S->allocated[i]) {
      for (j=0; j<E->maxdays; j++) {
	if (S->allocated[i][j]) {
	  free(S->gridobsind[i][j]);
	}
      }
      free(S->griddaycounts[i]);
      free(S->allocated[i]);
    }
  }
  for (i=0; i<E->maxdays; i++) {
    free(S->obs_grid_mask[i]);
  }
  free(S->obs_grid_mask);
  free(S->obs_grid_mask_counts);

  /* Free the pointers to the arrays etc. */
  free(S->x);
  free(S->y);
  free(S->z);
  free(S->t);
  //free(S->C);
  free(S->dv);
  free(S->dv_unc);
  if (E->use_wind) {
      free(S->u);
      free(S->v);
    }
  free(S);
}

void teardown_cholstruct(struct cholstruct *C) {

  /* Nothing should have been allocated if nobs=0 */
  if (!C->nobs) {return;}
  free(C->K);
  free(C->dv);
  free(C->Kxx);
  free(C->x);
  free(C->obsind);
  free(C);
}

void teardown_covfunconfig(struct covfunconfig *cfconf) {
  if (cfconf->kernels) {
    free(cfconf->kernels);
  }
  if (cfconf) {
    free(cfconf);
    cfconf = NULL;
  }
}

void teardown_config(struct config *E) {
  int i;

  free(E->x);
  free(E->y);
  free(E->z);
  free(E->lat);
  free(E->lon);
  free(E->t);
  /* If mode != 0, E has been used just as a list of points , and
     these are not allocated. */
  if (E->dv) {
    free(E->dv);
    free(E->dv_unc);
  }
  if (E->cfc) { /* Same thing for cfc */

    if (E->cfc->covftype < 0) {
      for (i=0; i<-E->cfc->covftype; i++) {
	if (E->cfc->kernels[i]) {
	  teardown_covfunconfig(E->cfc->kernels[i]);
	}
      }
    }
    teardown_covfunconfig(E->cfc);
  }
  free(E->mfs);
}

struct boxmullerstruct {
  float x;
  uint counter;
};

struct boxmullerstruct *initialize_boxmullerstruct() {
  struct boxmullerstruct *BM = malloc(sizeof(struct boxmullerstruct));
  BM->counter = 0;
  return BM;
}
