remove also NY21 as in R code
authorBenjamin Auder <benjamin.auder@somewhere>
Mon, 20 Mar 2017 03:18:01 +0000 (04:18 +0100)
committerBenjamin Auder <benjamin.auder@somewhere>
Mon, 20 Mar 2017 03:18:01 +0000 (04:18 +0100)
pkg/src/sources/EMGLLF.c

index e41fe3c..a028919 100644 (file)
@@ -54,7 +54,6 @@ void EMGLLF_core(
        const Real EPS = 1e-15;
        // Additional (not at this place, in R file)
        Real* gam2 = (Real*)malloc(k*sizeof(Real));
-       Real* nY21 = (Real*)malloc(n*m*k*sizeof(Real));
        Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
        gsl_matrix* matrix = gsl_matrix_alloc(m, m);
        gsl_permutation* permutation = gsl_permutation_alloc(m);
@@ -216,18 +215,17 @@ void EMGLLF_core(
                                                dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)];
                                        //ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
                                        ps1[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * dotProduct;
-                                       nY21[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)];
                                }
                                //ps[mm,r] = sum(ps1[,mm,r])
                                Real sumPs1 = 0.;
                                for (int u=0; u<n; u++)
                                        sumPs1 += ps1[ai(u,mm,r,n,m,k)];
                                ps[mi(mm,r,m,k)] = sumPs1;
-                               //nY2[mm,r] = sum(nY21[,mm,r])
-                               Real sumNy21 = 0.;
+                               //nY2[mm,r] = sum(nY2[,mm,r])
+                               Real sumNy2 = 0.;
                                for (int u=0; u<n; u++)
-                                       sumNy21 += nY21[ai(u,mm,r,n,m,k)];
-                               nY2[mi(mm,r,m,k)] = sumNy21;
+                                       sumNy2 += nY2[ai(u,mm,r,n,m,k)];
+                               nY2[mi(mm,r,m,k)] = sumNy2;
                                //rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
                                rho[ai(mm,mm,r,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
                                        + 4*nY2[mi(mm,r,m,k)] * gam2[r] ) ) / (2*nY2[mi(mm,r,m,k)]);
@@ -380,6 +378,8 @@ void EMGLLF_core(
                ite++;
        }
 
+       //TODO: affec = apply(gam, 1,which.max) à traduire, et retourner affec aussi.
+
        //free memory
        free(b);
        free(gam);
@@ -390,7 +390,6 @@ void EMGLLF_core(
        free(ps);
        free(nY2);
        free(ps1);
-       free(nY21);
        free(Gram2);
        free(ps2);
        gsl_matrix_free(matrix);