-
Notifications
You must be signed in to change notification settings - Fork 3
/
DEB_IBM.R
1592 lines (1325 loc) · 76.3 KB
/
DEB_IBM.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
## DEB IBM
# check bottom of page for diagnostics for running netlogo in rstudio
# version 1.6 (19-12-19)
# biocontrol, biocontrol with productivity
# DEB_INF_GUTS_IBM_1.1.nlogo with user-defined create-snails
# IndividualModel_IBM3.c (generates .o and .so files)
# FullStarve_shrink_dilute_damage3.Rd
# version 1.5 (17-8-19)
# new yRP, new DEB params from starvation exp, nocontrol event
# DEB_INF_GUTS_IBM_1.1.nlogo with user-defined create-snails
# IndividualModel_IBM3.c (generates .o and .so files)
# FullStarve_shrink_dilute_damage3.Rd
# version 1.4 (be event, me event, detr event)
# DEB_INF_GUTS_IBM_1.1.nlogo with user-defined create-snails
# IndividualModel_IBM3.c (generates .o and .so files)
# FullStarve_shrink_dilute_damage3.Rd
# version 1.3 (post sicb ver with molluscicide)
# DEB_INF_GUTS_IBM_1.1.nlogo
# IndividualModel_IBM3.c (generates .o and .so files)
# FullStarve_shrink_dilute_damage3.Rda
# version 1.2 (25-1-19)
# DEB_INF_GUTS_IBM_1.1.nlogo
# IndividualModel_IBM3.c (generates .o and .so files)
# FullStarve_shrink_dilute_damage3.Rda
# version 1.1
# DEB_INF_GUTS_IBM_1.1.nlogo
# IndividualModel_IBM2.c
# ILL_shrink_damageA5.Rda
# ver updates -------------------------------------------------------------
# 4-2-20
# set rg and rg_pars to 0 in detr sims for productivity (biocontrol)
# 2-2-20
# added rg0/ folder (biocontrol)
# 13-1-20
# added pred_ps = c(0,0.01,0.05,seq(0.1,1.5,0.1),2,2.5,3,3.5,4,4.5,5,10,15)
# added lower and upper snack size limits for heatmap (biocontrol)
# added highpred/ folder (biocontrol)
# 19-12-19
# removed old biocontrol toggles (be_event == 1)
# changed file handle from rg0.75_rep1.R to prod0.75_rep1.R
# 16-12-19
# added productivity to biocontrol (productivity == 1 and productivity fh)
# 9-12-19
# removed sicb fig2 and fig 3 user inputs and cleaned up user input
# 18-11-19
# added exponential predation rate (snail_snack_window = 0) (biocontrol)
# 23-10-19
# changed pred_a to 0.1 (biocontrol)
# added snail size dependent attack rate by predators (L3) for (biocontrol)
# 17-10-19
# changed pred_a to 0.6
# added lower pred_ps values to param range
# 17-9-19
# added snail_snack_min_vec and snail_snack_max_vec (snail_snack_window = 1)
# 16-9-19
# changed pred_ps > 1
# converted results in global_output to equal lengths
# 9-9-19
# set pred_p vector
# changed hb output from single value to vector for MoreArgs
# 6-9-19
# - changed hb to hb = as.numeric(pars["hb"]) + (pred_a * pred_p) / (1 + pred_a * pred_h * N.snails)*(snail.stats$L >= snail_snack_min & snail.stats$L <= snail_snack_max)
# - added snail_snack_min and snail_snack_max
# - added biocontrol = 0 and turned hb into vector
# - removed snail_snacks
# - added hb to snail.update matrix
# - removed hb from MoreArgs
# 3-9-19
# fixed ggplot gspl (error with color levels)
# added bio_list for saving biocontrol outputs
# added be_event 5 (non-selective feeding) to biocontrol
# 16-8-19 (v. 1.5)
# improved user input for M events
# added hailmary fh
# add nocontrol toggle and fh
# improved snailcontrol toggles
# 19-7-19
# added hailmary toggle
# 15-7-19
# new deb starvation params (samps_new)
# 10-7-19
# cleaned user input
# added pred vars for biocontrol sims
# 27-6-19
# created individual me_days for day30, monthly, bimonthly, skip30, 60, 90, 120, and hailmary events
# writes above outputs to dir
# 12-6-19
# added biocontrol (be_events)
# added init_host_pop to NL model v. 1.1
# added init_host_pop and init_host_pop_vec
# 6-6-19
# added rep loop in model input for me_days and me_im
# 31-5-19
# updated model input for detr_impact and me_days
# changed detritus impact to 50%
# final nih sims
# 30-5-19
# detr doesn't rebound after detr_impact
# 29-5-19
# removed Env_G[day] <- max(0, sum(Eggs),na.rm=T) after killing off eggs in env from me event
# added me_days for choosing independent me_days (multiple me events per sim)
# added detr_impact and detr_impact_days
# 25-5-19 (v. 1.4)
# molluscicide event and impact
# set hb to hb_pars[1] (0.001) on non snail control days
# added loop to automate save of sim outputs/plots
# 24-5-19
# added infected host length, summed host eggs, and mean host eggs
# 14-5-19
# mollusciciding kills eggs for all previous 10 days based on me_im [prob = exp(-me_im)]
# fixed empty agent set error with Env_G[day] <- 0 for no hosts (error: Error in NLCommand("ask patch 0 0 [set F", Env_F, "set M", Env_M, "set Z", )
# turned NULLs into NAs to get numeric values for hl_list and pmass_list (14-5-19. error: cannot coerce double)
# 13-5-19
# initial me and me_im sims
# 2-5-19
# added rep number for sicb fig 3
# made sim input section user-defined
# 1-5-19
# fixed create-snails error for me days
# 22-4-19 (v. 1.3)
# NL agentset error (v. 1.3)
## if snails > 0
## also kill eggs with me event
## if no snails
## no miracidia loss due to infection
## No cerceriae shedding
## no snail eggs
## food grows without feeding
# added molluscicide impact (me_im and me_im_pars)
# added legend to plots (sim_type). need to fix host biomass output
# 18-4-19 (v. 1.2)
# fixed snail control
# 15-4-19
# added chi param and host biomass, mean host length, and summed parasite mass to master output
# 9-4-19
# added molluscicide events
# changed yRP to 0.824 to match 2 hour hangover period
# 2-4-19
# changed hb pars to match sicb paper
# added file name (global_output_fh) to before sim
# 25-3-19
# added 0.04 to hb param range
# 5-3-19
# updated hb rates for sicb paper
# 11-2-19
# added harwell script
# 4-2-19
# added multi-panel plots for sim outputs
# 25-1-19 (v.1.2)
# MCMC DEB params (full starve model Rda)
# new detritus input
# new food growth in C file
# new mcmc deb estimates
# 23-1-19
# fixed individual snail deb update when using detritus: Food=environment[1]*(snail.update[,2]^2)/sum(snail.update[,2]^2),
# added detritus supply as food source
# added detritus param: pars["d_Z"] = 0.2 # detritus mg L-1 day-1
# removed resource cycle toggle from within sim loop to define param section (now alpha and rho of 0 = no cycle)
# 22-1-19
# added host length and parasite biomass to sim outputs
# 15-1-19
# set LHS parameter space
# 31-12-18
# added molluscide events (me_pars) for 95% host adult and eggs in env mortality rate (per day)
# 21-12-18
# added adult, juv, and infected host pop that sheds to sim results output
# fixed install packages section for windows
# 20-12-18
# added host length and parasite biomass to model outputs
# removed .so .o and .dll files from github and added to .gitignore and sensitive files dir
#18-12-18
# changed p to rho
# added rgrowth (pars rg) (rg)
# added master lists for total and infected hosts
# changed resource wave eq to r = r + alpha * r * sin(2 * pi * t/rho)
#16-12-18
# new damage density deb params (IndividualModel_IBM2.c, ILL_shrink_damageA5.Rda)
# added alpha and periodicity (p) param space to resource dynamics
# 11-12-18
# changed overdamped and periodic food dynamics to cyclical resource dynamics
# fixed cyclical resource dynamics
# 28-11-18
# added overdamped and periodic food dynamics to nl loop
# 27-11-18
# fixed NAs in 'create snails' command (Env_G param)
# list to check if 'create snails' command isn't producing NAs from Env_G
# set pop density outputs in NL loop to integer to pass into Env_G and rbinom func
# 23-11-18
# all user inputs at beginning of doc
# 22-11-19
# added debfunction.txt and pars.txt for defining params
# 19-11-18
# added "DEB_INF_GUTS_IBM_1.1.nlogo" as test model
###### TO DO ######
# for new snail deb feeding
## in new C file
### redefine dFdf (dfdt = -iM * f * sum(L^2) + rF(1-K/F) + Det)
# LHC sampling of r, alpha, and rho for different host sizes (for cyclical food)
# get ratio of infected hosts and shedding infected hosts
# change r and food
# - plot food against r as heatmap
# define
# - starvation-related hazard rate
# - pars and debfunction
# find papers on periodicity in resource loads in pops (Nisbet, Gurney, daphnia, etc)
# E.g. Nisbet Gurney Population dynamics in a periodically varying environment
# verify what volume of food density is reasonable (F)
# - sin wave of resource change
# - step function of resources (on/off season) (most unrealistic)
# - non-regenerative detritus (event-based)
# R = algae supply rate, don't vary K
# heat map of where peaks or resources and peaks of cercariae occur
# fix days parameter in NL
# output NL plots to R
# OUTPUTS
# survival and shell length results from DEBstep
# plot of rP vs. P/V (parasite biomass outcome)
#############################################################################################
#################################### Mac OSX test ###########################################
#############################################################################################
# Search "@netlogo" for netlogo code in file
### Files required
# "DEB_IBM.R"
# "DEB_INF_GUTS_IBM.nlogo"
# "Starvation_parameters.Rda"
# "FullStarve_shrink_dilute_damage3.Rda"
# "IndividualModel_IBM3.c"
# "IndividualModel_IBM3.so" # Mac OSX. generated from C
# "IndividualModel_IBM3.o" # Mac OSX. generated from C
# "IndividualModel_IBM.dll" # Windows. generated from C
test.java <- 0 # 1 = run java diagnostics
# run java test
install.packages("RCurl")
require(RCurl)
if(test.java==1){
require(RCurl)
script <- getURL("https://raw.githubusercontent.com/darwinanddavis/SchistoIBM/master/mac/java_test.R", ssl.verifypeer = FALSE)
eval(parse(text = script))
# check rJava version
.jinit()
.jcall("java/lang/System", "S", "getProperty", "java.runtime.version")
# get latest Java/Oracle version: https://www.oracle.com/technetwork/java/javase/downloads/index-jsp-138363.html
}
# :three: [GCC compiler in R (unconfirmed)](https://stackoverflow.com/questions/1616983/building-r-packages-using-alternate-gcc)
# [Running Netlogo 6.0.+](https://github.com/NetLogo/NetLogo/issues/1282)
################################# Running NetLogo in Mac ##################################
# if using Mac OSX El Capitan+ and not already in JQR, download and open JGR
mac <- 0
if(mac==1){
install.packages('JGR',,'http://www.rforge.net/')
library(JGR)
JGR::JGR()
}
#############################################################################################
############################# Windows or JGR onwards ########################################
#############################################################################################
#################################### set user inputs #######################################
# isolate sensitive data:
# "FullStarve_shrink_dilute_damage3.Rda"
# set user outputs
snab <- 1 # 1 = use remote access (snab comp), 0 = run model on your comp
mac <- 0 # mac or windows system? 1 = mac, 0 = windows
gui <- 0 # display the gui? 1 = yes, 0 = no
old_install <- 0 # 1 = use non pacman package install method
pck <- 1 # if not already, install rnetlogo and rjava from source? 1 = yes, 0 = already installed
save_to_file <- 0 # 1 = save simulation outputs to local dir, 0 = plot in current R session
mcmcplot <- 0 # 1 = save mcmc plots to dir
traceplot <- 0 # 1 = include traceplot in mcmc plots? intensive!!!
options(scipen=10) # set exponent to whole numbers
# set dir paths (for "/" for both Windows and Mac)
if(snab==1){
# set dir paths (for "/" for both Windows and Mac)
wd <- "R:/CivitelloLab/matt/schisto_ibm" # set working directory
ver_nl <-"6.0.4"# type in Netlogo version. found in local dir.
ver_gcc <-"4.6.3" # NULL # type in gcc version (if known). leave as "NULL" if unknown
nl.path <- "C:/Program Files" # set path to Netlogo program location
}else{
wd <- "/Users/malishev/Documents/Emory/research/schisto_ibm/SchistoIBM" # set working directory
ver_nl <-"6.0.4" # type in Netlogo version. found in local dir.
ver_gcc <-"4.6.3" # NULL # type in gcc version (if known). leave as "NULL" if unknown
nl.path <- "/Users/malishev/Documents/Melbourne Uni/Programs" # set path to Netlogo program location
}
# define starting conditions for simulation model @netlogo
n.ticks <- 150 # set number of days to simulate
day <- 1 # number of days to run simulation
resources <- "event" # set resources: "cyclical" or "event"
resource_type <- "algae" # set resource type as "algae" or "detritus"
#################################### set model paths #######################################
# load files
deb_samps <- "FullStarve_shrink_dilute_damage3.Rda" # old param estimates (15-7-19)
deb_samps_new <-"Starvation_parameters.Rda"
deb_compile <- "IndividualModel_IBM3"
setwd(wd)
nl.model <- list.files(pattern="*1.1.nlogo") ;nl.model # Netlogo model
if(mac==1){
nl.path <- paste0(nl.path,"/NetLogo ",ver_nl,"/Java/"); cat("Mac path:",nl.path)
}else{
nl.path <- paste0(nl.path,"/NetLogo ",ver_nl,"/app/"); cat("Windows path:",nl.path)
}
model.path <- paste0(wd,"/"); model.path # set path to Netlogo model
#################################### load packages #######################################
# new installation method (11-9-19)
# try pacman method again without concurrent r session running (15-10-19)
install.packages("pacman")
require(pacman)
packages <- c("Matrix","deSolve","mvtnorm","LaplacesDemon","coda","adaptMCMC","sp","RNetLogo","ggplot2","RCurl","RColorBrewer","Interpol.T","lubridate","ggExtra","tidyr","ggthemes","reshape2","pse","sensitivity","beepr")
pacman::p_load(Matrix,deSolve,mvtnorm,LaplacesDemon,coda,adaptMCMC,sp,RNetLogo,ggplot2,RCurl,RColorBrewer,Interpol.T,lubridate,ggExtra,tidyr,ggthemes,reshape2,pse,sensitivity,beepr)
p_load(packages)
# install dev versions of rjava and rnetlogo
install.packages("rJava", repos = "https://cran.r-project.org/", type="source"); library(rJava)
install.packages("RNetLogo", repos = "https://cran.r-project.org/", type="source"); library(RNetLogo)
require(rJava,RNetLogo)
#TS for RCurl package on Windows
remove.packages("httr")
remove.packages("rentrez")
remove.packages("curl")
install.packages("RCurl")
library(RCurl)
ppp <- lapply(packages,require,character.only=T)
names(ppp) <- packages
if(any(ppp==F)){cbind(packages,ppp);cat("\n\n\n ---> Check packages are loaded properly <--- \n\n\n")}
# old installation method --------------------------------------------------------
if(old_install==1){
# if already loaded, uninstall RNetlogo and rJava
if(pck==1){
p<-c("rJava", "RNetLogo"); remove.packages(p)
# then install rJava and RNetLogo from source
if(mac==1){
install.packages("rJava", repos = "https://cran.r-project.org/", type="source"); library(rJava)
install.packages("RNetLogo", repos = "https://cran.r-project.org/", type="source"); library(RNetLogo)
require(rJava,RNetLogo)
# check pck versions
installed.packages()["RNetLogo","Version"]
installed.packages()["rJava","Version"]
}
# install relevant packages
packages <- c("Matrix","deSolve","mvtnorm","LaplacesDemon","coda","adaptMCMC","sp","RNetLogo","ggplot2","RCurl","RColorBrewer","Interpol.T","lubridate","ggExtra","tidyr","ggthemes","reshape2","pse","sensitivity","beepr")
if(require(packages)){
install.packages(packages,dependencies = T)
}
# load annoying packages manually because they're stubborn
install.packages("RNetLogo")
install.packages("RCurl")
install.packages("Interpol.T")
install.packages("lubridate")
install.packages("tidyr")
install.packages("ggthemes")
install.packages("ggExtra")
install.packages("beepr")
install.packages("Rcpp")
install.packages("rlang")
}
ppp <- lapply(packages,require,character.only=T)
names(ppp) <- packages
if(any(ppp==F)){cbind(packages,ppp);cat("\n\n\n ---> Check packages are loaded properly <--- \n\n\n")}
} # old install method
cs <- list() # diagnostics list for checking NAs in create snails command
# load plot function
script <- getURL("https://raw.githubusercontent.com/darwinanddavis/plot_it/master/plot_it.R", ssl.verifypeer = FALSE)
eval(parse(text = script))
display.brewer.all()
# Set global plotting parameters
cat("plot_it( \n0 for presentation, 1 for manuscript, \nset colour for background, \nset colour palette 1. use 'display.brewer.all()', \nset colour palette 2. use 'display.brewer.all()', \nset alpha for colour transperancy, \nset font style \n)")
plot_it(0,"blue","YlOrRd","Greens",1,"mono") # set plot function params
plot_it_gg("white","black") # same as above for ggplot
# load harwell script
script <- getURL("https://raw.githubusercontent.com/darwinanddavis/harwell/master/harwell.R", ssl.verifypeer = FALSE)
eval(parse(text = script))
################################ compile packages and load files ###################################
### Install rtools and gcc for using C code and coda package
#### https://cran.r-project.org/bin/macosx/tools/
# define paths for gcc compiler
if(mac==1){ #### Mac OSX
rtools <- "/usr/local/clang6/bin"
gcc <- paste0("usr/local/clang6/gcc-",ver_gcc,"/bin")
# Mac OSX
}else{ #### Windows
rtools <- "C:\\Rtools\\bin"
gcc <- paste0("C:\\Rtools\\gcc-",ver_gcc,"\\bin")
}
#### point to path on comp to access rtools and gcc for C compiler
path <- strsplit(Sys.getenv("PATH"), ";")[[1]]
new_path <- c(rtools, gcc, path)
new_path <- new_path[!duplicated(tolower(new_path))]
Sys.setenv(PATH = paste(new_path, collapse = ";"))
if(mac==1){
# dyn.unload("IndividualModel_IBM3.so") # unLoad .so (Mac OSX
system(paste0("R CMD SHLIB ",deb_compile,".c")) # generates .o and .so files
dyn.load(paste0(deb_compile,".so")) # Load .so (Mac OSX)
}else{
# compile model from C definition
# dyn.unload(paste0(deb_compile,".dll")) # unload dll (Windows only)
system(paste0("R CMD SHLIB ",deb_compile,".c"))
dyn.load(paste0(deb_compile,".dll"))# Load dll (Windows only)
}
#################################### load deb params #######################################
# load DEB starvation model parameters and create mcmc (and convert mcmc chain to coda format)
# old starvation param estimates (use gammaH, gammaP, and lpost to add to new samps)
samps = readRDS(deb_samps)
samps <- as.mcmc(samps[, c("iM", "k", "M", "EM", "Fh", "muD", "DR", "fe", "yRP",
"ph", "yPE", "iPM", "eh", "mP", "alpha", "yEF", "LM",
"kd", "z", "kk", "hb", "theta", "mR", "yVE", "yEF2",
"sd.LI1", "sd.LU1", "sd.EI1", "sd.EU1", "sd.W1", "sd.LI2",
"sd.LU2", "sd.EI2", "sd.EU2", "sd.W2", "gammaH", "gammaP", "lpost")])
samps_new = readRDS(deb_samps_new) # new starvation param estimates (15-7-19)
# ---------- summarise and plot estimated params
svar <- "M" # select variable
sampsvar <- samps[,svar] # pull from mcmc
summary(sampsvar) # get mean, sd, se, and quantiles for each input variable
den <- density(sampsvar) # get AUC
densplot(sampsvar, show.obs = F,type="n") # density estimate of each variable
polygon(den, col=adjustcolor(colv,alpha=0.5),border=colv) # fill AUC
plot(sampsvar,trace=T,density=T,col=colv) # traceplot (below) and density plot (above)
# intensive
traceplot(sampsvar,smooth=T,type="l",lwd=0.3,xlim=c(0,length(sampsvar)),col=colv[2],xlab=paste0("Iterations"),ylab=paste0("Sampled values"),main=paste0("Sampled values over iterations for ",svar)) # iterations vs sampled valued per variable
if(mcmcplot==1){
par(mfrow=c(1,1))
plotlist <- list()
pdf("mcmc_vars.pdf",onefile = T,paper="a4")
for(i in colnames(samps)){
par(bty="n", las = 1)
if(traceplot==1){
traceplot(sampsvar,smooth=T,type="l",xlim=c(0,length(sampsvar)),col=colv[2],xlab=paste0("Iterations"),ylab=paste0("Sampled values"),main=paste0("Sampled values over iterations for ",svar)) # iterations vs sampled valued per variable
}
svar <- i # select variable
sampsvar <- samps[,svar] # pull from mcmc
den <- density(sampsvar) # get AUC
densplot(sampsvar, show.obs = F,type="n",main=paste0("Density estimate of ",i)) # density estimate of each variable
polygon(den, col=adjustcolor(colv,alpha=0.5),border=colv) # fill AUC
}
dev.off()
} # end mcmcplot
# ----------
# get the best fit DEB parameters to match the data (using mcmc)
read.csv("pars.txt",header=T,sep="/",fill=T,flush=T,strip.white=T,row.names=NULL)
pars = as.vector(data.frame(samps)[max(which(data.frame(samps)$lpost >= max(data.frame(samps)$lpost) -0.001)),]) # old best parameter estimates using mcmc (pre 15-7-19)
pars <- c(samps_new,pars[c("gammaH", "gammaP","lpost")]) # add gammaH, gammaP, and lpost from original pars to new best parameter estimates
pars <- pars %>% as.data.frame()
pars["Fh"] = 2 # f_scaled (for v.1.1)
pars["ENV"] = 500 # Units: L
pars["r"] = 1 # Units: day-1
pars["step"] = 1 # Units: day
pars["epsilon"] = 20 # Units: L host-1, day-1 (Rounded estimate from Civitello and Rohr)
pars["sigma"] = 0.5
pars["m_M"] = 1 # Units: day-1
pars["m_Z"] = 1 # Units: day-1
pars["M_in"] = 10
pars["K"] = 5
pars["Det"] = 0.1 # Units mg C/L-1 d-1 (detritus)
pars["yRP"] = pars["yRP"]*17.5 # new yRP = 0.824 * 17.5 (1-8-19)( # old yRP = 0.0471 * (1/0.4) * 7 [expected shedding output per week * (snails shed 40% of their total cercariae during 9-11 AM.) * 7 days]
#################################### solve deb eqs #######################################
# display list of param definitions
read.csv("debfunction.txt",header=T,sep="/",fill=T,flush=T,strip.white=T,row.names=NULL)
DEB = function(step, Food, L, e, D, RH, P, RP, DAM, HAZ, iM, k, M, EM,
Fh, muD, DR, yRP, ph, yPE, iPM, eh, mP, alpha, yEF,
LM, kd, z, kk, hb, theta, mR, yVE, ENV, Lp, SAtotal, r, K, Det){
# starting conditions
initials = c(Food=Food, L=L, e=e, D=D, RH=RH, P=P, RP=RP, DAM=DAM, HAZ=HAZ)
# deb parameters
parameters = c(iM, k, M, EM, Fh, muD, DR, yRP, ph, yPE, iPM,
eh, mP, alpha, yEF, LM, kd, z, kk, hb, theta, mR, yVE, ENV, Lp, SAtotal, r, K, Det)
# estimate starting deb conditions using fitted params by solving ode's
## return survival and host shell length
DEBstep <- lsoda(initials, c(0,step), func = "derivs", dllname = deb_compile,
initfunc = "initmod", nout=2, outnames=c("Survival", "LG"), maxsteps=500000,
as.numeric(parameters), rtol=1e-6, atol=1e-6, hmax=1)
DEBstep[2, 2:12] # 12 = survival
} # end deb model
### deb output for each timestep
result = DEB(step=1, Food=5, L=10, e=0.9, D=as.numeric(pars["DR"]), RH=0, P=0, RP=0, DAM=0, HAZ=0, iM=pars["iM"], k=pars["k"], M=pars["M"], EM=pars["EM"],
Fh=pars["Fh"], muD=pars["muD"], DR=pars["DR"], yRP=pars["yRP"], ph=pars["ph"], yPE=pars["yPE"], iPM=pars["iPM"], eh=pars["eh"],
mP=pars["mP"], alpha=pars["alpha"], yEF=pars["yEF"], LM=pars["LM"], kd=pars["kd"], z=pars["z"], kk=pars["kk"], hb=pars["hb"],
theta=pars["theta"], mR=pars["mR"], yVE=pars["yVE"], ENV=pars["ENV"], Lp=10,SAtotal=7007.822, r=pars["r"], K=pars["K"], Det=pars["Det"])
### Exposure submodel
# pass the deb state vars into infection model
Infection = function(snail.stats, miracidia, parameters){
# Parameters
epsilon = as.numeric(parameters["epsilon"])
sigma = as.numeric(parameters["sigma"])
ENV = as.numeric(parameters["ENV"])
m_M = as.numeric(parameters["m_M"])
step = as.numeric(parameters["step"])
# Later calculations depend on exposure probabilities
exp.rates = epsilon/ENV*(snail.stats[,"L"]>0) # This is just to get uniform exposure rates
sum.exp.rates = sum(exp.rates)
# Probabilities for fate of miracidia
## Still in water
P.left.in.water = exp(-(m_M+sum(exp.rates))*step)
## Infect a snail
P.infects.this.snail = (1 - P.left.in.water)*(sigma*exp.rates/(m_M+sum.exp.rates))
## Die in water or fail to infect
P.dead = (1 - P.left.in.water)*(m_M/(m_M+sum.exp.rates)) + sum((1 - P.left.in.water)*((1-sigma)*exp.rates/(m_M+sum.exp.rates)))
prob.vector = c(P.infects.this.snail, P.left.in.water, P.dead)
# Multinomial outcome from number of miracidia in env based on their survival probability
rmultinom(n=1, size=miracidia, prob=prob.vector)
#sum(P.left.in.water, P.invades.this.snail, P.dead)
} # end infection model
### update all the snails @netlogo
update.snails = function(who, new.L, new.e, new.D, new.RH, new.P, new.RP, new.DAM, new.HAZ, new.LG){
paste("ask snail", who,
"[set L", new.L,
"set ee", new.e,
"set D", new.D,
"set RH", new.RH,
"set P", new.P,
"set RPP", new.RP,
"set DAM", new.DAM,
"set HAZ", new.HAZ,
"set LG", new.LG, # new max length
"]")
} # end host update
#Example update
#paste(mapply(update.snails, who=snail.stats[,"who"], new.L=L, new.e=e, new.D=D, new.RH=RH, new.P=P, new.RP=RP, new.DAM=DAM, new.HAZ=HAZ), collapse=" ")
geterrmessage() # check if there were any error messages
###########################################################################################
#################################### load netlogo ########################################
###########################################################################################
# @netlogo
# working NLStart in RStudio. works with gui=F (2018/09/24)
if(gui==0){
NLStart(nl.path,gui=F,nl.jarname = paste0("netlogo-",ver_nl,".jar")) # open netlogo without a gui
}else{
NLStart(nl.path,nl.jarname = paste0("netlogo-",ver_nl,".jar")) # open netlogo
}
NLLoadModel(paste0(model.path,nl.model),nl.obj=NULL) # load model
# if java.lang error persists on Mac, try copying all .jar files from the 'Java' folder where Netlogo is installed into the main Netlogo folder
################################################################################################
#################################### start netlogo sim ########################################
################################################################################################
# WHO sims
# every 30 days
# 0.99, 0.95, 0.90, 0.75, 0.5
# 3 reps
#grant proposal sims
# rg = 0.25, det = 0.1
# infected length, egg density, mean eggs
# me sims (14-5-19)
# season = 150 days
# every 2 weeks at different me_im (4 sims) # me_event 1
# every month at different me_im (4 sims) # me_event 2
# every 2 months at different me_im (4 sims) # me_event 3
# 0.25 for algae and detritus
# 4 reps each
################################# 14-15-19 sims
# ~~~~~~ user input --------------------------------------------------------------
# for hailmary
# snail_control <- 1
# hailmary <- 1
# for(me_im_event in 4){ ... }
# for(me_event in 8){ ... }
# enable for(me in me_pars){ ... } (line 613) and closing bracket
# disable me_pars loop in sim model and closing bracket
# enable three @hailmary instances in sim model
# for all snail control events
# snail_control <- 1
# hailmary <- 0
# for(me_im_event in 1:5){ ... }
# for(me_event in 1:7){ ... }
# disable for(me in me_pars){ ... } (line 613) and closing bracket
# enable me_pars loop in sim model and closing bracket
# disable three @hailmary instances in sim model
#run 13-1-20
# 1.2
snail_snack_min_vec <- c(0) # min size host to eat
snail_snack_max_vec <- c(2.5) # max size host to eat
# 0.6 and 1.2
snail_snack_min_vec <- c(0) # min size host to eat
snail_snack_max_vec <- c(seq(5,20,2.5),60) # max size host to eat
# detritus
# 0_50 and pred_ps = c(0.6,0.7,0.8,0.9,1.1,1.2,1.3,1.4)
# 19-12-19
# to run
# algae and det 5 reps 0-15 mm pred_p c(0,1,1.5,2,2.5,3,3.5,4,4.5,5,10,15)
# algae and det 5 reps 5+ mm
# pred_ps sims 10-01-2020
pred_ps <- c(0,0.01,0.05,seq(0.1,1.5,0.1),2,2.5,3,3.5,4,4.5,5,10,15)
# per L
pred_ps / 50
# per m^2
pred_ps / 50 * 500
# -------------------------------------------------------------------------
resource_type="detritus" # detritus # set resource type
# these can both be > 0 (both are toggled on/off in sim code)
# algae params
rg_pars <- c(0,0.01,0.05,0.1,0.2,0.25) # resource growth rates (r)
# rg_pars <- 0.25
# detritus params
detr_pars <- c(0,0.001,0.005,0.01,0.05,0.1,0.2,0.25) # detritus input (mg L^-1 day^-1)
# detr_pars <- 0.25
# mortality params
hb_pars <- 0.001
snail_control <- 0 # run molluscicide sims?
hailmary <- 0 # run hailmary sims separately to other molluscicide sims
detr_impact = 0 # run detritus impact?
biocontrol = 1 # run biocontrol?
snail_snack_window = 0 # 1 = run biocontrol for sliding window of snail size classes, 0 = exp predation rate
productivity <- 0 # 1 = run biocontrol with varying productivity levels (rg and detr)
no_control <- 0 # run normal no control sims
# snail control
me_events <- 151 # 1 = on day 30, every month, skip a month. 8 = hailmary
me_im_events <- 1 # 1 = 0.69, 2 = 1.39, 3 = 2.30, 4 = 3.0, 5 = 4.60 for me_im
me_im <- 0 # create me_im vec
# detritus control
detr_im <- c(0.1,0.25,0.5,0.9,0.99) # percent reduction in detritus for detritus impact
detr_impact_days = c(5,15,30,60) # days for each detritus impact sim
# biocontrol
init_host_pop = 50 # initial population size
init_host_pop_vec <- c(50)
# functional pred feeding response: fN = aNP / 1 + ahN
pred_a <- 0.1 # predator attack rate: 50L sweeped / day or 8.33 * 6 L tanks (from sokolow etal 2014 acta tropica)
pred_h <- 0.1 # pred handling time # 0.1 = 10 snails per day
# pred_p <- 0.05 # predator population density
# pred_ps <- c(seq(0.1,1.5,0.1),2,2.5,3)
# pred_ps <- c(0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,10,15)
pred_ps <- c(0,0.5,1,1.5,2,5,10,15,20,30,50,100)
fh_buff = 10 # buffer to convert pred_ps digits into integer for saving file handle
snail_snack_min_vec <- c(0) # min size host to eat
snail_snack_max_vec <- c(0) # max size host to eat
# . 0-5 mm
# . 0-10
# . 0-15
# . 5-20
# . 10-20
# . 15+
rep_num <- 5 # number of replications
# for productivity == 1
### need to comment out rg and detr in next loop
# for(rg in rg_pars){ # loop through rgs (food growth rates)
for(detr in detr_pars){ # loop through det (food growth rates)
# for(snail_snack_min in snail_snack_min_vec){
# for(snail_snack_max in snail_snack_max_vec){
# file handle for bio
ifelse(snail_snack_window==1, be_fh <- paste0(snail_snack_min,"_",snail_snack_max), be_fh <- "exp")
if(hailmary==1){me_pars <- seq(10,140,10); me_event <- 8}else{me_pars <- n.ticks + 1} # set hailmary file handle
# save multiple sims to dir ---------------------------------------
for(init_host_pop in init_host_pop_vec){
for(pred_p in pred_ps){
# for(detr_impact_days in detr_impact_days){
for(rn in 1:rep_num){
# for(detr_im in detr_im){
for(me_im_event in me_im_events){ # 1:5 # me impact (0.69, 1.39, 2.3 ...)
for(me_event in me_events){ # 1:7 # run all me event scenarios and save to file.
cat("\npredp=",pred_p)
cat("\nprod=",rg)
cat("\nrep=",rn)
# for hailmary, uncomment below and comment out me_pars loop in next loop
# for(me in me_pars){ # loop through mes (molluscicide events) for saving one me event per sim to dir (hailmary)
resources="event" # set resource cycles
alpha_pars = 0
rho_pars = 1
# set type of resource input @netlogo
set_resource_type<-function(resource_type){ # set resource input in env
if(resource_type == "detritus"){NLCommand("set resource_type \"detritus\" ")}else{NLCommand("set resource_type \"algae\" ")}}
set_resource_type(resource_type) # set resource type: "detritus" or "algae" @netlogo
# set type of resource dynamics @netlogo
set_resources<-function(resources){ # set resource input in env
if (resources == "cyclical"){NLCommand("set resources \"cyclical\" ")}else{NLCommand("set resources \"event\" ")}}
set_resources(resources) # set resources: "cyclical" or "event" @netlogo
cat("\nResource type = ",resource_type,"\nResources = ",resources)
# set initial pop size @netlogo
NLCommand("set init_host_pop", init_host_pop)
if(save_to_file==1){pdf(paste0(wd,"/master_sim.pdf"),onefile=T,paper="a4")}
if(resource_type=="detritus"){
detr_pars <- detr_pars; alpha_pars <- 0; rho_pars <- 1; rg <- 0; rg_pars <- 0;cat("\ndetritus input = ",detr_pars,"\nrg = ",rg_pars,"\nhb = ", hb_pars)
}else{detr <- 0; detr_pars <- 0;cat("detritus input = ", detr_pars,"\nrg = ",rg_pars)}
# set resource to cycle or be constant
if(resource_type=="algae"){
if(resources=="cyclical"){
rg_pars <- rg_pars # resource growth rates (rs)
alpha_pars <- c(0,0.25,0.5,0.75,1) # amplitude of resources (alphas)
rho_pars <- c(1,seq(10,n.ticks,10)) # periodicity of resources (rhos)
cat("\nalphas = ",alpha_pars,"\nrhos = ",rho_pars,"\nrgs = ",rg_pars)
}else{alpha_pars <- 0; rho_pars <- 1; rg_pars <- rg_pars;cat("\nalphas = ",alpha_pars,"\nrhos = ",rho_pars,"\nrg = ",rg_pars,"\nhb = ", hb_pars)}
}
# # define param sample space with LHS
# require(sp)
# require(pse)
# require(lhs)
# lhsmodel <- function(params){
# params <- factors_space[[2]]*factors_space[[3]]*factors_space[[4]]
# }
# factors <- c("alpha","rho","rg","me") # name of params
# factors_space <- list(alpha_pars,rho_pars,rg_pars,me_pars)
# q <- rep("qnorm",length(factors)) # apply the dist to be used
# q.arg <- list(list(alpha_pars),list(rho_pars),list(rg_pars),list(me_pars)) # inputs for dist q
# # list(list(mean=1.7, sd=0.3), list(mean=40, sd=1),list(min=1, max=50) )
# N <- prod(as.numeric(summary(factors_space)[,1]))
# lhs_model <- LHS(model=lhsmodel,factors=factors,N=N,q=q,q.arg=q.arg,nboot=100)
# lhs_data <- get.data(lhs_model) # param space from LHS
# lhs_results <- get.results(lhs_model)
# get.N(lhs_model) # get the number of output points in hypercube
# lhs_data
Env_G = numeric() # create empty environment vector
day <- 1 # reset sim days
# individual outputs
cerc_list <- list() # cercariae
food_list <- list() # food in env
juv_list <- list() # juvenile hosts
adult_list <- list() # adult hosts
infec_list <- list() # infected hosts
infec_shed_list <- list() # infected shedding hosts
hl_list <- list() # host length
pmass_list <- list() # parasite biomass
host_biomass_list <- list() # host biomass list
egg_list <- list() # egg density in env
egg_mean_list <- list() # mean host eggs list
infec_shed_length_list <- list() # length of infected shedding hosts
# master outputs
cerc_master <- list() # master list for cerc density (Env_Z)
food_master <- list() # master list for food dynamics (Env_F)
juv_master <- list() # master list for total host pop ()
adult_master <- list() # master list for total host pop ()
infec_master <- list() # master list for infected host pop ()
infec_shed_master <- list() # master list for infected shedding host pop
hl_master <- list() # master list for host length
pmass_master <- list() # master list for parasite biomass
host_biomass_master <- list() # master list for host biomass
egg_master <- list() # master list for egg density in env
egg_mean_master <- list() # master list for mean host eggs
infec_shed_length_master <- list() # master list for length of infected shedding hosts
# define plot window
plot.matrix <- matrix(c(length(rg_pars),length(detr_pars)))
par(mfrow=plot.matrix)
# ~~~~~~ snail control -----------------------------------------------------------
if(snail_control == 1){
me_im_pars <- c(0.69, 1.39, 2.3, 3, 4.6) # 2.3 = 90% snail mortality from molluscicide event (per day)
if(me_event==1){me_days = c(60,120); me_fh = "bimonthly"}
if(me_event==2){me_days = c(60,90,120); me_fh = "skip30"}
if(me_event==3){me_days = c(30,90,120); me_fh = "skip60"}
if(me_event==4){me_days = c(30,60,120); me_fh = "skip90"}
if(me_event==5){me_days = c(30,60,90); me_fh = "skip120"}
if(me_event==6){me_days = 30; me_fh = "day30"}
if(me_event==7){me_days = c(30,60,90,120); me_fh = "monthly"}
if(me_event==8){me_days = seq(10,140,10);me_fh = paste0("hailmary",me)}
if(me_event==151){me_days <- n.ticks + 1; me_fh = ""}
if(me_im_event==1){me_im_pars = me_im_pars[1]}
if(me_im_event==2){me_im_pars = me_im_pars[2]}
if(me_im_event==3){me_im_pars = me_im_pars[3]}
if(me_im_event==4){me_im_pars = me_im_pars[4]}
if(me_im_event==5){me_im_pars = me_im_pars[5]}
hb_pars <- 0.001
if(resource_type=="algae"){detr_pars <- 0; algae <- rg_pars}else{detr_pars <- detr_pars; rg_pars <- 0}
cat("\nalgae:",rg_pars,"\ndetritus:",detr_pars,"\nrho:",0,"\nalpha:",alpha_pars,"\nmortality (if not mollusciciding):",hb_pars,"\nmolluscicide days:",me_pars, "\nmolluscicide impact: ",me_im_pars)
}else{ # dont run snail control
me_pars <- n.ticks+1
me_days <- n.ticks + 1
me_im_pars <- 0
me_im <- 0
hb_pars <- 0.001#hb_pars
cat("\nSnail control will occur every ",max(me_pars)/length(me_pars)," days \n Mortality is ",hb_pars)
}
cat("Mollusciding on day", me_pars)
me_im_pars
# if(resource_type=="algae"){detr =0; rg = 0.25;detr_impact=0;detr_impact_days=n.ticks+1}else{detr=0.25;rg=0; alpha_pars=1; rho_pars = 1}
# @detr_impact
if(detr_impact==1){
detr_impact_days = detr_impact_days # set detritus impact days
}else{
detr_impact_days = n.ticks+1 # no detritus impact
}
# file handles ----------------------------------------------------------
if(snail_control==1){
if(hailmary==1){ # @hailmary
fhh = paste0(resource_type,"_",me_fh,"_rep",rn);fhh # use for success/failure plot (fig 2) in plos_one sim
# fhh = paste0(resource_type,"_",me_fh,"_meim",me_im_event,"_rep",rn);fhh # use for hailmary120 for me_im_event = 1:5. plosone fig 3 (28-7-19))
global_output_fh = paste0(wd,"/plos_sims/hailmary/",fhh,".R")
}else{
fhh = paste0(resource_type,"_",me_fh,"_meim",me_im_event,"_rep",rn);fhh
global_output_fh = paste0(wd,"/plos_sims/",me_fh,"/",fhh,".R")
}
}
if(detr_impact==1){
fhh = paste0("detday_",detr_impact_days[1],"_detim_",detr_im*1000);fhh
global_output_fh = paste0(wd,"/detr_impact_sims/",fhh,".R")
}
if(biocontrol==1){
if(productivity==1){ # for altering rg and det values
ifelse(resource_type=="algae", resource_fh <- rg, resource_fh <- detr)
fhh = paste0(resource_type,"_",be_fh,"_hostpop",init_host_pop,"_predpop",pred_p*fh_buff,"_prod",resource_fh,"_rep",rn);fhh
global_output_fh = paste0(wd,"/biocontrol_sims/",init_host_pop,"/new_new_pred_a/productivity/productivity_0_25/pred_05_100/rg0/",fhh,".R")
}else{
fhh = paste0(resource_type,"_",be_fh,"_hostpop",init_host_pop,"_predpop",pred_p*fh_buff,"_rep",rn);fhh
global_output_fh = paste0(wd,"/biocontrol_sims/",init_host_pop,"/new_new_pred_a/highpred/",fhh,".R")
}
}
if(no_control==1){
# no control scenario
fhh = paste0(resource_type,"_meim",me_im_event,"_rep",rn);fhh
global_output_fh = paste0(wd,"/plos_sims/nocontrol/",fhh,".R")
}
if(resource_type=="algae"){detr_pars <- 0; algae <- rg_pars}else{detr_pars <- detr_pars; rg_pars <- 0}
if(length(me_days)== 1 | 2 | 3){me_days = rep(me_days,4)} # set vector same length as in sim model
if(length(detr_impact_days)==1){detr_impact_days = rep(detr_impact_days,2)} # set vector same length as in sim model
cat("\nalgae:",rg_pars,"\ndetritus:",detr_pars,"\nrho:",0,"\nalpha:",alpha_pars,"\nmortality (if not mollusciciding):",hb_pars,"\nmolluscicide days:",unique(me_days), "\nmolluscicide impact: ",me_im_pars)
cat("\nM days = ",unique(me_days)); cat("\nM impact = ",me_im); cat("\nDet impact days = ",detr_impact_days)
global_output_fh
#################################### start netlogo sim ########################################
for(hb in hb_pars){
# for(detr in detr_pars){ # loop through detritus inputs
for(alpha in alpha_pars){ # loop through alphas (amplitude in food cycle)
for(rho in rho_pars){ # loop through rhos (periodicity of food cycle)
# for(rg in rg_pars){ # loop through rgs (food growth rates)
for(me in me_pars){ # loop through me (molluscicide day events) comment out for @hailmary
for(me_im in me_im_pars){ # loop through me_im (molluscicide impact events)
NLCommand("setup")
day <- 1 # reset days
Env_G[is.na(Env_G)] <- 0 # turn NAs to 0 to feed into rbinom function
for(t in 1:n.ticks){ # start nl sim @netlogo
# define food dynamics for cyclical algal (logistic food growth equation) or detritus food sources
alpha <- alpha # amplitude of resources
rho <- rho # periodicity (time range of resource cycles)
rg <- rg # resource growth rate
rg_t <- rg + alpha * rg * sin(2 * pi * t/rho) # equilibrium cyclical resource dynamics (19-12-18)
pars["r"] <- rg_t # set resource growth rate
pars["Det"] <- detr # Units mg C/L-1 d-1 (detritus)
# @detr_impact -------------------------------------------------------------------