fix arrays reading in C, add type Real for double or float
[valse.git] / R / main.R
1 #' @useDynLib valse
2
3 Valse = setRefClass(
4 Class = "Valse",
5
6 fields = c(
7 # User defined
8
9 # regression data (size n*p, where n is the number of observations,
10 # and p is the number of regressors)
11 X = "matrix",
12 # response data (size n*m, where n is the number of observations,
13 # and m is the number of responses)
14 Y = "matrix",
15
16 # Optionally user defined (some default values)
17
18 # power in the penalty
19 gamma = "numeric",
20 # minimum number of iterations for EM algorithm
21 mini = "integer",
22 # maximum number of iterations for EM algorithm
23 maxi = "integer",
24 # threshold for stopping EM algorithm
25 eps = "numeric",
26 # minimum number of components in the mixture
27 kmin = "integer",
28 # maximum number of components in the mixture
29 kmax = "integer",
30 rangmin = "integer",
31 rangmax = "integer",
32
33 # Computed through the workflow
34
35 # initialisation for the reparametrized conditional mean parameter
36 phiInit = "numeric",
37 # initialisation for the reparametrized variance parameter
38 rhoInit = "numeric",
39 # initialisation for the proportions
40 piInit = "numeric",
41 # initialisation for the allocations probabilities in each component
42 tauInit = "numeric",
43 # values for the regularization parameter grid
44 gridLambda = "numeric",
45 # je ne crois pas vraiment qu'il faille les mettre en sortie, d'autant plus qu'on construit
46 # une matrice A1 et A2 pour chaque k, et elles sont grandes, donc ca coute un peu cher ...
47 A1 = "integer",
48 A2 = "integer",
49 # collection of estimations for the reparametrized conditional mean parameters
50 Phi = "numeric",
51 # collection of estimations for the reparametrized variance parameters
52 Rho = "numeric",
53 # collection of estimations for the proportions parameters
54 Pi = "numeric",
55
56 #immutable (TODO:?)
57 seuil = "numeric"
58 ),
59
60 methods = list(
61 #######################
62 #initialize main object
63 #######################
64 initialize = function(X,Y,...)
65 {
66 "Initialize Valse object"
67
68 callSuper(...)
69
70 X <<- X
71 Y <<- Y
72 gamma <<- ifelse (hasArg("gamma"), gamma, 1.)
73 mini <<- ifelse (hasArg("mini"), mini, as.integer(5))
74 maxi <<- ifelse (hasArg("maxi"), maxi, as.integer(10))
75 eps <<- ifelse (hasArg("eps"), eps, 1e-6)
76 kmin <<- ifelse (hasArg("kmin"), kmin, as.integer(2))
77 kmax <<- ifelse (hasArg("kmax"), kmax, as.integer(3))
78 rangmin <<- ifelse (hasArg("rangmin"), rangmin, as.integer(2))
79 rangmax <<- ifelse (hasArg("rangmax"), rangmax, as.integer(3))
80 seuil <<- 1e-15 #immutable (TODO:?)
81 },
82
83 ##################################
84 #core workflow: compute all models
85 ##################################
86
87 initParameters = function(k)
88 {
89 "Parameters initialization"
90
91 #smallEM initializes parameters by k-means and regression model in each component,
92 #doing this 20 times, and keeping the values maximizing the likelihood after 10
93 #iterations of the EM algorithm.
94 init = initSmallEM(k,X,Y,eps)
95 phiInit <<- init$phi0
96 rhoInit <<- init$rho0
97 piInit <<- init$pi0
98 tauInit <<- init$tau0
99 },
100
101 computeGridLambda = function()
102 {
103 "computation of the regularization grid"
104 #(according to explicit formula given by EM algorithm)
105
106 gridLambda <<- gridLambda(phiInit,rhoInit,piInit,tauInit,X,Y,gamma,mini,maxi,eps)
107 },
108
109 computeRelevantParameters = function()
110 {
111 "Compute relevant parameters"
112
113 #select variables according to each regularization parameter
114 #from the grid: A1 corresponding to selected variables, and
115 #A2 corresponding to unselected variables.
116 params = selectiontotale(
117 phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps)
118 A1 <<- params$A1
119 A2 <<- params$A2
120 Rho <<- params$Rho
121 Pi <<- params$Pi
122 },
123
124 runProcedure1 = function()
125 {
126 "Run procedure 1 [EMGLLF]"
127
128 #compute parameter estimations, with the Maximum Likelihood
129 #Estimator, restricted on selected variables.
130 return ( constructionModelesLassoMLE(
131 phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps,A1,A2) )
132 },
133
134 runProcedure2 = function()
135 {
136 "Run procedure 2 [EMGrank]"
137
138 #compute parameter estimations, with the Low Rank
139 #Estimator, restricted on selected variables.
140 return ( constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,eps,
141 A1,rangmin,rangmax) )
142 },
143
144 run = function()
145 {
146 "main loop: over all k and all lambda"
147
148 # Run the all procedure, 1 with the
149 #maximum likelihood refitting, and 2 with the Low Rank refitting.
150 p = dim(phiInit)[1]
151 m = dim(phiInit)[2]
152 for (k in kmin:kmax)
153 {
154 print(k)
155 initParameters(k)
156 computeGridLambda()
157 computeRelevantParameters()
158 if (procedure == 1)
159 {
160 r1 = runProcedure1()
161 Phi2 = Phi
162 Rho2 = Rho
163 Pi2 = Pi
164 p = ncol(X)
165 m = ncol(Y)
166 if (is.null(dim(Phi2))) #test was: size(Phi2) == 0
167 {
168 Phi[,,1:k] <<- r1$phi
169 Rho[,,1:k] <<- r1$rho
170 Pi[1:k,] <<- r1$pi
171 } else
172 {
173 Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(r1$phi)[4]))
174 Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2
175 Phi[,,1:k,dim(Phi2)[4]+1] <<- r1$phi
176 Rho <<- array(0., dim=c(m,m,kmax,dim(Rho2)[4]+dim(r1$rho)[4]))
177 Rho[,,1:(dim(Rho2)[3]),1:(dim(Rho2)[4])] <<- Rho2
178 Rho[,,1:k,dim(Rho2)[4]+1] <<- r1$rho
179 Pi <<- array(0., dim=c(kmax,dim(Pi2)[2]+dim(r1$pi)[2]))
180 Pi[1:nrow(Pi2),1:ncol(Pi2)] <<- Pi2
181 Pi[1:k,ncol(Pi2)+1] <<- r1$pi
182 }
183 } else
184 {
185 phi = runProcedure2()$phi
186 Phi2 = Phi
187 if (dim(Phi2)[1] == 0)
188 {
189 Phi[,,1:k,] <<- phi
190 } else
191 {
192 Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(phi)[4]))
193 Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2
194 Phi[,,1:k,-(1:(dim(Phi2)[4]))] <<- phi
195 }
196 }
197 }
198 }
199
200 ##################################################
201 #TODO: pruning: select only one (or a few best ?!) model
202 ##################################################
203 #
204 # function[model] selectModel(
205 # #TODO
206 # #model = odel(...)
207 # end
208
209 )
210 )