X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=src%2Fsources%2FEMGLLF.c;h=087116be27573bc679b8b41c0d264ae0ece95700;hp=14394164b62b8450cf06d4e2bcc843ded95a66f0;hb=9ff729fb6afff5bed327fa8619138fd3b6f6f13b;hpb=e8bb47649b9f4e625ef050aae096d23cfc058163 diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 1439416..087116b 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -5,23 +5,23 @@ // TODO: don't recompute indexes every time...... void EMGLLF_core( // IN parameters - const float* phiInit, // parametre initial de moyenne renormalisé - const float* rhoInit, // parametre initial de variance renormalisé - const float* piInit, // parametre initial des proportions - const float* 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 - float gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif - float lambda, // valeur du paramètre de régularisation du Lasso - const float* X, // régresseurs - const float* Y, // réponse - float 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) - float* phi, // parametre de moyenne renormalisé, calculé par l'EM - float* rho, // parametre de variance renormalisé, calculé par l'EM - float* pi, // parametre des proportions renormalisé, calculé par l'EM - float* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres - float* 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_core( //Other local variables //NOTE: variables order is always [maxi],n,p,m,k - float* gam = (float*)malloc(n*k*sizeof(float)); + Real* gam = (Real*)malloc(n*k*sizeof(Real)); copyArray(gamInit, gam, n*k); - float* b = (float*)malloc(k*sizeof(float)); - float* Phi = (float*)malloc(p*m*k*sizeof(float)); - float* Rho = (float*)malloc(m*m*k*sizeof(float)); - float* Pi = (float*)malloc(k*sizeof(float)); - float* gam2 = (float*)malloc(k*sizeof(float)); - float* pi2 = (float*)malloc(k*sizeof(float)); - float* Gram2 = (float*)malloc(p*p*k*sizeof(float)); - float* ps = (float*)malloc(m*k*sizeof(float)); - float* nY2 = (float*)malloc(m*k*sizeof(float)); - float* ps1 = (float*)malloc(n*m*k*sizeof(float)); - float* ps2 = (float*)malloc(p*m*k*sizeof(float)); - float* nY21 = (float*)malloc(n*m*k*sizeof(float)); - float* Gam = (float*)malloc(n*k*sizeof(float)); - float* X2 = (float*)malloc(n*p*k*sizeof(float)); - float* Y2 = (float*)malloc(n*m*k*sizeof(float)); + 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); - float* YiRhoR = (float*)malloc(m*sizeof(float)); - float* XiPhiR = (float*)malloc(m*sizeof(float)); - float dist = 0.; - float 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; - float EPS = 1e-15; - float* dotProducts = (float*)malloc(k*sizeof(float)); + 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_core( //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r); for (int u=0; u - float 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); - float 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_core( } //sum(pen(ite,:)) - float sumPen = 0.0; + Real sumPen = 0.0; for (int r=0; r Dist1) Dist1 = tmpDist; @@ -353,14 +353,14 @@ void EMGLLF_core( } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - float Dist2 = 0.0; + Real Dist2 = 0.0; for (int u=0; u Dist2) Dist2 = tmpDist; @@ -368,12 +368,12 @@ void EMGLLF_core( } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - float Dist3 = 0.0; + Real Dist3 = 0.0; for (int u=0; u Dist3) Dist3 = tmpDist; }