text:
/*
Copyright (c) 2003-2008 Rudi Cilibrasi, Rulers of the RHouse
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
3. Neither the name of the University nor the names of its contributors
   may be used to endorse or promote products derived from this software
   without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE RULERS AND CONTRIBUTORS ``AS IS'' AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED.  IN NO EVENT SHALL THE RULERS OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH DAMAGE.
*/
#include <glib.h>
#include <libintl.h>
#include <locale.h>
#include "complearn/newcomplearn.h"
#include "complearn/clconfig.h"

#define _(O) gettext(O)

#define clStringstackPush(a,v) {void*_u=(void*)v;g_array_append(a,_u);}while(0)

GString *complearn_matrix_prettyprint_text(LabeledMatrix *result) {
  int i, j;
  int html = complearn_ncd_get_html_output(complearn_ncd_top());
  char *rowstart = "", *rowend = "\n", *numstart = "", *numend = "", *labstart = "", *labend = " ", *lastcellstart="", *lastcellend ="";
  GString *toWrite = g_string_new("");
  CompLearnNcd *top = complearn_ncd_top();
  gsl_matrix *goodmat =gsl_matrix_alloc(result->mat->size1, result->mat->size2);
  gsl_matrix_memcpy(goodmat, result->mat);
  if (html) {
    g_string_append(toWrite, "
\n");
      rowstart = "
\n"; rowend = "
\n";
      lastcellstart = ""; lastcellend = "	";
      numstart = ""; numend = "	";
      labstart = ""; labend = "	";
    }
    for (j = 0; j < goodmat->size2; j += 1) {
      g_string_append(toWrite, rowstart);
      if (complearn_ncd_get_show_labels(top)) {
        g_string_append(toWrite, labstart);
        if (result->labels2[j] == NULL)
          g_error("bad labels2 %d with mat size2 %d", j, goodmat->size2);
        g_string_append(toWrite, result->labels2[j]);
        g_string_append(toWrite, labend);
      }
      for (i = 0; i < goodmat->size1; i += 1) {
        char buf[256];
        sprintf(buf, "%f ", gsl_matrix_get(goodmat, i, j));
        g_string_append(toWrite, i == goodmat->size1-1 ? lastcellstart : numstart);
        g_string_append(toWrite, buf);
        g_string_append(toWrite, i == goodmat->size1-1 ? lastcellend : numend);
      }
      g_string_append(toWrite, rowend);
    }
  if (html)
    g_string_append(toWrite, "
\n");
  gsl_matrix_free(goodmat);
  return toWrite;
}

#include <complearn/complearn.h>
#include <sys/types.h>
#include <unistd.h>
#include <assert.h>
#include <string.h>

#include <libxml/parser.h>
#include <libxml/tree.h>
#include <libxml/xmlwriter.h>
#include <libxml/encoding.h>

#define MY_ENCODING "ISO-8859-1"

#include <gsl/gsl_linalg.h>

#define PWD_RDY 1
#define UTS_RDY 1

#ifndef __MINGW32__

#if PWD_RDY
#include <pwd.h>
#endif

#if UTS_RDY
#include <sys/utsname.h>
#endif

#endif

GArray *clDefaultLabels(int i);

static void complearn_handle_string_list(xmlDocPtr doc, xmlNodePtr node, GArray *ss,
                             const char *tagname) {
  node = node->xmlChildrenNode;
  while (node != NULL) {
    if (strcmp(tagname, (char *) node->name) == 0) {
      guchar *v = xmlNodeListGetString(doc, node->xmlChildrenNode,1);
      g_array_append_val(ss, v);
    }
    node = node->next;
  }
}

static void handleDM(xmlDocPtr doc, xmlNodePtr node, GArray *ss,
struct CLDistMatrix *result) {
  result->title = (char *) xmlGetProp(node, (unsigned char *) "title");
  result->creationTime = (char *) xmlGetProp(node, (unsigned char *) "creationtime");
  node = node->xmlChildrenNode;
  GArray *numa, *lab1a, *lab2a;
  numa = result->numa;
  lab1a = result->labels1a;
  lab2a = result->labels2a;
  while (node != NULL) {
    do {
        if (strcmp("entries", (char *) node->name) == 0) {
          complearn_handle_string_list(doc, node, numa, "number");
          continue;
        }
        if (strcmp((char *) node->name, "axis1") == 0) {
          complearn_handle_string_list(doc, node, lab1a, "name");
          continue;
        }
        if (strcmp((char *) node->name, "axis2") == 0) {
          complearn_handle_string_list(doc, node, lab2a, "name");
          continue;
        }
      } while(0);
    node = node->next;
  }
}

static void complearn_write_string_list(xmlTextWriterPtr tw, char **ss, const char *topname, const char *itemname) {
  int rc;
  int i;
  rc = xmlTextWriterStartElement(tw, (unsigned char *) topname);
  for (i = 0; i < complearn_count_strings((const char * const * )ss); i += 1)
    rc = xmlTextWriterWriteElement(tw, (unsigned char *) itemname, (unsigned char *) ss[i]);
  rc = xmlTextWriterEndElement(tw);
}

static char *getTimestrNow(void)
{
  static char buf[32];
  sprintf(buf, "%lu", (unsigned long) time(NULL));
  return buf;
}

static struct CLDistMatrix *fill_in_blanks(struct CLDistMatrix *clb)
{
  if (clb->m == NULL) {
    g_error(_("Cannot write NULL gsl_matrix, exitting."));
  }
  if (clb->m->mat->size1 < 1 || clb->m->mat->size2 < 1) {
    g_error(_("Invalid gsl_matrix size, cannot write.  Exitting."));
  }
  if (clb->fileverstr == NULL)
    clb->fileverstr = g_strdup("1.0");
  if (clb->cllibver == NULL)
    clb->cllibver = g_strdup(complearn_package_version);
  if (clb->username == NULL)
    clb->username = g_strdup(g_get_user_name());
  if (clb->hostname == NULL)
    clb->hostname = g_strdup(complearn_get_hostname());
  if (clb->title == NULL)
    clb->title = g_strdup(_("untitled"));
  if (clb->compressor == NULL)
    clb->compressor = g_strdup(_("unknown"));
  if (clb->creationTime == NULL)
    clb->creationTime = g_strdup(getTimestrNow());
  return clb;
}

struct CLDistMatrix *complearn_clone_cldm(struct CLDistMatrix *clb)
{
  struct CLDistMatrix *r = calloc(sizeof(*r), 1);
  if (clb->m == NULL) {
    g_error(_("Cannot write NULL gsl_matrix, exitting."));
  }
#define DUPSTR(s) if (clb->s) r->s = g_strdup(clb->s)
  DUPSTR(cllibver);
  DUPSTR(username);
  DUPSTR(hostname);
  DUPSTR(title);
  DUPSTR(compressor);
  DUPSTR(creationTime);
  r->m = calloc(sizeof(*(r->m)), 1);
  r->m->mat = gsl_matrix_alloc(clb->m->mat->size1, clb->m->mat->size2);
  gsl_matrix_memcpy(r->m->mat, clb->m->mat);
  r->m->labels1 = complearn_dupe_strings((const char * const *)clb->m->labels1);
  r->m->labels2 = complearn_dupe_strings((const char * const *)clb->m->labels2);
  r->cmds = complearn_dupe_strings((const char * const *)clb->cmds);
  r->cmdtimes = complearn_dupe_strings((const char * const *)clb->cmdtimes);
  return r;
}

void complearn_free_cldm(struct CLDistMatrix *clb)
{
  if (clb->cllibver) g_free(clb->cllibver);
  if (clb->username) g_free(clb->username);
  if (clb->hostname) g_free(clb->hostname);
  if (clb->title) g_free(clb->title);
  if (clb->compressor) g_free(clb->compressor);
  if (clb->creationTime) g_free(clb->creationTime);
  if (clb->m->mat) gsl_matrix_free(clb->m->mat);
  g_strfreev(clb->m->labels1);
  g_strfreev(clb->m->labels2);
  g_strfreev(clb->cmdtimes);
  g_strfreev(clb->cmds);
  g_free(clb->m);
  memset(clb, 0, sizeof(*clb));
  g_free(clb);
}

static char *clXMLQuoteStr(char *inp)
{
  int i, c, j;
  static char *outstr;
  if (outstr != NULL)
    g_free(outstr);
  outstr = calloc(strlen(inp) * 10, 1);
  j = 0;
  for (i = 0; inp[i]; i += 1) {
    c = inp[i];
    switch(c) {
      case '&':
        outstr[j++] = '&'; outstr[j++] = 'a'; outstr[j++] = 'm';
        outstr[j++] = 'p'; outstr[j++] = ';'; break;
      case '>':
        outstr[j++] = '&'; outstr[j++] = 'g'; outstr[j++] = 't';
        outstr[j++] = ';'; break;
      case '<':
        outstr[j++] = '&'; outstr[j++] = 'l'; outstr[j++] = 't';
        outstr[j++] = ';'; break;
      case '"':
        outstr[j++] = '&'; outstr[j++] = 'q'; outstr[j++] = 'o';
        outstr[j++] = 'o'; outstr[j++] = 't'; outstr[j++] = ';'; break;
      case '\'':
        outstr[j++] = '&'; outstr[j++] = 'a'; outstr[j++] = 'p';
        outstr[j++] = 'o'; outstr[j++] = 's'; outstr[j++] = ';'; break;
//      case '\\':
 //     case '#':
//        outstr[j++] = '\\'; outstr[j++] = c; break;
      default:
        outstr[j++] = c;
    }
  }
  outstr[j++] = 0;
  return outstr;
}

static GString *clRealWriteCLBDistMatrix(struct CLDistMatrix *clb)
{
  int rc;
  GString *result, *cres;
  int dim1, dim2, i, j;
  gsl_matrix *m = clb->m->mat;
  xmlBufferPtr b;
  xmlTextWriterPtr tw;
  clb = fill_in_blanks(complearn_clone_cldm(clb));
  b = xmlBufferCreate();
  tw =xmlNewTextWriterMemory(b, 0);
  rc = xmlTextWriterStartDocument(tw, NULL, MY_ENCODING, NULL);
  rc = xmlTextWriterStartElement(tw, (unsigned char *) "clb");
  rc = xmlTextWriterWriteAttribute(tw, (unsigned char *) "version", (unsigned char *) clb->fileverstr);

    rc = xmlTextWriterWriteElement(tw, (unsigned char *) "cllibver", (unsigned char *) clb->cllibver);
    rc = xmlTextWriterWriteElement(tw, (unsigned char *) "username", (unsigned char *) clb->username);
    rc = xmlTextWriterWriteElement(tw, (unsigned char *) "hostname", (unsigned char *) clb->hostname);
    rc = xmlTextWriterWriteElement(tw, (unsigned char *) "compressor", (unsigned char *) clb->compressor);
    rc = xmlTextWriterStartElement(tw,(unsigned char *)  "distmatrix");
    rc = xmlTextWriterWriteAttribute(tw,(unsigned char *)  "creationtime", (unsigned char *) clb->creationTime);
    rc = xmlTextWriterWriteAttribute(tw, (unsigned char *) "title",(unsigned char *)  clb->title);
    complearn_write_string_list(tw, clb->m->labels1,  "axis1",  "name");
    complearn_write_string_list(tw, clb->m->labels2,  "axis2",  "name");
    dim1 = complearn_count_strings((const char * const *) clb->m->labels1);
    dim2 = complearn_count_strings((const char * const *) clb->m->labels2);
      rc = xmlTextWriterStartElement(tw,(unsigned char *)  "entries");
      for (i = 0; i < dim1; i += 1) {
        for (j = 0; j < dim2; j += 1) {
          double g = gsl_matrix_get(m, i, j);
            xmlTextWriterWriteFormatElement(tw, (unsigned char *) "number", (char *) "%f", g);
        }
      }
      rc = xmlTextWriterEndElement(tw);
    rc = xmlTextWriterEndElement(tw);
    if (clb->cmds) {
      char *str;
      rc = xmlTextWriterStartElement(tw,(unsigned char *)  "commands");
      for (i = 0; i < complearn_count_strings((const char * const *) clb->cmds); i += 1) {
        char *t = clb->cmdtimes ? (clb->cmdtimes[i]) : NULL;
        rc = xmlTextWriterStartElement(tw,(unsigned char *)  "cmdstring");
        if (t)
          rc = xmlTextWriterWriteAttribute(tw,(unsigned char *)  "creationtime", (unsigned char *) t);
        str = clb->cmds[i];
        xmlTextWriterWriteRaw(tw, (unsigned char *)  clXMLQuoteStr(str));
        rc = xmlTextWriterEndElement(tw);
      }
      rc = xmlTextWriterEndElement(tw);
    }
  rc = xmlTextWriterEndElement(tw);
  rc = xmlTextWriterEndDocument(tw);
  xmlFreeTextWriter(tw);
  result = g_string_new((char *) b->content);
  xmlBufferFree(b);
  /* TODO: remove more mem leaks from above func */
  complearn_free_cldm(clb);
  CompLearnRealCompressor *bz = COMPLEARN_REAL_COMPRESSOR(complearn_environment_get_nameable("bzlib"));
  if (bz == NULL)
    g_error(_("Cannot load bzlib"));
  cres = real_compressor_compress(bz, result);
  g_string_free(result, TRUE);
  return cres;
}

gsl_matrix *complearn_clb_dist_matrix(const GString *db)
{
  return complearn_read_clb_dist_matrix(db)->m->mat;
}

struct CLDistMatrix *complearn_read_clb_dist_matrix(const GString *udb)
{
  struct CLDistMatrix *result = calloc(sizeof(struct CLDistMatrix), 1);
  GArray *ents = g_array_new(FALSE, TRUE, sizeof(gpointer));
  GString *db;
  int dim1, dim2;
  xmlDocPtr doc;
  xmlNodePtr node;
  CompLearnRealCompressor *bz = COMPLEARN_REAL_COMPRESSOR(complearn_environment_get_nameable("bzlib"));
  if (bz == NULL)
    g_error(_("Cannot load bzlib compressor."));
  if (real_compressor_is_decompressible(bz, udb))
    db = real_compressor_decompress(bz, udb);
  else
    db = g_string_new_len(udb->str, udb->len);
  if (db == NULL || db->len < 1)
    return NULL;
  if (db->str[0] != '<')
    return NULL;
  result->fileverstr = "";
  result->username = "";
  result->title = "";
  result->creationTime = "";
  result->cmdsa = g_array_new(FALSE, TRUE, sizeof(gpointer));
  result->cmdtimesa = g_array_new(FALSE, TRUE, sizeof(gpointer));
  result->labels1a = g_array_new(FALSE, TRUE, sizeof(gpointer));
  result->labels2a = g_array_new(FALSE, TRUE, sizeof(gpointer));
  result->numa = g_array_new(FALSE, TRUE, sizeof(gpointer));
  result->m = calloc(sizeof(*result->m), 1);
  doc = xmlReadMemory((char *)db->str, db->len, "noname.xml",
                      NULL, 0);
  if (doc == NULL) {
    g_error(_("Failed to parse document."));
    g_string_free(db, TRUE);
    free(result);
    return NULL;
  }
  node = doc->children;
  if (strcmp((char *) node->name, "clb") != 0) {
    free(result);
    g_string_free(db, TRUE);
    return NULL;
  }
  result->fileverstr = (char *) xmlGetProp(node, (unsigned char *) "version");
  node = node->xmlChildrenNode;
 while (node != NULL) {
    do {
      if (strcmp((char *) node->name, "commands") == 0) {
        complearn_handle_string_list(doc, node, result->cmdsa, "cmdstring");
        continue;
      }
      if (strcmp((char *) node->name, "distmatrix") == 0) {
        handleDM(doc, node, ents, result);
  result->m->mat = gsl_matrix_alloc(result->labels1a->len, result->labels2a->len);
  int i, x=0, y=0;
  for (i = 0; i < result->numa->len; i += 1) {
    gsl_matrix_set(result->m->mat,x,y,atof(g_array_index(result->numa, char *, i)));
    g_assert(y < result->m->mat->size2);
    x += 1;
    if (x == result->m->mat->size1) {
      x = 0; y += 1;
    }
  }
        continue;
      }
      if (strcmp((char *) node->name, "cllibver") == 0) {
        result->cllibver = (char *) g_strdup((char *) xmlNodeListGetString(doc, node->xmlChildrenNode, 1));
        continue;
      }
      if (strcmp((char *) node->name, "compressor") == 0) {
        result->compressor = (char *) g_strdup((char *) xmlNodeListGetString(doc, node->xmlChildrenNode, 1));
        continue;
      }
      if (strcmp((char *) node->name, "username") == 0) {
        result->username = (char *) g_strdup((char *) xmlNodeListGetString(doc, node->xmlChildrenNode, 1));
        continue;
      }
    } while (0);
    node = node->next;
    }

  dim1 = result->labels1a->len;
  if (dim1 <= 0) {
    g_error(_("Error, no labels for dimension 1."));
    exit(1);
  }
  dim2 = result->labels2a->len;
  if (dim2 <= 0) {
    g_error(_("Error, no labels for dimension 2."));
    exit(1);
  }
  g_array_free(ents, TRUE);
  g_string_free(db, TRUE);
  result->m->labels1 = calloc(result->labels1a->len+1, sizeof(gpointer));
  result->m->labels2 = calloc(result->labels2a->len+1, sizeof(gpointer));
  memcpy(result->m->labels1, result->labels1a->data, sizeof(gpointer)*(result->labels1a->len));
  memcpy(result->m->labels2, result->labels2a->data, sizeof(gpointer)*(result->labels2a->len));
  return result;
}

const char *complearn_complearn_get_hostname(void)
{
  static char hostname[1024];
#ifdef __MINGW32__
	strcpy(hostname, "localhost");
#else
  gethostname(hostname, 1024);
#endif
  return hostname;
}

const char *complearn_get_uts_name(void)
{
  static char utsname[1024];
#if UTS_RDY
#ifdef __MINGW32__
	strcpy(utsname, "mingw32");
#else
  struct utsname uts;
  uname(&uts);
  sprintf(utsname, "[%s (%s), %s]", uts.sysname, uts.release, uts.machine);
#endif
#else
  strcpy(utsname, "[unknown]");
#endif
  return utsname;
}

GString *complearn_matrix_prettyprint_clb(LabeledMatrix *dm) {
  struct CLDistMatrix inp;
  memset(&inp, 0, sizeof(inp));
  inp.m = dm;
  char *old_locale = setlocale(LC_NUMERIC, NULL);
  setlocale(LC_NUMERIC, "C");
  GString *result = clRealWriteCLBDistMatrix(&inp);
  setlocale(LC_NUMERIC, old_locale);
  return result;
}

LabeledMatrix *complearn_load_any_matrix(const GString *inp)
{
  complearn_environment_top();
  if (complearn_is_nexus_file(inp))
    return complearn_load_nexus_matrix(inp);
  struct CLDistMatrix *cld;
  cld = complearn_read_clb_dist_matrix(inp);
  if (cld) {
    return cld->m;
  }
  if (complearn_is_text_matrix(inp)) {
    return complearn_load_text_matrix(inp);
  }
  return NULL;
}

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
#ifndef SPC_UTILS_H
#define SPC_UTILS_H

void init_search (const char *string);  /* Pbmsrch.C      */
char *strsearch (const char *string);   /* Pbmsrch.C      */

#include <stddef.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <limits.h>
#include <math.h>

#include <gsl/gsl_version.h>
#include <gsl/gsl_nan.h>
#include <gsl/gsl_sys.h>
#include <gsl/gsl_vector.h>

#include "aXe_errors.h"
#include "aXe_grism.h"

#define LAMBDACENTRAL 8000.0
#define COLNAMELENGTH 12
#define MAXCOLS       200

/**
 * Structure: sexcol
 * Structure for a column header
 */
typedef struct sexcol
{
  char name[COLNAMELENGTH];
  int  number;
}
sexcol;

/**
 * Structure: colinfo
 * Structure for the table header
 */
typedef struct colinfo
{
  int numcols;
  sexcol columns[MAXCOLS];
}
colinfo;

extern int
is_valid_entry(double mag);

extern int
get_valid_entries(const gsl_vector *magnitudes);

extern int
check_worldcoo_input(const colinfo * actcatinfo, const int thsky);

extern int
check_imagecoo_input(const colinfo * actcatinfo);

extern void
make_GOL_header(FILE *fout, const colinfo * actcatinfo,
                const gsl_vector * waves, const gsl_vector * cnums,
                const px_point backwin_cols, const px_point modinfo_cols);

extern colinfo *
get_sex_col_descr (char *filename);

extern double
get_col_value (const colinfo * actcatinfo, const char key[],
               gsl_vector * v, int fatal);

extern double
get_col_value2 (const colinfo * actcatinfo, const char key[],
                gsl_vector * v, int fatal);

extern int
line_is_valid (const colinfo * actcatinfo, char line[]);

extern int
get_magauto_col(const gsl_vector *wavelength, const gsl_vector *colnums,
                const double lambda_mark);

extern int
has_magnitudes(const colinfo * actcatinfo);

extern px_point
has_backwindow(const colinfo * actcatinfo);

extern px_point
has_modelinfo(const colinfo * actcatinfo);

extern int
get_magcols(const colinfo * actcatinfo, gsl_vector *wavelength,
            gsl_vector *colnums);

extern int
resolve_colname(const char colname[]);

extern void
get_columname(const colinfo * actcatinfo, const int colnum, char colname[]);

extern int
get_columnumber(const char cname[], const colinfo * actcatinfo);

extern char *
rmlead (char *str);

extern char *
stptok (const char *s, char *tok, size_t toklen, char *brk);

#define      strMove(d,s) memmove(d,s,strlen(s)+1)

extern void
lv1ws (char *str);

extern int
isnum2(char *string);

extern gsl_vector *
string_to_gsl_array (char *str);

extern void
check_libraries (void);

extern observation *
load_dummy_observation (void);
#endif

ext: h

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#include "NativeDataStream.h"
#include "PerFrameValue.h"
#include "ShaderCompiler.h"
#include "VertexArray.h"

#include <Babylon/JsRuntime.h>
#include <Babylon/JsRuntimeScheduler.h>

#include <GraphicsImpl.h>
#include <BgfxCallback.h>
#include <FrameBuffer.h>

#include <napi/napi.h>

#include <bgfx/bgfx.h>
#include <bgfx/platform.h>
#include <bimg/bimg.h>
#include <bx/allocator.h>

#include <gsl/gsl>

#include <arcana/threading/cancellation.h>
#include <unordered_map>

namespace Babylon
{
    struct UniformInfo final
    {
        UniformInfo(uint8_t stage, bgfx::UniformHandle handle) :
            Stage{stage},
            Handle{handle}
        {
        }

        uint8_t Stage{};
        bgfx::UniformHandle Handle{bgfx::kInvalidHandle};
    };

    struct ProgramData final
    {
        ProgramData() = default;
        ProgramData(ProgramData&& other) = delete;
        ProgramData(const ProgramData&) = delete;
        ProgramData& operator=(ProgramData&& other) = delete;
        ProgramData& operator=(const ProgramData& other) = delete;

        ~ProgramData()
        {
            Dispose();
        }

        void Dispose()
        {
            if (!Disposed && bgfx::isValid(Handle))
            {
                bgfx::destroy(Handle);
            }
            Disposed = true;
        }

        std::unordered_map<std::string, uint32_t> VertexAttributeLocations{};
        std::unordered_map<std::string, UniformInfo> UniformInfos{};

        bgfx::ProgramHandle Handle{bgfx::kInvalidHandle};
        bool Disposed{false};

        struct UniformValue
        {
            std::vector<float> Data{};
            uint16_t ElementLength{};
        };

        std::unordered_map<uint16_t, UniformValue> Uniforms{};

        void SetUniform(bgfx::UniformHandle handle, gsl::span<const float> data, size_t elementLength = 1)
        {
            UniformValue& value = Uniforms[handle.idx];
            value.Data.assign(data.begin(), data.end());
            value.ElementLength = static_cast<uint16_t>(elementLength);
        }
    };

    class NativeEngine final : public Napi::ObjectWrap<NativeEngine>
    {
        static constexpr auto JS_CLASS_NAME = "_NativeEngine";
        static constexpr auto JS_CONSTRUCTOR_NAME = "Engine";

    public:
        NativeEngine(const Napi::CallbackInfo& info);
        NativeEngine(const Napi::CallbackInfo& info, JsRuntime& runtime);
        ~NativeEngine();

        static void Initialize(Napi::Env env);

    private:
        void Dispose();

        void Dispose(const Napi::CallbackInfo& info);
        void RequestAnimationFrame(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexArray(const Napi::CallbackInfo& info);
        void DeleteVertexArray(NativeDataStream::Reader& data);
        void BindVertexArray(NativeDataStream::Reader& data);
        Napi::Value CreateIndexBuffer(const Napi::CallbackInfo& info);
        void DeleteIndexBuffer(NativeDataStream::Reader& data);
        void RecordIndexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicIndexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexBuffer(const Napi::CallbackInfo& info);
        void DeleteVertexBuffer(NativeDataStream::Reader& data);
        void RecordVertexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicVertexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateProgram(const Napi::CallbackInfo& info);
        Napi::Value GetUniforms(const Napi::CallbackInfo& info);
        Napi::Value GetAttributes(const Napi::CallbackInfo& info);
        void SetProgram(NativeDataStream::Reader& data);
        void DeleteProgram(NativeDataStream::Reader& data);
        void SetState(NativeDataStream::Reader& data);
        void SetZOffset(NativeDataStream::Reader& data);
        void SetZOffsetUnits(NativeDataStream::Reader& data);
        void SetDepthTest(NativeDataStream::Reader& data);
        void SetDepthWrite(NativeDataStream::Reader& data);
        void SetColorWrite(NativeDataStream::Reader& data);
        void SetBlendMode(NativeDataStream::Reader& data);
        void SetMatrix(NativeDataStream::Reader& data);
        void SetInt(NativeDataStream::Reader& data);
        void SetIntArray(NativeDataStream::Reader& data);
        void SetIntArray2(NativeDataStream::Reader& data);
        void SetIntArray3(NativeDataStream::Reader& data);
        void SetIntArray4(NativeDataStream::Reader& data);
        void SetFloatArray(NativeDataStream::Reader& data);
        void SetFloatArray2(NativeDataStream::Reader& data);
        void SetFloatArray3(NativeDataStream::Reader& data);
        void SetFloatArray4(NativeDataStream::Reader& data);
        void SetMatrices(NativeDataStream::Reader& data);
        void SetMatrix3x3(NativeDataStream::Reader& data);
        void SetMatrix2x2(NativeDataStream::Reader& data);
        void SetFloat(NativeDataStream::Reader& data);
        void SetFloat2(NativeDataStream::Reader& data);
        void SetFloat3(NativeDataStream::Reader& data);
        void SetFloat4(NativeDataStream::Reader& data);
        Napi::Value CreateTexture(const Napi::CallbackInfo& info);
        void LoadTexture(const Napi::CallbackInfo& info);
        void CopyTexture(const Napi::CallbackInfo& info);
        void LoadRawTexture(const Napi::CallbackInfo& info);
        void LoadCubeTexture(const Napi::CallbackInfo& info);
        void LoadCubeTextureWithMips(const Napi::CallbackInfo& info);
        Napi::Value GetTextureWidth(const Napi::CallbackInfo& info);
        Napi::Value GetTextureHeight(const Napi::CallbackInfo& info);
        void SetTextureSampling(NativeDataStream::Reader& data);
        void SetTextureWrapMode(NativeDataStream::Reader& data);
        void SetTextureAnisotropicLevel(NativeDataStream::Reader& data);
        void SetTexture(NativeDataStream::Reader& data);
        void DeleteTexture(const Napi::CallbackInfo& info);
        Napi::Value CreateFrameBuffer(const Napi::CallbackInfo& info);
        void DeleteFrameBuffer(NativeDataStream::Reader& data);
        void BindFrameBuffer(NativeDataStream::Reader& data);
        void UnbindFrameBuffer(NativeDataStream::Reader& data);
        void DrawIndexed(NativeDataStream::Reader& data);
        void Draw(NativeDataStream::Reader& data);
        void Clear(NativeDataStream::Reader& data);
        Napi::Value GetRenderWidth(const Napi::CallbackInfo& info);
        Napi::Value GetRenderHeight(const Napi::CallbackInfo& info);
        void SetViewPort(const Napi::CallbackInfo& info);
        Napi::Value GetHardwareScalingLevel(const Napi::CallbackInfo& info);
        void SetHardwareScalingLevel(const Napi::CallbackInfo& info);
        Napi::Value CreateImageBitmap(const Napi::CallbackInfo& info);
        Napi::Value ResizeImageBitmap(const Napi::CallbackInfo& info);
        void GetFrameBufferData(const Napi::CallbackInfo& info);
        void SetStencil(NativeDataStream::Reader& data);
        void SetCommandDataStream(const Napi::CallbackInfo& info);
        void SubmitCommands(const Napi::CallbackInfo& info);
        void DrawInternal(bgfx::Encoder* encoder, uint32_t fillMode);

        std::string ProcessShaderCoordinates(const std::string& vertexSource);

        GraphicsImpl::UpdateToken& GetUpdateToken();
        FrameBuffer& GetBoundFrameBuffer(bgfx::Encoder& encoder);

        std::shared_ptr<arcana::cancellation_source> m_cancellationSource{};

        ShaderCompiler m_shaderCompiler{};

        ProgramData* m_currentProgram{nullptr};

        JsRuntime& m_runtime;
        GraphicsImpl& m_graphicsImpl;

        JsRuntimeScheduler m_runtimeScheduler;

        std::optional<GraphicsImpl::UpdateToken> m_updateToken{};

        void ScheduleRequestAnimationFrameCallbacks();
        bool m_requestAnimationFrameCallbacksScheduled{};

        bx::DefaultAllocator m_allocator{};
        uint64_t m_engineState{BGFX_STATE_DEFAULT};
        uint32_t m_stencilState{BGFX_STENCIL_TEST_ALWAYS | BGFX_STENCIL_FUNC_REF(0) | BGFX_STENCIL_FUNC_RMASK(0xFF) | BGFX_STENCIL_OP_FAIL_S_KEEP | BGFX_STENCIL_OP_FAIL_Z_KEEP | BGFX_STENCIL_OP_PASS_Z_REPLACE};

        template<int size, typename arrayType>
        void SetTypeArrayN(const UniformInfo& uniformInfo, const uint32_t elementLength, const arrayType& array);

        template<int size>
        void SetIntArrayN(NativeDataStream::Reader& data);

        template<int size>
        void SetFloatArrayN(NativeDataStream::Reader& data);

        template<int size>
        void SetFloatN(NativeDataStream::Reader& data);

        template<int size>
        void SetMatrixN(NativeDataStream::Reader& data);

        // Scratch vector used for data alignment.
        std::vector<float> m_scratch{};

        std::vector<Napi::FunctionReference> m_requestAnimationFrameCallbacks{};

        VertexArray* m_boundVertexArray{};
        FrameBuffer m_defaultFrameBuffer;
        FrameBuffer* m_boundFrameBuffer{};
        PerFrameValue<bool> m_boundFrameBufferNeedsRebinding;

        // TODO: This should be changed to a non-owning ref once multi-update is available.
        NativeDataStream* m_commandStream{};
    };
}

max_stars_repo_name: olidum/Babylon

max_issues_repo_name: olidum/Babylon

max_forks_repo_name: olidum/Babylon


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_core.h>
#include <gbpCosmo_NFW_etc.h>
#include <gsl/gsl_sf_expint.h>

double V2_circ(double M, double r) {
    return (G_NEWTON * M / r);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
// Copyright © Matt Jones and Contributors. Licensed under the MIT Licence (MIT). See LICENCE.md in the repository root
// for more information.

#ifndef NOVELRT_GRAPHICS_H
#define NOVELRT_GRAPHICS_H

// Graphics dependencies
#include "NovelRT/EngineConfig.h"
#include "NovelRT/Maths/Maths.h"
#include "NovelRT/ResourceManagement/ResourceManagement.h"
#include "NovelRT/SceneGraph/SceneGraph.h"
#include "NovelRT/Threading/Threading.h"
#include "NovelRT/Utilities/Event.h"
#include "NovelRT/Utilities/Lazy.h"
#include "NovelRT/Utilities/Misc.h"
#include "RGBAColour.h"
#include <chrono>
#include <cstdint>
#include <filesystem>
#include <gsl/span>
#include <list>
#include <memory>
#include <mutex>
#include <optional>
#include <string>
#include <typeindex>
#include <utility>
#include <vector>

/**
 * @brief The experimental Graphics plugin API. Comes with built-in support for the ECS.
 */
namespace NovelRT::Graphics
{
    enum class ShaderProgramKind : uint32_t;
    enum class GraphicsResourceAccess : uint32_t;
    enum class GraphicsSurfaceKind : uint32_t;
    enum class GraphicsTextureAddressMode : uint32_t;
    enum class GraphicsPipelineBlendFactor : uint32_t;
    enum class GraphicsPipelineInputElementKind : uint32_t;
    enum class GraphicsPipelineResourceKind : uint32_t;
    enum class ShaderProgramVisibility : uint32_t;
    enum class GraphicsTextureKind : uint32_t;
    enum class GraphicsMemoryRegionAllocationFlags : uint32_t;
    enum class GraphicsMemoryRegionAllocationFlags : uint32_t;
    enum class GraphicsBufferKind : uint32_t;
    enum class TexelFormat : uint32_t;
    struct GraphicsMemoryAllocatorSettings;
    class GraphicsDeviceObject;
    class IGraphicsSurface;
    class GraphicsAdapter;
    class GraphicsDevice;
    class GraphicsResource;
    class GraphicsBuffer;
    class GraphicsTexture;
    class ShaderProgram;
    class GraphicsPipeline;
    class GraphicsPipelineSignature;
    class GraphicsPipelineInput;
    class GraphicsPipelineResource;
    class GraphicsPipelineInputElement;
    class GraphicsContext;
    class GraphicsFence;
    class GraphicsPrimitive;
    class GraphicsProvider;
    class GraphicsMemoryAllocator;
    class IGraphicsAdapterSelector;
    class GraphicsMemoryBlockCollection;
    class GraphicsMemoryBlock;
    class GraphicsMemoryBudget;
    class GraphicsSurfaceContext;
    class GraphicsResourceManager;
}

// Graphics types
// clang-format off

#include "ShaderProgramKind.h"
#include "GraphicsAdapter.h"
#include "GraphicsDeviceObject.h"
#include "GraphicsContext.h"
#include "GraphicsFence.h"
#include "GraphicsMemoryAllocatorSettings.h"
#include "GraphicsMemoryRegion.h"
#include "IGraphicsMemoryRegionCollection.h"
#include "GraphicsMemoryRegionAllocationFlags.h"
#include "TexelFormat.h"
#include "GraphicsMemoryAllocator.h"
#include "GraphicsMemoryBlockCollection.h"
#include "GraphicsMemoryBudget.h"
#include "GraphicsMemoryBlock.h"
#include "GraphicsResourceAccess.h"
#include "GraphicsSurfaceKind.h"
#include "GraphicsTextureKind.h"
#include "IGraphicsSurface.h"
#include "GraphicsSurfaceContext.h"
#include "GraphicsDevice.h"
#include "GraphicsResource.h"
#include "GraphicsBufferKind.h"
#include "GraphicsBuffer.h"
#include "GraphicsTextureAddressMode.h"
#include "GraphicsTexture.h"
#include "IGraphicsAdapterSelector.h"
#include "ShaderProgram.h"
#include "GraphicsPipelineBlendFactor.h"
#include "GraphicsPipeline.h"
#include "GraphicsPipelineSignature.h"
#include "GraphicsPrimitive.h"
#include "GraphicsProvider.h"
#include "GraphicsPipelineInput.h"
#include "GraphicsPipelineInputElement.h"
#include "GraphicsPipelineInputElementKind.h"
#include "GraphicsPipelineResource.h"
#include "GraphicsPipelineResourceKind.h"
#include "ShaderProgramVisibility.h"
#include "GraphicsResourceManager.h"

// clang-format on

#endif // !NOVELRT_GRAPHICS_H

max_stars_repo_name: Shidesu/NovelRT

max_issues_repo_name: Shidesu/NovelRT

max_forks_repo_name: Shidesu/NovelRT


text:
/*
 * Copyright (c) Facebook, Inc. and its affiliates.
 * All rights reserved.
 * This source code is licensed under the BSD-style license found in the
 * LICENSE file in the root directory of this source tree.
 */
#pragma once
#include <chrono>
#include <functional>
#include <vector>

#include <immintrin.h>

#ifdef USE_BLAS
#if __APPLE__
// not sure whether need to differentiate TARGET_OS_MAC or TARGET_OS_IPHONE,
// etc.
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif
#endif

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

#ifdef USE_MKL
#include <mkl.h>
#endif

#include "./AlignedVec.h"
#include "fbgemm/FbgemmBuild.h"
#include "fbgemm/FbgemmPackMatrixB.h"
#include "src/RefImplementations.h"

namespace fbgemm {

template <typename T>
void randFill(aligned_vector<T>& vec, T low, T high);

void llc_flush(std::vector<char>& llc);

// Same as omp_get_max_threads() when OpenMP is available, otherwise 1
int fbgemm_get_max_threads();
// Same as omp_get_num_threads() when OpenMP is available, otherwise 1
int fbgemm_get_num_threads();
// Same as omp_get_thread_num() when OpenMP is available, otherwise 0
int fbgemm_get_thread_num();

template <typename T>
NOINLINE float cache_evict(const T& vec) {
  auto const size = vec.size();
  auto const elemSize = sizeof(typename T::value_type);
  auto const dataSize = size * elemSize;

  const char* data = reinterpret_cast<const char*>(vec.data());
  constexpr int CACHE_LINE_SIZE = 64;
  // Not having this dummy computation significantly slows down the computation
  // that follows.
  float dummy = 0.0f;
  for (std::size_t i = 0; i < dataSize; i += CACHE_LINE_SIZE) {
    dummy += data[i] * 1.0f;
    _mm_mfence();
#ifndef _MSC_VER
    asm volatile("" ::: "memory");
#endif
    _mm_clflush(&data[i]);
  }

  return dummy;
}

/**
 * Parse application command line arguments
 *
 */
int parseArgumentInt(
    int argc,
    const char* argv[],
    const char* arg,
    int non_exist_val,
    int def_val);
bool parseArgumentBool(
    int argc,
    const char* argv[],
    const char* arg,
    bool def_val);

namespace {
struct empty_flush {
  void operator()() const {}
};
} // namespace

/**
 * @param Fn functor to execute
 * @param Fe data eviction functor
 */
template <class Fn, class Fe = std::function<void()>>
double measureWithWarmup(
    Fn&& fn,
    int warmupIterations,
    int measuredIterations,
    const Fe& fe = empty_flush(),
    bool useOpenMP = false) {
  for (int i = 0; i < warmupIterations; ++i) {
    // Evict data first
    fe();
    fn();
  }

  double ttot = 0.0;

#ifdef _OPENMP
#pragma omp parallel if (useOpenMP)
  {
#endif
    for (int i = 0; i < measuredIterations; ++i) {
      int thread_id = 0;
      std::chrono::time_point<std::chrono::high_resolution_clock> start, end;

#ifdef _OPENMP
      if (useOpenMP) {
        thread_id = omp_get_thread_num();
      }
#endif

      if (thread_id == 0) {
        fe();
      }

#ifdef _OPENMP
      if (useOpenMP) {
#pragma omp barrier
      }
#endif
      start = std::chrono::high_resolution_clock::now();

      fn();

#ifdef _OPENMP
      if (useOpenMP) {
#pragma omp barrier
      }
#endif

      end = std::chrono::high_resolution_clock::now();
      auto dur =
          std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);

      if (thread_id == 0) {
        // TODO: measure load imbalance
        ttot += dur.count();
      }
    }

#ifdef _OPENMP
  }
#endif
  return ttot / 1e9 / measuredIterations;
}

/*
 * @brief Out-of-place transposition for M*N matrix ref.
 * @param M number of rows in input
 * @param K number of columns in input
 */
template <typename T>
void transpose_matrix(
    int M,
    int N,
    const T* src,
    int ld_src,
    T* dst,
    int ld_dst) {
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < M; ++j) {
      dst[i * ld_dst + j] = src[i + j * ld_src];
    }
  } // for each output row
}

/*
 * @brief In-place transposition for nxk matrix ref.
 * @param n number of rows in input (number of columns in output)
 * @param k number of columns in input (number of rows in output)
 */
template <typename T>
void transpose_matrix(T* ref, int n, int k) {
  std::vector<T> local(n * k);
  transpose_matrix(n, k, ref, k, local.data(), n);
  memcpy(ref, local.data(), n * k * sizeof(T));
}

#if defined(USE_MKL)
void test_xerbla(char* srname, const int* info, int);
#endif

#define dataset 1

template <typename btype>
void performance_test(
    int num_instances,
    bool flush,
    int repetitions,
    bool is_mkl) {
#if defined(USE_MKL)
  mkl_set_xerbla((XerblaEntry)test_xerbla);
#endif

  float alpha = 1.f, beta = 1.f;
  matrix_op_t btran = matrix_op_t::Transpose;

#if dataset == 1
  const int NITER = (flush) ? 10 : 100;
  std::vector<std::vector<int>> shapes;
  for (auto m = 1; m < 120; m++) {
    // shapes.push_back({m, 128, 512});
    shapes.push_back({m, 512, 512});
  }

#elif dataset == 2
  const int NITER = (flush) ? 10 : 100;
#include "shapes_dataset.h"

#else
  flush = false;
  constexpr int NITER = 1;
  std::vector<std::vector<int>> shapes;
  std::random_device r;
  std::default_random_engine generator(r());
  std::uniform_int_distribution<int> dm(1, 100);
  std::uniform_int_distribution<int> dnk(1, 1024);
  for (int i = 0; i < 1000; i++) {
    int m = dm(generator);
    int n = dnk(generator);
    int k = dnk(generator);
    shapes.push_back({m, n, k});
  }
#endif

  std::string type;
  double gflops, gbs, ttot;
  for (auto s : shapes) {
    int m = s[0];
    int n = s[1];
    int k = s[2];

    // initialize with small numbers
    aligned_vector<int> Aint(m * k);
    randFill(Aint, 0, 4);
    std::vector<aligned_vector<float>> A;
    for (int i = 0; i < num_instances; ++i) {
      A.push_back(aligned_vector<float>(Aint.begin(), Aint.end()));
    }

    aligned_vector<int> Bint(k * n);
    randFill(Bint, 0, 4);
    aligned_vector<float> B(Bint.begin(), Bint.end());
    std::vector<std::unique_ptr<PackedGemmMatrixB<btype>>> Bp;
    for (int i = 0; i < num_instances; ++i) {
      Bp.emplace_back(std::unique_ptr<PackedGemmMatrixB<btype>>(
          new PackedGemmMatrixB<btype>(btran, k, n, alpha, B.data())));
    }
    auto kAligned = ((k * sizeof(float) + 64) & ~63) / sizeof(float);
    auto nAligned = ((n * sizeof(float) + 64) & ~63) / sizeof(float);
    std::vector<aligned_vector<float>> Bt(num_instances);
    auto& Bt_ref = Bt[0];

    if (btran == matrix_op_t::Transpose) {
      Bt_ref.resize(k * nAligned);
      for (auto row = 0; row < k; ++row) {
        for (auto col = 0; col < n; ++col) {
          Bt_ref[row * nAligned + col] = alpha * B[col * k + row];
        }
      }
    } else {
      Bt_ref.resize(kAligned * n);
      for (auto row = 0; row < k; ++row) {
        for (auto col = 0; col < n; ++col) {
          Bt_ref[col * kAligned + row] = alpha * B[col * k + row];
        }
      }
    }

    for (auto i = 1; i < num_instances; ++i) {
      Bt[i] = Bt_ref;
    }

    std::vector<aligned_vector<float>> C_ref;
    std::vector<aligned_vector<float>> C_fb;
    if (beta != 0.0f) {
      aligned_vector<int> Cint(m * n);
      randFill(Cint, 0, 4);
      for (int i = 0; i < num_instances; ++i) {
        C_ref.push_back(aligned_vector<float>(Cint.begin(), Cint.end()));
        C_fb.push_back(aligned_vector<float>(Cint.begin(), Cint.end()));
      }
    } else {
      for (int i = 0; i < num_instances; ++i) {
        C_ref.push_back(aligned_vector<float>(m * n, 1.f));
        C_fb.push_back(aligned_vector<float>(m * n, NAN));
      }
    }

    double nflops = 2.0 * m * n * k;
    double nbytes = 4.0 * m * k + sizeof(btype) * 1.0 * k * n + 4.0 * m * n;

    // warm up MKL and fbgemm
    // check correctness at the same time
    for (auto w = 0; w < 3; w++) {
#if defined(USE_MKL) || defined(USE_BLAS)
      cblas_sgemm(
          CblasRowMajor,
          CblasNoTrans,
          CblasNoTrans, // B is pretransposed, if required by operation
          m,
          n,
          k,
          1.0, // Mutliplication by Alpha is done during transpose of B
          A[0].data(),
          k,
          Bt[0].data(),
          btran == matrix_op_t::NoTranspose ? kAligned : nAligned,
          beta,
          C_ref[0].data(),
          n);
#else
      cblas_sgemm_ref(
          matrix_op_t::NoTranspose,
          matrix_op_t::NoTranspose,
          m,
          n,
          k,
          1.0,
          A[0].data(),
          k,
          Bt[0].data(),
          (btran == matrix_op_t::NoTranspose) ? kAligned : nAligned,
          beta,
          C_ref[0].data(),
          n);
#endif
#ifdef _OPENMP
#pragma omp parallel if (num_instances == 1)
#endif
      {
        int num_threads = num_instances == 1 ? fbgemm_get_num_threads() : 1;
        int tid = num_instances == 1 ? fbgemm_get_thread_num() : 0;
        cblas_gemm_compute(
            matrix_op_t::NoTranspose,
            m,
            A[0].data(),
            *Bp[0],
            beta,
            C_fb[0].data(),
            tid,
            num_threads);
      }

#if defined(USE_MKL) || defined(USE_BLAS)
      // Compare results
      for (auto i = 0; i < C_ref[0].size(); i++) {
        if (std::abs(C_ref[0][i] - C_fb[0][i]) > 1e-3) {
          fprintf(
              stderr,
              "Error: too high diff between fp32 ref %f and fp16 %f at %d\n",
              C_ref[0][i],
              C_fb[0][i],
              i);
          return;
        }
      }
#endif
    }

#if defined(USE_MKL)
    if (is_mkl) {
      // Gold via MKL sgemm
      type = "MKL_FP32";
#elif defined(USE_BLAS)
    type = "BLAS_FP32";
#else
    type = "REF_FP32";
#endif

      ttot = measureWithWarmup(
          [&]() {
            int copy = num_instances == 1 ? 0 : fbgemm_get_thread_num();
            for (int i = 0; i < repetitions; ++i) {
#if defined(USE_MKL) || defined(USE_BLAS)
              cblas_sgemm(
                  CblasRowMajor,
                  CblasNoTrans,
                  CblasNoTrans,
                  m,
                  n,
                  k,
                  1.0,
                  A[copy].data(),
                  k,
                  Bt[copy].data(),
                  btran == matrix_op_t::NoTranspose ? kAligned : nAligned,
                  beta,
                  C_ref[copy].data(),
                  n);
#else
            cblas_sgemm_ref(
                matrix_op_t::NoTranspose,
                matrix_op_t::NoTranspose,
                m,
                n,
                k,
                1.0,
                A[copy].data(),
                k,
                Bt[copy].data(),
                (btran == matrix_op_t::NoTranspose) ? kAligned : nAligned,
                beta,
                C_ref[copy].data(),
                n);
#endif
            }
          },
          3,
          NITER,
          [&]() {
            if (flush) {
              int copy = num_instances == 1 ? 0 : fbgemm_get_thread_num();
              cache_evict(A[copy]);
              cache_evict(Bt[copy]);
              cache_evict(C_ref[copy]);
            }
          },
          // Use OpenMP if num instances > 1
          num_instances > 1);

      gflops = nflops / ttot / 1e9;
      gbs = nbytes / ttot / 1e9;
      printf(
          "\n%30s m = %5d n = %5d k = %5d Gflops = %8.4lf GBytes = %8.4lf\n",
          type.c_str(),
          m,
          n,
          k,
          gflops * repetitions,
          gbs * repetitions);
#ifdef USE_MKL
    }
#endif
    type = "FBP_" + std::string(typeid(btype).name());

    ttot = measureWithWarmup(
        [&]() {
          // When executing in data decomposition (single-instance) mode
          // Different threads will access different regions of the same
          // matrices. Thus, copy to be used is always 0. The numbers of
          // threads would be the as number of threads in the parallel
          // region.
          // When running in functional decomposition (multi-instance) mode
          // different matrices are used. The copy to be used selected by
          // thread_id (thread_num), and the number of threads performance
          // the compute of the same instance is 1.
          int copy = num_instances == 1 ? 0 : fbgemm_get_thread_num();
          int num_threads = num_instances == 1 ? fbgemm_get_num_threads() : 1;
          int tid = num_instances == 1 ? fbgemm_get_thread_num() : 0;

          for (int i = 0; i < repetitions; ++i) {
            cblas_gemm_compute(
                matrix_op_t::NoTranspose,
                m,
                A[copy].data(),
                *Bp[copy],
                beta,
                C_fb[copy].data(),
                tid,
                num_threads);
          }
        },
        3,
        NITER,
        [&]() {
          if (flush) {
            int copy = num_instances == 1 ? 0 : fbgemm_get_thread_num();
            cache_evict(A[copy]);
            cache_evict(*Bp[copy]);
            cache_evict(C_fb[copy]);
          }
        },
        true /*useOpenMP*/);

    gflops = nflops / ttot / 1e9;
    gbs = nbytes / ttot / 1e9;
    printf(
        "%30s m = %5d n = %5d k = %5d Gflops = %8.4lf GBytes = %8.4lf\n",
        type.c_str(),
        m,
        n,
        k,
        gflops * repetitions,
        gbs * repetitions);
  }
}

} // namespace fbgemm


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

#ifndef __GSL_MATRIX_LONG_DOUBLE_H__
#define __GSL_MATRIX_LONG_DOUBLE_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_long_double.h>
#include <gsl/gsl_blas_types.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_long_double * block;
  int owner;
} gsl_matrix_long_double;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_view;

typedef _gsl_matrix_long_double_view gsl_matrix_long_double_view;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_const_view;

typedef const _gsl_matrix_long_double_const_view gsl_matrix_long_double_const_view;

/* Allocation */

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_block (gsl_block_long_double * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_matrix (gsl_matrix_long_double * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_long_double *
gsl_vector_long_double_alloc_row_from_matrix (gsl_matrix_long_double * m,
                                        const size_t i);

GSL_FUN gsl_vector_long_double *
gsl_vector_long_double_alloc_col_from_matrix (gsl_matrix_long_double * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_long_double_free (gsl_matrix_long_double * m);

/* Views */

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_submatrix (gsl_matrix_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_row (gsl_matrix_long_double * m, const size_t i);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_column (gsl_matrix_long_double * m, const size_t j);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_diagonal (gsl_matrix_long_double * m);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_superdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subrow (gsl_matrix_long_double * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subcolumn (gsl_matrix_long_double * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector (gsl_vector_long_double * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector_with_tda (gsl_vector_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_submatrix (const gsl_matrix_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_row (const gsl_matrix_long_double * m,
                            const size_t i);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_column (const gsl_matrix_long_double * m,
                               const size_t j);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_diagonal (const gsl_matrix_long_double * m);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subdiagonal (const gsl_matrix_long_double * m,
                                    const size_t k);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_superdiagonal (const gsl_matrix_long_double * m,
                                      const size_t k);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subrow (const gsl_matrix_long_double * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subcolumn (const gsl_matrix_long_double * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector (const gsl_vector_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector_with_tda (const gsl_vector_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_long_double_set_zero (gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_set_identity (gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_set_all (gsl_matrix_long_double * m, long double x);

GSL_FUN int gsl_matrix_long_double_fread (FILE * stream, gsl_matrix_long_double * m) ;
GSL_FUN int gsl_matrix_long_double_fwrite (FILE * stream, const gsl_matrix_long_double * m) ;
GSL_FUN int gsl_matrix_long_double_fscanf (FILE * stream, gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_fprintf (FILE * stream, const gsl_matrix_long_double * m, const char * format);

GSL_FUN int gsl_matrix_long_double_memcpy(gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);
GSL_FUN int gsl_matrix_long_double_swap(gsl_matrix_long_double * m1, gsl_matrix_long_double * m2);
GSL_FUN int gsl_matrix_long_double_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);

GSL_FUN int gsl_matrix_long_double_swap_rows(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_swap_columns(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_swap_rowcol(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_transpose (gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_transpose_memcpy (gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);
GSL_FUN int gsl_matrix_long_double_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);

GSL_FUN long double gsl_matrix_long_double_max (const gsl_matrix_long_double * m);
GSL_FUN long double gsl_matrix_long_double_min (const gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_minmax (const gsl_matrix_long_double * m, long double * min_out, long double * max_out);

GSL_FUN void gsl_matrix_long_double_max_index (const gsl_matrix_long_double * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_long_double_min_index (const gsl_matrix_long_double * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_long_double_minmax_index (const gsl_matrix_long_double * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_long_double_equal (const gsl_matrix_long_double * a, const gsl_matrix_long_double * b);

GSL_FUN int gsl_matrix_long_double_isnull (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_ispos (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_isneg (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_isnonneg (const gsl_matrix_long_double * m);

GSL_FUN int gsl_matrix_long_double_add (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_sub (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_mul_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_div_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_scale (gsl_matrix_long_double * a, const double x);
GSL_FUN int gsl_matrix_long_double_scale_rows (gsl_matrix_long_double * a, const gsl_vector_long_double * x);
GSL_FUN int gsl_matrix_long_double_scale_columns (gsl_matrix_long_double * a, const gsl_vector_long_double * x);
GSL_FUN int gsl_matrix_long_double_add_constant (gsl_matrix_long_double * a, const double x);
GSL_FUN int gsl_matrix_long_double_add_diagonal (gsl_matrix_long_double * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_long_double_get_row(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t i);
GSL_FUN int gsl_matrix_long_double_get_col(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t j);
GSL_FUN int gsl_matrix_long_double_set_row(gsl_matrix_long_double * m, const size_t i, const gsl_vector_long_double * v);
GSL_FUN int gsl_matrix_long_double_set_col(gsl_matrix_long_double * m, const size_t j, const gsl_vector_long_double * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL long double   gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x);
GSL_FUN INLINE_DECL long double * gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const long double * gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j);

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

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

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

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

#endif

__END_DECLS

#endif /* __GSL_MATRIX_LONG_DOUBLE_H__ */


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions are met: */

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

/* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND */
/* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED */
/* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
/* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE */
/* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL */
/* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR */
/* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, */
/* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE */
/* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */

#ifndef UKF_PARAMETER_NDIM_H
#define UKF_PARAMETER_NDIM_H

#include <gsl/gsl_linalg.h> // For the Cholesky decomposition
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>

#include "ukf_types.h"

namespace ukf
{
    namespace parameter
    {

        /**
          * @short Allocation of the vectors/matrices and initialization
          *
          */
        void ukf_init(ukf_param &p, ukf_state &s)
        {
            // Init the lambda
            p.lambda = p.alpha * p.alpha * (p.n + p.kpa) - p.n;
            p.gamma = sqrt(p.n + p.lambda);
            p.nbSamples = 2 * p.n + 1;

            // Init the matrices used to iterate
            s.Kk = gsl_matrix_alloc(p.n,p.no); // Kalman gain
            gsl_matrix_set_zero(s.Kk);

            s.Kk_T = gsl_matrix_alloc(p.no,p.n);
            gsl_matrix_set_zero(s.Kk_T);

            s.Pwdk = gsl_matrix_alloc(p.n,p.no);
            gsl_matrix_set_zero(s.Pwdk);

            // Whatever the type of evolution noise, its covariance is set to evolution_noise
            s.Prrk = gsl_matrix_alloc(p.n,p.n);
            p.evolution_noise->init(p,s);

            // Whatever the type of observation noise, its covariance is set to observation_noise
            s.Peek = gsl_matrix_alloc(p.no,p.no);
            gsl_matrix_set_identity(s.Peek);
            gsl_matrix_scale(s.Peek, p.observation_noise);

            s.Pddk = gsl_matrix_alloc(p.no, p.no); // Covariance of the output
            gsl_matrix_set_zero(s.Pddk);

            s.w = gsl_vector_alloc(p.n); // Parameter vector
            gsl_vector_set_zero(s.w);

            s.wk = gsl_vector_alloc(p.n); // Vector holding one sigma point
            gsl_vector_set_zero(s.wk);

            s.Pk = gsl_matrix_alloc(p.n,p.n); // Covariance matrix
            gsl_matrix_set_identity(s.Pk);
            gsl_matrix_scale(s.Pk,p.prior_pi);

            s.Sk = gsl_matrix_alloc(p.n,p.n); // Matrix holding the cholesky decomposition of Pk
            // Initialize Sk to the cholesky decomposition of Pk
            gsl_matrix_memcpy(s.Sk, s.Pk);
            gsl_linalg_cholesky_decomp(s.Sk);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; k++)
                for(int j = 0 ; j < k ; j++)
                    gsl_matrix_set(s.Sk,j,k,0.0);

            s.cSk = gsl_vector_alloc(p.n); // Vector holding one column of Lpi
            gsl_vector_set_zero(s.cSk);

            s.wm  = gsl_vector_alloc(p.nbSamples); // Weights used to compute the mean of the sigma points images
            s.wc = gsl_vector_alloc(p.nbSamples); // Weights used to update the covariance matrices

            // Set the weights
            gsl_vector_set(s.wm, 0, p.lambda / (p.n + p.lambda));
            gsl_vector_set(s.wc, 0, p.lambda / (p.n + p.lambda) + (1.0 - p.alpha*p.alpha + p.beta));
            for(int j = 1 ; j < p.nbSamples; j ++)
            {
                gsl_vector_set(s.wm, j, 1.0 / (2.0 * (p.n + p.lambda)));
                gsl_vector_set(s.wc, j, 1.0 / (2.0 * (p.n + p.lambda)));
            }

            s.dk = gsl_matrix_alloc(p.no, p.nbSamples); // Holds the image of the sigma points
            gsl_matrix_set_zero(s.dk);

            s.ino_dk = gsl_vector_alloc(p.no); // Holds the inovation
            gsl_vector_set_zero(s.ino_dk);

            s.d_mean = gsl_vector_alloc(p.no); // Holds the mean of the sigma points images
            gsl_vector_set_zero(s.d_mean);

            s.sigmaPoints = gsl_matrix_alloc(p.n,p.nbSamples); // Holds the sigma points in the columns
            gsl_matrix_set_zero(s.sigmaPoints);

            // Temporary vectors/matrices
            s.vec_temp_n = gsl_vector_alloc(p.n);
            s.vec_temp_output = gsl_vector_alloc(p.no);

            s.mat_temp_n_1 = gsl_matrix_alloc(p.n,1);
            s.mat_temp_n_output = gsl_matrix_alloc(p.n, p.no);
            s.mat_temp_output_n = gsl_matrix_alloc(p.no, p.n);
            s.mat_temp_1_output = gsl_matrix_alloc(1,p.no);
            s.mat_temp_output_1 = gsl_matrix_alloc(p.no, 1);
            s.mat_temp_output_output = gsl_matrix_alloc(p.no, p.no);
            s.mat_temp_n_n = gsl_matrix_alloc(p.n, p.n);
        }

        /**
          * @short Free of memory allocation
          *
          */
        void ukf_free(ukf_param &p, ukf_state &s)
        {
            gsl_matrix_free(s.Kk);
            gsl_matrix_free(s.Kk_T);
            gsl_matrix_free(s.Pwdk);
            gsl_matrix_free(s.Pddk);
            gsl_matrix_free(s.Peek);
            gsl_matrix_free(s.Prrk);

            gsl_vector_free(s.w);
            gsl_vector_free(s.wk);

            gsl_matrix_free(s.Pk);
            gsl_matrix_free(s.Sk);
            gsl_vector_free(s.cSk);

            gsl_vector_free(s.wm);
            gsl_vector_free(s.wc);

            gsl_matrix_free(s.dk);
            gsl_vector_free(s.ino_dk);
            gsl_vector_free(s.d_mean);

            gsl_matrix_free(s.sigmaPoints);

            gsl_vector_free(s.vec_temp_n);
            gsl_vector_free(s.vec_temp_output);

            gsl_matrix_free(s.mat_temp_n_1);
            gsl_matrix_free(s.mat_temp_n_output);
            gsl_matrix_free(s.mat_temp_output_n);
            gsl_matrix_free(s.mat_temp_1_output);
            gsl_matrix_free(s.mat_temp_output_1);

            gsl_matrix_free(s.mat_temp_output_output);
            gsl_matrix_free(s.mat_temp_n_n);
        }

        /**
          * @short Iteration for the statistical linearization
          *
          */
	template<typename GFUNC>
	void ukf_iterate(ukf_param &p, ukf_state &s,
			 GFUNC g,
			 gsl_vector * xk, gsl_vector* dk)
        {

            // Here, we implement the UKF for parameter estimation in the vectorial case
            // The notations follow p93 of the PhD thesis of Van Der Merwe, "Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models"

            // ************************************************** //
            // ************ Time update equations    ************ //
            // ************************************************** //
            // Add the evolution noise to the parameter covariance Eq 3.137
            gsl_matrix_add(s.Pk, s.Prrk);

            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //
            // Equations 3.138
            // w_k^j = w_(k-1)  <-- this is here denoted s.w
            // w_k^j = w_(k-1) + gamma Sk_j for 1 <= j <= n
            // w_k^j = w_(k-1) - gamma Sk_j for n+1 <= j <= 2n

            // Perform a cholesky decomposition of Pk
            gsl_matrix_memcpy(s.Sk, s.Pk);
            gsl_linalg_cholesky_decomp(s.Sk);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; ++k)
                for(int j = 0 ; j < k ; ++j)
                    gsl_matrix_set(s.Sk,j,k,0.0);

            gsl_matrix_set_col(s.sigmaPoints,0, s.w);
            for(int j = 1 ; j < p.n+1 ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) + p.gamma * gsl_matrix_get(s.Sk,i,j-1));

            for(int j = p.n+1 ; j < p.nbSamples ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) - p.gamma * gsl_matrix_get(s.Sk,i,j-(p.n+1)));

            /**************************************************/
            /***** Compute the images of the sigma points *****/
            /**************************************************/

            // Compute the images of the sigma points
            // and the mean of the dk
            gsl_vector_set_zero(s.d_mean);
            for(int j = 0 ; j < p.nbSamples ; j++)
            {
                // Equation 3.129
                gsl_matrix_get_col(s.wk, s.sigmaPoints,j);
                g(s.wk,xk, s.vec_temp_output);
                gsl_matrix_set_col(s.dk, j, s.vec_temp_output);

                // Equation 3.140
                // Update the mean : y_mean = sum_[j=0..2n] w_j y_j
                gsl_vector_scale(s.vec_temp_output, gsl_vector_get(s.wm,j));
                gsl_vector_add(s.d_mean, s.vec_temp_output);
            }

            /**************************************************/
            /************** Update the statistics *************/
            /**************************************************/

            gsl_matrix_set_zero(s.Pwdk);
            gsl_matrix_memcpy(s.Pddk, s.Peek); // Add R^e_k to Pddk, Eq 3.142
            for(int j = 0 ; j < p.nbSamples ; ++j)
            {
                // Update of Pwdk
                // (wk - w)
                gsl_matrix_get_col(s.wk, s.sigmaPoints,j);
                gsl_vector_sub(s.wk, s.w);
                gsl_matrix_set_col(s.mat_temp_n_1, 0, s.wk);

                // (dk - d_mean)
                gsl_matrix_get_col(s.vec_temp_output, s.dk, j);
                gsl_vector_sub(s.vec_temp_output, s.d_mean);
                gsl_matrix_set_col(s.mat_temp_output_1, 0, s.vec_temp_output);

                // compute wc_j . (wk - w_mean) * (dk - d_mean)^T
                gsl_blas_dgemm(CblasNoTrans, CblasTrans, gsl_vector_get(s.wc,j) , s.mat_temp_n_1, s.mat_temp_output_1, 0.0, s.mat_temp_n_output);

                // Equation 3.142
                // And add it to Pwdk
                gsl_matrix_add(s.Pwdk, s.mat_temp_n_output);

                // Equation 3.143
                // Update of Pddk
                gsl_matrix_get_col(s.vec_temp_output, s.dk, j);
                gsl_vector_sub(s.vec_temp_output, s.d_mean);
                gsl_matrix_set_col(s.mat_temp_output_1, 0, s.vec_temp_output);
                gsl_blas_dgemm(CblasNoTrans, CblasTrans, gsl_vector_get(s.wc,j) , s.mat_temp_output_1, s.mat_temp_output_1, 0.0, s.mat_temp_output_output);

                gsl_matrix_add(s.Pddk, s.mat_temp_output_output);
            }

            // ************************************************** //
            // ******* Kalman gain and parameters update ******** //
            // ************************************************** //

            //*** Ki = Pwdk Pddk^-1
            // Compute the inverse of Pddk
            gsl_matrix_memcpy(s.mat_temp_output_output, s.Pddk);
            gsl_linalg_cholesky_decomp(s.mat_temp_output_output);
            gsl_linalg_cholesky_invert(s.mat_temp_output_output);

            // Compute the product : Pwdk . Pddk^-1
            // Equation 3.144
            gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Pwdk, s.mat_temp_output_output, 0.0, s.Kk);

            // Update of the parameters
            // wk = w_(k-1) + Kk * (dk - d_mean)
            // Equation 3.145

            // Set the inovations
            /*for(int i = 0 ; i < p.no; ++i)
                s.ino_dk->data[i] = dk->data[i] - s.d_mean->data[i];
            gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Kk, &gsl_matrix_view_array(s.ino_dk->data,p.no,1).matrix, 0.0, s.mat_temp_n_1);
            gsl_matrix_get_col(s.vec_temp_n, s.mat_temp_n_1, 0);
            gsl_vector_add(s.w, s.vec_temp_n);*/

            for(int i = 0 ; i < p.no; ++i)
                s.ino_dk->data[i] = dk->data[i] - s.d_mean->data[i];
            gsl_blas_dgemv(CblasNoTrans, 1.0, s.Kk, s.ino_dk, 1.0, s.w);

            // Update of the parameter covariance
            // Pk = P_(k-1) - Kk Pddk Kk^T
            // Equation 3.146
            gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Kk, s.Pddk, 0.0, s.mat_temp_n_output);
            gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, s.mat_temp_n_output , s.Kk, 0.0, s.mat_temp_n_n);
            gsl_matrix_sub(s.Pk, s.mat_temp_n_n);

            // Update of the evolution noise
            p.evolution_noise->updateEvolutionNoise(p, s);
        }

        /**
          * @short Evaluation of the output from the sigma points
          *
          */
      void ukf_evaluate(ukf_param &p, ukf_state &s,
			void(*g)(gsl_vector*, gsl_vector*, gsl_vector*),
			gsl_vector * xk, gsl_vector * dk)
        {
            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //
            // Equations 3.138
            // w_k^j = w_(k-1)  <-- this is here denoted s.w
            // w_k^j = w_(k-1) + gamma Sk_j for 1 <= j <= n
            // w_k^j = w_(k-1) - gamma Sk_j for n+1 <= j <= 2n

            // Perform a cholesky decomposition of Pk
            gsl_matrix_memcpy(s.mat_temp_n_n, s.Pk);
            gsl_linalg_cholesky_decomp(s.mat_temp_n_n);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; ++k)
                for(int j = 0 ; j < k ; ++j)
                    gsl_matrix_set(s.mat_temp_n_n,j,k,0.0);

            gsl_matrix_set_col(s.sigmaPoints,0, s.w);
            for(int j = 1 ; j < p.n+1 ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) + p.gamma * gsl_matrix_get(s.mat_temp_n_n,i,j-1));

            for(int j = p.n+1 ; j < p.nbSamples ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) - p.gamma * gsl_matrix_get(s.mat_temp_n_n,i,j-(p.n+1)));

            /**************************************************/
            /***** Compute the images of the sigma points *****/
            /**************************************************/

            // Compute the images of the sigma points
            // and the mean of the dk
            gsl_vector_set_zero(dk);
            for(int j = 0 ; j < p.nbSamples ; j++)
            {
                // Equation 3.129
                gsl_matrix_get_col(s.wk, s.sigmaPoints,j);
                g(s.wk,xk, s.vec_temp_output);
                gsl_matrix_set_col(s.dk, j, s.vec_temp_output);

                // Equation 3.140
                // Update the mean : y_mean = sum_[j=0..2n] w_j y_j
                gsl_vector_scale(s.vec_temp_output, gsl_vector_get(s.wm,j));
                gsl_vector_add(dk, s.vec_temp_output);
            }
        }

        /**
          * @short Returns a set of sigma points
          */
        void getSigmaPoints(ukf_param &p, ukf_state &s, gsl_matrix * sigmaPoints)
        {
            //gsl_matrix * sigmaPoints = gsl_matrix_alloc(p.n, p.nbSamples);

            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //
            // Equations 3.138
            // w_k^j = w_(k-1)  <-- this is here denoted s.w
            // w_k^j = w_(k-1) + gamma Sk_j for 1 <= j <= n
            // w_k^j = w_(k-1) - gamma Sk_j for n+1 <= j <= 2n

            // Perform a cholesky decomposition of Pk
            gsl_matrix_memcpy(s.mat_temp_n_n, s.Pk);
            gsl_linalg_cholesky_decomp(s.mat_temp_n_n);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; ++k)
                for(int j = 0 ; j < k ; ++j)
                    gsl_matrix_set(s.mat_temp_n_n,j,k,0.0);

            gsl_matrix_set_col(sigmaPoints,0, s.w);
            for(int j = 1 ; j < p.n+1 ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(sigmaPoints, i, j, gsl_vector_get(s.w, i) + p.gamma * gsl_matrix_get(s.mat_temp_n_n,i,j-1));

            for(int j = p.n+1 ; j < p.nbSamples ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(sigmaPoints, i, j, gsl_vector_get(s.w, i) - p.gamma * gsl_matrix_get(s.mat_temp_n_n,i,j-(p.n+1)));

        }

    } // parameter
} // ukf

#endif // SL_PARAMETER_NDIM_H


text:
/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version __VERSION__, see http://ghmm.org
*
*       Filename: ghmm/ghmm/matrixop.c
*       Authors:  Christoph Hafemeister
*
*       Copyright (C) 2007-2008 Alexander Schliep
*	Copyright (C) 2007-2008 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: schliep@ghmm.org
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library 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
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*******************************************************************************/

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif

#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>

#include "matrixop.h"

#ifdef DO_WITH_GSL
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#endif

#ifdef DO_WITH_ATLAS
#include <clapack.h>
#endif

/*============================================================================*/
int ighmm_invert_det(double *sigmainv, double *det, int length, double *cov)
{
#define CUR_PROC "invert_det"
#ifdef DO_WITH_GSL
  int i, j, s;
  gsl_matrix *tmp;
  gsl_matrix *inv;
  tmp = gsl_matrix_alloc(length, length);
  inv = gsl_matrix_alloc(length, length);
  gsl_permutation *permutation = gsl_permutation_calloc(length);

  for (i=0; i<length; ++i) {
    for (j=0; j<length; ++j) {
#ifdef DO_WITH_GSL_DIAGONAL_HACK
        if (i == j){
          gsl_matrix_set(tmp, i, j, cov[i*length+j]);
        }else{
          gsl_matrix_set(tmp, i, j, 0.0);
        }
#else
        gsl_matrix_set(tmp, i, j, cov[i*length+j]);
#endif
    }
  }

  gsl_linalg_LU_decomp(tmp, permutation, &s);
  gsl_linalg_LU_invert(tmp, permutation, inv);
  *det = gsl_linalg_LU_det(tmp, s);
  gsl_matrix_free(tmp);
  gsl_permutation_free(permutation);

  for (i=0; i<length; ++i) {
    for (j=0; j<length; ++j) {
       sigmainv[i*length+j] = gsl_matrix_get(inv, i, j);
    }
  }

  gsl_matrix_free(inv);
#elif defined HAVE_CLAPACK_DGETRF && HAVE_CLAPACK_DGETRI
  char sign;
  int info, i;
  int *ipiv;
  double det_tmp;

  ipiv = malloc(length * sizeof *ipiv);

  /* copy cov. matrix entries to result matrix, the rest is done in-place */
  memcpy(sigmainv, cov, length * length * sizeof *cov);

  /* perform in-place LU factorization of covariance matrix */
  info = clapack_dgetrf(CblasRowMajor, length, length, sigmainv, length, ipiv);

  /* determinant */
  sign = 1;
  for( i=0; i<length; ++i)
    if( ipiv[i]!=i )
      sign *= -1;
  det_tmp = sigmainv[0];
  for( i=length+1; i<(length*length); i+=length+1 )
    det_tmp *= sigmainv[i];
  *det = det_tmp * sign;

  /* use the LU factorization to get inverse */
  info = clapack_dgetri(CblasRowMajor, length, sigmainv, length, ipiv);

  free(ipiv);
#else
  *det = ighmm_determinant(cov, length);
  ighmm_inverse(cov, length, *det, sigmainv);
#endif

  return 0;
#undef CUR_PROC
}

/*============================================================================*/
/* calculate determinant of a square matrix */
double ighmm_determinant(double *cov, int n)
{
  int j1, i, j, jm;
  double * m;
  double det;

  if (n == 1)
    return cov[0];
  if (n == 2)
    return cov[0] * cov[3] - cov[1] * cov[2];
  /* matrix dimension is bigger than 2 - we have to do some work */
  det = 0;
  for (j1=0;j1<n;j1++) {
    m = malloc((n-1)*(n-1)*sizeof(double));
    for (i=1;i<n;i++) {
      jm = 0;
      for (j=0;j<n;j++) {
        if (j == j1)
          continue;
        m[(i-1)*(n-1)+jm] = cov[i*n+j];
        jm++;
      }
    }
    det += pow(-1.0,1.0+j1+1.0) * cov[j1] * ighmm_determinant(m,n-1);
    free(m);
  }
  return det;
}

/*============================================================================*/
/*
  The inverse of a square matrix A with a non zero determinant is the adjoint
  matrix divided by the determinant.
  The adjoint matrix is the transpose of the cofactor matrix.
  The cofactor matrix is the matrix of determinants of the minors Aij
  multiplied by -1^(i+j).
  The i,j'th minor of A is the matrix A without the i'th column or
  the j'th row.
*/
int ighmm_inverse(double *cov, int n, double det, double *inv)
{
  int i, j, ic, jc, actrow, actcol;
  double * m;

  if (n == 1) {
    inv[0] = 1/cov[0];
    return 0;
  }
  if (n == 2) {
    inv[0] = cov[3] / (cov[0]*cov[3] - cov[1]*cov[2]);
    inv[1] = - cov[1] / (cov[0]*cov[3] - cov[1]*cov[2]);
    inv[2] = - cov[2] / (cov[0]*cov[3] - cov[1]*cov[2]);
    inv[3] = cov[0] / (cov[0]*cov[3] - cov[1]*cov[2]);
    return 0;
  }
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
      /* calculate minor i,j */
      m = malloc((n-1)*(n-1)*sizeof(double));
      actrow = 0;
      for (ic=0;ic<n;ic++) {
        if (ic == i)
          continue;
        actcol = 0;
        for (jc=0;jc<n;jc++) {
          if (jc == j)
            continue;
          m[actrow*(n-1)+actcol] = cov[ic*n+jc];
          actcol++;
        }
        actrow++;
      }
      /* cofactor i,j is determinant of m times -1^(i+j) */
      inv[j*n+i] = pow(-1.0,i+j+2.0) * ighmm_determinant(m,n-1) / det;
      free(m);
    }
  }
  return 0;
}

/*============================================================================*/
int ighmm_cholesky_decomposition (double *sigmacd, int dim, double *cov) {
#ifdef DO_WITH_GSL
  int i, j;
  gsl_matrix *tmp = gsl_matrix_alloc(dim, dim);

  /* temp matrix */
  for (i=0;i<dim;i++) {
    for (j=0;j<dim;j++) {
      gsl_matrix_set(tmp, i, j, cov[i*dim+j]);
    }
  }
  /* cholesky decomposition of temp matrix */
  gsl_linalg_cholesky_decomp(tmp);
  /* copy to sigmacd and set upper triangular to 0 */
  for (i=0;i<dim;i++) {
    for (j=0;j<dim;j++) {
      if (j>i)
        sigmacd[i*dim+j] = 0;
      else
        sigmacd[i*dim+j] = gsl_matrix_get(tmp, i, j);
    }
  }
  gsl_matrix_free(tmp);

#elif defined HAVE_CLAPACK_DPOTRF
  int info;
  /* copy cov. matrix entries to result matrix, the rest is done in-place */
  memcpy(sigmacd, cov, dim * dim * sizeof *cov);
  info = clapack_dpotrf(CblasRowMajor, CblasUpper, dim, sigmacd, dim);
#else
  int row, j, k;
  double sum;

  /* copy cov to sigmacd */
  for (row=0; row<dim; row++) {
    for (j=0; j<dim; j++) {
      sigmacd[row*dim+j] = cov[row*dim+j];
    }
  }

  for (row=0; row<dim; row++) {
    /* First compute U[row][row] */
    sum = cov[row*dim+row];
    for (j=0; j<(row-1); j++) sum -= sigmacd[j*dim+row]*sigmacd[j*dim+row];
    if (sum > DBL_MIN) {
      sigmacd[row*dim+row] = sqrt(sum);
      /* Now find elements sigmacd[row*dim+k], k > row. */
      for (k=(row+1); k<dim; k++) {
        sum = cov[row*dim+k];
        for (j=0; j<(row-1); j++)
          sum -= sigmacd[j*dim+row]*sigmacd[j*dim+k];
        sigmacd[row*dim+k] = sum/sigmacd[row*dim+row];
      }
    }
    else {  /* blast off the entire row. */
      for (k=row; k<dim; k++) sigmacd[row*dim+k] = 0.0;
    }
  }
#endif

  return 0;
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
/*
 * This file is part of MXE.
 * See index.html for further information.
 */

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int main(int argc, char *argv[])
{
    double x, y;

    (void)argc;
    (void)argv;

    x = 5.0;
    y = gsl_sf_bessel_J0 (x);
    printf ("J0(%g) = %.18e\n", x, y);
    return 0;
}

max_stars_repo_name: ChristianFrisson/mxe

max_issues_repo_name: ChristianFrisson/mxe

max_forks_repo_name: ChristianFrisson/mxe


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for structures representing a waveform as a list of modes in amplitude/phase form.
 *
 */

#ifndef _STRUCT_H
#define _STRUCT_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

/***************************************************/
/*************** Type definitions ******************/

/* Type for real functions */
typedef double (*RealFunctionPtr)(double);
typedef double (*RealObjectFunctionPtr)(const void *, double);

/* Type for real functions that reference an object */
typedef struct tagObjectFunction
{
  const void * object;
  RealObjectFunctionPtr function;
} ObjectFunction ;
double ObjectFunctionCall(const ObjectFunction*,double);

/* Complex frequency series in amplitude and phase representation (for one mode) */
typedef struct tagCAmpPhaseFrequencySeries
{
  gsl_vector* freq;
  gsl_vector* amp_real; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_vector* amp_imag; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_vector* phase;
} CAmpPhaseFrequencySeries;
/* GSL splines for complex amplitude and phase representation (for one mode) */
typedef struct tagCAmpPhaseGSLSpline
{
  gsl_vector* freq;
  gsl_spline* spline_amp_real; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_spline* spline_amp_imag; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_spline* spline_phase;
  gsl_interp_accel* accel_amp_real;
  gsl_interp_accel* accel_amp_imag;
  gsl_interp_accel* accel_phase;
} CAmpPhaseGSLSpline;
/* Splines in matrix form for complex amplitude and phase representation (for one mode) */
typedef struct tagCAmpPhaseSpline
{
  gsl_matrix* spline_amp_real; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_matrix* spline_amp_imag; /* We authorize complex amplitudes - will be used for the LISA response */
  gsl_matrix* quadspline_phase;
} CAmpPhaseSpline;

/* Complex frequency series in real/imaginary part representation (for one mode, or their sum) */
typedef struct tagReImFrequencySeries
{
  gsl_vector* freq;
  gsl_vector* h_real;
  gsl_vector* h_imag;
} ReImFrequencySeries;

/* Complex frequency series in real/imaginary part representation (for one mode, or their sum) */
/* NOTE: for now, exact duplicata of ReImFrequencySeries - differentiated for readability of the code */
typedef struct tagReImTimeSeries
{
  gsl_vector* times;
  gsl_vector* h_real;
  gsl_vector* h_imag;
} ReImTimeSeries;

/* Complex frequency series in amplitude/phase representation (representing one mode) */
/* NOTE: for now, exact duplicata of ReImFrequencySeries - differentiated for readability of the code */
typedef struct tagAmpPhaseTimeSeries
{
  gsl_vector* times;
  gsl_vector* h_amp;
  gsl_vector* h_phase;
} AmpPhaseTimeSeries;

/* Real time series */
/* NOTE: could change the h to something more general, like values - also used for TD 22 amplitude, for instance  */
typedef struct tagRealTimeSeries
{
  gsl_vector* times;
  gsl_vector* h;
} RealTimeSeries;

/* List structure, for a list of modes, each in amplitude and phase form */
typedef struct tagListmodesCAmpPhaseFrequencySeries
{
  CAmpPhaseFrequencySeries*                      freqseries; /* The frequencies series with amplitude and phase */
  int                                            l; /* Node mode l  */
  int                                            m; /* Node submode m  */
  struct tagListmodesCAmpPhaseFrequencySeries*    next; /* Next pointer */
} ListmodesCAmpPhaseFrequencySeries;

/* List structure, for a list of modes, each with interpolated splines in amplitude and phase form */
typedef struct tagListmodesCAmpPhaseSpline
{
  CAmpPhaseSpline*                       splines; /* The frequencies series with amplitude and phase */
  int                                    l;       /* Node mode l  */
  int                                    m;       /* Node submode m  */
  struct tagListmodesCAmpPhaseSpline*    next;    /* Next pointer */
} ListmodesCAmpPhaseSpline;

/**************************************************************/
/* Functions computing the max and min between two int */
int max (int a, int b);
int min (int a, int b);

/**************************************************************/
/************** GSL error handling and I/O ********************/

/* GSL error handler */
void Err_Handler(const char *reason, const char *file, int line, int gsl_errno);

/* Functions to read/write data from files */
int Read_Vector(const char dir[], const char fname[], gsl_vector *v);
int Read_Matrix(const char dir[], const char fname[], gsl_matrix *m);
int Read_Text_Vector(const char dir[], const char fname[], gsl_vector *v);
int Read_Text_Matrix(const char dir[], const char fname[], gsl_matrix *m);
int Write_Vector(const char dir[], const char fname[], gsl_vector *v);
int Write_Matrix(const char dir[], const char fname[], gsl_matrix *m);
int Write_Text_Vector(const char dir[], const char fname[], gsl_vector *v);
int Write_Text_Matrix(const char dir[], const char fname[], gsl_matrix *m);

/**********************************************************/
/**************** Internal functions **********************/

/* Functions for list manipulations */
ListmodesCAmpPhaseFrequencySeries* ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(
	   ListmodesCAmpPhaseFrequencySeries* appended,  /* List structure to prepend to */
	   CAmpPhaseFrequencySeries* freqseries,  /* data to contain */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
ListmodesCAmpPhaseFrequencySeries* ListmodesCAmpPhaseFrequencySeries_GetMode(
	   ListmodesCAmpPhaseFrequencySeries* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
void ListmodesCAmpPhaseFrequencySeries_Destroy(
	   ListmodesCAmpPhaseFrequencySeries* list  /* List structure to destroy; notice that the data is destroyed too */
);
ListmodesCAmpPhaseSpline* ListmodesCAmpPhaseSpline_AddModeNoCopy(
	   ListmodesCAmpPhaseSpline* appended,  /* List structure to prepend to */
	   CAmpPhaseSpline* freqseries,  /* data to contain */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
ListmodesCAmpPhaseSpline* ListmodesCAmpPhaseSpline_GetMode(
	   ListmodesCAmpPhaseSpline* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
void ListmodesCAmpPhaseSpline_Destroy(
	   ListmodesCAmpPhaseSpline* list  /* List structure to destroy; notice that the data is destroyed too */
);

/* Functions to initialize and clean up data structure */
void CAmpPhaseFrequencySeries_Init(
	 CAmpPhaseFrequencySeries** freqseries, /* double pointer for initialization */
	 const int n );                         /* length of the frequency series */
void CAmpPhaseFrequencySeries_Cleanup(CAmpPhaseFrequencySeries* freqseries);
void CAmpPhaseSpline_Init(
	 CAmpPhaseSpline** splines,             /* double pointer for initialization */
	 const int n );                         /* length of the frequency series setting the splines */
void CAmpPhaseSpline_Cleanup(CAmpPhaseSpline* splines);
void CAmpPhaseGSLSpline_Init(
	 CAmpPhaseGSLSpline** splines,          /* double pointer for initialization */
	 const int n );                         /* length of the frequency series setting the splines */
void CAmpPhaseGSLSpline_Cleanup(CAmpPhaseGSLSpline* splines);
void ReImFrequencySeries_Init(
	 ReImFrequencySeries** freqseries,      /* double pointer for initialization */
	 const int n );                         /* length of the frequency series */
void ReImFrequencySeries_Cleanup(ReImFrequencySeries* freqseries);
void ReImTimeSeries_Init(
	 ReImTimeSeries** timeseries,           /* double pointer for initialization */
	 const int n );                         /* length of the time series */
void ReImTimeSeries_Cleanup(ReImTimeSeries* timeseries);
void AmpPhaseTimeSeries_Init(
	 AmpPhaseTimeSeries** timeseries,       /* double pointer for initialization */
	 const int n );                         /* length of the time series */
void AmpPhaseTimeSeries_Cleanup(AmpPhaseTimeSeries* timeseries);
void RealTimeSeries_Init(
	 RealTimeSeries** timeseries,      /* double pointer for initialization */
	 const int n );                /* length of the time series */
void RealTimeSeries_Cleanup(RealTimeSeries* timeseries);

/***********************************************************************/
/**************** I/O functions for internal structures ****************/

/* Note: at the moment, requires external input for the number of lines in the data */
int Read_RealTimeSeries(RealTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary);
int Read_AmpPhaseTimeSeries(AmpPhaseTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary);
int Read_ReImTimeSeries(ReImTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary);
int Write_ReImFrequencySeries(const char dir[], const char file[], ReImFrequencySeries* freqseries, const int binary);
int Write_RealTimeSeries(const char dir[], const char file[], RealTimeSeries* timeseries, int binary);
int Write_AmpPhaseTimeSeries(const char dir[], const char file[], AmpPhaseTimeSeries* timeseries, int binary);
int Write_ReImTimeSeries(const char dir[], const char file[], ReImTimeSeries* timeseries, int binary);

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _STRUCT_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare