text:
/*********************************************************************
 * Program: psimci.c
 * Author:  Mauricio Caceres Bravo <caceres@nber.org>
 * Created: Sun Feb 12 19:28:43 EST 2017
 * Updated: Tue May 30 18:11:50 EDT 2017
 * Purpose: Stata plugin to simulate a CI under H0: b = 0 for a
 *          treatment effect given a regression specification.
 * Note:    See stata.com/plugins for more on Stata plugins
 *********************************************************************/

/**
 * @file psimci.c
 * @author Mauricio Caceres bravo
 * @date 30 May 2017
 * @brief Stata plugin to simulate a CI for a placebo treatment.
 *
 * See the documentation for simci.ado (e.g. help simci from Stata)
 *
 * @see http://www.stata.com/plugins
 */

#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_vector.h>
#include "psimci.h"
#include "stplugin.h"
#include "stutils.c"

/**
 * @brief Main function call will execute as Stata plugin
 *
 * The function takes the first variable in @argc as the dependent
 * variable and the next k - 1 variables are covariates. @argc contains
 * the comma options passed by Stata. Currently just the proportion
 * randomzed and the number of simulations to run. See documentation for
 * simci.ado or sim_ci below for moew.
 *
 * @param argc List of variables to use
 * @param argv Comma options from Stata
 * @return Modified variables in Stata
 * @see Documentation for simci.ado
 * @warning This is meant to be run from simci.ado and not by itself
 */
STDLL stata_call(int argc, char *argv[])
{
    // Simple call to see if the plugin loaded
    char buffer[16]; strcpy (buffer, argv[0]);
    if ( strcmp(buffer, "check") == 0 ) return (0);

    // Initialize the variables to use
    ST_int      i, j ;
    ST_double   z ;
    ST_retcode  rc ;

    // Get P and number of reps. Note the 0-based indexing! So the
    // functiona ssumes P and reps were the 1st and 3nd argument.
    double P    = strtod (argv[0], NULL);
    int    reps = strtod (argv[1], NULL);

    const size_t n = SF_in2();
    const int    k = SF_nvars();

    // If too few variables (at least 2 for regressio), exit
    if (k < 2) {
        return (102) ;
    }

    // Initialize GSL elements where to store data
    gsl_matrix *X  = gsl_matrix_alloc (n, k + 1);
    gsl_vector *y  = gsl_vector_alloc (n);

    // Not sure if there is another way to read data vs the double loop.
    // Note: Careful with the 0-based indexing!
    for (i = SF_in1(); i <= SF_in2(); i++) {
        if (SF_ifobs(i)) {

            // Variables 2 through k are covariates
            for (j = 2; j <= k; j++) {
                // Note we leave the first column empty
                if ( (rc = SF_vdata(j, i, &z)) ) return(rc);
                gsl_matrix_set (X, i - 1, j - 1, z);
            }

            // Note we add the constant
            gsl_matrix_set (X,  i - 1, k, 1.0);

            // Variable 1 is the dependent variable
            if ( (rc = SF_vdata(1, i, &z)) ) return(rc);
            gsl_vector_set (y,  i - 1, z);
        }
    }

    // Now we call the simulation function and output the results into b, mu
    gsl_vector *b   = gsl_vector_alloc (reps);
    gsl_vector *mu  = gsl_vector_alloc (reps);
    sim_ci (X, y, P, reps, b, mu);

    // Not sure that there is a good way to output this into Stata So I
    // write to a file and read it back.
    char outb[64], outmu[64];
    strcpy (outb,  argv[2]);
    strcpy (outmu, argv[2]);
    strcat (outb,  "b");
    strcat (outmu, "mu");
    FILE *fb  = fopen (outb, "wb");
    FILE *fmu = fopen (outmu, "wb");
    gsl_vector_fprintf (fb,  b,  "%15.9f");
    gsl_vector_fprintf (fmu, mu, "%15.9f");
    fclose (fb);
    fclose (fmu);

    // Cleanup
    gsl_matrix_free (X);
    gsl_vector_free (y);
    gsl_vector_free (b);
    gsl_vector_free (mu);

    return (0);
}

/**
 * @brief Simulate a confidence interval given X, y
 *
 * The idea is to simulate a non-parametric CI based on placebo
 * assignments of a treatment variable. The program assigns
 * treatment at random, hence a null effect, to individuals or
 * clusters, optionally stratifying by any number of variables (or
 * the means thereof, in the case of clusters). Consider
 *
 *     Y_ij = a + b T_j + g X_ij + e_ij
 *
 * There are C = J choose PJ ways to treat the clusters (or C =
 * P choose PN in the case of individuals). If we computed b_ols
 * for c = 1, ..., C we would know the exact distribution of our
 * estimator, conditional on the data being representative of the
 * study data. C is typically intractably large, hence we simulate
 * K draws with sum(T_jk = PJ) and run
 *
 *     Y_ij = a + b_k T_jk + g X_ij + e_ij
 *
 * Let Fhat be the empirical cdf of b_k; a valid 1 - alpha CI for
 * b is given by
 *
 *     CI(1 - a) = [Fhat^-1(a / 2), Fhat^-1(1 - a / 2)]
 *
 * The function takes the @X as the covariate matrix, which
 * must have k + 1 columns with the first column free, @y as
 * the dependent variable, and outputs the results to @b, @mu
 *
 * @param X Covariate matrix with first column blank
 * @param y Dependent variable
 * @param P Proportion in treatment
 * @param reps Number of reprtitions
 * @param b Vector of length @reps; will output coefficients here
 * @param mu Vector of length @reps; will output control means here
 * @return Modified @b, @mu with coefficients and means
 * @see Documentation for simci.ado
 */
int sim_ci (const gsl_matrix * X,
            const gsl_vector * y,
            const double P,
            const int reps,
            gsl_vector * b,
            gsl_vector * mu)
{

    const size_t n  = X->size1;
    const int k     = X->size2;
    const int np    = ceil(n * P);
    const int nc    = n - np;
    double *sy      = malloc (sizeof(double));

    gsl_vector *ones = gsl_vector_alloc (n);
    gsl_vector_set_all (ones, 1.0);
    gsl_blas_ddot (ones, y, sy);

    // Set the random seed based on the time of day (seconds)
    srand (time(NULL));
    gsl_rng *rng = gsl_rng_alloc (gsl_rng_default);
    gsl_rng_set (rng, rand());

    // Get vector of 1s and 0s
    gsl_vector *T = gsl_vector_alloc (n);
    gsl_vector_set_zero (T);
    for (int i = 0; i < np; i++) {
        gsl_vector_set (T, i, 1.0);
    }

    // Initialize elements for parallel loop
    gsl_vector *Tp ;
    gsl_matrix *Xp ;
    int nloops ;
    double *sty ;

    // Get the number of threads available to OMP
    sf_printf("Parallelizing simulation; %d threads found:\n",
              get_omp_num_threads());

    // Parallelize execution: Note We need a copy of Xp and Tp for each
    // thread since they will be modified at each iteration y does not
    // change, so it's shared.
    #pragma omp parallel private(Xp, Tp, nloops, sty) shared(y, b, sy)
    {
        nloops = 0;

        // Allocate to each therad their own copy
        Tp  = gsl_vector_alloc (n);
        Xp  = gsl_matrix_alloc (n, k);
        sty = malloc (sizeof(double));

        gsl_vector_memcpy (Tp, T);
        gsl_matrix_memcpy (Xp, X);

        // Parallel for loop through simulation
        #pragma omp for
        for (int r = 0; r < reps; r++) {
            // 1. Shuffle treatment
            // 2. Set as first column of covariate matrix
            // 3. Get mean of y over controls
            // 4. Store coefficient/mean
            // 5. Repeat 1-4
            // 6. ...
            // 7. Profit?
            gsl_ran_shuffle (rng, Tp->data, n, sizeof(size_t));
            gsl_matrix_set_col (Xp, 0, Tp);
            gsl_vector_set (b, r, sim_ols(Xp, y));
            gsl_blas_ddot (Tp, y, sty);
            gsl_vector_set (mu, r, (*sy - *sty) / nc);
            ++nloops;
        }

        // I want to print a pretty message saying how many iterations
        // each thread completed. Since threads finish on their own,
        // messages would be print at disparate times. However, one can
        // specify "critical" code which is executed only after all
        // threads are done running.
        #pragma omp critical
        {
            sf_printf("\tThread %d performed %d simulations.\n",
                      omp_get_thread_num(), nloops);
        }

        // Cleanup
        gsl_matrix_free (Xp);
        gsl_vector_free (Tp);
    }

    // Cleanup
    gsl_vector_free (T);
    gsl_rng_free (rng);

    return (0);
}

/**
 * @brief Number of threads available to OMP
 *
 * Short wrapper to get number of threads available to OMP
 *
 * @return Number of threads available to OMP
 */
int get_omp_num_threads()
{
    int thread_id;
    int nthreads = 0;

    #pragma omp parallel private(thread_id) shared(nthreads)
    {
        thread_id = omp_get_thread_num();
        #pragma omp critical
        {
            nthreads = thread_id > nthreads? thread_id: nthreads;
        }
    }
    nthreads++;

    return (nthreads);
}

/**
 * @brief Wrapper to run a linear regression
 *
 * All I want is the first coefficient of a linear regression. For
 *
 *     Y = X beta
 *
 * I want (X' X)^-1 X' Y. GSL has solvers for a system of the form
 *
 *     Ax = b
 *
 * Where A is a symmetric matrix. Take A = X' X and b = X' y, then
 * we can use any number of routines to find x (especially since A
 * is now symmetric).
 *
 * @param X A n by k gsl matrix containing covariates.
 * @param y A n by 1 gsl vector containing the dependent variable
 * @return The first coefficient of a linear regression.
 * @warning This is meant to be run within the main loop of stata_call
 */
double sim_ols(const gsl_matrix * X, const gsl_vector * y)
{

    // Allocate memory to express the system as Ax = b
    gsl_matrix *A = gsl_matrix_alloc (X->size2, X->size2);
    gsl_vector *b = gsl_vector_alloc (X->size2);
    gsl_vector *x = gsl_vector_alloc (X->size2);

    // Set A = X' X and b = X' y
    gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X, X, 0.0, A);
    gsl_blas_dgemv (CblasTrans, 1.0, X, y, 0.0, b);

    // Cholesky decomposition
    gsl_linalg_cholesky_decomp1 (A);
    gsl_linalg_cholesky_solve (A, b, x);

    // You don't have to use Cholesky; a number of methods are available
    //
    // int s;
    // gsl_permutation * P = gsl_permutation_alloc (X->size2);
    // gsl_vector * tau    = gsl_vector_alloc (X->size2);
    //
    // Householder
    // gsl_linalg_HH_solve (A, b, x);
    //
    // LU decomposition
    // gsl_linalg_LU_decomp (A, P, &s);
    // gsl_linalg_LU_solve (A, P, b, x);
    // gsl_permutation_free (P);
    //
    // QR decomposition
    // gsl_linalg_QR_decomp (A, tau);
    // gsl_linalg_QR_solve (A, tau, b, x);
    // gsl_vector_free (tau);

    // Free up space
    gsl_matrix_free (A);
    gsl_vector_free (b);

    return (gsl_vector_get(x, 0));
}

/**
 * @brief Get pctile of a function
 *
 * Basic wrapper to get the @pctile percentile of a function.
 *
 * @param x n by 1 gsl vector whose percentile we want.
 * @param pctile Percentile
 * @return @pctile percentile of x
 */
double sim_pctile(gsl_vector * x, double pctile)
{
    gsl_sort_vector (x);
    int n = x->size;
    int i = floor(n * pctile);
    double qq = gsl_vector_get (x, i);
    if (i / n == pctile) {
        qq = (qq + gsl_vector_get (x, i + 1)) / 2;
    }
    return (qq);
}

max_stars_repo_name: mcaceresb/stata-power

max_issues_repo_name: arlionn/stata-power

max_forks_repo_name: arlionn/stata-power


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

#define CHKERR(t1, t2, tol, name, str)\
{\
   if (fabs(t1-t2) > tol){\
       printf("%s: %s %e %e\n", name, str, t1, t2);\
       return EXIT_FAILURE;\
   }\
}

int check_lam2nualpha(void);
int check_Udetcheck(void);
int check_Uorth(void);
int check_lam2lune(void);
int check_CMTdecom(void);
int check_fangleSigned(void);
int check_CMT2faultpar(void);
int check_CMT2TT(void);
int check_TT2CMT(void);
int check_normal2strdip(void);
int check_CMT2omega(void);

int main(void)
{
    int ierr;

    ierr = check_lam2nualpha();
    if (ierr != 0){printf("failed lam2nualpha\n"); return EXIT_FAILURE;}
    printf("lam2nualpha was successful\n");

    ierr = check_fangleSigned();
    if (ierr != 0){printf("failed fangleSigned\n"); return EXIT_FAILURE;}
    printf("fangleSigned was successful\n");

    ierr = check_Udetcheck();
    if (ierr != 0){printf("failed Udetcheck\n"); return EXIT_FAILURE;}
    printf("udetcheck was successful\n");

    ierr = check_Uorth();
    if (ierr != 0){printf("failed Uorth\n"); return EXIT_FAILURE;}
    printf("Uorth was successful\n");

    ierr = check_lam2lune();
    if (ierr != 0){printf("failed lam2lune\n"); return EXIT_FAILURE;}
    printf("lam2lune was succesful\n");

    ierr = check_CMTdecom();
    if (ierr != 0){printf("failed CMTdecom\n"); return EXIT_FAILURE;}
    printf("CMTdecom was successful\n");

    ierr = check_CMT2TT();
    if (ierr != 0){printf("failed CMT2TT\n"); return EXIT_FAILURE;}
    printf("CMT2TT was successful\n");

    ierr = check_TT2CMT();
    if (ierr != 0){printf("failed TT2CMT\n"); return EXIT_FAILURE;}
    printf("TT2CMT was successful\n");

    ierr = check_CMT2faultpar();
    if (ierr != 0){printf("failed CMT2faultpar\n"); return EXIT_FAILURE;}
    printf("CMT2faultpar was successful\n");

    ierr = check_normal2strdip();
    if (ierr != 0){printf("failed noraml2strdip\n"); return EXIT_FAILURE;}
    printf("normal2strdip was successful\n");

    ierr = check_CMT2omega();
    if (ierr != 0){printf("failed CMT2Omega\n"); return EXIT_FAILURE;}
    printf("CMT2Omega was successful\n");

    return EXIT_SUCCESS;
}
//============================================================================//
int check_fangleSigned(void)
{
    const char *fcnm = "check_fangleSigned\0";
    const double va[3] = {1, 0, 0};
    const double vb[3] = {0, 0, 3};
    const double vnor[3] = {0, -2, 0};
    int ierr;
    double stheta;
    stheta = compearth_eulerUtil_fangleSigned(3, va,vb,vnor, &ierr);
    if (fabs(stheta - 90.0) > 1.e-14 || ierr != 0)
    {
        printf("%s: failed test 1 %f %d\n", fcnm, stheta, ierr);
        return EXIT_FAILURE;
    }
    const double vnor2[3] =  {0, 2, 0};
    stheta = compearth_eulerUtil_fangleSigned(3, va,vb,vnor2, &ierr);
    if (fabs(stheta + 90.0) > 1.e-14 || ierr != 0)
    {
        printf("%s: failed test 2 %f %d\n", fcnm, stheta, ierr);
        return EXIT_FAILURE;
    }
    return EXIT_SUCCESS;
}
//============================================================================//
int check_Udetcheck(void)
{
    const char *fcnm = "check_Udetcheck\0";
    const double U[9] = {0.6948,   0.3171,   0.9502,  // col 1
                         0.0344,   0.4387,   0.3816,  // col 2
                         0.7655,   0.7952,   0.1869}; // col 3
    const double Ur[9] = {0.6948,  0.3171,   0.9502,  // col 1
                         -0.0344, -0.4387,  -0.3816,  // col 2 negated
                          0.7655,  0.7952,   0.1869};
    double Un[9], Un2[9];
    int i, ierr;
    ierr = compearth_Udetcheck(1, U, Un);
    if (ierr != 0)
    {
        printf("%s: Error in Udetcheck1\n", fcnm);
        return EXIT_FAILURE;
    }
    ierr = compearth_Udetcheck(1, Un, Un2);
    if (ierr != 0)
    {
        printf("%s: Error in Udetcheck2\n", fcnm);
        return EXIT_FAILURE;
    }
    for (i=0; i<9; i++)
    {
        if (fabs(Ur[i] - Un[i]) > 1.e-12)
        {
            printf("%s: Udetcheck failed; %f %f\n", fcnm, Ur[i], Un[i]);
            return EXIT_FAILURE;
        }
        if (fabs(Un[i] - Un2[i]) > 1.e-12)
        {
            printf("%s: Udetcheck2 failed; %f %f\n", fcnm, Un[i], Un2[i]);
            return EXIT_FAILURE;
        }
    }
    return EXIT_SUCCESS;
}
//============================================================================//
int check_lam2nualpha(void)
{
    const char *fcnm = "check_lam2nualpha\0";
    // example from TT2013, App A.
    double lam[3] = { 8.802, 2.584, -1.851};
    double lamcheck[3], nu, alpha, mag;
    const double nu1 = 0.371745072651417;
    const double alpha1 = 80.365019327257144;
    int ierr;
    ierr = compearth_lam2nualpha(1, lam, &nu, &alpha);
    if (ierr != 0)
    {
        printf("%s: error calling lam2nualpha\n", fcnm);
        return EXIT_FAILURE;
    }
    CHKERR(nu, nu1, 1.e-10, fcnm, "error computing nu 1");
    CHKERR(alpha, alpha1, 1.e-10, fcnm, "error computing alpha 1");
    ierr = compearth_nualpha2lam(1, &nu, &alpha, lamcheck);
    if (ierr != 0)
    {
        printf("%s: error caling nualpha2lam\n", fcnm);
        return EXIT_FAILURE;
    }
    mag = cblas_dnrm2(3, lam, 1);
    lam[0] = lam[0]/mag;
    lam[1] = lam[1]/mag;
    lam[2] = lam[2]/mag;
    CHKERR(lam[0], lamcheck[0], 1.e-10, fcnm, "error computing lam 1");
    CHKERR(lam[1], lamcheck[1], 1.e-10, fcnm, "error computing lam 2");
    CHKERR(lam[2], lamcheck[2], 1.e-10, fcnm, "error computing lam 3");
    //lam / norm(lam)
    return EXIT_SUCCESS;
}
//============================================================================//
int check_CMT2faultpar(void)
{
    const double M[6] = {3.108304932835845, 3.044425632830430, 3.382269434333724,
                     -4.855033301709626,-1.949280336439431, 1.110527600460120};
    double nu, alpha, N1[3], N2[3], lam[3];
    const double nuRef = 0.371745072651417;
    const double alphaRef = 80.365019327257144;
    const double N1ref[3] = {-0.050351975426048, 0.930119048707825, 0.363790095799137};
    const double N2ref[3] = {-0.977452142939121, 0.046467470073259, 0.205980781843140};
    const double lamref[3] = { 8.802000000000001, 2.584000000000001, -1.851000000000002};
    int ierr;
    int nmt = 1;
    ierr = compearth_CMT2faultpar(nmt, M,
                                  &nu, &alpha, N1, N2, lam);
    if (ierr != 0)
    {
        fprintf(stderr, "Error calling CMT2faultpar\n");
        return EXIT_FAILURE;
    }
    CHKERR(nu,    nuRef,    1.e-10, __func__, "error computing nu");
    CHKERR(alpha, alphaRef, 1.e-10, __func__, "error computing alpha");
    CHKERR(N1[0], N1ref[0], 1.e-10, __func__, "error checking N1[0]");
    CHKERR(N1[1], N1ref[1], 1.e-10, __func__, "error checking N1[1]");
    CHKERR(N1[2], N1ref[2], 1.e-10, __func__, "error checking N1[2]");
    CHKERR(N2[0], N2ref[0], 1.e-10, __func__, "error checking N2[0]");
    CHKERR(N2[1], N2ref[1], 1.e-10, __func__, "error checking N2[1]");
    CHKERR(N2[2], N2ref[2], 1.e-10, __func__, "error checking N2[2]");
    CHKERR(lam[0], lamref[0], 1.e-10, __func__, "error checking lam[0]");
    CHKERR(lam[1], lamref[1], 1.e-10, __func__, "error checking lam[1]");
    CHKERR(lam[2], lamref[2], 1.e-10, __func__, "error checking lam[2]");
    return EXIT_SUCCESS;
}
//============================================================================//
int check_TT2CMT(void)
{
    const double M01[1] = {1.0};
    const double delta1[1] = {0.0};
    const double gamma1[1] = {0.0};
    const double kappa1[1] = {320.0};
    const double theta1[1] = {10.0};
    const double sigma1[1] = {20.0};
    const double lam1Ref[3] = {1.0, 0.0, -1.0};
    const double M61Ref[6] = { 0.116977778440511,
                               0.112364502228240,
                              -0.229342280668750,
                              -0.502322271868946,
                              -0.841048248646006,
                               0.029265111955970};
    const double U1Ref[9] = {-0.434841577578901, -0.515496946820204,
                              0.738360142632131,  0.856848940622339,
                             -0.489063917059259,  0.163175911166535,
                              0.276988619555150,  0.703618776646631,
                              0.654368338007907};
    double M61[6], lam1[3], U91[9];
    int i, ierr, j;
    memset(U91, 0, 9*sizeof(double));
    ierr = compearth_TT2CMT(1, gamma1, delta1, M01, kappa1, theta1, sigma1,
                            M61, lam1, U91);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: error calling TT2CMT\n", __func__);
        return EXIT_FAILURE;
    }
    CHKERR(M61[0], M61Ref[0], 1.e-10, __func__, "error checking M61[0]");
    CHKERR(M61[1], M61Ref[1], 1.e-10, __func__, "error checking M61[1]");
    CHKERR(M61[2], M61Ref[2], 1.e-10, __func__, "error checking M61[2]");
    CHKERR(M61[3], M61Ref[3], 1.e-10, __func__, "error checking M61[3]");
    CHKERR(M61[4], M61Ref[4], 1.e-10, __func__, "error checking M61[4]");
    CHKERR(M61[5], M61Ref[5], 1.e-10, __func__, "error checking M61[5]");

    CHKERR(lam1[0], lam1Ref[0], 1.e-10, __func__, "error checking lam[0]");
    CHKERR(lam1[1], lam1Ref[1], 1.e-10, __func__, "error checking lam[1]");
    CHKERR(lam1[2], lam1Ref[2], 1.e-10, __func__, "error checking lam[2]");

    CHKERR(U91[0], U1Ref[0], 1.e-10, __func__, "error checking U1[0]");
    CHKERR(U91[1], U1Ref[1], 1.e-10, __func__, "error checking U1[1]");
    CHKERR(U91[2], U1Ref[2], 1.e-10, __func__, "error checking U1[2]");
    CHKERR(U91[3], U1Ref[3], 1.e-10, __func__, "error checking U1[3]");
    CHKERR(U91[4], U1Ref[4], 1.e-10, __func__, "error checking U1[4]");
    CHKERR(U91[5], U1Ref[5], 1.e-10, __func__, "error checking U1[5]");
    CHKERR(U91[6], U1Ref[6], 1.e-10, __func__, "error checking U1[6]");
    CHKERR(U91[7], U1Ref[7], 1.e-10, __func__, "error checking U1[7]");
    CHKERR(U91[8], U1Ref[8], 1.e-10, __func__, "error checking U1[8]");
    // At this stage in the game CMT2TT works.  So let's use it to verify
    // its inverse operator.
    int ng = 3;
    double gamLoc[3] = {-29.0, 0.0, 29.0};
    int nd = 3;
    double deltaLoc[3] = {-89, 0.0, 89.0};
    int nk = 5;
    double kappaLoc[5] = {1.0, 90.0, 180.0, 270.0, 359.};
    int ns = 3;
    double sigmaLoc[3] = {-179.0, 1.0, 179.0};
    int nt = 3;
    double thetaLoc[3] = {10.0, 45.0, 80.0};
    int nm = 3;
    double M0loc[3] = {1.0, 2.0, 3.0};
    int nmt = ng*nd*nk*ns*nt*nm;
    double *gamma, *delta, *kappa, *sigma, *theta, *M0, *M, *U, *lam;
    gamma = (double *) calloc((size_t) nmt, sizeof(double));
    delta = (double *) calloc((size_t) nmt, sizeof(double));
    kappa = (double *) calloc((size_t) nmt, sizeof(double));
    sigma = (double *) calloc((size_t) nmt, sizeof(double));
    theta = (double *) calloc((size_t) nmt, sizeof(double));
    M0 = (double *) calloc((size_t) nmt, sizeof(double));
    int ig, id, ik, is, im, imt, it;
    imt = 0;
    for (im=0; im<nm; im++)
    {
        for (ig=0; ig<ng; ig++)
        {
            for (id=0; id<nd; id++)
            {
                for (ik=0; ik<nk; ik++)
                {
                    for (it=0; it<nt; it++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            gamma[imt] = gamLoc[ig];
                            delta[imt] = deltaLoc[id];
                            kappa[imt] = kappaLoc[ik];
                            sigma[imt] = sigmaLoc[is];
                            theta[imt] = thetaLoc[it];
                            M0[imt] = M0loc[im];
                            imt = imt + 1;
                        }
                    }
                }
            }
        }
    }
    M = (double *) calloc((size_t) (6*nmt), sizeof(double));
    lam = (double *) calloc((size_t) (3*nmt), sizeof(double));
    U = (double *) calloc((size_t) (9*nmt), sizeof(double));
    ierr = compearth_TT2CMT(nmt, gamma, delta, M0, kappa, theta, sigma,
                            M, lam, U);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: error calling TT2CMT\n", __func__);
        return EXIT_FAILURE;
    }
/*
printf("\n");
for (i=0; i<nmt; i++)
{
printf("%d\n", 6*i);
 printf("%e %e %e %e %e %e\n", M[6*i], M[6*i+1], M[6*i+2], M[6*i+3], M[6*i+4], M[6*i+5]);
}
getchar();
*/
    bool ldisplay = false;
    double *gammaNew = (double *) calloc((size_t) nmt, sizeof(double));
    double *deltaNew = (double *) calloc((size_t) nmt, sizeof(double));
    double *M0New = (double *) calloc((size_t) nmt, sizeof(double));
    double *kappaNew = (double *) calloc((size_t) nmt, sizeof(double));
    double *thetaNew = (double *) calloc((size_t) nmt, sizeof(double));
    double *sigmaNew = (double *) calloc((size_t) nmt, sizeof(double));
    double *K = (double *) calloc((size_t) (3*nmt), sizeof(double));
    double *N = (double *) calloc((size_t) (3*nmt), sizeof(double));
    double *S = (double *) calloc((size_t) (3*nmt), sizeof(double));
    double *thetadc = (double *) calloc((size_t) nmt, sizeof(double));
    double *lamNew = (double *) calloc((size_t) (3*nmt), sizeof(double));
    double *UNew = (double *) calloc((size_t) (9*nmt), sizeof(double));
    ierr = compearth_CMT2TT(nmt, M, ldisplay,
                            gammaNew, deltaNew, M0New,
                            kappaNew, thetaNew, sigmaNew,
                            K, N, S, thetadc,
                            lamNew, UNew);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling CMT2TT\n", __func__);
        return EXIT_FAILURE;
    }
    double kappa2, theta2, sigma2;
    for (i=0; i<nmt; i++)
    {
        ierr = compearth_auxiliaryPlane(1,
                                        &kappaNew[i], &thetaNew[i], &sigmaNew[i],
                                        &kappa2, &theta2, &sigma2);
        if (fabs(gammaNew[i] - gamma[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute gamma %e %e %d\n",
                    __func__, gamma[i], gammaNew[i], i);
            return -1;
        }
        if (fabs(deltaNew[i] - delta[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute delta %e %e %d\n",
                    __func__, delta[i], deltaNew[i], i);
            return -1;
        }
        if (fabs(M0[i] - M0New[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute M0 %e %e %d\n",
                    __func__, M0[i], M0New[i], i);
            return -1;
        }
        if (fabs(kappaNew[i] - kappa[i]) > 1.e-10 &&
            fabs(kappa2 - kappa[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute kappa %e %e %e %d\n",
                    __func__, kappa[i], kappaNew[i], kappa2, i);
            return -1;
        }
        if (fabs(thetaNew[i] - theta[i]) > 1.e-10 &&
            fabs(theta2 - theta[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute theta %e %e %d\n",
                    __func__, theta[i], thetaNew[i], i);
            return -1;
        }
        if (fabs(sigmaNew[i] - sigma[i]) > 1.e-10 &&
            fabs(sigma2 - sigma[i]) > 1.e-10)
        {
            fprintf(stderr, "%s: Failed to compute sigma %e %e %d\n",
                    __func__, sigma[i], sigmaNew[i], i);
            return -1;
        }
        for (j=0; j<3; j++)
        {
            if (fabs(lamNew[3*i+j] - lam[3*i+j]) > 1.e-10)
            {
                fprintf(stderr, "%s: Failed to compute lam %e %e %d\n",
                        __func__, lam[3*i+j], lamNew[3*i+j], i);
                return -1;
            }
        }
        for (j=0; j<9; j++)
        {
            // Can be off up to a sign factor
            if (fabs(UNew[9*i+j] - U[9*i+j]) > 1.e-10 &&
                fabs(UNew[9*i+j] + U[9*i+j]) > 1.e-10)
            {
                fprintf(stderr, "%s: Failed to compute U %e %e %d\n",
                        __func__, UNew[9*i+j], U[9*i+j], i);
            }
        }
    }
    // Free space
    free(gammaNew);
    free(deltaNew);
    free(M0New);
    free(kappaNew);
    free(thetaNew);
    free(sigmaNew);
    free(K);
    free(N);
    free(S);
    free(thetadc);
    free(lamNew);
    free(UNew);
    free(gamma);
    free(delta);
    free(kappa);
    free(sigma);
    free(theta);
    free(M0);
    free(M);
    free(U);
    free(lam);
    return EXIT_SUCCESS;
}

//============================================================================//
int check_CMT2TT(void)
{
    const char *fcnm = "check_CMT2TT\0";
    const double M[6] = {3.108304932835845, 3.044425632830430,
                         3.382269434333724,-4.855033301709626,
                        -1.949280336439431, 1.110527600460120};
    double U[9], lam[3];
    bool ldisplay;
    double gamma, delta, M0, kappa, sigma, theta, thetadc;
    double K[3], N[3], S[3];
    const double gamma1 = -5.519441015228609;
    const double delta1 = 36.032871051674995;
    const double M01 = 6.617343160211657;
    const double kappa1 = 69.492561954605137;
    const double theta1 = 88.144957350922951;
    const double sigma1 = -79.788954797632840;
    const double thetadc1 = 36.396475670817210;
    const double KR1[3] = {-0.350328975576259, 0.936626718000127, 0};
    const double NR1[3] = {0.936135854045508, 0.350145376429380, 0.032370945856051};
    const double SR1[3] = {-0.032265105849630, 0.177200864349715, -0.983645676358224};
    const double lam1[3] = {8.802000000000000, 2.584000000000003, -1.851000000000001};
    const double UR1[9] = {-0.639133135365464,-0.372890102888132, 0.672652812709492,
                            0.350155145207092,-0.919781533321278,-0.177181560118875,
                            0.684762885649414, 0.122290237260530, 0.718432243365964};
    // Test 2; problem child of horizontal fault
    const double M2[6] = {0, 0, 0, -sqrt(3)/2., .5, 0};
    const double gamma2 = 0.0; const double delta2 = 0.0;
    const double M02 = 1.0; const double kappa2 = 300.0;
    const double theta2 = 90.0; const double sigma2 = 90.0;
    const double thetadc2 = 1.207418269725733e-06;
    const double KR2[3] = {-0.500000000000000, -0.866025403784439, 0};
    const double NR2[3] = {-0.866025403784439, 0.500000000000000,  0};
    const double SR2[3] = {0, 0, 1};
    const double lam2[3] = {1, 0, -1};
    const double UR2[9] = {-0.612372435695794,0.353553390593274,0.707106781186547,
                            0.500000000000000,0.866025403784439,0,
                           -0.612372435695794,0.353553390593274,-0.707106781186547};

    int i, ierr;
    ldisplay = false;
    ierr = compearth_CMT2TT(1, M, ldisplay,
                            &gamma, &delta, &M0, &kappa, &theta, &sigma,
                            K, N, S, &thetadc,
                            lam, U);
    if (ierr != 0)
    {
        printf("%s: Error calling CMT2TT test 1\n", fcnm);
        return EXIT_FAILURE;
    }
    CHKERR(gamma, gamma1, 1.e-10, fcnm, "error computing gamma 1");
    CHKERR(delta, delta1, 1.e-10, fcnm, "error computing delta 1");
    CHKERR(M0, M01, 1.e-10, fcnm, "error computing M0 1");
    CHKERR(kappa, kappa1, 1.e-10, fcnm, "error computing kappa 1");
    CHKERR(theta, theta1, 1.e-10, fcnm, "error computing theta 1");
    CHKERR(sigma, sigma1, 1.e-10, fcnm, "error computing sigma 1");
    CHKERR(thetadc, thetadc1, 1.e-10, fcnm, "error computing thetadc 1");
    for (i=0; i<3; i++)
    {
        if (fabs(KR1[i] - K[i]) > 1.e-10)
        {
            printf("%s: Error computing K 1\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(NR1[i] - N[i]) > 1.e-10)
        {
            printf("%s: Error computing N 1\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(SR1[i] - S[i]) > 1.e-10)
        {
            printf("%s: Error computing S 1\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(lam1[i] - lam[i]) > 1.e-10)
        {
            printf("%s: Error computing lam 1\n", fcnm);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - UR1[i]) > 1.e-10)
        {
            printf("%s: Error computing U 1\n", fcnm);
            return EXIT_FAILURE;
        }
    }

    ierr = compearth_CMT2TT(1, M2, ldisplay,
                            &gamma, &delta, &M0, &kappa, &theta, &sigma,
                            K, N, S, &thetadc,
                            lam, U);
    if (ierr != 0)
    {
        printf("%s: Error calling CMT2TT horizontal fault\n", fcnm);
        return EXIT_FAILURE;
    }
    CHKERR(gamma, gamma2, 1.e-10, fcnm, "error computing gamma 2");
    CHKERR(delta, delta2, 1.e-10, fcnm, "error computing delta 2");
    CHKERR(M0, M02, 1.e-10, fcnm, "error computing M0 2");
    CHKERR(kappa, kappa2, 1.e-10, fcnm, "error computing kappa 2");
    CHKERR(theta, theta2, 1.e-10, fcnm, "error computing theta 2");
    CHKERR(sigma, sigma2, 1.e-10, fcnm, "error computing sigma 2");
    CHKERR(thetadc, thetadc2, 1.e-10, fcnm, "error computing thetadc 2");
    for (i=0; i<3; i++)
    {
        if (fabs(KR2[i] - K[i]) > 1.e-10)
        {
            printf("%s: Error computing K 2\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(NR2[i] - N[i]) > 1.e-10)
        {
            printf("%s: Error computing N 2\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(SR2[i] - S[i]) > 1.e-10)
        {
            printf("%s: Error computing S 2\n", fcnm);
            return EXIT_FAILURE;
        }
        if (fabs(lam2[i] - lam[i]) > 1.e-10)
        {
            printf("%s: Error computing lam 2\n", fcnm);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - UR2[i]) > 1.e-10)
        {
            printf("%s: Error computing U 2\n", fcnm);
            return EXIT_FAILURE;
        }
    }
    return EXIT_SUCCESS;
}
//============================================================================//
int check_CMTdecom(void)
{
    const char *fcnm = "check_CMTdecom\0";
    const double M[6] = {3.108304932835845, 3.044425632830430,
                         3.382269434333724,-4.855033301709626,
                        -1.949280336439431, 1.110527600460120};
    double U[9], lam[3];
    double lamRef1[3] = {8.8020000000e+00, 2.5840000000e+00, -1.8510000000e+00};
    double lamRef2[3] = {-1.851000000000002, 2.584000000000001, 8.802000000000001};
    double *lamRef3 = lamRef1;
    double *lamRef4 = lamRef2;
    double URef1[9] = {-0.672652812709492, 0.639133135365464, 0.372890102888132,
                        0.177181560118875,-0.350155145207092, 0.919781533321279,
                        0.718432243365964, 0.684762885649414, 0.122290237260530};
    double URef2[9] = { 0.718432243365964, 0.684762885649414, 0.122290237260530,
                       -0.177181560118875, 0.350155145207092,-0.919781533321279,
                       -0.672652812709492, 0.639133135365464, 0.372890102888132};
    double *URef3 = URef1;
    double *URef4 = URef2;
    int i, ierr;
    ierr = compearth_CMTdecom(1,  M, 1, lam, U);
    if (ierr != 0)
    {
        printf("%s: Error 1 calling CMTdecom\n", fcnm);
        return EXIT_FAILURE;
    }
    for (i=0; i<3; i++)
    {
        if (fabs(lam[i] - lamRef1[i]) > 1.e-13)
        {
            printf("%s: lam 1 is wrong %e %e\n", fcnm, lam[i], lamRef1[i]);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - URef1[i]) > 1.e-13)
        {
            printf("%s: U 1 is wrong %d %e %e\n", fcnm, i, U[i], URef1[i]);
            return EXIT_FAILURE;
        }
    }

    ierr = compearth_CMTdecom(1,  M, 2, lam, U);
    if (ierr != 0)
    {
        printf("%s: Error 2 calling CMTdecom\n", fcnm);
        return EXIT_FAILURE;
    }
    for (i=0; i<3; i++)
    {
        if (fabs(lam[i] - lamRef2[i]) > 1.e-13)
        {
            printf("%s: lam 2 is wrong %e %e\n", fcnm, lam[i], lamRef2[i]);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - URef2[i]) > 1.e-13)
        {
            printf("%s: U 2 is wrong %d %e %e\n", fcnm, i, U[i], URef2[i]);
            return EXIT_FAILURE;
        }
    }

    ierr = compearth_CMTdecom(1,  M, 3, lam, U);
    if (ierr != 0)
    {
        printf("%s: Error calling CMTdecom 3\n", fcnm);
        return EXIT_FAILURE;
    }
    for (i=0; i<3; i++)
    {
        if (fabs(lam[i] - lamRef3[i]) > 1.e-13)
        {
            printf("%s: lam 3 is wrong %e %e\n", fcnm, lam[i], lamRef3[i]);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - URef3[i]) > 1.e-13)
        {
            printf("%s: U 3 is wrong %d %e %e\n", fcnm, i, U[i], URef3[i]);
            return EXIT_FAILURE;
        }
    }

    ierr = compearth_CMTdecom(1,  M, 4, lam, U);
    if (ierr != 0)
    {
        printf("%s: Error calling CMTdecom 4\n", fcnm);
        return EXIT_FAILURE;
    }
    for (i=0; i<3; i++)
    {
        if (fabs(lam[i] - lamRef4[i]) > 1.e-13)
        {
            printf("%s: lam 4 is wrong %e %e\n", fcnm, lam[i], lamRef4[i]);
            return EXIT_FAILURE;
        }
    }
    for (i=0; i<9; i++)
    {
        if (fabs(U[i] - URef4[i]) > 1.e-13)
        {
            printf("%s: U 4 is wrong %d %e %e\n", fcnm, i, U[i], URef4[i]);
            return EXIT_FAILURE;
        }
    }
    return EXIT_SUCCESS;
}
//============================================================================//
int check_lam2lune(void)
{
    const char *fcnm = "check_lam2lune\0";
    const double lamRef1[6] = {
        -7.9621992143e+15,-8.2644837610e+15,-8.2644837610e+15,
        -7.5915701171e+15,-8.4370629711e+15,-8.4370629711e+15};
    double *gamma0, *gamma, *delta, *delta0, *lam, *M0, *M00;
    double *thetadc, *lamdev, *lamiso;
    int i, ib, ierr, ig, indx;
    const int ng = 100;
    const int nb = 100;
    int nmt = ng*nb;
    const double dg = (30.0 - -30.0)/(double) (ng - 1);
    const double db = (89.0 - -89.0)/(double) (nb - 1);
    gamma0 = (double *) calloc((size_t) nmt, sizeof(double));
    delta0 = (double *) calloc((size_t) nmt, sizeof(double));
    M00    = (double *) calloc((size_t) nmt, sizeof(double));
    gamma  = (double *) calloc((size_t) nmt, sizeof(double));
    delta  = (double *) calloc((size_t) nmt, sizeof(double));
    M0     = (double *) calloc((size_t) nmt, sizeof(double));
    lam  = (double *) calloc((size_t) (3*nmt), sizeof(double));
    for (ig=0; ig<ng; ig++)
    {
        for (ib=0; ib<nb; ib++)
        {
            indx = ig*nb + ib;
            gamma0[indx] =-30.0 + (double) ig*dg;
            delta0[indx] =-89.0 + (double) ib*db; // use latitude instead of colat
            M00[indx] = 1.e16;
        }
    }
    compearth_lune2lam(nmt, gamma0, delta0, M00, lam);
    // Verify the first two computations are right
    for (i=0; i<6; i++)
    {
        if (fabs(lam[i] - lamRef1[i])/1.e16 > 1.e-10)
        {
            printf("%s: error in lune2lam: %.10e %.10e\n",
                   fcnm, lam[i], lamRef1[i]);
            return EXIT_FAILURE;
        }
    }
    // Now perform the inverse operation
    thetadc = NULL;
    lamdev = NULL;
    lamiso = NULL;
    ierr = compearth_lam2lune(nmt, lam,
                              gamma, delta, M0,
                              thetadc, lamdev, lamiso);
    if (ierr != 0)
    {
        printf("%s: error calling lam2lune\n", fcnm);
        return EXIT_FAILURE;
    }
    // See how close i was
    for (i=0; i<nmt; i++)
    {
        if (fabs(gamma[i] - gamma0[i]) > 1.e-12)
        {
            printf("%s: gamma[%d] different %e %e %e\n",
                   __func__, i, gamma[i], gamma0[i],
                   fabs(gamma[i] - gamma0[i]));
            return EXIT_FAILURE;
        }
        if (fabs(delta[i] - delta0[i]) > 1.e-12)
        {
            printf("%s: delta[%d] different %e %e %e\n",
                   __func__, i, delta[i], delta0[i],
                   fabs(delta[i] - delta0[i]));
            return EXIT_FAILURE;
        }
        if (fabs(M0[i] - M00[i])/1.e15 > 1.e-12)
        {
            printf("%s: M0[%d] different %e %e %e\n",
                   __func__, i, M0[i], M00[i], fabs(M0[i] - M00[i]));
           return EXIT_FAILURE;
        }
    }
    free(gamma0);
    free(delta0);
    free(M00);
    free(gamma);
    free(delta);
    free(M0);
    free(lam);
    return EXIT_SUCCESS;
}

int check_CMT2omega(void)
{
    const double M1[6] = {1, 0, -1, 0, 0, 0};
    const double M2[6] = {1, 2,  3, 4, 5, 6};
    double *Mvec1, *Mvec2, *omega, *refOmega, om1;
    const double omRef1 = 96.263952719927232;
    int i, ierr, j, nmt1, nmt2;
    ierr = compearth_CMT2omega(1, M1, 1, M2, &om1);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling CMT2omega 1\n", __func__);
        return EXIT_FAILURE;
    }
    CHKERR(om1, omRef1, 1.e-10, __func__, "error computing CMT2Omega1");
    nmt2 = 80;
    Mvec1 = (double *) calloc((size_t) (6*nmt2), sizeof(double));
    Mvec2 = (double *) calloc((size_t) (6*nmt2), sizeof(double));
    omega = (double *) calloc((size_t) nmt2, sizeof(double));
    refOmega = (double *) calloc((size_t) nmt2, sizeof(double));
    for (i=0; i<nmt2; i++)
    {
        cblas_dcopy(6, M2, 1, &Mvec2[6*i], 1);
    }
    // one to many
    ierr = compearth_CMT2omega(1, M1, nmt2, Mvec2, omega);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling CMT2omega 2\n", __func__);
        return EXIT_FAILURE;
    }
    for (i=0; i<nmt2; i++)
    {
         CHKERR(omega[i], omRef1, 1.e-10, __func__,
                "error computing CMT2Omega2");
    }
    // flip roles
    ierr = compearth_CMT2omega(nmt2, Mvec2, 1, M1, omega);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling CMT2omega 3\n", __func__);
        return EXIT_FAILURE;
    }
    for (i=0; i<nmt2; i++)
    {
         CHKERR(omega[i], omRef1, 1.e-10, __func__,
                "error computing CMT2Omega3");
    }
    // many to many
    nmt1 = nmt2;
    srand(4093);
    for (i=0; i<nmt1; i++)
    {
        for (j=0; j<6; j++)
        {
            Mvec1[6*i+j] = ((double) (rand())/RAND_MAX - 0.5)*2.0;
            Mvec2[6*i+j] = ((double) (rand())/RAND_MAX - 0.5)*2.0;
        }
        compearth_CMT2omega(1, &Mvec1[6*i], 1, &Mvec2[6*i], &refOmega[i]);
    }
    ierr = compearth_CMT2omega(nmt1, Mvec1, nmt2, Mvec2, omega);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling CMT2omega 4\n", __func__);
        return EXIT_FAILURE;
    }
    for (i=0; i<nmt2; i++)
    {
         CHKERR(omega[i], refOmega[i], 1.e-10, __func__,
                "error computing CMT2Omega4");
    }

    free(Mvec1);
    free(Mvec2);
    free(omega);
    free(refOmega);
    return EXIT_SUCCESS;
}

int check_normal2strdip(void)
{
    const double Xin[3] = {-0.171010071662835, 0.969846310392954,
                            0.173648177666930};
    double Xout[2];
    compearth_normal2strdip(1, Xin, Xout);
    CHKERR(Xout[0], 350.0, 1.e-10, __func__, "error computing Xout[0]");
    CHKERR(Xout[1], 80.0,  1.e-10, __func__, "error computing Xout[1]");
    return EXIT_SUCCESS;
}

int check_Uorth(void)
{
    const double T[9] = {0.998800000000000, 0.047900000000000, 0.017300000000000,
             -0.048200000000000, 0.998800000000000, 0.015900000000000,
             -0.016500000000000,-0.016700000000000, 0.999800000000000};
    const double Tref[9] = { 0.998702105045258, 0.047907730990670, 0.017290306229078,
             -0.048183355871821, 0.998712063377127, 0.015892724182219,
             -0.016506653055635, -0.016705202073855,  0.999724195280165};
    const double detInRef1 = 1.000261915677000;
    const double detOutRef1 = 1.0;
    double Tout[9], detIn, detOut;
    int i, j, ierr;
    ierr = compearth_Uorth(1, CE_ORTH_SVD, T, Tout, &detIn, &detOut);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error calling Uorth_svd\n", __func__);
        return EXIT_FAILURE;
    }
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            if (fabs(Tout[3*j+i] - Tref[3*j+i]) > 1.e-14)
            {
                printf("%s: Failed calculating T(%d,%d)=%f ",
                       __func__, i+1,j+1, Tout[3*j+i] - Tref[3*j+i]);
                return EXIT_FAILURE;
            }
        }
    }
    CHKERR(detIn, detInRef1, 1.e-10, __func__, "error computing detIn1");
    CHKERR(detOut, detOutRef1, 1.e-10, __func__, "error computing detOUt1");
    return EXIT_SUCCESS;
}

ext: c


text:
#include <stdio.h>
#include <gsl/gsl_block.h>

int
main (void)
{
  gsl_block * b = gsl_block_alloc (100);

  printf ("length of block = %zu\n", b->size);
  printf ("block data address = %p\n", b->data);

  gsl_block_free (b);
  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include "rev/gl/Context.h"
#include "rev/gl/Resource.h"

#include <gsl/span>

namespace rev {

using Buffer = Resource<singleCreate<gl::genBuffers>, singleDestroy<gl::deleteBuffers>>;

} // namespace rev

ext: h

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
/// @file UniquePtr.h
/// @author Braxton Salyer <braxtonsalyer@gmail.com>
/// @brief This includes Hyperion's `constexpr` equivalents and extensions to the C++ standard
/// library's `std::unique_ptr<T, Deleter>`
/// @version 0.1
/// @date 2021-11-05
///
/// 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/Concepts.h>
#include <Hyperion/HyperionDef.h>
#include <Hyperion/memory/CompressedPair.h>
#include <gsl/gsl>

/// @ingroup memory
/// @{
///	@defgroup UniquePtr UniquePtr
/// Hyperion provides a `constexpr` equivalent to the C++ standard library's
/// `std::unique_ptr<T, Deleter>` in our `UniquePtr`, with additional functionality like allocator
/// aware factory functions (`allocate_unique`).
///
/// Example:
///
/// auto ptr = hyperion::make_unique<i32>(42);
/// ptr.reset(24);
///
/// auto ptr2 = hyperion::allocate_unique<i32>(std::allocator<i32>(), 42);
/// ptr2.reset(24);
/// @headerfile "Hyperion/memory/UniquePtr.h"
/// @}

namespace hyperion {

	/// @brief Default Deleter type for Hyperion smart pointers
	/// This is the default deleter type used by Hyperion smart pointers
	/// such as `hyperion::UniquePtr<T, Deleter>`
	///
	/// @tparam T The type to handle deletion for (eg, `i32`, `std::string`)
	///
	/// # Requirements
	/// - `concepts::Deletable<T>`: T must be a deletable type
	/// - `concepts::NotFunction<T>`: T must NOT be a function type
	/// @ingroup memory
	/// @headerfile "Hyperion/memory/UniquePtr.h"
	template<typename T>
	requires concepts::Deletable<T> && concepts::NotFunction<T>
	struct DefaultDeleter {

		/// @brief Default constructs a `DefaultDeleter`
		/// @ingroup memory
		constexpr DefaultDeleter() noexcept = default;
		/// @brief Copy-constructs a `DefaultDeleter` from the given one
		/// @ingroup memory
		constexpr DefaultDeleter(const DefaultDeleter&) noexcept = default;
		/// @brief Move-constructs a `DefaultDeleter` from the given one
		/// @ingroup memory
		constexpr DefaultDeleter(DefaultDeleter&&) noexcept = default;
		/// @brief Destructs a `DefaultDeleter`
		/// @ingroup memory
		constexpr ~DefaultDeleter() noexcept = default;
		/// @brief Constructs a `DefaultDeleter<T>` from a `DefaultDeleter<U>`
		///
		/// @tparam U - The type managed by the given deleter
		///
		/// @param deleter - The deleter to construct this from
		///
		/// # Requirements
		/// - `concepts::Convertible<U*, T*>`: The pointer type of `deleter` must be convertible two
		/// the pointer type of `DefaultDeleter<T>` in order to construct a `DefaultDeleter` from it
		/// @ingroup memory
		template<typename U>
		requires concepts::Convertible<U*, T*>
		explicit DefaultDeleter([[maybe_unused]] const DefaultDeleter<U>& deleter) noexcept {
		}

		/// @brief Deletes the data at the given pointer
		///
		/// @param ptr - The pointer to the data to delete
		/// @ingroup memory
		constexpr auto
		operator()(gsl::owner<T*> ptr) const noexcept(concepts::NoexceptDeletable<T>) {
			delete ptr;
		}

		/// @brief Copy-assigns this `DefaultDeleter` from the given one
		/// @ingroup memory
		constexpr auto operator=(const DefaultDeleter&) noexcept -> DefaultDeleter& = default;
		/// @brief Move-assigns this `DefaultDeleter` from the given one
		/// @ingroup memory
		constexpr auto operator=(DefaultDeleter&&) noexcept -> DefaultDeleter& = default;
	};

	template<typename T>
	requires concepts::Deletable<T[]> && concepts::NotFunction<T>
	struct DefaultDeleter<T[]> { // NOLINT
		constexpr DefaultDeleter() noexcept = default;
		constexpr DefaultDeleter(const DefaultDeleter&) noexcept = default;
		constexpr DefaultDeleter(DefaultDeleter&&) noexcept = default;
		constexpr ~DefaultDeleter() noexcept = default;
		template<typename U>
		requires concepts::Convertible<U (*)[], T (*)[]>							 // NOLINT
		explicit DefaultDeleter([[maybe_unused]] const DefaultDeleter<U[]>& deleter) // NOLINT
			noexcept {
		}

		constexpr auto
		operator()(gsl::owner<T*> ptr) const noexcept(concepts::NoexceptDeletable<T[]>) // NOLINT
		{
			delete[] ptr;
		}

		constexpr auto operator=(const DefaultDeleter&) noexcept -> DefaultDeleter& = default;
		constexpr auto operator=(DefaultDeleter&&) noexcept -> DefaultDeleter& = default;
	};

	/// @brief `UniquePtr<T, Deleter>` is Hyperion's `constexpr` equivalent to the standard
	/// library's `std::unique_ptr<T, Deleter>`
	///
	/// `UniquePtr<T, Deleter>` is a `constexpr` implementation of a unique owning pointer, with
	/// almost complete semantic equivalence with `std::unique_ptr<T, Deleter>`.
	/// It is a drop-in replacement in almost all situations(1).
	/// Hyperion also provides extended functionality for `UniquePtr<T, Deleter>`
	/// over `std::unique_ptr<T, Deleter>`, eg: an allocator-aware factory function set in
	/// `allocate_unique`
	///
	/// @tparam T - The type to store in the `UniquePtr`
	/// @tparam Deleter - The deleter type to handle deletion of the stored `T`, by default, this is
	/// `hyperion::DefaultDeleter<T>`. This type only needs to be explicitly provided when custom
	/// deletion or resource-freeing strategies are required.
	///
	/// # Requirements
	/// - `!std::is_rvalue_reference_v<Deleter>`: `Deleter` cannot be an r-value reference type
	///
	/// @note(1): The exception is in situations where destructor ordering is required and code was
	/// compiled with Clang. `UniquePtr<T, Deleter>` makes use of Clang's `trivial_abi` attribute,
	/// which allows `UniquePtr<T, Deleter>` to be passed via register in function calls
	/// (and thus removes its overhead compared to a raw pointer), but moves the call to its
	/// destructor from the calling code to the callee. The result is that the destructor for
	/// `UniquePtr<T, Deleter>` will not nest in the same way as `std::unique_ptr<T, Deleter>`, and
	/// code relying on that behavior may experience issues. This is identical behavior to
	/// `std::unique_ptr<T, Deleter>` when `_LIBCPP_ABI_ENABLE_UNIQUE_PTR_TRIVIAL_ABI` is defined
	/// @ingroup UniquePtr
	/// @headerfile "Hyperion/memory/UniquePtr.h"
	template<typename T, typename Deleter = DefaultDeleter<T>>
	requires concepts::NotRValueReference<Deleter>
	class HYPERION_TRIVIAL_ABI UniquePtr {
	  public:
		/// @brief The element type allocated in the `UniquePtr`
		/// @ingroup UniquePtr
		using element_type = T;
		/// @brief The deleter type used by the `UniquePtr`
		/// @ingroup UniquePtr
		using deleter_type = Deleter;
		/// @brief The pointer type stored in the `UniquePtr`
		/// @ingroup UniquePtr
		using pointer = T*;
		/// @brief The pointer type stored in the `UniquePtr`, as a pointer to const
		/// @ingroup UniquePtr
		using pointer_to_const = const T*;

		/// @brief Constructs a default `UniquePtr<T, Deleter>`, managing no pointer
		///
		/// # Requirements
		/// - `concepts::NoexceptDefaultConstructible<deleter_type>`: `deleter_type` must be
		/// noexcept default constructible in order to default construct a `UniquePtr`
		/// - `concepts::NotPointer<deleter_type>`: `deleter_type` must not be a pointer
		/// in order to default construct a `UniquePtr`
		/// @ingroup UniquePtr
		constexpr UniquePtr() noexcept requires concepts::NoexceptDefaultConstructible<
			deleter_type> && concepts::NotPointer<deleter_type>
			: m_ptr(pointer(), DefaultInitTag<deleter_type>()) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing no pointer
		///
		/// @param ptr - The `nullptr` signaling this `UniquePtr` should manage no pointer
		///
		/// # Requirements
		/// - `concepts::NoexceptDefaultConstructible<deleter_type>`: `deleter_type` must be
		/// noexcept default constructible in order to construct a `UniquePtr` from only a `nullptr`
		/// - `concepts::NotPointer<deleter_type>`: `deleter_type` must not be a pointer
		/// in order to construct a `UniquePtr` from only a `nullptr`
		/// @ingroup UniquePtr
		constexpr UniquePtr(std::nullptr_t ptr) noexcept // NOLINT
			requires concepts::NoexceptDefaultConstructible<deleter_type> && concepts::NotPointer<
				deleter_type> : m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		///
		/// # Requirements
		/// - `concepts::NoexceptDefaultConstructible<deleter_type>`: `deleter_type` must be
		/// noexcept default constructible in order to construct a `UniquePtr` from only a `pointer`
		/// - `concepts::NotPointer<deleter_type>`: `deleter_type` must not be a pointer
		/// in order to construct a `UniquePtr` from only a `pointer`
		/// @ingroup UniquePtr
		explicit constexpr UniquePtr(pointer ptr) noexcept requires concepts::
			NoexceptDefaultConstructible<deleter_type> && concepts::NotPointer<deleter_type>
			: m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @tparam U - The type of the pointer to manage in this `UniquePtr`
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		///
		/// # Requirements
		/// - `concepts::Convertible<U, pointer>`: `U` must be convertible to `pointer` to construct
		/// a `UniquePtr` from it
		/// - `concepts::NoexceptDefaultConstructible<deleter_type>`: `deleter_type` must be
		/// noexcept default constructible in order to construct a `UniquePtr` from only a `pointer`
		/// - `concepts::NotPointer<deleter_type>`: `deleter_type` must not be a pointer
		/// in order to construct a `UniquePtr` from only a `pointer`
		/// @ingroup UniquePtr
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NoexceptDefaultConstructible<
			deleter_type> && concepts::NotPointer<deleter_type>
		explicit constexpr UniquePtr(U ptr) noexcept : m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		/// @param deleter - The deleter to delete `ptr` with when the destructor of this
		/// `UniquePtr` is called
		///
		/// @ingroup UniquePtr
		constexpr UniquePtr(pointer ptr, const Deleter& deleter) noexcept requires
			concepts::NoexceptCopyConstructible<Deleter> : m_ptr(ptr, deleter) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		/// @param deleter - The deleter to delete `ptr` with when the destructor of this
		/// `UniquePtr` is called
		///
		/// @ingroup UniquePtr
		constexpr UniquePtr(pointer ptr, Deleter&& deleter) noexcept requires
			concepts::NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
			: m_ptr(ptr, std::move(deleter)) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @tparam U - The type of the pointer to manage in this `UniquePtr`
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		/// @param deleter - The deleter to delete `ptr` with when the destructor of this
		/// `UniquePtr` is called
		///
		/// # Requirements
		/// - `concepts::Convertible<U, pointer>`: `U` must be convertible to `pointer` to construct
		/// a `UniquePtr` from it
		/// @ingroup UniquePtr
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NotSame<U, pointer> && concepts::
			NoexceptCopyConstructible<Deleter>
		constexpr UniquePtr(U ptr, const Deleter& deleter) noexcept : m_ptr(ptr, deleter) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing the given pointer
		///
		/// @tparam U - The type of the pointer to manage in this `UniquePtr`
		///
		/// @param ptr - The pointer to manage with this `UniquePtr`
		/// @param deleter - The deleter to delete `ptr` with when the destructor of this
		/// `UniquePtr` is called
		///
		/// # Requirements
		/// - `concepts::Convertible<U, pointer>`: `U` must be convertible to `pointer` to construct
		/// a `UniquePtr` from it
		/// @ingroup UniquePtr
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NotSame<U, pointer> && concepts::
			NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
		constexpr UniquePtr(U ptr, Deleter&& deleter) noexcept : m_ptr(ptr, std::move(deleter)) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing no pointer
		///
		/// @param ptr - The `nullptr` signaling this `UniquePtr` should manage no pointer
		/// @param deleter - The deleter to delete a managed pointer with when the destructor of
		/// this `UniquePtr` is called
		///
		/// @ingroup UniquePtr
		constexpr UniquePtr(std::nullptr_t ptr, const Deleter& deleter) noexcept requires
			concepts::NoexceptCopyConstructible<Deleter> : m_ptr(ptr, deleter) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` managing no pointer
		///
		/// @param ptr - The `nullptr` signaling this `UniquePtr` should manage no pointer
		/// @param deleter - The deleter to delete a managed pointer with when the destructor of
		/// this `UniquePtr` is called
		///
		/// @ingroup UniquePtr
		constexpr UniquePtr(std::nullptr_t ptr, Deleter&& deleter) noexcept requires
			concepts::NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
			: m_ptr(ptr, std::move(deleter)) {
		}
		/// @brief Constructs a `UniquePtr<T, Deleter>` from the given moved `UniquePtr<U, D>`
		///
		/// @param ptr - The `UniquePtr` to construct this one from
		///
		/// # Requirements
		/// - `concepts::NoexceptConstructibleFrom<deleter_type,
		/// decltype((std::forward<D>(ptr.get_deleter())))`: `deleter_type` must be noexcept
		/// constructible from the deleter type returned by `ptr.get_deleter()`
		/// - `concepts::Convertible<typename UniquePtr<U, D>::pointer, pointer>`: The pointer type
		/// of `UniquePtr<U, D>` must be convertible to `pointer`
		/// - `concepts::Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type>`: The
		/// deleter type of `UniquePtr<U, D>` must be convertible to `deleter_type`
		/// @ingroup UniquePtr
		template<typename U, typename D>
		explicit constexpr UniquePtr(UniquePtr<U, D>&& ptr) noexcept requires concepts::
			NoexceptConstructibleFrom<deleter_type,
									  decltype((std::forward<D>(ptr.get_deleter())))> && concepts::
				Convertible<typename UniquePtr<U, D>::pointer, pointer> && concepts::
					Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type>
			: m_ptr(ptr.release(), std::forward<D>(ptr.get_deleter())) {
		}
		/// @brief `UniquePtr` cannot be copied
		/// @ingroup UniquePtr
		UniquePtr(const UniquePtr&) = delete;
		/// @brief Move Constructs a `UniquePtr` from the given one
		///
		/// @param ptr - The `UniquePtr` to move
		///
		/// # Requirements
		/// - `concepts::NoexceptMoveConstructible<deleter_type>`: `deleter_type` must be noexcept
		/// move constructible in order to move construct a `UniquePtr`
		/// @ingroup UniquePtr
		constexpr UniquePtr(
			UniquePtr&& ptr) noexcept requires concepts::NoexceptMoveConstructible<deleter_type>
			: m_ptr(ptr.release(), std::forward<deleter_type>(ptr.get_deleter())) {
		}
		/// @brief `UniquePtr` destructor
		///
		/// # Requirements
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		/// @ingroup UniquePtr
		constexpr ~UniquePtr() noexcept
			requires(noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset();
		}

		/// @brief Releases ownership of the managed pointer and returns it
		///
		/// @return The managed pointer
		/// @ingroup UniquePtr
		[[nodiscard]] inline constexpr auto release() noexcept -> pointer {
			auto* ptr = m_ptr.first();
			m_ptr.first() = pointer();
			return ptr;
		}

		/// @brief Deletes the currently managed pointer, if any, and begins managing the given one
		///
		/// @param ptr - The new pointer to manage with this `UniquePtr`
		///
		/// # Requirements
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		inline constexpr auto reset(pointer ptr = pointer()) noexcept -> void requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			gsl::owner<pointer> tmp = m_ptr.first(); // NOLINT
			m_ptr.first() = ptr;
			m_ptr.second()(tmp);
		}

		/// @brief Deletes the currently managed pointer, if any, and begins managing the given one
		///
		/// @tparam U - The type of the new pointer to manage
		///
		/// @param ptr - The new pointer to manage with this `UniquePtr`
		///
		/// # Requirements
		/// - `concepts::Convertible<U, pointer>`: `U` must be convertible to `pointer` in order for
		/// this `UniquePtr` to manage a `U`
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		/// @ingroup UniquePtr
		template<typename U>
		requires concepts::Convertible<U, pointer>
		inline constexpr auto reset(U ptr) noexcept -> void requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			gsl::owner<pointer> tmp = m_ptr.first(); // NOLINT
			m_ptr.first() = ptr;
			m_ptr.second()(tmp);
		}

		/// @brief Swaps the managed pointer and deleter of this `UniquePtr` with the given one
		///
		/// @param ptr - The `UniquePtr` to swap with
		///
		/// # Requirements
		/// - `std::is_nothrow_swappable_v<CompressedPair<pointer, deleter_type>>`:
		/// `CompressedPair<pointer, deleter_type> must be noexcept swappable in order to
		/// swap two `UniquePtr`s. (`CompressedPair` is used internally to store the managed pointer
		/// and deleter, to allow for
		/// [Empty Base Class Optimization](https://en.cppreference.com/w/cpp/language/ebo) )
		/// @ingroup UniquePtr
		inline constexpr auto swap(UniquePtr& ptr) noexcept
			-> void requires concepts::NoexceptSwappable<CompressedPair<pointer, deleter_type>> {
			m_ptr.swap(ptr.m_ptr);
		}

		/// @brief Returns the managed pointer
		///
		/// Returns the managed pointer. This does not release ownership of the pointer, it only
		/// provides unmanaged access to it. Use of the returned pointer after the lifetime of this
		/// `UniquePtr` has ended results in undefined behavior.
		///
		/// @return The managed pointer
		/// @ingroup UniquePtr
		[[nodiscard]] inline constexpr auto get() noexcept -> pointer {
			return m_ptr.first();
		}

		/// @brief Returns the managed pointer
		///
		/// Returns the managed pointer. This does not release ownership of the pointer, it only
		/// provides unmanaged access to it. Use of the returned pointer after the lifetime of this
		/// `UniquePtr` has ended results in undefined behavior.
		///
		/// @return The managed pointer
		/// @ingroup UniquePtr
		[[nodiscard]] inline constexpr auto get() const noexcept -> pointer_to_const {
			return m_ptr.first();
		}

		/// @brief Returns the associated deleter
		///
		/// @return The associated deleter
		/// @ingroup UniquePtr
		[[nodiscard]] inline constexpr auto get_deleter() const noexcept -> const deleter_type& {
			return m_ptr.second();
		}

		/// @brief Returns the associated deleter
		///
		/// @return The associated deleter
		/// @ingroup UniquePtr
		[[nodiscard]] inline constexpr auto get_deleter() noexcept -> deleter_type& {
			return m_ptr.second();
		}

		/// @brief Converts this `UniquePtr` to a `bool`
		///
		/// Returns `true` if the managed pointer is not null, `false` otherwise.
		/// @return this, as a `bool`
		/// @ingroup UniquePtr
		explicit constexpr operator bool() const noexcept {
			return m_ptr.first() != nullptr;
		}

		constexpr auto operator*() const -> typename std::add_lvalue_reference_t<T> {
			return *(m_ptr.first());
		}

		constexpr auto operator->() const noexcept -> pointer {
			return m_ptr.first();
		}

		inline constexpr auto operator==(const UniquePtr& ptr) const noexcept -> bool {
			return m_ptr.first() == ptr.first();
		}

		inline constexpr auto operator==(std::nullptr_t) const noexcept -> bool {
			return m_ptr.first() == nullptr;
		}

		auto operator=(const UniquePtr&) -> UniquePtr& = delete;
		/// @brief Move-assigns this `UniquePtr` from the given one
		///
		/// @param ptr - The `UniquePtr` to move into this one
		///
		/// # Requirements
		/// - `concepts::NoexceptMoveAssignable<deleter_type>`: `deleter_type` must be noexcept move
		/// assignable in order to move assign a `UniquePtr`
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		///
		/// @return this
		/// @ingroup UniquePtr
		constexpr auto operator=(UniquePtr&& ptr) noexcept
			-> UniquePtr& requires concepts::NoexceptMoveAssignable<deleter_type> &&(
				noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			if(this == &ptr) {
				return *this;
			}

			reset(ptr.release());
			m_ptr.second() = std::forward<deleter_type>(ptr.get_deleter());

			return *this;
		}

		/// @brief Move-assigns this `UniquePtr` from the given one
		///
		/// @tparam U - The managed type of the given `UniquePtr`
		/// @tparam D - The deleter type of the given `UniquePtr`
		///
		/// @param ptr - The `UniquePtr<U, D>` to move into this one
		///
		/// # Requirements
		/// - `concepts::NoexceptAssignable<deleter_type,
		/// decltype(std::forward<D>(ptr.get_deleter()))>`: `deleter_type` must be assignable with
		/// the deleter type returned by `ptr`'s `get_deleter()` member function
		/// - `concepts::Convertible<typename UniquePtr<U, D>::pointer, pointer>`: The pointer type
		/// of `ptr` must be convertible to the pointer type of `this` in order to assign from
		/// it
		/// - `concepts::Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type>`:
		/// The deleter type of `ptr` must be convertible to the pointer type of `this` in order to
		/// assign from it
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		///
		/// @return this
		/// @ingroup UniquePtr
		template<typename U, typename D>
		constexpr auto operator=(UniquePtr<U, D>&& ptr) noexcept -> UniquePtr& requires
			concepts::NoexceptAssignable<deleter_type,
										 decltype(std::forward<D>(ptr.get_deleter()))> && concepts::
				Convertible<typename UniquePtr<U, D>::pointer, pointer> && concepts::
					Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type> &&(
						noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset(ptr.release());
			m_ptr.second() = std::forward<D>(ptr.get_deleter());
			return *this;
		}

		/// @brief Assigns this `UniquePtr` with `nullptr`
		///
		/// Calls the deleter on the managed pointer and replaces it with `nullptr`
		///
		/// # Requirements
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		///
		/// @return this
		/// @ingroup UniquePtr
		constexpr auto operator=(std::nullptr_t) noexcept -> UniquePtr& requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset();
			return *this;
		}

		/// @brief Assigns this `UniquePtr` with `ptr`
		///
		/// Calls the deleter on the managed pointer and replaces it with `ptr`
		///
		/// # Requirements
		/// - `noexcept(std::declval<deleter_type>()(std::declval<pointer>()))`: Deleting the
		/// managed pointer via the associated `deleter_type`'s call operator must be noexcept
		///
		/// @return this
		/// @ingroup UniquePtr
		constexpr auto operator=(pointer ptr) noexcept -> UniquePtr& requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset(ptr);
			return *this;
		}

		// clang-format off
	  private:
		template<typename T_, typename Deleter_>
		using CompressedPair = CompressedPair<T_, Deleter_>;
		template<typename T_>
		using DefaultInitTag = DefaultInitTag<T_>;

		CompressedPair<pointer, deleter_type> m_ptr;
		// clang-format on
	};

	template<typename T, typename Deleter>
	requires concepts::NotRValueReference<Deleter>
	class HYPERION_TRIVIAL_ABI UniquePtr<T[], Deleter> { // NOLINT
	  public:
		using element_type = T;
		using deleter_type = Deleter;
		using pointer = std::add_pointer_t<element_type>;
		using pointer_to_const = std::add_pointer_t<std::add_const_t<element_type>>;

		constexpr UniquePtr() noexcept requires concepts::NoexceptDefaultConstructible<
			deleter_type> && concepts::NotPointer<deleter_type>
			: m_ptr(pointer(), DefaultInitTag<deleter_type>()) {
		}

		constexpr UniquePtr(std::nullptr_t ptr) noexcept // NOLINT
			requires concepts::NoexceptDefaultConstructible<deleter_type> && concepts::NotPointer<
				deleter_type> : m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		explicit constexpr UniquePtr(pointer ptr) noexcept requires concepts::
			NoexceptDefaultConstructible<deleter_type> && concepts::NotPointer<deleter_type>
			: m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NoexceptDefaultConstructible<
			deleter_type> && concepts::NotPointer<deleter_type> && concepts::NotSame<U, pointer>
		explicit constexpr UniquePtr(U ptr) noexcept : m_ptr(ptr, DefaultInitTag<deleter_type>()) {
		}
		constexpr UniquePtr(pointer ptr, const Deleter& deleter) noexcept requires
			concepts::NoexceptCopyConstructible<Deleter> : m_ptr(ptr, deleter) {
		}
		constexpr UniquePtr(pointer ptr, Deleter&& deleter) noexcept requires
			concepts::NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
			: m_ptr(ptr, std::move(deleter)) {
		}
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NotSame<U, pointer> && concepts::
			NoexceptCopyConstructible<Deleter>
		constexpr UniquePtr(U ptr, const Deleter& deleter) noexcept : m_ptr(ptr, deleter) {
		}
		template<typename U>
		requires concepts::Convertible<U, pointer> && concepts::NotSame<U, pointer> && concepts::
			NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
		constexpr UniquePtr(U ptr, Deleter&& deleter) noexcept : m_ptr(ptr, std::move(deleter)) {
		}
		constexpr UniquePtr(std::nullptr_t ptr, const Deleter& deleter) noexcept requires
			concepts::NoexceptCopyConstructible<Deleter> : m_ptr(ptr, deleter) {
		}
		constexpr UniquePtr(std::nullptr_t ptr, Deleter&& deleter) noexcept requires
			concepts::NoexceptMoveConstructible<Deleter> && concepts::NotRValueReference<Deleter>
			: m_ptr(ptr, std::move(deleter)) {
		}
		// clang-format off
		template<typename U, typename D>
		requires std::is_array_v<U>
		explicit constexpr UniquePtr(UniquePtr<U, D>&& ptr)
			noexcept(concepts::NoexceptConstructibleFrom<deleter_type,
														 decltype((std::forward<D>(ptr.get_deleter())))>)
			requires concepts::Convertible<typename UniquePtr<U, D>::element_type(*)[], // NOLINT
										   element_type (*)[] > 						// NOLINT
					 && concepts::Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type>
			: m_ptr(ptr.release(), std::forward<D>(ptr.get_deleter())) {
		}
		// clang-format on
		UniquePtr(const UniquePtr&) = delete;
		constexpr UniquePtr(
			UniquePtr&& ptr) noexcept requires concepts::NoexceptMoveConstructible<deleter_type>
			: m_ptr(ptr.release(), std::forward<deleter_type>(ptr.get_deleter())) {
		}
		constexpr ~UniquePtr() noexcept
			requires(noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset();
		}

		[[nodiscard]] inline constexpr auto release() noexcept -> pointer {
			auto* ptr = m_ptr.first();
			m_ptr.first() = pointer();
			return ptr;
		}

		inline constexpr auto reset(pointer ptr = pointer()) noexcept -> void requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			gsl::owner<pointer> tmp = m_ptr.first(); // NOLINT
			m_ptr.first() = ptr;
			m_ptr.second()(tmp);
		}

		template<typename U>
		requires concepts::Convertible<U, pointer>
		inline constexpr auto reset(U ptr) noexcept -> void requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			gsl::owner<pointer> tmp = m_ptr.first(); // NOLINT
			m_ptr.first() = ptr;
			m_ptr.second()(tmp);
		}

		inline constexpr auto swap(UniquePtr& ptr) noexcept
			-> void requires concepts::NoexceptSwappable<CompressedPair<pointer, deleter_type>> {
			m_ptr.swap(ptr.m_ptr);
		}

		[[nodiscard]] inline constexpr auto get() noexcept -> pointer {
			return m_ptr.first();
		}

		[[nodiscard]] inline constexpr auto get() const noexcept -> pointer_to_const {
			return m_ptr.first();
		}

		[[nodiscard]] inline constexpr auto get_deleter() const noexcept -> const deleter_type& {
			return m_ptr.second();
		}

		[[nodiscard]] inline constexpr auto get_deleter() noexcept -> deleter_type& {
			return m_ptr.second();
		}

		explicit constexpr operator bool() const noexcept {
			return m_ptr.first() != nullptr;
		}

		inline constexpr auto operator==(const UniquePtr& ptr) const noexcept -> bool {
			return m_ptr.first() == ptr.first();
		}

		inline constexpr auto operator==(std::nullptr_t) const noexcept -> bool {
			return m_ptr.first() == nullptr;
		}

		// constexpr inline auto operator[](concepts::Integral auto i) const
		//	noexcept -> std::add_const_t<std::add_lvalue_reference_t<element_type>> {
		//	return m_ptr.first()[i];
		// }

		inline constexpr auto operator[](concepts::Integral auto index) const noexcept
			-> std::add_lvalue_reference_t<element_type> {
			return m_ptr.first()[index];
		}

		auto operator=(const UniquePtr&) -> UniquePtr& = delete;
		constexpr auto operator=(UniquePtr&& ptr) noexcept
			-> UniquePtr& requires concepts::NoexceptMoveAssignable<deleter_type> &&(
				noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			if(this == &ptr) {
				return *this;
			}

			reset(ptr.release());
			m_ptr.second() = std::forward<deleter_type>(ptr.get_deleter());

			return *this;
		}
		// clang-format off

		template<typename U, typename D>
		requires std::is_array_v<U>
		constexpr auto operator=(UniquePtr<U, D>&& ptr) noexcept -> UniquePtr&
			requires concepts::NoexceptAssignable<deleter_type,
												  decltype(std::forward<D>(ptr.get_deleter()))>
					 && concepts::Convertible<typename UniquePtr<U, D>::element_type(*)	[], // NOLINT
										   	  element_type(*)[]> 							// NOLINT
					 && concepts::Convertible<typename UniquePtr<U, D>::deleter_type, deleter_type>
					 && (noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset(ptr.release());
			m_ptr.second() = std::forward<D>(ptr.get_deleter());
			return *this;
		}
		// clang-format on

		constexpr auto operator=(std::nullptr_t) noexcept -> UniquePtr& requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset();
			return *this;
		}

		constexpr auto operator=(pointer ptr) noexcept -> UniquePtr& requires(
			noexcept(std::declval<deleter_type>()(std::declval<pointer>()))) // NOLINT
		{
			reset(ptr);
			return *this;
		}

		// clang-format off
	  private:
		template<typename T_, typename Deleter_>
		using CompressedPair = CompressedPair<T_, Deleter_>;
		template<typename T_>
		using DefaultInitTag = DefaultInitTag<T_>;

		CompressedPair<pointer, deleter_type> m_ptr;
		// clang-format on
	};
	template<typename T, typename Deleter>
	UniquePtr(T*, Deleter) -> UniquePtr<T, Deleter>;

	/// @brief Swaps the managed pointer and deleter of the two `UniquePtr`s
	///
	/// @tparam T - The type held by the `UniquePtr`s
	/// @tparam Deleter - The deleter type used by the `UniquePtr`s
	///
	/// @param first - The `UniquePtr` to swap to
	/// @param second - The `UniquePtr` to swap from
	/// @ingroup UniquePtr
	template<typename T, typename Deleter>
	requires concepts::Swappable<Deleter>
	inline constexpr auto
	swap(UniquePtr<T, Deleter>& first, UniquePtr<T, Deleter>& second) noexcept -> void {
		first.swap(second);
	}
	// clang-format off

	IGNORE_UNUSED_TEMPLATES_START

	/// @brief Constructs a `UniquePtr<T, Deleter>` from the given arguments
	///
	/// Constructs a `UniquePtr<T, Deleter>` from the given arguments, passing them directly
	/// to `T`'s constructor
	///
	/// @tparam T - The type to hold in the `UniquePtr`
	/// @tparam Deleter - The deleter type to free the resources associated with `T`
	/// @tparam Args - The types of the arguments to pass to `T`'s constructor
	///
	/// @param args - The arguments to pass to `T`'s constructor
	///
	/// # Requirements
	/// - `concepts::ConstructibleFrom<T, Args...>`: `T` must be constructible from the given
	/// arguments in order to make a `UniquePtr` from them
	/// - `!std::is_array_v<T>`: Cannot make a `UniquePtr` holding an array type from a set of
	/// constructor arguments. Use the overload for arrays to make a `UniquePtr<T[]>`.
	///
	/// @return a `UniquePtr<T, Deleter>` with the `T` constructed from the given arguments
	/// @ingroup UniquePtr
	template<typename T, typename Deleter = DefaultDeleter<T>, typename... Args>
	requires concepts::ConstructibleFrom<T, Args...> && (!std::is_unbounded_array_v<T>)
	inline static constexpr auto make_unique(Args&&... args)
		noexcept(concepts::NoexceptConstructibleFrom<T, Args...>)
		-> UniquePtr<T, Deleter>
	{
		// NOLINTNEXTLINE(modernize-use-auto,hicpp-use-auto,bugprone-unhandled-exception-at-new)
		gsl::owner<T*> ptr = new T(std::forward<Args>(args)...);
		// NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
		return UniquePtr<T, Deleter>(ptr);
	}
	// clang-format on

	/// @brief Constructs a `UniquePtr` holding array type `T` of size `N`
	///
	/// Constructs a `UniquePtr<T, Deleter>`, where `T` is an unbounded array type with initial
	/// size of `N`.
	///
	/// @tparam T - The array type to hold in the `UniquePtr`
	/// @tparam Deleter - The deleter type to free the resources associated with `T`
	/// @tparam ElementType - The element type of the array
	///
	/// @param N - The initial size of the array managed by the `UniquePtr`
	///
	/// # Requirements
	/// - `std::is_unbounded_array_v<T>`: `T` must be an unbounded array type to make a `UniquePtr`
	/// managing an array.
	///
	/// @return a `UniquePtr<T, Deleter>` managing an array of `N` `ElementType`s
	/// @note The array managed by the returned `UniquePtr` will be uninitialized, and thus each
	/// element must be appropriately initialized before use. Otherwise, the result is undefined
	/// behavior.
	/// @ingroup UniquePtr
	template<typename T,
			 typename Deleter = DefaultDeleter<T>,
			 typename ElementType = std::remove_extent_t<T>> // NOLINT
	requires std::is_unbounded_array_v<T>
	inline static constexpr auto make_unique(usize N) noexcept -> UniquePtr<T, Deleter> {
		// NOLINTNEXTLINE(modernize-use-auto,hicpp-use-auto,bugprone-unhandled-exception-at-new)
		gsl::owner<ElementType*> ptr = new ElementType[N];
		// NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
		return UniquePtr<T, Deleter>(ptr);
	}
	IGNORE_UNUSED_TEMPLATES_STOP

	IGNORE_PADDING_START
	/// @brief An allocator-aware deleter type for Hyperion smart pointers
	///
	/// `AllocatorAwareDeleter` is useful for using smart pointers in situations that require custom
	/// allocation strategies. It allows use of smart pointers in an allocator-aware fashion.
	///
	/// @tparam T - The type to manage deletion for
	/// @tparam Allocator - The allocator type allocation was performed with
	/// @tparam ElementType - The element type of the managed deletion, in the case that `T` is an
	/// array
	///
	/// # Requirements
	/// - `concepts::Allocatable<ElementType, Allocator>`: `AllocatorAwareDeleter` can't be
	/// instantiated for an `ElementType`-`Allocator` pair where the element type `ElementType`
	/// isn't allocatable by the associated allocator type, `Allocator`
	/// @note In the case that `T` is an array type, `AllocatorAwareDeleter` requires that the
	/// managed array is fully-initialized and the same size as when the `AllocatorAwareDeleter`
	/// was constructed with it (by a call to one of the factory functions, `allocate_unique`, or
	/// `allocate_shared`) when its call operator is used to free the resources associated with the
	/// array.
	/// @note If `T` is an array type and the array is resized, the deleter associated with it in
	/// the owning smart pointer must be set to a new one to match the changed state.
	/// @note If `T` is an array type and the members of the array are not all initialized when the
	/// call operator of this deleter is used to free the resources associated with it, the result
	/// is undefined behavior.
	/// @ingroup memory
	template<typename T,
			 typename Allocator = std::allocator<T>,
			 typename ElementType = std::remove_cv_t<std::remove_all_extents_t<T>>>
	requires concepts::Allocatable<ElementType, Allocator>
	class AllocatorAwareDeleter {
	  public:
		/// @brief The rebound allocator type for this deleter
		/// @ingroup memory
		using Alloc = typename std::allocator_traits<Allocator>::template rebind_alloc<ElementType>;
		/// @brief The `std::allocator_traits` type for this deleter
		/// @ingroup memory
		using Traits = std::allocator_traits<Alloc>;
		/// @brief The pointer type for the associated allocator traits for this deleter
		/// @ingroup memory
		using pointer = typename Traits::pointer;

		/// @brief Default-Constructs an `AllocatorAwareDeleter`
		///
		/// # Requirements
		/// - `concepts::NoexceptDefaultConstructible<Alloc>`: The associated allocator type must be
		/// noexcept default constructible in order to construct an `AllocatorAwareDeleter` with one
		/// @ingroup memory
		constexpr AllocatorAwareDeleter() noexcept requires
			concepts::NoexceptDefaultConstructible<Alloc> : m_allocator() {
		}
		/// @brief Constructs an `AllocatorAwareDeleter` from the given allocator
		///
		/// @param alloc - The allocator to construct the allocator associated with this from
		///
		/// # Requirements
		/// - `concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)>`: The associated
		/// allocator type must be noexcept constructible from the given allocator in order to
		/// construct an `AllocatorAwareDeleter` from it
		/// @ingroup memory
		explicit constexpr AllocatorAwareDeleter(const Allocator& alloc) noexcept requires
			concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)> : m_allocator(alloc) {
		}
		/// @brief Copy-Constructs an `AllocatorAwareDeleter` from the given one
		///
		/// # Requirements
		/// - `concepts::NoexceptCopyConstructible<Alloc>`: The associated allocator type must be
		/// noexcept copy constructible in order to copy construct an `AllocatorAwareDeleter`
		/// @ingroup memory
		constexpr AllocatorAwareDeleter(const AllocatorAwareDeleter&) noexcept requires
			concepts::NoexceptCopyConstructible<Alloc>
		= default;
		/// @brief Move-Constructs an `AllocatorAwareDeleter` from the given one
		///
		/// # Requirements
		/// - `concepts::NoexceptMoveConstructible<Alloc>`: The associated allocator type must be
		/// noexcept move constructible in order to copy construct an `AllocatorAwareDeleter`
		/// @ingroup memory
		constexpr AllocatorAwareDeleter(
			AllocatorAwareDeleter&&) noexcept requires concepts::NoexceptMoveConstructible<Alloc>
		= default;

		/// @brief Destructs this `AllocatorAwareDeleter`
		/// @ingroup memory
		constexpr ~AllocatorAwareDeleter() noexcept = default;

		/// @brief Frees the resources associated with the pointed-to object
		///
		/// @param p - Pointer to the object to free associated resources of
		/// @ingroup memory
		inline constexpr auto operator()(pointer p) const noexcept {
			Alloc allocator = m_allocator;

			Traits::destroy(allocator, std::addressof(*p));

			Traits::deallocate(allocator, p, 1);
		}

		/// @brief Copy-Assigns this `AllocatorAwareDeleter` from the given one
		///
		/// # Requirements
		/// - `concepts::NoexceptCopyAssignable<Alloc>`: The associated allocator type must be
		/// noexcept copy assignable in order to copy assign this `AllocatorAwareDeleter`
		/// @ingroup memory
		constexpr auto operator=(const AllocatorAwareDeleter&) noexcept
			-> AllocatorAwareDeleter& requires concepts::NoexceptCopyAssignable<Alloc>
		= default;
		/// @brief Move-Assigns this `AllocatorAwareDeleter` from the given one
		///
		/// # Requirements
		/// - `concepts::NoexceptMoveAssignable<Alloc>`: The associated allocator type must be
		/// noexcept move assignable in order to move assign this `AllocatorAwareDeleter`
		/// @ingroup memory
		constexpr auto operator=(AllocatorAwareDeleter&&) noexcept
			-> AllocatorAwareDeleter& requires concepts::NoexceptMoveAssignable<Alloc>
		= default;

	  private:
		Alloc m_allocator;
	};

	template<typename T, typename Allocator, typename ElementType>
	requires concepts::Allocatable<ElementType, Allocator>
	class AllocatorAwareDeleter<T[], Allocator, ElementType> { // NOLINT
	  public:
		using Alloc = typename std::allocator_traits<Allocator>::template rebind_alloc<ElementType>;
		using Traits = std::allocator_traits<Alloc>;
		using pointer = typename Traits::pointer;

		constexpr AllocatorAwareDeleter() noexcept requires
			concepts::NoexceptDefaultConstructible<Alloc> : m_allocator() {
		}
		explicit constexpr AllocatorAwareDeleter(const Allocator& alloc,
												 usize num_elements) noexcept requires
			concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)>
			: m_allocator(alloc),
			  m_num_elements(num_elements) {
		}
		constexpr AllocatorAwareDeleter(const AllocatorAwareDeleter&) noexcept requires
			concepts::NoexceptCopyConstructible<Alloc>
		= default;
		constexpr AllocatorAwareDeleter(
			AllocatorAwareDeleter&&) noexcept requires concepts::NoexceptMoveConstructible<Alloc>
		= default;
		constexpr ~AllocatorAwareDeleter() noexcept = default;

		inline constexpr auto operator()(pointer p) const noexcept {
			Alloc allocator = m_allocator;

			auto i = m_num_elements >= 1_usize ? m_num_elements - 1_usize : 0_usize;
			// we need to use do-while to ensure we destroy the 0th element
			// in a 1 element array
			do {
				Traits::destroy(allocator, std::addressof(*p) + i); // NOLINT
				--i;
			} while(i != 0_usize);

			Traits::deallocate(allocator, p, m_num_elements);
		}

		constexpr auto operator=(const AllocatorAwareDeleter&) noexcept
			-> AllocatorAwareDeleter& requires concepts::NoexceptCopyAssignable<Alloc>
		= default;
		constexpr auto operator=(AllocatorAwareDeleter&&) noexcept
			-> AllocatorAwareDeleter& requires concepts::NoexceptMoveAssignable<Alloc>
		= default;

	  private:
		Alloc m_allocator;
		usize m_num_elements = 1_usize;
	};
	IGNORE_PADDING_STOP

	/// @brief Constructs an allocator-aware `UniquePtr`
	///
	/// Constructs a `UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>`. The managed `T` will be
	/// constructed from the given arguments, `args`, passing them directly to `T`'s constructor.
	/// Allocation of the managed `T` will be performed by the given allocator.
	///
	/// @tparam T - The type to manage in the `UniquePtr`
	/// @tparam Allocator - The type of the allocator to use for allocation and deletion of the `T`
	/// @tparam ElementType - The element type of the array, in the case that `T` is an unbounded
	/// array type
	/// @tparam Alloc - The allocator type `Allocator` rebound to the element type, `ElementType`
	///
	/// @param alloc - The allocator to perform allocation and deallocation with
	/// @param args - The arguments to pass to `T`'s constructor
	///
	/// # Requirements
	/// - `concepts::NoexceptConstructibleFrom<ElementType, Args...>`: `ElementType` must be
	/// noexcept constructible from `args` in order to make a `UniquePtr` from them
	/// - `concepts::Allocatable<ElementType, Alloc>`: `AllocatorAwareDeleter` can't be
	/// instantiated for an `ElementType`-`Alloc` pair where the element type `ElementType`
	/// isn't allocatable by the associated allocator type, `Alloc`
	/// - `!std::is_array_v<T>`: Cannot make a `UniquePtr` holding an array type from a set of
	/// constructor arguments. Use the overload for arrays to make a `UniquePtr<T[]>`.
	/// - `concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)>`: The associated allocator
	/// type, `Alloc`, must be noexcept constructible from the given allocator in order to create
	/// a `UniquePtr` using it
	/// @ingroup UniquePtr
	template<typename T,
			 typename Allocator = std::allocator<T>,
			 typename ElementType = std::remove_cv_t<std::remove_all_extents_t<T>>,
			 typename Alloc =
				 typename std::allocator_traits<Allocator>::template rebind_alloc<ElementType>,
			 typename... Args>
	requires concepts::NoexceptConstructibleFrom<ElementType, Args...> && concepts::
		Allocatable<ElementType, Alloc> &&(!std::is_unbounded_array_v<T>)
			[[nodiscard]] inline constexpr auto allocate_unique(const Allocator& alloc,
																Args&&... args) noexcept
		-> UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>
	requires concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)> {
		using Traits = std::allocator_traits<Alloc>;
		using Deleter = AllocatorAwareDeleter<T, Alloc>;

		Alloc allocator(alloc);
		auto* p = Traits::allocate(allocator, 1);
		Traits::construct(allocator, std::addressof(*p), std::forward<Args>(args)...);

		// NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
		return UniquePtr<T, Deleter>(p, Deleter(allocator));
	}

	/// @brief Constructs an allocator-aware `UniquePtr`
	///
	/// Constructs a `UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>`, where `T` is an unbounded
	/// array type. The managed array will have its elements default constructed in place.
	/// Allocation of the managed array will be performed by the given allocator.
	///
	/// @tparam T - The type to manage in the `UniquePtr`
	/// @tparam Allocator - The type of the allocator to use for allocation and deletion of the `T`
	/// @tparam ElementType - The element type of the array, in the case that `T` is an unbounded
	/// array type
	/// @tparam Alloc - The allocator type `Allocator` rebound to the element type, `ElementType`
	///
	/// @param alloc - The allocator to perform allocation and deallocation with
	/// @param N - The number of elements to allocate in the array
	///
	/// # Requirements
	/// - `concepts::NoexceptDefaultConstructible<ElementType>`: `ElementType` must be
	/// noexcept default-constructible in order to make a `UniquePtr` managing an array of them
	/// - `concepts::Allocatable<ElementType, Alloc>`: `AllocatorAwareDeleter` can't be
	/// instantiated for an `ElementType`-`Alloc` pair where the element type `ElementType`
	/// isn't allocatable by the associated allocator type, `Alloc`
	/// - `std::is_unbounded_array_v<T>`: `T` must be an unbounded array type to make a `UniquePtr`
	/// managing an array.
	/// - `concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)>`: The associated allocator
	/// type, `Alloc`, must be noexcept constructible from the given allocator in order to create
	/// a `UniquePtr` using it
	/// @ingroup UniquePtr
	template<typename T,
			 typename Allocator = std::allocator<T>,
			 typename ElementType = std::remove_cv_t<std::remove_all_extents_t<T>>,
			 typename Alloc =
				 typename std::allocator_traits<Allocator>::template rebind_alloc<ElementType>>
	requires concepts::NoexceptDefaultConstructible<
		ElementType> && concepts::Allocatable<ElementType, Alloc> && std::is_unbounded_array_v<T>
	[[nodiscard]] inline constexpr auto allocate_unique(const Allocator& alloc, usize N) noexcept
		-> UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>
	requires concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)> {
		using Traits = std::allocator_traits<Alloc>;
		using Deleter = AllocatorAwareDeleter<ElementType[], Alloc>; // NOLINT (c arrays)

		Alloc allocator(alloc);
		auto* p = Traits::allocate(allocator, N);
		// we need to use do-while to make sure we construct the 0th element in a one element array
		auto i = 0_usize;
		do {
			Traits::construct(allocator, std::addressof(*p) + i); // NOLINT
			++i;
		} while(i < N);

		// NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
		return UniquePtr<ElementType[], Deleter>(p, Deleter(allocator, N)); // NOLINT ( c arrays)
	}

	/// @brief Constructs an allocator-aware `UniquePtr`
	///
	/// Constructs a `UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>`, where `T` is an unbounded
	/// array type. The managed array will have its elements default constructed in place.
	/// Allocation of the managed array will be performed by the given allocator.
	///
	/// @tparam T - The type to manage in the `UniquePtr`
	/// @tparam Allocator - The type of the allocator to use for allocation and deletion of the `T`
	/// @tparam ElementType - The element type of the array, in the case that `T` is an unbounded
	/// array type
	/// @tparam Alloc - The allocator type `Allocator` rebound to the element type, `ElementType`
	/// @tparam Args - The types of the arguments to use to default-construct the elements of the
	/// array
	///
	/// @param alloc - The allocator to perform allocation and deallocation with
	/// @param N - The number of elements to allocate in the array
	///
	/// # Requirements
	/// - `concepts::NoexceptDefaultConstructible<ElementType>`: `ElementType` must be
	/// noexcept default-constructible in order to make a `UniquePtr` managing an array of them
	/// - `concepts::Allocatable<ElementType, Alloc>`: `AllocatorAwareDeleter` can't be
	/// instantiated for an `ElementType`-`Alloc` pair where the element type `ElementType`
	/// isn't allocatable by the associated allocator type, `Alloc`
	/// - `std::is_unbounded_array_v<T>`: `T` must be an unbounded array type to make a `UniquePtr`
	/// managing an array.
	/// - `concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)>`: The associated allocator
	/// type, `Alloc`, must be noexcept constructible from the given allocator in order to create
	/// a `UniquePtr` using it
	/// @ingroup UniquePtr
	template<typename T,
			 typename Allocator = std::allocator<T>,
			 typename ElementType = std::remove_cv_t<std::remove_all_extents_t<T>>,
			 typename Alloc =
				 typename std::allocator_traits<Allocator>::template rebind_alloc<ElementType>,
			 typename... Args>
	requires concepts::NoexceptConstructibleFrom<ElementType, Args...> && mpl::for_all_types_v<
		std::is_nothrow_copy_constructible,
		std::true_type,
		mpl::list<Args...>> && concepts::Allocatable<ElementType,
													 Alloc> && std::is_unbounded_array_v<T>
	[[nodiscard]] inline constexpr auto
	allocate_unique(const Allocator& alloc, usize N, Args&&... args) noexcept
		-> UniquePtr<T, AllocatorAwareDeleter<T, Alloc>>
	requires concepts::NoexceptConstructibleFrom<Alloc, decltype(alloc)> {
		using Traits = std::allocator_traits<Alloc>;
		using Deleter = AllocatorAwareDeleter<ElementType[], Alloc>; // NOLINT (c arrays)

		Alloc allocator(alloc);
		auto* p = Traits::allocate(allocator, N);
		auto tuple = std::make_tuple(std::forward<Args>(args)...);
		// we need to use do-while to make sure we construct the 0th element in a one element array
		auto i = 0_usize;
		do {
			auto t = std::tuple(tuple);
			Traits::construct(allocator,
							  std::addressof(*p) + i,
							  std::make_from_tuple<ElementType>(std::move(t))); // NOLINT
			++i;
		} while(i < N);

		return UniquePtr<ElementType[], Deleter>(p, Deleter(allocator, N)); // NOLINT (c arrays)
	}

#if HYPERION_DEFINE_TESTS

	// NOLINTNEXTLINE(modernize-use-trailing-return-type)
	TEST_SUITE("UniquePtr") {
		TEST_CASE("Constructor") {
			auto ptr1 = UniquePtr<i32>();
			// NOLINTNEXTLINE
			auto ptr2 = UniquePtr<i32>(new i32(3_i32));
			auto ptr3 = hyperion::make_unique<i32>(2_i32);

			CHECK_EQ(ptr1, nullptr);
			CHECK_NE(ptr2, nullptr);
			CHECK_NE(ptr3, nullptr);
			CHECK_EQ(*ptr2, 3_i32);
			CHECK_EQ(*ptr3, 2_i32);

			SUBCASE("move") {
				auto ptr4 = std::move(ptr3);
				// NOLINTNEXTLINE(bugprone-use-after-move,hicpp-invalid-access-moved)
				CHECK_EQ(ptr3, nullptr);
				CHECK_NE(ptr4, nullptr);
				CHECK_EQ(*ptr4, 2_i32);
			}

			SUBCASE("accessors_and_modifiers") {
				CHECK(ptr3);
				CHECK(static_cast<bool>(ptr3));
				CHECK_NE(ptr3.get(), nullptr);
				CHECK_EQ(*(ptr3.get()), 2_i32);

				auto* ptr4 = ptr3.release();
				CHECK_EQ(ptr3, nullptr);
				CHECK_NE(ptr4, nullptr);
				CHECK_EQ(*ptr4, 2_i32);

				*ptr4 = 4_i32;
				ptr3.reset(ptr4);
				CHECK_NE(ptr3, nullptr);
				CHECK(ptr3);
				CHECK_EQ(*ptr3, 4_i32);
				CHECK_EQ(*(ptr3.get()), 4_i32);
			}
		}
	}
#endif // HYPERION_DEFINE_TESTS
} // namespace hyperion