text:
module App
using COAST
using COAST.Zoned
using JuMP
using JSON3
using Genie
using Ipopt
using Genie.Router
using Genie.Requests
using Genie.Renderer.Html
include("test_app_running.jl")
export launchServer

function run_rdaam_dict(n_iter, TSeriesRun, tSeriesRun, data)
  L_dist = data["Letch"]
  U238 = data["U238"]
  Th232 = data["Th232"]
  Ea = data["Ea"]
  rad = data["rad"]
  D0 = data["D0"]
  rmr0 = data["rmr0"]
  alpha = data["alpha"]
  c0Value = data["c0Value"]
  c1Value = data["c1Value"]
  c2Value = data["c2Value"]
  c3Value = data["c3Value"]
  eta_q = data["etaq"]
  psi = data["psi"]
  omega = data["omega"]
  etrap = data["etrap"]

  (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol
  )= conc_to_atoms_per_volume(U238*1e-6,Th232*1e-6)
  L2 = 60
  log10D0L_a2_rdaam = log10(exp(D0)*data["rad"]^2/L2^2)

  he_est = COAST.rdaam_forward_diffusion(alpha,c0Value,c1Value,c2Value,c3Value,rmr0,eta_q,L_dist*1e-4,psi,omega,etrap*1e3,COAST.Rjoules,Ea*1e3,log10D0L_a2_rdaam,n_iter,
     U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,rad*1e-6,reverse(tSeriesRun)...,reverse(TSeriesRun[1:end-1])...)
  # initial guess
  t_val = [1e15]
  for i=1:10
      f0 = (8*(U238_mol*exp(t_val[1]/COAST.τ38)-U238_mol)+
      7*(U235_mol*exp(t_val[1]/COAST.τ35)-U235_mol)) - he_est
      f0_prime = 8*(U238_mol*exp(t_val[1]/COAST.τ38)/COAST.τ38)+7*(U235_mol*exp(t_val[1]/COAST.τ35)/COAST.τ35)
      t_val[1] = t_val[1] - f0/f0_prime
  end
  t_Ma = t_val[1]/(1e6*COAST.sec_in_yrs)

  return t_Ma

end

# string to say hello to users
html_coast_introduction = """
    \n
COAST is a program for thermochronology sensitivity analysis and
     modeling.
 \n
More information can be found on the official
    Github page.

"""

Genie.config.cors_headers["Access-Control-Allow-Credentials"] = "true"
Genie.config.cors_headers["Access-Control-Allow-Methods"] = "GET,HEAD,OPTIONS,POST,PUT"
Genie.config.cors_headers["Access-Control-Allow-Headers"] = "*"
Genie.config.cors_headers["Access-Control-Request-Headers"] = "*"
Genie.config.cors_headers["Access-Control-Allow-Origin"] = "*"
progress_test = [0.0]
run_list = []

function launchServer(progress_test,port)
  # holds run info
  run_queue = Dict("ids"=>[],"progress_percent"=>[],"results"=>[])

   # configuration to make website happy
   Genie.config.run_as_server = true
   Genie.config.server_host = "0.0.0.0"
   Genie.config.server_port = port

  route("/model", method = POST) do
    # extract payload from front-end
    payload = jsonpayload()

    # run appropriate function(s)
    output_string = parse_and_run_payload!(run_queue,payload)

    # make output a string to send back to front-end
    return Genie.Renderer.Json.json(output_string)
  end

  # simple summary for website visits/GET requests
  route("/model", method = GET) do
      print("test model get script \n")
      html_coast_introduction
  end

  # return progress of individual run
  route("/progress", method = POST) do
    payload = jsonpayload()
    id = payload["id"]
    id2find = run_queue["ids"].== id
    progress = run_queue["progress_percent"][id2find]
    return "$(progress)"
  end

  # create other page for testing purposes
  route_test_page_get_post!(progress_test,html_coast_introduction)

  Genie.AppServer.startup()
end

# function to run in COAST
function parse_and_run_payload!(queue,payload)

  # update stats
  push!(queue["progress_percent"],0)
  push!(queue["results"],Dict())
  function_to_COAST = payload["function_to_run"]

  if function_to_COAST=="zonation"
    print("Running zonation script! \n")

    Ea = payload["Ea"]
    Ea = convert(Float64,Ea)
    D0 = payload["D0"]
    D0 = convert(Float64,D0)
    Nt = payload["Nt"]
    Nt = parse(Int64,Nt)
    t_beg = payload["t_beg"]
    t_beg = parse(Float64,t_beg)
    t_end = payload["t_end"]
    t_end = parse(Float64,t_end)
    rad = payload["radius"]
    rad = parse(Float64,rad)
    data = payload["data"]
    vec_grains = convert(Vector{Vector{Vector{Float64}}},data)
    vec_grains = [collect(reduce(hcat,grain_data)') for grain_data in vec_grains]
    outstring = zonation_sensitivity(Ea,D0,rad, vec_grains, t_beg,t_end, Nt)
    return outstring
  elseif function_to_COAST=="store_params"
    userIP = payload["userIP"]
    has_IP = 0
    for (ind, dict_val) in enumerate(run_list)
        if haskey(dict_val,userIP)==1
          has_IP = 1
          run_list[ind][userIP] = payload
        end
    end

    if has_IP==0
      push!(run_list,Dict(userIP=>payload))
    end

    outstring = "params stored"
    return outstring

  ####
  # LOCAL SENSITIVITY/ X-Y-(Z) plot
  # NO ZONATION
  elseif function_to_COAST=="xy"
    # grab data
    tVaried = payload["tData"] # data points that user chooses to vary tT
    TVaried = payload["TData"]
    checkedList = payload["checkedList"]
    diffusion_model = payload["model"]

    # convert to correct units
    sec_in_Ma = 3.1558e13
    c_to_kelvin = 273.15
    tSeries = payload["tSeries"]  # all base tT data points wout variation of tT
    TSeries = payload["TSeries"]
    tSeries = tSeries*sec_in_Ma
    TSeries = TSeries.+c_to_kelvin

    # deal with concurrent users by storing their data separately and accessing here
    userIP = payload["userIP"]
    ind2run = findall(x -> haskey(x,userIP),run_list)
    data = run_list[ind2run[1]][userIP]
    data=JSON3dict_to_dict(data) # necessary

    # get number datapoints X-Y-(Z)
    numberX = parse(Int64,data["numberX"])
    if isempty(data["numberZ"])
      numberZ = 1
    else
      numberZ = parse(Int64,data["numberZ"])
    end

    # CHECK if time or temperature are varied or constant
    names_out = [] # to return to user for plotting
    t_var = [] # time array to (potentially) fill
    T_var = [] # temperature array to (potentially) fill
    t_ind = [] # time point index array to (potentially) fill
    T_ind = [] # temperature point index array to (potentially) fill

    # get any time variation that's been checked by user
    try
      t_var = [t["val"] for t in tVaried if t["check"]==true]
      t_ind = [t["ind"] for t in tVaried if t["check"]==true]
      t_var = [Dict("min"=>t_point[1]["x"]*sec_in_Ma,"main"=>t_point[2]["x"]*sec_in_Ma,"max"=>t_point[3]["x"]*sec_in_Ma)
         for t_point in t_var]
    catch e
      println("no time variation defined!")
    end

    # get any temperature variation that's been checked by user
    try
      T_var = [T["val"] for T in TVaried if T["check"]==true]
      T_ind = [T["ind"] for T in TVaried if T["check"]==true]

      T_var = [Dict("min"=>T_point[1]["y"]+c_to_kelvin,"main"=>T_point[2]["y"]+c_to_kelvin,"max"=>T_point[3]["y"]+c_to_kelvin)
         for T_point in T_var]
    catch e
      print("no temperature variation defined!")
    end
    print(t_var)
    print(T_var)
    # check TYPE of X-Y-(Z) data (t/T/diffusion parameter)
    if length(checkedList)==2
      tTchecked = "onlydiffparams"
      names_out = checkedList
      x = data[checkedList[1]]
      z = data[checkedList[2]]

    elseif (isempty(t_var)) & !(isempty(T_var)) & (length(checkedList)==1)
      tTchecked = "diffparam&T"
      names_out = [checkedList[1],string("temperature at point ",T_ind[1]+1," (°C)")]
      x = data[checkedList[1]]
      z = T_var[1]

    elseif !(isempty(t_var)) & (isempty(T_var)) & length(checkedList)==1
      tTchecked = "t&diffparam"
      names_out = [string("time at point ",t_ind[1]+1," (Ma)"),checkedList[1]]
      x = t_var[1]
      z = data[checkedList[1]]

    elseif (length(t_var)>=2) & (isempty(T_var)) & (length(checkedList)==0)
      tTchecked = "2 t"
      names_out = [string("time at point ",t_ind[1]+1," (Ma)"),string("time at point ",t_ind[2]+1," (Ma)")]
      x = t_var[1]
      z = t_var[2]

    elseif (length(T_var)>=2) & (length(t_var)==0) & (length(checkedList)==0)
      tTchecked = "2 T"
      names_out = [string("temperature at point ",T_ind[1]+1," (Ma)"),string("temperature at point ",T_ind[2]+1," (Ma)")]
      x = T_var[1]
      z = T_var[2]

    elseif (length(checkedList)==0) & !(isempty(t_var)) & !(isempty(T_var))
      tTchecked = "1 t&1 T"
      names_out = [string("time at point ",t_ind[1]+1," (Ma)"),string("temperature at point ",T_ind[1]+1," (°C)")]
      x = t_var[1]
      z = T_var[1]

    # cases where only X-Y plotted (not Z)
    elseif (length(checkedList)==1) & isempty(t_var) & isempty(T_var)
      tTchecked = "1diffparam"
      names_out = checkedList
      x = data[checkedList[1]]
      z = Dict{String,Any}("run"=>["val"])

    elseif (length(checkedList)==0) & (length(t_var)==1) & isempty(T_var)
      tTchecked = "1t"
      names_out = [string("time at point ",t_ind[1]+1," (Ma)")]
      x = t_var[1]
      z = Dict{String,Any}("run"=>["val"])

    elseif (length(checkedList)==0) & (length(T_var)==1) & isempty(t_var)
      tTchecked = "1T"
      names_out = [string("temperature at point ",T_ind[1]+1," (°C)")]
      x = T_var[1]
      z = Dict{String,Any}("run"=>["val"])
    end
    # END OF t/T/diff'n param CHECKING

    # convert from JSON3 to standard Dict
    x = Dict{String,Any}("min"=>x["min"],"main"=>x["main"],"max"=>x["max"])
    if haskey(z,"min")
      z = Dict{String,Any}("min"=>z["min"],"main"=>z["main"],"max"=>z["max"])
    end

    # convert types in case not converted on javascript side - harder to enforce there
    if typeof(x["min"])==String
      x["max"] = parse(Float64,x["max"])
      x["min"] = parse(Float64,x["min"])
    end

    if haskey(z,"min")
      if (typeof(z["min"])==String)
        z["max"] = parse(Float64,z["max"])
        z["min"] = parse(Float64,z["min"])
      end
    end

    # INITIALIZE PLOTTING LOOP VALS
    x["run"] = range(x["min"], stop=x["max"], length=round(Int64,numberX))
    if (haskey(z,"min"))
      z["run"] = range(z["min"], stop=z["max"], length=round(Int64,numberZ))
    end

    lentSeries = length(tSeries)
    n_iter = 50 # for diffusion solver
    n_time_seg_subsegs = 20
    t_Ma = zeros((round(Int64,numberX), round(Int64,numberZ)))

    for (i,xi) in enumerate(x["run"])
        for (j,zi) in enumerate(z["run"])
          print(tTchecked)
          if (tTchecked)=="onlydiffparams"
              data[checkedList[1]]["main"] = string(xi)
              data[checkedList[2]]["main"] = string(zi)
          elseif (tTchecked)=="diffparam&T"
            data[checkedList[1]]["main"] = xi
            TSeries[TVaried[1]["ind"]+1] = zi # convert to kelvin; correct for Javascript indexing
          elseif (tTchecked)=="t&diffparam"
            tSeries[tVaried[1]["ind"]+1] = xi
            data[checkedList[1]]["main"] = zi
          elseif (tTchecked)=="2 t"
              tSeries[tVaried[1]["ind"]+1] = xi
              tSeries[tVaried[2]["ind"]+1] = zi
          elseif (tTchecked)=="2 T"
            print("Tseries is: $TSeries\n")
            Tindex2 = T_ind
            print("T_ind is: $Tindex2\n")
              TSeries[T_ind[1]+1] = xi
              TSeries[T_ind[2]+1] = zi
          elseif (tTchecked)=="1 t&1 T"
              tSeries[tVaried[1]["ind"]+1] = xi
              TSeries[TVaried[1]["ind"]+1] = zi
          elseif (tTchecked)=="1diffparam"
              data[checkedList[1]]["main"] = xi
          elseif (tTchecked)=="1t"
              tSeries[tVaried[1]["ind"]+1] = xi
          elseif (tTchecked)=="1T"
              TSeries[TVaried[1]["ind"]+1] = xi
        end

        # javascript arrays are not necessarily sorted
        idSort = sortperm(tSeries)
        tSeries = tSeries[idSort]
        TSeries = TSeries[idSort]

        # subsample t-T path
        tSeriesRun = Float64[]
        TSeriesRun = Float64[]
        for ind in 1:lentSeries
          if ind<lentSeries
            if isempty(tSeriesRun)
              tSeriesRun = vcat(tSeriesRun,LinRange(tSeries[ind], tSeries[ind+1], n_time_seg_subsegs))
              TSeriesRun = vcat(TSeriesRun,LinRange(TSeries[ind], TSeries[ind+1], n_time_seg_subsegs))
            else
              tSectionNew = LinRange(tSeries[ind], tSeries[ind+1], n_time_seg_subsegs)
              TSectionNew = LinRange(TSeries[ind], TSeries[ind+1], n_time_seg_subsegs)
              tSeriesRun = vcat(tSeriesRun,tSectionNew[2:end])
              TSeriesRun = vcat(TSeriesRun,TSectionNew[2:end])
            end
          end
        end

        print("t is: $tSeries \n")
        print("T is: $TSeries for i = $i j= $j \n")

        # general params conversion
        L_dist = typeof(data["Letch"]["main"])==String ? parse(Float64, data["Letch"]["main"]) : convert(Float64, data["Letch"]["main"])
        U238 = typeof(data["U238"]["main"])==String ? parse(Float64, data["U238"]["main"]) : convert(Float64, data["U238"]["main"])
        Th232 = typeof(data["Th232"]["main"])==String ? parse(Float64, data["Th232"]["main"]) : convert(Float64, data["Th232"]["main"])
        Ea = typeof(data["Ea"]["main"])==String ? parse(Float64, data["Ea"]["main"]) : convert(Float64, data["Ea"]["main"])
        L = typeof(data["rad"]["main"])==String ? parse(Float64, data["rad"]["main"]) : convert(Float64, data["rad"]["main"])
        D0 = typeof(data["D0"]["main"])==String ? parse(Float64, data["D0"]["main"]) : convert(Float64, data["D0"]["main"])
        print("U238 concentration is: $U238 \n")
        print("Th232 concentration is: $Th232 \n")
        L2 = 60 # microns
        log10D0L_a2_rdaam = log10(exp(D0)*L^2/L2^2)

        # handle specific diffusion model param conversion
        if (diffusion_model=="flowers09")
          (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol
          )= conc_to_atoms_per_volume(U238*1e-6,Th232*1e-6)
          rmr0 = typeof(data["rmr0"]["main"])==String ? parse(Float64, data["rmr0"]["main"]) : convert(Float64, data["rmr0"]["main"])
          alpha = typeof(data["alpha"]["main"])==String ? parse(Float64, data["alpha"]["main"]) : convert(Float64, data["alpha"]["main"])
          c0 = typeof(data["c0Value"]["main"])==String ? parse(Float64, data["c0Value"]["main"]) : convert(Float64, data["c0Value"]["main"])
          c1 = typeof(data["c1Value"]["main"])==String ? parse(Float64, data["c1Value"]["main"]) : convert(Float64, data["c1Value"]["main"])
          c2 = typeof(data["c2Value"]["main"])==String ? parse(Float64, data["c2Value"]["main"]) : convert(Float64, data["c2Value"]["main"])
          c3 = typeof(data["c3Value"]["main"])==String ? parse(Float64, data["c3Value"]["main"]) : convert(Float64, data["c3Value"]["main"])
          eta_q = typeof(data["etaq"]["main"])==String ? parse(Float64, data["etaq"]["main"]) : convert(Float64, data["etaq"]["main"])
          psi = typeof(data["psi"]["main"])==String ? parse(Float64, data["psi"]["main"]) : convert(Float64, data["psi"]["main"])
          omega = typeof(data["omega"]["main"])==String ? parse(Float64, data["omega"]["main"]) : convert(Float64, data["omega"]["main"])
          etrap = typeof(data["etrap"]["main"])==String ? parse(Float64, data["etrap"]["main"]) : convert(Float64, data["etrap"]["main"])

        # solve diffusion
          he_est = COAST.rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist*1e-4,psi,omega,etrap*1e3,COAST.Rjoules,Ea*1e3,log10D0L_a2_rdaam,n_iter,
          U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L*1e-6,reverse(tSeriesRun)...,reverse(TSeriesRun[1:end-1])...)
        end
        # find date via newton iterations
        t_val = he_conc_to_date(U238_mol,U235_mol,Th232_mol,he_est)
        t_Ma[i,j] = t_val/(1e6*COAST.sec_in_yrs)
      end
    end

    print(names_out)
    # convert time values so results won't be in seconds
    if occursin("time",names_out[1])
      x["run"] /= sec_in_Ma
    end
    if (haskey(z,"min"))
      if occursin("time",names_out[2])
        z["run"] /= sec_in_Ma
      end
    end

    if occursin("temp",names_out[1])
      x["run"] = x["run"].-c_to_kelvin
    end
    if (haskey(z,"min"))
      if occursin("temp",names_out[2])
        z["run"] = z["run"].-c_to_kelvin
      end
    end

    # return depending on whether X-Y or X-Y-Z plot
    if (haskey(z,"min"))
      outstring = [t_Ma, x["run"], z["run"], names_out]
    else
      outstring = [t_Ma, x["run"], names_out]
    end
    return outstring

  ####
  # GLOBAL SENSITIVITY
  # NO ZONATION
  elseif function_to_COAST=="global_sensitivity"
    samples = payload["samples"]
    n_vars2vary = length(samples)
    data = payload["data"]
    diffusion_model = payload["model"]

    # convert data from weird JSON3 dictionary format to standard dict
    data=JSON3dict_to_dict(data)

    n_run = typeof(data["numberX"])==String ? parse(Int64, data["numberX"]) : convert(Int64, data["numberX"])
    tSeries = data["tSeries"]
    TSeries = data["TSeries"]
    print("TSeries is: $TSeries\n")
    print("tSeries is: $tSeries\n")
    tSeries = tSeries*3.1558e7*1e6
    TSeries = TSeries.+ 273.15
    lentSeries = length(tSeries)

    data_to_run = Dict{String, Float64}()

    data_to_run["U238"] = typeof(data["U238"]["main"])==String ? parse(Float64, data["U238"]["main"]) : convert(Float64, data["U238"]["main"])
    data_to_run["Th232"] = typeof(data["Th232"]["main"])==String ? parse(Float64, data["Th232"]["main"]) : convert(Float64, data["Th232"]["main"])
    data_to_run["Ea"] = typeof(data["Ea"]["main"])==String ? parse(Float64, data["Ea"]["main"]) : convert(Float64, data["Ea"]["main"])
    data_to_run["rad"] = typeof(data["rad"]["main"])==String ? parse(Float64, data["rad"]["main"]) : convert(Float64, data["rad"]["main"])
    data_to_run["D0"] = typeof(data["D0"]["main"])==String ? parse(Float64, data["D0"]["main"]) : convert(Float64, data["D0"]["main"])

    # diffusion model-specific conversion
    if (diffusion_model=="flowers09")
      data_to_run["Letch"] = typeof(data["Letch"]["main"])==String ? parse(Float64, data["Letch"]["main"]) : convert(Float64, data["Letch"]["main"])
      data_to_run["rmr0"] = typeof(data["rmr0"]["main"])==String ? parse(Float64, data["rmr0"]["main"]) : convert(Float64, data["rmr0"]["main"])
      data_to_run["alpha"] = typeof(data["alpha"]["main"])==String ? parse(Float64, data["alpha"]["main"]) : convert(Float64, data["alpha"]["main"])
      data_to_run["c0Value"] = typeof(data["c0Value"]["main"])==String ? parse(Float64, data["c0Value"]["main"]) : convert(Float64, data["c0Value"]["main"])
      data_to_run["c1Value"] = typeof(data["c1Value"]["main"])==String ? parse(Float64, data["c1Value"]["main"]) : convert(Float64, data["c1Value"]["main"])
      data_to_run["c2Value"] = typeof(data["c2Value"]["main"])==String ? parse(Float64, data["c2Value"]["main"]) : convert(Float64, data["c2Value"]["main"])
      data_to_run["c3Value"] = typeof(data["c3Value"]["main"])==String ? parse(Float64, data["c3Value"]["main"]) : convert(Float64, data["c3Value"]["main"])
      data_to_run["etaq"] = typeof(data["etaq"]["main"])==String ? parse(Float64, data["etaq"]["main"]) : convert(Float64, data["etaq"]["main"])
      data_to_run["psi"] = typeof(data["psi"]["main"])==String ? parse(Float64, data["psi"]["main"]) : convert(Float64, data["psi"]["main"])
      data_to_run["omega"] = typeof(data["omega"]["main"])==String ? parse(Float64, data["omega"]["main"]) : convert(Float64, data["omega"]["main"])
      data_to_run["etrap"] = typeof(data["etrap"]["main"])==String ? parse(Float64, data["etrap"]["main"]) : convert(Float64, data["etrap"]["main"])
    end

    n_time_seg_subsegs = 20
    n_iter = 50 # for diffusion model
    date_Ma = zeros(Float64,n_run)
    for i=1:n_run
      for j=1:n_vars2vary
        key = payload["key_list"][j]
        if occursin("tData", key)
          tData_ind = parse(Int64,key[6])+1
          tSeries[tData_ind] = samples[j][i]*3.1558e13 # correct from javascript indexing
        elseif occursin("TData", key)
          TData_ind = parse(Int64,key[6])+1 # jscript is zero indexed
          TSeries[TData_ind] = samples[j][i].+273.15
        else
          data[key]["main"] = samples[j][i]
        end
      end
      # javascript arrays are not necessarily sorted
      print("pre-sorted t is: $tSeries \n")
      print("pre-sorted T is: $TSeries for i = $i\n")

      idSort = sortperm(tSeries)
      tSeries = tSeries[idSort]
      TSeries = TSeries[idSort]

      tSeriesRun = Float64[]
      TSeriesRun = Float64[]

      for ind in 1:lentSeries
        if ind<lentSeries
          if isempty(tSeriesRun)
            tSeriesRun = vcat(tSeriesRun,LinRange(tSeries[ind], tSeries[ind+1], n_time_seg_subsegs))
            TSeriesRun = vcat(TSeriesRun,LinRange(TSeries[ind], TSeries[ind+1], n_time_seg_subsegs))
          else
            tSectionNew = LinRange(tSeries[ind], tSeries[ind+1], n_time_seg_subsegs)
            TSectionNew = LinRange(TSeries[ind], TSeries[ind+1], n_time_seg_subsegs)
            tSeriesRun = vcat(tSeriesRun,tSectionNew[2:end])
            TSeriesRun = vcat(TSeriesRun,TSectionNew[2:end])
          end
        end
      end

      print("t is: $tSeries \n")
      print("T is: $TSeries for i = $i\n")

      # conversion
      U238 = typeof(data["U238"]["main"])==String ? parse(Float64, data["U238"]["main"]) : convert(Float64, data["U238"]["main"])
      Th232 = typeof(data["Th232"]["main"])==String ? parse(Float64, data["Th232"]["main"]) : convert(Float64, data["Th232"]["main"])
      Ea = typeof(data["Ea"]["main"])==String ? parse(Float64, data["Ea"]["main"]) : convert(Float64, data["Ea"]["main"])
      L = typeof(data["rad"]["main"])==String ? parse(Float64, data["rad"]["main"]) : convert(Float64, data["rad"]["main"])
      D0 = typeof(data["D0"]["main"])==String ? parse(Float64, data["D0"]["main"]) : convert(Float64, data["D0"]["main"])

      L2 = 60

      print("U238 concentration is: $U238 \n")
      print("Th232 concentration is: $Th232 \n")
      log10D0L_a2_rdaam = log10(exp(D0)*L^2/L2^2)

      # diffusion model-specific conversion
      if (diffusion_model=="flowers09")
        (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol
        )= conc_to_atoms_per_volume(U238*1e-6,Th232*1e-6)
        L_dist = typeof(data["Letch"]["main"])==String ? parse(Float64, data["Letch"]["main"]) : convert(Float64, data["Letch"]["main"])
        rmr0 = typeof(data["rmr0"]["main"])==String ? parse(Float64, data["rmr0"]["main"]) : convert(Float64, data["rmr0"]["main"])
        alpha = typeof(data["alpha"]["main"])==String ? parse(Float64, data["alpha"]["main"]) : convert(Float64, data["alpha"]["main"])
        c0 = typeof(data["c0Value"]["main"])==String ? parse(Float64, data["c0Value"]["main"]) : convert(Float64, data["c0Value"]["main"])
        c1 = typeof(data["c1Value"]["main"])==String ? parse(Float64, data["c1Value"]["main"]) : convert(Float64, data["c1Value"]["main"])
        c2 = typeof(data["c2Value"]["main"])==String ? parse(Float64, data["c2Value"]["main"]) : convert(Float64, data["c2Value"]["main"])
        c3 = typeof(data["c3Value"]["main"])==String ? parse(Float64, data["c3Value"]["main"]) : convert(Float64, data["c3Value"]["main"])
        eta_q = typeof(data["etaq"]["main"])==String ? parse(Float64, data["etaq"]["main"]) : convert(Float64, data["etaq"]["main"])
        psi = typeof(data["psi"]["main"])==String ? parse(Float64, data["psi"]["main"]) : convert(Float64, data["psi"]["main"])
        omega = typeof(data["omega"]["main"])==String ? parse(Float64, data["omega"]["main"]) : convert(Float64, data["omega"]["main"])
        etrap = typeof(data["etrap"]["main"])==String ? parse(Float64, data["etrap"]["main"]) : convert(Float64, data["etrap"]["main"])

        he_est = COAST.rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist*1e-4,psi,omega,etrap*1e3,COAST.Rjoules,Ea*1e3,log10D0L_a2_rdaam,n_iter,
        U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L*1e-6,reverse(tSeriesRun)...,reverse(TSeriesRun[1:end-1])...)
      end
      t_val = he_conc_to_date(U238_mol,U235_mol,Th232_mol,he_est)
      date_Ma[i] = t_val/(1e6*COAST.sec_in_yrs)

    end
    print(date_Ma)
    return date_Ma

  # single load run
  elseif function_to_COAST=="single_grain"
    Etrap = payload["Etrap"]
    Etrap = parse(Float64,Etrap)
    alpha = payload["alpha"]
    alpha = parse(Float64, alpha)
    c0 = payload["c0"]
    c0 = parse(Float64, c0)
    c1 = payload["c1"]
    c1 = parse(Float64, c1)
    c2 = payload["c2"]
    c2 = parse(Float64, c2)
    c3 = payload["c3"]
    c3 = parse(Float64, c3)
    rmr0=payload["rmr0"]
    rmr0 = parse(Float64, rmr0)
    eta_q = payload["eta_q"]
    eta_q = parse(Float64, eta_q)
    L_dist = payload["L_dist"]
    L_dist = parse(Float64, L_dist)
    psi = payload["psi"]
    psi = parse(Float64, psi)
    omega = payload["omega"]
    omega = parse(Float64, omega)
    E_L = payload["E_L"]
    E_L = parse(Float64, E_L)
    D0L_a2 = payload["D0L_a2"]
    D0L_a2 = parse(Float64, D0L_a2)
    rad = payload["rad"]
    rad = parse(Float64, rad)
    u38 = payload["u38"]
    u38 = parse(Float64, u38)
    th32 = payload["th32"]
    th32 = parse(Float64, th32)
    raw_tT = JSON3.read(payload["tT"])

    times = [data_pt["x"]*COAST.sec_in_yrs for data_pt in raw_tT]

    T = [data_pt["z"]+273.15 for data_pt in raw_tT]

    # times = payload["times"]
    # times = convert(Vector{Float64},times)
    # T = payload["T"]
    # T = convert(Vector{Float64},T)
    if (
      isnan(rmr0) ||
      isnan(c0) ||
      isnan(c1) ||
      isnan(c2) ||
      isnan(c3) ||
      isnan(alpha) ||
      isnan(eta_q) ||
      isnan(L_dist)||
      isnan(psi) ||
      isnan(omega) ||
      isnan(Etrap)
        )
      rmr0 = 1.0
      c0 = 1.0
      c1 = 1.0
      c2 = 1.0
      c3 = 1.0
      alpha = 1.0
      eta_q = 1.0
      L_dist = 1.0
      psi = 0.0
      omega = 0.0
      Etrap = 0.0
    end

    outstring = single_grain(rmr0,c0,c1,c2,c3,alpha,eta_q,L_dist,psi,omega,Etrap,rad,u38,th32,E_L,D0L_a2,times,T)
    return outstring
  end

  return "Passing JSON payload to COAST successful!"
end

# TYPE CONVERSION HELPER FUNCS
function JSON3dict_to_dict(data)
  data_new = Dict()
  for (key, value) in data
    if (typeof(value)==JSON3.Object{Base.CodeUnits{UInt8, String}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}})
      value_new = Dict{String,Any}()
      for (nest_key, nest_value) in value
        value_new[string(nest_key)] = nest_value
      end
      value = deepcopy(value_new)
    end
    data_new[string(key)] = value
  end
  data = deepcopy(data_new)
  return data
end

# convert concentration to date using Newton iters
function he_conc_to_date(U238_mol,U235_mol,Th232_mol,he_est)
  t_val = [1e15] # initial value
  for i=1:20 # tested this val (10 usually sufficient)
    f0 = (8*(U238_mol*exp(t_val[1]/COAST.τ38)-U238_mol)+
    7*(U235_mol*exp(t_val[1]/COAST.τ35)-U235_mol)+
    6*(Th232_mol*exp(t_val[1]/COAST.τ32)-Th232_mol)) - he_est
    f0_prime = 8*(U238_mol*exp(t_val[1]/COAST.τ38)/COAST.τ38)+7*(U235_mol*exp(t_val[1]/COAST.τ35)/COAST.τ35)+6*(Th232_mol*exp(t_val[1]/COAST.τ32)/COAST.τ32)
    t_val[1] = t_val[1] - f0/f0_prime
  end
  return t_val[1]
end

function zonation_n_times_forward(distance,L,n_t,Ea,D0,U38_Pb06,sigU38_Pb06,dr,
  tmax,tmin,Tbound_top,Tbound_bot)
  print("running_zonation_n_times_forward \n")
  tmax = parse(Float64,tmax) # in yrs
  tmin = parse(Float64,tmin) # in yrs
  n_t  = parse(Int64,n_t)
  Ea = parse(Float64,Ea)
  D0 = parse(Float64,D0)
  dr = parse(Float64,dr)
  L = parse(Float64,L)
  distance = JSON3.read(distance, Vector{Float64})./1e6
  sigU38_Pb06 = JSON3.read(sigU38_Pb06, Vector{Vector{Float64}})
  Tbound_top = JSON3.read(Tbound_top, Vector{Float64})
  Tbound_bot = JSON3.read(Tbound_bot, Vector{Float64})
  U38_Pb06 = JSON3.read(U38_Pb06, Vector{Vector{Float64}})
  # rows are permutations of sample sets
  n_constraints = length(U38_Pb06[1])

  t = collect(LinRange(tmax*COAST.sec_in_yrs,tmin*COAST.sec_in_yrs,n_t))
  upper_bound_arr = 420.0*ones(n_t)
  upper_bound_arr[Tbound_top.>-1.0] = Tbound_top[Tbound_top.>-1.0]
  lower_bound_arr = 273.0*ones(n_t)
  lower_bound_arr[Tbound_bot.>-1.0] = Tbound_bot[Tbound_bot.>-1.0]
  model = Model(Ipopt.Optimizer)
  JuMP.set_optimizer_attributes(model,"print_level"=>2,"tol"=>1e2,"linear_solver"=>"mumps","print_timing_statistics"=>"no")
  register_objective_function!(n_t,model)
  @NLparameter(model, pa[i=1:n_constraints] == U38_Pb06[1][i])
  start_points = 0.5.*ones(Float64,n_t)
  (T,set_dev) = define_variables!(n_t,n_constraints,model,start_points,upper_bound=upper_bound_arr,lower_bound=lower_bound_arr)
  register_forward_model_zonation!(n_t,model)

  return lower_bound_arr
end

function zonation_sensitivity(Ea,D0,rad,vec_grains,t_beg,t_end,Nt)

  ngrains = length(vec_grains)
  data = Dict((1:ngrains).=>vec_grains)
  U238 = 1.0

  # convert to meters
  for val in values(data)
    val[:,1]/=1e6
  end
  rad /= 1e6
  Nx = zeros(Int64,ngrains) # preallocate for number of nodes in grain
  for (key, value) in data
    dr = diff(value[1:2,1])[1]
    r_added_min = maximum(value[:,1])
    r_added = collect(r_added_min:dr:(rad+dr))
    n_added = length(r_added)
    # fill in all other values as well
    Pb06U38_grain_center=fill(value[end,2],n_added)
    sigPb06U38_grain_center=fill(value[end,3],n_added)
    data_grain_center = hcat(r_added,Pb06U38_grain_center,sigPb06U38_grain_center)
    data[key] = vcat(value,data_grain_center)
    data[key] = vcat([0.0 0.0 0.0],data[key])
    Nx[key] = size(data[key],1)
  end

  # correct model to model end date
  end_model_Ma = t_end
  Pb_corr_end_model = 1.0*(exp(COAST.lambda_38*COAST.sec_in_yrs*end_model_Ma*1e6))-1.0
  # pull out individual components from data for clarity
  U238_new = Dict{Int64,Vector{Float64}}()
  r = Dict{Int64,Vector{Float64}}()
  Pbmeas = Dict{Int64,Vector{Float64}}()
  sigmeas = Dict{Int64,Vector{Float64}}()
  for (key,value) in data
      value[2:end,2] = value[2:end,2].-Pb_corr_end_model
      push!(r,key=>value[:,1])
      push!(U238_new,key=>U238*ones(length(value[:,1])))
      push!(Pbmeas,key=>reverse(value[:,2]))
      push!(sigmeas,key=>reverse(value[:,3]))
  end

  t = collect(LinRange(t_beg,end_model_Ma,Nt).*3.1558e7*1e6)
  dt = -maximum(diff(t))/(3.1558e7*1e6)
  T = collect(LinRange(800.0,673.0,Nt).+273.15)

  upper_bound_arr = 1100.0*ones(Nt)
  lower_bound_arr = 673.0*ones(Nt)
  lower_bound_arr[1] = 1100.0

  empty_dict = Dict(keys(U238_new).=>0.0.*values(U238_new))
  zon_loglike_objective = Zoned.create_zon_loglike_objective(Ea,COAST.Rjoules,D0,U238_new,empty_dict,empty_dict,Nt,Nx,t,r,Pbmeas,sigmeas,ngrains)
  grad_loglike_objective = Zoned.create_grad_loglike_objective(Ea,COAST.Rjoules,D0,U238_new,empty_dict,empty_dict,Nt,Nx,t,r,Pbmeas,sigmeas,ngrains)

  n_constr = Nt-1
  upper_bound_constr = 100*ones(n_constr)
  lower_bound_constr = 0.0*ones(n_constr)

  prob2 = createProblem(
      Nt, # should match length T_L, T_U otherwise julia may crash
      lower_bound_arr, # lower bound on temperature
      upper_bound_arr, # upper bound
      n_constr,
      lower_bound_constr, # lower bound of Pb concentration; i.e. constraints
      upper_bound_constr,
      Nt*n_constr, # number elements in constraint jacobian
      Nt*n_constr, # number elements in constraint hessian; not needed in approximation
      zon_loglike_objective, # objective function; sum smoothness (multiobjective not yet possible)
      Zoned.decrease_constraint,
      grad_loglike_objective,
      Zoned.jac_decrease_constraint,
  )
  addOption(prob2, "hessian_approximation", "limited-memory") # hessian rarely (if ever) called; approximation is ok
  addOption(prob2, "tol", 1e-4) # smoothness not as much of an issue here
  addOption(prob2,"print_level",5)
  prob2.x = copy(T)
  print(r)
  status=solveProblem(prob2)

  constraint_zon = create_constraint_zon(Ea,COAST.Rjoules,D0,U238_new[1],0.0,0.0,Nt,Nx[1],t,r[1])
  Pbs = zeros(Nx[1]) # preallocate
  constraint_zon(prob2.x,Pbs)

  constraint_zon = create_constraint_zon(Ea,COAST.Rjoules,D0,U238_new[2],0.0,0.0,Nt,Nx[2],t,r[2])
  Pbs2 = zeros(Nx[2]) # preallocate
  constraint_zon(prob2.x,Pbs2)
  if status==0
    out = prob2.x
  else
    out = "COAST could not find optimal solution!"
  end

  return out
end

function single_grain(rmr0,c0,c1,c2,c3,alpha,eta_q,L_dist,psi,omega,Etrap,L,u38,th32,E_L,D0L_a2,times,T)
  # preprocess concentrations
  (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = ppm_to_atoms_per_volume(u38,th32,density=3.20)
  L2 = 60*1e-4
  n_iter = 60
  log10D0L_a2_rdaam = log10(D0L_a2*L^2/L2^2)
  mass_he_t1 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,psi,omega,Etrap,COAST.Rjoules,E_L,log10D0L_a2_rdaam,n_iter,
                   U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L,times...,T...)

  outstring = string(mass_he_t1)
  #outstring = "ran successfully"
  return outstring
end

end # end module


text:
defaulteigalg(Ψ::UniformCMPS) = Arnoldi(; krylovdim = min(256, virtualdim(Ψ)^2))
defaultlinalg(Ψ::UniformCMPS) = GMRES(; krylovdim = min(256, virtualdim(Ψ)^2))

defaulteigalg(Ψ::FourierCMPS) =
    Arnoldi(; krylovdim = min(256, (10*nummodes(Ψ.Q)+1)*virtualdim(Ψ)^2))
defaultlinalg(Ψ::FourierCMPS) =
    GMRES(; krylovdim = min(256, (10*nummodes(Ψ.Q)+1)*virtualdim(Ψ)^2))

# CMPS environements
function leftenv(Ψ::InfiniteCMPS, ρ₀ = one(Ψ.Q);
                    eigalg = defaulteigalg(Ψ),
                    linalg = nothing, # ignored
                    kwargs...)

    eigsort = EigSorter(x->(abs(div(imag(x), pi/period(Ψ))), -real(x)))
    let TL = LeftTransfer(Ψ)
        _, ρs, λs, info = schursolve(ρ₀, 1, eigsort, eigalg) do x
            y = -∂(x) + TL(x; tol = eigalg.tol/10, kwargs...)
            return truncate!(y; tol = eigalg.tol/10, kwargs...)
        end
        λ, ρ = λs[1]/2, ρs[1]
        imag(λ) <= max(info.normres[1], defaulttol(λ)) ||
            @warn "Largest eigenvalue of transfer matrix not real? $λ"
        ρ = ρ + ρ'
        ρ = rmul!(ρ, 1/(norm(ρ)*sign(tr(ρ[0]))))
        ρ = truncate!(ρ; tol = eigalg.tol/100, kwargs...)
        res = -∂(ρ) + TL(ρ) - (2*λ)*ρ
        newinfo = ConvergenceInfo(info.converged, res, norm(res), info.numiter, info.numops)
        return ρ, real(λ), newinfo
    end
end

function leftenv!(Ψ::InfiniteCMPS, ρ₀ = one(Ψ.Q); kwargs...)
    ρ, λ, info = leftenv(Ψ, ρ₀; kwargs...)
    Q0 = Ψ.Q[0]
    view(Q0, diagind(Q0)) .-= λ
    return ρ, info
end

function rightenv(Ψ::InfiniteCMPS, ρ₀ = one(Ψ.Q);
                    eigalg = defaulteigalg(Ψ),
                    linalg = nothing, # ignored
                    kwargs...)

    eigsort = EigSorter(x->(abs(div(imag(x), pi/period(Ψ))), -real(x)))
    let TR = RightTransfer(Ψ)
        _, ρs, λs, info = schursolve(ρ₀, 1, eigsort, eigalg) do x
                y = ∂(x) + TR(x; kwargs...)
                return truncate!(y; tol = eigalg.tol/100, kwargs...)
            end
        λ, ρ = λs[1]/2, ρs[1]
        imag(λ) <= max(info.normres[1], defaulttol(λ)) ||
            @warn "Largest eigenvalue of transfer matrix not real? $λ"
        ρ = ρ + ρ'
        ρ = rmul!(ρ, 1/(norm(ρ)*sign(tr(ρ[0]))))
        ρ = truncate!(ρ; tol = eigalg.tol/100, kwargs...)
        res = ∂(ρ) + TR(ρ) - (2*real(λ))*ρ
        newinfo = ConvergenceInfo(info.converged, res, norm(res), info.numiter, info.numops)
        return ρ, real(λ), newinfo
    end
end

function rightenv!(Ψ::InfiniteCMPS, ρ₀ = one(Ψ.Q); kwargs...)
    ρ, λ, info = rightenv(Ψ, ρ₀; kwargs...)
    Q0 = Ψ.Q[0]
    view(Q0, diagind(Q0)) .-= λ
    return ρ, info
end

function environments!(Ψ::InfiniteCMPS, ρL₀ = one(Ψ.Q), ρR₀ = one(Ψ.Q); kwargs...)
    ρL, λL, infoL = leftenv(Ψ, ρL₀; kwargs...)
    ρR, λR, infoR = rightenv(Ψ, ρR₀; kwargs...)
    isapprox(λL, λR; atol = defaulttol(λL), rtol = defaulttol(λL)) ||
        @warn "Different left and right normalization: $λL versus $λR"

    Q0 = Ψ.Q[0]
    Q0d = view(Q0, diagind(Q0))
    Q0d .-= (λL+λR)/2

    Z = dot(ρL, ρR)
    imag(Z) < defaulttol(Z) && real(Z) > 0 ||
        @warn "Incompatible fixed point environments: Z = ⟨ρL|ρR⟩ = $Z"

    ρL = rmul!(ρL, 1/sqrt(Z))
    ρR = rmul!(ρR, 1/sqrt(Z))

    return ρL, ρR, infoL, infoR
end

# Environments of Hamiltonian with CMPS
const InfiniteCMPSData = Tuple{InfiniteCMPS,PeriodicMatrixFunction,PeriodicMatrixFunction}

function leftenv!(H::LocalHamiltonian, Ψ::InfiniteCMPS, HL₀ = zero(Ψ.Q); kwargs...)
    domain(H) == domain(Ψ) || throw(DomainMismatch())
    ρL, ρR, infoL, infoR = environments!(Ψ; kwargs...)
    return leftenv(H, (Ψ,ρL,ρR), HL₀; kwargs...)
end

function rightenv!(H::LocalHamiltonian, Ψ::InfiniteCMPS, HR₀ = zero(Ψ.Q); kwargs...)
    domain(H) == domain(Ψ) || throw(DomainMismatch())
    ρL, ρR, infoL, infoR = environments!(Ψ; kwargs...)
    return rightenv(H, (Ψ,ρL,ρR), HR₀; kwargs...)
end

# assumes Ψ is normalized and ⟨ρL|ρR⟩ = 1
function leftenv(H::LocalHamiltonian, Ψρs::InfiniteCMPSData, HL₀ = nothing;
                    eigalg = nothing, # ignored
                    linalg = defaultlinalg(Ψρs[1]),
                    kwargs...)

    Ψ, ρL, ρR = Ψρs
    domain(H) == domain(Ψ) || throw(DomainMismatch())

    hL = leftreducedoperator(H.h, Ψ, ρL)
    eL = real(localdot(hL, ρR))
    EL = real(dot(hL, ρR))
    hL = axpy!(-EL, ρL, hL)

    if isnothing(HL₀)
        HL₀ = zero(Ψρs[1].Q)
    else
        HL₀ = HL₀ - ρL * dot(HL₀, ρR)
    end
    let TL = LeftTransfer(Ψ)
        HL, infoL = linsolve(hL/norm(hL), HL₀, linalg) do x
            y = ∂(x) - TL(x; kwargs...)
            y = axpy!(dot(ρR, x), ρL, y)
            truncate!(y; tol = linalg.tol/100, kwargs...)
        end
        HL = rmul!(HL + HL', 0.5*norm(hL))
        # truncate!(HL; tol = linalg.tol/100, kwargs...)
        res = hL - (∂(HL)-TL(HL))
        infoL = ConvergenceInfo(infoL.converged, res, norm(res), infoL.numiter, infoL.numops)
        return HL, EL, eL, hL, infoL
    end
end

function rightenv(H::LocalHamiltonian, Ψρs::InfiniteCMPSData, HR₀ = nothing;
                    eigalg = nothing, # ignored
                    linalg = defaultlinalg(Ψρs[1]),
                    kwargs...)

    Ψ, ρL, ρR = Ψρs
    domain(H) == domain(Ψ) || throw(DomainMismatch())

    hR = rightreducedoperator(H.h, Ψ, ρR)
    eR = real(localdot(ρL, hR))
    ER = real(dot(ρL, hR))
    hR = axpy!(-ER, ρR, hR)

    if isnothing(HR₀)
        HR₀ = zero(Ψρs[1].Q)
    else
        HR₀ = HR₀ - ρR * dot(ρL, HR₀)
    end
    let TR = RightTransfer(Ψ)
        HR, infoR = linsolve(hR/norm(hR), HR₀, linalg) do x
            y = -∂(x) - TR(x; kwargs...)
            y = axpy!(dot(ρL, x), ρR, y)
            truncate!(y; tol = linalg.tol/100, kwargs...)
        end
        HR = rmul!(HR + HR', 0.5*norm(hR))
        # res = hR - (-∂(HR)-TR(HR))
        # truncate!(HR; tol = linalg.tol/100, kwargs...)
        res = hR - (-∂(HR)-TR(HR))
        infoR = ConvergenceInfo(infoR.converged, res, norm(res), infoR.numiter, infoR.numops)
        return HR, ER, eR, hR, infoR
    end
end

function environments!(H::LocalHamiltonian, Ψ::InfiniteCMPS; kwargs...)
    ρL, ρR, infoρL, infoρR = environments!(Ψ; kwargs...)
    HL, HR, hL, hR, E, e, infoHL, infoHR = environments(H,(Ψ,ρL,ρR); kwargs...)
    return HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR
end

function environments(H::LocalHamiltonian, Ψρs::InfiniteCMPSData,
                        HL₀ = nothing, HR₀ = nothing; kwargs...)
    HL, EL, eL, hL, infoHL = leftenv(H, Ψρs, HL₀; kwargs...)
    HR, ER, eR, hR, infoHR = rightenv(H, Ψρs, HR₀; kwargs...)
    eL ≈ eR ||
        @warn "non-matching energy from left and right environments"
    E = (EL+ER/2)
    e = rmul!(eL + eR, 1//2)
    return HL, HR, hL, hR, E, e, infoHL, infoHR
end

ext: jl


text:
using Reddit
using Test

creds = Credentials("id", "secret", "agent", "name", "password")

@testset "Credentials" begin
    @test creds.id == "id"
    @test creds.secret == "secret"
    @test creds.useragent == "agent"
    @test creds.username == "name"
    @test creds.password == "password"
end


text:
# functions that extract the mininum data necessary from the mesh to do a
# DG-type shock capturing scheme.

#TODO: is this unused?
"""
  This function gets the list of elements (and other required data) to do
  shock capturing with a DG-type dissipation.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * sensor

  **Inputs/Outputs**

   * shockmesh: a [`ShockedElement`](@ref) object.  Will be reset before
                writing a new set of shocked elements
"""
function getShockCapturingMesh(mesh::AbstractMesh, sbp::AbstractOperator,
                               eqn::EulerData, opts, sensor::AbstractShockSensor,
                               shockmesh::ShockedElements)

  reset(shockmesh)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)

    add_el = isShockElement(eqn.params, sbp, sensor, q_i, jac_i)

    if add_el
      push!(shockmesh, i)
    end
  end

  completeShockElements(shockmesh)

  return nothing
end

"""
  This function resets the `ShockElements` object to its (more-or-less) initial
  state while keeping the sizes of the existing arrays.  After calling reset
  a new set of shocked elements can be recorded, hopfully with fewer array
  resizes that the first time.

  **Inputs**

   * data: `ShockedElements`
"""
function reset(data::ShockedElements)

  fill!(data.elnums_mesh, 0)
  fill!(data.elnums_all, 0)
  data.numShock = 0
  data.numNeighbor = 0
  data.numInterfaces = 0
  data.numEl = 0

  data.idx_all = 1
  data.idx_if = 1
  data.idx_b = 1
  fill!(data.idx_sf, 1)

  return nothing
end

import Base.push!

"""
  Adds a new element to the list of elements that have a shock in them

  **Inputs**

   * data: `ShockedElements`
   * elnum: the element number
"""
@inline function push!(data::ShockedElements, elnum::Int)

  idx = data.idx_all
  sz = data.sz_all
  if idx > sz
    sz = max(sz*2, 8)
    resize!(data.elnums_all, sz)
    data.sz_all = sz
  end

  data.elnums_all[idx] = elnum
  data.elnums_mesh[elnum] = idx

  data.idx_all += 1
end

"""
  Adds a new element to the list of neighboring elements.  Internal function

  **Inputs**

   * data: `ShockedElements`
   * elnum: the new element number
   * idx: the index to insert the new element number
"""
@inline function push_neighbor(data::ShockedElements, elnum::Integer)

  idx = data.idx_all
  sz = data.sz_all

  if idx > sz
    sz = max(2*sz, 8)
    resize!(data.elnums_all, sz)
    data.sz_all = sz
  end

  data.elnums_all[idx] = elnum
  data.elnums_mesh[elnum] = idx
  data.idx_all += 1

  return nothing
end

"""
  Internal function for adding element to list of shared elements (ie.
  elements that live on other processes that share a face with an element
  on this process that has a shock in it.

  **Inputs**

   * data: `ShockedElements`
   * elnum: new element number
"""
@inline function push_shared(data::ShockedElements, elnum::Integer)

  idx = data.idx_all
  sz = data.sz_all

  if idx > sz
    sz = max(2*sz, 8)
    resize!(data.elnums_all, sz)
    data.sz_all = sz
  end

  data.elnums_all[idx] = elnum
  data.elnums_mesh[elnum] = idx
  data.idx_all += 1
  data.numShared[end] += 1

  return nothing
end

"""
  Internal function for pushing a new interface

  **Inputs**

   * data: `ShockElements`
   * iface: a `RelativeInterface`
   * idx: the current index in data.ifaces
   * sz: the current size of data.ifaces
"""
@inline function push_iface(data::ShockedElements, iface::RelativeInterface)

  idx = data.idx_if
  sz = data.sz_if

  if idx > sz
    sz = max(2*sz, 8)
    resize!(data.ifaces, sz)
    data.sz_if = sz
  end

  data.ifaces[idx] = iface
  data.idx_if += 1

  return nothing
end

@inline function push_sharediface(data::ShockedElements, iface::RelativeInterface)

  idx = data.idx_sf[end]
  sz = data.sz_sf[end]
  if idx > sz
    sz = max(2*sz, 8)
    resize!(data.shared_interfaces[end], sz)
    data.sz_sf[end] = sz
  end

  data.shared_interfaces[end][idx] = iface
  data.numSharedInterfaces[end] += 1
  data.idx_sf[end] += 1

  return nothing
end

"""
  Internal function for pushing new boundary

  **Inputs**

   * data: `ShockElements`
   * bndry: a `RelativeBoundary`

"""
@inline function push_bndry(data::ShockedElements, bndry::RelativeBoundary)

  idx = data.idx_b
  sz = data.sz_b

  if idx > sz
    sz = max(2*sz, 8)
    resize!(data.bndryfaces, sz)
    data.sz_b = sz
  end

  data.bndryfaces[idx] = bndry
  data.idx_b += 1

  return nothing
end

"""
  Replaces the element numbers for an interface

  **Inputs**

   * iface: an `Interface`
   * elementL: the new left element number
   * elementR: the new right elenent number

  **Outputs**

   * iface: a new `Interface` object
"""
function replace_interface(iface::Interface, elementL::Integer, elementR::Integer)
  return Interface(elementL, elementR, iface.faceL, iface.faceR, iface.orient)
end

"""
  Similar to [`replace_interface`](@ref), but for `Boundary` objects.
"""
function replace_boundary(bndry::Boundary, elnum::Integer)
  return Boundary(elnum, bndry.face)
end

#TODO: make idx and sz fields of the type
#      make the push_* functions update the counts
#      get rid of shocked and neighbors arrays, use only elnums_all

"""
  After all the elements that have shocks in them have been added, finished
  constructing the mesh data structure.

  **Inputs**

   * mesh
   * shockmesh
"""
function completeShockElements(mesh::AbstractMesh, data::ShockedElements)

  # Use mesh.interfaces to get both data.interfaces and the neighbor elements

  # The purpose of this is to define a reduced element and interface numbering
  # scheme that only numbers the entities needed to compute the shock capturing
  # terms.
  # The element numbers in data.elnums_all are assigned 1:n for the
  # elements with positive shock indicator, and (n+1):m for the neighboring
  # elements

  data.numShock = data.idx_all - 1
  # neighbor stuff
  # iface stuff
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    elementL = data.elnums_mesh[iface_i.elementL]
    elementR = data.elnums_mesh[iface_i.elementR]

    # if an element has not been seen before and its neighbor has a shock in
    # it, add to the list of neighbor elements
    if (elementL > 0 && elementL <= data.numShock) ||
       (elementR > 0 && elementR <= data.numShock)

      if elementL == 0
        elnum_full = iface_i.elementL  # element numbering in the full numbering
        elementL = data.idx_all  # element number in the reduced numbering

        push_neighbor(data, elnum_full)

      elseif elementR == 0  # it shouldn't be possible to see the same pair of
                          # elements more than once (even on periodic meshes)
        elnum_full = iface_i.elementR
        elementR = data.idx_all

        push_neighbor(data, elnum_full)
      end

    #if (elementL > 0 && elementL <= data.numShock) &&
    #   (elementR > 0 && elementR <= data.numShock)
      # record the new interface
      iface_new = RelativeInterface(replace_interface(iface_i, elementL, elementR), i)
      push_iface(data, iface_new)
    end  # end if
  end  # end for

  data.numNeighbor = data.idx_all - 1 - data.numShock
  data.numInterfaces = data.idx_if - 1

  # get shared interfaces
  setupShockmeshParallel(mesh, data)

#=
 # get the list of boundary faces
  for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    elnum = data.elnums_mesh[bndry_i.element]
    if elnum > 0 && elnum <= data.numShock
      bndry_new = RelativeBoundary(replace_boundary(bndry_i, elnum), i)
      push_bndry(data, bndry_new)
    end
  end
  data.numBoundaryFaces = data.idx_b - 1
=#
  # get boundary faces and Dirichlet info
  fill!(data.bndry_offsets, 0); data.bndry_offsets[1] = 1
  for i=1:mesh.numBC
    bc_range = mesh.bndry_offsets[i]:(mesh.bndry_offsets[i+1]-1)
    nfaces = 0
    for j in bc_range
      bndry_j = mesh.bndryfaces[j]
      elnum = data.elnums_mesh[bndry_j.element]
      # check if the element is in the data
      if elnum > 0 && elnum <= data.numShock
        bndry_new = RelativeBoundary(replace_boundary(bndry_j, elnum), j)
        push_bndry(data, bndry_new)
        nfaces += 1
      end
    end  # end j
    data.bndry_offsets[i+1] = data.bndry_offsets[i] + nfaces
  end  # end i
#=
  # add all boundaries that are not on the boundary of the original domain.
  # Don't include them in data.bndry_offsets
  @assert mesh.coord_order == 1
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    elementL = data.elnums_mesh[iface_i.elementL]
    elementR = data.elnums_mesh[iface_i.elementR]

    # if an element has not been seen before and its neighbor has a shock in
    # it, add to the list of neighbor elements
    elL_nonzero = elementL > 0 && elementL <= data.numShock
    elR_nonzero = elementR > 0 && elementR <= data.numShock

    # only one element non zero
    if elL_nonzero || elR_nonzero && !(elL_nonzero && elR_nonzero)

      if elementL == 0
        # add face of elementR as boundary
        bndry_new = RelativeBoundary(Boundary(elementR, iface_i.faceR), i, -1)
        push_bndry(data, bndry_new)

      elseif elementR == 0
        # add face of elementL as boundary
        #TODO: record orientation info needed for normal vector
        bndry_new = RelativeBoundary(Boundary(elementL, iface_i.faceL), i)
        push_bndry(data, bndry_new)
      end
    end
  end
=#
  data.numBoundaryFaces = data.idx_b - 1

  # compute final counts that will be useful for the user
  numEl = data.numShock + data.numNeighbor
  for i=1:data.npeers
    numEl += data.numShared[i]
  end
  data.numEl = numEl

  # setup ranges
  data.local_els = 1:data.numShock
  data.neighbor_els = (data.numShock+1):(data.numShock + data.numNeighbor)
  startidx = data.numShock + data.numNeighbor + 1
  data.shared_els = Vector{UnitRange{Int}}(data.npeers)
  for i=1:data.npeers
    data.shared_els[i] = startidx:(startidx + data.numShared[i] - 1)
    startidx = data.shared_els[i][end] + 1
  end

  return nothing
end

function setupShockmeshParallel(mesh::AbstractMesh, data::ShockedElements)

  # Its possible to re-use the arrays from last iteration, but it would be
  # tricky when the set of peer processes changes from one iteration to the next
  data.peer_indices = Array{Int}(0)
  data.shared_interfaces = Vector{Vector{RelativeInterface}}(0)
  data.numSharedInterfaces = Array{Int}(0)
  data.numShared = Array{Int}(0)
  data.idx_sf = Array{Int}(0)
  data.sz_sf = Array{Int}(0)
  data.npeers = 0
  const INITIAL_SIZE = 8
  found_peer = false
  for peer=1:mesh.npeers
    for i=1:length(mesh.shared_interfaces[peer])
      iface_i = mesh.shared_interfaces[peer][i]

      # elementL is always the local one, so only need to check that one
      elementL = data.elnums_mesh[iface_i.elementL]
      elementR = data.elnums_mesh[iface_i.elementR]

      if (elementL > 0) && (elementL <= data.numShock)

        # add the peer
        # need the found_peer flag because peer_indices[0] is an error, so
        # don't evaluate unless length(peer_indices) > 0
        if !found_peer || data.peer_indices[end] != peer
          found_peer = true
          data.npeers += 1
          push!(data.peer_indices, peer)
          # make the new vector non-zero size to save re-allocating it on the
          # first push
          push!(data.shared_interfaces, Vector{RelativeInterface}(INITIAL_SIZE))
          push!(data.numSharedInterfaces, 0)
          push!(data.numShared, 0)
          push!(data.idx_sf, 1)
          push!(data.sz_sf, INITIAL_SIZE)
        end

        # add elementR to the list of elements
        if elementR == 0
          elnum_full = iface_i.elementR
          elementR = data.idx_all

          push_shared(data, elnum_full)
        end

        # add the interface
        iface_new = RelativeInterface(replace_interface(iface_i, elementL, elementR), i)
        push_sharediface(data, iface_new)

      end  # end if elementL
    end  # end i
  end  # end peer

  return nothing
end
#TODO: have push_* functions return the new element number

#TODO: make one of these for regular meshes
"""
  Helper function to get the index of a shared element in the `RemoteMetrics`
  arrays.

  **Inputs**

   * data: a `ShockedElements`
   * mesh: the original mesh
   * peeridx: the peer index on the `ShockedElements`
   * elnum: the element number on the `ShockedElements`

  **Outputs**

   * idx: the index of the shared element in, for example, `RemoteMetrics`
          arrays
"""
function getSharedElementIndex(data::ShockedElements, mesh::AbstractMesh,
                               peeridx::Integer, elnum::Integer)

  peeridx2 = data.peer_indices[peeridx]
  firstnum = mesh.shared_element_offsets[peeridx2]
  elnum2 = data.elnums_all[elnum]
  return elnum2 - firstnum + 1
end


text:
module ZenFab

#using Blink
using Compat
using Reactive
using Mux
using JSON
using Markdown
using Gadfly
using Escher
using Patchwork
using Hiccup

include(Pkg.dir("Escher", "src", "cli", "serve.jl"))

end # module


text:
using DataStructures
using Statistics
using JLD

raw_data_dir  = "/home/user/Downloads/raw_data";
data_dir = "/home/user/thesis/data/scaling";
t_0_dir   = joinpath(raw_data_dir, "cluster_size_distribution/t_0");
t_1_dir   = joinpath(raw_data_dir, "cluster_size_distribution/t_1");
t_0_files = joinpath.(t_0_dir, readdir(t_0_dir));
t_1_files = joinpath.(t_1_dir, readdir(t_1_dir));

# %%

function μ(x)
	return sum(x) / length(x)
end

function σ(x)
	return sqrt(sum((x .- μ(x)).^2 / length(x)))
end

file_dict = Dict(
	"t_0" => joinpath.(t_0_dir, readdir(t_0_dir)),
	"t_1" => joinpath.(t_1_dir, readdir(t_1_dir))
);

function compute_distributions(file_dict)
	t_dicts = Dict()
	for (t, files) in file_dict
		t_dicts[t] = Dict()

		for file in files
			for (i, line) in enumerate(eachline(file))
				i ∈ (1, 2, 3) ? continue : nothing

				s, n_s = parse.(Int, split(line, ","))
				if s ∈ keys(t_dicts[t])
					push!(t_dicts[t][s], n_s)
				else
					t_dicts[t][s] = [n_s]
				end
			end
		end
	end

	t_dicts_out = Dict()
	for (t, distribution_arrays) in t_dicts
		s      = []
		counts = []
		n_s    = []
		σ_n_s  = []
		for (s_, distribution_array) in distribution_arrays
			push!(s, s_)
			push!(counts, length(n_s))
			push!(n_s, μ(distribution_array))
			push!(σ_n_s, std(distribution_array))
		end
		t_dicts_out[t] = Dict(
			"s"     => s,
			"counts" => counts,
			"n_s"   => n_s,
			"σ_n_s" => σ_n_s
		)
	end

	return t_dicts_out
end

t_dicts = compute_distributions(file_dict);
t_dicts["t_0"]
t_dicts["t_1"]

# %% Save and load data

save_file = joinpath(data_dir, "cluster_size_distribution_data.jld");
save(save_file, "distribution_avgs", t_dicts)

using PyCall
np = pyimport("numpy")

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
@testset "_mult!" begin
    a = collect(0:5)
    b = collect(1:6)
    _mult!(a, b, 6)
    @test a == 0:7:35
end

@testset "findcell" begin
    hrs = exampledata("hrs")
    @test_throws ArgumentError findcell((), hrs)
    @test_throws ArgumentError findcell((:wave,), hrs, falses(size(hrs, 1)))

    rows = findcell((:wave,), hrs)
    @test length(rows) == 5
    @test rows[UInt32(5)] == findall(x->x==7, hrs.wave)

    df = DataFrame(hrs)
    df.wave = PooledArray(df.wave)
    @test findcell((:wave,), df) == rows

    esample = hrs.wave.!=11
    rows = findcell((:wave, :wave_hosp), hrs, esample)
    @test length(rows) == 16
    @test rows[one(UInt32)] == intersect(findall(x->x==10, view(hrs.wave, esample)), findall(x->x==10, view(hrs.wave_hosp, esample)))

    rows = findcell((:wave, :wave_hosp, :male), hrs)
    @test length(rows) == 40
end

@testset "cellrows" begin
    hrs = exampledata("hrs")
    cols0 = subcolumns(hrs, (:wave, :wave_hosp), falses(size(hrs, 1)))
    cols = subcolumns(hrs, (:wave, :wave_hosp))
    rows_dict0 = IdDict{UInt32, Vector{Int}}()
    @test_throws ArgumentError cellrows(cols, rows_dict0)
    rows_dict = findcell(cols)
    @test_throws ArgumentError cellrows(cols0, rows_dict)

    cells, rows = cellrows(cols, rows_dict)
    @test length(cells[1]) == length(rows) == 20
    @test Tables.matrix(cells) ==
        sortslices(unique(hcat(hrs.wave, hrs.wave_hosp), dims=1), dims=1)
    @test propertynames(cells) == [:wave, :wave_hosp]
    @test rows[1] == intersect(findall(x->x==7, hrs.wave), findall(x->x==8, hrs.wave_hosp))

    df = DataFrame(hrs)
    df.wave = ScaledArray(hrs.wave, 1)
    df.wave_hosp = ScaledArray(hrs.wave_hosp, 1)
    cols1 = subcolumns(df, (:wave, :wave_hosp))
    cells1, rows1 = cellrows(cols1, rows_dict)
    @test rows1 == rows
    @test cells1 == cells
    @test cells1.wave isa ScaledArray
    @test cells1.wave_hosp isa ScaledArray

    rot = ones(size(df, 1))
    df.wave = RotatingTimeArray(rot, hrs.wave)
    df.wave_hosp = RotatingTimeArray(rot, hrs.wave_hosp)
    cols2 = subcolumns(df, (:wave, :wave_hosp))
    cells2, rows2 = cellrows(cols2, rows_dict)
    @test rows2 == rows
    @test cells2 == cells
    @test cells2.wave isa RotatingTimeArray
    @test cells2.wave_hosp isa RotatingTimeArray
end

@testset "settime aligntime" begin
    hrs = exampledata("hrs")
    t = settime(hrs.wave, 2, reftype=Int16)
    @test eltype(refarray(t)) == Int16
    @test sort!(unique(refarray(t))) == 1:3
    t = settime(hrs.wave, 0.5)
    @test eltype(t) == Float64

    rot = isodd.(hrs.wave)
    t = settime(hrs.wave, rotation=rot)
    @test eltype(t) == RotatingTimeValue{Bool, eltype(hrs.wave)}
    @test refpool(t) === nothing
    @test eltype(refarray(t.time)) == Int32
    @test t.rotation == rot
    @test t.time == hrs.wave

    df = DataFrame(hrs)
    df.wave1 = Date.(df.wave)
    t = settime(df.wave1, Year(1))
    @test t == df.wave1
    @test t.pool == Date(7):Year(1):Date(11)

    t = settime(df.wave1, Year(1), rotation=rot)
    @test eltype(t) == RotatingTimeValue{Bool, eltype(df.wave1)}
    @test eltype(refarray(t.time)) == Int32
    @test t.time isa ScaledArray

    df.t1 = settime(df.wave, 1)
    df.t2 = settime(df.wave, 2)
    t = aligntime(df, :t2, :t1)
    @test t == df.t2
    @test t.pool == df.t1.pool
    @test t.invpool == df.t1.invpool
    df.t2 = settime(df.wave, start=0, stop=20)
    t = aligntime(df, :t2, :t1)
    @test t == df.t2
    @test t.pool == df.t1.pool

    t1 = settime(df.wave1, Year(1), rotation=rot)
    t2 = settime(df.wave1, Year(1), start=Date(5))
    t = aligntime(t2, t1)
    @test t == t1
    @test t.time.pool == t1.time.pool
    t2 = settime(df.wave1, Year(1), start=Date(5), rotation=rot)
    t = aligntime(t2, t1)
    @test t == t1
    @test t.time.pool == t1.time.pool
end

@testset "PanelStructure" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    panel = setpanel(hrs, :hhidpn, :wave)
    @test length(unique(panel.refs)) == N
    inidbounds = view(hrs.hhidpn, 2:N) .== view(hrs.hhidpn, 1:N-1)
    @test view(diff(panel.refs), inidbounds) == view(diff(hrs.wave), inidbounds)
    @test length(panel.idpool) == 656
    @test panel.timepool == 7:11

    panel1 = setpanel(hrs, :hhidpn, :wave, step=0.5, reftype=Int)
    @test view(diff(panel1.refs), inidbounds) == 2 .* view(diff(hrs.wave), inidbounds)
    @test panel1.timepool == 7.0:0.5:11.0
    @test eltype(panel1.refs) == Int

    @test_throws ArgumentError setpanel(hrs, :hhidpn, :oop_spend)
    @test_throws DimensionMismatch setpanel(hrs.hhidpn, 1:100)

    @test sprint(show, panel) == "Panel Structure"
    t = VERSION < v"1.6.0" ? "Array{Int64,1}" : " Vector{Int64}"
    @test sprint(show, MIME("text/plain"), panel) == """
        Panel Structure:
          idpool:   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  647, 648, 649, 650, 651, 652, 653, 654, 655, 656]
          timepool:   7:1:11
          laginds:  Dict{Int64,$t}()"""

    lags = findlag!(panel)
    @test sprint(show, MIME("text/plain"), panel) == """
        Panel Structure:
          idpool:   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  647, 648, 649, 650, 651, 652, 653, 654, 655, 656]
          timepool:   7:1:11
          laginds:  Dict(1 => [4, 5, 1, 2, 0, 7, 0, 10, 8, 6  …  0, 3275, 3271, 3273, 3274, 3279, 3280, 3277, 3278, 0])"""

    leads = findlead!(panel, -1)
    @test leads == leads
    @test_throws ArgumentError findlag!(panel, 5)

    @test ilag!(panel) === panel.laginds[1]
    @test ilag!(panel, 2) === panel.laginds[2]
    @test ilead!(panel) === panel.laginds[-1]

    l1 = lag(panel, hrs.wave)
    @test view(l1, hrs.wave.!=7) == view(hrs.wave, hrs.wave.!=7) .- 1
    @test all(ismissing, view(l1, hrs.wave.==7))
    l2 = lag(panel, hrs.wave, 2, default=-1)
    @test view(l2, hrs.wave.>=9) == view(hrs.wave, hrs.wave.>=9) .- 2
    @test all(x->x==-1, view(l2, hrs.wave.<9))
    @test eltype(l2) == Int
    l_2 = lead(panel, hrs.wave, 2, default=-1.0)
    @test view(l_2, hrs.wave.<=9) == view(hrs.wave, hrs.wave.<=9) .+ 2
    @test all(x->x==-1, view(l_2, hrs.wave.>9))
    @test eltype(l_2) == Int

    @test_throws DimensionMismatch lag(panel, 1:10)

    d1 = diff(panel, hrs.wave)
    @test all(x->x==1, view(d1, hrs.wave.!=7))
    @test all(ismissing, view(d1, hrs.wave.==7))
    d2 = diff(panel, hrs.wave, l=2, default=-1)
    @test all(x->x==2, view(d2, hrs.wave.>=9))
    @test all(x->x==-1, view(d2, hrs.wave.<9))
    @test eltype(d2) == Int
    d_2 = diff(panel, hrs.wave, l=-2, default=-1.0)
    @test all(x->x==-2, view(d_2, hrs.wave.<=9))
    @test all(x->x==-1, view(d_2, hrs.wave.>9))
    @test eltype(d_2) == Int

    d2 = diff(panel, hrs.wave, order=2)
    @test all(x->x==0, view(d2, hrs.wave.>=9))
    @test all(ismissing, view(d2, hrs.wave.<9))
    d2 = diff(panel, hrs.wave, order=2, l=2)
    @test all(x->x==0, view(d2, hrs.wave.==11))
    @test all(ismissing, view(d2, hrs.wave.!=11))

    @test_throws ArgumentError diff(panel, hrs.wave, l=5)
    @test_throws ArgumentError diff(panel, hrs.wave, order=5)
end


text:
function ArrayRealVector()
    return ArrayRealVector(())
end

function ArrayRealVector(arg0::ArrayRealVector)
    return ArrayRealVector((ArrayRealVector,), arg0)
end

function ArrayRealVector(arg0::ArrayRealVector, arg1::ArrayRealVector)
    return ArrayRealVector((ArrayRealVector, ArrayRealVector), arg0, arg1)
end

function ArrayRealVector(arg0::ArrayRealVector, arg1::RealVector)
    return ArrayRealVector((ArrayRealVector, RealVector), arg0, arg1)
end

function ArrayRealVector(arg0::ArrayRealVector, arg1::Vector{jdouble})
    return ArrayRealVector((ArrayRealVector, Vector{jdouble}), arg0, arg1)
end

function ArrayRealVector(arg0::ArrayRealVector, arg1::jboolean)
    return ArrayRealVector((ArrayRealVector, jboolean), arg0, arg1)
end

function ArrayRealVector(arg0::RealVector)
    return ArrayRealVector((RealVector,), arg0)
end

function ArrayRealVector(arg0::RealVector, arg1::ArrayRealVector)
    return ArrayRealVector((RealVector, ArrayRealVector), arg0, arg1)
end

function ArrayRealVector(arg0::Vector{Double})
    return ArrayRealVector((Vector{Double},), arg0)
end

function ArrayRealVector(arg0::Vector{Double}, arg1::jint, arg2::jint)
    return ArrayRealVector((Vector{Double}, jint, jint), arg0, arg1, arg2)
end

function ArrayRealVector(arg0::Vector{jdouble})
    return ArrayRealVector((Vector{jdouble},), arg0)
end

function ArrayRealVector(arg0::Vector{jdouble}, arg1::ArrayRealVector)
    return ArrayRealVector((Vector{jdouble}, ArrayRealVector), arg0, arg1)
end

function ArrayRealVector(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return ArrayRealVector((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function ArrayRealVector(arg0::Vector{jdouble}, arg1::jboolean)
    return ArrayRealVector((Vector{jdouble}, jboolean), arg0, arg1)
end

function ArrayRealVector(arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return ArrayRealVector((Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function ArrayRealVector(arg0::jint)
    return ArrayRealVector((jint,), arg0)
end

function ArrayRealVector(arg0::jint, arg1::jdouble)
    return ArrayRealVector((jint, jdouble), arg0, arg1)
end

function add(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "add", ArrayRealVector, (RealVector,), arg0)
end

function add_to_entry(obj::ArrayRealVector, arg0::jint, arg1::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jdouble), arg0, arg1)
end

function append(obj::ArrayRealVector, arg0::ArrayRealVector)
    return jcall(obj, "append", ArrayRealVector, (ArrayRealVector,), arg0)
end

function append(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "append", RealVector, (RealVector,), arg0)
end

function append(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "append", RealVector, (jdouble,), arg0)
end

function combine(obj::ArrayRealVector, arg0::jdouble, arg1::jdouble, arg2::RealVector)
    return jcall(obj, "combine", ArrayRealVector, (jdouble, jdouble, RealVector), arg0, arg1, arg2)
end

function combine_to_self(obj::ArrayRealVector, arg0::jdouble, arg1::jdouble, arg2::RealVector)
    return jcall(obj, "combineToSelf", ArrayRealVector, (jdouble, jdouble, RealVector), arg0, arg1, arg2)
end

function copy(obj::ArrayRealVector)
    return jcall(obj, "copy", ArrayRealVector, ())
end

function cosine(obj::RealVector, arg0::RealVector)
    return jcall(obj, "cosine", jdouble, (RealVector,), arg0)
end

function dot_product(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "dotProduct", jdouble, (RealVector,), arg0)
end

function ebe_divide(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "ebeDivide", ArrayRealVector, (RealVector,), arg0)
end

function ebe_multiply(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "ebeMultiply", ArrayRealVector, (RealVector,), arg0)
end

function equals(obj::ArrayRealVector, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_data_ref(obj::ArrayRealVector)
    return jcall(obj, "getDataRef", Vector{jdouble}, ())
end

function get_dimension(obj::ArrayRealVector)
    return jcall(obj, "getDimension", jint, ())
end

function get_distance(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "getDistance", jdouble, (RealVector,), arg0)
end

function get_entry(obj::ArrayRealVector, arg0::jint)
    return jcall(obj, "getEntry", jdouble, (jint,), arg0)
end

function get_l1_distance(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "getL1Distance", jdouble, (RealVector,), arg0)
end

function get_l1_norm(obj::ArrayRealVector)
    return jcall(obj, "getL1Norm", jdouble, ())
end

function get_l_inf_distance(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "getLInfDistance", jdouble, (RealVector,), arg0)
end

function get_l_inf_norm(obj::ArrayRealVector)
    return jcall(obj, "getLInfNorm", jdouble, ())
end

function get_max_index(obj::RealVector)
    return jcall(obj, "getMaxIndex", jint, ())
end

function get_max_value(obj::RealVector)
    return jcall(obj, "getMaxValue", jdouble, ())
end

function get_min_index(obj::RealVector)
    return jcall(obj, "getMinIndex", jint, ())
end

function get_min_value(obj::RealVector)
    return jcall(obj, "getMinValue", jdouble, ())
end

function get_norm(obj::ArrayRealVector)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_sub_vector(obj::ArrayRealVector, arg0::jint, arg1::jint)
    return jcall(obj, "getSubVector", RealVector, (jint, jint), arg0, arg1)
end

function hash_code(obj::ArrayRealVector)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::ArrayRealVector)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::ArrayRealVector)
    return jcall(obj, "isNaN", jboolean, ())
end

function iterator(obj::RealVector)
    return jcall(obj, "iterator", JIterator, ())
end

function map(obj::ArrayRealVector, arg0::UnivariateFunction)
    return jcall(obj, "map", ArrayRealVector, (UnivariateFunction,), arg0)
end

function map_add(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapAdd", RealVector, (jdouble,), arg0)
end

function map_add_to_self(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "mapAddToSelf", RealVector, (jdouble,), arg0)
end

function map_divide(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapDivide", RealVector, (jdouble,), arg0)
end

function map_divide_to_self(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "mapDivideToSelf", RealVector, (jdouble,), arg0)
end

function map_multiply(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapMultiply", RealVector, (jdouble,), arg0)
end

function map_multiply_to_self(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "mapMultiplyToSelf", RealVector, (jdouble,), arg0)
end

function map_subtract(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapSubtract", RealVector, (jdouble,), arg0)
end

function map_subtract_to_self(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "mapSubtractToSelf", RealVector, (jdouble,), arg0)
end

function map_to_self(obj::ArrayRealVector, arg0::UnivariateFunction)
    return jcall(obj, "mapToSelf", ArrayRealVector, (UnivariateFunction,), arg0)
end

function outer_product(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "outerProduct", RealMatrix, (RealVector,), arg0)
end

function projection(obj::RealVector, arg0::RealVector)
    return jcall(obj, "projection", RealVector, (RealVector,), arg0)
end

function set(obj::ArrayRealVector, arg0::jdouble)
    return jcall(obj, "set", void, (jdouble,), arg0)
end

function set_entry(obj::ArrayRealVector, arg0::jint, arg1::jdouble)
    return jcall(obj, "setEntry", void, (jint, jdouble), arg0, arg1)
end

function set_sub_vector(obj::ArrayRealVector, arg0::jint, arg1::RealVector)
    return jcall(obj, "setSubVector", void, (jint, RealVector), arg0, arg1)
end

function set_sub_vector(obj::ArrayRealVector, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setSubVector", void, (jint, Vector{jdouble}), arg0, arg1)
end

function sparse_iterator(obj::RealVector)
    return jcall(obj, "sparseIterator", JIterator, ())
end

function subtract(obj::ArrayRealVector, arg0::RealVector)
    return jcall(obj, "subtract", ArrayRealVector, (RealVector,), arg0)
end

function to_array(obj::ArrayRealVector)
    return jcall(obj, "toArray", Vector{jdouble}, ())
end

function to_string(obj::ArrayRealVector)
    return jcall(obj, "toString", JString, ())
end

function unit_vector(obj::RealVector)
    return jcall(obj, "unitVector", RealVector, ())
end

function unitize(obj::RealVector)
    return jcall(obj, "unitize", void, ())
end

function unmodifiable_real_vector(::Type{RealVector}, arg0::RealVector)
    return jcall(RealVector, "unmodifiableRealVector", RealVector, (RealVector,), arg0)
end

function walk_in_default_order(obj::ArrayRealVector, arg0::RealVectorChangingVisitor)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorChangingVisitor,), arg0)
end

function walk_in_default_order(obj::ArrayRealVector, arg0::RealVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_default_order(obj::ArrayRealVector, arg0::RealVectorPreservingVisitor)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorPreservingVisitor,), arg0)
end

function walk_in_default_order(obj::ArrayRealVector, arg0::RealVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::ArrayRealVector, arg0::RealVectorChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::ArrayRealVector, arg0::RealVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::ArrayRealVector, arg0::RealVectorPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::ArrayRealVector, arg0::RealVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end


text:
import Base.Random: rand, rand!

abstract VSLContinuousDistribution <: VSLDistribution

macro vsl_distribution_continuous(name, methods, tmp, properties...)
    t2 = :(Union{})
    method_symbols = [Symbol(string("VSL_RNG_METHOD_", method)) for method in methods.args]
    for method in method_symbols
        push!(t2.args, method)
    end
    code = quote
        type $name{T1<:Union{Cfloat, Cdouble}, T2<:$t2} <: VSLContinuousDistribution
            brng::BasicRandomNumberGenerator
            ii::Int
        end
    end
    fields = code.args[2].args[end].args
    push!(fields, tmp)
    for property in properties
        push!(fields, property)
    end
    push!(fields, :(method::Type{T2}))

    push!(code.args, :(
        function $name{T1<:Union{Cfloat, Cdouble}, T2<:$t2}(brng::BasicRandomNumberGenerator, $(properties...), method::Type{T2})
            tmp = $(tmp.args[2].args[1] == :Vector ? :(Vector{T1}(BUFFER_LENGTH)) : :(Matrix{T1}(dimen, BUFFER_LENGTH)))
            $name(brng, BUFFER_LENGTH, tmp, $([property.args[1] for property in properties]...), method)
        end
    ))

    push!(code.args, :(
        $name{T1<:Union{Cfloat, Cdouble}}(brng::BasicRandomNumberGenerator, $(properties...)) =
        $name(brng, $([property.args[1] for property in properties]...), $(method_symbols[1]))
    ))

    esc(code)
end

macro register_rand_functions_continuous(name, methods, method_constants, types, arguments)
    code = quote end
    method_symbols = [Symbol(string("VSL_RNG_METHOD_", method)) for method in methods.args]
    for ttype in (:Cfloat, :Cdouble)
        ctype = ttype == :Cfloat ? 's' : 'd'
        argtypes = [atype == :T1 ? ttype : atype for atype in types.args]
        function_name = "v$(ctype)Rng$(name)"
        for (method, constant) in zip(method_symbols, method_constants.args)
            push!(
                code.args,
                :(function rand(d::$name{$ttype, $method}, ::Type{$ttype}=$ttype)
                    d.ii += 1
                    if d.ii > BUFFER_LENGTH
                        ccall(($function_name, libmkl), Cint, (Int, Ptr{Void}, Int, Ptr{$ttype}, $(argtypes...)),
                            $constant, d.brng.stream_state[1], BUFFER_LENGTH, d.tmp, $(arguments.args...))
                        d.ii = 1
                    end
                    d.tmp[d.ii]
                end)
            )
            push!(
                code.args,
                :(function rand!(d::$name{$ttype, $method}, A::Array{$ttype})
                    n = length(A)
                    t = BUFFER_LENGTH - d.ii
                    if n <= t
                        copy!(A, 1, d.tmp, d.ii + 1, n)
                        d.ii = d.ii + n
                        return A
                    end
                    copy!(A, 1, d.tmp, d.ii + 1, t)
                    d.ii = BUFFER_LENGTH
                    ccall(($function_name, libmkl), Cint, (Int, Ptr{Void}, Int, Ptr{$ttype}, $(argtypes...)),
                        $constant, d.brng.stream_state[1], n - t, view(A, t+1:n), $(arguments.args...))
                    A
                end)
            )
            push!(
                code.args,
                :(rand(d::$name{$ttype, $method}, dims::Dims) = rand!(d, Array($ttype, dims)))
            )
        end
    end
    esc(code)
end

macro register_rand_functions_continuous_multivariate(name, methods, method_constants, types, arguments)
    code = quote end
    method_symbols = [Symbol(string("VSL_RNG_METHOD_", method)) for method in methods.args]
    for ttype in (:Cfloat, :Cdouble)
        ctype = ttype == :Cfloat ? 's' : 'd'
        argtypes = [parse(replace("$atype", "T1", "$ttype")) for atype in types.args]
        function_name = "v$(ctype)Rng$(name)"
        for (method, constant) in zip(method_symbols, method_constants.args)
            push!(
                code.args,
                :(function rand(d::$name{$ttype, $method}, ::Type{Vector{$ttype}}=Vector{$ttype})
                    d.ii += 1
                    if d.ii > BUFFER_LENGTH
                        ccall(($function_name, libmkl), Cint, (Int, Ptr{Void}, Int, Ptr{$ttype}, Int, $(argtypes...)),
                            $constant, d.brng.stream_state[1], BUFFER_LENGTH, d.tmp, d.dimen, $(arguments.args...))
                        d.ii = 1
                    end
                    d.tmp[:, d.ii]
                end)
            )
        end
    end
    esc(code)
end

@vsl_distribution_continuous(
    Uniform,
    [STD, STD_ACCURATE],
    tmp::Vector{T1},
    a::T1,
    b::T1
)

@register_rand_functions_continuous(
    Uniform,
    [STD, STD_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1],
    [d.a, d.b]
)

@vsl_distribution_continuous(
    Gaussian,
    [BOXMULLER, BOXMULLER2, ICDF],
    tmp::Vector{T1},
    a::T1,
    σ::T1
)

@register_rand_functions_continuous(
    Gaussian,
    [BOXMULLER, BOXMULLER2, ICDF],
    [0x00000000, 0x00000001, 0x00000002],
    [T1, T1],
    [d.a, d.σ]
)

@enum(MatrixStorageType,
    VSL_MATRIX_STORAGE_FULL     = 0,
    VSL_MATRIX_STORAGE_PACKED   = 1,
    VSL_MATRIX_STORAGE_DIAGONAL = 2
)

@vsl_distribution_continuous(
    GaussianMV,
    [BOXMULLER, BOXMULLER2, ICDF],
    tmp::Matrix{T1},
    dimen::Int,
    mstorage::MatrixStorageType,
    a::Vector{T1},
    t::Matrix{T1}
)

@register_rand_functions_continuous_multivariate(
    GaussianMV,
    [BOXMULLER, BOXMULLER2, ICDF],
    [0x00000000, 0x00000001, 0x00000002],
    [Int, Ptr{T1}, Ptr{T1}],
    [Int(d.mstorage), d.a, d.t]
)

@vsl_distribution_continuous(
    Exponential,
    [ICDF, ICDF_ACCURATE],
    tmp::Vector{T1},
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Exponential,
    [ICDF, ICDF_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1],
    [d.a, d.β]
)

@vsl_distribution_continuous(
    Laplace,
    [ICDF],
    tmp::Vector{T1},
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Laplace,
    [ICDF],
    [0x00000000],
    [T1, T1],
    [d.a, d.β]
)

@vsl_distribution_continuous(
    Weibull,
    [ICDF, ICDF_ACCURATE],
    tmp::Vector{T1},
    α::T1,
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Weibull,
    [ICDF, ICDF_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1, T1],
    [d.α, d.a, d.β]
)

@vsl_distribution_continuous(
    Cauchy,
    [ICDF],
    tmp::Vector{T1},
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Cauchy,
    [ICDF],
    [0x00000000],
    [T1, T1],
    [d.a, d.β]
)

@vsl_distribution_continuous(
    Rayleigh,
    [ICDF, ICDF_ACCURATE],
    tmp::Vector{T1},
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Rayleigh,
    [ICDF, ICDF_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1],
    [d.a, d.β]
)

@vsl_distribution_continuous(
    Lognormal,
    [BOXMULLER2, BOXMULLER2_ACCURATE, ICDF, ICDF_ACCURATE],
    tmp::Vector{T1},
    a::T1,
    σ::T1,
    b::T1,
    β::T1
)

@register_rand_functions_continuous(
    Lognormal,
    [BOXMULLER2, BOXMULLER2_ACCURATE, ICDF, ICDF_ACCURATE],
    [0x00000000, 0x40000000, 0x00000001, 0x40000001],
    [T1, T1, T1, T1],
    [d.a, d.σ, d.b, d.β]
)

@vsl_distribution_continuous(
    Gumbel,
    [ICDF],
    tmp::Vector{T1},
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Gumbel,
    [ICDF],
    [0x00000000],
    [T1, T1],
    [d.a, d.β]
)

@vsl_distribution_continuous(
    Gamma,
    [GNORM, GNORM_ACCURATE],
    tmp::Vector{T1},
    α::T1,
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Gamma,
    [GNORM, GNORM_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1, T1],
    [d.α, d.a, d.β]
)

@vsl_distribution_continuous(
    Beta,
    [CJA, CJA_ACCURATE],
    tmp::Vector{T1},
    p::T1,
    q::T1,
    a::T1,
    β::T1
)

@register_rand_functions_continuous(
    Beta,
    [CJA, CJA_ACCURATE],
    [0x00000000, 0x40000000],
    [T1, T1, T1, T1],
    [d.p, d.q, d.a, d.β]
)


text:
module IO

using JSON2
using ClimateMARGO.Models

export export_state, import_state,
    export_parameters, import_parameters,
    included_configurations

include("json_io.jl")

end


text:
module TestTuning

# using Revise
using Test
using MLJ
using MLJBase
import Random.seed!
seed!(1234)

@load KNNRegressor

include("foobarmodel.jl")
include("simple_composite_model.jl")

x1 = rand(100);
x2 = rand(100);
x3 = rand(100);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

@testset "2-parameter tune, no nesting" begin

    sel = FeatureSelector()
    stand = UnivariateStandardizer()
    ridge = FooBarRegressor()
    composite = SimpleDeterministicCompositeModel(transformer=sel, model=ridge)

    features_ = range(composite, :(transformer.features), values=[[:x1], [:x1, :x2], [:x2, :x3], [:x1, :x2, :x3]])
    lambda_ = range(composite, :(model.lambda), lower=1e-6, upper=1e-1, scale=:log10)

    ranges = [features_, lambda_]

    holdout = Holdout(fraction_train=0.8)
    grid = Grid(resolution=10)

    tuned_model = TunedModel(model=composite, tuning=grid,
                             resampling=holdout, measure=rms,
                             ranges=ranges, full_report=false)

    MLJBase.info_dict(tuned_model)

    tuned = machine(tuned_model, X, y)

    fit!(tuned)
    r = report(tuned)
    @test r.best_report isa NamedTuple{(:machines, :reports)}
    tuned_model.full_report=true
    fit!(tuned)
    report(tuned)
    fp = fitted_params(tuned)
    @test fp.best_fitted_params isa NamedTuple{(:machines, :fitted_params)}
    b = fp.best_model
    @test b isa SimpleDeterministicCompositeModel

    measurements = tuned.report.measurements
    # should be all different:
    @test length(unique(measurements)) == length(measurements)

    @test length(b.transformer.features) == 3
    @test abs(b.model.lambda - 0.027825) < 1e-6

    # get the training error of the tuned_model:
    e = rms(y, predict(tuned, X))

    # check this error has same order of magnitude as best measurement
    # during tuning:
    r = e/tuned.report.best_measurement
    @test r < 10 && r > 0.1

    # test weights:
    tuned_model.weights = rand(length(y))
    fit!(tuned)
    @test tuned.report.measurements[1] != measurements[1]

end

#@testset "multiple resolutions" begin

    forest = EnsembleModel(atom=@load(KNNRegressor))
    pipe = @pipeline MyPipe(sel=FeatureSelector(), forest=forest)

    X = MLJ.table(rand(150, 2), names=[:x1, :x2])
    y = rand(150)

    r1 = range(pipe, :(sel.features), values=[[:x1,], [:x1, :x2]])
    r2 = range(pipe, :(forest.atom.K), lower = 1, upper=50)
    r3 = range(pipe, :(forest.bagging_fraction), lower=0.5, upper=1.0)
    ranges = [r1, r2, r3]

    holdout = Holdout(fraction_train=0.8)
    grid = Grid(resolution=5)

    tuned = TunedModel(model=pipe, tuning=grid,
                       resampling=holdout, measure=rms,
                       ranges=ranges)

    mach = machine(tuned, X, y)

    function num_grid_points(resolution)
        grid.resolution = resolution
        fit!(mach)
        return report(mach).measurements |> length
    end

    @test num_grid_points(6) == 72
    @test num_grid_points([:(forest.atom.K)=>7,
                           :(forest.bagging_fraction)=>4]) == 56
    @test num_grid_points([:(forest.atom.K)=>6, ]) == 60

    @testset "one parameter tune" begin
        ridge = FooBarRegressor()
        r = range(ridge, :lambda, lower=1e-7, upper=1e-3, scale=:log)
        tuned_model = TunedModel(model=ridge, ranges=r)
        tuned = machine(tuned_model, X, y)
        fit!(tuned)
        report(tuned)
    end

#end

@testset "nested parameter tune" begin
    tree_model = KNNRegressor()
    forest_model = EnsembleModel(atom=tree_model)
    r1 = range(forest_model, :(atom.K), lower=3, upper=4)
    r2 = range(forest_model, :bagging_fraction, lower=0.4, upper=1.0);
    self_tuning_forest_model = TunedModel(model=forest_model,
                                          tuning=Grid(resolution=3),
                                          resampling=Holdout(),
                                          ranges=[r1, r2],
                                          measure=rms)
    self_tuning_forest = machine(self_tuning_forest_model, X, y)
    fit!(self_tuning_forest, verbosity=2)
    r = report(self_tuning_forest)
    @test length(unique(r.measurements)) == 6
end

@testset "basic tuning with training weights" begin

    N = 100
    X = (x = rand(3N), );
    y = categorical(rand("abc", 3N));
    model = @load KNNClassifier
    r = range(model, :K, lower=2, upper=N)
    tuned_model = TunedModel(model=model,
                             measure=MLJ.BrierScore(),
                             resampling=Holdout(fraction_train=2/3),
                             range=r)

    # no weights:
    tuned = machine(tuned_model, X, y)
    fit!(tuned)
    best1 = fitted_params(tuned).best_model
    posterior1 = average([predict(tuned, X)...])

    # uniform weights:
    tuned = machine(tuned_model, X, y, fill(1, 3N))
    fit!(tuned)
    best2 = fitted_params(tuned).best_model
    posterior2 = average([predict(tuned, X)...])

    @test best1 == best2
    @test all([pdf(posterior1, c) ≈ pdf(posterior2, c) for c in levels(y)])

    # skewed weights:
    w = map(y) do η
        if η == 'a'
            return 2
        elseif η == 'b'
            return 4
        else
            return 1
        end
    end
    tuned = machine(tuned_model, X, y, w)
    fit!(tuned)
    best3 = fitted_params(tuned).best_model
    posterior3 = average([predict(tuned, X)...])

    # different tuning outcome:
    @test best1.K != best3.K

    # "posterior" is skewed appropriately in weighted case:
    @test abs(pdf(posterior3, 'b')/(2*pdf(posterior3, 'a'))  - 1) < 0.15
    @test abs(pdf(posterior3, 'b')/(4*pdf(posterior3, 'c'))  - 1) < 0.15

end

## LEARNING CURVE

@testset "learning curves" begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=50, rng=1)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda),
                     lower=0.0001, upper=0.1, scale=:log10)
    curve = MLJ.learning_curve!(mach; range=r_lambda)
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curve2 = MLJ.learning_curve!(mach; range=r_n)
    curve3 = learning_curve(ensemble, X, y; range=r_n)
    @test curve2.measurements ≈ curve3.measurements
end
end # module
true