From 9a6881ed8a16c31a3dbe995e3b1af76c1db6b5a0 Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Wed, 11 Dec 2019 16:44:28 +0100
Subject: [PATCH] Add comments, TODOs...

---
 pkg/R/optimParams.R | 15 +++++++++++----
 pkg/src/functions.c |  4 ++++
 2 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index 505b665..ecfae7f 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -115,13 +115,21 @@ setRefClass(
       c(L$p[1:(K-1)], as.double(L$β), L$b)
     },
 
+    #TODO: compare with R version?
+    #D <- diag(d) #matrix of ej vectors
+    #Y * X
+    #Y * ( t( apply(X, 1, function(row) row %o% row) ) - Reduce('+', lapply(1:d, function(j) as.double(D[j,] %o% D[j,])), rep(0, d*d)))
+    #Y * ( t( apply(X, 1, function(row) row %o% row %*% row) ) - Reduce('+', lapply(1:d, function(j) ), rep(0, d*d*d)))
     computeW = function(θ)
     {
       #require(MASS)
       dd <- d + d^2 + d^3
-      W <<- MASS::ginv( matrix( .C("Compute_Omega",
-        X=as.double(X), Y=Y, M=Moments(θ), pn=as.integer(n), pd=as.integer(d),
-        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) )
+      M <- Moments(θ)
+      Omega <- matrix( .C("Compute_Omega",
+        X=as.double(X), Y=as.double(Y), M=as.double(M),
+        pn=as.integer(n), pd=as.integer(d),
+        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
+      W <<- MASS::ginv(Omega, tol=1e-4)
       NULL #avoid returning W
     },
 
@@ -248,7 +256,6 @@ setRefClass(
         θ0$b = rep(0, K)
       else if (any(is.na(θ0$b)))
         stop("θ0$b cannot have missing values")
-
       # TODO: stopping condition? N iterations? Delta <= epsilon ?
       for (loop in 1:10)
       {
diff --git a/pkg/src/functions.c b/pkg/src/functions.c
index e534c57..f251eb3 100644
--- a/pkg/src/functions.c
+++ b/pkg/src/functions.c
@@ -54,6 +54,8 @@ void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3)
   }
 }
 
+#include <stdio.h>
+
 // W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) )
 // with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments
 void Compute_Omega(double* X, double* Y, double* M, int* pn, int* pd, double* W)
@@ -97,6 +99,8 @@ void Compute_Omega(double* X, double* Y, double* M, int* pn, int* pd, double* W)
         g[j] -= Y[i] * X[mi(i,idx1,n,d)];
       g[j] += Y[i] * X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)]*X[mi(i,idx3,n,d)] - M[j];
     }
+
+    // TODO: 1/n des gj empirique doit tendre vers 0
     // Add 1/n t(gi) %*% gi to W
     for (int j=0; j<dim; j++)
     {
-- 
2.44.0