X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=src%2Fsources%2FEMGLLF.c;h=087116be27573bc679b8b41c0d264ae0ece95700;hb=9ff729fb6afff5bed327fa8619138fd3b6f6f13b;hp=6966e9c7c15c6c6686a6239b1551cdaed6e3d4a8;hpb=8e92c49c15bdacebf46190e7c8279bd227873928;p=valse.git diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 6966e9c..087116b 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -3,25 +3,25 @@ #include // TODO: don't recompute indexes every time...... -void EMGLLF( +void EMGLLF_core( // IN parameters - const double* phiInit, // parametre initial de moyenne renormalisé - const double* rhoInit, // parametre initial de variance renormalisé - const double* piInit, // parametre initial des proportions - const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon + const Real* phiInit, // parametre initial de moyenne renormalisé + const Real* rhoInit, // parametre initial de variance renormalisé + const Real* piInit, // parametre initial des proportions + const Real* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon int mini, // nombre minimal d'itérations dans l'algorithme EM int maxi, // nombre maximal d'itérations dans l'algorithme EM - double gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif - double lambda, // valeur du paramètre de régularisation du Lasso - const double* X, // régresseurs - const double* Y, // réponse - double tau, // seuil pour accepter la convergence + Real gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif + Real lambda, // valeur du paramètre de régularisation du Lasso + const Real* X, // régresseurs + const Real* Y, // réponse + Real tau, // seuil pour accepter la convergence // OUT parameters (all pointers, to be modified) - double* phi, // parametre de moyenne renormalisé, calculé par l'EM - double* rho, // parametre de variance renormalisé, calculé par l'EM - double* pi, // parametre des proportions renormalisé, calculé par l'EM - double* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres - double* S, + Real* phi, // parametre de moyenne renormalisé, calculé par l'EM + Real* rho, // parametre de variance renormalisé, calculé par l'EM + Real* pi, // parametre des proportions renormalisé, calculé par l'EM + Real* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres + Real* S, // additional size parameters int n, // nombre d'echantillons int p, // nombre de covariables @@ -37,32 +37,32 @@ void EMGLLF( //Other local variables //NOTE: variables order is always [maxi],n,p,m,k - double* gam = (double*)malloc(n*k*sizeof(double)); + Real* gam = (Real*)malloc(n*k*sizeof(Real)); copyArray(gamInit, gam, n*k); - double* b = (double*)malloc(k*sizeof(double)); - double* Phi = (double*)malloc(p*m*k*sizeof(double)); - double* Rho = (double*)malloc(m*m*k*sizeof(double)); - double* Pi = (double*)malloc(k*sizeof(double)); - double* gam2 = (double*)malloc(k*sizeof(double)); - double* pi2 = (double*)malloc(k*sizeof(double)); - double* Gram2 = (double*)malloc(p*p*k*sizeof(double)); - double* ps = (double*)malloc(m*k*sizeof(double)); - double* nY2 = (double*)malloc(m*k*sizeof(double)); - double* ps1 = (double*)malloc(n*m*k*sizeof(double)); - double* ps2 = (double*)malloc(p*m*k*sizeof(double)); - double* nY21 = (double*)malloc(n*m*k*sizeof(double)); - double* Gam = (double*)malloc(n*k*sizeof(double)); - double* X2 = (double*)malloc(n*p*k*sizeof(double)); - double* Y2 = (double*)malloc(n*m*k*sizeof(double)); + Real* b = (Real*)malloc(k*sizeof(Real)); + Real* Phi = (Real*)malloc(p*m*k*sizeof(Real)); + Real* Rho = (Real*)malloc(m*m*k*sizeof(Real)); + Real* Pi = (Real*)malloc(k*sizeof(Real)); + Real* gam2 = (Real*)malloc(k*sizeof(Real)); + Real* pi2 = (Real*)malloc(k*sizeof(Real)); + Real* Gram2 = (Real*)malloc(p*p*k*sizeof(Real)); + Real* ps = (Real*)malloc(m*k*sizeof(Real)); + Real* nY2 = (Real*)malloc(m*k*sizeof(Real)); + Real* ps1 = (Real*)malloc(n*m*k*sizeof(Real)); + Real* ps2 = (Real*)malloc(p*m*k*sizeof(Real)); + Real* nY21 = (Real*)malloc(n*m*k*sizeof(Real)); + Real* Gam = (Real*)malloc(n*k*sizeof(Real)); + Real* X2 = (Real*)malloc(n*p*k*sizeof(Real)); + Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real)); gsl_matrix* matrix = gsl_matrix_alloc(m, m); gsl_permutation* permutation = gsl_permutation_alloc(m); - double* YiRhoR = (double*)malloc(m*sizeof(double)); - double* XiPhiR = (double*)malloc(m*sizeof(double)); - double dist = 0.; - double dist2 = 0.; + Real* YiRhoR = (Real*)malloc(m*sizeof(Real)); + Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); + Real dist = 0.; + Real dist2 = 0.; int ite = 0; - double EPS = 1e-15; - double* dotProducts = (double*)malloc(k*sizeof(double)); + Real EPS = 1e-15; + Real* dotProducts = (Real*)malloc(k*sizeof(Real)); while (ite < mini || (ite < maxi && (dist >= tau || dist2 >= sqrt(tau)))) { @@ -90,7 +90,7 @@ void EMGLLF( //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r); for (int u=0; u - double dotProduct = 0.0; + Real dotProduct = 0.0; for (int u=0; udata[u*m+v] = rho[ai(u,v,r,m,m,k)]; } gsl_linalg_LU_decomp(matrix, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrix, signum); + Real detRhoR = gsl_linalg_LU_det(matrix, signum); Gam[mi(i,r,n,k)] = pi[r] * detRhoR * exp(-0.5*dotProducts[r] + shift); sumLLF1 += Gam[mi(i,r,n,k)] / pow(2*M_PI,m/2.0); @@ -327,7 +327,7 @@ void EMGLLF( } //sum(pen(ite,:)) - double sumPen = 0.0; + Real sumPen = 0.0; for (int r=0; r Dist1) Dist1 = tmpDist; @@ -353,14 +353,14 @@ void EMGLLF( } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - double Dist2 = 0.0; + Real Dist2 = 0.0; for (int u=0; u Dist2) Dist2 = tmpDist; @@ -368,12 +368,12 @@ void EMGLLF( } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - double Dist3 = 0.0; + Real Dist3 = 0.0; for (int u=0; u Dist3) Dist3 = tmpDist; }