text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 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, 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.

 $Id: parse_exp.c 13079 2015-02-18 16:06:42Z dstrubbe $
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

#include "liboct_parser.h"
#include "symbols.h"

static char *par_string;
static int par_pos;
parse_result par_res;

static int oct_parser_lex(void);
static int oct_parser_error (char *s)  /* Called by oct_parser_parse on error */
{
  /* Do nothing */
  /* printf("%s\n", s); */
  return 0;
}

/* include the parser */
#include "grammar.c"

int parse_exp(char *exp, parse_result *r)
{
  int o;

  par_string = exp;
  par_pos = 0;

  o = oct_parser_parse();
  if(o == 0){
    r->type = par_res.type;
    if(r->type == PR_CMPLX)
      r->value.c = par_res.value.c;
    else
      r->value.s = par_res.value.s;
  }
  return o;
}

int get_real(char *s, double *d)
{
  int n=0;
  sscanf(s, "%lg", d);
  while(*s && (isdigit(*s) || *s=='.' || *s=='e' || *s=='E')){
    if((*s=='e' || *s=='E') && (*(s+1)=='+' || *(s+1)=='-')) {s++; n++;}
    s++; n++;
  }
  return n;
}

static int oct_parser_lex (){
  int c;
  char *symbuf = 0, *symbuf2 = 0;
  int length = 0;

  /* Ignore whitespace, get first nonwhite character.  */
  while ((c = par_string[par_pos++]) == ' ' || c == '\t');

  if (c == '\0')
    return '\n';

  /* Char starts a number => parse the number.         */
  if (c == '.' || isdigit (c)){
    par_pos--;
    par_pos += get_real(&par_string[par_pos], &GSL_REAL(yylval.val));
    return NUM;
  }

  /* get the logical operators */
  if (c == '<' && par_string[par_pos] == '='){
    par_pos++;
    return LE;
  }

  if (c == '>' && par_string[par_pos] == '='){
    par_pos++;
    return GE;
  }

  if (c == '=' && par_string[par_pos] == '='){
    par_pos++;
    return EQUAL;
  }

  /*
  if (c == ':' && par_string[par_pos] == '='){
    par_pos++;
    return SET;
  }
  */

  if (c == '&' && par_string[par_pos] == '&'){
    par_pos++;
    return LAND;
  }

  if (c == '|' && par_string[par_pos] == '|'){
    par_pos++;
    return LOR;
  }

  /* Char starts an identifier => read the name.       */
  if (isalpha (c) || c == '\'' || c == '\"'){
    symrec *s;
    char startc = (char)c;
    int i;

    /* Initially make the buffer long enough
       for a 40-character symbol name.  */
    if (length == 0)
      length = 40, symbuf = (char *)malloc (length + 1);

    if(startc == '\'' || startc == '\"')
      c = par_string[par_pos++];
    else
      startc = 0; /* false */

    i = 0;
    do{
      /* If buffer is full, make it larger.        */
      if (i == length){
	length *= 2;
	symbuf2 = (char *)realloc (symbuf, length + 1);
	if (symbuf2 == NULL) {
	  fprintf(stderr, "Parser error: failed to reallocate buffer to size %i.\n", length);
	  exit(1);
	} else {
	  symbuf = symbuf2;
	}
      }
      /* Add this character to the buffer.         */
      symbuf[i++] = (char)c;
      /* Get another character.                    */
      c = par_string[par_pos++];
    }while (c != '\0' &&
	    ((startc && c!=startc) || (!startc && (isalnum(c) || c == '_' ))));

    if(!startc) par_pos--;
    symbuf[i] = '\0';

    if(!startc){
      s = getsym (symbuf);
      if (s == 0){
	int jj;
	for (jj = 0; reserved_symbols[jj] != 0; jj++){
	  if(strcmp(symbuf, reserved_symbols[jj]) == 0){
	    fprintf(stderr, "Parser error: trying to redefine reserved symbol '%s'.\n", symbuf);
	    exit(1);
	  }
	}

	s = putsym (symbuf, S_CMPLX);
      }
      yylval.tptr = s;

      free(symbuf);
      if(s->type == S_CMPLX)
	return VAR;
      else
	return FNCT;
    }else{
      yylval.str = strdup(symbuf);

      free(symbuf);
      return STR;
    }
  }

  /* Any other character is a token by itself.        */
  return c;
}

max_stars_repo_name: gimunu/octopus-metric

max_issues_repo_name: gimunu/octopus-metric

max_forks_repo_name: gimunu/octopus-metric


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Mon Mar  8 17:44:42 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -notwiddle 15 */

/*
 * This function contains 156 FP additions, 56 FP multiplications,
 * (or, 128 additions, 28 multiplications, 28 fused multiply/add),
 * 62 stack variables, and 60 memory accesses
 */
static const fftw_real K587785252 = FFTW_KONST(+0.587785252292473129168705954639072768597652438);
static const fftw_real K951056516 = FFTW_KONST(+0.951056516295153572116439333379382143405698634);
static const fftw_real K250000000 = FFTW_KONST(+0.250000000000000000000000000000000000000000000);
static const fftw_real K559016994 = FFTW_KONST(+0.559016994374947424102293417182819058860154590);
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.36 1999/02/19 17:22:11 athena Exp $
 * $Id: fft.ml,v 1.41 1999/02/19 17:22:13 athena Exp $
 * $Id: to_c.ml,v 1.24 1999/02/19 17:22:17 athena Exp $
 */

void fftw_no_twiddle_15(const fftw_complex *input, fftw_complex *output, int istride, int ostride)
{
     fftw_real tmp5;
     fftw_real tmp33;
     fftw_real tmp57;
     fftw_real tmp145;
     fftw_real tmp100;
     fftw_real tmp124;
     fftw_real tmp21;
     fftw_real tmp26;
     fftw_real tmp27;
     fftw_real tmp49;
     fftw_real tmp54;
     fftw_real tmp55;
     fftw_real tmp136;
     fftw_real tmp137;
     fftw_real tmp147;
     fftw_real tmp61;
     fftw_real tmp62;
     fftw_real tmp63;
     fftw_real tmp112;
     fftw_real tmp113;
     fftw_real tmp126;
     fftw_real tmp83;
     fftw_real tmp88;
     fftw_real tmp94;
     fftw_real tmp10;
     fftw_real tmp15;
     fftw_real tmp16;
     fftw_real tmp38;
     fftw_real tmp43;
     fftw_real tmp44;
     fftw_real tmp139;
     fftw_real tmp140;
     fftw_real tmp146;
     fftw_real tmp58;
     fftw_real tmp59;
     fftw_real tmp60;
     fftw_real tmp115;
     fftw_real tmp116;
     fftw_real tmp125;
     fftw_real tmp72;
     fftw_real tmp77;
     fftw_real tmp93;
     ASSERT_ALIGNED_DOUBLE();
     {
	  fftw_real tmp1;
	  fftw_real tmp97;
	  fftw_real tmp4;
	  fftw_real tmp96;
	  fftw_real tmp32;
	  fftw_real tmp98;
	  fftw_real tmp29;
	  fftw_real tmp99;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp1 = c_re(input[0]);
	  tmp97 = c_im(input[0]);
	  {
	       fftw_real tmp2;
	       fftw_real tmp3;
	       fftw_real tmp30;
	       fftw_real tmp31;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp2 = c_re(input[5 * istride]);
	       tmp3 = c_re(input[10 * istride]);
	       tmp4 = tmp2 + tmp3;
	       tmp96 = K866025403 * (tmp3 - tmp2);
	       tmp30 = c_im(input[5 * istride]);
	       tmp31 = c_im(input[10 * istride]);
	       tmp32 = K866025403 * (tmp30 - tmp31);
	       tmp98 = tmp30 + tmp31;
	  }
	  tmp5 = tmp1 + tmp4;
	  tmp29 = tmp1 - (K500000000 * tmp4);
	  tmp33 = tmp29 - tmp32;
	  tmp57 = tmp29 + tmp32;
	  tmp145 = tmp97 + tmp98;
	  tmp99 = tmp97 - (K500000000 * tmp98);
	  tmp100 = tmp96 + tmp99;
	  tmp124 = tmp99 - tmp96;
     }
     {
	  fftw_real tmp17;
	  fftw_real tmp20;
	  fftw_real tmp45;
	  fftw_real tmp79;
	  fftw_real tmp80;
	  fftw_real tmp81;
	  fftw_real tmp48;
	  fftw_real tmp82;
	  fftw_real tmp22;
	  fftw_real tmp25;
	  fftw_real tmp50;
	  fftw_real tmp84;
	  fftw_real tmp85;
	  fftw_real tmp86;
	  fftw_real tmp53;
	  fftw_real tmp87;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp18;
	       fftw_real tmp19;
	       fftw_real tmp46;
	       fftw_real tmp47;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp17 = c_re(input[6 * istride]);
	       tmp18 = c_re(input[11 * istride]);
	       tmp19 = c_re(input[istride]);
	       tmp20 = tmp18 + tmp19;
	       tmp45 = tmp17 - (K500000000 * tmp20);
	       tmp79 = K866025403 * (tmp19 - tmp18);
	       tmp80 = c_im(input[6 * istride]);
	       tmp46 = c_im(input[11 * istride]);
	       tmp47 = c_im(input[istride]);
	       tmp81 = tmp46 + tmp47;
	       tmp48 = K866025403 * (tmp46 - tmp47);
	       tmp82 = tmp80 - (K500000000 * tmp81);
	  }
	  {
	       fftw_real tmp23;
	       fftw_real tmp24;
	       fftw_real tmp51;
	       fftw_real tmp52;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp22 = c_re(input[9 * istride]);
	       tmp23 = c_re(input[14 * istride]);
	       tmp24 = c_re(input[4 * istride]);
	       tmp25 = tmp23 + tmp24;
	       tmp50 = tmp22 - (K500000000 * tmp25);
	       tmp84 = K866025403 * (tmp24 - tmp23);
	       tmp85 = c_im(input[9 * istride]);
	       tmp51 = c_im(input[14 * istride]);
	       tmp52 = c_im(input[4 * istride]);
	       tmp86 = tmp51 + tmp52;
	       tmp53 = K866025403 * (tmp51 - tmp52);
	       tmp87 = tmp85 - (K500000000 * tmp86);
	  }
	  tmp21 = tmp17 + tmp20;
	  tmp26 = tmp22 + tmp25;
	  tmp27 = tmp21 + tmp26;
	  tmp49 = tmp45 - tmp48;
	  tmp54 = tmp50 - tmp53;
	  tmp55 = tmp49 + tmp54;
	  tmp136 = tmp80 + tmp81;
	  tmp137 = tmp85 + tmp86;
	  tmp147 = tmp136 + tmp137;
	  tmp61 = tmp45 + tmp48;
	  tmp62 = tmp50 + tmp53;
	  tmp63 = tmp61 + tmp62;
	  tmp112 = tmp82 - tmp79;
	  tmp113 = tmp87 - tmp84;
	  tmp126 = tmp112 + tmp113;
	  tmp83 = tmp79 + tmp82;
	  tmp88 = tmp84 + tmp87;
	  tmp94 = tmp83 + tmp88;
     }
     {
	  fftw_real tmp6;
	  fftw_real tmp9;
	  fftw_real tmp34;
	  fftw_real tmp68;
	  fftw_real tmp69;
	  fftw_real tmp70;
	  fftw_real tmp37;
	  fftw_real tmp71;
	  fftw_real tmp11;
	  fftw_real tmp14;
	  fftw_real tmp39;
	  fftw_real tmp73;
	  fftw_real tmp74;
	  fftw_real tmp75;
	  fftw_real tmp42;
	  fftw_real tmp76;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp7;
	       fftw_real tmp8;
	       fftw_real tmp35;
	       fftw_real tmp36;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp6 = c_re(input[3 * istride]);
	       tmp7 = c_re(input[8 * istride]);
	       tmp8 = c_re(input[13 * istride]);
	       tmp9 = tmp7 + tmp8;
	       tmp34 = tmp6 - (K500000000 * tmp9);
	       tmp68 = K866025403 * (tmp8 - tmp7);
	       tmp69 = c_im(input[3 * istride]);
	       tmp35 = c_im(input[8 * istride]);
	       tmp36 = c_im(input[13 * istride]);
	       tmp70 = tmp35 + tmp36;
	       tmp37 = K866025403 * (tmp35 - tmp36);
	       tmp71 = tmp69 - (K500000000 * tmp70);
	  }
	  {
	       fftw_real tmp12;
	       fftw_real tmp13;
	       fftw_real tmp40;
	       fftw_real tmp41;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp11 = c_re(input[12 * istride]);
	       tmp12 = c_re(input[2 * istride]);
	       tmp13 = c_re(input[7 * istride]);
	       tmp14 = tmp12 + tmp13;
	       tmp39 = tmp11 - (K500000000 * tmp14);
	       tmp73 = K866025403 * (tmp13 - tmp12);
	       tmp74 = c_im(input[12 * istride]);
	       tmp40 = c_im(input[2 * istride]);
	       tmp41 = c_im(input[7 * istride]);
	       tmp75 = tmp40 + tmp41;
	       tmp42 = K866025403 * (tmp40 - tmp41);
	       tmp76 = tmp74 - (K500000000 * tmp75);
	  }
	  tmp10 = tmp6 + tmp9;
	  tmp15 = tmp11 + tmp14;
	  tmp16 = tmp10 + tmp15;
	  tmp38 = tmp34 - tmp37;
	  tmp43 = tmp39 - tmp42;
	  tmp44 = tmp38 + tmp43;
	  tmp139 = tmp69 + tmp70;
	  tmp140 = tmp74 + tmp75;
	  tmp146 = tmp139 + tmp140;
	  tmp58 = tmp34 + tmp37;
	  tmp59 = tmp39 + tmp42;
	  tmp60 = tmp58 + tmp59;
	  tmp115 = tmp71 - tmp68;
	  tmp116 = tmp76 - tmp73;
	  tmp125 = tmp115 + tmp116;
	  tmp72 = tmp68 + tmp71;
	  tmp77 = tmp73 + tmp76;
	  tmp93 = tmp72 + tmp77;
     }
     {
	  fftw_real tmp134;
	  fftw_real tmp28;
	  fftw_real tmp133;
	  fftw_real tmp142;
	  fftw_real tmp144;
	  fftw_real tmp138;
	  fftw_real tmp141;
	  fftw_real tmp143;
	  fftw_real tmp135;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp134 = K559016994 * (tmp16 - tmp27);
	  tmp28 = tmp16 + tmp27;
	  tmp133 = tmp5 - (K250000000 * tmp28);
	  tmp138 = tmp136 - tmp137;
	  tmp141 = tmp139 - tmp140;
	  tmp142 = (K951056516 * tmp138) - (K587785252 * tmp141);
	  tmp144 = (K951056516 * tmp141) + (K587785252 * tmp138);
	  c_re(output[0]) = tmp5 + tmp28;
	  tmp143 = tmp134 + tmp133;
	  c_re(output[9 * ostride]) = tmp143 - tmp144;
	  c_re(output[6 * ostride]) = tmp143 + tmp144;
	  tmp135 = tmp133 - tmp134;
	  c_re(output[12 * ostride]) = tmp135 - tmp142;
	  c_re(output[3 * ostride]) = tmp135 + tmp142;
     }
     {
	  fftw_real tmp110;
	  fftw_real tmp56;
	  fftw_real tmp109;
	  fftw_real tmp118;
	  fftw_real tmp120;
	  fftw_real tmp114;
	  fftw_real tmp117;
	  fftw_real tmp119;
	  fftw_real tmp111;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp110 = K559016994 * (tmp44 - tmp55);
	  tmp56 = tmp44 + tmp55;
	  tmp109 = tmp33 - (K250000000 * tmp56);
	  tmp114 = tmp112 - tmp113;
	  tmp117 = tmp115 - tmp116;
	  tmp118 = (K951056516 * tmp114) - (K587785252 * tmp117);
	  tmp120 = (K951056516 * tmp117) + (K587785252 * tmp114);
	  c_re(output[5 * ostride]) = tmp33 + tmp56;
	  tmp119 = tmp110 + tmp109;
	  c_re(output[14 * ostride]) = tmp119 - tmp120;
	  c_re(output[11 * ostride]) = tmp119 + tmp120;
	  tmp111 = tmp109 - tmp110;
	  c_re(output[2 * ostride]) = tmp111 - tmp118;
	  c_re(output[8 * ostride]) = tmp111 + tmp118;
     }
     {
	  fftw_real tmp150;
	  fftw_real tmp148;
	  fftw_real tmp149;
	  fftw_real tmp154;
	  fftw_real tmp156;
	  fftw_real tmp152;
	  fftw_real tmp153;
	  fftw_real tmp155;
	  fftw_real tmp151;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp150 = K559016994 * (tmp146 - tmp147);
	  tmp148 = tmp146 + tmp147;
	  tmp149 = tmp145 - (K250000000 * tmp148);
	  tmp152 = tmp21 - tmp26;
	  tmp153 = tmp10 - tmp15;
	  tmp154 = (K951056516 * tmp152) - (K587785252 * tmp153);
	  tmp156 = (K951056516 * tmp153) + (K587785252 * tmp152);
	  c_im(output[0]) = tmp145 + tmp148;
	  tmp155 = tmp150 + tmp149;
	  c_im(output[6 * ostride]) = tmp155 - tmp156;
	  c_im(output[9 * ostride]) = tmp156 + tmp155;
	  tmp151 = tmp149 - tmp150;
	  c_im(output[3 * ostride]) = tmp151 - tmp154;
	  c_im(output[12 * ostride]) = tmp154 + tmp151;
     }
     {
	  fftw_real tmp129;
	  fftw_real tmp127;
	  fftw_real tmp128;
	  fftw_real tmp123;
	  fftw_real tmp132;
	  fftw_real tmp121;
	  fftw_real tmp122;
	  fftw_real tmp131;
	  fftw_real tmp130;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp129 = K559016994 * (tmp125 - tmp126);
	  tmp127 = tmp125 + tmp126;
	  tmp128 = tmp124 - (K250000000 * tmp127);
	  tmp121 = tmp49 - tmp54;
	  tmp122 = tmp38 - tmp43;
	  tmp123 = (K951056516 * tmp121) - (K587785252 * tmp122);
	  tmp132 = (K951056516 * tmp122) + (K587785252 * tmp121);
	  c_im(output[5 * ostride]) = tmp124 + tmp127;
	  tmp131 = tmp129 + tmp128;
	  c_im(output[11 * ostride]) = tmp131 - tmp132;
	  c_im(output[14 * ostride]) = tmp132 + tmp131;
	  tmp130 = tmp128 - tmp129;
	  c_im(output[2 * ostride]) = tmp123 + tmp130;
	  c_im(output[8 * ostride]) = tmp130 - tmp123;
     }
     {
	  fftw_real tmp95;
	  fftw_real tmp101;
	  fftw_real tmp102;
	  fftw_real tmp106;
	  fftw_real tmp107;
	  fftw_real tmp104;
	  fftw_real tmp105;
	  fftw_real tmp108;
	  fftw_real tmp103;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp95 = K559016994 * (tmp93 - tmp94);
	  tmp101 = tmp93 + tmp94;
	  tmp102 = tmp100 - (K250000000 * tmp101);
	  tmp104 = tmp58 - tmp59;
	  tmp105 = tmp61 - tmp62;
	  tmp106 = (K951056516 * tmp104) + (K587785252 * tmp105);
	  tmp107 = (K951056516 * tmp105) - (K587785252 * tmp104);
	  c_im(output[10 * ostride]) = tmp100 + tmp101;
	  tmp108 = tmp102 - tmp95;
	  c_im(output[7 * ostride]) = tmp107 + tmp108;
	  c_im(output[13 * ostride]) = tmp108 - tmp107;
	  tmp103 = tmp95 + tmp102;
	  c_im(output[ostride]) = tmp103 - tmp106;
	  c_im(output[4 * ostride]) = tmp106 + tmp103;
     }
     {
	  fftw_real tmp65;
	  fftw_real tmp64;
	  fftw_real tmp66;
	  fftw_real tmp90;
	  fftw_real tmp92;
	  fftw_real tmp78;
	  fftw_real tmp89;
	  fftw_real tmp91;
	  fftw_real tmp67;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp65 = K559016994 * (tmp60 - tmp63);
	  tmp64 = tmp60 + tmp63;
	  tmp66 = tmp57 - (K250000000 * tmp64);
	  tmp78 = tmp72 - tmp77;
	  tmp89 = tmp83 - tmp88;
	  tmp90 = (K951056516 * tmp78) + (K587785252 * tmp89);
	  tmp92 = (K951056516 * tmp89) - (K587785252 * tmp78);
	  c_re(output[10 * ostride]) = tmp57 + tmp64;
	  tmp91 = tmp66 - tmp65;
	  c_re(output[7 * ostride]) = tmp91 - tmp92;
	  c_re(output[13 * ostride]) = tmp91 + tmp92;
	  tmp67 = tmp65 + tmp66;
	  c_re(output[4 * ostride]) = tmp67 - tmp90;
	  c_re(output[ostride]) = tmp67 + tmp90;
     }
}

fftw_codelet_desc fftw_no_twiddle_15_desc =
{
     "fftw_no_twiddle_15",
     (void (*)()) fftw_no_twiddle_15,
     15,
     FFTW_FORWARD,
     FFTW_NOTW,
     241,
     0,
     (const int *) 0,
};

ext: c

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void free_MCMC_arrays(MCMC_info *MCMC) {
    int           i_DS;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;

    if(MCMC->n_M != NULL) {
        SID_log("Freeing MCMC target arrays for %d dataset(s)...", SID_LOG_OPEN, MCMC->n_DS);

        // Free chain arrays
        i_DS       = 0;
        current_DS = MCMC->DS;
        while(current_DS != NULL) {
            next_DS = current_DS->next;
            SID_free(SID_FARG MCMC->M_new[i_DS]);
            SID_free(SID_FARG MCMC->M_last[i_DS]);
            current_DS = next_DS;
            i_DS++;
        }
        SID_free(SID_FARG MCMC->n_M);
        SID_free(SID_FARG MCMC->M_new);
        SID_free(SID_FARG MCMC->M_last);
        MCMC->n_M_total = 0;

        // Free results arrays
        SID_free(SID_FARG MCMC->P_min);
        SID_free(SID_FARG MCMC->P_max);
        SID_free(SID_FARG MCMC->P_avg);
        SID_free(SID_FARG MCMC->dP_avg);
        SID_free(SID_FARG MCMC->P_best);
        SID_free(SID_FARG MCMC->P_peak);
        SID_free(SID_FARG MCMC->P_lo_68);
        SID_free(SID_FARG MCMC->P_hi_68);
        SID_free(SID_FARG MCMC->P_lo_95);
        SID_free(SID_FARG MCMC->P_hi_95);

        SID_free(SID_FARG MCMC->ln_likelihood_DS);
        SID_free(SID_FARG MCMC->ln_likelihood_DS_best);
        SID_free(SID_FARG MCMC->ln_likelihood_DS_peak);
        SID_free(SID_FARG MCMC->n_DoF_DS);
        SID_free(SID_FARG MCMC->n_DoF_DS_best);
        SID_free(SID_FARG MCMC->n_DoF_DS_peak);

        // These only need to be deallocated if we are using MCMC_MODE_MINIMIZE_IO
        if(SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_MINIMIZE_IO)) {
            SID_free(SID_FARG MCMC->flag_success_buffer);
            SID_free(SID_FARG MCMC->ln_likelihood_new_buffer);
            SID_free(SID_FARG MCMC->P_new_buffer);
            SID_free(SID_FARG MCMC->M_new_buffer);
        }

        SID_log("Done.", SID_LOG_CLOSE);
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


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

#ifndef SiconosBlas_H
#define SiconosBlas_H
#include "SiconosConfig.h"

#if defined(__cplusplus)
extern "C"
{
#endif

// tells include-what-you-use to keep this file
// and not to suggest cblas.h or alike.
// IWYU pragma: begin_exports
#if defined(HAS_MKL_CBLAS)
#include <mkl_cblas.h>
#elif defined(HAS_MATLAB_BLAS)
#include <blas.h>
#define cblas_daxpy daxpy
#define cblas_dcopy dcopy
#define cblas_ddot ddot
#define cblas_dgemm dgemm
#define cblas_dgemv dgemv
#define cblas_dnrm2 dnrm2
#define cblas_dscal dscal
#else
#include <cblas.h>
#endif
// IWYU pragma: end_exports

#ifdef __cplusplus
}
#undef restrict
#define restrict __restrict
#endif

static inline double* NMD_row_rmajor(double* restrict mat, unsigned ncols, unsigned rindx)
{
  return &mat[rindx*ncols];
}

static inline void NMD_copycol_rmajor(int nrows, double* col, double* restrict mat, int ncols, unsigned cindx)
{
  cblas_dcopy(nrows, col, 1, &mat[cindx], ncols);
}

static inline void NMD_dense_gemv(int nrows, int ncols, double alpha, double* restrict mat, double* restrict y, double beta, double* restrict x)
{
  cblas_dgemv(CblasColMajor, CblasTrans, ncols, nrows, alpha, mat, ncols, y, 1, beta, x, 1);
}

#endif // SiconosBlas_H

ext: h

max_stars_repo_name: BuildJet/siconos

max_issues_repo_name: BuildJet/siconos

max_forks_repo_name: BuildJet/siconos


text:
/******************************************************************************
 * myrng.h: header for abstract random number generator                       *
 * Author: Yan Y. Liu <yanliu@illinois.edu>                                   *
 * Date: 2014/08/17                                                           *
 * Copyright and license of this source file are specified in LICENSE.TXT     *
 * under the root directory of this software package.                         *
 ******************************************************************************/
#ifndef MYRNG_H
#define MYRNG_H

// macro: random number generator
#ifdef GSL_SPRNG
#include <gsl/gsl_rng.h>
#include "gsl-sprng.h"
#include <gsl/gsl_randist.h>
#define MYRANDI(kkk) ((int)(gsl_rng_uniform_int(gsl_sprng_r, (unsigned long int)(kkk))))
#define MYRANDF() (gsl_rng_uniform(gsl_sprng_r))
#else
#ifdef SPRNG
#include "mysprng.h"
#define MYRANDI(kkk) ((int)(mysprng() * (kkk)))
#define MYRANDF() (mysprng())
#else
#define MYRANDI(kkk) ((int)(drand48() * (kkk)))
#define MYRANDF() (drand48())
#endif
#endif

#endif


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

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>

#include "gsl_linalg.h"

#define REAL double

/* Factorise a general N x N matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.
 *
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation.
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (p->size != A->size1)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i, j, k;

      *signum = 1;
      gsl_permutation_init (p);

      for (j = 0; j < N - 1; j++)
	{
	  /* Find maximum in the j-th column */

	  REAL ajj, max = fabs (gsl_matrix_get (A, j, j));
	  size_t i_pivot = j;

	  for (i = j + 1; i < N; i++)
	    {
	      REAL aij = fabs (gsl_matrix_get (A, i, j));

	      if (aij > max)
		{
		  max = aij;
		  i_pivot = i;
		}
	    }

	  if (i_pivot != j)
	    {
	      gsl_matrix_swap_rows (A, j, i_pivot);
	      gsl_permutation_swap (p, j, i_pivot);
	      *signum = -(*signum);
	    }

	  ajj = gsl_matrix_get (A, j, j);

	  if (ajj != 0.0)
	    {
	      for (i = j + 1; i < N; i++)
		{
		  REAL aij = gsl_matrix_get (A, i, j) / ajj;
		  gsl_matrix_set (A, i, j, aij);

		  for (k = j + 1; k < N; k++)
		    {
		      REAL aik = gsl_matrix_get (A, i, k);
		      REAL ajk = gsl_matrix_get (A, j, k);
		      gsl_matrix_set (A, i, k, aik - aij * ajk);
		    }
		}
	    }
	}

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve for x */

      gsl_linalg_LU_svx (LU, p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* Apply permutation to RHS */

      gsl_permute_vector (p, x);

      /* Solve for c using forward-substitution, L c = P b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_refine (const gsl_matrix * A, const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Compute residual, residual = (A * x  - b) */

      gsl_vector_memcpy (residual, b);
      gsl_blas_dgemv (CblasNoTrans, 1.0, A, x, -1.0, residual);

      /* Find correction, delta = - (A^-1) * residual, and apply it */

      gsl_linalg_LU_svx (LU, p, residual);
      gsl_blas_daxpy (-1.0, residual, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)
{
  size_t i, n = LU->size1;

  int status = GSL_SUCCESS;

  gsl_matrix_set_identity (inverse);

  for (i = 0; i < n; i++)
    {
      gsl_vector_view c = gsl_matrix_column (inverse, i);
      int status_i = gsl_linalg_LU_svx (LU, p, &(c.vector));

      if (status_i)
	status = status_i;
    }

  return status;
}

double
gsl_linalg_LU_det (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  double det = (double) signum;

  for (i = 0; i < n; i++)
    {
      det *= gsl_matrix_get (LU, i, i);
    }

  return det;
}

double
gsl_linalg_LU_lndet (gsl_matrix * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      lndet += log (fabs (gsl_matrix_get (LU, i, i)));
    }

  return lndet;
}

int
gsl_linalg_LU_sgndet (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  int s = signum;

  for (i = 0; i < n; i++)
    {
      double u = gsl_matrix_get (LU, i, i);

      if (u < 0)
	{
	  s *= -1;
	}
      else if (u == 0)
	{
	  s = 0;
	  break;
	}
    }

  return s;
}


text:
#include "stdio.h"
#include "gsl/gsl_integration.h"
#include "gsl/gsl_errno.h"
#include "limber.h"
#include <gsl/gsl_interp2d.h>
//#include <gsl/gsl_spline2d.h>
//#include <gsl/gsl_spline.h>

// This is a workspace size for the gsl integrator
#define LIMBER_FIXED_TABLE_SIZE 4096

// data that is passed into the integrator
// This is everything we need to compute the
// integrand
typedef struct IntegrandData{
	double chimin;
	double chimax;
	double ell;
	gsl_spline * WbX; //K1
	gsl_spline * WfX; //K2
	gsl_spline * WmX; //K3
	gsl_spline * WbY; //K4
	gsl_spline * WfY; //K5
	gsl_spline * WmY; //K6
	Interpolator2D * P;
	Interpolator2D * f; //for the f(k,z(x))
	Interpolator2D * D; //for the D(k,z(x))
	Interpolator2D * BB; //for the b(k,z(x))
	gsl_interp_accel * accelerator_xb;
	gsl_interp_accel * accelerator_yb;
	gsl_interp_accel * accelerator_xf;
	gsl_interp_accel * accelerator_yf;
	gsl_interp_accel * accelerator_xm;
	gsl_interp_accel * accelerator_ym;
} IntegrandData;

// the integrand for all
static double integrand(double chi, void * data_void)
{
	IntegrandData * data = (IntegrandData*) data_void;
	// Return 0 if outside range, for convenience.
	// Important to ensure that ranges are wide enough.
	if(chi < data->chimin || chi > data->chimax) return 0.0;

//exactly the same notation of paperIII

	double wx0 = gsl_spline_eval(data->WbX,chi,data->accelerator_xb);
	double wy0 = gsl_spline_eval(data->WbY,chi,data->accelerator_yb);

	double wx0m = gsl_spline_eval(data->WmX,chi,data->accelerator_xm);
	double wy0m = gsl_spline_eval(data->WmY,chi,data->accelerator_ym);

	double wx1 = gsl_spline_eval(data->WfX,chi,data->accelerator_xf);
	double wy1 = gsl_spline_eval(data->WfY,chi,data->accelerator_yf);

	double wx2 = gsl_spline_eval(data->WfX,chi*(2.0*data->ell-3.0)/(2.0*data->ell+1.0),data->accelerator_xf);
	double wy2 = gsl_spline_eval(data->WfY,chi*(2.0*data->ell-3.0)/(2.0*data->ell+1.0),data->accelerator_yf);

	double chi5 = chi*(2.0*data->ell+5.0)/(2.0*data->ell+1.0) < data->chimax ? chi*(2.0*data->ell+5.0)/(2.0*data->ell+1.0) : 0.0;
	double wx3 = gsl_spline_eval(data->WfX, chi5, data->accelerator_xf);
	double wy3 = gsl_spline_eval(data->WfY, chi5, data->accelerator_yf);

        double c1 = (2.*data->ell*data->ell+2.*data->ell-1.)/((2.*data->ell-1.)*(2.*data->ell+3.));
        double c2 = -((data->ell*(data->ell-1.))/((2.*data->ell-1.)*sqrt((2.*data->ell+1.)*(2.*data->ell-3.))));
        double c3 = -(((data->ell+1.)*(data->ell+2.))/((2.*data->ell+3.)*sqrt((2.*data->ell+1.)*(2.*data->ell+5.))));

	// Get P(k,0) using k=ell/chi.
	// The interp_2d interpolator returns 0 if either
	// parameter is outside its range
	double ka = (data->ell+0.5) / chi;
	double pa = interp_2d(ka, chi, data->P);

        //for the scale dependent bias b(k,z)
        double b0 = interp_2d(ka, chi, data->BB);

        //for the scale dependent growth factor D(k,z)
        double d0 = interp_2d(ka, chi, data->D);
        double d1 = d0;
	double d2 = interp_2d(ka, chi*(2.0*data->ell-3.0)/(2.0*data->ell+1.0), data->D);
	double d3 = interp_2d(ka, chi5, data->D);

        //for the scale dependent growth rate f(k,z)
        double f1 = interp_2d(ka, chi, data->f);
	double f2 = interp_2d(ka, chi*(2.0*data->ell-3.0)/(2.0*data->ell+1.0), data->f);
	double f3 = interp_2d(ka, chi5, data->f);

        double WWX= (wx0*b0*d0 + wx0m*d0 + c1*wx1*f1*d1 + c2*wx2*f2*d2 + c3*wx3*f3*d3);
        double WWY= (wy0*b0*d0 + wy0m*d0 + c1*wy1*f1*d1 + c2*wy2*f2*d2 + c3*wy3*f3*d3);

	double result = WWX * WWY * pa / chi / chi;

	return result;

}

// These two convenience functions
// peer into the internals of the gsl_spline.
// This is probably a bit naughty, since they could
// in theory change the internals.
static double inline limber_gsl_spline_min_x(gsl_spline * s)
{
	return s->x[0];
}
static double inline limber_gsl_spline_max_x(gsl_spline * s)
{
	return s->x[s->size-1];
}

double get_kernel_peak(gsl_spline * WbX, gsl_spline * WbY, int n_chi)
{
  double chimin_x = limber_gsl_spline_min_x(WbX);
  double chimin_y = limber_gsl_spline_min_x(WbY);
  double chimax_x = limber_gsl_spline_max_x(WbX);
  double chimax_y = limber_gsl_spline_max_x(WbY);
  double chimin = chimin_x>chimin_y ? chimin_x : chimin_y;
  double chimax = chimax_x<chimax_y ? chimax_x : chimax_y;
  double dchi = (chimax - chimin)/n_chi;
  double chi_peak = chimin;
  double chi;
  double kernel_val=0.;
  double kernel_peak=0.;
  for (int i_chi=0; i_chi<=n_chi; i_chi++){
    chi=chimin+i_chi*dchi;
    kernel_val = gsl_spline_eval(WbX,chi,NULL) * gsl_spline_eval(WbY,chi,NULL) / chi / chi;
    if (kernel_val>kernel_peak){
      kernel_peak = kernel_val;
      chi_peak = chi;
    }
  }
  // printf("chi_peak = %f\n",chi_peak);
  return chi_peak;
}

gsl_integration_workspace * W = NULL;
gsl_integration_glfixed_table *table = NULL;

void setup_integration_workspaces(){
	if (W==NULL){
		W = gsl_integration_workspace_alloc(LIMBER_FIXED_TABLE_SIZE);
	}
	if (table==NULL){
		table = gsl_integration_glfixed_table_alloc((size_t) LIMBER_FIXED_TABLE_SIZE);
	}
}

static
void limber_gsl_fallback_integrator(gsl_function * F, double chimin, double chimax,
	double abstol, double reltol, double * c_ell, double * error){

	// Only one warning per process
	static int fallback_warning_given = 0;

	// Deactivate error handling - if this one fails we will fall back to a more reliable but slower integrator
	gsl_error_handler_t * old_handler = gsl_set_error_handler_off();

	// Try the fast but flaky integrator.
	int status = gsl_integration_qag(F, chimin, chimax, abstol, reltol, LIMBER_FIXED_TABLE_SIZE, GSL_INTEG_GAUSS61, W, c_ell, error);

	// Restore the old error handler
	gsl_set_error_handler(old_handler);

	// If the fast integrator failed fall back to the old one.
	if (status){
		IntegrandData * data = (IntegrandData*) F->params;
		double ell = data->ell;
		if (fallback_warning_given==0){
			fprintf(stderr, "Falling back to the old integrator for ell=%lf (status=%d)\n", ell,status);
			fallback_warning_given=1;
		}
		*c_ell = gsl_integration_glfixed(F,chimin,chimax,table);
	}

}

// The only function in this little library callable from the outside
// world.  The limber_config structure is defined in limber.h but is fairly
// obvious.  The splines and the interpolator need to be functions of
// chi NOT z.
gsl_spline * limber_integral(limber_config * config, gsl_spline * WbX,  gsl_spline * WfX,  gsl_spline * WmX, gsl_spline * WbY,
	                 gsl_spline * WfY, gsl_spline * WmY, Interpolator2D * P, Interpolator2D * f, Interpolator2D * D, Interpolator2D * BB)

{

    config->status = LIMBER_STATUS_ERROR;
    int any_parameter_error=0;
    if (WbX==NULL){
        fprintf(stderr, "NULL WbX parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WbY==NULL){
        fprintf(stderr, "NULL WbY parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WfX==NULL){
        fprintf(stderr, "NULL WfX parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WfY==NULL){
        fprintf(stderr, "NULL WfY parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WmX==NULL){
        fprintf(stderr, "NULL WmX parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WmY==NULL){
        fprintf(stderr, "NULL WmY parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (P==NULL){
        fprintf(stderr, "NULL P parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (f==NULL){
        fprintf(stderr, "NULL f parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (D==NULL){
        fprintf(stderr, "NULL D parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (config->n_ell<0){
        fprintf(stderr, "Negative n_ell parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (config->n_ell==0){
        fprintf(stderr, "Error: n_ell=0 in limber calculation.\n");
        any_parameter_error = 1;
    }
    if (any_parameter_error){
        return NULL;
    }

	config->status = LIMBER_STATUS_OK;

	// Get the appropriate ranges over which to integrate
	// It is assumed that (at least one of) the kernel
	// splines should go to zero in some kind of reasonable
	// place, so we just use the range they specify
	IntegrandData data;
	double chimin_xb = limber_gsl_spline_min_x(WbX);
	double chimin_yb = limber_gsl_spline_min_x(WbY);
	double chimax_xb = limber_gsl_spline_max_x(WbX);
	double chimax_yb = limber_gsl_spline_max_x(WbY);

	double chimin_xf = limber_gsl_spline_min_x(WfX);
	double chimin_yf = limber_gsl_spline_min_x(WfY);
	double chimax_xf = limber_gsl_spline_max_x(WfX);
	double chimax_yf = limber_gsl_spline_max_x(WfY);

	double chimin_xm = limber_gsl_spline_min_x(WmX);
	double chimin_ym = limber_gsl_spline_min_x(WmY);
	double chimax_xm = limber_gsl_spline_max_x(WmX);
	double chimax_ym = limber_gsl_spline_max_x(WmY);

        //my ifs here
        double chimin_x = chimin_xb > chimin_xf ? chimin_xb : chimin_xf;
        double chimin_y = chimin_yb > chimin_yf ? chimin_yb : chimin_yf;
        double chimax_x = chimax_xb < chimax_xf ? chimax_xb : chimax_xf;
        double chimax_y = chimax_yb < chimax_yf ? chimax_yb : chimax_yf;
	double c_ell, error;

	// Workspaces for the main and falback integrators.
	// Static, so only allocated once as it has a fixed size.
	setup_integration_workspaces();

	double reltol = config->relative_tolerance;
	double abstol = config->absolute_tolerance;
	// double reltol = 0.001;
	// double abstol = 0.00001;
	// printf("TOLS: %le %le\n",reltol,abstol);

	// Take the smallest range since we want both the
	// splines to be valid there.
	// This range as well as all the data needed to compute
	// the integrand is put into a struct to be passed
	// through the integrator to the function above.
	data.chimin = chimin_x>chimin_y ? chimin_x : chimin_y;
	data.chimax = chimax_x<chimax_y ? chimax_x : chimax_y;
	data.WbX = WbX;
	data.WfX = WfX;
	data.WmX = WmX;
	data.WbY = WbY;
	data.WfY = WfY;
	data.WmY = WmY;
	data.P = P;
	data.f = f;//new
	data.D = D;//new
	data.BB = BB;//new
	data.accelerator_xb = gsl_interp_accel_alloc();
	data.accelerator_xf = gsl_interp_accel_alloc();
	data.accelerator_xm = gsl_interp_accel_alloc();
	data.accelerator_yb = gsl_interp_accel_alloc();
	data.accelerator_yf = gsl_interp_accel_alloc();
	data.accelerator_ym = gsl_interp_accel_alloc();

	// Set up the workspace and inputs to the integrator.
	// Not entirely sure what the table is.
	gsl_function F;
	F.function = integrand;
	F.params = &data;

	//gsl_integration_workspace * workspace = gsl_integration_workspace_alloc(2048);

	// results of the integration go into these arrays.
	double c_ell_vector[config->n_ell];
	double ell_vector[config->n_ell];

	// loop through ell values according to the input configuration
	for (int i_ell = 0; i_ell<config->n_ell; i_ell++){
		double ell = config->ell[i_ell];
		data.ell=ell;

		// Perform the main integration.
		// This particular function is used because that's what Matt Becker
		// found to work best.
		//c_ell = gsl_integration_glfixed(&F,data.chimin,data.chimax,table);
		// New function still attributable to the legacy of Matt Becker's integrator wisdom.
		// gsl_integration_qag(&F, data.chimin, data.chimax, abstol, reltol, LIMBER_FIXED_TABLE_SIZE, GSL_INTEG_GAUSS61, W, table, &c_ell, &error);
		//printf("%d %f %f\n",i_ell,c_ell_old,c_ell);
		limber_gsl_fallback_integrator(&F, data.chimin, data.chimax,
			abstol, reltol, &c_ell, &error);

		//Include the prefactor scaling
		c_ell *= config->prefactor;

		// Record the results into arrays
		c_ell_vector[i_ell] = c_ell;
		ell_vector[i_ell] = ell;
	}

		// It is often useful to interpolate into the logs of the functions
		// This is optional in the config. We move this outside the main loop
		// since we may have all zeros in the output
		if (config->xlog) {
			for (int i_ell = 0; i_ell<config->n_ell; i_ell++){
				ell_vector[i_ell] = log(ell_vector[i_ell]);
			}
		}
		if (config->ylog){
			for (int i_ell = 0; i_ell<config->n_ell; i_ell++){
				if (c_ell_vector[i_ell]<0){
					config->status = LIMBER_STATUS_NEGATIVE;
				}
					// negative is worse than zero so only set to zero it not already negative
				else if ((c_ell_vector[i_ell]==0) && (config->status<LIMBER_STATUS_ZERO)){
					config->status = LIMBER_STATUS_ZERO;
				}
			}
			// If none of the values are <= 0 then we are okay to go ahead and take the logs.
			if (config->status == LIMBER_STATUS_OK){
				for (int i_ell = 0; i_ell<config->n_ell; i_ell++) c_ell_vector[i_ell] = log(c_ell_vector[i_ell]);
			}

		}

	// Create a spline of the arrays as the output
	gsl_spline * output = gsl_spline_alloc(gsl_interp_akima, (size_t) config->n_ell);
	gsl_spline_init(output, ell_vector, c_ell_vector, (size_t) config->n_ell);

	// Tidy up
	gsl_interp_accel_free(data.accelerator_xb);
	gsl_interp_accel_free(data.accelerator_yb);
	gsl_interp_accel_free(data.accelerator_xf);
	gsl_interp_accel_free(data.accelerator_yf);
	gsl_interp_accel_free(data.accelerator_xm);
	gsl_interp_accel_free(data.accelerator_ym);

	// These two are not deallocated because they are static and only initialized once.
	// gsl_integration_glfixed_table_free(table);
	// gsl_integration_workspace_free(W);

	// And that's it
	return output;
}


text:
/* ################################################################################## */
/* ###                                                                            ### */
/* ###                                 Gadgetmp2                                  ### */
/* ###                                                                            ### */
/* ###   Original: Gadget2 in the version used in Amuse                           ### */
/* ###   Author: Gadget2 and Amuse contributors                                   ### */
/* ###                                                                            ### */
/* ###   Modified: July 2020                                                      ### */
/* ###   Author: Thomas Schano                                                    ### */
/* ###                                                                            ### */
/* ###   Changes are intended to enable precise calculations in                   ### */
/* ###   non periodic small domain simulations in which comoving parts            ### */
/* ###   are simulated in std precision                                           ### */
/* ###                                                                            ### */
/* ################################################################################## */
#include "src/tags.hpp"
#include <gsl/gsl_rng.h>
#include "src/proto.hpp"

typedef struct {
    my_float mass;                                        /// mass
    my_float x, y, z;                                     /// position
    my_float vx, vy, vz;                                  /// velocity
    my_float radius;                                      /// radius
} dynamics_state;

typedef struct {
    my_float mass;                                        /// mass
    my_float x, y, z;                                     /// position
    my_float vx, vy, vz;                                  /// velocity
    my_float u;                                           /// entropy
#ifdef MORRIS97VISC
    my_float alpha, dalphadt;				///viscosity
#endif
} sph_state;

//void   gadgetmp2::begrun(void);
//double gadgetmp2::second(void);

int found_particle(int index_of_the_particle, int *local_index);
void update_particle_map(void);

/*void gadgetmp2::hydro_state_at_point(double pos[3], double vel[3], double *h_out,
  double *ngb_out, double *dhsml_out, double *rho_out, double *rhov_out,
  double *rhov2_out, double *rhoe_out);
*/

max_stars_repo_name: GFTwrt/amuse

max_issues_repo_name: GFTwrt/amuse

max_forks_repo_name: GFTwrt/amuse


text:
//! \file heat_driver.h
//! Driver for Magnolia's heat transfer solver
#ifndef ENRICO_SURROGATE_HEAT_DRIVER_H
#define ENRICO_SURROGATE_HEAT_DRIVER_H

#include "enrico/geom.h"
#include "enrico/heat_fluids_driver.h"

#include <gsl/gsl>
#include <mpi.h>
#include <pugixml.hpp>
#include <xtensor/xtensor.hpp>

#include <cstddef>

namespace enrico {

//! Struct containing geometric information for a flow channel
struct Channel {
  //! Channel index
  int index_;

  //! Channel flow area
  double area_;

  //! Vector of rod IDs connected to this channel, all with a fractional perimeter
  //! in contact with the channel equal to 0.25
  std::vector<std::size_t> rod_ids_;
};

//! Struct containing geometry information for a cylindrical solid rod
struct Rod {
  //! Rod index
  int index_;

  //! Rod cladding outer radius
  double clad_outer_radius_;

  //! Rod cladding inner radius
  double clad_inner_radius_;

  //! Rod pellet radius
  double pellet_radius_;

  //! Vector of channel IDs connected to this rod, all with a fractional perimeter
  //! in contact with the rod equal to 0.25
  std::vector<std::size_t> channel_ids_;
};

//! Class to construct flow channels for a Cartesian lattice of pins
class ChannelFactory {
public:
  ChannelFactory(double pitch, double rod_radius)
    : pitch_(pitch)
    , radius_(rod_radius)
    , interior_area_(pitch_ * pitch_ - M_PI * radius_ * radius_)
  {}

  //! Make a corner subchannel connected to given rods
  Channel make_corner(const std::vector<std::size_t>& rods) const
  {
    Channel c;
    c.index_ = index_++;
    c.area_ = 0.25 * interior_area_;
    c.rod_ids_ = rods;
    return c;
  }

  //! Make an edge subchannel connected to given rods
  Channel make_edge(const std::vector<std::size_t>& rods) const
  {
    Channel c;
    c.index_ = index_++;
    c.area_ = 0.5 * interior_area_;
    c.rod_ids_ = rods;
    return c;
  }

  //! Make an interior subchannel connected to given rods
  Channel make_interior(const std::vector<std::size_t>& rods) const
  {
    Channel c;
    c.index_ = index_++;
    c.area_ = interior_area_;
    c.rod_ids_ = rods;
    return c;
  }

private:
  //! rod pitch
  double pitch_;

  //! rod outer radius
  double radius_;

  //! interior channel flow area, which is proportional to flow areas for all
  //! other channel types
  double interior_area_;

  //! index of constructed channel
  static int index_;
};

class RodFactory {
public:
  RodFactory(double clad_OR, double clad_IR, double pellet_OR)
    : clad_outer_r_(clad_OR)
    , clad_inner_r_(clad_IR)
    , pellet_outer_r_(pellet_OR)
  {}

  //! Make a rod connected to given channels
  Rod make_rod(const std::vector<std::size_t>& channels) const
  {
    Rod r;
    r.index_ = index_++;
    r.clad_inner_radius_ = clad_inner_r_;
    r.clad_outer_radius_ = clad_outer_r_;
    r.pellet_radius_ = pellet_outer_r_;
    r.channel_ids_ = channels;
    return r;
  }

private:
  //! Cladding outer radius
  double clad_outer_r_;

  //! Cladding inner radius
  double clad_inner_r_;

  //! Pellet outer radius
  double pellet_outer_r_;

  //! Index of constructed rod
  static int index_;
};

/**
 * Class providing surrogate thermal-hydraulic solution for a Cartesian
 * bundle of rods with upwards-flowing coolant. A conduction model is used
 * for the solid phase, with axial conduction neglected. The solid phase is
 * linked to the fluid phase by conjugate heat transfer, which is treated
 * here with a pseudo-steady-state approach where the power entering the fluid
 * matches the power in the rod at that axial elevation. It is assumed that
 * there is zero thermal resistance between the rod and the fluid.
 *
 * The fluid solution is obtained with a very simplified "subchannel" method
 * that neglects all crossflow terms between channels such that the method
 * is more akin to a single-pin analysis in a coolant-centered basis. The
 * enthalpy is solved by simply axial energy balance, while the axial momentum
 * equation is solved for pressure (the mass flow rate in each channel being
 * fixed) while neglecting friction effects.
 */
class SurrogateHeatDriver : public HeatFluidsDriver {
public:
  //! Initializes heat-fluids surrogate with the given MPI communicator.
  //!
  //! \param comm  The MPI communicator used to initialze the surrogate
  //! \param node  XML node containing settings for surrogate
  explicit SurrogateHeatDriver(MPI_Comm comm, pugi::xml_node node);

  //! Verbosity options for printing simulation results
  enum class verbose { NONE, LOW, HIGH };

  bool has_coupling_data() const final { return comm_.rank == 0; }

  //! Get the number of local mesh elements
  //! \return Number of local mesh elements
  int n_local_elem() const override;

  //! Get the number of global mesh elements
  //! \return Number of global mesh elements
  std::size_t n_global_elem() const override;

  //! Set the heat source for a given local element
  //!
  //! \param local_elem A local element ID
  //! \param heat A heat source term
  //! \return Error code
  int set_heat_source_at(int32_t local_elem, double heat) override;

  //! Solves the heat-fluids surrogate solver
  void solve_step() final;

  void solve_heat();

  void solve_fluid();

  //! Returns Number of rings in fuel and clad
  std::size_t n_rings() const { return n_fuel_rings_ + n_clad_rings_; }

  //! Returns cladding inner radius
  double clad_inner_radius() const { return clad_inner_radius_; }

  //! Returns cladding outer radius
  double clad_outer_radius() const { return clad_outer_radius_; }

  //! Returns pellet outer radius
  double pellet_radius() const { return pellet_radius_; }

  //! Returns number of fuel rings
  std::size_t n_fuel_rings() const { return n_fuel_rings_; }

  //! Returns number of clad rings
  std::size_t n_clad_rings() const { return n_clad_rings_; }

  //! Returns number of pins in x-direction
  std::size_t n_pins_x() const { return n_pins_x_; }

  //! Returns number of pins in y-direction
  std::size_t n_pins_y() const { return n_pins_y_; }

  //! Returns pin pitch
  double pin_pitch() const { return pin_pitch_; }

  //! Returns inlet temperature boundary condition in [K]
  double inlet_temperature() const { return inlet_temperature_; }

  //! Returns inlet mass flowrate boundary condition in [kg/s]
  double mass_flowrate() const { return mass_flowrate_; }

  //! Returns maximum number of subchannel iterations
  std::size_t max_subchannel_its() const { return max_subchannel_its_; }

  //! Returns subchannel convergence tolerance for enthalpy
  double subchannel_tol_h() const { return subchannel_tol_h_; }

  //! Returns subchannel convergence tolerance for pressure
  double subchannel_tol_p() const { return subchannel_tol_p_; }

  //! Returns convergence tolerance for solid energy equation
  double heat_tol() const { return heat_tol_; }

  //! Write data to VTK
  void write_step(int timestep, int iteration) final;

  //! Returns solid temperature in [K] for given region
  double solid_temperature(std::size_t pin, std::size_t axial, std::size_t ring) const;

  //! Returns fluid density in [g/cm^3] for given region
  double fluid_density(std::size_t pin, std::size_t axial) const;

  //! Returns fluid temperature in [K] for given region
  double fluid_temperature(std::size_t pin, std::size_t axial) const;

  // Data on fuel pins
  xt::xtensor<double, 2> pin_centers_; //!< (x,y) values for center of fuel pins
  xt::xtensor<double, 1> z_;           //!< Bounding z-values for axial segments
  std::size_t n_axial_;                //!< number of axial segments
  std::size_t n_azimuthal_{4};         //!< number of azimuthal segments

  //! Total number of pins
  std::size_t n_pins_;

  // Dimensions for a single fuel pin axial segment
  double clad_outer_radius_;     //!< clad outer radius in [cm]
  double clad_inner_radius_;     //!< clad inner radius in [cm]
  double pellet_radius_;         //!< fuel pellet radius in [cm]
  std::size_t n_fuel_rings_{20}; //!< number of fuel rings
  std::size_t n_clad_rings_{2};  //!< number of clad rings

  //!< Channels in the domain
  std::vector<Channel> channels_;

  //!< Rods in the domain
  std::vector<Rod> rods_;

  //! Mass flowrate for coolant-centered channels; this is determine by distributing
  //! a total inlet mass flowrate among the channels based on the fractional flow area.
  xt::xtensor<double, 1> channel_flowrates_;

  // solver variables and settings
  xt::xtensor<double, 4>
    source_; //!< heat source for each (pin, axial segment, ring, azimuthal segment)
  xt::xtensor<double, 1> r_grid_clad_; //!< radii of each clad ring in [cm]
  xt::xtensor<double, 1> r_grid_fuel_; //!< radii of each fuel ring in [cm]

  //! Cross-sectional areas of rings in fuel and cladding
  xt::xtensor<double, 1> solid_areas_;

  // visualization
  std::string viz_basename_{
    "heat_surrogate"}; //!< base filename for visualization files (default: magnolia)
  std::string viz_iterations_{
    "none"};                    //!< visualization iterations to write (none, all, final)
  std::string viz_data_{"all"}; //!< visualization data to write
  std::string viz_regions_{"all"}; //!< visualization regions to write
  size_t vtk_radial_res_{20};      //!< radial resolution of resulting vtk files

private:
  //! Get temperature of local mesh elements
  //! \return Temperature of local mesh elements in [K]
  std::vector<double> temperature_local() const override;

  //! Get density of local mesh elements
  //! \return Density of local mesh elements in [g/cm^3]
  std::vector<double> density_local() const override;

  //! States whether each local region is in fluid
  //! \return For each local region, 1 if region is in fluid and 0 otherwise
  std::vector<int> fluid_mask_local() const override;

  //! Get centroids of local mesh elements
  //! \return Centroids of local mesh elements
  std::vector<Position> centroid_local() const override;

  //! Get volumes of local mesh elements
  //! \return Volumes of local mesh elements
  std::vector<double> volume_local() const override;

  //! Create internal arrays used for heat equation solver
  void generate_arrays();

  //! Channel index in terms of row, column index
  int channel_index(int row, int col) const { return row * (n_pins_x_ + 1) + col; }

  //! Rod power at a given node in a given pin, computed by integrating the heat source
  //! (assumed constant in each ring) over the pin.
  //! \param pin   pin index
  //! \param axial axial index
  double rod_axial_node_power(const int pin, const int axial) const;

  //! Diagnostic function to assess whether the mass is conserved by the subchannel
  //! solver by comparing the mass flowrate in each axial plane (at cell-centered
  //! positions) to the specified inlet mass flowrate.
  //! \param rho density in a cell-centered basis
  //! \param u   axial velocity in a face-centered basis
  bool is_mass_conserved(const xt::xtensor<double, 2>& rho,
                         const xt::xtensor<double, 2>& u) const;

  //! Diagnostic function to assess whether the energy is conserved by the subchannel
  //! solver by comparing the energy deposition in each channel in each axial plane
  //! (at cell-centered positions) to the powers of the rods connected to that channel.
  //! \param rho density in a cell-centered basis
  //! \param u   axial velocity in a face-centered basis
  //! \param h   enthalpy in a face-centered basis
  //! \param q   powers in each channel in a cell-centered basis
  bool is_energy_conserved(const xt::xtensor<double, 2>& rho,
                           const xt::xtensor<double, 2>& u,
                           const xt::xtensor<double, 2>& h,
                           const xt::xtensor<double, 2>& q) const;

  //!< solid temperature in [K] for each (pin, axial segment, ring)
  xt::xtensor<double, 3> solid_temperature_;

  //! Flow areas for coolant-centered channels
  xt::xtensor<double, 1> channel_areas_;

  //! Fluid temperature in a rod-centered basis indexed by rod ID and axial ID
  xt::xtensor<double, 2> fluid_temperature_;

  //! Fluid density in [g/cm^3] in a rod-centered basis indexed by rod ID and axial ID
  xt::xtensor<double, 2> fluid_density_;

  //! Number of pins in the x-direction in a Cartesian grid
  std::size_t n_pins_x_;

  //! Number of pins in the y-direction in a Cartesian grid
  std::size_t n_pins_y_;

  //! Pin pitch, assumed the same for the x and y directions
  double pin_pitch_;

  //! Inlet fluid temperature [K]
  double inlet_temperature_;

  //! Mass flowrate of fluid into the domain [kg/s]
  double mass_flowrate_;

  //! Number of channels
  std::size_t n_channels_;

  //! Maximum number of iterations for subchannel solution, set to a default value
  //! of 100 if not set by the user
  int max_subchannel_its_ = 100;

  //! Convergence tolerance on enthalpy for the subchannel solution for use in
  //! convergence based on the L-1 norm, set to a default value of 1e-2
  double subchannel_tol_h_ = 1e-2;

  //! Convergence tolerance on pressure for the subchannel solution for use in
  //! convergence based on the L-1 norm, set to a default value of 1e-2
  double subchannel_tol_p_ = 1e-2;

  //! Convergence tolerance for solid temperature solution, set to a default value
  //! of 1e-4
  double heat_tol_ = 1e-4;

  //! Gravitational acceleration
  const double g_ = 9.81;

  //! Verbosity setting for printing simulation results; defaults to NONE
  verbose verbosity_ = verbose::NONE;

}; // end SurrogateHeatDriver

} // namespace enrico

#endif // ENRICO_SURROGATE_HEAT_DRIVER_H

max_stars_repo_name: pshriwise/enrico

max_issues_repo_name: pshriwise/enrico

max_forks_repo_name: lebuller/enrico


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <alloca.h>

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_statistics_float.h>

#include "rank.h"
#include "stattest.h"
#include "num.h"

typedef float con_t;

struct ConCovars {

	/**
	  * Maximum number of samples
	  */
	int SAMPLE_CAPACITY;

	/**
	  * Count of samples pushed since last con_clear.
	  */
	int sample_count;

	/**
	  * The total amount of malloc'ed space (for the purposes
	  * of fast clearing with memset).
	  */
	size_t SIZEOF_BUFFERS;

	/**
	  * Two buffers
	  */
	con_t *l, *r;

	void *rank_scratch;
};

#if defined(_UNITTEST_NUM_)
static void dbg_dump( struct ConCovars *co, FILE *fp ) {
	for(unsigned int i = 0; i < co->sample_count; i++ )
		fprintf( fp, "%f\t%f\n", co->l[i], co->r[i] );
}
#endif

void con_destroy( void *pv ) {

	if( pv ) {
		struct ConCovars *co = (struct ConCovars *)pv;
		if( co->rank_scratch )
			rank_free( co->rank_scratch );
		if( co->l )
			free( co->l );
		free( pv );
	}
}

/**
  * Pre-allocate a set of working buffers large enough for all anticipated
  * calculations (max feature length) and a struct to wrap them.
  */
void *con_create( unsigned int cap ) {
	struct ConCovars *co
		= calloc( 1, sizeof(struct ConCovars) );
	if( co ) {
		co->SAMPLE_CAPACITY = cap;
		co->SIZEOF_BUFFERS  = 2*cap*sizeof(con_t);
		// Allocate one large buffer and partition it up.
		co->l = calloc( co->SIZEOF_BUFFERS, sizeof(char) );
		co->r = co->l + cap;
		co->rank_scratch = rank_alloc( cap );
		// If -anything- failed clean up any successes.
		if( (NULL == co->l) ||
			(NULL == co->rank_scratch) ) {
			con_destroy( co );
			return NULL;
		}
		return co;
	}
	return NULL;
}

void con_clear( void *pv ) {
	struct ConCovars *co = (struct ConCovars *)pv;
	co->sample_count = 0;
	memset( co->l, 0, co->SIZEOF_BUFFERS );
}

void con_push( void *pv, float n1, float n2 ) {
	struct ConCovars *co = (struct ConCovars *)pv;
	const int i = co->sample_count++;
	assert( i <= co->SAMPLE_CAPACITY );
	co->l[i] = n1;
	co->r[i] = n2;
}

size_t con_size( void *pv ) {
	return ((struct ConCovars *)pv)->sample_count;
}

bool con_complete( void *pv ) {
	return ((struct ConCovars *)pv)->sample_count > 2;
	// ...otherwise p-value computation will fail.
}

/**
 */
int con_spearman_correlation( void *pv, struct Statistic *result ) {

	struct ConCovars *co = (struct ConCovars *)pv;
	const int N = co->sample_count;

	assert( N > 2 );
	assert( NULL != co->rank_scratch );

	const int rinfo1
		= rank_floats( co->l, N, 0, co->rank_scratch );
	const int rinfo2
		= rank_floats( co->r, N, 0, co->rank_scratch );

	if( RANK_STATUS_CONST & rinfo1 ) // vectors were in fact constant!
		result->extra_value[0] = N-1;
	if( RANK_STATUS_CONST & rinfo2 )
		result->extra_value[1] = N-1;

	{
		const double rho
			= gsl_stats_float_correlation( co->l, 1, co->r, 1, N );

		/**
		 * P-value computation for the correlation.
		 */

#ifdef HAVE_FISHER_TRANSFORM
		const double FisherTransform
			//= 0.5 * log( (1.+rho) / (1.-rho) );
			= atanh(rho);
		// ...absolute value to simplify CDF use below, since the
		// transformation is symmetric.
		const double z
			= sqrt( (N - 3.0) / 1.06 ) * FisherTransform;
		// ...z ~ N(0,1) under null hyp of statistical independence.
		result->name
			= "Spearman_rho,Fisher_transform";
		result->probability = gsl_cdf_ugaussian_Q( fabs(z) );
#else
		const double t
			=  fabs( rho*sqrt((N-2.0)/(1.0-rho*rho)) );
		// ...abs so that I can always test the upper tail.
		// x2 below to make it a two-tailed test. (t-distribution
		// is symmetric).
		result->name
			= "Spearman_rho,t-distribution";
		result->probability = 2* gsl_cdf_tdist_Q(t,N-2.0);
#endif
		result->value = rho;
	}
	result->sample_count = N;

	return 0;
}

#ifdef HAVE_SCALAR_PEARSON
int con_pearson_correlation( void *pv, struct Statistic *result ) {
	struct ConCovars *co = (struct ConCovars *)pv;
	const int N = co->sample_count;
	ps->rho = gsl_stats_float_correlation( co->l, 1, co->r, 1, N );
	return 0;
}
#endif

#ifdef _UNITTEST_NUM_

/**
 * This is intended to be exercised with the following R script
 * that generates a small table of grouped floating-point values
 * executes a Kruskal-Wallis test on it, and dumps the table
 * to a tab-delimited file with the test results as a comment on
 * the first line.
 *
 * x <- data.frame(
 * 		cat=as.integer(gl(3,4))-1,
 * 		num=c( runif(4)*10, runif(4)*20, 10+runif(4)*10 ) );
 * k <- with( x, kruskal.test( num, cat ) );
 * cat( sprintf( "# K=%f p-value=%f\n", k$statistic, k$p.value ), file="foo.tab" );
 * write.table( x, 'foo.tab', quote=F, sep='\t', row.names=F, col.names=F, append=TRUE );
 */
#include <err.h>

int main( int argc, char *argv[] ) {

	if( argc >= 2 ) {

		struct Statistic result;
		struct ConCovars *accum
			= con_create( atoi(argv[1]) );
		FILE *fp
			= argc > 2
			? fopen( argv[2], "r" )
			: stdin;
		char *line = NULL;
		size_t n = 0;

		con_clear( accum );

		while( getline( &line, &n, fp ) > 0 ) {
			float l, r;
			if( line[0] == '#' ) {
				fputs( line, stdout );
				continue;
			}
			if( 2 == sscanf( line, "%f\t%f\n", &l, &r ) ) {
				con_push( accum, l, r );
			} else {
				fprintf( stderr, "Failure parsing line: %s", line );
			}
		}
		free( line );
		fclose( fp );

		memset( &result, 0, sizeof(struct Statistic) );
#ifdef HAVE_SCALAR_PEARSON
		con_pearson_correlation( accum, &result );
		printf( "pearson=%f\n", result.value );
#endif
		memset( &result, 0, sizeof(struct Statistic) );
		if( con_spearman_correlation( accum, &result ) == 0 ) {
			printf( "p-value=%f, spearman=%f\n", result.probability, result.value );
		} else
			printf( "error\n" );

		con_destroy( accum );

	} else
		err( -1, "%s <sample count> [ <input file> ]", argv[0] );

	return 0;
}
#endif

max_stars_repo_name: IlyaLab/kramtools

max_issues_repo_name: IlyaLab/kramtools

max_forks_repo_name: IlyaLab/kramtools


text:
/** \file
 * This file contains an implementation of a MPMC queue by Dmitry Vyukov. It is
 * largely the same as the original source except for: minor formatting
 * changes, movement to a namespace, use of alignas(X) instead of pad variables,
 * and use of dyn_array from the guideline support library instead of
 * `new char[size]`.
 */

/*  Multi-producer/multi-consumer bounded queue.
 *  http://www.1024cores.net/home/lock-free-algorithms/queues/bounded-mpmc-queue
 *
 *  Copyright (c) 2010-2011, Dmitry Vyukov. All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *     1. Redistributions of source code must retain the above copyright notice,
 *        this list of conditions and the following disclaimer.
 *     2. Redistributions in binary form must reproduce the above copyright
 *        notice, this list of conditions and the following disclaimer in the
 *        documentation and/or other materials provided with the distribution.
 *
 *  THIS SOFTWARE IS PROVIDED BY DMITRY VYUKOV "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 DMITRY VYUKOV 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.
 *
 *  The views and conclusions contained in the software and documentation are
 *  those of the authors and should not be interpreted as representing official
 *  policies, either expressed or implied, of Dmitry Vyukov.
 */

#include <cassert>
#include <atomic>

#include <gsl_p/dyn_array.h>

namespace dvyukov {

    /**
     * A 'lockless' bounded multi-producer, multi-consumer queue
     *
     * Has the caveat that the queue can *appear* empty even if there are
     * returned items within it as a single thread can block progression
     * of the queue.
     */
    template<typename T>
    class mpmc_bounded_queue {
    public:
        /**
         * Constructs a bounded multi-producer, multi-consumer queue
         *
         * Note: Due to the algorithm used, buffer_size must be a power
         *       of two and must be greater than or equal to two.
         *
         * @param buffer_size Number of spaces available in the queue.
         */
        mpmc_bounded_queue(size_t buffer_size)
                : buffer_(buffer_size),
                  buffer_mask_(buffer_size - 1) {

            assert((buffer_size >= 2) &&
                   ((buffer_size & (buffer_size - 1)) == 0));
            for (size_t i = 0; i != buffer_size; i += 1) {
                buffer_[i].sequence_.store(i, std::memory_order_relaxed);
            }
            enqueue_pos_.store(0, std::memory_order_relaxed);
            dequeue_pos_.store(0, std::memory_order_relaxed);
        }

        /**
         * Enqueues an item into the queue
         *
         * @param data Argument to place into the array
         * @return false if the queue was full (and enqueing failed),
         *         true otherwise
         */
        bool enqueue(const T& data) {
            cell_t *cell;
            size_t pos = enqueue_pos_.load(std::memory_order_relaxed);
            for (; ;) {
                cell = &buffer_[pos & buffer_mask_];
                size_t seq = cell->sequence_.load(std::memory_order_acquire);
                intptr_t dif = (intptr_t) seq - (intptr_t) pos;
                if (dif == 0) {
                    if (enqueue_pos_.compare_exchange_weak(pos, pos + 1,
                                                           std::memory_order_relaxed)) {
                        break;
                    }
                } else if (dif < 0) {
                    return false;
                } else {
                    pos = enqueue_pos_.load(std::memory_order_relaxed);
                }
            }

            cell->data_ = data;
            cell->sequence_.store(pos + 1, std::memory_order_release);

            return true;
        }

        /**
         * Dequeues an item from the queue
         *
         * @param[out] data Reference to place item into
         * @return false if the queue was empty (and dequeuing failed),
         *         true if successful
         */
        bool dequeue(T& data) {
            cell_t *cell;
            size_t pos = dequeue_pos_.load(std::memory_order_relaxed);
            for (; ;) {
                cell = &buffer_[pos & buffer_mask_];
                size_t seq = cell->sequence_.load(std::memory_order_acquire);
                intptr_t dif = (intptr_t) seq - (intptr_t) (pos + 1);
                if (dif == 0) {
                    if (dequeue_pos_.compare_exchange_weak(pos, pos + 1,
                                                           std::memory_order_relaxed)) {
                        break;
                    }
                } else if (dif < 0) {
                    return false;
                } else {
                    pos = dequeue_pos_.load(std::memory_order_relaxed);
                }
            }

            data = cell->data_;
            cell->sequence_.store(pos + buffer_mask_ + 1,
                                  std::memory_order_release);

            return true;
        }

    private:
        struct cell_t {
            std::atomic<size_t> sequence_;
            T data_;
        };

        static const size_t cacheline_size = 64;
        typedef char cacheline_pad_t [cacheline_size];

        cacheline_pad_t pad0_;
        gsl_p::dyn_array<cell_t> buffer_;
        const size_t buffer_mask_;
        cacheline_pad_t pad1_;
        std::atomic<size_t> enqueue_pos_;
        cacheline_pad_t pad2_;
        std::atomic<size_t> dequeue_pos_;
        cacheline_pad_t pad3_;

        mpmc_bounded_queue(mpmc_bounded_queue const &);

        void operator=(mpmc_bounded_queue const &);
    };
}

max_stars_repo_name: Chippiewill/Phosphor

max_issues_repo_name: Chippiewill/Phosphor

max_forks_repo_name: Chippiewill/Phosphor


text:
#pragma once

#include <gsl/gsl>

#include <memory>
#include <atlcomcli.h>

#include "textures.h"

struct ID3D11Texture2D;
struct ID3D11DeviceContext;
struct ID3D11RenderTargetView;
struct ID3D11ShaderResourceView;
struct ID3D11DepthStencilView;

namespace gfx {

	enum class BufferFormat;

	class DynamicTexture : public Texture {
		friend class RenderingDevice;
	public:
		DynamicTexture(ID3D11DeviceContext *context,
			ID3D11Texture2D *texture,
			ID3D11ShaderResourceView *resourceView,
			const Size &size,
			uint32_t bytesPerPixel);

		int GetId() const override;

		const std::string& GetName() const;

		const ContentRect& GetContentRect() const override {
			return mContentRect;
		}

		const Size& GetSize() const override {
			return mSize;
		}

		TextureType GetType() const override {
			return TextureType::Dynamic;
		}

		// Unloads the device texture (does't prevent it from being loaded again later)
		void FreeDeviceTexture() override;

		ID3D11ShaderResourceView* GetResourceView() override {
			return mResourceView;
		}

		uint32_t GetBytesPerPixel() const {
			return mBytesPerPixel;
		}

		template<typename T>
		void Update(gsl::span<T> data) {
			UpdateRaw(
				gsl::span(reinterpret_cast<uint8_t*>(&data[0]), data.size_bytes()),
				mSize.width * sizeof(T)
			);
		}

	private:

		void UpdateRaw(gsl::span<uint8_t> data, size_t pitch);

		CComPtr<ID3D11DeviceContext> mContext;
		CComPtr<ID3D11Texture2D> mTexture;
		CComPtr<ID3D11ShaderResourceView> mResourceView;
		Size mSize;
		size_t mBytesPerPixel;
		ContentRect mContentRect;

	};

	using DynamicTexturePtr = std::shared_ptr<DynamicTexture>;

	class RenderTargetTexture : public Texture {
	friend class RenderingDevice;
	public:

		RenderTargetTexture(ID3D11Texture2D *texture,
			ID3D11RenderTargetView *rtView,
			ID3D11Texture2D *resolvedTexture,
			ID3D11ShaderResourceView *resourceView,
			const Size &size,
			bool multiSampled);

		int GetId() const override;

		const std::string& GetName() const;

		const ContentRect& GetContentRect() const override {
			return mContentRect;
		}

		const Size& GetSize() const override {
			return mSize;
		}

		TextureType GetType() const override {
			return TextureType::RenderTarget;
		}

		// Unloads the device texture (does't prevent it from being loaded again later)
		void FreeDeviceTexture() override;

		/**
		 * Only valid for MSAA targets. Will return the texture used to resolve the multisampling
		 * so it can be used as a shader resource.
		 */
		ID3D11Texture2D* GetResolvedTexture() const {
			return mResolvedTexture;
		}

		/**
		 * For MSAA targets, this will return the resolvedTexture, while for normal targets,
		 * this will just be the texture itself.
		 */
		ID3D11ShaderResourceView* GetResourceView() override {
			return mResourceView;
		}

		bool IsMultiSampled() const {
			return mMultiSampled;
		}

		gfx::BufferFormat GetFormat() const;

	private:
		CComPtr<ID3D11RenderTargetView> mRtView;
		CComPtr<ID3D11Texture2D> mTexture;
		CComPtr<ID3D11Texture2D> mResolvedTexture;
		CComPtr<ID3D11ShaderResourceView> mResourceView;
		Size mSize;
		ContentRect mContentRect;
		bool mMultiSampled;

	};

	using RenderTargetTexturePtr = std::shared_ptr<RenderTargetTexture>;

	class RenderTargetDepthStencil {
		friend class RenderingDevice;
	public:

		RenderTargetDepthStencil(ID3D11Texture2D *textureNew,
			ID3D11DepthStencilView *dsView,
			const Size &size);

		const Size& GetSize() const {
			return mSize;
		}

	private:

		CComPtr<ID3D11DepthStencilView> mDsView;
		CComPtr<ID3D11Texture2D> mTextureNew;
		Size mSize;

	};

	using RenderTargetDepthStencilPtr = std::shared_ptr<RenderTargetDepthStencil>;

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus