text:
//This does CELL (~soma) stage of GRU (gated recurrent unit) model.
//This requires each neuron to have 3 input time series, X, Xr, Xz,
//where X is the usual input and Xr, Xz the inputs for the reset and update gates,
//stacked into X = [X; Xr; Xz] for dim=0, or X = [X Xr Xz] for dim=1.

//For dim=0, R[:,t] = sig{Xr[:,t] + Ur*Y[:,t-1]} \n";
//           Z[:,t] = sig{Xz[:,t] + Uz*Y[:,t-1]} \n";
//           H[:,t] = R[:,t].*Y[:,t-1] \n";
//           Y[:,t] = Z[:,t].*Y[:,t-1] + (1-Z[:,t]).*tanh{X[:,t] + U*H[:,t]} \n";
//with sizes X, Xr, Xz: N x T \n";
//           U, Ur, Uz: N x N \n";
//           Y        : N x T \n";
//
//For dim=1, R[t,:] = sig{Xr[t,:] + Y[t-1,:]*Ur} \n";
//           Z[t,:] = sig{Xz[t,:] + Y[t-1,:]*Uz} \n";
//           H[t,:] = R[t,:].*Y[t-1,:] \n";
//           Y[t,:] = Z[t,:].*Y[t-1,:] + (1-Z[t,:]).*tanh{X[t,:] + H[t,:]*U} \n";
//with sizes X, Xr, Xz: T x N \n";
//           U, Ur, Uz: N x N \n";
//           Y        : T x N \n";
//
//where sig is the logistic (sigmoid) nonlinearity = 1/(1+exp(-x)),
//R is the reset gate, Z is the update gate, H is an intermediate (hidden) vector,
//U, Ur, Uz are NxN matrices, and Y is the output.

//Note that, the neurons of a layer are independent only if U, Ur, Uz are diagonal matrices.
//This is only really a CELL (~soma) stage in that case.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int gru_s (float *Y, const float *X, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru_d (double *Y, const double *X, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru_inplace_s (float *X, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru_inplace_d (double *X, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru_s (float *Y, const float *X, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t N2 = 2*N, NT = N*T, NT2 = 2*N*T;
    size_t nT, tN, tN3;

    float *R, *Z, *H;
    if (!(R=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                //R[n] = 1.0f / (1.0f+expf(-X[N+n]));
                Z[n] = 1.0f / (1.0f+expf(-X[N2+n]));
                Y[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_scopy((int)N,&X[tN3+N],1,R,1); cblas_scopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    //R[n] = 1.0f / (1.0f+expf(-R[n]));
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    //H[n] = R[n] * Y[tN-N+n];
                    H[n] = Y[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[tN3],1,&Y[tN],1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0f-Z[n])*tanhf(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-X[NT2+nT]));
                Y[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,R,1); cblas_scopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0f-Z[n])*tanhf(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-X[NT2+nT]));
                Y[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,R,1); cblas_scopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0f-Z[n])*tanhf(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0f / (1.0f+expf(-X[N2+n]));
                Y[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_scopy((int)N,&X[tN3+N],1,R,1); cblas_scopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0f+expf(-R[n]));
                }
                cblas_scopy((int)N,&X[tN3],1,&Y[tN],1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0f-Z[n])*tanhf(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_d (double *Y, const double *X, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t N2 = 2*N, NT = N*T, NT2 = 2*N*T;
    size_t nT, tN, tN3;

    double *R, *Z, *H;
    if (!(R=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-X[N2+n]));
                Y[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_dcopy((int)N,&X[tN3+N],1,R,1); cblas_dcopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[tN3],1,&Y[tN],1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0-Z[n])*tanh(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-X[NT2+nT]));
                Y[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,R,1); cblas_dcopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0-Z[n])*tanh(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-X[NT2+nT]));
                Y[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,R,1); cblas_dcopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = Z[n]*Y[t-1+nT] + (1.0-Z[n])*tanh(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-X[N2+n]));
                Y[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_dcopy((int)N,&X[tN3+N],1,R,1); cblas_dcopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&Y[tN-N],1,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&Y[tN-N],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = Y[tN-N+n] / (1.0+exp(-R[n]));
                }
                cblas_dcopy((int)N,&X[tN3],1,&Y[tN],1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = Z[n]*Y[tN-N+n] + (1.0-Z[n])*tanh(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_inplace_s (float *X, const float *U, const float *Ur, const float *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t N2 = 2*N, N3 = 3*N, NT = N*T, NT2 = 2*N*T;
    size_t nT, tN, tN3;

    float *R, *Z, *H;
    if (!(R=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_inplace_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_inplace_s: problem with malloc. "); perror("malloc"); return 1; }

    //clock_gettime(CLOCK_REALTIME,&tic);
    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0f / (1.0f+expf(-X[N2+n]));
                X[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_scopy((int)N,&X[tN3+N],1,R,1); cblas_scopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN3-N3],1,o,R,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN3-N3],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[tN3-N3+n] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN3],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN3+n] = Z[n]*X[tN-N3+n] + (1.0f-Z[n])*tanhf(X[tN3+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-X[NT2+nT]));
                X[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,R,1); cblas_scopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0f-Z[n])*tanhf(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0f / (1.0f+expf(-X[NT2+nT]));
                X[nT] = (1.0f-Z[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,R,1); cblas_scopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0f-Z[n])*tanhf(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0f / (1.0f+expf(-X[N2+n]));
                X[n] = (1.0f-Z[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_scopy((int)N,&X[tN3+N],1,R,1); cblas_scopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN3-N3],1,o,R,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN3-N3],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0f / (1.0f+expf(-Z[n]));
                    H[n] = X[tN3-N3+n] / (1.0f+expf(-R[n]));
                }
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN3],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN3+n] = Z[n]*X[tN-N3+n] + (1.0f-Z[n])*tanhf(X[tN3+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_inplace_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_inplace_d (double *X, const double *U, const double *Ur, const double *Uz, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t N2 = 2*N, N3 = 3*N, NT = N*T, NT2 = 2*N*T;
    size_t nT, tN, tN3;

    double *R, *Z, *H;
    if (!(R=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(Z=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_inplace_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_inplace_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-X[N2+n]));
                X[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_dcopy((int)N,&X[tN3+N],1,R,1); cblas_dcopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN3-N3],1,o,R,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN3-N3],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[tN3-N3+n] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN3],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN3+n] = Z[n]*X[tN-N3+n] + (1.0-Z[n])*tanh(X[tN3+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-X[NT2+nT]));
                X[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,R,1); cblas_dcopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uz,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0-Z[n])*tanh(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                Z[n] = 1.0 / (1.0+exp(-X[NT2+nT]));
                X[nT] = (1.0-Z[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,R,1); cblas_dcopy((int)N,&X[NT2+t],(int)T,Z,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[t-1],(int)T,o,R,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[t-1],(int)T,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[t-1+n*T] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = Z[n]*X[t-1+nT] + (1.0-Z[n])*tanh(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                Z[n] = 1.0 / (1.0+exp(-X[N2+n]));
                X[n] = (1.0-Z[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN3 = 3*tN;
                cblas_dcopy((int)N,&X[tN3+N],1,R,1); cblas_dcopy((int)N,&X[tN3+N2],1,Z,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Ur,(int)N,&X[tN3-N3],1,o,R,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uz,(int)N,&X[tN3-N3],1,o,Z,1);
                for (size_t n=0u; n<N; ++n)
                {
                    Z[n] = 1.0 / (1.0+exp(-Z[n]));
                    H[n] = X[tN3-N3+n] / (1.0+exp(-R[n]));
                }
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN3],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN3+n] = Z[n]*X[tN-N3+n] + (1.0-Z[n])*tanh(X[tN3+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_inplace_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
/* randist/geometric.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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 <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Geometric distribution (bernoulli trial with probability p)

   prob(k) =  p (1 - p)^(k-1) for n = 1, 2, 3, ...

   It gives the distribution of "waiting times" for an event that
   occurs with probability p. */

unsigned int
gsl_ran_geometric (const gsl_rng * r, const double p)
{
  double u = gsl_rng_uniform_pos (r);

  unsigned int k;

  if (p == 1)
    {
      k = 1;
    }
  else
    {
      k = log (u) / log (1 - p) + 1;
    }

  return k;
}

double
gsl_ran_geometric_pdf (const unsigned int k, const double p)
{
  if (k == 0)
    {
      return 0 ;
    }
  else if (k == 1)
    {
      return p ;
    }
  else
    {
      double P = p * pow (1 - p, k - 1.0);
      return P;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl/gsl_vector.h>

double vs_gsl_Gradient_ForwardDiff (
    const gsl_vector * x, /**< [in] Point at which the gradient is to be evaluated */
    void * data,          /**< [in] Optional parameters passed directly to the function func_f */
    double (*func_f)(const gsl_vector * x, void *data), /**< [in] User-supplied routine that returns the value of the function at x */
    gsl_vector * J,       /**< [out] Gradient vector (same length as x) */
    double dh             /**< [in] Increment in variable for numerical differentiation */
    );

ext: h


text: #include <lapacke.h>


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

#ifndef __GSL_STATISTICS_CHAR_H__
#define __GSL_STATISTICS_CHAR_H__

#include <stddef.h>
#include <gsl/gsl_types.h>

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

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_char_mean (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_variance (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_sd (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_variance_with_fixed_mean (const char data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_char_sd_with_fixed_mean (const char data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_char_absdev (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_skew (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_kurtosis (const char data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_char_lag1_autocorrelation (const char data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_char_covariance (const char data1[], const size_t stride1,const char data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_char_variance_m (const char data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_char_sd_m (const char data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_char_absdev_m (const char data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_char_skew_m_sd (const char data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_char_kurtosis_m_sd (const char data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_char_lag1_autocorrelation_m (const char data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_char_covariance_m (const char data1[], const size_t stride1,const char data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

GSL_EXPORT double gsl_stats_char_pvariance (const char data1[], const size_t stride1, const size_t n1, const char data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_char_ttest (const char data1[], const size_t stride1, const size_t n1, const char data2[], const size_t stride2, const size_t n2);

GSL_EXPORT char gsl_stats_char_max (const char data[], const size_t stride, const size_t n);
GSL_EXPORT char gsl_stats_char_min (const char data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_char_minmax (char * min, char * max, const char data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_char_max_index (const char data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_char_min_index (const char data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_char_minmax_index (size_t * min_index, size_t * max_index, const char data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_char_median_from_sorted_data (const char sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_char_quantile_from_sorted_data (const char sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_CHAR_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/**
 * ex_gnuplot.h
 *
 * utility functions for plotting with gnuplot (v 4.6)
 *
 *  - minimizes boilerplate needed to use
 *  - popen -> process piping available, less temp files
 *  - idea based on myexamples/gnuplot (syntax) and gnuplot_i by N. Devillard (pipes), except a bit simpler
 *
 * Aaro Salosensaari 2016
 *
 * Original gnuplot_i was written by N. Devillard 1998 (version 2.10 2003) and is in public domain.
 */

#ifndef EX_GNUPLOT
#define EX_GNUPLOT

#include <stdarg.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include "ex_util.h"

#define GP_MAX_TMP_FILES 32

typedef struct {
    /* Pipe to gnuplot process */
    FILE *gnucmd ;

    /* Number of currently active plots */
    int nplots ;

    /* Pointer to table of names of temporary files */
    char *tmp_filename_tbl[GP_MAX_TMP_FILES] ;

    /* Number of temporary files */
    int ntmp ;
} gnuplot_ctrl ;

/**
 * Public interface function declarations
 */

gnuplot_ctrl *gnuplot_init(void);

void gnuplot_close(gnuplot_ctrl *handle);

/*
 * gnuplot_cmd
 *
 * This sends a string to an active gnuplot session, to be executed.
 * There is strictly no way to know if the command has been
 * successfully executed or not.
 * The command syntax is the same as printf.
 */
void gnuplot_cmd(gnuplot_ctrl *handle, char const *cmd, ...);

/*
 * ex_plot_xf can plot a variable number n of x, f(x) plots.
 *
 * n: number of plots
 * xlow:
 * xhigh:
 * xstep: (self-evident)
 * f1: function to plot
 * title1: title
 * style1: gnuplot style spef
 * (repeat fi, titlei, stylei until i=n)
 */
void ex_plot_xf(gnuplot_ctrl *handle, const int n, double xlow, double xhigh, double xstep, db_fx f1, const char *title1, const char *style1, ...);

/*
 * ex_plot_xys can plot a variable number n of x, y plots.
 *
 * n: number of plots
 * xdata: x points
 * ydata: y points
 * title1: title
 * style1: gnuplot style spef
 * (repeat xdatai, ydatai, titlei, stylei until i=n)
 */
void ex_plot_xys(gnuplot_ctrl *handle, const int n, gsl_vector *xdata1, gsl_vector *ydata1, const char *title1, const char *style1, ...);

void gnuplot_xy(gnuplot_ctrl *handle, gsl_vector *xdata, gsl_vector *ydata, const char *title1, const char *style1);

#endif /* ifndef EX_GNUPLOT */


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fitsio.h>
#include <gsl/gsl_statistics.h>

#ifdef SPECTRIM
#include "FBSpectrim.h"
#else
#include "fbus.h"
#endif

typedef struct _Box {
  double  x;
  double  y;
  double  fwhm;
  double  cenx;
  double  ceny;
  double  counts;
  double  background;
  double  noise;
  double  sigmaxythresh;
  double  sigmaxy;
  double  sigmafwhmthresh;
  double  sigmafwhm;
  int     r;
} Box;

typedef struct _Back {
  double x;
  double y;
  double background;
  double sigma;
  int r;
  int width;
} Back;

Box box[4];
Back back;
long nelements, naxes[2], fpixel;
int boxsize = 13;

double stardist(int i, int j) {
  return( sqrt( (box[i].cenx-box[j].cenx)*(box[i].cenx-box[j].cenx) +
		(box[i].ceny-box[j].ceny)*(box[i].ceny-box[j].ceny) ) );
}

double seeing(double var, double d, double r) {
  double lambda;
  double b, K, seeing;

  lambda = 0.65e-6;
  b = r/d;

  /* pixel scale of 0.78 "/pixel, convert to rad */
  var = var*pow(0.78/206265.0,2);

  K = 0.364*(1.0 - 0.532*pow(b, -1/3) - 0.024*pow(b, -7/3));

  seeing = 206265.0*0.98*pow(d/lambda, 0.2)*pow(var/K, 0.6);

  return seeing;
}

double old_seeing(double var, double d, double r) {
  double lambda;
  double r0;

  lambda = 0.55e-6;

  /* pixel scale of 0.78 "/pixel, convert to rad */
  var = var*pow(0.78/206265.0,2);

  r0 = pow(2.0*(lambda*lambda)*(
				( 0.1790*pow(d, (-1.0/3.0)) -
				  0.0968*pow(r, (-1.0/3.0)) )/var
				), 0.6);

  return 206265.0*0.98*lambda/r0;
}

int background(char *image) {

  int i, j, backpix;
  int low_y, up_y, low_x, up_x;
  double dist, sum, sumsq;

  backpix = 0;
  sum = 0.0;
  sumsq = 0.0;

  low_y = back.y - back.r - back.width;
  up_y  = back.y + back.r + back.width;
  low_x = back.x - back.r - back.width;
  up_x  = back.x + back.r + back.width;
  if (low_y < 0) {
    low_y = 0;
  }
  if (up_y >= 480) {
    up_y = 479;
  }
  if (low_x < 0) {
    low_x = 0;
  }
  if (up_x >= 640) {
    up_x = 639;
  }

  for (i=low_y; i<up_y; i++) {
    for (j=low_x; j<up_x; j++) {
      dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
      if (dist >= back.r && dist <= back.r+back.width) {
	sum += image[i*naxes[0]+j];
	backpix++;
      }
    }
  }

  back.background = sum/backpix;

  for (i=low_y; i<up_y; i++) {
    for (j=low_x; j<up_x; j++) {
      dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
      if (dist >= back.r && dist <= back.r+back.width) {
	sumsq += (image[i*naxes[0]+j]-back.background)*
	  (image[i*naxes[0]+j]-back.background);
      }
    }
  }

  back.sigma = sqrt(sumsq/backpix);

  return 1;

}

int centroid(char *image, int num) {

  int i, j;
  double  sum    = 0.0;
  double  sumx   = 0.0;
  double  sumxx  = 0.0;
  double  sumy   = 0.0;
  double  sumyy  = 0.0;
  double  val = 0.0;
  double  rmom;
  double  dist;
  int low_y, up_y, low_x, up_x;
  int sourcepix = 0;

  low_y = box[num].y - box[num].r;
  up_y  = box[num].y + box[num].r;
  low_x = box[num].x - box[num].r;
  up_x  = box[num].x + box[num].r;
  if (low_y < 0) {
    low_y = 0;
  }
  if (up_y >= 480) {
    up_y = 479;
  }
  if (low_x < 0) {
    low_x = 0;
  }
  if (up_x >= 640) {
    up_x = 639;
  }

  for (i=low_y; i<up_y; i++) {
    for (j=low_x; j<up_x; j++) {

      dist = sqrt(pow(box[num].x-j, 2) + pow(box[num].y-i, 2));
      if (dist <= box[num].r) {
	val = image[i*naxes[0]+j] - back.background;
	sum   += val;
	sumx  += val*j;
	sumxx += val*j*j;
	sumy  += val*i;
	sumyy += val*i*i;
	sourcepix++;
      }

    }
  }

  if ( sum <= 0.0 ) {
    box[num].sigmaxy = -1.0;
    box[num].sigmafwhm = -1.0;
    box[num].cenx = 0.0;
    box[num].ceny = 0.0;
    box[num].fwhm = -1.0;
  } else {
    rmom = ( sumxx - sumx * sumx / sum + sumyy - sumy * sumy / sum ) / sum;

    if ( rmom <= 0 ) {
      box[num].fwhm = -1.0;
    } else {
      box[num].fwhm = sqrt(rmom)  * 2.354 / sqrt(2.0);
    }

    box[num].counts = sum;
    box[num].cenx   = sumx / sum;
    box[num].ceny   = sumy / sum;
    box[num].x = box[num].cenx;
    box[num].y = box[num].ceny;
    box[num].sigmaxy= box[num].noise * sourcepix / box[num].counts / sqrt(6.0);
    box[num].sigmafwhm = box[num].noise * pow(sourcepix,1.5) / 10.
      / box[num].fwhm / box[num].counts
      * 2.354 * 2.354 / 2.0;
  }

  return 1;

}

int main() {

  char * buffer;
  fitsfile *fptr;
  int i, j, f, status, nimages, anynul, nboxes;
  char fitsfile[256];
  char *froot;
  FILE *init, *out;
  float xx = 0.0, yy = 0.0, xsum = 0.0, ysum = 0.0;
  double dist01[900], dist02[900], dist03[900];
  double dist12[900], dist13[900], dist23[900];
  double sig01[900], sig02[900], sig03[900];
  double sig12[900], sig13[900], sig23[900];
  double mean01, mean02, mean03, mean12, mean13, mean23;
  double var01, var02, var03, var12, var13, var23, rad;
  double avesig01, avesig02, avesig03, avesig12, avesig13, avesig23;
  double d = 0.05;
  double r_01, r_02, r_03, r_12, r_13, r_23, seeing_ave;
  r_01 = 2.0*d;
  r_02 = 2.0*d;
  r_03 = 3.0*d;
  r_12 = 3.0*d;
  r_13 = 2.0*d;
  r_23 = 2.0*d;

  out = fopen("/dev/shm/video/seeing.dat", "a");
  init = fopen("/dev/shm/video/init_cen_all", "r");
  i = 0;
  while (fscanf(init, "%f %f\n", &xx, &yy) != EOF) {
    box[i].x = xx;
    box[i].cenx = xx;
    box[i].y = yy;
    box[i].ceny = yy;
    box[i].r = boxsize/2.0;
    i++;
  }
  fclose(init);

  back.r = 40;
  back.width = 5;

  status = 0;
  anynul = 0;
  naxes[0] = 640;
  naxes[1] = 480;
  fpixel = 1;

  nelements = naxes[0]*naxes[1];

  /* init card */
  if (FB_Init() == RET_ERROR) {
    printf("FlashBus failed init\n");
    return(-1);
  }

  /* configure video capture */
  FB_SetVideoConfig(TYPE_COMPOSITE, STANDARD_NTSC, 1, 0);

  /* set up offscreen video capture */
  if (FB_VideoOffscreen(naxes[0], naxes[1], 8, TRUE) != RET_SUCCESS) {
    printf("Couldn't set offscreen capture\n");
    return(-1);
  }

  /* allocate the buffer */
  if (!(buffer = malloc(nelements*sizeof(char)))) {
    printf("Couldn't Allocate Image Buffer\n");
    exit(-1);
  }

  FB_VideoLive(TRUE, ALIGN_ANY);
  FB_WaitVSync(12);

  nimages = 900;
  froot = "/dev/shm/video/seeing.fits";

  for (f=0; f<nimages; f++) {
    FB_WaitVSync(1);
    FB_CopyVGARect(0,
		   0,
		   naxes[0],
		   naxes[1],
		   buffer,
		   -1,
		   COPYDIR_FROMOFFSCREEN);
    xsum = 0.0;
    ysum = 0.0;
    for (i=0; i<4; i++) {
      xsum += box[i].cenx;
      ysum += box[i].ceny;
    }
    back.x = xsum/4.0;
    back.y = ysum/4.0;
    background(buffer);

    //printf("file: %s  \t  background = %f noise = %f\n",
    //	   fitsfile, back.background, back.sigma);

    for (i=0; i<4; i++) {
      box[i].noise = back.sigma;
      //box[i].r = boxsize/2.0;
      centroid(buffer, i);
      //box[i].r = boxsize/3.0;
      centroid(buffer, i);
      //box[i].r = boxsize/4.0;
      centroid(buffer, i);

      /* printf("\t star = %d, counts = %f, \n\t x = %f, y = %f, sigmaxy = %f\n",
	     i,
	     box[i].counts,
	     box[i].cenx+1,
	     box[i].ceny+1,
	     box[i].sigmaxy);
      */
    }

    dist01[f] = stardist(0, 1);
    sig01[f] = box[0].sigmaxy*box[0].sigmaxy + box[1].sigmaxy*box[1].sigmaxy;

    dist02[f] = stardist(0, 2);
    sig02[f] = box[0].sigmaxy*box[0].sigmaxy + box[2].sigmaxy*box[2].sigmaxy;

    dist03[f] = stardist(0, 3);
    sig03[f] = box[0].sigmaxy*box[0].sigmaxy + box[3].sigmaxy*box[3].sigmaxy;

    dist12[f] = stardist(1, 2);
    sig12[f] = box[1].sigmaxy*box[1].sigmaxy + box[2].sigmaxy*box[2].sigmaxy;

    dist13[f] = stardist(1, 3);
    sig13[f] = box[1].sigmaxy*box[1].sigmaxy + box[3].sigmaxy*box[3].sigmaxy;

    dist23[f] = stardist(2, 3);
    sig23[f] = box[3].sigmaxy*box[3].sigmaxy + box[2].sigmaxy*box[2].sigmaxy;

  }

  sprintf(fitsfile, "!%s", froot);
  fits_create_file(&fptr, fitsfile, &status);
  fits_create_img(fptr, BYTE_IMG, 2, naxes, &status);
  fits_write_img(fptr, TBYTE, fpixel, nelements, buffer, &status);
  fits_close_file(fptr, &status);
  fits_report_error(stderr, status);

  mean01 = gsl_stats_mean(dist01, 1, nimages);
  avesig01 = gsl_stats_mean(sig01, 1, nimages);
  printf("mean01 = %f, avesig01 = %f\n", mean01, avesig01);

  mean02 = gsl_stats_mean(dist02, 1, nimages);
  avesig02 = gsl_stats_mean(sig02, 1, nimages);
  printf("mean02 = %f, avesig02 = %f\n", mean02, avesig02);

  mean03 = gsl_stats_mean(dist03, 1, nimages);
  avesig03 = gsl_stats_mean(sig03, 1, nimages);
  printf("mean03 = %f, avesig03 = %f\n", mean03, avesig03);

  mean12 = gsl_stats_mean(dist12, 1, nimages);
  avesig12 = gsl_stats_mean(sig12, 1, nimages);
  printf("mean12 = %f, avesig12 = %f\n", mean12, avesig12);

  mean13 = gsl_stats_mean(dist13, 1, nimages);
  avesig13 = gsl_stats_mean(sig13, 1, nimages);
  printf("mean13 = %f, avesig13 = %f\n", mean13, avesig13);

  mean23 = gsl_stats_mean(dist23, 1, nimages);
  avesig23 = gsl_stats_mean(sig23, 1, nimages);
  printf("mean23 = %f, avesig23 = %f\n", mean23, avesig23);

  printf("\n");

  var01 = gsl_stats_variance_m(dist01, 1, nimages, mean01);
  var01 = var01 - avesig01;
  printf("sigma01 = %f, seeing01 = %f\n", sqrt(var01), seeing(var01, d, r_01));

  var02 = gsl_stats_variance_m(dist02, 1, nimages, mean02);
  var02 = var02 - avesig02;
  printf("sigma02 = %f, seeing02 = %f\n", sqrt(var02), seeing(var02, d, r_02));

  var03 = gsl_stats_variance_m(dist03, 1, nimages, mean03);
  var03 = var03 - avesig03;
  printf("sigma03 = %f, seeing03 = %f\n", sqrt(var03), seeing(var03, d, r_03));

  var12 = gsl_stats_variance_m(dist12, 1, nimages, mean12);
  var12 = var12 - avesig12;
  printf("sigma12 = %f, seeing12 = %f\n", sqrt(var12), seeing(var12, d, r_12));

  var13 = gsl_stats_variance_m(dist13, 1, nimages, mean13);
  var13 = var13 - avesig13;
  printf("sigma13 = %f, seeing13 = %f\n", sqrt(var13), seeing(var13, d, r_13));

  var23 = gsl_stats_variance_m(dist23, 1, nimages, mean23);
  var23 = var23 - avesig23;
  printf("sigma23 = %f, seeing23 = %f\n\n", sqrt(var23), seeing(var23, d, r_23));

  seeing_ave = (seeing(var01, d, r_01) +
		seeing(var02, d, r_02) +
		seeing(var03, d, r_03) +
		seeing(var12, d, r_12) +
		seeing(var13, d, r_13) +
		seeing(var23, d, r_23))/6.0;

  printf("Ave. seeing = %4.2f\"\n", seeing_ave);

  fprintf(out, "%f %f %f %f %f %f %f\n",
	  seeing(var01, d, r_01),
	  seeing(var02, d, r_02),
	  seeing(var03, d, r_03),
	  seeing(var12, d, r_12),
	  seeing(var13, d, r_13),
	  seeing(var23, d, r_23),
	  seeing_ave);

  init = fopen("/dev/shm/video/init_cen_all", "w");
  for (i=0; i<4; i++) {
    fprintf(init, "%f %f\n", box[i].cenx, box[i].ceny);
  }
  fclose(init);

  return (status);
}

max_stars_repo_name: marissakotze/timDIMM

max_issues_repo_name: marissakotze/timDIMM

max_forks_repo_name: marissakotze/timDIMM


text:
#ifndef MESH_H
#define MESH_H
#define EIGEN_SUPERLU_SUPPORT

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <chrono>
#include <armadillo>
#include <petsc.h>
#include "../TPLs/yaml-cpp/include/yaml-cpp/yaml.h"
#include "../TPLs/eigen-git-mirror/Eigen/Eigen"
#include "../TPLs/eigen-git-mirror/Eigen/IterativeLinearSolvers"
#include "../TPLs/eigen-git-mirror/unsupported/Eigen/IterativeSolvers"
#include "../TPLs/eigen-git-mirror/Eigen/SuperLUSupport"

using namespace std;

class WriteData; // forward declaration

class qdCell
{

	public:
	qdCell();
	//number of ordinates on this quadrature level
  	int cFIndex=0,nFIndex=0,sFIndex=0,eFIndex=0,wFIndex=0;
        int nCIndex=0,sCIndex=0,eCIndex=0,wCIndex=0;

};

class quadLevel
{
	public:
	quadLevel(vector< vector<double> > myQuad,\
  		vector<double> myAlpha,\
  		vector<double> myTau,\
                int myStartIndex);
	//number of ordinates on this quadrature level
  	int nOrd;
	//quadrature set
  	vector< vector<double> > quad;
        //difference coefficients
  	vector<double> alpha;
	//factor for linear interpolation of half angles
        vector<double> tau;
	//index of each ordinate in angular flux matrices
	vector<int> ordIdx;

};

class Mesh
{
	public:
	Mesh(YAML::Node * myInput);
  	int n,nAngles,nR,nZ;
        int state = 1;
        double dz,dr,drCorner,dzCorner,Z,R,dt,T,totalWeight;
        bool verbose = false, petsc = false, verbose_keff_only = false;
  	vector< vector<double> > quadSet;
  	vector< vector<double> > alpha;
        vector< vector<double> > tau;
        vector< vector<double> > cellVol;
        vector< vector<double> > cellVSA;
        vector<double> dts;
        vector<double> ts;
        vector<bool> outputOnStep;
        arma::rowvec dzs,dzsCorner;
        arma::rowvec drs,drsCorner;
	arma::rowvec rEdge;
	arma::rowvec zEdge;
        arma::rowvec rCent;
        arma::rowvec zCent;
	arma::rowvec rCornerEdge;
	arma::rowvec zCornerEdge;
        arma::rowvec rCornerCent;
        arma::rowvec zCornerCent;
        arma::rowvec rVWCornerCent;
        vector<quadLevel> quadrature;
	vector<qdCell> qdCells;
        Eigen::MatrixXd volume;
        string outputDir = "mesh/";
        WriteData * output;

        // Recirculation loop parameters
        double dzCornerRecirc,recircZ;
        int nZrecirc;
	arma::rowvec zCornerEdgeRecirc,dzsCornerRecirc;
        arma::rowvec zCornerCentRecirc;

        // Functions
  	vector<int> getQDCellIndices(int iR, int iZ);
  	vector<double> getGeoParams(int iR, int iZ);
  	vector<double> getRecircGeoParams(int iR, int iZ);
  	void advanceOneTimeStep();
  	void calcQuadSet();
        void writeVars();
	void printQuadSet();
        void checkOptionalParams();

        private:
	vector<double> mu;
  	vector< vector<double> > ordinates;
  	void calcMu();
	void calcAlpha();
        void calcTau();
        void calcSpatialMesh();
        void calcQDCellIndices(int nCornersR,int nCornersZ);
        void addLevels();
 	void calcNumAnglesTotalWeight();
 	void calcTimeMesh();
  	void calcRecircMesh();
	int quad_index(int p,int q);
	int low_quad_index(int p,int q);
	int getQDCellIndex(int iR,int iZ);
	YAML::Node * input;
};

typedef Eigen::Array<bool,Eigen::Dynamic,1> VectorXb;

#endif

max_stars_repo_name: khurrumsaleem/QuasiMolto

max_issues_repo_name: khurrumsaleem/QuasiMolto

max_forks_repo_name: khurrumsaleem/QuasiMolto


text:
/* linalg/cod.c
 *
 * Copyright (C) 2016, 2017 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 <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

/*
 * This module contains routines for factoring an M-by-N matrix A as:
 *
 * A P = Q R Z^T
 *
 * known as the Complete Orthogonal Decomposition, where:
 *
 * P is a N-by-N permutation matrix
 * Q is M-by-M orthogonal
 * R has an r-by-r upper triangular block
 * Z is N-by-N orthogonal
 *
 * When A is full rank, Z = I and this becomes the QR decomposition
 * with column pivoting. When A is rank deficient, then
 *
 * R = [ R11 0 ] where R11 is r-by-r and r = rank(A)
 *     [  0  0 ]
 */

static int cod_RZ(gsl_matrix * A, gsl_vector * tau);
static double cod_householder_transform(double *alpha, gsl_vector * v);
static int cod_householder_mh(const double tau, const gsl_vector * v,
                              gsl_matrix * A, gsl_vector * work);
static int cod_householder_hv(const double tau, const gsl_vector * v, gsl_vector * w);
static int cod_householder_Zvec(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
                                gsl_vector * v);
static int cod_trireg_solve(const gsl_matrix * R, const double lambda, const gsl_vector * b,
                            gsl_matrix * S, gsl_vector * x, gsl_vector * work);

int
gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
                        gsl_permutation * p, double tol, size_t * rank, gsl_vector * work)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau_Q->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau_Q must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (tau_Z->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau_Z must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
    }
  else if (work->size != N)
    {
      GSL_ERROR ("work size must be N", GSL_EBADLEN);
    }
  else
    {
      int status, signum;
      size_t r;

      /* decompose: A P = Q R */
      status = gsl_linalg_QRPT_decomp(A, tau_Q, p, &signum, work);
      if (status)
        return status;

      /* estimate rank of A */
      r = gsl_linalg_QRPT_rank(A, tol);

      if (r < N)
        {
          /*
           * matrix is rank-deficient, so that the R factor is
           *
           * R = [ R11 R12 ] =~ [ R11 R12 ]
           *     [  0  R22 ]    [  0   0  ]
           *
           * compute RZ decomposition of upper trapezoidal matrix
           * [ R11 R12 ] = [ R11~ 0 ] Z
           */
          gsl_matrix_view R_upper = gsl_matrix_submatrix(A, 0, 0, r, N);
          gsl_vector_view t = gsl_vector_subvector(tau_Z, 0, r);

          cod_RZ(&R_upper.matrix, &t.vector);
        }

      *rank = r;

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
                      gsl_permutation * p, size_t * rank, gsl_vector * work)
{
  return gsl_linalg_COD_decomp_e(A, tau_Q, tau_Z, p, -1.0, rank, work);
}

/*
gsl_linalg_COD_lssolve()
  Find the least squares solution to the overdetermined system

min ||b - A x||^2

for M >= N using the COD factorization A P = Q R Z

Inputs: QRZT     - matrix A, in COD compressed format, M-by-N
        tau_Q    - Householder scalars for Q, length min(M,N)
        tau_Z    - Householder scalars for Z, length min(M,N)
        perm     - permutation matrix
        rank     - rank of A
        b        - rhs vector, length M
        x        - (output) solution vector, length N
        residual - (output) residual vector, b - A x, length M
*/

int
gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
                        const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
                        gsl_vector * x, gsl_vector * residual)
{
  const size_t M = QRZT->size1;
  const size_t N = QRZT->size2;

  if (M < N)
    {
      GSL_ERROR ("QRZT matrix must have M>=N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (rank > GSL_MIN (M, N))
    {
      GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (M != residual->size)
    {
      GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R11 = gsl_matrix_const_submatrix (QRZT, 0, 0, rank, rank);
      gsl_vector_view QTb1 = gsl_vector_subvector(residual, 0, rank);
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, rank);

      gsl_vector_set_zero(x);

      /* compute residual = Q^T b = [ c1 ; c2 ] */
      gsl_vector_memcpy(residual, b);
      gsl_linalg_QR_QTvec (QRZT, tau_Q, residual);

      /* solve x1 := R11^{-1} (Q^T b)(1:r) */
      gsl_vector_memcpy(&(x1.vector), &(QTb1.vector));
      gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &(R11.matrix), &(x1.vector));

      /* compute Z ( R11^{-1} x1; 0 ) */
      cod_householder_Zvec(QRZT, tau_Z, rank, x);

      /* compute x = P Z^T ( R11^{-1} x1; 0 ) */
      gsl_permute_vector_inverse(perm, x);

      /* compute residual = b - A x = Q (Q^T b - R [ R11^{-1} x1; 0 ]) = Q [ 0 ; c2 ] */
      gsl_vector_set_zero(&(QTb1.vector));
      gsl_linalg_QR_Qvec(QRZT, tau_Q, residual);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_COD_lssolve2()
  Find the least squares solution to the Tikhonov regularized
system in standard form:

min ||b - A x||^2 + lambda^2 ||x||^2

for M >= N using the COD factorization A P = Q R Z

Inputs: lambda   - parameter
        QRZT     - matrix A, in COD compressed format, M-by-N
        tau_Q    - Householder scalars for Q, length min(M,N)
        tau_Z    - Householder scalars for Z, length min(M,N)
        perm     - permutation matrix
        rank     - rank of A
        b        - rhs vector, length M
        x        - (output) solution vector, length N
        residual - (output) residual vector, b - A x, length M
        S        - workspace, rank-by-rank
        work     - workspace, length rank
*/

int
gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
                         const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
                         gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work)
{
  const size_t M = QRZT->size1;
  const size_t N = QRZT->size2;

  if (M < N)
    {
      GSL_ERROR ("QRZT matrix must have M>=N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (rank > GSL_MIN (M, N))
    {
      GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (M != residual->size)
    {
      GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
    }
  else if (S->size1 != rank || S->size2 != rank)
    {
      GSL_ERROR ("S must be rank-by-rank", GSL_EBADLEN);
    }
  else if (work->size != rank)
    {
      GSL_ERROR ("work must be length rank", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R11 = gsl_matrix_const_submatrix (QRZT, 0, 0, rank, rank);
      gsl_vector_view c1 = gsl_vector_subvector(residual, 0, rank);
      gsl_vector_view y1 = gsl_vector_subvector(x, 0, rank);

      gsl_vector_set_zero(x);

      /* compute residual = Q^T b = [ c1 ; c2 ]*/
      gsl_vector_memcpy(residual, b);
      gsl_linalg_QR_QTvec (QRZT, tau_Q, residual);

      /* solve [ R11 ; lambda*I ] y1 = [ (Q^T b)(1:r) ; 0 ] */
      cod_trireg_solve(&(R11.matrix), lambda, &(c1.vector), S, &(y1.vector), work);

      /* save y1 for later residual calculation */
      gsl_vector_memcpy(work, &(y1.vector));

      /* compute Z [ y1; 0 ] */
      cod_householder_Zvec(QRZT, tau_Z, rank, x);

      /* compute x = P Z^T ( y1; 0 ) */
      gsl_permute_vector_inverse(perm, x);

      /* compute residual = b - A x = Q (Q^T b - [ R11 y1; 0 ]) = Q [ c1 - R11*y1 ; c2 ] */

      /* work = R11*y1 */
      gsl_blas_dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, &(R11.matrix), work);

      gsl_vector_sub(&(c1.vector), work);
      gsl_linalg_QR_Qvec(QRZT, tau_Q, residual);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_COD_unpack()
  Unpack encoded COD decomposition into the matrices Q,R,Z,P

Inputs: QRZT  - encoded COD decomposition
        tau_Q - Householder scalars for Q
        tau_Z - Householder scalars for Z
        rank  - rank of matrix (as determined from gsl_linalg_COD_decomp)
        Q     - (output) M-by-M matrix Q
        R     - (output) M-by-N matrix R
        Z     - (output) N-by-N matrix Z
*/

int
gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,
                      const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,
                      gsl_matrix * R, gsl_matrix * Z)
{
  const size_t M = QRZT->size1;
  const size_t N = QRZT->size2;

  if (tau_Q->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau_Q must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (tau_Z->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau_Z must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (rank > GSL_MIN (M, N))
    {
      GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
    }
  else if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q must by M-by-M", GSL_EBADLEN);
    }
  else if (R->size1 != M || R->size2 != N)
    {
      GSL_ERROR ("R must by M-by-N", GSL_EBADLEN);
    }
  else if (Z->size1 != N || Z->size2 != N)
    {
      GSL_ERROR ("Z must by N-by-N", GSL_EBADLEN);
    }
  else
    {
      size_t i;
      gsl_matrix_view R11 = gsl_matrix_submatrix(R, 0, 0, rank, rank);
      gsl_matrix_const_view QRZT11 = gsl_matrix_const_submatrix(QRZT, 0, 0, rank, rank);

      /* form Q matrix */

      gsl_matrix_set_identity(Q);

      for (i = GSL_MIN (M, N); i-- > 0;)
        {
          gsl_vector_const_view h = gsl_matrix_const_subcolumn (QRZT, i, i, M - i);
          gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
          double ti = gsl_vector_get (tau_Q, i);
          gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
        }

      /* form Z matrix */
      gsl_matrix_set_identity(Z);

      if (rank < N)
        {
          gsl_vector_view work = gsl_matrix_row(R, 0); /* temporary workspace, size N */

          /* multiply I by Z from the right */
          gsl_linalg_COD_matZ(QRZT, tau_Z, rank, Z, &work.vector);
        }

      /* copy rank-by-rank upper triangle of QRZT into R and zero the rest */
      gsl_matrix_set_zero(R);
      gsl_matrix_tricpy('U', 1, &R11.matrix, &QRZT11.matrix);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_COD_matZ
  Multiply an M-by-N matrix A on the right by Z (N-by-N)

Inputs: QRZT  - encoded COD matrix
        tau_Z - Householder scalars for Z
        rank  - matrix rank
        A     - on input, M-by-N matrix
                on output, A * Z
        work  - workspace of length M
*/

int
gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
                    gsl_matrix * A, gsl_vector * work)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau_Z->size != GSL_MIN (QRZT->size1, QRZT->size2))
    {
      GSL_ERROR("tau_Z must be GSL_MIN(M,N)", GSL_EBADLEN);
    }
  else if (QRZT->size2 != N)
    {
      GSL_ERROR("QRZT must have N columns", GSL_EBADLEN);
    }
  else if (work->size != M)
    {
      GSL_ERROR("workspace must be length M", GSL_EBADLEN);
    }
  else
    {
      /* if rank == N, then Z = I and there is nothing to do */
      if (rank < N)
        {
          size_t i;

          for (i = rank; i > 0 && i--; )
            {
              gsl_vector_const_view h = gsl_matrix_const_subrow (QRZT, i, rank, N - rank);
              gsl_matrix_view m = gsl_matrix_submatrix (A, 0, i, M, N - i);
              double ti = gsl_vector_get (tau_Z, i);
              cod_householder_mh (ti, &h.vector, &m.matrix, work);
            }
        }

      return GSL_SUCCESS;
    }
}

/*********************************************
 * INTERNAL ROUTINES                         *
 *********************************************/

/*
cod_RZ()
  Perform RZ decomposition of an upper trapezoidal matrix,

A = [ A11 A12 ] = [ R 0 ] Z

where A is M-by-N with N >= M, A11 is M-by-M upper triangular,
and A12 is M-by-(N-M). On output, Z is stored as Householder
reflectors in the A12 portion of A,

Z = Z(1) Z(2) ... Z(M)

Inputs: A   - M-by-N matrix with N >= M
              On input, upper trapezoidal matrix [ A11 A12 ]
              On output, A11 is overwritten by R (subdiagonal elements
              are not touched), and A12 is overwritten by Z in packed storage
        tau - (output) Householder scalars, size M
*/

static int
cod_RZ(gsl_matrix * A, gsl_vector * tau)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau->size != M)
    {
      GSL_ERROR("tau has wrong size", GSL_EBADLEN);
    }
  else if (N < M)
    {
      GSL_ERROR("N must be >= M", GSL_EINVAL);
    }
  else if (M == N)
    {
      /* quick return */
      gsl_vector_set_all(tau, 0.0);
      return GSL_SUCCESS;
    }
  else
    {
      size_t k;

      for (k = M; k > 0 && k--; )
        {
          double *alpha = gsl_matrix_ptr(A, k, k);
          gsl_vector_view z = gsl_matrix_subrow(A, k, M, N - M);
          double tauk;

          /* compute Householder reflection to zero [ A(k,k) A(k,M+1:N) ] */
          tauk = cod_householder_transform(alpha, &z.vector);
          gsl_vector_set(tau, k, tauk);

          if ((tauk != 0) && (k > 0))
            {
              gsl_vector_view w = gsl_vector_subvector(tau, 0, k);
              gsl_matrix_view B = gsl_matrix_submatrix(A, 0, k, k, N - k);

              cod_householder_mh(tauk, &z.vector, &B.matrix, &w.vector);
            }
        }

      return GSL_SUCCESS;
    }
}

static double
cod_householder_transform(double *alpha, gsl_vector * v)
{
  double beta, tau;
  double xnorm = gsl_blas_dnrm2(v);

  if (xnorm == 0)
    {
      return 0.0; /* tau = 0 */
    }

  beta = - (*alpha >= 0.0 ? +1.0 : -1.0) * gsl_hypot(*alpha, xnorm);
  tau = (beta - *alpha) / beta;

  {
    double s = (*alpha - beta);

    if (fabs(s) > GSL_DBL_MIN)
      {
        gsl_blas_dscal (1.0 / s, v);
      }
    else
      {
        gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
        gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
      }

    *alpha = beta;
  }

  return tau;
}

/*
cod_householder_hv
  Apply Householder reflection H = (I - tau*v*v') to vector v from the left,

w' = H * w

Inputs: tau  - Householder scalar
        v    - Householder vector, size M
        w    - on input, w vector, size M
               on output, H * w

Notes:
1) Based on LAPACK routine DLARZ
*/

static int
cod_householder_hv(const double tau, const gsl_vector * v, gsl_vector * w)
{
  if (tau == 0)
    {
      return GSL_SUCCESS; /* H = I */
    }
  else
    {
      const size_t M = w->size;
      const size_t L = v->size;
      double w0 = gsl_vector_get(w, 0);
      gsl_vector_view w1 = gsl_vector_subvector(w, M - L, L);
      double d1, d;

      /* d1 := v . w(M-L:M) */
      gsl_blas_ddot(v, &w1.vector, &d1);

      /* d := w(1) + v . w(M-L:M) */
      d = w0 + d1;

      /* w(1) = w(1) - tau * d */
      gsl_vector_set(w, 0, w0 - tau * d);

      /* w(M-L:M) = w(M-L:M) - tau * d * v */
      gsl_blas_daxpy(-tau * d, v, &w1.vector);

      return GSL_SUCCESS;
    }
}

/*
cod_householder_mh
  Apply Householder reflection H = (I - tau*v*v') to matrix A from the right

Inputs: tau  - Householder scalar
        v    - Householder vector, size N-M
        A    - matrix, size M-by-N
        work - workspace, size M

Notes:
1) Based on LAPACK routine DLARZ
*/

static int
cod_householder_mh(const double tau, const gsl_vector * v, gsl_matrix * A,
                   gsl_vector * work)
{
  if (tau == 0)
    {
      return GSL_SUCCESS; /* H = I */
    }
  else
    {
      const size_t M = A->size1;
      const size_t N = A->size2;
      const size_t L = v->size;
      gsl_vector_view A1 = gsl_matrix_subcolumn(A, 0, 0, M);
      gsl_matrix_view C = gsl_matrix_submatrix(A, 0, N - L, M, L);

      /* work(1:M) = A(1:M,1) */
      gsl_vector_memcpy(work, &A1.vector);

      /* work(1:M) = work(1:M) + A(1:M,M+1:N) * v(1:N-M) */
      gsl_blas_dgemv(CblasNoTrans, 1.0, &C.matrix, v, 1.0, work);

      /* A(1:M,1) = A(1:M,1) - tau * work(1:M) */
      gsl_blas_daxpy(-tau, work, &A1.vector);

      /* A(1:M,M+1:N) = A(1:M,M+1:N) - tau * work(1:M) * v(1:N-M)' */
      gsl_blas_dger(-tau, work, v, &C.matrix);

      return GSL_SUCCESS;
    }
}

/*
cod_householder_Zvec
  Multiply a vector by Z

Inputs: QRZT  - encoded COD matrix
        tau_Z - Householder scalars for Z
        rank  - matrix rank
        v     - on input, vector of length N
                on output, Z^T * v
*/

static int
cod_householder_Zvec(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
                     gsl_vector * v)
{
  const size_t M = QRZT->size1;
  const size_t N = QRZT->size2;

  if (tau_Z->size != GSL_MIN (M, N))
    {
      GSL_ERROR("tau_Z must be GSL_MIN(M,N)", GSL_EBADLEN);
    }
  else if (v->size != N)
    {
      GSL_ERROR("v must be length N", GSL_EBADLEN);
    }
  else
    {
      if (rank < N)
        {
          size_t i;

          for (i = 0; i < rank; ++i)
            {
              gsl_vector_const_view h = gsl_matrix_const_subrow (QRZT, i, rank, N - rank);
              gsl_vector_view w = gsl_vector_subvector (v, i, N - i);
              double ti = gsl_vector_get (tau_Z, i);
              cod_householder_hv(ti, &h.vector, &w.vector);
            }
        }

      return GSL_SUCCESS;
    }
}

/*
cod_trireg_solve()

  This function computes the solution to the least squares system

  [    R     ] x = [ b ]
  [ lambda*I ]     [ 0 ]

where R is an N-by-N upper triangular matrix, lambda is a scalar parameter,
and b is a vector of length N. This is done by computing the QR factorization

[    R     ] = W S^T
[ lambda*I ]

where S^T is upper triangular, and solving

S^T x = W^T [ b ]
            [ 0 ]

Inputs: R      - full rank upper triangular matrix; the diagonal
                 elements are modified but restored on output
        lambda - scalar parameter lambda
        b      - right hand side vector b
        S      - workspace, N-by-N
        x      - (output) least squares solution of the system
        work   - workspace of length N
*/

static int
cod_trireg_solve (const gsl_matrix * R, const double lambda, const gsl_vector * b,
                  gsl_matrix * S, gsl_vector * x, gsl_vector * work)
{
  const size_t N = R->size2;
  gsl_vector_const_view diag = gsl_matrix_const_diagonal(R);
  size_t i, j, k;

  if (lambda <= 0.0)
    {
      GSL_ERROR("lambda must be positive", GSL_EINVAL);
    }

  /* copy R and b to preserve input and initialise S; store diag(R) in work */
  gsl_matrix_transpose_tricpy('U', 0, S, R);
  gsl_vector_memcpy(work, &diag.vector);
  gsl_vector_memcpy(x, b);

  /* eliminate the diagonal matrix lambda*I using Givens rotations */

  for (j = 0; j < N; j++)
    {
      double bj = 0.0;

      gsl_matrix_set (S, j, j, lambda);

      for (k = j + 1; k < N; k++)
        {
          gsl_matrix_set (S, k, k, 0.0);
        }

      /* the transformations to eliminate the row of lambda*I modify only a
         single element of b beyond the first n, which is initially
         zero */

      for (k = j; k < N; k++)
        {
          /* determine a Givens rotation which eliminates the
             appropriate element in the current row of lambda*I */

          double sine, cosine;

          double xk = gsl_vector_get (x, k);
          double rkk = gsl_vector_get (work, k);
          double skk = gsl_matrix_get (S, k, k);

          if (skk == 0)
            {
              continue;
            }

          if (fabs (rkk) < fabs (skk))
            {
              double cotangent = rkk / skk;
              sine = 0.5 / sqrt (0.25 + 0.25 * cotangent * cotangent);
              cosine = sine * cotangent;
            }
          else
            {
              double tangent = skk / rkk;
              cosine = 0.5 / sqrt (0.25 + 0.25 * tangent * tangent);
              sine = cosine * tangent;
            }

          /* Compute the modified diagonal element of r and the
             modified element of [b,0] */

          {
            double new_rkk = cosine * rkk + sine * skk;
            double new_xk = cosine * xk + sine * bj;

            bj = -sine * xk + cosine * bj;

            gsl_vector_set(work, k, new_rkk);
            gsl_matrix_set(S, k, k, new_rkk);
            gsl_vector_set(x, k, new_xk);
          }

          /* Accumulate the transformation in the row of s */

          for (i = k + 1; i < N; i++)
            {
              double sik = gsl_matrix_get (S, i, k);
              double sii = gsl_matrix_get (S, i, i);

              double new_sik = cosine * sik + sine * sii;
              double new_sii = -sine * sik + cosine * sii;

              gsl_matrix_set(S, i, k, new_sik);
              gsl_matrix_set(S, i, i, new_sii);
            }
        }
    }

  /* solve: S^T x = rhs in place */
  gsl_blas_dtrsv(CblasLower, CblasTrans, CblasNonUnit, S, x);

  return GSL_SUCCESS;
}


text:
#ifndef PROB_FUNCTIONS_H
#define PROB_FUNCTIONS_H

#include <gsl/gsl_vector.h>
#include "data.h"

double my_f (const gsl_vector *x, void *params);
void my_df (const gsl_vector *x, void *params, gsl_vector *g);
void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *g);
void gradientQ (Dataset *data, double *dQdAlpha, double *dQdBeta);
double computeLikelihood (Dataset *data);
double computeQ (Dataset *data);
double logProbL (int l, int z, double alphaI, double betaJ);
double prob(double alpha, double beta);
void EStep (Dataset *data);
void MStep (Dataset *data);

#endif

max_stars_repo_name: ShreyaR/WorkerPoolSelection

max_issues_repo_name: ShreyaR/WorkerPoolSelection

max_forks_repo_name: ShreyaR/WorkerPoolSelection


text:
//This computes the IDFT (inverse discrete Fourier transformation) along dim of matrix X.
//This uses a CBLAS matrix multiplication by the IDFT matrix.

//The real-valued case is complicated (i.e., where only real part is output).
//This will only exactly invert with DFT if the DFT output all F=ndft/2+1 positive freqs.
//That is, use of this function to invert a specific DFT requires knowledge of how that DFT was run.

//Profile note: I tried splitting the matrix multiply into two halves,
//for positive and negative freqs, so that X would not need to be copied into Xc,
//but it was definitely slower.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int idft_cblas_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc);
int idft_cblas_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc);
int idft_cblas_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc);
int idft_cblas_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc);

int idft_cblas_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idft_cblas_s: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in idft_cblas_s: ndft must be >= Lx (length of vecs in X)\n"); return 1; }

    //Scaling
    const float s = sc ? 2.0f*sqrtf(0.5f*(float)ndft)/(float)ndft : 1.0f/(float)ndft;

    if (N==0u || ndft==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Init IDFT matrix
        const size_t NN = ndft * ndft;
        const float P2_N = (float)(2.0*M_PI/(double)ndft);
        float *IDFT;
        IDFT = (float *)aligned_alloc(sizeof(float),2u*NN*sizeof(float));
        if (!IDFT) { fprintf(stderr,"error in idft_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<ndft; ++l) { *IDFT++ = s; *IDFT++ = 0.0f; }
        for (size_t n=1u; n<ndft; ++n)
        {
            for (size_t l=0u; l<ndft; ++l)
            {
                *IDFT++ = s * cosf(P2_N*(float)n*(float)l);
                *IDFT++ = s * sinf(P2_N*(float)n*(float)l);
            }
        }
        IDFT -= 2u*NN;

        //Init complex matrix multiply
        const float z[2] = {0.0f,0.0f}, o[2] = {1.0f,0.0f};
        float *Xc, *Yc;
        Xc = (float *)aligned_alloc(sizeof(float),2u*ndft*sizeof(float));
        Yc = (float *)aligned_alloc(sizeof(float),2u*ndft*sizeof(float));
        if (!Xc) { fprintf(stderr,"error in idft_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        if (!Yc) { fprintf(stderr,"error in idft_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<2u*ndft; ++n, ++Xc) { *Xc = 0.0f; }
        Xc -= 2u*ndft;

        if (Lx==N)
        {
            //Copy into Xc
            for (size_t nf=Lx; nf>0u; --nf) { *Xc++ = *X++; *Xc++ = *X++; }
            if (ndft%2u==0u) { X -= 2; }
            for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2; *Xc++ = *X; *Xc++ = -*(X+1); }
            Xc -= 2u*ndft;

            //Matrix multiply
            cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

            //Output real part only
            cblas_scopy((int)ndft,Yc,2,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=ndft)
                {
                    //Copy into Xc
                    for (size_t nf=Lx; nf>0u; --nf) { *Xc++ = *X++; *Xc++ = *X++; }
                    if (ndft%2u==0u) { X -= 2; }
                    for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2; *Xc++ = *X; *Xc++ = -*(X+1); }
                    Xc -= 2u*ndft;
                    X += 2u*(ndft-ndft/2u-ndft%2u);

                    //Matrix multiply
                    cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

                    //Output real part only
                    cblas_scopy((int)ndft,Yc,2,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=B*(ndft-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2u*K*(ndft-ndft/2u-ndft%2u-Lx)+2u, ++Y)
                    {
                        //Copy into Xc
                        for (size_t nf=Lx; nf>0u; --nf, X+=2u*K) { *Xc++ = *X; *Xc++ = *(X+1); }
                        if (ndft%2u==0u) { X -= 2u*K; }
                        for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2u*K; *Xc++ = *X; *Xc++ = -*(X+1); }
                        Xc -= 2u*ndft;

                        //Matrix multiply
                        cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

                        //Output real part only
                        cblas_scopy((int)ndft,Yc,2,Y,(int)K);
                    }
                }
            }
        }
        free(IDFT); free(Xc); free(Yc);
    }

    return 0;
}

int idft_cblas_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idft_cblas_d: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in idft_cblas_d: ndft must be >= Lx (length of vecs in X)\n"); return 1; }

    //Scaling
    const double s = sc ? 2.0*sqrt(0.5*(double)ndft)/(double)ndft : 1.0/(double)ndft;

    if (N==0u || ndft==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Init IDFT matrix
        const size_t NN = ndft * ndft;
        const double P2_N = 2.0*M_PI/(double)ndft;
        double *IDFT;
        IDFT = (double *)aligned_alloc(sizeof(double),2u*NN*sizeof(double));
        if (!IDFT) { fprintf(stderr,"error in idft_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<ndft; ++l) { *IDFT++ = s; *IDFT++ = 0.0; }
        for (size_t n=1u; n<ndft; ++n)
        {
            for (size_t l=0u; l<ndft; ++l)
            {
                *IDFT++ = s * cos(P2_N*(double)n*(double)l);
                *IDFT++ = s * sin(P2_N*(double)n*(double)l);
            }
        }
        IDFT -= 2u*NN;

        //Init complex matrix multiply
        const double z[2] = {0.0f,0.0f}, o[2] = {1.0f,0.0f};
        double *Xc, *Yc;
        Xc = (double *)aligned_alloc(sizeof(double),2u*ndft*sizeof(double));
        Yc = (double *)aligned_alloc(sizeof(double),2u*ndft*sizeof(double));
        if (!Xc) { fprintf(stderr,"error in idft_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        if (!Yc) { fprintf(stderr,"error in idft_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<2u*ndft; ++n, ++Xc) { *Xc = 0.0; }
        Xc -= 2u*ndft;

        if (Lx==N)
        {
            //Copy into Xc
            for (size_t nf=Lx; nf>0u; --nf) { *Xc++ = *X++; *Xc++ = *X++; }
            if (ndft%2u==0u) { X -= 2; }
            for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2; *Xc++ = *X; *Xc++ = -*(X+1); }
            Xc -= 2u*ndft;

            //Matrix multiply
            cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

            //Output real part only
            cblas_dcopy((int)ndft,Yc,2,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=ndft)
                {
                    //Copy into Xc
                    for (size_t nf=Lx; nf>0u; --nf) { *Xc++ = *X++; *Xc++ = *X++; }
                    if (ndft%2u==0u) { X -= 2; }
                    for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2; *Xc++ = *X; *Xc++ = -*(X+1); }
                    Xc -= 2u*ndft;
                    X += 2u*(ndft-ndft/2u-ndft%2u);

                    //Matrix multiply
                    cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

                    //Output real part only
                    cblas_dcopy((int)ndft,Yc,2,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=B*(ndft-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2u*K*(ndft-ndft/2u-ndft%2u-Lx)+2u, ++Y)
                    {
                        //Copy into Xc
                        for (size_t nf=Lx; nf>0u; --nf, X+=2u*K) { *Xc++ = *X; *Xc++ = *(X+1); }
                        if (ndft%2u==0u) { X -= 2u*K; }
                        for (size_t nf=ndft-ndft/2u; nf>1u; --nf) { X-=2u*K; *Xc++ = *X; *Xc++ = -*(X+1); }
                        Xc -= 2u*ndft;

                        //Matrix multiply
                        cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)ndft,o,IDFT,(int)ndft,Xc,1,z,Yc,1);

                        //Output real part only
                        cblas_dcopy((int)ndft,Yc,2,Y,(int)K);
                    }
                }
            }
        }
        free(IDFT); free(Xc); free(Yc);
    }

    return 0;
}

int idft_cblas_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idft_cblas_c: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in idft_cblas_c: ndft must be >= Lx (length of vecs in X)\n"); return 1; }

    //Scaling
    const float s = sc ? 2.0f*sqrtf(0.5f*(float)ndft)/(float)ndft : 1.0f/(float)ndft;

    if (N==0u || ndft==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X * s; }
    }
    else
    {
        //Init IDFT matrix
        const size_t LN = Lx * ndft;
        const float P2_N = (float)(2.0*M_PI/(double)ndft);
        float *IDFT;
        IDFT = (float *)aligned_alloc(sizeof(float),2u*LN*sizeof(float));
        if (!IDFT) { fprintf(stderr,"error in idft_cblas_c: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l) { *IDFT++ = s; *IDFT++ = 0.0f; }
        for (size_t n=1u; n<ndft; ++n)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *IDFT++ = s * cosf(P2_N*(float)n*(float)l);
                *IDFT++ = s * sinf(P2_N*(float)n*(float)l);
            }
        }
        IDFT -= 2u*LN;

        //Init complex matrix multiply
        const float z[2] = {0.0f,0.0f}, o[2] = {1.0f,0.0f};

        if (Lx==N)
        {
            //Matrix multiply
            cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,1,z,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*ndft)
                {
                    //Matrix multiply
                    cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,1,z,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndft-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2, Y+=2)
                    {
                        //Matrix multiply
                        cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,(int)K,z,Y,(int)K);
                    }
                }
            }
        }
        free(IDFT);
    }

    return 0;
}

int idft_cblas_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idft_cblas_z: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in idft_cblas_z: ndft must be >= Lx (length of vecs in X)\n"); return 1; }

    //Scaling
    const double s = sc ? 2.0*sqrt(0.5*(double)ndft)/(double)ndft : 1.0/(double)ndft;

    if (N==0u || ndft==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X * s; }
    }
    else
    {
        //Init IDFT matrix
        const size_t LN = Lx * ndft;
        const double P2_N = 2.0*M_PI/(double)ndft;
        double *IDFT;
        IDFT = (double *)aligned_alloc(sizeof(double),2u*LN*sizeof(double));
        if (!IDFT) { fprintf(stderr,"error in idft_cblas_z: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l) { *IDFT++ = s; *IDFT++ = 0.0; }
        for (size_t n=1u; n<ndft; ++n)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *IDFT++ = s * cos(P2_N*(double)n*(double)l);
                *IDFT++ = s * sin(P2_N*(double)n*(double)l);
            }
        }
        IDFT -= 2u*LN;

        //Init complex matrix multiply
        const double z[2] = {0.0,0.0}, o[2] = {1.0,0.0};

        if (Lx==N)
        {
            //Matrix multiply
            cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,1,z,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*ndft)
                {
                    //Matrix multiply
                    cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,1,z,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndft-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2, Y+=2)
                    {
                        //Matrix multiply
                        cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)ndft,(int)Lx,o,IDFT,(int)Lx,X,(int)K,z,Y,(int)K);
                    }
                }
            }
        }
        free(IDFT);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#ifndef RENDERING_ENTITY_STEREO_VIEWERCENTRICENTITY_H
#define RENDERING_ENTITY_STEREO_VIEWERCENTRICENTITY_H

#include "stereoimageentity.h"

#include <gsl/gsl>

class QPaintEvent;
class QPainter;

namespace s3d {
struct ViewerContext;
}  // namespace s3d

class ViewerCentricEntity : public StereoImageEntity {
 public:
  ViewerCentricEntity();

  void init() override;
  void setHorizontalShift(float shift) override;
  void setAspectRatio(float ratio) override;
  void setTextureLeft(QOpenGLTexture* texture) override;
  void setTextureRight(QOpenGLTexture* texture) override;
  void setViewerContext(gsl::not_null<s3d::ViewerContext*> viewerContext);
  void setDepthRangeMeters(float min, float max);
  void setDisplayZoom(float displayZoom);

  void draw(QPaintDevice* paintDevice) override;

  QPoint worldToWidget(QPointF pos, int deviceWidth, int deviceHeight);
  float worldYToWidgetX(float y, int deviceWidth, int deviceHeight);
  float worldZToWidgetY(float z, int deviceHeight);
  float getPixelToWorldRatio(float deviceHeight);

 private:
  void drawScreen(QPainter* painter);
  void drawViewer(QPainter* painter);

 private:
  s3d::ViewerContext* m_viewerContext{};
  float m_displayZoom{1.0f};

  float m_minZ{-1.0f};
  float m_maxZ{3.0f};
};

#endif  // RENDERING_ENTITY_STEREO_VIEWERCENTRICENTITY_H


text:
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_statistics_double.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_stack_bwmorph.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_stack_draw.h"
#include "tz_voxel_graphics.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];

  sprintf(file_path, "../data/%s.tif", neuron_name);
  Stack *canvas = Read_Stack(file_path);

  Translate_Stack(canvas, COLOR, 1);
  Print_Stack_Info(canvas);

  sprintf(file_path, "../data/%s/grow_bundle.tif", neuron_name);
  Stack *bundle = Read_Stack(file_path);
  Stack *bundle_boundary = Stack_Perimeter(bundle, NULL, 4);
  Kill_Stack(bundle);

  sprintf(file_path, "../data/%s/bundle_seed.tif", neuron_name);
  Stack *seed = Read_Stack(file_path);

  Rgb_Color color;
  Set_Color(&color, 0, 255, 0);
  Stack_Label_Bwc(canvas, bundle_boundary, color);
  Set_Color(&color, 255, 0, 0);
  Stack_Label_Bwc(canvas, seed, color);

  sprintf(file_path, "../data/%s/bundle_label.tif", neuron_name);
  Write_Stack(file_path, canvas);

  Kill_Stack(canvas);
  Kill_Stack(bundle_boundary);
  Kill_Stack(seed);

  return 0;
}


text:
/* permutation/gsl_permutation.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_PERMUTATION_H__
#define __GSL_PERMUTATION_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>

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

__BEGIN_DECLS

struct gsl_permutation_struct
{
  size_t size;
  size_t *data;
};

typedef struct gsl_permutation_struct gsl_permutation;

gsl_permutation *gsl_permutation_alloc (const size_t n);
gsl_permutation *gsl_permutation_calloc (const size_t n);
void gsl_permutation_init (gsl_permutation * p);
void gsl_permutation_free (gsl_permutation * p);
int gsl_permutation_memcpy (gsl_permutation * dest, const gsl_permutation * src);

int gsl_permutation_fread (FILE * stream, gsl_permutation * p);
int gsl_permutation_fwrite (FILE * stream, const gsl_permutation * p);
int gsl_permutation_fscanf (FILE * stream, gsl_permutation * p);
int gsl_permutation_fprintf (FILE * stream, const gsl_permutation * p, const char *format);

size_t gsl_permutation_size (const gsl_permutation * p);
size_t * gsl_permutation_data (const gsl_permutation * p);

int gsl_permutation_swap (gsl_permutation * p, const size_t i, const size_t j);

int gsl_permutation_valid (const gsl_permutation * p);
void gsl_permutation_reverse (gsl_permutation * p);
int gsl_permutation_inverse (gsl_permutation * inv, const gsl_permutation * p);
int gsl_permutation_next (gsl_permutation * p);
int gsl_permutation_prev (gsl_permutation * p);
int gsl_permutation_mul (gsl_permutation * p, const gsl_permutation * pa, const gsl_permutation * pb);

int gsl_permutation_linear_to_canonical (gsl_permutation * q, const gsl_permutation * p);
int gsl_permutation_canonical_to_linear (gsl_permutation * p, const gsl_permutation * q);

size_t gsl_permutation_inversions (const gsl_permutation * p);
size_t gsl_permutation_linear_cycles (const gsl_permutation * p);
size_t gsl_permutation_canonical_cycles (const gsl_permutation * q);

INLINE_DECL size_t gsl_permutation_get (const gsl_permutation * p, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
size_t
gsl_permutation_get (const gsl_permutation * p, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= p->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return p->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_PERMUTATION_H__ */

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#pragma once

#include <filesystem>
#include <unordered_map>
#include <optional>
#include <variant>
#include <string>
#include <vector>
#include <gsl/gsl>
#include <DirectML.h>
#include "BucketAllocator.h"

class Model
{
public:
    // When binding a buffer to an operator it is possible to use a subregion of
    // the buffer by specifying an elementOffset, elementCount, and elementSizeInBytes.
    // Additionally, an optional format specifier dictates how to interpret the buffer
    // contents; when omitted the buffer will be interpreted using the same data type used
    // to initialize it.
    struct BufferBindingSource
    {
        std::string name;
        uint64_t elementCount;
        uint64_t elementSizeInBytes;
        uint64_t elementOffset;
        std::optional<DXGI_FORMAT> format;

        // For Append/Consume buffers only:
        std::optional<std::string> counterName;
        uint64_t counterOffsetBytes;
    };

    using Bindings = std::unordered_map<std::string, std::vector<BufferBindingSource>>;

    // RESOURCES
    // ------------------------------------------------------------------------

    struct BufferDesc
    {
        uint64_t sizeInBytes;
        std::vector<std::byte> initialValues;
        DML_TENSOR_DATA_TYPE initialValuesDataType;
        uint64_t initialValuesOffsetInBytes;
    };

    struct ResourceDesc
    {
        std::string name;
        std::variant<BufferDesc> value;
    };

    // DISPATCHABLES
    // ------------------------------------------------------------------------

    struct DmlDispatchableDesc
    {
        struct BindPoint
        {
            std::string name;
            uint32_t resourceCount;
            bool required;
        };

        struct BindPoints
        {
            std::vector<BindPoint> inputs;
            std::vector<BindPoint> outputs;
        };

        DML_OPERATOR_DESC* desc;
        BindPoints bindPoints;
        DML_EXECUTION_FLAGS executionFlags;
        Bindings initBindings;
    };

    struct HlslDispatchableDesc
    {
        enum class Compiler
        {
            DXC
        };

        std::filesystem::path sourcePath;
        Compiler compiler;
        std::vector<std::string> compilerArgs;
    };

    struct DispatchableDesc
    {
        std::string name;
        std::variant<DmlDispatchableDesc, HlslDispatchableDesc> value;
    };

    // COMMANDS
    // ------------------------------------------------------------------------

    struct DispatchCommand
    {
        std::string dispatchableName;
        Bindings bindings;
        std::array<uint32_t, 3> threadGroupCount;
    };

    struct PrintCommand
    {
        std::string resourceName;
    };

    using Command = std::variant<DispatchCommand, PrintCommand>;

    Model() = default;

    Model(
        std::vector<ResourceDesc>&& resourceDescs,
        std::vector<DispatchableDesc>&& dispatchableDescs,
        std::vector<Command>&& commands,
        BucketAllocator&& allocator);

    Model(const Model&) = delete;
    Model& operator=(const Model&) = delete;
    Model(Model&&) = default;
    Model& operator=(Model&&) = default;

    gsl::span<const ResourceDesc> GetResourceDescs() const { return m_resourceDescs; }
    gsl::span<const DispatchableDesc> GetDispatchableDescs() const { return m_dispatchableDescs; }
    gsl::span<const Command> GetCommands() const { return m_commands; }

    const ResourceDesc& GetResource(std::string_view name) const { return *m_resourceDescsByName.find(name.data())->second; }
    const DispatchableDesc& GetDispatchable(std::string_view name) const { return *m_dispatchableDescsByName.find(name.data())->second; }

private:
    std::vector<ResourceDesc> m_resourceDescs;
    std::vector<DispatchableDesc> m_dispatchableDescs;
    std::vector<Command> m_commands;
    BucketAllocator m_allocator;
    std::unordered_map<std::string, ResourceDesc*> m_resourceDescsByName;
    std::unordered_map<std::string, DispatchableDesc*> m_dispatchableDescsByName;
};

max_stars_repo_name: miaobin/DirectML

max_issues_repo_name: miaobin/DirectML

max_forks_repo_name: miaobin/DirectML