text:
#include <model.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

#define P_inv_gamma gsl_cdf_gamma_Pinv // invverse cumulative gamma dist
#define P_inc_gamma gsl_sf_gamma_inc_P // regularised incomplete Gamma function
#define lngamma gsl_sf_lngamma       // ln Gamma function

#define V0 9
#define NUMBER_SHEEP 30
#define NUM_DONORS 10
#define SHEEP_FIRST_TEST_WEEK 1

// month defined as 365/12 days
// week defined as 7 days
# define WEEKS_PER_MONTH ((365/12)/7)

#include "Texel_clean_data.h"

typedef struct {
  double a, b;
  double *lambda;
  double *Lambda;
  double **W4;
  double *W2;
} parameters_t;

double u(int, double*);
double lngexpan(double, double, int*);
double ln_inc_gamma_star(double, double);
double w(int, int, int, parameters_t*);
double Prob_seroconvert(int, parameters_t*);
extern double R0(double, double, double, int, double*);

void Variables(struct variables *V)
{
  Var(V, 0, "week");
  Var(V, 1, "field infected");
  Var(V, 2, "new positives");
  Var(V, 3, "cumulative new positives");
  Var(V, 4, "seroconversion rate");
  Var(V, 5, "number infectious");
  Var(V, 6, "number infectious");
  Var(V, 6, "housed infected");
  Var(V, 7, "infected");
  Var(V, 8, "infected");
  Var(V, V0, "log_p_pos");
}

void Parameters(int trt, int ind, double *p, struct parameterset *Q)
{
  Par(Q,  0, Real, "beta_{field}", Exponential, 0.0008);
  Par(Q,  1, Real, "beta_{housed}", Gamma, 2.0, 0.1);
  Par(Q,  2, Real, "a", Exponential, 3.2);
  Par(Q,  3, Real, "k1", Gamma, 2.0, 4.0);
  // Par(Q,  0, Real, "beta_{field}", Exponential, 0.0082/WEEKS_PER_MONTH);
  // Par(Q,  1, Real, "beta_{housed}", Gamma, 2.0, 0.01);
  // Par(Q,  2, Real, "mu", Gamma, 2.0, 1.5);
  // Par(Q,  3, Real, "sigma", Gamma, 2.0, 1.5);

  // Par(Q,  8, Real, "a", Der);
  // Par(Q,  9, Real, "b", Der);
  // Par(Q, 10, Real, "beta_{housed}/beta_{field}", Der);
  // Par(Q, 12, Real, "E[inf_{field}]", Der);
  // Par(Q, 13, Real, "E[inf_{housed}]", Der);
  // Par(Q, 14, Real, "E[inf]", Der);
  // Par(Q, 15, Real, "s_{low}", Der);
  // Par(Q, 16, Real, "s_{high}", Der);
}

int Model(const uint trt, const uint ind, const uint nweeks, double *p, double **V, TimePoints *Data)
{
  int i, j, m, *num_infected;
  double **p_seropos = V, *Einf_field, *Einf_housed;
  parameters_t pars;

  double beta_field = p[0]/WEEKS_PER_MONTH;
  double beta_housed = p[1]/WEEKS_PER_MONTH;
  double S_a = p[2];
  double S_b = S_a/p[3]/WEEKS_PER_MONTH;
  int seroconversion = lrint(p[3]);

  for (m = 0; m < nweeks; m++)
    V[m][0] = m;

  int *num_infectious = integervector(nweeks);
  pars.Lambda = doublevector(nweeks);
  pars.lambda = doublevector(nweeks);
  pars.W4 = doublematrix(nweeks, nweeks);
  pars.W2 = doublevector(nweeks);
  pars.a = S_a;
  pars.b = S_b;

  // ewe to ewe transmission
  // the 10 donor sheep are infectious at start of experiment
  for (j = 0; j < NUM_DONORS; j++)
    for (m = 0; m <= sheep_last_week[j]; m++)
        num_infectious[m]++;

  for (m = 0; m < nweeks; m++)
    if (m == 20 || m == 21)
      pars.lambda[m] = num_infectious[m]*beta_housed;
    else
      pars.lambda[m] = num_infectious[m]*beta_field;

  // cumulative sum of lambda
  pars.Lambda[0] = 0;
  for (m = 1; m < nweeks; m++)
    pars.Lambda[m] = pars.Lambda[m-1]+pars.lambda[m-1];

  for (m = 0; m < nweeks; m++) {
    for (i = 0; i < m; i++)
      pars.W4[i][m] = w(i, i+1, m+1, &pars)
                    - w(i, i+1, m,   &pars)
                    - w(i, i,   m+1, &pars)
                    + w(i, i,   m,   &pars);
    pars.W2[m] = w(m, m, m+1, &pars);
  }

  // for each testing date calculate the probability that a sheep tests positive since the last test
  for (m = 0; m < nweeks; m++)
    p_seropos[m][V0] = Prob_seroconvert(m, &pars);

  /******************** OUTPUT ***********************/
  if (Data->mode == OUTPUT) {
    for (m = 0; m < nweeks; m++)
       V[m][0] = m;

    // expected number of new positives on each test date
    for (j = NUM_DONORS; j < NUMBER_SHEEP; j++)
      for (i = 0; i < sheep_last_test_week[j]-1; i++)
        V[i+1][2] += p_seropos[i][V0];

    // cumulative expected number of new positives
    V[0][3] = V[0][2];
    for (i = 0; i < nweeks-1; i++)
      V[i+1][3] = V[i][3]+V[i+1][2];

    // seroconversion rate (new positives per week)
    // V[0][4] = 0;
    // for (i = 1; i < Data->n; i++)
    //   V[i][4] = V[i][2]/(test_weeks[i]-test_weeks[i-1]);

    // // number of infectious animals
    // for (i = 0; i < Data->n; i++)
    //   V[i][5] = num_infectious[test_weeks[i]];

    // Einf_field = doublevector(nweeks);
    // Einf_housed = doublevector(nweeks);
    // for (j = NUM_DONORS; j < NUMBER_SHEEP; j++)
    //   for (m = 0; m <= SHEEP_LAST_WEEK; m++)
    //     if (test_weeks[m] == 19 || test_weeks[m] == 20)
    //       Einf_housed[m] += (1.-exp(-pars.lambda[m]))*u(m, pars.Lambda);
    //     else
    //       Einf_field[m] += (1.-exp(-pars.lambda[m]))*u(m, pars.Lambda);

    // for (m = 1; m < nweeks; m++) {
    //   Einf_field[m] += Einf_field[m-1];
    //   Einf_housed[m] += Einf_housed[m-1];
    // }
    // for (i = 0; i < Data->n; i++) {
    //   V[i][1] = Einf_field[test_weeks[i]];
    //   V[i][6] = Einf_housed[test_weeks[i]];
    //   V[i][7] = V[i][1]+V[i][6];
    // }
    // free(Einf_field);
    // free(Einf_housed);

    // number infected
    // num_infected = integervector(nweeks);
    // for (j = 0; j < NUMBER_SHEEP; j++)
    //   if (sheep_pos_week[j] != -1)
    //     for (m = max(0, sheep_pos_week[j]-seroconversion); m <= SHEEP_LAST_WEEK; m++)
    //       num_infected[m]++;

    // for (i = 0; i < Data->n; i++)
    //   V[i][8] = num_infected[test_weeks[i]];
    // free(num_infected);
  }

  free(num_infectious);
  free(pars.lambda);
  free(pars.Lambda);
  free(pars.W4[0]); free(pars.W4);
  free(pars.W2);
  return SUCCESS;
}

double P_pos(int s1, int s2, double **V)
{
  int m;
  double sum = 0;
  for (m = s1; m < s2; m++)
    sum += V[m][V0];
  return sum;
}

double logLikelihood(const uint trt, const uint ind, const double *p, double **V, const TimePoints *Data)
{
  int j, pos_idx;
  double lnL = 0;

  for (j = NUM_DONORS; j < NUMBER_SHEEP; j++)
    if (sheep_pos_index[j] != 0) {
      pos_idx = sheep_pos_index[j];
      if (pos_idx > 0) {
        do {
          pos_idx--;
        } while (Data->t[pos_idx][j] == 0);
        // ewe is first seropositive at start of month m, therefore it seroconverted
        // sometime in the months from the last test to the month prior to seroconverting
        lnL += log(P_pos(lrint(Data->t[pos_idx][0]), sheep_pos_week[j], V));
      }
      else
        // ewe never seroconverted
        lnL += log(1. - P_pos(SHEEP_FIRST_TEST_WEEK, sheep_last_test_week[j], V));
    }
  return lnL;
}

void OutputModel(int trt, int ind, double *Output, double *V)
{
  int i;
  for (i = 1; i < V0; i++)
    Output[i] = V[i];
}

void OutputData(int trt, int ind, double *Output, double *Var, double *Data, uint *value, int index)
{
  Output[2] = Data[31];
  Output[3] = Data[32];
  Output[4] = Data[33];
  Output[5] = Data[34];
}

int f0(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 2 month housed
  int i, n = 2;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = 100*i+1;
    y[i] = R0(x[i], 2., 0.1, 5, p);
  }
  return n;
}

int f1(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 0 months housed
  int i, n = 2;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = 100*i+1;
    y[i] = R0(x[i], 0., 0.1, 5, p);
  }
  return n;
}

int f2(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 for 25 sheep housed for different durations
  int i, n = 14;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = i;
    y[i] = R0(25., x[i]*12./365., 0.1, 5, p);
  }
  return n;
}

int f3(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // CDF of normal seroconversion period
  int i, n = 1000;
  if (x == NULL) return n;
  double S_a = p[2];
  double S_k1 = p[3];

  for (i = 0; i < n; i++) {
    x[i] = (double)i/20.;
    // y[i] = gsl_ran_gamma_pdf(x[i], a, 1./b);
    y[i] = gsl_cdf_gamma_P(x[i], S_a, S_k1);
  }
  return n;
}

void function_list(functions_t *F)
{
  F->n_func = 0;
  // F->function_list = (function_ptr*) malloc(F->n_func*sizeof(function_ptr));
  // F->function_list[0] = f0;
  // F->function_list[1] = f1;
  // F->function_list[2] = f2;
  // F->function_list[3] = f3;
}

double u(int m, double *Lambda)
{
  // The probability a sheep is uninfected at the start of week m
  return exp(-Lambda[m]);
}

double lngseries(double a, double x, int *ierr)
{
  double giant = HUGE_VAL/1000., eps = 1e-15;
  double t = 1./a, v = t;
  double p, q, lnigam;
  int k = 0;
  *ierr = 0;

  while ((fabs(t/v) > eps) && *ierr == 0) {
    p = (a+k)/(a+k+1);
    q = k+1;
    t *= -x*p/q;
    v += t;
    k += 1;
    if (t > giant)
      *ierr = 1;
  }
  if (*ierr == 0)
    if (lngamma(a) < log(giant))
      lnigam = log(v)-lngamma(a);
    else {
      lnigam = 0;
      *ierr = 1;
    }
  else {
    lnigam = 0;
    *ierr = 1;
  }
  return lnigam;
}

double lngexpan(double a, double x, int *ierr)
{
  double giant = HUGE_VAL/1000., eps = 1e-15;
  double t = 1, v = t;
  double p, lnigam;
  int k = 0;
  *ierr = 0;

  if (x > log(giant)) {
      *ierr = 1;
      return 0;
  }
  while ((fabs(t/v) > eps) && *ierr == 0) {
    p = 1-a+k;
    t *= p/x;
    v += t;
    k += 1;
    if (t > giant)
      *ierr = 1;
  }
  if (*ierr == 0)
    if (lngamma(a) < log(giant))
      lnigam = log(v)+x-log(x)-lngamma(a);
    else {
      lnigam = 0;
      *ierr = 1;
    }
  else {
    lnigam = 0;
    *ierr = 1;
  }
  return lnigam;
}

double ln_inc_gamma_star(double a, double z)
{
  int ierr;
  double lnigam;

  if (z < -50) {
    lnigam = lngexpan(a, -z, &ierr);
    if (ierr != 0)
      printf("error1\n");
  }
  else {
    lnigam = lngseries(a, z, &ierr);
    if (ierr != 0)
      printf("error2\n");
  }
  return lnigam;
}

double w(int n, int t, int s, parameters_t *p)
{
  double a = p->a;
  double l = p->lambda[n];
  double b = p->b;
  double r = (b-l)/b;
  double Z, z, v;

  Z = b*(s-t);
  if (Z == 0)
    return 0;

  z = r*Z;
  if (z <= 0)
    v = exp(a*log(Z) - l*Z/b + ln_inc_gamma_star(a, z));
  else
    v = exp(-a*log(r) - l*Z/b + log(P_inc_gamma(a, z)));

  return exp(-l*(t-n))*(v - P_inc_gamma(a, Z));
}

double Prob_seroconvert(int m, parameters_t *p)
{
  int n;
  double P_m = 0;

  for (n = 0; n < m; n++)
    P_m += u(n, p->Lambda)*p->W4[n][m];
  P_m -= u(m, p->Lambda)*p->W2[m];

  return min(max(0, P_m), 1);
}

void DerivedParameters(int trt, int ind, int nmonths, double *p, double **V, TimePoints *Data)
{
  // int i, j, m;
  // double a, b;
  // double beta_field = p[0];
  // double beta_housed = p[1];
  // double mu = p[2];
  // double sd = p[3];
  // int seroconversion = lrint(mu);
  // int latency = lrint(p[4]);
  // int nweeks = lrint(Data->t[Data->n-1][0])+1;
  // int *test_weeks = integervector(Data->n);
  // int *num_infectious = integervector(nweeks);
  // double *Lambda = doublevector(nweeks);
  // double *lambda = doublevector(nweeks);

  // p[8]  = a = pow(mu/sd, 2.);
  // p[9]  = b = a/mu;
  // p[10] = beta_housed/beta_field;
  // p[15] = P_inv_gamma(0.025, a, 1./b);
  // p[16] = P_inv_gamma(0.975, a, 1./b);

  // for (i = 0; i < Data->n; i++)
  //   test_weeks[i] = lrint(Data->t[i][0]);

  // for (j = 0; j < NUM_DONORS; j++)
  //   for (m = 0; m <= SHEEP_LAST_WEEK; m++)
  //       num_infectious[m]++;

  // for (j = NUM_DONORS; j < NUMBER_SHEEP; j++)
  //   if (sheep_pos_week[j] > 0)
  //     for (m = max(0, sheep_pos_week[j]-seroconversion+latency); m <= SHEEP_LAST_WEEK; m++)
  //       num_infectious[m]++;

  // for (m = 0; m < nweeks; m++)
  //   if (test_weeks[m] == 19 || test_weeks[m] == 20)
  //     lambda[m] = num_infectious[m]*beta_housed;
  //   else
  //     lambda[m] = num_infectious[m]*beta_field;

  // Lambda[0] = 0;
  // for (m = 1; m < nweeks; m++)
  //   Lambda[m] = Lambda[m-1]+lambda[m-1];

  // p[12] = 0;
  // p[13] = 0;
  // for (j = NUM_DONORS; j < NUMBER_SHEEP; j++)
  //   for (m = 0; m <= SHEEP_LAST_WEEK; m++)
  //     if (test_weeks[m] == 19 || test_weeks[m] == 20)
  //       p[13] += (1.-exp(-lambda[m]))*u(m, Lambda);
  //     else
  //       p[12] += (1.-exp(-lambda[m]))*u(m, Lambda);
  // p[14] = p[12]+p[13];

  // free(lambda);
  // free(Lambda);
  // free(num_infectious);
  // free(test_weeks);
}

void WAIC(int trt, int ind, double **lnL, double **V, double *p, const TimePoints *Data) {}
void PredictData(int trt, int ind, double *Output, double *V, double *p, gsl_rng *stream) {}
void SimulateData(int trt, int ind, double *Output, double *V, double *p, double *Data, uint *value, gsl_rng *stream) {}
double timestep(void) {return 1;}
void GlobalParameters() {}
double UserPDF(double x) {return 0;}
void HyperParameters(Treatments T, Hyperparameters *H) {}
void SaturatedModel(int trt, int ind, double **V, double *p, const TimePoints *Data) {}
void Residual(int trt, int ind, double *R, double *V, double *p, double *Data, uint *value) {}

max_stars_repo_name: nicksavill/maedi-visna-epidemiology

max_issues_repo_name: nicksavill/maedi-visna-epidemiology

max_forks_repo_name: nicksavill/maedi-visna-epidemiology


text:
/* eigen/gsl_eigen.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_EIGEN_H__
#define __GSL_EIGEN_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct {
  size_t size;
  double * d;
  double * sd;
} gsl_eigen_symm_workspace;

gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * gc;
  double * gs;
} gsl_eigen_symmv_workspace;

gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
} gsl_eigen_herm_workspace;

gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
                         gsl_eigen_herm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
  double * gc;
  double * gs;
} gsl_eigen_hermv_workspace;

gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
                           gsl_matrix_complex * evec,
                           gsl_eigen_hermv_workspace * w);

typedef enum {
  GSL_EIGEN_SORT_VAL_ASC,
  GSL_EIGEN_SORT_VAL_DESC,
  GSL_EIGEN_SORT_ABS_ASC,
  GSL_EIGEN_SORT_ABS_DESC
}
gsl_eigen_sort_t;

/* Sort eigensystem results based on eigenvalues.
 * Sorts in order of increasing value or increasing
 * absolute value.
 *
 * exceptions: GSL_EBADLEN
 */

int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
                         gsl_eigen_sort_t sort_type);

int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

/* The following functions are obsolete: */

/* Eigensolve by Jacobi Method
 *
 * The data in the matrix input is destroyed.
 *
 * exceptions:
 */
int
gsl_eigen_jacobi(gsl_matrix * matrix,
                      gsl_vector * eval,
                      gsl_matrix * evec,
                      unsigned int max_rot,
                      unsigned int * nrot);

/* Invert by Jacobi Method
 *
 * exceptions:
 */
int
gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
                             gsl_matrix * ainv,
                             unsigned int max_rot);

__END_DECLS

#endif /* __GSL_EIGEN_H__ */

ext: h


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

#include <stdarg.h>

#include "compat.h"
#include "gslutils.h"

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

int gslutils_solve_leastsquares_v(gsl_matrix* A, int NB, ...) {
    int i, res;
    gsl_vector**  B = malloc(NB * sizeof(gsl_vector*));
    // Whoa, three-star programming!
    gsl_vector*** X = malloc(NB * sizeof(gsl_vector**));
    gsl_vector*** R = malloc(NB * sizeof(gsl_vector**));

    gsl_vector** Xtmp = malloc(NB * sizeof(gsl_vector*));
    gsl_vector** Rtmp = malloc(NB * sizeof(gsl_vector*));

    va_list va;
    va_start(va, NB);
    for (i=0; i<NB; i++) {
        B[i] = va_arg(va, gsl_vector*);
        X[i] = va_arg(va, gsl_vector**);
        R[i] = va_arg(va, gsl_vector**);
    }
    va_end(va);

    res = gslutils_solve_leastsquares(A, B, Xtmp, Rtmp, NB);
    for (i=0; i<NB; i++) {
        if (X[i])
            *(X[i]) = Xtmp[i];
        else
            gsl_vector_free(Xtmp[i]);
        if (R[i])
            *(R[i]) = Rtmp[i];
        else
            gsl_vector_free(Rtmp[i]);
    }
    free(Xtmp);
    free(Rtmp);
    free(X);
    free(R);
    free(B);
    return res;
}

int gslutils_solve_leastsquares(gsl_matrix* A, gsl_vector** B,
                                gsl_vector** X, gsl_vector** resids,
                                int NB) {
    int i;
    gsl_vector *tau, *resid = NULL;
    int ret;
    size_t M, N;

    M = A->size1;
    N = A->size2;

    for (i=0; i<NB; i++) {
        assert(B[i]);
        assert(B[i]->size == M);
    }

    tau = gsl_vector_alloc(MIN(M, N));
    assert(tau);

    ret = gsl_linalg_QR_decomp(A, tau);
    assert(ret == 0);
    // A,tau now contains a packed version of Q,R.

    for (i=0; i<NB; i++) {
        if (!resid) {
            resid = gsl_vector_alloc(M);
            assert(resid);
        }
        X[i] = gsl_vector_alloc(N);
        assert(X[i]);
        ret = gsl_linalg_QR_lssolve(A, tau, B[i], X[i], resid);
        assert(ret == 0);
        if (resids) {
            resids[i] = resid;
            resid = NULL;
        }
    }

    gsl_vector_free(tau);
    if (resid)
        gsl_vector_free(resid);

    return 0;
}

ext: c

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: mgreter/astrometrylib

max_forks_repo_name: mgreter/astrometrylib


text:
/// @file Span.h
/// @author Braxton Salyer <braxtonsalyer@gmail.com>
/// @brief A thin wrapper around `gsl::span`, providing `at` as a member function instead of a free
/// function
/// @version 0.1
/// @date 2021-10-15
///
/// MIT License
/// @copyright Copyright (c) 2021 Braxton Salyer <braxtonsalyer@gmail.com>
///
/// Permission is hereby granted, free of charge, to any person obtaining a copy
/// of this software and associated documentation files (the "Software"), to deal
/// in the Software without restriction, including without limitation the rights
/// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
/// copies of the Software, and to permit persons to whom the Software is
/// furnished to do so, subject to the following conditions:
///
/// The above copyright notice and this permission notice shall be included in all
/// copies or substantial portions of the Software.
///
/// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
/// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
/// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
/// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
/// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
/// SOFTWARE.
#pragma once

#include <Hyperion/BasicTypes.h>
#include <Hyperion/HyperionDef.h>
#include <array>
#include <gsl/gsl>
#include <gsl/span>

namespace hyperion {

	/// @brief Thin wrapper around `gsl::span`
	///
	/// @tparam T - The type contained in the `Span`
	/// @tparam Size - The number of elements for constexpr sizes, or gsl::dynamic_extent (default)
	/// for run-time sizes
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename T, usize Size = gsl::dynamic_extent>
	class Span {
	  public:
		using Iterator = typename gsl::span<T, Size>::iterator;
		using ReverseIterator = typename gsl::span<T, Size>::reverse_iterator;
		template<usize Offset, usize Count>
		using SubSpan = typename gsl::details::calculate_subspan_type<T, Size, Offset, Count>::type;

		Span() noexcept = default;

		/// @brief Constructs a `Span` from a `gsl::span`
		///
		/// @param span - The `gsl::span` to wrap
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		explicit constexpr Span(gsl::span<T, Size> span) noexcept : m_span_internal(span) {
		}

		/// @brief Copy constructs a `Span` from the given one
		///
		/// @param span - The `Span` to copy
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		constexpr Span(const Span<T, Size>& span) noexcept = default;

		/// @brief Move constructs the given `Span`
		///
		/// @param span - The `Span` to move
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		constexpr Span(Span<T, Size>&& span) noexcept = default;
		~Span() noexcept = default;

		/// @brief Returns the element at the given `index`
		///
		/// @param index - The index of the desired element
		/// @return - The element at `index`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto at(usize index) noexcept -> T& {
			return gsl::at(m_span_internal, static_cast<gsl::index>(index));
		}

		/// @brief Returns the first `Count` elements in the `Span`
		///
		/// @tparam Count - The number of elements to get
		/// @tparam size - The size of the `Span`, this should always be left defaulted
		/// @return - The first `Count` elements in the `Span`, as a `Span`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		template<usize Count, usize size = Size>
		requires(size != gsl::dynamic_extent)
			[[nodiscard]] inline constexpr auto first() const noexcept -> Span<T, Count> {
			return Span(m_span_internal.template first<Count>());
		}

		/// @brief Returns the first `count` elements in the `Span`
		///
		/// @tparam size - The size of the `Span`, this should always be left defaulted
		/// @param count - The number of elements to get
		/// @return - The first `count` elements in the `Span`, as a `Span`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		template<usize size = Size>
		requires(size == gsl::dynamic_extent)
			[[nodiscard]] inline constexpr auto first(usize count) const noexcept -> Span<T, Size> {
			return Span(m_span_internal.first(count));
		}

		/// @brief Returns the last `Count` elements in the `Span`
		///
		/// @tparam Count - The number of elements to get
		/// @tparam size - The size of the `Span`, this should always be left defaulted
		/// @return - The last `Count` elements in the `Span`, as a `Span`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		template<usize Count, usize size = Size>
		requires(size != gsl::dynamic_extent)
			[[nodiscard]] inline constexpr auto last() const noexcept -> Span<T, Count> {
			return Span(m_span_internal.template last<Count>());
		}

		/// @brief Returns the last `count` elements in the `Span`
		///
		/// @tparam size - The size of the `Span`, this should always be left defaulted
		/// @param count - The number of elements to get
		/// @return- The last `count` elements in the `Span` as a `Span`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		template<usize size = Size>
		requires(size == gsl::dynamic_extent)
			[[nodiscard]] inline constexpr auto last(usize count) const noexcept -> Span<T, Size> {
			return Span(m_span_internal.last(count));
		}

		/// @brief Returns a subspan starting at the given `Offset` with `Count` elements
		///
		/// @tparam Offset - The offset to start the subspan at
		/// @tparam Count - The number of elements to get in the subspan
		/// @return - The subspan starting at `Offset` of size `Count`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		template<usize Offset, usize Count = gsl::dynamic_extent>
		[[nodiscard]] inline constexpr auto subspan() const noexcept -> SubSpan<Offset, Count> {
			return Span(m_span_internal.template subspan<Offset, Count>());
		}

		/// @brief Returns a subspan starting at the given `offset` with `count` elements
		///
		/// @param offset - The offset to start the subspan at
		/// @param count - The number of elements to get in the subspan
		/// @return - The subspan starting at `offset` of size `count`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto
		subspan(usize offset, usize count = gsl::dynamic_extent) const noexcept
			-> Span<T, gsl::dynamic_extent> {
			return Span(m_span_internal.subspan(offset, count));
		}

		/// @brief Returns the number of elements in the `Span`
		///
		/// @return - The number of elements in the `Span`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto size() const noexcept -> usize {
			return m_span_internal.size();
		}

		/// @brief Returns the size of the `Span` in bytes
		///
		/// @return - The size of the `Span` in bytes
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto size_bytes() const noexcept -> usize {
			return m_span_internal.size_bytes();
		}

		/// @brief Returns whether the `Span` is empty
		///
		/// @return - If the `Span` is empty
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto is_empty() const noexcept -> bool {
			return size() == 0;
		}

		/// @brief A pointer to the data contained in the `Span`
		///
		/// @return - A pointer to the data
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto data() const noexcept -> T* {
			return m_span_internal.data();
		}

		/// @brief Returns the first element in the `Span`
		///
		/// @return - The first element
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto front() const noexcept -> T& {
			return m_span_internal.front();
		}

		/// @brief Returns the last element in the `Span`
		///
		/// @return - The last element
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto back() const noexcept -> T& {
			return m_span_internal.back();
		}

		/// @brief Returns an iterator at the beginning of the span
		///
		/// @return an iterator at the beginning of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto begin() noexcept -> Iterator {
			return m_span_internal.begin();
		}

		/// @brief Returns an iterator at the beginning of the span
		///
		/// @return an iterator at the beginning of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto begin() const noexcept -> Iterator {
			return m_span_internal.begin();
		}

		/// @brief Returns an iterator at the end of the span
		///
		/// @return an iterator at the end of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto end() noexcept -> Iterator {
			return m_span_internal.end();
		}

		/// @brief Returns an iterator at the end of the span
		///
		/// @return an iterator at the end of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto end() const noexcept -> Iterator {
			return m_span_internal.end();
		}

		/// @brief Returns an iterator at the beginning of the reversed iteration of the span
		///
		/// @return an iterator at the beginning of the reversed iteration of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto rbegin() noexcept -> ReverseIterator {
			return m_span_internal.rbegin();
		}

		/// @brief Returns an iterator at the beginning of the reversed iteration of the span
		///
		/// @return an iterator at the beginning of the reversed iteration of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto rbegin() const noexcept -> ReverseIterator {
			return m_span_internal.rbegin();
		}

		/// @brief Returns an iterator at the end of the reversed iteration of the span
		///
		/// @return an iterator at the end of the reversed iteration of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto rend() noexcept -> ReverseIterator {
			return m_span_internal.rend();
		}

		/// @brief Returns an iterator at the end of the reversed iteration of the span
		///
		/// @return an iterator at the end of the reversed iteration of the span
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto rend() const noexcept -> ReverseIterator {
			return m_span_internal.rend();
		}

#if HYPERION_PLATFORM_COMPILER_MSVC
		[[nodiscard]] inline constexpr auto _Unchecked_begin() const noexcept -> T* {
			return m_span_internal._Unchecked_begin();
		}

		[[nodiscard]] inline constexpr auto _Unchecked_end() const noexcept -> T* {
			return m_span_internal._Unchecked_end();
		}
#endif // HYPERION_PLATFORM_COMPILER_MSVC

		/// @brief Returns the element at the given `index`
		///
		/// @param index - The index of the desired element
		/// @return - The element at `index`
		/// @ingroup utils
		/// @headerfile "Hyperion/Span.h"
		[[nodiscard]] inline constexpr auto operator[](usize index) noexcept -> T& {
			return this->at(index);
		}

		constexpr auto operator=(const Span<T, Size>& span) noexcept -> Span<T, Size>& = default;
		constexpr auto operator=(Span<T, Size>&& span) noexcept -> Span<T, Size>& = default;

	  private:
		gsl::span<T, Size> m_span_internal = gsl::span<T, Size>();
	};

	IGNORE_UNUSED_TEMPLATES_START

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam T - The type stored in the array
	///
	/// @param array - The array to get a `Span` over
	/// @param size - The size of the array
	///
	/// @return a `Span` over the given array
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename T>
	[[nodiscard]] inline static constexpr auto make_span(T* array, usize size) noexcept -> Span<T> {
		return Span(gsl::make_span(array, static_cast<typename gsl::span<T>::size_type>(size)));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam T - The type stored in the array
	///
	/// @param first - Pointer to the first element in the array
	/// @param last - Pointer to the last element in the array
	///
	/// @return a `Span` over the given array
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename T>
	[[nodiscard]] inline static constexpr auto make_span(T* first, T* last) noexcept -> Span<T> {
		return Span(gsl::make_span(first, last));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam T - The type stored in the array
	/// @tparam Size - The size of the array
	///
	/// @param array - The array to get a `Span` over
	///
	/// @return a `Span` over the given array
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename T, usize Size>
	[[nodiscard]] inline static constexpr auto
	make_span(T (&array)[Size]) noexcept -> Span<T, Size> { // NOLINT
		return Span(gsl::make_span(array));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam Collection - The collection type to get a `Span` over
	///
	/// @param collection - The collection to get a `Span` over
	///
	/// @return a `Span` over the given collection
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename Collection>
	[[nodiscard]] inline static constexpr auto
	make_span(Collection& collection) noexcept -> Span<typename Collection::value_type> {
		return Span(gsl::make_span(collection));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam Collection - The collection type to get a `Span` over
	///
	/// @param collection - The collection to get a `Span` over
	///
	/// @return a `Span` over the given collection
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename Collection>
	[[nodiscard]] inline static constexpr auto make_span(const Collection& collection) noexcept
		-> Span<const typename Collection::value_type> {
		return Span(gsl::make_span(collection));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam SmartPtr - The smart pointer type wrapping the array
	///
	/// @param array - The smart pointer wrapping the array to get a `Span` over
	/// @param size - The size of the array
	///
	/// @return a `Span` over the given array
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename SmartPtr>
	[[nodiscard]] inline static constexpr auto
	make_span(SmartPtr& array, usize size) noexcept -> Span<typename SmartPtr::element_type> {
		return Span(gsl::make_span(array, size));
	}

	/// @brief Creates a `Span` from the given arguments
	///
	/// @tparam T - The type stored in the array
	/// @tparam Size - The size of the array
	///
	/// @param array - The array to get a `Span` over
	///
	/// @return a `Span` over the given array
	/// @ingroup utils
	/// @headerfile "Hyperion/Span.h"
	template<typename T, usize Size>
	[[nodiscard]] inline static constexpr auto
	make_span(std::array<T, Size> array) noexcept -> Span<T, Size> {
		return Span(gsl::make_span(array));
	}

	IGNORE_UNUSED_TEMPLATES_STOP
} // namespace hyperion