https://github.com/jsollari/popABC
Tip revision: 3ebe7b76043299bc4b6a541b5be0d7012895ca1c authored by Joao Sollari Lopes on 21 November 2017, 15:04:40 UTC
remove binaries folder
remove binaries folder
Tip revision: 3ebe7b7
mylib.c
/*
@author: joao lopes
@workplace: Reading University
@date: 29th May 2009
NBB - based on Mark's myutil.h
*/
#include "mylib.h"
#define repeat for(;;)
int jindic,
rand_table[98];
void opengfsr(char*path){
static int rand[] = {
-899909286,-1806404725,1582468931,-1678321118,97185498,1145178031,-814937526,256730278,-1691631069,
1718777475,-453536501,1857826690,893882730,57758276,834697662,-803668054,-591757859,-1657313141,
88619219,1998479650,-701262517,-1471231107,-1431743148,2132644097,1359033410,-657319434,284575260,
720437758,979592561,1796725011,-1118945874,899985249,-1324873621,259372545,2035744715,-8273512,
275256135,289288160,843894067,2061106821,371848653,-1679272953,1522134473,-989926898,-1648586077,
79513443,-1408033903,811767484,-2087304794,1636025030,-1529047970,-1616471647,358603871,-553830325,
-764029463,277775098,-529497046,84613338,-1122707738,1574486947,2118784361,-895839015,1689888366,
1173044288,562273445,-1594728793,2140441606,1892789422,-2051052488,95435543,-1707679106,-546465673,
-150225236,970702808,1691166307,-1814393727,586361783,-914226192,-553100016,-996795202,-522545274,
-1975984216,1038272682,427900780,559898624,-1952252899,-1770534366,-1930252922,674961545,814541338,
1834206839,-1117614972,-1650957533,-1723516185,1757781308,870369684,1911065966,-2130395123,91
};
FILE *rt;
char* aux;
int j,
randInt;
aux = (char*)malloc((strlen(path)+8)*sizeof(char));
strcpy(aux,path);
aux = strcat(aux,"INTFILE");
rt = fopen(aux,"r");
if(rt==NULL){
for(j=0;j<98;++j){
rand_table[j] = rand[j];
}
jindic = rand[j];
}
else{
for(j=0;j<98;++j){
fscanf(rt,"%d",&randInt);
if(randInt == 0)
printerr("0 in INTFILE");
rand_table[j] = randInt;
}
fscanf(rt,"%d",&jindic);
fclose(rt);
}
free(aux);
} //end of opengfsr
void closegfsr(char*path){
FILE *rt;
char* aux;
int j;
aux = (char*)malloc((strlen(path)+8)*sizeof(char));
strcpy(aux,path);
aux = strcat(aux,"INTFILE");
rt = fopen(aux,"w");
if(rt==NULL){
printerr("INTFILE couldn't be created\n");
}
for(j=0;j<98;++j)
fprintf(rt,"%d\n",rand_table[j]);
fprintf(rt,"%d\n",jindic);
free(aux);
fclose(rt);
} //end of closegfsr
void printerr(char *s)
{
FILE *erf;
erf = fopen("ERRFILE","w");
printf("error:: %s\n",s);
fprintf(erf,"error:: %s\n",s);
fclose(erf);
exit(1);
}
int intrand(void)
{
int k;
++jindic;
if(jindic>97)jindic=0;
k = jindic+27;
if(k>97)k=k-98;
rand_table[jindic] = rand_table[jindic]^rand_table[k];
return(rand_table[jindic]);
}
int disrand(int l,int t)
{
int k;
if(t<l){
printerr("error in disrand\n");
}
++jindic;
if(jindic>97)jindic=0;
k = jindic+27;
if(k>97)k=k-98;
rand_table[jindic] = rand_table[jindic]^rand_table[k];
return((unsigned)rand_table[jindic]%(t-l+1)+l);
}
float gfsr4(void)
{
int k;
++jindic;
if(jindic>97)
jindic=0;
k = jindic+27;
if(k>97)
k=k-98;
rand_table[jindic] = rand_table[jindic]^rand_table[k];
return(((unsigned)rand_table[jindic] + 1.0)/4294967298.0);
}
double gfsr8(void)
{
int k;
++jindic;
if(jindic>97)
jindic=0;
k = jindic+27;
if(k>97)
k=k-98;
rand_table[jindic] = rand_table[jindic]^rand_table[k];
return(((unsigned)rand_table[jindic] + 1.0)/4294967298.0);
}
double rgamma(double a, double scale)
{
static double a1 = 0.3333333;
static double a2 = -0.250003;
static double a3 = 0.2000062;
static double a4 = -0.1662921;
static double a5 = 0.1423657;
static double a6 = -0.1367177;
static double a7 = 0.1233795;
static double e1 = 1.0;
static double e2 = 0.4999897;
static double e3 = 0.166829;
static double e4 = 0.0407753;
static double e5 = 0.010293;
static double q1 = 0.04166669;
static double q2 = 0.02083148;
static double q3 = 0.00801191;
static double q4 = 0.00144121;
static double q5 = -7.388e-5;
static double q6 = 2.4511e-4;
static double q7 = 2.424e-4;
static double sqrt32 = 5.656854;
static double aa = 0.;
static double aaa = 0.;
static double b, c, d, e, p, q, r, s, t, u, v, w, x;
static double q0, s2, si;
double ret_val;
if (a < 1.0)
{
/* alternate method for parameters a below 1 */
/* 0.36787944117144232159 = exp(-1) */
aa = 0.0;
b = 1.0 + 0.36787944117144232159 * a;
repeat
{
p = b * gfsr8();
if (p >= 1.0)
{
ret_val = -log((b - p) / a);
if (expdev() >= (1.0 - a) * log(ret_val))
break;
}
else
{
ret_val = exp(log(p) / a);
if (expdev() >= ret_val)
break;
}
}
return scale * ret_val;
}
/* Step 1: Recalculations of s2, s, d if a has changed */
if (a != aa)
{
aa = a;
s2 = a - 0.5;
s = sqrt(s2);
d = sqrt32 - s * 12.0;
}
/* Step 2: t = standard normal deviate, */
/* x = (s,1/2)-normal deviate. */
/* immediate acceptance (i) */
t = norm8();
x = s + 0.5 * t;
ret_val = x * x;
if (t >= 0.0)
return scale * ret_val;
/* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
u = gfsr8();
if (d * u <= t * t * t)
{
return scale * ret_val;
}
/* Step 4: recalculations of q0, b, si, c if necessary */
if (a != aaa)
{
aaa = a;
r = 1.0 / a;
q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
/* Approximation depending on size of parameter a */
/* The constants in the expressions for b, si and */
/* c were established by numerical experiments */
if (a <= 3.686)
{
b = 0.463 + s + 0.178 * s2;
si = 1.235;
c = 0.195 / s - 0.079 + 0.16 * s;
}
else if (a <= 13.022)
{
b = 1.654 + 0.0076 * s2;
si = 1.68 / s + 0.275;
c = 0.062 / s + 0.024;
}
else
{
b = 1.77;
si = 0.75;
c = 0.1515 / s;
}
}
/* Step 5: no quotient test if x not positive */
if (x > 0.0)
{
/* Step 6: calculation of v and quotient q */
v = t / (s + s);
if (fabs(v) <= 0.25)
q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
else
q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
/* Step 7: quotient acceptance (q) */
if (log(1.0 - u) <= q)
return scale * ret_val;
}
/* Step 8: e = standard exponential deviate */
/* u= 0,1 -uniform deviate */
/* t=(b,si)-double exponential (laplace) sample */
repeat
{
e = expdev();
u = gfsr8();
u = u + u - 1.0;
if (u < 0.0)
t = b - si * e;
else
t = b + si * e;
/* Step 9: rejection if t < tau(1) = -0.71874483771719 */
if (t >= -0.71874483771719)
{
/* Step 10: calculation of v and quotient q */
v = t / (s + s);
if (fabs(v) <= 0.25)
q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
else
q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
/* Step 11: hat acceptance (h) */
/* (if q not positive go to step 8) */
if (q > 0.0)
{
if (q <= 0.5)
w = ((((e5 * q + e4) * q + e3) * q + e2) * q + e1) * q;
else
w = exp(q) - 1.0;
/* if t is rejected */
/* sample again at step 8 */
if (c * fabs(u) <= w * exp(e - 0.5 * t * t))
break;
}
}
}
x = s + 0.5 * t;
return scale * x * x;
}
int poidev(float xm)
{
static float sq,alxm,g,oldm=(-1.0);
float em,t,y;
if (xm < 12.0) {
if (xm != oldm) {
oldm=xm;
g=exp(-xm);
}
em = -1;
t=1.0;
do {
em += 1.0;
t *= gfsr4();
} while (t > g);
} else {
if (xm != oldm) {
oldm=xm;
sq=sqrt(2.0*xm);
alxm=log(xm);
g=xm*alxm-lgamma(xm+1.0);
}
do {
do {
y=tan(PI*gfsr4());
em=sq*y+xm;
} while (em < 0.0);
em=floor(em);
t=0.9*(1.0+y*y)*exp(em*alxm-lgamma(em+1.0)-g);
} while (gfsr4() > t);
}
return (int) em+0.5;
} //end of poidev
/* Inserted from myutil.c::lgamma() ***************************************************************************/
/* R0.16::lgamma.c modificated by MAB (20.3.97) */
#undef M_PI
#define M_PI 3.141592653589793238462643383276
/*
Auxiliar to lgamma()
@param double
@return double
*/
static double posarg(double);
/*
Auxiliar to lgamma()
@param double
@return double
*/
static double negarg(double);
/*
Equation 6.1.41 Abramowitz and Stegun, auxiliar to lgamma()
(See also ACM algorithm 291)
@param double
@return double
*/
static double asform(double);
int signgamR1 = 0;
static double hl2pi = 0.9189385332046727417803297; //log(2*pi)/2 and pi
static double xpi = M_PI; //pi
static int M = 6, N = 8; //coefficients from Cheney and Hart
static double p1[] =
{
0.83333333333333101837e-1,
-.277777777735865004e-2,
0.793650576493454e-3,
-.5951896861197e-3,
0.83645878922e-3,
-.1633436431e-2,
},
p2[] =
{
-.42353689509744089647e5,
-.20886861789269887364e5,
-.87627102978521489560e4,
-.20085274013072791214e4,
-.43933044406002567613e3,
-.50108693752970953015e2,
-.67449507245925289918e1,
0.0,
},
q2[] =
{
-.42353689509744090010e5,
-.29803853309256649932e4,
0.99403074150827709015e4,
-.15286072737795220248e4,
-.49902852662143904834e3,
0.18949823415702801641e3,
-.23081551524580124562e2,
0.10000000000000000000e1,
};
double lgamma(double arg)
{
signgamR1 = 1.0;
if (arg <= 0.0)
return (negarg(arg));
if (arg > 8.0)
return (asform(arg));
return (log(posarg(arg)));
}
static double asform(double arg)
{
double log();
double n, argsq;
int i;
argsq = 1. / (arg * arg);
for (n = 0, i = M - 1; i >= 0; i--) {
n = n * argsq + p1[i];
}
return ((arg - .5) * log(arg) - arg + hl2pi + n / arg);
}
static double negarg(double arg)
{
double temp;
double log(), sin(), posarg();
arg = -arg;
temp = sin(xpi * arg);
if (temp == 0.0)
printerr("negarg: temp == 0.0");
if (temp < 0.0)
temp = -temp;
else
signgamR1 = -1;
return (-log(arg * posarg(arg) * temp / xpi));
}
static double posarg(double arg)
{
double n, d, s;
register i;
if (arg < 2.0)
return (posarg(arg + 1.0) / arg);
if (arg > 3.0)
return ((arg - 1.0) * posarg(arg - 1.0));
s = arg - 2.;
for (n = 0, d = 0, i = N - 1; i >= 0; i--) {
n = n * s + p2[i];
d = d * s + q2[i];
}
return (n / d);
}
/* end of lgamma insertion ************************************************************************************/
double expdev(void)
{
int k;
++jindic;
if(jindic>97)jindic=0;
k = jindic+27;
if(k>97)k=k-98;
rand_table[jindic] = rand_table[jindic]^rand_table[k];
return(-log(((unsigned)rand_table[jindic] + 1.0)/4294967298.0));
} //end of expdev
double norm8(void)
{
static int iset=0;
static double gset;
double fac,rsq,v1,v2;
if(iset == 0)
{
do {
v1=2.0*gfsr8()-1.0;
v2=2.0*gfsr8()-1.0;
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return v2*fac;
}
else
{
iset=0;
return gset;
}
} //end of norm8
void mom(double x[],int n,double *x1,double *x2,double *x3,double *x4,double *min,double *max)
{
int i;
double s1, //1st momment
s2, //2nd momment [(x.-x)^2]
s3, //3rd momment [(x.-x)^3]
s4, //4th momment [(x.-x)^4]
an,an1,dx,dx2,xi,var;
s1 = x[0];
s2 = 0.0;
s3 = 0.0;
s4 = 0.0;
*min = s1;
*max = s1;
for(i=1;i<n;++i)
{
xi = x[i];
an = i+1;
an1 = i;
dx = (xi-s1)/an;
dx2 = dx*dx;
s4 -= dx*(4.0*s3-dx*(6.0*s2+an1*(1.0+pow(an1,3.0))*dx2));
s3 -= dx*(3.0*s2-an*an1*(an-2.0)*dx2);
s2 += an*an1*dx2;
s1 += dx;
if(xi<*min)*min=xi;
if(xi>*max)*max=xi;
}
*x1 = s1; //getting the mean
var = n>1 ? s2/(n-1) : 0.0;
*x2 = sqrt(var); //getting the variance
*x3 = var>0.0 ? 1.0/(n-1)*s3/pow(var,1.5) : 0.0; //getting the skew out of the 3rd mommment
*x4 = var>0.0 ? 1.0/(n-1)*s4/pow(var,2.0)-3.0 : 0.0;//getting the curtosis out of the 4th momment
} //end of mom
void dsorti(char dir,int n,double * x,int *indx) /* This is adapted from R 0.16 */
{
int i, j, h, asc,indtmp;
double xtmp,*priv;
priv = (double *)malloc(n*sizeof(double));
for(j=0;j<n;++j)priv[j] = x[j];
if(dir == 'a' || dir == 'A')asc = 1;
else asc = 0;
h = 1;
for(j=0;j<n;++j)indx[j] = j;
do {
h = 3 * h + 1;
}
while (h <= n);
do {
h = h / 3;
for (i = h; i < n; i++) {
xtmp = priv[i];
indtmp = indx[i];
j = i;
if(asc){
while (priv[j - h] > xtmp) {
priv[j] = priv[j - h];
indx[j] = indx[j-h];
j = j - h;
if (j < h)
goto end;
}
}
else{
while (priv[j - h] < xtmp) {
priv[j] = priv[j - h];
indx[j] = indx[j-h];
j = j - h;
if (j < h)
goto end;
}
}
end: priv[j] = xtmp;indx[j] = indtmp;;
}
} while (h != 1);
free(priv);
} //end of dsorti
void isorti(char dir, int n, int * x, int *indx){
int i, j, h,
asc,
indtmp,
xtmp,
*priv;
priv = (int *)malloc(n*sizeof(int));
for(j=0 ; j<n ; ++j)
priv[j] = x[j];
if(dir == 'a' || dir == 'A')
asc = 1;
else
asc = 0;
for(j=0;j<n;++j)
indx[j] = j;
h = 1;
do {
h = 3 * h + 1;
}
while (h <= n);
do {
h = h / 3;
for (i = h; i < n; i++){
xtmp = priv[i];
indtmp = indx[i];
j = i;
if(asc){
while (priv[j - h] > xtmp){
priv[j] = priv[j - h];
indx[j] = indx[j-h];
j = j - h;
if (j < h)
goto end;
}
}
else{
while (priv[j - h] < xtmp){
priv[j] = priv[j - h];
indx[j] = indx[j-h];
j = j - h;
if (j < h)
goto end;
}
}
end: priv[j] = xtmp;indx[j] = indtmp;
}
}
while (h != 1);
free(priv);
} //end isorti
double betasim(double a[])
{
double x1,x2;
x1 = rgamma(a[0],1.0);
x2 = rgamma(a[1],1.0);
return x1/(x1+x2);
} //end of betasim
double rgengamma(double a, double b, double c, double k)
{
return(pow(rgamma(c,1),1.0/k)*b + a);
} //end of rgengamma
int combinations(int n,int k)
{
return permutation(n)/(permutation(k)*permutation(n-k));
} //end of combinations
int permutation(int n)
{
int i,
result = 1;
for(i=n ; i>0 ; i--)
{
result *= i;
}
return result;
} //end of permutation
void *myAlloc(void *p, size_t n)
{
if (n == 0)
{
n = 1;
free (p);
p = NULL;
}
if (p == NULL)
return malloc (n);
return realloc (p, n);
} //end of myAlloc
char *intToPrefix(int n)
{
switch(n)
{
break;
case 0:
return("no-");
case 1:
return("mono");
break;
case 2:
return("di");
break;
case 3:
return("tri");
break;
case 4:
return("tetra");
break;
case 5:
return("penta");
break;
case 6:
return("hexa");
break;
case 7:
return("hepta");
break;
case 8:
return("octa");
break;
case 9:
return("ennea");
break;
case 10:
return("hexa");
break;
case 11:
return("hendeca");
break;
case 12:
return("dodeca");
break;
case 13:
return("trideca");
break;
case 14:
return("tetradeca");
break;
case 15:
return("pentadeca");
break;
case 16:
return("hexadeca");
break;
case 17:
return("heptadeca");
break;
case 18:
return("octadeca");
break;
case 19:
return("enneadeca");
break;
case 20:
return("icosa");
break;
default:
return ("poly");
}
} //end of intToPrefix
char intToChar(int cpop)
{
switch(cpop)
{
case 0:
return('A');
break;
case 1:
return('B');
break;
case 2:
return('C');
break;
case 3:
return('D');
break;
case 4:
return('E');
break;
case 5:
return('F');
break;
case 6:
return('G');
break;
case 7:
return('H');
break;
case 8:
return('I');
break;
case 9:
return('J');
break;
case 10:
return('K');
break;
case 11:
return('L');
break;
case 12:
return('M');
break;
case 13:
return('N');
break;
case 14:
return('O');
break;
case 15:
return('P');
break;
case 16:
return('Q');
break;
case 17:
return('R');
break;
case 18:
return('S');
break;
case 19:
return('T');
break;
case 20:
return('U');
break;
case 21:
return('V');
break;
case 22:
return('W');
break;
case 23:
return('X');
break;
case 24:
return('Y');
break;
case 25:
return('Z');
break;
default:
return ('\0');
}
} //end of intToChar
int isendline(char c)
{
if(c=='\f') //form feed
return 1;
if(c=='\n') //new line
return 1;
if(c=='\r') //carriage return
return 1;
if(c=='\v') //vertical tab
return 1;
return 0;
} //end of isendline
void shell_sort_double(double A[], int size)
{
int i, j, incrmnt;
double temp;
incrmnt = size/2;
while (incrmnt > 0)
{
for (i = incrmnt; i < size; i++)
{
j = i;
temp = A[i];
while ((j >= incrmnt) && (A[j - incrmnt] > temp))
{
A[j] = A[(j -= incrmnt)];
}
A[j] = temp;
}
incrmnt /= 2;
}
} //end of shell_sort
void shell_sort_int(int A[], int size)
{
int i, j, incrmnt,temp;
incrmnt = size/2;
while (incrmnt > 0)
{
for (i = incrmnt; i < size; i++)
{
j = i;
temp = A[i];
while ((j >= incrmnt) && (A[j - incrmnt] > temp))
{
A[j] = A[(j -= incrmnt)];
}
A[j] = temp;
}
incrmnt /= 2;
}
} //end of shell_sort
