diff --git a/CUtils/CUtils.xs b/CUtils/CUtils.xs
index 9c977ddfe7594d824f9806f36ed15e30bcfbdfb1..5c4f2da61b0628a1956196642f06eeb288dba519 100644
--- a/CUtils/CUtils.xs
+++ b/CUtils/CUtils.xs
@@ -128,7 +128,7 @@ DoublePermutation(nb_sample, nb_chi2, data)
############################################################
HV *
-ResamplingChi2(leaf_refs, leaf_depth, leaf_parent, nleaf_parent, max_depth, prolonge, nb_permutations)
+ResamplingChi2(leaf_refs, leaf_depth, leaf_parent, nleaf_parent, max_depth, prolonge, nb_permutations, parallel)
AV * leaf_refs
AV * leaf_depth
AV * leaf_parent
@@ -136,6 +136,7 @@ ResamplingChi2(leaf_refs, leaf_depth, leaf_parent, nleaf_parent, max_depth, prol
int max_depth
int prolonge
int nb_permutations
+ SV * parallel
INIT:
struct cc *lcc;
struct tree tree;
@@ -144,6 +145,7 @@ ResamplingChi2(leaf_refs, leaf_depth, leaf_parent, nleaf_parent, max_depth, prol
int i;
AV * ra;
+ int cparallel;
CODE:
tree.nb_leaves=av_len(leaf_refs);
tree.nb_nodes=av_len(nleaf_parent);
@@ -185,9 +187,14 @@ ResamplingChi2(leaf_refs, leaf_depth, leaf_parent, nleaf_parent, max_depth, prol
for (i=0; i<tree.nb_nodes; i++) {
tree.np[i]=SvNV(*av_fetch(nleaf_parent, i, 0));
}
+ if (!SvOK(parallel) || !SvIOK(parallel)) {
+ cparallel=0;
+ } else {
+ cparallel=SvIV(parallel);
+ }
int res=resampling_chi2(&tree, lcc, prolonge, nb_permutations,
- results);
+ results, cparallel);
RETVAL = newHV();
sv_2mortal((SV *)RETVAL);
diff --git a/CUtils/Makefile.PL b/CUtils/Makefile.PL
index d630637b030252d5bd3ae91be3e9b2f9aa1af200..8d5da2d45f460d3f2d3f1e7676756e72c67586bc 100644
--- a/CUtils/Makefile.PL
+++ b/CUtils/Makefile.PL
@@ -9,7 +9,7 @@ WriteMakefile(
($] >= 5.005 ? ## Add these new keywords supported since 5.005
(ABSTRACT_FROM => 'lib/ALTree/CUtils.pm', # retrieve abstract from module
AUTHOR => 'Claire Bardel <Claire.Bardel@univ-lyon1.fr>') : ()),
- LIBS => ['-lm '], # e.g., '-lm'
+ LIBS => ['-lm ', '-lpthread'], # e.g., '-lm'
DEFINE => '', # e.g., '-DHAVE_SOMETHING'
INC => '-I. -Ic_sources', # e.g., '-I. -I/usr/include/other'
# Un-comment this if you add C files to link with later:
diff --git a/CUtils/c_sources/mt19937ar.c b/CUtils/c_sources/mt19937ar.c
new file mode 100644
index 0000000000000000000000000000000000000000..3c99f71e619566f7a02e8f07cdc8276b44d19ee5
--- /dev/null
+++ b/CUtils/c_sources/mt19937ar.c
@@ -0,0 +1,214 @@
+/*
+ A C-program for MT19937, with initialization improved 2002/1/26.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+
+ Before using, initialize the state by using init_genrand(seed)
+ or init_by_array(init_key, key_length).
+
+ Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ 3. The names of its contributors may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 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.
+
+
+ Any feedback is very welcome.
+ http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+ email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+*/
+
+#include <stdio.h>
+#include "mt19937ar.h"
+
+/* Period parameters */
+#define N MT19937AR_N /* 624 */
+#define M 397
+#define MATRIX_A 0x9908b0dfUL /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+static mt19937ar_t local = { {0, }, /* the array for the state vector */
+ N+1, /* mti==N+1 means mt[N] is not initialized */
+};
+
+/* initializes mt[N] with a seed */
+void init_genrand_mt(mt19937ar_t* p, unsigned long s)
+{
+ int mti;
+ unsigned long *mt=&p->mt[0];
+ mt[0]= s & 0xffffffffUL;
+ for (mti=1; mti<N; mti++) {
+ mt[mti] =
+ (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array mt[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ mt[mti] &= 0xffffffffUL;
+ /* for >32 bit machines */
+ }
+ p->mti=mti;
+}
+
+void init_genrand(unsigned long s)
+{
+ init_genrand_mt(&local, s);
+}
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+/* slight change for C++, 2004/2/26 */
+void init_by_array_mt(mt19937ar_t* p, unsigned long init_key[], int key_length)
+{
+ int i, j, k;
+ unsigned long *mt=&p->mt[0];
+ init_genrand_mt(p, 19650218UL);
+ i=1; j=0;
+ k = (N>key_length ? N : key_length);
+ for (; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ + init_key[j] + j; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++; j++;
+ if (i>=N) { mt[0] = mt[N-1]; i=1; }
+ if (j>=key_length) j=0;
+ }
+ for (k=N-1; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+ - i; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ if (i>=N) { mt[0] = mt[N-1]; i=1; }
+ }
+
+ mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
+}
+void init_by_array(unsigned long init_key[], int key_length)
+{
+ init_by_array_mt(&local, init_key, key_length);
+}
+
+/* generates a random number on [0,0xffffffff]-interval */
+unsigned long genrand_int32_mt(mt19937ar_t* p)
+{
+ unsigned long *mt=&p->mt[0];
+ unsigned long y;
+ static const unsigned long mag01[2]={0x0UL, MATRIX_A};
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ if (p->mti >= N) { /* generate N words at one time */
+ int kk;
+
+ for (kk=0;kk<N-M;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ for (;kk<N-1;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+ mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+ p->mti = 0;
+ }
+
+ y = mt[p->mti++];
+
+ /* Tempering */
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ y ^= (y >> 18);
+
+ return y;
+}
+
+unsigned long genrand_int32(void)
+{
+ if (local.mti == N+1) /* if init_genrand() has not been called, */
+ init_genrand(5489UL); /* a default initial seed is used */
+
+ return genrand_int32_mt(&local);
+}
+
+/* generates a random number on [0,0x7fffffff]-interval */
+long genrand_int31(void)
+{
+ return (long)(genrand_int32()>>1);
+}
+
+/* generates a random number on [0,1]-real-interval */
+double genrand_real1(void)
+{
+ return genrand_int32()*(1.0/4294967295.0);
+ /* divided by 2^32-1 */
+}
+
+/* generates a random number on [0,1)-real-interval */
+double genrand_real2(void)
+{
+ return genrand_int32()*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/* generates a random number on (0,1)-real-interval */
+double genrand_real3(void)
+{
+ return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/* generates a random number on [0,1) with 53-bit resolution*/
+double genrand_res53(void)
+{
+ unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
+ return(a*67108864.0+b)*(1.0/9007199254740992.0);
+}
+/* These real versions are due to Isaku Wada, 2002/01/09 added */
+
+#if 0
+int main(void)
+{
+ int i;
+ unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
+ init_by_array(init, length);
+ printf("1000 outputs of genrand_int32()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10lu ", genrand_int32());
+ if (i%5==4) printf("\n");
+ }
+ printf("\n1000 outputs of genrand_real2()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10.8f ", genrand_real2());
+ if (i%5==4) printf("\n");
+ }
+ return 0;
+}
+#endif
+
diff --git a/CUtils/c_sources/mt19937ar.h b/CUtils/c_sources/mt19937ar.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ec8623f5f83bd70a53939449537cf11886bd2a9
--- /dev/null
+++ b/CUtils/c_sources/mt19937ar.h
@@ -0,0 +1,86 @@
+/*
+ A C-program for MT19937, with initialization improved 2002/1/26.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+
+ Before using, initialize the state by using init_genrand(seed)
+ or init_by_array(init_key, key_length).
+
+ Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ 3. The names of its contributors may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 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.
+
+
+ Any feedback is very welcome.
+ http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+ email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+#define MT19937AR_N 624
+ unsigned long mt[MT19937AR_N]; /* the array for the state vector */
+ int mti; /* mti==N+1 means mt[N] is not initialized */
+} mt19937ar_t;
+
+/* initializes mt[N] with a seed */
+void init_genrand(unsigned long s);
+void init_genrand_mt(mt19937ar_t* priv, unsigned long s);
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+/* slight change for C++, 2004/2/26 */
+void init_by_array(unsigned long init_key[], int key_length);
+void init_by_array_mt(mt19937ar_t* priv, unsigned long init_key[], int key_length);
+
+/* generates a random number on [0,0xffffffff]-interval */
+unsigned long genrand_int32(void);
+unsigned long genrand_int32_mt(mt19937ar_t* priv);
+
+/* generates a random number on [0,0x7fffffff]-interval */
+long genrand_int31(void);
+
+/* generates a random number on [0,1]-real-interval */
+double genrand_real1(void);
+
+/* generates a random number on [0,1)-real-interval */
+double genrand_real2(void);
+
+/* generates a random number on (0,1)-real-interval */
+double genrand_real3(void);
+
+/* generates a random number on [0,1) with 53-bit resolution*/
+double genrand_res53(void);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/CUtils/c_sources/myrand.c b/CUtils/c_sources/myrand.c
index 812e2dc1427d1d6e973bda6591ea92be951906e7..4590e9af829446a4f91571a904082e20dea6c68e 100644
--- a/CUtils/c_sources/myrand.c
+++ b/CUtils/c_sources/myrand.c
@@ -5,22 +5,48 @@
#include "debug.h"
-int myrand(int up) {
+#define MT
- static int init=0;
- if (!init) {
- srand(getpid());
- init=1;
- //debug("init rand done");
- }
+#ifdef MT
+# include "mt19937ar.h"
+
+static __thread mt19937ar_t rdata;
+# define MYRAND_MAX (0xffffffffUL)
+#else
+# define MYRAND_MAX RAND_MAX
+#endif
+
+
+int myrand_init(unsigned long value) {
+#ifdef MT
+ init_genrand_mt(&rdata, value);
+#endif
+}
+
+static void init() __attribute__((constructor));
+static void init() {
+#ifdef MT
+ myrand_init(getpid());
+#else
+ srand(getpid());
+#endif
+ debug("init rand done");
+}
+
+
+int myrand(int up) {
for(;;) {
+#ifdef MT
+ unsigned long r=genrand_int32_mt(&rdata);
+#else
int r=rand();
+#endif
- if (r<=RAND_MAX-up) {
+ if (r<=MYRAND_MAX-up) {
return r%up;
}
- if (r<RAND_MAX-(RAND_MAX%up)) {
+ if (r<MYRAND_MAX-(MYRAND_MAX%up)) {
return r%up;
}
}
diff --git a/CUtils/c_sources/myrand.h b/CUtils/c_sources/myrand.h
index bfad08d40333e757224cb1d88c45e880919539f2..fc859e34a0e8b2ebc34839e6152354485359a564 100644
--- a/CUtils/c_sources/myrand.h
+++ b/CUtils/c_sources/myrand.h
@@ -1,6 +1,7 @@
#ifndef _MYRAND_H
#define _MYRAND_H
+int myrand_init(unsigned long value);
int myrand(int);
#endif
diff --git a/CUtils/c_sources/resampling.c b/CUtils/c_sources/resampling.c
index fbe1b0f4c1041000fdc3c5e1e20828b5b80f17cc..d9d5b8861471bcbe75a0db2921055e065d8fd2ad 100644
--- a/CUtils/c_sources/resampling.c
+++ b/CUtils/c_sources/resampling.c
@@ -2,7 +2,9 @@
#include "debug.h"
#include "resampling.h"
#include <strings.h>
-
+#include <stdlib.h>
+#include <pthread.h>
+#include "myrand.h"
struct chi2s {
int nb_leaves;
@@ -83,28 +85,103 @@ static void compute_chi2s(const struct tree *tree, const struct cc *lcc,
}
}
-int resampling_chi2(struct tree *tree, struct cc *lcc, int prolonge,
- int nb_permutations, datatype_t *results)
-{
+struct paractl {
+ volatile int permutation;
+ int nb_permutations;
+ const struct tree *tree;
+ const struct cc *lcc;
+ int cases;
+ int controls;
+ int prolonge;
+ datatype_t *results;
+};
- struct cc rand_lcc[tree->nb_leaves];
+struct parainfo {
+ struct paractl *ctl;
+ int permutation;
+};
- struct cc leaves[tree->nb_leaves];
- struct cc th[tree->nb_leaves];
- struct cc nodes[tree->nb_nodes];
- int stashed_nodes[tree->nb_nodes];
- int next_nodes[tree->nb_nodes];
- struct chi2s temp={
+struct memory {
+ struct cc *rand_lcc;
+ struct chi2s temp;
+};
+
+static struct memory *mem_alloc(const struct tree *tree)
+{
+ struct memory *m=malloc(sizeof(struct memory));
+
+ m->rand_lcc=malloc(tree->nb_leaves*sizeof(struct cc));
+ m->temp=(struct chi2s){
.nb_leaves=tree->nb_leaves,
- .leaves=leaves,
- .th=th,
+ .leaves=malloc(tree->nb_leaves*sizeof(struct cc)),
+ .th=malloc(tree->nb_leaves*sizeof(struct cc)),
.nb_nodes=tree->nb_nodes,
- .stashed_nodes=stashed_nodes,
- .next_nodes=next_nodes,
- .nodes=nodes,
+ .stashed_nodes=malloc(tree->nb_nodes*sizeof(int)),
+ .next_nodes=malloc(tree->nb_nodes*sizeof(int)),
+ .nodes=malloc(tree->nb_nodes*sizeof(struct cc)),
};
-
- compute_chi2s(tree, lcc, &temp, prolonge, results);
+
+ return m;
+}
+
+static void free_alloc(struct memory *m)
+{
+ free(m->temp.leaves);
+ free(m->temp.th);
+ free(m->temp.stashed_nodes);
+ free(m->temp.next_nodes);
+ free(m->temp.nodes);
+ free(m->rand_lcc);
+}
+
+static void *resampling_worker(void* param) {
+ struct parainfo *p=(struct parainfo *)param;
+
+ /* Copy values for ro access, to avoid ping-pong with cache lines
+ when updating permutation */
+ struct paractl data= *(p->ctl);
+ int permutation=p->permutation;
+
+ struct memory *mem=mem_alloc(data.tree);
+
+ myrand_init(pthread_self()+getpid());
+
+ do {
+ //debug("thread %i handle perm %i", p->permutation, permutation);
+ random_clades(data.tree->nb_leaves, data.lcc,
+ data.cases, data.controls, mem->rand_lcc);
+ compute_chi2s(data.tree, mem->rand_lcc, &mem->temp, data.prolonge,
+ data.results+(permutation * data.tree->max_depth));
+ permutation=__sync_fetch_and_add(&p->ctl->permutation, 1);
+ } while (permutation < data.nb_permutations);
+
+ free_alloc(mem);
+
+ return NULL;
+}
+
+int resampling_chi2(const struct tree *tree, const struct cc *lcc, int prolonge,
+ int nb_permutations, datatype_t *results, int parallel)
+{
+
+ FILE* dump=fopen("/tmp/dump", "w+");
+ fwrite(tree, sizeof(struct tree), 1, dump);
+ fwrite(tree->ld, sizeof(int), tree->nb_leaves, dump);
+ fwrite(tree->lp, sizeof(int), tree->nb_leaves, dump);
+ fwrite(tree->np, sizeof(int), tree->nb_nodes, dump);
+ fwrite(lcc, sizeof(struct cc), tree->nb_leaves, dump);
+ fwrite(&prolonge, sizeof(int), 1, dump);
+ fwrite(&nb_permutations, sizeof(int), 1, dump);
+ fclose(dump);
+
+ char* envvar=getenv("ALTREE_PARALLEL");
+ if (envvar) {
+ parallel=atoi(envvar);
+ }
+ debug("parallel=%i", parallel);
+
+ struct memory *mem=mem_alloc(tree);
+ compute_chi2s(tree, lcc, &mem->temp, prolonge, results);
int i;
int cases=0;
int controls=0;
@@ -112,11 +189,105 @@ int resampling_chi2(struct tree *tree, struct cc *lcc, int prolonge,
cases += lcc[i].cases;
controls += lcc[i].controls;
}
- for(i=0; i<nb_permutations; i++) {
- results += tree->max_depth;
- random_clades(tree->nb_leaves, lcc, cases, controls, rand_lcc);
- compute_chi2s(tree, rand_lcc, &temp, prolonge, results);
+ if (!parallel) {
+ for(i=0; i<nb_permutations; i++) {
+ results += tree->max_depth;
+ random_clades(tree->nb_leaves, lcc, cases, controls, mem->rand_lcc);
+ compute_chi2s(tree, mem->rand_lcc, &mem->temp, prolonge, results);
+ }
+ } else {
+ struct paractl ctl = {
+ .permutation=parallel,
+ .nb_permutations=nb_permutations,
+ .tree=tree,
+ .lcc=lcc,
+ .cases=cases,
+ .controls=controls,
+ .prolonge=prolonge,
+ .results=results+tree->max_depth,
+ };
+
+ struct parainfo infos[parallel];
+ pthread_t tids[parallel];
+
+ for(i=0; i<parallel; i++) {
+ infos[i].ctl=&ctl;
+ infos[i].permutation=i;
+ pthread_create(&tids[i], NULL, &resampling_worker, &infos[i]);
+ }
+ for(i=0; i<parallel; i++) {
+ pthread_join(tids[i], NULL);
+ }
+ }
+ free_alloc(mem);
+
+ return 0;
+}
+
+#ifdef MAIN_PROG
+int main(int argc, char *argv[])
+{
+ int resampling_chi2(const struct tree *tree, const struct cc *lcc, int prolonge,
+ int nb_permutations, datatype_t *results, int parallel);
+ struct tree tree;
+
+ FILE* dump=fopen("/tmp/dump.read", "r");
+ fread(&tree, sizeof(struct tree), 1, dump);
+
+ int ld[tree.nb_leaves];
+ int lp[tree.nb_leaves];
+ int np[tree.nb_nodes];
+ struct cc lcc[tree.nb_leaves];
+ int prolonge;
+ int nb_permutations;
+
+ tree.ld=ld;
+ tree.lp=lp;
+ tree.np=np;
+
+ fread(tree.ld, sizeof(int), tree.nb_leaves, dump);
+ fread(tree.lp, sizeof(int), tree.nb_leaves, dump);
+ fread(tree.np, sizeof(int), tree.nb_nodes, dump);
+ fread(lcc, sizeof(struct cc), tree.nb_leaves, dump);
+ fread(&prolonge, sizeof(int), 1, dump);
+ fread(&nb_permutations, sizeof(int), 1, dump);
+ nb_permutations=10;
+ fclose(dump);
+
+ datatype_t results[(nb_permutations+1)*tree.max_depth];
+ //bzero(results, (nb_permutations+1)*tree.max_depth*sizeof(datatype_t));
+
+ resampling_chi2(&tree, lcc, prolonge,
+ nb_permutations, results, 0);
+
+ int i,j;
+ datatype_t *r=results;
+ for(i=0; i<=nb_permutations; i++) {
+ for(j=0; j<tree.max_depth; j++) {
+ printf("\t"CONV, *(r++));
+ }
+ printf("\n");
}
+ /* ensemble_t ens_min_pval; */
+
+ /* ens_min_pval=alloc_ensemble(nb_sample); */
+
+ /* min=double_permutation(nb_sample, nb_chi2, mat, rep, ens_min_pval); */
+
+ /* free_ensemble(ens_min_pval); */
+
+ /* for (j=0; j<nb_chi2; j++) { */
+ /* printf("chi2 niveau %i, pval nc "CONV"\n", j+1, rep[j]); */
+ /* } */
+ /* for (i=0; i<nb_sample; i++) { */
+ /* printf("sample %i, pval min "CONV"\n", i, ens_min_pval[i]); */
+ /* } */
+ /* printf("pmin corrigé: "CONV"\n", min); */
+
+ /* free_matrice(mat, nb_sample, nb_chi2); */
+ /* free_replicat(rep); */
+
return 0;
}
+#endif
diff --git a/CUtils/c_sources/resampling.h b/CUtils/c_sources/resampling.h
index 998435b56d1a8747cab347734b9c6e9509eba66b..c84e02680d5d89296b7ad2996312865ed3aa6678 100644
--- a/CUtils/c_sources/resampling.h
+++ b/CUtils/c_sources/resampling.h
@@ -13,7 +13,8 @@ struct tree {
int *np;
};
-int resampling_chi2(struct tree *tree, struct cc *lcc, int prolonge,
- int nb_permutations, datatype_t *results);
+int resampling_chi2(const struct tree *tree, const struct cc *lcc,
+ int prolonge, int nb_permutations,
+ datatype_t *results, int parallel);
#endif
diff --git a/progs/altree b/progs/altree
index d02dfb8caed25f44b75d6de7adfac7df7907432e..834100e3513f4ce108821832f5a132e88d8b4c28 100755
--- a/progs/altree
+++ b/progs/altree
@@ -712,7 +712,7 @@ sub RepeatAssociation
}
my $resamp=ALTree::CUtils::ResamplingChi2(
$leaf_refs, $leaf_depth, $leaf_parent, $nleaf_parent,
- $max_depth, $prolonge, $nb_permutation
+ $max_depth, $prolonge, $nb_permutation, 0
);
$ligne_stats=$resamp->{"chi2s"};