## ----message = FALSE----------------------------------------------------------

library(broom)
library(dplyr)
library(photosynthesis)

# Read in your data
dat = system.file("extdata", "A_Ci_Q_data_1.csv", package = "photosynthesis") |>
  read.csv() |>
  # Set grouping variable
  mutate(group = round(CO2_s, digits = 0)) |>
  # For this example, round sequentially due to CO2_s set points
  mutate(group = as.factor(round(group, digits = -1)))

# Fit one light-response curve
fit = fit_photosynthesis(
  .data = filter(dat, group == 600),
  .photo_fun = "aq_response",
  .vars = list(.A = A, .Q = Qabs),
)

# The 'fit' object inherits class 'nls' and many methods can be used

## Model summary:
summary(fit)

## Estimated parameters:
coef(fit)

## 95% confidence intervals:
confint(fit)

## Tidy summary table using 'broom::tidy()'
tidy(fit, conf.int = TRUE, conf.level = 0.95)

## Calculate light compensation point
coef(fit) |>
  t() |>
  as.data.frame() |>
  mutate(LCP = ((Rd) * (Rd * theta_J - k_sat) / (phi_J * (Rd - k_sat)))) |>

## Calculate residual sum-of-squares
sum(resid(fit) ^ 2)


## -----------------------------------------------------------------------------
library(ggplot2)

b = coef(fit)

df_predict = data.frame(Qabs = seq(0, 0.84 * 1500, length.out = 100)) |>
  mutate(
    A = marshall_biscoe_1980(
      Q_abs = Qabs,
      k_sat = b["k_sat"],
      b["phi_J"],
      b["theta_J"]
    ) - b["Rd"]
  )

ggplot(mapping = aes(Qabs, A)) +
  geom_line(data = df_predict) +
  geom_point(data = filter(dat, group == 600)) +
  labs(
    x = expression("Irradiance (" * mu * mol ~ m^{-2} ~ s^{-1} * ")"),
    y = expression(A[net] ~ "(" * mu * mol ~ m^{-2} ~ s^{-1} * ")")
  ) +
  theme_bw()


## -----------------------------------------------------------------------------
library(purrr)

fits = dat |>
  split(~ group) |>
  map(fit_photosynthesis, .photo_fun = "aq_response", .vars = list(.A = A, .Q = Qabs))

## Estimated parameters:
fits |>
  map(coef) |>
  map(t) |>
  map(as.data.frame) |>
  imap_dfr(~ mutate(.x, CO2_s = .y))


## ----eval = FALSE-------------------------------------------------------------
# 
# fit = fit_photosynthesis(
#   .data = filter(dat, group == 600),
#   .photo_fun = "aq_response",
#   .vars = list(.A = A, .Q = Qabs),
#   .method = "brms",
#   brm_options = list(chains = 1)
# )
# 
# # The 'fit' object inherits class 'brmsfit' and many methods can be used
# summary(fit)
# 

## ----message = FALSE----------------------------------------------------------
library(broom)
library(dplyr)
library(tibble)
library(photosynthesis)

# Simulate data
set.seed(202411123)

## Parameters 
k_sat = 30
phi_J = 0.05
theta_J = 0.8
Rd = 3
b_inh = 0.003

## Vector of light values
Q = c(0, 25, 50, 75, 100, 125, 150, 200, 300, 400, 500, 600, 700, 800, 1000, 1200, 1400, 1600, 1800, 2000)

dat = tibble(
  Qabs = 0.84 * Q,
  A = photoinhibition(Qabs, k_sat, phi_J, theta_J, b_inh) - Rd + rnorm(length(Qabs), 0, 0.1)
)

# Fit one light-response curve
fit = fit_photosynthesis(
  .data = dat,
  .photo_fun = "aq_response",
  .model = "photoinhibition",
  .vars = list(.A = A, .Q = Qabs),
)

# The 'fit' object inherits class 'nls' and many methods can be used

## Model summary:
summary(fit)

## Estimated parameters:
coef(fit)

## 95% confidence intervals:
confint(fit)

## Tidy summary table using 'broom::tidy()'
tidy(fit, conf.int = TRUE, conf.level = 0.95)


## -----------------------------------------------------------------------------

library(ggplot2)

b = coef(fit)

df_predict = data.frame(Qabs = seq(0, 0.84 * 2000, length.out = 100)) |>
  mutate(
    A = photoinhibition(
      Q_abs = Qabs,
      k_sat = b["k_sat"],
      b["phi_J"],
      b["theta_J"],
      b["b_inh"]
    ) - b["Rd"]
  )

ggplot(mapping = aes(Qabs, A)) +
  geom_line(data = df_predict) +
  geom_point(data = dat) +
  labs(
    x = expression("Irradiance (" * mu * mol ~ m^{-2} ~ s^{-1} * ")"),
    y = expression(A[net] ~ "(" * mu * mol ~ m^{-2} ~ s^{-1} * ")")
  ) +
  theme_bw()


## ----message = FALSE----------------------------------------------------------

library(dplyr)
library(photosynthesis)
# Read in your data
dat = system.file("extdata", "A_Ci_Q_data_1.csv", package = "photosynthesis") |>
  read.csv() |>
  # Set grouping variable
  mutate(group = round(CO2_s, digits = 0)) |>
  # For this example, round sequentially due to CO2_s setpoints
  mutate(group = as.factor(round(group, digits = -1))) |>
  rename(A_net = A, PPFD = Qin)

# To fit one AQ curve
fit = fit_aq_response(filter(dat, group == 600))

# Print model summary
summary(fit[[1]])

# Print fitted parameters
fit[[2]]

# Print graph
fit[[3]]

# Fit many curves
fits = fit_many(
  data = dat,
  funct = fit_aq_response,
  group = "group",
  progress = FALSE
)

# Look at model summary for a given fit
# First set of double parentheses selects an individual group value
# Second set selects an element of the sublist
summary(fits[[3]][[1]])

# Print the parameters
fits[[2]][[2]]

# Print the graph
fits[[3]][[3]]

#Compile graphs into a list for plotting
fits_graphs = compile_data(fits, list_element = 3)

# Print graphs to jpeg
# print_graphs(data = fits_graphs, path = tempdir(), output_type = "jpeg")

#Compile parameters into data.frame for analysis
fits_pars = compile_data(fits, output_type = "dataframe", list_element = 2)