HugoR<-function(F,N,w0,w1,t_max){ # F ... the mxn matrix of regressors of LRM, where # m is the dimension of the parameter, n is the size of the design space. # N ... the number of trials # w0 ... the nx1 vector representing the design to be augmented. # w1 ... the nx1 vector representing the initial design. # t_max ... upper limit on the time taken by the computation. back_max<-16 n_round<-9 V_max<-500000 eps<-1/1000000 m<-dim(F)[1] n<-dim(F)[2] u<-rep(1,n) ############subfunctions############################ Phi<-function(w){ # D-optimality det(F%*%diag(w)%*%t(F)+eps*diag(m)) } ############### explore_up<-function(){ wei<-w val<-(-u) for (i in 1:n){ if (d[i]>eps){ wei[i]<-w[i]+1 atribei<-signif(Phi(wei),n_round) if (!is.element(atribei,V)){ dei<-(N-sum(wei))*u val[i]<-Phi(wei+1/n*dei) } wei[i]<-w[i] } } up_index<-which.max(val) max_val<-val[up_index] if (max_val<(-eps)){ up_index<-0 } up_index } ############### explore_down<-function(){ wei<-w val<-(-u) for (i in 1:n){ if (w[i]>w0[i]){ wei[i]<-w[i]-1 atribei<-signif(Phi(wei),n_round) if (!is.element(atribei,V)){ dei<-(N-sum(wei))*u val[i]<-Phi(wei+1/n*dei) } wei[i]<-w[i] } } down_index<-which.max(val) max_val<-val[down_index] if (max_val<(-eps)){ down_index<-0 } down_index } ################ random_step<-function(){ bloc<-0 wloc<-w success<-FALSE while(!success){ dir<-sample(c(0,1),1) if(dir){ up_index<-sample(1:n,1) if (d[up_index]>0){ wloc[up_index]<-wloc[up_index]+1 success<-TRUE } } else{ down_index<-sample(1:n,1) if (w[down_index]>w0[down_index]){ wloc[down_index]<-wloc[down_index]-1 success<-TRUE bloc<-1 } } } list(wloc=wloc,bloc=bloc) } ############MAIN BODY################################# start <- proc.time () if((length(w1)!=n)||(min(w1-w0)<0)||(sum(w1)>N)){ w<-w0 d<-(N-sum(w))*u maximal<-FALSE while(!maximal){ val<-(d>0) if(sum(val)==0){ maximal<-TRUE } else{ i<-sample(1:n,1,prob=val) w[i]<-w[i]+1 d<-(N-sum(w))*u } } atrib<-signif(Phi(w),n_round) } else{ w<-w1 d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) } w_best<-w d_best<-d atrib_best<-atrib V<-NULL back_no<-0 finish<-FALSE Phi_best<-0 while(!finish){ if(!is.element(atrib,V)){ V<-union(atrib,V) up_index<-explore_up() if(up_index){ w[up_index]<-w[up_index]+1 d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) } else{ if (sum(w)==N){ Phi_w_cmp<-Phi(w) if (Phi_w_cmp>Phi_best){ w_best<-w d_best<-d atrib_best<-atrib Phi_best<-Phi_w_cmp back_no<-0 } } down_index<-explore_down() if(down_index){ w[down_index]<-w[down_index]-1 d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) back_no<-back_no+1 } else{ res<-random_step() w<-res$wloc; back_no<-back_no+res$bloc d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) } } } else{ down_index<-explore_down() if(down_index){ w[down_index]<-w[down_index]-1 d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) back_no<-back_no+1 } else{ up_index<-explore_up() if(up_index){ w[up_index]<-w[up_index]+1 d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) } else{ res<-random_step() w<-res$wloc; back_no<-back_no+res$bloc d<-(N-sum(w))*u atrib<-signif(Phi(w),n_round) } } } if (back_no>back_max){ w<-w_best d<-d_best atrib<-atrib_best back_no<-0 } ptm<-proc.time () - start if((ptm>t_max)||(length(V)>=V_max)){ finish<-TRUE } } Phi_best<-(det(F%*%diag(w_best)%*%t(F)))^(1/(m-1)) out<-list(design=w_best, value=Phi_best, Vl=length(V)) out }