text:
#include "asf.h"
#include "asf_meta.h"
#include <unistd.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

struct refine_offset_params {
        meta_parameters *meta;
        double x_off, y_off;
};

static double err_at_pixel(struct refine_offset_params *p, double off_t,
                           double off_x, double line, double samp)
{
    const double huge_err = 999999.9;

    meta_parameters *meta = p->meta;
    double line_out, samp_out;
    double lat, lon;
    int bad=0;

    double saved_timeOffset = meta->sar->time_shift;
    double saved_slant = meta->sar->slant_shift;

    // converting the pixel's location back to lat/lon is done before
    // applying the prospective slant and range shifts.
    bad = meta_get_latLon(meta, line, samp, 0.0, &lat, &lon);
    if (bad) { printf("bad: meta_get_latLon\n"); return huge_err; }

    meta->sar->time_shift += off_t;
    meta->sar->slant_shift += off_x;

    bad = meta_get_lineSamp(meta, lat, lon, 0.0, &line_out, &samp_out);

    meta->sar->time_shift = saved_timeOffset;
    meta->sar->slant_shift = saved_slant;

    // meta_get_lineSamp returns a non-zero value if, for the specified
    // lat&lon, we weren't able to determine a line & sample.  In this case,
    // we just return a giant error value.
    if (bad) { printf("bad: meta_get_lineSamp\n"); return huge_err; }

    //printf("trying %f,%f\n",off_t,off_x);
    //printf("at %7.3f,%7.3f,%7.3f,%7.3f,%7.3f,%7.3f,%7.3f,%7.3f -> %7.3f\n",
    //       off_t, off_x, line, samp, line_out, samp_out, p->y_off, p->x_off,
    //       fabs(line-line_out-p->y_off)+fabs(samp-samp_out-p->x_off));
    double xdiff = line_out - line - p->y_off;
    double ydiff = samp_out - samp - p->x_off;
    return sqrt(xdiff*xdiff + ydiff*ydiff);
}

static int
getObjective(const gsl_vector *x, void *params, gsl_vector *f)
{
    double off_t = gsl_vector_get(x,0);
    double off_x = gsl_vector_get(x,1);

    if (!meta_is_valid_double(off_x) ||
        !meta_is_valid_double(off_t)) {
        // This does happen sometimes, when we've already found the root
	return GSL_FAILURE;
    }

    struct refine_offset_params *p = (struct refine_offset_params *)params;
    meta_parameters *meta = p->meta;

    int nl = meta->general->line_count;
    int ns = meta->general->sample_count;

    double err = 0.0;

    err += err_at_pixel(p, off_t, off_x, nl/8, ns/8);
    err += err_at_pixel(p, off_t, off_x, 7*nl/8, ns/8);
    err += err_at_pixel(p, off_t, off_x, nl/8, 7*ns/8);
    err += err_at_pixel(p, off_t, off_x, 7*nl/8, 7*ns/8);

    err += err_at_pixel(p, off_t, off_x, nl/8, ns/2);
    err += err_at_pixel(p, off_t, off_x, 7*nl/8, ns/2);
    err += err_at_pixel(p, off_t, off_x, nl/2, 7*ns/8);
    err += err_at_pixel(p, off_t, off_x, nl/2, ns/8);

    err += err_at_pixel(p, off_t, off_x, nl/2, ns/2);

    gsl_vector_set(f,0,err);
    gsl_vector_set(f,1,err);
    return GSL_SUCCESS;
}

/*
static void print_state(int iter, gsl_multiroot_fsolver *s)
{
    printf("iter = %3d   x = (%.3f %.3f)   f(x) = %.8f\n", iter,
           gsl_vector_get(s->x, 0), gsl_vector_get(s->x, 1),
           gsl_vector_get(s->f, 0));
}
*/

static void coarse_search(double t_extent_min, double t_extent_max,
                          double x_extent_min, double x_extent_max,
                          double *t_min, double *x_min,
                          struct refine_offset_params *params)
{
    double the_min = 9999999;
    double min_t=99, min_x=99;
    int i,j,k=6;
    double t_extent = t_extent_max - t_extent_min;
    double x_extent = x_extent_max - x_extent_min;
    gsl_vector *v = gsl_vector_alloc(2);
    gsl_vector *u = gsl_vector_alloc(2);
    //printf("           ");
    //for (j = 0; j <= k; ++j) {
    //    double x = x_extent_min + ((double)j)/k*x_extent;
    //    printf("%9.3f ", x);
    //}
    //printf("\n           ");
    //for (j = 0; j <= k; ++j)
    //    printf("--------- ");
    //printf("\n");
    for (i = 0; i <= k; ++i) {
        double t = t_extent_min + ((double)i)/k*t_extent;
        //printf("%9.3f | ", t);

        for (j = 0; j <= k; ++j) {
            double x = x_extent_min + ((double)j)/k*x_extent;

            gsl_vector_set(v, 0, t);
            gsl_vector_set(v, 1, x);
            getObjective(v,(void*)params, u);
            double n = gsl_vector_get(u,0);
            //printf("%9.3f ", n);
            if (n<the_min) {
                the_min=n;
                min_t=gsl_vector_get(v,0);
                min_x=gsl_vector_get(v,1);
            }
        }
        //printf("\n");
    }

    *t_min = min_t;
    *x_min = min_x;

    gsl_vector_free(v);
    gsl_vector_free(u);
}

static void generate_start(struct refine_offset_params *params,
                           double *start_t, double *start_x)
{
    int i;

    double extent_t_min = -10;
    double extent_t_max = 10;

    double extent_x_min = -5000;
    double extent_x_max = 5000;

    double t_range = extent_t_max - extent_t_min;
    double x_range = extent_x_max - extent_x_min;

    for (i = 0; i < 12; ++i)
    {
        coarse_search(extent_t_min, extent_t_max, extent_x_min, extent_x_max,
                      start_t, start_x, params);

        t_range /= 3;
        x_range /= 3;

        extent_t_min = *start_t - t_range/2;
        extent_t_max = *start_t + t_range/2;

        extent_x_min = *start_x - x_range/2;
        extent_x_max = *start_x + x_range/2;

        //printf("refining search to region: (%9.3f,%9.3f)\n"
        //       "                           (%9.3f,%9.3f)\n",
        //       extent_t_min, extent_t_max,
        //       extent_x_min, extent_x_max);
    }

}

void refine_offset(double x_off, double y_off, meta_parameters *meta,
                   double *out_t, double *out_x)
{
    int status;
    int iter = 0, max_iter = 1000;
    const gsl_multiroot_fsolver_type *T;
    gsl_multiroot_fsolver *s;
    gsl_error_handler_t *prev;
    struct refine_offset_params params;
    const size_t n = 2;

    params.meta = meta;
    params.x_off = x_off;
    params.y_off = y_off;

    gsl_multiroot_function F = {&getObjective, n, &params};
    gsl_vector *x = gsl_vector_alloc(n);

    double start_t, start_x;
    generate_start(&params, &start_t, &start_x);
    //printf("Starting point at (%f,%f)\n", start_t, start_x);

    gsl_vector_set (x, 0, start_t);
    gsl_vector_set (x, 1, start_x);

    T = gsl_multiroot_fsolver_hybrid;
    s = gsl_multiroot_fsolver_alloc(T, n);
    gsl_multiroot_fsolver_set(s, &F, x);

    prev = gsl_set_error_handler_off();

    do {
        ++iter;
        status = gsl_multiroot_fsolver_iterate(s);

        //print_state(iter, s);

        // abort if stuck
        if (status) break;

        status = gsl_multiroot_test_residual (s->f, 1e-8);
    } while (status == GSL_CONTINUE && iter < max_iter);

    *out_t = gsl_vector_get(s->x, 0);
    *out_x = gsl_vector_get(s->x, 1);

    gsl_vector *retrofit = gsl_vector_alloc(n);
    gsl_vector_set(retrofit, 0, *out_t);
    gsl_vector_set(retrofit, 1, *out_x);
    gsl_vector *output = gsl_vector_alloc(n);
    getObjective(retrofit, (void*)&params, output);
    //double val= gsl_vector_get(output,0);
    //printf("GSL Result: %f at (%f,%f)\n", val, *out_t, *out_x);
    gsl_vector_free(retrofit);
    gsl_vector_free(output);

    gsl_multiroot_fsolver_free(s);
    gsl_vector_free(x);
    gsl_set_error_handler(prev);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef NOMPI
#include <mpi.h>
#endif

#include <fcntl.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <sys/file.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>

#include "allvars.h"
#include "proto.h"

/*! \file restart.c
 *  \brief Code for reading and writing restart files
 */

static FILE *fd;

static void in(int *x, int modus);
static void byten(void *x, size_t n, int modus);

/*! This function reads or writes the restart files.  Each processor writes
 *  its own restart file, with the I/O being done in parallel. To avoid
 *  congestion of the disks you can tell the program to restrict the number of
 *  files that are simultaneously written to NumFilesWrittenInParallel.
 *
 *  If modus>0 the restart()-routine reads, if modus==0 it writes a restart
 *  file.
 */
void restart(int modus)
{
  char buf[200], buf_bak[200], buf_mv[500];
  double save_PartAllocFactor, save_TreeAllocFactor;
  int i, nprocgroup, masterTask, groupTask, old_MaxPart, old_MaxNodes;
  struct global_data_all_processes all_task0;

  sprintf(buf, "%s%s.%d", All.OutputDir, All.RestartFile, ThisTask);
  sprintf(buf_bak, "%s%s.%d.bak", All.OutputDir, All.RestartFile, ThisTask);
  sprintf(buf_mv, "mv %s %s", buf, buf_bak);

  if((NTask < All.NumFilesWrittenInParallel))
    {
      printf
	("Fatal error.\nNumber of processors must be a smaller or equal than `NumFilesWrittenInParallel'.\n");
      endrun(2131);
    }

  nprocgroup = NTask / All.NumFilesWrittenInParallel;

  if((NTask % All.NumFilesWrittenInParallel))
    {
      nprocgroup++;
    }

  masterTask = (ThisTask / nprocgroup) * nprocgroup;

  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))	/* ok, it's this processor's turn */
	{
	  if(modus)
	    {
	      if(!(fd = fopen(buf, "r")))
		{
		  printf("Restart file '%s' not found.\n", buf);
		  endrun(7870);
		}
	    }
	  else
	    {
	      system(buf_mv);	/* move old restart files to .bak files */

	      if(!(fd = fopen(buf, "w")))
		{
		  printf("Restart file '%s' cannot be opened.\n", buf);
		  endrun(7878);
		}
	    }

	  save_PartAllocFactor = All.PartAllocFactor;
	  save_TreeAllocFactor = All.TreeAllocFactor;

	  /* common data  */
	  byten(&All, sizeof(struct global_data_all_processes), modus);

	  if(ThisTask == 0 && modus > 0)
	    all_task0 = All;

#ifndef NOMPI
	  if(modus > 0 && groupTask == 0)	/* read */
	    {
	      MPI_Bcast(&all_task0, sizeof(struct global_data_all_processes), MPI_BYTE, 0, GADGET_WORLD);
	    }
#endif
	  old_MaxPart = All.MaxPart;
	  old_MaxNodes = All.TreeAllocFactor * All.MaxPart;

	  if(modus)		/* read */
	    {
	      if(All.PartAllocFactor != save_PartAllocFactor)
		{
		  All.PartAllocFactor = save_PartAllocFactor;
		  All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask);
		  All.MaxPartSph = All.PartAllocFactor * (All.TotN_gas / NTask);
		  save_PartAllocFactor = -1;
		}

	      if(All.TreeAllocFactor != save_TreeAllocFactor)
		{
		  All.TreeAllocFactor = save_TreeAllocFactor;
		  save_TreeAllocFactor = -1;
		}

	      if(all_task0.Time != All.Time)
		{
		  printf("The restart file on task=%d is not consistent with the one on task=0\n", ThisTask);
		  fflush(stdout);
		  endrun(16);
		}

	      allocate_memory();
	    }

	  in(&NumPart, modus);

	  if(NumPart > All.MaxPart)
	    {
	      printf
		("it seems you have reduced(!) 'PartAllocFactor' below the value of %g needed to load the restart file.\n",
		 NumPart / (((double) All.TotNumPart) / NTask));
	      printf("fatal error\n");
	      endrun(22);
	    }

	  /* Particle data  */
	  byten(&P[0], NumPart * sizeof(struct particle_data), modus);

	  in(&N_gas, modus);

	  if(N_gas > 0)
	    {
	      if(N_gas > All.MaxPartSph)
		{
		  printf
		    ("SPH: it seems you have reduced(!) 'PartAllocFactor' below the value of %g needed to load the restart file.\n",
		     N_gas / (((double) All.TotN_gas) / NTask));
		  printf("fatal error\n");
		  endrun(222);
		}
	      /* Sph-Particle data  */
	      byten(&SphP[0], N_gas * sizeof(struct sph_particle_data), modus);
	    }

	  /* write state of random number generator */
	  byten(gsl_rng_state(random_generator), gsl_rng_size(random_generator), modus);

	  /* now store relevant data for tree */

	  if(modus)		/* read */
	    {
	      ngb_treeallocate(MAX_NGB);

	      force_treeallocate(All.TreeAllocFactor * All.MaxPart, All.MaxPart);
	    }

	  in(&Numnodestree, modus);

	  if(Numnodestree > MaxNodes)
	    {
	      printf
		("Tree storage: it seems you have reduced(!) 'PartAllocFactor' below the value needed to load the restart file (task=%d). "
		 "Numnodestree=%d  MaxNodes=%d\n", ThisTask, Numnodestree, MaxNodes);
	      endrun(221);
	    }

	  byten(Nodes_base, Numnodestree * sizeof(struct NODE), modus);
	  byten(Extnodes_base, Numnodestree * sizeof(struct extNODE), modus);

	  byten(Father, NumPart * sizeof(int), modus);

	  byten(Nextnode, NumPart * sizeof(int), modus);
	  byten(Nextnode + All.MaxPart, MAXTOPNODES * sizeof(int), modus);

	  byten(DomainStartList, NTask * sizeof(int), modus);
	  byten(DomainEndList, NTask * sizeof(int), modus);
	  byten(DomainTask, MAXTOPNODES * sizeof(int), modus);
	  byten(DomainNodeIndex, MAXTOPNODES * sizeof(int), modus);
	  byten(DomainTreeNodeLen, MAXTOPNODES * sizeof(FLOAT), modus);
	  byten(DomainHmax, MAXTOPNODES * sizeof(FLOAT), modus);
	  byten(DomainMoment, MAXTOPNODES * sizeof(struct DomainNODE), modus);

	  byten(DomainCorner, 3 * sizeof(double), modus);
	  byten(DomainCenter, 3 * sizeof(double), modus);
	  byten(&DomainLen, sizeof(double), modus);
	  byten(&DomainFac, sizeof(double), modus);
	  byten(&DomainMyStart, sizeof(int), modus);
	  byten(&DomainMyLast, sizeof(int), modus);

	  if(modus)		/* read */
	    if(All.PartAllocFactor != save_PartAllocFactor || All.TreeAllocFactor != save_TreeAllocFactor)
	      {
		for(i = 0; i < NumPart; i++)
		  Father[i] += (All.MaxPart - old_MaxPart);

		for(i = 0; i < NumPart; i++)
		  if(Nextnode[i] >= old_MaxPart)
		    {
		      if(Nextnode[i] >= old_MaxPart + old_MaxNodes)
			Nextnode[i] += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxPart);
		      else
			Nextnode[i] += (All.MaxPart - old_MaxPart);
		    }

		for(i = 0; i < Numnodestree; i++)
		  {
		    if(Nodes_base[i].u.d.sibling >= old_MaxPart)
		      {
			if(Nodes_base[i].u.d.sibling >= old_MaxPart + old_MaxNodes)
			  Nodes_base[i].u.d.sibling +=
			    (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			else
			  Nodes_base[i].u.d.sibling += (All.MaxPart - old_MaxPart);
		      }

		    if(Nodes_base[i].u.d.father >= old_MaxPart)
		      {
			if(Nodes_base[i].u.d.father >= old_MaxPart + old_MaxNodes)
			  Nodes_base[i].u.d.father += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			else
			  Nodes_base[i].u.d.father += (All.MaxPart - old_MaxPart);
		      }

		    if(Nodes_base[i].u.d.nextnode >= old_MaxPart)
		      {
			if(Nodes_base[i].u.d.nextnode >= old_MaxPart + old_MaxNodes)
			  Nodes_base[i].u.d.nextnode +=
			    (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			else
			  Nodes_base[i].u.d.nextnode += (All.MaxPart - old_MaxPart);
		      }
		  }

		for(i = 0; i < MAXTOPNODES; i++)
		  if(Nextnode[i + All.MaxPart] >= old_MaxPart)
		    {
		      if(Nextnode[i + All.MaxPart] >= old_MaxPart + old_MaxNodes)
			Nextnode[i + All.MaxPart] += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
		      else
			Nextnode[i + All.MaxPart] += (All.MaxPart - old_MaxPart);
		    }

		for(i = 0; i < MAXTOPNODES; i++)
		  if(DomainNodeIndex[i] >= old_MaxPart)
		    {
		      if(DomainNodeIndex[i] >= old_MaxPart + old_MaxNodes)
			DomainNodeIndex[i] += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
		      else
			DomainNodeIndex[i] += (All.MaxPart - old_MaxPart);
		    }
	      }

	  fclose(fd);
	}
      else			/* wait inside the group */
	{
#ifndef NOMPI
	  if(modus > 0 && groupTask == 0)	/* read */
	    {
	      MPI_Bcast(&all_task0, sizeof(struct global_data_all_processes), MPI_BYTE, 0, GADGET_WORLD);
	    }
#endif
	}

#ifndef NOMPI
	      MPI_Barrier(GADGET_WORLD);
#endif
    }
}

/*! reads/writes n bytes in restart routine
 */
void byten(void *x, size_t n, int modus)
{
  if(modus)
    my_fread(x, n, 1, fd);
  else
    my_fwrite(x, n, 1, fd);
}

/*! reads/writes one `int' variable in restart routine
 */
void in(int *x, int modus)
{
  if(modus)
    my_fread(x, 1, sizeof(int), fd);
  else
    my_fwrite(x, 1, sizeof(int), fd);
}

ext: c

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


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

int main()
{
	printf("PI: %lf\n", M_PI);
}

max_stars_repo_name: quchunguang/test


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

#pragma once

#ifndef xml_f13c2b07_0230_4967_8a16_0c276e19b8a2_h
#define xml_f13c2b07_0230_4967_8a16_0c276e19b8a2_h

#include <gslib/std.h>
#include <gslib/tree.h>

__gslib_begin__

struct xml_kvpair { string key, value; };
typedef xml_kvpair xml_attr;

struct xml_key:
    public xml_kvpair
{
    xml_key(const string& k) { key = k; }
    xml_key(const gchar* c) { key.assign(c); }
};

struct xml_kvp_hash
{
    size_t operator()(const xml_kvpair& kvp) const
    { return string_hash(kvp.key); }
};

struct xml_kvp_equalto
{
    bool operator()(const xml_kvpair& p1, const xml_kvpair& p2) const
    { return p1.key == p2.key; }
};

template<class _kvpvsl>
struct xml_kvp_search
{
    static string* find(_kvpvsl& v, const string& k)
    {
        auto i = v.find(xml_key(k));
        return i != v.end() ? &(const_cast<string&>(i->value)) : nullptr;
    }
    static const string* find(const _kvpvsl& v, const string& k)
    {
        auto i = v.find(xml_key(k));
        return i != v.end() ? &(i->value) : nullptr;
    }
};

struct __gs_novtable xml_node abstract
{
    typedef unordered_set<xml_kvpair, xml_kvp_hash, xml_kvp_equalto> kvplist;
    typedef kvplist::iterator iterator;
    typedef kvplist::const_iterator const_iterator;

public:
    virtual bool is_value() const = 0;
    virtual void add_attribute(const xml_attr& attr) = 0;
    virtual string* get_attribute(const string& k) = 0;
    virtual const string* get_attribute(const string& k) const = 0;
    virtual void set_attribute(const gchar* k, int len1, const gchar* v, int len2) = 0;
    virtual void set_attribute(const string& k, const string& v) = 0;
    virtual int get_attribute_count() const = 0;
    virtual const string& get_name() const = 0;
    virtual void set_name(const gchar* str, int len) = 0;
    virtual string to_string() const = 0;   /* for debug */

public:
    string* get_attribute(const gchar* k) { return get_attribute(string(k)); }
    const string* get_attribute(const gchar* k) const { return get_attribute(string(k)); }
    string* get_attribute(const gchar* k, int len) { return get_attribute(string(k, len)); }
    const string* get_attribute(const gchar* k, int len) const { return get_attribute(string(k, len)); }
    bool has_attributes() const { return get_attribute_count() != 0; }
    void set_name(const string& str) { set_name(str.c_str(), str.length()); }
};

struct xml_element:
    public xml_node
{
    string      key;
    kvplist     attrlist;

public:
    bool is_value() const override { return false; }
    void add_attribute(const xml_attr& attr) override { attrlist.insert(attr); }
    string* get_attribute(const string& k) override { return xml_kvp_search<kvplist>::find(attrlist, k); }
    const string* get_attribute(const string& k) const override { return xml_kvp_search<kvplist>::find(attrlist, k); }
    void set_attribute(const gchar* k, int len1, const gchar* v, int len2) override
    {
        assert(k && v);
        string* vstr = get_attribute(string(k, len1));
        if(vstr) {
            vstr->assign(v, len2);
            return;
        }
        xml_attr attr;
        attr.key.assign(k, len1);
        attr.value.assign(v, len2);
        add_attribute(attr);
    }
    void set_attribute(const string& k, const string& v) override
    {
        string* vstr = get_attribute(k);
        if(vstr) {
            vstr->assign(v);
            return;
        }
        xml_attr attr;
        attr.key = k;
        attr.value = v;
        add_attribute(attr);
    }
    int get_attribute_count() const override { return (int)attrlist.size(); }
    const string& get_name() const override { return key; }
    void set_name(const gchar* str, int len) override { key.assign(str, len); }
    void set_value(const string& v) { __super::set_name(v); }
    kvplist& get_attributes() { return attrlist; }
    const kvplist& const_attributes() const { return attrlist; }
    string to_string() const override;
};

struct xml_value:
    public xml_node
{
    string      value;

public:
    bool is_value() const override { return true; }
    void add_attribute(const xml_attr& attr) override { assert(!"error!"); }
    string* get_attribute(const string& k) override { return nullptr; }
    const string* get_attribute(const string& k) const override { return nullptr; }
    void set_attribute(const gchar* k, int len1, const gchar* v, int len2) override { assert(!"error!"); }
    void set_attribute(const string& k, const string& v) override { assert(!"error!"); }
    int get_attribute_count() const override { return 0; }
    const string& get_name() const override { return value; }
    void set_name(const gchar* str, int len) override { value.assign(str, len); }
    string to_string() const override;
};

struct xml_enum
{
    typedef xml_node::iterator iterator;
    typedef xml_node::const_iterator const_iterator;

    xml_element*  node;
    iterator    iter;

    xml_enum(xml_node* n)
    {
        assert(n && !n->is_value());
        node = static_cast<xml_element*>(n);
        iter = node->attrlist.begin();
    }
    xml_enum(xml_node* n, iterator i)
    {
        assert(n && !n->is_value());
        node = static_cast<xml_element*>(n);
        iter = i;
    }
    void rewind(iterator i) { iter = i; }
    const xml_attr& get_attribute() const { return *iter; }
    bool next()
    {
        if(iter == node->attrlist.end())
            return false;
        return ++ iter != node->attrlist.end();
    }
    template<class _lambda>
    void for_each(_lambda lam) { for_each(node->attrlist.end(), lam); }
    template<class _lambda>
    void for_each(iterator end, _lambda lam)
    {
        iterator i = iter;
        for( ; i != end; ++ i)
            lam(*i);
    }
    template<class _lambda>
    void const_for_each(_lambda lam) const { const_for_each(node->attrlist.end(), lam); }
    template<class _lambda>
    void const_for_each(const_iterator end, _lambda lam) const
    {
        const_iterator i = iter;
        for( ; i != end; ++ i)
            lam(*i);
    }
};

struct xml_const_enum
{
    typedef xml_node::iterator iterator;
    typedef xml_node::const_iterator const_iterator;

    const xml_element* node;
    const_iterator  iter;

    xml_const_enum(const xml_node* n)
    {
        assert(n && !n->is_value());
        node = static_cast<const xml_element*>(n);
        iter = node->attrlist.begin();
    }
    xml_const_enum(const xml_node* n, const_iterator i)
    {
        assert(n && !n->is_value());
        node = static_cast<const xml_element*>(n);
        iter = i;
    }
    void rewind(const_iterator i) { iter = i; }
    const xml_attr& get_attribute() const { return *iter; }
    bool next()
    {
        if(iter == node->attrlist.end())
            return false;
        return ++ iter != node->attrlist.end();
    }
    template<class _lambda>
    void const_for_each(_lambda lam) const { const_for_each(node->attrlist.end(), lam); }
    template<class _lambda>
    void const_for_each(const_iterator end, _lambda lam) const
    {
        const_iterator i = iter;
        for( ; i != end; ++ i)
            lam(*i);
    }
};

typedef _treenode_wrapper<xml_node> xml_wrapper;

class xmltree:
    public tree<xml_node, xml_wrapper>
{
public:
    enum encodesys
    {
        encodesys_unknown,
        encodesys_ascii,
        encodesys_mbc,          /* multi-byte compatible */
        encodesys_wstr,         /* wide string */
    };
    enum encode
    {
        encode_unknown,
        encode_ansi,
        encode_gb2312,
        encode_utf8,
        encode_utf16,
    };
    encodesys       _encodesys;
    encode          _encode;

protected:
    static const gchar* skip(const gchar* str);
    static const gchar* skip(const gchar* str, const gchar* skp);
    static void replace(string& str, const gchar* s1, const gchar* s2);
    static void filter_entity_reference(string& s);
    static void recover_entity_reference(string& s);
    static void write_attribute(string& str, const gchar* prefix, const xml_attr& attr, const gchar* postfix);
    static void close_write_item(string& str, const gchar* prefix, const xml_node* node, const gchar* postfix);
    static void begin_write_item(string& str, const gchar* prefix, const xml_node* node, const gchar* postfix);
    static void end_write_item(string& str, const gchar* prefix, const xml_node* node, const gchar* postfix);

protected:
    const gchar* check_header(const gchar* str, string& version, string& encoding);
    const gchar* read_attribute(const gchar* str, string& k, string& v);
    const gchar* read_item(const gchar* str, iterator p);
    void setup_encoding_info(const string& enc);
    void write_header(string& str) const;
    void write_item(string& str, const_iterator i) const;

public:
    xmltree();
    bool load(const gchar* filename);
    bool parse(const gchar* src);
    void output(string& str) const;
    void save(const gchar* filename) const;
    void set_encode(encode enc) { _encode = enc; }

public:
    /*
     * Unique path locater
     * It was designed to access the DOM easier by a specified string path.
     * Usage:
     * 1.basically, we use the string "xx/yy/zz" to specify a path, which means in current position, find a sub node named "xx",
     *   go to the sub node; then find "yy", then "zz", and so on.
     *   the separator could be "/", or backslash "\", as well.
     * 2.about the name selector "$".
     *   you might also write the path as "$xx/$yy/$zz", it was the same as above; the selector "$" was skippable.
     * 3.about the counter "@"
     *   you could add a counter into the path, like "$xx@3/yy/zz", in this case, the selector "$" was unskippable.
     *   you could make a nameless selection, like "$@3/yy/zz", which means in the first jump, it would always choose the third node.
     *   or simply wrote "$/yy/zz", which means it would always choose the first node in the first jump.
     * 4.about the condition specifier "#"
     *   you could add a condition content into the selector;
     *   a condition was defined like: "$xx#aa,bb=cc,!dd|!ee,ff~=gg",
     *   there were two major binary connectors of conditions, "," refer to logic "AND", and "|" refer to logic "OR";
     *   the priority of "OR" was higher than "AND".
     *   if you need a negative condition, use "!" to reverse the logic, and the "!" should always follow "|"
     *   if a condition content was like the above "aa", which means the selected node must have the attribute named "aa"
     *   if a condition content was wrote like "bb=cc", which means the selected node must have an attribute named "bb" and valued "cc"
     *   if a condition content was wrote like "ff~=gg", which means the selected node must have an attribute named "ff" and valued "gg",
     *   the comparison was case insensitive.
     * 5.the "$", "@", "#" could be used as a combo, and the announcement should always be in the specified order, $#@
     */
    iterator unique_path_locater(const gchar* ctlstr) { return unique_path_locater(get_root(), ctlstr); }
    static iterator unique_path_locater(iterator i, const gchar* ctlstr)
    {
        assert(ctlstr);
        string cpyctlstr(ctlstr);
        ctlstr = cpyctlstr.c_str(); /* for trap */
        ctlstr = upl_run_over(i, ctlstr);
        if(ctlstr && !ctlstr[0])
            return i;
        return iterator(nullptr);
    }

public:
    static const gchar* upl_run_once(iterator& i, const gchar* ctlstr);
    static const gchar* upl_run_over(iterator& i, const gchar* ctlstr)
    {
        while(ctlstr && ctlstr[0]) {
            const gchar* ctlret = upl_run_once(i, ctlstr);
            if(!ctlret)
                return ctlstr;
            ctlstr = ctlret;
        }
        return ctlstr;
    }
};

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>

#if defined(HAVE_IRIX_IEEE_INTERFACE)
/* don't test denormals on IRIX */
#else
#define TEST_DENORMAL 1
#endif

int
main (void)
{
  float zerof = 0.0f, minus_onef = -1.0f ;
  double zero = 0.0, minus_one = -1.0 ;

  /* Check for +ZERO (float) */

  {
    float f = 0.0f;
    const char mantissa[] = "00000000000000000000000";
    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = 0, sign is +");
    gsl_test_int (r.exponent, -127, "float x = 0, exponent is -127");
    gsl_test_str (r.mantissa, mantissa, "float x = 0, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "float x = 0, type is ZERO");
  }

  /* Check for -ZERO (float) */

  {
    float f = minus_onef;
    const char mantissa[] = "00000000000000000000000";
    gsl_ieee_float_rep r;

    while (f < 0) {
      f *= 0.1f;
    }

    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 1, "float x = -1*0, sign is -");
    gsl_test_int (r.exponent, -127, "float x = -1*0, exponent is -127");
    gsl_test_str (r.mantissa, mantissa, "float x = -1*0, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "float x = -1*0, type is ZERO");
  }

  /* Check for a positive NORMAL number (e.g. 2.1) (float) */

  {
    float f = 2.1f;
    const char mantissa[] = "00001100110011001100110";

    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = 2.1, sign is +");
    gsl_test_int (r.exponent, 1, "float x = 2.1, exponent is 1");
    gsl_test_str (r.mantissa, mantissa, "float x = 2.1, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = 2.1, type is NORMAL");
  }

  /* Check for a negative NORMAL number (e.g. -1.3304...) (float) */

  {
    float f = -1.3303577090924210f ;
    const char mantissa[] = "01010100100100100101001";

    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 1, "float x = -1.3304..., sign is -");
    gsl_test_int (r.exponent, 0, "float x = -1.3304..., exponent is 0");
    gsl_test_str (r.mantissa, mantissa, "float x = -1.3304..., mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "float x = -1.3304..., type is NORMAL");
  }

  /* Check for a large positive NORMAL number (e.g. 3.37e31) (float) */

  {
    float f = 3.37e31f;
    const char mantissa[] = "10101001010110101001001";
    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = 3.37e31, sign is +");
    gsl_test_int (r.exponent, 104, "float x = 3.37e31, exponent is 104");
    gsl_test_str (r.mantissa, mantissa, "float x = 3.37e31, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = 3.37e31, type is NORMAL");
  }

  /* Check for a small positive NORMAL number (e.g. 3.37e-31) (float) */

  {
    float f = 3.37e-31f;
    const char mantissa[] = "10110101011100110111011";

    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = 3.37e-31, sign is +");
    gsl_test_int (r.exponent, -102, "float x = 3.37e-31, exponent is -102");
    gsl_test_str (r.mantissa, mantissa, "float x = 3.37e-31, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "float x = 3.37e-31, type is NORMAL");
  }

  /* Check for FLT_MIN (smallest possible number that is not denormal) */

  {
    float f = 1.17549435e-38f;	/* FLT_MIN (float) */
    const char mantissa[] = "00000000000000000000000";
    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = FLT_MIN, sign is +");
    gsl_test_int (r.exponent, -126, "float x = FLT_MIN, exponent is -126");
    gsl_test_str (r.mantissa, mantissa, "float x = FLT_MIN, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = FLT_MIN, type is NORMAL");
  }

  /* Check for FLT_MAX (largest possible number that is not Inf) */

  {
    float f = 3.40282347e+38f;	/* FLT_MAX */
    const char mantissa[] = "11111111111111111111111";

    gsl_ieee_float_rep r;
    gsl_ieee_float_to_rep (&f, &r);

    gsl_test_int (r.sign, 0, "float x = FLT_MAX, sign is +");
    gsl_test_int (r.exponent, 127, "float x = FLT_MAX, exponent is 127");
    gsl_test_str (r.mantissa, mantissa, "float x = FLT_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = FLT_MAX, type is NORMAL");
  }

  /* Check for DENORMAL numbers (e.g. FLT_MIN/2^n) */

#ifdef TEST_DENORMAL
  {
    float f = 1.17549435e-38f;	/* FLT_MIN */
    char mantissa[] = "10000000000000000000000";

    int i;
    gsl_ieee_float_rep r;

    for (i = 0; i < 23; i++)
      {
	float x = f / (float)pow (2.0, 1 + (float) i);
	mantissa[i] = '1';
	gsl_ieee_float_to_rep (&x, &r);

	gsl_test_int (r.sign, 0, "float x = FLT_MIN/2^%d, sign is +", i + 1);
	gsl_test_int (r.exponent, -127,
		      "float x = FLT_MIN/2^%d, exponent is -127", i + 1);
	gsl_test_str (r.mantissa, mantissa,
		      "float x = FLT_MIN/2^%d, mantissa", i + 1);
	gsl_test_int (r.type, GSL_IEEE_TYPE_DENORMAL,
		      "float x = FLT_MIN/2^%d, type is DENORMAL", i + 1);
	mantissa[i] = '0';
      }
  }
#endif

  /* Check for positive INFINITY (e.g. 2*FLT_MAX) */

  {
    float f = 3.40282347e+38f;	/* FLT_MAX */
    const char mantissa[] = "00000000000000000000000";

    gsl_ieee_float_rep r;

    float x;
    x = 2 * f;
    gsl_ieee_float_to_rep (&x, &r);

    gsl_test_int (r.sign, 0, "float x = 2*FLT_MAX, sign is +");
    gsl_test_int (r.exponent, 128, "float x = 2*FLT_MAX, exponent is 128");
    gsl_test_str (r.mantissa, mantissa, "float x = 2*FLT_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "float x = -2*FLT_MAX, type is INF");
  }

  /* Check for negative INFINITY (e.g. -2*FLT_MAX) */

  {
    float f = 3.40282347e+38f;	/* FLT_MAX */
    const char mantissa[] = "00000000000000000000000";

    gsl_ieee_float_rep r;

    float x;
    x = -2 * f;
    gsl_ieee_float_to_rep (&x, &r);

    gsl_test_int (r.sign, 1, "float x = -2*FLT_MAX, sign is -");
    gsl_test_int (r.exponent, 128, "float x = -2*FLT_MAX, exponent is 128");
    gsl_test_str (r.mantissa, mantissa, "float x = -2*FLT_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "float x = -2*FLT_MAX, type is INF");
  }

  /* Check for NAN (e.g. Inf - Inf) (float) */

  {
    gsl_ieee_float_rep r;
    float x = 1.0f, y = 2.0f, z = zerof;

    x = x / z;
    y = y / z;
    z = y - x;

    gsl_ieee_float_to_rep (&z, &r);

    /* We don't check the sign and we don't check the mantissa because
       they could be anything for a NaN */

    gsl_test_int (r.exponent, 128, "float x = NaN, exponent is 128");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NAN, "float x = NaN, type is NAN");
  }

  /* Check for +ZERO */

  {
    double d = 0.0;
    const char mantissa[]
      = "0000000000000000000000000000000000000000000000000000";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = 0, sign is +");
    gsl_test_int (r.exponent, -1023, "double x = 0, exponent is -1023");
    gsl_test_str (r.mantissa, mantissa, "double x = 0, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "double x = 0, type is ZERO");
  }

  /* Check for -ZERO */

  {
    double d =  minus_one;
    const char mantissa[]
      = "0000000000000000000000000000000000000000000000000000";
    gsl_ieee_double_rep r;

    while (d < 0) {
      d *= 0.1;
    }

    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 1, "double x = -1*0, sign is -");
    gsl_test_int (r.exponent, -1023, "double x = -1*0, exponent is -1023");
    gsl_test_str (r.mantissa, mantissa, "double x = -1*0, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "double x = -1*0, type is ZERO");
  }

  /* Check for a positive NORMAL number (e.g. 2.1) */

  {
    double d = 2.1;
    const char mantissa[]
      = "0000110011001100110011001100110011001100110011001101";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = 2.1, sign is +");
    gsl_test_int (r.exponent, 1, "double x = 2.1, exponent is 1");
    gsl_test_str (r.mantissa, mantissa, "double x = 2.1, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "double x = 2.1, type is NORMAL");
  }

  /* Check for a negative NORMAL number (e.g. -1.3304...) */

  {
    double d = -1.3303577090924210146738460025517269968986511230468750;
    const char mantissa[]
      = "0101010010010010010100101010010010001000100011101110";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 1, "double x = -1.3304..., sign is -");
    gsl_test_int (r.exponent, 0, "double x = -1.3304..., exponent is 0");
    gsl_test_str (r.mantissa, mantissa, "double x = -1.3304..., mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "double x = -1.3304..., type is NORMAL");
  }

  /* Check for a large positive NORMAL number (e.g. 3.37e297) */

  {
    double d = 3.37e297;
    const char mantissa[]
      = "0100100111001001100101111001100000100110011101000100";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = 3.37e297, sign is +");
    gsl_test_int (r.exponent, 988, "double x = 3.37e297, exponent is 998");
    gsl_test_str (r.mantissa, mantissa, "double x = 3.37e297, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "double x = 3.37e297, type is NORMAL");
  }

  /* Check for a small positive NORMAL number (e.g. 3.37e-297) */

  {
    double d = 3.37e-297;
    const char mantissa[]
    = "0001101000011011101011100001110010100001001100110111";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = 3.37e-297, sign is +");
    gsl_test_int (r.exponent, -985, "double x = 3.37e-297, exponent is -985");
    gsl_test_str (r.mantissa, mantissa, "double x = 3.37e-297, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "double x = 3.37e-297, type is NORMAL");
  }

  /* Check for DBL_MIN (smallest possible number that is not denormal) */

  {
    double d = 2.2250738585072014e-308;		/* DBL_MIN */
    const char mantissa[]
      = "0000000000000000000000000000000000000000000000000000";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = DBL_MIN, sign is +");
    gsl_test_int (r.exponent, -1022, "double x = DBL_MIN, exponent is -1022");
    gsl_test_str (r.mantissa, mantissa, "double x = DBL_MIN, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "double x = DBL_MIN, type is NORMAL");
  }

  /* Check for DBL_MAX (largest possible number that is not Inf) */

  {
    double d = 1.7976931348623157e+308;		/* DBL_MAX */
    const char mantissa[]
    = "1111111111111111111111111111111111111111111111111111";
    gsl_ieee_double_rep r;
    gsl_ieee_double_to_rep (&d, &r);

    gsl_test_int (r.sign, 0, "double x = DBL_MAX, sign is +");
    gsl_test_int (r.exponent, 1023, "double x = DBL_MAX, exponent is 1023");
    gsl_test_str (r.mantissa, mantissa, "double x = DBL_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
		  "double x = DBL_MAX, type is NORMAL");
  }

  /* Check for DENORMAL numbers (e.g. DBL_MIN/2^n) */

#ifdef TEST_DENORMAL
  {
    double d = 2.2250738585072014e-308;		/* DBL_MIN */
    char mantissa[]
      = "1000000000000000000000000000000000000000000000000000";
    int i;
    gsl_ieee_double_rep r;

    for (i = 0; i < 52; i++)
      {
	double x = d / pow (2.0, 1 + (double) i);
	mantissa[i] = '1';
	gsl_ieee_double_to_rep (&x, &r);

	gsl_test_int (r.sign, 0, "double x = DBL_MIN/2^%d, sign is +", i + 1);
	gsl_test_int (r.exponent, -1023,
		      "double x = DBL_MIN/2^%d, exponent", i + 1);
	gsl_test_str (r.mantissa, mantissa,
		      "double x = DBL_MIN/2^%d, mantissa", i + 1);
	gsl_test_int (r.type, GSL_IEEE_TYPE_DENORMAL,
		      "double x = DBL_MIN/2^%d, type is DENORMAL", i + 1);
	mantissa[i] = '0';
      }
  }
#endif

  /* Check for positive INFINITY (e.g. 2*DBL_MAX) */

  {
    double d = 1.7976931348623157e+308;		/* DBL_MAX */
    const char mantissa[]
      = "0000000000000000000000000000000000000000000000000000";
    gsl_ieee_double_rep r;

    double x;
    x = 2.0 * d;
    gsl_ieee_double_to_rep (&x, &r);

    gsl_test_int (r.sign, 0, "double x = 2*DBL_MAX, sign is +");
    gsl_test_int (r.exponent, 1024, "double x = 2*DBL_MAX, exponent is 1024");
    gsl_test_str (r.mantissa, mantissa, "double x = 2*DBL_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "double x = 2*DBL_MAX, type is INF");
  }

  /* Check for negative INFINITY (e.g. -2*DBL_MAX) */

  {
    double d = 1.7976931348623157e+308;		/* DBL_MAX */
    const char mantissa[]
      = "0000000000000000000000000000000000000000000000000000";
    gsl_ieee_double_rep r;

    double x;
    x = -2.0 * d;
    gsl_ieee_double_to_rep (&x, &r);

    gsl_test_int (r.sign, 1, "double x = -2*DBL_MAX, sign is -");
    gsl_test_int (r.exponent, 1024, "double x = -2*DBL_MAX, exponent is 1024");
    gsl_test_str (r.mantissa, mantissa, "double x = -2*DBL_MAX, mantissa");
    gsl_test_int (r.type, GSL_IEEE_TYPE_INF,"double x = -2*DBL_MAX, type is INF");
  }

  /* Check for NAN (e.g. Inf - Inf) */

  {
    gsl_ieee_double_rep r;
    double x = 1.0, y = 2.0, z = zero;

    x = x / z;
    y = y / z;
    z = y - x;

    gsl_ieee_double_to_rep (&z, &r);

    /* We don't check the sign and we don't check the mantissa because
       they could be anything for a NaN */

    gsl_test_int (r.exponent, 1024, "double x = NaN, exponent is 1024");
    gsl_test_int (r.type, GSL_IEEE_TYPE_NAN, "double x = NaN, type is NAN");
  }

  exit (gsl_test_summary ());
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpHalos.h>
#include <gsl/gsl_sf_zeta.h>

#define READ_MATCHES_GOODNESS_A 1.4421650782960167E+01
#define READ_MATCHES_GOODNESS_B 2.7010890377990435E+01
#define READ_MATCHES_GOODNESS_C ONE_THIRD
#define READ_MATCHES_GOODNESS_D -2.4490239433470165E+00

#define MAXIMUM_SCORE_COEFF 2.61      // Only if MATCH_SCORE_RANK_INDEX=-1.5.  This comes from Wolfram Alpha.  Query "sum k^(-1.5) from k = 1 to N"
#define EULER_MASCHERONI_CONST 0.5772 // Used if MATCH_SCORE_RANK_INDEX=-1

float maximum_match_score(double n_particles) {
    if(n_particles < 1)
        return (0);
    if(MATCH_SCORE_RANK_INDEX == (-1.5))
        return (MAXIMUM_SCORE_COEFF - gsl_sf_hzeta(-MATCH_SCORE_RANK_INDEX, n_particles));
    else if(MATCH_SCORE_RANK_INDEX == (-1.))
        return (EULER_MASCHERONI_CONST + take_ln(n_particles));
    else if(MATCH_SCORE_RANK_INDEX == (-TWO_THIRDS))
        return ((float)pow(READ_MATCHES_GOODNESS_A + READ_MATCHES_GOODNESS_B * n_particles, READ_MATCHES_GOODNESS_C) + READ_MATCHES_GOODNESS_D);
    else
        SID_exit_error(
                "Invalid MATCH_SCORE_RANK_INDEX passed to maximum_match_score().  Please update the code to accomodate MATCH_SCORE_RANK_INDEX=%lf.",
                SID_ERROR_LOGIC, MATCH_SCORE_RANK_INDEX);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <string.h>

//** Init Multinomial
gsl_rng* init_gsl_rng(){
		gsl_rng* r;
 		const gsl_rng_type* T2;

    srand(time(NULL));
    unsigned int seed;
    seed= rand();

    gsl_rng_env_setup();

    T2 = gsl_rng_default;
    r = gsl_rng_alloc (T2);
    gsl_rng_set (r, seed);

		return r;
}

unsigned int substract(unsigned int val){
	if(val>0)
		return val-1;
	else
		return 0;
}

//** Data must be dictionated
//** Doc-Word-topic
void dt2b_inference(int** data , unsigned int** uKu, unsigned  int** pKp, unsigned  int** KuKp, const int it, const int Ku,const int Kp, const int nusers, const int nplaces, const int nrecs){

	float l,prev_l,prob=0;

 	//** Allocate Dynamic Tables
	unsigned int zu[Ku];
	memset( zu, 0, Ku*sizeof(unsigned int) );
	unsigned int zp[Kp];
	memset( zp, 0, Kp*sizeof(unsigned int) );
	double p_zuzp[Ku][Kp];

  //** Initialization
	gsl_rng * r;
	r=init_gsl_rng();

	//** Initialize Dynamic Tables
	unsigned int oneutopic,oneptopic;
	size_t rowid;
	for(rowid=0;rowid<nrecs;rowid++){
			oneutopic=data[rowid][2];
			oneptopic=data[rowid][3];
			zu[(size_t)oneutopic]++;
			zp[(size_t)oneptopic]++;
	}

	//** Iterate
	size_t i,s,kp,ku;
	int curr_u,curr_p,curr_ku,curr_kp;
	for(i=0;i<it;i++){
		l=0;
		//printf("iteration %zu\n",i);
		for(rowid=0;rowid<nrecs;rowid++){

			//** Discount
			curr_u=data[rowid][0];
			curr_p=data[rowid][1];
			curr_ku=data[rowid][2];
			curr_kp=data[rowid][3];

			uKu[curr_u][curr_ku]=substract(uKu[curr_u][curr_ku]);
			pKp[curr_p][curr_kp]=substract(pKp[curr_p][curr_kp]);
			KuKp[curr_ku][curr_kp]=substract(KuKp[curr_ku][curr_kp]);
			zu[curr_ku]=substract(zu[curr_ku]);
			zp[curr_kp]=substract(zp[curr_kp]);

			//** Obtain Posterior
			prob=0;
			for(ku=0;ku<Ku;ku++){
				for(kp=0;kp<Kp;kp++){
					p_zuzp[ku][kp]=((KuKp[ku][kp]+0.5)*(uKu[curr_u][ku]+0.5)*(pKp[curr_p][kp]+0.5))/((zu[ku]+0.5*nusers)*(zp[kp]+0.5*nplaces));
					prob+=p_zuzp[ku][kp];
				}
			}
			l+=log(prob);

			//** Obtain one sample
    	unsigned int samples[Ku*Kp];
    	double* p_z = (double*) p_zuzp;

   		gsl_ran_multinomial(r, Ku*Kp, 1, p_z, samples);

 			size_t chosen_ku,chosen_kp;
 			for(s=0;s<Ku*Kp;s++){

 				if (samples[s]>0){

 					//get two indexes from one
 					chosen_ku=(unsigned int) (s/Kp);
 					chosen_kp=(unsigned int) (s%Kp);

 					break;
 				}
 			}

 			//** Update
 			data[rowid][2]=chosen_ku;
 			data[rowid][3]=chosen_kp;
 			uKu[curr_u][chosen_ku]+=1;
 			pKp[curr_p][chosen_kp]+=1;
 			KuKp[chosen_ku][chosen_kp]+=1;
 			zu[chosen_ku]+=1;
 			zp[chosen_kp]+=1;

		}

	if (i==0)
		prev_l=l;

	//** Check convergence
		if (l<prev_l && i>50){
			printf("Stopped at iteration %zu",i);
			break;
		}

		printf("Likelihood at iteration %zu: %f \n",i,l);

	}

	gsl_rng_free (r);
}

max_stars_repo_name: danrugeles/YLDA

max_issues_repo_name: danrugeles/YLDA

max_forks_repo_name: danrugeles/YLDA


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <nubspline.h>
#include "SimulationSteps.h"
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

void eval_NUBspline_1d_d(NUBspline_1d_d * restrict spline, double x,
	double* restrict val);

void eval_NUBspline_1d_d_vg(NUBspline_1d_d * restrict spline, double x,
	double* restrict val, double* restrict grad);

void eval_NUBspline_1d_d_vgl(NUBspline_1d_d * restrict spline, double x,
	double* restrict val, double* restrict grad, double* restrict lapl);

void InterpolateWithSplines(int NumPointsToInterpolate, int NumPointsToEvaluate, double* x, double* NodalX, double *func, double *NodalFunc, double* NodalDeriv)
{

	NUgrid* myGrid = create_general_grid(x, NumPointsToInterpolate);
	BCtype_d myBC;
	myBC.lCode = NATURAL;
	myBC.rCode = NATURAL;
	NUBspline_1d_d *mySpline = create_NUBspline_1d_d(myGrid, myBC, func);
	for(int i =0; i < NumPointsToEvaluate; i++)
	{
		if (NodalDeriv != NULL)
			eval_NUBspline_1d_d_vg(mySpline, NodalX[i], &NodalFunc[i], &NodalDeriv[i]);
		else
			eval_NUBspline_1d_d(mySpline, NodalX[i], &NodalFunc[i]);
	}

}

void InterpolateWithGSLSplines(int NumPointsToInterpolate, int NumPointsToEvaluate, double* x, double* NodalX, double* func, double* NodalFunc, double* NodalDeriv)
{
	gsl_interp_accel *acc = gsl_interp_accel_alloc();
	gsl_spline* spline = gsl_spline_alloc(gsl_interp_cspline, NumPointsToInterpolate);
	gsl_spline_init(spline, x, func, NumPointsToInterpolate);

	for(int j = 0; j < NumPointsToEvaluate; j++)
	{
		NodalFunc[j] = gsl_spline_eval(spline, NodalX[j], acc);
		if (NodalDeriv != NULL)
			NodalDeriv[j] = gsl_spline_eval_deriv(spline, NodalX[j], acc);
	}

}

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


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

#ifndef __GSL_RNG_H__
#define __GSL_RNG_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif
#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>

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

__BEGIN_DECLS

typedef struct
  {
    const char *name;
    unsigned long int max;
    unsigned long int min;
    size_t size;
    void (*set) (void *state, unsigned long int seed);
    unsigned long int (*get) (void *state);
    double (*get_double) (void *state);
  }
gsl_rng_type;

typedef struct
  {
    const gsl_rng_type * type;
    void *state;
  }
gsl_rng;

/* These structs also need to appear in default.c so you can select
   them via the environment variable GSL_RNG_TYPE */

GSL_VAR const gsl_rng_type *gsl_rng_borosh13;
GSL_VAR const gsl_rng_type *gsl_rng_coveyou;
GSL_VAR const gsl_rng_type *gsl_rng_cmrg;
GSL_VAR const gsl_rng_type *gsl_rng_fishman18;
GSL_VAR const gsl_rng_type *gsl_rng_fishman20;
GSL_VAR const gsl_rng_type *gsl_rng_fishman2x;
GSL_VAR const gsl_rng_type *gsl_rng_gfsr4;
GSL_VAR const gsl_rng_type *gsl_rng_knuthran;
GSL_VAR const gsl_rng_type *gsl_rng_knuthran2;
GSL_VAR const gsl_rng_type *gsl_rng_knuthran2002;
GSL_VAR const gsl_rng_type *gsl_rng_lecuyer21;
GSL_VAR const gsl_rng_type *gsl_rng_minstd;
GSL_VAR const gsl_rng_type *gsl_rng_mrg;
GSL_VAR const gsl_rng_type *gsl_rng_mt19937;
GSL_VAR const gsl_rng_type *gsl_rng_mt19937_1999;
GSL_VAR const gsl_rng_type *gsl_rng_mt19937_1998;
GSL_VAR const gsl_rng_type *gsl_rng_r250;
GSL_VAR const gsl_rng_type *gsl_rng_ran0;
GSL_VAR const gsl_rng_type *gsl_rng_ran1;
GSL_VAR const gsl_rng_type *gsl_rng_ran2;
GSL_VAR const gsl_rng_type *gsl_rng_ran3;
GSL_VAR const gsl_rng_type *gsl_rng_rand;
GSL_VAR const gsl_rng_type *gsl_rng_rand48;
GSL_VAR const gsl_rng_type *gsl_rng_random128_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random128_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random128_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_random256_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random256_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random256_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_random32_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random32_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random32_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_random64_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random64_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random64_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_random8_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random8_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random8_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_random_bsd;
GSL_VAR const gsl_rng_type *gsl_rng_random_glibc2;
GSL_VAR const gsl_rng_type *gsl_rng_random_libc5;
GSL_VAR const gsl_rng_type *gsl_rng_randu;
GSL_VAR const gsl_rng_type *gsl_rng_ranf;
GSL_VAR const gsl_rng_type *gsl_rng_ranlux;
GSL_VAR const gsl_rng_type *gsl_rng_ranlux389;
GSL_VAR const gsl_rng_type *gsl_rng_ranlxd1;
GSL_VAR const gsl_rng_type *gsl_rng_ranlxd2;
GSL_VAR const gsl_rng_type *gsl_rng_ranlxs0;
GSL_VAR const gsl_rng_type *gsl_rng_ranlxs1;
GSL_VAR const gsl_rng_type *gsl_rng_ranlxs2;
GSL_VAR const gsl_rng_type *gsl_rng_ranmar;
GSL_VAR const gsl_rng_type *gsl_rng_slatec;
GSL_VAR const gsl_rng_type *gsl_rng_taus;
GSL_VAR const gsl_rng_type *gsl_rng_taus2;
GSL_VAR const gsl_rng_type *gsl_rng_taus113;
GSL_VAR const gsl_rng_type *gsl_rng_transputer;
GSL_VAR const gsl_rng_type *gsl_rng_tt800;
GSL_VAR const gsl_rng_type *gsl_rng_uni;
GSL_VAR const gsl_rng_type *gsl_rng_uni32;
GSL_VAR const gsl_rng_type *gsl_rng_vax;
GSL_VAR const gsl_rng_type *gsl_rng_waterman14;
GSL_VAR const gsl_rng_type *gsl_rng_zuf;

GSL_FUN const gsl_rng_type ** gsl_rng_types_setup(void);

GSL_VAR const gsl_rng_type *gsl_rng_default;
GSL_VAR unsigned long int gsl_rng_default_seed;

GSL_FUN gsl_rng *gsl_rng_alloc (const gsl_rng_type * T);
GSL_FUN int gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src);
GSL_FUN gsl_rng *gsl_rng_clone (const gsl_rng * r);

GSL_FUN void gsl_rng_free (gsl_rng * r);

GSL_FUN void gsl_rng_set (const gsl_rng * r, unsigned long int seed);
GSL_FUN unsigned long int gsl_rng_max (const gsl_rng * r);
GSL_FUN unsigned long int gsl_rng_min (const gsl_rng * r);
GSL_FUN const char *gsl_rng_name (const gsl_rng * r);

GSL_FUN int gsl_rng_fread (FILE * stream, gsl_rng * r);
GSL_FUN int gsl_rng_fwrite (FILE * stream, const gsl_rng * r);

GSL_FUN size_t gsl_rng_size (const gsl_rng * r);
GSL_FUN void * gsl_rng_state (const gsl_rng * r);

GSL_FUN void gsl_rng_print_state (const gsl_rng * r);

GSL_FUN const gsl_rng_type * gsl_rng_env_setup (void);

GSL_FUN INLINE_DECL unsigned long int gsl_rng_get (const gsl_rng * r);
GSL_FUN INLINE_DECL double gsl_rng_uniform (const gsl_rng * r);
GSL_FUN INLINE_DECL double gsl_rng_uniform_pos (const gsl_rng * r);
GSL_FUN INLINE_DECL unsigned long int gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n);

#ifdef HAVE_INLINE

INLINE_FUN unsigned long int
gsl_rng_get (const gsl_rng * r)
{
  return (r->type->get) (r->state);
}

INLINE_FUN double
gsl_rng_uniform (const gsl_rng * r)
{
  return (r->type->get_double) (r->state);
}

INLINE_FUN double
gsl_rng_uniform_pos (const gsl_rng * r)
{
  double x ;
  do
    {
      x = (r->type->get_double) (r->state) ;
    }
  while (x == 0) ;

  return x ;
}

/* Note: to avoid integer overflow in (range+1) we work with scale =
   range/n = (max-min)/n rather than scale=(max-min+1)/n, this reduces
   efficiency slightly but avoids having to check for the out of range
   value.  Note that range is typically O(2^32) so the addition of 1
   is negligible in most usage. */

INLINE_FUN unsigned long int
gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
{
  unsigned long int offset = r->type->min;
  unsigned long int range = r->type->max - offset;
  unsigned long int scale;
  unsigned long int k;

  if (n > range || n == 0)
    {
      GSL_ERROR_VAL ("invalid n, either 0 or exceeds maximum value of generator",
                     GSL_EINVAL, 0) ;
    }

  scale = range / n;

  do
    {
      k = (((r->type->get) (r->state)) - offset) / scale;
    }
  while (k >= n);

  return k;
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_RNG_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* very simple approximation */
double st_gamma(double x)
{
  return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
}

#define A 12
double sp_gamma(double z)
{
  const int a = A;
  static double c_space[A];
  static double *c = NULL;
  int k;
  double accm;

  if ( c == NULL ) {
    double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/
    c = c_space;
    c[0] = sqrt(2.0*M_PI);
    for(k=1; k < a; k++) {
      c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
	  k1_factrl *= -k;
    }
  }
  accm = c[0];
  for(k=1; k < a; k++) {
    accm += c[k] / ( z + k );
  }
  accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */
  return accm/z;
}

int main()
{
  double x;

  printf("%15s%15s%15s%15s\n", "Stirling", "Spouge", "GSL", "libm");
  for(x=1.0; x <= 10.0; x+=1.0) {
    printf("%15.8lf%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0),
	   gsl_sf_gamma(x/3.0), tgamma(x/3.0));
  }
  return 0;
}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
/* linalg/invtri_complex.c
 *
 * Copyright (C) 2019 Patrick Alken
 *
 * This 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, or (at your option) any
 * later version.
 *
 * This source 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.
 *
 * This module contains code to invert complex triangular matrices
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

#include "recurse.h"

static int complex_tri_invert_L2(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T);
static int complex_tri_invert_L3(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T);
static int triangular_singular(const gsl_matrix_complex * T);

/*
gsl_linalg_complex_tri_invert()
  Invert a complex triangular matrix T using Level 3 BLAS

Inputs: Uplo - CblasUpper or CblasLower
        Diag - unit triangular?
        T    - on output the upper (or lower) part of T
               is replaced by its inverse

Return: success/error
*/

int
gsl_linalg_complex_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      int status;

      status = triangular_singular(T);
      if (status)
        return status;

      return complex_tri_invert_L3(Uplo, Diag, T);
    }
}

/*
complex_tri_invert_L2()
  Invert a complex triangular matrix T

Inputs: Uplo - CblasUpper or CblasLower
        Diag - unit triangular?
        T    - on output the upper (or lower) part of T
               is replaced by its inverse

Return: success/error

Notes:
1) Based on LAPACK routine ZTRTI2 using Level 2 BLAS
*/

static int
complex_tri_invert_L2(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t i;

      if (Uplo == CblasUpper)
        {
          for (i = 0; i < N; ++i)
            {
              gsl_complex * Tii = gsl_matrix_complex_ptr(T, i, i);
              gsl_complex aii;

              if (Diag == CblasNonUnit)
                {
                  *Tii = gsl_complex_inverse(*Tii);
                  GSL_REAL(aii) = -GSL_REAL(*Tii);
                  GSL_IMAG(aii) = -GSL_IMAG(*Tii);
                }
              else
                aii = GSL_COMPLEX_NEGONE;

              if (i > 0)
                {
                  gsl_matrix_complex_view m = gsl_matrix_complex_submatrix(T, 0, 0, i, i);
                  gsl_vector_complex_view v = gsl_matrix_complex_subcolumn(T, i, 0, i);

                  gsl_blas_ztrmv(CblasUpper, CblasNoTrans, Diag, &m.matrix, &v.vector);
                  gsl_blas_zscal(aii, &v.vector);
                }
            }
        }
      else
        {
          for (i = 0; i < N; ++i)
            {
              size_t j = N - i - 1;
              gsl_complex * Tjj = gsl_matrix_complex_ptr(T, j, j);
              gsl_complex ajj;

              if (Diag == CblasNonUnit)
                {
                  *Tjj = gsl_complex_inverse(*Tjj);
                  GSL_REAL(ajj) = -GSL_REAL(*Tjj);
                  GSL_IMAG(ajj) = -GSL_IMAG(*Tjj);
                }
              else
                ajj = GSL_COMPLEX_NEGONE;

              if (j < N - 1)
                {
                  gsl_matrix_complex_view m = gsl_matrix_complex_submatrix(T, j + 1, j + 1, N - j - 1, N - j - 1);
                  gsl_vector_complex_view v = gsl_matrix_complex_subcolumn(T, j, j + 1, N - j - 1);

                  gsl_blas_ztrmv(CblasLower, CblasNoTrans, Diag, &m.matrix, &v.vector);
                  gsl_blas_zscal(ajj, &v.vector);
                }
            }
        }

      return GSL_SUCCESS;
    }
}

/*
complex_tri_invert_L3()
  Invert a complex triangular matrix T

Inputs: Uplo - CblasUpper or CblasLower
        Diag - unit triangular?
        T    - on output the upper (or lower) part of T
               is replaced by its inverse

Return: success/error

Notes:
1) Based on ReLAPACK using Level 3 BLAS
*/

static int
complex_tri_invert_L3(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_INVTRI)
    {
      /* use Level 2 algorithm */
      return complex_tri_invert_L2(Uplo, Diag, T);
    }
  else
    {
      /*
       * partition matrix:
       *
       * T11 T12
       * T21 T22
       *
       * where T11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT_COMPLEX(N);
      const size_t N2 = N - N1;
      gsl_matrix_complex_view T11 = gsl_matrix_complex_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_complex_view T12 = gsl_matrix_complex_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_complex_view T21 = gsl_matrix_complex_submatrix(T, N1, 0, N2, N1);
      gsl_matrix_complex_view T22 = gsl_matrix_complex_submatrix(T, N1, N1, N2, N2);

      /* recursion on T11 */
      status = complex_tri_invert_L3(Uplo, Diag, &T11.matrix);
      if (status)
        return status;

      if (Uplo == CblasLower)
        {
          /* T21 = - T21 * T11 */
          gsl_blas_ztrmm(CblasRight, Uplo, CblasNoTrans, Diag, GSL_COMPLEX_NEGONE, &T11.matrix, &T21.matrix);

          /* T21 = T22 * T21^{-1} */
          gsl_blas_ztrsm(CblasLeft, Uplo, CblasNoTrans, Diag, GSL_COMPLEX_ONE, &T22.matrix, &T21.matrix);
        }
      else
        {
          /* T12 = - T11 * T12 */
          gsl_blas_ztrmm(CblasLeft, Uplo, CblasNoTrans, Diag, GSL_COMPLEX_NEGONE, &T11.matrix, &T12.matrix);

          /* T12 = T12 * T22^{-1} */
          gsl_blas_ztrsm(CblasRight, Uplo, CblasNoTrans, Diag, GSL_COMPLEX_ONE, &T22.matrix, &T12.matrix);
        }

      /* recursion on T22 */
      status = complex_tri_invert_L3(Uplo, Diag, &T22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

static int
triangular_singular(const gsl_matrix_complex * T)
{
  size_t i;

  for (i = 0; i < T->size1; ++i)
    {
      gsl_complex z = gsl_matrix_complex_get(T, i, i);
      if (GSL_REAL(z) == 0.0 && GSL_IMAG(z) == 0.0)
        return GSL_ESING;
    }

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_pow_int.h"
#include "gsl_sf_elljac.h"

/* See [Thompson, Atlas for Computing Mathematical Functions] */

int
gsl_sf_elljac_e(double u, double m, double * sn, double * cn, double * dn)
{
  if(fabs(m) > 1.0) {
    *sn = 0.0;
    *cn = 0.0;
    *dn = 0.0;
    GSL_ERROR ("|m| > 1.0", GSL_EDOM);
  }
  else if(fabs(m) < 2.0*GSL_DBL_EPSILON) {
    *sn = sin(u);
    *cn = cos(u);
    *dn = 1.0;
    return GSL_SUCCESS;
  }
  else if(fabs(m - 1.0) < 2.0*GSL_DBL_EPSILON) {
    *sn = tanh(u);
    *cn = 1.0/cosh(u);
    *dn = *cn;
    return GSL_SUCCESS;
  }
  else {
    int status = GSL_SUCCESS;
    const int N = 16;
    double   a[16];
    double   b[16];
    double   c[16];
    double phi[16];
    double psi[16]; /* psi[i] := phi[i] - Pi 2^{i-1} */
    double two_N;
    int n = 0;

    a[0] = 1.0;
    b[0] = sqrt(1.0 - m);
    c[0] = sqrt(m);

    while( fabs(c[n]) > 4.0 * GSL_DBL_EPSILON) {
      a[n+1] = 0.5 * (a[n] + b[n]);
      b[n+1] = sqrt(a[n] * b[n]);
      c[n+1] = 0.5 * (a[n] - b[n]);
      if(n >= N - 2) {
        status = GSL_EMAXITER;
	c[N-1] = 0.0;
	break;
      }
      ++n;
    }

    --n;
    two_N = (double)(1 << n ); /* 2^n */  /* gsl_sf_pow_int(2.0, n); */
    phi[n] = two_N * a[n] * u;
    psi[n] = two_N * (a[n]*u - 0.5*M_PI);

    while(n > 0) {
      const double psi_sgn = ( n == 1 ? -1.0 : 1.0 );
      const double phi_asin_arg = c[n] * sin(phi[n])/a[n];
      const double psi_asin_arg = c[n]/a[n] * psi_sgn * sin(psi[n]);
      const double phi_asin = asin(phi_asin_arg);
      const double psi_asin = asin(psi_asin_arg);
      phi[n-1] = 0.5 * (phi[n] + phi_asin);
      psi[n-1] = 0.5 * (psi[n] + psi_asin);
      --n;
    }

    *sn = sin(phi[0]);
    *cn = cos(phi[0]);
    {
      /* const double dn_method_1 = *cn / cos(phi[1] - phi[0]); */
      const double dn_method_2 = sin(psi[0])/sin(psi[1] - psi[0]);
      *dn = dn_method_2;
      /* printf("%18.16g  %18.16g\n", dn_method_1, dn_method_2); */
    }

    return status;
  }
}


text:
#pragma once

#include <gsl/gsl>
#include <iostream>

#include "../utils/Concepts.h"
#include "General.h"
#include "Random.h"

namespace math {
	using gsl::narrow_cast;
	using utils::concepts::SignedNumeric, utils::concepts::Integral,
		utils::concepts::SignedIntegral, utils::concepts::FloatingPoint;

	enum class Vec2Idx : size_t
	{
		X = 0ULL,
		Y
	};

	template<SignedNumeric T = float>
	class Vec2 {
	  public:
		/// @brief Creates a default `Vec2`
		constexpr Vec2() noexcept = default;

		/// @brief Creates a new `Vec2` with the given x and y components
		///
		/// @param x - The x component
		/// @param y - The y component
		constexpr Vec2(T x, T y) noexcept : elements{x, y} {
		}
		constexpr Vec2(const Vec2& vec) noexcept = default;
		constexpr Vec2(Vec2&& vec) noexcept = default;
		constexpr ~Vec2() noexcept = default;

		/// @brief Returns the x component
		///
		/// @return a const ref to the x component
		[[nodiscard]] inline constexpr auto x() const noexcept -> const T& {
			return elements[X];
		}

		/// @brief Returns the x component
		///
		/// @return a mutable (ie: non-const) ref to the x component
		inline constexpr auto x() noexcept -> T& {
			return elements[X];
		}

		/// @brief Returns the y component
		///
		/// @return a const ref to the y component
		[[nodiscard]] inline constexpr auto y() const noexcept -> const T& {
			return elements[Y];
		}

		/// @brief Returns the y component
		///
		/// @return a mutable (ie: non-const) ref to the y component
		inline constexpr auto y() noexcept -> T& {
			return elements[Y];
		}

		/// @brief Returns the magnitude (length) of the vector
		///
		/// @return The magnitude
		template<FloatingPoint TT = float>
		[[nodiscard]] inline constexpr auto magnitude() const noexcept -> TT {
			return General::sqrt(narrow_cast<TT>(magnitude_squared()));
		}

		/// @brief Returns the dot product of this and `vec`
		///
		/// @param vec - The vector to perform the dot product with
		///
		/// @return The dot product
		template<FloatingPoint TT = float>
		[[nodiscard]] inline constexpr auto dot_prod(const Vec2<TT>& vec) const noexcept -> TT {
			return narrow_cast<TT>(x()) * vec.x() + narrow_cast<TT>(y()) * vec.y();
		}

		/// @brief Returns the dot product of this and `vec`
		///
		/// @param vec - The vector to perform the dot product with
		///
		/// @return The dot product
		template<SignedIntegral TT = int>
		[[nodiscard]] inline constexpr auto dot_prod(const Vec2<TT>& vec) const noexcept -> T {
			return x() * narrow_cast<T>(vec.x()) + y() * narrow_cast<T>(vec.y());
		}

		/// @brief Performs the 2D equivalent to the cross product
		/// @note This is equivalent to the z component of a 3D cross product
		///
		/// @param vec The vector to perform the cross product with
		///
		/// @return The 2D cross product
		template<FloatingPoint TT = float>
		[[nodiscard]] inline constexpr auto cross_prod(const Vec2<TT>& vec) const noexcept -> TT {
			return narrow_cast<TT>(x()) * vec.y() - narrow_cast<TT>(y()) * vec.x();
		}

		/// @brief Performs the 2D equivalent to the cross product
		/// @note This is equivalent to the z component of a 3D cross product
		///
		/// @param vec The vector to perform the cross product with
		///
		/// @return The 2D cross product
		template<SignedIntegral TT = int>
		[[nodiscard]] inline constexpr auto cross_prod(const Vec2<TT>& vec) const noexcept -> T {
			return x() * narrow_cast<T>(vec.y()) - y() * narrow_cast<T>(vec.x());
		}

		/// @brief Returns **a** vector normal to this one
		/// @note this is not necessarily the **only** vector normal to this one
		///
		/// @return a vector normal to this
		template<FloatingPoint TT = float>
		[[nodiscard]] inline constexpr auto normal() const noexcept -> Vec2<TT> {
			constexpr auto _x = narrow_cast<TT>(1.0);
			const auto _y = (-_x * narrow_cast<TT>(x())) / narrow_cast<TT>(y());
			return {_x, _y};
		}

		/// @brief Returns this vector with normalized magnitude
		///
		/// @return this vector, normalized
		template<FloatingPoint TT = float>
		[[nodiscard]] inline constexpr auto normalized() const noexcept -> Vec2<TT> {
			return std::move(*this / magnitude<TT>());
		}

		template<SignedNumeric TT = float>
		[[nodiscard]] inline static constexpr auto random() noexcept -> Vec2<TT> {
			return {random_value<TT>(), random_value<TT>()};
		}

		template<SignedNumeric TT = float>
		[[nodiscard]] inline static constexpr auto random(TT min, TT max) noexcept -> Vec2<TT> {
			return {random_value<TT>(min, max), random_value<TT>(min, max)};
		}

		template<FloatingPoint TT = float>
		[[nodiscard]] inline static constexpr auto random_in_unit_circle() noexcept -> Vec2<TT> {
			do {
				auto val = random(narrow_cast<TT>(-1), narrow_cast<TT>(1));
				if(val.magnitude_squared() < narrow_cast<TT>(1)) {
					return val;
				}
			} while(true);
		}

		constexpr auto operator=(const Vec2& vec) noexcept -> Vec2& = default;
		constexpr auto operator=(Vec2&& vec) noexcept -> Vec2& = default;

		template<FloatingPoint TT = float>
		inline constexpr auto operator==(const Vec2<TT>& vec) const noexcept -> bool {
			const auto xEqual
				= General::abs<TT>(narrow_cast<TT>(x()) - vec.x()) < narrow_cast<TT>(0.01);
			const auto yEqual
				= General::abs<TT>(narrow_cast<TT>(y()) - vec.y()) < narrow_cast<TT>(0.01);
			return xEqual && yEqual;
		}

		template<SignedIntegral TT = int>
		inline constexpr auto operator==(const Vec2<TT>& vec) const noexcept -> bool {
			if constexpr(FloatingPoint<T>) {
				const auto xEqual
					= General::abs<T>(x() - narrow_cast<T>(vec.x())) < narrow_cast<TT>(0.01);
				const auto yEqual
					= General::abs<T>(y() - narrow_cast<T>(vec.y())) < narrow_cast<TT>(0.01);
				return xEqual && yEqual;
			}
			else {
				return x() == narrow_cast<T>(vec.x()) && y() == narrow_cast<T>(vec.y());
			}
		}

		template<SignedNumeric TT = T>
		inline constexpr auto operator!=(const Vec2<TT>& vec) const noexcept -> bool {
			return !(*this == vec);
		}

		inline constexpr auto operator-() const noexcept -> Vec2 {
			return {-x(), -y()};
		}

		inline constexpr auto operator[](Vec2Idx i) const noexcept -> T {
			const auto index = static_cast<size_t>(i);
			return elements[index]; // NOLINT
		}
		inline constexpr auto operator[](Vec2Idx i) noexcept -> T& {
			const auto index = static_cast<size_t>(i);
			return elements[index]; // NOLINT
		}

		template<FloatingPoint TT = float>
		inline constexpr auto operator+(const Vec2<TT>& vec) const noexcept -> Vec2<TT> {
			return {narrow_cast<TT>(x()) + vec.x(), narrow_cast<TT>(y()) + vec.y()};
		}

		template<SignedIntegral TT = int>
		inline constexpr auto operator+(const Vec2<TT>& vec) const noexcept -> Vec2<T> {
			return {x() + narrow_cast<T>(vec.x()), y() + narrow_cast<T>(vec.y())};
		}

		template<SignedNumeric TT = T>
		inline constexpr auto operator+=(const Vec2<TT>& vec) noexcept -> Vec2& {
			x() += narrow_cast<T>(vec.x());
			y() += narrow_cast<T>(vec.y());
			return *this;
		}

		template<FloatingPoint TT = float>
		inline constexpr auto operator-(const Vec2<TT>& vec) const noexcept -> Vec2<TT> {
			return {narrow_cast<TT>(x()) - vec.x(), narrow_cast<TT>(y()) - vec.y()};
		}

		template<SignedIntegral TT = int>
		inline constexpr auto operator-(const Vec2<TT>& vec) const noexcept -> Vec2<T> {
			return {x() - narrow_cast<T>(vec.x()), y() - narrow_cast<T>(vec.y())};
		}

		template<SignedNumeric TT = T>
		inline constexpr auto operator-=(const Vec2<TT>& vec) noexcept -> Vec2& {
			x() -= narrow_cast<T>(vec.x());
			y() -= narrow_cast<T>(vec.y());
			return *this;
		}

		inline constexpr auto operator*(FloatingPoint auto s) const noexcept -> Vec2<decltype(s)> {
			using TT = decltype(s);
			return {narrow_cast<TT>(x()) * s, narrow_cast<TT>(y()) * s};
		}

		inline constexpr auto operator*(SignedIntegral auto s) const noexcept -> Vec2 {
			auto scalar = narrow_cast<T>(s);
			return {x() * scalar, y() * scalar};
		}

		inline constexpr auto operator*=(FloatingPoint auto s) noexcept -> Vec2& {
			using TT = decltype(s);
			x() = narrow_cast<T>(narrow_cast<TT>(x()) * s);
			y() = narrow_cast<T>(narrow_cast<TT>(y()) * s);
			return *this;
		}

		inline constexpr auto operator*=(SignedIntegral auto s) noexcept -> Vec2& {
			auto scalar = narrow_cast<T>(s);
			x() *= scalar;
			y() *= scalar;
			return *this;
		}

		friend inline constexpr auto
		operator*(FloatingPoint auto lhs, const Vec2& rhs) noexcept -> Vec2<decltype(lhs)> {
			return rhs * lhs;
		}

		friend inline constexpr auto
		operator*(SignedIntegral auto lhs, const Vec2& rhs) noexcept -> Vec2 {
			return rhs * lhs;
		}

		inline constexpr auto operator/(FloatingPoint auto s) const noexcept -> Vec2<decltype(s)> {
			using TT = decltype(s);
			return {narrow_cast<TT>(x()) / s, narrow_cast<TT>(y()) / s};
		}

		inline constexpr auto operator/(SignedIntegral auto s) const noexcept -> Vec2 {
			auto scalar = narrow_cast<T>(s);
			return {x() / scalar, y() / scalar};
		}

		friend inline constexpr auto
		operator/(FloatingPoint auto lhs, const Vec2& rhs) noexcept -> Vec2<decltype(lhs)> {
			return rhs / lhs;
		}

		friend inline constexpr auto
		operator/(SignedIntegral auto lhs, const Vec2& rhs) noexcept -> Vec2 {
			return rhs / lhs;
		}

		inline constexpr auto operator/=(FloatingPoint auto s) noexcept -> Vec2 {
			using TT = decltype(s);
			x() = narrow_cast<T>(narrow_cast<TT>(x()) / s);
			y() = narrow_cast<T>(narrow_cast<TT>(y()) / s);
			return *this;
		}

		inline constexpr auto operator/=(SignedIntegral auto s) noexcept -> Vec2 {
			auto scalar = narrow_cast<T>(s);
			x() /= scalar;
			y() /= scalar;
			return *this;
		}

		friend inline constexpr auto
		operator<<(std::ostream& out, const Vec2& vec) noexcept -> std::ostream& {
			return out << vec.x() << ' ' << vec.y();
		}

	  private:
		static constexpr size_t NUM_ELEMENTS = static_cast<size_t>(Vec2Idx::Y) + 1;
		T elements[NUM_ELEMENTS] = {narrow_cast<T>(0), narrow_cast<T>(0)}; // NOLINT

		/// @brief Calculates the magnitude squared of this vector
		///
		/// @return The magnitude squared
		[[nodiscard]] inline constexpr auto magnitude_squared() const noexcept -> T {
			return x() * x() + y() * y();
		}

		/// Index for x component
		static constexpr size_t X = static_cast<size_t>(Vec2Idx::X);
		/// Index for y component
		static constexpr size_t Y = static_cast<size_t>(Vec2Idx::Y);
	};

	template<FloatingPoint T>
	explicit Vec2(T, T) -> Vec2<T>;

	template<SignedIntegral T>
	explicit Vec2(T, T) -> Vec2<T>;
} // namespace math


text:
#ifndef TSNE_H
#define TSNE_H

/* Cmake will define ${LIBRARY_NAME}_EXPORTS on Windows when it
configures to build a shared library. If you are going to use
another build system on windows or create the visual studio
projects by hand you need to define ${LIBRARY_NAME}_EXPORTS when
building a DLL on windows.
*/
// We are using the Visual Studio Compiler and building Shared libraries

#if defined (_WIN32)
#if defined(examplelib_EXPORTS)
#define DLLEXPORT __declspec(dllexport)
#else
#define DLLEXPORT __declspec(dllimport)
#endif
#else
#define DLLEXPORT
#endif
#ifndef SWIG
#include <math.h>
#include <iostream>    // std::cout, std::fixed
#include <iomanip>      // std::setprecision
#include <stdlib.h>
#include <cstring>
#include <time.h>
#include <stdlib.h>
#include <algorithm>
#include <vector>
#include <stdio.h>
#include <queue>
#include <limits>
#include <iostream>
#include <cfloat>
#include <limits>
#include <armadillo>
using namespace std;
using namespace arma;
#include "quadtree.h"
#include "vptree.h"

extern "C" {
#include <cblas.h>
}
#else
#endif
/*
 *  tsne.h
 *  Header file for t-SNE.
 *
 *  Created by Laurens van der Maaten.
 *  Copyright 2012, Delft University of Technology. All rights reserved.
 *
 */

static inline double sign(double x)
{
    return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0));
}

class DLLEXPORT TSNE {
public:
    void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta);
    bool load_data(const char * filename,double** data, int* n, int* d, double* theta, double* perplexity);
    void save_data(const char* filename,double* data, int* landmarks, double* costs, int n, int d);

    void symmetrizeMatrix(int** row_P, int** col_P, double** val_P, int N); // should be static?!

private:
    void computeGradient(double* P, int* inp_row_P, int* inp_col_P, double* inp_val_P, double* Y, int N, int D, double* dC, double theta);
    void computeExactGradient(double* P, double* Y, int N, int D, double* dC);
    double evaluateError(double* P, double* Y, int N);
    double evaluateError(int* row_P, int* col_P, double* val_P, double* Y, int N, double theta);
    void zeroMean(double* X, int N, int D);
    void computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity);
    void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, int K);
    void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, double threshold);
    void computeSquaredEuclideanDistance(double* X, int N, int D, double* DD);
    double randn();
};
/* the following Math  class is used to test the python wrapper work for basic pointer manipulation*/
class DLLEXPORT Math4py {
public:
    Math4py();
    int pi() const;
    void pi(int pi);
private:
    int _pi;
};
int DLLEXPORT fact(int n);

#endif

max_stars_repo_name: lyase/barnes-hut-sne

max_issues_repo_name: lyase/barnes-hut-sne

max_forks_repo_name: lyase/barnes-hut-sne


text:
#ifndef WALK_H_
#define WALK_H_

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <gsl/gsl_rng.h>

#include "Lattice.h"

/* Random number generator parameters used across this program */
namespace WalkRNG {
    /* GSL random number generator, if using */
    extern gsl_rng *rng;

    /* Whether or not to use the GSL random number generator */
    extern bool use_GSL;

    /* Whether or not we've generated a new random seed since program start */
    extern bool random_seeded;
};

/* Class to perform a random walk on an N-dimensional lattice */
template<unsigned int N>
class Walk : public std::vector<Vector<N> > {
 private:
    Lattice<N> lattice;

 public:
 Walk(Lattice<N> lattice) : lattice(lattice) {
	/* If we haven't seeded the RNG yet... */
	if (!WalkRNG::random_seeded) {
	    if (!WalkRNG::use_GSL)
		/* If not using GSL's random number generator, seed with time */
		std::srand(std::time(0));
	    else {
		/* else set up GSL random number generator */
		const gsl_rng_type *T;

		gsl_rng_env_setup();

		T = gsl_rng_default;
		WalkRNG::rng = gsl_rng_alloc(T);
	    }
	    WalkRNG::random_seeded = true;
	}
    }

    /**
     * Generate random walk on the lattice of length `length` modifying in place
     * with non-basis transformed vectors and also returning the walk.
     * Clears the current walk.
     */
    Walk<N> &generate(int length) {
	this->clear();

	/* Get the new translation set for this lattice */
	std::vector<Vector<N> > translation_set = this->lattice.getTranslationSet();

	for (int i = 0; i < length; ++i) {
	    /* Get random index for the translation_set */
	    size_t random_index;
	    if (WalkRNG::use_GSL)
		random_index = gsl_rng_get(WalkRNG::rng) % translation_set.size();
	    else
		random_index = std::rand() % translation_set.size();

	    /* Get the element */
	    Vector<N> random_element = translation_set.at(random_index);

	    /* Add it to the walk */
	    this->push_back(random_element);
	}

	return *this;
    }

    /**
     * Step walk for one iteration, returning single vector and adding to Walk.
     * Same as generate(int length) with a length == 1, and doesn't clear the walk.
     * Returns the generated Vector<N>
     */
    Vector<N> step() {
	std::vector<Vector<N> > translation_set = this->lattice.getTranslationSet();

	size_t random_index;
	if (WalkRNG::use_GSL)
	    random_index = gsl_rng_get(WalkRNG::rng) % translation_set.size();
	else
	    random_index = std::rand() % translation_set.size();

	Vector<N> random_element = translation_set.at(random_index);

	this->push_back(random_element);

	return random_element;
    }

    /**
     * Apply this->lattice's basis to each vector in the walk
     */
    Walk &applyBasis() {
	for (int i = 0; i < this->size(); ++i) {
	    this->at(i) = this->lattice.applyBasis(this->at(i));
	}

	return *this;
    }

    /**
     * Get distance between start and end points of walk
     */
    double getDistance() {
	return this->getDistanceBetween(0, this->size());
    }

    /**
     * Get distance between inclusive point `start` and
     * exclusive point `end` of walk
     */
    double getDistanceBetween(size_t start, size_t end) {
	Vector<N> res;

	for (int i = start; i < end; ++i) {
	    res += this->at(i);
	}

	return res.getMagnitude();
    }

    /**
     * Return new walk with each step added, so each step corresponds to the
     * current position on the lattice.
     */
    Walk<N> accumulateVectors()
    {
	Walk<N> walk(this->lattice);
	Vector<N> running_total = this->at(0);

	// Initial value
	walk.push_back(running_total);

	/* Accumulate each vector to the running total, and push it
	 * to the walk */
	for (int i = 1; i < this->size(); ++i) {
	    running_total += this->at(i);
	    walk.push_back(running_total);
	}

	return walk;
    }

    /**
     * Return as string of newline-seperated comma seperated rows
     * with first column as x-component, second column as y-component, etc.
     */
    std::string toCSV() {
	std::stringstream ss;
	ss << *this;
	return ss.str();
    }

    // Display each vector in walk, one on each line
    friend std::ostream& operator<<(std::ostream& os, const Walk<N> &walk) {
	for (int i = 0; i < walk.size(); ++i)
	    os << walk.at(i) << std::endl;

	return os;
    }
};

#endif /* WALK_H_ */

max_stars_repo_name: joebentley/walk-gen

max_issues_repo_name: joebentley/walk-gen

max_forks_repo_name: joebentley/walk-gen