Mercurial > repos > greg > bmsb

comparison bmsb.R @ 35:7c40c2b303f1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 16 Dec 2016 10:16:00 -0500
parents 390ed5192839
children d3a6b291e096
comparison
equal deleted inserted replaced
34:7c43bc0ce69c 35:7c40c2b303f1
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"), 6 make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"),
7 make_option(c("-b", "--adult_nymph_accum"), action="store", dest="adult_nymph_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"), 7 make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"),
8 make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"), 8 make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"),
9 make_option(c("-d", "--latitude"), action="store", dest="latitude", type="double", help="Latitude of selected location"), 9 make_option(c("-d", "--latitude"), action="store", dest="latitude", type="double", help="Latitude of selected location"),
10 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"), 10 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
11 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 11 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
12 make_option(c("-g", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 12 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
13 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"), 13 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),
14 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"), 14 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),
15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
204 mort.prob=temperature * 0.0005 + 0.02 204 mort.prob=temperature * 0.0005 + 0.02
205 } 205 }
206 return=mort.prob 206 return=mort.prob
207 return 207 return
208 } 208 }
209
210 cat("Replications: ", opt$replications, "\n")
211 cat("Photoperiod: ", opt$photoperiod, "\n")
212 cat("Oviposition rate: ", opt$oviposition, "\n")
213 cat("Egg mortality rate: ", opt$egg_mort, "\n")
214 cat("Nymph mortality rate: ", opt$nymph_mort, "\n")
215 cat("Adult mortality rate: ", opt$adult_mort, "\n")
216 cat("Min clutch size: ", opt$min_clutch_size, "\n")
217 cat("Max clutch size: ", opt$max_clutch_size, "\n")
218 cat("(egg->young nymph): ", opt$young_nymph_accum, "\n")
219 cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n")
220 cat("(old nymph->adult): ", opt$adult_accum, "\n")
209 221
210 # Read in the input temperature datafile 222 # Read in the input temperature datafile
211 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset) 223 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset)
212 224
213 # Initialize matrix for results from all replications 225 # Initialize matrix for results from all replications
617 # SE for adults 629 # SE for adults
618 sa.se <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications) 630 sa.se <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications)
619 # SE for nymphs 631 # SE for nymphs
620 sn.se <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd) 632 sn.se <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd)
621 # SE for eggs 633 # SE for eggs
622 se.se <- apply(S0.rep ,1 ,sd) / sqrt(opt$replications) 634 se.se <- apply(S0.rep, 1, sd) / sqrt(opt$replications)
623 # SE value for P 635 # SE value for P
624 g0.se <- apply(g0.rep, 1, sd) / sqrt(opt$replications) 636 g0.se <- apply(g0.rep, 1, sd) / sqrt(opt$replications)
625 # SE for F1 637 # SE for F1
626 g1.se <- apply(g1.rep, 1, sd) / sqrt(opt$replications) 638 g1.se <- apply(g1.rep, 1, sd) / sqrt(opt$replications)
627 # SE for F2 639 # SE for F2
631 # SE for F1 adult 643 # SE for F1 adult
632 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) 644 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
633 # SE for F2 adult 645 # SE for F2 adult
634 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) 646 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
635 647
636 dev.new(width = 9, height = 9) 648 dev.new(width=20, height=20)
637 649
638 # Start PDF device driver to save charts to output. 650 # Start PDF device driver to save charts to output.
639 pdf(file=opt$output, height=20, width=20, bg="white") 651 pdf(file=opt$output, height=20, width=20, bg="white")
640 652
641 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1)) 653 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))