text:
#include <stdio.h>
#include "sweeny_bfs.h"
#include "../src/queue_2.h"
#include <ctype.h>
#include <strings.h>
#include <unistd.h>
#include <math.h>
#include <gsl/gsl_rng.h>

// Define Macro SEQUENTIAL to use sequential BFS; Default is interleaved BFS

#define MIN(a,b) a<=b ? a: b
#define MAX(a,b) a>=b? a:b

static char setup=0; // 0 not setup; 1 ibfs; 2 sbfs
static char verbose=0;
static double rcWeights[4]; // array with precalculated mc weights
static double p_min_del,p_max_del, p_min_ins,p_max_ins;
static __s8 dN,dB;
static char equilibration=1;
static __u32 DX;
static __u32 N;
static __u32 seed;
static double q;
static double coupling;
static double beta;
static double v;
static double K;
static __u32 steps;
static __s8 *bonds;
static gsl_rng * r;
static __u32 edge[2];
static __s32 adj1[4];
static __s32 adj2[4];
static __u32 offset_1;
static __u32 offset_2;
static __u32 cs2;
static struct queue_2_node *todo_pool;
static __u32 *visited;
static __u32 cs1;
static __u32 activeEdges=0;
static __u32 cutoff;
static __u32 *num_bonds , *num_cluster, *size_giant;
static __u64 *sec_cs_moment, *four_cs_moment;

/******************************************************************************
 *****************************************************************************/
static char init(void)
{
    K=coupling*beta;
    v = exp(K) - 1.0;
    N = DX*DX;
    // Deletion accept. ratios
    rcWeights[0] = pow(v,-1); //db==-1,dN==0
    rcWeights[1] = rcWeights[0]*q; //db==-1,dN=1
    // Insertion accept. ratios
    rcWeights[2] = v; //db==1,dN==0
    rcWeights[3] = v*pow(q,-1); //db==1,dN==-1

    p_min_del = MIN(rcWeights[0],rcWeights[1]);
    p_max_del = MAX(rcWeights[0],rcWeights[1]);
    p_min_ins = MIN(rcWeights[2],rcWeights[3]);
    p_max_ins = MAX(rcWeights[2],rcWeights[3]);

    r = gsl_rng_alloc (gsl_rng_mt19937);
    if(!r)
        return 0;
    gsl_rng_set(r,seed);
    __u32 i;
    bonds = malloc(sizeof(*bonds)*2*N);
    if(!bonds) {
        return 0;
    }
    for(i=0;i<2*N;i++) bonds[i] = -1;
    visited = calloc(N,sizeof(*visited));
    if(!visited) {
        return 0;
    }
    todo_pool = malloc(N*sizeof(struct queue_2_node));
    if(!todo_pool) {
        return 0;
    }
    return 1;

}
/******************************************************************************
 *****************************************************************************/
static void destroy(void)
{
    if(!setup)return;
    free(todo_pool);
    free(bonds);
    gsl_rng_free(r);
    free(visited);
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcXnext(__u32 idx)
{
    if((idx+1) % DX)
        return idx+1;
    else
        return idx +1- DX ;
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcXprev(__u32 idx)
{
    if((idx%DX))
        return idx-1;
    else
        return idx + DX  -1;
}
/******************************************************************************
 *****************************************************************************/
 static inline __u32 ltcYnext(__u32 idx)
{

    if(idx <= N -1 && idx >= N - DX)
        return idx +DX- N ;
    else
        return idx+DX;
}
/******************************************************************************
 *****************************************************************************/
static inline __u32 ltcYprev(__u32 idx)
{
    if(idx <= DX - 1)
        return idx + N - DX;
    else
        return idx - DX;

}
/******************************************************************************
 *****************************************************************************/
static void Adjacent(__u32 bidx)
{

  if(bidx%2) {
    bidx = (bidx -1)/2;
    edge[0] = bidx;
    edge[1] = ltcXnext(bidx);
  }
  else {
    bidx = bidx/2;
    edge[0] = bidx;
    edge[1] = ltcYprev(bidx);
  }

}
/******************************************************************************
 *****************************************************************************/
static void neighbours(__u32 idx,__u8 a) {
  if(a == 1)
  {
    if(bonds[2*idx] == 1)  adj1[0] = ltcYprev(idx);
    else  adj1[0] = -1;
    if(bonds[(2*idx)+1] == 1)  adj1[1] = ltcXnext(idx);
    else adj1[1] = -1;
    if(bonds[2*ltcYnext(idx)] == 1)	adj1[2] = ltcYnext(idx);
    else adj1[2] = -1;
    if(bonds[(2*ltcXprev(idx))+1] ==1)adj1[3] = ltcXprev(idx);
    else adj1[3] = -1;
  }
  else {
    if(bonds[2*idx] == 1)adj2[0] = ltcYprev(idx);
    else adj2[0] = -1;
    if(bonds[(2*idx)+1] == 1)adj2[1] = ltcXnext(idx);
    else adj2[1] = -1;
    if(bonds[2*ltcYnext(idx)] == 1)	adj2[2] = ltcYnext(idx);
    else adj2[2] = -1;
    if(bonds[(2*ltcXprev(idx))+1] == 1)adj2[3] = ltcXprev(idx);
    else adj2[3] = -1;
  }

}
/******************************************************************************
 *****************************************************************************/
static __u8 breadthFirstSearch_s(__u32 start, __u32 goal)
{
  static struct queue_2 todo1;
  __u32 i=0,activeP1=0;
  init_queue_2(&todo1); //todo_pool);
  cs1=0;
  enqueue_2(&todo1,start);
  cs1++; // Increase cluster size of bfs 1
  visited[start] = offset_1; // mark starting point as visited
    while(!queue_2_empty_p(&todo1)) {
      dequeue_2(&todo1,&activeP1);
      neighbours(activeP1,1);
      for(i=0;i<4;i++) {
        if(adj1[i] != -1) {
            if(visited[adj1[i]] == offset_1) continue;
            if((__u32)adj1[i] == goal) {
                    while(!queue_2_empty_p(&todo1)){ dequeue_2(&todo1,&activeP1);}
                    return 1;
                }
            enqueue_2(&todo1,adj1[i]);
            visited[adj1[i]] = offset_1;
            cs1++;
        }
      }
  }
  return 0;
}
/******************************************************************************
 *****************************************************************************/
static __u8 breadthFirstSearch(__u32 start, __u32 goal)
{
  static struct queue_2 todo1,todo2;
  __u32 i=0,activeP1=0,activeP2=0;
  init_queue_2(&todo1);
  init_queue_2(&todo2);
  cs1=0;
  cs2=0;
  enqueue_2(&todo1,start); // Put starting point onto queue_2 1
  cs1++; // Increase cluster size of bfs 1
  visited[start] = offset_1; // mark starting point as visited
  enqueue_2(&todo2,goal); // Put goal point onto queue_2 2
  visited[goal] = offset_2; // mark goal point as visited
  cs2++; // increase cluster size of bfs 2
    while(!queue_2_empty_p(&todo1) && !queue_2_empty_p(&todo2)) {
      dequeue_2(&todo2,&activeP2); // get next vertex of bfs 2
      neighbours(activeP2,2);	 // get all adjacent vertices of current vertex
      for(i=0;i<4;i++) {

        if(adj2[i] != -1) { // if accessable (edge exists)
	  if(visited[adj2[i]] == offset_2) continue; // already visited
	  if((__u32)adj2[i] == start || visited[adj2[i]] == offset_1) { // reconnected

            while(!queue_2_empty_p(&todo1)){ dequeue_2(&todo1,&activeP1);}
            while(!queue_2_empty_p(&todo2)){ dequeue_2(&todo2,&activeP2);}
	    return 1; // 1 indicates success
	  }
	  enqueue_2(&todo2,adj2[i]);
	  visited[adj2[i]] = offset_2; // mark as visited
	  cs2++; // increase cluster size
	  }
      }
      dequeue_2(&todo1,&activeP1);
      neighbours(activeP1,1);
      for(i=0;i<4;i++) {
	if(adj1[i] != -1) {
	  if(visited[adj1[i]] == offset_1) continue;
	  if((__u32)adj1[i] == goal || visited[adj1[i]] == offset_2) {
            while(!queue_2_empty_p(&todo1)){ dequeue_2(&todo1,&activeP1);}
            while(!queue_2_empty_p(&todo2)){ dequeue_2(&todo2,&activeP2);}
			return 1;
	  }
	  enqueue_2(&todo1,adj1[i]);
	  visited[adj1[i]] = offset_1;
	  cs1++;
	}
      }
  }
  while(!queue_2_empty_p(&todo1)){dequeue_2(&todo1,&activeP1);}
  while(!queue_2_empty_p(&todo2)){dequeue_2(&todo2,&activeP2);}
  return 0;
}
/******************************************************************************
 *****************************************************************************/
static inline __u8 connected(__u32 start, __u32 goal) {
    offset_1+=2;
    offset_2+=2;
    return breadthFirstSearch(start,goal);
}

/******************************************************************************
 *****************************************************************************/
__u8 static inline connected_s(__u32 start, __u32 goal) {
    ++offset_1;
    return breadthFirstSearch_s(start,goal);
}
/******************************************************************************
 *****************************************************************************/
static inline void   sweep_ibfs(void)
{
    static __u32 bidx;
    static double rnd_num;
    __u32 i=0;

    for(;i<2*N ;++i) {
        bidx = gsl_rng_uniform_int(r,2*N);
        rnd_num  = gsl_rng_uniform(r);
        Adjacent(bidx);
        if(bonds[bidx] == 1) { //edge is active hence delete it
            dB = -1;
            if(rnd_num < p_min_del) { // accept bond deletion
                bonds[bidx]*=-1;
                activeEdges--;
            }
            else {
                if(rnd_num < p_max_del) {
                    bonds[bidx]*=-1;
                    dN = (connected(edge[0],edge[1])? 0 : 1);
                    if(rnd_num >= rcWeights[dB == -1? dN : 2-dN]) bonds[bidx]*=-1; // undo bond deletion
                    else activeEdges--; // accept bond deletion
                }

            }

        }
        else {
            // Insert cedge. If adjacent vertices are already connected, the new
            // edge will be a non-tree edge, i.e. dN = 0
            dB = 1;
            if(rnd_num < p_min_ins) { // accept bond insertion
                bonds[bidx] *=-1;
                activeEdges++;
            }
            else {
                if(rnd_num < p_max_ins) {
                    dN = ( connected(edge[0],edge[1]) ?0 : -1);
                    if(rnd_num < rcWeights[dB == -1? dN : 2-dN]) {
                        bonds[bidx]*=-1;
                        activeEdges++;
                    }
                }

            }
        }
    }
}
/******************************************************************************
 *****************************************************************************/
static inline void   sweep_sbfs(void)
{
    static __u32 bidx;
    static double rnd_num;
    __u32 i=0;

    for(;i<2*N ;++i) {
        bidx = gsl_rng_uniform_int(r,2*N);
        rnd_num  = gsl_rng_uniform(r);
        Adjacent(bidx);
        if(bonds[bidx] == 1) { //edge is active hence delete it
            dB = -1;
            if(rnd_num < p_min_del) { // accept bond deletion
                bonds[bidx]*=-1;
                activeEdges--;
            }
            else {
                if(rnd_num < p_max_del) {
                    bonds[bidx]*=-1;
                    dN = (connected_s(edge[0],edge[1])? 0 : 1);
                    if(rnd_num >= rcWeights[dB == -1? dN : 2-dN]) bonds[bidx]*=-1; // undo bond deletion
                    else activeEdges--; // accept bond deletion
                }

            }

        }
        else {
            // Insert cedge. If adjacent vertices are already connected, the new
            // edge will be a non-tree edge, i.e. dN = 0
            dB = 1;
            if(rnd_num < p_min_ins) { // accept bond insertion
                bonds[bidx] *=-1;
                activeEdges++;
            }
            else {
                if(rnd_num < p_max_ins) {
                    dN = ( connected_s(edge[0],edge[1]) ?0 : -1);
                    if(rnd_num < rcWeights[dB == -1? dN : 2-dN]) {
                        bonds[bidx]*=-1;
                        activeEdges++;
                    }
                }

            }
        }
    }
}
/******************************************************************************
 *****************************************************************************/
static void extract_observables(__u32 i) {
    static __u32 j=0;
    static __u32 clust_cnt=0;
    static __u64 sum=0,sum_2;
    static __u32 maxc=0,chksum=0;
    num_bonds[i] = activeEdges;
    offset_1=0;
    for(j=0;j<N;++j) visited[j] = offset_1;
    offset_1=1;
    for(j=0;j<N&&chksum<=N;++j) {
        if(!visited[j]){
            breadthFirstSearch_s(j,N);
            clust_cnt++;
            if(cs1 > maxc) maxc = cs1;
            sum+=cs1*cs1;
            sum_2+=pow(cs1,4);
            chksum+=cs1;
        }

    }
    num_cluster[i] = clust_cnt;
    size_giant[i] = maxc;
    sec_cs_moment[i] = sum;
    four_cs_moment[i] = sum_2;
    maxc=clust_cnt=sum=sum_2=chksum=0;
    offset_1=3;
    offset_2=2;
}
/******************************************************************************
 *****************************************************************************/
static void generateTimeSeries(void)
{

  __u32 i;
  if(setup==1) {
    for(i=0;i<cutoff;++i)sweep_ibfs();
    equilibration=0;
    if(verbose)
        printf("Equilibration done!\n");
    for(i=0;i<steps;++i) {
        sweep_ibfs();
        extract_observables(i);
    }
  }
  if(setup==2) {
    for(i=0;i<cutoff;++i)sweep_sbfs();
    equilibration=0;
    if(verbose)
        printf("Equilibration done!\n");
    for(i=0;i<steps;++i) {
        sweep_sbfs();
        extract_observables(i);
    }

  }

}
/******************************************************************************
 *****************************************************************************/
char init_sweeny_ibfs(double _q,unsigned int _l,double _beta,double _coupl,
        unsigned int _cutoff,unsigned _tslength,unsigned int rng_seed,
        void *ts_0,void *ts_1,void * ts_2,void *ts_3,void *ts_4) {
    setup=1;
    offset_1 =1;
    offset_2 = 2;
    activeEdges =0;
    q = _q;
    DX = _l;
    beta = _beta;
    coupling = _coupl;
    cutoff = _cutoff;
    steps = _tslength;
    seed = rng_seed;
    num_bonds = (__u32 *)ts_0;
    num_cluster = (__u32 *)ts_1;
    size_giant = (__u32 *)ts_2;
    sec_cs_moment = (__u64 *)ts_3;
    four_cs_moment = (__u64 *)ts_4;
    if(!init())
        setup=0;
    return setup;

}
/******************************************************************************
 *****************************************************************************/
char init_sweeny_sbfs(double _q,unsigned int _l,double _beta,double _coupl,
        unsigned int _cutoff,unsigned _tslength,unsigned int rng_seed,
        void *ts_0,void *ts_1,void * ts_2,void *ts_3, void *ts_4) {
    setup=init_sweeny_ibfs(_q,_l,_beta,_coupl,_cutoff,_tslength,rng_seed,
            ts_0,ts_1,ts_2,ts_3,ts_4);
    if(setup)
        return setup=2;
    else
        return setup;
}
/******************************************************************************
 *****************************************************************************/
void destroy_sweeny_ibfs(void) {
    if(setup) {
        destroy();
    }
    setup=0;
}
/******************************************************************************
 *****************************************************************************/
void destroy_sweeny_sbfs(void) {
    destroy_sweeny_ibfs();
}
/******************************************************************************
 *****************************************************************************/
char simulate_sweeny_ibfs(void) {
    if(setup!=1)
        return 0;
    else
        generateTimeSeries();
    return 1;
}
/******************************************************************************
 *****************************************************************************/
char simulate_sweeny_sbfs(void) {
    if(setup!=2)
        return 0;
    else
        generateTimeSeries();
    return 1;
}

ext: c

max_stars_repo_name: ernmeel/sweeny

max_issues_repo_name: ernmeel/sweeny


text:
#include <math.h>
#include <string.h>

#include <gsl/gsl_matrix.h>

#include <gpc/gpc.h>
#include <egsl/egsl_macros.h>

#include "../csm_all.h"

#include "icp.h"

void sm_journal_open(const char* file) {
	file = 0;
/*	journal_open(file);*/
}

void ld_invalid_if_outside(LDP ld, double min_reading, double max_reading) {
	int i;
	for(i=0;i<ld->nrays;i++) {
		if(!ld_valid_ray(ld, i)) continue;
		double r = ld->readings[i];
		if( r <= min_reading || r > max_reading)
			ld->valid[i] = 0;
	}
}

void sm_icp(struct sm_params*params, struct sm_result*res) {
	res->valid = 0;

	LDP laser_ref  = params->laser_ref;
	LDP laser_sens = params->laser_sens;

	if(!ld_valid_fields(laser_ref) ||
	   !ld_valid_fields(laser_sens)) {
		return;
	}

	sm_debug("sm_icp: laser_sens has %d/%d; laser_ref has %d/%d rays valid\n",
		count_equal(laser_sens->valid, laser_sens->nrays, 1), laser_sens->nrays,
		count_equal(laser_ref->valid, laser_ref->nrays, 1), laser_ref->nrays);

	/** Mark as invalid the rays outside of (min_reading, max_reading] */
	ld_invalid_if_outside(laser_ref, params->min_reading, params->max_reading);
	ld_invalid_if_outside(laser_sens, params->min_reading, params->max_reading);

	sm_debug("sm_icp:  laser_sens has %d/%d; laser_ref has %d/%d rays valid (after removing outside interval [%f, %f])\n",
		count_equal(laser_sens->valid, laser_sens->nrays, 1), laser_sens->nrays,
		count_equal(laser_ref->valid, laser_ref->nrays, 1), laser_ref->nrays,
   	   params->min_reading, params->max_reading);

	if(JJ) jj_context_enter("sm_icp");

	egsl_push_named("sm_icp");

	if(params->use_corr_tricks || params->debug_verify_tricks)
		ld_create_jump_tables(laser_ref);

	ld_compute_cartesian(laser_ref);
	ld_compute_cartesian(laser_sens);

	if(params->do_alpha_test) {
		ld_simple_clustering(laser_ref, params->clustering_threshold);
		ld_compute_orientation(laser_ref, params->orientation_neighbourhood, params->sigma);
		ld_simple_clustering(laser_sens, params->clustering_threshold);
		ld_compute_orientation(laser_sens, params->orientation_neighbourhood, params->sigma);
	}

	if(JJ) jj_add("laser_ref",  ld_to_json(laser_ref));
	if(JJ) jj_add("laser_sens", ld_to_json(laser_sens));

	gsl_vector * x_new = gsl_vector_alloc(3);
	gsl_vector * x_old = vector_from_array(3, params->first_guess);

	if(params->do_visibility_test) {
		sm_debug("laser_ref:\n");
		visibilityTest(laser_ref, x_old);

		sm_debug("laser_sens:\n");
		gsl_vector * minus_x_old = gsl_vector_alloc(3);
		ominus(x_old,minus_x_old);
		visibilityTest(laser_sens, minus_x_old);
		gsl_vector_free(minus_x_old);
	}

	double error;
	int iterations;
	int nvalid;
	if(!icp_loop(params, x_old->data, x_new->data, &error, &nvalid, &iterations)) {
		sm_error("icp: ICP failed for some reason. \n");
		res->valid = 0;
		res->iterations = iterations;
		res->nvalid = 0;

	} else {
		/* It was succesfull */

		int restarted = 0;
		double best_error = error;
		gsl_vector * best_x = gsl_vector_alloc(3);
		gsl_vector_memcpy(best_x, x_new);

		if(params->restart &&
			(error/nvalid)>(params->restart_threshold_mean_error) ) {
			sm_debug("Restarting: %f > %f \n",(error/nvalid),(params->restart_threshold_mean_error));
			restarted = 1;
			double dt  = params->restart_dt;
			double dth = params->restart_dtheta;
			sm_debug("icp_loop: dt = %f dtheta= %f deg\n",dt,rad2deg(dth));

			double perturb[6][3] = {
				{dt,0,0}, {-dt,0,0},
				{0,dt,0}, {0,-dt,0},
				{0,0,dth}, {0,0,-dth}
			};

			int a; for(a=0;a<6;a++){
				sm_debug("-- Restarting with perturbation #%d\n", a);
				struct sm_params my_params = *params;
				gsl_vector * start = gsl_vector_alloc(3);
					gvs(start, 0, gvg(x_new,0)+perturb[a][0]);
					gvs(start, 1, gvg(x_new,1)+perturb[a][1]);
					gvs(start, 2, gvg(x_new,2)+perturb[a][2]);
				gsl_vector * x_a = gsl_vector_alloc(3);
				double my_error; int my_valid; int my_iterations;
				if(!icp_loop(&my_params, start->data, x_a->data, &my_error, &my_valid, &my_iterations)){
					sm_error("Error during restart #%d/%d. \n", a, 6);
					break;
				}
				iterations+=my_iterations;

				if(my_error < best_error) {
					sm_debug("--Perturbation #%d resulted in error %f < %f\n", a,my_error,best_error);
					gsl_vector_memcpy(best_x, x_a);
					best_error = my_error;
				}
				gsl_vector_free(x_a); gsl_vector_free(start);
			}
		}

		/* At last, we did it. */
		res->valid = 1;
		vector_to_array(best_x, res->x);
		sm_debug("icp: final x =  %s  \n", gsl_friendly_pose(best_x));

		if (restarted) { // recompute correspondences in case of restarts
			ld_compute_world_coords(laser_sens, res->x);
			if(params->use_corr_tricks)
				find_correspondences_tricks(params);
			else
				find_correspondences(params);
		}

		if(params->do_compute_covariance)  {

			val cov0_x, dx_dy1, dx_dy2;
			compute_covariance_exact(
				laser_ref, laser_sens, best_x,
				&cov0_x, &dx_dy1, &dx_dy2);

			val cov_x = sc(square(params->sigma), cov0_x);
/*			egsl_v2da(cov_x, res->cov_x); */

			res->cov_x_m = egsl_v2gslm(cov_x);
			res->dx_dy1_m = egsl_v2gslm(dx_dy1);
			res->dx_dy2_m = egsl_v2gslm(dx_dy2);

			if(0) {
				egsl_print("cov0_x", cov0_x);
				egsl_print_spectrum("cov0_x", cov0_x);

				val fim = ld_fisher0(laser_ref);
				val ifim = inv(fim);
				egsl_print("fim", fim);
				egsl_print_spectrum("ifim", ifim);
			}
		}

		res->error = best_error;
		res->iterations = iterations;
		res->nvalid = nvalid;

		gsl_vector_free(best_x);
	}
	gsl_vector_free(x_new);
	gsl_vector_free(x_old);

	egsl_pop_named("sm_icp");

	if(JJ) jj_context_exit();
}

void sm_icp_xy(struct sm_params*params, struct sm_result*res)
{
	res->valid = 0;

	LDP laser_ref  = params->laser_ref;
	LDP laser_sens = params->laser_sens;

	if(!ld_valid_fields(laser_ref) ||
	   !ld_valid_fields(laser_sens)) {
		return;
	}

/*
	sm_debug("sm_icp: laser_sens has %d/%d; laser_ref has %d/%d rays valid\n",
		count_equal(laser_sens->valid, laser_sens->nrays, 1), laser_sens->nrays,
		count_equal(laser_ref->valid, laser_ref->nrays, 1), laser_ref->nrays);
	*/

	/** Mark as invalid the rays outside of (min_reading, max_reading] */
	ld_invalid_if_outside(laser_ref, params->min_reading, params->max_reading);
	ld_invalid_if_outside(laser_sens, params->min_reading, params->max_reading);

/*
	sm_debug("sm_icp:  laser_sens has %d/%d; laser_ref has %d/%d rays valid (after removing outside interval [%f, %f])\n",
		count_equal(laser_sens->valid, laser_sens->nrays, 1), laser_sens->nrays,
		count_equal(laser_ref->valid, laser_ref->nrays, 1), laser_ref->nrays,
   	   params->min_reading, params->max_reading);

	if(JJ) jj_context_enter("sm_icp");

	egsl_push_named("sm_icp");
	*/

	if(params->use_corr_tricks || params->debug_verify_tricks)
		ld_create_jump_tables(laser_ref);
/*
	ld_compute_cartesian(laser_ref);
	ld_compute_cartesian(laser_sens);
*/

	if(params->do_alpha_test) {
		ld_simple_clustering(laser_ref, params->clustering_threshold);
		ld_compute_orientation(laser_ref, params->orientation_neighbourhood, params->sigma);
		ld_simple_clustering(laser_sens, params->clustering_threshold);
		ld_compute_orientation(laser_sens, params->orientation_neighbourhood, params->sigma);
	}

	if(JJ) jj_add("laser_ref",  ld_to_json(laser_ref));
	if(JJ) jj_add("laser_sens", ld_to_json(laser_sens));

	gsl_vector * x_new = gsl_vector_alloc(3);
	gsl_vector * x_old = vector_from_array(3, params->first_guess);

	if(params->do_visibility_test) {
		sm_debug("laser_ref:\n");
		visibilityTest(laser_ref, x_old);

		sm_debug("laser_sens:\n");
		gsl_vector * minus_x_old = gsl_vector_alloc(3);
		ominus(x_old,minus_x_old);
		visibilityTest(laser_sens, minus_x_old);
		gsl_vector_free(minus_x_old);
	}

	double error;
	int iterations;
	int nvalid;
	if(!icp_loop(params, x_old->data, x_new->data, &error, &nvalid, &iterations)) {
		sm_error("icp: ICP failed for some reason. \n");
		res->valid = 0;
		res->iterations = iterations;
		res->nvalid = 0;

	} else {
		/* It was succesfull */

		double best_error = error;
		gsl_vector * best_x = gsl_vector_alloc(3);
		gsl_vector_memcpy(best_x, x_new);

		if(params->restart &&
			(error/nvalid)>(params->restart_threshold_mean_error) ) {
			sm_debug("Restarting: %f > %f \n",(error/nvalid),(params->restart_threshold_mean_error));
			double dt  = params->restart_dt;
			double dth = params->restart_dtheta;
			sm_debug("icp_loop: dt = %f dtheta= %f deg\n",dt,rad2deg(dth));

			double perturb[6][3] = {
				{dt,0,0}, {-dt,0,0},
				{0,dt,0}, {0,-dt,0},
				{0,0,dth}, {0,0,-dth}
			};

			int a; for(a=0;a<6;a++){
				sm_debug("-- Restarting with perturbation #%d\n", a);
				struct sm_params my_params = *params;
				gsl_vector * start = gsl_vector_alloc(3);
					gvs(start, 0, gvg(x_new,0)+perturb[a][0]);
					gvs(start, 1, gvg(x_new,1)+perturb[a][1]);
					gvs(start, 2, gvg(x_new,2)+perturb[a][2]);
				gsl_vector * x_a = gsl_vector_alloc(3);
				double my_error; int my_valid; int my_iterations;
				if(!icp_loop(&my_params, start->data, x_a->data, &my_error, &my_valid, &my_iterations)){
					sm_error("Error during restart #%d/%d. \n", a, 6);
					break;
				}
				iterations+=my_iterations;

				if(my_error < best_error) {
					sm_debug("--Perturbation #%d resulted in error %f < %f\n", a,my_error,best_error);
					gsl_vector_memcpy(best_x, x_a);
					best_error = my_error;
				}
				gsl_vector_free(x_a); gsl_vector_free(start);
			}
		}

		/* At last, we did it. */
		res->valid = 1;
		vector_to_array(best_x, res->x);
		sm_debug("icp: final x =  %s  \n", gsl_friendly_pose(best_x));

		if(params->do_compute_covariance)  {

			val cov0_x, dx_dy1, dx_dy2;
			compute_covariance_exact(
				laser_ref, laser_sens, best_x,
				&cov0_x, &dx_dy1, &dx_dy2);

			val cov_x = sc(square(params->sigma), cov0_x);
//			egsl_v2da(cov_x, res->cov_x);

			res->cov_x_m = egsl_v2gslm(cov_x);
			res->dx_dy1_m = egsl_v2gslm(dx_dy1);
			res->dx_dy2_m = egsl_v2gslm(dx_dy2);

			if(0) {
				egsl_print("cov0_x", cov0_x);
				egsl_print_spectrum("cov0_x", cov0_x);

				val fim = ld_fisher0(laser_ref);
				val ifim = inv(fim);
				egsl_print("fim", fim);
				egsl_print_spectrum("ifim", ifim);
			}
		}

		res->error = best_error;
		res->iterations = iterations;
		res->nvalid = nvalid;

		gsl_vector_free(x_new);
		gsl_vector_free(x_old);
		gsl_vector_free(best_x);
	}
/*
	egsl_pop_named("sm_icp");

	if(JJ) jj_context_exit();
*/
}


text:
#include <gsl/gsl_matrix.h>
#include "BandEnergy.h"

int main(int argc, char *argv[]) {
    HTightBinding *Hrs = ExtractHTightBinding("test_data/Fe_soc/Fe_soc_hr.dat");

    gsl_matrix *R = gsl_matrix_calloc(3, 3); // R -> all zeros
    // Overall scale for R doesn't matter.
    gsl_matrix_set(R, 0, 1, 1.0);
    gsl_matrix_set(R, 0, 2, 1.0);
    gsl_matrix_set(R, 1, 0, 1.0);
    gsl_matrix_set(R, 1, 2, 1.0);
    gsl_matrix_set(R, 2, 0, 1.0);
    gsl_matrix_set(R, 2, 1, 1.0);
    double num_electrons = 8.0;
    int na = 8;
    int nb = 8;
    int nc = 8;
    bool use_cache = true;
    double E_Fermi = 0.0;

    double energy = BandEnergy(&E_Fermi, Hrs, R, num_electrons, na, nb, nc, use_cache);

    printf("energy = %f\n", energy);
    printf("E_Fermi = %f\n", E_Fermi);

    return 0;
}

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
#include <math.h>
#include <stdlib.h>
#include <malloc.h>
#include <stdio.h>
#include <assert.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include "../../theory/basics.c"
#include "../../emu13/emu.c"
#include "../../theory/parameters.c"
#include "../../theory/cosmo3D.c"
#include "../../theory/redshift.c"

typedef struct {
  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
  gsl_spline2d *spline;
  gsl_interp_accel *xacc;
  gsl_interp_accel *yacc;
  double *grid_values;
} spline_interpolator;
spline_interpolator global_baryon_spline;

void init_baryon_spline_interpolation_to_TNG100(){

  const size_t N = 100;             /* number of points to interpolate */
  const double lnk_values[] = { 0.0, 1.0 };
  const double z_values[] = { 3.71, 3.49, 3.28, 2.90, 2.44, 2.10, 1.74, 1.41, 1.04, 0.7, 0.35, 0.18, 0.0 }; /* define unit square */
  const size_t nk = sizeof(lnk_values) / sizeof(double); /* y grid points */
  const size_t nz = sizeof(z_values) / sizeof(double); /* x grid points */
  global_baryon_spline.grid_values = malloc(nk * nz * sizeof(double));
  global_baryon_spline.spline = gsl_spline2d_alloc(global_baryon_spline.T, nk, nz);
  global_baryon_spline.xacc = gsl_interp_accel_alloc();
  global_baryon_spline.yacc = gsl_interp_accel_alloc();
  size_t i, j;

  for(int i = 0; i < nk; i++){
    for(int j = 0; j < nz; j++){
      gsl_spline2d_set(spline, global_baryon_spline.grid_values, i, j, TNG100[i][j]);
    }
  }

  /* initialize interpolation */
  gsl_spline2d_init(global_baryon_spline.spline, lnk_values, z_values, Pk_ratio, nk, nz);

}

void free_baryon_spline_interpolation_to_TNG100(){
  gsl_spline2d_free(global_baryon_spline.spline);
  gsl_interp_accel_free(global_baryon_spline.xacc);
  gsl_interp_accel_free(global_baryon_spline.yacc);
  free(global_baryon_spline.grid_values);
}

double fraction_Pdelta_baryon_from_sims_DES(double k,double z)
{
  return gsl_spline2d_eval(global_baryon_spline.spline, k, z, global_baryon_spline.xacc, global_baryon_spline.yacc);
}

double fraction_Pdelta_baryon_from_spline_DES(double k,double z)
{
  int i;
  double val;
  static int READ_TABLE=0;
  static double **table;
  FILE *F;
  int baryon_zbin=13;
  int baryon_kbin=325;
  double logk_min = -3.3010299956639813;
  double logk_max = 3.1760912590556813;
  double dz=(redshift.shear_zdistrpar_zmax)/(10.0);
  static double dlogk=(logk_max-logk_min)/(double(baryon_kbin-1));
  val = interpol2d(AGN_DESdepth, baryon_kbin, logk_min, logk_max, dk, log(k), baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0);
  return val;
}

//   if (READ_TABLE==0){
//     if (strcmp(pdeltaparams.baryons,"AGN_DESdepth")==0) P_type = 0;
//     if (strcmp(pdeltaparams.baryons,"NOSN_DESdepth")==0) P_type = 1;
//     if (strcmp(pdeltaparams.baryons,"NOSN_NOZCOOL_DESdepth")==0) P_type = 2;
//     if (strcmp(pdeltaparams.baryons,"NOZCOOL_DESdepth")==0) P_type = 3;
//     if (strcmp(pdeltaparams.baryons,"REF_DESdepth")==0) P_type = 4;
//     if (strcmp(pdeltaparams.baryons,"WDENS_DESdepth") ==0) P_type = 5;
//     if (strcmp(pdeltaparams.baryons,"DBLIMFV1618_DESdepth")  ==0) P_type = 6;
//     if (strcmp(pdeltaparams.baryons,"WML4_DESdepth")==0) P_type = 7;
//     if (strcmp(pdeltaparams.baryons,"WML1V848_DESdepth")==0) P_type = 8;
//     if (strcmp(pdeltaparams.baryons,"AD_DESdepth")==0) P_type = 9;
//     if (strcmp(pdeltaparams.baryons,"CX_DESdepth")==0) P_type = 10;
//     if (strcmp(pdeltaparams.baryons,"CW_DESdepth")==0) P_type = 11;
//     if (strcmp(pdeltaparams.baryons,"A_DESdepth")==0) P_type = 12;
//     if (strcmp(pdeltaparams.baryons,"CSF_DESdepth")==0) P_type = 13;
//   }
//   switch (P_type){
//     case 0:  val = interpol2d(AGN_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 1:  val = interpol2d(NOSN_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 2:  val = interpol2d(NOSN_NOZCOOL_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 3:  val = interpol2d(NOZCOOL_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 4:  val = interpol2d(REF_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 5:  val = interpol2d(WDENS_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 6:  val = interpol2d(DBLIMFV1618_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 7:  val = interpol2d(WML4_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 8:  val = interpol2d(WML1V848_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 9:  val = interpol2d(AD_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 10:  val = interpol2d(CX_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 11:  val = interpol2d(CW_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 12:  val = interpol2d(A_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;
//     case 13:  val = interpol2d(CSF_DESdepth, baryon_kbin, 0.3, 10.0, dk, kintern, baryon_zbin, 0.0, 2.0, dz, z, 0.0, 0.0); break;

//     default:
//             printf("baryons.c: %s baryonic scenario not defined\n",pdeltaparams.baryons);

//             break;
//     }
//     kintern=k/cosmology.coverH0;

//  return val;
// }

//    if (strcmp(pdeltaparams.baryons,"AGN_LSSTdepth")==0) P_type = 14;
//     if (strcmp(pdeltaparams.baryons,"NOSN_LSSTdepth")==0) P_type = 15;
//     if (strcmp(pdeltaparams.baryons,"NOSN_NOZCOOL_LSSTdepth")==0) P_type = 16;
//     if (strcmp(pdeltaparams.baryons,"NOZCOOL_LSSTdepth")==0) P_type = 17;
//     if (strcmp(pdeltaparams.baryons,"REF_LSSTdepth")==0) P_type = 18;
//     if (strcmp(pdeltaparams.baryons,"WDENS_LSSTdepth") ==0) P_type = 19;
//     if (strcmp(pdeltaparams.baryons,"DBLIMFV1618_LSSTdepth")  ==0) P_type = 20;
//     if (strcmp(pdeltaparams.baryons,"WML4_LSSTdepth")==0) P_type = 21;
//     if (strcmp(pdeltaparams.baryons,"WML1V848_LSSTdepth")==0) P_type = 22;
//     if (strcmp(pdeltaparams.baryons,"AD_LSSTdepth")==0) P_type = 23;
//     if (strcmp(pdeltaparams.baryons,"CX_LSSTdepth")==0) P_type = 24;
//     if (strcmp(pdeltaparams.baryons,"CW_LSSTdepth")==0) P_type = 25;
//   if (table!=0) free_double_matrix(table, 0, OWLS_kbin-1, 0, OWLS_zbin-1);
//     table   = create_double_matrix(0, OWLS_kbin-1, 0, OWLS_zbin-1);
//     printf("%s\n",file.DATA_FILE);
//     F=fopen(file.DATA_FILE,"r");
//       for(i=0;i<OWLS_kbin;i++){
// 	fscanf(F,"%le %le %le %le %le %le %le %le %le %le %le\n",&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11);
// //	printf("%le %le %le %le %le %le %le %le %le %le %le\n",a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11);
// 	table[i][0]=a1;
// 	table[i][1]=a2;
// 	table[i][2]=a3;
// 	table[i][3]=a4;
// 	table[i][4]=a5;
// 	table[i][5]=a6;
// 	table[i][6]=a7;
// 	table[i][7]=a8;
// 	table[i][8]=a9;
// 	table[i][9]=a10;
// 	table[i][10]=a11;
//       }
//       READ_TABLE=1;
//   }
//  kintern=k/cosmology.coverH0;
//  val = interpol2d(table, OWLS_kbin, 0.3, 10.0, dk, kintern,OWLS_zbin, 0.0, 2.0, dz, z, 0.0, 0.0);
//  return val;
// }


text:
//IN method.
//Affine transformation (weights and biases) of Ni inputs to No outputs.
//This version uses CBLAS.

//Input X has Ni neurons and output Y has No neurons.
//Each output neuron has a bias term, so B is a vector of length No.

//The vecs of length Ni are always contiguous in memory, such that:

//If col-major: Y[:,l] = W' * X[:,l] + B
//where:
//X has size Ni x L
//Y has size No x L
//W has size Ni x No
//B has size No x 1

//If row-major: Y[l,:] = X[l,:] * W' + B
//X has size L x Ni
//Y has size L x No
//W has size No x Ni
//B has size 1 x No

//For a different set-up that allows affine transformation of vecs in
//any orientation, use the affine function from math.

#include <cblas.h>

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

int affine_cblas_s (float *Y, const float *X, const float *W, const float *B, const size_t Ni, const size_t No, const size_t L);
int affine_cblas_d (double *Y, const double *X, const double *W, const double *B, const size_t Ni, const size_t No, const size_t L);
int affine_cblas_c (float *Y, const float *X, const float *W, const float *B, const size_t Ni, const size_t No, const size_t L);
int affine_cblas_z (double *Y, const double *X, const double *W, const double *B, const size_t Ni, const size_t No, const size_t L);

int affine_cblas_s (float *Y, const float *X, const float *W, const float *B, const size_t Ni, const size_t No, const size_t L)
{
    for (size_t l=L; l>0u; --l, X+=Ni, Y+=No)
    {
        cblas_scopy((int)No,B,1,Y,1);
        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,1.0f,W,(int)Ni,X,1,1.0f,Y,1);
    }

    return 0;
}

int affine_cblas_d (double *Y, const double *X, const double *W, const double *B, const size_t Ni, const size_t No, const size_t L)
{
    for (size_t l=L; l>0u; --l, X+=Ni, Y+=No)
    {
        cblas_dcopy((int)No,B,1,Y,1);
        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,1.0,W,(int)Ni,X,1,1.0,Y,1);
    }

    return 0;
}

int affine_cblas_c (float *Y, const float *X, const float *W, const float *B, const size_t Ni, const size_t No, const size_t L)
{
    const float o[2] = {1.0f,0.0f};

    for (size_t l=L; l>0u; --l, X+=2u*Ni, Y+=2u*No)
    {
        cblas_ccopy((int)No,B,1,Y,1);
        cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,o,W,(int)Ni,X,1,o,Y,1);
    }

    return 0;
}

int affine_cblas_z (double *Y, const double *X, const double *W, const double *B, const size_t Ni, const size_t No, const size_t L)
{
    const double o[2] = {1.0,0.0};

    for (size_t l=L; l>0u; --l, X+=2u*Ni, Y+=2u*No)
    {
        cblas_zcopy((int)No,B,1,Y,1);
        cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,o,W,(int)Ni,X,1,o,Y,1);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#ifndef __MATH_GSL_VECTOR__
#define __MATH_GSL_VECTOR__

#include <assert.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_cblas.h>

#include <gsl/gsl_matrix.h>
#include "math/vectorops.h"

class GslVectorItem {
 public:
  GslVectorItem(gsl_vector* ptr, size_t index) :
    ptr_(ptr),
    index_(index) { }

  operator const double() {
    return gsl_vector_get(ptr_, index_);
  }

  double operator =(const double v) {
    gsl_vector_set(ptr_, index_, v);
    return v;
  }

  double operator +=(const double v) {
    double oldv = gsl_vector_get(ptr_, index_);
    gsl_vector_set(ptr_, index_, oldv + v);
    return oldv + v;
  }

  double operator -=(const double v) {
    double oldv = gsl_vector_get(ptr_, index_);
    gsl_vector_set(ptr_, index_, oldv - v);
    return oldv - v;
  }
 private:
  gsl_vector* ptr_;
  size_t index_;
};

class GslVectorBase {
 public:
  /*
  double operator[](const size_t index) const {
    assert(ptr_ != NULL);
    return gsl_vector_get(ptr_, index);
  }
  */

  GslVectorItem operator[](const size_t index) const {
    assert(ptr_ != NULL);
    return GslVectorItem(ptr_, index);
  }

  GslVectorBase& operator+=(const gsl_vector* x) {
    assert(ptr_ != NULL);
    gsl_vector_add(ptr_, x);
    return *this;
  }

  GslVectorBase& operator-=(const gsl_vector* x) {
    assert(ptr_ != NULL);
    gsl_vector_sub(ptr_, x);
    return *this;
  }

  GslVectorBase& operator+=(const GslVectorBase& x) {
    assert(ptr_ != NULL);
    gsl_vector_add(ptr_, x.ptr());
    return *this;
  }

  GslVectorBase& operator-=(const GslVectorBase& x) {
    assert(ptr_ != NULL);
    gsl_vector_sub(ptr_, x.ptr());
    return *this;
  }

  GslVectorBase& operator*=(const gsl_vector* x) {
    assert(ptr_ != NULL);
    gsl_vector_mul(ptr_, x);
    return *this;
  }

  GslVectorBase& operator*=(const GslVectorBase& x) {
    assert(ptr_ != NULL);
    gsl_vector_mul(ptr_, x.ptr());
    return *this;
  }

  GslVectorBase& operator/=(const GslVectorBase& x) {
    assert(ptr_ != NULL);
    gsl_vector_div(ptr_, x.ptr());
    return *this;
  }

  double Sum() const {
    return gsl_blas_dsum(ptr_);
  }

  double L2Norm() const {
    return gsl_blas_dnrm2(ptr_);
  }

  void Normalize() const {
    assert(ptr_ != NULL);
    double s = Sum();
    gsl_vector_scale(ptr_, 1. / s);
  }

  size_t size() const {
    return ptr_->size;
  }

  GslVectorBase& operator*=(double v) {
    assert(ptr_ != NULL);
    gsl_vector_scale(ptr_, v);
    return *this;
  }

  GslVectorBase& operator/=(double v) {
    assert(ptr_ != NULL);
    gsl_vector_scale(ptr_, 1. / v);
    return *this;
  }

  // Note that the standalone product is a dot product!
  const double operator*(const gsl_vector* x) const {
    double result;
    assert(ptr_ != NULL);
    gsl_blas_ddot(ptr_, x, &result);
    return result;
  }

  const double operator*(const GslVectorBase& x) const {
    double result;
    assert(ptr_ != NULL);
    gsl_blas_ddot(ptr_, x.ptr(), &result);
    return result;
  }

  GslVectorBase& operator+=(const double x) {
    for (size_t ii = 0; ii < ptr_->size; ++ii) {
      gsl_vector_set(ptr_, ii, gsl_vector_get(ptr_, ii) + x);
    }
    return *this;
  }

  GslVectorBase& operator-=(const double x) {
    for (size_t ii = 0; ii < ptr_->size; ++ii) {
      gsl_vector_set(ptr_, ii, gsl_vector_get(ptr_, ii) - x);
    }
    return *this;
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    assert(ptr_ != NULL);
    gsl_vector_memcpy(ptr_, x);

    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    assert(ptr_ != NULL);
    return *this = x.ptr();
  }

  GslVectorBase& operator=(const double v) {
    if (v == 0.0) {
      SetZero();
    } else {
      SetAll(v);
    }
    return *this;
  }

  void SetZero() {
    assert(ptr_ != NULL);
    gsl_vector_set_zero(ptr_);
  }

  void SetAll(const double v) {
    assert(ptr_ != NULL);
    gsl_vector_set_all(ptr_, v);
  }

  int Fprintf(FILE* stream, const char* format) const {
    assert(ptr_ != NULL);
    return gsl_vector_fprintf(stream, ptr_, format);
  }

  int Fscanf(FILE* stream) {
    assert(ptr_ != NULL);
    return gsl_vector_fscanf(stream, ptr_);
  }

  const gsl_vector* ptr() const { return ptr_; }
  gsl_vector* mutable_ptr() { return ptr_; }

 protected:
  GslVectorBase() : ptr_(NULL) { }
  gsl_vector* ptr_;

 private:
  GslVectorBase(const GslVectorBase&) { }
};

class GslVector : public GslVectorBase {
 public:
  GslVector(const size_t size) : GslVectorBase() {
    Allocate(size);
  }

  GslVector() : GslVectorBase() {
  }

  ~GslVector() {
    if (ptr_ != NULL) {
      gsl_vector_free(ptr_);
    }
  }

  void Reset(gsl_vector* val) {
    if (ptr_ != NULL) {
      gsl_vector_free(ptr_);
    }
    ptr_ = val;
  }

  void Allocate(const size_t size) {
    assert(ptr_ == NULL);
    ptr_ = gsl_vector_alloc(size);
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const double v) {
    GslVectorBase::operator=(v);
    return *this;
  }
 private:
  GslVector(const GslVector&) { }
};

class GslMatrixRow : public GslVectorBase {
 public:
  GslMatrixRow(GslMatrix& matrix, const size_t row) :
   view_(gsl_matrix_row(matrix.mutable_ptr(), row)) {
     ptr_ = &view_.vector;
  }

  GslMatrixRow(gsl_matrix* matrix, const size_t row) :
   view_(gsl_matrix_row(matrix, row)) {
     ptr_ = &view_.vector;
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const double v) {
    GslVectorBase::operator=(v);
    return *this;
  }
 private:
  gsl_vector_view view_;
  GslMatrixRow(const GslMatrixRow&) { }
};

class GslMatrixColumn : public GslVectorBase {
 public:
  GslMatrixColumn(GslMatrix& matrix, const size_t col) :
   view_(gsl_matrix_column(matrix.mutable_ptr(), col)) {
     ptr_ = &view_.vector;
  }

  GslMatrixColumn(gsl_matrix* matrix, const size_t col) :
   view_(gsl_matrix_column(matrix, col)) {
     ptr_ = &view_.vector;
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const double v) {
    GslVectorBase::operator=(v);
    return *this;
  }
 private:
  gsl_vector_view view_;
  GslMatrixColumn(const GslMatrixColumn&) { }
};

class GslMatrixDiagonal : public GslVectorBase {
 public:
  GslMatrixDiagonal(GslMatrix& matrix) :
   view_(gsl_matrix_diagonal(matrix.mutable_ptr())) {
     ptr_ = &view_.vector;
  }

  GslMatrixDiagonal(gsl_matrix* matrix) :
   view_(gsl_matrix_diagonal(matrix)) {
     ptr_ = &view_.vector;
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const double v) {
    GslVectorBase::operator=(v);
    return *this;
  }
 private:
  gsl_vector_view view_;
  GslMatrixDiagonal(const GslMatrixDiagonal&) { }
};

class GslSubvector : public GslVectorBase {
 public:
 GslSubvector(GslVectorBase& vector, size_t i, size_t n) :
  view_(gsl_vector_subvector(vector.mutable_ptr(), i, n)) {
     ptr_ = &view_.vector;
  }

  GslVectorBase& operator=(const gsl_vector* x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const GslVectorBase& x) {
    GslVectorBase::operator=(x);
    return *this;
  }

  GslVectorBase& operator=(const double v) {
    GslVectorBase::operator=(v);
    return *this;
  }
 private:
  gsl_vector_view view_;
  GslSubvector(const GslSubvector&) { }
};

#endif  // __MATH_GSL_VECTOR__

ext: h

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri


text:
/* ******************************************************************************
 *
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://www.apache.org/licenses/LICENSE-2.0.
 *
 *  See the NOTICE file distributed with this work for additional
 *  information regarding copyright ownership.
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations
 * under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
// @author raver119@gmail.com
//

#ifndef LIBND4J_GEMM_H
#define LIBND4J_GEMM_H

#include <cblas.h>
#include <math/templatemath.h>
#include <system/op_boilerplate.h>

namespace sd {
     namespace blas {
         template <typename T>
         static void * transpose(int orderSource, int orderTarget, int rows, int cols, void *source);

         static inline int linearIndexC(int rows, int cols, int r, int c);
         static inline int linearIndexF(int rows, int cols, int r, int c);

         template <typename X, typename Y, typename Z>
         class GEMM {
         protected:
         public:
             static void op(int Order, int TransA, int TransB, int M, int N, int K, double alpha, void *A, int lda, void *B, int ldb, double beta, void *C, int ldc);
         };

         template <typename X, typename Y, typename Z>
         class GEMV : public sd::blas::GEMM<X, Y, Z>{
         public:
             static void op(int TRANS, int M, int N, double alpha, void* vA, int lda, void* vX, int incx, double beta, void* vY, int incy );
         };

         int FORCEINLINE linearIndexC(int rows, int cols, int r, int c) {
             return (r * cols + c);
         }

         int FORCEINLINE linearIndexF(int rows, int cols, int r, int c) {
             return (c * rows + r);
         }

    }
}

#endif //LIBND4J_GEMM_H


text:
/* blas/gsl_blas.h
 *
 * 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
 */
#ifndef __GSL_BLAS_H__
#define __GSL_BLAS_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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

/* ========================================================================
 * Level 1
 * ========================================================================
 */

int gsl_blas_sdsdot (float alpha,
                     const gsl_vector_float * X,
                     const gsl_vector_float * Y,
                     float * result
                     );

int gsl_blas_dsdot (const gsl_vector_float * X,
                    const gsl_vector_float * Y,
                    double * result
                    );

int gsl_blas_sdot (const gsl_vector_float * X,
                   const gsl_vector_float * Y,
                   float * result
                   );

int gsl_blas_ddot (const gsl_vector * X,
                   const gsl_vector * Y,
                   double * result
                   );

int  gsl_blas_cdotu (const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_complex_float * dotu);

int  gsl_blas_cdotc (const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_complex_float * dotc);

int  gsl_blas_zdotu (const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_complex * dotu);

int  gsl_blas_zdotc (const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_complex * dotc);

float  gsl_blas_snrm2  (const gsl_vector_float * X);
float  gsl_blas_sasum  (const gsl_vector_float * X);
double gsl_blas_dnrm2  (const gsl_vector * X);
double gsl_blas_dasum  (const gsl_vector * X);
float  gsl_blas_scnrm2 (const gsl_vector_complex_float * X);
float  gsl_blas_scasum (const gsl_vector_complex_float * X);
double gsl_blas_dznrm2 (const gsl_vector_complex * X);
double gsl_blas_dzasum (const gsl_vector_complex * X);

CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * X);
CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * X);
CBLAS_INDEX_t gsl_blas_icamax (const gsl_vector_complex_float * X);
CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * X);

int  gsl_blas_sswap (gsl_vector_float * X,
                     gsl_vector_float * Y);

int  gsl_blas_scopy (const gsl_vector_float * X,
                     gsl_vector_float * Y);

int  gsl_blas_saxpy (float alpha,
                     const gsl_vector_float * X,
                     gsl_vector_float * Y);

int  gsl_blas_dswap (gsl_vector * X,
                     gsl_vector * Y);

int  gsl_blas_dcopy (const gsl_vector * X,
                     gsl_vector * Y);

int  gsl_blas_daxpy (double alpha,
                     const gsl_vector * X,
                     gsl_vector * Y);

int  gsl_blas_cswap (gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

int  gsl_blas_ccopy (const gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

int  gsl_blas_caxpy (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

int  gsl_blas_zswap (gsl_vector_complex * X,
                     gsl_vector_complex * Y);

int  gsl_blas_zcopy (const gsl_vector_complex * X,
                     gsl_vector_complex * Y);

int  gsl_blas_zaxpy (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     gsl_vector_complex * Y);

int  gsl_blas_srotg (float a[], float b[], float c[], float s[]);

int  gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[]);

int  gsl_blas_srot (gsl_vector_float * X,
                    gsl_vector_float * Y,
                    float c, float s);

int  gsl_blas_srotm (gsl_vector_float * X,
                     gsl_vector_float * Y,
                     const float P[]);

int  gsl_blas_drotg (double a[], double b[], double c[], double s[]);

int  gsl_blas_drotmg (double d1[], double d2[], double b1[],
                      double b2, double P[]);

int  gsl_blas_drot (gsl_vector * X,
                    gsl_vector * Y,
                    const double c, const double s);

int  gsl_blas_drotm (gsl_vector * X,
                     gsl_vector * Y,
                     const double P[]);

void gsl_blas_sscal  (float  alpha, gsl_vector_float * X);
void gsl_blas_dscal  (double alpha, gsl_vector * X);
void gsl_blas_cscal  (const gsl_complex_float alpha, gsl_vector_complex_float * X);
void gsl_blas_zscal  (const gsl_complex alpha, gsl_vector_complex * X);
void gsl_blas_csscal (float  alpha, gsl_vector_complex_float * X);
void gsl_blas_zdscal (double alpha, gsl_vector_complex * X);

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
int  gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_vector_float * X,
                     float beta,
                     gsl_vector_float * Y);

int  gsl_blas_strmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_float * A,
                     gsl_vector_float * X);

int  gsl_blas_strsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_float * A,
                     gsl_vector_float * X);

int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_vector * X,
                     double beta,
                     gsl_vector * Y);

int  gsl_blas_dtrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix * A,
                     gsl_vector * X);

int  gsl_blas_dtrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix * A,
                     gsl_vector * X);

int  gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_vector_complex_float * X,
                     const gsl_complex_float beta,
                     gsl_vector_complex_float * Y);

int  gsl_blas_ctrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex_float * A,
                     gsl_vector_complex_float * X);

int  gsl_blas_ctrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex_float * A,
                     gsl_vector_complex_float * X);

int  gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_vector_complex * X,
                     const gsl_complex beta,
                     gsl_vector_complex * Y);

int  gsl_blas_ztrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex * A,
                     gsl_vector_complex * X);

int  gsl_blas_ztrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex * A,
                     gsl_vector_complex *X);

/*
 * Routines with S and D prefixes only
 */
int  gsl_blas_ssymv (CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_vector_float * X,
                     float beta,
                     gsl_vector_float * Y);

int  gsl_blas_sger (float alpha,
                    const gsl_vector_float * X,
                    const gsl_vector_float * Y,
                    gsl_matrix_float * A);

int  gsl_blas_ssyr (CBLAS_UPLO_t Uplo,
                    float alpha,
                    const gsl_vector_float * X,
                    gsl_matrix_float * A);

int  gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_vector_float * X,
                     const gsl_vector_float * Y,
                     gsl_matrix_float * A);

int  gsl_blas_dsymv (CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_vector * X,
                     double beta,
                     gsl_vector * Y);
int  gsl_blas_dger (double alpha,
                    const gsl_vector * X,
                    const gsl_vector * Y,
                    gsl_matrix * A);

int  gsl_blas_dsyr (CBLAS_UPLO_t Uplo,
                    double alpha,
                    const gsl_vector * X,
                    gsl_matrix * A);

int  gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_vector * X,
                     const gsl_vector * Y,
                     gsl_matrix * A);

/*
 * Routines with C and Z prefixes only
 */

int  gsl_blas_chemv (CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_vector_complex_float * X,
                     const gsl_complex_float beta,
                     gsl_vector_complex_float * Y);

int  gsl_blas_cgeru (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

int  gsl_blas_cgerc (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

int  gsl_blas_cher (CBLAS_UPLO_t Uplo,
                    float alpha,
                    const gsl_vector_complex_float * X,
                    gsl_matrix_complex_float * A);

int  gsl_blas_cher2 (CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

int  gsl_blas_zhemv (CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_vector_complex * X,
                     const gsl_complex beta,
                     gsl_vector_complex * Y);

int  gsl_blas_zgeru (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

int  gsl_blas_zgerc (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

int  gsl_blas_zher (CBLAS_UPLO_t Uplo,
                    double alpha,
                    const gsl_vector_complex * X,
                    gsl_matrix_complex * A);

int  gsl_blas_zher2 (CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
int  gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_matrix_float * B,
                     float beta,
                     gsl_matrix_float * C);

int  gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_matrix_float * B,
                     float beta,
                     gsl_matrix_float * C);

int  gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                     float alpha,
                     const gsl_matrix_float * A,
                     float beta,
                     gsl_matrix_float * C);

int  gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                      float alpha,
                      const gsl_matrix_float * A,
                      const gsl_matrix_float * B,
                      float beta,
                      gsl_matrix_float * C);

int  gsl_blas_strmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     float alpha,
                     const gsl_matrix_float * A,
                     gsl_matrix_float * B);

int  gsl_blas_strsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     float alpha,
                     const gsl_matrix_float * A,
                     gsl_matrix_float * B);

int  gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_matrix * B,
                     double beta,
                     gsl_matrix * C);

int  gsl_blas_dsymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_matrix * B,
                     double beta,
                     gsl_matrix * C);

int  gsl_blas_dsyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     double alpha,
                     const gsl_matrix * A,
                     double beta,
                     gsl_matrix * C);

int  gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      double alpha,
                      const  gsl_matrix * A,
                      const  gsl_matrix * B,
                      double beta,
                      gsl_matrix * C);

int  gsl_blas_dtrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     double alpha,
                     const gsl_matrix * A,
                     gsl_matrix * B);

int  gsl_blas_dtrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     double alpha,
                     const gsl_matrix * A,
                     gsl_matrix * B);

int  gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

int  gsl_blas_csymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

int  gsl_blas_csyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

int  gsl_blas_csyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex_float alpha,
                      const gsl_matrix_complex_float * A,
                      const gsl_matrix_complex_float * B,
                      const gsl_complex_float beta,
                      gsl_matrix_complex_float * C);

int  gsl_blas_ctrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     gsl_matrix_complex_float * B);

int  gsl_blas_ctrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     gsl_matrix_complex_float * B);

int  gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

int  gsl_blas_zsymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

int  gsl_blas_zsyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

int  gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex alpha,
                      const gsl_matrix_complex * A,
                      const gsl_matrix_complex * B,
                      const gsl_complex beta,
                      gsl_matrix_complex *C);

int  gsl_blas_ztrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     gsl_matrix_complex * B);

int  gsl_blas_ztrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     gsl_matrix_complex * B);

/*
 * Routines with prefixes C and Z only
 */
int  gsl_blas_chemm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

int  gsl_blas_cherk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     float alpha,
                     const gsl_matrix_complex_float * A,
                     float beta,
                     gsl_matrix_complex_float * C);

int  gsl_blas_cher2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex_float alpha,
                      const gsl_matrix_complex_float * A,
                      const gsl_matrix_complex_float * B,
                      float beta,
                      gsl_matrix_complex_float * C);

int  gsl_blas_zhemm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

int  gsl_blas_zherk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     double alpha,
                     const gsl_matrix_complex * A,
                     double beta,
                     gsl_matrix_complex * C);

int  gsl_blas_zher2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex alpha,
                      const gsl_matrix_complex * A,
                      const gsl_matrix_complex * B,
                      double beta,
                      gsl_matrix_complex * C);

__END_DECLS

#endif /* __GSL_BLAS_H__ */


text:
/*

 */

#include "kern.h"

#include <cblas.h>
#include <err.h>
#include <fcntl.h>
#include <float.h>
#include <math.h>
#include <stdbool.h>
#include <stdlib.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <time.h>
#include <unistd.h>

#ifdef WITH_THREADS
#include <omp.h>
#define MP_LOOP() PRAGMA(omp for simd)
#else
#define MP_LOOP()
#endif

#define PRAGMA(X) _Pragma(#X)

/*
 * Transform the input vector to the output stream.
 * The m x n matrix is stored with column arrays.
 *
 * y[k * n + j] = x[k * m + i] * w[m * j + i]
 * y(l, n) = x(l,m) * w(m, n)
 *
 */
void trans(uint32_t batch_len, uint32_t m, uint32_t n, const float *w,
           const float *x, float *y) {
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, batch_len, n, m,
                1.0f, x, batch_len, w, m, 0.0f, y, batch_len);
}

/*
 * Original:
 * w(m,n) -= N * x(l,m)^T * dy(l,n)
 * w[m * $j + i] -= N * x[m * $k + i] * dy[n * k + $j]
 *
 * Optimized:
 * w(m,n)^T -= N * dy(l,n)^T * x(l,m)
 * w[m * j + $i] -= N * dy[n * $k + j] * x[m * k + $i]
 */
void train_sgd(uint32_t batch_len, uint32_t m, uint32_t n,
               const float x[m * batch_len], const float y[n * batch_len],
               float rate, float w[m * n]) {
    cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, n, m, batch_len, -rate,
                y, n, x, m, 1.0f, w, m);
}

/*
 * adam optimizer for the weight matrix
 */
float train_adam(uint32_t batch_len, uint32_t m, uint32_t n,
                 const float x[restrict m * batch_len],
                 const float dy[restrict const n * batch_len], float counter,
                 float N, float beta1, float beta2, float epsilon,
                 float w[m * n], float mom[n * m], float veloc[n * m]) {
    float *restrict wr;
    float *restrict mr;
    float *restrict vr;
    const float *restrict xr;
    float dr;
    uint32_t k, i;
    for (uint32_t j = 0; j < n; j++) {
        for (k = 0, wr = &w[m * j], vr = &veloc[j * m], mr = &mom[j * m];
             k < batch_len; k++) {
            for (i = 0, xr = &x[k * m], dr = dy[k * n + j]; i < m; i++) {
                const float g      = dr * xr[i];
                mr[i]              = beta1 * mr[i] + ((1 - beta1) * g);
                const float mr_hat = mr[i] / (1 - powf(beta1, counter));
                vr[i] = beta2 * vr[i] + ((1 - beta2) * (powf(g, 2)));
                const float vr_hat = vr[i] / (1 - powf(beta2, counter));
                wr[i] -= N * mr_hat / (sqrtf(vr_hat + epsilon));
            }
        }
    }
    return counter + 1.0f;
}

/*
 * Original:
 * dx(l,m) = dy(l,n) * w(m,n)^T
 * dx[m * k + $i] = dy[n * k + $j] * w[m * $j + i]
 */
void loss(uint32_t batch_len, uint32_t m, uint32_t n, const float *w,
          const float *dy, float *dx) {
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, batch_len, m, n, 1.0f,
                dy, n, w, m, .0f, dx, m);
}

/*
 * Calculate the difference between two vectors of size 'size'.
 *
 * size:   The size of the vectors.
 * vec1:   The first vector
 * vec2:   The second vector
 *
 * Return the mean difference
 */
double vec_delta(uint32_t size, const float *vec1, const float *vec2,
                 float *deltas) {
    double error = 0.0;
    for (size_t i = 0; i < size; i++) {
        deltas[i] = vec1[i] - vec2[i];
        error += pow(deltas[i], 2.0);
    }
    return error / (float)size;
}

bool vec_is_equal_f32(uint32_t n, const float a[n], const float b[n],
                      float epsilon) {
    for (uint32_t i = 0; i < n; i++) {
        if (fabsf(a[i] - b[i]) >= epsilon) return false;
        // better:
        // if(fabs(a[ii]-b[ii]) < 1e-10 * (fabs(a[ii]) + fabs(b[ii]))) {
        // with the appropriate tolerance
    }
    return true;
}

/*
 * In each iteration of the while loop two normal random variables are
 * generated. On the first call of the function, two normal random variables are
 * generated.
 * On the after call, the second generated number will be returned.
 */
float rand_normal(float mu, float sigma) {
    float U1, U2, W, mult;
    static float X1, X2;
    static bool call = false;

    if (call == true) {
        call = !call;
        return mu + sigma * X2;
    }

    do {
        U1 = -1 + ((float)random() / (float)(RAND_MAX)) * 2;
        U2 = -1 + ((float)random() / (float)(RAND_MAX)) * 2;
        W  = powf(U1, 2) + powf(U2, 2);
    } while (W >= 1 || W == 0);

    mult = sqrtf((-2 * logf(W)) / W);
    X1   = U1 * mult;
    X2   = U2 * mult;

    call = !call;
    return (mu + sigma * X1);
}

void weights_norm_init(uint32_t in_size, uint32_t out_size, float *weights) {
    srandom(time(NULL));
    for (unsigned long long i = 0; i < (in_size * out_size); i++) {
        weights[i] = rand_normal(0.0f, sqrtf(2.0f / in_size));
    }
}

uint32_t argmax(uint32_t len, const float x[len], float *max) {
    *max             = x[0];
    uint32_t max_pos = 0;
    for (uint32_t i = 0; i < len; i++) {
        if (x[i] > *max) {
            *max    = x[i];
            max_pos = i;
        }
    }
    return max_pos;
}

void softmax(uint32_t len, const float x[len], float xs[len]) {
    float max;
    argmax(len, x, &max);
    float sum = .0f;
    for (uint32_t i = 0; i < len; i++) {
        xs[i] = expf(x[i] - max);
        sum += xs[i];
    }
    for (uint32_t i = 0; i < len; i++) {
        xs[i] = xs[i] / sum;
    }
}

// List of activation functions
// Convention: derived activation functions are prefixed with label 'derived'
/*
 * The sigmoid function
 */
float fsigmoid(float x) { return 1.0f / (1.0f + expf(-x)); }

/*
 * The derived of the sigmoid function
 */
float derived_fsigmoid(float y) { return (y * (1.0f - y)); }

float frelu(float x) { return (float)(x > .0) * x; }

float derived_frelu(float x) { return (float)(x > .0); }

float ftanh(float x) { return tanhf(x); }

float derived_ftanh(float x) { return 1.0f - powf(ftanh(x), 2.0f); }

static void vec_func_f32(uint32_t len, float result[len], float (*f)(float)) {
    for (uint32_t d = 0; d < len; d++) {
        result[d] = f(result[d]);
    }
}

void sigmoid(uint32_t len, float *result) {
    vec_func_f32(len, result, fsigmoid);
}

void relu(uint32_t len, float *result) { vec_func_f32(len, result, frelu); }

void tanhg(uint32_t len, float *result) { vec_func_f32(len, result, ftanh); }

static void vec_derived_f32(uint32_t len, const float result[len],
                            float delta[len], float (*f)(float)) {
    for (uint32_t d = 0; d < len; d++) {
        delta[d] = f(result[d]) * delta[d];
    }
}

void relu_derived(uint32_t len, const float *result, float *delta) {
    vec_derived_f32(len, result, delta, derived_frelu);
}

void tanhg_derived(uint32_t len, const float *result, float *delta) {
    vec_derived_f32(len, result, delta, derived_ftanh);
}

void sigmoid_derived(uint32_t len, const float *result, float *delta) {
    vec_derived_f32(len, result, delta, derived_fsigmoid);
}

/*
 * Allocate memory for a float matrix.
 */
float *matrix_alloc(uint32_t m, uint32_t n) {
    return reallocarray(NULL, m * n, sizeof(float));
}

void matrix_init(uint32_t m, uint32_t n, float matrix[m * n]) {
    for (uint32_t i = 0; i < m * n; i++) {
        matrix[i] = 0.0f;
    }
}

float *weights_create_or_load(const char *filename, uint32_t input_len,
                              uint32_t output_len) {
    float *weight;
    if (filename == NULL) {
        weight = matrix_alloc(input_len, output_len);
        weights_norm_init(input_len, output_len, weight);
    } else {
        struct stat statbuf;
        size_t size_expected = input_len * output_len * sizeof(float);

        int fd = open(filename, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR);
        if (fd < 0) {
            /* failure */
            err(EXIT_FAILURE, "open file '%s'", filename);
        }
        /* find size of input file */
        if (fstat(fd, &statbuf) < 0) {
            err(EXIT_FAILURE, "fstat error");
        } else if (statbuf.st_size == 0) {
            // file does not exist before - create it with the required size
            if (ftruncate(fd, size_expected)) {
                err(EXIT_FAILURE, "file truncate");
            }
            weight = mmap(0, size_expected, PROT_WRITE, MAP_SHARED, fd, 0);
            close(fd);
            weights_norm_init(input_len, output_len, weight);
            return weight;
        } else if (statbuf.st_size > 0 &&
                   (unsigned long)statbuf.st_size < size_expected) {
            errx(EXIT_FAILURE, "invalid data size. Expected: %lu; given: %ld",
                 size_expected, statbuf.st_size);
        }
        weight = mmap(0, size_expected, PROT_WRITE, MAP_SHARED, fd, 0);
        close(fd);
    }
    return weight;
}

static inline float fbernoulli(float p /* must between [0,1] */) {
    return (float)(p < (((float)random() / (float)(RAND_MAX))));
}

void dropout(uint32_t len, const float vec[len], float p, float result[len]) {
    for (uint32_t i = 0; i < len; i++) {
        result[i] = vec[i] * fbernoulli(p);
    }
}

max_stars_repo_name: geisten/gstnn

max_issues_repo_name: geisten/gstnn

max_forks_repo_name: geisten/gstnn


text:
#ifndef DART_CPP14_SHIM_H
#define DART_CPP14_SHIM_H

// Figure out what compiler we have.
#if defined(__clang__)
#define DART_USING_CLANG 1
#elif defined(__GNUC__) || defined(__GNUG__)
#define DART_USING_GCC 1
#elif defined(_MSC_VER)
#define DART_USING_MSVC 1
#endif

#ifdef DART_USING_MSVC
#define _CRT_SECURE_NO_WARNINGS 1
#endif

#include <ctime>
#include <cstring>

#if DART_USING_CLANG && __clang_major__ >= 5 && __clang_major__ <= 7
// Side-step a disagreement between clang (5/6) and GNU std::variant.
#define DART_USE_MPARK_VARIANT 1
#elif defined(__APPLE__)
// Side-step AppleClang misreporting compiler capabilities on macos 10.13 and below.
#include <availability.h>
#ifndef __MAC_10_14
#define DART_USE_MPARK_VARIANT 1
#endif
#endif

// MSVC doesn't have a signed size type, obviously
#if DART_USING_MSVC
#include <BaseTsd.h>
using ssize_t = SSIZE_T;
#endif

// Make sure we have a fallback for compilers that don't support attributes at all.
#ifndef __has_cpp_attribute
#define __has_cpp_attribute(name) 0
#endif

// Figure out how to declare things [[nodiscard]]
#if __has_cpp_attribute(gnu::warn_unused_result)
#define DART_NODISCARD [[gnu::warn_unused_result]]
#elif __has_cpp_attribute(nodiscard)
#define DART_NODISCARD [[nodiscard]]
#else
#define DART_NODISCARD
#endif

#define DART_STRINGIFY_IMPL(x) #x
#define DART_STRINGIFY(x) DART_STRINGIFY_IMPL(x)

#if DART_USING_MSVC
#define DART_UNLIKELY(x) !!(x)
#else
#define DART_UNLIKELY(x) __builtin_expect(!!(x), 0)
#endif

#ifndef NDEBUG

#if DART_USING_MSVC
#include <io.h>
#define DART_WRITE(fd, ptr, bytes) _write(fd, ptr, static_cast<unsigned int>(bytes))
#define DART_STDERR_FILENO _fileno(stderr)
#else
#include <unistd.h>
#define DART_WRITE(fd, ptr, bytes) write(fd, ptr, bytes)
#define DART_STDERR_FILENO STDERR_FILENO
#endif

/**
 *  @brief
 *  Macro customizes functionality usually provided by assert().
 *
 *  @details
 *  Not strictly necessary, but tries to provide a bit more context and
 *  information as to why I just murdered the user's program (in production, no doubt).
 *
 *  @remarks
 *  Don't actually know if Doxygen lets you document macros, guess we'll see.
 */
#define DART_ASSERT(cond)                                                                                \
  if (DART_UNLIKELY(!(cond))) {                                                                               \
    auto& msg = "dart::packet has detected fatal memory corruption and cannot continue execution.\n"          \
      "\"" DART_STRINGIFY(cond) "\" violated.\nSee " __FILE__ ":" DART_STRINGIFY(__LINE__) "\n";              \
    int errfd = DART_STDERR_FILENO;                                                                           \
    ssize_t spins {0}, written {0}, total {sizeof(msg)};                                                      \
    do {                                                                                \
      ssize_t ret = DART_WRITE(errfd, msg + written, total - written);                                        \
      if (ret >= 0) written += ret;                                                                           \
    } while (written != total && spins++ < 16);                                                               \
    std::abort();                                                                                \
  }
#else
#define DART_ASSERT(cond)
#endif

// Conditionally include different implementations of different data structures
// depending on what standard we have access to.
#if __cplusplus >= 201703L && !DART_USE_MPARK_VARIANT
#define DART_HAS_CPP17 1
#include <variant>
#include <optional>
#include <string_view>
#else
#define DART_HAS_CPP14 1
#include "support/variant.h"
#include "support/optional.h"
#include "support/string_view.h"
#endif

// We support two versions of GSL, but unfortunately they don't agree about
// the template signature of gsl::span (why would we expect two production-ready
// implementations of the same exact specification to agree about something like
// that), so we have to declare some of our partial specializations differently
// depending on what we included.
#include <gsl/gsl>
#ifndef gsl_lite_VERSION
#define DART_USING_GSL
#else
#define DART_USING_GSL_LITE
#endif

// Conditionally pull each of those types into our namespace.
namespace dart {
  namespace shim {

#if DART_USING_MSVC
    inline int aligned_alloc(void** memptr, size_t alignment, size_t size) {
      void* ret = _aligned_malloc(size, alignment);
      if (ret) {
        *memptr = ret;
        return 0;
      } else {
        return -1;
      }
    }

    inline void aligned_free(void* ptr) {
      _aligned_free(ptr);
    }

    inline void gmtime(time_t const* src, std::tm* out) {
      gmtime_s(out, src);
    }

    inline time_t timegm(std::tm* src) {
      return _mkgmtime(src);
    }
#else
    inline int aligned_alloc(void** memptr, size_t alignment, size_t size) {
      return posix_memalign(memptr, alignment, size);
    }

    inline void aligned_free(void* ptr) {
      free(ptr);
    }

    inline void gmtime(time_t const* src, std::tm* out) {
      gmtime_r(src, out);
    }

    inline time_t timegm(std::tm* src) {
      return ::timegm(src);
    }
#endif

#ifdef DART_HAS_CPP17

    // Pull in names of types.
    using std::optional;
    using std::nullopt_t;
    using std::variant;
    using std::monostate;
    using std::string_view;
    using std::basic_string_view;

    // Pull in constants.
    static inline constexpr auto nullopt = std::nullopt;

    // Pull in non-member helpers.
    using std::get;
    using std::visit;
    using std::get_if;
    using std::launder;
    using std::holds_alternative;
    using std::variant_alternative;
    using std::variant_alternative_t;

    // Define a way to compose lambdas.
    template <class... Ls>
    struct compose : Ls... {
      using Ls::operator ()...;
    };
    template <class... Ls>
    compose(Ls...) -> compose<Ls...>;

    template <class... Ls>
    auto compose_together(Ls&&... lambdas) {
      return compose {std::forward<Ls>(lambdas)...};
    }
#else

    // Pull in names of types.
    using dart::optional;
    using dart::nullopt_t;
    using dart::string_view;
    using dart::basic_string_view;
    using mpark::variant;
    using mpark::monostate;

    // Pull in non-member helpers.
    using mpark::get;
    using mpark::visit;
    using mpark::get_if;
    using mpark::holds_alternative;
    using mpark::variant_alternative;
    using mpark::variant_alternative_t;

    // Pull in constants.
    static constexpr auto nullopt = dart::nullopt;

    // Define a way to compose lambdas.
    template <class... Ls>
    struct compose;
    template <class L, class... Ls>
    struct compose<L, Ls...> : L, compose<Ls...> {
      compose(L l, Ls... the_rest) : L(std::move(l)), compose<Ls...>(std::move(the_rest)...) {}
      using L::operator ();
      using compose<Ls...>::operator ();
    };
    template <class L>
    struct compose<L> : L {
      compose(L l) : L(std::move(l)) {}
      using L::operator ();
    };

    template <class... Ls>
    auto compose_together(Ls&&... lambdas) {
      return compose<std::decay_t<Ls>...> {std::forward<Ls>(lambdas)...};
    }
#endif
  }
}

#endif


text:
/* combination/test.c
 * based on permutation/test.c 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_combination.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

size_t c63[20][3] = {
  { 0, 1, 2 },  { 0, 1, 3 },  { 0, 1, 4 },  { 0, 1, 5 },
  { 0, 2, 3 },  { 0, 2, 4 },  { 0, 2, 5 },  { 0, 3, 4 },
  { 0, 3, 5 },  { 0, 4, 5 },  { 1, 2, 3 },  { 1, 2, 4 },
  { 1, 2, 5 },  { 1, 3, 4 },  { 1, 3, 5 },  { 1, 4, 5 },
  { 2, 3, 4 },  { 2, 3, 5 },  { 2, 4, 5 },  { 3, 4, 5 }
} ;

void my_error_handler (const char *reason, const char *file, int line, int err);

int
main (void)
{
  size_t i, j;
  int status = 0, s;
  gsl_combination * c ;

  gsl_ieee_env_setup ();

  c = gsl_combination_alloc (6,3);

  /* Test combinations in forward order */

  gsl_combination_init_first (c);

  i = 0;

  do
    {
      if ( i >= 20 )
        {
          status = 1;
          break;
        }
      for (j = 0; j < 3; j++)
        {
          status |= (c->data[j] != c63[i][j]);
        }

      {
        int s1 = gsl_combination_valid (c);
        gsl_test (s1, "gsl_combination_valid (%u)", i);
      }

      i++;
    }
  while (gsl_combination_next(c) == GSL_SUCCESS);

  gsl_test(status, "gsl_combination_next, 6 choose 3 combination, 20 steps");

  gsl_combination_next(c);
  gsl_combination_next(c);
  gsl_combination_next(c);
  for (j = 0; j < 3; j++)
    {
      status |= (c->data[j] != c63[19][j]);
    }
  gsl_test(status, "gsl_combination_next on the last combination");

  {
    int s1 = gsl_combination_valid (c);
    gsl_test (s1, "gsl_combination_valid on the last combination");
  }

  {
    gsl_combination * d = gsl_combination_alloc (6,3);
    gsl_combination_memcpy (d, c);

    status = 0;

    for (j = 0; j < 3; j++)
      {
        status |= (d->data[j] != c->data[j]);
      }

    gsl_test (status, "gsl_combination_memcpy, 6 choose 3 combination");
    gsl_combination_free(d);
  }

  /* Now test combinations in reverse order */

  gsl_combination_init_last (c);

  i = 20;
  do
    {
      if ( i == 0 )
        {
          status = 1;
          break;
        }

      i--;

      for (j = 0; j < 3; j++)
        {
          status |= (c->data[j] != c63[i][j]);
        }

      {
        int s1 = gsl_combination_valid (c);
        gsl_test (s1, "gsl_combination_valid (%u)", i);
      }
    }
  while (gsl_combination_prev(c) == GSL_SUCCESS);

  gsl_test(status, "gsl_combination_prev, 6 choose 3 combination, 20 steps");

  gsl_combination_prev(c);
  gsl_combination_prev(c);
  gsl_combination_prev(c);
  for (j = 0; j < 3; j++)
    {
      status |= (c->data[j] != c63[0][j]);
    }
  gsl_test(status, "gsl_combination_prev on the first combination");

  {
    int s1 = gsl_combination_valid (c);
    gsl_test (s1, "gsl_combination_valid on the first combination");
  }

  {
    gsl_combination * d = gsl_combination_alloc (6,3);
    gsl_combination_memcpy (d, c);

    status = 0;

    for (j = 0; j < 3; j++)
      {
        status |= (d->data[j] != c->data[j]);
      }

    gsl_test (status, "gsl_combination_memcpy, 6 choose 3 combination");
    gsl_combination_free(d);
  }

  gsl_combination_free (c);

  c = gsl_combination_calloc(7, 0);
  /* should return GSL_FAILURE every time */
  status |= (gsl_combination_next(c) != GSL_FAILURE);
  status |= (gsl_combination_next(c) != GSL_FAILURE);
  status |= (gsl_combination_prev(c) != GSL_FAILURE);
  status |= (gsl_combination_prev(c) != GSL_FAILURE);
  gsl_test(status, "gsl_combination 7 choose 0");
  gsl_combination_free (c);

  c = gsl_combination_calloc(7, 7);
  /* should return GSL_FAILURE every time */
  for(j = 0; j < 7; j++)
  {
    status |= (gsl_combination_get(c, j) != j);
  }
  status |= (gsl_combination_next(c) != GSL_FAILURE);
  for(j = 0; j < 7; j++)
  {
    status |= (gsl_combination_get(c, j) != j);
  }
  status |= (gsl_combination_next(c) != GSL_FAILURE);
  for(j = 0; j < 7; j++)
  {
    status |= (gsl_combination_get(c, j) != j);
  }
  status |= (gsl_combination_prev(c) != GSL_FAILURE);
  for(j = 0; j < 7; j++)
  {
    status |= (gsl_combination_get(c, j) != j);
  }
  status |= (gsl_combination_prev(c) != GSL_FAILURE);
  for(j = 0; j < 7; j++)
  {
    status |= (gsl_combination_get(c, j) != j);
  }
  gsl_test(status, "gsl_combination 7 choose 7");
  gsl_combination_free (c);

  c = gsl_combination_calloc(6, 3);

  gsl_set_error_handler (&my_error_handler);

  c->data[0] = 1;
  c->data[1] = 1;
  c->data[2] = 2;
  s = gsl_combination_valid (c);
  gsl_test (!s, "gsl_combination_valid on an invalid combination (1,1,2)");

  c->data[0] = 2;
  c->data[1] = 1;
  c->data[2] = 0;
  s = gsl_combination_valid (c);
  gsl_test (!s, "gsl_combination_valid on an invalid combination (2,1,0)");

  c->data[0] = 1;
  c->data[1] = 2;
  c->data[2] = 0;
  s = gsl_combination_valid (c);
  gsl_test (!s, "gsl_combination_valid on an invalid combination (1,2,0)");

  {
    gsl_combination * d = gsl_combination_alloc (6,4);
    int s = gsl_combination_memcpy (d, c);
    gsl_test (!s, "gsl_combination_memcpy, (6,4) vs (6,3)");
    gsl_combination_free(d);
  }

  {
    gsl_combination * d = gsl_combination_alloc (7,3);
    int s = gsl_combination_memcpy (d, c);
    gsl_test (!s, "gsl_combination_memcpy, (7,3) vs (6,3)");
    gsl_combination_free(d);
  }

  {
    gsl_combination * d = gsl_combination_alloc (7,2);
    int s = gsl_combination_memcpy (d, c);
    gsl_test (!s, "gsl_combination_memcpy, (7,2) vs (6,3)");
    gsl_combination_free(d);
  }

  gsl_combination_free (c);

  exit (gsl_test_summary());
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc