text:
// Implementation of interface described in viewer_image.h.

#include <math.h>
#include <unistd.h>

#include <gsl/gsl_statistics_double.h>

#include "meta_read_wrapper.h"
#include "utilities.h"
#include "viewer_image.h"

// We use a seperate thread to fill in the stuff that takes some time
// to create.
static gpointer
create_pyramid (ViewerImage *self)
{
  g_assert (self->pixmaps != NULL);

  if ( self->pyramid_cache != NULL ) {
    g_assert (self->scratch_dir == NULL);
    // We never unref this pyramid.  The reason is that we don't want
    // to keep this thread alive just for that purpose, and since the
    // unreffing the pyramid triggers a pyramid_cache_release, it has
    // to be done by its creator.  It will get reclaimed anyway by the
    // OS.
    self->pyramid
      = pyramid_new_using_cache (self->base_name->str, self->pixmaps,
				 self->pyramid_cache);
  }
  else {
    g_assert (self->scratch_dir != NULL);
    g_assert (self->pyramid_cache == NULL);
    self->pyramid = pyramid_new (self->base_name->str, self->pixmaps,
				 self->scratch_dir->str);
  }

  // Uncomment to get debugging snapshots of the pyramid layers.
  // pyramid_export_as_float_blobs (self->pyramid, "pyramid_snapshot");
  // pyramid_export_as_jpeg_files (self->pyramid, "pyramid_snapshot");

  g_mutex_lock (self->pyramid_mutex);
  self->pyramid_ready = TRUE;
  g_mutex_unlock (self->pyramid_mutex);

  // This thread isn't even joinable, but we keep the compiler happy
  // with a return statement anyway.
  return NULL;
}

ViewerImage *
viewer_image_new (const char *base_name, PixMapsSpec *pixmaps, gdouble sigmas,
		  size_t x_offset, size_t y_offset,
		  PyramidCache *pyramid_cache)
{
  ViewerImage *self = g_new (ViewerImage, 1);

  self->base_name = g_string_new (base_name);

  GString *mn = g_string_append (g_string_new (base_name), ".meta");
  GString *dn = g_string_append (g_string_new (base_name), ".img");

  self->meta = meta_read_wrapper (mn->str);
  meta_parameters *md = self->meta; // Convenience alias.

  self->size_x = md->general->sample_count;
  self->size_y = md->general->line_count;

  self->x_offset = x_offset;
  self->y_offset = y_offset;

  self->pixmaps = pix_maps_spec_ref (pixmaps);

  self->sigmas = sigmas;

  g_assert (self->size_x <= INT_MAX);
  gint sx = self->size_x;
  g_assert (self->size_y <= INT_MAX);
  gint sy = self->size_y;

  FILE *df = fopen (dn->str, "r");
  g_assert (df != NULL);

  // We will form the quick look by sampling every pixel_strid'th
  // pixel.  This number reflects a critical tradeoff -- bigger means
  // that we won't have a quick view of sufficiently high resolution
  // ready when the user tries to zoom in on an image, and interactive
  // performance will suffer when we create one on the fly.  On the
  // other hand, setting this number too small may run some machines
  // out of RAM and into virtual memory, which is also likely to
  // degrade interactive performance significantly.  Both these
  // problems go away once we have the pyramid fully biult, but of
  // course this takes a few seconds to do.
  const size_t pixel_stride = 8;

  self->quick_view_stride = pixel_stride;

  // Quick view dimensions.
  self->qvw = ceil ((double) sx / pixel_stride);
  self->qvh = ceil ((double) sy / pixel_stride);

  // Memory for quick_view.
  self->quick_view = g_new (float, self->qvw * self->qvh);

  // Storage for one image row.
  float *row_buffer = g_new (float, sx);

  size_t ii, jj;
  size_t cqvr = 0, cqvc = 0;		// Current quick_view row and column.
  for ( jj = 0 ; jj < sy ; jj += pixel_stride, cqvr++ ) {
    cqvc = 0;
    size_t read_count = fread (row_buffer, sizeof (float), sx, df);
    g_assert (read_count == sx);
    int return_code = fseek (df, sizeof (float)  * sx * (pixel_stride - 1),
			     SEEK_CUR);
    g_assert (return_code == 0);
    for ( ii = 0 ; ii < sx ; ii += pixel_stride, cqvc++ ) {
      float pv = row_buffer[ii];
      if ( G_BYTE_ORDER == G_LITTLE_ENDIAN ) {
	swap_bytes_32 ((unsigned char *) &pv);
      }
      (self->quick_view)[cqvr * self->qvw + cqvc] = pv;
    }
  }

  size_t qvpc = self->qvw * self->qvh;   // Quick view pixel count.

  // Apply pixmaps.
  for ( ii = 0 ; ii < qvpc ; ii++ ) {
    (self->quick_view)[ii]
      = pix_maps_spec_map (self->pixmaps, (self->quick_view)[ii]);
  }

  // Double precision version of quickview data so we can use GSL
  // library functions for statistics.
  double *qvad = g_new (double, qvpc);
  for ( ii = 0 ; ii < qvpc ; ii++ ) {
    qvad[ii] = (self->quick_view)[ii];
  }

  // Compute statistics for OpenGL to use for scaling.  IMPROVEME: it
  // might be possible to use OpenGL to do this or near enough using
  // the faster functions from the OpenGL imaging subset.
  gsl_stats_minmax (&(self->min), &(self->max), qvad, 1, qvpc);
  self->mean = gsl_stats_mean (qvad, 1, qvpc);
  self->sdev = gsl_stats_sd_m (qvad, 1, qvpc, self->mean);

  g_free (qvad);

  g_free (row_buffer);

  if ( pyramid_cache == NULL ) {
    self->pyramid_cache = NULL;
    self->scratch_dir = g_string_new ("/tmp/");
  }
  else {
    self->pyramid_cache = pyramid_cache_ref (pyramid_cache);
    self->scratch_dir = NULL;
  }

  if ( !g_thread_supported () ) {
    g_thread_init (NULL);
  }

  self->pyramid_ready = FALSE;
  self->pyramid_mutex = g_mutex_new ();
  self->pyramid = NULL;

  GThread *pyramid_thread
    = g_thread_create ((GThreadFunc) create_pyramid, self, FALSE, NULL);
  g_assert (pyramid_thread != NULL);
  // We try to keep things snappy for the use during the time before
  // the pyramids are finished.
  g_thread_set_priority (pyramid_thread, G_THREAD_PRIORITY_LOW);

  g_string_free (mn, TRUE);
  g_string_free (dn, TRUE);

  self->reference_count = 1;

  return self;
}

Pyramid *
viewer_image_pyramid (ViewerImage *self)
{
  Pyramid *result = NULL;
  g_mutex_lock (self->pyramid_mutex);
  if ( self->pyramid_ready ) {
    result = self->pyramid;
  }
  g_mutex_unlock (self->pyramid_mutex);

  return result;
}

ViewerImage *
viewer_image_ref (ViewerImage *self)
{
  self->reference_count++;

  return self;
}

void
viewer_image_unref (ViewerImage *self)
{
  self->reference_count--;

  if ( self->reference_count == 0 ) {
    while ( viewer_image_pyramid (self) == NULL ) {
      // oo icky poo a polling loop in a destructor!
      const guint sleep_time = 1;
      sleep (sleep_time);
    }
    viewer_image_free (self);
  }
}

gboolean
viewer_image_free (ViewerImage *self)
{
  if ( viewer_image_pyramid (self) == NULL ) {
    return FALSE;
  }

  g_string_free (self->base_name, TRUE);
  meta_free (self->meta);
  pix_maps_spec_unref (self->pixmaps);
  g_free (self->quick_view);
  if ( self->pyramid_cache != NULL ) {
    g_assert (self->scratch_dir == NULL);
    pyramid_cache_unref (self->pyramid_cache);
    // Note that we deliberately don't unref the pyramid here; see
    // comments in create_pyramid().
  }
  else {
    g_assert (self->scratch_dir != NULL);
    g_assert (self->pyramid_cache == NULL);
    g_string_free (self->scratch_dir, TRUE);
    pyramid_unref (self->pyramid);
  }
  g_mutex_free (self->pyramid_mutex);
  g_free (self);

  return TRUE;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/*
   Copyright [2019,2021] [IBM Corporation]
   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at
       http://www.apache.org/licenses/LICENSE-2.0
   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.
*/

#ifndef __CCPM_CRASH_CONSISTENT_ALLOCATOR_H__
#define __CCPM_CRASH_CONSISTENT_ALLOCATOR_H__

#include <ccpm/interfaces.h>
#include <common/exceptions.h>
#include <gsl/pointers>
#include <iosfwd>
#include <memory>
#include <string>
#include <vector>

namespace ccpm
{

	struct area_top;

	struct cca
		: public IHeap_expandable
	{
		using byte_span = common::byte_span;
		using persist_type = gsl::not_null<persister *>;

		/*
		 * An objection arose to requiring multiple arguments for the cca constructor.
		 * This struct allows construction from a single argument in all cases.
		 * cca.h may not be the best place for this, as cca does not itsenf use the struct.
		 */
		struct ctor_args
		{
			persist_type persister;
			region_span regions;
			bool force_init;
			ownership_callback_type callee_owns;
		};

	private:
		using top_vec_t = std::vector<std::unique_ptr<area_top>>;
		top_vec_t _top;
		top_vec_t::difference_type _last_top_allocate;
		top_vec_t::difference_type _last_top_free;
		persist_type _persist;

		void init(
			region_span regions
			, ownership_callback_type resolver
			, bool force_init
		);
	public:
		explicit cca(persist_type persist, region_span regions, ownership_callback_type resolver);

		explicit cca(persist_type persist, region_span regions);

		explicit cca(persist_type persist);

		~cca();

		cca(const cca &) = delete;
		const cca& operator=(const cca &) = delete;

		bool reconstitute(
			region_span regions
			, ownership_callback_type resolver
			, bool force_init
		) override;

		status_t allocate(
			void * & ptr_
			, std::size_t bytes_
			, std::size_t alignment_
		) override;

    void * allocate(std::size_t bytes_, std::size_t alignment_ = 1) {
      void * ptr = nullptr;
      if(allocate(ptr, bytes_, alignment_) != 0)
        throw General_exception("ccpm::cca::allocate failed");
      return ptr;
    }

    void * allocate_root(std::size_t bytes_, std::size_t alignment_ = 8) {
      void * ptr = nullptr;
      if(allocate(ptr, bytes_, alignment_) != 0)
        throw General_exception("ccpm::cca::allocate failed");
      set_root(common::make_byte_span(ptr, bytes_));
      return ptr;
    }

		status_t free(
			void * & ptr_
			, std::size_t bytes_
		) override;

		void add_regions(
			region_span
		) override;

		bool includes(
			const void *addr
		) const override;

		status_t remaining(
			std::size_t & out_size_
		) const override;

		region_vector_t get_regions() const override;

    void set_root(
      byte_span
    );

    byte_span get_root() const;

		void print(std::ostream &, const std::string &title = "cca") const;
	};
}

#endif

ext: h

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
#include "vector.h"

#include "math.h"
#include "mem.h"
#include "string_utils.h"
#include "ut_utils.h"

typedef struct {
    char *name;
    char age;
} nage_t;

void free_nage_t(void *pnage)
{
    nage_t *nage = pnage;
    ufree(nage->name);
    ufree(nage);
}

void *copy_nage_t(const void *pnage)
{
    const nage_t *nage = pnage;

    nage_t *n = umalloc(sizeof(*n));
    UASSERT(n);

    n->name = ustring_dup(nage->name);
    n->age = nage->age;

    return n;
}

void test_uvector_copy(bool verbose)
{
    (void)verbose;

    uvector_t *v = uvector_create();
    uvector_set_void_destroyer(v, free_nage_t);
    uvector_set_void_copier(v, copy_nage_t);

    nage_t *n = umalloc(sizeof(*n));

    *n = (nage_t){ustring_dup("john"), 34};
    uvector_append(v, G_PTR(n)); // v owns *n data

    uvector_t *v2 = uvector_copy(v); // v2 has shallow copy of *n2
    uvector_t *v3 = uvector_deep_copy(v); // v3 has deep copy of *n2, must set own

    uvector_destroy(v);
    uvector_destroy(v2);
    uvector_destroy(v3);
}

void test_uvector_bsearch(void)
{
    uvector_t *v = uvector_create();
    uvector_append(v, G_INT(1));
    uvector_append(v, G_INT(2));
    uvector_append(v, G_INT(4));
    uvector_append(v, G_INT(3));

    uvector_sort(v); // [1, 2, 3, 4]]
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(1)), 0);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(2)), 1);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(3)), 2);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(4)), 3);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(400)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(-400)), SIZE_MAX);

    uvector_remove_at(v, 1); // [1, 3, 4]
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(1)), 0);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(2)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(3)), 1);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(4)), 2);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(400)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(-400)), SIZE_MAX);

    uvector_pop_back(v); // [1, 3]
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(1)), 0);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(2)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(3)), 1);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(4)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(400)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(-400)), SIZE_MAX);

    uvector_pop_at(v, 0); // [3]
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(1)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(2)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(3)), 0);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(4)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(400)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(-400)), SIZE_MAX);

    uvector_pop_at(v, 0);
    UASSERT_INT_EQ(uvector_get_size(v), 0);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(1)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(2)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(3)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(4)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(400)), SIZE_MAX);
    UASSERT_INT_EQ(uvector_bsearch(v, G_INT(-400)), SIZE_MAX);

    uvector_destroy(v);
}

void test_uvector_next_permutation(bool verbose)
{
    uvector_t *v = uvector_create();
    uvector_append(v, G_INT(4));
    uvector_append(v, G_INT(3));
    uvector_append(v, G_INT(1));
    uvector_append(v, G_INT(2));

    do {
        if (verbose)
            uvector_print(v);
    } while (uvector_next_permutation(v));

    uvector_destroy(v);
}

void test_uvector_serialization(bool verbose)
{
    char *vs;
    uvector_t *v = uvector_create();

    vs = uvector_as_str(v);
    UASSERT_STR_EQ(vs, "[]");
    ufree(vs);

    uvector_append(v, G_INT(11));
    uvector_append(v, G_INT(22));
    uvector_append(v, G_INT(33));
    uvector_append(v, G_STR(ustring_dup("44")));
    uvector_append(v, G_CSTR("5"));

    uvector_append(v, G_VECTOR(uvector_deep_copy(v)));
    uvector_append(v, G_VECTOR(uvector_deep_copy(v)));

    if (verbose)
        uvector_print(v);
    vs = uvector_as_str(v);
    UASSERT_STR_EQ(vs,
                  "[11, 22, 33, \"44\", \"5\", [11, 22, 33, \"44\", \"5\"], [11, 22, 33, \"44\", \"5\", [11, 22, 33, \"44\", \"5\"]]]");
    ufree(vs);
    uvector_destroy(v);
}

void test_uvector_api()
{
    uvector_t *v = uvector_create();
    uvector_append(v, G_INT(0));
    UASSERT(uvector_get_capacity(v) >= VECTOR_INITIAL_CAPACITY);
    uvector_reserve_capacity(v, 1000);
    UASSERT_INT_EQ(uvector_get_capacity(v), 1000);
    uvector_destroy(v);

    v = uvector_create_with_size(3, G_NULL());
    UASSERT_INT_EQ(uvector_get_size(v), 3);

    uvector_append(v, G_INT(11));
    uvector_append(v, G_INT(22));
    uvector_append(v, G_INT(33));
    UASSERT_INT_EQ(uvector_get_size(v), 6);
    UASSERT(G_IS_NULL(uvector_get_at(v, 0)));
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 3)), 11);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 4)), 22);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 5)), 33);

    uvector_swap(v, 4, 5);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 4)), 33);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 5)), 22);

    uvector_resize(v, 0, G_NULL());
    UASSERT(uvector_is_empty(v));

    uvector_append(v, G_INT(44));
    UASSERT_INT_EQ(uvector_get_size(v), 1);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 0)), 44);

    uvector_set_at(v, 0, G_INT(2222));
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 0)), 2222);

    ugeneric_t *a = uvector_get_cells(v);
    UASSERT_INT_EQ(G_AS_INT(a[0]), 2222);
    uvector_destroy(v);

    v = uvector_create();
    UASSERT_INT_EQ(uvector_get_size(v), 0);
    UASSERT(uvector_is_empty(v));
    uvector_append(v, G_INT(11));
    UASSERT(!uvector_is_empty(v));
    UASSERT_INT_EQ(uvector_get_size(v), 1);
    UASSERT_INT_EQ(G_AS_INT(uvector_pop_back(v)), 11);
    UASSERT(uvector_is_empty(v));
    UASSERT_INT_EQ(uvector_get_size(v), 0);

    uvector_resize(v, 0, G_NULL());
    uvector_append(v, G_INT(200));
    uvector_append(v, G_INT(100));
    UASSERT_INT_EQ(uvector_get_size(v), 2);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 0)), 200);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 1)), 100);
    uvector_insert_at(v, 1, G_INT(300));
    UASSERT_INT_EQ(uvector_get_size(v), 3);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 0)), 200);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 1)), 300);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 2)), 100);
    UASSERT_INT_EQ(G_AS_INT(uvector_pop_at(v, 0)), 200);
    UASSERT_INT_EQ(uvector_get_size(v), 2);
    UASSERT_INT_EQ(G_AS_INT(uvector_pop_at(v, 0)), 300);
    UASSERT_INT_EQ(uvector_get_size(v), 1);
    UASSERT_INT_EQ(G_AS_INT(uvector_pop_at(v, 0)), 100);
    UASSERT_INT_EQ(uvector_get_size(v), 0);
    uvector_resize(v, 1000, G_NULL());
    UASSERT_INT_EQ(uvector_get_size(v), 1000);
    uvector_destroy(v);

    v = uvector_create();
    uvector_append(v, G_STR(ustring_dup("1")));
    uvector_append(v, G_STR(ustring_dup("2")));
    uvector_append(v, G_STR(ustring_dup("3")));
    uvector_append(v, G_STR(ustring_dup("4")));
    uvector_resize(v, 2, G_NULL());
    uvector_resize(v, 1, G_NULL());
    uvector_clear(v);
    uvector_append(v, G_STR(ustring_dup("5")));
    UASSERT_INT_EQ(uvector_get_size(v), 1);
    uvector_destroy(v);

    v = uvector_create();
    uvector_append(v, G_STR(ustring_dup("string")));
    uvector_set_at(v, 0, G_STR(ustring_dup("string2")));
    uvector_destroy(v);

    v = uvector_create();
    uvector_append(v, G_INT(123));
    UASSERT_INT_EQ(G_AS_INT(uvector_get_at(v, 0)), 123);
    UASSERT_INT_EQ(G_AS_INT(uvector_get_back(v)), 123);
    uvector_t *vcopy = uvector_copy(v);
    UASSERT(uvector_compare(v, vcopy) == 0);
    uvector_destroy(v);
    uvector_destroy(vcopy);

    v = uvector_create();
    uvector_reserve_capacity(v, 1024);
    uvector_shrink_to_size(v);
    uvector_destroy(v);

    v = uvector_create();
    uvector_append(v, G_STR(ustring_dup("1")));
    uvector_append(v, G_INT(42));
    uvector_append(v, G_BOOL(false));
    UASSERT(uvector_contains(v, G_STR("1")));
    UASSERT(uvector_contains(v, G_CSTR("1")));
    UASSERT(uvector_contains(v, G_BOOL(false)));
    UASSERT(uvector_contains(v, G_INT(42)));
    UASSERT(!uvector_contains(v, G_STR("2")));
    UASSERT(!uvector_contains(v, G_BOOL(true)));
    UASSERT(!uvector_contains(v, G_INT(100)));
    uvector_destroy(v);

#if defined(__unix__) && defined(ENABLE_UASSERT_INPUT)
    // Check asserts.
    UASSERT_ABORTS(uvector_append(NULL, G_INT(0)));
    UASSERT_ABORTS(uvector_pop_back(NULL));
    UASSERT_ABORTS(uvector_get_cells(NULL));
    UASSERT_ABORTS(uvector_get_size(NULL));
    UASSERT_ABORTS(uvector_resize(NULL, 100, G_NULL()));
    UASSERT_ABORTS(uvector_reserve_capacity(NULL, 100));
    UASSERT_ABORTS(uvector_get_capacity(NULL));
    UASSERT_ABORTS(uvector_is_empty(NULL));
    UASSERT_ABORTS(uvector_swap(NULL, 1, 2));
    UASSERT_ABORTS(uvector_insert_at(NULL, 1, G_INT(100)));
    UASSERT_ABORTS(uvector_get_at(NULL, 1));
    UASSERT_ABORTS(uvector_set_at(NULL, 1, G_INT(0)));
    uvector_destroy(NULL);

    v = uvector_create_with_size(0, G_NULL());
    UASSERT_ABORTS(uvector_get_at(v, 1));
    UASSERT_ABORTS(uvector_set_at(v, 1, G_INT(0)));

    uvector_resize(v, 100, G_NULL());
    UASSERT_ABORTS(uvector_swap(v, 150, 250));
    UASSERT_ABORTS(uvector_swap(v, 50, 250));
    UASSERT_ABORTS(uvector_insert_at(v, 500, G_INT(400)));
    uvector_destroy(v);
#endif

}

int cmp(const char *s1, const char *s2)
{
    ugeneric_t gv1 = ugeneric_parse(s1);
    ugeneric_t gv2 = ugeneric_parse(s2);
    UASSERT_NO_ERROR(gv1);
    UASSERT_NO_ERROR(gv2);
    uvector_t *v1 = G_AS_PTR(gv1);
    uvector_t *v2 = G_AS_PTR(gv2);
    int diff = uvector_compare(v1, v2);
    uvector_destroy(v1);
    uvector_destroy(v2);
    return diff;
}

void test_uvector_compare(void)
{
    uvector_t *v1 = uvector_create();
    uvector_t *v2 = uvector_create();
    uvector_append(v1, G_STR(ustring_dup("s1")));
    uvector_append(v2, G_STR(ustring_dup("s1")));
    UASSERT(uvector_compare(v1, v2) == 0);
    uvector_set_at(v1, 0, G_STR(ustring_dup("s3")));
    UASSERT(uvector_compare(v1, v2) != 0);
    uvector_set_at(v2, 0, G_STR(ustring_dup("s4")));
    UASSERT(uvector_compare(v1, v2) < 0);
    uvector_destroy(v1);
    uvector_destroy(v2);

    UASSERT(cmp("[]", "[]") == 0);
    UASSERT(cmp("[1]", "[1]") == 0);
    UASSERT(cmp("[1, 2, 3]", "[1, 2, 3]") == 0);
    UASSERT(cmp("[1, 2]", "[1, 2, 3]") < 0);
    UASSERT(cmp("[1, 2, 3]", "[1, 2]") > 0);
    UASSERT(cmp("[3, 2, 3]", "[1, 2, 3]") > 0);
}

/*
double f(double x, void *p)
{
    uvector_t *v = p;
    return G_AS_REAL(uvector_get_at(v, x));
}

void test_gnuplot(void)
{
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>

    gsl_function F;
    double result, abserr;
    #define N 10

    gnuplot_attrs_t a = {
        .xlabel = "xlabel",
        .ylabel = "ylabel",
        .data_label = "samples",
        .title = "chart",
    };

    uvector_t *v = uvector_create();
    for (size_t i = 0; i < N; i++)
    {
        uvector_append(v, G_REAL(sin(i/50.0)));
    }
    uvector_dump_to_gnuplot(v, &a, stdout);

    uvector_t *dv = uvector_create();
    F.function = &f;
    F.params = v;

    for (size_t i = 0; i < N; i++)
    {
        UASSERT(gsl_deriv_central(&F, (double)i, 1e-8, &result, &abserr) == 0);
        uvector_append(dv, G_REAL(result));
    }
    uvector_dump_to_gnuplot(dv, &a, stdout);

    uvector_destroy(v);
    uvector_destroy(dv);
}
*/

void _check_slice(const uvector_t *v, size_t b, size_t e, size_t s, char *exp)
{
    uvector_t *slice = uvector_get_slice(v, b, e, s);
    char *str = uvector_as_str(slice);
    UASSERT_STR_EQ(str, exp);
    ufree(str);
    uvector_destroy(slice);
}

void test_uvector_slice(void)
{
    ugeneric_t g = ugeneric_parse("[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]");
    UASSERT_NO_ERROR(g);
    uvector_t *v = G_AS_PTR(g);

    _check_slice(v, 0, 10, 1, "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]");
    _check_slice(v, 0, 10, 2, "[1, 3, 5, 7, 9]");
    _check_slice(v, 1, 10, 2, "[2, 4, 6, 8, 10]");
    _check_slice(v, 0, 10, 3, "[1, 4, 7, 10]");
    _check_slice(v, 1, 10, 3, "[2, 5, 8]");
    _check_slice(v, 0, 10, 4, "[1, 5, 9]");
    _check_slice(v, 1, 10, 4, "[2, 6, 10]");
    _check_slice(v, 0, 10, 5, "[1, 6]");
    _check_slice(v, 1, 10, 5, "[2, 7]");
    _check_slice(v, 0, 10, 6, "[1, 7]");
    _check_slice(v, 1, 10, 6, "[2, 8]");
    _check_slice(v, 0, 10, 9, "[1, 10]");
    _check_slice(v, 1, 10, 9, "[2]");
    _check_slice(v, 1, 1, 1, "[]");
    _check_slice(v, 1, 1, 5, "[]");
    _check_slice(v, 0, 10, 25, "[1]");

    uvector_destroy(v);
}

void test_uvector_data_ownership(void)
{
    uvector_t *v = uvector_create();
    UASSERT(uvector_is_data_owner(v));
    //uvector_append(v, G_CSTR("s1"));
    //uvector_append(v, G_CSTR("s1"));

    // shallow copy
    uvector_t *copy1 = uvector_copy(v);
    UASSERT(!uvector_is_data_owner(copy1));

    // deep copy
    uvector_t *copy2 = uvector_deep_copy(v);
    UASSERT(uvector_is_data_owner(copy2));
    uvector_drop_data_ownership(copy2);
    UASSERT(!uvector_is_data_owner(copy2));

    // kick the bucket
    uvector_destroy(v);
    uvector_destroy(copy1);
    uvector_destroy(copy2);
}

void _check_reverse(const char *in, const char *rev)
{
    ugeneric_t g = ugeneric_parse(in);
    UASSERT_NO_ERROR(g);
    uvector_t *v = G_AS_PTR(g);
    uvector_reverse(v);
    char *t = uvector_as_str(v);
    UASSERT_STR_EQ(t, rev);
    ufree(t);
    uvector_destroy(v);
}

void test_uvector_reverse(void)
{
    _check_reverse("[]", "[]");
    _check_reverse("[1]", "[1]");
    _check_reverse("[1, 2]", "[2, 1]");
    _check_reverse("[1, 2, 3]", "[3, 2, 1]");
    _check_reverse("[1, 2, 3, 4]", "[4, 3, 2, 1]");
    _check_reverse("[1, 2, 3, 4, 5]", "[5, 4, 3, 2, 1]");
}

int main(int argc, char **argv)
{
//    test_gnuplot();

    (void)argv;
    test_uvector_api();
    test_uvector_serialization(argc > 1);
    test_uvector_next_permutation(argc > 1);
    test_uvector_copy(argc > 1);
    test_uvector_bsearch();
    test_uvector_compare();
    test_uvector_slice();
    test_uvector_data_ownership();
    test_uvector_reverse();

    return EXIT_SUCCESS;
}

ext: c

max_stars_repo_name: YauheniKaratsevich/ugeneric

max_issues_repo_name: YauheniKaratsevich/ugeneric

max_forks_repo_name: YauheniKaratsevich/ugeneric


text:
#pragma once

#include "internal/Buffer.h"
#include <gsl/span>

namespace ge::gl4
{

using namespace gl45;

class BufferView
{
public:
	BufferView(Buffer &buffer, GLenum type, gsl::span<uint8_t> range);

	// Mapping Functions
	auto map(BufferAccessMask accessFlags) -> void *;
	void flush();
	void unmap();

	void bindToIndex(GLuint index);
	static void unbind(GLenum target, GLuint index);

	auto getBuffer() const -> Buffer &;

	const GLenum type;
	const gsl::span<uint8_t> range;
private:
	Buffer *buffer_;
};

}

max_stars_repo_name: mmha/gravity

max_issues_repo_name: mmha/gravity

max_forks_repo_name: mmha/gravity


text:
/**
 *
 * @file core_slaswp.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

#define A(m, n) BLKADDR(descA, float, m, n)

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_slaswp performs a series of row interchanges on the matrix A.
 *  One row interchange is initiated for each of rows I1 through I2 of A.
 *
 *******************************************************************************
 *
 *  @param[in] N
 *          The number of columns in the matrix A. N >= 0.
 *
 *  @param[in,out] A
 *         On entry, the matrix of column dimension N to which the row
 *         interchanges will be applied.
 *         On exit, the permuted matrix.
 *
 *  @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,max(IPIV[I1..I2])).
 *
 *  @param[in] I1
 *          The first element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] I2
 *          The last element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] IPIV
 *          The pivot indices; Only the element in position i1 to i2
 *          are accessed. The pivot are offset by A.i.
 *
 *  @param[in] INC
 *          The increment between successive values of IPIV.  If IPIV
 *          is negative, the pivots are applied in reverse order.
 *
 *******************************************************************************
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp = PCORE_slaswp
#define CORE_slaswp PCORE_slaswp
#endif
void CORE_slaswp(int N, float *A, int LDA, int I1, int I2, const int *IPIV, int INC)
{
    LAPACKE_slaswp_work( LAPACK_COL_MAJOR, N, A, LDA, I1, I2, IPIV, INC );
}

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_slaswp_ontile apply the slaswp function on a matrix stored in
 *  tile layout
 *
 *******************************************************************************
 *
 *  @param[in,out] descA
 *          The descriptor of the matrix A to permute.
 *
 *  @param[in] i1
 *          The first element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] i2
 *          The last element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] ipiv
 *          The pivot indices; Only the element in position i1 to i2
 *          are accessed. The pivot are offset by A.i.
 *
 *  @param[in] inc
 *          The increment between successive values of IPIV.  If IPIV
 *          is negative, the pivots are applied in reverse order.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 *******************************************************************************
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp_ontile = PCORE_slaswp_ontile
#define CORE_slaswp_ontile PCORE_slaswp_ontile
#endif
int CORE_slaswp_ontile(PLASMA_desc descA, int i1, int i2, const int *ipiv, int inc)
{
    int i, j, ip, it;
    float *A1;
    int lda1, lda2;

    /* Change i1 to C notation */
    i1--;

    /* Check parameters */
    if ( descA.nt > 1 ) {
        coreblas_error(1, "Illegal value of descA.nt");
        return -1;
    }
    if ( i1 < 0 ) {
        coreblas_error(2, "Illegal value of i1");
        return -2;
    }
    if ( (i2 <= i1) || (i2 > descA.m) ) {
        coreblas_error(3, "Illegal value of i2");
        return -3;
    }
    if ( ! ( (i2 - i1 - i1%descA.mb -1) < descA.mb ) ) {
        coreblas_error(2, "Illegal value of i1,i2. They have to be part of the same block.");
        return -3;
    }

    if (inc > 0) {
        it = i1 / descA.mb;
        A1 = A(it, 0);
        lda1 = BLKLDD(descA, 0);

        for (j = i1; j < i2; ++j, ipiv+=inc) {
            ip = (*ipiv) - descA.i - 1;
            if ( ip != j )
            {
                it = ip / descA.mb;
                i  = ip % descA.mb;
                lda2 = BLKLDD(descA, it);
                cblas_sswap(descA.n, A1       + j, lda1,
                                     A(it, 0) + i, lda2 );
            }
        }
    }
    else
    {
        it = (i2-1) / descA.mb;
        A1 = A(it, 0);
        lda1 = BLKLDD(descA, it);

        i1--;
        ipiv = &ipiv[(1-i2)*inc];
        for (j = i2-1; j > i1; --j, ipiv+=inc) {
            ip = (*ipiv) - descA.i - 1;
            if ( ip != j )
            {
                it = ip / descA.mb;
                i  = ip % descA.mb;
                lda2 = BLKLDD(descA, it);
                cblas_sswap(descA.n, A1       + j, lda1,
                                     A(it, 0) + i, lda2 );
            }
        }
    }

    return PLASMA_SUCCESS;
}

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_sswptr_ontile apply the slaswp function on a matrix stored in
 *  tile layout, followed by a strsm on the first tile of the panel.
 *
 *******************************************************************************
 *
 *  @param[in,out] descA
 *          The descriptor of the matrix A to permute.
 *
 *  @param[in] i1
 *          The first element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] i2
 *          The last element of IPIV for which a row interchange will
 *          be done.
 *
 *  @param[in] ipiv
 *          The pivot indices; Only the element in position i1 to i2
 *          are accessed. The pivot are offset by A.i.
 *
 *  @param[in] inc
 *          The increment between successive values of IPIV.  If IPIV
 *          is negative, the pivots are applied in reverse order.
 *
 *  @param[in] Akk
 *          The triangular matrix Akk. The leading descA.nb-by-descA.nb
 *          lower triangular part of the array Akk contains the lower triangular matrix, and the
 *          strictly upper triangular part of A is not referenced. The
 *          diagonal elements of A are also not referenced and are assumed to be 1.
 *
 * @param[in] ldak
 *          The leading dimension of the array Akk. ldak >= max(1,descA.nb).
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 *******************************************************************************
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sswptr_ontile = PCORE_sswptr_ontile
#define CORE_sswptr_ontile PCORE_sswptr_ontile
#endif
int CORE_sswptr_ontile(PLASMA_desc descA, int i1, int i2, const int *ipiv, int inc,
                       const float *Akk, int ldak)
{
    float zone  = 1.0;
    int lda;
    int m = descA.mt == 1 ? descA.m : descA.mb;

    if ( descA.nt > 1 ) {
        coreblas_error(1, "Illegal value of descA.nt");
        return -1;
    }
    if ( i1 < 1 ) {
        coreblas_error(2, "Illegal value of i1");
        return -2;
    }
    if ( (i2 < i1) || (i2 > m) ) {
        coreblas_error(3, "Illegal value of i2");
        return -3;
    }

    CORE_slaswp_ontile(descA, i1, i2, ipiv, inc);

    lda = BLKLDD(descA, 0);
    cblas_strsm( CblasColMajor, CblasLeft, CblasLower,
                 CblasNoTrans, CblasUnit,
                 m, descA.n, (zone),
                 Akk,     ldak,
                 A(0, 0), lda );

    return PLASMA_SUCCESS;
}

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_slaswpc_ontile apply the slaswp function on a matrix stored in
 *  tile layout
 *
 *******************************************************************************
 *
 *  @param[in,out] descA
 *          The descriptor of the matrix A to permute.
 *
 *  @param[in] i1
 *          The first element of IPIV for which a column interchange will
 *          be done.
 *
 *  @param[in] i2
 *          The last element of IPIV for which a column interchange will
 *          be done.
 *
 *  @param[in] ipiv
 *          The pivot indices; Only the element in position i1 to i2
 *          are accessed. The pivot are offset by A.i.
 *
 *  @param[in] inc
 *          The increment between successive values of IPIV.  If IPIV
 *          is negative, the pivots are applied in reverse order.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 *******************************************************************************
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswpc_ontile = PCORE_slaswpc_ontile
#define CORE_slaswpc_ontile PCORE_slaswpc_ontile
#endif
int CORE_slaswpc_ontile(PLASMA_desc descA, int i1, int i2, const int *ipiv, int inc)
{
    int i, j, ip, it;
    float *A1;
    int lda;

    /* Change i1 to C notation */
    i1--;

    /* Check parameters */
    if ( descA.mt > 1 ) {
        coreblas_error(1, "Illegal value of descA.mt");
        return -1;
    }
    if ( i1 < 0 ) {
        coreblas_error(2, "Illegal value of i1");
        return -2;
    }
    if ( (i2 <= i1) || (i2 > descA.n) ) {
        coreblas_error(3, "Illegal value of i2");
        return -3;
    }
    if ( ! ( (i2 - i1 - i1%descA.nb -1) < descA.nb ) ) {
        coreblas_error(2, "Illegal value of i1,i2. They have to be part of the same block.");
        return -3;
    }

    lda = BLKLDD(descA, 0);

    if (inc > 0) {
        it = i1 / descA.nb;
        A1 = A(0, it);

        for (j = i1-1; j < i2; ++j, ipiv+=inc) {
            ip = (*ipiv) - descA.j - 1;
            if ( ip != j )
            {
                it = ip / descA.nb;
                i  = ip % descA.nb;
                cblas_sswap(descA.m, A1       + j*lda, 1,
                                     A(0, it) + i*lda, 1 );
            }
        }
    }
    else
    {
        it = (i2-1) / descA.mb;
        A1 = A(0, it);

        i1--;
        ipiv = &ipiv[(1-i2)*inc];
        for (j = i2-1; j > i1; --j, ipiv+=inc) {
            ip = (*ipiv) - descA.j - 1;
            if ( ip != j )
            {
                it = ip / descA.nb;
                i  = ip % descA.nb;
                cblas_sswap(descA.m, A1       + j*lda, 1,
                                     A(0, it) + i*lda, 1 );
            }
        }
    }

    return PLASMA_SUCCESS;
}


text:
//#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

static double fun(double x)
{
    return x + cos(x * x);
}

int main(void)
{
    const double a = 1.0;
    const double b = 10.0;
    const int steps = 10;
    double xi, yi, x[100], y[100], dx;
    FILE *input, *output;
    int i;

    input = fopen("wartosci.txt", "w");
    output = fopen("inter-cspline.txt", "w");

    dx = (b - a) / (double)steps;

    for (i = 0; i <= steps; ++i)
    {
        x[i] = a + (double)i * dx + 0.5 * sin((double)i * dx);
        y[i] = fun(x[i]);
        fprintf(input, "%g %g\n", x[i], y[i]);
    }

    {
        gsl_interp_accel *acc = gsl_interp_accel_alloc();

        gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, steps + 1);

        gsl_spline_init(spline, x, y, steps + 1);

        for (xi = a; xi <= b; xi += 0.01)
        {
            yi = gsl_spline_eval(spline, xi, acc);
            fprintf(output, "%g %g\n", xi, yi);
        }

        gsl_spline_free(spline);
        gsl_interp_accel_free(acc);
    }
    return 0;
}


text:
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the LICENSE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// 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.

#pragma once

#include <boost/multiprecision/cpp_int.hpp>
#include <gsl/gsl>
#include <spdlog/spdlog.h>

using uint256_t = boost::multiprecision::checked_uint256_t;
using byte = std::byte;

#include "inc/exception.h"
#include "inc/log.h"
#include "inc/string.h"

#define BETTER_ENUMS_DEFAULT_CONSTRUCTOR(Enum)                                 \
public:                                                                        \
  Enum() = default;

#define ENUM BETTER_ENUM

max_stars_repo_name: bandprotocol/bandprotocol

max_issues_repo_name: bandprotocol/bandchain-legacy

max_forks_repo_name: bandprotocol/bandchain-legacy


text:
/* specfunc/erfc.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:  J. Theiler (modifications by G. Jungman) */

/*
 * See Hart et al, Computer Approximations, John Wiley and Sons, New York (1968)
 * (This applies only to the erfc8 stuff, which is the part
 *  of the original code that survives. I have replaced much of
 *  the other stuff with Chebyshev fits. These are simpler and
 *  more precise than the original approximations. [GJ])
 */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_erf.h"

#include "check.h"

#include "chebyshev.h"
#include "cheb_eval.c"

#define LogRootPi_  0.57236494292470008706

static double erfc8_sum(double x)
{
  /* estimates erfc(x) valid for 8 < x < 100 */
  /* This is based on index 5725 in Hart et al */

  static double P[] = {
      2.97886562639399288862,
      7.409740605964741794425,
      6.1602098531096305440906,
      5.019049726784267463450058,
      1.275366644729965952479585264,
      0.5641895835477550741253201704
  };
  static double Q[] = {
      3.3690752069827527677,
      9.608965327192787870698,
      17.08144074746600431571095,
      12.0489519278551290360340491,
      9.396034016235054150430579648,
      2.260528520767326969591866945,
      1.0
  };
  double num=0.0, den=0.0;
  int i;

  num = P[5];
  for (i=4; i>=0; --i) {
      num = x*num + P[i];
  }
  den = Q[6];
  for (i=5; i>=0; --i) {
      den = x*den + Q[i];
  }

  return num/den;
}

inline
static double erfc8(double x)
{
  double e;
  e = erfc8_sum(x);
  e *= exp(-x*x);
  return e;
}

inline
static double log_erfc8(double x)
{
  double e;
  e = erfc8_sum(x);
  e = log(e) - x*x;
  return e;
}

#if 0
/* Abramowitz+Stegun, 7.2.14 */
static double erfcasympsum(double x)
{
  int i;
  double e = 1.;
  double coef = 1.;
  for (i=1; i<5; ++i) {
    /* coef *= -(2*i-1)/(2*x*x); ??? [GJ] */
    coef *= -(2*i+1)/(i*(4*x*x*x*x));
    e += coef;
    /*
    if (fabs(coef) < 1.0e-15) break;
    if (fabs(coef) > 1.0e10) break;

    [GJ]: These tests are not useful. This function is only
    used below. Took them out; they gum up the pipeline.
    */
  }
  return e;
}
#endif /* 0 */

/* Abramowitz+Stegun, 7.1.5 */
static int erfseries(double x, gsl_sf_result * result)
{
  double coef = x;
  double e    = coef;
  double del;
  int k;
  for (k=1; k<30; ++k) {
    coef *= -x*x/k;
    del   = coef/(2.0*k+1.0);
    e += del;
  }
  result->val = 2.0 / M_SQRTPI * e;
  result->err = 2.0 / M_SQRTPI * (fabs(del) + GSL_DBL_EPSILON);
  return GSL_SUCCESS;
}

/* Chebyshev fit for erfc((t+1)/2), -1 < t < 1
 */
static double erfc_xlt1_data[20] = {
  1.06073416421769980345174155056,
 -0.42582445804381043569204735291,
  0.04955262679620434040357683080,
  0.00449293488768382749558001242,
 -0.00129194104658496953494224761,
 -0.00001836389292149396270416979,
  0.00002211114704099526291538556,
 -5.23337485234257134673693179020e-7,
 -2.78184788833537885382530989578e-7,
  1.41158092748813114560316684249e-8,
  2.72571296330561699984539141865e-9,
 -2.06343904872070629406401492476e-10,
 -2.14273991996785367924201401812e-11,
  2.22990255539358204580285098119e-12,
  1.36250074650698280575807934155e-13,
 -1.95144010922293091898995913038e-14,
 -6.85627169231704599442806370690e-16,
  1.44506492869699938239521607493e-16,
  2.45935306460536488037576200030e-18,
 -9.29599561220523396007359328540e-19
};
static cheb_series erfc_xlt1_cs = {
  erfc_xlt1_data,
  19,
  -1, 1,
  12
};

/* Chebyshev fit for erfc(x) exp(x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1
 */
static double erfc_x15_data[25] = {
  0.44045832024338111077637466616,
 -0.143958836762168335790826895326,
  0.044786499817939267247056666937,
 -0.013343124200271211203618353102,
  0.003824682739750469767692372556,
 -0.001058699227195126547306482530,
  0.000283859419210073742736310108,
 -0.000073906170662206760483959432,
  0.000018725312521489179015872934,
 -4.62530981164919445131297264430e-6,
  1.11558657244432857487884006422e-6,
 -2.63098662650834130067808832725e-7,
  6.07462122724551777372119408710e-8,
 -1.37460865539865444777251011793e-8,
  3.05157051905475145520096717210e-9,
 -6.65174789720310713757307724790e-10,
  1.42483346273207784489792999706e-10,
 -3.00141127395323902092018744545e-11,
  6.22171792645348091472914001250e-12,
 -1.26994639225668496876152836555e-12,
  2.55385883033257575402681845385e-13,
 -5.06258237507038698392265499770e-14,
  9.89705409478327321641264227110e-15,
 -1.90685978789192181051961024995e-15,
  3.50826648032737849245113757340e-16
};
static cheb_series erfc_x15_cs = {
  erfc_x15_data,
  24,
  -1, 1,
  16
};

/* Chebyshev fit for erfc(x) x exp(x^2), 5 < x < 10, x = (5t + 15)/2, -1 < t < 1
 */
static double erfc_x510_data[20] = {
  1.11684990123545698684297865808,
  0.003736240359381998520654927536,
 -0.000916623948045470238763619870,
  0.000199094325044940833965078819,
 -0.000040276384918650072591781859,
  7.76515264697061049477127605790e-6,
 -1.44464794206689070402099225301e-6,
  2.61311930343463958393485241947e-7,
 -4.61833026634844152345304095560e-8,
  8.00253111512943601598732144340e-9,
 -1.36291114862793031395712122089e-9,
  2.28570483090160869607683087722e-10,
 -3.78022521563251805044056974560e-11,
  6.17253683874528285729910462130e-12,
 -9.96019290955316888445830597430e-13,
  1.58953143706980770269506726000e-13,
 -2.51045971047162509999527428316e-14,
  3.92607828989125810013581287560e-15,
 -6.07970619384160374392535453420e-16,
  9.12600607264794717315507477670e-17
};
static cheb_series erfc_x510_cs = {
  erfc_x510_data,
  19,
  -1, 1,
  12
};

#if 0
inline
static double
erfc_asymptotic(double x)
{
  return exp(-x*x)/x * erfcasympsum(x) / M_SQRTPI;
}
inline
static double
log_erfc_asymptotic(double x)
{
  return log(erfcasympsum(x)/x) - x*x - LogRootPi_;
}
#endif /* 0 */

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_erfc_e(double x, gsl_sf_result * result)
{
  const double ax = fabs(x);
  double e_val, e_err;

  /* CHECK_POINTER(result) */

  if(ax <= 1.0) {
    double t = 2.0*ax - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&erfc_xlt1_cs, t, &c);
    e_val = c.val;
    e_err = c.err;
  }
  else if(ax <= 5.0) {
    double ex2 = exp(-x*x);
    double t = 0.5*(ax-3.0);
    gsl_sf_result c;
    cheb_eval_e(&erfc_x15_cs, t, &c);
    e_val = ex2 * c.val;
    e_err = ex2 * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON);
  }
  else if(ax < 10.0) {
    double exterm = exp(-x*x) / ax;
    double t = (2.0*ax - 15.0)/5.0;
    gsl_sf_result c;
    cheb_eval_e(&erfc_x510_cs, t, &c);
    e_val = exterm * c.val;
    e_err = exterm * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON + GSL_DBL_EPSILON);
  }
  else {
    e_val = erfc8(ax);
    e_err = (x*x + 1.0) * GSL_DBL_EPSILON * fabs(e_val);
  }

  if(x < 0.0) {
    result->val  = 2.0 - e_val;
    result->err  = e_err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  }
  else {
    result->val  = e_val;
    result->err  = e_err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  }

  return GSL_SUCCESS;
}

int gsl_sf_log_erfc_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x*x < 10.0*GSL_ROOT6_DBL_EPSILON) {
    const double y = x / M_SQRTPI;
    /* series for -1/2 Log[Erfc[Sqrt[Pi] y]] */
    const double c3 = (4.0 - M_PI)/3.0;
    const double c4 = 2.0*(1.0 - M_PI/3.0);
    const double c5 = -0.001829764677455021;  /* (96.0 - 40.0*M_PI + 3.0*M_PI*M_PI)/30.0  */
    const double c6 =  0.02629651521057465;   /* 2.0*(120.0 - 60.0*M_PI + 7.0*M_PI*M_PI)/45.0 */
    const double c7 = -0.01621575378835404;
    const double c8 =  0.00125993961762116;
    const double c9 =  0.00556964649138;
    const double c10 = -0.0045563339802;
    const double c11 =  0.0009461589032;
    const double c12 =  0.0013200243174;
    const double c13 = -0.00142906;
    const double c14 =  0.00048204;
    double series = c8 + y*(c9 + y*(c10 + y*(c11 + y*(c12 + y*(c13 + c14*y)))));
    series = y*(1.0 + y*(1.0 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*(c7 + y*series)))))));
    result->val = -2.0 * series;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  /*
  don't like use of log1p(); added above series stuff for small x instead, should be ok [GJ]
  else if (fabs(x) < 1.0) {
    gsl_sf_result result_erf;
    gsl_sf_erf_e(x, &result_erf);
    result->val  = log1p(-result_erf.val);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  */
  else if(x > 8.0) {
    result->val = log_erfc8(x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result result_erfc;
    gsl_sf_erfc_e(x, &result_erfc);
    result->val  = log(result_erfc.val);
    result->err  = fabs(result_erfc.err / result_erfc.val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_erf_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(fabs(x) < 1.0) {
    return erfseries(x, result);
  }
  else {
    gsl_sf_result result_erfc;
    gsl_sf_erfc_e(x, &result_erfc);
    result->val  = 1.0 - result_erfc.val;
    result->err  = result_erfc.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_erf_Z_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    const double ex2 = exp(-x*x/2.0);
    result->val  = ex2 / (M_SQRT2 * M_SQRTPI);
    result->err  = fabs(x * result->val) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_erf_Q_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    gsl_sf_result result_erfc;
    int stat = gsl_sf_erfc_e(x/M_SQRT2, &result_erfc);
    result->val  = 0.5 * result_erfc.val;
    result->err  = 0.5 * result_erfc.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_erfc(double x)
{
  EVAL_RESULT(gsl_sf_erfc_e(x, &result));
}

double gsl_sf_log_erfc(double x)
{
  EVAL_RESULT(gsl_sf_log_erfc_e(x, &result));
}

double gsl_sf_erf(double x)
{
  EVAL_RESULT(gsl_sf_erf_e(x, &result));
}

double gsl_sf_erf_Z(double x)
{
  EVAL_RESULT(gsl_sf_erf_Z_e(x, &result));
}

double gsl_sf_erf_Q(double x)
{
  EVAL_RESULT(gsl_sf_erf_Q_e(x, &result));
}


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "matrix_operations.c"
#include "ellipse/parametric_function_roots.c"
#include "polynomial.c"

void characteristic_ellipse_matrix(double *X, double *R, double phi, double exponent)
{

    // rotation matrix
    double Q[2][2] = {{ 0 }};
    double Qt[2][2] = {{ 0 }};
    Q[0][0] = cos(phi);
    Q[0][1] = sin(phi);
    Q[1][0] = -sin(phi);
    Q[1][1] = cos(phi);
    matrix_transpose(&Qt[0][0], &Q[0][0], 2, 2);

    // radii matrix
    double O[2][2] = {{ 0 }};
    double diag_vals[2];
    diag_vals[0] = pow(R[0], exponent * -2.0);
    diag_vals[1] = pow(R[1], exponent * -2.0);
    set_diagonal(&O[0][0], diag_vals, 2, 2);

    // characteristic ellipse matrix
    double X_temp[2][2] = {{ 0 }};
    matrix_multiply(&X_temp[0][0], &O[0][0], &Q[0][0], 2, 2, 2);
    matrix_multiply(X, &Qt[0][0], &X_temp[0][0], 2, 2, 2);

}

double ellipse_overlap(double *rA, double *radiiA, double phiA, double *rB, double *radiiB, double phiB)
{

    // find XA^(-1) and XB^(1/2)
    double XA[2][2] = {{ 0 }};
    double XB[2][2] = {{ 0 }};
    characteristic_ellipse_matrix(&XA[0][0], &radiiA[0], phiA, -1.0);
    characteristic_ellipse_matrix(&XB[0][0], &radiiB[0], phiB, 0.5);

    // find A_AB
    double A_AB[2][2] = {{ 0 }};
    double A_temp[2][2] = {{ 0 }};
    matrix_multiply(&A_temp[0][0], &XA[0][0], &XB[0][0], 2, 2, 2);
    matrix_multiply(&A_AB[0][0], &XB[0][0], &A_temp[0][0], 2, 2, 2);

    // find r_AB
    double rAB[2];
    rAB[0] = rB[0] - rA[0];
    rAB[1] = rB[1] - rA[1];

    // find a_AB
    double a_AB[2];
    matrix_multiply(&a_AB[0], &XB[0][0], &rAB[0], 2, 2, 1);

    // extract elements of the matrix A_AB and a_AB
    double a11, a12, a21, a22;
    double b1, b2;
    a11 = A_AB[0][0];
    a12 = A_AB[0][1];
    a21 = A_AB[1][0];
    a22 = A_AB[1][1];
    b1 = a_AB[0];
    b2 = a_AB[1];

    // find coefficients for the parametric polynomial derivative used to find max
    double h[5];
    double z[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
    size_t n = 5;
    h[0] = h0(a11, a12, a21, a22, b1, b2);
    h[1] = h1(a11, a12, a21, a22, b1, b2);
    h[2] = h2(a11, a12, a21, a22, b1, b2);
    h[3] = h3(a11, a12, a21, a22, b1, b2);
    h[4] = h4(a11, a12, a21, a22, b1, b2);

    // find roots
    size_t m;
    m = find_roots(&z[0], &h[0], n);

    double F;
    int i;

    if (m > 1)
    {
        if (f(0, a11, a12, a21, a22, b1, b2) > f(1, a11, a12, a21, a22, b1, b2))
        {

           F = f(0, a11, a12, a21, a22, b1, b2);

        }
        else
        {

           F = f(1, a11, a12, a21, a22, b1, b2);

        }

        for(i=0; i<=m-2; i++)
        {

           if (( z[2 * i + 1] == 0 ) & (z[2 * i] > 0) & (z[2 * i] < 1))
           {

              F = f(z[2 * i], a11, a12, a21, a22, b1, b2);

           }

        }
    }
    else
    {
        F = 0.;
    }

    return F;

}

double container_square_overlap_potential(double *rA, double *radiiA, double phiA)
{

    double rB[2];
    double radiiB[2];
    double phiB;

    double top = 0;
    double bottom = 0;
    double left = 0;
    double right = 0;

    // top
    rB[0] = 0.5;
    rB[1] = 2.0;
    radiiB[0] = INFINITY;
    radiiB[1] = 1.0;
    phiB = 0.0;
    top = ellipse_overlap(&rA[0], &radiiA[0], phiA, &rB[0], &radiiB[0], phiB);

    // bottom
    rB[0] = 0.5;
    rB[1] = -1.0;
    radiiB[0] = INFINITY;
    radiiB[1] = 1.0;
    phiB = 0.0;
    bottom = ellipse_overlap(&rA[0], &radiiA[0], phiA, &rB[0], &radiiB[0], phiB);

    // left
    rB[0] = -1.0;
    rB[1] = 0.5;
    radiiB[0] = 1.0;
    radiiB[1] = INFINITY;
    phiB = 0.0;
    left = ellipse_overlap(&rA[0], &radiiA[0], phiA, &rB[0], &radiiB[0], phiB);

    // right
    rB[0] = 2.0;
    rB[1] = 0.5;
    radiiB[0] = 1.0;
    radiiB[1] = INFINITY;
    phiB = 0.0;
    right = ellipse_overlap(&rA[0], &radiiA[0], phiA, &rB[0], &radiiB[0], phiB);

    return fminf(top, fminf(bottom, fminf(left, right)));

}

size_t rsa_align_square(double *x, double *y,
    size_t npoints, double *radius, double phi, int step_limit,
    unsigned long randSeed)
{

    // Setup GSL random number generator
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    // Set the seed
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    // arrays for overlap_potential functions
    double rA[2];
    double rB[2];
    double radiiA[2];
    double radiiB[2];
    double F;
    double xn = 0.;
    double yn = 0.;

    radiiA[0] = radius[0];
    radiiA[1] = radius[1];

    // Set the initial position
    double C;
    C = fmaxf(radius[0], radius[1]);

    F = 0;
    while (F < 1.)
    {

        xn = gsl_rng_uniform (r);
        yn = gsl_rng_uniform (r);

        rA[0] = xn;
        rA[1] = yn;

        F = container_square_overlap_potential(&rA[0], &radiiA[0], phi);

    }

    x[0] = xn;
    y[0] = yn;

    size_t valid_pts;
    int k, flag, step;

    step = 0;
    valid_pts = 1;
    while ((valid_pts < npoints) & (step < step_limit))
    {

        // Generate new ellipse inside the container
        F = 0;
        while (F < 1.)
        {

            xn = gsl_rng_uniform (r);
            yn = gsl_rng_uniform (r);

            rA[0] = xn;
            rA[1] = yn;

            F = container_square_overlap_potential(&rA[0], &radiiA[0], phi);

        }

        // Determine if new ellipse overlaps with existing ellipses
        flag = 1;
        for (k = 0; k < valid_pts; k++)
        {

            rA[0] = x[k];
            rA[1] = y[k];
            rB[0] = xn;
            rB[1] = yn;
            radiiA[0] = radius[0];
            radiiA[1] = radius[1];
            radiiB[0] = radius[0];
            radiiB[1] = radius[1];

            F = ellipse_overlap(&rA[0], &radiiA[0], phi, &rB[0], &radiiB[0], phi);

            if (F < 1.0)
            {

                flag = 0;
                break;

            }
        }
        if (flag == 1)
        {

           x[valid_pts] = xn;
           y[valid_pts] = yn;
           valid_pts += 1;

        }

        step += 1;

    }

    gsl_rng_free (r);

    return valid_pts;

}

size_t rsa_square(double *x, double *y,
    size_t npoints, double *radius, double *phi, int step_limit,
    unsigned long randSeed)
{

    // Setup GSL random number generator
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    // Set the seed
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    // arrays for overlap_potential functions
    double rA[2];
    double rB[2];
    double radiiA[2];
    double radiiB[2];
    double phiA;
    double phiB;
    double F;
    double xn = 0;
    double yn = 0;
    double phin = 0;

    radiiA[0] = radius[0];
    radiiA[1] = radius[1];

    // Set the initial position
    double C;
    C = fmaxf(radius[0], radius[1]);

    F = 0;
    while (F < 1.)
    {

        xn = gsl_rng_uniform (r);
        yn = gsl_rng_uniform (r);

        rA[0] = xn;
        rA[1] = yn;
        phiA = 2 * M_PI * gsl_rng_uniform (r);

        F = container_square_overlap_potential(&rA[0], &radiiA[0], phiA);

    }

    x[0] = xn;
    y[0] = yn;
    phi[0] = phiA;

    size_t valid_pts;
    int k, flag, step;

    step = 0;
    valid_pts = 1;
    while ((valid_pts < npoints) & (step < step_limit))
    {

        // Generate new ellipse inside the container
        F = 0;
        while (F < 1.)
        {

            xn = gsl_rng_uniform (r);
            yn = gsl_rng_uniform (r);
            phin = 2 * M_PI * gsl_rng_uniform (r);

            rA[0] = xn;
            rA[1] = yn;
            phiA = phin;

            F = container_square_overlap_potential(&rA[0], &radiiA[0], phiA);

        }

        // Determine if new ellipse overlaps with existing ellipses
        flag = 1;
        for (k = 0; k < valid_pts; k++)
        {

            rA[0] = x[k];
            rA[1] = y[k];
            phiA = phi[k];

            rB[0] = xn;
            rB[1] = yn;
            phiB = phin;

            radiiA[0] = radius[0];
            radiiA[1] = radius[1];
            radiiB[0] = radius[0];
            radiiB[1] = radius[1];

            F = ellipse_overlap(&rA[0], &radiiA[0], phiA, &rB[0], &radiiB[0], phiB);

            if (F < 1.0)
            {

                flag = 0;
                break;

            }
        }
        if (flag == 1)
        {

           x[valid_pts] = xn;
           y[valid_pts] = yn;
           phi[valid_pts] = phin;
           valid_pts += 1;

        }

        step += 1;

    }

    gsl_rng_free (r);

    return valid_pts;

}


text:
/* hmm class header */
// Andrew Kern 8/24/05
//

#include "adkGSL.h"
#include <math.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>

typedef struct {
  int nstates;
  gsl_matrix *transition_matrix;
  gsl_matrix *transition_score_matrix; /* entries are logs of entries
                                        * in transition matrix */
  gsl_vector *piStart, *piEnd, *piStart_scores, *piEnd_scores;
  gsl_vector *otherData; //this is meant to hold distances between obs
} HMM;

HMM* newHMM(gsl_matrix *mat, gsl_vector *piStart, gsl_vector *piEnd, gsl_vector *otherData);
void freeHMM(HMM *hmm);
void printHMM(HMM *hmm);
void initHMM(HMM *hmm);
double forwardSumPrevious(HMM *hmm, gsl_matrix *alphas, gsl_matrix *transition_score_matrix_i, int state, int obsIndex);
double viterbiMaxPrevious(HMM *hmm, gsl_matrix *deltas, gsl_matrix *transition_score_matrix_i, int state, int obsIndex, int *point);
double backwardSumNext(HMM *hmm, gsl_matrix *alphas,  gsl_matrix *transition_score_matrix_i, gsl_matrix *emissions, int state, int obsIndex);
double forwardAlg(HMM *hmm, gsl_matrix *alphas, gsl_matrix *emissions, int length, int transitionFlag);
double backwardAlg(HMM *hmm, gsl_matrix *betas, gsl_matrix *emissions, int length,int transitionPowerFlag);
void posteriorProbs(HMM *hmm, gsl_matrix *postProbMat, gsl_matrix *emissions, int length, int transitionPowerFlag);
void posteriorProbsReduced(HMM *hmm, gsl_matrix *postProbMat, double pFor, gsl_matrix *alphas, gsl_matrix *betas, gsl_matrix *emissions, int length, int transitionPowerFlag);
void reestimateTransitions(HMM *hmm, gsl_matrix *postProbMat, gsl_matrix **jpMatArray, double pFor, gsl_matrix *alphas, \
			   gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag);
void reestimateTransitions_vanilla(HMM *hmm, double pFor, gsl_matrix *alphas, \
	gsl_matrix *betas, gsl_matrix *emisLogs, int length);

double viterbi(HMM *hmm, gsl_vector *viterbiPath, gsl_matrix *emissions,  int length, int transitionPowerFlag);
void hmm_logify_transitions(HMM *hmm);
void hmm_normalize_transitions(HMM *hmm);
double forwardAlgOffset(HMM *hmm, gsl_matrix *alphas, gsl_matrix *emissions, int length, int offset, int transitionPowerFlag);
double hmm_subpath_score(HMM *hmm, gsl_matrix *emisLogs, gsl_vector *stateInclusion, int begin, int length,  int transitionPowerFlag);
double hmm_subpath_log_odds(HMM *hmm, gsl_matrix *emisLogs, gsl_vector *stateInclusion1,gsl_vector *stateInclusion2,int begin, int length, int powerTransitionFlag);
gsl_vector *simulatePath(HMM *h, void *r, int obsNumber);
gsl_vector *simulatePathSpaced(HMM *h, void *r, int obsNumber, gsl_vector *locs, double mean);
void posteriorDecode(HMM *hmm, gsl_vector *posteriorPath, gsl_matrix *posts, int length);
void printTransitions(HMM *h, char *outfileName);

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
#ifndef XYCO_RUNTIME_ASYNC_H_
#define XYCO_RUNTIME_ASYNC_H_

#include <gsl/pointers>

#include "future.h"
#include "poll.h"
#include "runtime.h"

namespace xyco::runtime {
template <typename Return>
class AsyncFuture : public Future<Return> {
 public:
  [[nodiscard]] auto poll(Handle<void> self) -> Poll<Return> override {
    if (!ready_) {
      ready_ = true;
      auto res = RuntimeCtx::get_ctx()->blocking_handle()->Register(event_);
      return Pending();
    }

    gsl::owner<Return *> ret = gsl::owner<Return *>(
        std::get<AsyncFutureExtra>(event_.extra_).after_extra_);
    auto result = std::move(*ret);
    delete ret;

    return Ready<Return>{std::move(result)};
  }

  template <typename Fn>
  explicit AsyncFuture(Fn &&f) requires(std::is_invocable_r_v<Return, Fn>)
      : Future<Return>(nullptr),
        f_([&]() {
          auto extra = AsyncFutureExtra{.after_extra_ = new Return(f())};
          event_.extra_ = extra;
        }),
        ready_(false),
        event_(xyco::runtime::Event{
            .future_ = this, .extra_ = AsyncFutureExtra{.before_extra_ = f_}}) {
  }

  AsyncFuture(const AsyncFuture<Return> &) = delete;

  AsyncFuture(AsyncFuture<Return> &&) = delete;

  auto operator=(const AsyncFuture<Return> &) -> AsyncFuture<Return> & = delete;

  auto operator=(AsyncFuture<Return> &&) -> AsyncFuture<Return> & = delete;

  ~AsyncFuture() = default;

 private:
  bool ready_;
  std::function<void()> f_;
  Event event_;
};
}  // namespace xyco::runtime

#endif  // XYCO_RUNTIME_ASYNC_H_


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_multifit.h>

//DEFINES
#define SUCCESS 1
#define FAILURE -1
#define PROGRESS "===================="
#define PRG_WIDTH 20

//FUNCTIONS DECLARATION
int dcca_comp(double *, double *, double *, int, int, int, int, double *, int);
int rows_number(char *);
double mean(double *, int, int);
void cumsum(double *, double *, int);
void slice_vec(double *, double *, int, int);
void polynomialFit(int, int, double *, double *, double *);

//MAIN
int main(int argc, char **argv)
{
    //inputs
    int input_args = 7;
    if(argc < input_args){
        printf("Not enough input arguments!\n");
        return FAILURE;
    }else if(argc > input_args){
        printf("Too many input arguments!\n");
        return FAILURE;
    }
    char file_name1[255]; //path+name for the first file
    memset(file_name1, 0x00, sizeof(file_name1));
    sprintf(file_name1, "%s", argv[1]);
    char file_name2[255]; //path+name for the second file
    memset(file_name2, 0x00, sizeof(file_name2));
    sprintf(file_name2, "%s", argv[2]);
    int smallest_win_size = atoi(argv[3]); //minimum size for the sliding window
    int biggest_win_num = atoi(argv[4]); //number of windows in the last step
    int pol_ord = atoi(argv[5]); //polynomial order for fit
    char path_out[255]; //path+name for the output file
    memset(path_out, 0x00, sizeof(path_out));
    sprintf(path_out, "%s/rho_mw=%d_Mw=%d_pol=%d.txt", argv[6], smallest_win_size, biggest_win_num, pol_ord);
    //files length
    int N1 = rows_number(file_name1);
    if(N1 == FAILURE){
        printf("Cannot open input file (1).\n");
        return FAILURE;
    }
    int N2 = rows_number(file_name2);
    if(N2 == FAILURE){
        printf("Cannot open input file (2).\n");
        return FAILURE;
    }
    //time series
    double *ts1, *ts2;
    ts1 = calloc(N1, sizeof(double));
    if(!ts1){
        printf("MALLOC ERROR (ts1)\n");
        return FAILURE;
    }
    ts2 = calloc(N2, sizeof(double));
    if(!ts2){
        printf("MALLOC ERROR (ts2)\n");
        return FAILURE;
    }
    FILE *f1, *f2;
    f1 = fopen(file_name1, "r");
    if(!f1){
        printf("Cannot open input file (1).\n");
        return FAILURE;
    }
    for(int i = 0; i < N1; i++)
        fscanf(f1, "%lf", ts1+i);
    fclose(f1);
    f2 = fopen(file_name2, "r");
    if(!f2){
        printf("Cannot open input file (2).\n");
        return FAILURE;
    }
    for(int i = 0; i < N2; i++)
        fscanf(f2, "%lf", ts2+i);
    fclose(f2);
    //time series must have the same length
    int L = N1;
    if(N1 < N2){
        printf("Series 2 is longer than series 1. Data length reduced to: <%d>\n", L);
        slice_vec(ts2, ts2, 0, N1-1);
    }else if(N1 > N2){
        slice_vec(ts1, ts1, 0, N2-1);
        L = N2;
        printf("Series 1 is longer than series 2. Data length reduced to: <%d>\n", L);
    }else{
        printf("Data length: <%d>\n", L);
    }
    //time vector
    double *time;
    time = calloc(L, sizeof(double));
    if(!time){
        printf("MALLOC ERROR (time)\n");
        return FAILURE;
    }
    for(int i = 0; i < L; i++)
        time[i] = (double)(i+1);
    //time series minus its mean
    double ts1_ave = mean(ts1, L, L);
    double ts2_ave = mean(ts2, L, L);
    double *ts1_nomean, *ts2_nomean;
    ts1_nomean = calloc(L, sizeof(double));
    if(!ts1_nomean){
        printf("MALLOC ERROR (ts1_nomean)\n");
        return FAILURE;
    }
    ts2_nomean = calloc(L, sizeof(double));
    if(!ts2_nomean){
        printf("MALLOC ERROR (ts2_nomean)\n");
        return FAILURE;
    }
    for(int i = 0; i < L; i++){
        ts1_nomean[i] = ts1[i] - ts1_ave;
        ts2_nomean[i] = ts2[i] - ts2_ave;
    }
    //cumulative sum
    double *ts1_cum, *ts2_cum;
    ts1_cum = calloc(L, sizeof(double));
    if(!ts1_cum){
        printf("MALLOC ERROR (ts1_cum)\n");
        return FAILURE;
    }
    ts2_cum = calloc(L, sizeof(double));
    if(!ts2_cum){
        printf("MALLOC ERROR (ts2_cum)\n");
        return FAILURE;
    }
    printf("Data integration...");
    cumsum(ts1_nomean, ts1_cum, L);
    cumsum(ts2_nomean, ts2_cum, L);
    printf("done!\n");
    //dcca computation
    int biggest_win_size = L / biggest_win_num;
    int win_range = biggest_win_size - smallest_win_size + 1;
    double *F_DCCA, *F_DFA_x, *F_DFA_y, *rho;
    F_DCCA = calloc(win_range, sizeof(double));
    F_DFA_x = calloc(win_range, sizeof(double));
    F_DFA_y = calloc(win_range, sizeof(double));
    rho = calloc(win_range, sizeof(double));
    printf("Computing DCCA between series 1 and series 2...\n");
    int ret = dcca_comp(time, ts1_cum, ts2_cum, L, smallest_win_size, biggest_win_size, pol_ord, F_DCCA, win_range);
    if(ret == FAILURE)
        return FAILURE;
    printf("Computing DCCA between series 1 and series 1...\n");
    ret = dcca_comp(time, ts1_cum, ts1_cum, L, smallest_win_size, biggest_win_size, pol_ord, F_DFA_x, win_range);
    if(ret == FAILURE)
        return FAILURE;
    printf("Computing DCCA between series 2 and series 2...\n");
    ret = dcca_comp(time, ts2_cum, ts2_cum, L, smallest_win_size, biggest_win_size, pol_ord, F_DFA_y, win_range);
    if(ret == FAILURE)
        return FAILURE;
    FILE *fOut;
    fOut = fopen(path_out, "w");
    if(!fOut){
        printf("Cannot open output file.\n");
        return FAILURE;
    }
    printf("Computing rho...");
    for(int i = 0; i < win_range; i++){
        F_DFA_x[i] = sqrt(F_DFA_x[i]);
        F_DFA_y[i] = sqrt(F_DFA_y[i]);
        rho[i] = F_DCCA[i] / (double)(F_DFA_x[i] * F_DFA_y[i]);
        fprintf(fOut, "%lf %lf %lf %lf\n", F_DCCA[i], F_DFA_x[i], F_DFA_y[i], rho[i]);
    }
    printf("done\n");
    printf("Output file...");
    fclose(fOut);
    printf("done\n");
    //free memory
    free(ts1); free(ts2); free(time);
    free(ts1_nomean); free(ts2_nomean);
    free(ts1_cum); free(ts2_cum);
    free(F_DCCA); free(F_DFA_x); free(F_DFA_y); free(rho);
    return 0;
}

//FUNCTIONS
int dcca_comp(double *t, double *y1, double *y2, int N, int min_win_size, int last_win_len, int ord, double *F, int range_dcca)
{
    //fluctuations vector and other arrays
    int F_len = N - min_win_size;
    double *F_nu, *t_fit, *y_fit1, *y_fit2, *diff_vec, *fit_coeffs1, *fit_coeffs2;
    F_nu = calloc(F_len, sizeof(double));
    if(!F_nu){
        printf("MALLOC ERROR (F_nu)\n");
        return FAILURE;
    }
    t_fit = calloc((last_win_len+1), sizeof(double));
    if(!t_fit){
        printf("MALLOC ERROR (t_fit)\n");
        return FAILURE;
    }
    y_fit1 = calloc((last_win_len+1), sizeof(double));
    if(!y_fit1){
        printf("MALLOC ERROR (y_fit1)\n");
        return FAILURE;
    }
    y_fit2 = calloc((last_win_len+1), sizeof(double));
    if(!y_fit2){
        printf("MALLOC ERROR (y_fit2)\n");
        return FAILURE;
    }
    diff_vec = calloc((last_win_len+1), sizeof(double));
    if(!diff_vec){
        printf("MALLOC ERROR (diff_vec)\n");
        return FAILURE;
    }
    fit_coeffs1 = calloc(ord+1, sizeof(double));
    if(!fit_coeffs1){
        printf("MALLOC ERROR (fit_coeffs1)\n");
        return FAILURE;
    }
    fit_coeffs2 = calloc(ord+1, sizeof(double));
    if(!fit_coeffs2){
        printf("MALLOC ERROR (fit_coeffs2)\n");
        return FAILURE;
    }
    //computation
    int s, N_s, start_lim, end_lim;
    for(int i = 0; i < range_dcca; i++){
        double perc = i * 100 / (double)range_dcca;
        int prg = (i * PRG_WIDTH) / range_dcca;
        printf(" [%.*s%*s] %.2lf%%\r", prg, PROGRESS, PRG_WIDTH-prg, "", perc);
        fflush(stdout);
        s = i + min_win_size;
        N_s = N - s;
        for(int v = 0; v < N_s; v++){
            start_lim = v;
            end_lim = v + s;
            slice_vec(t, t_fit, start_lim, end_lim);
            slice_vec(y1, y_fit1, start_lim, end_lim);
            slice_vec(y2, y_fit2, start_lim, end_lim);
            polynomialFit(s+1, ord+1, t_fit, y_fit1, fit_coeffs1);
            polynomialFit(s+1, ord+1, t_fit, y_fit2, fit_coeffs2);
            for(int j = 0; j <= s; j++){
                for(int k = 0; k < ord+1; k++){
                    y_fit1[j] -= fit_coeffs1[k] * pow(t_fit[j], k);
                    y_fit2[j] -= fit_coeffs2[k] * pow(t_fit[j], k);
                }
                diff_vec[j] = y_fit1[j] * y_fit2[j];
            }
            F_nu[v] = mean(diff_vec, s+1, s-1);
        }
        F[i] = mean(F_nu, N_s, N_s);
    }
    printf(" [%s] 100.00%%\r", PROGRESS);
    fflush(stdout);
    printf("\n");
    free(F_nu); free(t_fit); free(y_fit1); free(y_fit2);
    free(diff_vec); free(fit_coeffs1); free(fit_coeffs2);
    return SUCCESS;
}

int rows_number(char *file_name)
{
    FILE *f;
    int stop;
    int lines = 0;
    f = fopen(file_name,"r");
    if(!f){
        printf("Cannot open file %s\n", file_name);
        return FAILURE;
    }
    while(!feof(f)){
        stop = fgetc(f);
        if(stop == '\n')
            lines++;
    }
    fclose(f);
    return lines;
}

double mean(double *vec, int vecL, int L)
{
    double avg = 0.0;
    for(int i = 0; i < vecL; i++)
        avg += vec[i];
    avg /= (double)L;
    return avg;
}

void cumsum(double *vec, double *sum_vec, int L)
{
    sum_vec[0] = vec[0];
    for(int i = 1; i < L; i++)
        sum_vec[i] = sum_vec[i-1] + vec[i];
}

void slice_vec(double *all_vec, double *sliced_vec, int start, int end)
{
    for(int i = 0; i <= (end-start); i++)
        sliced_vec[i] = all_vec[start+i];
}

void polynomialFit(int obs, int degree, double *dx, double *dy, double *store)
{
    gsl_multifit_linear_workspace *ws;
    gsl_matrix *cov, *X;
    gsl_vector *y, *c;
    double chisq;
    int i, j;
    //alloc
    X = gsl_matrix_alloc(obs, degree);
    y = gsl_vector_alloc(obs);
    c = gsl_vector_alloc(degree);
    cov = gsl_matrix_alloc(degree, degree);
    //computation
    for(i = 0; i < obs; i++){
        for(j = 0; j < degree; j++)
            gsl_matrix_set(X, i, j, pow(dx[i], j));
        gsl_vector_set(y, i, dy[i]);
    }
    ws = gsl_multifit_linear_alloc(obs, degree);
    gsl_multifit_linear(X, y, c, cov, &chisq, ws);
    for(i = 0; i < degree; i++)
        store[i] = gsl_vector_get(c, i);
    //free memory
    gsl_multifit_linear_free(ws);
    gsl_matrix_free(X);
    gsl_matrix_free(cov);
    gsl_vector_free(y);
    gsl_vector_free(c);
}