comparison insect_phenology_model.R @ 45:315c5e1bc44a draft

Uploaded

44:c61d3d9d44db

    make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
    make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
    make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
    make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
    make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
    make_option(c("--location"), action="store", dest="location", help="Selected location"),
    make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
    make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
    make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", default=NULL, type="integer", help="Total number of days in the year-to-date temperature dataset"),
    make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
    make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),

        mortality.probability = temperature * 0.0008 + 0.03;
    }
    return(mortality.probability);
}

parse_input_data = function(input_ytd, input_norm, num_days_ytd) {
    if (is.null(input_ytd)) {
        # We're analysing only the 30 year normals data, so create an empty
        # data frame for containing temperature data after it is converted
        # from the 30 year normals format to the year-to-date format.
        temperature_data_frame = data.frame(matrix(ncol=6, nrow=0));

    # All normals data includes Feb 29 which is row 60 in
    # the data, so delete that row if we're not in a leap year.
    if (!is_leap_year) {
        norm_data_frame = norm_data_frame[-c(60),];
    }
    if (is.null(input_ytd)) {
        # Convert the 30 year normals data to the year-to-date format.
        for (i in 1:total_days) {
            temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
        }

            temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
        }
    }
    # Add a column containing the daylight length for each day.
    temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
    return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days));
}

render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
    replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
    life_stages_adult=NULL, life_stages_nymph=NULL) {

}
# Display the total number of days in the Galaxy history item blurb.
cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");

# Parse the inputs.
data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd);
temperature_data_frame = data_list[[1]];
# Information needed for plots.
start_date = data_list[[2]];
end_date = data_list[[3]];
start_doy_ytd = data_list[[4]];
end_doy_ytd = data_list[[5]];
is_leap_year = data_list[[6]];
total_days = data_list[[7]];
total_days_vector = c(1:total_days);

# Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
if (plot_generations_separately) {
    temperature_data_frame_P = data.frame(temperature_data_frame);
    temperature_data_frame_F1 = data.frame(temperature_data_frame);

            file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Egg population size by generation.
            maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
            render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error,
                group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        } else if (life_stage == "Nymph") {

                    group2 = F1_total_nymphs;
                    group2_std_error = F1_total_nymphs.std_error;
                    group3 = F2_total_nymphs;
                    group3_std_error = F2_total_nymphs.std_error;
                }
                render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
                # Turn off device driver to flush output.
                dev.off();
            }

                    group2 = F1_total_adults;
                    group2_std_error = F1_total_adults.std_error;
                    group3 = F2_total_adults;
                    group3_std_error = F2_total_adults.std_error;
                }
                render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
                # Turn off device driver to flush output.
                dev.off();
            }

            file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Total population size by generation.
            maxval = max(P+F1+F2) + 100;
            render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error,
                group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
            # Turn off device driver to flush output.
            dev.off();
        }

            file_path = get_file_path(life_stage, "egg_pop.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Egg population size.
            maxval = max(eggs+eggs.std_error) + 100;
            render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        } else if (life_stage == "Nymph") {
            for (life_stage_nymph in life_stages_nymph) {

                    # Old nymph population size.
                    group = old_nymphs;
                    group_std_error = old_nymphs.std_error;
                }
                maxval = max(group+group_std_error) + 100;
                render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    life_stages_nymph=life_stage_nymph);
                # Turn off device driver to flush output.
                dev.off();
            }

                    # Diapausing adult population size.
                    group = diapausing_adults;
                    group_std_error = diapausing_adults.std_error
                }
                maxval = max(group+group_std_error) + 100;
                render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    life_stages_adult=life_stage_adult);
                # Turn off device driver to flush output.
                dev.off();
            }

            file_path = get_file_path(life_stage, "total_pop.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Total population size.
            maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
            render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error,
                group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        }

45:315c5e1bc44a

    make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
    make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
    make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
    make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
    make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
    make_option(c("--location"), action="store", dest="location", default=NULL, help="Selected location"),
    make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
    make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
    make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", default=NULL, type="integer", help="Total number of days in the year-to-date temperature dataset"),
    make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
    make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),

        mortality.probability = temperature * 0.0008 + 0.03;
    }
    return(mortality.probability);
}

parse_input_data = function(input_ytd, input_norm, num_days_ytd, location) {
    if (is.null(input_ytd)) {
        # We're analysing only the 30 year normals data, so create an empty
        # data frame for containing temperature data after it is converted
        # from the 30 year normals format to the year-to-date format.
        temperature_data_frame = data.frame(matrix(ncol=6, nrow=0));

    # All normals data includes Feb 29 which is row 60 in
    # the data, so delete that row if we're not in a leap year.
    if (!is_leap_year) {
        norm_data_frame = norm_data_frame[-c(60),];
    }
    # Set the location to be the station name if the user elected no to enter it.
    if (is.null(location)) {
        location = norm_data_frame$NAME[1];
    }
    if (is.null(input_ytd)) {
        # Convert the 30 year normals data to the year-to-date format.
        for (i in 1:total_days) {
            temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
        }

            temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
        }
    }
    # Add a column containing the daylight length for each day.
    temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
    return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days, location));
}

render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
    replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
    life_stages_adult=NULL, life_stages_nymph=NULL) {

}
# Display the total number of days in the Galaxy history item blurb.
cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");

# Parse the inputs.
data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd, opt$location);
temperature_data_frame = data_list[[1]];
# Information needed for plots.
start_date = data_list[[2]];
end_date = data_list[[3]];
start_doy_ytd = data_list[[4]];
end_doy_ytd = data_list[[5]];
is_leap_year = data_list[[6]];
total_days = data_list[[7]];
total_days_vector = c(1:total_days);
location =  data_list[[8]];

# Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
if (plot_generations_separately) {
    temperature_data_frame_P = data.frame(temperature_data_frame);
    temperature_data_frame_F1 = data.frame(temperature_data_frame);

            file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Egg population size by generation.
            maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
            render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error,
                group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        } else if (life_stage == "Nymph") {

                    group2 = F1_total_nymphs;
                    group2_std_error = F1_total_nymphs.std_error;
                    group3 = F2_total_nymphs;
                    group3_std_error = F2_total_nymphs.std_error;
                }
                render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
                # Turn off device driver to flush output.
                dev.off();
            }

                    group2 = F1_total_adults;
                    group2_std_error = F1_total_adults.std_error;
                    group3 = F2_total_adults;
                    group3_std_error = F2_total_adults.std_error;
                }
                render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
                # Turn off device driver to flush output.
                dev.off();
            }

            file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Total population size by generation.
            maxval = max(P+F1+F2) + 100;
            render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error,
                group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
            # Turn off device driver to flush output.
            dev.off();
        }

            file_path = get_file_path(life_stage, "egg_pop.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Egg population size.
            maxval = max(eggs+eggs.std_error) + 100;
            render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        } else if (life_stage == "Nymph") {
            for (life_stage_nymph in life_stages_nymph) {

                    # Old nymph population size.
                    group = old_nymphs;
                    group_std_error = old_nymphs.std_error;
                }
                maxval = max(group+group_std_error) + 100;
                render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    life_stages_nymph=life_stage_nymph);
                # Turn off device driver to flush output.
                dev.off();
            }

                    # Diapausing adult population size.
                    group = diapausing_adults;
                    group_std_error = diapausing_adults.std_error
                }
                maxval = max(group+group_std_error) + 100;
                render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
                    start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
                    life_stages_adult=life_stage_adult);
                # Turn off device driver to flush output.
                dev.off();
            }

            file_path = get_file_path(life_stage, "total_pop.pdf")
            pdf(file=file_path, width=20, height=30, bg="white");
            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
            # Total population size.
            maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
            render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
                start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error,
                group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
            # Turn off device driver to flush output.
            dev.off();
        }