Open-Systems-Pharmacology / OSPSuite.ParameterIdentification

R package for parameter identification for OSPS models
http://www.open-systems-pharmacology.org/OSPSuite.ParameterIdentification/
GNU General Public License v2.0
3 stars 2 forks source link

`plotOFVProfiles()` fails in the vignette example #91

Open PavelBal opened 4 months ago

PavelBal commented 4 months ago

In the current user guide, when trying to call plotOFVProfiles() after running the Aciclovir PI:

> profile <- piTask$calculateOFVProfiles()
Error in `tibble::as_tibble()` at OSPSuite.ParameterIdentification/R/ParameterIdentification.R:875:7:
! Column name `Neighborhoods|Kidney_pls_Kidney_ur|Aciclovir|Renal Clearances-TS|TSspec` must not be duplicated.
Use `.name_repair` to specify repair.
Caused by error in `repaired_names()`:
! Names must be unique.
✖ These names are duplicated:
  * "Neighborhoods|Kidney_pls_Kidney_ur|Aciclovir|Renal Clearances-TS|TSspec" at locations 2 and 3.
Backtrace:
 1. piTask$calculateOFVProfiles()
 3. tibble:::as_tibble.list(gridList)
 4. tibble:::lst_to_tibble(x, .rows, .name_repair, col_lengths(x))
 5. tibble:::set_repaired_names(...)
 6. tibble:::repaired_names(...)

The functionality of plotOFVProfiles() will be removed until fixed.

PavelBal commented 4 months ago

Removed documentation:

Around the global minimum, the objective function is expected to be convex and parabolic. We can check if we the solution is around a minimum (global or local) by plotting the objective function value (OFV) profile around the tentative optimal point. For example, after running the parameter identification task, run to plot the OFV profile:

{r include=TRUE, eval=TRUE} profile <- piTask$calculateOFVProfiles() plotOFVProfiles(profile)[[1]]

{r include=TRUE, echo=FALSE, eval=TRUE, fig.width=6, fig.height=4, fig.align='center'}
plotOFVProfiles(cachedAciclovirProfile)[[1]]

We observe that the tentative lipophilicity value of r signif(taskResults$par[[1]], 3) is close to at least a local (and potentially, a global) minimum. By passing explicit parameters to the calculateOFVProfiles function, we can check the whole parameter space, at least for 1-dimensional parameter estimation problems:

{r include=TRUE, eval=FALSE}
profile <- task$calculateOFVProfiles(lower = -10, upper = 10, totalEvaluations = 200)
plotOFVProfiles(profile)[[1]]
{r include=TRUE, echo=FALSE, eval=TRUE, fig.width=6, fig.height=4, fig.align='center'}
plotOFVProfiles(cachedAciclovirFullProfile)[[1]]
PavelBal commented 4 months ago

Code removed from ParameterIdentification:

    #' @description
    #' Calculates the values of the objective function on all orthogonal lines
    #' passing through a given point in the parameter space.
    #' @param par A vector of parameter values, with the same length as the number of parameters.
    #' If not supplied, the current parameter values are used.
    #' @param lower A vector of lower bounds for parameters, with the same length as the number of parameters.
    #' By default, uses 0.9 of the current parameter value.
    #' @param upper A vector of upper bounds for parameters, with the same length as the number of parameters.
    #' By default, uses 1.1 of the current parameter value.
    #' @param totalEvaluations An integer number. The combined profiles will not contain more than `totalEvaluations`
    #' points. If not supplied, 21 points per parameter are plotted to cover a uniform grid from 0.9 to 1.1.
    #' @return A list of tibbles, one tibble per parameter, with one column for parameter values
    #' and one column for the matching objective function values.
    calculateOFVProfiles = function(par = NULL, lower = NULL, upper = NULL, totalEvaluations = NULL) {
      # If the batches have not been initialized yet (i.e., no run has been
      # performed), this must be done prior to plotting
      private$.batchInitialization()

      nrOfParameters <- length(private$.piParameters)

      # if par is not supplied, we use the current parameter values
      if (missing(par)) {
        par <- unlist(lapply(private$.piParameters, function(x) {
          x$currValue
        }), use.names = FALSE)
      }

      # if lower and upper are not supplied, we calculate them as 0.9 and 1.1
      # of the current parameter values
      if (missing(lower)) {
        lower <- 0.9 * par
      }
      if (missing(upper)) {
        upper <- 1.1 * par
      }

      # calculate the grid for each parameter separately
      if (missing(totalEvaluations)) {
        gridSize <- 21
        # creates a grid with values at 0.9, 0.91, 0.92 .. 1.0 .. 1.09, 1.1
        # of the current parameter values
      } else {
        gridSize <- floor(totalEvaluations / nrOfParameters)
      }
      gridList <- vector(mode = "list", length = nrOfParameters)
      for (idx in seq_along(private$.piParameters)) {
        gridList[[idx]] <- rep(par[[idx]], gridSize)
        names(gridList)[[idx]] <- private$.piParameters[[idx]]$parameters[[1]]$path
      }
      defaultGrid <- tibble::as_tibble(gridList)

      profileList <- vector(mode = "list", length = nrOfParameters)
      for (idx in seq_along(private$.piParameters)) {
        # the names of the parameters are extracted from the first available path
        parameterName <- private$.piParameters[[idx]]$parameters[[1]]$path
        grid <- seq(from = lower[[idx]], to = upper[[idx]], length.out = gridSize)
        currentGrid <- defaultGrid
        currentGrid[[parameterName]] <- grid
        # creates a tibble with the column name from the `parameterName` variable
        profileList[[idx]] <- tibble::tibble(!!parameterName := grid)
        profileList[[idx]][["ofv"]] <- purrr::pmap_dbl(currentGrid, function(...) {
          private$.targetFunction(c(...))$model
        })
        names(profileList)[[idx]] <- parameterName
      }

      if (!is.null(private$.savedSimulationState)) {
        .restoreSimulationState(private$.simulations, private$.savedSimulationState)
      }

      return(profileList)
    },
PavelBal commented 4 months ago

Please add more explanation for the calculateOFVProfiles method. Examples with lower/upper. What is the output of the method, it's structure?

rengelke commented 4 months ago

Need to review and recover the global-algorithm-example once issue #91 and #92 are resolved