# This is a wrapper that comes with the PopJSON library (written in NodeJS)
# It exists in $POPJSON_WPATH or node_modules/popjson/wrappers/
# Here, it is symbolically linked to the working directory.
#
# sys.call(sprintf("cp /opt/conda/lib/node_modules/popjson/wrappers/population.* ./"))
#
source("population.R")

# This is a convenience script to parse PopJSON models into C (which implements the Population model)
# and compile as a "dylib" for access from R
# 
prepare <- function(filename) {
    out <- system(sprintf("popjson %s", filename), intern = TRUE)
    print(out)
    #
    if (exists("model") && inherits(model, "PopulationModel")) model$destroy()
    #
    model <- PopulationModel$new(sprintf("%s.dylib", filename))
    #
    return(model)
}

model <- prepare("models/example")

ret <- model$sim(ftime = 20, 
                 y0 = list(larva = 100.0))

plot(ret$ret[1, , model$popids[["larva"]]], t = "l")

df <- read.csv("data/data_semifield.csv")
df_time <- as.POSIXct(df$Date, format="%d/%m/%Y")

tm <- read.csv("data/datalogger.csv")
tm_time <- as.POSIXct(tm$Time, format="%Y-%m-%d %H:%M:%S")

plot(tm_time,tm$Tmean,t='l',ylim=c(0,40))
lines(df_time,df$L,t='b',col='blue')
lines(df_time,df$P,t='b',col='red')

dsig <-function(T,L,l,R,r,M,X) {
    return (((M) + ((X) / ((1.0 + exp(exp((l)) * ((L) - (T)))) * (1.0 + exp(exp((r)) * ((T) - ((L) + (R)))))))))
}

temp <- seq(0, 50, length.out = 1000)

plot(temp, dsig(temp,10,0,25,0,24,2000) / 24, t='l', col='blue')
lines(temp, dsig(temp,10,0,25,0,1000,-500) / 24, t='l', col='green')

model <- prepare("models/model_larva")

temp <- seq(0,50,length.out=1000)

param <- model$param

ret <- model$sim(
  ftime = length(temp)+1,
  envir = list(
    temp = temp
  ),
  pr = param,
  rep = -1
)

plot(temp, ret$iret[1, , model$intids[["n2m"]]] / 24.0, t='l', col='blue')
lines(temp, ret$iret[1, , model$intids[["d2m"]]] / 24.0, t='l', col='green')

temp <- tm$Tmean

param <- model$param

ret <- model$sim(
  ftime = length(temp),
  envir = list(
    temp = temp
  ),
  pr = param,
  y0 = list(
    larva = 20.0
  )
)

plot(ret$ret[1, , model$popids[["larva"]]], t = "l")
lines(cumsum(ret$iret[1, , model$intids[["larva_to_pupa"]]]), t = "l")

plot(tm_time,tm$Tmean,t='l',ylim=c(0,40))
lines(df_time,df$L,t='b',col='blue')
lines(df_time,df$P,t='b',col='red')
lines(tm_time,ret$ret[1, , model$popids[["larva"]]], t = "l")
lines(tm_time[-length(tm_time)],cumsum(ret$iret[1, , model$intids[["larva_to_pupa"]]]), t = "l")

obs <- df
obs['hour'] <- (1:nrow(obs))*12

simulate_system <- function(param) {
  temp <- tm$Tmean
  hours <- length(temp)

  ret <- model$sim(
    ftime = hours,
    envir = list(
      temp = temp
    ),
    pr = param,
    y0 = list(
      larva = 20.0
    )
  )

  return(data.frame(hour = (1:hours)-1, 
                    larva = ret$ret[1, , model$popids[["larva"]]],
                    pupa = c(cumsum(ret$iret[1, , model$intids[["larva_to_pupa"]]]), NA)))
}

objective_function <- function(param) {
  sim <- simulate_system(param)
  # Extract simulation only at obs times
  sim_at_obs <- sim[tm_time %in% df_time, ]
  # sum of squared errors (larva + pupa)
  SSE <- sum((obs$L - sim_at_obs$larva)^2 + (obs$P - sim_at_obs$pupa)^2)
  return(SSE)
}

objective_function(param)

lower_bounds <- model$parmin
upper_bounds <- model$parmax

fit <- optim(
  par = model$param,
  fn = objective_function,
  method = "L-BFGS-B",
  lower = lower_bounds,
  upper = upper_bounds
)

print(fit)

best_param <- fit$par
best_sim <- simulate_system(best_param)

plot(df_time,df$L,t='b',col='blue')
lines(df_time,df$P,t='b',col='red')

lines(tm_time, best_sim$larva, t='l')
lines(tm_time, best_sim$pupa, t='l')

best_sim_at_obs <- best_sim[tm_time %in% df_time, ]

lines(df_time, best_sim_at_obs$pupa, t='p')

# This works when the standard deviation is 1% of the mean
best_param <- c(-10.0, 6.0, 50.0, 6.0, 655.094187, 2631.398751, 26.608980, -6.0, -3.672721, -1.440932, 540.172220, 1045.062222)

best_sim <- simulate_system(best_param)
best_sim_at_obs <- best_sim[tm_time %in% df_time, ]

plot(df_time,df$L,t='b',col='blue')
lines(df_time,df$P,t='b',col='red')
lines(tm_time, best_sim$larva, t='l')
lines(tm_time, best_sim$pupa, t='l')
lines(df_time, best_sim_at_obs$pupa, t='p')

temp <- seq(0,50,length.out=1000)

ret <- model$sim(
  ftime = length(temp)+1,
  envir = list(
    temp = temp
  ),
  pr = best_param,
  rep = -1
)

plot(temp, ret$iret[1, , model$intids[["n2m"]]] / 24.0, t='l', col='blue')
lines(temp, ret$iret[1, , model$intids[["d2m"]]] / 24.0, t='l', col='green')

model <- prepare("models/model_aegypti")

print(model$popnames)
print(model$parnames)
print(model$param)

dat <- read.csv("https://veclim.com/files/workshops/September2025/clim/MissionI/clim_MissionI_Larnaca_V2025-09-11.csv")

source("aux_micro.R")

dat[['time']] <- as.Date(dat[['time']])

temp_raw <- dat[['temp']]
temp_mild <- micro_optimum(temp_raw,2.5)

plot(dat[['time']],temp_raw,t='l',col='black')
lines(dat[['time']],temp_mild,t='l',col='red')
legend("topright",
       legend = c("Raw temperature", "Slightly milder"),
       col = c("black", "red"),
       lty = 1, bty = "n")

pr <- c(10,-0.88795,25,-0.146617,400,1,10,-6.00179,-1.82367,16.5314,12877.2,1,240.108,0.543329,14.9951,0.9,12,0.6,22,0.929549,1000,0.5,0,2051.56,0.1,300,-4.8,1,0.1,17,0.6,17,0.929549,1000,0.5,36,2051.56,0.15,0.9,0.5,0.5,0.000923007,0.00205535,20.0499,100,100,26,-0.5,45,1700,1,0.1,16.5,1,15,-1.2,870,0.001,0.1)

sim <- model$sim(
    ftime = length(dat[['time']]),
    envir = list(
                 temp = temp_mild,
                 prec = dat[['prec']],
                 evap = dat[['evap']],
                 rehum = dat[['rhum']],
                 pdens = dat[['popd']],
                 experiment = c(0),
                 configure = c(1)
                ),
    pr = pr,
    y0 = list(pop_Eh = 1000.0)
    )

plot(dat[['time']], sim$ret[1, , model$popids[['pop_Af']]], t='l', xlab="Date", ylab="Number of females per grid cell")