text:
function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
                constant_speed::Val{false}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  for element in eachelement(dg, cache)
    max_λ1 = max_λ2 = max_λ3 = zero(max_scaled_speed)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)
      λ1, λ2, λ3 = max_abs_speeds(u_node, equations)
      max_λ1 = max(max_λ1, λ1)
      max_λ2 = max(max_λ2, λ2)
      max_λ3 = max(max_λ3, λ3)
    end
    inv_jacobian = cache.elements.inverse_jacobian[element]
    max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2 + max_λ3))
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end

function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
                constant_speed::Val{true}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  for element in eachelement(dg, cache)
    max_λ1, max_λ2, max_λ3 = max_abs_speeds(equations)
    inv_jacobian = cache.elements.inverse_jacobian[element]
    max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2 + max_λ3))
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end


text:
# Autogenerated wrapper script for Sundials_jll for x86_64-w64-mingw32-libgfortran4
export libsundials_arkode, libsundials_cvode, libsundials_cvodes, libsundials_ida, libsundials_idas, libsundials_kinsol, libsundials_nvecmanyvector, libsundials_nvecserial, libsundials_sunlinsolband, libsundials_sunlinsoldense, libsundials_sunlinsolklu, libsundials_sunlinsollapackband, libsundials_sunlinsollapackdense, libsundials_sunlinsolpcg, libsundials_sunlinsolspbcgs, libsundials_sunlinsolspfgmr, libsundials_sunlinsolspgmr, libsundials_sunlinsolsptfqmr, libsundials_sunmatrixband, libsundials_sunmatrixdense, libsundials_sunmatrixsparse, libsundials_sunnonlinsolfixedpoint, libsundials_sunnonlinsolnewton

using CompilerSupportLibraries_jll
using OpenBLAS_jll
using SuiteSparse_jll
JLLWrappers.@generate_wrapper_header("Sundials")
JLLWrappers.@declare_library_product(libsundials_arkode, "libsundials_arkode.dll")
JLLWrappers.@declare_library_product(libsundials_cvode, "libsundials_cvode.dll")
JLLWrappers.@declare_library_product(libsundials_cvodes, "libsundials_cvodes.dll")
JLLWrappers.@declare_library_product(libsundials_ida, "libsundials_ida.dll")
JLLWrappers.@declare_library_product(libsundials_idas, "libsundials_idas.dll")
JLLWrappers.@declare_library_product(libsundials_kinsol, "libsundials_kinsol.dll")
JLLWrappers.@declare_library_product(libsundials_nvecmanyvector, "libsundials_nvecmanyvector.dll")
JLLWrappers.@declare_library_product(libsundials_nvecserial, "libsundials_nvecserial.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolband, "libsundials_sunlinsolband.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsoldense, "libsundials_sunlinsoldense.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolklu, "libsundials_sunlinsolklu.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackband, "libsundials_sunlinsollapackband.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackdense, "libsundials_sunlinsollapackdense.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolpcg, "libsundials_sunlinsolpcg.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspbcgs, "libsundials_sunlinsolspbcgs.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspfgmr, "libsundials_sunlinsolspfgmr.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspgmr, "libsundials_sunlinsolspgmr.dll")
JLLWrappers.@declare_library_product(libsundials_sunlinsolsptfqmr, "libsundials_sunlinsolsptfqmr.dll")
JLLWrappers.@declare_library_product(libsundials_sunmatrixband, "libsundials_sunmatrixband.dll")
JLLWrappers.@declare_library_product(libsundials_sunmatrixdense, "libsundials_sunmatrixdense.dll")
JLLWrappers.@declare_library_product(libsundials_sunmatrixsparse, "libsundials_sunmatrixsparse.dll")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolfixedpoint, "libsundials_sunnonlinsolfixedpoint.dll")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolnewton, "libsundials_sunnonlinsolnewton.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS_jll, SuiteSparse_jll)
    JLLWrappers.@init_library_product(
        libsundials_arkode,
        "bin\\libsundials_arkode.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvode,
        "bin\\libsundials_cvode.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvodes,
        "bin\\libsundials_cvodes.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_ida,
        "bin\\libsundials_ida.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_idas,
        "bin\\libsundials_idas.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_kinsol,
        "bin\\libsundials_kinsol.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecmanyvector,
        "bin\\libsundials_nvecmanyvector.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecserial,
        "bin\\libsundials_nvecserial.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolband,
        "bin\\libsundials_sunlinsolband.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsoldense,
        "bin\\libsundials_sunlinsoldense.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolklu,
        "bin\\libsundials_sunlinsolklu.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackband,
        "bin\\libsundials_sunlinsollapackband.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackdense,
        "bin\\libsundials_sunlinsollapackdense.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolpcg,
        "bin\\libsundials_sunlinsolpcg.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspbcgs,
        "bin\\libsundials_sunlinsolspbcgs.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspfgmr,
        "bin\\libsundials_sunlinsolspfgmr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspgmr,
        "bin\\libsundials_sunlinsolspgmr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolsptfqmr,
        "bin\\libsundials_sunlinsolsptfqmr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixband,
        "bin\\libsundials_sunmatrixband.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixdense,
        "bin\\libsundials_sunmatrixdense.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixsparse,
        "bin\\libsundials_sunmatrixsparse.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolfixedpoint,
        "bin\\libsundials_sunnonlinsolfixedpoint.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolnewton,
        "bin\\libsundials_sunnonlinsolnewton.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# Matrix Multiplication Driver
#
# This is a driver program to test various ways of computing
# the product:
#                 C = A * B
#
# A and B are constant matrices, square and the order is
# set as a constant, ORDER (see definitions.py). This is so
# we can make a quick test of the multiplication result.
#
# History:   C++ version written by Tim Mattson, August 2010
#            Modified by Simon McIntosh-Smith, September 2011
#            Modified by Tom Deakin and Simon McIntosh-Smith, October 2012
#            Ported to Python by Tom Deakin, July 2013
#            Ported to Julia by Jake Bolewski, Nov 2013

using OpenCL

const kernel_source = """
__kernel void mmul(
	const int Mdim,
	const int Ndim,
	const int Pdim,
	__global float* A,
	__global float* B,
	__global float* C)
{
	int k;
	int i = get_global_id(0);
	int j = get_global_id(1);
	float tmp;
	if ((i < Ndim) && (j < Mdim))
	{
		tmp = 0.0f;
		for (k = 0; k < Pdim; k++)
			tmp += A[i*Ndim+k] * B[k*Pdim+j];
		C[i*Ndim+j] = tmp;
	}
}
"""

#### Definitions ###

# Order of the square matrices A, B and C
ORDER = 512

# A elemetns are constant and equal to AVAL
AVAL = 3.0

# B elemetns are constant and equal to BVAL
BVAL = 5.0

# tolerance used in floating point comparisons
TOL = 0.001

# Max dim for NDRange
DIM = 2

# number of times to do each multiplication
COUNT = 1

# Helper functions
include("helper.jl")

# A[N,P], B[P M], C[N,M]
Ndim = ORDER
Pdim = ORDER
Mdim = ORDER

# Number of elements in the matrix
sizeA = Ndim * Pdim
sizeB = Pdim * Mdim
sizeC = Ndim * Mdim

# Number of elements in the matrix
h_A = fill(Float32(AVAL), sizeA)
h_B = fill(Float32(BVAL), sizeB)
h_C = Vector{Float32}(sizeC)

# %20 improvment using @inbounds
function seq_mat_mul_sdot{T}(Mdim::Int, Ndim::Int, Pdim::Int,
                             A::Array{T}, B::Array{T}, C::Array{T})
    for i in 1:Ndim
        for j in 1:Mdim
            tmp = zero(Float32)
            for k in 1:Pdim
                @inbounds tmp += A[(i-1)*Ndim+k] * B[(k-1)*Pdim+j]
            end
            @inbounds C[(i-1)*Ndim+j] = tmp
        end
    end
end

info("=== Julia, matix mult (dot prod), order $ORDER ===")

# force compilation
seq_mat_mul_sdot(Mdim, Ndim, Pdim, h_A, h_B, h_C)

for i in 1:COUNT
    fill!(h_C, 0.0)
    t1 = time()
    #seq_mat_mul_sdot(Mdim, Ndim, Pdim, h_A, h_B, h_C)
    t2 = time()
    results(Mdim, Ndim, Pdim, h_C, t2 - t1)
end

# set up OpenCL
ctx = cl.create_some_context()

# You can enable profiling events on the queue
# by calling the constructor with the :profile flag
queue = cl.CmdQueue(ctx, :profile)

# create OpenCL Buffers
d_a = cl.Buffer(Float32, ctx, (:r,:copy), hostbuf=h_A)
d_b = cl.Buffer(Float32, ctx, (:r,:copy), hostbuf=h_B)
d_c = cl.Buffer(Float32, ctx, :w, length(h_C))

prg  = cl.Program(ctx, source=kernel_source) |> cl.build!
mmul = cl.Kernel(prg, "mmul")

info("=== OpenCL, matrix mult, C(i, j) per work item, order $Ndim ====")

for i in 1:COUNT
    fill!(h_C, 0.0)

    global_range = (Ndim, Mdim)
    mmul_ocl = mmul[queue, global_range]

    evt = mmul_ocl(Int32(Mdim), Int32(Ndim), Int32(Pdim),
                           d_a, d_b, d_c)

    # profiling events are measured in ns
    run_time = evt[:profile_duration] / 1e9
    cl.copy!(queue, h_C, d_c)
    results(Mdim, Ndim, Pdim, h_C, run_time)
end


text:
using DelimitedFiles

mutable struct Model
    diffusion_coeff::Array{Float64,2}       # Diffusion coefficient: depends on the sine of latitude
    heat_capacity::Array{Float64,3}         # Heat capacity: depends on the geography (land, ocean, ice, etc.). Position [:,:,1]: Previous time step, position [:,:,2]: current time step.
    albedo::Array{Float64,2}                # Albedo coefficient: depends on the geography (land, ocean, ice, etc.)
    solar_forcing::Array{Float64,3}         # Time-dependent incoming solar radiation: depends on the orbital parameters and the albedo
    radiative_cooling_co2::Array{Float64,1} # Outgoing long-wave radiation: depends on the CO2 concentration. Position 1: Previous time step, position 2: current time step.
    radiative_cooling_feedback::Float64     # Outgoing long-wave radiation (feedback effects): models the water vapor cyces, lapse rate and cloud cover
    co_albedo::Array{Float64,2}             # Absorbed radiatian coefficient: depends on the albedo
    compute_albedo::Bool                    # Attribute for deciding whether the albedo should be calculated instead of being read from a file
    geography::Array{Int8,2}                # Monthly land-sea-snow/sea ice distribution: depends on initial geography input if compute_sea_ice_extent is true
    solar_irradiance::Array{Float64,2}      # Solar irradiance: depends on the orbital parameters
    compute_sea_ice_extent::Bool            # Attribute for deciding whether sea ice extent should be calculated for the geography or read from given maps
    year::Int64
  end

  function Model(mesh, num_steps_year; co2_concentration = 315.0, compute_albedo = false, compute_sea_ice_extent = false, sea_ice_extent_year = 1979) # co2_concentration in [ppm], default value from year 1950
    @unpack nx,ny = mesh

    # Constants
    radiative_cooling_feedback = 2.15   #[W/m^2/°C]: sensitivity of the seasonal cycle and annual change in the forcing agents

    radiative_cooling_co2 = zeros(Float64,2)
    radiative_cooling_co2[1] = calc_radiative_cooling_co2(co2_concentration) # Initialize both entries of radiative_cooling_co2 to the same value
    radiative_cooling_co2[2] = radiative_cooling_co2[1]

    # Read parameters

    geography = read_geography(joinpath(@__DIR__, "..", "input", "world", string("The_World", string(nx, "x", ny),".dat")),nx,ny)

    compute_sea_ice_extent = compute_sea_ice_extent
    year = sea_ice_extent_year

    compute_albedo = compute_albedo

    if compute_albedo == false
      albedo = read_albedo(joinpath(@__DIR__, "..", "input", "albedo", string("albedo", string(nx, "x", ny), ".dat")),nx,ny)
    else

      albedo = calc_albedo(geography,nx,ny)
    end

    diffusion_coeff = calc_diffusion_coefficients(geography,nx,ny)

    # Compute heat capacity (we initialize the heat capacity in the previous time step as the current heat capacity)
    heat_capacity = zeros(nx,ny,2)
    calc_heat_capacity!(view(heat_capacity,:,:,1),geography,radiative_cooling_feedback) # TODO: Remove unused variables
    heat_capacity[:,:,2] = heat_capacity[:,:,1]

    co_albedo = 1.0 .- albedo

    ecc = eccentricity(1950)
    ob = obliquity(1950)
    per = perihelion(1950)

    solar_irradiance = calc_solar_forcing_params(nlongitude=nx, nlatitude=ny, ntimesteps=num_steps_year, ecc=ecc, ob=ob, per=per)
    solar_forcing = calc_solar_forcing(co_albedo, solar_irradiance, nx, ny, num_steps_year)

    return Model(diffusion_coeff, heat_capacity, albedo, solar_forcing, radiative_cooling_co2, radiative_cooling_feedback,co_albedo,compute_albedo,geography,solar_irradiance,compute_sea_ice_extent,year)
  end

  function set_co2_concentration!(model, co2_concentration)
    model.radiative_cooling_co2[1] = model.radiative_cooling_co2[2]
    model.radiative_cooling_co2[2] = calc_radiative_cooling_co2(co2_concentration)
  end

 """
  update_model()
  Gives the possibility to interactively decide which parameters should be updated in a simulation
 """
 function update_model(model, mesh, time_step, update_geography=true, update_albedo=true, update_diffusion=true,update_heat_capacity=true, update_solar_forcing=true)
  @unpack nx,ny = mesh

  if update_geography == true
    set_geography!(model,mesh, time_step)
  end

  if update_albedo == true
    set_albedo!(model,mesh)
    set_co_albedo!(model)
  end

  if update_diffusion == true
    set_diffusion_coeff!(model, mesh)
  end
  if update_heat_capacity == true
    set_heat_capacity!(model)
  end
  if update_solar_forcing == true
    set_solar_forcing!(model, mesh)
  end
  end

  # TODO: this routine works only for NT = 48

  function set_geography!(model, mesh, time_step)
  @unpack nx,ny = mesh

  # monthly geography can be updated by reading the input maps or by computing the sea ice extent
  # first time step is the last week of March as the simulation starts at the vernal equinox
   if model.compute_sea_ice_extent == false
      if time_step in (1,46,47,48)
        model.geography = read_geography(joinpath(@__DIR__, "..", "input","world", "year2004", string("The_World_from_image", nx, "x", ny,"_1")),nx,ny)

      elseif mod(time_step,4) == 2 && time_step in (2:45)
       month::Int64 = 1
        for nt in 1:time_step
          if mod(nt,4) == 1
           month +=1
          end
        end
        model.geography = read_geography(joinpath(@__DIR__, "..", "input","world", "year2004", string("The_World_from_image", nx, "x", ny,"_",month)),nx,ny)
      end
   else
     annual_sea_ice_extent = read_sea_ice_extent_N(model.year)
     calc_geography_per_month_extent!(model.geography,annual_sea_ice_extent, time_step, nx, ny)
   end
  end

  function set_albedo!(model,mesh)
  @unpack nx,ny = mesh
  if model.compute_albedo == false
    model.albedo = read_albedo(joinpath(@__DIR__, "..", "input", "albedo.dat"),nx,ny)
  else
    model.albedo = calc_albedo(model.geography,nx,ny)
  end
  end

  function set_co_albedo!(model)
  model.co_albedo = 1.0 .- model.albedo
  end

  function set_diffusion_coeff!(model, mesh)
  @unpack nx,ny = mesh
  model.diffusion_coeff = calc_diffusion_coefficients(model.geography,nx,ny)
  end

  function set_heat_capacity!(model)
  model.heat_capacity[:,:,1] = model.heat_capacity[:,:,2]
  calc_heat_capacity!(view(model.heat_capacity,:,:,2), model.geography, model.radiative_cooling_feedback)
  end

  function set_solar_forcing!(model, mesh)
  @unpack nx,ny = mesh
  solar_irradiance = model.solar_irradiance
  model.solar_forcing = calc_solar_forcing(model.co_albedo, solar_irradiance, nx, ny)
  end

  Base.size(model::Model) = size(model.heat_capacity)

  function Base.show(io::IO, model::Model)
  nx, ny = size(model)
  print(io, "Model() with ", nx, "×", ny, " degrees of freedom")
  end

  """
  calc_radiative_cooling_co2()
  Computes the CO2 parameter depending on the CO2 concentration

  * Default CO2 concentration is 315 ppm (equivalent to year 1950)
  """
  function calc_radiative_cooling_co2(co2_concentration=315.0)

  # Define base values for co2_concentration and radiative_cooling_co2
  co2_concentration_base = 315.0
  radiative_cooling_co2_base = 210.3

  # Doesn't change as long as CP2ppm doesn't change
  radiative_cooling_co2=radiative_cooling_co2_base-5.35*log(co2_concentration/co2_concentration_base)

  return radiative_cooling_co2
  end

  """
  _calc_insolation()
  auxiliar function of calc_solar_forcing
  """
  function _calc_insolation(dt, ob, ecc, per, nt, nlat, siny, cosy, tany, s0)

  eccfac = 1.0 - ecc^2
  rzero  = (2.0*pi)/eccfac^1.5

  #  Solve the orbital equation for lambda as a function of time with
  #  a fourth-order Runge-Kutta method
  lambda = zeros(Float64,nt+1)
  solar = zeros(Float64,nlat,nt)

  for n in 2:nt+1
    nu = lambda[n-1] - per
    t1 = dt*(rzero*(1.0 - ecc * cos(nu))^2)
    t2 = dt*(rzero*(1.0 - ecc * cos(nu+0.5*t1))^2)
    t3 = dt*(rzero*(1.0 - ecc * cos(nu+0.5*t2))^2)
    t4 = dt*(rzero*(1.0 - ecc * cos(nu + t3))^2)
    lambda[n] = lambda[n-1] + (t1 + 2.0*t2 + 2.0*t3 + t4)/6.0
  end

  #  Compute the average daily solar irradiance as a function of
  #  latitude and longitude (time)

  for n in 1:nt
    nu = lambda[n] - per
    rhofac = ((1.0- ecc*cos(nu))/eccfac)^2
    sindec = sin(ob)*sin(lambda[n])
    cosdec = sqrt(1.0-sindec^2)
    tandec = sindec/cosdec
    for j in 1:nlat
      z = -tany[j]*tandec
      if z >= 1.0    # polar latitudes when there is no sunrise (winter)
        solar[j,n] = 0.0
      else
        if z <= -1.0                # when there is no sunset (summer)
          solar[j,n] = rhofac * s0 * siny[j] * sindec
        else
          h_zero = acos(z)
          solar[j,n] = rhofac/pi*s0*(h_zero*siny[j]*sindec+cosy[j]*cosdec*sin(h_zero))
        end
      end
    end
  end

  return lambda, solar
  end

  """
  calc_solar_forcing_params()

  * Default s0 is 1371.685 [W/m²] (Current solar constant)
  * Default orbital parameters of correspond to year 1950 AD:
    ecc = 0.016740
    ob  = 0.409253
    per = 1.783037
  """
  function calc_solar_forcing_params(yr=0; nlongitude=nx, nlatitude=ny, ntimesteps=48, solar_cycle=false, s0=1371.685, orbital=false, ecc=0.016740, ob=0.409253, per=1.783037)
  # Calculate the sin, cos, and tan of the latitudes of Earth from the
  # colatitudes, calculate the insolation

  # nlatitude   = ny
  # nlongitude  = nx
  # ntimesteps  = num_steps_year

  dy = pi/(nlatitude-1.0)
  dt = 1.0 / ntimesteps

  siny = zeros(Float64,nlatitude)
  cosy = zeros(Float64,nlatitude)
  tany = zeros(Float64,nlatitude)

  if yr == 0
    for j in 1:nlatitude
      lat = pi/2.0 - dy*(j-1)    # latitude in radians
      siny[j] = sin(lat)
      if j == 1
        cosy[j] = 0.0
        tany[j] = 1000.0
      elseif j == nlatitude
        cosy[j] = 0.0
        tany[j] = -1000.0
      else
        cosy[j] = cos(lat)
        tany[j] = tan(lat)
      end
    end
    lambda, solar = _calc_insolation(dt, ob, ecc, per, ntimesteps, nlatitude, siny, cosy, tany, s0)
  elseif yr > 0 && (solar_cycle || orbital)
    lambda, solar = _calc_insolation(dt, ob, ecc, per, ntimesteps, nlatitude, siny, cosy, tany, s0)
  end

  return solar

  end

  function calc_solar_forcing(co_albedo, solar, nlongitude=128, nlatitude=65, ntimesteps=48)

  solar_forcing = zeros(Float64,nlongitude,nlatitude,ntimesteps)
  # calcualte the seasonal forcing
  for ts in 1:ntimesteps
    for j in 1:nlatitude
      for i in 1:nlongitude
        solar_forcing[i,j,ts] = solar[j,ts]*co_albedo[i,j]
      end
    end
  end
  return solar_forcing

  end

  function calc_heat_capacity!(heatcap,geography,radiative_cooling_feedback=2.15)

  # Depths
  depth_atmos = 5000.         # meters #TODO: not used???
  depth_mixed_layer = 70.0    # meters
  depth_soil = 2.0            # meters
  depth_seaice = 2.5          # meters
  depth_snow = 2.0            # meters
  layer_depth = 0.5           # kilometers

  # Physical properties of atmosphere
  rho_atmos = 1.293           # kg m^-3  dry air (STP)
  csp_atmos = 1005.0          # J kg^-1 K^-1 (STP)
  scale_height = 7.6          # kilometers

  # Physical properties of water
  rho_water = 1000.0          # kg m^-3
  csp_water = 4186.0          # J kg^-1 K^-1

  # Physical properties of soil
  rho_soil = 1100.0           # kg m^-3
  csp_soil = 850.0            # J kg^-1 K^-1

  # Physical properties of sea ice
  rho_sea_ice = 917.0         # kg m^-3
  csp_sea_ice = 2106.0        # J kg^-1 K^-1

  # Physical properties of snow covered surface
  rho_snow = 400.0            # kg m^-3
  csp_snow = 1900.0           # J kg^-1 K^-1

  # Other constants
  sec_per_yr = 3.15576e7      # seconds per year
  days_per_yr = 365.2422      # days per year

  # atmosphere with exponentially decaying density
  sum = 0.0
  for n in 1:10
    z = (0.25 + layer_depth*real(n-1))/scale_height
    sum = sum + exp(-z)
  end

  c_atmos  	= csp_atmos*layer_depth*1000.0*rho_atmos*sum/sec_per_yr
  c_soil   	= depth_soil*rho_soil*csp_soil/sec_per_yr
  c_seaice 	= depth_seaice*rho_sea_ice*csp_sea_ice/sec_per_yr
  c_snow   	= depth_snow * rho_snow * csp_snow/sec_per_yr
  c_mixed_layer = depth_mixed_layer*rho_water*csp_water/sec_per_yr

  # Assign the correct value of the heat capacity of the columns
  for j in 1:size(geography,2)
    for i in 1:size(geography,1)
      geo  = geography[i,j]
      if geo == 1                            # land
        heatcap[i,j] = c_soil + c_atmos
      elseif geo == 2                        # perennial sea ice
        heatcap[i,j] = c_seaice + c_atmos
      elseif geo == 3                        # permanent snow cover
        heatcap[i,j] = c_snow + c_atmos
      elseif geo == 4                        # lakes, inland seas
        heatcap[i,j] = c_mixed_layer/3.0 + c_atmos
      elseif geo in (5, 6, 7, 8)
        # 5: Pacific ocean
        # 6: Atlantic ocean
        # 7: Indian ocean
        # 8: Mediterranean
        heatcap[i,j] = c_mixed_layer + c_atmos
      end
    end
  end

  end

  """
  calc_diffusion_coefficients()
  Calculate the diffusion coefficients at finest grid level.
  """
  function calc_diffusion_coefficients(geography,nlongitude=128,nlatitude=65)

  coeff_eq     = 0.65 # coefficinet for diffusion at equator
  coeff_ocean  = 0.40 # coefficient for ocean diffusion
  coeff_land   = 0.65 # coefficinet for land diffusion
  coeff_landNP = 0.28 # coefficinet for land diffusion (north pole)
  coeff_landSP = 0.20 # coefficinet for land diffusion (south pole)

  diffusion = zeros(Float64,nlongitude,nlatitude)

  j_equator = div(nlatitude,2) + 1

  for j = 1:nlatitude
      theta = pi*real(j-1)/real(nlatitude-1)
      colat = sin(theta)^5

      for i = 1:nlongitude
          let geo = geography[i,j]
              if geo >= 5 && geo <= 7 # oceans
                  diffusion[i,j] = (coeff_eq-coeff_ocean)*colat + coeff_ocean
              else # land, sea ice, etc
                  if j <= j_equator # northern hemisphere
                      diffusion[i,j] = (coeff_land-coeff_landNP)*colat + coeff_landNP
                  else # southern hemisphere
                      diffusion[i,j] = (coeff_land-coeff_landSP)*colat + coeff_landSP
                  end
                end
          end
      end
  end

  return diffusion
end

function read_albedo(filepath=joinpath(@__DIR__,"..","input","albedo","albedo128x65.dat"),nlongitude=128,nlatitude=65)

  result = zeros(Float64,nlongitude,nlatitude)
  open(filepath) do fh
      for lat = 1:nlatitude
          if eof(fh) break end
         result[:,lat] = parse.(Float64,split(strip(readline(fh) ),r"\s+"))
     end
  end
  return result
  end

  #Calculate the albedo depending on the geography

  function calc_albedo(geography,nlongitude=128,nlatitude=65)
  result = zeros(Float64,nlongitude,nlatitude)
  dtheta= pi/(nlatitude-1.0)

  for j in 1:nlatitude
    theta = 0.5*pi-dtheta*(j-1)
    sintheta = sin(theta)
    legendrepol = (3*sintheta^2-1)*0.5
    for i in 1:nlongitude
      let geo = geography[i,j]
        if geo == 1 # land
          result[i,j] = 0.3+0.12*legendrepol
        elseif geo == 2 # sea ice
          result[i,j] = 0.6
        elseif geo == 3 # snow cover
          result[i,j] = 0.75
        elseif geo in (5,6,7,8) # oceans
          result[i,j] = 0.29+0.12*legendrepol
        end
      end
    end
  end
  return result
end

function read_geography(filepath=joinpath(@__DIR__,"..","input","world","The_World128x65.dat"),nlongitude=128,nlatitude=65)
  result = zeros(Int8,nlongitude,nlatitude)
  open(filepath) do fh
      for lat = 1:nlatitude
          if eof(fh) break end
          result[:,lat] = parse.(Int8,split(strip(readline(fh) ),r""))
      end
  end
  return result
  end

  # Table 4 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  # Higher accuracy values from Gary Russel, https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html
  const eccentricity_parameters_amplitude = [0.01860798, 0.01627522, -0.01300660, 0.00988829,
                                           -0.00336700, 0.00333077, -0.00235400, 0.00140015,
                                           0.00100700, 0.00085700, 0.00064990, 0.00059900,
                                           0.00037800, -0.00033700, 0.00027600, 0.00018200,
                                           -0.00017400, -0.00012400, 0.00001250]
  const eccentricity_parameters_mean_rate = [4.2072050, 7.3460910, 17.8572630, 17.2205460, 16.8467330,
                                           5.1990790, 18.2310760, 26.2167580, 6.3591690, 16.2100160,
                                           3.0651810, 16.5838290, 18.4939800, 6.1909530, 18.8677930,
                                           17.4255670, 6.1860010, 18.4174410, 0.6678630]
  const eccentricity_parameters_phase = [28.620089, 193.788772, 308.307024, 320.199637, 279.376984,
                                       87.195000, 349.129677, 128.443387, 154.143880, 291.269597,
                                       114.860583, 332.092251, 296.414411, 145.769910, 337.237063,
                                       152.092288, 126.839891, 210.667199, 72.108838]

  # Table 1 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  # Higher accuracy values from Gary Russel, https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html
  const obliquity_parameters_amplitude = [-2462.2214466, -857.3232075, -629.3231835, -414.2804924,
                                        -311.7632587, 308.9408604, -162.5533601, -116.1077911,
                                        101.1189923, -67.6856209, 24.9079067, 22.5811241,
                                        -21.1648355, -15.6549876, 15.3936813, 14.6660938,
                                        -11.7273029, 10.2742696, 6.4914588, 5.8539148, -5.4872205,
                                        -5.4290191, 5.160957, 5.0786314, -4.0735782, 3.7227167,
                                        3.3971932, -2.8347004, -2.6550721, -2.5717867, -2.4712188,
                                        2.462541, 2.2464112, -2.0755511, -1.9713669, -1.8813061,
                                        -1.8468785, 1.8186742, 1.7601888, -1.5428851, 1.4738838,
                                        -1.4593669, 1.4192259, -1.181898, 1.1756474, -1.1316126,
                                        1.0896928]
  const obliquity_parameters_mean_rate = [31.609974, 32.620504, 24.172203, 31.983787, 44.828336,
                                        30.973257, 43.668246, 32.246691, 30.599444, 42.681324,
                                        43.836462, 47.439436, 63.219948, 64.230478, 1.01053,
                                        7.437771, 55.782177, 0.373813, 13.218362, 62.583231,
                                        63.593761, 76.43831, 45.815258, 8.448301, 56.792707,
                                        49.747842, 12.058272, 75.27822, 65.241008, 64.604291,
                                        1.647247, 7.811584, 12.207832, 63.856665, 56.15599,
                                        77.44884, 6.801054, 62.209418, 20.656133, 48.344406,
                                        55.14546, 69.000539, 11.07135, 74.291298, 11.047742,
                                        0.636717, 12.844549]
  const obliquity_parameters_phase = [251.9025, 280.8325, 128.3057, 292.7252, 15.3747, 263.7951,
                                    308.4258, 240.0099, 222.9725, 268.7809, 316.7998, 319.6024,
                                    143.805, 172.7351, 28.93, 123.5968, 20.2082, 40.8226, 123.4722,
                                    155.6977, 184.6277, 267.2772, 55.0196, 152.5268, 49.1382,
                                    204.6609, 56.5233, 200.3284, 201.6651, 213.5577, 17.0374,
                                    164.4194, 94.5422, 131.9124, 61.0309, 296.2073, 135.4894,
                                    114.875, 247.0691, 256.6114, 32.1008, 143.6804, 16.8784,
                                    160.6835, 27.5932, 348.1074, 82.6496]

  # Table 5 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  # Higher accuracy values from Gary Russel, https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html
  const perihelion_parameters_amplitude = [7391.022589, 2555.1526947, 2022.7629188, -1973.6517951,
                                         1240.2321818, 953.8679112, -931.7537108, 872.3795383,
                                         606.3544732, -496.0274038, 456.9608039, 346.946232,
                                         -305.8412902, 249.6173246, -199.10272, 191.0560889,
                                         -175.2936572, 165.9068833, 161.1285917, 139.7878093,
                                         -133.5228399, 117.0673811, 104.6907281, 95.3227476,
                                         86.7824524, 86.0857729, 70.5893698, -69.9719343,
                                         -62.5817473, 61.5450059, -57.9364011, 57.1899832,
                                         -57.0236109, -54.2119253, 53.2834147, 52.1223575,
                                         -49.0059908, -48.3118757, -45.4191685, -42.235792,
                                         -34.7971099, 34.4623613, -33.8356643, 33.6689362,
                                         -31.2521586, -30.8798701, 28.4640769, -27.1960802,
                                         27.0860736, -26.3437456, 24.725374, 24.6732126, 24.4272733,
                                         24.0127327, 21.7150294, -21.5375347, 18.1148363,
                                         -16.9603104, -16.1765215, 15.5567653, 15.4846529,
                                         15.2150632, 14.5047426, -14.3873316, 13.1351419,
                                         12.8776311, 11.9867234, 11.9385578, 11.7030822, 11.6018181,
                                         -11.2617293, -10.4664199, 10.433397, -10.2377466,
                                         10.1934446, -10.1280191, 10.0289441, -10.0034259]
  const perihelion_parameters_mean_rate = [31.609974, 32.620504, 24.172203, 0.636717, 31.983787,
                                         3.138886, 30.973257, 44.828336, 0.991874, 0.373813,
                                         43.668246, 32.246691, 30.599444, 2.147012, 10.511172,
                                         42.681324, 13.650058, 0.986922, 9.874455, 13.013341,
                                         0.262904, 0.004952, 1.142024, 63.219948, 0.205021,
                                         2.151964, 64.230478, 43.836462, 47.439436, 1.384343,
                                         7.437771, 18.829299, 9.500642, 0.431696, 1.16009,
                                         55.782177, 12.639528, 1.155138, 0.168216, 1.647247,
                                         10.884985, 5.610937, 12.658184, 1.01053, 1.983748,
                                         14.023871, 0.560178, 1.273434, 12.021467, 62.583231,
                                         63.593761, 76.43831, 4.28091, 13.218362, 17.818769,
                                         8.359495, 56.792707, 8.448301, 1.978796, 8.863925,
                                         0.186365, 8.996212, 6.771027, 45.815258, 12.002811,
                                         75.27822, 65.241008, 18.870667, 22.009553, 64.604291,
                                         11.498094, 0.578834, 9.237738, 49.747842, 2.147012,
                                         1.196895, 2.133898, 0.173168]
  const perihelion_parameters_phase = [251.9025, 280.8325, 128.3057, 348.1074, 292.7252, 165.1686,
                                     263.7951, 15.3747, 58.5749, 40.8226, 308.4258, 240.0099,
                                     222.9725, 106.5937, 114.5182, 268.7809, 279.6869, 39.6448,
                                     126.4108, 291.5795, 307.2848, 18.93, 273.7596, 143.805,
                                     191.8927, 125.5237, 172.7351, 316.7998, 319.6024, 69.7526,
                                     123.5968, 217.6432, 85.5882, 156.2147, 66.9489, 20.2082,
                                     250.7568, 48.0188, 8.3739, 17.0374, 155.3409, 94.1709, 221.112,
                                     28.93, 117.1498, 320.5095, 262.3602, 336.2148, 233.0046,
                                     155.6977, 184.6277, 267.2772, 78.9281, 123.4722, 188.7132,
                                     180.1364, 49.1382, 152.5268, 98.2198, 97.4808, 221.5376,
                                     168.2438, 161.1199, 55.0196, 262.6495, 200.3284, 201.6651,
                                     294.6547, 99.8233, 213.5577, 154.1631, 232.7153, 138.3034,
                                     204.6609, 106.5938, 250.4676, 332.3345, 27.3039]

  # Eqn. 4 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  function e_sincos_pi(year)
  M    = eccentricity_parameters_amplitude
  g    = eccentricity_parameters_mean_rate
  beta = eccentricity_parameters_phase

  arg = deg2rad.(g/3600) * (year-1950) .+ deg2rad.(beta)

  return sum(M .* sin.(arg)), sum(M .* cos.(arg))
  end

  eccentricity(year, e_sin_pi, e_cos_pi) = sqrt(e_sin_pi^2 + e_cos_pi^2)

  eccentricity(year) = eccentricity(year, e_sincos_pi(year)...)

  # Eqn. 1 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  function obliquity(year)
  A     = obliquity_parameters_amplitude
  f     = obliquity_parameters_mean_rate
  delta = obliquity_parameters_phase

  epsilon_star = 23.320556

  return deg2rad(epsilon_star + sum(A/3600 .* cos.(deg2rad.(f/3600) * (year-1950) .+ deg2rad.(delta))))
  end

  # Eqn. 6 from Berger 1978, Long-Term Variations of daily Insolation and Quaternary Climatic Changes
  function perihelion(year)
  F     = perihelion_parameters_amplitude
  f     = perihelion_parameters_mean_rate
  delta = perihelion_parameters_phase

  pi_ = atan(e_sincos_pi(year)...)
  psi_tilde = 50.439273
  zeta = 3.392506

  psi = psi_tilde/3600 * (year-1950) + zeta + sum(F/3600 .* sin.(deg2rad.(f/3600) * (year-1950) .+ deg2rad.(delta)))

  return (pi_ + deg2rad(psi)) % (2*pi)
  end

"""
read_sea_ice_extent_N()
* Default year: 1979
* Possible years to pick annual sea ice extent for the Northern Hemisphere : 1979 to 2020
* returns a vector with sea ice extent per month from March to February

"""
function read_sea_ice_extent_N(year=1979)

  matrix = readdlm(joinpath(@__DIR__,"..","input/sea_ice_data/","N_ALL_extent_v3.0.dat"),comments=true)
  yearindex = findfirst(isequal(year),matrix[:,1])
  annual_sea_ice_extent = matrix[yearindex,2:13]

  return annual_sea_ice_extent
end

"""
calc_geography_per_month_extent!()
* This routine a new geography matrix by determining sea ice extent with a monthly rate of change
  and mapping extent or reduction to adjacent sea ice cells and ocean cells on the Northern Hemisphere.

* This routine assumes that the sea ice extent of the geography file corresponds to
  the measured sea ice at the vernal equinox.

* There is a bias in the order sea ice cells are added/removed.
  sea ice extent: starts from top left to bottom right
  sea ice reduction: starts from bottom left to top right

"""
function calc_geography_per_month_extent!(geography,annual_sea_ice_extent, time_step, nlongitude=128, nlatitude =65)

  month::Int64 = 0

 # Gets the right month when the simulation starts at the vernal equinox
 # TODO: this part works only for NT = 48

 if time_step in (1,46,47,48)
  month = 1
elseif mod(time_step,4) == 2
    month = 1
    for nt in 1:time_step
      if mod(nt,4) == 1
        month +=1
      end
    end
 end

  nlat::Int64 = (nlatitude-1)*0.5
  geography_start = read_geography(joinpath(@__DIR__, "..", "input","world", "The_World128x65.dat"),128,65)
  geo = geography_start[:,1:nlat]

  if month == 0 || month == 1
    geography = geography_start
    return geography
  else

  sea_ice_index_unit = size(findall(isequal(2),geo),1)/annual_sea_ice_extent[1]

   # collect submatrix and indices for the Northern Hemisphere

   G = geography[:,1:nlat]
   N = reverse(G, dims = 2)
   IndG = CartesianIndices(G)
   IndN = CartesianIndices(N)
   Gfirst, Glast = first(IndG), last(IndG)
   Ifirst, Ilast = first(IndN), last(IndN)
   G1 = oneunit(Gfirst)
   I1 = oneunit(Ifirst)

   # get all sea ice cells for the Northern Hemispehre and determine how many cells will be affected by the extent/reduction

   sea_ice_index = findall(isequal(2),G)
   old_counter = size(sea_ice_index,1)
   new_counter = sea_ice_index_unit*annual_sea_ice_extent[month]
   sea_ice_cells_delta = abs(floor(Int,new_counter - old_counter))

   # case for an extent of sea ice
   # adjacent ocean cells (geo=5,6,7,8) will turn to sea ice cells
   if sea_ice_cells_delta == 0
     return geography
   end
   if new_counter > old_counter
      for I in IndG
          if G[I[1], I[2]] == 2
              if sea_ice_cells_delta <= 0
                  break
              else
                  for J in max(Gfirst, I-G1):min(Glast, I+G1)
                      if G[J] in (5,6,7,8)
                          G[J[1], J[2]] = 2
                          sea_ice_cells_delta -= 1
                         if sea_ice_cells_delta <= 0
                             break
                         end
                      end
                  end
              end
          end
      end
      geography[:,1:nlat] = G

   # case for a reduction of sea ice
   # adjacent sea ice cells will turn to ocean cells
   elseif new_counter < old_counter
      for I in IndN
          if N[I[1], I[2]] == 2 && sea_ice_cells_delta > 0
              N[I[1], I[2]] = 5
              sea_ice_cells_delta -= 1
              if sea_ice_cells_delta <= 0
                  break
              else
                  for J in max(Ifirst, I-I1):min(Ilast, I+I1)
                      if N[J] == 2
                          N[J[1], J[2]] = 5
                          sea_ice_cells_delta -= 1
                          if sea_ice_cells_delta <= 0
                              break
                          end
                      end
                  end

              end

          end
      end
      geography[:,1:nlat] = reverse(N, dims = 2)
   end

   return geography
  end

end


text:
# module TestCombinations

using Test, Dates, Temporal

X1 = copy(X)
X2 = TS(Y.values, X1.index+Week(1), Y.fields)
x1 = X1[:,1]
x2 = X1[:,1]
@testset "Combining" begin
    @testset "Joins" begin
        # outer joins
        Z = ojoin(X1, X2)
        @test size(Z,1) == length(union(X1.index, X2.index))
        @test size(Z,2) == size(X1,2) + size(X2,2)
        @test dropnan(Z) == dropnan(merge(X1, X2, join='o'))
        # inner joins
        Z = ijoin(X1, X2)
        @test size(Z,2) == size(X1,2) + size(X2,2)
        @test Z == merge(X1, X2, join='i')
        # right joins
        Z = rjoin(X1, X2)
        @test size(Z,1) == size(X1,1)
        @test dropnan(Z) == dropnan(merge(X1, X2, join='r'))
        # left joins
        Z = ljoin(X1,X2)
        @test size(Z,1) == size(X2,1)
        @test dropnan(Z) == dropnan(merge(X1, X2, join='l'))
    end
    @testset "Horizontal Concatenation" begin
        # two ts objects
        Z = [X1 X2]
        @test size(Z,1) == length(union(X1.index, X2.index))
        @test size(Z,2) == size(X2,2) + size(X2,2)
        # ts object with number
        Z = [X1 1.0]
        @test size(Z,1) == size(X1,1)
        @test size(Z,2) == size(X1,2) + 1
        # ts object with array
        Z = [X1 ones(size(X1))]
        @test size(Z,1) == size(X1,1)
        @test size(Z,2) == size(X1,2)*2
        # multiple ts objects
        Z = [X1 X2 x1 x2]
        @test size(Z,1) == length(union(X1.index, X2.index))
        @test size(Z,2) == size(X1,2) + size(X2,2) + size(x1,2) + size(x2,2)
        # ts object with multiple arrays
        Z = [X1 ones(N) 2.0.+zeros(N,2)]
        @test size(Z,1) == size(X1,1)
        @test size(Z,2) == size(X1,2) + 3
        @test all(Z[:,5].values .== 1.0)
        @test all(Z[:,6:end].values .== 2.0)
        # ts object with multiple numbers
        Z = [X1 1 2.0]
        @test size(Z,1) == size(X1,1)
        @test size(Z,2) == size(X1,2) + 2
        @test all(Z[:,5].values .== 1.0)
        @test all(Z[:,6].values .== 2.0)
    end
end

# end


text:
#=
P17 (*) Split a list into two parts; the length of the first part is given.
Do not use any predefined predicates.

Example:
* (split '(a b c d e f g h i k) 3)
( (A B C) (D E F G H I K))
=#

function my_split(coll, index)
    vcat(Vector{}[coll[1:index]], Vector{}[coll[index + 1:end]])
end

splitted = my_split([:a, :b, :c, :d, :e, :f, :g, :h, :i, :k], 3)
@assert splitted[1] == [:a, :b, :c]
@assert splitted[2] == [:d, :e, :f, :g, :h, :i, :k]

println("Tests passed: JL-17.jl")


text:
module sortedindices

"""
6, [1, 4] => [2, 3, 5, 6]
"""
function complement(N::Int, indices::Vector{Int})
    L = length(indices)
    x = Vector{Int}(undef, N - L)
    i_ = 1 # Position in the x vector
    j = 1 # Position in indices vector
    for i=1:N
        if j > L || indices[j]!=i
            x[i_] = i
            i_ += 1
        else
            j += 1
        end
    end
    x
end

"""
[1, 4, 5], [2, 4, 7] => [1, 2, 4, 5, 7]
"""
function union(ind1::Vector{Int}, ind2::Vector{Int})
    i1 = 1
    i2 = 1
    N1 = length(ind1)
    N2 = length(ind2)
    xvec = Vector{Int}()
    while true
        if i1 > N1
            for j=i2:N2
                push!(xvec, ind2[j])
            end
            return xvec
        elseif i2 > N2
            for j=i1:N1
                push!(xvec, ind1[j])
            end
            return xvec
        end
        x1 = ind1[i1]
        x2 = ind2[i2]
        if x1 == x2
            i1 += 1
            i2 += 1
            push!(xvec, x1)
        elseif x1 < x2
            i1 += 1
            push!(xvec, x1)
        else
            i2 += 1
            push!(xvec, x2)
        end
    end
end

"""
[1, 4, 5], [2, 4, 7] => [1, 5]
"""
function remove(ind1::Vector{Int}, ind2::Vector{Int})
    x = Int[]
    for i in ind1
        if i ∉ ind2
            push!(x, i)
        end
    end
    x
end

"""
[1, 4, 5], [2, 4, 7] => [1, 3]
"""
function shiftremove(ind1::Vector{Int}, ind2::Vector{Int})
    x = Int[]
    for i in ind1
        if i ∉ ind2
            counter = 0
            for i2 in ind2
                if i2 < i
                    counter += 1
                else
                    break
                end
            end
            push!(x, i-counter)
        end
    end
    x
end

"""
[2, 3, 6], [1, 3, 4, 6, 7] => [3, 6]
"""
function intersect(ind1::Vector{Int}, ind2::Vector{Int})
    i1 = 1
    i2 = 1
    N1 = length(ind1)
    N2 = length(ind2)
    xvec = Vector{Int}()
    if i1 > N1 || i2 > N2
        return xvec
    end
    x1 = ind1[i1]
    x2 = ind2[i2]
    while true
        if x1 == x2
            i1 += 1
            i2 += 1
            push!(xvec, x1)
            if i1 > N1 || i2 > N2
                return xvec
            end
            x1 = ind1[i1]
            x2 = ind2[i2]
        elseif x1 < x2
            i1 += 1
            if i1 > N1
                return xvec
            end
            x1 = ind1[i1]
        else
            i2 += 1
            if i2 > N2
                return xvec
            end
            x2 = ind2[i2]
        end
    end
end

function reducedindices(I_::Vector{Int}, I::Vector{Int})
    N = length(I_)
    x = Vector{Int}(undef, N)
    for n in 1:N
        x[n] = findfirst(isequal(I_[n]), I)
    end
    x
end

function reducedindices!(I_::Vector{Int}, I::Vector{Int})
    for n in 1:length(I_)
        I_[n] = findfirst(isequal(I_[n]), I)
    end
end

"""
Check if all indices are unique and smaller than or equal to imax.
"""
function check_indices(imax::Int, indices::Vector{Int})
    N = length(indices)
    for n=1:N
        i = indices[n]
        @assert 0 < i <= imax
        for m in n+1:N
            @assert i != indices[m]
        end
    end
end

"""
Check if the indices are sorted, unique and smaller than or equal to imax.
"""
function check_sortedindices(imax::Int, indices::Vector{Int})
    N = length(indices)
    if N == 0
        return nothing
    end
    i_ = indices[1]
    @assert 0 < i_ <= imax
    for i in indices[2:end]
        @assert 0 < i <= imax
        @assert i > i_
    end
end

end # module


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: cloudhsm

using AWS.Compat
using AWS.UUIDs
"""
    AddTagsToResource()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Adds or overwrites one or more tags for the specified AWS CloudHSM resource. Each tag consists of a key and a value. Tag keys must be unique to each resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the AWS CloudHSM resource to tag.
- `TagList`: One or more tags.

"""

add_tags_to_resource(ResourceArn, TagList; aws_config::AWSConfig=global_aws_config()) = cloudhsm("AddTagsToResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagList"=>TagList); aws_config=aws_config)
add_tags_to_resource(ResourceArn, TagList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("AddTagsToResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagList"=>TagList), args)); aws_config=aws_config)

"""
    CreateHapg()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Creates a high-availability partition group. A high-availability partition group is a group of partitions that spans multiple physical HSMs.

# Required Parameters
- `Label`: The label of the new high-availability partition group.

"""

create_hapg(Label; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateHapg", Dict{String, Any}("Label"=>Label); aws_config=aws_config)
create_hapg(Label, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateHapg", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Label"=>Label), args)); aws_config=aws_config)

"""
    CreateHsm()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Creates an uninitialized HSM instance. There is an upfront fee charged for each HSM instance that you create with the CreateHsm operation. If you accidentally provision an HSM and want to request a refund, delete the instance using the DeleteHsm operation, go to the AWS Support Center, create a new case, and select Account and Billing Support.  It can take up to 20 minutes to create and provision an HSM. You can monitor the status of the HSM with the DescribeHsm operation. The HSM is ready to be initialized when the status changes to RUNNING.

# Required Parameters
- `IamRoleArn`: The ARN of an IAM role to enable the AWS CloudHSM service to allocate an ENI on your behalf.
- `SshKey`: The SSH public key to install on the HSM.
- `SubnetId`: The identifier of the subnet in your VPC in which to place the HSM.
- `SubscriptionType`:

# Optional Parameters
- `ClientToken`: A user-defined token to ensure idempotence. Subsequent calls to this operation with the same token will be ignored.
- `EniIp`: The IP address to assign to the HSM's ENI. If an IP address is not specified, an IP address will be randomly chosen from the CIDR range of the subnet.
- `ExternalId`: The external ID from IamRoleArn, if present.
- `SyslogIp`: The IP address for the syslog monitoring server. The AWS CloudHSM service only supports one syslog monitoring server.
"""

create_hsm(IamRoleArn, SshKey, SubnetId, SubscriptionType; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateHsm", Dict{String, Any}("IamRoleArn"=>IamRoleArn, "SshKey"=>SshKey, "SubnetId"=>SubnetId, "SubscriptionType"=>SubscriptionType); aws_config=aws_config)
create_hsm(IamRoleArn, SshKey, SubnetId, SubscriptionType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateHsm", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IamRoleArn"=>IamRoleArn, "SshKey"=>SshKey, "SubnetId"=>SubnetId, "SubscriptionType"=>SubscriptionType), args)); aws_config=aws_config)

"""
    CreateLunaClient()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Creates an HSM client.

# Required Parameters
- `Certificate`: The contents of a Base64-Encoded X.509 v3 certificate to be installed on the HSMs used by this client.

# Optional Parameters
- `Label`: The label for the client.
"""

create_luna_client(Certificate; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateLunaClient", Dict{String, Any}("Certificate"=>Certificate); aws_config=aws_config)
create_luna_client(Certificate, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("CreateLunaClient", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Certificate"=>Certificate), args)); aws_config=aws_config)

"""
    DeleteHapg()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Deletes a high-availability partition group.

# Required Parameters
- `HapgArn`: The ARN of the high-availability partition group to delete.

"""

delete_hapg(HapgArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteHapg", Dict{String, Any}("HapgArn"=>HapgArn); aws_config=aws_config)
delete_hapg(HapgArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteHapg", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HapgArn"=>HapgArn), args)); aws_config=aws_config)

"""
    DeleteHsm()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Deletes an HSM. After completion, this operation cannot be undone and your key material cannot be recovered.

# Required Parameters
- `HsmArn`: The ARN of the HSM to delete.

"""

delete_hsm(HsmArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteHsm", Dict{String, Any}("HsmArn"=>HsmArn); aws_config=aws_config)
delete_hsm(HsmArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteHsm", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HsmArn"=>HsmArn), args)); aws_config=aws_config)

"""
    DeleteLunaClient()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Deletes a client.

# Required Parameters
- `ClientArn`: The ARN of the client to delete.

"""

delete_luna_client(ClientArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteLunaClient", Dict{String, Any}("ClientArn"=>ClientArn); aws_config=aws_config)
delete_luna_client(ClientArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DeleteLunaClient", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientArn"=>ClientArn), args)); aws_config=aws_config)

"""
    DescribeHapg()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Retrieves information about a high-availability partition group.

# Required Parameters
- `HapgArn`: The ARN of the high-availability partition group to describe.

"""

describe_hapg(HapgArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeHapg", Dict{String, Any}("HapgArn"=>HapgArn); aws_config=aws_config)
describe_hapg(HapgArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeHapg", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HapgArn"=>HapgArn), args)); aws_config=aws_config)

"""
    DescribeHsm()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Retrieves information about an HSM. You can identify the HSM by its ARN or its serial number.

# Optional Parameters
- `HsmArn`: The ARN of the HSM. Either the HsmArn or the SerialNumber parameter must be specified.
- `HsmSerialNumber`: The serial number of the HSM. Either the HsmArn or the HsmSerialNumber parameter must be specified.
"""

describe_hsm(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeHsm"; aws_config=aws_config)
describe_hsm(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeHsm", args; aws_config=aws_config)

"""
    DescribeLunaClient()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Retrieves information about an HSM client.

# Optional Parameters
- `CertificateFingerprint`: The certificate fingerprint.
- `ClientArn`: The ARN of the client.
"""

describe_luna_client(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeLunaClient"; aws_config=aws_config)
describe_luna_client(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("DescribeLunaClient", args; aws_config=aws_config)

"""
    GetConfig()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Gets the configuration files necessary to connect to all high availability partition groups the client is associated with.

# Required Parameters
- `ClientArn`: The ARN of the client.
- `ClientVersion`: The client version.
- `HapgList`: A list of ARNs that identify the high-availability partition groups that are associated with the client.

"""

get_config(ClientArn, ClientVersion, HapgList; aws_config::AWSConfig=global_aws_config()) = cloudhsm("GetConfig", Dict{String, Any}("ClientArn"=>ClientArn, "ClientVersion"=>ClientVersion, "HapgList"=>HapgList); aws_config=aws_config)
get_config(ClientArn, ClientVersion, HapgList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("GetConfig", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientArn"=>ClientArn, "ClientVersion"=>ClientVersion, "HapgList"=>HapgList), args)); aws_config=aws_config)

"""
    ListAvailableZones()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Lists the Availability Zones that have available AWS CloudHSM capacity.

"""

list_available_zones(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListAvailableZones"; aws_config=aws_config)
list_available_zones(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListAvailableZones", args; aws_config=aws_config)

"""
    ListHapgs()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Lists the high-availability partition groups for the account. This operation supports pagination with the use of the NextToken member. If more results are available, the NextToken member of the response contains a token that you pass in the next call to ListHapgs to retrieve the next set of items.

# Optional Parameters
- `NextToken`: The NextToken value from a previous call to ListHapgs. Pass null if this is the first call.
"""

list_hapgs(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListHapgs"; aws_config=aws_config)
list_hapgs(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListHapgs", args; aws_config=aws_config)

"""
    ListHsms()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Retrieves the identifiers of all of the HSMs provisioned for the current customer. This operation supports pagination with the use of the NextToken member. If more results are available, the NextToken member of the response contains a token that you pass in the next call to ListHsms to retrieve the next set of items.

# Optional Parameters
- `NextToken`: The NextToken value from a previous call to ListHsms. Pass null if this is the first call.
"""

list_hsms(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListHsms"; aws_config=aws_config)
list_hsms(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListHsms", args; aws_config=aws_config)

"""
    ListLunaClients()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Lists all of the clients. This operation supports pagination with the use of the NextToken member. If more results are available, the NextToken member of the response contains a token that you pass in the next call to ListLunaClients to retrieve the next set of items.

# Optional Parameters
- `NextToken`: The NextToken value from a previous call to ListLunaClients. Pass null if this is the first call.
"""

list_luna_clients(; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListLunaClients"; aws_config=aws_config)
list_luna_clients(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListLunaClients", args; aws_config=aws_config)

"""
    ListTagsForResource()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Returns a list of all tags for the specified AWS CloudHSM resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the AWS CloudHSM resource.

"""

list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListTagsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    ModifyHapg()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Modifies an existing high-availability partition group.

# Required Parameters
- `HapgArn`: The ARN of the high-availability partition group to modify.

# Optional Parameters
- `Label`: The new label for the high-availability partition group.
- `PartitionSerialList`: The list of partition serial numbers to make members of the high-availability partition group.
"""

modify_hapg(HapgArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyHapg", Dict{String, Any}("HapgArn"=>HapgArn); aws_config=aws_config)
modify_hapg(HapgArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyHapg", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HapgArn"=>HapgArn), args)); aws_config=aws_config)

"""
    ModifyHsm()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Modifies an HSM.  This operation can result in the HSM being offline for up to 15 minutes while the AWS CloudHSM service is reconfigured. If you are modifying a production HSM, you should ensure that your AWS CloudHSM service is configured for high availability, and consider executing this operation during a maintenance window.

# Required Parameters
- `HsmArn`: The ARN of the HSM to modify.

# Optional Parameters
- `EniIp`: The new IP address for the elastic network interface (ENI) attached to the HSM. If the HSM is moved to a different subnet, and an IP address is not specified, an IP address will be randomly chosen from the CIDR range of the new subnet.
- `ExternalId`: The new external ID.
- `IamRoleArn`: The new IAM role ARN.
- `SubnetId`: The new identifier of the subnet that the HSM is in. The new subnet must be in the same Availability Zone as the current subnet.
- `SyslogIp`: The new IP address for the syslog monitoring server. The AWS CloudHSM service only supports one syslog monitoring server.
"""

modify_hsm(HsmArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyHsm", Dict{String, Any}("HsmArn"=>HsmArn); aws_config=aws_config)
modify_hsm(HsmArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyHsm", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HsmArn"=>HsmArn), args)); aws_config=aws_config)

"""
    ModifyLunaClient()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Modifies the certificate used by the client. This action can potentially start a workflow to install the new certificate on the client's HSMs.

# Required Parameters
- `Certificate`: The new certificate for the client.
- `ClientArn`: The ARN of the client.

"""

modify_luna_client(Certificate, ClientArn; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyLunaClient", Dict{String, Any}("Certificate"=>Certificate, "ClientArn"=>ClientArn); aws_config=aws_config)
modify_luna_client(Certificate, ClientArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("ModifyLunaClient", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Certificate"=>Certificate, "ClientArn"=>ClientArn), args)); aws_config=aws_config)

"""
    RemoveTagsFromResource()

This is documentation for AWS CloudHSM Classic. For more information, see AWS CloudHSM Classic FAQs, the AWS CloudHSM Classic User Guide, and the AWS CloudHSM Classic API Reference.  For information about the current version of AWS CloudHSM, see AWS CloudHSM, the AWS CloudHSM User Guide, and the AWS CloudHSM API Reference. Removes one or more tags from the specified AWS CloudHSM resource. To remove a tag, specify only the tag key to remove (not the value). To overwrite the value for an existing tag, use AddTagsToResource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the AWS CloudHSM resource.
- `TagKeyList`: The tag key or keys to remove. Specify only the tag key to remove (not the value). To overwrite the value for an existing tag, use AddTagsToResource.

"""

remove_tags_from_resource(ResourceArn, TagKeyList; aws_config::AWSConfig=global_aws_config()) = cloudhsm("RemoveTagsFromResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeyList"=>TagKeyList); aws_config=aws_config)
remove_tags_from_resource(ResourceArn, TagKeyList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudhsm("RemoveTagsFromResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeyList"=>TagKeyList), args)); aws_config=aws_config)


text:
function build_ensemble_stochastic(
    A::Annealing,
    tf::Real,
    output_func,
    reduction;
    tspan=(0.0, tf),
    initializer=DEFAULT_INITIALIZER,
    save_positions=(false, false),
    kwargs...
)
    u0 = build_u0(A.u0, :v)
    initializer =
        initializer == DEFAULT_INITIALIZER ? (x, y) -> rand([-1, 1], x, y) :
        initializer

    flist = OpenQuantumBase.fluctuator_from_interactions(A.interactions)
    if length(flist) == 1
        cb = FluctuatorCallback(flist[1], initializer, save_positions)
    else
        cb = CallbackSet([FluctuatorCallback(f, initializer, save_positions) for f in flist]...)
    end
    ff = DiffEqLiouvillian(A.H, [], flist, size(A.H, 1))
    p = ODEParams(ff, float(tf), A.annealing_parameter)

    update_func = function (cache, u, p, t)
        update_cache!(cache, p.L, p, t)
    end

    cache = get_cache(A.H)
    diff_op = DiffEqArrayOperator(cache, update_func=update_func)
    jac_cache = similar(cache)
    jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
    ff = ODEFunction(diff_op, jac_prototype=jac_op)
    prob = ODEProblem{true}(ff, u0, tspan, p, callback=cb)

    ensemble_prob =
        EnsembleProblem(prob; output_func=output_func, reduction=reduction, kwargs...)
end


text:
"""
    pinmask(N::Integer) -> mask
Generates a `mask` that can be supplied to `likwid pin -s <mask>` to pin `N` Julia threads.

Taken from https://discourse.julialang.org/t/thread-affinitization-pinning-julia-threads-to-cores/58069/8.
"""
function pinmask(N::Integer)
    mask = UInt(0)
    for i in 1:N
        mask |= 1<<i
    end
    return _uint_to_hexmask(~mask) # Invert the mask to only pin Julia threads
end

_uint_to_hexmask(mask::UInt) = "0x" * string(mask, pad = sizeof(mask)<<1, base = 16)

"""
Set the verbosity level of the LIKWID library. Returns `true` on success.

Options are:
  * `LIKWID.LibLikwid.DEBUGLEV_ONLY_ERROR` or `0`

  * `LIKWID.LibLikwid.DEBUGLEV_INFO` or `1`

  * `LIKWID.LibLikwid.DEBUGLEV_DETAIL` or `2`

  * `LIKWID.LibLikwid.DEBUGLEV_DEVELOP` or `3`
"""
function setverbosity(verbosity::Integer)
    if verbosity ≥ LibLikwid.DEBUGLEV_ONLY_ERROR && verbosity ≤ LibLikwid.DEBUGLEV_DEVELOP
        LibLikwid.perfmon_setVerbosity(verbosity)
        return true
    end
    return false
end

"""
Returns whether LIKWID has been compiled with GPU support
(i.e. has been compiled with `NVIDIA_INTERFACE=true`).
"""
function gpusupport()
    # Note: Should be replaced by a less "hacky" implementation when
    # https://github.com/RRZE-HPC/likwid/issues/430 has been addressed.
    if isnothing(likwid_gpusupport[])
        likwid_gpusupport[] = _check_likwid_gpusupport()
    end
    return likwid_gpusupport[]
end

function _check_likwid_gpusupport()
    dlopen(liblikwid) do handle
        return !isnothing(dlsym(handle, :likwid_gpuMarkerInit; throw_error=false))
    end
end

function _check_likwid_gpusupport_alternative()
    x = _execute(`likwid-topology`)
    # if there is a gpu section in the output, return true
    contains(lowercase(x[1]), "gpu") && return true
    # if there is a cuda related error, return true
    contains(lowercase(x[2]), "cuda") && return true
    # else
    return false
end

"Run a Cmd object, returning the stdout & stderr contents plus the exit code"
function _execute(cmd::Cmd)
    out = Pipe()
    err = Pipe()

    process = run(pipeline(ignorestatus(cmd); stdout=out, stderr=err))
    close(out.in)
    close(err.in)

    out = (stdout=String(read(out)), stderr=String(read(err)), exitcode=process.exitcode)
    return out
end

function _execute_test(cmd::Cmd; print_only_on_fail=true)
    out = _execute(cmd)
    hasfailed = !iszero(out[:exitcode]) || (out[:stderr] != "")
    if !print_only_on_fail || hasfailed
        @info("Command: $cmd")
        @info("stdout:\n" * out[:stdout])
        @warn("stderr:\n" * out[:stderr])
    end
    return out[:exitcode] == 0
end

# working around https://github.com/JuliaLang/julia/issues/12711
function capture_stderr!(f::Function, io::IO)
    old_stderr = stderr
    rd, = redirect_stderr()
    task = @async write(io, rd)
    try
        ret = f()
        Libc.flush_cstdio()
        flush(stderr)
        return ret
    finally
        Base.close(rd)
        redirect_stderr(old_stderr)
        wait(task)
    end
end

function capture_stdout!(f::Function, io::IO)
    old_stdout = stdout
    rd, = redirect_stdout()
    task = @async write(io, rd)
    try
        ret = f()
        Libc.flush_cstdio()
        flush(stdout)
        return ret
    finally
        Base.close(rd)
        redirect_stdout(old_stdout)
        wait(task)
    end
end

function capture_stderr(f::Function)
    mktemp() do path, io
        redirect_stderr(io) do
            f()
        end
        flush(io)
        s = read(path, String)
    end
end

"""
List the values of `LIKWID_*` environment variables.
"""
env() = filter(x->startswith(x[1], "LIKWID"), ENV)

"""
Unset all `LIKWID_*` environment variables (for the current session).
"""
function clearenv()
    for k in keys(env())
        delete!(ENV, k)
    end
    return nothing
end

"""
Enables the overwriting of counters that are detected to be in-use.
The environment variable is similar to the `-f`/`--force` command line switch for `likwid-perfctr`.
"""
function LIKWID_FORCE end

"""
The execution does not require the access layer (access to hardware counters).
For example, this variable is set by `likwid-topology` or `likwid-pin`.
"""
function LIKWID_NO_ACCESS end

"""
The comma-separated list contains the CPUs the application threads should be pinned to.
Careful, the first CPU in the cpuset must be the last entry because the application is pinned to this CPU per default.
"""
function LIKWID_PIN end

"""
Disable stdout output caused by the library and the scripts.
Some scripts provide the `-q`/`--quiet` command line switch which provides the same functionality.
"""
function LIKWID_SILENT end

"""
Variable content must be a hexmask.
This hexmask describes which threads should be skipped while pinning.
This function is required to avoid pinning the shepherd threads used by some OpenMP and MPI implementations.
The version 4.3.1 introduced an automatic detection of the shepherd threads.
In most cases the detection works, but if not, the hexmask overwrites the automatic detection.
"""
function LIKWID_SKIP end

"""
Verbosity settings for the LIKWID library.
"""
function LIKWID_DEBUG end

"""
Always setup all counters in `setupCounters(...)` and don't respect the previous configuration.
Without this environment variable set, LIKWID writes the configuration to the register only if
the configuration has changed (compared to the last `setupCounters(...)` call)
"""
function LIKWID_FORCE_SETUP end

"""
LIKWID respects the CPUset of the calling process.
If you want to measure/run outside of this CPUset, use this environment variable.
It will not ignore the CPUset but create a new CPUset internally which contains `sysconf(_SC_NPROCESSORS_CONF)` hardware threads.
"""
function LIKWID_IGNORE_CPUSET end

"""
Filepath for the result file of the MarkerAPI.
"""
function LIKWID_FILEPATH end

"""
Access mode for MarkerAPI. `1` is the code for the access daemon.
"""
function LIKWID_MODE end

"""
Event string or performance group name.
Multiple event strings or performance group names can be separated by `|`.
"""
function LIKWID_EVENTS end

"""
The CPUs LIKWID is configured to run on (comma-separated list).
"""
function LIKWID_THREADS end

"""
Connection method for Intel MPI. Default is `ssh`, see option `-r` of `mpdboot` or similar.
"""
function LIKWID_MPI_CONNECT end

LIKWID_FORCE() = get(ENV, "LIKWID_FORCE", "")
LIKWID_NO_ACCESS() = get(ENV, "LIKWID_NO_ACCESS", "")
LIKWID_PIN() = get(ENV, "LIKWID_PIN", "")
LIKWID_SILENT() = get(ENV, "LIKWID_SILENT", "")
LIKWID_SKIP() = get(ENV, "LIKWID_SKIP", "")
LIKWID_DEBUG() = get(ENV, "LIKWID_DEBUG", "")
LIKWID_IGNORE_CPUSET() = get(ENV, "LIKWID_IGNORE_CPUSET", "")
LIKWID_FILEPATH() = get(ENV, "LIKWID_FILEPATH", "")
LIKWID_MODE() = get(ENV, "LIKWID_MODE", "")
LIKWID_EVENTS() = get(ENV, "LIKWID_EVENTS", "")
LIKWID_THREADS() = get(ENV, "LIKWID_THREADS", "")
LIKWID_MPI_CONNECT() = get(ENV, "LIKWID_MPI_CONNECT", "")

LIKWID_FORCE(v::Bool) = ENV["LIKWID_FORCE"] = Int(v);
LIKWID_NO_ACCESS(v::Bool) = ENV["LIKWID_NO_ACCESS"] = Int(v);
LIKWID_PIN(cpustr::AbstractString) = ENV["LIKWID_PIN"] = cpustr;
LIKWID_SILENT(v::Bool) = ENV["LIKWID_SILENT"] = Int(v);
LIKWID_SKIP(hexmask::AbstractString) = ENV["LIKWID_SKIP"] = hexmask;
LIKWID_DEBUG(v::Int) = ENV["LIKWID_DEBUG"] = v;
LIKWID_IGNORE_CPUSET(v::Bool) = ENV["LIKWID_IGNORE_CPUSET"] = Int(v);
LIKWID_FILEPATH(path::AbstractString) = ENV["LIKWID_FILEPATH"] = path;
LIKWID_MODE(mode) = ENV["LIKWID_MODE"] = mode;
LIKWID_EVENTS(eventstr::AbstractString) = ENV["LIKWID_EVENTS"] = eventstr;
LIKWID_THREADS(cpustr::AbstractString) = ENV["LIKWID_THREADS"] = cpustr;
LIKWID_MPI_CONNECT(x::AbstractString) = ENV["LIKWID_MPI_CONNECT"] = x;


text:
######## CmdStan diagnose example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoulli = "
  data {
    int<lower=0> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
      y ~ bernoulli(theta);
  }
  "

  bernoullidata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
  global stanmode, rc, diags, cnames
  stanmodel = Stanmodel(Diagnose(CmdStan.Gradient(epsilon=1e-6)),
    output_format=:array, name="bernoulli", model=bernoulli);

  rc, diags, cnames = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    println()
    display(diags)
    println()
    tmp = diags[:error][1]
    println("Test round.(diags[:error], digits=6) ≈ 0.0")
    @test round.(tmp, digits=6) ≈ 0.0
  end
end # cd


text:
using Plots
using StartUpDG

# Brusselator
B = 3
f(y1,y2) = [1+y1^2*y2-(B+1)*y1, B*y1-y1^2*y2]
f(Q) = f(Q[1],Q[2])
p,q = 1.5, 3.0
Q = [p;q]

dt = .1
FinalTime = 20

res = zero.(Q) # init RK residual
rk4a,rk4b,rk4c = ck45()
Nsteps = ceil(FinalTime/dt)
dt = FinalTime/Nsteps

plot() # init plot
for i = 1:Nsteps
    global res # yes, I know...this is just for simplicty
    for INTRK=1:5
        rhsQ = f(Q)
        @. res = rk4a[INTRK]*res + dt*rhsQ # i = RK stage
        @. Q =  Q + rk4b[INTRK]*res
    end
    scatter!([i*dt;i*dt],Q)
end
display(plot!(leg=false))


text:
function backward_rand_kernel(x, Gf, Gpred, Φ) where T
    xf, Pf = meancov(Gf)
    xpred, Ppred = meancov(Gpred)

    J = Pf*Φ'/Ppred
    xs = xf +  J*(x - xpred)
    Ps = Pf - J*Ppred*J'
    rand(Gaussian(xs, Ps))
end

struct BackwardRand{T,S}
    iter::T # iterate backwards through
    M::S
end

function dyniterate(B::BackwardRand, (x, state))
    ϕ = dyniterate(B.iter, state)
    ϕ === nothing && return nothing
    (n, U), state = ϕ
    Gf, Gpred = U[1], U[2]
    Gs = backward_rand_kernel(x, Gf, Gpred, Phi(RTS.M.sys, (n, U)))
    n=>x, (x, state)
end

function backward_sampler(M, prior, Y)
    P = filter(Y, M)

    ϕ = dyniterate(P, Start(prior))
    ϕ === nothing && error("no observations")
    (t, u), state = ϕ

    Xf = trajectory((t => (u[1], u[2]),))
    while true
        ϕ = dyniterate(P, state)
        ϕ === nothing && break
        (t, u), state = ϕ
        push!(Xf, t => (u[1], u[2])) # safe only what is needed
    end

    n = length(Xf)
    U = Xf.x[end]
    Gf, Gpred = U[1], U[2]

    x = rand(Gf)
    xs = fill(NaN*x, n)
    xs[n] = x
    for i in n-1:-1:1
        U = Xf.x[i]
        Gf, Gpredᵒ = U[1], U[2]
        x = backward_rand_kernel(x, Gf, Gpred, Phi(M.sys, i => U))
        xs[i] = x
        Gpred = Gpredᵒ
    end
    Trajectory(Xf.t, xs)
end


text:
using ABM
using Test
using DataStructures
using Distributed
using DistributedArrays

const testdir = dirname(@__FILE__)

#addprocs(4)

tests = [
    "priority",
    "schedule",
    "field2d",
    "agent",
    "real2d"
]

@testset "ABM" begin
    for t in tests
        tp = joinpath(testdir, "test-$(t).jl")
        include(tp)
    end
end


text:
"""
    initialize_random(target :: AbstractArray[, shape])

Initialize a system for the annealing procedure. The system is to be
filled with random data in such a way that the porosity of the
`target` is preserved.
"""
function initialize_random(target :: AbstractArray{T, N},
                           shape  :: NTuple{N, Int} = size(target)) where {T, N}
    # The porosity may be already calculated as $S^1_2(0)$ or
    # $L^1_2(0)$, but what if we do not use these functions for
    # reconstruction? Better to recalculate it again, it's fast.
    porosity = count(x -> x == 1, target) / length(target)

    # Only two-phase systems are supported by now. So to begin the
    # reconstruction we start with a system with all elements
    # belonging to the phase "zero".
    array = zeros(T, shape)

    # All indices in array
    indices = CartesianIndices(array)

    # Determine the number of elements with phase "one" to keep the
    # porosity of the original
    n = porosity * prod(shape)

    # Populate the array
    while n > 0
        idx = rand(indices)
        # If we have a free place, put phase "one" to it
        if array[idx] == 0
            array[idx] = 1
            n -= 1
        end
    end

    return array
end

"""
    initialize_spheres(target :: AbstractArray[, shape])

Initialize a system for the annealing procedure. The system is to be
filled with random spheres so its two-point correlation function for
the void phase is as close as possible to the `target`'s one.
"""
function initialize_spheres(target :: AbstractArray{T, N};
                            shape  :: NTuple{N, Int}  = size(target),
                            R0     :: Float64         = 10.0,
                            λ0     :: Float64         = 1e-3) where {T, N}
    R, λ = spheres_parameters(target, R0, λ0)
    return generate_spheres(shape, T, R, λ)
end


text:
using DynamicAxisWarping, DSP, LPVSpectral, Distances, SpectralDistances

# Create two signals
N = 48_000
g(x,N) = exp(-10*(x-N/2)^2/N^2)
t = 1:N
f = range(0.01, stop=1, length=N)
y = sin.(t .* f) .* g.(t, N)
y1 = [y; 0y] .+ 0.1 .* randn.()
y2 = [0y; y] .+ 0.1 .* randn.()

plot([y1 y2])
#
M1,M2 = melspectrogram.((y1,y2), 2048)
plot(plot(M1), plot(M2))
# Calculate the Dynamic Time-Warping cost between them (I do this in √ domain since that appears to produce the nicest looking plots)
@btime dtw(.√(M1.power), .√(M2.power))[1] evals=1 samples=5
# Calculating this distance is rather fast.
# Visualize the coupling matrix. The green line should cross diagonally (since the two signals are the same) starting halfway along the axes (since one signal is shifted by 50% of the length copared to the other).
dtwplot(.√(M1.power), .√(M2.power), linecolor=:green)
#

# The code above used the squared Euclidean distance between spectra for each time point. We can replace that distance with a transport-based distance instead.
# We thus combine dynamic time warping with a transport-based cost along the frequency axis
n,m = size(M1.power)
dist = DiscreteGridTransportDistance(Cityblock(), n, n)
dtw(.√(M1.power), .√(M2.power), dist)[1]
#
dtwplot(.√(M1.power), .√(M2.power), dist, linecolor=:green)
# We now see that the valley through which the path can move is slightly wider. This happens since there is now an additional degree of freedom due to transport along the frequency axis. Since there is no need for that transport in this case, the resulting path remains almost the same.
# We can create the need for some transport along the frequency axis

f1 = range(0.01, stop=1, length=N)
f2 = range(0.1, stop=1, length=N).^2
y10 = sin.(t .* f1) .* g.(t, N)
y20 = sin.(t .* f2) .* g.(t, N)
y1 = [y10; 0y10] .+ 0.1 .* randn.()
y2 = [0y20; y20] .+ 0.1 .* randn.()
M1,M2 = melspectrogram.((y1,y2), 2048)
plot(plot(M1), plot(M2))
# We now redo the same stuff as above
dtw(.√(M1.power), .√(M2.power))[1]
#
dtwplot(.√(M1.power), .√(M2.power), linecolor=:green)
#
dtw(.√(M1.power), .√(M2.power), dist)[1]
#
dtwplot(.√(M1.power), .√(M2.power), dist, linecolor=:green)
# We expect to see less time warping (straighter green line) now when it's allowed to fudge the distance a bit by shuffling mass along the frequency axis.
# How long time does it take to compute such a distance? (on a laptop, expect about 3x improvement on desktop)
@btime dtw($(M1.power), $(M2.power), $dist)
@btime dtw($(M1.power), $(M2.power), radiuslimits(50,m,m)..., $dist)
# The size of the spectrogram
size(M1.power)
# The complexity is almost linear in the number of frequency bins, but quadratic in the number of time steps.


text:
using Documenter, BitTorrent

if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

makedocs(
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "GITHUB_ACTIONS"),
        canonical = "https://fredrikekre.github.io/BitTorrent.jl/v1",
    ),
    modules = [BitTorrent],
    sitename = "BitTorrent.jl",
    pages = Any[
        "index.md",
    ]
)

deploydocs(
    repo = "github.com/fredrikekre/BitTorrent.jl.git",
    push_preview=true,
)