diff --git a/CUtils/CUtils.xs b/CUtils/CUtils.xs
index 2cbfd78f7217e43ca6f6b9f75ccdf56264346310..fc4dba327d1451441c71572781a0590fcac66ab1 100644
--- a/CUtils/CUtils.xs
+++ b/CUtils/CUtils.xs
@@ -45,32 +45,30 @@ right(a, b, c, d)
double c
double d
-SV *
+HV *
double_permutation(nb_sample, nb_chi2, data)
int nb_sample
int nb_chi2
- SV * data
+ AV * data
INIT:
matrice_t mat;
replicat_t rep;
+ ensemble_t ens;
int i,j,n;
datatype_t min;
- HV * rh;
AV * ra;
CODE:
//fprintf(stderr, "\nSTART(%i, %i)\n", nb_sample, nb_chi2);
if ((nb_sample <= 0)
|| (nb_chi2 <= 0)
- || (!SvROK(data))
- || (SvTYPE(SvRV(data)) != SVt_PVAV)
- || (av_len((AV *)SvRV(data)) != nb_sample*nb_chi2-1))
+ || (av_len(data) != nb_sample*nb_chi2-1))
{
XSRETURN_UNDEF;
}
- rh = (HV *)sv_2mortal((SV *)newHV());
mat=alloc_matrice(nb_sample, nb_chi2);
rep=alloc_replicat(nb_chi2);
+ ens=alloc_ensemble(nb_sample);
n=0;
for (i=0; i<nb_sample; i++) {
@@ -78,26 +76,34 @@ double_permutation(nb_sample, nb_chi2, data)
//fprintf(stderr, "[%i][%i](%i)...\n", i, j, n);
/* Attention: on place un réplicat par colonne
* (et pas par ligne) */
- mat[j][i]=SvNV(*av_fetch((AV *)SvRV(data), n, 0));
+ mat[j][i]=SvNV(*av_fetch(data, n, 0));
//fprintf(stderr, "[%i][%i](%i)=%lf\n", i, j, n,mat[j][i]);
n++;
}
}
- min=calcul(nb_sample, nb_chi2, mat, rep);
+ min=calcul_distrib_pmin(nb_sample, nb_chi2, mat, rep, ens);
- hv_store(rh, "pmin", 4, newSVnv(min), 0);
+ RETVAL = newHV();
+ sv_2mortal((SV *)RETVAL);
+
+ hv_store(RETVAL, "pmin", 4, newSVnv(min), 0);
ra = (AV *)sv_2mortal((SV *)newAV());
for (j=0; j<nb_chi2; j++) {
av_push(ra, newSVnv(rep[j]));
//fprintf(stderr, "chi2[%i]=%lf\n", j ,rep[j]);
- }
-
- hv_store(rh, "chi2", 4, newRV((SV *)ra), 0);
+ }
+ hv_store(RETVAL, "chi2", 4, newRV((SV *)ra), 0);
+
+ ra = (AV *)sv_2mortal((SV *)newAV());
+ for (j=0; j<nb_sample; j++) {
+ av_push(ra, newSVnv(ens[j]));
+ //fprintf(stderr, "pmin[%i]=%lf\n", j ,rep[j]);
+ }
+ hv_store(RETVAL, "distrib_pmin", 12, newRV((SV *)ra), 0);
free_matrice(mat, nb_sample, nb_chi2);
free_replicat(rep);
-
- RETVAL = newRV((SV *)rh);
+ free_ensemble(ens);
OUTPUT:
RETVAL
diff --git a/CUtils/c_sources/double_permutation.c b/CUtils/c_sources/double_permutation.c
index e99875635aa0114369ba249fda49a362bc5a2c2a..66cec1b6228f86b0f264b63a7a698ff9560b7bc3 100644
--- a/CUtils/c_sources/double_permutation.c
+++ b/CUtils/c_sources/double_permutation.c
@@ -6,7 +6,7 @@
#include "double_permutation.h"
#define CALC_PVAL(count, nb_sample) \
- ((datatype_t)(count-1))/* On s'enlève soi-même*/ \
+ ((datatype_t)(count-1))/* On s'enlève soi-même*/ \
/nb_sample
@@ -19,10 +19,10 @@ int read_matrice(matrice_t mat, int nb_sample, int nb_chi2)
for (j=0; j<nb_chi2; j++) {
res=scanf(CONV, &d);
if (res!=1) {
- fprintf(stderr, "Erreur de lecture. Probablement pas assez de données\n");
+ fprintf(stderr, "Erreur de lecture. Probablement pas assez de données\n");
exit(1);
}
- /* Attention: on place un réplicat par colonne
+ /* Attention: on place un réplicat par colonne
* (et pas par ligne) */
mat[j][i]=d;
}
@@ -39,7 +39,7 @@ ensemble_t alloc_replicat(int nb_chi2)
}
return rep;
err:
- fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
+ fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
exit(1);
}
@@ -57,7 +57,7 @@ ensemble_t alloc_ensemble(int nb_sample)
}
return ens;
err:
- fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
+ fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
exit(1);
}
@@ -79,7 +79,7 @@ matrice_t alloc_matrice(int nb_sample, int nb_chi2)
return mat;
err:
- fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
+ fprintf(stderr, "Erreur d'allocation mémoire. Aborting\n");
exit(1);
}
@@ -93,7 +93,8 @@ void free_matrice(matrice_t mat, int nb_sample, int nb_chi2)
free(mat);
}
-inline static int count_superieur(ensemble_t ens, datatype_t val_ref, int nb_sample)
+inline static int count_superieur(ensemble_t ens, datatype_t val_ref,
+ int nb_sample)
{
int i, count=0;
for (i=0; i< nb_sample; i++) {
@@ -105,6 +106,19 @@ inline static int count_superieur(ensemble_t ens, datatype_t val_ref, int nb_sam
return count;
}
+inline static int count_inferieur(ensemble_t ens, datatype_t val_ref,
+ int nb_sample)
+{
+ int i, count=0;
+ for (i=0; i< nb_sample; i++) {
+ if (ens[i]<=val_ref) {
+ count++;
+ }
+ }
+ //printf( "count=%i\n", count);
+ return count;
+}
+
inline static datatype_t pval_min(replicat_t rep, int nb_chi2)
{
int i;
@@ -119,38 +133,50 @@ inline static datatype_t pval_min(replicat_t rep, int nb_chi2)
datatype_t calcul(int nb_sample, int nb_chi2, matrice_t mat, replicat_t rep)
{
- int i, j;
- ensemble_t ens_min_pval;
datatype_t min;
- datatype_t local[nb_chi2];
+ ensemble_t ens_min_pval;
ens_min_pval=alloc_ensemble(nb_sample);
+
+ min=calcul_distrib_pmin(nb_sample, nb_chi2, mat, rep, ens_min_pval);
+
+ free_ensemble(ens_min_pval);
+ return min;
+}
+datatype_t calcul_distrib_pmin(int nb_sample, int nb_chi2, matrice_t mat,
+ replicat_t rep, ensemble_t ens_min_pval)
+{
+ int i, j;
+ datatype_t min;
+ datatype_t local[nb_chi2];
+
i=0;
for (j=0; j<nb_chi2; j++) {
rep[j]=CALC_PVAL(count_superieur(mat[j], mat[j][i], nb_sample),
nb_sample);
//printf("cal rep[%i]=%lf\n", j, rep[j]);
}
- ens_min_pval[i]=-pval_min(rep, nb_chi2);
+ /* i is still 0 here */
+ ens_min_pval[i]=pval_min(rep, nb_chi2);
//printf("pmin for sample %i: "CONV"\n", i, ens_min_pval[i]);
for (i=1; i<nb_sample; i++) {
for (j=0; j<nb_chi2; j++) {
local[j]=CALC_PVAL(count_superieur(mat[j], mat[j][i],
- nb_sample),
- nb_sample);
+ nb_sample),
+ nb_sample);
}
- ens_min_pval[i]=-pval_min(local, nb_chi2);
+ ens_min_pval[i]=pval_min(local, nb_chi2);
//printf("pmin for sample %i: "CONV"\n", i, ens_min_pval[i]);
}
- min=CALC_PVAL(count_superieur(ens_min_pval, ens_min_pval[0],
+ min=CALC_PVAL(count_inferieur(ens_min_pval, ens_min_pval[0],
nb_sample),
nb_sample);
-
return min;
}
+
#ifdef MAIN_PROG
int main(int argc, char *argv[])
{
@@ -174,7 +200,7 @@ int main(int argc, char *argv[])
for (j=0; j<nb_chi2; j++) {
printf("chi2 niveau %i, pval nc "CONV"\n", j+1, rep[j]);
}
- printf("pmin corrigé: "CONV"\n", min);
+ printf("pmin corrigé: "CONV"\n", min);
free_matrice(mat, nb_sample, nb_chi2);
free_replicat(rep);
diff --git a/CUtils/c_sources/double_permutation.h b/CUtils/c_sources/double_permutation.h
index 669089543f54828be7d28968c3db160eba09c7cd..0231b2f43c237d506c008b97432eb899e9d92f74 100644
--- a/CUtils/c_sources/double_permutation.h
+++ b/CUtils/c_sources/double_permutation.h
@@ -11,7 +11,12 @@ typedef datatype_t *replicat_t;
ensemble_t alloc_replicat(int nb_chi2);
void free_replicat(ensemble_t rep);
+ensemble_t alloc_ensemble(int nb_sample);
+void free_ensemble(ensemble_t ens);
+
matrice_t alloc_matrice(int nb_sample, int nb_chi2);
void free_matrice(matrice_t mat, int nb_sample, int nb_chi2);
datatype_t calcul(int nb_sample, int nb_chi2, matrice_t mat, replicat_t rep);
+datatype_t calcul_distrib_pmin(int nb_sample, int nb_chi2, matrice_t mat,
+ replicat_t rep, ensemble_t ens_min_pval);
diff --git a/CUtils/t/double_permutation.t b/CUtils/t/double_permutation.t
index 7d0f4be137f4ba890d07c67fcdc3ddae3c74ebce..ca52835fd04cd2d7a0e393d37c9c519c48123ea0 100644
--- a/CUtils/t/double_permutation.t
+++ b/CUtils/t/double_permutation.t
@@ -22,5 +22,8 @@ $res=ALTree::CUtils::double_permutation(10, 2,
);
is_deeply($res, { "pmin" => 0.2,
- "chi2" => [ 0.8, 0.1] }, "double_permutation");
+ "chi2" => [ 0.8, 0.1],
+ "distrib_pmin" =>
+ [0.1, 0.1, 0, 0.2, 0.2, 0.8, 0.6, 0.5, 0.4, 0.3]
+ }, "double_permutation");