text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include "parmt_postProcess.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

static double *diff(const int n, const double *__restrict__ x, int *ierr);
static int computePadding8i(const int n);

int marginal_getOptimum(const int nloc, const int nm, const int nb,
                        const int ng, const int nk, const int ns, const int nt,
                        const double *__restrict__ phi,
                        int *jloc, int *jm, int *jb, int *jg,
                        int *jk, int *js, int *jt)
{
    int ib, ig, ik, iloc, imax, im, imt, indx, is, it, nmt;
    enum isclError_enum ierr;
    nmt = nm*nb*ng*nk*ns*nt;
    imax = array_argmax64f(nmt*nloc, phi, &ierr);
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk; ik++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            for (it=0; it<nt; it++)
                            {
                                imt = im*nb*ng*nk*ns*nt
                                    + ib*ng*nk*ns*nt
                                    + ig*nk*ns*nt
                                    + ik*ns*nt
                                    + is*nt
                                    + it;
                                indx = iloc*nmt + imt;
                                if (indx == imax)
                                {
                                    *jloc = iloc;
                                    *jm = im;
                                    *jb = ib;
                                    *jg = ig;
                                    *jk = ik;
                                    *js = is;
                                    *jt = it;
                                }
                            }
                        }
                    }
                }
            }
         }
     }
     return 0;
}
/*!
 * @brief Computes the marginal PDF at lune points using Eqn 48 of
 *        Tape and Tape 2015.
 *
 * @param[in] nloc       number of locations in grid-search.
 * @param[in] nm         number of moment magnitudes in grid search
 * @param[in] nb         number of colatitudes
 * @param[in] betas      colatitudes (radians) of lune points [nb]
 * @param[in] ng         number of longitudes
 * @param[in] gammas     longitudes (radians) of lune points [ng]
 * @param[in] nk         number of strike angles
 * @param[in] kappas     fault strike angles (radians) [nk]
 * @param[in] ns         number of slip angles
 * @param[in] sigmas     fault slip angles (radians) [ns]
 * @param[in] nh         number of dip angles
 * @param[in] h          related to dip angles (radians) [nh]
 * @param[in] phi        objective function computed throughout grid-search
 *                       space to marginalize [nloc*nm*nb*ng*nk*ns*nh]
 * @param[out] luneMPDF  marginalized lune [nb x ng].  the ib, ig'th
 *                       colatitude and longitude is given by ib*ng + ig.
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 * @bug There is inadequate handling of moment tensor scaling and locs
  *     is a misnomer.  Technically this can only loop in 1D and I would
 *      require the cell volumes for position to loop on locations.
 */
int marginal_computeLuneMPDF(const int nloc,
                             const int nm,
                             const int nb, const double *__restrict__ betas,
                             const int ng, const double *__restrict__ gammas,
                             const int nk, const double *__restrict__ kappas,
                             const int ns, const double *__restrict__ sigmas,
                             const int nh, const double *__restrict__ h,
                             const double *__restrict__ phi,
                             double *__restrict__ luneMPDF)
{
    const char *fcnm = "marginal_computeLuneMPDF\0";
    int ib, ierr, ig, igIb, ih, ik, iloc, im, imt, indx, is, nmt;
    double *cos3g, *dh, *dk, *ds, *twoSinB3;
    double sinb, sinb3, dkdsdh, dV;
    nmt = nm*nb*ng*nk*ns*nh;
    memset(luneMPDF, 0, (size_t) (ng*nb)*sizeof(double));
    if (nk < 2 || ns < 2 || nh < 2)
    {
        if (nk < 2){printf("%s: Error nk must be > 1\n", fcnm);}
        if (ns < 2){printf("%s: Error ns must be > 1\n", fcnm);}
        if (nh < 2){printf("%s: Error nh must be > 1\n", fcnm);}
    }
    // Differentiate weights for Riemann quadrature
    dk = diff(nk, kappas, &ierr);
    ds = diff(ns, sigmas, &ierr);
    dh = diff(nh, h, &ierr);
    // theta goes from 0 to pi/2 as h goes from 1 to 0 hence it's
    // orientation is reversed from our ini file (recall that it asks
    // for theta and not h) and hence -dh is not conducive to quadrature.
    cblas_dscal(nh-1, -1.0, dh, 1);
    // Precompute the trigonemetric terms
    twoSinB3 = memory_calloc64f(nb);
    for (ib=0; ib<nb; ib++)
    {
        twoSinB3[ib] = 2.0*pow(sin(betas[ib]), 3);
    }
    cos3g = memory_calloc64f(ng);
    for (ig=0; ig<ng; ig++)
    {
        cos3g[ig] = cos(3.0*gammas[ig]);
    }
    // Marginalize the lune
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk-1; ik++)
                    {
                        for (is=0; is<ns-1; is++)
                        {
                            for (ih=0; ih<nh-1; ih++)
                            {
                                imt = im*nb*ng*nk*ns*nh
                                    + ib*ng*nk*ns*nh
                                    + ig*nk*ns*nh
                                    + ik*ns*nh
                                    + is*nh
                                    + ih;
                                igIb = ib*ng + ig;
                                indx = iloc*nmt + imt;
                                //sinb = sin(betas[ib]);
                                //sinb3 = (sinb*sinb)*sinb;
                                //twoSinb3 = 2.0*sinb3;
                                //cos3g = cos(3.0*gammas[ig]);
                                dkdsdh = (dk[ik]*ds[is])*dh[ih];
                                dV = twoSinB3[ib]*cos3g[ig]*dkdsdh;
                                luneMPDF[igIb] = luneMPDF[igIb] + phi[indx]*dV;
                            }
                        }
                    }
                }
            }
        }
    }
    memory_free64f(&dh);
    memory_free64f(&dk);
    memory_free64f(&ds);
    memory_free64f(&twoSinB3);
    memory_free64f(&cos3g);
    return 0;
}

int marginal_computeLuneUVMPDF(const int nloc,
                               const int nm,
                               const int nu, const double *__restrict__ us,
                               const int nv, const double *__restrict__ vs,
                               const int nk, const double *__restrict__ kappas,
                               const int ns, const double *__restrict__ sigmas,
                               const int nh, const double *__restrict__ hs,
                               const double *__restrict__ phi,
                               double *__restrict__ luneMPDF)
{
    const char *fcnm = "marginal_computeLuneMPDF\0";
    int ierr, ih, ik, iloc, im, imt, indx, is, iu, iuIv, iv, nmt;
    double *dh, *dk, *ds, *du, *dv;
    double dudv, dkdsdh, dV;
    const double sqrt8 = sqrt(8.0);
    nmt = nm*nu*nv*nk*ns*nh;
    memset(luneMPDF, 0, (size_t) (nu*nv)*sizeof(double));
    if (nk < 2 || ns < 2 || nh < 2)
    {
        if (nk < 2){printf("%s: Error nk must be > 1\n", fcnm);}
        if (ns < 2){printf("%s: Error ns must be > 1\n", fcnm);}
        if (nh < 2){printf("%s: Error nh must be > 1\n", fcnm);}
    }
    // Differentiate weights for Riemann quadrature
    dk = diff(nk, kappas, &ierr);
    ds = diff(ns, sigmas, &ierr);
    dh = diff(nh, hs,     &ierr);
/*
    du = diff(nu, us,     &ierr);
    dv = diff(nv, vs,     &ierr);
for (int i=0; i<nu; i++)
{
  printf("du: %f: \n", du[i]);
}
for (int i=0; i<nv; i++)
{
 printf("dv: %f: \n", dv[i]);
}
*/
    // theta goes from 0 to pi/2 as h goes from 1 to 0 hence it's
    // orientation is reversed from our ini file (recall that it asks
    // for theta and not h) and hence -dh is not conducive to quadrature.
    cblas_dscal(nh-1, -1.0, dh, 1);
    // Marginalize the lune
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (iu=0; iu<nu; iu++)
            {
                for (iv=0; iv<nv; iv++)
                {
                    for (ik=0; ik<nk-1; ik++)
                    {
                        for (is=0; is<ns-1; is++)
                        {
                            for (ih=0; ih<nh-1; ih++)
                            {
                                imt = im*nu*nv*nk*ns*nh
                                    + iu*nv*nk*ns*nh
                                    + iv*nk*ns*nh
                                    + ik*ns*nh
                                    + is*nh
                                    + ih;
                                iuIv = iv*nu + iu;
                                indx = iloc*nmt + imt;
                                //sinb = sin(betas[ib]);
                                //sinb3 = (sinb*sinb)*sinb;
                                //twoSinb3 = 2.0*sinb3;
                                //cos3g = cos(3.0*gammas[ig]);
                                dudv = sqrt8; //du[iu]*dv[iv];
                                dkdsdh = (dk[ik]*ds[is])*dh[ih];
                                dV = dudv*dkdsdh;
                                luneMPDF[iuIv] = luneMPDF[iuIv] + phi[indx]*dV;
                            }
                        }
                    }
                }
            }
        }
    }
/*
    memory_free64f(&du);
    memory_free64f(&dv);
*/
    memory_free64f(&dk);
    memory_free64f(&ds);
    memory_free64f(&dh);
    return 0;
}
//============================================================================//
int marginal_computeDepthMPDF(const int nloc,
                              const int nm, const double *__restrict__ M0s,
                              const int nb, const double *__restrict__ betas,
                              const int ng, const double *__restrict__ gammas,
                              const int nk, const double *__restrict__ kappas,
                              const int ns, const double *__restrict__ sigmas,
                              const int nt, const double *__restrict__ thetas,
                              const double *__restrict__ phi,
                              double *__restrict__ depMagMPDF,
                              double *__restrict__ depMPDF)
{
    const char *fcnm = "marginal_computeDepthMPDF\0";
    double *cos3g, *db,  *dg, *dk, *dm, *ds, *dt, *twoSinB4, *sint,
           arg, dbdg, dkdsdt, dV5, geom, sum;
    int ib, ierr, ig, ik, iloc, im, imt, indx, is, it, jloc;
    ierr = 0;
    memset(depMPDF, 0, (size_t) nloc*sizeof(double));
    memset(depMagMPDF, 0, (size_t) (nloc*nm)*sizeof(double));
    // Cell spacings for for Riemann quadrature
    dm = memory_calloc64f(nm);
    db = memory_calloc64f(nb);
    dg = memory_calloc64f(ng);
    dk = memory_calloc64f(nk);
    ds = memory_calloc64f(ns);
    dt = memory_calloc64f(nt);
    ierr += postprocess_computeBetaCellSpacing(nb,  betas, db);
    ierr += postprocess_computeGammaCellSpacing(ng, gammas, dg);
    ierr += postprocess_computeKappaCellSpacing(nk, kappas, dk);
    ierr += postprocess_computeSigmaCellSpacing(ns, sigmas, ds);
    ierr += postprocess_computeThetaCellSpacing(nt, thetas, dt);
    ierr += postprocess_computeM0CellSpacing(nm, M0s, dm);
printf("%s: for now i'm setting dm = 1 and dl = 1\n", fcnm);
array_set64f_work(nm, 1.0, dm);
    if (ierr != 0)
    {
        printf("%s: Failed to compute cell spacing\n", fcnm);
        return -1;
    }
    if (array_min64f(nm, dm, &ierr) <= 0.0 ||
        array_min64f(nb, db, &ierr) <= 0.0 ||
        array_min64f(ng, dg, &ierr) <= 0.0 ||
        array_min64f(nk, dk, &ierr) <= 0.0 ||
        array_min64f(ns, ds, &ierr) <= 0.0 ||
        array_min64f(nt, dt, &ierr) <= 0.0)
    {
        printf("%s: Negative jacobian\n", fcnm);
        return -1;
    }
    // Compute the geometric factors
    twoSinB4 = memory_calloc64f(nb);
    for (ib=0; ib<nb; ib++)
    {
        twoSinB4[ib] = 2.0*pow(sin(betas[ib]), 4);
    }
    sint = memory_calloc64f(nt);
    for (it=0; it<nt; it++)
    {
        sint[it] = sin(thetas[it]);
    }
    cos3g = memory_calloc64f(ng);
    for (ig=0; ig<ng; ig++)
    {
        cos3g[ig] = cos(3.0*gammas[ig]);
    }
    if (array_min64f(nb, twoSinB4, &ierr) <= 0.0 ||
        array_min64f(nt, sint, &ierr) <= 0.0 ||
        array_min64f(ng, cos3g, &ierr) <= 0.0)
    {
        printf("%s: Warning negative jacobian from geometric factors\n", fcnm);
    }
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk; ik++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            for (it=0; it<nt; it++)
                            {
                                imt = iloc*nm*nb*ng*nk*ns*nt
                                    +      im*nb*ng*nk*ns*nt
                                    +         ib*ng*nk*ns*nt
                                    +            ig*nk*ns*nt
                                    +               ik*ns*nt
                                    +                  is*nt
                                    +                     it;
                                dkdsdt = (dk[ik]*ds[is])*dt[it];
                                dbdg = db[ib]*dg[ig];
                                geom = (twoSinB4[ib]*cos3g[ig])*sint[it];
                                dV5 = geom*(dbdg*dkdsdt);
                                //if (dV5 < 0.0){printf("error\n");}
                                arg = phi[imt]*dV5;
                                jloc = im*nloc + iloc;
                                depMagMPDF[jloc] = depMagMPDF[jloc] + arg;
                            }
                        }
                    }
                }
            }
        }
    }
    // Integrate out the depMagMPDF
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            jloc = im*nloc + iloc;
            depMPDF[iloc] = depMPDF[iloc] +  depMagMPDF[jloc]*dm[im];
        }
    }
    memory_free64f(&cos3g);
    memory_free64f(&sint);
    memory_free64f(&twoSinB4);
    memory_free64f(&db);
    memory_free64f(&dg);
    memory_free64f(&dk);
    memory_free64f(&ds);
    memory_free64f(&dt);
    memory_free64f(&dm);
    return ierr;
}
//============================================================================//
double marginal_computeNormalization(
    const int nloc, const double *__restrict__ deps,
    const int nm, const double *__restrict__ M0s,
    const int nb, const double *__restrict__ betas,
    const int ng, const double *__restrict__ gammas,
    const int nk, const double *__restrict__ kappas,
    const int ns, const double *__restrict__ sigmas,
    const int nt, const double *__restrict__ thetas,
    const double *__restrict__ phi,
    int *ierr)
{
    const char *fcnm = "marginal_computeNormalization\0";
    double *cos3g, *db,  *dg, *dk, *dl, *dm, *ds, *dt, *twoSinB4, *sint,
           dbdg, dkdsdt, dV2, dV5, dV7, geom, sum;
    int ib, ig, ik, iloc, im, imt, is, it, nmt;
    *ierr = 0;
    sum = 0.0;
    nmt = nm*nb*ng*nk*ns*nt;
    // Cell spacings for for Riemann quadrature
    dl = memory_calloc64f(nloc);
    dm = memory_calloc64f(nm);
    db = memory_calloc64f(nb);
    dg = memory_calloc64f(ng);
    dk = memory_calloc64f(nk);
    ds = memory_calloc64f(ns);
    dt = memory_calloc64f(nt);
    *ierr += postprocess_computeBetaCellSpacing(nb,  betas, db);
    *ierr += postprocess_computeGammaCellSpacing(ng, gammas, dg);
    *ierr += postprocess_computeKappaCellSpacing(nk, kappas, dk);
    *ierr += postprocess_computeSigmaCellSpacing(ns, sigmas, ds);
    *ierr += postprocess_computeThetaCellSpacing(nt, thetas, dt);
    *ierr += postprocess_computeM0CellSpacing(nm, M0s, dm);
    *ierr += postprocess_computeM0CellSpacing(nloc, deps, dl);
printf("%s: for now i'm setting dm = 1 and dl = 1\n", fcnm);
array_set64f_work(nm, 1.0, dm);
array_set64f_work(nloc, 1.0, dl);
    if (*ierr != 0)
    {
        printf("%s; Error computing normalization\n", fcnm);
        return sum;
    }
    if (array_min64f(nloc, dl, ierr) <= 0.0 ||
        array_min64f(nm, dm, ierr) <= 0.0 ||
        array_min64f(nb, db, ierr) <= 0.0 ||
        array_min64f(ng, dg, ierr) <= 0.0 ||
        array_min64f(nk, dk, ierr) <= 0.0 ||
        array_min64f(ns, ds, ierr) <= 0.0 ||
        array_min64f(nt, dt, ierr) <= 0.0)
    {
        printf("%s: Negative jacobian\n", fcnm);
        *ierr = 1;
        return sum;
    }
    // Compute the geometric factors
    twoSinB4 = memory_calloc64f(nb);
    for (ib=0; ib<nb; ib++)
    {
        twoSinB4[ib] = 2.0*pow(sin(betas[ib]), 4);
    }
    sint = memory_calloc64f(nt);
    for (it=0; it<nt; it++)
    {
        sint[it] = sin(thetas[it]);
    }
    cos3g = memory_calloc64f(ng);
    for (ig=0; ig<ng; ig++)
    {
        cos3g[ig] = cos(3.0*gammas[ig]);
    }
    sum = 0.0;
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk; ik++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            for (it=0; it<nt; it++)
                            {
                                imt = iloc*nm*nb*ng*nk*ns*nt
                                    +      im*nb*ng*nk*ns*nt
                                    +         ib*ng*nk*ns*nt
                                    +            ig*nk*ns*nt
                                    +               ik*ns*nt
                                    +                  is*nt
                                    +                     it;
                                dkdsdt = (dk[ik]*ds[is])*dt[it];
                                dbdg = db[ib]*dg[ig];
                                geom = (twoSinB4[ib]*cos3g[ig])*sint[it];
                                dV5 = geom*(dbdg*dkdsdt);
                                dV2 = dl[iloc]*dm[im];
                                dV7 = dV5*dV2;
                                sum = sum + phi[imt]*dV7;
                            }
                        }
                    }
                }
            }
        }
    }
    memory_free64f(&cos3g);
    memory_free64f(&sint);
    memory_free64f(&twoSinB4);
    memory_free64f(&db);
    memory_free64f(&dg);
    memory_free64f(&dk);
    memory_free64f(&ds);
    memory_free64f(&dt);
    memory_free64f(&dl);
    memory_free64f(&dm);
    return sum;
}
//============================================================================//
int marginal_computeMarginalBeachball(
    const int nloc,
    const int nm,
    const int nb, const double *__restrict__ betas,
    const int ng, const double *__restrict__ gammas,
    const int nk, const double *__restrict__ kappas,
    const int ns, const double *__restrict__ sigmas,
    const int nt, const double *__restrict__ thetas,
    const double *__restrict__ phi
    )
{
    const char *fcnm = "marginal_computeMarginalBeachball\0";
    const double M0loc = 1.0/sqrt(2.0);
    double pAxis[3], nAxis[3], tAxis[3];
    double *bb, *bbAvg, *bbWt, *cos3g, *db,  *dg, *dk, *ds, *dt,
           *twoSinB4, *sint, *xw1, *yw1, dbdg, dkdsdt, dV5, geom, xscal;
    int ib, ig, ik, ierr, iloc, im, imt, indx, is, it, j, jmt, jndx,
        ldi, nmt, nmtBB;
    int8_t *bmap;
    size_t nwork;
    const int nxp = 51;
    const int nyp = nxp;
    const double xc = 1.5;
    const double yc = xc;
    const double rad = 1.0;
    nmt = nm*nb*ng*nk*ns*nt;
    nmtBB = nb*ng*nk*ns*nt;
    ldi = nxp*nyp + computePadding8i(nxp*nyp);
printf("%d %d\n", ldi, nxp*nyp);
    nwork = (size_t)(nmtBB)*(size_t) (ldi);
    if (nwork > INT_MAX)
    {
        printf("%s: Insufficient space for bb\n", fcnm);
        return -1;
    }
    xw1 = memory_calloc64f(nxp*nyp);
    yw1 = memory_calloc64f(nxp*nyp);
    bmap = (int8_t *) calloc(nwork, sizeof(int8_t));
    bbAvg = memory_calloc64f(nxp*nyp);
    bbWt = memory_calloc64f(nxp*nyp);
    bb = memory_calloc64f(nxp*nyp);
    printf("%s: Drawing...\n", fcnm);
    for (ib=0; ib<nb; ib++)
    {
        printf("%s: Drawing beta: %d\n", fcnm, ib+1);
        for (ig=0; ig<ng; ig++)
        {
            for (ik=0; ik<nk; ik++)
            {
                for (is=0; is<ns; is++)
                {
                    for (it=0; it<nt; it++)
                    {
                        postprocess_tt2tnp(betas[ib], gammas[ig], kappas[ik],
                                           sigmas[is], thetas[it],
                                           pAxis, nAxis, tAxis);
                        jmt = ib*ng*nk*ns*nt
                            + ig*nk*ns*nt
                            + ik*ns*nt
                            + is*nt
                            + it;
                        jndx = jmt*ldi;
                        postprocess_tnp2beachballPolarity(nxp, xc, yc, rad,
                                                       pAxis, nAxis, tAxis,
                                                       xw1, yw1, &bmap[jndx]);
                    }
                }
            }
        }
    }
    // Differentiate weights for Riemann quadrature
    db = diff(nb, betas,  &ierr);
    dg = diff(ng, gammas, &ierr);
    dk = diff(nk, kappas, &ierr);
    ds = diff(ns, sigmas, &ierr);
    dt = diff(nt, thetas, &ierr);
    // Compute the geometric factors
    twoSinB4 = memory_calloc64f(nb);
    for (ib=0; ib<nb; ib++)
    {
        twoSinB4[ib] = 2.0*pow(sin(betas[ib]), 4);
    }
    sint = memory_calloc64f(nt);
    for (it=0; it<nt; it++)
    {
        sint[it] = sin(thetas[it]);
    }
    cos3g = memory_calloc64f(ng);
    for (ig=0; ig<ng; ig++)
    {
        cos3g[ig] = cos(3.0*gammas[ig]);
    }
printf("integrating\n");
    // Stack the result
    for (iloc=0; iloc<nloc; iloc++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk; ik++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            for (it=0; it<nt; it++)
                            {
                                imt = im*nb*ng*nk*ns*nt
                                    + ib*ng*nk*ns*nt
                                    + ig*nk*ns*nt
                                    + ik*ns*nt
                                    + is*nt
                                    + it;
                                indx = iloc*nmt + imt;
                                dkdsdt = (dk[ik]*ds[is])*dt[it];
                                dbdg = db[ib]*dg[ig];
                                geom = (twoSinB4[ib]*cos3g[ig])*sint[it];
                                dV5 = geom*(dbdg*dkdsdt);
//if (dV5 < 0.0){printf("error %f %f %f %f %f %f %f\n", dV5, dbdg, dkdsdt, geom, twoSinB4[ib], cos3g[ig], sint[it]);}
                                xscal = (phi[indx] - 1);//*dV5;
                                jmt = ib*ng*nk*ns*nt
                                    + ig*nk*ns*nt
                                    + ik*ns*nt
                                    + is*nt
                                    + it;
                                jndx = jmt*ldi;
                                for (j=0; j<nxp*nxp; j++)
                                {
                                    bb[j] = bb[j] + xscal*(double) bmap[jndx+j];
//                                    bbAvg[j] = bbAvg[j] + (double) bmap[jndx+j];
//                                    bbWt[j] = bbWt[j] + (double) bmap[jndx+j]*dV5;
                                }
                            }
                        }
                    }
                }
            }
        }
    }
printf("dumping result\n");
int ix, iy;
    FILE *fwork;
    fwork = fopen("depmag/beachball.txt", "w");
    for (iy=0; iy<nyp; iy++)
    {
        for (ix=0; ix<nxp; ix++)
        {
            int k = iy*nxp + ix;
            fprintf(fwork, "%e %e %e %e %e %d\n", xw1[k], yw1[k], bb[k], bbAvg[k], bbWt[k], bmap[nmtBB/2*ldi+k]);
        }
        fprintf(fwork, "\n");
    }
    fclose(fwork);
    memory_free64f(&bbAvg);
    memory_free64f(&bbWt);
    memory_free64f(&xw1);
    memory_free64f(&yw1);
    memory_free64f(&cos3g);
    memory_free64f(&sint);
    memory_free64f(&twoSinB4);
    memory_free64f(&db);
    memory_free64f(&dg);
    memory_free64f(&dk);
    memory_free64f(&ds);
    memory_free64f(&dt);
    free(bmap);
    return 0;
}

static double *diff(const int n, const double *__restrict__ x, int *ierr)
{
    double *d;
    int i;
    d = array_set64f(n, 1.0, ierr);
    for (i=0; i<n; i++)
    {
        if (i < n - 1)
        {
            d[i] = x[i+1] - x[i];
        }
        else
        {
            d[i] = x[i] - x[i-1];
        }
    }
    return d;
}

/*! Related to beta */
/*
static double g11(void)
{
    return 1;
}
*/

/*! Related to gamma */
/*
static double g22(const double sinBeta)
{
    return sinBeta*sinBeta;
}
*/

/*! Related to kappa */
/*
static double g33( )
{
    const double twoSqrt3 = 3.4641016151377544;
    coss2 = coss*coss;
    g33 = 0.5*((4.0 + (-2.0 + 3.0*(1.0 - cos2t))*coss2)*cos2g
        + twoSqrt3*sin2g*sins*sin2t)*sin2b;
}
*/

/*! Related to sigma */
/*
static void g44( )
{
    for (ib=0; ib<nb; ib++)
    {
        for (ig=0; ig<ng; ig++)
        {
            g33 = (2.0 - cos(2.0*gamma))*sinb2;
            for (ik=0; ik<nk; ik++)
            {
                for (is=0; is<ns; is++)
                {
                    for (it=0; it<nt; it++)
                    {

                    }
                }
            }
        }
    }
    return;
}
*/

/*!
 * @brief Computes the determinant of the upper 3 x 3 matrix
 */
/*
static void computeDet123(const int nmt,
                          const double *__restrict__ g11,
                          const double *__restrict__ g22,
                          const double *__restrict__ g33,
                          double *__restrict__ det)
{
    int i;
    for (i=0; i<nmt; i++)
    {
        det[i] = (g11[i]*g22[i])*g33;
    }
    return
}

static void computeDet124(const int nmt,
                          const double *__restrict__ g11,
                          const double *__restrict__ g22,
                          const double *__restrict__ g44,
                          double *__restrict__ det)
{
    computeDet123(nmt, g11, g22, g44, det);
    return;
}

static void computeDet125(const int nmt,
                          const double *__restrict__ g11,
                          const double *__restrict__ g22,
                          const double *__restrict__ g55,
                          double *__restrict__ det)
{
    computeDet123(nmt, g11, g22, g55, det);
    return;
}
*/

int marginal_write1DToGnuplot(const char *fname,
                              const int nx, const double *__restrict__ xlocs,
                              const double *__restrict__ vals)
{
    FILE *fout;
    int ix;
    fout = fopen(fname, "w");
    for (ix=0; ix<nx; ix++)
    {
        fprintf(fout, "%.8e %.10e\n", xlocs[ix], vals[ix]);
    }
    fclose(fout);
    return 0;
}

int marginal_write2DToGnuplot(const char *fname,
                              const int nx, const double *__restrict__ xlocs,
                              const int ny, const double *__restrict__ ylocs,
                              const double *__restrict__ vals)
{
    FILE *fout;
    int ix, ixy, iy;
    fout = fopen(fname, "w");
    for (iy=0; iy<ny; iy++)
    {
        for (ix=0; ix<nx; ix++)
        {
            ixy = iy*nx + ix;
            fprintf(fout, "%.8e %.8e %.10e\n", xlocs[ix], ylocs[iy], vals[ixy]);
        }
        fprintf(fout, "\n");
    }
    fclose(fout);
    return 0;
}
//============================================================================//
static int computePadding8i(const int n)
{
    size_t mod, pad;
    int ipad;
    // Set space and make G matrix
    pad = 0;
    mod = ((size_t) n*sizeof(int8_t))%64;
    if (mod != 0)
    {
        pad = (64 - mod)/sizeof(int8_t);
    }
    ipad = (int) pad;
    return ipad;
}


text:
/* linalg/test_common.c
 *
 * Copyright (C) 2017, 2018, 2019, 2020 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>

static int create_random_vector(gsl_vector * v, gsl_rng * r);
static int create_random_matrix(gsl_matrix * m, gsl_rng * r);
static int create_posdef_matrix(gsl_matrix * m, gsl_rng * r);
static int create_hilbert_matrix2(gsl_matrix * m);

static int
create_random_vector(gsl_vector * v, gsl_rng * r)
{
  const size_t N = v->size;
  size_t i;

  for (i = 0; i < N; ++i)
    {
      double vi = gsl_rng_uniform(r);
      gsl_vector_set(v, i, vi);
    }

  return GSL_SUCCESS;
}

static int
create_random_complex_vector(gsl_vector_complex * v, gsl_rng * r)
{
  const size_t N = v->size;
  size_t i;

  for (i = 0; i < N; ++i)
    {
      gsl_complex vi;
      GSL_REAL(vi) = gsl_rng_uniform(r);
      GSL_IMAG(vi) = gsl_rng_uniform(r);
      gsl_vector_complex_set(v, i, vi);
    }

  return GSL_SUCCESS;
}

static int
create_random_matrix(gsl_matrix * m, gsl_rng * r)
{
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double mij = gsl_rng_uniform(r);
          gsl_matrix_set(m, i, j, mij);
        }
    }

  return GSL_SUCCESS;
}

static int
create_random_complex_matrix(gsl_matrix_complex * m, gsl_rng * r)
{
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex mij;

          GSL_REAL(mij) = gsl_rng_uniform(r);
          GSL_IMAG(mij) = gsl_rng_uniform(r);

          gsl_matrix_complex_set(m, i, j, mij);
        }
    }

  return GSL_SUCCESS;
}

static int
create_symm_matrix(gsl_matrix * m, gsl_rng * r)
{
  const size_t N = m->size1;
  size_t i, j;

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j <= i; ++j)
        {
          double mij = gsl_rng_uniform(r);
          gsl_matrix_set(m, i, j, mij);
        }
    }

  /* copy lower triangle to upper */
  gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, m, m);

  return GSL_SUCCESS;
}

static int
create_herm_matrix(gsl_matrix_complex * m, gsl_rng * r)
{
  const size_t N = m->size1;
  size_t i, j;

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j <= i; ++j)
        {
          double re = gsl_rng_uniform(r);
          double im = (i != j) ? gsl_rng_uniform(r) : 0.0;
          gsl_complex z = gsl_complex_rect(re, im);

          gsl_matrix_complex_set(m, i, j, z);

          if (i != j)
            gsl_matrix_complex_set(m, j, i, gsl_complex_conjugate(z));
        }
    }

  return GSL_SUCCESS;
}

/* create symmetric banded matrix with p sub/super-diagonals */
static int
create_symm_band_matrix(const size_t p, gsl_matrix * m, gsl_rng * r)
{
  size_t i;

  gsl_matrix_set_zero(m);

  for (i = 0; i < p + 1; ++i)
    {
      gsl_vector_view subdiag = gsl_matrix_subdiagonal(m, i);
      create_random_vector(&subdiag.vector, r);

      if (i > 0)
        {
          gsl_vector_view superdiag = gsl_matrix_superdiagonal(m, i);
          gsl_vector_memcpy(&superdiag.vector, &subdiag.vector);
        }
    }

  return GSL_SUCCESS;
}

/* create (p,q) banded matrix */
static int
create_band_matrix(const size_t p, const size_t q, gsl_matrix * m, gsl_rng * r)
{
  size_t i;

  gsl_matrix_set_zero(m);

  for (i = 0; i <= p; ++i)
    {
      gsl_vector_view v = gsl_matrix_subdiagonal(m, i);
      create_random_vector(&v.vector, r);
    }

  for (i = 1; i <= q; ++i)
    {
      gsl_vector_view v = gsl_matrix_superdiagonal(m, i);
      create_random_vector(&v.vector, r);
    }

  return GSL_SUCCESS;
}

static int
create_posdef_matrix(gsl_matrix * m, gsl_rng * r)
{
  const size_t N = m->size1;
  const double alpha = 10.0 * N;
  size_t i;

  /* The idea is to make a symmetric diagonally dominant
   * matrix. Make a symmetric matrix and add alpha*I to
   * its diagonal */

  create_symm_matrix(m, r);

  for (i = 0; i < N; ++i)
    {
      double mii = gsl_matrix_get(m, i, i);
      gsl_matrix_set(m, i, i, mii + alpha);
    }

  return GSL_SUCCESS;
}

static int
create_posdef_complex_matrix(gsl_matrix_complex *m, gsl_rng *r)
{
  const size_t N = m->size1;
  const double alpha = 10.0 * N;
  size_t i;

  create_herm_matrix(m, r);

  for (i = 0; i < N; ++i)
    {
      gsl_complex * mii = gsl_matrix_complex_ptr(m, i, i);
      GSL_REAL(*mii) += alpha;
    }

  return GSL_SUCCESS;
}

static int
create_hilbert_matrix2(gsl_matrix * m)
{
  const size_t N = m->size1;
  size_t i, j;

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
        }
    }

  return GSL_SUCCESS;
}

static int
create_posdef_band_matrix(const size_t p, gsl_matrix * m, gsl_rng * r)
{
  const size_t N = m->size1;
  const double alpha = 10.0 * N;
  size_t i;

  /* The idea is to make a symmetric diagonally dominant
   * matrix. Make a symmetric matrix and add alpha*I to
   * its diagonal */

  create_symm_band_matrix(p, m, r);

  for (i = 0; i < N; ++i)
    {
      double *mii = gsl_matrix_ptr(m, i, i);
      *mii += alpha;
    }

  return GSL_SUCCESS;
}

/* transform dense symmetric banded matrix to compact form, with bandwidth p */
static int
symm2band_matrix(const size_t p, const gsl_matrix * m, gsl_matrix * bm)
{
  const size_t N = m->size1;

  if (bm->size1 != N)
    {
      GSL_ERROR("banded matrix requires N rows", GSL_EBADLEN);
    }
  else if (bm->size2 != p + 1)
    {
      GSL_ERROR("banded matrix requires p + 1 columns", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      gsl_matrix_set_zero(bm);

      for (i = 0; i < p + 1; ++i)
        {
          gsl_vector_const_view diag = gsl_matrix_const_subdiagonal(m, i);
          gsl_vector_view v = gsl_matrix_subcolumn(bm, i, 0, N - i);

          gsl_vector_memcpy(&v.vector, &diag.vector);
        }

      return GSL_SUCCESS;
    }
}

/* transform general dense (p,q) banded matrix to compact banded form */
static int
gen2band_matrix(const size_t p, const size_t q, const gsl_matrix * A, gsl_matrix * AB)
{
  const size_t N = A->size2;

  if (AB->size1 != N)
    {
      GSL_ERROR("banded matrix requires N rows", GSL_EBADLEN);
    }
  else if (AB->size2 != 2*p + q + 1)
    {
      GSL_ERROR("banded matrix requires 2*p + q + 1 columns", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      gsl_matrix_set_zero(AB);

      /* copy diagonal and subdiagonals */
      for (i = 0; i <= p; ++i)
        {
          gsl_vector_const_view v = gsl_matrix_const_subdiagonal(A, i);
          gsl_vector_view w = gsl_matrix_subcolumn(AB, p + q + i, 0, v.vector.size);
          gsl_vector_memcpy(&w.vector, &v.vector);
        }

      /* copy superdiagonals */
      for (i = 1; i <= q; ++i)
        {
          gsl_vector_const_view v = gsl_matrix_const_superdiagonal(A, i);
          gsl_vector_view w = gsl_matrix_subcolumn(AB, p + q - i, i, v.vector.size);
          gsl_vector_memcpy(&w.vector, &v.vector);
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
// =============================================================================
// == IMG.h
// == --------------------------------------------------------------------------
// == An image class to be used with temporal superpixels.
// ==
// == All work using this code should cite:
// == J. Chang, D. Wei, and J. W. Fisher III. A Video Representation Using
// ==    Temporal Superpixels. CVPR 2013.
// == --------------------------------------------------------------------------
// == Written by Jason Chang and Donglai Wei 06-20-2013
// =============================================================================

#ifndef _IMG_H_INCLUDED_
#define _IMG_H_INCLUDED_
#include <iostream>
#include <fstream>
#include "utils.h"
#include "SP.h"
#include "NormalD.h"
#include "topology.h"
#include "array.h"
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

#include "helperMEX.h"

using namespace std;

class IMG
{
private:
   int debug;
   int xdim;
   int ydim;
   int N;

   arr(double) data;
   arr(int) label;
   int K;
   //NIW new_pos;
   NormalD new_pos;
   NormalD new_app;
   arr(SP*) SP_arr;
   SP* SP_new; // a pointer to an exemplar "new" SP
   arr(bool) SP_old; // indicates the SPs that are from a prev frame
   arr(bool) SP_changed; // indicates which SPs changed in the previous iteration

   arr(linkedListNode<int>*) border_ptr;
   arr(linkedListNode<int>*) pixel_ptr;

   arr(double) T4Table;

   double log_alpha; // size parameter
   double log_beta;
   double area;
   double area_var;
   double log_area_var;
   unsigned long max_UID;

   double lambda; // flow parameters
   double lambda_sigma;

   arr(int) prev_label;
   arr(double) prev_app_mean;
   arr(double) prev_pos_mean;
   int prev_K;
   arr(double) prev_covariance;
   arr(double) prev_precision;
   linkedList< std::pair<int, int> > prev_neighbors;
   linkedList< double > prev_neighbors_scale;
   bool alive_dead_changed;

   double dummy_log_prob;
   arr(bool) boundary_mask;

   // flow stuff
   gsl_vector *vgsl;
   gsl_vector *ugsl;
   arr(double) Sxy;
   arr(double) Syy;
   arr(double) SxySyy;
   arr(double) obs_u;
   arr(double) obs_v;
   arr(double) temp_Syy;
   arr(double) temp_d;
   arr(bool) temp_still_alive;
   int num_alive;
   int num_dead;

   arr(int) all2alive;
   arr(int) alive2all;
   arr(int) all2dead;
   arr(int) dead2all;

public:
   IMG();
   virtual ~IMG();

   void ReadIMG();
   void mxReadIMG(const mxArray *mstruct);
   void mxWriteIMG(mxArray* plhs[], const mxArray* oldstruct);

   // --------------------------------------------------------------------------
   // -- move_local
   // --   finds the optimal local joint move of labels and parameters. Chooses
   // -- a super pixel at random and loops through and updates its borders.
   // --------------------------------------------------------------------------
   bool move_local_IMG();

   bool move_switch_IMG();

   // --------------------------------------------------------------------------
   // -- move_merge
   // --   finds the optimal merge between all pairs of super pixels
   // --------------------------------------------------------------------------
   void move_merge_IMG();
   void move_split_IMG();
   void move_IMG();

   void Flow_QP2();
   void find_SxySyy(arr(int) all2alive, arr(int) alive2all, arr(int) all2dead, arr(int) dead2all);
private:
   //
   // -- 1. I/O
   void U_initialize();
   void U_relabel_SP(bool final=false);

   double U_calc_energy();
   double U_calc_model_order(int N, bool is_old);

   // -- 3. border bookkeeping
   bool U_check_border_pix(int index);
   bool U_check_border_pix(int index, int cur_label);
   void U_find_border_SP(int k,arr(bool) neighbors,linkedList<int> &neighborsLL);
   // --------------------------------------------------------------------------
   // -- U_update_border_changed
   // --   Updates all border linked lists and pointers for the neighbors of a
   // -- changed pixel at index.
   // --
   // --   parameters:
   // --     - index : the index of the recently changed pixel
   // --------------------------------------------------------------------------
   void U_update_border_changed(int index);
   // --------------------------------------------------------------------------
   // -- U_update_border_changed_pixel
   // --   Updates all border linked lists and pointers at index
   // --
   // --   parameters:
   // --     - index : the index of the recently changed pixel
   // --------------------------------------------------------------------------
   void U_update_border_changed_pixel(int index);

   // --------------------------------------------------------------------------
   // -- U_update_neighbor_list
   // --   Updates super pixel neighbors used in merge
   // --
   // --   parameters:
   // --     - neighbors : a boolean array indicating which labels are added
   // --     - neighborsLL : a linked list of the neighbor indices
   // --     - index : the index of the point to consider adding
   // --------------------------------------------------------------------------
   void U_update_neighbor_list(arr(bool) neighbors, linkedList<int> &neighborsLL, int index);

   void U_fix_neighbors_self(int k);
   void U_fix_neighbors_neighbors(int k);
   void U_fix_neighbors_neighbors(int k, int kignore);
   void U_print_neighbors(int k);

   // --------------------------------------------------------------------------
   // -- U_update_neighbors_rem
   // --   Updates all neighbor lists when a pixel is "removed". A subsequent
   // -- call to update_neighbors_add should be completed right after this one.
   // -- The neighboring label should be changed before calling this function.
   // --
   // --   parameters:
   // --     - old_label : the label of the pixel before it was changed
   // --     - index : the index bordering the removed pixel
   // --------------------------------------------------------------------------
   void U_update_neighbors_rem(int old_label, int index);

   // --------------------------------------------------------------------------
   // -- U_update_neighbors_add
   // --   Updates all neighbor lists when a pixel is "added". A previous
   // -- call to update_neighbors_rem should be completed right before this one.
   // -- The neighboring label should be changed before calling this function.
   // --
   // --   parameters:
   // --     - index : the index bordering the removed pixel
   // --------------------------------------------------------------------------
   void U_update_neighbors_add(int index);

   // --------------------------------------------------------------------------
   // -- U_update_neighbors_merge
   // --   Updates the neighbor lists for merging two super pixels.
   // -- All labels and SP_arr things should be updated *before* calling this.
   // --
   // --   parameters:
   // --     - index : the index bordering the removed pixel
   // --------------------------------------------------------------------------
   void U_update_neighbors_merge(int new_label, int old_label);

   // --------------------------------------------------------------------------
   // -- U_update_neighbors_split
   // --   Updates the neighbor lists for merging two super pixels.
   // -- All labels and SP_arr things should be updated *before* calling this.
   // --
   // --   parameters:
   // --     - index : the index bordering the removed pixel
   // --------------------------------------------------------------------------
   void U_update_neighbors_split(int label1, int label2);

   // -- 4. Modularized Move
   // -- 4.0 Change of Energy

   // --------------------------------------------------------------------------
   // -- link_cost
   // --   calculates the energy of linking SP_arr[ko] (old) to SP_arr[kn] (new)
   // --
   // --   parameters:
   // --     - ko : the old k
   // --     - kn : the new k
   // --------------------------------------------------------------------------
   double link_cost(int ko, int kn);

   double move_local_calc_delta_MM(int index, int new_k, double& max_prob, int& max_k);

   // --------------------------------------------------------------------------
   // -- move_local_calc_neigbor
   // --   calculates the probability of assigning the pixel at index to the
   // -- cluster of nindex (neighbor index). Updates the max_prob and max_k if
   // -- the new_app value is greater than the old one
   // --
   // --   parameters:
   // --     - index : the new point to add
   // --     - new_k : the neighbor to add this point to
   // --     - max_prob : the maximum probability of all neighbors
   // --     - max_k : the super pixel index of the maximum probability
   // --------------------------------------------------------------------------
   double move_local_calc_delta(int index, int new_k, bool add, double& max_prob, int& max_k);

   double move_switch_calc_delta(SP* &oldSP, SP* &newSP);

   // --------------------------------------------------------------------------
   // -- move_merge_calc_delta
   // --   calculates the probability of assigning the pixel at index to the
   // -- cluster of nindex (neighbor index).
   // --
   // --   parameters:
   // --     - index : the new point to add
   // --     - nindex : the neighbor to add this point to
   // --------------------------------------------------------------------------
   double move_merge_calc_delta(int k, int merge_k);

   // --------------------------------------------------------------------------
   // -- move_split_calc_delta
   // --   calculates the change in energy for
   // -- (k1 U new_k1) && (k2 U new_k2) - (k1 U new_k1 U new_k) && (k2)
   // --
   // --   parameters:
   // --     - SP1 : the SP that originates the splitting
   // --     - SP2 : the SP to split to
   // --     - new_SP1 : temporary SP that contains pixels that will go in k1
   // --     - new_SP2 : temporary SP that contains pixels that will go in k2
   // --     - SP1_old : indicates if SP1 is an old SP
   // --     - SP2_old : indicates if SP2 is an old SP
   // --------------------------------------------------------------------------
   double move_split_calc_delta(SP* SP1, SP* SP2, SP* new_SP1, SP* new_SP2, bool SP1_old, bool SP2_old);

   // -- 4.1 move local
   // need to be done
   int move_local_SP_region(int k,linkedList<int> &check_labels);
   int move_local_SP(int k);
   // need to be done
   int move_local_pix(int index);
   //-- 4.2 move merge
   void move_merge_SP_propose(int k,arr(bool) neighbors,double &max_E,int &max_k);
   void move_merge_SP_propose_region(int k,arr(bool) neighbors, linkedList<int> &check_labels,double &max_E,int &max_k);
   void move_merge_SP(int k,arr(bool) neighbors);

   // -- 4.3 move split
   void move_split_SP(int k);
   bool U_Kmeans_plusplus(int k, int *bbox, int num_SP, int numiter, int index2=-1);
   void U_connect_newSP(int *bbox, int num_SP);
   double U_dist(int index1 , double* center);
   double U_dist(int index1 , int index2);

   void move_split_SP_propose(int index, int num_SP, int option, double &max_E,int &ksplit,int* new_ks);

};

#endif

ext: h


text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file genetic.c
 * \brief Source file to define genetic algorithm main function.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#if HAVE_MPI
#include <mpi.h>
#endif
#include "bits.h"
#include "entity.h"
#include "population.h"
#include "reproduction.h"
#include "selection.h"
#include "evolution.h"
#include "genetic.h"

#define DEBUG_GENETIC 0         ///< Macro to debug the genetic functions.

static Population genetic_population[1];
///< Population of the genetic algorithm.
static double (*genetic_simulation) (Entity *);
///< Pointer to the function to perform a simulation.

/**
 * Function to get a variable encoded in the genome of an entity.
 *
 * \return Variable value.
 */
double
genetic_get_variable (Entity * entity,  ///< Entity struct.
                      GeneticVariable * variable)       ///< Variable data.
{
  double x;
#if DEBUG_GENETIC
  fprintf (stderr, "get variable: start\n");
#endif
  x = variable->minimum
    + bit_get_value (entity->genome, variable->location, variable->nbits)
    * (variable->maximum - variable->minimum)
    / ((unsigned long long int) 1L << variable->nbits);
#if DEBUG_GENETIC
  fprintf (stderr, "get variable: value=%lg\n", x);
  fprintf (stderr, "get variable: end\n");
#endif
  return x;
}

/**
 * Funtion to apply the evolution of a population.
 */
static inline void
genetic_evolution (Population * population,     ///< Population
                   gsl_rng * rng)       ///< GSL random numbers generator.
{
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_evolution: start\n");
#endif
  evolution_mutation (population, rng);
  evolution_reproduction (population, rng);
  evolution_adaptation (population, rng);
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_evolution: end\n");
#endif
}

/**
 * Funtion to perform the simulations on a thread.
 */
static void
genetic_simulation_thread (GeneticThreadData * data)    ///< Thread data.
{
  unsigned int i;
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_thread: start\n");
  fprintf (stderr, "genetic_simulation_thread: nmin=%u nmax=%u\n",
           data->nmin, data->nmax);
#endif
  for (i = data->nmin; i < data->nmax && !genetic_population->stop; ++i)
    {
      genetic_population->objective[i]
        = genetic_simulation (genetic_population->entity + i);
      if (genetic_population->objective[i] < genetic_population->threshold)
        {
          g_mutex_lock (mutex);
          genetic_population->stop = 1;
          g_mutex_unlock (mutex);
          break;
        }
    }
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_thread: end\n");
#endif
}

/**
 * Function to perform the simulations on a task.
 */
static void
genetic_simulation_master (unsigned int nsurvival)
                           ///< Number of survival entities already simulated.
{
  unsigned int j, nsimulate, nmin, nmax;
  GThread *thread[nthreads];
  GeneticThreadData thread_data[nthreads];
  Population *population;
#if HAVE_MPI
  unsigned int i;
  unsigned int n[ntasks + 1];
  unsigned int stop[ntasks];
  char *genome_array;
  MPI_Status mpi_status;
#endif

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_master: start\n");
#endif

  population = genetic_population;
  nmax = population->nentities;
  nsimulate = nmax - nsurvival;
  nmin = nsurvival;
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_master: nmax=%u nmin=%u nsimulate=%u\n",
           nmax, nmin, nsimulate);
#endif

#if HAVE_MPI
  nmax = nmin + nsimulate / ntasks;

  // Send genome information to the slaves
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_master: nmax=%u nmin=%u nsimulate=%u\n",
           nmax, nmin, nsimulate);
#endif

  genome_array = (char *) g_malloc (nsimulate * population->genome_nbytes);
  for (j = 0; j < nsimulate; ++j)
    memcpy (genome_array + j * population->genome_nbytes,
            population->entity[nmin + j].genome, population->genome_nbytes);
  n[0] = 0;
  for (i = 0; (int) ++i < ntasks;)
    {
      n[i] = i * nsimulate / ntasks;
      j = (i + 1) * nsimulate / ntasks - n[i];
#if DEBUG_GENETIC
      fprintf (stderr,
               "genetic_simulation_master: send %u bytes on %u to task %u\n",
               j * population->genome_nbytes, nsurvival + n[i], i);
#endif
      MPI_Send (genome_array + n[i] * population->genome_nbytes,
                j * population->genome_nbytes, MPI_CHAR, i, 1, MPI_COMM_WORLD);
    }
  n[i] = nsimulate;

#if DEBUG_GENETIC
  for (j = 0; j <= ntasks; ++j)
    fprintf (stderr, "genetic_simulation_master: n[%u]=%u\n", j, n[j]);
#endif
  g_free (genome_array);

  nsimulate = nmax - nmin;

#endif

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_master: nsimulate=%u\n", nsimulate);
  fprintf (stderr, "genetic_simulation_master: performing simulations\n");
#endif

  thread_data[0].nmin = nmin;
  for (j = 0; ++j < nthreads;)
    thread_data[j - 1].nmax = thread_data[j].nmin
      = nmin + j * nsimulate / nthreads;
  thread_data[j - 1].nmax = nmax;
  for (j = 0; j < nthreads; ++j)
    thread[j] = g_thread_new
      (NULL, (GThreadFunc) (void (*)(void)) genetic_simulation_thread,
       thread_data + j);
  for (j = 0; j < nthreads; ++j)
    g_thread_join (thread[j]);

#if HAVE_MPI
  // Receive objective function resuts from the slaves
  for (j = 0; (int) ++j < ntasks;)
    {
#if DEBUG_GENETIC
      fprintf (stderr,
               "genetic_simulation_master: "
               "receive %u reals from task %u on %u\n",
               n[j + 1] - n[j], j, nsurvival + n[j]);
#endif
      MPI_Recv (population->objective + nsurvival + n[j], n[j + 1] - n[j],
                MPI_DOUBLE, j, 1, MPI_COMM_WORLD, &mpi_status);
#if DEBUG_GENETIC
      fprintf (stderr,
               "genetic_simulation_master: receive one integer from task %u\n",
               j);
#endif
      MPI_Recv (stop + j, j, MPI_UNSIGNED, j, 1, MPI_COMM_WORLD, &mpi_status);
      if (stop[j])
        genetic_population->stop = 1;
    }
  // Sending stop instruction to the slaves
  for (j = 0; (int) ++j < ntasks;)
    {
#if DEBUG_GENETIC
      fprintf (stderr,
               "genetic_simulation_master: sending one integer to task %u\n",
               j);
#endif
      MPI_Send (&genetic_population->stop, 1, MPI_UNSIGNED, j, 1,
                MPI_COMM_WORLD);
    }
#endif

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_master: end\n");
#endif
}

#if HAVE_MPI

/**
 * Function to perform the simulations on a task.
 */
static void
genetic_simulation_slave (unsigned int nsurvival,
                          ///< Number of survival entities already simulated.
                          int rank)     ///< Number of task.
{
  unsigned int j, nsimulate, nmin, nmax, stop;
  GThread *thread[nthreads];
  GeneticThreadData thread_data[nthreads];
  Population *population;
  char *genome_array;
  MPI_Status mpi_status;

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_slave: rank=%d start\n", rank);
#endif

  population = genetic_population;
  nmax = population->nentities;
  nsimulate = nmax - nsurvival;
  nmin = nsurvival;
  nmax = nmin + (rank + 1) * nsimulate / ntasks;
  nmin += rank * nsimulate / ntasks;

  // Receive genome information from the master
  nsimulate = nmax - nmin;
  genome_array = (char *) g_malloc (nsimulate * population->genome_nbytes);
#if DEBUG_GENETIC
  fprintf (stderr,
           "genetic_simulation_slave: rank=%d receive %u bytes from master\n",
           rank, nsimulate * population->genome_nbytes);
#endif
  MPI_Recv (genome_array, nsimulate * population->genome_nbytes, MPI_CHAR, 0,
            1, MPI_COMM_WORLD, &mpi_status);
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_slave: rank=%d nmin=%u nsimulate=%u\n",
           rank, nmin, nsimulate);
#endif
  for (j = 0; j < nsimulate; ++j)
    memcpy (population->entity[nmin + j].genome,
            genome_array + j * population->genome_nbytes,
            population->genome_nbytes);
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_slave: rank=%d freeing\n", rank);
#endif
  g_free (genome_array);

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_slave: rank=%d performing simulations\n",
           rank);
#endif
  thread_data[0].nmin = nmin;
  for (j = 0; ++j < nthreads;)
    thread_data[j - 1].nmax = thread_data[j].nmin
      = nmin + j * nsimulate / nthreads;
  thread_data[j - 1].nmax = nmax;
  for (j = 0; j < nthreads; ++j)
    thread[j] = g_thread_new
      (NULL, (GThreadFunc) (void (*)(void)) genetic_simulation_thread,
       thread_data + j);
  for (j = 0; j < nthreads; ++j)
    g_thread_join (thread[j]);

#if DEBUG_GENETIC
  fprintf (stderr,
           "genetic_simulation_slave: rank=%d send %u reals on %u to master\n",
           rank, nsimulate, nmin);
#endif
  // Send the objective function resuts to the master
  MPI_Send (population->objective + nmin, nsimulate, MPI_DOUBLE, 0, 1,
            MPI_COMM_WORLD);
#if DEBUG_GENETIC
  fprintf (stderr,
           "genetic_simulation_slave: rank=%d send one integer to master\n",
           rank);
#endif
  // Send the stop variable from the master
  MPI_Send (&genetic_population->stop, 1, MPI_UNSIGNED, 0, 1, MPI_COMM_WORLD);
#if DEBUG_GENETIC
  fprintf (stderr,
           "genetic_simulation_slave: "
           "rank=%d receive one integer from master\n", rank);
#endif
  // Receive the stop variable from the master
  MPI_Recv (&stop, 1, MPI_UNSIGNED, 0, 1, MPI_COMM_WORLD, &mpi_status);
  if (stop)
    genetic_population->stop = 1;

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_simulation_slave: rank=%d end\n", rank);
#endif
}

#endif

/**
 * Function to create the data of the genetic algorithm.
 *
 * \return 1 on succes, 0 on error.
 */
int
genetic_new (unsigned int nvariables,   ///< Number of variables.
             GeneticVariable * variable,        ///< Array of variables data.
             unsigned int nentities,
             ///< Number of entities in each generation.
             double mutation_ratio,     ///< Mutation ratio.
             double reproduction_ratio, ///< Reproduction ratio.
             double adaptation_ratio,   ///< Adaptation ratio.
             double threshold)  ///< Threshold to finish the simulations.
{
  unsigned int i, genome_nbits, nprocesses;

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_new: start\n");
#endif

  // Checking variables number
  if (!nvariables)
    {
      fprintf (stderr, "ERROR: no variables\n");
      return 0;
    }

  // Checking variable bits number
  for (i = genome_nbits = 0; i < nvariables; ++i)
    {
      if (!variable[i].nbits || variable[i].nbits > 32)
        {
          fprintf (stderr, "ERROR: bad bits number in variable %u\n", i + 1);
          return 0;
        }
      variable[i].location = genome_nbits;
      genome_nbits += variable[i].nbits;
    }

  // Checking processes number
  nprocesses = ntasks * nthreads;
  if (!nprocesses)
    {
      fprintf (stderr, "ERROR: no processes\n");
      return 0;
    }

  // Init the population
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_new: init the population\n");
#endif
  if (!population_new (genetic_population, variable, nvariables, genome_nbits,
                       nentities, mutation_ratio, reproduction_ratio,
                       adaptation_ratio, threshold))
    return 0;

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_new: end\n");
#endif
  return 1;
}

/**
 * Function to perform the genetic algorithm.
 *
 * \return 1 on succes, 0 on error.
 */
int
genetic_algorithm (unsigned int nvariables,     ///< Number of variables.
                   GeneticVariable * variable,  ///< Array of variables data.
                   unsigned int nentities,
                   ///< Number of entities in each generation.
                   unsigned int ngenerations,   ///< Number of generations.
                   double mutation_ratio,       ///< Mutation ratio.
                   double reproduction_ratio,   ///< Reproduction ratio.
                   double adaptation_ratio,     ///< Adaptation ratio.
                   const gsl_rng_type * type_random,
                   ///< Type of GSL random numbers generator algorithm.
                   unsigned long random_seed,
                   ///< Seed of the GSL random numbers generator.
                   unsigned int type_reproduction,
                   ///< Type of reproduction algorithm.
                   unsigned int type_selection_mutation,
                   ///< Type of mutation selection algorithm.
                   unsigned int type_selection_reproduction,
                   ///< Type of reproduction selection algorithm.
                   unsigned int type_selection_adaptation,
                   ///< Type of adaptation selection algorithm.
                   double threshold,
                   ///< Threshold to finish the simulations.
                   double (*simulate_entity) (Entity *),
///< Pointer to the function to perform a simulation of an entity.
                   char **best_genome,  ///< Best genome.
                   double **best_variables,     ///< Best array of variables.
                   double *best_objective)
                   ///< Best objective function value.
{
  unsigned int i;
  double *bv;
  gsl_rng *rng;
  Entity *best_entity;
#if HAVE_MPI
  int rank;
#endif

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_algorithm: start\n");
#endif

  // Init the data
  if (!genetic_new (nvariables, variable, nentities,
                    mutation_ratio, reproduction_ratio, adaptation_ratio,
                    threshold))
    return 0;

  // Init the evaluation function
  genetic_simulation = simulate_entity;

  // Get the rank
#if HAVE_MPI
  MPI_Comm_rank (MPI_COMM_WORLD, &rank);
#endif

  // Variables to init only by master task
#if HAVE_MPI
  if (rank == 0)
    {
#endif
      // Init the GSL random numbers generator
      rng = gsl_rng_alloc (type_random);
      if (random_seed)
        gsl_rng_set (rng, random_seed);

      // Init genomes
      population_init_genomes (genetic_population, rng);

      // Init selection, mutation and reproduction algorithms
      reproduction_init (type_reproduction);
      selection_init (type_selection_mutation, type_selection_reproduction,
                      type_selection_adaptation);

#if HAVE_MPI
    }
  if (rank == 0)
    {
#endif

      // First simulation of all entities
#if DEBUG_GENETIC
      fprintf (stderr, "genetic_algorithm: first simulation of all entities\n");
#endif
      genetic_simulation_master (0);

      // Sorting by objective function results
#if DEBUG_GENETIC
      fprintf (stderr,
               "genetic_algorithm: sorting by objective function results\n");
#endif
      evolution_sort (genetic_population);

      // Population generations
#if DEBUG_GENETIC
      fprintf (stderr, "genetic_algorithm: ngenerations=%u\n", ngenerations);
#endif
      for (i = 1; i < ngenerations && !genetic_population->stop; ++i)
        {
          // Evolution
#if DEBUG_GENETIC
          fprintf (stderr, "genetic_algorithm: evolution of the population\n");
#endif
          genetic_evolution (genetic_population, rng);

          // Simulation of the new entities
#if DEBUG_GENETIC
          fprintf (stderr,
                   "genetic_algorithm: simulation of the new entities\n");
#endif
          genetic_simulation_master (genetic_population->nsurvival);

          // Sorting by objective function results
#if DEBUG_GENETIC
          fprintf (stderr,
                   "genetic_algorithm: sorting by objective function results\n");
#endif
          evolution_sort (genetic_population);
        }

      // Saving the best
#if DEBUG_GENETIC
      fprintf (stderr, "genetic_algorithm: saving the best\n");
#endif
      best_entity = genetic_population->entity;
      *best_objective = genetic_population->objective[0];
      *best_genome = (char *) g_malloc (genetic_population->genome_nbytes);
      memcpy (*best_genome, best_entity->genome,
              genetic_population->genome_nbytes);
      *best_variables = bv = (double *) g_malloc (nvariables * sizeof (double));
      for (i = 0; i < nvariables; ++i)
        bv[i] = genetic_get_variable (best_entity, variable + i);
      gsl_rng_free (rng);

#if HAVE_MPI

    }
  else
    {
      // First simulation of all entities
#if DEBUG_GENETIC
      fprintf (stderr, "genetic_algorithm: first simulation of all entities\n");
#endif
      genetic_simulation_slave (0, rank);

      // Population generations
      for (i = 1; i < ngenerations && !genetic_population->stop; ++i)
        {
          // Simulation of the new entities
#if DEBUG_GENETIC
          fprintf (stderr,
                   "genetic_algorithm: simulation of the new entities\n");
#endif
          genetic_simulation_slave (genetic_population->nsurvival, rank);
        }
    }

#endif

  // Freeing memory
#if DEBUG_GENETIC
  fprintf (stderr, "genetic_algorithm: freeing memory\n");
#endif
  population_free (genetic_population);

#if DEBUG_GENETIC
  fprintf (stderr, "genetic_algorithm: rank=%d\n", rank);
  fprintf (stderr, "genetic_algorithm: end\n");
#endif
  return 1;
}

/**
 * Function to perform the genetic algorithm with default random and evolution
 *   algorithms.
 *
 * \return 1 on succes, 0 on error.
 */
int
genetic_algorithm_default (unsigned int nvariables,
                           ///< Number of variables.
                           GeneticVariable * variable,
                           ///< Array of variables data.
                           unsigned int nentities,
                           ///< Number of entities in each generation.
                           unsigned int ngenerations,
                           ///< Number of generations.
                           double mutation_ratio,       ///< Mutation ratio.
                           double reproduction_ratio,   ///< Reproduction ratio.
                           double adaptation_ratio,     ///< Adaptation ratio.
                           unsigned long random_seed,
                           ///< Seed of the GSL random numbers generator.
                           double threshold,
                           ///< Threshold to finish the simulations.
                           double (*simulate_entity) (Entity *),
///< Pointer to the function to perform a simulation of an entity.
                           char **best_genome,  ///< Best genome.
                           double **best_variables,
                           ///< Best array of variables.
                           double *best_objective)
                           ///< Best objective function value.
{
  return genetic_algorithm (nvariables,
                            variable,
                            nentities,
                            ngenerations,
                            mutation_ratio,
                            reproduction_ratio,
                            adaptation_ratio,
                            gsl_rng_mt19937,
                            random_seed,
                            0,
                            0,
                            0,
                            0,
                            threshold,
                            simulate_entity,
                            best_genome, best_variables, best_objective);
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


text:
#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <ctype.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
//#include <pthread.h>
#include <omp.h>
#include <complex.h>
#include <fftw3.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include "21CMMC.h"

void ComputeInitialConditions(struct UserParams user_params, struct CosmoParams cosmo_params, struct InitialConditions boxes) {

    /*
     Generates the initial conditions: gaussian random density field (DIM^3) as well as the equal or lower resolution velocity fields, and smoothed density field (HII_DIM^3).
     See INIT_PARAMS.H and ANAL_PARAMS.H to set the appropriate parameters.
     Output is written to ../Boxes

     Author: Andrei Mesinger
     Date: 9/29/06
     */

    fftwf_plan plan;

    unsigned long long ct;
    int n_x, n_y, n_z, i, j, k, ii;
    float k_x, k_y, k_z, k_mag, p, a, b, k_sq;
    double pixel_deltax;

    float f_pixel_factor;

    gsl_rng * r;

    /************  INITIALIZATION **********************/

    // Removed all references to threads as 21CMMC is always a single core implementation

    printf("%d\n",cosmo_params.RANDOM_SEED);
    // seed the random number generators
//    r = gsl_rng_alloc(gsl_rng_mt19937);
//    gsl_rng_set(r, cosmo_params.RANDOM_SEED);

/*
    // allocate array for the k-space and real-space boxes
    HIRES_box = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*user_params.KSPACE_NUM_PIXELS);
    HIRES_box_saved = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*user_params.KSPACE_NUM_PIXELS);

    // now allocate memory for the lower-resolution box
    // use HII_DIM from ANAL_PARAMS
    LOWRES_density = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
    LOWRES_vx = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
    LOWRES_vy= (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
    LOWRES_vz = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);

    if(SECOND_ORDER_LPT_CORRECTIONS){
        LOWRES_vx_2LPT = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
        LOWRES_vy_2LPT = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
        LOWRES_vz_2LPT = (float *) malloc(sizeof(float)*HII_TOT_NUM_PIXELS);
    }

    // find factor of HII pixel size / deltax pixel size
    f_pixel_factor = DIM/(float)HII_DIM;
*/

    /************  END INITIALIZATION ******************/

    /************ CREATE K-SPACE GAUSSIAN RANDOM FIELD ***********/
/*
    for (n_x=0; n_x<DIM; n_x++){
        // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
        if (n_x>MIDDLE)
            k_x =(n_x-DIM) * DELTA_K;  // wrap around for FFT convention
        else
            k_x = n_x * DELTA_K;

        for (n_y=0; n_y<DIM; n_y++){
            // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
            if (n_y>MIDDLE)
                k_y =(n_y-DIM) * DELTA_K;
            else
                k_y = n_y * DELTA_K;

            // since physical space field is real, only half contains independent modes
            for (n_z=0; n_z<=MIDDLE; n_z++){
                // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
                k_z = n_z * DELTA_K;

                // now get the power spectrum; remember, only the magnitude of k counts (due to issotropy)
                // this could be used to speed-up later maybe
                k_mag = sqrt(k_x*k_x + k_y*k_y + k_z*k_z);
                p = power_in_k(k_mag);

                // ok, now we can draw the values of the real and imaginary part
                // of our k entry from a Gaussian distribution
                a = gsl_ran_ugaussian(r);
                b = gsl_ran_ugaussian(r);
                HIRES_box[C_INDEX(n_x, n_y, n_z)] = sqrt(VOLUME*p/2.0) * (a + b*I);
            }
        }
    }
*/
    /*****  Adjust the complex conjugate relations for a real array  *****/

//    adj_complex_conj(HIRES_box);

    /*** Let's also create a lower-resolution version of the density field  ***/

/*    memcpy(HIRES_box_saved, HIRES_box, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

    if (DIM != HII_DIM)
        filter(HIRES_box, 0, L_FACTOR*BOX_LEN/(HII_DIM+0.0));
    // FFT back to real space
    plan = fftwf_plan_dft_c2r_3d(DIM, DIM, DIM, (fftwf_complex *)HIRES_box, (float *)HIRES_box, FFTW_ESTIMATE);
    fftwf_execute(plan);
    // now sample the filtered box
    for (i=0; i<HII_DIM; i++){
        for (j=0; j<HII_DIM; j++){
            for (k=0; k<HII_DIM; k++){
                LOWRES_density[HII_R_INDEX(i,j,k)] =
                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                   (unsigned long long)(j*f_pixel_factor+0.5),
                                                   (unsigned long long)(k*f_pixel_factor+0.5)))/VOLUME;
            }
        }
    }
*/

    /******* PERFORM INVERSE FOURIER TRANSFORM *****************/
    // add the 1/VOLUME factor when converting from k space to real space
/*
    memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

    for (ct=0; ct<KSPACE_NUM_PIXELS; ct++){
        HIRES_box[ct] /= VOLUME;
    }
    plan = fftwf_plan_dft_c2r_3d(DIM, DIM, DIM, (fftwf_complex *)HIRES_box, (float *)HIRES_box, FFTW_ESTIMATE);
    fftwf_execute(plan);
    fftwf_destroy_plan(plan);
    fftwf_cleanup();

    for (i=0; i<DIM; i++){
        for (j=0; j<DIM; j++){
            for (k=0; k<DIM; k++){
                *((float *)HIRES_density + R_FFT_INDEX(i,j,k)) = *((float *)HIRES_box + R_FFT_INDEX(i,j,k));
            }
        }
    }

    for(ii=0;ii<3;ii++) {

        memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
        // Now let's set the velocity field/dD/dt (in comoving Mpc)

        for (n_x=0; n_x<DIM; n_x++){
            if (n_x>MIDDLE)
                k_x =(n_x-DIM) * DELTA_K;  // wrap around for FFT convention
            else
                k_x = n_x * DELTA_K;

            for (n_y=0; n_y<DIM; n_y++){
                if (n_y>MIDDLE)
                    k_y =(n_y-DIM) * DELTA_K;
                else
                    k_y = n_y * DELTA_K;

                for (n_z=0; n_z<=MIDDLE; n_z++){
                    k_z = n_z * DELTA_K;

                    k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                    // now set the velocities
                    if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                        HIRES_box[0] = 0;
                    }
                    else{
                        if(ii==0) {
                            HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_x*I/k_sq/VOLUME;
                        }
                        if(ii==1) {
                            HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_y*I/k_sq/VOLUME;
                        }
                        if(ii==2) {
                            HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_z*I/k_sq/VOLUME;
                        }
                        // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                    }
                }
            }
        }

        if (DIM != HII_DIM)
            filter(HIRES_box, 0, L_FACTOR*BOX_LEN/(HII_DIM+0.0));

        plan = fftwf_plan_dft_c2r_3d(DIM, DIM, DIM, (fftwf_complex *)HIRES_box, (float *)HIRES_box, FFTW_ESTIMATE);
        fftwf_execute(plan);
        // now sample to lower res
        // now sample the filtered box
        for (i=0; i<HII_DIM; i++){
            for (j=0; j<HII_DIM; j++){
                for (k=0; k<HII_DIM; k++){
                    if(ii==0) {
                        LOWRES_vx[HII_R_INDEX(i,j,k)] =
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                           (unsigned long long)(j*f_pixel_factor+0.5),
                                                           (unsigned long long)(k*f_pixel_factor+0.5)));
                    }
                    if(ii==1) {
                        LOWRES_vy[HII_R_INDEX(i,j,k)] =
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                           (unsigned long long)(j*f_pixel_factor+0.5),
                                                           (unsigned long long)(k*f_pixel_factor+0.5)));
                    }
                    if(ii==2) {
                        LOWRES_vz[HII_R_INDEX(i,j,k)] =
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                           (unsigned long long)(j*f_pixel_factor+0.5),
                                                           (unsigned long long)(k*f_pixel_factor+0.5)));
                    }
                }
            }
        }
    }
    // write out file
*/
    /* *************************************************** *
     *              BEGIN 2LPT PART                        *
     * *************************************************** */

    // Generation of the second order Lagrangian perturbation theory (2LPT) corrections to the ZA
    // reference: Scoccimarro R., 1998, MNRAS, 299, 1097-1118 Appendix D
/*
    // Parameter set in ANAL_PARAMS.H
    if(SECOND_ORDER_LPT_CORRECTIONS){
        // use six supplementary boxes to store the gradients of phi_1 (eq. D13b)
        // Allocating the boxes
#define PHI_INDEX(i, j) ((int) ((i) - (j)) + 3*((j)) - ((int)(j))/2  )
        // ij -> INDEX
        // 00 -> 0
        // 11 -> 3
        // 22 -> 5
        // 10 -> 1
        // 20 -> 2
        // 21 -> 4

        fftwf_complex *phi_1[6];

        for(i = 0; i < 3; ++i){
            for(j = 0; j <= i; ++j){
                phi_1[PHI_INDEX(i, j)] = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
            }
        }

        for(i = 0; i < 3; ++i){
            for(j = 0; j <= i; ++j){

                // read in the box
                memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

                // generate the phi_1 boxes in Fourier transform
                for (n_x=0; n_x<DIM; n_x++){
                    if (n_x>MIDDLE)
                        k_x =(n_x-DIM) * DELTA_K;  // wrap around for FFT convention
                    else
                        k_x = n_x * DELTA_K;

                    for (n_y=0; n_y<DIM; n_y++){
                        if (n_y>MIDDLE)
                            k_y =(n_y-DIM) * DELTA_K;
                        else
                            k_y = n_y * DELTA_K;

                        for (n_z=0; n_z<=MIDDLE; n_z++){
                            k_z = n_z * DELTA_K;

                            k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                            float k[] = {k_x, k_y, k_z};
                            // now set the velocities
                            if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                                phi_1[PHI_INDEX(i, j)][0] = 0;
                            }
                            else{
                                phi_1[PHI_INDEX(i, j)][C_INDEX(n_x,n_y,n_z)] = -k[i]*k[j]*HIRES_box[C_INDEX(n_x, n_y, n_z)]/k_sq/VOLUME;
                                // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                            }
                        }
                    }
                }
                // Now we can generate the real phi_1[i,j]
                plan = fftwf_plan_dft_c2r_3d(DIM, DIM, DIM, (fftwf_complex *)phi_1[PHI_INDEX(i, j)], (float *)phi_1[PHI_INDEX(i, j)], FFTW_ESTIMATE);
                fftwf_execute(plan);
            }
        }

        // Then we will have the laplacian of phi_2 (eq. D13b)
        // After that we have to return in Fourier space and generate the Fourier transform of phi_2
        int m, l;
        for (i=0; i<DIM; i++){
            for (j=0; j<DIM; j++){
                for (k=0; k<DIM; k++){
                    *( (float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i), (unsigned long long)(j), (unsigned long long)(k) )) = 0.0;
                    for(m = 0; m < 3; ++m){
                        for(l = m+1; l < 3; ++l){
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k)) ) += ( *((float *)(phi_1[PHI_INDEX(l, l)]) + R_FFT_INDEX((unsigned long long) (i),(unsigned long long) (j),(unsigned long long) (k)))  ) * (  *((float *)(phi_1[PHI_INDEX(m, m)]) + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k)))  );
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k)) ) -= ( *((float *)(phi_1[PHI_INDEX(l, m)]) + R_FFT_INDEX((unsigned long long)(i),(unsigned long long) (j),(unsigned long long)(k) ) )  ) * (  *((float *)(phi_1[PHI_INDEX(l, m)]) + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k) ))  );
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k)) ) /= TOT_NUM_PIXELS;
                        }
                    }
                }
            }
        }

        plan = fftwf_plan_dft_r2c_3d(DIM, DIM, DIM, (float *)HIRES_box, (fftwf_complex *)HIRES_box, FFTW_ESTIMATE);
        fftwf_execute(plan);

        memcpy(HIRES_box_saved, HIRES_box, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

        // Now we can store the content of box in a back-up file
        // Then we can generate the gradients of phi_2 (eq. D13b and D9)
*/
        /***** Write out back-up k-box RHS eq. D13b *****/
/*
        // For each component, we generate the velocity field (same as the ZA part)

        // Now let's set the velocity field/dD/dt (in comoving Mpc)

        // read in the box
        // TODO correct free of phi_1

        for(ii=0;ii<3;ii++) {

            if(ii>0) {
                memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
            }
            // set velocities/dD/dt
            for (n_x=0; n_x<DIM; n_x++){
                if (n_x>MIDDLE)
                    k_x =(n_x-DIM) * DELTA_K;  // wrap around for FFT convention
                else
                    k_x = n_x * DELTA_K;

                for (n_y=0; n_y<DIM; n_y++){
                    if (n_y>MIDDLE)
                        k_y =(n_y-DIM) * DELTA_K;
                    else
                        k_y = n_y * DELTA_K;

                    for (n_z=0; n_z<=MIDDLE; n_z++){
                        k_z = n_z * DELTA_K;

                        k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                        // now set the velocities
                        if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                            HIRES_box[0] = 0;
                        }
                        else{
                            if(ii==0) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_x*I/k_sq;
                            }
                            if(ii==1) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_y*I/k_sq;
                            }
                            if(ii==2) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_z*I/k_sq;
                            }
                            // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                        }
                    }
                }
            }

            if (DIM != HII_DIM)
                filter(HIRES_box, 0, L_FACTOR*BOX_LEN/(HII_DIM+0.0));

            plan = fftwf_plan_dft_c2r_3d(DIM, DIM, DIM, (fftwf_complex *)HIRES_box, (float *)HIRES_box, FFTW_ESTIMATE);
            fftwf_execute(plan);
            // now sample to lower res
            // now sample the filtered box
            for (i=0; i<HII_DIM; i++){
                for (j=0; j<HII_DIM; j++){
                    for (k=0; k<HII_DIM; k++){
                        if(ii==0) {
                            LOWRES_vx_2LPT[HII_R_INDEX(i,j,k)] =
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                               (unsigned long long)(j*f_pixel_factor+0.5),
                                                               (unsigned long long)(k*f_pixel_factor+0.5)));
                        }
                        if(ii==1) {
                            LOWRES_vy_2LPT[HII_R_INDEX(i,j,k)] =
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                               (unsigned long long)(j*f_pixel_factor+0.5),
                                                               (unsigned long long)(k*f_pixel_factor+0.5)));
                        }
                        if(ii==2) {
                            LOWRES_vz_2LPT[HII_R_INDEX(i,j,k)] =
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                               (unsigned long long)(j*f_pixel_factor+0.5),
                                                               (unsigned long long)(k*f_pixel_factor+0.5)));
                        }
                    }
                }
            }
        }
        // deallocate the supplementary boxes
        for(i = 0; i < 3; ++i){
            for(j = 0; j <= i; ++j){
                fftwf_free(phi_1[PHI_INDEX(i,j)]);
            }
        }
    }
*/
    /* *********************************************** *
     *               END 2LPT PART                     *
     * *********************************************** */

    // deallocate
//    fftwf_free(HIRES_box);
//    fftwf_free(HIRES_box_saved);
}

/*****  Adjust the complex conjugate relations for a real array  *****/
/*
void adj_complex_conj(){
    int i, j, k;

    // corners
    HIRES_box[C_INDEX(0,0,0)] = 0;
    HIRES_box[C_INDEX(0,0,MIDDLE)] = crealf(HIRES_box[C_INDEX(0,0,MIDDLE)]);
    HIRES_box[C_INDEX(0,MIDDLE,0)] = crealf(HIRES_box[C_INDEX(0,MIDDLE,0)]);
    HIRES_box[C_INDEX(0,MIDDLE,MIDDLE)] = crealf(HIRES_box[C_INDEX(0,MIDDLE,MIDDLE)]);
    HIRES_box[C_INDEX(MIDDLE,0,0)] = crealf(HIRES_box[C_INDEX(MIDDLE,0,0)]);
    HIRES_box[C_INDEX(MIDDLE,0,MIDDLE)] = crealf(HIRES_box[C_INDEX(MIDDLE,0,MIDDLE)]);
    HIRES_box[C_INDEX(MIDDLE,MIDDLE,0)] = crealf(HIRES_box[C_INDEX(MIDDLE,MIDDLE,0)]);
    HIRES_box[C_INDEX(MIDDLE,MIDDLE,MIDDLE)] = crealf(HIRES_box[C_INDEX(MIDDLE,MIDDLE,MIDDLE)]);

    // do entire i except corners
    for (i=1; i<MIDDLE; i++){
        // just j corners
        for (j=0; j<=MIDDLE; j+=MIDDLE){
            for (k=0; k<=MIDDLE; k+=MIDDLE){
                HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX(DIM-i,j,k)]);
            }
        }

        // all of j
        for (j=1; j<MIDDLE; j++){
            for (k=0; k<=MIDDLE; k+=MIDDLE){
                HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX(DIM-i,DIM-j,k)]);
                HIRES_box[C_INDEX(i,DIM-j,k)] = conjf(HIRES_box[C_INDEX(DIM-i,j,k)]);
            }
        }
    } // end loop over i

    // now the i corners
    for (i=0; i<=MIDDLE; i+=MIDDLE){
        for (j=1; j<MIDDLE; j++){
            for (k=0; k<=MIDDLE; k+=MIDDLE){
                HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX(i,DIM-j,k)]);
            }
        }
    } // end loop over remaining j

}
*/


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include "configuration.h"
#include "utility/geometry.h"
#include "utility/SDL++.h"
#include <SDL.h>
#include <gsl/gsl>
#include <mutex>
#include <vector>
#include <memory>
class File_Browser;
class File_Entry;
class Metadata_Display;
class Level_Meter;
class Modal_Box;
class Player;
struct Player_State;
class Main_Layout;
struct Midi_Output;

class Application
{
public:
    static const Point size_;

    Application();
    ~Application();

    SDL_Window *init_window();
    SDL_Renderer *init_renderer();
    void exec();

    void set_scale_factor(SDL_Window *win, unsigned sf);
    void paint(SDL_Renderer *rr, int paint);
    void paint_scaled(SDL_Renderer *rr, int paint, unsigned scale);
    void paint_cached_background(SDL_Renderer *rr);
    bool handle_key_pressed(const SDL_KeyboardEvent &event);
    bool handle_key_released(const SDL_KeyboardEvent &event);
    bool handle_mouse_pressed(const SDL_MouseButtonEvent &event);
    bool handle_mouse_released(const SDL_MouseButtonEvent &event);
    bool handle_text_input(const SDL_TextInputEvent &event);

    void play_file(const std::string &dir, const File_Entry *entries, size_t index, size_t count);
    void play_random(const std::string &dir, const File_Entry &entry);
    void set_current_path(const std::string &path);
    static bool filter_file_name(const std::string &name);
    static bool filter_file_entry(const File_Entry &ent);

    void request_update();
    void update_modals();
    void open_help_dialog();
    void open_fx_dialog();
    void choose_midi_output(bool ask, gsl::cstring_span choice);
    void choose_synth(bool ask, gsl::cstring_span choice);
    void get_midi_outputs(std::vector<Midi_Output> &outputs);

    void choose_theme(gsl::cstring_span choice);
    void get_themes(std::vector<std::string> &themes);

    void load_theme(gsl::cstring_span theme);
    void load_default_theme();
    void load_theme_configuration(const CSimpleIniA &ini);
    void get_fx_parameters(const CSimpleIniA &ini, int *values) const;

    void engage_shutdown();
    void engage_shutdown_if_esc_key();
    void advance_shutdown();
    bool should_quit() const;

private:
    std::unique_ptr<CSimpleIniA> initialize_config();

private:
    void receive_state_in_other_thread(const Player_State &ps);

private:
    SDLpp_Window_u window_;
    SDLpp_Renderer_u renderer_;

    SDL_TimerID update_timer_ = 0;

    std::unique_ptr<Main_Layout> layout_;
    std::vector<std::unique_ptr<Modal_Box>> modal_;
    SDLpp_Texture_u cached_background_;

    SDLpp_Surface_u logo_image_;
    SDLpp_Surface_u wallpaper_image_;

    bool fadeout_engaged_ = false;
    int fadeout_time_ = 0;

    uint32_t esc_key_timer_ = 0;

    unsigned scale_factor_ = 1;
    std::unique_ptr<File_Browser> file_browser_;
    std::unique_ptr<Metadata_Display> metadata_display_;
    std::unique_ptr<Level_Meter> level_meter_[10];
    std::unique_ptr<Player> player_;

    std::string last_midi_output_choice_;
    std::string last_synth_choice_;
    std::string last_theme_choice_;

    std::unique_ptr<Player_State> ps_;
    std::mutex ps_mutex_;

    enum Info_Mode {
        Info_File,
        Info_Metadata,
        Info_Mode_Count,
    };
    Info_Mode info_mode_ = Info_File;
};

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
#include "lah.h"

#ifdef HAVE_LAPACK /* Use a LAPACK package to do the heavy work */
/*#include <cblas.h>*/
#include <lapacke.h>
lah_Return lah_forwardSub(lah_mat *B, lah_mat const *L)
{
    if ( L == NULL ||  B == NULL || L->nR != B->nR
        || B->nR != L->nC)
    {
        return lahReturnParameterError;
    }

    /*
    matrix_layout = LAH_LAPACK_LAYOUT;
    if (matrix_layout == LAPACK_COL_MAJOR)
    {
        lda = L->nR;
        ldb = B->nR;
    }
    else
    {
        lda = L->nC;
        ldb = B->nC;
    }
    */

    return (!POTRS(LAH_LAPACK_LAYOUT, 'L', L->nR, B->nC,
                   L->data, LAH_LEADING_DIM(L), B->data, LAH_LEADING_DIM(B))) ? lahReturnOk : lahReturnExternError;
}

#else /* fallback to primitive linear algebra subroutines */
#include <stdlib.h>
#include <math.h>
lah_Return lah_forwardSub(lah_mat *B, const lah_mat  *L)
{
    lah_index i, j, k;
    lah_value *diag, *row, *Aij, *bi, *bj;

    if ( L == NULL ||  B == NULL || L->nR != B->nR
        || B->nR != L->nC)
    {
        return lahReturnParameterError;
    }

    /* loop over cols of x and b */
    for (k = 0; k < B->nC; k++)
    {
        /* initialize pointers for i-loop*/
        bi = B->data + k * B->incCol;
        diag = L->data;
        row = L->data;
        /* loop over rows of x*/
        for (i = 0; i < B->nR; i++)
        {
            /* initialize pointers for j-loop */
            bj = B->data + k * B->incCol;
            Aij = row;

            /* substitute up to (excluding) current x */
            for (j = 0; j < i; j++)
            {
                *bi -= *Aij * *bj;

                /* End of Loop: Increment Pointers */
                bj += B->incRow;
                Aij += L->incCol;
            }
            /* divide xi by diagonal element of L */
            if (fabs(*diag) > 0.00001)
                *bi /= *diag;
            else
                return lahReturnMathError;

            /* End of Loop: Increment Pointers */
            diag += L->incRow + L->incCol;
            bi += B->incRow;
            row += L->incRow;
        }
    }
    return lahReturnOk;
}

#endif /* endif primitive implementation */