text:
//This computes "the" IDCT (inverse discrete cosine transformation) along dim of matrix X.
//This uses a CBLAS matrix multiplication by the DCT-III matrix.

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

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

#ifndef M_SQRT1_2
    #define M_SQRT1_2 0.707106781186547524401
#endif

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

int idct_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 ndct, const int sc);
int idct_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 ndct, const int sc);
int idct_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 ndct, const int sc);
int idct_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 ndct, const int sc);

int idct_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 ndct, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idct_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 (ndct<Lx) { fprintf(stderr,"error in idct_cblas_s: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const float xsc = (sc) ? (float)(2.0*M_SQRT1_2) : 1.0f;
        const float ysc = (sc) ? 2.0f/sqrtf((float)(2u*ndct)) : 2.0f/(float)(2u*ndct);
        const float dcsc = 0.5f * xsc * ysc;

        //Initialize DCT-III matrix
        const size_t LN = Lx * ndct;
        const float P_N = (float)(M_PI/(double)ndct);
        float *DCT, x0;
        DCT = (float *)aligned_alloc(sizeof(float),LN*sizeof(float));
        if (!DCT) { fprintf(stderr,"error in idct_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<ndct; ++n)
        {
            *DCT++ = xsc * ysc;
            for (size_t l=1u; l<Lx; ++l, ++DCT)
            {
                *DCT = ysc * cosf(P_N*(0.5f+(float)n)*(float)l);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            x0 = dcsc * *X;
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,1,0.0f,Y,1);
            for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0; }
        }
        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=V; v>0u; --v, X+=Lx)
                {
                    x0 = dcsc * *X;
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,1,0.0f,Y,1);
                    for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, Y-=K*ndct-1u)
                    {
                        x0 = dcsc * *X;
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,(int)K,0.0f,Y,(int)K);
                        for (size_t n=ndct; n>0u; --n, Y+=K) { *Y -= x0; }
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

int idct_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 ndct, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idct_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 (ndct<Lx) { fprintf(stderr,"error in idct_cblas_d: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
     else
    {
        //Scaling
        const double xsc = (sc) ? 2.0*M_SQRT1_2 : 1.0;
        const double ysc = (sc) ? 2.0/sqrt((double)(2u*ndct)) : 2.0/(double)(2u*ndct);
        const double dcsc = 0.5 * xsc * ysc;

        //Initialize DCT-III matrix
        const size_t LN = Lx * ndct;
        const double P_N = M_PI/(double)ndct;
        double *DCT, x0;
        DCT = (double *)aligned_alloc(sizeof(double),LN*sizeof(double));
        if (!DCT) { fprintf(stderr,"error in idct_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<ndct; ++n)
        {
            *DCT++ = xsc * ysc;
            for (size_t l=1u; l<Lx; ++l, ++DCT)
            {
                *DCT = ysc * cos(P_N*(0.5+(double)n)*(double)l);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            x0 = dcsc * *X;
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,1,0.0,Y,1);
            for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0; }
        }
        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=V; v>0u; --v, X+=Lx)
                {
                    x0 = dcsc * *X;
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,1,0.0,Y,1);
                    for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, Y-=K*ndct-1u)
                    {
                        x0 = dcsc * *X;
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,(int)K,0.0,Y,(int)K);
                        for (size_t n=ndct; n>0u; --n, Y+=K) { *Y -= x0; }
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

int idct_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 ndct, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idct_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 (ndct<Lx) { fprintf(stderr,"error in idct_cblas_c: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const float xsc = (sc) ? (float)(2.0*M_SQRT1_2) : 1.0f;
        const float ysc = (sc) ? 2.0f/sqrtf((float)(2u*ndct)) : 2.0f/(float)(2u*ndct);
        const float dcsc = 0.5f * xsc * ysc;

        //Initialize DCT-III matrix
        const size_t LN = Lx * ndct;
        const float P_N = (float)(M_PI/(double)ndct);
        float *DCT, x0r, x0i;
        DCT = (float *)aligned_alloc(sizeof(float),LN*sizeof(float));
        if (!DCT) { fprintf(stderr,"error in idct_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<ndct; ++n)
        {
            *DCT++ = xsc * ysc;
            for (size_t l=1u; l<Lx; ++l, ++DCT)
            {
                *DCT = ysc * cosf(P_N*(0.5f+(float)n)*(float)l);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            x0r = dcsc * *X; x0i = dcsc * *(X+1);
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2,0.0f,Y,2);
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X+1,2,0.0f,Y+1,2);
            for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0r; *++Y -= x0i; }
        }
        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, X+=2u*Lx)
                {
                    x0r = dcsc * *X; x0i = dcsc * *(X+1);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2,0.0f,Y,2);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X+1,2,0.0f,Y+1,2);
                    for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0r; *++Y -= x0i; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2u, Y-=2u*K*ndct-2u)
                    {
                        x0r = dcsc * *X; x0i = dcsc * *(X+1);
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2*(int)K,0.0f,Y,2*(int)K);
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X+1,2*(int)K,0.0f,Y+1,2*(int)K);
                        for (size_t n=ndct; n>0u; --n, Y+=2u*K) { *Y -= x0r; *(Y+1) -= x0i; }
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

int idct_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 ndct, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in idct_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 (ndct<Lx) { fprintf(stderr,"error in idct_cblas_z: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const double xsc = (sc) ? 2.0*M_SQRT1_2 : 1.0;
        const double ysc = (sc) ? 2.0/sqrt((double)(2u*ndct)) : 2.0/(double)(2u*ndct);
        const double dcsc = 0.5 * xsc * ysc;

        //Initialize DCT-III matrix
        const size_t LN = Lx * ndct;
        const double P_N = M_PI/(double)ndct;
        double *DCT, x0r, x0i;
        DCT = (double *)aligned_alloc(sizeof(double),LN*sizeof(double));
        if (!DCT) { fprintf(stderr,"error in idct_cblas_z: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t n=0u; n<ndct; ++n)
        {
            *DCT++ = xsc * ysc;
            for (size_t l=1u; l<Lx; ++l, ++DCT)
            {
                *DCT = ysc * cos(P_N*(0.5+(double)n)*(double)l);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            x0r = dcsc * *X; x0i = dcsc * *(X+1);
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2,0.0,Y,2);
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X+1,2,0.0,Y+1,2);
            for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0r; *++Y -= x0i; }
        }
        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, X+=2u*Lx)
                {
                    x0r = dcsc * *X; x0i = dcsc * *(X+1);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2,0.0,Y,2);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X+1,2,0.0,Y+1,2);
                    for (size_t n=ndct; n>0u; --n, ++Y) { *Y -= x0r; *++Y -= x0i; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2u, Y-=2u*K*ndct-2u)
                    {
                        x0r = dcsc * *X; x0i = dcsc * *(X+1);
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2*(int)K,0.0,Y,2*(int)K);
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X+1,2*(int)K,0.0,Y+1,2*(int)K);
                        for (size_t n=ndct; n>0u; --n, Y+=2u*K) { *Y -= x0r; *(Y+1) -= x0i; }
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sf_gamma.h>

int
main(int argc, char *argv[]){
	gsl_permutation *P;
	size_t l;

	if(2!=argc){
		return 1;
	}

	l=strlen(argv[1]);

	if(GSL_SF_FACT_NMAX<l){
		printf("String too long\n");
		return 1;
	}

	P=gsl_permutation_calloc(l);
	assert(NULL!=P);

	do{
		size_t *d=gsl_permutation_data(P),i;

		for(i=0;i<l;i++)
			printf("%c",argv[1][d[i]]);

		printf("\n");
	}while(GSL_SUCCESS==gsl_permutation_next(P));

	return 0;
}

ext: c


text:
#pragma once

#include <csignal>
#include <iostream>
#include <sstream>
#include <cassert>
#include <vector>
#include <algorithm>
#include <array>
#include <memory>
#include <iomanip>
#include <fstream>
#include <filesystem>
#include <cstdarg>
#include <unordered_set>
#include <unordered_map>
#include <bitset>
#include <numeric>
#include <cmath>

#ifdef _MSC_VER
#	pragma  warning(disable:4267)
#endif

#include <GLTFAsset.h>
#include <GLTFScene.h>
#include <GLTFBuffer.h>
#include <GLTFBufferView.h>
#include <GLTFAccessor.h>
#include <GLTFMesh.h>
#include <GLTFPrimitive.h>

#ifdef _MSC_VER
#	pragma  warning(default:4267)
#endif

#include "rapidjson/document.h"
#include "rapidjson/prettywriter.h"
#include "rapidjson/stringbuffer.h"
#include "rapidjson/writer.h"

#include <gsl/span>

#include <coveo/linq.h>
#include <coveo/enumerable.h>

#include <maya/MPxCommand.h>
#include <maya/MFnPlugin.h>
#include <maya/MIOStream.h>
#include <maya/MGlobal.h>
#include <maya/MArgList.h>
#include <maya/MArgDatabase.h>
#include <maya/MSyntax.h>
#include <maya/MStreamUtils.h>
#include <maya/MFileObject.h>
#include <maya/MFileIO.h>
#include <maya/MSelectionList.h>
#include <maya/MFnMesh.h>
#include <maya/MItMeshPolygon.h>
#include <maya/MFloatPointArray.h>
#include <maya/MFloatVectorArray.h>
#include <maya/MDagPath.h>
#include <maya/MPointArray.h>
#include <maya/MFnAttribute.h>
#include <maya/MFnStringArrayData.h>
#include <maya/MFnMessageAttribute.h>
#include <maya/MFnNumericAttribute.h>
#include <maya/MFnTypedAttribute.h>
#include <maya/MDagPathArray.h>
#include <maya/MFnSet.h>
#include <maya/MFnSingleIndexedComponent.h>
#include <maya/MFnComponentListData.h>
#include <maya/MDagModifier.h>
#include <maya/MMatrix.h>
#include <maya/MFnMatrixData.h>
#include <maya/MFnTransform.h>
#include <maya/MItDependencyGraph.h>
#include <maya/MFnBlendShapeDeformer.h>
#include <maya/MFnPhongShader.h>
#include <maya/MFnLambertShader.h>
#include <maya/MFnBlinnShader.h>
#include <maya/MUuid.h>
#include <maya/MImage.h>
#include <maya/MFloatMatrix.h>
#include <maya/MTime.h>
#include <maya/MAnimControl.h>
#include <maya/M3dView.h>
#include <maya/MFnSkinCluster.h>
#include <maya/MItGeometry.h>
#include <maya/MItDependencyNodes.h>
#include <maya/MItMeshFaceVertex.h>

ext: h


text:
/*
 * Copyright (C) 2015  University of Oregon
 *
 * You may distribute under the terms of either the GNU General Public
 * License or the Apache License, as specified in the LICENSE file.
 *
 * For more information, see the LICENSE file.
 */

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

void vnmr_svd(double **a, int irow, int jcol, double w[], double **v)
{
   gsl_matrix *gA;
   gsl_matrix *gV;
   gsl_vector *gS;
   int i,j;
   int res;
   gA = gsl_matrix_alloc(irow,jcol);
   for (i=0; i<irow; i++)
      for (j=0; j<jcol; j++)
         gsl_matrix_set(gA, i, j, a[i+1][j+1]);

   gV = gsl_matrix_alloc(jcol,jcol);
   gS = gsl_vector_alloc(jcol);

   res = gsl_linalg_SV_decomp_jacobi(gA, gV, gS);

   for (i=0; i<irow; i++)
      for (j=0; j<jcol; j++)
         a[i+1][j+1] = gsl_matrix_get( gA, i, j);
   for (i=0; i<jcol; i++)
      for (j=0; j<jcol; j++)
         v[i+1][j+1] = gsl_matrix_get( gV, i, j);
   for (i=0; i<jcol; i++)
      w[i+1] = gsl_vector_get( gS, i);

   gsl_matrix_free(gA);
   gsl_matrix_free(gV);
   gsl_vector_free(gS);
}

void vnmr_svd_solve(double **u, double w[], double **v, int irow, int jcol, double b[], double x[])
{
   gsl_matrix *gU;
   gsl_matrix *gV;
   gsl_vector *gS;
   gsl_vector *gB;
   gsl_vector *gX;
   int i,j;
   int res;
   gU = gsl_matrix_alloc(irow,jcol);
   for (i=0; i<irow; i++)
      for (j=0; j<jcol; j++)
         gsl_matrix_set(gU, i, j, u[i+1][j+1]);

   gV = gsl_matrix_alloc(jcol,jcol);
   for (i=0; i<jcol; i++)
      for (j=0; j<jcol; j++)
         gsl_matrix_set( gV, i, j, v[i+1][j+1]);

   gS = gsl_vector_alloc(jcol);
   for (i=0; i<jcol; i++)
      gsl_vector_set( gS, i, w[i+1]);

   gB = gsl_vector_alloc(irow);
   for (i=0; i<irow; i++)
      gsl_vector_set( gB, i, b[i+1]);

   gX = gsl_vector_alloc(jcol);
   res = gsl_linalg_SV_solve(gU, gV, gS, gB, gX);

   for (i=0; i<jcol; i++)
      x[i+1] = gsl_vector_get( gX, i);

   gsl_matrix_free(gU);
   gsl_matrix_free(gV);
   gsl_vector_free(gS);
   gsl_vector_free(gB);
   gsl_vector_free(gX);
fprintf(stderr,"solve m=%d n=%d\n",irow, jcol);
for (i=1; i<=jcol; i++)
     fprintf(stderr," w[%d]: %g\n",i,w[i]);
for (i=0; i<jcol; i++)
     fprintf(stderr," x[%d]: %g\n",i+1,x[i+1]);
}

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "SpotLight.h"

namespace Library
{
	class ProxyModel;
}

namespace Rendering
{
	class SpotLightMaterial;

	class SpotLightDemo final : public Library::DrawableGameComponent
	{
	public:
		SpotLightDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		SpotLightDemo(const SpotLightDemo&) = delete;
		SpotLightDemo(SpotLightDemo&&) = default;
		SpotLightDemo& operator=(const SpotLightDemo&) = default;
		SpotLightDemo& operator=(SpotLightDemo&&) = default;
		~SpotLightDemo();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float SpotLightIntensity() const;
		void SetSpotLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightPosition() const;
		const DirectX::XMVECTOR LightPositionVector() const;
		void SetLightPosition(const DirectX::XMFLOAT3& position);
		void SetLightPosition(DirectX::FXMVECTOR position);

		const DirectX::XMFLOAT3& LightLookAt() const;
		void RotateSpotLight(const DirectX::XMFLOAT2& amount);

		float LightRadius() const;
		void SetLightRadius(float radius);

		float SpecularIntensity() const;
		void SetSpecularIntensity(float intensity);

		float SpecularPower() const;
		void SetSpecularPower(float power);

		float SpotLightInnerAngle() const;
		void SetSpotLightInnerAngle(float angle);

		float SpotLightOuterAngle() const;
		void SetSpotLightOuterAngle(float angle);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		std::shared_ptr<SpotLightMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		std::uint32_t mVertexCount{ 0 };
		Library::SpotLight mSpotLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		float mModelRotationAngle{ 0.0f };
		bool mUpdateMaterial{ true };
	};
}


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

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_blas_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

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  double * data;
  gsl_block * block;
  int owner;
} gsl_matrix;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_view;

typedef _gsl_matrix_view gsl_matrix_view;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_const_view;

typedef const _gsl_matrix_const_view gsl_matrix_const_view;

/* Allocation */

GSL_FUN gsl_matrix *
gsl_matrix_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix *
gsl_matrix_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix *
gsl_matrix_alloc_from_block (gsl_block * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix *
gsl_matrix_alloc_from_matrix (gsl_matrix * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector *
gsl_vector_alloc_row_from_matrix (gsl_matrix * m,
                                        const size_t i);

GSL_FUN gsl_vector *
gsl_vector_alloc_col_from_matrix (gsl_matrix * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_free (gsl_matrix * m);

/* Views */

GSL_FUN _gsl_matrix_view
gsl_matrix_submatrix (gsl_matrix * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_view
gsl_matrix_row (gsl_matrix * m, const size_t i);

GSL_FUN _gsl_vector_view
gsl_matrix_column (gsl_matrix * m, const size_t j);

GSL_FUN _gsl_vector_view
gsl_matrix_diagonal (gsl_matrix * m);

GSL_FUN _gsl_vector_view
gsl_matrix_subdiagonal (gsl_matrix * m, const size_t k);

GSL_FUN _gsl_vector_view
gsl_matrix_superdiagonal (gsl_matrix * m, const size_t k);

GSL_FUN _gsl_vector_view
gsl_matrix_subrow (gsl_matrix * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_view
gsl_matrix_subcolumn (gsl_matrix * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_array (double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_array_with_tda (double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_vector (gsl_vector * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_vector_with_tda (gsl_vector * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_submatrix (const gsl_matrix * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_row (const gsl_matrix * m,
                            const size_t i);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_column (const gsl_matrix * m,
                               const size_t j);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_diagonal (const gsl_matrix * m);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subdiagonal (const gsl_matrix * m,
                                    const size_t k);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_superdiagonal (const gsl_matrix * m,
                                      const size_t k);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subrow (const gsl_matrix * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subcolumn (const gsl_matrix * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_array (const double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_array_with_tda (const double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_vector (const gsl_vector * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_vector_with_tda (const gsl_vector * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_set_zero (gsl_matrix * m);
GSL_FUN void gsl_matrix_set_identity (gsl_matrix * m);
GSL_FUN void gsl_matrix_set_all (gsl_matrix * m, double x);

GSL_FUN int gsl_matrix_fread (FILE * stream, gsl_matrix * m) ;
GSL_FUN int gsl_matrix_fwrite (FILE * stream, const gsl_matrix * m) ;
GSL_FUN int gsl_matrix_fscanf (FILE * stream, gsl_matrix * m);
GSL_FUN int gsl_matrix_fprintf (FILE * stream, const gsl_matrix * m, const char * format);

GSL_FUN int gsl_matrix_memcpy(gsl_matrix * dest, const gsl_matrix * src);
GSL_FUN int gsl_matrix_swap(gsl_matrix * m1, gsl_matrix * m2);
GSL_FUN int gsl_matrix_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * dest, const gsl_matrix * src);

GSL_FUN int gsl_matrix_swap_rows(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_swap_rowcol(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_transpose (gsl_matrix * m);
GSL_FUN int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);
GSL_FUN int gsl_matrix_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix * dest, const gsl_matrix * src);

GSL_FUN double gsl_matrix_max (const gsl_matrix * m);
GSL_FUN double gsl_matrix_min (const gsl_matrix * m);
GSL_FUN void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out);

GSL_FUN void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_min_index (const gsl_matrix * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_equal (const gsl_matrix * a, const gsl_matrix * b);

GSL_FUN int gsl_matrix_isnull (const gsl_matrix * m);
GSL_FUN int gsl_matrix_ispos (const gsl_matrix * m);
GSL_FUN int gsl_matrix_isneg (const gsl_matrix * m);
GSL_FUN int gsl_matrix_isnonneg (const gsl_matrix * m);

GSL_FUN double gsl_matrix_norm1 (const gsl_matrix * m);

GSL_FUN int gsl_matrix_add (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_sub (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_mul_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_div_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_scale (gsl_matrix * a, const double x);
GSL_FUN int gsl_matrix_scale_rows (gsl_matrix * a, const gsl_vector * x);
GSL_FUN int gsl_matrix_scale_columns (gsl_matrix * a, const gsl_vector * x);
GSL_FUN int gsl_matrix_add_constant (gsl_matrix * a, const double x);
GSL_FUN int gsl_matrix_add_diagonal (gsl_matrix * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_get_row(gsl_vector * v, const gsl_matrix * m, const size_t i);
GSL_FUN int gsl_matrix_get_col(gsl_vector * v, const gsl_matrix * m, const size_t j);
GSL_FUN int gsl_matrix_set_row(gsl_matrix * m, const size_t i, const gsl_vector * v);
GSL_FUN int gsl_matrix_set_col(gsl_matrix * m, const size_t j, const gsl_vector * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL double   gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x);
GSL_FUN INLINE_DECL double * gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const double * gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
double
gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
double *
gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (double *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const double *
gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_DOUBLE_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#ifndef VKST_PLAT_LOG_H
#define VKST_PLAT_LOG_H

#include <plat/core.h>
#include <plat/filesystem.h>
#include <gsl.h>
#include <string>

namespace plat {

// Initialize logging with output to a given path.
// Must be called before logging output will occur.
void init_logging(plat::filesystem::path logfile, std::error_code& ec) noexcept;

enum class log_severities : uint8_t {
  trace = 1,
  debug,
  info,
  warn,
  error,
  fatal,
  none
};

inline constexpr gsl::czstring to_string(log_severities severity) noexcept {
  switch(severity) {
  case log_severities::trace: return "TRACE";
  case log_severities::debug: return "DEBUG";
  case log_severities::info: return "INFO";
  case log_severities::warn: return "WARN";
  case log_severities::error: return "ERROR";
  case log_severities::fatal: return "FATAL";
  case log_severities::none: return "NONE";
  }
  PLAT_MARK_UNREACHABLE;
}

void log(log_severities severity, gsl::czstring fmt, ...) noexcept;

#define LOG_TRACE(...) ::plat::log(::plat::log_severities::trace, __VA_ARGS__)
#define LOG_DEBUG(...) ::plat::log(::plat::log_severities::debug, __VA_ARGS__)
#define LOG_INFO(...) ::plat::log(::plat::log_severities::info, __VA_ARGS__)
#define LOG_WARN(...) ::plat::log(::plat::log_severities::warn, __VA_ARGS__)
#define LOG_ERROR(...) ::plat::log(::plat::log_severities::error, __VA_ARGS__)
#define LOG_FATAL(...) ::plat::log(::plat::log_severities::fatal, __VA_ARGS__)

#define LOG_ENTER                                                              \
  ::plat::log(::plat::log_severities::trace, "enter: %s (%s:%d)", __func__,    \
              __FILE__, __LINE__)
#define LOG_LEAVE                                                              \
  ::plat::log(::plat::log_severities::trace, "leave: %s (%s:%d)", __func__,    \
              __FILE__, __LINE__)

} // namespace plat

#endif // VKST_PLAT_LOG_H

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst


text:
#pragma once

#include "text_parser.h"
#include "width_cache.h"

#include <string>
#include <vector>

#include <gsl/span>

namespace terminal_editor {

enum class GraphemeKind {
    NORMAL,      ///< Normal, displayable characters.
    INVALID,     ///< Invalid characters for given encoding.
    REPLACEMENT, ///< Replacement representation of valid (possibly control) characters (such as 4 spaces for tabs, or [NUL] for 0x00).
};

/// Each grapheme represents one logical 'image' on the screen.
/// It can consist of many actual characters (either because we have combining characters, or because we replaced bytes with some special string, like "[CR]").
/// It is produced from one or more bytes in the underlying byte data.
/// It can have a width of 0 or more (either because actual characters take two columns, or because replacement string is longer, like "[x66]" or "[TAB]").
struct Grapheme {
    GraphemeKind kind;    ///< Kind of grapheme: normal grapheme, invalid bytes or replacement.
    std::string rendered; ///< Valid UTF-8 string to display on the screen.
    std::string info;     ///< Valid UTF-8 string. Can contain arbitrary informational text about the grapheme. Will contain error description in case of invalid graphemes.
    int width;            ///< Width (in terminal cells) of the 'rendered' string, once it will be displayed on the terminal.
    gsl::span<const char> consumedInput; ///< Span of input data that was rendered into this grapheme.
};

/// Converts given span of CodePoinInfos into a Grapheme, by concatenating their representations.
/// @param codePoinInfos    A span of CodePoinInfos.
///                         @note If codePoinInfos is empty a NORMAL, but zero width grapheme is returned.
Grapheme renderGrapheme(gsl::span<const CodePointInfo> codePointInfos);

/// CodePointWidthCache used by renderLine().
/// @todo This should not be a gloabl variable.
extern CodePointWidthCache textRendererWidthCache;

/// Renders data into Graphemes.
/// Each byte of invalid CodePointInfos are rendered as separate graphemes.
/// Valid CodePointInfos are grouped into maximal chunks where only first CodePointInfo can have non-zero wcwidth() (after processing replacements). One grapheme is created for
/// each such group.
/// @note Some resulting graphemes can have zero-width.
/// @note Result can be empty.
/// @note Width of characters is computed using global CodePointWidthCache.
/// @param data     Input string to render. It is assumed to be in UTF-8, but can contain invalid characters (which will be rendered as special Graphemes).
///                 @note Any control characters (including new line characters) will be rendered as a replacement string (i.e. [LF]).
std::vector<Grapheme> renderLine(gsl::span<CodePointInfo> codePointInfos);

/// Returns concatenation of rendered property of all graphemes.
/// @param graphemes    Span of graphemes to concatenate.
/// @param useBrackets  If true all invalid and replacement sequences will be enclosed with brackets.
std::string renderGraphemes(gsl::span<const Grapheme> graphemes, bool useBrackets);

/// Returns width of given text after rendering on screen.
/// @note This functions takes into consideration replacement strings, and thus differs from wcswidth().
/// @note This function is very slow, as it needs to convert text to graphemes first.
/// @param text     UTF-8 string, can be invalid.
int getRenderedWidth(gsl::span<const char> text);

/// Returns width of given line of graphemes after rendering on screen.
/// @param graphemes    A span of graphemes.
int getRenderedWidth(gsl::span<const Grapheme> graphemes);

} // namespace terminal_editor

max_stars_repo_name: Zbyl/terminal-editor

max_issues_repo_name: Zbyl/terminal-editor

max_forks_repo_name: Zbyl/terminal-editor


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#define DIM 1
#define NPTS 100

#define abserr 1e-5
#define relerr 1e-5
#define mem 10000
typedef double (* function) (double x, void * params);
double integrate(function f, void * params, double lower_bound, double upper_bound);

//int func (double t, const double y[], double f[], void *params);
//void gsl_ode(gsl_odeiv_system sys, double Xo, double maxtime);

max_stars_repo_name: cboettig/AdaptiveDynamics

max_issues_repo_name: cboettig/AdaptiveDynamics

max_forks_repo_name: cboettig/AdaptiveDynamics


text:
//CELL operation: does output side of Elman RNN layer.
//There is no strict definition of Elman RNN in the literature.
//For example, the Wikipedia article doesn't specify the dims of the hidden or other vectors.
//Here I implement an interpretation where there are N neurons and thus N inputs/outputs.
//The inputs here are from the IN stage, so reduced to N driving input time-series in X.

//Again, this is not an "Elman network", rather a layer of N neurons
//that I have named "Elman" neurons due to their great similarity to an Elman RNN.

//To do: should I allow other output activations other than logistic?

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

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

int elman_s (float *Y, const float *X, const float *U, float *H, const float *W, const float *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int elman_d (double *Y, const double *X, const double *U, double *H, const double *W, const double *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int elman_s (float *Y, const float *X, const float *U, float *H, const float *W, const float *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;

    if (N==1u)
    {
        for (size_t t=0u; t<T; ++t)
        {
            H[0] = 1.0f / (1.0f+expf(-X[t]-U[0]*H[0]));
            Y[t] = 1.0f / (1.0f+expf(-B[0]-W[0]*H[0]));
        }
    }
    else
    {
        float *tmp;
        if (!(tmp=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in elman_s: problem with malloc. "); perror("malloc"); return 1; }

        if (dim==0u)
        {
            if (iscolmajor)
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t*N],1,tmp,1);
                    cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-tmp[n])); }
                    //cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t*N],1);
                    //for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-X[t*N+n])); }
                    //for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0f/(1.0f+expf(-cblas_sdot((int)N,&W[n],(int)N,H,1)-B[n])); }
                    cblas_scopy((int)N,B,1,tmp,1);
                    cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0f/(1.0f+expf(-tmp[n])); }
                }
            }
            else
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t],(int)T,tmp,1);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-tmp[n])); }
                    cblas_scopy((int)N,B,1,tmp,1);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0f/(1.0f+expf(-tmp[n])); }
                }
            }
        }
        else if (dim==1u)
        {
            if (iscolmajor)
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t],(int)T,tmp,1);
                    cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-tmp[n])); }
                    cblas_scopy((int)N,B,1,tmp,1);
                    cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0f/(1.0f+expf(-tmp[n])); }
                }
            }
            else
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t*N],1,tmp,1);
                    cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-tmp[n])); }
                    cblas_scopy((int)N,B,1,tmp,1);
                    cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0f/(1.0f+expf(-tmp[n])); }
                }
            }
        }
        else
        {
            fprintf(stderr,"error in elman_s: dim must be 0 or 1.\n"); return 1;
        }
    }

    return 0;
}

int elman_d (double *Y, const double *X, const double *U, double *H, const double *W, const double *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;

    if (N==1u)
    {
        for (size_t t=0u; t<T; ++t)
        {
            H[0] = 1.0 / (1.0+exp(-X[t]-U[0]*H[0]));
            Y[t] = 1.0 / (1.0+exp(-B[0]-W[0]*H[0]));
        }
    }
    else
    {
        double *tmp;
        if (!(tmp=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in elman_d: problem with malloc. "); perror("malloc"); return 1; }

        if (dim==0u)
        {
            if (iscolmajor)
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t*N],1,tmp,1);
                    cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-tmp[n])); }
                    cblas_dcopy((int)N,B,1,tmp,1);
                    cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0/(1.0+exp(-tmp[n])); }
                }
            }
            else
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t],(int)T,tmp,1);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-tmp[n])); }
                    cblas_dcopy((int)N,B,1,tmp,1);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0/(1.0+exp(-tmp[n])); }
                }
            }
        }
        else if (dim==1u)
        {
            if (iscolmajor)
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t],(int)T,tmp,1);
                    cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-tmp[n])); }
                    cblas_dcopy((int)N,B,1,tmp,1);
                    cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0/(1.0+exp(-tmp[n])); }
                }
            }
            else
            {
                for (size_t t=0u; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t*N],1,tmp,1);
                    cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-tmp[n])); }
                    cblas_dcopy((int)N,B,1,tmp,1);
                    cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,tmp,1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0/(1.0+exp(-tmp[n])); }
                }
            }
        }
        else
        {
            fprintf(stderr,"error in elman_d: dim must be 0 or 1.\n"); return 1;
        }
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


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

/* Author:  G. Jungman, modified by C. Krueger */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/* Factorise a general complex-valued M x N matrix A into
 *
 *   A = Q R
 *
 * where Q is unitary (M x M) and R is upper triangular (M x N).
 *
 * Q is stored as a packed set of Householder transformations in the
 * strict lower triangular part of the input matrix.
 *
 * R is stored in the diagonal and upper triangle of the input matrix.
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_k .. Q_2 Q_1
 *
 * where k = MIN(M,N) and
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
 *
 * This storage scheme is the same as in LAPACK.  */

int
gsl_linalg_complex_QR_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau->size != N)
    {
      GSL_ERROR ("size of tau must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          /* Compute the Householder transformation to reduce the j-th
             column of the matrix to a multiple of the j-th unit vector */

          gsl_vector_complex_view c = gsl_matrix_complex_subcolumn (A, i, i, M - i);
          gsl_complex tau_i = gsl_linalg_complex_householder_transform (&(c.vector));

          gsl_vector_complex_set (tau, i, tau_i);

          /* apply the transformation to the remaining columns and update the norms */

          if (i + 1 < N)
            {
              gsl_matrix_complex_view m = gsl_matrix_complex_submatrix (A, i, i + 1, M - i, N - i - 1);
              gsl_complex tau_i_conj = gsl_complex_conjugate(tau_i);
              gsl_vector_complex_view work = gsl_vector_complex_subvector(tau, i + 1, N - i - 1);

              gsl_linalg_complex_householder_left(tau_i_conj, &(c.vector), &(m.matrix), &(work.vector));
            }
        }

      return GSL_SUCCESS;
    }
}

/* Solves the system A x = b using the QR factorisation,

 *  R x = Q^H b
 *
 * to obtain x.
 */

int
gsl_linalg_complex_QR_solve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
                             const gsl_vector_complex * b, gsl_vector_complex * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* copy x <- b */
      gsl_vector_complex_memcpy (x, b);

      /* solve for x */
      gsl_linalg_complex_QR_svx (QR, tau, x);

      return GSL_SUCCESS;
    }
}

/* Solves the system A x = b in place using the QR factorisation,

 *  R x = Q^H b
 *
 * to obtain x.
 */

int
gsl_linalg_complex_QR_svx (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * x)
{

  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* compute rhs = Q^H b */
      gsl_linalg_complex_QR_QHvec (QR, tau, x);

      /* solve R x = rhs, storing x in-place */
      gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      return GSL_SUCCESS;
    }
}

/* Find the least squares solution to the overdetermined system
 *
 *   A x = b
 *
 * for M >= N using the QR factorization A = Q R.
 */

int
gsl_linalg_complex_QR_lssolve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
                               const gsl_vector_complex * b, gsl_vector_complex * x,
                               gsl_vector_complex * residual)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", 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_complex_const_view R = gsl_matrix_complex_const_submatrix (QR, 0, 0, N, N);
      gsl_vector_complex_view c = gsl_vector_complex_subvector(residual, 0, N);

      gsl_vector_complex_memcpy(residual, b);

      /* compute rhs = Q^H b */
      gsl_linalg_complex_QR_QHvec (QR, tau, residual);

      /* solve R x = rhs */
      gsl_vector_complex_memcpy(x, &(c.vector));
      gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);

      /* compute residual = b - A x = Q (Q^H b - R x) */
      gsl_vector_complex_set_zero(&(c.vector));
      gsl_linalg_complex_QR_Qvec(QR, tau, residual);

      return GSL_SUCCESS;
    }
}

/* form the product v := Q^H v from a QR factorized matrix */

int
gsl_linalg_complex_QR_QHvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (tau->size != N)
    {
      GSL_ERROR ("size of tau must be N", GSL_EBADLEN);
    }
  else if (v->size != M)
    {
      GSL_ERROR ("vector size must be M", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q^H v */

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          gsl_vector_complex_const_view h = gsl_matrix_complex_const_subcolumn (QR, i, i, M - i);
          gsl_vector_complex_view w = gsl_vector_complex_subvector (v, i, M - i);
          gsl_complex ti = gsl_vector_complex_get (tau, i);
          gsl_complex ti_conj = gsl_complex_conjugate(ti);
          gsl_linalg_complex_householder_hv (ti_conj, &(h.vector), &(w.vector));
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_QR_Qvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (v->size != M)
    {
      GSL_ERROR ("vector size must be M", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q v */

      for (i = GSL_MIN (M, N); i-- > 0;)
        {
          gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i);
          gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&(c.vector),
                                                                i, M - i);
          gsl_vector_complex_view w = gsl_vector_complex_subvector (v, i, M - i);
          gsl_complex ti = gsl_vector_complex_get (tau, i);
          /* we do not need the conjugate of ti here */
          gsl_linalg_complex_householder_hv (ti, &h.vector, &w.vector);
        }
      return GSL_SUCCESS;
    }
}

/* form the unitary matrix Q from the packed QR matrix */
int
gsl_linalg_complex_QR_unpack (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
                              gsl_matrix_complex * Q, gsl_matrix_complex * R)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
    }
  else if (R->size1 != M || R->size2 != N)
    {
      GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
    }
  else if (tau->size != N)
    {
      GSL_ERROR ("size of tau must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      /* initialize Q to the identity */
      gsl_matrix_complex_set_identity (Q);

      for (i = GSL_MIN (M, N); i-- > 0;)
        {
          gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i);
          gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&c.vector, i, M - i);
          gsl_matrix_complex_view m = gsl_matrix_complex_submatrix (Q, i, i, M - i, M - i);
          gsl_complex ti = gsl_vector_complex_get (tau, i);
          gsl_vector_complex_view work = gsl_matrix_complex_subcolumn(R, 0, 0, M - i);

          /* we do not need the conjugate of ti here */
          gsl_linalg_complex_householder_left (ti, &h.vector, &m.matrix, &work.vector);
        }

      /* form the right triangular matrix R from a packed QR matrix */

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < i && j < N; j++)
            gsl_matrix_complex_set (R, i, j, GSL_COMPLEX_ZERO);

          for (j = i; j < N; j++)
            gsl_matrix_complex_set (R, i, j, gsl_matrix_complex_get (QR, i, j));
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* combination/gsl_combination.h
 * based on permutation/gsl_permutation.h by Brian Gough
 *
 * Copyright (C) 2001 Szymon Jaroszewicz
 *
 * 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_COMBINATION_H__
#define __GSL_COMBINATION_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_types.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_combination_struct
{
  size_t n;
  size_t k;
  size_t *data;
};

typedef struct gsl_combination_struct gsl_combination;

gsl_combination *gsl_combination_alloc (const size_t n, const size_t k);
gsl_combination *gsl_combination_calloc (const size_t n, const size_t k);
void gsl_combination_init_first (gsl_combination * c);
void gsl_combination_init_last (gsl_combination * c);
void gsl_combination_free (gsl_combination * c);
int gsl_combination_memcpy (gsl_combination * dest, const gsl_combination * src);

int gsl_combination_fread (FILE * stream, gsl_combination * c);
int gsl_combination_fwrite (FILE * stream, const gsl_combination * c);
int gsl_combination_fscanf (FILE * stream, gsl_combination * c);
int gsl_combination_fprintf (FILE * stream, const gsl_combination * c, const char *format);

size_t gsl_combination_n (const gsl_combination * c);
size_t gsl_combination_k (const gsl_combination * c);
size_t * gsl_combination_data (const gsl_combination * c);

size_t gsl_combination_get (const gsl_combination * c, const size_t i);

int gsl_combination_valid (gsl_combination * c);
int gsl_combination_next (gsl_combination * c);
int gsl_combination_prev (gsl_combination * c);

#ifdef HAVE_INLINE

extern inline
size_t
gsl_combination_get (const gsl_combination * c, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= c->k)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return c->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_COMBINATION_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _GLOBAL_SETTING_
#define _GLOBAL_SETTING_

#define MATRIX_WIDTH 128
#define MATRIX_HEIGHT 128
#define MATRIX_CHANNEL 64
#define MATRIX_DEPTH 1

#define KERNEL_WIDTH 3
#define KERNEL_HEIGHT 3
#define KERNEL_CHANNEL MATRIX_CHANNEL
#define KERNEL_DEPTH 128

#define FEATURE_WIDTH (MATRIX_WIDTH - KERNEL_WIDTH + 1)
#define FEATURE_HEIGHT (MATRIX_HEIGHT - KERNEL_HEIGHT + 1)
#define FEATURE_CHANNEL 1
#define FEATURE_DEPTH 128

#ifdef CBLAS
    #include <cblas.h>
#endif  // #ifdef CBLAS

#ifdef OMP
	#include <omp.h>
#endif	// #ifdef OMP

#ifdef OPENCL
	#ifdef __APPLE__
		#include <OpenCL/opencl.h>
	#elif linux
		#include <CL/cl.h>
	#elif __WIN32__
		#include <cl.h>
	#else
		#include <CL/cl.h>
	#endif	//	#ifdef __APPLE__
#endif	//	#ifdef OPENCL

#endif	//	_GLOBAL_SETTING_


text:
/* bclsqr.c
   $Revision: 273 $ $Date: 2006-09-04 15:59:04 -0700 (Mon, 04 Sep 2006) $

   ----------------------------------------------------------------------
   This file is part of BCLS (Bound-Constrained Least Squares).

   Copyright (C) 2006 Michael P. Friedlander, Department of Computer
   Science, University of British Columbia, Canada. All rights
   reserved. E-mail: <mpf@cs.ubc.ca>.

   BCLS is free software; you can redistribute it and/or modify it
   under the terms of the GNU Lesser General Public License as
   published by the Free Software Foundation; either version 2.1 of the
   License, or (at your option) any later version.

   BCLS 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 Lesser General
   Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with BCLS; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
   USA
   ----------------------------------------------------------------------
*/
/*!
   \file
   Interface to LSQR routine.  Used to compute a Newton step.
*/

#include <cblas.h>
#include <string.h>
#include <stdio.h>

#include "bcls.h"
#include "bclib.h"
#include "bclsqr.h"
#include "lsqr.h"

/*!

  \brief Mat-vec routine called by LSQR.

  LSQR calls this routine, which in turn calls bcls_aprod: This
  routine is declared "static" so that it won't be confused with the
  user's own Aprod routine.

  - If mode = 1,
    - y(1:m)  <-  y(1:m)  +  A(:,ix) * dxFree.
    - y(m+1:) <-  y(m+1:) +     damp * dxFree

  - If mode = 2,
    - dxFree <- dxFree + A'* y(1:m) + damp y(m+1:).

  Note that mSubProb is the number of rows in the subproblem.  This
  may or may not be equal to m, which is the number of rows in the
  original problem.  If the subproblem is damped (because the user has
  either provided a linear term or an explicit damping parameter),
  then mSubProb = m + nFree.  Otherwise, mSubProb = m.

  \param[in]      mode      Determines which producte with A is required.

  \param[in]      mSubProb  Number of rows in the matrix seen by LSQR.  Also:
                            - length of y.
  \param[in]      nFree     Number of columns in A(:,ix).  Also:
                            - length of ix
                            - length of dxFree.
  \param[in,out]  dxFree    Primal variables
  \param[in,out]  y         Dual variables.
  \param[in,out]  UsrWrk    Transit pointer to the BCLS problem context.

*/
static void
aprod_free_lsqr( const int mode, const int mSubProb, const int nFree,
                 double dxFree[], double y[], void *UsrWrk)
{
    int j;
    BCLS *ls = (BCLS *)UsrWrk;  // Reclaim access to the BCLS workspace.
    const int    m      = ls->m;
    const int    preconditioned = ls->Usolve != NULL;
    const int    damped = ls->damp_actual > 0.0;
    const double damp   = ls->damp_actual;
    int          *ix    = ls->ix;
    double       *dx    = ls->dx;    // Used as workspace.
    double       *dy    = ls->wrk_u; // ...

    if (mode == 1) {

        // Solve U dy = dxFree.  U is nFree-by-nFree, and so the
        // relevant part of dy has length nFree.
        if (preconditioned)
            bcls_usolve( ls, BCLS_PRECON_U, nFree, ix, dy, dxFree );
        else
            cblas_dcopy( nFree, dxFree, 1, dy, 1 );

        // y2 <- y2 + damp * dy(1:nFree).
        if (damped)
            cblas_daxpy( nFree, damp, dy, 1, &y[m], 1 );

        // Scatter dy into dx: dx(ix) <- dy(1:nFree).
        for (j = 0; j < nFree; j++)  dx[ ix[j] ] = dy[j];

        // dy <- A(:,ix) * dx(ix).
        bcls_aprod( ls, BCLS_PROD_A, nFree, ix, dx, dy );

        // y1 <- y1 + dy.
        cblas_daxpy( m, 1.0, dy, 1, y, 1 );

    }
    else { // mode == 2

        // dx <- A' * y1.
        bcls_aprod( ls, BCLS_PROD_At, nFree, ix, dx, y );

        // Gather dx(ix) into dy: dy(1:nFree) <- dx(ix).
        for (j = 0; j < nFree; j++) dy[j] = dx[ ix[j] ];

        // dy <- dy + damp * y2.
        if (damped)
            cblas_daxpy( nFree, damp, &y[m], 1, dy, 1 );

        // Solve U' dx = dy.
        if (preconditioned)
            bcls_usolve( ls, BCLS_PRECON_Ut, nFree, ix, dy, dx );
        else
            cblas_dcopy( nFree, dy, 1, dx, 1 );

        // dxFree <- dxFree + dx.
        cblas_daxpy( nFree, 1.0, dx, 1, dxFree, 1 );

    }
    return;
}

/*!

  \brief Compute a Newton step using LSQR.

  \see  bcls_newton_step_cgls.

  \param[in,out] ls      BCLS problem context.
  \param[in]     m       Number of rows in A.
  \param[in]     nFree   Number of columns in A(:,ix).  Also:
                         - length of ix and dxFree
  \param[in]     ix      Index of free variables.
  \param[in]     damp    Regularization parameter.
  \param[in]     itnLim  Iteration limit on current LSQR call.
  \param[in]     tol     LSQR's atol and btol.
  \param[in,out] dxFree  Search direction on free variables.
  \param[in,out] x       Current point.
  \param[in]     c       Linear term.
  \param[in,out] r       Residual. Used as RHS for LSQR.
  \param[out]    itns    Number of LSQR iterations on current subproblem.
  \param[out]    opt     Optimality achieved by LSQR on current subproblem.

  \return
  - 0: Required accurace was achieved.
  - 1: The iteration limit (itnLim) was reached.
  - 2: A(:,ix) is excessively ill-conditioned.

*/
int
bcls_newton_step_lsqr( BCLS *ls, int m, int nFree, int ix[], double damp,
		       int itnLim, double tol, double dxFree[], double x[],
		       double c[], double r[], int *itns, double *opt )
{
    int j, k;                             // Misc. counters.
    int mpn;                              // No. of rows in [ N; damp I ].
    int          unscale_dxFree   = 0;
    const int    rescaling_method = 1;
    const int    linear   = c != NULL;
    const int    damped   = damp  > 0.0;
    const double damp_min = ls->damp_min;
    const double damp2    = damp * damp;
    const double zero     = 0.0;
    double damp_actual;
    double beta1, beta2;

    // LSQR outputs
    int    istop;      // Termination flag.
    double anorm;      // Estimate of Frobenious norm of Abar.
    double acond;      // Estimate of condition no. of Abar.
    double rnorm;      // Estimate of the final value of norm(rbar).
    double xnorm;      // Estimate of the norm of the final solution dx.

    // Set r(m+1:) <- - 1/beta c + damp^2/beta x, where
    //     beta = max(min_damp, damp).
    // Note that r is declared length (m+n), so there is always enough
    // space.  The chosen damping parameter is stored in
    // ls->damp_actual so that it can be used in bcls_aprod_free.
    if (!damped  &&  !linear) {
	mpn = m;
	ls->damp_actual = 0.0;
    }
    else {

	mpn = m + nFree;

	//--------------------------------------------------------------
	// Rescale the RHS.  Will need to unscale LSQR's solution later.
	//--------------------------------------------------------------
	if (rescaling_method) {

	    if (linear) {

		unscale_dxFree  = 1;
		damp_actual     = fmax( damp, damp_min );
		ls->damp_actual = damp_actual;

		// r(1:m) = damp_actual * r(1:m)
		cblas_dscal( m, damp_actual, r, 1 );

		// r(m+1:) = - c
		for (j = 0; j < nFree; j++ )
		    r[m+j] = - c[ ix[j] ];

		// r(m+1:) = - c - damp^2 * x
		if (damped)
		    for (j = 0; j < nFree; j++)
			r[m+j] -= damp2 * x[ ix[j] ];
	    }
	    else {

	        ls->damp_actual = damp;
		for (j = 0; j < nFree; j++)
		    r[m+j] = - damp * x[ ix[j] ];
	    }
	}
	//--------------------------------------------------------------
	// No rescaling.
	//--------------------------------------------------------------
	else {
	    if (damped  &&  linear) {

		damp_actual     =  fmax( damp, damp_min );
		beta1           =  -1.0 / damp_actual;
		beta2           =  -damp * damp / damp_actual;
		ls->damp_actual =  damp_actual;

		for (j = 0; j < nFree; j++) {
		    k = ix[j];
		    r[m+j] = beta1 * c[k] + beta2 * x[k];
		}
	    }
	    else if ( damped  &&  !linear ) {

		beta2 = - damp;
		ls->damp_actual = damp;
		for (j = 0; j < nFree; j++)
		    r[m+j] = beta2 * x[ ix[j] ];
	    }
	    else if (!damped  &&   linear ) {

		beta1 = - 1.0 / damp_min;
		ls->damp_actual = damp_min;
		for (j = 0; j < nFree; j++)
		    r[m+j] = beta1 * c[ ix[j] ];
	    }
	}
    }

    // -----------------------------------------------------------------
    // Solve the subproblem with LSQR.
    // -----------------------------------------------------------------
    bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR]), BCLS_TIMER_START );

    lsqr( mpn, nFree, aprod_free_lsqr, zero, (void *)ls,
          r, ls->wrk_v, ls->wrk_w, dxFree, NULL,
          tol, tol, ls->conlim, itnLim, ls->minor_file,
          &istop, itns, &anorm, &acond,
          &rnorm, opt, &xnorm );

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR]), BCLS_TIMER_STOP );

    // -----------------------------------------------------------------
    // Cleanup and exit.  First unscale dxFree if needed.
    // -----------------------------------------------------------------
    if (rescaling_method  &&  unscale_dxFree)
	cblas_dscal( nFree, 1.0/damp_actual, dxFree, 1 );

    if (istop <= 3)
	return 0;
    else if (istop == 5)
	return 1;
    else
	return 2;
}

max_stars_repo_name: echristakopoulou/glslim

max_issues_repo_name: echristakopoulou/glslim

max_forks_repo_name: echristakopoulou/glslim


text:
#include "asf_sar.h"
#include <asf.h>
#include <asf_meta.h>
#include <asf_raster.h>

#include <glib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_statistics_double.h>
#include <assert.h>

static double min2(double a, double b)
{
    return a<b ? a : b;
}

static double min4(double a, double b, double c, double d)
{
    return min2(min2(a,b), min2(c,d));
}

static double max2(double a, double b)
{
    return a>b ? a : b;
}

static double max4(double a, double b, double c, double d)
{
    return max2(max2(a,b), max2(c,d));
}

static void ts2ls(meta_parameters *meta, double time, double slant,
                  double *pLine, double *pSamp)
{
  double lat, lon;

  meta_timeSlantDop2latLon(meta, time, slant, 0, 0, &lat, &lon);
  meta_get_lineSamp(meta, lat, lon, 0, pLine, pSamp);

  //printf("ts2ls: %f,%f -> %f,%f -> %f,%f\n",
  //       time,slant, lat,lon, *pLine,*pSamp);
}

static int
proj_to_sr(const char *infile, const char *outfile, double pixel_size)
{
    int ii, jj, kk;
    const float_image_sample_method_t sampling_method =
        FLOAT_IMAGE_SAMPLE_METHOD_BILINEAR;

    // overall algorithm:
    // 1. find extents in time/slant space
    // 2. for each pixel in output, resample in input space

    meta_parameters *inMeta = meta_read(infile);
    int nl = inMeta->general->line_count;
    int ns = inMeta->general->sample_count;

    if (!inMeta->projection && !inMeta->transform)
        asfPrintError("Expected a projection/transform block!\n");
    if (!inMeta->state_vectors)
        asfPrintError("Input data does not have state vectors!\n");

    //asfPrintStatus("Converting %s to slant range...\n", infile);

    // first, find extents in time/slant space
    // do this by projecting image corners to time/slant
    int tl_x=0,    tl_y=0;
    int tr_x=ns-1, tr_y=0;
    int bl_x=0,    bl_y=nl-1;
    int br_x=ns-1, br_y=nl-1;

    // we have to find the "real" corners of the image
    // do this using the first band of the input image as a reference
    if (inMeta->general->band_count == 1)
        asfPrintStatus("Tiling the input image...\n");
    else
        asfPrintStatus("Tiling the reference band of the input image...\n");

    FloatImage *in = float_image_new_from_metadata(inMeta, infile);

    // find top left pixel -- TOP-most non-no-data pixel in the image
    for (ii=0; ii<nl; ++ii)
        for (jj=0; jj<ns; ++jj) {
            double val = float_image_get_pixel(in, jj, ii);
            if (val != inMeta->general->no_data && val != 0.0) {
                tl_x = jj; tl_y = ii;
                goto found_tl;
            }
        }

    asfPrintError("Couldn't find top-left pixel! Entire image no data?\n");

  found_tl:

    // find top right pixel -- RIGHT-most non-no-data pixel in the image
    for (jj=ns-1; jj>=0; --jj)
        for (ii=0; ii<nl; ++ii) {
            double val = float_image_get_pixel(in, jj, ii);
            if (val != inMeta->general->no_data && val != 0.0) {
                tr_x = jj; tr_y = ii;
                goto found_tr;
            }
        }

    asfPrintError("Couldn't find top-right pixel! Entire image no data?\n");

  found_tr:

    // find bottom left pixel -- LEFT-most non-no-data pixel in the image
    for (jj=0; jj<ns; ++jj)
        for (ii=nl-1; ii>=0; --ii) {
            double val = float_image_get_pixel(in, jj, ii);
            if (val != inMeta->general->no_data && val != 0.0) {
                bl_x = jj; bl_y = ii;
                goto found_bl;
            }
        }

    asfPrintError("Couldn't find bottom-left pixel! Entire image no data?\n");

  found_bl:

    // find bottom right pixel -- BOTTOM-most non-no-data pixel in the image
    for (ii=nl-1; ii>=0; --ii)
        for (jj=ns-1; jj>=0; --jj) {
            double val = float_image_get_pixel(in, jj, ii);
            if (val != inMeta->general->no_data && val != 0.0) {
                br_x = jj; br_y = ii;
                goto found_br;
            }
        }

    asfPrintError("Couldn't find bottom-right pixel! Entire image no data?\n");

  found_br:

    asfPrintStatus("Determining image extents in time/slant coordinates.\n");
    //asfPrintStatus("Corners are at: TL (%d,%d)\n", tl_y, tl_x);
    //asfPrintStatus(" (line,sample)  TR (%d,%d)\n", tr_y, tr_x);
    //asfPrintStatus("                BL (%d,%d)\n", bl_y, bl_x);
    //asfPrintStatus("                BR (%d,%d)\n", br_y, br_x);

    double tl_time, tl_slant;
    double tr_time, tr_slant;
    double bl_time, bl_slant;
    double br_time, br_slant;

    meta_get_timeSlantDop(inMeta, tl_y, tl_x, &tl_time, &tl_slant, NULL);
    meta_get_timeSlantDop(inMeta, tr_y, tr_x, &tr_time, &tr_slant, NULL);
    meta_get_timeSlantDop(inMeta, bl_y, bl_x, &bl_time, &bl_slant, NULL);
    meta_get_timeSlantDop(inMeta, br_y, br_x, &br_time, &br_slant, NULL);

    //asfPrintStatus("Corners are at: TL (%f,%f)\n", tl_time, tl_slant);
    //asfPrintStatus(" (time,slant)   TR (%f,%f)\n", tr_time, tr_slant);
    //asfPrintStatus("                BL (%f,%f)\n", bl_time, bl_slant);
    //asfPrintStatus("                BR (%f,%f)\n", br_time, br_slant);

    double slant_start = min4(tl_slant, tr_slant, bl_slant, br_slant);
    double slant_end = max4(tl_slant, tr_slant, bl_slant, br_slant);

    double time_min = min4(tl_time, tr_time, bl_time, br_time);
    double time_max = max4(tl_time, tr_time, bl_time, br_time);

    double slant_incr;
    double time_start, time_end, time_incr;
    int onl, ons;

    if (pixel_size > 0) {
        slant_incr = pixel_size;
        ons = (slant_end - slant_start) / slant_incr;

        if (inMeta->sar) {
            // in this case, the original data has a SAR block, we will use the
            // same azimuth time per pixel.
            time_incr = inMeta->sar->azimuth_time_per_pixel;

            // we always want to be DECREASING in time
            // latest time is on top (line 1), earliest on bottom (line ONL)
            if (time_incr > 0) {
                time_incr = -time_incr;
                inMeta->sar->azimuth_time_per_pixel =
                    -inMeta->sar->azimuth_time_per_pixel;
            }

            time_start = time_max;
            time_end = time_min;

            onl = (time_end - time_start) / time_incr;
        }
        else {
            // here, no sar block in the original data, just make a square
            // image with decreasing time
            onl = ons;
            time_incr = (time_min - time_max) / (double)onl;
            time_start = time_max;
            time_end = time_min;
        }
    }
    else {
        // not provided a slant range pixel size, we'll figure something out

        if (inMeta->sar) {
            // use the same azimuth time per pixel.
            time_incr = inMeta->sar->azimuth_time_per_pixel;

            // we always want to be DECREASING in time
            // latest time is on top (line 1), earliest on bottom (line ONL)
            if (time_incr > 0) {
                time_incr = -time_incr;
                inMeta->sar->azimuth_time_per_pixel =
                    -inMeta->sar->azimuth_time_per_pixel;
            }

            time_start = time_max;
            time_end = time_min;

            onl = (time_end - time_start) / time_incr;
        }
        else {
            // no info... determine azimuth time per pixel by keeping
            // the height the same as in the original image
            onl = nl;
            time_incr = (time_min - time_max) / (double)onl;
            time_start = time_max;
            time_end = time_min;
        }

        // make it square, to get the slant range pixel size
        ons = onl;
        pixel_size = slant_incr = (slant_end - slant_start) / (double)ons;
    }

    asfRequire(onl > 0, "Internal Error: Invalid output line count: %d\n", onl);
    asfRequire(ons > 0, "Internal Error: Invalid output sample count: %d\n", ons);

    asfPrintStatus("  Slant range values: %f -> %f\n", slant_start, slant_end);
    asfPrintStatus("  Slant range pixel size: %f\n", pixel_size);
    asfPrintStatus("  Time values: %f -> %f\n", time_start, time_end);
    asfPrintStatus("  Output Image will be %5d x %5d LxS\n", onl, ons);
    asfPrintStatus("      (Input Image was %5d x %5d LxS)\n", nl, ns);

    // generate a grid over the image, to generate our splines
    // this grid size seems to work pretty well...
    int n = 120;

    asfPrintStatus("Creating %dx%d mapping grid...\n", n, n);

    double time_grid_incr = ((double)(onl))/((double)(n)) * time_incr;
    double slant_grid_incr = ((double)(ons))/((double)(n)) * slant_incr;

    // allocating memory for the splines, and the arrays to generate them
    gsl_interp_accel **samp_accels = MALLOC(sizeof(gsl_interp_accel *) * n);
    gsl_spline **samp_splines = MALLOC(sizeof(gsl_spline *) * n);

    gsl_interp_accel **line_accels = MALLOC(sizeof(gsl_interp_accel *) * n);
    gsl_spline **line_splines = MALLOC(sizeof(gsl_spline *) * n);

    double *slant_in = MALLOC(sizeof(double)*n);
    double *line_out = MALLOC(sizeof(double)*n);
    double *samp_out = MALLOC(sizeof(double)*n);

    // an alias -- use the same array (to save memory -- these are not used
    // at the same time), but create an alias for it, so it is not so confusing
    double *time_in = slant_in;
    //double max_err = 0;

    // set up the vertical splines
    for (jj=0; jj<n; ++jj) {
        double slant = slant_start + jj*slant_grid_incr;

        for (ii=0; ii<n; ++ii) {
            // splines need strictly increasing range variables
            if (time_grid_incr > 0)
                time_in[ii] = time_start + ii*time_grid_incr;
            else
                time_in[ii] = time_end - ii*time_grid_incr;

            ts2ls(inMeta, time_in[ii], slant, &line_out[ii], &samp_out[ii]);
        }

        samp_accels[jj] = gsl_interp_accel_alloc();
        samp_splines[jj] = gsl_spline_alloc(gsl_interp_cspline, n);
        gsl_spline_init(samp_splines[jj], time_in, samp_out, n);

        line_accels[jj] = gsl_interp_accel_alloc();
        line_splines[jj] = gsl_spline_alloc(gsl_interp_cspline, n);
        gsl_spline_init(line_splines[jj], time_in, line_out, n);
    }

    // now, we're on to the resampling stage.. loop through output pixels
    asfPrintStatus("Generating slant range image...\n");

    double no_data_value = meta_is_valid_double(inMeta->general->no_data) ?
        inMeta->general->no_data : 0;

    // keep track of error sizes
    double max_error = 0;
    double avg_error = 0;
    int count = 0;

    // these stride values allow us to track when we're in between grid points
    int ii_n = onl/n;
    int jj_n = ons/n;
    int ii_n2 = ii_n/2;
    int jj_n2 = jj_n/2;

    // set up output metadata
    meta_parameters *outMeta = meta_read(infile);

    if (outMeta->transform) {
        FREE(outMeta->transform);
        outMeta->transform = NULL;
    }

    if (outMeta->projection) {
        FREE(outMeta->projection);
        outMeta->projection = NULL;
    }

    outMeta->general->line_count = onl;
    outMeta->general->sample_count = ons;

    if (!outMeta->sar)
        outMeta->sar = meta_sar_init();

    outMeta->sar->image_type = 'S';

    outMeta->sar->azimuth_time_per_pixel = time_incr;
    assert(outMeta->sar->azimuth_time_per_pixel < 0);

    outMeta->sar->time_shift = time_start;

    outMeta->general->y_pixel_size = inMeta->sar->azimuth_time_per_pixel /
                                     time_incr * inMeta->general->y_pixel_size;
    assert(outMeta->general->y_pixel_size > 0);

    outMeta->sar->slant_range_first_pixel = slant_start;
    outMeta->general->x_pixel_size = slant_incr;

    outMeta->sar->line_increment = outMeta->sar->sample_increment = 1;
    outMeta->general->start_sample = outMeta->general->start_line = 0;

    outMeta->general->no_data = no_data_value;

    char **band_name = extract_band_names(inMeta->general->bands,
                                          inMeta->general->band_count);

    // now generate output image
    char *img_file = appendExt(outfile, ".img");
    float *out = MALLOC(sizeof(float) * ons);

    for (kk=0; kk<inMeta->general->band_count; ++kk) {
        if (inMeta->general->band_count != 1)
            asfPrintStatus("Working on band: %s\n", band_name[kk]);

        // for the 2nd and higher bands, free the band from the previous iteration,
        // and read in the next band from the input image
        if (kk>0) {
            float_image_free(in);
            asfPrintStatus("Loading input...\n");
            in = float_image_band_new_from_metadata(inMeta, kk, infile);
        }

        FILE *ofp = FOPEN(img_file, kk==0 ? "wb" : "ab");
        asfPrintStatus("Generating output...\n");
        for (ii=0; ii<onl; ++ii) {
            asfLineMeter(ii,onl);
            double time = time_start + ii * time_incr;

            // set up horizontal splines for this row
            gsl_interp_accel *samp_accel = gsl_interp_accel_alloc();
            gsl_spline *samp_spline = gsl_spline_alloc(gsl_interp_cspline, n);

            gsl_interp_accel *line_accel = gsl_interp_accel_alloc();
            gsl_spline *line_spline = gsl_spline_alloc(gsl_interp_cspline, n);

            //printf("time: %f slant: %f\n", time, slant_start);
            for (jj=0; jj<n; ++jj) {
                slant_in[jj] = slant_start + jj * slant_grid_incr;
                //printf("time: %f slant: %f\n", time, slant_in[jj]);
                samp_out[jj] = gsl_spline_eval(samp_splines[jj], time,
                                               samp_accels[jj]);
                line_out[jj] = gsl_spline_eval(line_splines[jj], time,
                                               line_accels[jj]);
                //printf("samp_out: %f line_out: %f\n", samp_out[jj], line_out[jj]);
            }

            gsl_spline_init(samp_spline, slant_in, samp_out, n);
            gsl_spline_init(line_spline, slant_in, line_out, n);

            // use the splines to produce output pixels
            for (jj=0; jj<ons; ++jj) {
                double slant = slant_start + jj * slant_incr;
                double samp = gsl_spline_eval(samp_spline, slant, samp_accel);
                double line = gsl_spline_eval(line_spline, slant, line_accel);

                // check the spline every so often (halfway between grid points)
                // only do this on band #1 (the reference band)
                if (kk==0 && ii%ii_n2==0 &&
                    ii%ii_n!=0 && jj%jj_n2==0 && jj%jj_n!=0)
                {
                    double samp_real, line_real;
                    ts2ls(inMeta, time, slant, &line_real, &samp_real);

                    double err = (line-line_real)*(line-line_real) +
                                 (samp-samp_real)*(samp-samp_real);

                    //printf("(%d,%d) -- Actual: (%f,%f) Splined: (%f,%f)\n",
                    //       ii, jj, line_real, samp_real, line, samp);

                    if (err > max_error) max_error = err;
                    avg_error += err;
                    ++count;
                }

                // now interpolate within the original image
                // if we are outside, use "no_data" from metadata
                double val = no_data_value;
                if (line > 0 && line < nl-1 && samp > 0 && samp < ns-1)
                    val = float_image_sample(in, samp, line, sampling_method);

                out[jj] = (float)val;
            }

            gsl_interp_accel_free(samp_accel);
            gsl_spline_free(samp_spline);

            gsl_interp_accel_free(line_accel);
            gsl_spline_free(line_spline);

            put_float_line(ofp, outMeta, ii, out);
        }

        fclose(ofp);
    }

    // free the last band of the input
    float_image_free(in);

    FREE(slant_in);
    FREE(line_out);
    FREE(samp_out);

    for (ii=0; ii<n; ++ii) {
        gsl_interp_accel_free(samp_accels[ii]);
        gsl_spline_free(samp_splines[ii]);
        gsl_interp_accel_free(line_accels[ii]);
        gsl_spline_free(line_splines[ii]);
    }

    FREE(samp_accels);
    FREE(samp_splines);
    FREE(line_accels);
    FREE(line_splines);

    FREE(out);

    for (kk=0; kk<inMeta->general->band_count; ++kk)
        FREE(band_name[kk]);
    FREE(band_name);

    // see how bad our errors were
    avg_error /= (double)count;
    asfPrintStatus("Model max error: %f, avg: %f\n",
                   max_error, avg_error);

    double thresh = 0.1;
    if (max_error > 100*thresh)
        asfPrintError("Maximum error exceeded threshold: %f > %f\n",
                      max_error, 100*thresh);
    else if (avg_error > 10*thresh)
        asfPrintError("Average error exceeded threshold: %f > %f\n",
                      avg_error, 10*thresh);
    if (max_error > 10*thresh)
        asfPrintWarning("Maximum error exceeds threshold: %f > %f\n",
                        max_error, 10*thresh);
    if (avg_error > thresh)
        asfPrintWarning("Average error exceeds threshold: %f > %f\n",
                        avg_error, thresh);

    char *meta_file = appendExt(outfile, ".meta");
    asfPrintStatus("Writing %s\n", meta_file);
    meta_write(outMeta, meta_file);
    free(meta_file);
    free(img_file);

    meta_free(outMeta);
    meta_free(inMeta);

    return 0; //success
}

int to_sr_pixsiz(const char *infile, const char *outfile, double pixel_size)
{
   meta_parameters *inMeta = meta_read(infile);
   int ret;

   if (inMeta->sar && inMeta->sar->image_type == 'G') {
       // ground range image
       ret = gr2sr_pixsiz(infile, outfile, pixel_size);
   }
   else if (inMeta->sar && inMeta->sar->image_type == 'S') {
       // already a slant range image, just copy it
       copyImgAndMeta(infile, outfile);
       ret = 0; // success
   }
   else if ((inMeta->sar && inMeta->sar->image_type == 'P') ||
            (inMeta->projection)) {
       // projected image
       ret = proj_to_sr(infile, outfile, pixel_size);
   }
   else if (inMeta->sar && inMeta->transform &&
            inMeta->sar->image_type == 'R') {
       // georeferenced, ALOS most likely
       ret = proj_to_sr(infile, outfile, pixel_size);
   }
   else {
       asfPrintError("Couldn't figure out what kind of image this is.\n");
       ret = 1;
   }

   meta_free(inMeta);
   return ret;
}

int to_sr(const char *infile, const char *outfile)
{
    char *img_in = MALLOC(sizeof(char)*(strlen(infile)+10));
    char *img_out = MALLOC(sizeof(char)*(strlen(outfile)+10));
    create_name(img_in, infile, ".img");
    create_name(img_out, outfile, ".img");
    int ret = to_sr_pixsiz(img_in, img_out, -1);
    FREE(img_in);
    FREE(img_out);
    return ret;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#pragma once

#include "icode_generator.h"
#include <halley/data_structures/hash_map.h>
#include <gsl/gsl>
#include "halley/file/path.h"

namespace YAML
{
	class Node;
}

namespace Halley
{
	class ECSData;
	class ComponentSchema;
	class SystemSchema;
	class MessageSchema;
	class CustomTypeSchema;

	class Codegen
	{
		struct Stats
		{
			int written = 0;
			int skipped = 0;
			Vector<Path> files;
		};

	public:
		constexpr static int currentCodegenVersion = 105;

		using ProgressReporter = std::function<bool(float, String)>;

		static void run(Path inDir, Path outDir);
		static Vector<Path> generateCode(const ECSData& data, Path directory);

	private:
		static bool writeFile(const Path& path, gsl::span<const char> data, bool stub);
		static void writeFiles(const Path& directory, const CodeGenResult& files, Stats& stats);
		static int getHeaderVersion(gsl::span<const char> data, size_t& endOfHeader);
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


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

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_interp.h>

#ifndef HIDE_INLINE_STATIC
size_t
gsl_interp_bsearch (
  const double x_array[], double x,
  size_t index_lo,
  size_t index_hi
  )
{
  size_t ilo = index_lo;
  size_t ihi = index_hi;
  while (ihi > ilo + 1)
    {
      size_t i = (ihi + ilo) / 2;
      if (x_array[i] > x)
        ihi = i;
      else
        ilo = i;
    }

  return ilo;
}
#endif