### R code from vignette source 'MEIGOR-vignette.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: preliminaries ################################################### options(width=70, useFancyQuotes="UTF-8", prompt=" ", continue=" ") set.seed(123L) ################################################### ### code chunk number 2: installMEIGO2 (eval = FALSE) ################################################### ## if (!requireNamespace("BiocManager", quietly=TRUE)) ## install.packages("BiocManager") ## BiocManager::install("CellNOptR") ################################################### ### code chunk number 3: load_package ################################################### library(MEIGOR) ################################################### ### code chunk number 4: quickProblem (eval = FALSE) ################################################### ## library(MEIGOR) ## problem<-list(f=objective, x_L=rep(-1,2), x_U=rep(1,2)) ## opts<-list(maxeval=500, local_solver='DHC') ################################################### ### code chunk number 5: quickProblemSolve ################################################### ex3<-function(x,k1,k2,k3,k4){ f=-x[4]; #Equality constraints (declare them before the inequality ones) g<-rep(0,5); g[1]=x[4]-x[3]+x[2]-x[1]+k4*x[4]*x[6]; g[2]=x[1]-1+k1*x[1]*x[5]; g[3]=x[2]-x[1]+k2*x[2]*x[6]; g[4]=x[3]+x[1]-1+k3*x[3]*x[5]; #Inequality constraint g[5]=x[5]^0.5+x[6]^0.5; return(list(f=f,g=g)); } ################################################### ### code chunk number 6: multiStart (eval = FALSE) ################################################### ## Results_multistart<-MEIGO(problem,opts,algorithm="MULTISTART",p1,p2,...) ################################################### ### code chunk number 7: unConstrainedF ################################################### ex1 <- function(x){ y<-4*x[1]*x[1]-2.1*x[1]^4+1/3*x[1]^6+x[1]*x[2]-4*x[2]*x[2]+4*x[2]^4; return(y) } ################################################### ### code chunk number 8: unConstrainedSolve ################################################### #========================= #PROBLEM SPECIFICATIONS # ========================= problem<-list(f="ex1",x_L=rep(-1,2),x_U=rep(1,2)) opts<-list(maxeval=500, ndiverse=40, local_solver='DHC', local_finish='LBFGSB', local_iterprint=1) #========================= # END OF PROBLEM SPECIFICATIONS #========================= Results<-MEIGO(problem,opts,algorithm="ESS"); ################################################### ### code chunk number 9: ex2F ################################################### ex2<-function(x){ F=-x[1]-x[2]; g<-rep(0,2); g[1]<-x[2]-2*x[1]^4+8*x[1]^3-8*x[1]^2; g[2]<-x[2]-4*x[1]^4+32*x[1]^3-88*x[1]^2+96*x[1]; return(list(F=F,g=g)) } ################################################### ### code chunk number 10: ex2Solve ################################################### #========================= #PROBLEM SPECIFICATIONS #========================= problem<-list(f="ex2",x_L=rep(0,2),x_U=c(3,4), c_L=rep(-Inf,2), c_U=c(2,36)) opts<-list(maxeval=750, local_solver="DHC", local_n1=2, local_n2=3) #======================== #END OF PROBLEM SPECIFICATIONS # ======================== Results<-MEIGO(problem,opts,algorithm="ESS"); ################################################### ### code chunk number 11: ex3f ################################################### ex3<-function(x,k1,k2,k3,k4){ f=-x[4]; #Equality constraints g<-rep(0,5); g[1]=x[4]-x[3]+x[2]-x[1]+k4*x[4]*x[6]; g[2]=x[1]-1+k1*x[1]*x[5]; g[3]=x[2]-x[1]+k2*x[2]*x[6]; g[4]=x[3]+x[1]-1+k3*x[3]*x[5]; #Inequality constraint g[5]=x[5]^0.5+x[6]^0.5; return(list(f=f,g=g)); } ################################################### ### code chunk number 12: ex2Solve3 ################################################### #========================= #PROBLEM SPECIFICATIONS #========================= problem<-list(f="ex3",x_L=rep(0,6),x_U=c(rep(1,4),16,16), neq=4, c_L=-Inf, c_U=4) opts<-list(maxtime=5, local_solver='solnp') #========================= #END OF PROBLEM SPECIFICATIONS #========================= k1=0.09755988; k3=0.0391908; k2=0.99*k1; k4=0.9*k3; Results<-MEIGO(problem,opts,algorithm="ESS",k1,k2,k3,k4); ################################################### ### code chunk number 13: ex4f ################################################### ex4<-function(x){ F = x[2]^2 + x[3]^2 + 2.0*x[1]^2 + x[4]^2 - 5.0*x[2] - 5.0*x[3] - 21.0*x[1] + 7.0*x[4]; g<-rep(0,3); g[1] = x[2]^2 + x[3]^2 + x[1]^2 + x[4]^2 + x[2] - x[3] + x[1] - x[4]; g[2] = x[2]^2 + 2.0*x[3]^2 + x[1]^2 + 2.0*x[4]^2 - x[2] - x[4]; g[3] = 2.0*x[2]^2 + x[3]^2 + x[1]^2 + 2.0*x[2] - x[3] - x[4]; return(list(F=F, g=g)); } ################################################### ### code chunk number 14: ex4f ################################################### #========================= PROBLEM SPECIFICATIONS =========================== problem<-list(f="ex4", x_L=rep(0,4), x_U=rep(10,4), x_0=c(3,4,5,1), int_var=3, c_L=rep(-Inf,3), c_U=c(8,10,5)) opts<-list(maxtime=2) #========================= END OF PROBLEM SPECIFICATIONS ===================== Results<-MEIGO(problem,opts,algorithm="ESS"); ################################################### ### code chunk number 15: ex5f ################################################### ex5<-function(x,texp,yexp){ yini<-c(100, 0, 0, 0 ,0); times<-texp; out <- lsodes(yini, times, ex5_dynamics, parms = x) tout<-out[,1]; yout<-out[,-1]; J <- sum((yout-yexp)^2); g<-0; residuals<-(yout-yexp); return(list(J,g,residuals)) } #*************************************************** #Function of the dynamic system ex5_dynamics<-function(t,y,p){ dy<-rep(0,5) dy[1]<--(p[1]+p[2])*y[1]; dy[2]<-p[1]*y[1]; dy[3]<-p[2]*y[1]-(p[3]+p[4])*y[3]+p[5]*y[5]; dy[4]<-p[3]*y[3]; dy[5]<-p[4]*y[3]-p[5]*y[5]; return(list(dy)) } #*************************************************** ################################################### ### code chunk number 16: ex5Solve ################################################### #======================= #PROBLEM SPECIFICATIONS #======================= problem<-list(f="ex5", x_L=rep(0.01,5), x_U=rep(1,5), x_0=rep(0.1,5)) opts<-list(maxeval=1000, log_var=1:5, local_solver='NL2SOL', inter_save=1) #======================== #END OF PROBLEM SPECIFICATIONS #======================= #time intervals texp<-c(0, 1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420) #Distribution of species concentration # y(1) y(2) y(3) y(4) y(5) yexp<-rbind(c(100.0, 0.0, 0.0, 0.0, 0.0), c(88.35, 7.3, 2.3, 0.4, 1.75), c( 76.4 , 15.6, 4.5 , 0.7 , 2.8), c(65.1 , 23.1 , 5.3 , 1.1 , 5.8), c(50.4 , 32.9 , 6.0 , 1.5 , 9.3), c(37.5 , 42.7 , 6.0 , 1.9 , 12.0), c(25.9 , 49.1 , 5.9 , 2.2 , 17.0), c( 14.0 , 57.4 , 5.1 , 2.6 , 21.0), c(4.5 , 63.1 , 3.8 , 2.9 , 25.7)); Results<-MEIGO(problem,opts,algorithm="ESS",texp,yexp); ################################################### ### code chunk number 17: ex3solve ################################################### #========================= #PROBLEM SPECIFICATIONS #========================= problem<-list(f="ex3",x_L=rep(0,6),x_U=c(rep(1,4),16,16), neq=4, c_L=-Inf, c_U=4) opts<-list(ndiverse=2, local_solver='SOLNP', local_tol=3) #========================= #END OF PROBLEM SPECIFICATIONS #========================= k1=0.09755988; k3=0.0391908; k2=0.99*k1; k4=0.9*k3; Results<-MEIGO(problem,opts,algorithm="multistart",k1,k2,k3,k4); ################################################### ### code chunk number 18: sourceRosen10 ################################################### rosen10<-function(x){ f<-0; n=length(x); for (i in 1:(n-1)){ f <- f + 100*(x[i]^2 - x[i+1])^2 + (x[i]-1)^2; } return(f) } nvar<-10; ################################################### ### code chunk number 19: rosen10settings ################################################### problem<-list(f="rosen10", x_L=rep(-5,nvar), x_U=rep(1,nvar)) ################################################### ### code chunk number 20: algVNS ################################################### algorithm="VNS"; ################################################### ### code chunk number 21: vnsOPts ################################################### opts<-list(maxeval=2000, maxtime=3600*69, use_local=1, aggr=0, local_search_type=1, decomp=1, maxdist=0.5) ################################################### ### code chunk number 22: vnsCallMEIGO ################################################### Results<-MEIGO(problem,opts,algorithm); ################################################### ### code chunk number 23: options (eval = FALSE) ################################################### ## ## opts=list(); ## opts[[1]]$maxeval<-value1 ## opts[[2]]$maxeval<-value2 ## opts[[3]]$maxeval<-value3 ################################################### ### code chunk number 24: ceSSexample ################################################### rosen10<-function(x){ f<-0; n=length(x); for (i in 1:(n-1)){ f <- f + 100*(x[i]^2 - x[i+1])^2 + (x[i]-1)^2; } return(f) } nvar=20; problem<-list(f=rosen10, x_L=rep(-1000,nvar), x_U=rep(1000,nvar)); #Set 1 nodes and 2 cpu's per node n_nodes=1; n_cpus_per_node=3; #Set different values for dim_refset, bal and n2 #for each of the 10 cpu's to be used dim1 = 23; bal1 = 0; n2_1 = 0; dim2 = 33; bal2 = 0; n2_2 = 0; dim3 = 46; bal3 = 0; n2_3 = 2; dim4 = 56; bal4 = 0; n2_4 = 4; dim5 = 72; bal5 = 0.25; n2_5 = 7; dim6 = 72; bal6 = 0.25; n2_6 = 10; dim7 = 88; bal7 = 0.25; n2_7 = 15; dim8 = 101; bal8 = 0.5; n2_8 = 20; dim9 = 111; bal9 = 0.25; n2_9 = 50; dim10 = 123; bal10 = 0.25; n2_10 = 100; opts_dim=c(dim1,dim2,dim3,dim4,dim5,dim6,dim7,dim8,dim9,dim10); opts_bal=c(bal1,bal2,bal3,bal4,bal5,bal6,bal7,bal8,bal9,bal10); opts_n2=c(n2_1,n2_2,n2_3,n2_4,n2_5,n2_6,n2_7,n2_8,n2_9,n2_10); D=10; #Initialize counter and options counter=0; opts=list(); hosts=c(); for(i in 1:n_nodes){ for(j in 1:n_cpus_per_node){ counter=counter+1; #Set the name of every thread if(i<10)hosts=c(hosts,paste('node0',i,sep="")); if(i>=10 && i<100)hosts=c(hosts,paste('node',i,sep="")); opts[[counter]]=list(); #Set specific options for each thread opts[[counter]]$local_balance = opts_bal[counter]; opts[[counter]]$dim_refset = opts_dim[counter]; opts[[counter]]$local_n2 = opts_n2[counter]; #Set common options for each thread opts[[counter]]$maxeval = 10000; opts[[counter]]$local_solver = "dhc"; #Options not set will take default values for every thread } } #Set the address of each machine, defined inside the 'for' loop opts$hosts=c('localhost','localhost','localhost'); #Do not define the additional options for cooperative methods (e.g., ce_maxtime, ce_isparallel, etc..) #They will take their default values opts$ce_niter=5; opts$ce_type="SOCKS"; opts$ce_isparallel=TRUE; #Call the solver Results<-MEIGO(problem,opts,algorithm="CeSSR") ################################################### ### code chunk number 25: ceVNSexample ################################################### rosen10<-function(x){ f<-0; n=length(x); for (i in 1:(n-1)){ f <- f + 100*(x[i]^2 - x[i+1])^2 + (x[i]-1)^2; } return(f) } nvar=20; problem<-list(f=rosen10, x_L=rep(-1000,nvar), x_U=rep(1000,nvar)) opts=list(); opts[[1]]=list(use_local=1,aggr=1,local_search=1,decomp=1,maxdist=0.8,maxeval=2000); opts[[2]]=list(use_local=1,aggr=0,local_search=2,decomp=0,maxdist=0.5,maxeval=2000); opts[[3]]=list(use_local=1,aggr=0,local_search=2,decomp=0,maxdist=0.5,maxeval=2000); opts[[4]]=list(use_local=1,aggr=0,local_search=2,decomp=0,maxdist=0.5,maxeval=2000); opts$hosts=c('localhost','localhost','localhost','localhost'); opts$ce_niter=10; opts$ce_type="SOCKS"; opts$ce_isparallel= TRUE; Results=MEIGO(problem,opts, algorithm="CeVNSR"); ################################################### ### code chunk number 26: CNORodeExample_fig (eval = FALSE) ################################################### ## library(MEIGOR); ## library(CNORode); ## data("CNORodeExample",package="MEIGOR"); ## plotCNOlist(cnolist); ################################################### ### code chunk number 27: CNORodeExample_fig_plot ################################################### library(MEIGOR); library(CNORode); data("CNORodeExample",package="MEIGOR"); plotCNOlist(cnolist); ################################################### ### code chunk number 28: CNORodeExample_fig_model (eval = FALSE) ################################################### ## plotModel(model, cnolist); ################################################### ### code chunk number 29: CNORodeExample_fig_model_plot ################################################### plotModel(model, cnolist); ################################################### ### code chunk number 30: CNORodeExample_initial (eval = FALSE) ################################################### ## ## initial_pars=createLBodeContPars(model, LB_n = 1, LB_k = 0.09, ## LB_tau = 0.1, UB_n = 5, UB_k = 0.95, UB_tau = 10, random = TRUE); ## ## simData=plotLBodeFitness(cnolist, model,initial_pars, ## reltol = 1e-05, atol = 1e-03, maxStepSize = 0.01); ################################################### ### code chunk number 31: CNORodeExample_initial_plot ################################################### initial_pars=createLBodeContPars(model, LB_n = 1, LB_k = 0.09, LB_tau = 0.1, UB_n = 5, UB_k = 0.95, UB_tau = 10, random = TRUE); simData=plotLBodeFitness(cnolist, model,initial_pars, reltol = 1e-05, atol = 1e-03, maxStepSize = 0.01); ################################################### ### code chunk number 32: CNORodeExample_optimize ################################################### f_hepato<-getLBodeContObjFunction(cnolist, model, initial_pars, indices=NULL, time = 1, verbose = 0, transfer_function = 3, reltol = 1e-04, atol = 1e-03, maxStepSize =Inf, maxNumSteps = 1e4, maxErrTestsFails = 50, nan_fac = 5) ################################################### ### code chunk number 33: CNORodeExample_define_prob ################################################### n_pars=length(initial_pars$LB); problem<-list(f=f_hepato, x_L=initial_pars$LB[initial_pars$index_opt_pars], x_U=initial_pars$UB[initial_pars$index_opt_pars],x_0=initial_pars$LB[initial_pars$index_opt_pars]); ################################################### ### code chunk number 34: CNORodeExample_define_opts ################################################### opts<-list(maxeval=100, local_solver=0,ndiverse=10,dim_refset=6); ################################################### ### code chunk number 35: CNORodeExample_2 ################################################### Results<-MEIGO(problem,opts,algorithm="ESS") ################################################### ### code chunk number 36: CNORodeExample_opt (eval = FALSE) ################################################### ## opt_pars=initial_pars; ## opt_pars$parValues=Results$xbest; ## simData=plotLBodeFitness(cnolist, model,opt_pars, ## reltol = 1e-05, atol = 1e-03, maxStepSize = 0.01); ################################################### ### code chunk number 37: CNORodeExample_opt_plot ################################################### opt_pars=initial_pars; opt_pars$parValues=Results$xbest; simData=plotLBodeFitness(cnolist, model,opt_pars, reltol = 1e-05, atol = 1e-03, maxStepSize = 0.01); ################################################### ### code chunk number 38: CellNOptR_example_cnolist_fig (eval = FALSE) ################################################### ## library(MEIGOR); ## library(CellNOptR); ## data("CellNOptR_example",package="MEIGOR"); ## cnolist=CNOlist(cnolist); ## plotCNOlist(cnolist) ################################################### ### code chunk number 39: CellNOptR_example_cnolist ################################################### library(MEIGOR); library(CellNOptR); data("CellNOptR_example",package="MEIGOR"); cnolist=CNOlist(cnolist); plotCNOlist(cnolist); ################################################### ### code chunk number 40: CellNOptR_example_fig_model (eval = FALSE) ################################################### ## model <- preprocessing(cnolist, model, ## expansion=TRUE, compression=TRUE, verbose=FALSE); ## plotModel(model, cnolist); ################################################### ### code chunk number 41: CellNOptR_example_fig_model_plot ################################################### model <- preprocessing(cnolist, model, expansion=TRUE, compression=TRUE, verbose=FALSE); plotModel(model, cnolist); ################################################### ### code chunk number 42: CellNOptR_prep_fig_model ################################################### get_fobj<- function(cnolist,model){ f<-function(x,model1=model,cnolist1=cnolist){ simlist=prep4sim(model1) score=computeScoreT1(cnolist1, model1, x) return(score) } } fobj=get_fobj(cnolist,model) ################################################### ### code chunk number 43: CellNOptR_prep_3 ################################################### nvar=16; problem<-list(f=fobj, x_L=rep(0,nvar), x_U=rep(1,nvar)) opts<-list(maxeval=2000, maxtime=30, use_local=1, aggr=0, local_search_type=1, decomp=1, maxdist=0.5) ################################################### ### code chunk number 44: CellNOptR_4 ################################################### Results<-MEIGO(problem,opts,"VNS") ################################################### ### code chunk number 45: CellNOptR_example_plot_optModel (eval = FALSE) ################################################### ## optModel=cutModel(model,Results$xbest); ## plotModel(optModel,cnolist); ################################################### ### code chunk number 46: CellNOptR_example_plot_optModel_fig ################################################### optModel=cutModel(model,Results$xbest); plotModel(optModel,cnolist);