[WIP] Remove mean F value calculation for selectivity curve

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(generate_selectivity_curve)
 export(get_functional_groups)
 export(load_csv_environmental_data)
 export(load_csv_ewe)

R/fit_selectivity_helpers.R

@@ -0,0 +1,235 @@
+#' Fit logistic selectivity parameters
+#'
+#' @description
+#' Fit logistic curve parameters to a vector of normalized values. The logistic
+#' function is defined as: S(x) = 1 / (1 + exp(-slope * (x - inflection_point)))
+#' where x is the index position (1, 2, 3, ...).
+#'
+#' @param f_values A numeric vector of normalized values (0-1) for each
+#'   group, in the order they should appear on the x-axis. These values
+#'   represent relative selectivity, maturity, or any other quantity that
+#'   follows a logistic pattern.
+#'
+#' @return A tibble with one row per parameter, containing the standard FIMS
+#'   parameter format columns.
+#'
+#' @details
+#' The function uses nonlinear least squares (nls) to fit a logistic curve to
+#' the normalized values. The x-axis index is simply 1, 2, 3, ... based on
+#' the position in the values vector.
+#'
+#' If the optimization fails, the function falls back to a simple heuristic
+#' method to estimate the inflection point and slope.
+#'
+#' @keywords internal
+fit_logistic_selectivity <- function(f_values) {
+  # Create index (1, 2, 3, ...)
+  index_values <- seq_along(f_values)
+
+  # Define logistic function
+  logistic <- function(index, inflection_point, slope) {
+    1 / (1 + exp(-slope * (index - inflection_point)))
+  }
+
+  # Initial parameter guesses
+  # Inflection point: index where value is around 0.5
+  midpoint_index <- which.min(abs(f_values - 0.5))
+  initial_inflection <- if (length(midpoint_index) > 0) index_values[midpoint_index] else mean(index_values)
+
+  # Slope: positive value, start with 1
+  initial_slope <- 1
+
+  # Fit using nonlinear least squares
+  fit_result <- tryCatch(
+    {
+      stats::nls(
+        f_values ~ logistic(index_values, inflection_point, slope),
+        start = list(inflection_point = initial_inflection, slope = initial_slope),
+        control = stats::nls.control(maxiter = 100, warnOnly = TRUE)
+      )
+    },
+    error = function(e) {
+      # Fallback: use simple heuristic
+      NULL
+    }
+  )
+
+  if (!is.null(fit_result)) {
+    params <- stats::coef(fit_result)
+    inflection_point <- params["inflection_point"]
+    slope <- params["slope"]
+  } else {
+    # Fallback method: find inflection as midpoint of rise
+    # and estimate slope from the steepest part of the curve
+    inflection_point <- initial_inflection
+
+    # Estimate slope from differences
+    if (length(index_values) > 1) {
+      diffs <- diff(f_values) / diff(index_values)
+      slope <- max(abs(diffs), 0.1) * 4 # Scale up since max derivative of logistic is slope/4
+    } else {
+      slope <- 1
+    }
+  }
+
+  # Create output tibble in FIMS format (replace time with year/month)
+  tibble::tibble(
+    model_family = "catch_at_age",
+    module_name = "Selectivity",
+    module_type = "Logistic",
+    label = c("inflection_point", "slope"),
+    distribution_link = NA_character_,
+    age = NA_real_,
+    length = NA_real_,
+    year = NA_integer_,
+    month = NA_integer_,
+    value = c(inflection_point, slope),
+    estimation_type = "fixed_effects",
+    distribution_type = NA_character_,
+    distribution = NA_character_
+  )
+}
+
+
+#' Fit double logistic selectivity parameters
+#'
+#' @description
+#' Fit double logistic (dome-shaped) curve parameters to a vector of normalized
+#' values. The double logistic function combines ascending and descending logistic
+#' curves: S(x) = 1 / (1 + exp(-slope_asc * (x - inflection_asc))) *
+#'              (1 - 1 / (1 + exp(-slope_desc * (x - inflection_desc))))
+#'
+#' @param f_values A numeric vector of normalized values (0-1) for each
+#'   group, in the order they should appear on the x-axis. These values
+#'   represent relative selectivity, maturity, or any other quantity that
+#'   follows a dome-shaped pattern.
+#'
+#' @return A tibble with one row per parameter, containing the standard FIMS
+#'   parameter format columns.
+#'
+#' @details
+#' The function attempts to fit a dome-shaped curve. If the data
+#' appears monotonic (no dome shape), it will still fit the parameters but may
+#' produce parameters that effectively reduce to a simple logistic.
+#'
+#' The fitting uses nonlinear least squares with constraints to ensure:
+#' \itemize{
+#'   \item Ascending inflection point comes before descending
+#'   \item Slopes are positive
+#'   \item The curve shape is biologically reasonable
+#' }
+#'
+#' @keywords internal
+fit_double_logistic_selectivity <- function(f_values) {
+  # Create index (1, 2, 3, ...)
+  index_values <- seq_along(f_values)
+
+  # Define double logistic function
+  double_logistic <- function(index, inflection_point_ascending, slope_ascending, inflection_point_descending, slope_descending) {
+    ascending_component <- 1 / (1 + exp(-slope_ascending * (index - inflection_point_ascending)))
+    descending_component <- 1 - 1 / (1 + exp(-slope_descending * (index - inflection_point_descending)))
+    ascending_component * descending_component
+  }
+
+  # Initial parameter guesses
+  # Find the peak of the values
+  peak_idx <- which.max(f_values)
+  peak_index <- index_values[peak_idx]
+
+  # Ascending parameters: before peak
+  if (peak_idx > 1) {
+    ascending_f <- f_values[1:peak_idx]
+    ascending_indices <- index_values[1:peak_idx]
+    midpoint_ascending_index <- which.min(abs(ascending_f - 0.5))
+    initial_inflection_point_ascending <- if (length(midpoint_ascending_index) > 0) ascending_indices[midpoint_ascending_index] else peak_index - 1
+  } else {
+    initial_inflection_point_ascending <- peak_index - 1
+  }
+
+  # Descending parameters: after peak
+  if (peak_idx < length(index_values)) {
+    descending_f <- f_values[peak_idx:length(index_values)]
+    descending_indices <- index_values[peak_idx:length(index_values)]
+    # For descending, look for where it drops below peak/2
+    midpoint_descending_index <- which.min(abs(descending_f - max(f_values) / 2))
+    initial_inflection_point_descending <- if (length(midpoint_descending_index) > 0) descending_indices[midpoint_descending_index] else peak_index + 1
+  } else {
+    initial_inflection_point_descending <- peak_index + 1
+  }
+
+  # Ensure descending inflection is after ascending
+  if (initial_inflection_point_descending <= initial_inflection_point_ascending) {
+    initial_inflection_point_descending <- initial_inflection_point_ascending + 1
+  }
+
+  initial_slope_ascending <- 1
+  initial_slope_descending <- 1
+
+  # Fit using nonlinear least squares
+  fit_result <- tryCatch(
+    {
+      stats::nls(
+        f_values ~ double_logistic(index_values, inflection_point_ascending, slope_ascending, inflection_point_descending, slope_descending),
+        start = list(
+          inflection_point_ascending = initial_inflection_point_ascending,
+          slope_ascending = initial_slope_ascending,
+          inflection_point_descending = initial_inflection_point_descending,
+          slope_descending = initial_slope_descending
+        ),
+        control = stats::nls.control(maxiter = 200, warnOnly = TRUE),
+        algorithm = "port",
+        lower = c(
+          inflection_point_ascending = min(index_values),
+          slope_ascending = 0.01,
+          inflection_point_descending = initial_inflection_point_ascending + 0.5,
+          slope_descending = 0.01
+        ),
+        upper = c(
+          inflection_point_ascending = max(index_values),
+          slope_ascending = 10,
+          inflection_point_descending = max(index_values) + 10,
+          slope_descending = 10
+        )
+      )
+    },
+    error = function(e) {
+      NULL
+    }
+  )
+
+  if (!is.null(fit_result)) {
+    params <- stats::coef(fit_result)
+    inflection_point_ascending <- params["inflection_point_ascending"]
+    slope_ascending <- params["slope_ascending"]
+    inflection_point_descending <- params["inflection_point_descending"]
+    slope_descending <- params["slope_descending"]
+  } else {
+    # Fallback to initial guesses
+    inflection_point_ascending <- initial_inflection_point_ascending
+    slope_ascending <- initial_slope_ascending
+    inflection_point_descending <- initial_inflection_point_descending
+    slope_descending <- initial_slope_descending
+  }
+
+  # Create output tibble in FIMS format (replace time with year/month)
+  tibble::tibble(
+    model_family = "catch_at_age",
+    module_name = "Selectivity",
+    module_type = "DoubleLogistic",
+    label = c(
+      "inflection_point_asc",
+      "slope_asc",
+      "inflection_point_desc",
+      "slope_desc"
+    ),
+    distribution_link = NA_character_,
+    age = NA_real_,
+    length = NA_real_,
+    year = NA_integer_,
+    month = NA_integer_,
+    value = c(inflection_point_ascending, slope_ascending, inflection_point_descending, slope_descending),
+    estimation_type = "fixed_effects",
+    distribution_type = NA_character_,
+    distribution = NA_character_
+  )
+}