diff --git a/.Rbuildignore b/.Rbuildignore index 78ea75d2..38dc18f6 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,5 +14,6 @@ ^data-raw$ ^docs$ ^man-roxygen$ +^organisation.yml$ ^pkgdown$ ^wercker\.yml$ diff --git a/.github/workflows/check_on_branch.yml b/.github/workflows/check_on_branch.yml index 80bf4d8f..182dd948 100644 --- a/.github/workflows/check_on_branch.yml +++ b/.github/workflows/check_on_branch.yml @@ -2,7 +2,6 @@ on: push: branches-ignore: - main - - master - ghpages name: "check package with checklist" @@ -14,11 +13,11 @@ jobs: env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + N2KBUCKET: ${{ secrets.N2KBUCKET }} AWS_ACCESS_KEY_ID: ${{ secrets.AWS_ACCESS_KEY_ID }} AWS_DEFAULT_REGION: ${{ secrets.AWS_DEFAULT_REGION }} AWS_SECRET_ACCESS_KEY: ${{ secrets.AWS_SECRET_ACCESS_KEY }} - N2KBUCKET: ${{ secrets.N2KBUCKET }} permissions: contents: read steps: - - uses: inbo/actions/check_pkg@checklist-0.3.6 + - uses: inbo/actions/check_pkg@checklist-0.4.1 diff --git a/.github/workflows/check_on_different_r_os.yml b/.github/workflows/check_on_different_r_os.yml index 2ac75562..bb243621 100644 --- a/.github/workflows/check_on_different_r_os.yml +++ b/.github/workflows/check_on_different_r_os.yml @@ -21,6 +21,7 @@ jobs: matrix: config: - {os: macOS-latest, r: 'release'} + - {os: windows-latest, r: 'release'} - {os: ubuntu-22.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/jammy/latest"} - {os: ubuntu-22.04, r: 'oldrel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/jammy/latest"} @@ -40,7 +41,7 @@ jobs: - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} - extra-repositories: https://inla.r-inla-download.org/R/testing + extra-repositories: https://inla.r-inla-download.org/R/stable - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/check_on_main.yml b/.github/workflows/check_on_main.yml index 57e2a642..3ffc91a0 100644 --- a/.github/workflows/check_on_main.yml +++ b/.github/workflows/check_on_main.yml @@ -15,9 +15,9 @@ jobs: env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + N2KBUCKET: ${{ secrets.N2KBUCKET }} AWS_ACCESS_KEY_ID: ${{ secrets.AWS_ACCESS_KEY_ID }} AWS_DEFAULT_REGION: ${{ secrets.AWS_DEFAULT_REGION }} AWS_SECRET_ACCESS_KEY: ${{ secrets.AWS_SECRET_ACCESS_KEY }} - N2KBUCKET: ${{ secrets.N2KBUCKET }} steps: - - uses: inbo/actions/check_pkg@checklist-0.3.6 + - uses: inbo/actions/check_pkg@checklist-0.4.1 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 8146b6a0..9b10b9f6 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -10,32 +10,21 @@ on: - completed jobs: - prepare: + publish: runs-on: ubuntu-latest - outputs: - tag: ${{ steps.gettag.outputs.tag }} - body: ${{ steps.gettag.outputs.body }} + permissions: + contents: write steps: - uses: actions/checkout@v3 - name: Get tag run: | git fetch --tags --force TAG=$(git tag --contains $(git rev-parse HEAD)) - echo "tag=$TAG" >> "$GITHUB_OUTPUT" TAG_BODY=$(git tag --contains $(git rev-parse HEAD) --format='%(contents)') - TAG_BODY="${TAG_BODY//'%'/'%25'}" - TAG_BODY="${TAG_BODY//$'\n'/'%0A'}" - TAG_BODY="${TAG_BODY//$'\r'/'%0D'}" - echo "body=$TAG_BODY" >> "$GITHUB_OUTPUT" - id: gettag - publish: - runs-on: ubuntu-latest - permissions: - contents: write - needs: prepare - steps: + echo "tag=$TAG" >> $GITHUB_ENV + echo "$TAG_BODY" > body.md - uses: ncipollo/release-action@v1 with: - name: Release ${{needs.prepare.outputs.tag}} - tag: ${{needs.prepare.outputs.tag}} - body: ${{needs.prepare.outputs.body}} + name: Release ${{ env.tag }} + tag: ${{ env.tag }} + bodyFile: body.md diff --git a/.gitignore b/.gitignore index 84dba30a..8de6414e 100644 --- a/.gitignore +++ b/.gitignore @@ -17,7 +17,7 @@ .Ruserdata .httr-oauth docs +inst/doc libs output renv/library -inst/doc diff --git a/.zenodo.json b/.zenodo.json index b1088cec..1d1a09aa 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,9 +1,9 @@ { "title": "n2kanalysis: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring", - "version": "0.3.2", + "version": "0.4.0", "license": "GPL-3.0", "upload_type": "software", - "description": "

All generic functions and classes for the analysis for the ‘Natura\n2000’ monitoring. The classes contain all required data and definitions\nto fit the model without the need to access other sources. Potentially\nthey might need access to one or more parent objects. An aggregation\nobject might for example need the result of an imputation object. The\nactual definition of the analysis, using these generic function and\nclasses, is defined in dedictated analysis R packages for every\nmonitoring scheme. For example ‘abvanalysis’ and\n‘watervogelanalysis’.<\/p>", + "description": "

All generic functions and classes for the analysis for the ‘Natura 2000’ monitoring. The classes contain all required data and definitions to fit the model without the need to access other sources. Potentially they might need access to one or more parent objects. An aggregation object might for example need the result of an imputation object. The actual definition of the analysis, using these generic function and classes, is defined in dedictated analysis R packages for every monitoring scheme. For example ‘abvanalysis’ and ‘watervogelanalysis’.<\/p>", "keywords": [ "analysis, reproducible research, natura 2000, monitoring" ], @@ -14,11 +14,11 @@ "name": "Onkelinx, Thierry", "affiliation": "Research Institute for Nature and Forest (INBO)", "orcid": "0000-0001-8804-4216", - "type": "ContactPerson" + "type": "contactperson" }, { "name": "Research Institute for Nature and Forest (INBO)", - "type": "RightsHolder" + "type": "rightsholder" } ], "creators": [ diff --git a/CITATION.cff b/CITATION.cff index dcd4a619..6f445561 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -27,4 +27,6 @@ abstract: "All generic functions and classes for the analysis for the 'Natura 20 identifiers: - type: doi value: 10.5281/zenodo.3576047 -version: 0.3.2 +- type: url + value: https://inbo.github.io/n2kanalysis/ +version: 0.4.0 diff --git a/DESCRIPTION b/DESCRIPTION index 173d8bfe..f5d110de 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: n2kanalysis Title: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring -Version: 0.3.2 +Version: 0.4.0 Authors@R: c( person("Thierry", "Onkelinx", , "thierry.onkelinx@inbo.be", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8804-4216", affiliation = "Research Institute for Nature and Forest (INBO)")), @@ -16,7 +16,8 @@ Description: All generic functions and classes for the analysis for the every monitoring scheme. For example 'abvanalysis' and 'watervogelanalysis'. License: GPL-3 -URL: https://doi.org/10.5281/zenodo.3576047 +URL: https://doi.org/10.5281/zenodo.3576047, + https://inbo.github.io/n2kanalysis/ BugReports: https://github.com/inbo/n2kanalysis/issues Depends: R (>= 4.2.0) @@ -26,22 +27,26 @@ Imports: digest (>= 0.6.23.2), dplyr, fs, + git2rdata (>= 0.5.0), MASS, methods, - multimput (>= 0.2.13), + multimput (>= 0.2.14), n2khelper (>= 0.5.0), purrr, rlang, RODBC, tibble, tidyr (>= 0.4.0), + tools, yaml Suggests: + fmesher, INLA (>= 23.04.24), knitr, Matrix, parallel, rmarkdown, + sf, sn, testthat (>= 2.0.1) VignetteBuilder: @@ -57,7 +62,7 @@ Config/checklist/keywords: analysis, reproducible research, natura 2000, Encoding: UTF-8 Language: en-GB Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Collate: 'n2k_parameter_class.R' 'n2k_anomaly_class.R' @@ -71,24 +76,27 @@ Collate: 'delete_model.R' 'display.R' 'n2k_model_class.R' + 'n2k_inla_comparison_class.R' 'n2k_inla_class.R' + 'n2k_aggregate_class.R' + 'n2k_model_imputed_class.R' 'extract.R' 'fit_every_model.R' 'fit_model.R' 'fit_model_character.R' - 'n2k_aggregate_class.R' 'fit_model_n2k_aggregate.R' 'n2k_composite_class.R' 'fit_model_n2k_composite.R' 'n2k_hurdle_imputed_class.R' 'fit_model_n2k_hurdle_imputed.R' 'fit_model_n2k_inla.R' - 'n2k_inla_comparison_class.R' 'fit_model_n2k_inla_comparison.R' 'n2k_manifest_class.R' 'fit_model_n2k_manifest.R' - 'n2k_model_imputed_class.R' 'fit_model_n2k_model_imputed.R' + 'spde_class.R' + 'n2k_spde_class.R' + 'fit_model_n2k_spde.R' 'fit_model_s3_object.R' 'get_analysis_date.R' 'get_analysis_version.R' @@ -96,6 +104,7 @@ Collate: 'get_anomaly_n2k_inla.R' 'get_anomaly_n2k_model.R' 'get_data.R' + 'get_datafield_id.R' 'get_file_fingerprint.R' 'get_formula.R' 'get_location_group_id.R' @@ -120,8 +129,12 @@ Collate: 'get_species_group_id.R' 'get_status_fingerprint.R' 'inla_inverse.R' + 'make_a.R' 'manifest_yaml_to_bash.R' 'mark_obsolete_dataset.R' + 'moving_average.R' + 'moving_difference.R' + 'moving_trend.R' 'n2k_aggregated.R' 'n2k_composite.R' 'n2k_hurdle_imputed.R' @@ -131,6 +144,7 @@ Collate: 'n2k_inla_comparison.R' 'n2k_manifest.R' 'n2k_model_imputed.R' + 'n2k_spde.R' 'parent_status.R' 'read_manifest.R' 'read_model.R' @@ -146,6 +160,9 @@ Collate: 'select_observed_range.R' 'session_package.R' 'sha1.R' + 'spde.R' + 'spde2matern.R' + 'spde2mesh.R' 'status.R' 'store_manifest.R' 'store_manifest_yaml.R' diff --git a/LICENSE.md b/LICENSE.md index 2fb2e74d..379c1b2d 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,4 +1,4 @@ -### GNU GENERAL PUBLIC LICENSE +# GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 diff --git a/NAMESPACE b/NAMESPACE index 6d3c9271..40a0f55b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,14 +5,21 @@ S3method(sha1,inla) S3method(sha1,rawImputed) export(display) export(fit_every_model) +export(get_datafield_id) export(get_parents) export(inla_inverse) export(mark_obsolete_dataset) +export(moving_average) +export(moving_difference) +export(moving_trend) +export(order_manifest) export(select_factor_count_strictly_positive) export(select_factor_threshold) export(select_observed_range) export(sha1) +export(spde) export(union) +exportClasses(Spde) exportClasses(n2kAggregate) exportClasses(n2kAnalysisMetadata) exportClasses(n2kAnalysisVersion) @@ -28,6 +35,7 @@ exportClasses(n2kModel) exportClasses(n2kModelImputed) exportClasses(n2kParameter) exportClasses(n2kResult) +exportClasses(n2kSpde) exportMethods("parent_status<-") exportMethods("status<-") exportMethods(combine) @@ -49,6 +57,7 @@ exportMethods(get_scheme_id) exportMethods(get_seed) exportMethods(get_species_group_id) exportMethods(get_status_fingerprint) +exportMethods(make_a) exportMethods(manifest_yaml_to_bash) exportMethods(n2k_aggregate) exportMethods(n2k_composite) @@ -58,6 +67,7 @@ exportMethods(n2k_inla) exportMethods(n2k_inla_comparison) exportMethods(n2k_manifest) exportMethods(n2k_model_imputed) +exportMethods(n2k_spde) exportMethods(parent_status) exportMethods(read_manifest) exportMethods(read_model) @@ -65,6 +75,8 @@ exportMethods(read_result) exportMethods(result_estimate) exportMethods(result_metadata) exportMethods(session_package) +exportMethods(spde2matern) +exportMethods(spde2mesh) exportMethods(status) exportMethods(store_manifest) exportMethods(store_manifest_yaml) @@ -131,6 +143,10 @@ importFrom(fs,file_delete) importFrom(fs,file_exists) importFrom(fs,path) importFrom(fs,path_abs) +importFrom(git2rdata,is_git2rmeta) +importFrom(git2rdata,update_metadata) +importFrom(git2rdata,verify_vc) +importFrom(git2rdata,write_vc) importFrom(methods,new) importFrom(methods,setClass) importFrom(methods,setClassUnion) @@ -160,6 +176,7 @@ importFrom(rlang,"!!") importFrom(rlang,.data) importFrom(stats,as.formula) importFrom(stats,coef) +importFrom(stats,median) importFrom(stats,na.fail) importFrom(stats,na.omit) importFrom(stats,qnorm) @@ -169,6 +186,7 @@ importFrom(stats,terms) importFrom(tibble,rownames_to_column) importFrom(tibble,tibble) importFrom(tidyr,pivot_longer) +importFrom(tools,R_user_dir) importFrom(utils,file_test) importFrom(utils,flush.console) importFrom(utils,head) diff --git a/NEWS.md b/NEWS.md index 028e2164..e139687a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# `n2kanalysis` 0.4.0 + +* Handle INLA models with an SPDE element. + # `n2kanalysis` 0.3.2 * Make `fit_model()` more efficient. diff --git a/R/extract.R b/R/extract.R index f791a819..6b8fe4b2 100644 --- a/R/extract.R +++ b/R/extract.R @@ -63,3 +63,17 @@ setMethod( extractor(object@Model) } ) + + +#' @rdname extract +#' @aliases extract,n2kModelImputed-methods +#' @importFrom methods setMethod new +#' @include n2k_model_imputed_class.R +setMethod( + f = "extract", + signature = signature(object = "n2kModelImputed"), + definition = function(extractor, object, base = NULL, project = NULL) { + assert_that(inherits(extractor, "function")) + extractor(object) + } +) diff --git a/R/fit_model_character.R b/R/fit_model_character.R index 9decb4a3..578de151 100644 --- a/R/fit_model_character.R +++ b/R/fit_model_character.R @@ -48,7 +48,11 @@ setMethod( fit_model(base = base, project = project, verbose = verbose, ...) return(invisible(NULL)) } - if (!has_name(list(...), "local") || !inherits(base, "s3_bucket")) { + dots <- list(...) + if ( + !has_name(dots, "local") || is.null(dots$local) || + !inherits(base, "s3_bucket") + ) { analysis <- read_model(hash, base = base, project = project) display(verbose, paste(status(analysis), "-> "), FALSE) analysis <- fit_model( @@ -56,9 +60,12 @@ setMethod( ) display(verbose, status(analysis)) store_model(analysis, base = base, project = project) + result <- data.frame( + fingerprint = get_file_fingerprint(analysis), status = status(analysis) + ) rm(analysis) gc(verbose = FALSE) - return(invisible(NULL)) + return(invisible(result)) } dots <- list(...) to_do <- object_status( @@ -67,7 +74,10 @@ setMethod( if (length(to_do) == 0) { display(verbose, "skipping") gc(verbose = FALSE) - return(invisible(NULL)) + result <- data.frame( + fingerprint = hash, status = "converged" + ) + return(invisible(result)) } download_model( hash = hash, base = base, local = dots$local, project = project, @@ -91,8 +101,11 @@ setMethod( hash = hash, local = base, base = dots$local, project = project, verbose = verbose ) + result <- data.frame( + fingerprint = get_file_fingerprint(analysis), status = status(analysis) + ) rm(analysis) gc(verbose = FALSE) - return(invisible(NULL)) + return(invisible(result)) } ) diff --git a/R/fit_model_n2k_composite.R b/R/fit_model_n2k_composite.R index cd03f46c..1170cd11 100644 --- a/R/fit_model_n2k_composite.R +++ b/R/fit_model_n2k_composite.R @@ -1,6 +1,6 @@ #' @rdname fit_model #' @importFrom methods setMethod new -#' @importFrom dplyr %>% filter group_by n summarise transmute mutate arrange +#' @importFrom dplyr arrange filter group_by mutate n summarise transmute #' @importFrom rlang .data #' @importFrom utils file_test #' @importFrom stats qnorm @@ -21,21 +21,19 @@ setMethod( status(x) <- "error" return(x) } - x@Parameter %>% - filter(!is.na(.data$estimate), !is.na(.data$variance)) %>% - group_by(.data$value) %>% + x@Parameter |> + filter(!is.na(.data$estimate), !is.na(.data$variance)) |> + group_by(.data$value) |> summarise( - estimate = mean(.data$estimate), - se = sqrt(sum(.data$variance)) / n() - ) %>% + estimate = mean(.data$estimate), se = sqrt(sum(.data$variance)) / n() + ) |> transmute( - .data$value, - .data$estimate, + .data$value, .data$estimate, lower_confidence_limit = qnorm(0.025, mean = .data$estimate, sd = .data$se), upper_confidence_limit = qnorm(0.975, mean = .data$estimate, sd = .data$se) - ) %>% + ) |> as.data.frame() -> x@Index status(x) <- "converged" return(x) @@ -43,8 +41,8 @@ setMethod( status(x) <- "waiting" parent_status <- parent_status(x) - parent_status %>% - filter(.data$parent_status %in% c("new", "waiting", status)) %>% + parent_status |> + filter(.data$parent_status %in% c("new", "waiting", status)) |> pull("parent_analysis") -> todo for (this_parent in todo) { @@ -67,12 +65,12 @@ setMethod( ] <- get_status_fingerprint(model) x@AnalysisRelation <- parent_status if (status(model) == "converged") { - extract(extractor = x@Extractor, object = model) %>% - mutate(parent = this_parent) %>% + extract(extractor = x@Extractor, object = model) |> + mutate(parent = this_parent) |> bind_rows( - x@Parameter %>% + x@Parameter |> filter(.data$parent != this_parent) - ) %>% + ) |> arrange(.data$parent, .data$value) -> x@Parameter } } diff --git a/R/fit_model_n2k_hurdle_imputed.R b/R/fit_model_n2k_hurdle_imputed.R index 8635c822..edaf2571 100644 --- a/R/fit_model_n2k_hurdle_imputed.R +++ b/R/fit_model_n2k_hurdle_imputed.R @@ -18,7 +18,7 @@ setMethod( return(x) } - if (status(x) != "new") { + if (is.null(x@Presence) || status(x) %in% c("converged", "error")) { presence <- read_model( x@AnalysisRelation$parent_analysis[1], base = base, project = project ) @@ -29,7 +29,26 @@ setMethod( x@AnalysisRelation$parentstatus[1] <- status(presence) rm(presence) gc() - + x@AnalysisMetadata$status <- ifelse( + all(x@AnalysisRelation$parentstatus == "converged"), "new", + ifelse( + any(!x@AnalysisRelation$parentstatus %in% + c("new", "waiting", "converged")), + "error", "waiting" + ) + ) + x@AnalysisMetadata$status_fingerprint <- sha1( + list( + get_file_fingerprint(x), x@AnalysisMetadata$status, + x@AnalysisVersion$fingerprint, x@AnalysisVersion, x@RPackage, + x@AnalysisVersionRPackage, x@AnalysisRelation, x@Presence, x@Count, + x@Hurdle + ), + digits = 6L + ) + store_model(x, base = base, project = project) + } + if (is.null(x@Count) || status(x) %in% c("converged", "error")) { count <- read_model( x@AnalysisRelation$parent_analysis[2], base = base, project = project ) @@ -58,9 +77,10 @@ setMethod( ), digits = 6L ) + store_model(x, base = base, project = project) } - if (status(x) != "new") { + if (!all(x@AnalysisRelation$parentstatus == "converged")) { return(x) } diff --git a/R/fit_model_n2k_inla.R b/R/fit_model_n2k_inla.R index 1e34d955..2fe03d1b 100644 --- a/R/fit_model_n2k_inla.R +++ b/R/fit_model_n2k_inla.R @@ -1,6 +1,6 @@ #' @rdname fit_model #' @importFrom methods setMethod new -#' @importFrom assertthat assert_that is.number +#' @importFrom assertthat assert_that #' @importMethodsFrom multimput impute #' @include n2k_inla_class.R #' @param timeout the optional number of second until the model will time out @@ -40,43 +40,14 @@ setMethod( fm <- terms(x@AnalysisFormula[[1]]) response <- all.vars(fm)[attr(fm, "response")] if (mean(is.na(data[[response]])) < 0.10) { - # directly fit model when less than 10% missing data - control$data <- data - control$lincomb <- lc - model <- try({ - if (!is.null(timeout)) { - assert_that(is.number(timeout), timeout > 0) - setTimeLimit(cpu = timeout, elapsed = timeout) - } - do.call(INLA::inla, control) - }, silent = TRUE) + model <- direct_fit( + control = control, data = data, lc = lc, timeout = timeout + ) } else { - # first fit model without missing data - control$data <- data[!is.na(data[[response]]), ] - m0 <- try({ - if (!is.null(timeout)) { - assert_that(is.number(timeout), timeout > 0) - setTimeLimit(cpu = timeout, elapsed = timeout) - } - do.call(INLA::inla, control) - }, silent = TRUE) - if (inherits(m0, "try-error")) { - status(x) <- ifelse( - grepl("time limit", m0), "time-out", "error" - ) - return(x) - } - # then refit with missing data - control$data <- data - control$lincomb <- lc - control$control.update <- list(result = m0) - model <- try({ - if (!is.null(timeout)) { - assert_that(is.number(timeout), timeout > 0) - setTimeLimit(cpu = timeout, elapsed = timeout) - } - do.call(INLA::inla, control) - }, silent = TRUE) + model <- indirect_fit( + response = response, control = control, data = data, lc = lc, + timeout = timeout + ) } # handle error in model fit @@ -126,3 +97,56 @@ model2lincomb <- function(lincomb) { } return(lc) } + +# directly fit model when less than 10% missing data +direct_fit <- function(control, data, lc, timeout = NULL) { + control$data <- data + control$lincomb <- lc + try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) +} + +#' @importFrom assertthat assert_that is.number +indirect_fit <- function(control, data, lc, response, timeout = NULL) { + # first fit model without missing data + control$data <- data[!is.na(data[[response]]), ] + m0 <- try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) + if (inherits(m0, "try-error") && "control.family" %in% names(control)) { + # when model failed to fit, try again without family + old_control_family <- control$control.family + control$control.family <- NULL + m0 <- try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) + control$control.family <- old_control_family + } + if (inherits(m0, "try-error")) { + return(m0) + } + # then refit with missing data + control$data <- data + control$lincomb <- lc + control$control.update <- list(result = m0) + try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) +} diff --git a/R/fit_model_n2k_manifest.R b/R/fit_model_n2k_manifest.R index 3b044477..c291f3e9 100644 --- a/R/fit_model_n2k_manifest.R +++ b/R/fit_model_n2k_manifest.R @@ -4,92 +4,65 @@ #' @importFrom methods setMethod new #' @importFrom purrr walk #' @importFrom stats na.omit +#' @importFrom tools R_user_dir #' @include n2k_manifest_class.R #' @param local A local folder into which objects from an AWS S3 bucket are #' downloaded. -#' @param first A logical. -#' `first = TRUE` implies to fit only the first object in the manifest with -#' matching status. -#' `first = FALSE` implies to fit all objects in the manifest with matching -#' status. -#' Defaults to `FALSE`. setMethod( f = "fit_model", signature = signature(x = "n2kManifest"), definition = function( x, base, project, status = c("new", "waiting"), verbose = TRUE, ..., - local = tempfile("fit_model"), first = FALSE + local = NULL ) { assert_that( is.string(project), noNA(project), is.character(status), noNA(status), length(status) >= 1 ) to_do <- order_manifest(x) - remaining <- length(to_do) - while (length(to_do) > 1 && first) { - head(to_do, 1) |> - hash_status(base = base, project = project) -> stat - if (stat %in% status) { - to_do <- head(to_do, 1) - } else { - to_do <- tail(to_do, -1) - } - } - if (length(to_do) == 0) { - return(invisible(0)) + R_user_dir("n2kanalysis", which = "cache") |> + file.path(x@Fingerprint) -> cache_file + dirname(cache_file) |> + dir.create(recursive = TRUE, showWarnings = FALSE) + if (file_test("-f", cache_file)) { + processed <- read.table(cache_file, header = TRUE, sep = "\t") + done <- processed$status %in% status + to_do <- to_do[!to_do %in% processed$fingerprint[!done]] + } else { + data.frame(fingerprint = character(0), status = character(0)) |> + write.table( + file = cache_file, sep = "\t", row.names = FALSE, quote = FALSE + ) } - if (inherits(base, "character")) { - walk( - to_do, fit_model, base = base, project = project, - status = status, verbose = verbose, ... + start_time <- Sys.time() + for (i in seq_along(to_do)) { + display( + verbose = verbose, + message = sprintf( + "Processing %i from %i (%.2f%%) %s ETA %s", i, length(to_do), + 100 * (i - 1) / length(to_do), + format(Sys.time(), "%Y-%m-%d %H:%M:%S"), + format( + start_time + (Sys.time() - start_time) * length(to_do) / (i - 1), + "%d %H:%M" + ) + ) ) - return(invisible(remaining)) + result <- try(fit_model( + x = to_do[i], base = base, project = project, status = status, + verbose = verbose, ..., local = local + )) + if (!inherits(result, "try-error")) { + write.table( + result, file = cache_file, append = TRUE, sep = "\t", + row.names = FALSE, quote = FALSE, col.names = FALSE + ) + } } - - display(verbose, "Downloading objects") - x@Manifest$parent[x@Manifest$fingerprint %in% to_do] |> - c(to_do) |> - unique() |> - na.omit() -> to_download - path(local, project) |> - dir_create() - path(local, project) |> - dir_ls(recurse = TRUE, type = "file") |> - basename() -> local_files - to_download[!paste0(to_download, ".rds") %in% local_files] |> - walk( - download_model, base = base, project = project, local = local, - verbose = verbose - ) - walk( - to_do, fit_model, base = local, project = project, - status = status, verbose = verbose, ... - ) - display(verbose, "Uploading objects") - walk( - to_do, download_model, base = local, project = project, local = base, - verbose = verbose - ) - return(invisible(remaining)) + return(invisible(NULL)) } ) -#' @importFrom aws.s3 get_bucket -#' @importFrom purrr map_chr -hash_status <- function(hash, base, project) { - if (inherits(base, "s3_bucket")) { - substr(hash, 1, 4) |> - sprintf(fmt = "%2$s/%1$s/", project) |> - get_bucket(bucket = base, max = Inf) |> - map_chr("Key") -> keys - keys[grepl(hash, keys)] |> - gsub(pattern = sprintf(".*/(.*)/%s\\.rds", hash), replacement = "\\1") |> - unname() -> output - return(output) - } - stop("hash status for ", class(base), " still do to") -} - download_model <- function(hash, base, local, project, verbose = FALSE) { display(verbose, paste("Moving", hash)) read_model(x = hash, base = base, project = project) |> diff --git a/R/fit_model_n2k_spde.R b/R/fit_model_n2k_spde.R new file mode 100644 index 00000000..2c0816aa --- /dev/null +++ b/R/fit_model_n2k_spde.R @@ -0,0 +1,158 @@ +#' @rdname fit_model +#' @importFrom methods setMethod new +#' @importFrom assertthat assert_that is.number +#' @importMethodsFrom multimput impute +#' @include n2k_spde_class.R +#' @param timeout the optional number of second until the model will time out +#' @inheritParams multimput::impute +setMethod( + f = "fit_model", + signature = signature(x = "n2kSpde"), + definition = function( + x, status = "new", ..., timeout = NULL, seed = get_seed(x), + num_threads = NULL, parallel_configs = TRUE + ) { + assert_that( + requireNamespace("INLA", quietly = TRUE), + msg = "INLA package required but not installed." + ) + validObject(x) + assert_that(is.character(status), length(status) >= 1) + + # don't fit model when status doesn't match + if (!(status(x) %in% status)) { + return(x) + } + + set.seed(seed) + + data <- get_data(x) + matern <- spde2matern(x@Spde) + index <- INLA::inla.spde.make.index("matern_spde", n.spde = matern$n.spde) + x@AnalysisMetadata$formula |> + c("f(matern_spde, model = matern)") |> + paste(collapse = " + ") |> + as.formula() -> model_formula + fm <- terms(model_formula) + response <- all.vars(fm)[attr(fm, "response")] + stack_observed <- INLA::inla.stack( + data = data[!is.na(data[[response]]), response, drop = FALSE], + A = list(1, make_a(object = x@Spde, data[!is.na(data[[response]]), ])), + tag = "observed", + effects = list( + data[ + !is.na(data[[response]]), colnames(data) != response, drop = FALSE + ], + index + ) + ) + INLA::inla.stack( + data = data[, response, drop = FALSE], + A = list(1, make_a(object = x@Spde, data)), + tag = "observed", + effects = list( + data[, colnames(data) != response, drop = FALSE], + index + ) + ) -> stack_total + # prepare linear combinations + lc <- model2lincomb(x@LinearCombination) + # prepare inla() arguments + control <- x@Control + control$formula <- model_formula + control$family <- x@Family + # fit model + if (mean(is.na(data[[response]])) < 0.10) { + # directly fit model when less than 10% missing data + control$data <- INLA::inla.stack.data(stack_total, spde = spde) + control$control.predictor <- list( + A = INLA::inla.stack.A(stack_total), link = 1 + ) + control$lincomb <- lc + model <- try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) + } else { + # first fit model without missing data + control$data <- INLA::inla.stack.data(stack_observed, spde = spde) + control$control.predictor <- list(A = INLA::inla.stack.A(stack_observed)) + m0 <- try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) + if (inherits(m0, "try-error")) { + status(x) <- ifelse( + grepl("time limit", m0), "time-out", "error" + ) + return(x) + } + # then refit with missing data + control$data <- INLA::inla.stack.data(stack_total, spde = spde) + control$control.predictor <- list( + A = INLA::inla.stack.A(stack_total), link = 1 + ) + control$lincomb <- lc + control$control.update <- list(result = m0) + model <- try({ + if (!is.null(timeout)) { + assert_that(is.number(timeout), timeout > 0) + setTimeLimit(cpu = timeout, elapsed = timeout) + } + do.call(INLA::inla, control) + }, silent = TRUE) + } + + # handle error in model fit + if (inherits(model, "try-error")) { + status(x) <- ifelse( + grepl("time limit", model), "time-out", "error" + ) + return(x) + } + # return fitted model when no imputation is required + if (x@ImputationSize == 0) { + return(n2k_spde(data = x, model_fit = model, status = "converged")) + } + + imputed <- try(impute( + model = model, n_imp = x@ImputationSize, minimum = x@Minimum, + seed = seed, num_threads = num_threads, extra = x@Extra, + parallel_configs = parallel_configs + )) + if (inherits(imputed, "try-error")) { + return(n2k_spde(data = x, model_fit = model, status = "error")) + } + # return fitted model with imputations + return(n2k_spde( + data = x, model_fit = model, status = "converged", raw_imputed = imputed + )) + } +) + +model2lincomb <- function(lincomb) { + if (is.null(lincomb)) { + return(NULL) + } + if (inherits(lincomb, "matrix")) { + lincomb |> + as.data.frame() |> + as.list() %>% + INLA::inla.make.lincombs() |> + setNames(rownames(lincomb)) -> lc + return(lc) + } + lc <- INLA::inla.make.lincombs(lincomb) + if (is.matrix(lincomb[[1]])) { + names(lc) <- rownames(lincomb[[1]]) + } else { + names(lc) <- names(lincomb[[1]]) + } + return(lc) +} diff --git a/R/get_datafield_id.R b/R/get_datafield_id.R new file mode 100644 index 00000000..6836271d --- /dev/null +++ b/R/get_datafield_id.R @@ -0,0 +1,49 @@ +#' Get the data field id +#' @param table The table name +#' @param field The field name +#' @param datasource The data source name +#' @inheritParams git2rdata::write_vc +#' @export +#' @importFrom assertthat assert_that is.string noNA +#' @importFrom git2rdata is_git2rmeta update_metadata verify_vc write_vc +get_datafield_id <- function(table, field, datasource, root, stage = FALSE) { + assert_that( + is.string(table), is.string(field), is.string(datasource), noNA(table), + noNA(field), noNA(datasource) + ) + if (!is_git2rmeta(file = "datafield", root = root)) { + data.frame(id = 1L, table = table, field = field, source = datasource) |> + write_vc(file = "datafield", root = root, sorting = "id", stage = stage) + update_metadata( + file = "datafield", root = root, stage = stage, name = "datafield", + title = "Pointer to external code identifiers", + description = + "This dataset describes the external code identifiers used in the data. +It points to the original source of the external code: which datasource, which +table in that datasource and which field in that table.", + field_description = c( + id = "unique identifier of the datafield", + table = "table name of the identifier", + field = "field name of the identifier", + source = "data source which stores the data" + ) + ) + return(1L) + } + datafield <- verify_vc( + file = "datafield", root = root, + variables = c("id", "table", "field", "source") + ) + which( + datafield$table == table & datafield$field == field & + datafield$source == datasource + ) -> relevant + stopifnot("multiple matching datafield id found" = length(relevant) <= 1) + if (length(relevant) == 1) { + return(datafield$id[relevant]) + } + new_id <- max(datafield$id) + 1L + data.frame(id = new_id, table = table, field = field, source = datasource) |> + write_vc(file = "datafield", root = root, stage = stage, append = TRUE) + return(new_id) +} diff --git a/R/make_a.R b/R/make_a.R new file mode 100644 index 00000000..b0216e1c --- /dev/null +++ b/R/make_a.R @@ -0,0 +1,36 @@ +#' Convert an `Spde` object to a A object +#' @param object The `Spde` object +#' @param data The `data.frame` to combine with the `Spde` object to create the +#' `A` matrix. +#' @name make_a +#' @rdname make_a +#' @exportMethod make_a +#' @docType methods +#' @importFrom methods setGeneric +setGeneric( + name = "make_a", + def = function(object, data) { + standardGeneric("make_a") # nocov + } +) + +#' @rdname make_a +#' @importFrom methods setMethod new +#' @include spde_class.R +setMethod( + f = "make_a", + signature = signature(object = "Spde"), + definition = function(object, data) { + assert_that( + inherits(data, "data.frame"), + all(colnames(object@Coordinates) %in% colnames(data)) + ) + stopifnot( + "INLA package required but not installed." = + requireNamespace("INLA", quietly = TRUE) + ) + data[colnames(object@Coordinates)] |> + as.matrix() |> + INLA::inla.spde.make.A(mesh = spde2mesh(object)) + } +) diff --git a/R/manifest_yaml_to_bash.R b/R/manifest_yaml_to_bash.R index c08298d1..fb40a147 100644 --- a/R/manifest_yaml_to_bash.R +++ b/R/manifest_yaml_to_bash.R @@ -67,42 +67,16 @@ setMethod( gsub("\\.manifest$", "", yaml$hash) |> read_manifest(base = base, project = project) -> manifest docker_hash <- get_file_fingerprint(manifest) - sprintf( - "RUN Rscript -e 'remotes::install_github(\\\"%s\\\"%s)'", yaml$github, - ", dependencies = FALSE, upgrade = \\\"never\\\"" - ) -> deps - sprintf( - "#!/bin/bash -export $(cat .env | xargs) -echo \"FROM %s -%s\" > Dockerfile -docker build --pull --tag rn2k:%s . -rm Dockerfile", - yaml$docker, paste(deps, collapse = "\n"), docker_hash - ) -> init + init <- create_docker_init(yaml = yaml, docker_hash = docker_hash) volume <- "/n2kanalysis:/n2kanalysis:rw" models <- order_manifest(manifest = manifest) to_do <- object_status(base = base, project = project, status = status) models <- models[models %in% to_do] - c( - "echo \"\n\nmodel %i of %i\n\n\"\ndate\n", - "timeout --kill-after=2m %ih docker run %s --name=%s -v %s rn2k:%s", - "./fit_model_aws.sh -b %s -p %s -m %s%s" - ) |> - paste(collapse = " ") |> - sprintf( - seq_along(models), length(models), timeout, - paste( - c( - "--rm", "--env AWS_ACCESS_KEY_ID=$AWS_ACCESS_KEY_ID", - "--env AWS_SECRET_ACCESS_KEY=$AWS_SECRET_ACCESS_KEY", - "--env AWS_DEFAULT_REGION=$AWS_DEFAULT_REGION", - "--cap-add NET_ADMIN"[limit], "--cpu-shares=512"[limit] - ), - collapse = " " - ), models, volume, docker_hash, get_bucketname(base), project, models, - ifelse(limit, " -s 1", "") - ) -> model_scripts + model_scripts <- create_docker_model_scripts( + models = models, base = base, timeout = timeout, + limit = limit, volume = volume, docker_hash = docker_hash, + project = project + ) vapply( seq_len(split), FUN.VALUE = character(1), project = project, init = init, split = split, shutdown = shutdown, base = base, @@ -121,6 +95,52 @@ rm Dockerfile", } ) +create_docker_init <- function(yaml, docker_hash) { + sprintf( + "RUN Rscript -e 'pak::pkg_install(\\\"%s\\\"%s)'", yaml$github, + ", dependencies = FALSE, upgrade = FALSE, ask = FALSE" + ) -> deps + sprintf( + "#!/bin/bash +export $(cat .env | xargs) +echo \"FROM %s +%s\" > Dockerfile +docker build --pull --tag rn2k:%s . +rm Dockerfile", + yaml$docker, paste(deps, collapse = "\n"), docker_hash + ) +} + +create_docker_model_scripts <- function( + models, base, timeout = 4, limit = FALSE, volume, docker_hash, project +) { + if (inherits(base, "character")) { + script <- "./fit_model_file.sh" + } else { + script <- "./fit_model_aws.sh" + base <- get_bucketname(base) + } + c( + "echo \"\n\nmodel %i of %i\n\n\"\ndate\n", + "timeout --kill-after=2m %ih docker run %s --name=%s -v %s rn2k:%s", + script, " -b %s -p %s -m %s%s" + ) |> + paste(collapse = " ") |> + sprintf( + seq_along(models), length(models), timeout, + paste( + c( + "--rm", "--env AWS_ACCESS_KEY_ID=$AWS_ACCESS_KEY_ID", + "--env AWS_SECRET_ACCESS_KEY=$AWS_SECRET_ACCESS_KEY", + "--env AWS_DEFAULT_REGION=$AWS_DEFAULT_REGION", + "--cap-add NET_ADMIN"[limit], "--cpu-shares=512"[limit] + ), + collapse = " " + ), models, volume, docker_hash, base, project, models, + ifelse(limit, " -s 1", "") + ) +} + #' @export #' @rdname manifest_yaml_to_bash #' @importFrom assertthat assert_that is.count is.flag is.string noNA @@ -134,7 +154,7 @@ setMethod( signature = signature(base = "character"), definition = function( base, project, hash, shutdown = FALSE, split = 1, - status = c("new", "waiting"), limit = FALSE + status = c("new", "waiting"), limit = FALSE, timeout = 4 ) { assert_that( is.string(base), noNA(base), file_test("-d", base), is.string(project), @@ -143,13 +163,11 @@ setMethod( ) assert_that(split == 1, msg = "`split > 1` to do on local file systems.") assert_that( - file_test("-d", path(base, project)), - msg = sprintf("`%s` is not a subdirectory of `%s`", project, base) - ) - assert_that( - file_test("-d", path(base, project, "yaml")), - msg = sprintf("`yaml` is not a subdirectory of `%s/%s`", base, project) + file_test("-d", base), + msg = sprintf("`%s` is not an existing folder", base) ) + path(base, project, "yaml") |> + dir.create(showWarnings = FALSE) if (missing(hash)) { path(base, project, "yaml") |> dir_info(type = "file", regexp = "\\.yaml$") |> @@ -171,31 +189,14 @@ setMethod( gsub("\\.manifest$", "", yaml$hash) |> read_manifest(base = base, project = project) -> manifest docker_hash <- get_file_fingerprint(manifest) - sprintf( - "Rscript -e 'remotes::install_github(\\\"%s\\\"%s)'", yaml$github, - ", dependencies = TRUE, upgrade = \\\"never\\\", keep_source = FALSE" - ) -> deps - sprintf( - "#!/bin/bash -echo \"FROM %s -RUN %s\" > Dockerfile -docker build --pull --tag rn2k:%s . -rm Dockerfile", - yaml$docker, paste(deps, collapse = " \\\n&& "), docker_hash - ) -> init + init <- create_docker_init(yaml = yaml, docker_hash = docker_hash) base <- normalizePath(base, winslash = "/") volume <- paste(base, base, "rw", sep = ":") models <- order_manifest(manifest = manifest) - sprintf( - "echo \"model %i of %i\" -docker run %s%s --name=%s -v %s rn2k:%s ./fit_model_file.sh -b %s -p %s -m %s -date -docker stop --time 14400 %s -date", - seq_along(models), length(models), "--rm -d", - ifelse(limit, "--cpu-shares=512", ""), models, volume, docker_hash, - base, project, models, models - ) -> model_scripts + model_scripts <- create_docker_model_scripts( + models = models, base = base, timeout = timeout, limit = limit, + volume = volume, docker_hash = docker_hash, project = project + ) path(base, project, "bash") |> dir_create() script <- path(base, project, sprintf("bash/%s.sh", docker_hash)) @@ -206,7 +207,15 @@ date", } ) +#' Convert a manifest in vector of analysis fingerprints +#' +#' The order of the analysis is determined by the parent-child relationship. +#' It starts with the analyses without parents. +#' Then it adds the analyses with parents that have already been added. +#' This process is repeated until all analyses have been added. +#' @param manifest the `n2kManifest` #' @importFrom assertthat assert_that +#' @export order_manifest <- function(manifest) { assert_that(inherits(manifest, "n2kManifest")) full <- slot(manifest, "Manifest") diff --git a/R/moving_average.R b/R/moving_average.R new file mode 100644 index 00000000..57f3d843 --- /dev/null +++ b/R/moving_average.R @@ -0,0 +1,32 @@ +#' Calculate moving trend coefficients +#' @inheritParams moving_trend +#' @return A matrix with moving average coefficients +#' One row for each window and and one column for each year in the data. +#' The format of the row names is `average_{window mid point}_{window length}`. +#' `average_2001.5_4` is the average for the years 2000 to 2003. +#' @export +#' @examples +#' moving_average(5, 3) +#' moving_average(5, 3, 2000) +#' @importFrom assertthat assert_that is.count is.number +#' @importFrom stats median +moving_average <- function(n_year, duration, first_year = 0) { + assert_that(is.count(n_year), is.count(duration), is.number(first_year)) + duration <- min(n_year, duration) + vapply( + seq_len(n_year - duration + 1) - 1, + FUN.VALUE = vector(mode = "numeric", length = n_year), + FUN = function(i, trend_coef, n_year) { + c(rep(0, i), trend_coef, rep(0, n_year - length(trend_coef) - i)) + }, trend_coef = rep(1 / duration, duration), n_year = n_year + ) |> + `colnames<-`( + sprintf( + "average_%.1f_%i", + seq_len(n_year - duration + 1) + median(seq_len(duration)) - 2 + + first_year, + duration + ) + ) |> + t() +} diff --git a/R/moving_difference.R b/R/moving_difference.R new file mode 100644 index 00000000..33b893d4 --- /dev/null +++ b/R/moving_difference.R @@ -0,0 +1,45 @@ +#' Calculate coefficients for a moving difference +#' @inheritParams moving_trend +#' @return A matrix with moving average coefficients +#' One row for each window and and one column for each year in the data. +#' The format of the row names is +#' `difference_{window mid point start}_{window mid point end}_{window length}`. +#' `difference_2001.5_2010.5_4` is the difference of the average for period for +#' the years 2009 to 2012 compared to the period from 2000 to 2003. +#' @export +#' @examples +#' moving_difference(6, 2) +#' moving_difference(6, 2, 2000) +#' @importFrom assertthat assert_that is.count is.number +#' @importFrom stats median +moving_difference <- function(n_year, duration, first_year = 1) { + assert_that(is.count(n_year), is.count(duration), is.number(first_year)) + duration <- min(floor(n_year / 2), duration) + list(seq_len(n_year - 2 * duration + 1) - 1) |> + rep(2) |> + expand.grid() -> extra_zero + extra_zero <- extra_zero[ + rowSums(extra_zero) <= n_year - 2 * duration, + ] + vapply( + seq_len(nrow(extra_zero)), + FUN.VALUE = vector(mode = "numeric", length = n_year), + FUN = function(i, trend_coef, n_year, extra_zero) { + c( + rep(0, extra_zero[i, 1]), -trend_coef, + rep(0, n_year - 2 * length(trend_coef) - sum(extra_zero[i, ])), + trend_coef, rep(0, extra_zero[i, 2]) + ) + }, trend_coef = rep(1 / duration, duration), n_year = n_year, + extra_zero = extra_zero + ) |> + `colnames<-`( + sprintf( + "difference_%.1f_%.1f_%i", + extra_zero[, 1] + median(seq_len(duration)) - 1 + first_year, + n_year - extra_zero[, 2] - median(seq_len(duration)) + first_year, + duration + ) + ) |> + t() +} diff --git a/R/moving_trend.R b/R/moving_trend.R new file mode 100644 index 00000000..8cd35779 --- /dev/null +++ b/R/moving_trend.R @@ -0,0 +1,39 @@ +#' Calculate moving trend coefficients +#' @param n_year Number of available years in the data. +#' @param duration Number of years in the moving window. +#' If the number of available years is less than this value, the trend will be +#' calculated for the available years. +#' @param first_year First year of the data. +#' Only used to name the rows. +#' @return A matrix with the moving trend coefficients. +#' One row for each window and and one column for each year in the data. +#' The format of the row names is `trend_{window mid point}_{window length}`. +#' `trend_2001.5_4` is the trend for the years 2000 to 2003. +#' @export +#' @examples +#' moving_trend(5, 3) +#' moving_trend(5, 4, 2000) +#' @export +#' @importFrom assertthat assert_that is.count is.number +#' @importFrom stats median +moving_trend <- function(n_year, duration, first_year = 0) { + assert_that(is.count(n_year), is.count(duration), is.number(first_year)) + duration <- min(n_year, duration) + trend_coef <- seq_len(duration) - (duration + 1) / 2 + vapply( + seq_len(n_year - duration + 1), + function(i, trend_coef) { + c(rep(0, i - 1), trend_coef, rep(0, n_year - duration - i + 1)) + }, + numeric(n_year), trend_coef = trend_coef / sum(trend_coef ^ 2) + ) |> + `colnames<-`( + sprintf( + "trend_%.1f_%i", + seq_len(n_year - duration + 1) + median(seq_len(duration)) - 2 + + first_year, + duration + ) + ) |> + t() +} diff --git a/R/n2k_hurdle_imputed_class.R b/R/n2k_hurdle_imputed_class.R index 56a85b1f..5ff80a00 100644 --- a/R/n2k_hurdle_imputed_class.R +++ b/R/n2k_hurdle_imputed_class.R @@ -10,6 +10,7 @@ #' @aliases n2kHurdleImputed-class #' @importFrom methods setClass #' @docType class +#' @include n2k_aggregate_class.R setClass( "n2kHurdleImputed", representation = representation( diff --git a/R/n2k_model_imputed_class.R b/R/n2k_model_imputed_class.R index 25973ca1..c5c5d905 100644 --- a/R/n2k_model_imputed_class.R +++ b/R/n2k_model_imputed_class.R @@ -31,7 +31,7 @@ setClassUnion("maybeFunction", c("function", "character")) #' @aliases n2kModelImputed-class #' @importFrom methods setClass #' @docType class -#' @include n2k_model_imputed_class.R +#' @include n2k_aggregate_class.R #' @include n2k_inla_comparison_class.R setClass( "n2kModelImputed", diff --git a/R/n2k_spde.R b/R/n2k_spde.R new file mode 100644 index 00000000..cc7b9609 --- /dev/null +++ b/R/n2k_spde.R @@ -0,0 +1,220 @@ +#' Create an `n2kSpde` object +#' @param data a `data.frame` with the data to analyse +#' @param model_fit The fitted model +#' @param ... other arguments. See below +#' @name n2k_spde +#' @rdname n2k_spde +#' @exportMethod n2k_spde +#' @docType methods +#' @importFrom methods setGeneric +setGeneric( + name = "n2k_spde", + def = function( + data, ..., model_fit + ) { + standardGeneric("n2k_spde") # nocov + } +) + +#' @description A new `n2kSpde` model is created when `data` is a `data.frame`. +#' @param spde the `Spde` object. See `[spde]`. +#' @rdname n2k_spde +#' @aliases n2k_spde,n2kSpde-methods +#' @importFrom methods setMethod new +#' @importFrom assertthat assert_that is.count is.string is.time +#' @importFrom digest sha1 +#' @importFrom stats as.formula +#' @importFrom utils sessionInfo +#' @include n2k_spde_class.R +#' @inheritParams n2k_inla_comparison +#' @inheritParams n2k_inla +#' @inheritParams multimput::impute +setMethod( + f = "n2k_spde", + signature = signature(data = "data.frame"), + definition = function( + data, status = "new", result_datasource_id, scheme_id, family = "poisson", + formula, species_group_id, location_group_id, model_type, spde, + first_imported_year, last_imported_year, duration, last_analysed_year, + analysis_date, lin_comb = NULL, minimum = "", imputation_size, + parent = character(0), seed, replicate_name = list(), control = list(), + parent_status = "converged", parent_statusfingerprint, extra, ..., model_fit + ) { + assert_that(is.string(status)) + assert_that(is.string(minimum)) + if (missing(seed)) { + seed <- sample(.Machine$integer.max, 1) + } + assert_that(is.count(seed)) + seed <- as.integer(seed) + if (missing(imputation_size)) { + imputation_size <- 0L + } else { + assert_that(is.count(imputation_size)) + imputation_size <- as.integer(imputation_size) + } + assert_that( + is.string(result_datasource_id), is.string(scheme_id), + is.string(species_group_id), is.string(location_group_id), + is.string(model_type), is.string(formula), is.count(first_imported_year), + is.count(last_imported_year) + ) + first_imported_year <- as.integer(first_imported_year) + last_imported_year <- as.integer(last_imported_year) + if (missing(duration)) { + duration <- last_imported_year - first_imported_year + 1L + } else { + assert_that(is.count(duration)) + duration <- as.integer(duration) + } + if (missing(last_analysed_year)) { + last_analysed_year <- last_imported_year + } + assert_that(is.count(last_analysed_year)) + last_analysed_year <- as.integer(last_analysed_year) + assert_that(is.time(analysis_date)) + assert_that( + is.null(lin_comb) || inherits(lin_comb, "list") || + (inherits(lin_comb, "matrix") && length(dim(lin_comb) == 2)), + msg = "lin_comb must be either a list or a matrix" + ) + assert_that(is.list(replicate_name)) + assert_that( + length(replicate_name) == 0 || !is.null(names(replicate_name)), + msg = "replicate_name must have names" + ) + assert_that(is.character(family), length(family) >= 1) + assert_that(is.list(control)) + control$control.compute$dic <- ifelse( + is.null(control$control.compute$dic), TRUE, control$control.compute$dic + ) + control$control.compute$waic <- ifelse( + is.null(control$control.compute$waic), TRUE, control$control.compute$waic + ) + control$control.compute$cpo <- ifelse( + is.null(control$control.compute$cpo), TRUE, control$control.compute$cpo + ) + control$control.compute$config <- ifelse( + is.null(control$control.compute$config), TRUE, + control$control.compute$config + ) + control$control.predictor$compute <- ifelse( + is.null(control$control.predictor$compute), TRUE, + control$control.predictor$compute + ) + if (is.null(control$control.predictor$link)) { + control$control.predictor$link <- 1 + } + control$control.fixed$prec.intercept <- ifelse( + is.null(control$control.fixed$prec.intercept), + 1, control$control.fixed$prec.intercept + ) + if (missing(extra)) { + extra <- data[0, ] + } + + file_fingerprint <- sha1( + list( + data, result_datasource_id, scheme_id, species_group_id, + location_group_id, family, model_type, formula, first_imported_year, + last_imported_year, duration, last_analysed_year, + format(analysis_date, tz = "UTC"), seed, parent, replicate_name, + lin_comb, imputation_size, minimum, control, extra + ) + ) + + if (length(parent) == 0) { + analysis_relation <- data.frame( + analysis = character(0), parent_analysis = character(0), + parentstatus_fingerprint = character(0), parent_status = character(0), + stringsAsFactors = FALSE + ) + } else { + assert_that(is.string(parent)) + assert_that(is.string(parent_status)) + if (missing(parent_statusfingerprint)) { + parent_statusfingerprint <- sha1(parent_status) + } else { + assert_that(is.string(parent_statusfingerprint)) + } + analysis_relation <- data.frame( + analysis = file_fingerprint, parent_analysis = parent, + parentstatus_fingerprint = parent_statusfingerprint, + parent_status = parent_status, stringsAsFactors = FALSE + ) + } + version <- get_analysis_version(sessionInfo()) + status_fingerprint <- sha1( + list( + file_fingerprint, status, NULL, version@AnalysisVersion$fingerprint, + version@AnalysisVersion, version@RPackage, + version@AnalysisVersionRPackage, analysis_relation, NULL + ), + digits = 6L + ) + + new( + "n2kSpde", + AnalysisVersion = version@AnalysisVersion, RPackage = version@RPackage, + AnalysisVersionRPackage = version@AnalysisVersionRPackage, + AnalysisMetadata = data.frame( + result_datasource_id = result_datasource_id, scheme_id = scheme_id, + species_group_id = species_group_id, + location_group_id = location_group_id, model_type = model_type, + formula = formula, first_imported_year = first_imported_year, + last_imported_year = last_imported_year, duration = duration, + last_analysed_year = last_analysed_year, analysis_date = analysis_date, + seed = seed, status = status, + analysis_version = version@AnalysisVersion$fingerprint, + file_fingerprint = file_fingerprint, + status_fingerprint = status_fingerprint, stringsAsFactors = FALSE + ), + AnalysisFormula = list(as.formula(formula)), LinearCombination = lin_comb, + AnalysisRelation = analysis_relation, Data = data, Model = NULL, + ReplicateName = replicate_name, Family = family, Control = control, + ImputationSize = imputation_size, Minimum = minimum, RawImputed = NULL, + Extra = extra, Spde = spde + ) + } +) + +#' @description In case `data` is an `n2kSpde` object, then only the model and +#' status are updated. +#' All other slots are unaffected. +#' @rdname n2k_spde +#' @aliases n2k_spde,n2kSpde-methods +#' @importFrom methods setMethod validObject new +#' @importFrom digest sha1 +#' @importFrom utils sessionInfo +#' @include n2k_spde_class.R +#' @param raw_imputed the optional `rawImputed` object +setMethod( + f = "n2k_spde", + signature = signature(data = "n2kSpde", model_fit = "inla"), + definition = function( + data, status, raw_imputed = NULL, ..., model_fit + ) { + assert_that(is.string(status)) + data@Model <- model_fit + data@AnalysisMetadata$status <- status + version <- get_analysis_version(sessionInfo()) + new_version <- union(data, version) + data@AnalysisVersion <- new_version$Union@AnalysisVersion + data@RPackage <- new_version$Union@RPackage + data@AnalysisVersionRPackage <- new_version$Union@AnalysisVersionRPackage + data@AnalysisMetadata$analysis_version <- new_version$Unionfingerprint + data@RawImputed <- raw_imputed + data@AnalysisMetadata$status_fingerprint <- sha1( + list( + data@AnalysisMetadata$file_fingerprint, data@AnalysisMetadata$status, + data@Model, data@AnalysisMetadata$analysis_version, + data@AnalysisVersion, data@RPackage, data@AnalysisVersionRPackage, + data@AnalysisRelation, data@RawImputed + ), + digits = 6L + ) + + validObject(data) + return(data) + } +) diff --git a/R/n2k_spde_class.R b/R/n2k_spde_class.R new file mode 100644 index 00000000..10b8d6e8 --- /dev/null +++ b/R/n2k_spde_class.R @@ -0,0 +1,26 @@ +#' The `n2kSpde` class +#' +#' It hold analysis data based on an INLA model with SPDE. +#' @slot Spde A `list` containing the information for the SPDE. +#' @name n2kSpde-class +#' @rdname n2kSpde-class +#' @exportClass n2kSpde +#' @aliases n2kSpde-class +#' @importFrom methods setClass +#' @docType class +#' @include n2k_inla_class.R +#' @include spde_class.R +setClass( + "n2kSpde", + representation = representation(Spde = "Spde"), + contains = "n2kInla" +) + +#' @importFrom methods setValidity +#' @importFrom digest sha1 +setValidity( + "n2kSpde", + function(object) { + return(TRUE) + } +) diff --git a/R/spde.R b/R/spde.R new file mode 100644 index 00000000..ce8738b6 --- /dev/null +++ b/R/spde.R @@ -0,0 +1,18 @@ +#' Create an `Spde` object +#' @param coordinates a `data.frame` of coordinates use to define the mesh. +#' @param range a numeric vector of length 2. +#' Will be used as the `prior.range` argument of `[INLA::inla.spde2.pcmatern]`. +#' @param sigma a numeric vector of length 2. +#' Will be used as the `prior.sigma` argument of `[INLA::inla.spde2.pcmatern]`. +#' @export +#' @importFrom assertthat assert_that noNA +#' @importFrom methods new +spde <- function(coordinates, range, sigma) { + assert_that( + inherits(coordinates, "data.frame"), is.numeric(range), is.numeric(sigma), + noNA(coordinates), noNA(range), noNA(sigma), ncol(coordinates) == 2, + length(range) == 2, length(sigma) == 2, all(range > 0), all(sigma > 0), + range[2] < 1, sigma[2] < 1 + ) + new("Spde", Coordinates = coordinates, Range = range, Sigma = sigma) +} diff --git a/R/spde2matern.R b/R/spde2matern.R new file mode 100644 index 00000000..7946eb58 --- /dev/null +++ b/R/spde2matern.R @@ -0,0 +1,31 @@ +#' Convert an `Spde` object to a `matern` object +#' @param object The `Spde` object +#' @name spde2matern +#' @rdname spde2matern +#' @exportMethod spde2matern +#' @docType methods +#' @importFrom methods setGeneric +setGeneric( + name = "spde2matern", + def = function(object) { + standardGeneric("spde2matern") # nocov + } +) + +#' @rdname spde2mesh +#' @importFrom methods setMethod new +#' @include spde_class.R +setMethod( + f = "spde2matern", + signature = signature(object = "Spde"), + definition = function(object) { + stopifnot( + "INLA package required but not installed." = + requireNamespace("INLA", quietly = TRUE) + ) + INLA::inla.spde2.pcmatern( + mesh = spde2mesh(object), prior.range = object@Range, + prior.sigma = object@Sigma + ) + } +) diff --git a/R/spde2mesh.R b/R/spde2mesh.R new file mode 100644 index 00000000..91f4136e --- /dev/null +++ b/R/spde2mesh.R @@ -0,0 +1,39 @@ +#' Convert an `Spde` object to a mesh object +#' @param object The `Spde` object +#' @name spde2mesh +#' @rdname spde2mesh +#' @exportMethod spde2mesh +#' @docType methods +#' @importFrom methods setGeneric +setGeneric( + name = "spde2mesh", + def = function(object) { + standardGeneric("spde2mesh") # nocov + } +) + +#' @rdname spde2mesh +#' @importFrom methods setMethod new +#' @include spde_class.R +setMethod( + f = "spde2mesh", + signature = signature(object = "Spde"), + definition = function(object) { + stopifnot( + "fmesher package required but not installed." = + requireNamespace("fmesher", quietly = TRUE), + "sf package required but not installed." = + requireNamespace("sf", quietly = TRUE) + ) + max_dist <- object@Range[1] + object@Coordinates |> + sf::st_as_sf(coords = colnames(object@Coordinates)) |> + sf::st_buffer(dist = max_dist) |> + sf::st_union() |> + sf::st_simplify(dTolerance = max_dist / 10) -> region + fmesher::fm_mesh_2d_inla( + boundary = region, max.edge = c(max_dist / 3, max_dist * 2), + cutoff = max_dist / 10 + ) + } +) diff --git a/R/spde_class.R b/R/spde_class.R new file mode 100644 index 00000000..68fe52a7 --- /dev/null +++ b/R/spde_class.R @@ -0,0 +1,35 @@ +#' The `spde` class +#' +#' It holds the coordinates, range and sigma parameters for the SPDE model. +#' @slot Coordinates a `data.frame` with the coordinates used for the mesh. +#' @slot Range a numeric vector of length 2. +#' Will be used as the `prior.range` argument of `[INLA::inla.spde2.pcmatern]`. +#' @slot Sigma a numeric vector of length 2. +#' Will be used as the `prior.sigma` argument of `[INLA::inla.spde2.pcmatern]`. +#' @name Spde-class +#' @rdname Spde-class +#' @exportClass Spde +#' @aliases Spde-class +#' @importFrom methods setClass +#' @docType class +setClass( + "Spde", + representation = representation( + Coordinates = "data.frame", Range = "numeric", Sigma = "numeric" + ) +) + +#' @importFrom assertthat assert_that noNA +#' @importFrom methods setValidity +setValidity( + "Spde", + function(object) { + assert_that( + noNA(object@Coordinates), noNA(object@Range), noNA(object@Sigma), + ncol(object@Coordinates) == 2, length(object@Range) == 2, + length(object@Sigma) == 2, all(object@Range > 0), all(object@Sigma > 0), + object@Range[2] < 1, object@Sigma[2] < 1 + ) + return(TRUE) + } +) diff --git a/_pkgdown.yml b/_pkgdown.yml index b98f97fe..844160cf 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,27 +1,32 @@ +url: https://inbo.github.io/n2kanalysis/ +template: + bootstrap: 5 + light-switch: false navbar: title: ~ type: default - left: - - text: NEWS - href: news/index.html - - text: Tutorials - href: articles/index.html - menu: - - text: Using n2kanalysis to analyse monitoring data - href: articles/workflow.html - - text: Functions - href: reference/index.html - - text: Contributing - href: CONTRIBUTING.html - right: - - icon: "fa fa-github" - href: https://github.com/inbo/n2kanalysis - - icon: "fa fa-twitter" - href: https://twitter.com/INBOVlaanderen - - icon: "fa fa-facebook" - href: https://www.facebook.com/pg/INBOVlaanderen + structure: + left: [intro, reference, news, tutorials, contributing] + right: [search, github] + components: + news: + text: Changelog + href: news/index.html + contributing: + text: Contributing + menu: + - text: Contributing + href: CONTRIBUTING.html + - text: Code of conduct + href: CODE_OF_CONDUCT.html + tutorials: + text: Tutorials + href: articles/index.html + menu: + - text: Using n2kanalysis to analyse monitoring data + href: articles/workflow.html authors: Research Institute for Nature and Forest (INBO): href: "https://www.vlaanderen.be/inbo/en-gb" - html: "" + html: "logo of the Research Institute for Nature and Forest (INBO)" diff --git a/inst/CITATION b/inst/CITATION index f0100d0f..8146bc88 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -2,12 +2,12 @@ citHeader("To cite `n2kanalysis` in publications please use:") # begin checklist entry bibentry( bibtype = "Manual", - title = "n2kanalysis: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring. Version 0.3.2", + title = "n2kanalysis: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring. Version 0.4.0", author = c( author = c(person(given = "Thierry", family = "Onkelinx"))), - year = 2024, - url = "https://github.com/inbo/n2kanalysis/", + year = 2025, + url = "https://inbo.github.io/n2kanalysis/", abstract = "All generic functions and classes for the analysis for the 'Natura 2000' monitoring. The classes contain all required data and definitions to fit the model without the need to access other sources. Potentially they might need access to one or more parent objects. An aggregation object might for example need the result of an imputation object. The actual definition of the analysis, using these generic function and classes, is defined in dedictated analysis R packages for every monitoring scheme. For example 'abvanalysis' and 'watervogelanalysis'.", - textVersion = "Onkelinx, Thierry (2024) n2kanalysis: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring. Version 0.3.2. https://github.com/inbo/n2kanalysis/", + textVersion = "Onkelinx, Thierry (2025) n2kanalysis: Generic Functions to Analyse Data from the 'Natura 2000' Monitoring. Version 0.4.0. https://inbo.github.io/n2kanalysis/", keywords = "analysis, reproducible research, natura 2000, monitoring", doi = "10.5281/zenodo.3576047", ) diff --git a/man/Spde-class.Rd b/man/Spde-class.Rd new file mode 100644 index 00000000..3a34d1df --- /dev/null +++ b/man/Spde-class.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spde_class.R +\docType{class} +\name{Spde-class} +\alias{Spde-class} +\title{The \code{spde} class} +\description{ +It holds the coordinates, range and sigma parameters for the SPDE model. +} +\section{Slots}{ + +\describe{ +\item{\code{Coordinates}}{a \code{data.frame} with the coordinates used for the mesh.} + +\item{\code{Range}}{a numeric vector of length 2. +Will be used as the \code{prior.range} argument of \verb{[INLA::inla.spde2.pcmatern]}.} + +\item{\code{Sigma}}{a numeric vector of length 2. +Will be used as the \code{prior.sigma} argument of \verb{[INLA::inla.spde2.pcmatern]}.} +}} + diff --git a/man/extract.Rd b/man/extract.Rd index 9fdce4d0..4bd2aadf 100644 --- a/man/extract.Rd +++ b/man/extract.Rd @@ -7,6 +7,8 @@ \alias{extract,character-methods} \alias{extract,ANY,n2kInla-method} \alias{extract,n2kInla-methods} +\alias{extract,ANY,n2kModelImputed-method} +\alias{extract,n2kModelImputed-methods} \title{Extract the relevant coefficients} \usage{ extract(extractor, object, base, project) @@ -14,6 +16,8 @@ extract(extractor, object, base, project) \S4method{extract}{ANY,character}(extractor, object, base, project) \S4method{extract}{ANY,n2kInla}(extractor, object, base = NULL, project = NULL) + +\S4method{extract}{ANY,n2kModelImputed}(extractor, object, base = NULL, project = NULL) } \arguments{ \item{extractor}{the extractor function} diff --git a/man/fit_model.Rd b/man/fit_model.Rd index beb447d6..a89abd58 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -3,7 +3,8 @@ % R/fit_model_n2k_aggregate.R, R/fit_model_n2k_composite.R, % R/fit_model_n2k_hurdle_imputed.R, R/fit_model_n2k_inla.R, % R/fit_model_n2k_inla_comparison.R, R/fit_model_n2k_manifest.R, -% R/fit_model_n2k_model_imputed.R, R/fit_model_s3_object.R +% R/fit_model_n2k_model_imputed.R, R/fit_model_n2k_spde.R, +% R/fit_model_s3_object.R \docType{methods} \name{fit_model} \alias{fit_model} @@ -15,6 +16,7 @@ \alias{fit_model,n2kInlaComparison-method} \alias{fit_model,n2kManifest-method} \alias{fit_model,n2kModelImputed-method} +\alias{fit_model,n2kSpde-method} \alias{fit_model,s3_object-method} \title{Fit an \code{n2kModel} object} \usage{ @@ -55,12 +57,21 @@ fit_model(x, ...) status = c("new", "waiting"), verbose = TRUE, ..., - local = tempfile("fit_model"), - first = FALSE + local = NULL ) \S4method{fit_model}{n2kModelImputed}(x, ...) +\S4method{fit_model}{n2kSpde}( + x, + status = "new", + ..., + timeout = NULL, + seed = get_seed(x), + num_threads = NULL, + parallel_configs = TRUE +) + \S4method{fit_model}{s3_object}(x, status = c("new", "waiting"), ...) } \arguments{ @@ -111,13 +122,6 @@ using \code{B:0} threads.} \item{local}{A local folder into which objects from an AWS S3 bucket are downloaded.} - -\item{first}{A logical. -\code{first = TRUE} implies to fit only the first object in the manifest with -matching status. -\code{first = FALSE} implies to fit all objects in the manifest with matching -status. -Defaults to \code{FALSE}.} } \description{ Fit an \code{n2kModel} object diff --git a/man/get_datafield_id.Rd b/man/get_datafield_id.Rd new file mode 100644 index 00000000..1bb2ad16 --- /dev/null +++ b/man/get_datafield_id.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_datafield_id.R +\name{get_datafield_id} +\alias{get_datafield_id} +\title{Get the data field id} +\usage{ +get_datafield_id(table, field, datasource, root, stage = FALSE) +} +\arguments{ +\item{table}{The table name} + +\item{field}{The field name} + +\item{datasource}{The data source name} + +\item{root}{The root of a project. Can be a file path or a \code{git-repository}. +Defaults to the current working directory (\code{"."}).} + +\item{stage}{Logical value indicating whether to stage the changes after +writing the data. Defaults to \code{FALSE}.} +} +\description{ +Get the data field id +} diff --git a/man/make_a.Rd b/man/make_a.Rd new file mode 100644 index 00000000..f761343b --- /dev/null +++ b/man/make_a.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/make_a.R +\docType{methods} +\name{make_a} +\alias{make_a} +\alias{make_a,Spde-method} +\title{Convert an \code{Spde} object to a A object} +\usage{ +make_a(object, data) + +\S4method{make_a}{Spde}(object, data) +} +\arguments{ +\item{object}{The \code{Spde} object} + +\item{data}{The \code{data.frame} to combine with the \code{Spde} object to create the +\code{A} matrix.} +} +\description{ +Convert an \code{Spde} object to a A object +} diff --git a/man/manifest_yaml_to_bash.Rd b/man/manifest_yaml_to_bash.Rd index 6ef9fb59..c7096eb7 100644 --- a/man/manifest_yaml_to_bash.Rd +++ b/man/manifest_yaml_to_bash.Rd @@ -36,7 +36,8 @@ manifest_yaml_to_bash( shutdown = FALSE, split = 1, status = c("new", "waiting"), - limit = FALSE + limit = FALSE, + timeout = 4 ) } \arguments{ diff --git a/man/moving_average.Rd b/man/moving_average.Rd new file mode 100644 index 00000000..520bfa32 --- /dev/null +++ b/man/moving_average.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/moving_average.R +\name{moving_average} +\alias{moving_average} +\title{Calculate moving trend coefficients} +\usage{ +moving_average(n_year, duration, first_year = 0) +} +\arguments{ +\item{n_year}{Number of available years in the data.} + +\item{duration}{Number of years in the moving window. +If the number of available years is less than this value, the trend will be +calculated for the available years.} + +\item{first_year}{First year of the data. +Only used to name the rows.} +} +\value{ +A matrix with moving average coefficients +One row for each window and and one column for each year in the data. +The format of the row names is \verb{average_\{window mid point\}_\{window length\}}. +\code{average_2001.5_4} is the average for the years 2000 to 2003. +} +\description{ +Calculate moving trend coefficients +} +\examples{ +moving_average(5, 3) +moving_average(5, 3, 2000) +} diff --git a/man/moving_difference.Rd b/man/moving_difference.Rd new file mode 100644 index 00000000..b9fa099f --- /dev/null +++ b/man/moving_difference.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/moving_difference.R +\name{moving_difference} +\alias{moving_difference} +\title{Calculate coefficients for a moving difference} +\usage{ +moving_difference(n_year, duration, first_year = 1) +} +\arguments{ +\item{n_year}{Number of available years in the data.} + +\item{duration}{Number of years in the moving window. +If the number of available years is less than this value, the trend will be +calculated for the available years.} + +\item{first_year}{First year of the data. +Only used to name the rows.} +} +\value{ +A matrix with moving average coefficients +One row for each window and and one column for each year in the data. +The format of the row names is +\verb{difference_\{window mid point start\}_\{window mid point end\}_\{window length\}}. +\code{difference_2001.5_2010.5_4} is the difference of the average for period for +the years 2009 to 2012 compared to the period from 2000 to 2003. +} +\description{ +Calculate coefficients for a moving difference +} +\examples{ +moving_difference(6, 2) +moving_difference(6, 2, 2000) +} diff --git a/man/moving_trend.Rd b/man/moving_trend.Rd new file mode 100644 index 00000000..e3d555c6 --- /dev/null +++ b/man/moving_trend.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/moving_trend.R +\name{moving_trend} +\alias{moving_trend} +\title{Calculate moving trend coefficients} +\usage{ +moving_trend(n_year, duration, first_year = 0) +} +\arguments{ +\item{n_year}{Number of available years in the data.} + +\item{duration}{Number of years in the moving window. +If the number of available years is less than this value, the trend will be +calculated for the available years.} + +\item{first_year}{First year of the data. +Only used to name the rows.} +} +\value{ +A matrix with the moving trend coefficients. +One row for each window and and one column for each year in the data. +The format of the row names is \verb{trend_\{window mid point\}_\{window length\}}. +\code{trend_2001.5_4} is the trend for the years 2000 to 2003. +} +\description{ +Calculate moving trend coefficients +} +\examples{ +moving_trend(5, 3) +moving_trend(5, 4, 2000) +} diff --git a/man/n2kSpde-class.Rd b/man/n2kSpde-class.Rd new file mode 100644 index 00000000..d9fa67d6 --- /dev/null +++ b/man/n2kSpde-class.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/n2k_spde_class.R +\docType{class} +\name{n2kSpde-class} +\alias{n2kSpde-class} +\title{The \code{n2kSpde} class} +\description{ +It hold analysis data based on an INLA model with SPDE. +} +\section{Slots}{ + +\describe{ +\item{\code{Spde}}{A \code{list} containing the information for the SPDE.} +}} + diff --git a/man/n2k_spde.Rd b/man/n2k_spde.Rd new file mode 100644 index 00000000..d9b7352b --- /dev/null +++ b/man/n2k_spde.Rd @@ -0,0 +1,124 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/n2k_spde.R +\docType{methods} +\name{n2k_spde} +\alias{n2k_spde} +\alias{n2k_spde,data.frame,ANY-method} +\alias{n2k_spde,n2kSpde-methods} +\alias{n2k_spde,n2kSpde,inla-method} +\title{Create an \code{n2kSpde} object} +\usage{ +n2k_spde(data, ..., model_fit) + +\S4method{n2k_spde}{data.frame,ANY}( + data, + status = "new", + result_datasource_id, + scheme_id, + family = "poisson", + formula, + species_group_id, + location_group_id, + model_type, + spde, + first_imported_year, + last_imported_year, + duration, + last_analysed_year, + analysis_date, + lin_comb = NULL, + minimum = "", + imputation_size, + parent = character(0), + seed, + replicate_name = list(), + control = list(), + parent_status = "converged", + parent_statusfingerprint, + extra, + ..., + model_fit +) + +\S4method{n2k_spde}{n2kSpde,inla}(data, status, raw_imputed = NULL, ..., model_fit) +} +\arguments{ +\item{data}{a \code{data.frame} with the data to analyse} + +\item{...}{other arguments. See below} + +\item{model_fit}{The fitted model} + +\item{status}{A single character indicating the status of the model. +Defaults to \code{"waiting"}.} + +\item{result_datasource_id}{A string identifying the data source.} + +\item{scheme_id}{A single integer holding the id of the scheme.} + +\item{family}{the family to use in the INLA model.} + +\item{formula}{A single character identifying the comparison.} + +\item{species_group_id}{A string identifying the species group.} + +\item{location_group_id}{A string identifying the location group.} + +\item{model_type}{The type of the models. +Must start with \code{"inla comparison:"}.} + +\item{spde}{the \code{Spde} object. See \verb{[spde]}.} + +\item{first_imported_year}{Oldest year considered in the data.} + +\item{last_imported_year}{Most recent year considered in the data.} + +\item{duration}{The width of the moving window. +Defaults to the \code{last_imported_year - first_imported_year + 1}.} + +\item{last_analysed_year}{Most recent year in the window. +Defaults to \code{last_imported_year}.} + +\item{analysis_date}{A \code{POSIXct} date indicating the date that the dataset +was imported.} + +\item{lin_comb}{A model matrix to calculate linear combinations.} + +\item{minimum}{The name of the variable which holds the minimum counts. +Only relevant in case of multiple imputation.} + +\item{imputation_size}{The required number of imputations defaults to 0.} + +\item{parent}{The file fingerprint of the optional parent analysis.} + +\item{seed}{A single integer used as a seed for all calculations. +A random seed will be inserted when missing.} + +\item{replicate_name}{A list with the names of replicates. +Defaults to an empty list. +Used in case of \code{f(X, ..., replicate = Z)}. +Should be a named list like e.g. \code{list(X = c("a", "b", "c"))}.} + +\item{control}{A named list passed to \code{\link[INLA:inla]{INLA::inla()}} when fitting +the model.} + +\item{parent_status}{A \code{data.frame} with columns +\code{parent_analysis} (the file fingerprint of the parent), +\code{parentstatus_fingerprint} (the status fingerprint of the parent), +and \code{parent_status} (the status of the parent).} + +\item{parent_statusfingerprint}{The status fingerprint of the parent +analysis.} + +\item{extra}{a \code{data.frame} with extra observations not used in the model. +They will be added in subsequent analyses.} + +\item{raw_imputed}{the optional \code{rawImputed} object} +} +\description{ +A new \code{n2kSpde} model is created when \code{data} is a \code{data.frame}. + +In case \code{data} is an \code{n2kSpde} object, then only the model and +status are updated. +All other slots are unaffected. +} diff --git a/man/order_manifest.Rd b/man/order_manifest.Rd new file mode 100644 index 00000000..d9178649 --- /dev/null +++ b/man/order_manifest.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/manifest_yaml_to_bash.R +\name{order_manifest} +\alias{order_manifest} +\title{Convert a manifest in vector of analysis fingerprints} +\usage{ +order_manifest(manifest) +} +\arguments{ +\item{manifest}{the \code{n2kManifest}} +} +\description{ +The order of the analysis is determined by the parent-child relationship. +It starts with the analyses without parents. +Then it adds the analyses with parents that have already been added. +This process is repeated until all analyses have been added. +} diff --git a/man/spde.Rd b/man/spde.Rd new file mode 100644 index 00000000..07aeecf7 --- /dev/null +++ b/man/spde.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spde.R +\name{spde} +\alias{spde} +\title{Create an \code{Spde} object} +\usage{ +spde(coordinates, range, sigma) +} +\arguments{ +\item{coordinates}{a \code{data.frame} of coordinates use to define the mesh.} + +\item{range}{a numeric vector of length 2. +Will be used as the \code{prior.range} argument of \verb{[INLA::inla.spde2.pcmatern]}.} + +\item{sigma}{a numeric vector of length 2. +Will be used as the \code{prior.sigma} argument of \verb{[INLA::inla.spde2.pcmatern]}.} +} +\description{ +Create an \code{Spde} object +} diff --git a/man/spde2matern.Rd b/man/spde2matern.Rd new file mode 100644 index 00000000..75450131 --- /dev/null +++ b/man/spde2matern.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spde2matern.R +\docType{methods} +\name{spde2matern} +\alias{spde2matern} +\title{Convert an \code{Spde} object to a \code{matern} object} +\usage{ +spde2matern(object) +} +\arguments{ +\item{object}{The \code{Spde} object} +} +\description{ +Convert an \code{Spde} object to a \code{matern} object +} diff --git a/man/spde2mesh.Rd b/man/spde2mesh.Rd new file mode 100644 index 00000000..a4c35b3b --- /dev/null +++ b/man/spde2mesh.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spde2matern.R, R/spde2mesh.R +\docType{methods} +\name{spde2matern,Spde-method} +\alias{spde2matern,Spde-method} +\alias{spde2mesh} +\alias{spde2mesh,Spde-method} +\title{Convert an \code{Spde} object to a mesh object} +\usage{ +\S4method{spde2matern}{Spde}(object) + +spde2mesh(object) + +\S4method{spde2mesh}{Spde}(object) +} +\arguments{ +\item{object}{The \code{Spde} object} +} +\description{ +Convert an \code{Spde} object to a mesh object +} diff --git a/n2kanalysis.Rproj b/n2kanalysis.Rproj index b067fd56..24461809 100644 --- a/n2kanalysis.Rproj +++ b/n2kanalysis.Rproj @@ -1,4 +1,5 @@ Version: 1.0 +ProjectId: 5ea8d182-fbfb-485c-8a11-cd3348bf81a0 RestoreWorkspace: No SaveWorkspace: No diff --git a/pkgdown/extra.css b/pkgdown/extra.css index 00938dd1..dfdfc0b2 100644 --- a/pkgdown/extra.css +++ b/pkgdown/extra.css @@ -18,13 +18,20 @@ a:hover { .navbar, .label-default, -.navbar-default .navbar-nav>.active>a, .navbar-default .navbar-nav>.active>a:hover, .navbar-default .navbar-nav>.active>a:focus { - background-color: #356196; +.navbar-default .navbar-nav>.active>a, .navbar-default .navbar-nav>.active>a:hover, .navbar-default .navbar-nav>.active>a:focus, +#toc>.nav a.nav-link.active { + background-color: #356196 !important; } -.navbar-default .navbar-link, -.navbar-default .navbar-nav>li>a { - color: #ffffff; +#toc>.nav a.nav-link { + background-color: #729BB7 !important; +} + +.navbar, +.navbar-brand, +.nav-link, +.nav-text.text-muted { + color: #ffffff !important; } .nav-pills li.active>a, .nav-pills li>a:hover { diff --git a/tests/testthat/test_aaa_get_datafield_id.R b/tests/testthat/test_aaa_get_datafield_id.R new file mode 100644 index 00000000..d6870d10 --- /dev/null +++ b/tests/testthat/test_aaa_get_datafield_id.R @@ -0,0 +1,24 @@ +test_that("get_datafield_id", { + root <- tempfile("get_datafield_id") + dir.create(root) + on.exit(unlink(root, recursive = TRUE)) + + expect_is( + get_datafield_id( + table = "test", field = "id", datasource = "database", root = root + ), + "integer" + ) + expect_is( + get_datafield_id( + table = "test", field = "id", datasource = "database", root = root + ), + "integer" + ) + expect_is( + get_datafield_id( + table = "test2", field = "id", datasource = "database", root = root + ), + "integer" + ) +}) diff --git a/tests/testthat/test_abb_n2k_spde.R b/tests/testthat/test_abb_n2k_spde.R new file mode 100644 index 00000000..a92af460 --- /dev/null +++ b/tests/testthat/test_abb_n2k_spde.R @@ -0,0 +1,40 @@ +this_result_datasource_id <- sha1(sample(letters)) +this_scheme_id <- sha1(sample(letters)) +this_species_group_id <- sha1(sample(letters)) +this_location_group_id <- sha1(sample(letters)) +this_analysis_date <- Sys.time() +this_model_type <- "inla poisson: A * (B + C) + C:D" +this_formula <- + "Count ~ A * (B + C) + C:D + + f(E, model = \"rw1\", replicate = as.integer(A)) + + f(G, model = \"iid\")" +this_first_imported_year <- 1990L +this_last_imported_year <- 2015L +this_last_analysed_year <- 2014L +this_duration <- 1L +dataset <- test_data() +this_lc <- dataset %>% + select("A", "B", "C", "D") %>% + filter(.data$C == max(.data$C), .data$D == max(.data$D)) %>% + distinct() %>% + model.matrix(object = ~A * (B + C) + C:D) +test_that("n2k_spde() creates the object", { + expect_s4_class( + spde_model <- spde( + dataset[, c("C", "D")], range = c(0.5, 0.05), sigma = c(0.5, 0.05) + ), + "Spde" + ) + expect_s4_class( + object <- n2k_spde( + result_datasource_id = this_result_datasource_id, + scheme_id = this_scheme_id, species_group_id = this_species_group_id, + location_group_id = this_location_group_id, model_type = this_model_type, + formula = this_formula, spde = spde_model, data = dataset, + first_imported_year = this_first_imported_year, + last_imported_year = this_last_imported_year, + analysis_date = this_analysis_date + ), + "n2kSpde" + ) +}) diff --git a/tests/testthat/test_caa_fit_model.R b/tests/testthat/test_caa_fit_model.R index 08a42518..a5839e1d 100644 --- a/tests/testthat/test_caa_fit_model.R +++ b/tests/testthat/test_caa_fit_model.R @@ -218,7 +218,8 @@ test_that("fit_model() works on n2kInlaComparison", { fit_model(filename2, verbose = FALSE) fit_model(filename3, verbose = FALSE) filename3 <- gsub("waiting", "converged", filename3) - expect_null(fit_model(filename3, verbose = FALSE)) + expect_invisible(output <- fit_model(filename3, verbose = FALSE)) + expect_is(output, "data.frame") }) test_that("fit_model() works on n2kInlaComposite", { @@ -294,7 +295,8 @@ test_that("fit_model() works on n2kInlaComposite", { fit_model(filename2, verbose = FALSE) fit_model(filename3, verbose = FALSE) filename3 <- gsub("waiting", "converged", filename3) - expect_null(fit_model(filename3, verbose = FALSE)) + expect_invisible(output <- fit_model(filename3, verbose = FALSE)) + expect_is(output, "data.frame") }) test_that("fit_model() works on n2kHurdleImputed", { @@ -337,18 +339,30 @@ test_that("fit_model() works on n2kHurdleImputed", { sha_count <- store_model(count, base = base, project = project) sha_presence <- store_model(presence, base = base, project = project) sha_hurdle <- store_model(hurdle, base = base, project = project) - expect_null(fit_model(basename(sha_hurdle), base = base, project = project)) + expect_invisible( + output <- fit_model(basename(sha_hurdle), base = base, project = project) + ) + expect_is(output, "data.frame") suppressWarnings( - expect_null(fit_model(basename(sha_count), base = base, project = project)) + expect_invisible( + output <- fit_model(basename(sha_count), base = base, project = project) + ) ) + expect_is(output, "data.frame") suppressWarnings( - expect_null( - fit_model(basename(sha_presence), base = base, project = project) + expect_invisible( + output <- fit_model( + basename(sha_presence), base = base, project = project + ) ) ) + expect_is(output, "data.frame") suppressWarnings( - expect_null(fit_model(basename(sha_hurdle), base = base, project = project)) + expect_invisible( + output <- fit_model(basename(sha_hurdle), base = base, project = project) + ) ) + expect_is(output, "data.frame") expect_s4_class( basename(sha_hurdle) |> read_model(base = base, project = project) |> diff --git a/tests/testthat/test_cba_fit_model_manifest.R b/tests/testthat/test_cba_fit_model_manifest.R index 526a8321..2f0f0ec8 100644 --- a/tests/testthat/test_cba_fit_model_manifest.R +++ b/tests/testthat/test_cba_fit_model_manifest.R @@ -38,20 +38,32 @@ test_that("it handles a manifest", { NA, get_file_fingerprint(object), get_file_fingerprint(object2) ), stringsAsFactors = FALSE - ) %>% + ) |> n2k_manifest() + hash <- store_manifest_yaml( + x = x, base = base, project = project, docker = "inbobmk/rn2k:dev-0.10", + dependencies = c("inbo/n2khelper@v0.5.0", "inbo/n2kanalysis@0.4.0") + ) + script <- manifest_yaml_to_bash( + base = base, project = project, hash = basename(hash) + ) expect_invisible(fit_model(x, base = base, project = project)) - x <- store_manifest(x, base, project) - expect_null(fit_model(x, base = base, project = project)) - expect_null(fit_model(x)) + y <- store_manifest(x, base, project) + expect_null(fit_model(y, base = base, project = project)) + expect_null(fit_model(y)) - file.path(base, project) %>% - list.files(recursive = TRUE, full.names = TRUE) %>% + file.path(base, project) |> + list.files(recursive = TRUE, full.names = TRUE) |> + c( + R_user_dir("n2kanalysis", which = "cache") |> + file.path(x@Fingerprint) + ) |> file.remove() + # works with an S3 bucket skip_if(Sys.getenv("AWS_SECRET_ACCESS_KEY") == "", message = "No AWS access") - aws_base <- get_bucket(Sys.getenv("N2KBUCKET")) + aws_base <- get_bucket(Sys.getenv("N2KBUCKET"), max = 1) store_model(object, base = aws_base, project = project) store_model(object2, base = aws_base, project = project) store_model(object3, base = aws_base, project = project) @@ -64,18 +76,38 @@ test_that("it handles a manifest", { NA, get_file_fingerprint(object), get_file_fingerprint(object2) ), stringsAsFactors = FALSE - ) %>% + ) |> n2k_manifest() + hash <- store_manifest_yaml( + x = x, base = aws_base, project = project, docker = "inbobmk/rn2k:dev-0.10", + dependencies = c("inbo/n2khelper@v0.5.0", "inbo/n2kanalysis@0.4.0") + ) + script <- manifest_yaml_to_bash( + base = aws_base, project = project, hash = basename(hash) + ) + expect_s3_class( + results <- get_result(x, base = aws_base, project = project), + "data.frame" + ) + expect_true(all(results$status == "new")) expect_invisible( fit_model(x, base = aws_base, project = project, verbose = TRUE) ) + expect_s3_class( + results <- get_result(x, base = aws_base, project = project), + "data.frame" + ) + expect_true(all(results$status == "converged")) - x <- store_manifest(x, base = aws_base, project = project) - expect_invisible(fit_model(x$Contents)) + y <- store_manifest(x, base = aws_base, project = project) + expect_invisible(fit_model(y$Contents)) - expect_null(fit_model(x$Contents$Key, base = aws_base, project = project)) + expect_null(fit_model(y$Contents$Key, base = aws_base, project = project)) available <- get_bucket(aws_base, prefix = project) %>% sapply("[[", "Key") expect_true(all(sapply(available, delete_object, bucket = aws_base))) + R_user_dir("n2kanalysis", which = "cache") |> + file.path(x@Fingerprint) |> + file.remove() })