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

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

#include <gsl/gsl_heapsort.h>

static struct peano_hilbert_data
{
  peanokey key;
  int index;
}
 *mp;

static int *Id;

void peano_hilbert_order(void)
{
  int i;

  if(ThisTask == 0)
    printf("begin Peano-Hilbert order...\n");

  if(N_gas)
    {
      mp = (struct peano_hilbert_data *) mymalloc(sizeof(struct peano_hilbert_data) * N_gas);
      Id = (int *) mymalloc(sizeof(int) * N_gas);

      for(i = 0; i < N_gas; i++)
	{
	  mp[i].index = i;
	  mp[i].key = Key[i];
	}

#ifdef MYSORT
      mysort_peano(mp, N_gas, sizeof(struct peano_hilbert_data), peano_compare_key);
#else
      qsort(mp, N_gas, sizeof(struct peano_hilbert_data), peano_compare_key);
#endif

      for(i = 0; i < N_gas; i++)
	Id[mp[i].index] = i;

      reorder_gas();

      myfree(Id);
      myfree(mp);
    }

  if(NumPart - N_gas > 0)
    {
      mp = (struct peano_hilbert_data *) mymalloc(sizeof(struct peano_hilbert_data) * (NumPart - N_gas));
      mp -= (N_gas);

      Id = (int *) mymalloc(sizeof(int) * (NumPart - N_gas));
      Id -= (N_gas);

      for(i = N_gas; i < NumPart; i++)
	{
	  mp[i].index = i;
	  mp[i].key = Key[i];
	}

#ifdef MYSORT
      mysort_peano(mp + N_gas, NumPart - N_gas, sizeof(struct peano_hilbert_data), peano_compare_key);
#else
      qsort(mp + N_gas, NumPart - N_gas, sizeof(struct peano_hilbert_data), peano_compare_key);
#endif

      for(i = N_gas; i < NumPart; i++)
	Id[mp[i].index] = i;

      reorder_particles();

      Id += N_gas;
      myfree(Id);
      mp += N_gas;
      myfree(mp);
    }

  if(ThisTask == 0)
    printf("Peano-Hilbert done.\n");
}

int peano_compare_key(const void *a, const void *b)
{
  if(((struct peano_hilbert_data *) a)->key < (((struct peano_hilbert_data *) b)->key))
    return -1;

  if(((struct peano_hilbert_data *) a)->key > (((struct peano_hilbert_data *) b)->key))
    return +1;

  return 0;
}

void reorder_gas(void)
{
  int i;
  struct particle_data Psave, Psource;
  struct sph_particle_data SphPsave, SphPsource;
  int idsource, idsave, dest;

  for(i = 0; i < N_gas; i++)
    {
      if(Id[i] != i)
	{
	  Psource = P[i];
	  SphPsource = SphP[i];

	  idsource = Id[i];
	  dest = Id[i];

	  do
	    {
	      Psave = P[dest];
	      SphPsave = SphP[dest];
	      idsave = Id[dest];

	      P[dest] = Psource;
	      SphP[dest] = SphPsource;
	      Id[dest] = idsource;

	      if(dest == i)
		break;

	      Psource = Psave;
	      SphPsource = SphPsave;
	      idsource = idsave;

	      dest = idsource;
	    }
	  while(1);
	}
    }
}

void reorder_particles(void)
{
  int i;
  struct particle_data Psave, Psource;
  int idsource, idsave, dest;

  for(i = N_gas; i < NumPart; i++)
    {
      if(Id[i] != i)
	{
	  Psource = P[i];
	  idsource = Id[i];

	  dest = Id[i];

	  do
	    {
	      Psave = P[dest];
	      idsave = Id[dest];

	      P[dest] = Psource;
	      Id[dest] = idsource;
#ifdef LT_STELLAREVOLUTION
	      if(P[dest].Type == 4)
		MetP[P[dest].MetID].PID = dest;
#endif
	      if(dest == i)
		break;

	      Psource = Psave;
	      idsource = idsave;

	      dest = idsource;
	    }
	  while(1);
	}
    }
}

/*  The following rewrite of the original function
 *  peano_hilbert_key_old() has been written by MARTIN REINECKE.
 *  It is about a factor 2.3 - 2.5 faster than Volker's old routine!
 */
const unsigned char rottable3[48][8] = {
  {36, 28, 25, 27, 10, 10, 25, 27},
  {29, 11, 24, 24, 37, 11, 26, 26},
  {8, 8, 25, 27, 30, 38, 25, 27},
  {9, 39, 24, 24, 9, 31, 26, 26},
  {40, 24, 44, 32, 40, 6, 44, 6},
  {25, 7, 33, 7, 41, 41, 45, 45},
  {4, 42, 4, 46, 26, 42, 34, 46},
  {43, 43, 47, 47, 5, 27, 5, 35},
  {33, 35, 36, 28, 33, 35, 2, 2},
  {32, 32, 29, 3, 34, 34, 37, 3},
  {33, 35, 0, 0, 33, 35, 30, 38},
  {32, 32, 1, 39, 34, 34, 1, 31},
  {24, 42, 32, 46, 14, 42, 14, 46},
  {43, 43, 47, 47, 25, 15, 33, 15},
  {40, 12, 44, 12, 40, 26, 44, 34},
  {13, 27, 13, 35, 41, 41, 45, 45},
  {28, 41, 28, 22, 38, 43, 38, 22},
  {42, 40, 23, 23, 29, 39, 29, 39},
  {41, 36, 20, 36, 43, 30, 20, 30},
  {37, 31, 37, 31, 42, 40, 21, 21},
  {28, 18, 28, 45, 38, 18, 38, 47},
  {19, 19, 46, 44, 29, 39, 29, 39},
  {16, 36, 45, 36, 16, 30, 47, 30},
  {37, 31, 37, 31, 17, 17, 46, 44},
  {12, 4, 1, 3, 34, 34, 1, 3},
  {5, 35, 0, 0, 13, 35, 2, 2},
  {32, 32, 1, 3, 6, 14, 1, 3},
  {33, 15, 0, 0, 33, 7, 2, 2},
  {16, 0, 20, 8, 16, 30, 20, 30},
  {1, 31, 9, 31, 17, 17, 21, 21},
  {28, 18, 28, 22, 2, 18, 10, 22},
  {19, 19, 23, 23, 29, 3, 29, 11},
  {9, 11, 12, 4, 9, 11, 26, 26},
  {8, 8, 5, 27, 10, 10, 13, 27},
  {9, 11, 24, 24, 9, 11, 6, 14},
  {8, 8, 25, 15, 10, 10, 25, 7},
  {0, 18, 8, 22, 38, 18, 38, 22},
  {19, 19, 23, 23, 1, 39, 9, 39},
  {16, 36, 20, 36, 16, 2, 20, 10},
  {37, 3, 37, 11, 17, 17, 21, 21},
  {4, 17, 4, 46, 14, 19, 14, 46},
  {18, 16, 47, 47, 5, 15, 5, 15},
  {17, 12, 44, 12, 19, 6, 44, 6},
  {13, 7, 13, 7, 18, 16, 45, 45},
  {4, 42, 4, 21, 14, 42, 14, 23},
  {43, 43, 22, 20, 5, 15, 5, 15},
  {40, 12, 21, 12, 40, 6, 23, 6},
  {13, 7, 13, 7, 41, 41, 22, 20}
};

const unsigned char subpix3[48][8] = {
  {0, 7, 1, 6, 3, 4, 2, 5},
  {7, 4, 6, 5, 0, 3, 1, 2},
  {4, 3, 5, 2, 7, 0, 6, 1},
  {3, 0, 2, 1, 4, 7, 5, 6},
  {1, 0, 6, 7, 2, 3, 5, 4},
  {0, 3, 7, 4, 1, 2, 6, 5},
  {3, 2, 4, 5, 0, 1, 7, 6},
  {2, 1, 5, 6, 3, 0, 4, 7},
  {6, 1, 7, 0, 5, 2, 4, 3},
  {1, 2, 0, 3, 6, 5, 7, 4},
  {2, 5, 3, 4, 1, 6, 0, 7},
  {5, 6, 4, 7, 2, 1, 3, 0},
  {7, 6, 0, 1, 4, 5, 3, 2},
  {6, 5, 1, 2, 7, 4, 0, 3},
  {5, 4, 2, 3, 6, 7, 1, 0},
  {4, 7, 3, 0, 5, 6, 2, 1},
  {6, 7, 5, 4, 1, 0, 2, 3},
  {7, 0, 4, 3, 6, 1, 5, 2},
  {0, 1, 3, 2, 7, 6, 4, 5},
  {1, 6, 2, 5, 0, 7, 3, 4},
  {2, 3, 1, 0, 5, 4, 6, 7},
  {3, 4, 0, 7, 2, 5, 1, 6},
  {4, 5, 7, 6, 3, 2, 0, 1},
  {5, 2, 6, 1, 4, 3, 7, 0},
  {7, 0, 6, 1, 4, 3, 5, 2},
  {0, 3, 1, 2, 7, 4, 6, 5},
  {3, 4, 2, 5, 0, 7, 1, 6},
  {4, 7, 5, 6, 3, 0, 2, 1},
  {6, 7, 1, 0, 5, 4, 2, 3},
  {7, 4, 0, 3, 6, 5, 1, 2},
  {4, 5, 3, 2, 7, 6, 0, 1},
  {5, 6, 2, 1, 4, 7, 3, 0},
  {1, 6, 0, 7, 2, 5, 3, 4},
  {6, 5, 7, 4, 1, 2, 0, 3},
  {5, 2, 4, 3, 6, 1, 7, 0},
  {2, 1, 3, 0, 5, 6, 4, 7},
  {0, 1, 7, 6, 3, 2, 4, 5},
  {1, 2, 6, 5, 0, 3, 7, 4},
  {2, 3, 5, 4, 1, 0, 6, 7},
  {3, 0, 4, 7, 2, 1, 5, 6},
  {1, 0, 2, 3, 6, 7, 5, 4},
  {0, 7, 3, 4, 1, 6, 2, 5},
  {7, 6, 4, 5, 0, 1, 3, 2},
  {6, 1, 5, 2, 7, 0, 4, 3},
  {5, 4, 6, 7, 2, 3, 1, 0},
  {4, 3, 7, 0, 5, 2, 6, 1},
  {3, 2, 0, 1, 4, 5, 7, 6},
  {2, 5, 1, 6, 3, 4, 0, 7}
};

/*! This function computes a Peano-Hilbert key for an integer triplet (x,y,z),
  *  with x,y,z in the range between 0 and 2^bits-1.
  */
peanokey peano_hilbert_key(int x, int y, int z, int bits)
{
  int mask;
  unsigned char rotation = 0;
  peanokey key = 0;

  for(mask = 1 << (bits - 1); mask > 0; mask >>= 1)
    {
      unsigned char pix = ((x & mask) ? 4 : 0) | ((y & mask) ? 2 : 0) | ((z & mask) ? 1 : 0);

      key <<= 3;
      key |= subpix3[rotation][pix];
      rotation = rottable3[rotation][pix];
    }

  return key;
}

peanokey morton_key(int x, int y, int z, int bits)
{
  int mask;
  peanokey morton = 0;

  for(mask = 1 << (bits - 1); mask > 0; mask >>= 1)
    {
      morton <<= 3;
      morton += ((z & mask) ? 4 : 0) + ((y & mask) ? 2 : 0) + ((x & mask) ? 1 : 0);
    }

  return morton;
}

peanokey peano_and_morton_key(int x, int y, int z, int bits, peanokey * morton_key)
{
  int mask;
  unsigned char rotation = 0;
  peanokey key = 0;
  peanokey morton = 0;

  for(mask = 1 << (bits - 1); mask > 0; mask >>= 1)
    {
      unsigned char pix = ((x & mask) ? 4 : 0) | ((y & mask) ? 2 : 0) | ((z & mask) ? 1 : 0);

      key <<= 3;
      key |= subpix3[rotation][pix];
      rotation = rottable3[rotation][pix];

      morton <<= 3;
      morton += ((z & mask) ? 4 : 0) + ((y & mask) ? 2 : 0) + ((x & mask) ? 1 : 0);
    }

  *morton_key = morton;

  return key;
}

static int quadrants[24][2][2][2] = {
  /* rotx=0, roty=0-3 */
  {{{0, 7}, {1, 6}}, {{3, 4}, {2, 5}}},
  {{{7, 4}, {6, 5}}, {{0, 3}, {1, 2}}},
  {{{4, 3}, {5, 2}}, {{7, 0}, {6, 1}}},
  {{{3, 0}, {2, 1}}, {{4, 7}, {5, 6}}},
  /* rotx=1, roty=0-3 */
  {{{1, 0}, {6, 7}}, {{2, 3}, {5, 4}}},
  {{{0, 3}, {7, 4}}, {{1, 2}, {6, 5}}},
  {{{3, 2}, {4, 5}}, {{0, 1}, {7, 6}}},
  {{{2, 1}, {5, 6}}, {{3, 0}, {4, 7}}},
  /* rotx=2, roty=0-3 */
  {{{6, 1}, {7, 0}}, {{5, 2}, {4, 3}}},
  {{{1, 2}, {0, 3}}, {{6, 5}, {7, 4}}},
  {{{2, 5}, {3, 4}}, {{1, 6}, {0, 7}}},
  {{{5, 6}, {4, 7}}, {{2, 1}, {3, 0}}},
  /* rotx=3, roty=0-3 */
  {{{7, 6}, {0, 1}}, {{4, 5}, {3, 2}}},
  {{{6, 5}, {1, 2}}, {{7, 4}, {0, 3}}},
  {{{5, 4}, {2, 3}}, {{6, 7}, {1, 0}}},
  {{{4, 7}, {3, 0}}, {{5, 6}, {2, 1}}},
  /* rotx=4, roty=0-3 */
  {{{6, 7}, {5, 4}}, {{1, 0}, {2, 3}}},
  {{{7, 0}, {4, 3}}, {{6, 1}, {5, 2}}},
  {{{0, 1}, {3, 2}}, {{7, 6}, {4, 5}}},
  {{{1, 6}, {2, 5}}, {{0, 7}, {3, 4}}},
  /* rotx=5, roty=0-3 */
  {{{2, 3}, {1, 0}}, {{5, 4}, {6, 7}}},
  {{{3, 4}, {0, 7}}, {{2, 5}, {1, 6}}},
  {{{4, 5}, {7, 6}}, {{3, 2}, {0, 1}}},
  {{{5, 2}, {6, 1}}, {{4, 3}, {7, 0}}}
};

static int rotxmap_table[24] = { 4, 5, 6, 7, 8, 9, 10, 11,
  12, 13, 14, 15, 0, 1, 2, 3, 17, 18, 19, 16, 23, 20, 21, 22
};

static int rotymap_table[24] = { 1, 2, 3, 0, 16, 17, 18, 19,
  11, 8, 9, 10, 22, 23, 20, 21, 14, 15, 12, 13, 4, 5, 6, 7
};

static int rotx_table[8] = { 3, 0, 0, 2, 2, 0, 0, 1 };
static int roty_table[8] = { 0, 1, 1, 2, 2, 3, 3, 0 };

static int sense_table[8] = { -1, -1, -1, +1, +1, -1, -1, -1 };

peanokey peano_hilbert_key_old(int x, int y, int z, int bits)
{
  int i, quad, bitx, bity, bitz;
  int mask, rotation, rotx, roty, sense;
  peanokey key;

  mask = 1 << (bits - 1);
  key = 0;
  rotation = 0;
  sense = 1;

  for(i = 0; i < bits; i++, mask >>= 1)
    {
      bitx = (x & mask) ? 1 : 0;
      bity = (y & mask) ? 1 : 0;
      bitz = (z & mask) ? 1 : 0;

      quad = quadrants[rotation][bitx][bity][bitz];

      key <<= 3;
      key += (sense == 1) ? (quad) : (7 - quad);

      rotx = rotx_table[quad];
      roty = roty_table[quad];
      sense *= sense_table[quad];

      while(rotx > 0)
	{
	  rotation = rotxmap_table[rotation];
	  rotx--;
	}

      while(roty > 0)
	{
	  rotation = rotymap_table[rotation];
	  roty--;
	}
    }

  return key;
}

peanokey peano_and_morton_key_old(int x, int y, int z, int bits, peanokey * morton_key)
{
  int i, quad, bitx, bity, bitz;
  int mask, rotation, rotx, roty, sense;
  peanokey key, morton;

  mask = 1 << (bits - 1);
  key = 0;
  rotation = 0;
  sense = 1;
  morton = 0;

  for(i = 0; i < bits; i++, mask >>= 1)
    {
      bitx = (x & mask) ? 1 : 0;
      bity = (y & mask) ? 1 : 0;
      bitz = (z & mask) ? 1 : 0;

      quad = quadrants[rotation][bitx][bity][bitz];

      key <<= 3;
      key += (sense == 1) ? (quad) : (7 - quad);

      rotx = rotx_table[quad];
      roty = roty_table[quad];
      sense *= sense_table[quad];

      while(rotx > 0)
	{
	  rotation = rotxmap_table[rotation];
	  rotx--;
	}

      while(roty > 0)
	{
	  rotation = rotymap_table[rotation];
	  roty--;
	}

      morton <<= 3;
      morton += (bitz << 2) + (bity << 1) + bitx;
    }

  *morton_key = morton;

  return key;
}

peanokey morton_key_old(int x, int y, int z, int bits)
{
  int i, bitx, bity, bitz;
  int mask;
  peanokey morton;

  mask = 1 << (bits - 1);
  morton = 0;

  for(i = 0; i < bits; i++, mask >>= 1)
    {
      bitx = (x & mask) ? 1 : 0;
      bity = (y & mask) ? 1 : 0;
      bitz = (z & mask) ? 1 : 0;

      morton <<= 3;
      morton += (bitz << 2) + (bity << 1) + bitx;
    }

  return morton;
}

static void msort_peano_with_tmp(struct peano_hilbert_data *b, size_t n, struct peano_hilbert_data *t)
{
  struct peano_hilbert_data *tmp;
  struct peano_hilbert_data *b1, *b2;
  size_t n1, n2;

  if(n <= 1)
    return;

  n1 = n / 2;
  n2 = n - n1;
  b1 = b;
  b2 = b + n1;

  msort_peano_with_tmp(b1, n1, t);
  msort_peano_with_tmp(b2, n2, t);

  tmp = t;

  while(n1 > 0 && n2 > 0)
    {
      if(b1->key <= b2->key)
	{
	  --n1;
	  *tmp++ = *b1++;
	}
      else
	{
	  --n2;
	  *tmp++ = *b2++;
	}
    }

  if(n1 > 0)
    memcpy(tmp, b1, n1 * sizeof(struct peano_hilbert_data));
  memcpy(b, t, (n - n2) * sizeof(struct peano_hilbert_data));
}

void mysort_peano(void *b, size_t n, size_t s, int (*cmp) (const void *, const void *))
{
  const size_t size = n * s;

  struct peano_hilbert_data *tmp = (struct peano_hilbert_data *) mymalloc(size);

  msort_peano_with_tmp((struct peano_hilbert_data *) b, n, tmp);

  myfree(tmp);
}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The lognormal distribution has the form

   p(x) dx = 1/(x * sqrt(2 pi sigma^2)) exp(-(ln(x) - zeta)^2/2 sigma^2) dx

   for x > 0. Lognormal random numbers are the exponentials of
   gaussian random numbers */

double
gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma)
{
  double u, v, r2, normal, z;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */

      u = -1 + 2 * gsl_rng_uniform (r);
      v = -1 + 2 * gsl_rng_uniform (r);

      /* see if it is in the unit circle */
      r2 = u * u + v * v;
    }
  while (r2 > 1.0 || r2 == 0);

  normal = u * sqrt (-2.0 * log (r2) / r2);

  z =  exp (sigma * normal + zeta);

  return z;
}

double
gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma)
{
  if (x <= 0)
    {
      return 0 ;
    }
  else
    {
      double u = (log (x) - zeta)/sigma;
      double p = 1 / (x * fabs(sigma) * sqrt (2 * M_PI)) * exp (-(u * u) /2);
      return p;
    }
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
 * Copyright (c) 2010-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include <lapacke.h>
#include "dplasma.h"
#include "dplasmatypes.h"

#include "map2.h"

static int
dplasma_zlacpy_trim_operator( parsec_execution_stream_t *es,
                              const parsec_tiled_matrix_dc_t *descA,
                              const parsec_tiled_matrix_dc_t *descB,
                              const void *_A, void *_B,
                              PLASMA_enum uplo, int m, int n,
                              void *args )
{
    int tempmm, tempnn, ldam, ldbm;
    const parsec_complex64_t *A = (const parsec_complex64_t*)_A;
    parsec_complex64_t       *B = (parsec_complex64_t*)_B;
    (void)es;
    (void)args;
    (void)uplo;

    tempmm = ((m)==((descB->mt)-1)) ? ((descB->m)-(m*(descB->mb))) : (descB->mb);
    tempnn = ((n)==((descB->nt)-1)) ? ((descB->n)-(n*(descB->nb))) : (descB->nb);
    ldam = BLKLDD( descA, m );
    ldbm = BLKLDD( descB, m );

    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR, lapack_const( PlasmaUpperLower ), tempmm, tempnn, A, ldam, B, ldbm);

    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy_New - Generates an object that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2_New() for further information.
 *
 * WARNING: The computations are not done by this call.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval NULL if incorrect parameters are given.
 *          \retval The PaRSEC object describing the operation that can be
 *          enqueued in the runtime. It, then, needs to be
 *          destroy with the corresponding _destruct function.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy_New
 * @sa dplasma_dlacpy_New
 * @sa dplasma_slacpy_New
 *
 ******************************************************************************/
parsec_taskpool_t*
dplasma_zlacpy_trim_New( PLASMA_enum uplo,
                         const parsec_tiled_matrix_dc_t *A,
                         parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t* tp;

    tp = dplasma_map2_New(uplo, PlasmaNoTrans, A, B,
                          dplasma_zlacpy_trim_operator, NULL );
    return tp;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 *  dplasma_zlacpy_Destruct - Free the data structure associated to an object
 *  created with dplasma_zlacpy_New().
 *
 *******************************************************************************
 *
 * @param[in,out] o
 *          On entry, the object to destroy.
 *          On exit, the object cannot be used anymore.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy
 *
 ******************************************************************************/
void
dplasma_zlacpy_trim_Destruct( parsec_taskpool_t *tp )
{
    parsec_taskpool_free(tp);
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy - Generates an object that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2() for further information.
 *
 *******************************************************************************
 *
 * @param[in,out] dague
 *          The dague context of the application that will run the operation.
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy
 * @sa dplasma_dlacpy
 * @sa dplasma_slacpy
 *
 ******************************************************************************/
int
dplasma_zlacpy_trim( parsec_context_t *parsec,
                     PLASMA_enum uplo,
                     const parsec_tiled_matrix_dc_t *A,
                     parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t *parsec_zlacpy_trim = NULL;

    if ((uplo != PlasmaUpperLower) &&
        (uplo != PlasmaUpper)      &&
        (uplo != PlasmaLower))
    {
        dplasma_error("dplasma_zlacpy_trim", "illegal value of uplo");
        return -2;
    }

    parsec_zlacpy_trim = dplasma_zlacpy_trim_New(uplo, A, B);

    if ( parsec_zlacpy_trim != NULL )
    {
        parsec_context_add_taskpool(parsec, parsec_zlacpy_trim);
        dplasma_wait_until_completion(parsec);
        dplasma_zlacpy_trim_Destruct( parsec_zlacpy_trim );
    }
    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
/***************************************************************************
 *   Copyright (C) 2008 by Regis Behmo,,,   *
 *   regis.behmo@ecp.fr   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 ***************************************************************************/

#ifndef _MPLSH_FIT_TUNE_H_INCLUDED_
#define _MPLSH_FIT_TUNE_H_INCLUDED_

#include <lshkit.h>
#include <lshkit/matrix.h>
#include <lshkit/tune.h>

#include <cstdlib>
#include <gsl/gsl_multifit.h>
#include <boost/program_options.hpp>
#include <boost/progress.hpp>

bool is_good_value (double v);
bool constraint (const lshkit::tune::Input &x);

double mplsh_tune(const double& param_M, const double& param_G,
                  const double& param_a_M, const double& param_b_M, const double& param_c_M,
                  const double& param_a_G, const double& param_b_G, const double& param_c_G,
                  int N, int K, int L, int T, int& M, float& W, float R );

template<class DATA, class METRIC>
void mplsh_fit( const DATA& data, const METRIC& metric, int NN, double& param_M, double& param_G,
                double& param_a_M, double& param_b_M, double& param_c_M,
                double& param_a_G, double& param_b_G, double& param_c_G )
{
	// Parameters
	unsigned N, P, Q, K, F;
	N = 0;// # points
	P = 50000;// number of pairs to sample
	Q = 1000;// number of queries to sample
	K = NN;// search for K nearest neighbors
	F = 10;// divide the sample to F folds

	std::vector<unsigned> idx(data.getSize());
	for (unsigned i = 0; i < idx.size(); ++i) idx[i] = i;
	random_shuffle(idx.begin(), idx.end());

	if (N > 0 && N < data.getSize()) idx.resize(N);

	lshkit::DefaultRng rng;
	rng.seed(0);//plant the same seed, always
	boost::variate_generator<lshkit::DefaultRng &, lshkit::UniformUnsigned> gen(rng, lshkit::UniformUnsigned(0, idx.size()-1));

	double gM = 0.0;
	double gG = 0.0;
	{
		// sample P pairs of points
		for (unsigned k = 0; k < P; ++k)
		{
			double dist, logdist;
			for (;;)
			{
				unsigned i = gen();
				unsigned j = gen();
				if (i == j) continue;
				dist = metric( data[idx[i]], data[idx[j]] );
				logdist = log(dist);
				if (is_good_value(logdist)) break;
			}
			gM += dist;
			gG += logdist;
		}
		gM /= P;
		gG /= P;
		gG = exp(gG);
	}

  // TODO check that
	// Custom
	Q = (idx.size() > 1000)? 1000 : idx.size();// Added

	// TODO check that
  // Supprime
  //if (Q > idx.size()) Q = idx.size();
  //if (K > idx.size() - Q) K = idx.size() - Q;

	/* sample query */
	std::vector<unsigned> qry(Q);

	lshkit::SampleQueries(&qry, idx.size(), rng);

	/* do the queries */
	std::vector<lshkit::Topk<unsigned> > topks(Q);
	for (unsigned i = 0; i < Q; ++i) topks[i].reset(K);

	/* ... */
	gsl_matrix *X = gsl_matrix_alloc(F * K, 3);
	gsl_vector *yM = gsl_vector_alloc(F * K);
	gsl_vector *yG = gsl_vector_alloc(F * K);
	gsl_vector *pM = gsl_vector_alloc(3);
	gsl_vector *pG = gsl_vector_alloc(3);
	gsl_matrix *cov = gsl_matrix_alloc(3,3);

	std::vector<double> M(K);
	std::vector<double> G(K);

	boost::progress_display progress(F, std::cerr);
	unsigned m = 0;
	for (unsigned l = 0; l < F; l++)
	{
		// Scan
		for (unsigned i = l; i< idx.size(); i += F)
		{
			for (unsigned j = 0; j < Q; j++)
			{
				int id = qry[j];
				if (i != id)
				{
					float d = metric( data[idx[id]], data[idx[i]] );
					if (is_good_value(log(double(d)))) topks[j] << lshkit::Topk<unsigned>::Element(i, d);
				}
			}
		}

		std::fill(M.begin(), M.end(), 0.0);
		std::fill(G.begin(), G.end(), 0.0);

		for (unsigned i = 0; i < Q; i++)
		{
			for (unsigned k = 0; k < K; k++)
			{
				M[k] += topks[i][k].dist;
				G[k] += log(topks[i][k].dist);
			}
		}

		for (unsigned k = 0; k < K; k++)
		{
			M[k] = log(M[k]/Q);
			G[k] /= Q;
			gsl_matrix_set(X, m, 0, 1.0);
			gsl_matrix_set(X, m, 1, log(double(data.getSize() * (l + 1)) / double(F)));
			gsl_matrix_set(X, m, 2, log(double(k + 1)));
			gsl_vector_set(yM, m, M[k]);
			gsl_vector_set(yG, m, G[k]);
			++m;
		}
        ++progress;
	}

	gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(F * K, 3);

	double chisq;

	gsl_multifit_linear(X, yM, pM, cov, &chisq, work);
	gsl_multifit_linear(X, yG, pG, cov, &chisq, work);

	param_M = gM;
	param_G = gG;
	param_a_M = gsl_vector_get(pM, 0);
	param_b_M = gsl_vector_get(pM, 1);
	param_c_M = gsl_vector_get(pM, 2);
	param_a_G = gsl_vector_get(pG, 0);
	param_b_G = gsl_vector_get(pG, 1);
	param_c_G = gsl_vector_get(pG, 2);

	gsl_matrix_free(X);
	gsl_matrix_free(cov);
	gsl_vector_free(pM);
	gsl_vector_free(pG);
	gsl_vector_free(yM);
	gsl_vector_free(yG);
}

template<class DATA, class METRIC>
double mplsh_fit_tune( const DATA& data, const METRIC& metric, const int& L, const int& T, int& M, float& W, float R = 0.8, int K = 1 )
{

	/** ***************************/
	/** ********** FIT  ***********/
	/** ***************************/
  std::cout << "Fitting..." << std::endl;
	double param_M, param_G, param_a_M, param_b_M, param_c_M, param_a_G, param_b_G, param_c_G;
	mplsh_fit( data, metric, K, param_M, param_G, param_a_M, param_b_M, param_c_M, param_a_G, param_b_G, param_c_G );

  std::cout << param_M << "    " << param_G << std::endl;
  std::cout << param_a_M << "    " << param_b_M << "    " << param_c_M << std::endl;
  std::cout << param_a_G << "    " << param_b_G << "    " << param_c_G << std::endl;

	/** ***************************/
	/** ********** TUNE ***********/
	/** ***************************/
  std::cout << "Tuning..." << std::endl;
	double cost = mplsh_tune(param_M, param_G, param_a_M, param_b_M, param_c_M, param_a_G, param_b_G, param_c_G,
                    data.getSize(), K, L, T, M, W, R );
  std::cout << "L = " << L;
  std::cout << "    " << "T = " << T;
  std::cout << "    " << "M = " << M;
  std::cout << "    " << "W = " << W;
  std::cout << "    " << "R = " << R;
  std::cout << "    " << "M = " << M;
  std::cout << std::endl;

  return cost;
}

#endif

ext: h

max_stars_repo_name: ocallaco/LuaSHkit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: ocallaco/LuaSHkit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: ocallaco/LuaSHkit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
/*
 *  Copyright 2017 NVIDIA Corporation
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

/* include the ATLAS headers */

extern "C"
{
#include <cblas.h>
}

/* choose precision to train and classify.  Only float and double are
 * currently suppored
 */

typedef float floatType_t;

/* macro to convert 2d coords to 1d offset */

#define INDX(row,col,ld) (((col) * (ld)) + (row))

inline float sigmoid_f( float z )
{
  return 1.0f / ( 1.0f + expf( -z ) );
} /* end sigmoid_f */

inline double sigmoid( double z )
{
  return 1.0 / ( 1.0 + exp( -z ) );
} /* end sigmoid */

inline float sigmoidGradient_f( float z )
{
  float temp = sigmoid_f( z );
  return temp * ( 1.0f - temp );
} /* end sigGrad_f */

inline double sigmoidGradient( double z )
{
  double temp = sigmoid( z );
  return temp * ( 1.0 - temp );
} /* end sigGrad_f */

/* hardcoded constants for training and test set size and feature
 * vector size
 */

#define FEATURE_VECTOR_SIZE (784)
#define TRAINING_SET_SIZE (60000)
#define TEST_SET_SIZE (10000)
#define HIDDEN_LAYER_SIZE (25)
#define NUM_OUTPUT_CLASSES (10)

/* CUDA debugging */

#ifdef DEBUG
#define CUDA_CALL(F)  if( (F) != cudaSuccess ) \
  {printf("Error %s at %s:%d\n", cudaGetErrorString(cudaGetLastError()), \
   __FILE__,__LINE__); exit(-1);}
#define CUDA_CHECK()  if( (cudaPeekAtLastError()) != cudaSuccess ) \
  {printf("Error %s at %s:%d\n", cudaGetErrorString(cudaGetLastError()), \
   __FILE__,__LINE__-1); exit(-1);}
#else
#define CUDA_CALL(F) (F)
#define CUDA_CHECK()
#endif

/* function defs */

void readMatrixFromFile( char *, float *, const int, const int, const int );
void readCommandLineArgs( int, char *[], float *, int *, int *, int *);

void calcGradient( floatType_t *X,
                   int const XRows,
                   int const XCols,
                   floatType_t const *theta1,
                   int         const theta1Rows,
                   int         const theta1Cols,
                   floatType_t const *theta2,
                   int         const theta2Rows,
                   int         const theta2Cols,
                   floatType_t const *Y,
                   floatType_t *cost,
                   floatType_t       *theta1Grad,
                   floatType_t       *theta2Grad,
                   floatType_t       *tempMatrix );
void predict( floatType_t *X,
                   int const XRows,
                   int const XCols,
                   floatType_t const *theta1,
                   int         const theta1Rows,
                   int         const theta1Cols,
                   floatType_t const *theta2,
                   int         const theta2Rows,
                   int         const theta2Cols,
                   int               *predictVector);
void trainNetwork( floatType_t       *X,
                   int         const XRows,
                   int         const XCols,
                   floatType_t       *theta1,
                   int         const theta1Rows,
                   int         const theta1Cols,
                   floatType_t       *theta2,
                   int         const theta2Rows,
                   int         const theta2Cols,
                   floatType_t const *Y,
                   float       const learningRate,
                   int         const iterations,
                   int         const batchSize );

max_stars_repo_name: jnbntz/gpu-edu-workshops


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

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.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``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 Javier Burguete Tolosa 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.
*/

/**
 * \file utils.c
 * \brief Source file with util functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"

gchar *error_message = NULL;    ///< error message string.

/**
 * Function to print an error message.
 */
void
show_error (const char *message)        ///< error message string.
{
  printf ("%s\n%s\n", _("ERROR!"), message);
}

/**
 * Function to get an integer number of a XML node property.
 *
 * \return Integer number value.
 */
int
xml_node_get_int (xmlNode * node,       ///< XML node.
                  const xmlChar * prop, ///< XML property.
                  int *error_code)      ///< Error code.
{
  int i = 0;
  xmlChar *buffer;
  buffer = xmlGetProp (node, prop);
  if (!buffer)
    *error_code = 1;
  else
    {
      if (sscanf ((char *) buffer, "%d", &i) != 1)
        *error_code = 2;
      else
        *error_code = 0;
      xmlFree (buffer);
    }
  return i;
}

/**
 * Function to get an unsigned integer number of a XML node property.
 *
 * \return Unsigned integer number value.
 */
unsigned int
xml_node_get_uint (xmlNode * node,      ///< XML node.
                   const xmlChar * prop,        ///< XML property.
                   int *error_code)     ///< Error code.
{
  unsigned int i = 0;
  xmlChar *buffer;
  buffer = xmlGetProp (node, prop);
  if (!buffer)
    *error_code = 1;
  else
    {
      if (sscanf ((char *) buffer, "%u", &i) != 1)
        *error_code = 2;
      else
        *error_code = 0;
      xmlFree (buffer);
    }
  return i;
}

/**
 * Function to get an unsigned integer number of a XML node property with a
 *   default value.
 *
 * \return Unsigned integer number value.
 */
unsigned int
xml_node_get_uint_with_default (xmlNode * node, ///< XML node.
                                const xmlChar * prop,   ///< XML property.
                                unsigned int default_value,
                                ///< default value.
                                int *error_code)        ///< Error code.
{
  unsigned int i;
  if (xmlHasProp (node, prop))
    i = xml_node_get_uint (node, prop, error_code);
  else
    {
      i = default_value;
      *error_code = 0;
    }
  return i;
}

/**
 * Function to get an long double number of a XML node property.
 *
 * \return Long double number value.
 */
long double
xml_node_get_float (xmlNode * node,     ///< XML node.
                    const xmlChar * prop,       ///< XML property.
                    int *error_code)    ///< Error code.
{
  long double x = 0.L;
  xmlChar *buffer;
  buffer = xmlGetProp (node, prop);
  if (!buffer)
    *error_code = 1;
  else
    {
      if (sscanf ((char *) buffer, "%Lg", &x) != 1)
        *error_code = 2;
      else
        *error_code = 0;
      xmlFree (buffer);
    }
  return x;
}

/**
 * Function to read the minimum, the interval and the random function type of a
 * variable on a XML node.
 *
 * \return 1 on success, 0 on error.
 */
int
read_variable (xmlNode * node,  ///< XML node.
               long double *minimum,    ///< minimum value.
               long double *interval,   ///< interval value.
               unsigned int *type,      ///< random function type.
               unsigned int n)  ///< variable number.
{
  char number[16];
  gchar *buffer;
  xmlChar *prop;
  int code;

  if (!node)
    {
      error_message = g_strdup (_("No XML node"));
      goto exit_on_error;
    }
  if (xmlStrcmp (node->name, XML_VARIABLE))
    {
      error_message = g_strdup (_("Bad XML node"));
      goto exit_on_error;
    }
  minimum[n] = xml_node_get_float (node, XML_MINIMUM, &code);
  if (code)
    {
      error_message = g_strdup (_("Bad minimum"));
      goto exit_on_error;
    }
  interval[n] = xml_node_get_float (node, XML_INTERVAL, &code);
  if (code)
    {
      error_message = g_strdup (_("Bad interval"));
      goto exit_on_error;
    }
  prop = xmlGetProp (node, XML_TYPE);
  if (!prop || !xmlStrcmp (prop, XML_RANDOM))
    type[n] = RANDOM_TYPE_UNIFORM;
  else if (!xmlStrcmp (prop, XML_BOTTOM))
    type[n] = RANDOM_TYPE_BOTTOM;
  else if (!xmlStrcmp (prop, XML_EXTREME))
    type[n] = RANDOM_TYPE_EXTREME;
  else if (!xmlStrcmp (prop, XML_TOP))
    type[n] = RANDOM_TYPE_TOP;
  else if (!xmlStrcmp (prop, XML_REGULAR))
    type[n] = RANDOM_TYPE_REGULAR;
  else if (!xmlStrcmp (prop, XML_ORTHOGONAL))
    type[n] = RANDOM_TYPE_ORTHOGONAL;
  else
    {
      xmlFree (prop);
      error_message = g_strdup (_("Bad random type function"));
      goto exit_on_error;
    }
  xmlFree (prop);
  return 1;

exit_on_error:
  snprintf (number, 16, "%u", n + 1);
  buffer = error_message;
  error_message
    = g_strconcat (_("Variable"), " ", number, ":\n", error_message, NULL);
  g_free (buffer);
  return 0;
}

/**
 * Function to set the precision in a maxima program file.
 */
void
print_maxima_precision (FILE * file)    ///< maxima program file.
{
  fprintf (file, "fpprec:%u;\nfpprintprec:fpprec;\n", MAXIMA_PRECISION);
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
// -*- mode:C++; tab-width:4; c-basic-offset:4; indent-tabs-mode:nil -*-

/*
 * Copyright (C) 2009 RobotCub Consortium
 * Author: Alessandro Scalzo alessandro.scalzo@iit.it
 * CopyPolicy: Released under the terms of the GNU GPL v2.0.
 *
 */

#include <gsl/gsl_math.h>
#include <memory.h>

#include <stdio.h>

#ifndef __ALE_TOUCHSENSOR_H__
#define __ALE_TOUCHSENSOR_H__

#define MAX_TAXELS 128

class TouchSensor
{
    bool calibrated_skin;

public:

    void setColor(unsigned char r, unsigned char g, unsigned char b)
    {
        R_MAX = r;
        G_MAX = g;
        B_MAX = b;
    }

    void setCalibrationFlag (bool use_calibrated_skin)
    {
        calibrated_skin=use_calibrated_skin;
    }

    void resize(int width,int height,int margin)
    {
        if (3*margin>=width || 3*margin>=height) margin=0;

        double scaleX=double(width -2*margin)/(dXmax-dXmin);
        double scaleY=double(height-2*margin)/(dYmax-dYmin);
        double scale=scaleX<scaleY?scaleX:scaleY;

        m_Radius=m_RadiusOrig*scale;

        m_maxRangeLight=int(m_Radius);

        double dXmid=0.5*(dXmin+dXmax);
        double dYmid=0.5*(dYmin+dYmax);

        int w2=width /2;
        int h2=height/2;

        for (int i=0; i<nTaxels; ++i)
        {
            x[i]=w2+int(scale*(dX[i]-dXmid));
            y[i]=h2+int(scale*(dY[i]-dYmid));
        }

        for (int i=0; i<nVerts; ++i)
        {
            xv[i]=w2+int(scale*(dXv[i]-dXmid));
            yv[i]=h2+int(scale*(dYv[i]-dYmid));
        }

        double sigma=0.5*5.55*scale;
        int maxRange=int(2.5*sigma);

        if (maxRange!=m_maxRange)
        {
            m_maxRange=maxRange;

            delete [] Exponential;
            Exponential=new double[maxRange];

            double k=-0.5/(sigma*sigma);
            for (int x=0; x<maxRange; ++x)
            {
                Exponential[x]=exp(k*double(x*x));
            }
        }

        xMin=w2+int(scale*(dXc-dXmid-15.0))-maxRange;
        xMax=w2+int(scale*(dXc-dXmid+15.0))+maxRange;
        yMin=h2+int(scale*(dYc-dYmid-15.0))-maxRange;
        yMax=h2+int(scale*(dYc-dYmid+15.0))+maxRange;

        if (xMin<0)      xMin=0;
        if (xMax>width)  xMax=width;
        if (yMin<0)      yMin=0;
        if (yMax>height) yMax=height;

        m_Width=width;
        m_Height=height;
    }

    virtual ~TouchSensor()
    {
        if (Exponential)
        {
            delete [] Exponential;
            Exponential=0;
        }
    }

    int Abs(int x)
    {
        return x>=0?x:-x;
    }

    int get_nTaxels ()
    {
        return nTaxels;
    }

    void eval_light(unsigned char *image)
    {
        int act;
        int dx,dy;
        int Y0,Y1;
        int dya,dyb,dxa,dxb;

        double remapped_activation[MAX_TAXELS];
        switch (ilayoutNum)
        {
            case 0:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[i]=activation[i];
                break;
            case 1:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[nTaxels-1-i]=activation[i];
                break;
            default:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[i]=activation[i];
                printf("WARN: unkwnown layout number.\n");
                break;
        }

        int maxRange2=m_maxRangeLight*m_maxRangeLight;

        for (int i=0; i<nTaxels; ++i) if (remapped_activation[i]>0.0)
        {
            act=int(dGain*remapped_activation[i]);
            Y0=(m_Height-y[i]-1)*m_Width+x[i];

            dya=(y[i]>=m_maxRangeLight)?-m_maxRangeLight:-y[i];
            dyb=(y[i]+m_maxRangeLight<m_Height)?m_maxRangeLight:m_Height-y[i]-1;

            dxa=(x[i]>=m_maxRangeLight)?-m_maxRangeLight:-x[i];
            dxb=(x[i]+m_maxRangeLight<m_Width)?m_maxRangeLight:m_Width-x[i]-1;

            for (dy=dya; dy<=dyb; ++dy)
            {
                Y1=Y0-dy*m_Width;

                for (dx=dxa; dx<=dxb; ++dx)
                {
                    if (dx*dx+dy*dy<=maxRange2)
                    {
                        image[(dx+Y1)*3]=act<255?act:255;
                    }
                }
            }
        }
    }

    void eval(unsigned char *image)
    {
        int act;
        int dx,dy;
        int Y0,Y1;
        int index;
        double k0,k1;
        int dya,dyb,dxa,dxb;
        double remapped_activation[MAX_TAXELS];
        switch (ilayoutNum)
        {
            case 0:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[i]=activation[i];
                break;
            case 1:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[nTaxels-1-i]=activation[i];
                break;
            default:
                for (int i=0; i<nTaxels; ++i)  remapped_activation[i]=activation[i];
                printf("WARN: unkwnown layout number.\n");
                break;
        }

        for (int i=0; i<nTaxels; ++i) if (remapped_activation[i]>0.0)
        {
            k0=dGain*remapped_activation[i];
            Y0=(m_Height-y[i]-1)*m_Width+x[i];

            dya=(y[i]>=m_maxRange)?-m_maxRange:-y[i];
            dyb=(y[i]+m_maxRange<m_Height)?m_maxRange:m_Height-y[i]-1;

            dxa=(x[i]>=m_maxRange)?-m_maxRange:-x[i];
            dxb=(x[i]+m_maxRange<m_Width)?m_maxRange:m_Width-x[i]-1;

            for (dy=dya; dy<=dyb; ++dy)
            {
                k1=k0*Exponential[Abs(dy)];
                Y1=Y0-dy*m_Width;

                for (dx=dxa; dx<=dxb; ++dx)
                {
                    index=(dx+Y1)*3;

                    if (image[index]<R_MAX || image[index+1]<G_MAX || image[index+2]<B_MAX)
                    {
                        act=int(k1*Exponential[Abs(dx)]);

                        int actR=image[index  ]+(act*R_MAX)/255;
                        int actG=image[index+1]+(act*G_MAX)/255;
                        int actB=image[index+2]+(act*B_MAX)/255;

                        image[index  ]=actR<R_MAX?actR:R_MAX;
                        image[index+1]=actG<G_MAX?actG:G_MAX;
                        image[index+2]=actB<B_MAX?actB:B_MAX;
                    }
                }
            }
        }
    }

    void setActivationFirst7(unsigned char* data)
    {
        for (int i=0; i<7; ++i)
        {
            activation[i]=data[i+1]<=244?double(244-data[i+1]):0.0;
        }
    }
    void setActivationLast5(unsigned char* data)
    {
        for (int i=1; i<=5; ++i)
        {
            activation[i+6]=data[i]<=244?double(244-data[i]):0.0;
        }
    }

    void setActivationFromPortData(double val, int id)
    {
        if (calibrated_skin)
        {
            if (val>244.0)
            {
                activation[id]=244.0;
            }
            else if (val<0.0)
            {
                activation[id]=0.0;
            }
            else
            {
                activation[id]=val;
            }
        }
        else
        {
            activation[id]=val<=244?double(244-val):0.0;
        }
    }

    virtual void draw(unsigned char *image)
    {
        for (int i=0; i<nVerts; ++i)
        {
            drawLine(image,xv[i],yv[i],xv[(i+1)%nVerts],yv[(i+1)%nVerts]);
        }

        for (int i=0; i<nTaxels; ++i)
        {
            drawCircle(image,x[i],y[i],m_Radius);
        }
    }

protected:
    void dither(int x,int y,unsigned char *image)
    {
        static const unsigned char R1=0x80,G1=0x50,B1=0x00;
        static const unsigned char R2=3*R1/4,G2=3*G1/4,B2=3*B1/4;
        static const unsigned char R4=3*R2/4,G4=3*G2/4,B4=3*B2/4;

        //y=m_Width-y-1;
        //int bytePos=(x+(y-1)*m_Width)*3;
        int bytePos=(x+(m_Height-y-2)*m_Width)*3;

        if (image[bytePos-3]<R4) image[bytePos-3]=R4;
        if (image[bytePos-2]<G4) image[bytePos-2]=G4;
        //if (image[bytePos-1]<B4) image[bytePos-1]=B4;

        if (image[bytePos  ]<R2) image[bytePos  ]=R2;
        if (image[bytePos+1]<G2) image[bytePos+1]=G2;
        //if (image[bytePos+2]<B2) image[bytePos+2]=B2;

        if (image[bytePos+3]<R4) image[bytePos+3]=R4;
        if (image[bytePos+4]<G4) image[bytePos+4]=G4;
        //if (image[bytePos+5]<B4) image[bytePos+5]=B4;

        bytePos+=m_Width*3;

        if (image[bytePos-3]<R2) image[bytePos-3]=R2;
        if (image[bytePos-2]<G2) image[bytePos-2]=G2;
        //if (image[bytePos-1]<B2) image[bytePos-1]=B2;

        if (image[bytePos  ]<R1) image[bytePos  ]=R1;
        if (image[bytePos+1]<G1) image[bytePos+1]=G1;
        //if (image[bytePos+2]<B1) image[bytePos+2]=B1;

        if (image[bytePos+3]<R2) image[bytePos+3]=R2;
        if (image[bytePos+4]<G2) image[bytePos+4]=G2;
        //if (image[bytePos+5]<B2) image[bytePos+5]=B2;

        bytePos+=m_Width*3;

        if (image[bytePos-3]<R4) image[bytePos-3]=R4;
        if (image[bytePos-2]<G4) image[bytePos-2]=G4;
        //if (image[bytePos-1]<B4) image[bytePos-1]=B4;

        if (image[bytePos  ]<R2) image[bytePos  ]=R2;
        if (image[bytePos+1]<G2) image[bytePos+1]=G2;
        //if (image[bytePos+2]<B2) image[bytePos+2]=B2;

        if (image[bytePos+3]<R4) image[bytePos+3]=R4;
        if (image[bytePos+4]<G4) image[bytePos+4]=G4;
        //if (image[bytePos+5]<B4) image[bytePos+5]=B4;
    }

    void drawLine(unsigned char *image,int x0,int y0,int x1,int y1)
    {
        if (x1==x0 && y1==y0) return;
        double Vx=double(x1-x0);
        double Vy=double(y1-y0);
        double dt=1.0/sqrt(Vx*Vx+Vy*Vy);

        for (double t=0.0; t<=1.0; t+=dt)
        {
            dither(x0+int(t*Vx),y0+int(t*Vy),image);
        }
    }

    void drawCircle(unsigned char *image,int cx,int cy,double radius)
    {
        double dt=1.0/(2*M_PI*radius);

        int dx,dy;

        double cs,sn;

        for (double t=0.0; t<=M_PI_4; t+=dt)
        {
            cs=cos(t);
            sn=sin(t);

            dx=int(radius*cs);
            dy=int(radius*sn);

            dither(cx+dx,cy+dy,image);
            dither(cx+dx,cy-dy,image);
            dither(cx-dx,cy-dy,image);
            dither(cx-dx,cy+dy,image);

            dither(cx+dy,cy+dx,image);
            dither(cx+dy,cy-dx,image);
            dither(cx-dy,cy-dx,image);
            dither(cx-dy,cy+dx,image);
        }
    }

    // original
    double dX[MAX_TAXELS],dY[MAX_TAXELS];
    static double dXmin,dXmax,dYmin,dYmax;
    double dXv[8],dYv[8];
    double dXc,dYc;
    double dGain;
    int    ilayoutNum;
    int    ilrMirror;

    double m_Radius,m_RadiusOrig;
    double activation[MAX_TAXELS];

    unsigned char R_MAX, G_MAX, B_MAX;

    int m_maxRangeLight;
    static int m_maxRange;
    static double *Exponential;

    // scaled
    int x[MAX_TAXELS],y[MAX_TAXELS];
    int xv[8],yv[8];

    int nVerts;
    int nTaxels;

    int xMin,xMax,yMin,yMax;

    int m_Width,m_Height;

    public:
    int min_tax;
    int max_tax;
};

#endif

max_stars_repo_name: robotology-legacy/icub-gtk-modules

max_issues_repo_name: robotology-legacy/icub-gtk-modules

max_forks_repo_name: robotology-legacy/icub-gtk-modules


text:
#ifndef OPS_REDUCERS_H
#define OPS_REDUCERS_H

#include <fundamental/tensor.h>
#include <nn/container.h>

#ifdef TENSOR_USE_CUBLAS
  #include <cublas_v2.h>
#endif

#ifdef TENSOR_USE_OPENBLAS
  #include <cblas.h>
#endif

#include <ops/transform.h>

namespace operators{

  template<typename xpu, index_t stream_id, typename T>
  void asum_backward(nn::TensorContainer<xpu, stream_id, T>& x,
                     nn::TensorContainer<xpu, stream_id, T>& y,
                     int offset, int inc, int n){

     auto bwd = [](nn::TensorContainer<xpu, stream_id, float> x,
                   nn::TensorContainer<xpu, stream_id, float> y,
                   int offset, int inc, int n)
     {

       auto strided_x = expression::transform::stride(x, offset, inc, n);

       core::Shape shape = strided_x.self().shape();

       //auto repeated_y = expression::transform::repeat(shape.size(), yexp);
       core::Constant<xpu, T, stream_id> repeated_y(y.gradient().scalar(),
                                                    core::Shape(x.shape().size()));

       expression::Exceturer<xpu>::backward(repeated_y, strided_x);

     };

     STREAM_BACKWARD(xpu, stream_id).put(bwd, x, y, offset, inc, n);

   }

  template<index_t stream_id>
  void asum(nn::TensorContainer<cpu, stream_id, float>& x,
            nn::TensorContainer<cpu, stream_id, float>& y,
            int offset, int inc, int n){

    //y.conditional_override(core::Shape(1));

    ASSERT_ERROR( ((n+offset) - (x.data().shape().size())) <= 0 , " Out of bound ");

    auto fwd = [](nn::TensorContainer<cpu, stream_id, float> x,
                  nn::TensorContainer<cpu, stream_id, float> y,
                  int offset, int inc, int n)
    {
      float* ptr_offseted = x.data().ptr() + offset;
      float result = cblas_sasum(n, ptr_offseted, inc);
      y.data().set(0, result);
    };

    STREAM_FORWARD(cpu, stream_id).put(fwd, x, y, offset, inc, n);

    if (x.require_grad() && y.require_grad())
      asum_backward(x, y, offset, inc, n);

  }

  #ifdef TENSOR_USE_CUBLAS

    template<index_t stream_id>
    void asum(nn::TensorContainer<gpu, stream_id, float>& x,
              nn::TensorContainer<gpu, stream_id, float>& y,
              int offset, int inc, int n){

      DEBUG_ASSERT(y.data().shape().size() == 1);
      ASSERT_ERROR( ((n+offset) - (x.data().shape().size())) <= 0 , " Out of bound ");

      auto fwd = [](nn::TensorContainer<gpu, stream_id, float> x,
                    nn::TensorContainer<gpu, stream_id, float> y,
                    int offset, int inc, int n)
      {
        auto cublasHandle = STREAM_FORWARD(gpu, stream_id).cublasHandle();
        float result;
        float* ptr_offseted = x.data().ptr() + offset;
        CUDBLAS_CHECK_ERROR_ENFORCE(cublasSasum(cublasHandle, n,
                                                ptr_offseted, inc,
                                                &result));
        y.data().set(0, result);

      };

      STREAM_FORWARD(gpu, stream_id).put(fwd, x, y, offset, inc, n);

      if (x.require_grad() && y.require_grad())
        asum_backward(x, y, offset, inc, n);

    }

  #endif

}

#ifdef TENSOR_USE_CUDNN

  namespace cudnn {

    template<cudnnReduceTensorOp_t OPERATION, index_t stream_id>
    void reduce_cudnn_i(core::Tensor<gpu, stream_id, float>& tensorA,
                        core::Tensor<gpu, stream_id, float>& tensorB,
                        core::Tensor<gpu, stream_id, indices_t>& indicies,
                        float alpha = 1, float beta = 0){

      DEBUG_ASSERT( indicies.shape() == tensorB.shape() );

      cudnnReduceTensorDescriptor_t reduce;
      cudnnCreateReduceTensorDescriptor(&reduce);

      CHECK_CUDNN( cudnnSetReduceTensorDescriptor(
                   reduce,
                   OPERATION,
                   CUDNN_DATA_FLOAT,
                   CUDNN_NOT_PROPAGATE_NAN,
                   CUDNN_REDUCE_TENSOR_FLATTENED_INDICES,
                   CUDNN_32BIT_INDICES) );

      size_t workspaces_sz;

      auto cudnnhandle = STREAM_FORWARD(gpu, stream_id).cudnnHandle();

      CHECK_CUDNN( cudnnGetReductionWorkspaceSize(cudnnhandle,
                                                  reduce,
                                                  tensorA.descriptor(),
                                                  tensorB.descriptor(),
                                                  &workspaces_sz) );

      auto workspace_ptr = cuda::generate_workspace<void>(workspaces_sz);

      CHECK_CUDNN( cudnnReduceTensor(
                                     cudnnhandle,
                                     reduce, // reduceTensorDesc
                                     indicies.ptr(), // indicies
                                     sizeof(indices_t) * indicies.shape().size(),  // indicesSizeInBytes,
                                     workspace_ptr, // workspace,
                                     workspaces_sz, //workspaceSizeInBytes,
                                     (void*) &alpha,
                                     tensorA.descriptor(),
                                     (void*) tensorA.ptr(),
                                     (void*) &beta,
                                     tensorB.descriptor(),
                                     (void*) tensorB.ptr()) ) ;

      for (int i = 0; i < indicies.shape()[0];  ++i)
        printf(" %f\n", tensorB.at(i) );

    }

    template<cudnnReduceTensorOp_t OPERATION, index_t stream_id>
    void reduce_cudnn(core::Tensor<gpu, stream_id, float>& tensorA,
                      core::Tensor<gpu, stream_id, float>& tensorB,
                      float alpha = 1, float beta = 0){

      cudnnReduceTensorDescriptor_t reduce;
      cudnnCreateReduceTensorDescriptor(&reduce);

      CHECK_CUDNN( cudnnSetReduceTensorDescriptor(
                   reduce,
                   OPERATION,
                   CUDNN_DATA_FLOAT,
                   CUDNN_NOT_PROPAGATE_NAN,
                   CUDNN_REDUCE_TENSOR_NO_INDICES,
                   CUDNN_32BIT_INDICES) );

      size_t workspaces_sz;

      auto cudnnhandle = STREAM_FORWARD(gpu, stream_id).cudnnHandle();

      CHECK_CUDNN( cudnnGetReductionWorkspaceSize(cudnnhandle,
                                                  reduce,
                                                  tensorA.descriptor(),
                                                  tensorB.descriptor(),
                                                  &workspaces_sz) );

      auto workspace_ptr = cuda::generate_workspace<void>(workspaces_sz);

      float dummy;

      CHECK_CUDNN( cudnnReduceTensor(
                                     cudnnhandle,
                                     reduce, // reduceTensorDesc
                                     &dummy, // indicies
                                     0,  // indicesSizeInBytes,
                                     workspace_ptr, // workspace,
                                     workspaces_sz, //workspaceSizeInBytes,
                                     (void*) &alpha,
                                     tensorA.descriptor(),
                                     (void*) tensorA.ptr(),
                                     (void*) &beta,
                                     tensorB.descriptor(),
                                     (void*) tensorB.ptr()) ) ;

    }

  }

#endif

namespace operators{

  #ifdef TENSOR_USE_CUDNN

    template<index_t stream_id>
    void reduce_sum(core::Tensor<gpu, stream_id, float>& tensorA,
                    core::Tensor<gpu, stream_id, float>& tensorB,
                    float alpha = 1, float beta = 0){
      cudnn::reduce_cudnn<CUDNN_REDUCE_TENSOR_ADD, stream_id>
              (tensorA, tensorB, alpha, beta);
    }

    template<index_t stream_id>
    void reduce_max(core::Tensor<gpu, stream_id, float>& tensorA,
                    core::Tensor<gpu, stream_id, float>& tensorB,
                    float alpha = 1, float beta = 0){
      cudnn::reduce_cudnn<CUDNN_REDUCE_TENSOR_MAX, stream_id>
              (tensorA, tensorB, alpha, beta);
    }

    template<index_t stream_id>
    void reduce_min(core::Tensor<gpu, stream_id, float>& tensorA,
                    core::Tensor<gpu, stream_id, float>& tensorB,
                    float alpha = 1, float beta = 0){
      cudnn::reduce_cudnn<CUDNN_REDUCE_TENSOR_MIN, stream_id>
              (tensorA, tensorB, alpha, beta);
    }

    template<index_t stream_id>
    void reduce_max_i(core::Tensor<gpu, stream_id, float>& tensorA,
                      core::Tensor<gpu, stream_id, float>& tensorB,
                      core::Tensor<gpu, stream_id, indices_t>& indices,
                      float alpha = 1, float beta = 0){
      cudnn::reduce_cudnn_i<CUDNN_REDUCE_TENSOR_MAX, stream_id>
              (tensorA, tensorB, indices, alpha, beta);
    }

    template<index_t stream_id>
    void reduce_min_i(core::Tensor<gpu, stream_id, float>& tensorA,
                      core::Tensor<gpu, stream_id, float>& tensorB,
                      core::Tensor<gpu, stream_id, indices_t>& indices,
                      float alpha = 1, float beta = 0){
      cudnn::reduce_cudnn_i<CUDNN_REDUCE_TENSOR_MIN, stream_id>
              (tensorA, tensorB, indices, alpha, beta);
    }

    template<index_t stream_id>
    TENSOR_INLINE_HOST void reduce_min(
        std::shared_ptr<nn::TensorContainer<gpu, stream_id, float>> A,
        std::shared_ptr<nn::TensorContainer<gpu, stream_id, float>> B,
        index_t axis){

        B->conditional_override(A->shape().reduce_axis(axis));

        if (B->require_grad()){
          // Create tensor for incicies
          core::Tensor<gpu, stream_id, float> I(B->shape()); I.allocate();

          auto fwd = [](core::Tensor<gpu, stream_id, float> A,
                        core::Tensor<gpu, stream_id, float> B,
                        core::Tensor<gpu, stream_id, float> I){
            reduce_min_i(A, B, I);
          };

          STREAM_FORWARD(gpu, stream_id).put(fwd, A->data(), B->data(), I);

        }else{

          auto fwd = [](core::Tensor<gpu, stream_id, float> A,
                        core::Tensor<gpu, stream_id, float> B){
            reduce_min(A, B);
          };

          STREAM_FORWARD(gpu, stream_id).put(fwd, A->data(), B->data());

        }

    }

  #endif

}

#endif //OPS_REDUCERS_H


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_sf.h>

BC_INLINE1(GSL_MODE_PREC,gsl_mode_t,unsigned)


text:
#ifndef __GSL_BLOCK_H__
#define __GSL_BLOCK_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 <gsl/gsl_block_complex_long_double.h>
#include <gsl/gsl_block_complex_double.h>
#include <gsl/gsl_block_complex_float.h>

#include <gsl/gsl_block_long_double.h>
#include <gsl/gsl_block_double.h>
#include <gsl/gsl_block_float.h>

#include <gsl/gsl_block_ulong.h>
#include <gsl/gsl_block_long.h>

#include <gsl/gsl_block_uint.h>
#include <gsl/gsl_block_int.h>

#include <gsl/gsl_block_ushort.h>
#include <gsl/gsl_block_short.h>

#include <gsl/gsl_block_uchar.h>
#include <gsl/gsl_block_char.h>

#endif /* __GSL_BLOCK_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#ifndef DATA_H
#define DATA_H

#include <gsl/gsl_vector.h>

typedef struct {
  int imageIdx;
  int labelerId;
  int workflowId;
  int workerPoolId;
  int label;
  double trueLabel;
  double trueDifficulty;
} Label;

typedef struct {
  Label *labels;
  int numLabels;
  int numLabelers;
  int numImages;
  int numWorkflows;
  int numWorkerPools;
  double *priorAlpha;
  double *priorBeta;
  double *probZ1, *probZ0;
  double *alpha, *beta;
  double *priorZ1;
} Dataset;

void packX (gsl_vector *x, Dataset *data);
void unpackX (const gsl_vector *x, Dataset *data);
void readData (char *filename, Dataset *data, double*, double);
void outputResults (Dataset *data);

#endif

max_stars_repo_name: ShreyaR/WorkerPoolSelection

max_issues_repo_name: ShreyaR/WorkerPoolSelection

max_forks_repo_name: ShreyaR/WorkerPoolSelection


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for structures for EOBNRv2HM reduced order model (non-spinning version).
 * See CQG 31 195010, 2014, arXiv:1402.4146 for details on the reduced order method.
 * See arXiv:1106.1021 for the EOBNRv2HM model.
 *
 * Borrows from the SEOBNR ROM LAL code written by Michael Puerrer and John Veitch.
 *
 * Put the untared data in the directory designated by the environment variable ROM_DATA_PATH.
 *
 * Parameter range:
 *   q = 1-12 (almost)
 *   No spin
 *   Mtot >= 10Msun for fstart=8Hz
 *
 */

#ifndef _EOBNRV2HMROMSTRUCT_H
#define _EOBNRV2HMROMSTRUCT_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 ******************/

typedef struct tagSplineList {
    gsl_spline*            spline; /* The gsl spline */
    gsl_interp_accel*      accel; /* The gsl accelerator */
    int                    i; /* Index in the list  */
    struct tagSplineList*  next; /* Pointer to the next list element */
} SplineList;

typedef struct tagEOBNRHMROMdata
{
  gsl_vector* q;
  gsl_vector* freq;
  gsl_matrix* Camp;
  gsl_matrix* Cphi;
  gsl_matrix* Bamp;
  gsl_matrix* Bphi;
  gsl_vector* shifttime;
  gsl_vector* shiftphase;
} EOBNRHMROMdata;

typedef struct tagEOBNRHMROMdata_interp
{
  SplineList* Camp_interp; /* List of splines for the amp coefficients - SplineList, index of reduced basis */
  SplineList* Cphi_interp; /* List of splines for the phase coefficients - SplineList, index of reduced basis */
  SplineList* shifttime_interp; /* interpolated shift in time - SplineList with one element */
  SplineList* shiftphase_interp; /* interpolated shift in phase - SplineList with one element */
} EOBNRHMROMdata_interp;

typedef struct tagEOBNRHMROMdata_coeff
{
  gsl_vector* Camp_coeff;
  gsl_vector* Cphi_coeff;
  double*     shifttime_coeff;
  double*     shiftphase_coeff;
} EOBNRHMROMdata_coeff;

typedef struct tagListmodesEOBNRHMROMdata
{
    EOBNRHMROMdata*                     data; /* The ROM data. */
    int                                 l; /* Node mode l  */
    int                                 m; /* Node submode m  */
    struct tagListmodesEOBNRHMROMdata*  next; /* next pointer */
} ListmodesEOBNRHMROMdata;

typedef struct tagListmodesEOBNRHMROMdata_interp
{
    EOBNRHMROMdata_interp*                     data_interp; /* The splines built from the coefficients. */
    int                                        l; /* Node mode l  */
    int                                        m; /* Node submode m  */
    struct tagListmodesEOBNRHMROMdata_interp*  next; /* next pointer */
} ListmodesEOBNRHMROMdata_interp;

typedef struct tagListmodesEOBNRHMROMdata_coeff
{
    EOBNRHMROMdata_coeff*                     data_coeff; /* The data of coefficients. */
    int                                       l; /* Node mode l  */
    int                                       m; /* Node submode m  */
    struct tagListmodesEOBNRHMROMdata_coeff*  next; /* next pointer */
} ListmodesEOBNRHMROMdata_coeff;

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

/* Functions associated to list manipulations */
SplineList* SplineList_AddElementNoCopy(
	   SplineList* appended,  /* List structure to prepend to */
	   gsl_spline* spline,  /* spline to contain */
           gsl_interp_accel* accel,  /* accelerator to contain */
	   int  i /* index in the list */
);
SplineList* SplineList_GetElement(
	   SplineList* const splinelist,  /* List structure to get a particular mode from */
	   const int i /* index in the list */
);
void SplineList_Destroy(
	   SplineList* list  /* List structure to destroy; notice that the content is destroyed too */
);
ListmodesEOBNRHMROMdata* ListmodesEOBNRHMROMdata_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata* indata,  /* data to contain */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
ListmodesEOBNRHMROMdata* ListmodesEOBNRHMROMdata_GetMode(
	   ListmodesEOBNRHMROMdata* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
void ListmodesEOBNRHMROMdata_Destroy(
	   ListmodesEOBNRHMROMdata* list  /* List structure to destroy; notice that the data is destroyed too */
);
ListmodesEOBNRHMROMdata_interp* ListmodesEOBNRHMROMdata_interp_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata_interp* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata_interp* data,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */
);
ListmodesEOBNRHMROMdata_interp* ListmodesEOBNRHMROMdata_interp_GetMode(
	   ListmodesEOBNRHMROMdata_interp* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
void ListmodesEOBNRHMROMdata_interp_Destroy(
	   ListmodesEOBNRHMROMdata_interp* list  /* List structure to destroy; notice that the data is destroyed too */
);
ListmodesEOBNRHMROMdata_coeff* ListmodesEOBNRHMROMdata_coeff_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata_coeff* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata_coeff* data,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */
);
ListmodesEOBNRHMROMdata_coeff* ListmodesEOBNRHMROMdata_coeff_GetMode(
	   ListmodesEOBNRHMROMdata_coeff* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */
);
void ListmodesEOBNRHMROMdata_coeff_Destroy(
	   ListmodesEOBNRHMROMdata_coeff* list  /* List structure to destroy; notice that the data is destroyed too */
);

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

#endif /* _EOBNRV2HMROMSTRUCT_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_movstat.h>

/* compute moving mean by explicitely constructing window and computing mean */
int
slow_movmean(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y,
             const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double mean = gsl_stats_mean(window, 1, wsize);

      gsl_vector_set(y, i, mean);
    }

  free(window);

  return GSL_SUCCESS;
}

static double
func_mean(const size_t n, double x[], void * params)
{
  (void) params;
  return gsl_stats_mean(x, 1, n);
}

static void
test_mean_proc(const double tol, const size_t n, const size_t H, const size_t J,
               const gsl_movstat_end_t etype, gsl_rng * rng_p)
{
  gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
  gsl_vector * x = gsl_vector_alloc(n);
  gsl_vector * y = gsl_vector_alloc(n);
  gsl_vector * z = gsl_vector_alloc(n);
  gsl_movstat_function F;
  char buf[2048];

  F.function = func_mean;
  F.params = NULL;

  random_vector(x, rng_p);

  /* y = mean(x) with slow brute force method */
  slow_movmean(etype, x, y, H, J);

  /* y = mean(x) with fast method */
  gsl_movstat_mean(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u mean random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = mean(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_mean(etype, z, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u mean random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = mean(x) with user-defined function */
  gsl_movstat_apply(etype, &F, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u mean user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  gsl_movstat_free(w);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
}

static void
test_mean(gsl_rng * rng_p)
{
  const double eps = 1.0e-10;

  test_mean_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 1000, 0, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 1000, 5, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_mean_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_mean_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 1000, 0, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 1000, 5, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_mean_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_mean_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 1000, 0, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 1000, 5, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_mean_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* histogram/add2d.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>

#include "find2d.c"

int
gsl_histogram2d_increment (gsl_histogram2d * h, double x, double y)
{
  int status = gsl_histogram2d_accumulate (h, x, y, 1.0);
  return status;
}

int
gsl_histogram2d_accumulate (gsl_histogram2d * h,
                            double x, double y, double weight)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;

  size_t i = 0, j = 0;

  int status = find2d (h->nx, h->xrange, h->ny, h->yrange, x, y, &i, &j);

  if (status)
    {
      return GSL_EDOM;
    }

  if (i >= nx)
    {
      GSL_ERROR ("index lies outside valid range of 0 .. nx - 1",
                 GSL_ESANITY);
    }

  if (j >= ny)
    {
      GSL_ERROR ("index lies outside valid range of 0 .. ny - 1",
                 GSL_ESANITY);
    }

  h->bin[i * ny + j] += weight;

  return GSL_SUCCESS;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* rng/types.c
 *
 * Copyright (C) 2001, 2007 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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define N 100

const gsl_rng_type * gsl_rng_generator_types[N];

#define ADD(t) {if (i==N) abort(); gsl_rng_generator_types[i] = (t); i++; };

const gsl_rng_type **
gsl_rng_types_setup (void)
{
  int i = 0;

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

  return gsl_rng_generator_types;
}


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_dscal (const int N, const double alpha, double *X, const int incX)
{
#define BASE double
#include "source_scal_r.h"
#undef BASE
}


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "Data/Data.h"
#include <gsl/gsl>

// forward declare cv::KeyPoint
namespace cv
{
    class KeyPoint;
}

namespace mage
{
    struct InitializationData;
    struct FrameData;
    class AnalyzedImage;

    class Introspector
    {
    public:
        virtual ~Introspector() = default;

        virtual void Introspect(const InitializationData& /*data*/) {}

        virtual void IntrospectEstimatedPose(const mage::FrameId& /*frameId*/, const mage::Matrix& /*viewMatrix*/) {}

        virtual void IntrospectAnalyzedImage(const mage::FrameData& /*frame*/, const mage::AnalyzedImage& /* image */) {}
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
/**
 *
 * @file core_zgessm.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <cblas.h>
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_zgessm applies the factors L computed by CORE_zgetrf_incpiv to
 *  a complex M-by-N tile A.
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the tile A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile A.  N >= 0.
 *
 * @param[in] K
 *         The number of columns of the tile L. K >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in] IPIV
 *         The pivot indices array of size K as returned by
 *         CORE_zgetrf_incpiv.
 *
 * @param[in] L
 *         The M-by-K lower triangular tile.
 *
 * @param[in] LDL
 *         The leading dimension of the array L.  LDL >= max(1,M).
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile A.
 *         On exit, updated by the application of L.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgessm = PCORE_zgessm
#define CORE_zgessm PCORE_zgessm
#endif
int CORE_zgessm(int M, int N, int K, int IB,
                const int *IPIV,
                const PLASMA_Complex64_t *L, int LDL,
                PLASMA_Complex64_t *A, int LDA)
{
    static PLASMA_Complex64_t zone  =  1.0;
    static PLASMA_Complex64_t mzone = -1.0;
    static int                ione  =  1;

    int i, sb;
    int tmp, tmp2;

    /* Check input arguments */
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (K < 0) {
        coreblas_error(3, "Illegal value of K");
        return -3;
    }
    if (IB < 0) {
        coreblas_error(4, "Illegal value of IB");
        return -4;
    }
    if ((LDL < max(1,M)) && (M > 0)) {
        coreblas_error(7, "Illegal value of LDL");
        return -7;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(9, "Illegal value of LDA");
        return -9;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (K == 0) || (IB == 0))
        return PLASMA_SUCCESS;

    for(i = 0; i < K; i += IB) {
        sb = min(IB, K-i);
        /*
         * Apply interchanges to columns I*IB+1:IB*( I+1 )+1.
         */
        tmp  = i+1;
        tmp2 = i+sb;
        LAPACKE_zlaswp_work(LAPACK_COL_MAJOR, N, A, LDA, tmp, tmp2, IPIV, ione);
        /*
         * Compute block row of U.
         */
        cblas_ztrsm(
            CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasUnit,
            sb, N, CBLAS_SADDR(zone),
            &L[LDL*i+i], LDL,
            &A[i], LDA );

        if (i+sb < M) {
        /*
        * Update trailing submatrix.
        */
        cblas_zgemm(
            CblasColMajor, CblasNoTrans, CblasNoTrans,
            M-(i+sb), N, sb,
            CBLAS_SADDR(mzone), &L[LDL*i+(i+sb)], LDL,
            &A[i], LDA,
            CBLAS_SADDR(zone), &A[i+sb], LDA );
        }
    }
    return PLASMA_SUCCESS;
}


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef __COMPUTELISASNR_H__
#define __COMPUTELISASNR_H__ 1

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

#define _XOPEN_SOURCE 500

#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"
#include "struct.h"
#include "waveform.h"
#include "fft.h"
#include "LISAgeometry.h"
#include "LISAFDresponse.h"
#include "LISAutils.h"

/********************************** Structures ******************************************/

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagComputeLISASNRparams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double fRef;               /* reference frequency at which phiRef is set (Hz, default 0 which is interpreted as Mf=0.14) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double inclination;        /* inclination of source (rad, default pi/3) */
  double lambda;             /* first angle for the position in the sky (rad, default 0) */
  double beta;               /* second angle for the position in the sky (rad, default 0) */
  double polarization;       /* polarization angle (rad, default 0) */

  int nbmode;                /* number of modes to generate (starting with 22) - defaults to 5 (all modes) */
  double deltatobs;          /* Observation duration (years, default=2) */
  double minf;               /* Minimal frequency (Hz, default=0) - when set to 0, use the lowest frequency where the detector noise model is trusted __LISASimFD_Noise_fLow (set somewhat arbitrarily)*/
  double maxf;               /* Maximal frequency (Hz, default=0) - when set to 0, use the highest frequency where the detector noise model is trusted __LISASimFD_Noise_fHigh (set somewhat arbitrarily)*/
  int tagextpn;              /* Tag to allow PN extension of the waveform at low frequencies */
  double Mfmatch;            /* When PN extension allowed, geometric matching frequency: will use ROM above this value. If <=0, use ROM down to the lowest covered frequency */
  int tagtdi;                /* Tag selecting the desired TDI observables */
  int tagint;                /* Tag choosing the integrator: 0 for Fresnel (default), 1 for linear integration */
  int nbptsoverlap;          /* Number of points to use in loglinear overlaps (default 32768) */
  LISAconstellation *variant;  /* A structure defining the LISA constellation features */
  int frozenLISA;            /* Freeze the orbital configuration to the time of peak of the injection (default 0) */
  ResponseApproxtag responseapprox;    /* Approximation in the GAB and orb response - choices are full (full response, default), lowfL (keep orbital delay frequency-dependence but simplify constellation response) and lowf (simplify constellation and orbital response) - WARNING : at the moment noises are not consistent, and TDI combinations from the GAB are unchanged */
  int fromtditdfile;         /* Option for loading time series for TDI observables and FFTing */
  int nlinesinfile;          /* Number of lines of input file */
  char indir[256];           /* Input directory */
  char infile[256];          /* Input file */
  int loadparamsfile;        /* Option to load physical parameters from file and to output result to file (default 0) */
  int nlinesparams;          /* Number of lines in params file */
  char paramsdir[256];       /* Directory for the input/output file */
  char paramsfile[256];      /* Input file with the parameters */
  char outputfile[256];      /* Output file */
} ComputeLISASNRparams;

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

#endif /* _COMPUTELISASNR_H */

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
#include <vector>
#include <string>
#include <math.h>
#include <cmath>
#include <memory>
#include <iostream>

#include "hdf5.h"
#include "hdf5_hl.h"
#include "H5Gpublic.h"
#include "H5Fpublic.h"

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_integration.h>

#ifndef NUFATE_H
#define NUFATE_H

namespace nufate{

/// \struct Square_matrix_double
/// \brief Very simple matrix struct used for pybinding.
/// vec_ is the front address of a chank of double data,
/// which will be shaped in NxN 2dim matrix in python.
struct Square_matrix_double {
 public :
 unsigned int n_;
 std::shared_ptr<double> vec_;
};

///\class Result
///\brief Simple class to hold the results.
class Result {
  public :
   std::vector<double> eval;          // eigen values
   std::shared_ptr<double> evec;      // pointer to top address of eigen vectors
   std::vector<double> ci;            // coefficients
   std::vector<double> energy_nodes_; // energy nodes in GeV
   std::vector<double> phi_0_;        // initial flux * E^(pedestal_gamma)

   //
   // functions below are not used by main nuFATE program,
   // but may be used for pybinding.
   //

   /// \brief getter for eigenvalues
   /// @return 1D vector
   std::vector<double> get_eigenvalues() { return eval; }

   /// \brief getter for coefficients
   /// @return 1D vector
   std::vector<double> get_coefficients()   { return ci; }

   /// \brief getter for energy nodes
   /// @return 1D vector
   std::vector<double> get_energy_nodes() { return energy_nodes_; }

   /// \brief getter for initial flux * E^(pedestal_gamma)
   /// @return 1D vector
   std::vector<double> get_phi_0() { return phi_0_; }

   /// \brief converter from pointer of double to Square_matrix_double object
   /// @return Square_matrix_double
   Square_matrix_double get_eigenvec_matrix() {
      unsigned int n = energy_nodes_.size();
      Square_matrix_double smatrix;
      // evec is n x n data
      smatrix.n_ = n;
      smatrix.vec_ = evec;
      return smatrix;
   }

   /// \brief debug print function
   void Print(unsigned int target_index = 0) {
      unsigned int dim = energy_nodes_.size();

      std::cout << "***** Result print *****"  << std::endl;
      std::cout << "eigenvectors at row " << target_index << " ---"  << std::endl;
      for (unsigned int i=target_index; i<target_index+1; ++i) {
         for (unsigned int j=0; j<dim; ++j) {
             std::cout << *(evec.get() + i*dim + j) << " " ;
         }
         std::cout << std::endl;
      }

      std::cout << std::endl;
      std::cout << "eigenvalue at i = " << target_index << " ---"  << std::endl;
      std::cout << eval[target_index] << std::endl;

      std::cout << std::endl;
      std::cout << "coefficient at i = " << target_index << " ---"  << std::endl;
      std::cout << ci[target_index] << std::endl;

      std::cout << std::endl;
      std::cout << "energy_node at i = " << target_index << " ---"  << std::endl;
      std::cout << energy_nodes_[target_index] << std::endl;
   }
};

///\class nuFATE
///\brief nuFATE main class
class nuFATE {
  private:
    const double GF = 1.16e-5;
    const double hbarc = 1.97e-14;
    const double GW = 2.085;
    const double MW = 80.385e0;
    const double mmu = 0.106e0;
    const double me = 511.e-6;
    const double pi = M_PI;
    const double MZ = 91.18;
    const double s2t = 0.23;
    const double gL =  s2t-0.5;
    const double REarth = 6371.;
    const double gR = s2t;
  private:
    int newflavor_;
    double newgamma_;
    std::string newh5_filename_;
  public:
    /// \brief Constructor
    /// @param flavor position ID of the system (1:NuE, -1:NuEBar, 2:NuMu, -2:NuMuBar, 3:NuTau, -3:NuTauBar)
    /// @param gamma spectral index of input flux.
    /// @param h5_filename name of hdf5 file containing the cross sections.
    /// @param include_secondaries if true secondaries are added to the flux propagation.
    nuFATE(int flv, double gamma, std::string h5_filename, bool include_secondaries);
    /// \brief Constructor
    /// @param flavor position ID of the system (1:NuE, -1:NuEBar, 2:NuMu, -2:NuMuBar, 3:NuTau, -3:NuTauBar)
    /// @param gamma spectral index of input flux.
    /// @param energy_nodes energy nodes in GeV.
    /// @param sigma_array total cross section(CC+NC) at each energy node in cm^2.
    /// @param dsigma_dE square array of differential cross section (NC only) in cm^2/GeV.
    /// @param include_secondaries if true secondaries are added to the flux propagation.
    nuFATE(int flv, double gamma, std::vector<double> energy_nodes, std::vector<double> sigma_array, std::vector<std::vector<double>> dsigma_dE, bool include_secondaries);
    /// \brief Constructor
    /// @param flavor position ID of the system (1:NuE, -1:NuEBar, 2:NuMu, -2:NuMuBar, 3:NuTau, -3:NuTauBar)
    /// @param gamma spectral index of input flux.
    /// @param cc_sigma_file file path for charge current(CC) total cross section in input format of nuSQuIDS
    /// @param nc_sigma_file file path for neutral current(NC) total cross section in input format of nuSQuIDS
    /// @param nc_dsigma_file file path to NC dSigma/dE differential cross section in input format of nuSQuIDS
    /// @param include_secondaries if true secondaries are added to the flux propagation.
    nuFATE(int flv, double gamma, std::string cc_sigma_file, std::string nc_sigma_file, std::string nc_dsigma_file, bool include_secondaries);

    /// \brief Eigensystem calculator
    Result getEigensystem();
    /// \brief Function to get Earth column density
    /// @param theta Zenith angle in radians.
    double getEarthColumnDensity(double theta);
    /// \brief Function to get flavor
    int getFlavor() const;
    /// \brief Function to get input spectral index
    double getGamma() const;
    /// \brief Function to get input filename
    std::string getFilename() const;
    /// \brief Function to get number of energy nodes
    double getNumNodes() const;
    /// \brief Function to toggle secondaries
    void setAddSecondaries(bool opt) { add_secondary_term_ = opt;}

    /// \brief Function to get energy nodes
    std::vector<double> getEnergyNodes() { return energy_nodes_; }
    /// \brief Function to get total crosssections (CC + NC)
    std::vector<double> getTotalCrossSections() { return sigma_array_; }
    /// \brief Function to get dsigma_dE differential crossection (for neutral-current only)
    Square_matrix_double getNCDifferentialCrossSections() {
       Square_matrix_double smatrix;
       smatrix.n_ = NumNodes_;
       smatrix.vec_ = dxs_array_;
       return smatrix;
    }

    /// \brief Calculate Relative Attenuation
    /// @param total number of targets, X[g/cm^2]*Na(Avogadro number) for example
    /// @return attenuation factor (arrival_flux / initial_flux), 1D array in energy bins
    std::vector<double> getRelativeAttenuation(double number_of_targets);

  protected:
    void Init(double gamma, const std::vector<double> &enodes, const std::vector<double> &sigmas, const std::vector<std::vector<double> > &dsigma);
    void AddAdditionalTerms();
    void LoadCrossSectionFromHDF5();
    void SetCrossSectionsFromInput(const std::vector<double> &sigma, const std::vector<std::vector<double>> &dsigma_dE);
    void SetInitialFlux();
    void set_glashow_total();
    void set_glashow_partial();
    void set_RHS_matrices(std::shared_ptr<double> RMatrix_, std::shared_ptr<double> dxs_array);
    static double rho_earth(double theta, void * p); //returns density for a given angle in radians along x distance
    double readDoubleAttribute(hid_t, std::string) const;
    unsigned int readUIntAttribute(hid_t, std::string) const;
    std::vector<double> logspace(double min,double max,unsigned int samples) const;
  private:
    void AllocateMemoryForMembers(unsigned int num_nodes);
    void SetEnergyBinWidths();
    std::string grpdiff_;
    std::string grptot_;
    hid_t file_id_;
    hid_t root_id_;
    unsigned int NumNodes_;
    unsigned int rsize_;
    double Emax_;
    double Emin_;
    int dxsdim_[2];
    std::vector<double> energy_nodes_;
    std::vector<double> sigma_array_;
    std::vector<double> sigma_array_orig_;
    std::vector<double> DeltaE_;
    std::vector<double> phi_0_;
    std::vector<double> glashow_total_;
    std::vector<double> sig3_array_;
    std::shared_ptr<double> glashow_partial_;
    std::shared_ptr<double> dxs_array_;
    std::shared_ptr<double> tau_array_;
    std::shared_ptr<double> sec_array_;
    std::shared_ptr<double> regen_array_;
    std::shared_ptr<double> RHSMatrix_;
    std::shared_ptr<double> RHSMatrix1_;
    std::shared_ptr<double> RHSMatrix2_;
    std::shared_ptr<double> RHSMatrix3_;
    std::shared_ptr<double> RHSMatrix4_;
    std::shared_ptr<double> RHregen_;
    std::shared_ptr<double> Enuin_;
    std::shared_ptr<double> Enu_;
    std::shared_ptr<double> selectron_;
    std::shared_ptr<double> den_;
    std::shared_ptr<double> t1_;
    std::shared_ptr<double> t2_;
    std::shared_ptr<double> t3_;
  private:
    bool include_secondaries_ = false;
    bool memory_allocated_ = false;
    bool initial_flux_set_ = false;
    bool total_cross_section_set_ = false;
    bool differential_cross_section_set_ = false;
    bool RHS_set_ = false;
    bool add_tau_regeneration_ = true;
    bool add_glashow_term_= true;
    bool add_secondary_term_ = true;
};

}

#endif

max_stars_repo_name: kotoyo/nuFATE

max_issues_repo_name: kotoyo/nuFATE

max_forks_repo_name: kotoyo/nuFATE


text:
#pragma once
#include "halley/support/logger.h"
#include "halley/net/connection/network_message.h"
#include <gsl/gsl>

namespace Halley
{
	class Serializer;
	class String;
	class MessageQueue;

	namespace DevCon
	{
		void setupMessageQueue(MessageQueue& queue);

		enum class MessageType
		{
			Log,
			ReloadAssets
		};

		class DevConMessage : public NetworkMessage
		{
		public:
			virtual ~DevConMessage() = default;
			virtual MessageType getMessageType() const = 0;
		};

		class LogMsg : public DevConMessage
		{
		public:
			LogMsg(gsl::span<const gsl::byte> data);
			LogMsg(LoggerLevel level, const String& msg);

			void serialize(Serializer& s) const override;

			LoggerLevel getLevel() const;
			const String& getMessage() const;
			MessageType getMessageType() const override;

		private:
			LoggerLevel level;
			String msg;
		};

		class ReloadAssetsMsg : public DevConMessage
		{
		public:
			ReloadAssetsMsg(gsl::span<const gsl::byte> data);
			ReloadAssetsMsg(std::vector<String> ids);

			void serialize(Serializer& s) const override;

			std::vector<String> getIds() const;

			MessageType getMessageType() const override;

		private:
			std::vector<String> ids;
		};
	}
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley