text:
//cmc_topol.c --
// this uses stateSpace generated previously to
// create the topology matrix and the move type matrix
//defined by the markov chain

//this was only for testing purposes

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>

#include "AFS.h"
#include "adkGSL.h"
#include "cs.h"
#include "time.h"
#include "adkCSparse.h"
#include "AFS_pthreads.c"

#define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
#define NTHREADS 4

int n1, n2;
const char *filename;

afsStateSpace *stateSpace;
void getParameters(int argc,const char **argv);
void usage();

int main(int argc, const char * argv[]){
	int i,N,tmpS,win,t;
	int rc;
	void *v;
	int nnz;
	matrixThreadObject_sparse *tObj[NTHREADS];
 	pthread_t threads[NTHREADS];

	getParameters(argc,argv);
	stateSpace = afsStateSpaceImportFromFile(filename);
	//afsStateSpaceRemoveAbsorbing(stateSpace);
	N = stateSpace->nstates;

	//Multithreaded construction of original transition matrix
//	printf("initializing original transition matrix using %d threads....\n",NTHREADS);
	//setup threads
	tmpS=0;
	win =(int) ceil( (((float)stateSpace->nstates) / NTHREADS) );
//	printf("nstates: %d winSize: %d\n",stateSpace->nstates,win);
	for(t=0;t<NTHREADS;t++){
		tObj[t] = malloc(sizeof( matrixThreadObject_sparse));
		tObj[t]->nnz = 0;
		tObj[t]->dim1 = malloc(sizeof(int) * stateSpace->nstates * 3);
		tObj[t]->dim2 = malloc(sizeof(int) * stateSpace->nstates * 3);
		tObj[t]->moveType = malloc(sizeof(int) * stateSpace->nstates *  3);
		tObj[t]->topol = malloc(sizeof(double) * stateSpace->nstates *  3);
		tObj[t]->start = tmpS;
		tObj[t]->stop = MIN(tmpS+win,stateSpace->nstates);
		tObj[t]->im_params = malloc(sizeof(im_lik_params));
		tObj[t]->im_params->stateSpace = stateSpace;
//		printf("start: %d stop: %d\n",tmpS,tmpS+win);
		rc = pthread_create(&threads[t], NULL, threadDoRows_sparse, (void *)tObj[t]);
		tmpS += win;
	}
//	printf("threads alloced....\n");

	//harvest jobs
	nnz=0;
	for(t=0;t<NTHREADS;t++){
		rc = pthread_join(threads[t], &v);
		tObj[t] = ((matrixThreadObject_sparse *)v);
		nnz += tObj[t]->nnz;
	}
	//now combine outputs for nnz, dim1, and dim2
	printf("nnz: %d\n",nnz);
	for(t=0;t<NTHREADS;t++){
		for(i=0;i<tObj[t]->nnz;i++){
			printf("%lf\t%d\t%d\t%d\n",tObj[t]->topol[i],tObj[t]->moveType[i],tObj[t]->dim1[i],tObj[t]->dim2[i]);
		}
	}
	//currentParams->nnz -= 1;
	////////// Multithreading done!

	// dim1 = malloc(sizeof(int) * N * 10);
	// dim2 = malloc(sizeof(int) * N * 10);
	// moveA = malloc(sizeof(int) * N * 10);
	// topA = malloc(sizeof(double) * N * 10);
	// for(i=0;i< (N*10);i++){
	// 	moveA[i]=0;
	// 	topA[i]=0;
	// }
	// nnz = coalMarkovChainTopologyMatrix_sparse(stateSpace,topA, moveA, dim1, dim2);
	// printf("nnz: %d\n",nnz);
	// for(i=0;i<nnz;i++){
	// 	printf("%f\t%d\t%d\t%d\n",topA[i],moveA[i],dim1[i],dim2[i]);
	// }
	return(0);
}

void getParameters(int argc, const char **argv){
//	int args;
//	int i;

	if( argc <= 1){
		usage();
	}
	filename = argv[1];

}

void usage(){
	fprintf(stderr,"usage: cmc_topol stateSpaceFile\n");
//	fprintf(stderr,"parameters: \n");
	exit(1);
}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
/* Copyright 2013 Perttu Luukko

 * This file is part of libeemd.

 * libeemd 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.

 * libeemd 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 libeemd.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _EEMD_H_
#define _EEMD_H_

#ifndef EEMD_DEBUG
#define EEMD_DEBUG 0
#endif

#if EEMD_DEBUG == 0
#ifndef NDEBUG
#define NDEBUG
#endif
#endif

#include <assert.h>
#include <limits.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_poly.h>

#ifdef _OPENMP
#include <omp.h>
#endif

// Possible error codes returned by functions eemd, ceemdan and
// emd_evaluate_spline
typedef enum {
	EMD_SUCCESS = 0,
	// Errors from invalid parameters
	EMD_INVALID_ENSEMBLE_SIZE = 1,
	EMD_INVALID_NOISE_STRENGTH = 2,
	EMD_NOISE_ADDED_TO_EMD = 3,
	EMD_NO_NOISE_ADDED_TO_EEMD = 4,
	EMD_NO_CONVERGENCE_POSSIBLE = 5,
	EMD_NOT_ENOUGH_POINTS_FOR_SPLINE = 6,
	EMD_INVALID_SPLINE_POINTS = 7,
	// Other errors
	EMD_GSL_ERROR = 8
} libeemd_error_code;

// Helper functions to print an error message if an error occured
void emd_report_if_error(libeemd_error_code err);
void emd_report_to_file_if_error(FILE* file, libeemd_error_code err);

// Main EEMD decomposition routine as described in:
//   Z. Wu and N. Huang,
//   Ensemble Empirical Mode Decomposition: A Noise-Assisted Data Analysis
//   Method, Advances in Adaptive Data Analysis,
//   Vol. 1, No. 1 (2009) 1–41
//
// Parameters 'input' and 'N' denote the input data and its length,
// respectively. Output from the routine is written to array 'output', which
// needs to be able to store at least N*M doubles, where M is the number of
// Intrinsic Mode Functions (IMFs) to compute. If M is set to zero, a value of
// M = emd_num_imfs(N) will be used, which corresponds to a maximal number of
// IMFs. Note that the final residual is also counted as an IMF in this
// respect, so you most likely want at least num_imfs=2. The following
// parameters are the ensemble size and the relative noise standard deviation,
// respectively. These are followed by the parameters for the stopping
// criterion. The stopping parameter can be defined by a S-number (see the
// article for details) or a fixed number of siftings. If both are specified,
// the sifting ends when either criterion is fulfilled. The final parameter is
// the seed given to the random number generator. A value of zero denotes a
// RNG-specific default value.
libeemd_error_code eemd(double const* restrict input, size_t N,
		double* restrict output, size_t M,
		unsigned int ensemble_size, double noise_strength, unsigned int
		S_number, unsigned int num_siftings, unsigned long int rng_seed);

// A complete variant of EEMD as described in:
//   M. Torres et al,
//   A Complete Ensemble Empirical Mode Decomposition with Adaptive Noise
//   IEEE Int. Conf. on Acoust., Speech and Signal Proc. ICASSP-11,
//   (2011) 4144-4147
//
// Parameters are identical to routine eemd
libeemd_error_code ceemdan(double const* restrict input, size_t N,
		double* restrict output, size_t M,
		unsigned int ensemble_size, double noise_strength, unsigned int
		S_number, unsigned int num_siftings, unsigned long int rng_seed);

// A method for finding the local minima and maxima from input data specified
// with parameters x and N. The memory for storing the coordinates of the
// extrema and their number are passed as the rest of the parameters. The
// arrays for the coordinates must be at least size N. The method also counts
// the number of zero crossings in the data, and saves the results into the
// pointer given as num_zero_crossings_ptr.
void emd_find_extrema(double const* restrict x, size_t N,
		double* restrict maxx, double* restrict maxy, size_t* num_max_ptr,
		double* restrict minx, double* restrict miny, size_t* num_min_ptr,
		size_t* num_zero_crossings_ptr);

// Return the number of IMFs that can be extracted from input data of length N,
// including the final residual.
size_t emd_num_imfs(size_t N);

// This routine evaluates a cubic spline with nodes defined by the arrays x and
// y, each of length N. The spline is evaluated using the not-a-node end point
// conditions (same as Matlab). The y values of the spline curve will be
// evaluated at integer points from 0 to x[N-1], and these y values will be
// written to the array spline_y. The endpoint x[N-1] is assumed to be an
// integer, and the x values are assumed to be in ascending order, with x[0]
// equal to 0. The workspace required is 5*N-10 doubles, except that N==2
// requires no extra memory. For N<=3 the routine falls back to polynomial
// interpolation, same as Matlab.
//
// This routine is mainly exported so that it can be tested separately to
// produce identical results to the Matlab routine 'spline'.
libeemd_error_code emd_evaluate_spline(double const* restrict x, double const* restrict y,
		size_t N, double* restrict spline_y, double* spline_workspace);

#endif // _EEMD_H_

ext: h

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
#ifndef __calcAinvb__
#define __calcAinvb__
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_multiroots.h>
#define GSL 2
int GSLtag;
void least_gsl_error_handler( const char * reason,
			      const char * file,
			      int line,
			      int gsl_errno)
{
  fprintf(stderr,"%s:%d: %s (gsl_error %d)\n",  file, line, reason,gsl_errno);
  GSLtag=1;
}

#define ZERO (1e-10)
double VarY;
int calcAinvb1(double *x_data, double *A_data, double *b_data, int nx, int ndata)
{
  /*
    b= A x  --> x= Ainv b
     b  in R^{ndata x 1}
     A  in R^{ndata x nx}   (ndata>nx)
     x  in R^{nx    x 1}
   */
  int i;

  if(ndata<=0){
    x_data[0]=1;
    for(i=1;i<nx;i++) x_data[i]=1./nx;
  }
  else{
    gsl_matrix *A =gsl_matrix_alloc(ndata,nx);
    gsl_vector *b =gsl_vector_alloc(ndata);
    gsl_matrix *V = gsl_matrix_alloc (nx,nx);
    gsl_vector *S = gsl_vector_alloc (nx);
    gsl_vector *work = gsl_vector_alloc(nx);
    gsl_vector *x = gsl_vector_alloc(nx);
    A->data=A_data;
    b->data=b_data;
    x->data=x_data;
    gsl_set_error_handler( &least_gsl_error_handler );GSLtag=0;

#if GSL == 1
    gsl_linalg_SV_decomp_jacobi (A, V, S); //higer accuracy than Golub-Reinsh
#elif GSL == 2
    gsl_linalg_SV_decomp (A, V, S, work); //Golub-Reinsh
#elif GSL == 3
    {//faster
      gsl_matrix *X;
      X = gsl_matrix_alloc (nx,nx);
      gsl_linalg_SV_decomp_mod (A, X, V, S, work);
      gsl_matrix_free(X);
    }
#endif
    VarY=0;
    if(GSLtag==1){
      fprintf(stderr,"\nError:EigenValue=");
      for(i=0;i<nx;i++) fprintf(stderr,"%e ",S->data[i]); fprintf(stderr,"\n");
      for(i=0;i<nx;i++) x_data[i]=1./nx;
    }
    else{
      int i,j,l;
      double *Si=(double*)malloc(sizeof(double)*nx);
      double *Ai=(double*)malloc(sizeof(double)*nx*ndata);
      for(i=0;i<nx;i++) if(S->data[i]>1e-10) Si[i]=1./S->data[i]; else Si[i]=0;
      for(i=0;i<nx;i++){
	for(j=0;j<ndata;j++){
	  double *Ai_ij=&Ai[i*ndata+j];
	  Ai_ij[0]=0;
	  for(l=0;l<nx;l++) Ai_ij[0]+=V->data[i*nx+l]*Si[l]*A->data[j*nx+l];
	}
      }
      //      for(i=0;i<nx;i++) x->data[i]=0;
      x->data[0]=0;for(j=0;j<nx;j++) x->data[0]+=Ai[j]*b->data[j];
//      int check=0;
//      if(check){
//	gsl_linalg_SV_solve(A,V,S,b,x);
//	double xi;
//	for(i=0;i<nx;i++){
//	  xi=0;
//	  for(j=0;j<nx;j++) xi+= Ai[i*ndata+j]*b->data[j];
//	  fprintf(stderr,"x[%d]=%g =?= %g\n",i,xi,x->data[i]);
//	}
//	fprintf(stderr,"p:");
//	for(i=0;i<nx;i++){
//	  fprintf(stderr,"\n p(i%d,*):",i);
//	  for(j=0;j<nx;j++){
//	    fprintf(stderr,"%+.3f ",Ai[i*nx+j]);
//	  }
//	}
//	fprintf(stderr,"\n:");
//	fprintf(stderr,"check\n");
//      }
      free(Si);
      free(Ai);
    }
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(x);
    gsl_matrix_free(V);
    gsl_matrix_free(A);
    gsl_vector_free(b);
  }
  return 0;
}
int calc_Ainvb(double *x_data, double *A_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /*
    b= A x  --> x= Ainv b
     b  in R^{ndata x 1}
     A  in R^{ndata x nx}   (ndata>nx)
     x  in R^{nx    x 1}
   */
  int i;

  if(ndata<=0){
    x_data[0]=1;
    for(i=1;i<nx;i++) x_data[i]=1./nx;
  }
  else{
    gsl_matrix *A =gsl_matrix_alloc(ndata,nx);
    gsl_vector *b =gsl_vector_alloc(ndata);
    //A->data=A_data;
    //	b->data=b_data;
    {
      int n=ndata*nx;
	for(i=0;i<n;i++){
	  A->data[i]=A_data[i];
	}
	for(i=0;i<ndata;i++) b->data[i]=b_data[i];
    }

    gsl_matrix *V = gsl_matrix_alloc (nx,nx);
    gsl_vector *S = gsl_vector_alloc (nx);
    gsl_vector *work = gsl_vector_alloc(nx);
    gsl_vector *x = gsl_vector_alloc(nx);
    gsl_set_error_handler( &least_gsl_error_handler );GSLtag=0;

#if GSL == 1
    gsl_linalg_SV_decomp_jacobi (A, V, S); //higer accuracy than Golub-Reinsh
#elif GSL == 2
    gsl_linalg_SV_decomp (A, V, S, work); //Golub-Reinsh
#elif GSL == 3
    {//faster
      gsl_matrix *X;
      X = gsl_matrix_alloc (nx,nx);
      gsl_linalg_SV_decomp_mod (A, X, V, S, work);
      gsl_matrix_free(X);
    }
#endif
    VarY=0;
    if(GSLtag==1){
      fprintf(stderr,"\nError:EigenValue=");for(i=0;i<nx;i++) fprintf(stderr,"%e ",S->data[i]); fprintf(stderr,"\n");
      for(i=0;i<nx;i++) x_data[i]=1./nx;
    }
    else{
      if(nxmarume>=0){//Principal Component Analysis
	for(i=0;i<nxmarume;i++) VarY+=S->data[i];
	for(i=nxmarume;i<nx;i++) S->data[i]=0;
	gsl_linalg_SV_solve(A,V,S,b,x);
	for(i=0;i<nx;i++) x_data[i]=gsl_vector_get(x,i);
      }
      else{//Minor Component Analysis
	int ii=0;
	int im=-nxmarume;
	for(i=nx-1;i>=0;i--){
	  if(S->data[i]>1e-20){
	    ii++;
	    VarY+=S->data[i];
	    if(ii>=im) break;
	  }
	}
	int j;

	double M0=-V->data[i];
	for(j=1;j<nx;j++){
	  x_data[j-1]=V->data[j*nx+i]/M0;
	}
      }
    }
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(x);
    gsl_matrix_free(V);
    gsl_matrix_free(A);
    gsl_vector_free(b);
  }
  //  fprintf(stderr,"OK");
  return 0;
}
void calc_AtinvbMCA(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  int nx1=nx+1,i,j,jnx1;
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx1);
  for(j=0;j<ndata;j++){
    jnx1=j*nx1;
    A->data[jnx1]=b_data[j];
    for(i=0;i<nx;i++){
      A->data[jnx1+i+1]=At_data[i*ndata+j];
    }
  }
  calc_Ainvb(M, (double *)(A->data), b_data, nx1, ndata, -nxmarume);

  gsl_matrix_free(A);
}
void calc_Atinvb(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  gsl_matrix *At =gsl_matrix_alloc(nx,ndata);
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx);
  At->data=At_data;
  gsl_matrix_transpose_memcpy(A, At);
  calc_Ainvb(M, (double *)(A->data), b_data, nx, ndata, nxmarume);
  gsl_matrix_free(A);
  gsl_matrix_free(At);
}

#endif // __calcAinvb__

ext: c


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

int main()
{
	double bottom_tail = gsl_cdf_gaussian_P(-1.96,1);
	printf("Area between +-1.96: %g\n",1-2*bottom_tail);

	return 0;
}

max_stars_repo_name: jieyaren/hello-world

max_issues_repo_name: jieyaren/hello-world

max_forks_repo_name: jieyaren/hello-world


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef pool_81b70dbf_d70f_42d7_8770_678651c585c4_h
#define pool_81b70dbf_d70f_42d7_8770_678651c585c4_h

#include <assert.h>
#include <memory.h>
#include <stdlib.h>
#include <gslib/type.h>

__gslib_begin__

template<class _construct>
class factory
{
public:
    typedef _construct mycon;
    typedef factory<mycon> myref;

public:
    factory(void* ptr = nullptr) { _generated = (mycon*)ptr; }
    const mycon* get_ptr() const { return _generated; }
    mycon* get_ptr() { return _generated; }
    mycon* emerge()
    {
        if(!_generated)
            _generated = new mycon;
        assert(_generated);
        return _generated;
    }
    void destroy()
    {
        if(!_generated)
            return;
        delete _generated;
        _generated = nullptr;
    }

protected:
    mycon*              _generated;
};

class vessel
{
protected:
    void*               _buf;
    int                 _cap;
    int                 _cur;

public:
    vessel()
    {
        _buf = nullptr;
        _cap = 0;
        _cur = 0;
    }
    ~vessel()
    {
        destroy();
    }
    void destroy()
    {
        if(_buf) {
            free(_buf);
            _buf = nullptr;
        }
        _cap = _cur = 0;
    }
    int get_cap() const
    {
        return _cap;
    }
    int get_cur() const
    {
        return _cur;
    }
    void* get_ptr() const
    {
        return _buf;
    }
    int get_rest() const
    {
        return _cap - _cur;
    }
    void flex(int size)
    {
        _buf = realloc(_buf, size);
        _cap = size;
    }
    void expand(int size)
    {
        flex(get_cap() + size);
    }
    void fit()
    {
        _buf = realloc(_buf, _cur);
        _cap = _cur;
    }
    void occupy(int size)
    {
        _cur += size;
        assert(_cur <= get_cap());
    }
    void set_cur(int size)
    {
        _cur = size;
        assert(_cur <= get_cap());
    }
    int curaddr(int ofs = 0) const
    {
        return (int)_buf + _cur + ofs;
    }
    void store(const void* buf, int size)
    {
        if(_cur + size > get_cap())
            expand(get_cap() > size ? get_cap() : size);
        assert(get_cap() >= _cur + size);
        memcpy_s((void*)curaddr(), get_rest(), buf, size);
        occupy(size);
    }
    void attach(void* buf, int size)
    {
        flex(0);
        _buf = buf;
        _cap = size;
        _cur = size;
    }
    void attach(vessel* vsl)
    {
        assert(vsl);
        attach(vsl->get_ptr(), vsl->get_cap());
        _cur = vsl->get_cur();
        vsl->detach();
    }
    void detach()
    {
        _buf = nullptr;
        _cap = 0;
        _cur = 0;
    }
    template<class tpl>
    tpl& current(int ofs = 0)
    {
        return *((tpl*)((int)_buf + _cur + ofs));
    }
    template<class tpl>
    tpl& front(int ofs = 0)
    {
        return *((tpl*)((int)_buf + ofs));
    }
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

#include <GL/glut.h>

#define SLICE	360

int *jac;

int rhs (double t, const double y[], double f[], void *params_ptr);
//int jacobian (double t, const double y[], double *dfdy, double dfdt[], void *params_ptr);

void draw_differential_eq_plot(void);
void draw_lattice(void);

typedef struct _xy_plot xy_plot;
struct _xy_plot
{
	double x;
	double y;
};

xy_plot plot_res[1001];

void display(void)
{
	glClearColor(0.0, 0.0, 0.0, 1.0);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();

	glColor3f(1, 0, 0);

	glBegin(GL_LINE_LOOP);
	glVertex3f(100.0, 0.0, 0.0);
        glVertex3f(-100.0, 0.0, 0.0);
        glEnd();

        glColor3f(0.0, 1.0, 0.0);

        glBegin(GL_LINE_LOOP);
        glVertex3f(0.0, 100.0, 0.0);
        glVertex3f(0.0, -100.0, 0.0);
    glEnd();

	draw_differential_eq_plot();
	draw_lattice();
	glutSwapBuffers();
}

void reshape(int w, int h)
{
    GLfloat n_range = 100.0f;

    if(h == 0)
        h = 1;

    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

#if 0
    if(w <= h)
        glOrtho(-n_range, n_range, -n_range * h / w, n_range * h / w, -n_range, n_range);
    else
        glOrtho(-n_range * w / h, n_range * w / h, -n_range, n_range, -n_range, n_range);
#endif

	glOrtho(-10, 10, -0.5, 0.5, -1, 1);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

int calc_vertices_num(void)
{
	double tmin = 0.0, tmax = 10.0, delta_t = 0.01;
	return (tmax - tmin) / delta_t;
}

void calc_differential_eq(void)
{
    int dimension = 2;		/* number of differential equations */

    double eps_abs = 1.e-8;	/* absolute error requested */
    double eps_rel = 1.e-10;	/* relative error requested */

    const gsl_odeiv_step_type *type_ptr = gsl_odeiv_step_rkf45;

    gsl_odeiv_step *step_ptr = gsl_odeiv_step_alloc (type_ptr, dimension);
    gsl_odeiv_control *control_ptr = gsl_odeiv_control_y_new (eps_abs, eps_rel);
    gsl_odeiv_evolve *evolve_ptr = gsl_odeiv_evolve_alloc (dimension);

    gsl_odeiv_system my_system;	/* structure with the rhs function, etc. */

    double mu = 10;		/* parameter for the diffeq */
    double y[2];			/* current solution vector */

    double t, t_next;		/* current and next independent variable */
    double tmin, tmax, delta_t;	/* range of t and step size for output */

    double h = 1e-6;		/* starting step size for ode solver */

	int cnt = 0;
	int vertex_num = calc_vertices_num();
#if 0
	double res_x[vertex_num + 1] = {0};
	double res_y[vertex_num + 1] = {0};

	double *res_x = (double *)malloc(sizeof(double) * (vertex_num + 1));
	double *res_y = (double *)malloc(sizeof(double) * (vertex_num + 1));
#endif

    my_system.function = rhs;	/* the right-hand-side functions dy[i]/dt */
    my_system.jacobian = NULL;
    my_system.dimension = dimension;	/* number of diffeq's */
    my_system.params = NULL;

    tmin = 0.;			/* starting t value */
    tmax = 10.;			/* final t value */
    delta_t = 0.01;

    y[0] = 0.16;			/* initial x value */
    y[1] = 0.;			/* initial v value */

    t = tmin;             /* initialize t */
    printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);	/* initial values */

    for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)
    {
        while (t < t_next)	/* evolve from t to t_next */
        {
            gsl_odeiv_evolve_apply (evolve_ptr, control_ptr, step_ptr,
                                    &my_system, &t, t_next, &h, y);
        }

		plot_res[cnt].x = t;
		plot_res[cnt].y = y[0];
        printf ("%.5e %.5e %.5e\n", t, plot_res[cnt].x, plot_res[cnt].y); /* print at t=t_next */

		cnt++;
    }

    gsl_odeiv_evolve_free (evolve_ptr);
    gsl_odeiv_control_free (control_ptr);
    gsl_odeiv_step_free (step_ptr);
}

void draw_lattice(void)
{
	int i;

	glColor3f(1, 1, 0);
	glBegin(GL_LINES);
		for (i = 1; i <= 10; i++)
		{
			glVertex2f(i, -0.01);
			glVertex2f(i, 0.01);
		}
	glEnd();
}

void draw_differential_eq_plot(void)
{
	float t = -100.0, step = 0.01;
	float x = 0, x2 = 0, y2, cx, cy;
	float tmp;
	int cache = 0;
	int i;

	glBegin(GL_LINES);
	//for(; ; t += step)
	for (i = 0; i < 1000; i++)
	{
		x2 = plot_res[i].x;
		y2 = plot_res[i].y;

		if(cache)
		{
			glVertex2f(cx, cy);	// 이전값
			glVertex2f(x2, y2);	// 현재값
		}

		cache = 1;
		cx = x2;
		cy = y2;
		printf("t = %f, y2 = %f\n", x2, y2);
	}
	glEnd();

}

int main (int argc, char **argv)
{
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE);
	glutInitWindowSize(800, 800);
	glutInitWindowPosition(0, 0);
	glutCreateWindow("Digital Signal Processing");

	calc_differential_eq();

	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutMainLoop();

    return 0;
}

/*************************** rhs ****************************/
/*
   x' = v                  ==>  dy[0]/dt = f[0] = y[1]
   v' = -x + \mu v (1-x^2) ==>  dy[1]/dt = f[1] = -y[0] + mu*y[1]*(1-y[0]*y[0])

   x''(t) + x'(t) + x(t) = 0
   x''(t) = -x'(t) - x(t)

   x' = v            ===> dy[0]/dt = f[0] = y[1]
   v' = -v - x       ===> dy[1]/dt = f[1] = -y[1] - y[0]
*/

// https://www.wolframalpha.com/input/?i=y%27%27+%2B+y%27+%2B+y+%3D+0%2C+y%280%29+%3D+1%2C+y%27%280%29+%3D+0, y(0) = 1
// https://www.wolframalpha.com/input/?i=%28sqrt%283%29+*+sin%28sqrt%283%29+%2F+2%29+%2B+3+*+cos%28sqrt%283%29+%2F+2%29%29+%2F+%283+*+sqrt%28e%29%29, y(1)
// (sqrt(3) * sin(sqrt(3) / 2) + 3 * cos(sqrt(3) / 2)) / (3 * sqrt(e)), y(1) = 0.6597001...
int
rhs (double t, const double y[], double f[], void *params_ptr)
{
	/*
    f[0] = y[1];
    f[1] = -y[0] + mu * y[1] * (1. - y[0] * y[0]);
	*/

	// 10 * y'' + 10 * y' + 90 * y = 0
	// y'' + y' + 9 * y = 0
	// y'' = - y' - 9 * y
	// 파이썬의 dsolve가 자체적으로 표준형을 만들어주는 반면
	// GSL의 경우 성능은 더 좋지만 이와 같은 작업을 사용자에게 위임하고 있음
	f[0] = y[1];
	f[1] = -y[1] - 9 * y[0];

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}

/*************************** Jacobian ****************************/
/*
*/

#if 0
int
jacobian (double t, const double y[], double *dfdy,
          double dfdt[], void *params_ptr)
{
    /* get parameter(s) from params_ptr; here, just a double */
    double mu = *(double *) params_ptr;

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);

    gsl_matrix *m_ptr = &dfdy_mat.matrix;	/* m_ptr points to the matrix */

    /* fill the Jacobian matrix as shown */
    gsl_matrix_set (m_ptr, 0, 0, 0.0);	/* df[0]/dy[0] = 0 */
    gsl_matrix_set (m_ptr, 0, 1, 1.0);	/* df[0]/dy[1] = 1 */
    gsl_matrix_set (m_ptr, 1, 0, -2.0 * mu * y[0] * y[1] - 1.0); /* df[1]/dy[0] */
    gsl_matrix_set (m_ptr, 1, 1, -mu * (y[0] * y[0] - 1.0));     /* df[1]/dy[1] */

    /* set explicit t dependence of f[i] */
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}
#endif


text:
#ifndef VKST_RENDERER_H
#define VKST_RENDERER_H

#include <plat/filesystem.h>
#include <wsi/window.h>
#include <vk/result.h>
#include <gsl.h>
#include <filesystem>
#include <vector>

// Holds all of the data for a surface. This includes the swapchain,
// renderpass, and framebuffers.
// Surfaces must be resized when the window is resized.
class surface {
public:
  std::size_t num_images() const noexcept { return _color_images.size(); }

  VkRenderPass render_pass() const noexcept { return _render_pass; }

  VkFramebuffer framebuffer(std::size_t index) const noexcept {
    return _framebuffers[index];
  }

  VkSampleCountFlagBits samples() const noexcept { return _samples; }

  VkViewport& viewport() noexcept { return _viewport; }
  VkViewport const& viewport() const noexcept { return _viewport; }

  VkRect2D& scissor() noexcept { return _scissor; }
  VkRect2D const& scissor() const noexcept { return _scissor; }

  surface() noexcept {}
  surface(surface const&) = delete;
  surface(surface&& other) noexcept;
  surface& operator=(surface const&) = delete;
  surface& operator=(surface&& rhs) noexcept;
  ~surface() noexcept = default;

private:
  VkSurfaceKHR _surface{VK_NULL_HANDLE};

  VkSurfaceFormatKHR _color_format{};
  VkFormat _depth_format{};
  VkSampleCountFlagBits _samples{};
  VkPresentModeKHR _present_mode{};

  VkSemaphore _image_available{VK_NULL_HANDLE};
  VkSemaphore _render_finished{VK_NULL_HANDLE};
  VkRenderPass _render_pass{VK_NULL_HANDLE};

  VkSurfaceCapabilitiesKHR _capabilities{};
  VkExtent2D _extent{};
  VkViewport _viewport{};
  VkRect2D _scissor{};
  VkSwapchainKHR _swapchain{VK_NULL_HANDLE};

  constexpr static uint32_t MAX_IMAGES = 4;
  std::vector<VkImage> _color_images{};
  std::vector<VkImageView> _color_image_views{};

  VkImage _depth_image{VK_NULL_HANDLE};
  VkDeviceMemory _depth_image_memory{VK_NULL_HANDLE};
  VkImageView _depth_image_view{VK_NULL_HANDLE};

  VkImage _color_target{VK_NULL_HANDLE};
  VkDeviceMemory _color_target_memory{VK_NULL_HANDLE};
  VkImageView _color_target_view{VK_NULL_HANDLE};

  VkImage _depth_target{VK_NULL_HANDLE};
  VkDeviceMemory _depth_target_memory{VK_NULL_HANDLE};
  VkImageView _depth_target_view{VK_NULL_HANDLE};

  std::vector<VkFramebuffer> _framebuffers{};

  friend class renderer;
}; // class surface

// Convenience class to hold both the VkShaderModule for a shader as well as
// any error message from compiling the shader.
class shader {
public:
  enum class types : uint8_t {
    vertex = 0,
    fragment = 1,
  }; // enum class types

  operator VkShaderModule() const noexcept { return _module; }

  std::string const& error_message() const noexcept { return _error_message; }

  shader() noexcept {}
  shader(shader const&) = delete;
  shader(shader&& other) noexcept;
  shader& operator=(shader const&) = delete;
  shader& operator=(shader&& rhs) noexcept;
  ~shader() noexcept = default;

private:
  VkShaderModule _module{VK_NULL_HANDLE};
  std::string _error_message{};

  friend class renderer;
}; // class shader

enum class renderer_result {
  success = 0,
  no_device = 1,
  initialization_failed = 2,
  surface_not_supported = 3,
  no_memory_type = 4,
}; // class renderer_result

class renderer_result_category_impl : public std::error_category {
public:
  virtual char const* name() const noexcept override {
    return "renderer_result";
  }

  virtual std::string message(int ev) const override;
}; // class renderer_result_category_impl

std::error_category const& renderer_result_category();

inline std::error_code make_error_code(renderer_result e) noexcept {
  return {static_cast<int>(e), renderer_result_category()};
}

enum class renderer_options : uint8_t {
  none = 0,
  use_integrated_gpu = (1 << 1),
}; // renderer_options

// Holds all of the data for rendering. Also provides methods for creating
// surfaces, shaders, pipelines, and command buffers, as well as submitting
// command buffers for execution on the device.
class renderer {
public:
  // Create a new renderer. If ec is true, then an error occurred during
  // creation and the renderer object is in an invalid state.
  static renderer create(gsl::czstring application_name, renderer_options opts,
                         PFN_vkDebugReportCallbackEXT debug_report_callback,
                         uint32_t push_constant_size,
                         std::error_code& ec) noexcept;

  // Create a new surface. If ec is true, then an error occurred during
  // creation and the surface object is in an invalid state.
  surface create_surface(wsi::window const& window,
                         std::error_code& ec) noexcept;

  // Resize a surface. Must be called when the window that was passed for
  // surface creation is resized. This is not automatically done to allow
  // the render loop to determine when to perform the resize. If ec is true,
  // then an error occurred during the resize.
  void resize(surface& s, wsi::extent2d const& extent,
              std::error_code& ec) noexcept;

  // Acquire the next ready image in the swapchain for rendering to. This
  // must be called and the returned index passed to submit_present. If
  // ec is true, then an error occurred and the returned index is invalid.
  uint32_t acquire_next_image(surface& s, std::error_code& ec) noexcept;

  // Submit a set of command buffers for execution and then present the
  // previously acquired swapchain image. image_index must come from an
  // immediately preceding call to acquire_next_image. fence can be
  // VK_NULL_HANDLE which indicates no fence should be signaled when the
  // command buffers can be reused. If ec is true, then an error occurred
  // during either the submit or present.
  void submit_present(gsl::span<VkCommandBuffer> buffers, surface& s,
                      uint32_t image_index, VkFence fence,
                      std::error_code& ec) noexcept;

  void destroy(surface& s) noexcept;

private:
  void release(surface& s) noexcept;

public:
  // Allocate a set of command buffers. If ec is true, then an error occurred
  // and the vector is invalid.
  std::vector<VkCommandBuffer>
  allocate_command_buffers(uint32_t count, std::error_code& ec) noexcept;

  // Submit a set of command buffers. onetime indicates that the submit should
  // use the onetime fence and wait for the submit to complete before
  // continuing. If ec is true, then an error occurred.
  void submit(gsl::span<VkCommandBuffer> command_buffers, bool onetime,
              std::error_code& ec) noexcept;

  void free(std::vector<VkCommandBuffer>& command_buffers) noexcept;

  // Create a new shader from the given source code. path is expected to hold
  // GLSL source code which will be compiled before creating the shader. If ec
  // is true, then an error occurred and the shader is valid such that
  // shader::error_message can be called to get any compilation errors.
  shader create_shader(plat::filesystem::path const& path, shader::types type,
                       std::error_code& ec) noexcept;

  void destroy(shader& s) noexcept;

  // Create a new pipeline layout. If ec is true, then an error occurred and
  // the pipeline layout object is invalid.
  VkPipelineLayout create_pipeline_layout(
    gsl::span<VkDescriptorSetLayout> descriptor_set_layouts,
    gsl::span<VkPushConstantRange> push_constant_ranges,
    std::error_code& ec) noexcept;

  void destroy(VkPipelineLayout layout) noexcept;

  // Create a new set of pipelines. A single pipeline can also be created.
  // If ec is true, then an error occurred and the vector of pipelines is
  // invalid.
  std::vector<VkPipeline>
  create_pipelines(gsl::span<VkGraphicsPipelineCreateInfo> cinfos,
                   std::error_code& ec) noexcept;

  void destroy(gsl::span<VkPipeline> pipes) noexcept;

  // Create a new fence. If ec is true, then an error occurred and the fence
  // is invalid.
  VkFence create_fence(bool signaled, std::error_code& ec) noexcept;

  // Wait for a set of fences to be signaled. If wait_all is true, then all
  // fences must be signaled before the call will return, if wait_all is
  // false, then a single fence will cause the call to return. timeout
  // specifies how long to wait, 0 returns immediately, and UINT64_MAX will
  // wait indefinitely. If ec is true, and error occurred while waiting.
  void wait(gsl::span<VkFence> fences, bool wait_all, uint64_t timeout,
            std::error_code& ec) noexcept;

  // Reset a set of fences from the signaled state to unsignaled. If ec is
  // true then an error occured.
  void reset(gsl::span<VkFence> fences, std::error_code& ec) noexcept;

  void destroy(VkFence fence) noexcept;

  constexpr renderer() noexcept {};
  renderer(renderer const&) = delete;
  renderer(renderer&& other) noexcept;
  renderer& operator=(renderer const&) = delete;
  renderer& operator=(renderer&& rhs) noexcept;
  ~renderer() noexcept;

private:
  VkInstance _instance{VK_NULL_HANDLE};
  VkDebugReportCallbackEXT _callback{VK_NULL_HANDLE};

  VkPhysicalDevice _physical{VK_NULL_HANDLE};
  VkDevice _device{VK_NULL_HANDLE};

  uint32_t _graphics_queue_family_index{UINT32_MAX};
  VkQueue _graphics_queue{VK_NULL_HANDLE};
  VkCommandPool _graphics_command_pool{VK_NULL_HANDLE};
  VkFence _graphics_onetime_fence{VK_NULL_HANDLE};
}; // class renderer

inline constexpr auto operator|(renderer_options a,
                                renderer_options b) noexcept {
  using U = std::underlying_type_t<renderer_options>;
  return static_cast<renderer_options>(static_cast<U>(a) | static_cast<U>(b));
}

inline constexpr auto operator&(renderer_options a,
                                renderer_options b) noexcept {
  using U = std::underlying_type_t<renderer_options>;
  return static_cast<renderer_options>(static_cast<U>(a) & static_cast<U>(b));
}

#endif // VKST_RENDERER_H

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst