
begin
    using Markdown
    using Plots
end

md"""
Function that optimizes sinusoid parameters to fit $\delta^{18}O$ data.

The second core function of the ShellChron growth model. Loops
through all data windows and uses the `growth_model` function
to create d18O series that match the input data. This step is
iterated and optimized (minimizing the Sum of Squared Residuals)
through the SCEUA algorithm (by Duan et al., 1992) which finds
the optimal input parameters to the growth rate and Sea Surface
Temperature (SST) sinusoids to simulate d18O data.

[dat] Matrix containing the input data

[dynwindow] Information on the position and length of modeling
windows

[transfer_function] Transfer function used to convert d18Oc to temperature
data.

[$\delta^{18}O_w$] Either a single value (constant $\delta^{18}O_w$) or a vector of length
equal to the period in SST data (365 days by default) containing information
about seasonality in d18Ow. Defaults to constant d18Ow of 0 permille VSMOW
(the modern mean ocean value)

[T_per] Period of SST sinusoid (in days; default = 365)

[G_per] Period of growth rate sinusoid (in days; default = 365)

[t_int] Time interval (in days; default = 1)

[t_maxtemp] Timing of the warmest day of the year (in julian day; 
default = 182.5, or May 26th halfway through the year)

[SCEUApar] Parameters for SCEUA optimization (iniflg, ngs, maxn, kstop
pcento, peps). For details, refer to Duan et al. (1992) in references

[sinfit] Apply sinusoidal fitting to guess initial parameters for SCEUA
optimization? \code{true/FALSE}

[MC] Number of Monte Carlo simulations to apply for error propagation
Default = 1000

[plot] Should results of modeling be plotted? \code{true/FALSE}
@return A list containing the \code{resultarray} which contains the full
result of all simulations on each data window and the \code{parmat} listing
all optimized growth rate and SST parameters used to model d18O in each data
window

seealso Duan, Qingyun, Soroosh Sorooshian, and Vijai Gupta. "Effective and
efficient global optimization for conceptual rainfall runoff models." Water
resources research 28.4 (1992): 1015-1031. https://doi.org/10.1029/91WR02985
@references package dependencies: ggplot2 3.2.1; rtop 0.5.14
Function dependencies: sinreg, d18O_model, growth_model

\doi{10.1029/91WR02985}

"""
function run_model(  # Core function to run the entire model
            dat, # data to run the model on
            dynwindow; # The window vetor resulting from reading in the data 
            transfer_function = "KimONeil97",
            d18Ow = 0,
            T_per = 365, # Temperature sinusoid parameters
            G_per = 365, # Growth sinusoid parameters
            t_int = 1, # Default time interval = 1 day
            t_maxtemp = 182.5, # Default time (day) at which maximum temperature is reached is 182.5 (exactly halfway through the year, or 1st of June)
            SCEUApar = (1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
            sinfit = true, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (true/FALSE)
            MC = 1000, # If errors = true, give the number of iterations for Monte Carlo simulation used in error propagation (default = 1000, if MC = 0 no eror propagation is done)
            plot = FALSE # Should the progress of model fitting be plotted?
    )

    d18Oc = nothing
    d18Oc_err = nothing
    Omod = nothing # Predefine variables to circumvent global variable binding error

    # Prepare data arrays for storage of modeling results
    resultarray = [ # Create array to contain all modeling results of overlapping windows
        zeros(dat), # Replicate matrix of dat + extra columns to contain all variables
        dim = [length(dat.D), length(dynwindow.x) + length(dat[:, 1]), 9] # Six variables, being: Modeled d18O, residuals, Day of the Year, Growth between datapoints, Instantaneous growth rate at datapoint and Modeled temperature
    ]

    parmat = DataFrame(zeros(7, length(dynwindow.x), names=dynwindow.x) # DataFrame in which to store the modeling parameters
    parmat.index = ["T_amp", "T_pha", "T_av", "G_amp", "G_pha", "G_av", "G_skw"]
    
    # Prepare plot to show model progress
    if plot == true 
        p = plot(dat.D, dat.d18Oc, m=:c, yerr=dat.d18Oc_err, xerr=dat.D_err)
        fitplot = ggplot2::ggplot(dat, ggplot2::aes(D, d18Oc)) + # Create a plot showing the fit of each window on the original data, start with a plot of the original data
        ggplot2::geom_point() +
        ggplot2::geom_line() +
        ggplot2::geom_errorbar(ggplot2::aes(ymin = d18Oc - d18Oc_err, # Add error bars on model result (1 SD)
            ymax = d18Oc + d18Oc_err),
            width = dat$D_err) +
        ggplot2::geom_errorbarh(ggplot2::aes(xmin = D - D_err,
            xmax = D + D_err),
            height = 0.05) +
        ggplot2::ggtitle("Plot of d18Oc values vs. depth. Black = data, Red = model, Errorbars = 1 SD")
        plot(fitplot)
    end

    # Estimate growth rate variability and round up to nearest higher magnitude of 10 for conservative boundary
    GRavest = max(diff(dat[dat$YEARMARKER == 1,1])) / 365 # Estimate maximum growth rate from yearmarkers
    GRavmax = 10 ^ (ceiling(log(GRavest, 10))) # Round up to nearest higher magnitude of 10

    # Find tailored range of temperatures from data
    d18Oc_range = range(dat$d18Oc) # Find d18Oc range in data
    if(transfer_function == "KimONeil97"){ # Find temperature range (to be superseded with inverse d18O_model function in later updates)
        T_range = 18.03 * 1000 / (log((d18Oc_range - (0.97002 * rev(range(d18Ow)) - 29.98)) / 1000 + 1) * 1000 + 32.42) - 273.15 # Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
    }else if(transfer_function == "GrossmanKu86"){
        T_range =  20.6 - 4.34 * (d18Oc_range - rev(range(d18Ow)) - 0.2) # Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
    }else{
        print("ERROR: Supplied transfer function is not recognized")
    }
    T_max = max(T_range)
    T_amp_max = 2 * abs(diff(T_range))

    # Collate lower boundaries of parameters
    parl = c(
        T_amp = 0, # Minimum T amplitude in degrees C
        T_pha = 0, # Minimum phase in days
        T_av = -4, # Minimum average T in degrees C
        G_amp = 0, # Minimum seasonal GR range in um/d
        G_pha = 0, # Minimum GR phase in days
        G_av = -1 * GRavmax, # Minimum average GR in um/d.
        G_skw = 0 # Minimum skew factor
    )

    # Collate upper boundaries of parameters
    paru = c(
        T_amp = round(T_amp_max + 0.5, 0), # Maximum T amplitude in degrees C
        T_pha = 365, # Maximum phase in days
        T_av = round(T_max + 0.5, 0), # Maximum average T in degrees C
        G_amp = 2 * GRavmax, # Maximum seasonal GR range in um/d
        G_pha = 365, # Maximum GR phase in days
        G_av = GRavmax, # Maximum average GR in um/d (based on conservative boundaries of YEARMARKER indicators)
        G_skw = 100 # Maximum skew factor 
    )
    
    # Set parameters for SCEUA optimization
    iniflg = SCEUApar[1] # Flag for initial parameter array (default = 1; included)
    ngs = SCEUApar[2] # Number of complexes (sub-populations, default = 25)
    maxn = SCEUApar[3] # Maximum number of function evaluations allowed during optimization (default = 10000)
    kstop = SCEUApar[4] # Maximum number of evolution loops before convergency (default = 5)
    pcento = SCEUApar[5] # Percentage change allowed in function value criterion before stop (default = 0.01)
    peps = SCEUApar[6] # Convergence level for parameter set (difference between parameters required for stop; default = 0.01)

    # Run the model on all windows

    for(i in 1:length(dynwindow$x)){ # Loop over shell record
        print(paste("Processing Datawindow:", i, "of", length(dynwindow$x))) # Keep track of progress
        
        # Isolate year of d18O data based on window data
        Dsam = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 1]
        Osam = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 2]
        if(MC > 0){
            D_err = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 4] # Optional: include error on D
            O_err = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 5] # Optional: include error on d18Oc
        }else{
            D_err = rep(0, dynwindow$y)
            O_err = rep(0, dynwindow$y)
            MC = 0
        }

        if(sinfit){ # If sinusoidal fitting is enabled
            sinlist = sinreg(Dsam, Osam) # Run sinusoidal regression to find initial parameter values
            # Estimate starting parameters from regression results
            O_av_start = sinlist[[1]][1] # Export starting value for annual d18O average
            O_amp_start = sinlist[[1]][2] # Export starting value for d18O amplitude
            O_pha_start = sinlist[[1]][4] %% sinlist[[1]][3] # Estimate position (in depth of the first peak in d18O)
            O_per_start = sinlist[[1]][3] # Export starting value for period in distance domain
        }else{
            O_av_start = mean(Osam) # Estimate starting value for annual d18O average by mean of d18O in record
            O_amp_start = diff(range(Osam)) / 2 # Estimate starting value for d18O amplitude by half the difference between minimum and maximum d18Oc
            O_per_start = diff(range(Dsam)) # Estimate starting period as thickness of isolated year
            O_pha_start = 0.25 * O_per_start # Set starting phase to one quarter of a cycle if it cannot be estimated from sinusoidal regression 
        }

        if(transfer_function == "KimONeil97"){
            T_av_start = 18.03 * 1000 / (1000 * log((O_av_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15  # Estimate mean temperature. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
            T_amp_start = 18.03 * 1000 / (1000 * log((O_av_start - O_amp_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 - T_av_start # Estimate temperature amplitude. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
        }else if(transfer_function == "GrossmanKu86"){
            T_av_start = 20.6 - 4.34 * (O_av_start - mean(d18Ow) - 0.2) # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
            T_amp_start = 20.6 - 4.34 * (O_av_start - O_amp_start - mean(d18Ow) - 0.2) - T_av_start # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
        }else{
            print("ERROR: Supplied transfer function is not recognized")
        }

        O_peak = O_pha_start + Dsam[1] # Find position of d18O peak in distance domain
        T_pha_start = ((O_pha_start - 0.5 * O_per_start) %% O_per_start) / O_per_start * T_per # Estimate position of first peak in temperature (low in d18O) relative to annual cycle (days)
        G_av_start = O_per_start / G_per # Estimate average growth rate in distance/day

        years = 3 # Set default number of years to 3

        # Collate starting parameters
        par0 = c(
            T_amp = T_amp_start,
            T_pha = T_pha_start,
            T_av = T_av_start,

            G_amp = G_av_start / 2, # Start by estimating growth rate changes by half the average
            G_pha = T_pha_start, # Start by estimating that the peak in growth rate coincides with the peak in temperature
            G_av = G_av_start,
            G_skw = 50 # Start with a no skew
        )

        invisible(capture.output( # Suppress the details on converging SCEUA
            sceua_list = rtop::sceua(growth_model,
                par0,
                T_per = T_per,
                G_per = G_per,
                years = years,
                t_int = t_int,
                transfer_function = transfer_function,
                d18Ow = d18Ow,
                Dsam = Dsam,
                Osam = Osam,
                t_maxtemp = t_maxtemp,
                parl,
                paru,
                maxn,
                kstop,
                pcento,
                ngs,
                iniflg = iniflg,
                peps = peps,
                implicit = function(pars){sum(pars[4]/2, pars[6]) < 1} # Make sure that the cumulative GR curve is not located below 0 (if G_av < - G_amp / 2)
                )
            ))

        par1 = sceua_list[[1]] # Eport parameters of final model
        names(par1) = names(par0)

        result = growth_model(par1, T_per, G_per, years, t_int, transfer_function, d18Ow, Dsam, Osam, t_maxtemp, plot = FALSE, MC, D_err, O_err, return = "result") # Calculate the end result of the best fit
        
        if(plot == true){
            fitplot = fitplot + # Add the new model fit to the plot to track progress of the model
                ggplot2::geom_point(data = as.data.frame(result), ggplot2::aes(Dsam, Omod), colour = "red") +
                ggplot2::geom_line(data = as.data.frame(result), ggplot2::aes(Dsam, Omod), colour = "red") 
            print(fitplot)
        }
        
        resultarray[, i + length(dat[1, ]), ] = rbind(matrix(NA, nrow = i - 1, ncol = 9), result[, 3:11], matrix(NA, nrow = length(dat$D) - i - length(result[,1]) + 1, ncol = 9)) # Add results to result array
        parmat[, i] = par1 # Add parameters to parameter array
    }

    # Provide names for all dimensions in the result array
    dimnames(resultarray) = list(
        paste("sample", 1:length(resultarray[, 1, 3])),
        c(colnames(dat), paste("window", 1:length(dynwindow$x))),
        c("Modeled_d18O", "d18O_residuals", "Time_of_year", "Instantaneous_growth_rate", "Modeled temperature", "Modeled_d18O_SD", "Time_of_Year_SD", "Instantaneous_growth_rate_SD", "Modeled_temperature_SD")
    )

    colnames(parmat) = paste("window", 1:length(parmat[1,]))
    return(list(resultarray, parmat))
end