| 1 | #include <stdlib.h> |
| 2 | |
| 3 | // Index matrix (by columns) |
| 4 | int mi(int i, int j, int d1, int d2) |
| 5 | { |
| 6 | return j*d1 + i; |
| 7 | } |
| 8 | |
| 9 | // Index 3-tensor (by columns, matrices ordered by last dim) |
| 10 | int ti(int i, int j, int k, int d1, int d2, int d3) |
| 11 | { |
| 12 | return k*d1*d2 + j*d1 + i; |
| 13 | } |
| 14 | |
| 15 | // Empirical cross-moment of order 2 between X size nxd and Y size n |
| 16 | void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2) |
| 17 | { |
| 18 | int n=*pn, d=*pd; |
| 19 | //double* M2 = (double*)calloc(d*d,sizeof(double)); |
| 20 | |
| 21 | // M2 = E[Y*X^*2] - E[Y*e^*2] = E[Y (X^*2 - I)] |
| 22 | for (int j=0; j<d; j++) |
| 23 | { |
| 24 | for (int i=0; i<n; i++) |
| 25 | { |
| 26 | M2[mi(j,j,d,d)] -= Y[i] / n; |
| 27 | for (int k=0; k<d; k++) |
| 28 | M2[mi(j,k,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)] / n; |
| 29 | } |
| 30 | } |
| 31 | } |
| 32 | |
| 33 | // Empirical cross-moment of order 3 between X size nxd and Y size n |
| 34 | void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3) |
| 35 | { |
| 36 | int n=*pn, d=*pd; |
| 37 | //double* M3 = (double*)calloc(d*d*d,sizeof(double)); |
| 38 | |
| 39 | // M3 = E[Y*X^*3] - E[Y*e*X*e] - E[Y*e*e*X] - E[Y*X*e*e] |
| 40 | for (int j=0; j<d; j++) |
| 41 | { |
| 42 | for (int k=0; k<d; k++) |
| 43 | { |
| 44 | for (int i=0; i<n; i++) |
| 45 | { |
| 46 | double tensor_elt = Y[i]*X[mi(i,k,n,d)] / n; |
| 47 | M3[ti(j,k,j,d,d,d)] -= tensor_elt; |
| 48 | M3[ti(j,j,k,d,d,d)] -= tensor_elt; |
| 49 | M3[ti(k,j,j,d,d,d)] -= tensor_elt; |
| 50 | for (int o=0; o<d; o++) |
| 51 | M3[ti(j,k,o,d,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)]*X[mi(i,o,n,d)] / n; |
| 52 | } |
| 53 | } |
| 54 | } |
| 55 | } |
| 56 | |
| 57 | #include <stdio.h> |
| 58 | |
| 59 | // W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) ) |
| 60 | // with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments |
| 61 | //void Compute_Omega(double* X, int* Y, double* M, int* pn, int* pd, double* W) |
| 62 | //{ |
| 63 | // int n=*pn, d=*pd; |
| 64 | // int dim = d + d*d + d*d*d; |
| 65 | // //double* W = (double*)malloc(dim*dim*sizeof(double)); |
| 66 | // |
| 67 | // // (Re)Initialize W: |
| 68 | // for (int j=0; j<dim; j++) |
| 69 | // { |
| 70 | // for (int k=0; k<dim; k++) |
| 71 | // W[j*dim+k] = 0.0; |
| 72 | // } |
| 73 | // double* g = (double*)malloc(dim*sizeof(double)); |
| 74 | // for (int i=0; i<n; i++) |
| 75 | // { |
| 76 | // // g == gi: |
| 77 | // for (int j=0; j<d; j++) |
| 78 | // g[j] = Y[i] * X[mi(i,j,n,d)] - M[j]; |
| 79 | // for (int j=d; j<d+(d*d); j++) |
| 80 | // { |
| 81 | // int idx1 = (j-d) % d; //num row |
| 82 | // int idx2 = ((j-d) - idx1) / d; //num col |
| 83 | // g[j] = 0.0; |
| 84 | // if (idx1 == idx2) |
| 85 | // g[j] -= Y[i]; |
| 86 | // g[j] += Y[i] * X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)] - M[j]; |
| 87 | // } |
| 88 | // for (int j=d+d*d; j<dim; j++) |
| 89 | // { |
| 90 | // int idx1 = (j-d-d*d) % d; //num row |
| 91 | // int idx2 = ((j-d-d*d - idx1) / d) %d; //num col |
| 92 | // int idx3 = (((j-d-d*d - idx1) / d) - idx2) / d; //num "depth" |
| 93 | // g[j] = 0.0; |
| 94 | // if (idx1 == idx2) |
| 95 | // g[j] -= Y[i] * X[mi(i,idx3,n,d)]; |
| 96 | // if (idx1 == idx3) |
| 97 | // g[j] -= Y[i] * X[mi(i,idx2,n,d)]; |
| 98 | // if (idx2 == idx3) |
| 99 | // g[j] -= Y[i] * X[mi(i,idx1,n,d)]; |
| 100 | // g[j] += Y[i] * X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)]*X[mi(i,idx3,n,d)] - M[j]; |
| 101 | // } |
| 102 | // // Add 1/n t(gi) %*% gi to W |
| 103 | // for (int j=0; j<dim; j++) |
| 104 | // { |
| 105 | // for (int k=0; k<dim; k++) |
| 106 | // W[j*dim+k] += g[j] * g[k] / n; |
| 107 | // } |
| 108 | // } |
| 109 | // free(g); |
| 110 | //} |
| 111 | |
| 112 | // Optimisation attempt: |
| 113 | void Compute_Omega(double* X, int* Y, double* M, int* pn, int* pd, double* W) |
| 114 | { |
| 115 | int n=*pn, d=*pd; |
| 116 | int dim = d + d*d + d*d*d; |
| 117 | //double* W = (double*)malloc(dim*dim*sizeof(double)); |
| 118 | |
| 119 | // (Re)Initialize W: |
| 120 | for (int j=0; j<dim; j++) |
| 121 | { |
| 122 | for (int k=0; k<dim; k++) |
| 123 | W[j*dim+k] = 0.0; |
| 124 | } |
| 125 | double* g = (double*)malloc(dim*sizeof(double)); |
| 126 | for (int i=0; i<n; i++) |
| 127 | { |
| 128 | printf("i: %i\n",i); |
| 129 | // g == gi: |
| 130 | for (int j=0; j<d; j++) |
| 131 | g[j] = (Y[i] ? X[mi(i,j,n,d)] - M[j] : 0.0); |
| 132 | for (int j=d; j<d+(d*d); j++) |
| 133 | { |
| 134 | int idx1 = (j-d) % d; //num row |
| 135 | int idx2 = ((j-d) - idx1) / d; //num col |
| 136 | g[j] = 0.0; |
| 137 | if (Y[i]) |
| 138 | { |
| 139 | if (idx1 == idx2) |
| 140 | g[j]--; |
| 141 | g[j] += X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)] - M[j]; |
| 142 | } |
| 143 | } |
| 144 | for (int j=d+d*d; j<dim; j++) |
| 145 | { |
| 146 | int idx1 = (j-d-d*d) % d; //num row |
| 147 | int idx2 = ((j-d-d*d - idx1) / d) %d; //num col |
| 148 | int idx3 = (((j-d-d*d - idx1) / d) - idx2) / d; //num "depth" |
| 149 | g[j] = 0.0; |
| 150 | if (Y[i]) |
| 151 | { |
| 152 | if (idx1 == idx2) |
| 153 | g[j] -= X[mi(i,idx3,n,d)]; |
| 154 | if (idx1 == idx3) |
| 155 | g[j] -= X[mi(i,idx2,n,d)]; |
| 156 | if (idx2 == idx3) |
| 157 | g[j] -= X[mi(i,idx1,n,d)]; |
| 158 | g[j] += X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)]*X[mi(i,idx3,n,d)] - M[j]; |
| 159 | } |
| 160 | } |
| 161 | // Add 1/n t(gi) %*% gi to W |
| 162 | for (int j=0; j<dim; j++) |
| 163 | { |
| 164 | for (int k=0; k<dim; k++) |
| 165 | W[j*dim+k] += g[j] * g[k] / n; |
| 166 | } |
| 167 | } |
| 168 | free(g); |
| 169 | } |