Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 59:892cf703be62 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 21 Nov 2018 11:42:37 -0500 (2018-11-21)
parents 2194155309f4
children 393085589438
comparison
equal deleted inserted replaced
58:2194155309f4 59:892cf703be62
248 mortality.probability = temperature * 0.0005 + 0.02; 248 mortality.probability = temperature * 0.0005 + 0.02;
249 } 249 }
250 return(mortality.probability) 250 return(mortality.probability)
251 } 251 }
252 252
253 #mortality.egg = function(temperature, adj=0) { 253 mortality.egg = function(temperature, adj=0) {
254 # # If no input from adjustment, default 254 # If no input from adjustment, default
255 # # value is 0 (data from Nielsen, 2008). 255 # value is 0 (data from Nielsen, 2008).
256 # T.mortality = c(15, 17, 20, 25, 27, 30, 33, 35); 256 T.mortality = c(15, 17, 20, 25, 27, 30, 33, 35);
257 # egg.mortality = c(50, 2, 1, 0, 0, 0, 5, 100); 257 egg.mortality = c(50, 2, 1, 0, 0, 0, 5, 100);
258 # # Calculates slopes and intercepts for lines. 258 # Calculates slopes and intercepts for lines.
259 # slopes = NULL; 259 slopes = NULL;
260 # intercepts = NULL; 260 intercepts = NULL;
261 # for (i in 1:length(T.mortality)) { 261 for (i in 1:length(T.mortality)) {
262 # slopes[i] = (egg.mortality[i+1] - egg.mortality[i]) / (T.mortality[i+1] - T.mortality[i]); 262 slopes[i] = (egg.mortality[i+1] - egg.mortality[i]) / (T.mortality[i+1] - T.mortality[i]);
263 # intercepts[i] = -slopes[i] * T.mortality[i] + egg.mortality[i]; 263 intercepts[i] = -slopes[i] * T.mortality[i] + egg.mortality[i];
264 # } 264 }
265 # # Calculates mortality based on temperature. 265 # Calculates mortality based on temperature.
266 # mortality.probability = NULL; 266 mortality.probability = NULL;
267 # for (j in 1:length(temperature)) { 267 for (j in 1:length(temperature)) {
268 # mortality.probability[j] = if(temperature[j] <= T.mortality[2]) { 268 mortality.probability[j] = if(temperature[j] <= T.mortality[2]) {
269 # temperature[j] * slopes[1] + intercepts[1]; 269 temperature[j] * slopes[1] + intercepts[1];
270 # } else if (temperature[j] > T.mortality[2] && temperature[j] <= T.mortality[3]) { 270 } else if (temperature[j] > T.mortality[2] && temperature[j] <= T.mortality[3]) {
271 # temperature[j] * slopes[2] + intercepts[2]; 271 temperature[j] * slopes[2] + intercepts[2];
272 # } else if (temperature[j] > T.mortality[3] && temperature[j] <= T.mortality[4]) { 272 } else if (temperature[j] > T.mortality[3] && temperature[j] <= T.mortality[4]) {
273 # temperature[j] * slopes[3] + intercepts[3]; 273 temperature[j] * slopes[3] + intercepts[3];
274 # } else if (temperature[j] > T.mortality[4] && temperature[j] <= T.mortality[5]) { 274 } else if (temperature[j] > T.mortality[4] && temperature[j] <= T.mortality[5]) {
275 # temperature[j] * slopes[4] + intercepts[4]; 275 temperature[j] * slopes[4] + intercepts[4];
276 # } else if (temperature[j] > T.mortality[5] && temperature[j] <= T.mortality[6]) { 276 } else if (temperature[j] > T.mortality[5] && temperature[j] <= T.mortality[6]) {
277 # temperature[j] * slopes[5] + intercepts[5]; 277 temperature[j] * slopes[5] + intercepts[5];
278 # } else if (temperature[j] > T.mortality[6] && temperature[j] <= T.mortality[7]) { 278 } else if (temperature[j] > T.mortality[6] && temperature[j] <= T.mortality[7]) {
279 # temperature[j] * slopes[6] + intercepts[6]; 279 temperature[j] * slopes[6] + intercepts[6];
280 # } else if (temperature[j] > T.mortality[7]) { 280 } else if (temperature[j] > T.mortality[7]) {
281 # temperature[j] * slopes[7] + intercepts[7]; 281 temperature[j] * slopes[7] + intercepts[7];
282 # } 282 }
283 # # If mortality > 100, make it equal to 100. 283 # If mortality > 100, make it equal to 100.
284 # mortality.probability[mortality.probability>100] = 100; 284 mortality.probability[mortality.probability>100] = 100;
285 # # If mortality <0, make equal to 0. 285 # If mortality <0, make equal to 0.
286 # mortality.probability[mortality.probability<0] = 0; 286 mortality.probability[mortality.probability<0] = 0;
287 # } 287 }
288 # # Make mortality adjustments based on adj parameter. 288 # Make mortality adjustments based on adj parameter.
289 # mortality.probability = (100 - mortality.probability) * adj + mortality.probability; 289 mortality.probability = (100 - mortality.probability) * adj + mortality.probability;
290 # # if mortality > 100, make it equal to 100. 290 # if mortality > 100, make it equal to 100.
291 # mortality.probability[mortality.probability>100] = 100; 291 mortality.probability[mortality.probability>100] = 100;
292 # # If mortality <0, make equal to 0. 292 # If mortality <0, make equal to 0.
293 # mortality.probability[mortality.probability<0] = 0; 293 mortality.probability[mortality.probability<0] = 0;
294 # # Change percent to proportion. 294 # Change percent to proportion.
295 # mortality.probability = mortality.probability / 100; 295 mortality.probability = mortality.probability / 100;
296 # return(mortality.probability)
297 #}
298
299 mortality.egg = function(temperature) {
300 if (temperature < 12.7) {
301 mortality.probability = 0.8;
302 }
303 else {
304 mortality.probability = 0.8 - temperature / 40.0;
305 if (mortality.probability < 0) {
306 mortality.probability = 0.01;
307 }
308 }
309 return(mortality.probability) 296 return(mortality.probability)
297 }
310 298
311 mortality.nymph = function(temperature) { 299 mortality.nymph = function(temperature) {
312 if (temperature < 12.7) { 300 if (temperature < 12.7) {
313 mortality.probability = 0.03; 301 mortality.probability = 0.03;
314 } 302 }
612 num_days = 365; 600 num_days = 365;
613 } 601 }
614 cat("Number of days in year: ", num_days, "\n"); 602 cat("Number of days in year: ", num_days, "\n");
615 } 603 }
616 604
617 # Get the ticks date labels for plots.
618 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm);
619 ticks = c(unlist(ticks_and_labels[1]));
620 date_labels = c(unlist(ticks_and_labels[2]));
621 # All latitude values are the same, so get the value for plots from the first row. 605 # All latitude values are the same, so get the value for plots from the first row.
622 latitude = temperature_data_frame$LATITUDE[1]; 606 latitude = temperature_data_frame$LATITUDE[1];
623 607
624 # Determine the specified life stages for processing. 608 # Determine the specified life stages for processing.
625 # Split life_stages into a list of strings for plots. 609 # Split life_stages into a list of strings for plots.
754 } 738 }
755 } 739 }
756 # Total population. 740 # Total population.
757 population.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications); 741 population.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
758 742
743 doy_zero_insects = NULL;
759 # Process replications. 744 # Process replications.
760 for (current_replication in 1:opt$replications) { 745 for (current_replication in 1:opt$replications) {
761 # Start with the user-defined number of insects per replication. 746 # Start with the user-defined number of insects per replication.
762 num_insects = opt$insects_per_replication; 747 num_insects = opt$insects_per_replication;
763 # Generation, Stage, degree-days, T, Diapause. 748 # Generation, Stage, degree-days, T, Diapause.
854 # Trash bin for death. 839 # Trash bin for death.
855 death.vector = NULL; 840 death.vector = NULL;
856 # Newborn. 841 # Newborn.
857 birth.vector = NULL; 842 birth.vector = NULL;
858 # All individuals. 843 # All individuals.
859 for (i in 1:num_insects) { 844 if (num_insects > 0) {
860 # Individual record. 845 for (i in 1:num_insects) {
861 vector.individual = vector.matrix[i,]; 846 # Individual record.
862 # Adjustment for late season mortality rate (still alive?). 847 vector.individual = vector.matrix[i,];
863 if (latitude < 40.0) { 848 # Adjustment for late season mortality rate (still alive?).
864 post.mortality = 1; 849 if (latitude < 40.0) {
865 day.kill = 300; 850 post.mortality = 1;
866 } 851 day.kill = 300;
867 else {
868 post.mortality = 2;
869 day.kill = 250;
870 }
871 if (vector.individual[2] == 0) {
872 # Egg.
873 #death.probability = opt$egg_mortality * mortality.egg(mean.temp, adj=opt$egg_mortality);
874 death.probability = opt$egg_mortality * mortality.egg(mean.temp);
875 }
876 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
877 # Nymph.
878 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp);
879 }
880 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
881 # Adult.
882 if (doy < day.kill) {
883 death.probability = opt$adult_mortality * mortality.adult(mean.temp);
884 } 852 }
885 else { 853 else {
886 # Increase adult mortality after fall equinox. 854 post.mortality = 2;
887 death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp); 855 day.kill = 250;
888 } 856 }
889 } 857 if (vector.individual[2] == 0) {
890 # Dependent on temperature and life stage? 858 # Egg.
891 u.d = runif(1); 859 death.probability = opt$egg_mortality * mortality.egg(mean.temp, adj=opt$egg_mortality);
892 if (u.d < death.probability) { 860 }
893 death.vector = c(death.vector, i); 861 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
894 } 862 # Nymph.
895 else { 863 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp);
896 # End of diapause. 864 }
897 if (vector.individual[1] == 0 && vector.individual[2] == 3) { 865 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
898 # Overwintering adult (pre-vittelogenic). 866 # Adult.
899 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) { 867 if (doy < day.kill) {
900 # Add 68C to become fully reproductively matured. 868 death.probability = opt$adult_mortality * mortality.adult(mean.temp);
901 # Transfer to vittelogenic.
902 vector.individual = c(0, 4, 0, 0, 0);
903 vector.matrix[i,] = vector.individual;
904 } 869 }
905 else { 870 else {
871 # Increase adult mortality after fall equinox.
872 death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp);
873 }
874 }
875 # Dependent on temperature and life stage?
876 u.d = runif(1);
877 if (u.d < death.probability) {
878 death.vector = c(death.vector, i);
879 }
880 else {
881 # End of diapause.
882 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
883 # Overwintering adult (pre-vittelogenic).
884 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) {
885 # Add 68C to become fully reproductively matured.
886 # Transfer to vittelogenic.
887 vector.individual = c(0, 4, 0, 0, 0);
888 vector.matrix[i,] = vector.individual;
889 }
890 else {
891 # Add average temperature for current day.
892 vector.individual[3] = vector.individual[3] + averages.temp;
893 # Add 1 day in current stage.
894 vector.individual[4] = vector.individual[4] + 1;
895 vector.matrix[i,] = vector.individual;
896 }
897 }
898 if (vector.individual[1] != 0 && vector.individual[2] == 3) {
899 # Not overwintering adult (pre-vittelogenic).
900 current.gen = vector.individual[1];
901 if (vector.individual[3] > 68) {
902 # Add 68C to become fully reproductively matured.
903 # Transfer to vittelogenic.
904 vector.individual = c(current.gen, 4, 0, 0, 0);
905 vector.matrix[i,] = vector.individual;
906 }
907 else {
908 # Add average temperature for current day.
909 vector.individual[3] = vector.individual[3] + averages.temp;
910 # Add 1 day in current stage.
911 vector.individual[4] = vector.individual[4] + 1;
912 vector.matrix[i,] = vector.individual;
913 }
914 }
915 # Oviposition -- where population dynamics comes from.
916 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) {
917 # Vittelogenic stage, overwintering generation.
918 if (vector.individual[4] == 0) {
919 # Just turned in vittelogenic stage.
920 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size));
921 }
922 else {
923 # Daily probability of birth.
924 p.birth = opt$oviposition * 0.01;
925 u1 = runif(1);
926 if (u1 < p.birth) {
927 num_insects.birth = round(runif(1, 2, 8));
928 }
929 }
906 # Add average temperature for current day. 930 # Add average temperature for current day.
907 vector.individual[3] = vector.individual[3] + averages.temp; 931 vector.individual[3] = vector.individual[3] + averages.temp;
908 # Add 1 day in current stage. 932 # Add 1 day in current stage.
909 vector.individual[4] = vector.individual[4] + 1; 933 vector.individual[4] = vector.individual[4] + 1;
910 vector.matrix[i,] = vector.individual; 934 vector.matrix[i,] = vector.individual;
935 if (num_insects.birth > 0) {
936 # Add new birth -- might be in different generations.
937 new.gen = vector.individual[1] + 1;
938 # Egg profile.
939 new.individual = c(new.gen, 0, 0, 0, 0);
940 new.vector = rep(new.individual, num_insects.birth);
941 # Update batch of egg profile.
942 new.vector = t(matrix(new.vector, nrow=5));
943 # Group with total eggs laid in that day.
944 birth.vector = rbind(birth.vector, new.vector);
945 }
911 } 946 }
912 } 947 # Oviposition -- for generation 1.
913 if (vector.individual[1] != 0 && vector.individual[2] == 3) { 948 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
914 # Not overwintering adult (pre-vittelogenic). 949 # Vittelogenic stage, 1st generation
915 current.gen = vector.individual[1]; 950 if (vector.individual[4] == 0) {
916 if (vector.individual[3] > 68) { 951 # Just turned in vittelogenic stage.
917 # Add 68C to become fully reproductively matured. 952 num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size));
918 # Transfer to vittelogenic. 953 }
919 vector.individual = c(current.gen, 4, 0, 0, 0); 954 else {
920 vector.matrix[i,] = vector.individual; 955 # Daily probability of birth.
921 } 956 p.birth = opt$oviposition * 0.01;
922 else { 957 u1 = runif(1);
958 if (u1 < p.birth) {
959 num_insects.birth = round(runif(1, 2, 8));
960 }
961 }
923 # Add average temperature for current day. 962 # Add average temperature for current day.
924 vector.individual[3] = vector.individual[3] + averages.temp; 963 vector.individual[3] = vector.individual[3] + averages.temp;
925 # Add 1 day in current stage. 964 # Add 1 day in current stage.
926 vector.individual[4] = vector.individual[4] + 1; 965 vector.individual[4] = vector.individual[4] + 1;
927 vector.matrix[i,] = vector.individual; 966 vector.matrix[i,] = vector.individual;
928 } 967 if (num_insects.birth > 0) {
929 } 968 # Add new birth -- might be in different generations.
930 # Oviposition -- where population dynamics comes from. 969 new.gen = vector.individual[1] + 1;
931 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) { 970 # Egg profile.
932 # Vittelogenic stage, overwintering generation. 971 new.individual = c(new.gen, 0, 0, 0, 0);
933 if (vector.individual[4] == 0) { 972 new.vector = rep(new.individual, num_insects.birth);
934 # Just turned in vittelogenic stage. 973 # Update batch of egg profile.
935 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)); 974 new.vector = t(matrix(new.vector, nrow=5));
936 } 975 # Group with total eggs laid in that day.
937 else { 976 birth.vector = rbind(birth.vector, new.vector);
938 # Daily probability of birth.
939 p.birth = opt$oviposition * 0.01;
940 u1 = runif(1);
941 if (u1 < p.birth) {
942 num_insects.birth = round(runif(1, 2, 8));
943 } 977 }
944 } 978 }
945 # Add average temperature for current day. 979 # Egg to young nymph.
946 vector.individual[3] = vector.individual[3] + averages.temp; 980 if (vector.individual[2] == 0) {
947 # Add 1 day in current stage. 981 # Add average temperature for current day.
948 vector.individual[4] = vector.individual[4] + 1; 982 vector.individual[3] = vector.individual[3] + averages.temp;
949 vector.matrix[i,] = vector.individual; 983 if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) {
950 if (num_insects.birth > 0) { 984 # From egg to young nymph, degree-days requirement met.
951 # Add new birth -- might be in different generations. 985 current.gen = vector.individual[1];
952 new.gen = vector.individual[1] + 1; 986 # Transfer to young nymph stage.
953 # Egg profile. 987 vector.individual = c(current.gen, 1, 0, 0, 0);
954 new.individual = c(new.gen, 0, 0, 0, 0);
955 new.vector = rep(new.individual, num_insects.birth);
956 # Update batch of egg profile.
957 new.vector = t(matrix(new.vector, nrow=5));
958 # Group with total eggs laid in that day.
959 birth.vector = rbind(birth.vector, new.vector);
960 }
961 }
962 # Oviposition -- for generation 1.
963 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
964 # Vittelogenic stage, 1st generation
965 if (vector.individual[4] == 0) {
966 # Just turned in vittelogenic stage.
967 num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size));
968 }
969 else {
970 # Daily probability of birth.
971 p.birth = opt$oviposition * 0.01;
972 u1 = runif(1);
973 if (u1 < p.birth) {
974 num_insects.birth = round(runif(1, 2, 8));
975 }
976 }
977 # Add average temperature for current day.
978 vector.individual[3] = vector.individual[3] + averages.temp;
979 # Add 1 day in current stage.
980 vector.individual[4] = vector.individual[4] + 1;
981 vector.matrix[i,] = vector.individual;
982 if (num_insects.birth > 0) {
983 # Add new birth -- might be in different generations.
984 new.gen = vector.individual[1] + 1;
985 # Egg profile.
986 new.individual = c(new.gen, 0, 0, 0, 0);
987 new.vector = rep(new.individual, num_insects.birth);
988 # Update batch of egg profile.
989 new.vector = t(matrix(new.vector, nrow=5));
990 # Group with total eggs laid in that day.
991 birth.vector = rbind(birth.vector, new.vector);
992 }
993 }
994 # Egg to young nymph.
995 if (vector.individual[2] == 0) {
996 # Add average temperature for current day.
997 vector.individual[3] = vector.individual[3] + averages.temp;
998 if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) {
999 # From egg to young nymph, degree-days requirement met.
1000 current.gen = vector.individual[1];
1001 # Transfer to young nymph stage.
1002 vector.individual = c(current.gen, 1, 0, 0, 0);
1003 }
1004 else {
1005 # Add 1 day in current stage.
1006 vector.individual[4] = vector.individual[4] + 1;
1007 }
1008 vector.matrix[i,] = vector.individual;
1009 }
1010 # Young nymph to old nymph.
1011 if (vector.individual[2] == 1) {
1012 # Add average temperature for current day.
1013 vector.individual[3] = vector.individual[3] + averages.temp;
1014 if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) {
1015 # From young to old nymph, degree_days requirement met.
1016 current.gen = vector.individual[1];
1017 # Transfer to old nym stage.
1018 vector.individual = c(current.gen, 2, 0, 0, 0);
1019 if (photoperiod < opt$photoperiod && doy > 180) {
1020 vector.individual[5] = 1;
1021 } # Prepare for diapausing.
1022 }
1023 else {
1024 # Add 1 day in current stage.
1025 vector.individual[4] = vector.individual[4] + 1;
1026 }
1027 vector.matrix[i,] = vector.individual;
1028 }
1029 # Old nymph to adult: pre-vittelogenic or diapausing?
1030 if (vector.individual[2] == 2) {
1031 # Add average temperature for current day.
1032 vector.individual[3] = vector.individual[3] + averages.temp;
1033 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
1034 # From old to adult, degree_days requirement met.
1035 current.gen = vector.individual[1];
1036 if (vector.individual[5] == 0) {
1037 # Previttelogenic.
1038 vector.individual = c(current.gen, 3, 0, 0, 0);
1039 } 988 }
1040 else { 989 else {
1041 # Diapausing. 990 # Add 1 day in current stage.
1042 vector.individual = c(current.gen, 5, 0, 0, 1); 991 vector.individual[4] = vector.individual[4] + 1;
1043 } 992 }
993 vector.matrix[i,] = vector.individual;
1044 } 994 }
1045 else { 995 # Young nymph to old nymph.
1046 # Add 1 day in current stage. 996 if (vector.individual[2] == 1) {
997 # Add average temperature for current day.
998 vector.individual[3] = vector.individual[3] + averages.temp;
999 if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) {
1000 # From young to old nymph, degree_days requirement met.
1001 current.gen = vector.individual[1];
1002 # Transfer to old nym stage.
1003 vector.individual = c(current.gen, 2, 0, 0, 0);
1004 if (photoperiod < opt$photoperiod && doy > 180) {
1005 vector.individual[5] = 1;
1006 } # Prepare for diapausing.
1007 }
1008 else {
1009 # Add 1 day in current stage.
1010 vector.individual[4] = vector.individual[4] + 1;
1011 }
1012 vector.matrix[i,] = vector.individual;
1013 }
1014 # Old nymph to adult: pre-vittelogenic or diapausing?
1015 if (vector.individual[2] == 2) {
1016 # Add average temperature for current day.
1017 vector.individual[3] = vector.individual[3] + averages.temp;
1018 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
1019 # From old to adult, degree_days requirement met.
1020 current.gen = vector.individual[1];
1021 if (vector.individual[5] == 0) {
1022 # Previttelogenic.
1023 vector.individual = c(current.gen, 3, 0, 0, 0);
1024 }
1025 else {
1026 # Diapausing.
1027 vector.individual = c(current.gen, 5, 0, 0, 1);
1028 }
1029 }
1030 else {
1031 # Add 1 day in current stage.
1032 vector.individual[4] = vector.individual[4] + 1;
1033 }
1034 vector.matrix[i,] = vector.individual;
1035 }
1036 # Growing of diapausing adult (unimportant, but still necessary).
1037 if (vector.individual[2] == 5) {
1038 vector.individual[3] = vector.individual[3] + averages.temp;
1047 vector.individual[4] = vector.individual[4] + 1; 1039 vector.individual[4] = vector.individual[4] + 1;
1040 vector.matrix[i,] = vector.individual;
1048 } 1041 }
1049 vector.matrix[i,] = vector.individual; 1042 } # Else if it is still alive.
1050 } 1043 } # End of the individual bug loop.
1051 # Growing of diapausing adult (unimportant, but still necessary). 1044
1052 if (vector.individual[2] == 5) { 1045 # Number of deaths.
1053 vector.individual[3] = vector.individual[3] + averages.temp; 1046 num_insects.death = length(death.vector);
1054 vector.individual[4] = vector.individual[4] + 1; 1047 if (num_insects.death > 0) {
1055 vector.matrix[i,] = vector.individual; 1048 # Remove record of dead.
1056 } 1049 vector.matrix = vector.matrix[-death.vector,];
1057 } # Else if it is still alive. 1050 }
1058 } # End of the individual bug loop. 1051 # Number of births.
1059 1052 num_insects.newborn = length(birth.vector[,1]);
1060 # Number of deaths. 1053 vector.matrix = rbind(vector.matrix, birth.vector);
1061 num_insects.death = length(death.vector); 1054 # Update population size for the next day.
1062 if (num_insects.death > 0) { 1055 num_insects = num_insects - num_insects.death + num_insects.newborn;
1063 # Remove record of dead. 1056 } else {
1064 vector.matrix = vector.matrix[-death.vector,]; 1057 if (is.null(doy_zero_insects)) {
1065 } 1058 # Only set the doy for zero insects if
1066 # Number of births. 1059 # it has not yet been set.
1067 num_insects.newborn = length(birth.vector[,1]); 1060 doy_zero_insects = doy;
1068 vector.matrix = rbind(vector.matrix, birth.vector); 1061 }
1069 # Update population size for the next day. 1062 }
1070 num_insects = num_insects - num_insects.death + num_insects.newborn;
1071 1063
1072 # Aggregate results by day. Due to multiple transpose calls 1064 # Aggregate results by day. Due to multiple transpose calls
1073 # on vector.matrix above, the columns of vector.matrix 1065 # on vector.matrix above, the columns of vector.matrix
1074 # are now Generation, Stage, degree-days, T, Diapause, 1066 # are now Generation, Stage, degree-days, T, Diapause,
1075 if (process_eggs) { 1067 if (process_eggs) {
1545 # Save the analyzed data for generation F2. 1537 # Save the analyzed data for generation F2.
1546 file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/"); 1538 file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/");
1547 write.csv(temperature_data_frame_F2, file=file_path, row.names=F); 1539 write.csv(temperature_data_frame_F2, file=file_path, row.names=F);
1548 } 1540 }
1549 1541
1542 # Get the ticks date labels for plots.
1543 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm=prepend_end_doy_norm, append_start_doy_norm=append_start_doy_norm, date_interval=FALSE, doy_zero_insects=doy_zero_insects);
1544 ticks = c(unlist(ticks_and_labels[1]));
1545 date_labels = c(unlist(ticks_and_labels[2]));
1550 total_days_vector = c(1:dim(temperature_data_frame)[1]); 1546 total_days_vector = c(1:dim(temperature_data_frame)[1]);
1547
1551 if (plot_generations_separately) { 1548 if (plot_generations_separately) {
1552 for (life_stage in life_stages) { 1549 for (life_stage in life_stages) {
1553 if (life_stage == "Egg") { 1550 if (life_stage == "Egg") {
1554 # Start PDF device driver. 1551 # Start PDF device driver.
1555 dev.new(width=20, height=30); 1552 dev.new(width=20, height=30);