text:
# Alternative shebang:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#
using ArgParse
using Distributed
using Printf

VERSION = "1.0.0"

function parse_commandline()
    s = ArgParseSettings()
    s.epilog = "MentaLiST -- The MLST pipeline developed by the PathOGiST research group. https://github.com/WGS-TB/MentaLiST\n" *
      "To cite: Feijao P, Yao H, Fornika D, Gardy J, Hsiao W, Chauve C, Chindelevitch L. 10/01/2018. Microbial Genomics 4(2): doi:10.1099/mgen.0.000146\n"
    s.version = "MentaLiST $VERSION."
    s.allow_ambiguous_opts= true # to allow -1 and -2.
    s.preformatted_epilog = true
    @add_arg_table! s begin
      "call"
        help = "MLST caller, given a sample and a k-mer database."
        action = :command
      "build_db"
        help = "Build a MLST k-mer database, given a list of FASTA files."
        action = :command
      "db_info"
        help = "Extract information from an existing MentaLiST k-mer database"
        action = :command
      "list_pubmlst"
        help = "List all available MLST schemes from www.pubmlst.org. "
        action = :command
      "download_pubmlst"
        help = "Dowload a MLST scheme from pubmlst and build a MLST k-mer database."
        action = :command
      "list_cgmlst"
        help = "List all available cgMLST schemes from www.cgmlst.org."
        action = :command
      "download_cgmlst"
          help = "Dowload a MLST scheme from cgmlst.org and build a MLST k-mer database."
          action = :command
      "download_enterobase"
          help = "Dowload a MLST scheme from Enterobase (enterobase.warwick.ac.uk) and build a MLST k-mer database."
          action = :command
      "-v", "--version"
          action = :show_version
          help = "show version information and exit"
      end

    # Calling MLST options:
    s["call"].preformatted_epilog = true
    s["call"].epilog = "MentaLiST MLST calling function. Calls alleles on a given MLST database.\n" *
      "You can create a custom DB with 'create_db' or other MentaLiST functions that download schemes from pubmlst, cgmlst.org or Enterobase.\n\n" *
      "Examples:\n" *
      "mentalist call -o my_sample.mlst --db my_scheme.db -1 sample_1.fastq.gz -2 sample_2.fastq.gz # one paired-end sample.\n" *
      "mentalist call -o all_samples.mlst --db my_scheme.db -1 *.fastq.gz -2 *.fastq.gz # multiple paired-end samples."

    @add_arg_table! s["call"] begin
        "-o"
          help = "Output file with MLST call"
          arg_type = String
          required = true
        "--db"
          help = "Kmer database"
          required = true
          arg_type = String
        "-t", "--mutation_threshold"
          help = "Maximum number of mutations when looking for novel alleles."
          arg_type = Int
          default = 6
        "--kt"
          help = "Minimum # of times a kmer is seen to be considered present in the sample (solid)."
          arg_type = Int
          default = 10
        "--output_votes"
          help = "Outputs the results for the original voting algorithm."
          action = :store_true
        "--output_special"
          help = "Outputs a FASTA file with the alleles from 'special cases' such as incomplete coverage, novel, and multiple alleles. "
          action = :store_true
        "-i", "--sample_input_file"
          help = "Input TXT file for multiple samples. First column has the sample name, second the FASTQ file. Repeat the sample name for samples with more than one file (paired reads, f.i.)"
        "-1"
          nargs = '*'
          help = "FastQ input files, one per sample, forward reads (or unpaired reads)."
          arg_type = String
        "-2"
          nargs = '*'
          help = "FastQ input files, one per sample, reverse reads."
          arg_type = String
        "--fasta"
          help = "Input files are in FASTA format, instead of the default FASTQs."
          action = :store_true
    end

    ## Common option for all db building commands:
    s_db = ArgParseSettings()
    @add_arg_table! s_db begin
      "--db"
        help = "Output file (kmer database)"
        arg_type = String
        required = true
      "-k"
        help = "Kmer size"
        required = true
        arg_type = Int8
      "--threads" # disabled for julia 1.1, the current syntax does not seem to work.
        arg_type = Int
        default = 1
        help = "Number of threads used in parallel."
      "-c", "--allele_coverage"
       arg_type = Float64
       default = 1.0
       help = "Minimum percentage of allele coverage in number of kmers of each allele."
    end
    # Build DB from FASTA, options:
    import_settings!(s["build_db"], s_db)
    @add_arg_table! s["build_db"] begin
        "-d", "--database"
            arg_type = String
            help = "MLST Fasta Database Directory"
            required = true
        "-p", "--profile"
            arg_type = String
            help = "Profile file for known genotypes."
    end

    @add_arg_table! s["db_info"] begin
      "--db"
        help = "MentaLiST kmer database"
        arg_type = String
        required = true
    end

    # Listing functions, common options:
    s_list = ArgParseSettings()
    @add_arg_table! s_list begin
      "-p", "--prefix"
      help = "Only list schemes where the species name starts with this prefix."
      arg_type = String
    end

    # List pubmlst
    import_settings!(s["list_pubmlst"], s_list)

    # List cgmlst
    import_settings!(s["list_cgmlst"], s_list)

    # Common options for MLST download functions:
    s_db_download = ArgParseSettings()
    import_settings!(s_db_download, s_db) # import build_db common options
    @add_arg_table! s_db_download begin
      "-o", "--output"
        help = "Output folder for the scheme Fasta files."
        arg_type = String
        required = true
      "-s", "--scheme"
        help = "Species name or scheme ID."
        arg_type = String
        required = true
    end

    # Download pubmlst:
    import_settings!(s["download_pubmlst"], s_db_download) # import common options

    # Download cgmlst:
    import_settings!(s["download_cgmlst"], s_db_download)  # import common options

    # Download enterobase:
    import_settings!(s["download_enterobase"], s_db_download) # import common options
    # add (and override) specific enterobase options:
    s["download_enterobase"].error_on_conflict = false  # do not error-out when trying to override an option
    @add_arg_table! s["download_enterobase"] begin
      "-s", "--scheme"
        help = "Letter identifying which scheme: (S)almonella, (Y)ersinia, or (E)scherichia/Shigella."
        arg_type = String
        required = true
      "-t", "--type"
        help = "Choose the type: 'cg' or 'wg' for cgMLST or wgMLST scheme, respectively."
        arg_type = String
        required = true
    end

    return parse_args(s)
end

# Check if files exist:
function check_files(files)
  dont_exist = [file for file in files if !isfile(file)]
  if length(dont_exist) > 0
    exit_error("The following input file(s) could not be found: $(join(dont_exist,',')).")
  end
end

### Error FUNCTIONS
@everywhere function exit_error(msg)
  @error(msg)
  println("exiting ...")
  exit(-1)
end

#### Main COMMAND functions:
function call_mlst(args)
  run_calling_pipeline(args)
end

function list_pubmlst(args)
  list_pubmlst_schema(args["prefix"])
end

function download_pubmlst(args)
  loci_files, profile_file = download_pubmlst_scheme(args["scheme"], args["output"])
  @info "Building the k-mer database ..."
  args["database"] = dirname(loci_files[1])
  args["profile"] = profile_file
  build_db(args)
end

function list_cgmlst(args)
  list_cgmlst_schema(args["prefix"])
end

function download_cgmlst(args)
  loci_files = download_cgmlst_scheme(args["scheme"], args["output"])
  @info "Building the k-mer database ..."
  args["database"] = dirname(loci_files[1])
  args["profile"] = nothing
  build_db(args)
end

function download_enterobase(args)
  loci_files = download_enterobase_scheme(args["scheme"], args["type"], args["output"])
  @info "Building the k-mer database ..."
  args["database"] = dirname(loci_files[1])
  args["profile"] = nothing
  build_db(args)
end

function build_db(args, version=VERSION)

  # check if we need Gurobi:
  if args["allele_coverage"] < 1
    ok = include_gurobi()
    if !ok
      exit_error("Could not load libraries. To use the allele coverage functionality (option -c), a proper installation of Gurobi (www.gurobi.com) and the Gurobi and JuMP julia packages are required.\nCheck https://github.com/JuliaOpt/Gurobi.jl for installation instructions.")
    end
  end

  # get files
  fasta_files = readdir(args["database"])

  # get abs path
  for i in 1:length(fasta_files)
    fasta_files[i]=joinpath(args["database"], fasta_files[i])
  end

  # check if files exist:
  check_files(fasta_files)
  # get arguments and call the kmer db builder for each locus:
  k::Int8 = args["k"]
  db_file = args["db"]
  profile = args["profile"]
  @info "Opening FASTA files ... "
  results, loci = kmer_class_for_each_locus(DNAKmer{k}, fasta_files, args["allele_coverage"])
  # Combine results:
  @info "Combining results for each locus ..."
  kmer_classification = combine_loci_classification(DNAKmer{k}, results, loci)

  for (index, fasta_file) in enumerate(fasta_files)
    scheme_fasta_directory = pop!(split(dirname(fasta_file), "/"))
    if scheme_fasta_directory == ""
      scheme_fasta_directory = "."
    end
    fasta_files[index] = scheme_fasta_directory * "/" * basename(fasta_file)
  end
  @info "Saving DB ..."
  save_db(DNAKmer{k}, kmer_classification, loci, db_file, profile, args, version)
  @info "Done!"
end

function db_info(args)
  filename = args["db"]
  @info "Opening kmer database ... "
  # Compressed database, open and decompress/decode in memory:
  d = load("$filename")
  @info "Finished the JLD load."
  build_args = JSON.parse(d["args"])
  k = build_args["k"]
  loci = d["loci"]
  loci_list = Blosc.decompress(Int32, d["loci_list"])
  num_loci = length(loci_list)
  mentalist_version = try
    d["mentalist_version"]
  catch
    "unknown"
  end
  println("mentalist_version\t$mentalist_version")
  println("k\t$k")
  println("num_loci\t$num_loci")
end

## Gurobi ILP:

# Kmer coverage is based on Gurobi:
function has_gurobi()
  packages = keys(Pkg.installed())
  return (in("Gurobi", packages) & in("JuMP", packages))
end

function include_gurobi()
  if has_gurobi()
    try
      include("kmer_coverage.jl")
    catch e
      if isa(e, LoadError)
        @error("Error trying to load Gurobi package.")
        return false
      end
    end
    return true
  else
    return false
  end
end

##### Main function: just calls the appropriate commands, with arguments:

args = parse_commandline()
# determine command:
cmd = args["%COMMAND%"]
if cmd == "call"
  include("calling_functions.jl")
  call_mlst(args[cmd])

elseif cmd == "build_db"
  if args[cmd]["threads"] > 1
    addprocs(args[cmd]["threads"])
  end
  include("build_db_functions.jl")
  build_db(args[cmd])

elseif cmd == "db_info"
  import JLD: load, JSON, Blosc
  db_info(args[cmd])

elseif cmd == "list_pubmlst"
  include("mlst_download_functions.jl")
  list_pubmlst(args[cmd])

elseif cmd == "download_pubmlst"
  include("mlst_download_functions.jl")
  addprocs(args[cmd]["threads"])
  include("build_db_functions.jl")
  download_pubmlst(args[cmd])

elseif cmd == "list_cgmlst"
  include("mlst_download_functions.jl")
  list_cgmlst(args[cmd])

elseif cmd == "download_cgmlst"
  if args[cmd]["threads"] > 1
    addprocs(args[cmd]["threads"])
  end
  include("mlst_download_functions.jl")
  include("build_db_functions.jl")
  download_cgmlst(args[cmd])

elseif cmd == "download_enterobase"
  if args[cmd]["threads"] > 1
    addprocs(args[cmd]["threads"])
  end
  include("mlst_download_functions.jl")
  include("build_db_functions.jl")
  download_enterobase(args[cmd])

end

max_stars_repo_name: CNRResistanceAntibiotic/MentaLiST

max_issues_repo_name: CNRResistanceAntibiotic/MentaLiST

max_forks_repo_name: CNRResistanceAntibiotic/MentaLiST


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Atm data file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds

    global mask = ds["mask"][:] |> mreplace
    global lat = ds["yc"][1, :] |> mreplace

    global lat_weight = 2π * Re^2.0 * cos.(deg2rad.(lat))
    global ϕ          = deg2rad.(lat)

    global L_idx = (mask .== 1)
    global O_idx = (mask .== 0)

end

Dataset(parsed["data-file"], "r") do ds

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12
    rng = (:,:,beg_t:end_t)

    global A_flx_2D = (
        FFLX_TOA = ( - mean(ds["FSNT"][rng...], dims=(3, ))  + mean(ds["FLNT"][rng...],  dims=(3, )) )[:, :, 1],
        FFLX_SFC = ( - mean(ds["FSNS"][rng...], dims=(3, ))  + mean(ds["FLNS"][rng...],  dims=(3, )) )[:, :, 1],
        HFLX_SFC = (   mean(ds["SHFLX"][rng...], dims=(3, )) + mean(ds["LHFLX"][rng...], dims=(3, )) )[:, :, 1],
    )

    global (Nx, Ny) = size(A_flx_2D.FFLX_TOA)

    O_flx_2D = (
        FFLX_SFC = zeros(Nx, Ny),
        HFLX_SFC = zeros(Nx, Ny),
    )

    L_flx_2D = (
        FFLX_SFC = zeros(Nx, Ny),
        HFLX_SFC = zeros(Nx, Ny),
    )

    O_flx_2D.FFLX_SFC[O_idx] .= A_flx_2D.FFLX_SFC[O_idx]
    O_flx_2D.HFLX_SFC[O_idx] .= A_flx_2D.HFLX_SFC[O_idx]

    L_flx_2D.FFLX_SFC[L_idx] .= A_flx_2D.FFLX_SFC[L_idx]
    L_flx_2D.HFLX_SFC[L_idx] .= A_flx_2D.HFLX_SFC[L_idx]

    # zonal average
    global A_flx = (
        FFLX_TOA = mean(A_flx_2D.FFLX_TOA, dims=(1,))[1,:] .* lat_weight,
        FFLX_SFC = mean(A_flx_2D.FFLX_SFC, dims=(1,))[1,:] .* lat_weight,
        HFLX_SFC = mean(A_flx_2D.HFLX_SFC, dims=(1,))[1,:] .* lat_weight,
    )

    global O_flx = (
        FFLX_SFC = mean(O_flx_2D.FFLX_SFC, dims=(1,))[1,:] .* lat_weight,
        HFLX_SFC = mean(O_flx_2D.HFLX_SFC, dims=(1,))[1,:] .* lat_weight,
    )

    global L_flx = (
        FFLX_SFC = mean(L_flx_2D.FFLX_SFC, dims=(1,))[1,:] .* lat_weight,
        HFLX_SFC = mean(L_flx_2D.HFLX_SFC, dims=(1,))[1,:] .* lat_weight,
    )

end

A_EFLX_CONV = - ( A_flx.FFLX_TOA - A_flx.FFLX_SFC - A_flx.HFLX_SFC )
O_EFLX_CONV = - ( O_flx.FFLX_SFC + O_flx.HFLX_SFC )
L_EFLX_CONV = - ( L_flx.FFLX_SFC + L_flx.HFLX_SFC )

A_IET = integrate(ϕ, A_EFLX_CONV)
O_IET = integrate(ϕ, O_EFLX_CONV)
L_IET = integrate(ϕ, L_EFLX_CONV)

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Ny", Ny)

    for (varname, vardata, vardim, attrib) in [
        ("A_FFLX_TOA",  A_flx.FFLX_TOA,  ("Ny",), Dict()),
        ("A_FFLX_SFC",  A_flx.FFLX_SFC,  ("Ny",), Dict()),
        ("A_HFLX_SFC",  A_flx.HFLX_SFC,  ("Ny",), Dict()),

        ("O_FFLX_SFC",  O_flx.FFLX_SFC,  ("Ny",), Dict()),
        ("O_HFLX_SFC",  O_flx.HFLX_SFC,  ("Ny",), Dict()),

        ("L_FFLX_SFC",  L_flx.FFLX_SFC,  ("Ny",), Dict()),
        ("L_HFLX_SFC",  L_flx.HFLX_SFC,  ("Ny",), Dict()),

        ("A_IET",       A_IET,           ("Ny",), Dict()),
        ("O_IET",       O_IET,           ("Ny",), Dict()),
        ("L_IET",       L_IET,           ("Ny",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

ext: jl

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
## VISIBLE SECTION
#

function flock_basic!(B; tau=0.1)
  B0 = deepcopy(B)
  @simd for bi in B
    deli = neighbors(bi, B0)

    si = seperate(bi, deli)
    ai = align(bi, deli)
    ci = cohere(bi, deli)

    bi.vel += bi.ws*si + bi.wa*ai + bi.wc*ci
    bi.pos += tau*bi.vel
  end
end

function neighbors(bi, B)
  deli = [bj for bj in B if
    bi != bj && dis(bi, bj) < bi.dv && ang(bi, bj) < bi.av]
  return deli
end

function seperate(bi, deli)
  si = zeros(dim(bi))
  for bj in deli
    if dis(bi, bj) < bi.ds
      si -= bj.pos - bi.pos
    end
  end
  return si
end

function align(bi, deli)
  ai = zeros(dim(bi))
  if isempty(deli)
    return ai
  end
  for bj in deli
    ai += bj.vel
  end
  return ai/length(deli) - bi.vel
end

function cohere(bi, deli)
  ci = zeros(dim(bi))
  if isempty(deli)
    return ci
  end
  for bj in deli
    ci += bj.pos
  end
  return ci/length(deli) - bi.pos
end

## HIDDEN SECTION
#

function initialize(;
  dim=2,
  nb=100, dv=5., av=3.05, ds=dv/2., nc=4, de=3*dv, dp=3*dv,
  w0=1., ws=.1, wa=.05, wc=.01, ww=0., wb=.1, wp=1., we=1.,
  vmax=10., vmin=.4*vmax,
  bmin=-(200/dim^2)*ones(dim), bmax=(200/dim^2)*ones(dim)
)
  Random.seed!(0)
  pfun = x -> bmin + x .* (bmax .- bmin)
  vfun = x -> (x ./ norm(x)) .* vmax / 2
  B = [Boid(i, pfun(rand(dim)), vfun(rand(dim)*2 - ones(dim)),
    dv, av, ds, w0, ws, wa, wc, vmax, vmin, bmin, bmax, ww, wb, wp, we, nc, dp, de) for i = 1:nb]
  W = World()
  if dim == 2
    W.xy = Node(collect(Point2f0(b.pos) for b in B))
  else
    W.xy = Node(collect(Point3f0(b.pos) for b in B))
  end
  W.sc = scatter(W.xy)
  return (B, W)
end

function draw!(W, B; t, dim=2, fps=60)
  if dim == 2
    W.xy[] = collect(Point2f0(b.pos) for b in B)
  else
    W.xy[] = collect(Point3f0(b.pos) for b in B)
  end
  dt = time() - t
  if dt > 1/fps
    sleep(1/fps)
    return time()
  else
    sleep(1/fps - dt)
    return time()
  end
  #sleep(1/fps)
  #return time()
end

function stat!(W, B)
  # TODO: make relevant draw function
  n = dim(first(B))
  pgv = zeros(n)
  cg = zeros(n)
  for bi in B
    pgv += bi.vel/norm(bi.vel)
    cg += bi.pos
  end
  mgv = zeros(n, n)
  for bi in B
    mgv += (bi.pos - cg)*bi.vel'
  end
  push!(W.sc, norm(pgv)/n, norm(mgv)/n)
  for bi in B
  end
end

function run_basic(; nb=100, nmax=1000, dim=2)
  B, W = initialize(nb=nb, dim=dim)
  t = time()
  display(W.sc)

  try
    @fastmath for _ = 1:nmax
      t = draw!(W, B, t=t, dim=dim)
      flock_basic!(B)
    end
  catch InterruptException
    println("\nSimulation interrupted.")
  end
end


text:
function _RheoTimeData_explicit_nolog()
    a = Vector{RheoFloat}(1.0:1.0:3.0)
    b = Vector{RheoFloat}(4.0:1.0:6.0)
    c = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(a, b, c, nothing)

    data.σ==a && data.ϵ==b && data.t==c && isnothing(data.log)
end

function _operators_logs()
    d1=strainfunction(timeline(),t->exp(-t))
    d2=strainfunction(timeline(),t->1-exp(-t))
    d=2*d1 - (-d2) + d2

    d3=rheologrun(d.log)

    (d3.ϵ == d.ϵ) && all([ abs(e-2.)<=eps(RheoFloat) for e in d.ϵ ])
end
@test _operators_logs()

function _union_strain1stress2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_strain|imposed_stress

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_strain1stress2()

function _union_stress1strain2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_stress|imposed_strain

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_stress1strain2()

function _freeze_params()
    SLS2_mod = freeze_params( SLS2, (G₀=2,η₂=3.5))

    SLS2.G(1,[2,1,2,3,3.5])≈SLS2_mod.G(1,[1,2,3])
end
@test _freeze_params()


text:
"""
    `get_rew_rate(rewards, alpha, initialvalue)`

Use accumulate to update reward rate according to
a starting point 'initialvalue'
a sequence of outcomes 'rewards'
and a learning rate 'alpha'
"""

function get_rew_rate(rewards, alpha, initialvalue)
   return accumulate(rewards; init=initialvalue) do value, reward
       #=(1 - alpha) * value + alpha*reward=#
       value + alpha*(reward - value)
   end
end

"""
    `get_consecutive_omissions_rate(n, alpha, initialvalue)`

Calculate rew_rate update for n consecutive no rewards
"""
function get_consecutive_omissions_rate(n, alpha, initialvalue)
    initialvalue*(1-alpha)^n
end

"""
    `get_consecutive_rewards_rate(n, alpha, initialvalue)`

Calculate rew_rate update for n consecutive no rewards
"""

function get_consecutive_rewards_rate(n, alpha, initialvalue)
    initialvalue*((1-alpha)^n) + alpha*(ploynomial_alpha(1-alpha,n))
end

"""
    `ploynomial_alpha(n, alpha, initialvalue)`

    Helper function to calculate the value of the polynomial necessary to calculate consecutive rewards rate
"""
function ploynomial_alpha(oneminusalpha,n)
    n<1 && error("minimum one outcome to consider")
    sum(oneminusalpha.^(0:n-1))
end
"""
    `get_1bin_reward_rate(n, alpha, initialvalue)`
"""
function get_1bin_reward_rate(n, alpha, initialvalue)
    pre_rew_rate = initialvalue*(1-alpha)^(n-1) #omission until the last bin
    pre_rew_rate + alpha*(1-pre_rew_rate)
end

function get_rew_rate(rewards::AbstractVector,n::AbstractVector, initialvalue::Float64, alpha::Float64)
    x = collect(zip(rewards,n))
    return accumulate(x; init=initialvalue) do value, x
        #= x is a tuple containing
        the outcome val (1 = rew, 0 = omission): x[1]
        and the number of bins in that poke: x[2]
        =#
        # Bool(x[1]) ? get_consecutive_rewards_rate(x[2], alpha, value) : get_consecutive_omissions_rate(x[2], alpha, value)
        Bool(x[1]) ? get_1bin_reward_rate(x[2], alpha, value) : get_consecutive_omissions_rate(x[2], alpha, value)
    end
end

function env_get_rew_rate(rewards::AbstractVector,n::AbstractVector, initialvalue::AbstractVector, reset::AbstractVector, alpha)
    x = collect(zip(reset, initialvalue,rewards,n))
    #= x is a tuple containing
    x[1]: whether the reward rate has to be reset i.e. first poke in a new session for env_rew
    x[2]: the initial env_reward rate value (precomputed per session)
    x[3]: the outcome val (1 = rew, 0 = omission)
    x[4]: and the number of bins in that poke
    =#
    return accumulate(x; init=x[2]) do value, x
        if Bool(x[1])
            # reset env_rew or patch_rew if a new session or trial starts respectively
            return Bool(x[3]) ? get_1bin_reward_rate(x[4], alpha, x[2]) : get_consecutive_omissions_rate(x[4], alpha, x[2])
        else
            return Bool(x[3]) ? get_1bin_reward_rate(x[4], alpha, value) : get_consecutive_omissions_rate(x[4], alpha, value)
        end
    end
end

function patch_get_rew_rate(rewards::AbstractVector,n::AbstractVector, resetvalue::Float64, resetcases::AbstractVector, alpha::Float64)
    x = collect(zip(resetcases,rewards,n))
    #= x is a tuple containing
    x[1]: whether the reward rate has to be reset i.e. first poke in a new trial for patch_rew
    x[2]: the outcome val (1 = rew, 0 = omission)
    x[3]: and the number of bins in that poke
    =#
    return accumulate(x; init=resetvalue) do value, x
        if Bool(x[1])
            # reset env_rew or patch_rew if a new session or trial starts respectively
            return Bool(x[2]) ? get_1bin_reward_rate(x[3], alpha, resetvalue) : get_consecutive_omissions_rate(x[3], alpha, resetvalue)
        else
            return Bool(x[2]) ? get_1bin_reward_rate(x[3], alpha, value) : get_consecutive_omissions_rate(x[3], alpha, value)
        end
    end
end

# rew = [0,0,1,1,0,0,0,1,1,1]
# alpha_val = 0.1
# initialvalue = 0.5
# get_rew_rate(rew, alpha_val, initialvalue)
# get_rew_rate([0,1,0,1],[2,2,3,3],alpha_val,initialvalue)


text:
using TinyAutoGrad
using Test

@testset "TinyAutoGrad.jl" begin
    @testset "Simple Calculations" begin
        x = Dual(3)
        y = Dual(6)

        @test x^2 == Dual(9, 6) # Because dx^2/dx = 2x
        @test y^2 == Dual(36, 12)
    end

    @testset "Derivative" begin
        f(x) = x^2
        g(x) = 6x^3 + 2x^2
        df(x) = derive(f, x)
        dg(x) = derive(g, x)

        # Manually derive g(x) for comparison
        h(x) = 18x^2 + 4x

        @test df(3) == 2*3 # Because f'(x) = 2x
        @test dg(1) == h(1)
        @test dg(2) == h(2)
    end

    @testset "Gradients" begin
        f(x, y) = 3x + 10y
        g(x, y) = 3x^2 + 10y^2

        @test gradient(f, 1, 1) == [3, 10]
        @test gradient(f, 2, 3) == [3, 10]

        @test gradient(g, 1, 1) == [6, 20]
        @test gradient(g, 2, 3) == [12, 60]
    end

    @testset "Gradient of Params" begin
        f(x, y) = 3x + 10y

        ps = Params()

        g() = 3ps.x + 10ps.y

        ∇f = gradient(f, 1, 1)
        ps.x = 1
        ps.y = 1
        ∇g = gradient(g, ps)

        @test ∇f[1] == ∇g.x
        @test ∇f[2] == ∇g.y

        ∇f = gradient(f, 2, 3)
        ps.x = 2
        ps.y = 3
        ∇g = gradient(g, ps)

        @test ∇f[1] == ∇g.x
        @test ∇f[2] == ∇g.y
    end

end


text:
function restart(backcolor)
    empty_figure(filename = joinpath(@__DIR__, "test_24.png"))
    background(backcolor)
    sethue(PALETTE[8])
end
"""
    prline_1(vx, vy)

This is not for use in the final algorithm, just for debugging.
Plots a polyline until NaN or end of vectors.
"""
function prline_1(origo, vx, vy)
    move(origo + (vx[1], vy[1]))
    for (x, y) in zip(vx, vy)
        isnan(x) && break
        isnan(y) && break
        line(origo + (x, y))
    end
    do_action(:stroke)
end

"""
    euler_step_complex_1!(f, vx, vy, h, n)

Coordinates (vx[n + 1], vy[n + 1]) are updated with the estimated position,
along function f.
     f is a complex-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size factor
"""
function euler_step_complex_1!(f, vx, vy, h, n)
    @assert length(vx) == length(vy)
    @assert n < length(vx)
    f0 = f(vx[n], vy[n])
    vx[n + 1] = vx[n] + real(f0) * h
    vy[n + 1] = vy[n] + imag(f0) * h
end
"""
    euler_forwardsteps_1!(f, vx, vy, h)

Coordinates (vx[n > 1], vy[n > 1]) are updated with the estimated positions,
along function f.
     f is a complex-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size factor
"""
function euler_forwardsteps_1!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        euler_step_complex_1!(f, vx, vy, h, n)
    end
end

"""
    rk4_step_complex_1!(f, vx, vy, h, n)

Coordinates (vx[n + 1], vy[n + 1]) are updated with the estimated position,
along function f.
     f is a complex-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size quantity. If the integration variable is time,
     it defines the length of time betweeen points.
"""
function rk4_step_complex_1!(f, vx, vy, h, n)
    @assert length(vx) == length(vy)
    @assert n < length(vx)
    f0 = f(vx[n], vy[n])
    x1 = vx[n] + real(f0) * h * 0.5
    y1 = vy[n] + imag(f0) * h * 0.5
    f1 = f(x1, y1)
    x2 = vx[n] + real(f1) * h * 0.5
    y2 = vy[n] + imag(f1) * h * 0.5
    f2 = f(x2, y2)
    x3 = vx[n] + real(f2) * h
    y3 = vy[n] + imag(f2) * h
    f3 = f(x3, y3)
    vx[n + 1] = vx[n] + 1/6 * h * real( f0  + 2∙f1 + 2∙f2 + f3 )
    vy[n + 1] = vy[n] + 1/6 * h * imag( f0  + 2∙f1 + 2∙f2 + f3 )
end

"""
    rrk4_forwardsteps_1!(f, vx, vy, h)

Coordinates (vx[n > 1], vy[n > 1]) are updated with the estimated positions,
along function f.
     f(x,y) returns a complex-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size quantity. If the integration variable is time,
     it defines the length of time betweeen points.
"""
function rrk4_forwardsteps_1!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        rk4_step_complex_1!(f, vx, vy, h, n)
    end
end
@inline function euler_step_complex_2!(f, vx, vy, h, n)
    @assert length(vx) == length(vy)
    @assert n < length(vx)
    f0 = f(vx[n], vy[n])
    vx[n + 1] = vx[n] + real(f0) * h
    vy[n + 1] = vy[n] + imag(f0) * h
end
@inline function rk4_step_complex_2!(f, vx, vy, h, n)
    @assert length(vx) == length(vy)
    @assert n < length(vx)
    f0 = f(vx[n], vy[n])
    x1 = vx[n] + real(f0) * h * 0.5
    y1 = vy[n] + imag(f0) * h * 0.5
    f1 = f(x1, y1)
    x2 = vx[n] + real(f1) * h * 0.5
    y2 = vy[n] + imag(f1) * h * 0.5
    f2 = f(x2, y2)
    x3 = vx[n] + real(f2) * h
    y3 = vy[n] + imag(f2) * h
    f3 = f(x3, y3)
    vx[n + 1] = vx[n] + 1/6 * h * real( f0  + 2∙f1 + 2∙f2 + f3 )
    vy[n + 1] = vy[n] + 1/6 * h * imag( f0  + 2∙f1 + 2∙f2 + f3 )
end
function euler_forwardsteps_2!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        euler_step_complex_2!(f, vx, vy, h, n)
    end
end
function rrk4_forwardsteps_2!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        rk4_step_complex_2!(f, vx, vy, h, n)
    end
end
# Inconclusive inlining for both . Redefine the inner functions without assertions

function euler_step_complex_3!(f, vx, vy, h, n)
    f0 = f(vx[n], vy[n])
    vx[n + 1] = vx[n] + real(f0) * h
    vy[n + 1] = vy[n] + imag(f0) * h
end
function rk4_step_complex_3!(f, vx, vy, h, n)
    f0 = f(vx[n], vy[n])
    x1 = vx[n] + real(f0) * h * 0.5
    y1 = vy[n] + imag(f0) * h * 0.5
    f1 = f(x1, y1)
    x2 = vx[n] + real(f1) * h * 0.5
    y2 = vy[n] + imag(f1) * h * 0.5
    f2 = f(x2, y2)
    x3 = vx[n] + real(f2) * h
    y3 = vy[n] + imag(f2) * h
    f3 = f(x3, y3)
    vx[n + 1] = vx[n] + 1/6 * h * real( f0  + 2∙f1 + 2∙f2 + f3 )
    vy[n + 1] = vy[n] + 1/6 * h * imag( f0  + 2∙f1 + 2∙f2 + f3 )
end
function euler_forwardsteps_3!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        euler_step_complex_3!(f, vx, vy, h, n)
    end
end
function rrk4_forwardsteps_3!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    for n in 1:(length(vx) - 1)
        rk4_step_complex_3!(f, vx, vy, h, n)
    end
end
# Removing inner assertions was infavourable for Euler, favourable for RK4.

# Redefine the outer functions too without assertions

function euler_forwardsteps_4!(f, vx, vy, h)
    for n in 1:(length(vx) - 1)
        euler_step_complex_3!(f, vx, vy, h, n)
    end
end
function rrk4_forwardsteps_4!(f, vx, vy, h)
    for n in 1:(length(vx) - 1)
        rk4_step_complex_3!(f, vx, vy, h, n)
    end
end

# Removing outer assertions was favourable for Euler, not for RK4.

# Proceed with the optimal versions so far.
# Use @fastmath

@fastmath function euler_step_complex_4!(f, vx, vy, h, n)
    @assert length(vx) == length(vy)
    @assert n < length(vx)
    @fastmath f0 = f(vx[n], vy[n])
    @fastmath vx[n + 1] = vx[n] + real(f0) * h
    @fastmath vy[n + 1] = vy[n] + imag(f0) * h
end
@fastmath function rk4_step_complex_4!(f, vx, vy, h, n)
    @fastmath f0 = f(vx[n], vy[n])
    @fastmath x1 = vx[n] + real(f0) * h * 0.5
    @fastmath y1 = vy[n] + imag(f0) * h * 0.5
    @fastmath f1 = f(x1, y1)
    @fastmath  x2 = vx[n] + real(f1) * h * 0.5
    @fastmath y2 = vy[n] + imag(f1) * h * 0.5
    @fastmath f2 = f(x2, y2)
    @fastmath x3 = vx[n] + real(f2) * h
    @fastmath y3 = vy[n] + imag(f2) * h
    @fastmath f3 = f(x3, y3)
    @fastmath vx[n + 1] = vx[n] + 1/6 * h * real( f0  + 2∙f1 + 2∙f2 + f3 )
    @fastmath vy[n + 1] = vy[n] + 1/6 * h * imag( f0  + 2∙f1 + 2∙f2 + f3 )
end
@fastmath function euler_forwardsteps_5!(f, vx, vy, h)
    @fastmath for n in 1:(length(vx) - 1)
        euler_step_complex_4!(f, vx, vy, h, n)
    end
end
@fastmath function rrk4_forwardsteps_5!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    @fastmath for n in 1:(length(vx) - 1)
        rk4_step_complex_4!(f, vx, vy, h, n)
    end
end

# Below, we're reusing code to generate the same flow field as in test_23. The cutoff value is no longer NaN.
if !@isdefined m²
    @import_expand(m, s)
end
ϕ_vortex_23 =  generate_complex_potential_vortex(; pos = complex(0.0, 1.0)m, vorticity = 1.0m²/s / 2π)
ϕ_source_23 = generate_complex_potential_source(; pos = complex(3.0, 0.0)m, massflowout = 1.0m²/s)
ϕ_sink_23 = generate_complex_potential_source(; pos = -complex(3.0, 0.0)m, massflowout = -1.0m²/s)
ϕ_23(p) = ϕ_vortex_23(p) + ϕ_source_23(p) + ϕ_sink_23(p)

PHYSWIDTH_23 = 10.0m
HEIGHT_RELATIVE_WIDTH_23 = 0.4
PHYSHEIGHT_23 = PHYSWIDTH_23 * HEIGHT_RELATIVE_WIDTH_23
SCREEN_WIDTH_FRAC_23 = 2 / 3
CUTOFF_23 = 0.5m/s
set_scale_sketch(PHYSWIDTH_23, round(Int, SCREEN_WIDTH_FRAC_23 * WI))

"""
The test flow field from test_23 as a matrix of complex velocities with one element per pixel.
We're filling the above-cutoff values with the cutoff value instead of with NaN.
"""
function flowfield_23()
    unclamped = ∇_rectangle(ϕ_23,
        physwidth = PHYSWIDTH_23,
        physheight  = PHYSWIDTH_23 * HEIGHT_RELATIVE_WIDTH_23);
    map(unclamped) do u
        hypot(u) > CUTOFF_23 ? CUTOFF_23∙u / hypot(u) : u
    end
end


text:
# visualize 90 degree yaw with ROV as bimodal example.

using Caesar
using RoME
using TransformUtils
using IncrementalInference
using Distributions
# using Rotations
# using CoordinateTransformations

# start visualizer with some scene
vc = startdefaultvisualization(draworigin=false)
sleep(3.0)
# load necessary data
rovt = loadmodel(:rov)
rovt(vc)
# defaultscene01!(vc)
sc01 = loadmodel(:scene01)
sc01(vc)

# ground truth positions
gt = Dict{Symbol, Vector{Float64}}()
gt[:x1] = [0.0;0.0;0;0;0;pi/2]
gt[:x2] = [-4.0;4.0;0;0;0;0]
gt[:l1] = [0.0;4.0;0.7]

# am = Translation(0,0,0.0) ∘ LinearMap(Rotations.AngleAxis(pi/2,0,0,1.0))
# rovt(vc, am )

# covariance parameters
initCov = 0.05*eye(6); [initCov[i,i] = 0.001 for i in 4:6];
odoCov = 0.001*eye(6); [odoCov[i,i] = 0.001 for i in 4:6];
rangecov, bearingcov = 3e-2, 2e-2
tfx1 = SE3(gt[:x1][1:3],Euler(gt[:x1][4:6]) )
tfx2 = SE3(gt[:x2][1:3],Euler(gt[:x2][4:6]) )

# factor graph
N=200
fg = identitypose6fg(initpose=tfx1, initCov=initCov, N=N)

drawposepoints!(vc, fg, :x1)

addOdoFG!(fg, Pose3Pose3(MvNormal([4.0;4;0;  0;0;-pi/2], odoCov) ) )

visualizeallposes!(vc, fg, drawlandms=false)
addLinearArrayConstraint(fg, (4.23, 0.0), :x1, :l1, rangecov=rangecov,bearingcov=bearingcov)
visualizeDensityMesh!(vc, fg, :l1)

drawposepoints!(vc, fg, :x2)

# animate visualization
# @async begin
  as = ArcPointsRangeSolve(
        [0.0;0;0.0],
        [-4.0/sqrt(2);4-4.0/sqrt(2);0.0],
        [-4.0;4.0;0.0], 4.0 )
  findaxiscenter!(as)
  animatearc(vc, rovt, as, initrot=Rotations.AngleAxis(pi/2,0,0,1.0), delaytime=0.02)
# end

addLinearArrayConstraint(fg, (4.23, 0.0), :x2, :l1, rangecov=rangecov,bearingcov=bearingcov)
visualizeallposes!(vc, fg, drawlandms=false)

pose2prior = PriorPose3( MvNormal(veeEuler(tfx2), initCov) )
addFactor!(fg,[getVert(fg, :x2)], pose2prior)

solveTree!(fg)
visualize(fg, vc, densitymeshes=[:l1], N=N)

deletemeshes!(vc)

#


text:
using IntervalSets
import CompactBases: RightContinuous, find_interval, within_interval, within_support

"""
    within_interval_linear(x, interval)

Return the range of indices of elements of the sorted vector/range `x`
that fall within `interval`. This loops through all elements of `x`
(linear complexity), so it should not be used for large number of
elements.
"""
function within_interval_linear(x, interval)
    sel = findall(e -> e ∈ interval, x)
    isempty(sel) ? (1:0) : (minimum(sel):maximum(sel))
end

function test_within_interval(x, interval, expected=nothing; reversed=true)
    if reversed
        x = reverse(x)
    end
    N = length(x)

    linear_expected = within_interval_linear(x, interval)
    if isnothing(expected)
        expected = linear_expected
    else
        if reversed
            expected = (N-expected[end]+1):(N-expected[1]+1)
        end
        # This is mostly a sanity check to test if the expectations
        # are actually correct.
        expected ≠ linear_expected &&
            throw(ArgumentError("Invalid expectations: $(expected), actual: $(linear_expected)"))
    end

    result = :(within_interval($x, $interval))
    expr = :($result == $expected || isempty($result) && isempty($expected))
    if !(@eval $expr)
        println("Looking for elements of $x ∈ $interval, got $(@eval $result), expected $expected")
        length(x) < 30 && println("    x = ", collect(enumerate(x)), "\n")
    end
    @eval @test $expr
end

@testset "Knot sets" begin
    @testset "Simple tests" begin
        @test ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3) == [0.0, 1, 1, 3, 4, 6, 6, 6]
        @test LinearKnotSet(3, 0, 1, 2) == [0,0,0,0.5,1,1,1]
        @test LinearKnotSet(2, 0, 1, 2, 1, 1) == [0,0.5,1]
        @test ExpKnotSet(2, -4, 2, 7) == [0,0,0.0001,0.001,0.01,0.1,1,10,100,100]
    end

    @testset "k = $k" for k ∈ 1:6
        @testset "Full multiplicity" begin
            t = LinearKnotSet(k, 0, 1, 1)
            @test length(t) == 2k
            @test collect(t) == vcat(fill(0, k), fill(1, k))
            @test length(nonempty_intervals(t)) == 1
        end
        @testset "Simple endpoints" begin
            t = LinearKnotSet(k, 0, 1, k, 1, 1)
            @test length(t) == k+1
            @test collect(t) == range(0, stop=1, length=k+1)
            @test length(nonempty_intervals(t)) == k
        end
    end
    @test length(nonempty_intervals(ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3))) == 4

    @testset "Invalid order/multiplicities" begin
        @testset "Invalid order" begin
            @test_throws ArgumentError LinearKnotSet(0, 0, 1, 4)
            @test_throws ArgumentError LinearKnotSet(-1, 0, 1, 4)
        end
        @testset "Invalid multipicities" begin
            @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 0, 1)
            @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 1, 0)
            @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 1, 2)
        end
        @testset "Not enough interior points" begin
            @test_throws ArgumentError LinearKnotSet(2, 0, 1, 1, 1, 1)
        end
    end

    @testset "Properties" begin
        @testset "#intervals = $N" for N = 1:4
            @testset "k = $k" for k = 1:N
                @testset "Simple multiplicities" begin
                    t = LinearKnotSet(k, 0.0, 1.0, N, 1, 1)
                    @test numintervals(t) == N
                    @test numfunctions(t) == N - k + 1
                end
            end
        end
    end

    @testset "Compact support" begin
        a,b = 0,1
        x = range(a, stop=b, length=21)
        @testset "Interval coverage" begin
            @testset "Two intervals" begin
                @testset "Reversed: $reversed" for reversed in [false, true]
                    test_within_interval(x, 0..0.5, 1:11, reversed=reversed)
                end
                @testset "L=$L" for L=[:closed,:open]
                    @testset "R=$R" for R=[:closed,:open]
                        @testset "Reversed: $reversed" for reversed in [false, true]
                            test_within_interval(x, Interval{L,R}(0,0.5), reversed=reversed)
                            test_within_interval(x, Interval{L,R}(0.25,0.5), reversed=reversed)
                        end
                    end
                end
            end
            @testset "Three intervals" begin
                @testset "Reversed: $reversed" for reversed in [false, true]
                    test_within_interval(x, RightContinuous(0,1/3), 1:7, reversed=reversed)
                    test_within_interval(x, RightContinuous(1/3,2/3), 8:14, reversed=reversed)
                    test_within_interval(x, 2/3..1, 15:21, reversed=reversed)
                end
            end
            @testset "Open interval" begin
                @testset "Reversed: $reversed" for reversed in [false, true]
                    test_within_interval(x, OpenInterval(0.2,0.4), 6:8, reversed=reversed)
                end
            end
            @testset "Random intervals" begin
                @testset "L=$L" for L=[:closed,:open]
                    @testset "R=$R" for R=[:closed,:open]
                        @testset "Reversed: $reversed" for reversed in [false, true]
                            for i = 1:20
                                interval = Interval{L,R}(minmax(rand(),rand())...)
                                test_within_interval(x, interval, reversed=reversed)
                            end
                        end
                    end
                end
            end
        end

        @testset "Partially covered intervals" begin
            k = 3
            t = LinearKnotSet(k, 0, 1, 2)
            @testset "$name, x = $x" for (name,x) in [
                ("Outside left",range(-1,stop=-0.5,length=10)),
                ("Touching left",range(-1,stop=0,length=10)),
                ("Touching left-ϵ",range(-1,stop=0-eps(),length=10)),
                ("Touching left+ϵ",range(-1,stop=0+eps(),length=10)),

                ("Outside right",range(1.5,stop=2,length=10)),
                ("Touching right",range(1,stop=2,length=10)),
                ("Touching right-ϵ",range(1-eps(),stop=2,length=10)),
                ("Touching right+ϵ",range(1+eps(),stop=2,length=10)),

                ("Other right",range(0.5,stop=1,length=10)),
                ("Other right-ϵ",range(0.5-eps(0.5),stop=1,length=10)),
                ("Other right+ϵ",range(0.5+eps(0.5),stop=1,length=10)),

                ("Complete", range(0,stop=1,length=10)),
                ("Complete-ϵ", range(eps(),stop=1-eps(),length=10)),
                ("Complete+ϵ", range(-eps(),stop=1+eps(),length=10)),

                ("Left partial", range(-0.5,stop=0.6,length=10)),
                ("Left", range(-0.5,stop=1.0,length=10)),
                ("Right partial", range(0.5,stop=1.6,length=10)),
                ("Right", range(0,stop=1.6,length=10))]
                @testset "L=$L" for L=[:closed,:open]
                    @testset "R=$R" for R=[:closed,:open]
                        @testset "Reversed: $reversed" for reversed in [false, true]
                            for i in nonempty_intervals(t)
                                interval = Interval{L,R}(t[i], t[i+1])
                                test_within_interval(x, interval, reversed=reversed)
                            end
                        end
                    end
                end
            end
        end
    end

    @testset "Finding intervals" begin
        t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3)
        R = BSpline(t)
        # Test without specifying initial interval
        @test find_interval(t, 0.5) == 1
        for (i,xs) in enumerate([[0, 0.5, 1.0-eps()],
                                 [],
                                 [1.0, 2.0, 3.0-eps(3.0)],
                                 [3, 3.5, 4.0-eps(4.0)],
                                 [4, 5, 6.0-eps(6.0),6]])
            for j in [1,i] # Test starting from the beginning, and from the interval itself
                for x in xs
                    result = :(find_interval($t, $x, $j))
                    res = @eval $result
                    expr = :($result == $i)
                    (@eval $expr) ||
                        println("Expected to find $x ∈ interval #$i, got $(isnothing(res) ? "nothing" : res)")
                    @eval @test $expr
                    # @test find_interval(t, x) ==
                    find_interval(t, x)
                    within_support(x:x, t, i)[1][2]
                end
            end
        end
        @test isnothing(find_interval(t, 0.5, 2))
    end

    @testset "Support of Heavyside splines" begin
        a,b = 0,1
        x = range(a, stop=b, length=21)
        t = LinearKnotSet(1, a, b, 4)
        supports = [within_support(x, t, j)
                    for j = 1:numfunctions(t)]
        @test length(supports) == 4
        # Each basis function should cover one interval only (since order = 1).
        @test all(length.(supports) .== 1)
        # Test that all elements of x are covered by the basis
        # functions, and that none of the basis functions overlap.
        @test first.(first.(supports)) == [1:5, 6:10, 11:15, 16:21]
    end

    @testset "Support of linear splines" begin
        a,b = 0,1
        x = range(a, stop=b, length=21)
        @testset "Simple multiplicity" begin
            t = LinearKnotSet(2, a, b, 2, 1, 1)
            supports = [within_support(x, t, j)
                        for j = 1:numfunctions(t)]
            @test length(supports) == 1
            @test supports[1] == [(1:10,1), (11:21,2)]
        end
        @testset "Full multiplicity" begin
            @testset "#intervals = 2" begin
                t = LinearKnotSet(2, a, b, 2)
                supports = [within_support(x, t, j)
                            for j = 1:numfunctions(t)]
                @test length(supports) == 3
                @test supports[1] == [(1:10,2)]
                @test supports[2] == [(1:10,2), (11:21,3)]
                @test supports[3] == [(11:21,3)]
            end
            @testset "#intervals = 3" begin
                t = LinearKnotSet(2, a, b, 3)
                supports = [within_support(x, t, j)
                            for j = 1:numfunctions(t)]
                @test length(supports) == 4
                @test supports[1] == [(1:7,2)]
                @test supports[2] == [(1:7,2), (8:14,3)]
                @test supports[3] == [(8:14,3), (15:21,4)]
                @test supports[4] == [(15:21,4)]
            end
        end
    end
end

@testset "Quadrature" begin
    t = LinearKnotSet(1, 0, 1, 2)
    x,w = CompactBases.lgwt(t, 2)
    @test all(w .== 1/4)
    @test x ≈ [-1,1,-1,1]/(4*√3) + [1,1,3,3]/4 rtol=1e-13

    @test_logs (:warn, "N = 1 quadrature point not enough to calculate overlaps between polynomials of order k = 3") CompactBases.lgwt(LinearKnotSet(3, 0, 1, 3), 1)

    @test_logs (:warn, "N = 2 quadrature points not enough to calculate overlaps between polynomials of order k = 3") CompactBases.lgwt(LinearKnotSet(3, 0, 1, 3), 2)
end


text:
"""
DuckDB value
"""
mutable struct Value
    handle::duckdb_value

    function Value(handle::duckdb_value)
        result = new(handle)
        finalizer(_destroy_value, result)
        return result
    end
end

function _destroy_value(val::Value)
    if val.handle != C_NULL
        duckdb_destroy_value(val.handle)
    end
    val.handle = C_NULL
    return
end

GetValue(val::Value, ::Type{T}) where {T <: Int64} = duckdb_get_int64(val.handle)
function GetValue(val::Value, ::Type{T}) where {T <: String}
    ptr = duckdb_get_varchar(val.handle)
    result = unsafe_string(ptr)
    duckdb_free(ptr)
    return result
end
function GetValue(val::Value, ::Type{T}) where {T}
    throw(NotImplementedException("Unsupported type for GetValue"))
end

CreateValue(val::T) where {T <: Int64} = Value(duckdb_create_int64(val))
function CreateValue(val::T) where {T}
    throw(NotImplementedException("Unsupported type for GetValue"))
end


text:
function evaluate_decision(structure::DeterministicEquivalent{2}, decision::AbstractVector)
    # Update decisions (checks handled by structure model)
    take_decisions!(structure.model, all_decision_variables(structure.model, 1), decision)
    # Optimize model
    optimize!(structure.model)
    # Switch on return status
    status = termination_status(structure.model)
    result = if status in AcceptableTermination
        result = objective_value(structure.model)
    else
        result = if status == MOI.INFEASIBLE
            result = objective_sense(structure.model) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            result = objective_sense(structure.model) == MOI.MAX_SENSE ? Inf : -Inf
        else
            error("Deterministically equivalent model could not be solved, returned status: $status")
        end
    end
    # Revert back to untaken decisions
    untake_decisions!(structure.model, all_decision_variables(structure.model, 1))
    # Return evaluation result
    return result
end

function statistically_evaluate_decision(structure::DeterministicEquivalent{2}, decision::AbstractVector)
    # Update decisions (checks handled by structure model)
    take_decisions!(structure.model, all_decision_variables(structure.model, 1), decision)
    # Optimize model
    optimize!(structure.model)
    # Get sense-correted objective value
    status = termination_status(structure.model)
    Q̂ = if status in AcceptableTermination
        Q̂ = objective_value(structure.model)
    else
        Q̂ = if status == MOI.INFEASIBLE
            return objective_sense(structure.model) == MOI.MAX_SENSE ? (-Inf, 0) : (Inf, 0)
        elseif status == MOI.DUAL_INFEASIBLE
            return objective_sense(structure.model) == MOI.MAX_SENSE ? (Inf, 0) : (-Inf, 0)
        else
            error("Deterministically equivalent model could not be solved, returned status: $status")
        end
    end
    # Calculate subobjectives
    N = num_scenarios(structure)
    Q = Vector{Float64}(undef, N)
    obj_sense = objective_sense(structure.model)
    for (i, sub_objective) in enumerate(structure.sub_objectives[2])
        (sub_sense, sub_obj) = sub_objective
        Qᵢ = MOIU.eval_variables(sub_obj) do idx
            return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
        end
        if obj_sense == sub_sense
            Q[i] = Qᵢ
        else
            Q[i] = -Qᵢ
        end
    end
    probabilities = map(1:num_scenarios(structure, 2)) do i
        probability(structure, 2, i)
    end
    weights = Distributions.StatsBase.ProbabilityWeights(probabilities)
    σ = std(Q, weights, corrected = true)
    return Q̂, σ
end


text:
"""
Jacobienne des contraintes actives à precmpcc près
"""
function jac_actif(rlxmpcc :: RlxMPCC, x :: Vector, prec :: Float64)

  n = rlxmpcc.n
  ncc = rlxmpcc.ncc
  r, s, t, tb = rlxmpcc.r, rlxmpcc.s, rlxmpcc.t, rlxmpcc.tb
  mod = rlxmpcc.mod

  Il = findall(z->z<=prec,abs.(x - mod.meta.lvar))
  Iu = findall(z->z<=prec,abs.(x - mod.meta.uvar))
  jl = zeros(n);jl[Il].=1.0;Jl=diagm(0 => jl);
  ju = zeros(n);ju[Iu].=1.0;Ju=diagm(0 => ju);

  IG=[];IH=[];IPHI=[];Ig=[];Ih=[];

 if mod.meta.ncon+ncc ==0

  A=[]

 else

  if mod.meta.ncon > 0
   c = cons_nl(mod, x)
   J = jac_nl(mod, x)
  else
   c = Float64[]
   J = sparse(zeros(0,2))
  end

  Ig=findall(z->z<=prec,abs.(c- mod.meta.lcon))
  Ih=findall(z->z<=prec,abs.(c- mod.meta.ucon))

  Jg, Jh = zeros(mod.meta.ncon,n), zeros(mod.meta.ncon,n)

  Jg[Ig,1:n] = J[Ig,1:n]
  Jh[Ih,1:n] = J[Ih,1:n]

  if ncc>0

   G = consG(mod,x)
   H = consH(mod,x)
   IG   = findall(z->z<=prec,abs.(G - tb*ones(ncc) - mod.meta.lccG))
   IH   = findall(z->z<=prec,abs.(H - tb*ones(ncc) - mod.meta.lccH))
   IPHI = findall(z->z<=prec,abs.(phi(G, H, r, s, t)))

   JPHI = dphi(G,H,rlxmpcc.r,rlxmpcc.s,rlxmpcc.t)
   JPHIG, JPHIH = JPHI[1:ncc][IPHI], JPHI[ncc+1:2*ncc][IPHI]

   JG   = jacG(mod,x)
   JH   = jacH(mod,x)

   JBG, JBH = zeros(ncc, n), zeros(ncc, n)
   JBG[IG,1:n] = JG[IG,1:n]
   JBH[IH,1:n] = JH[IH,1:n]

   JP = zeros(ncc, n)
   JP[IPHI,1:n] = diagm(0 => JPHIG)*JG[IPHI,1:n] + diagm(0 => JPHIH)*JH[IPHI,1:n]

   A=[Jl;Ju;-Jg;Jh;-JBG;-JBH;JP]'

  else

   A=[Jl;Ju;-Jg;Jh]' #c'est à l'envers!

  end
 end

 return A, Il,Iu,Ig,Ih,IG,IH,IPHI
end

function jac_actif2(rlxmpcc :: RlxMPCC, x :: Vector, prec :: Float64)

  n = rlxmpcc.n
  ncc = rlxmpcc.ncc
  r, s, t, tb = rlxmpcc.r, rlxmpcc.s, rlxmpcc.t, rlxmpcc.tb
  mod = rlxmpcc.mod

  IG=[];IH=[];IPHI=[];Ig=[];Ih=[];

 if mod.meta.ncon+ncc ==0

  A=[]

 else

  if mod.meta.ncon > 0
   c = cons_nl(mod, x)
   J = jac_nl(mod, x)
  else
   c = Float64[]
   J = sparse(zeros(0,2))
  end

  Ig=findall(z->z<=prec,abs.(c- mod.meta.lcon))
  Ih=findall(z->z<=prec,abs.(c- mod.meta.ucon))

  Jg, Jh = zeros(mod.meta.ncon,n), zeros(mod.meta.ncon,n)

  Jg[Ig,1:n] = J[Ig,1:n]
  Jh[Ih,1:n] = J[Ih,1:n]

  if ncc>0

   G = consG(mod,x)
   H = consH(mod,x)

   JPHI = dphi(G,H,rlxmpcc.r,rlxmpcc.s,rlxmpcc.t)
   JPHIG, JPHIH = JPHI[1:ncc], JPHI[ncc+1:2*ncc]

   JG   = jacG(mod,x)
   JH   = jacH(mod,x)
   JP = diagm(0 => JPHIG)*JG + diagm(0 => JPHIH)*JH

   A=[J;JG;JH;JP]

  else

   A=J

  end
 end

 return A
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

# adapted from https://github.com/FluxML/Zygote.jl/blob/master/test/gradcheck.jl
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "gate gradients" begin

    @testset "single qubit gates" begin

        # fictitious gradients of cost function with respect to quantum gate
        Δ = randn(ComplexF64, 2, 2)

        for g in [RxGate, RyGate, RzGate]
            f = (θ) -> 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[]))))
            θ = 2π*rand()
            ngrad = ngradient(f, [θ])
            dg = Qaintessent.backward(g(θ), conj(Δ))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        for g in [XGate, YGate, ZGate, HadamardGate, SGate, SdagGate, TGate, TdagGate]
            dg = Qaintessent.backward(g(), conj(Δ))
            @test dg ≈ g()
        end

        begin
            f = (ϕ) -> 2*real(sum(Δ .* Qaintessent.sparse_matrix(PhaseShiftGate(ϕ[]))))
            ϕ = 2π*rand()
            ngrad = ngradient(f, [ϕ])
            dg = Qaintessent.backward(PhaseShiftGate(ϕ), conj(Δ))
            @test isapprox(dg.ϕ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        begin
            f = (nθ) -> 2*real(sum(Δ .* Qaintessent.sparse_matrix(RotationGate(nθ))))
            nθ = randn(3)
            ngrad = ngradient(f, nθ)
            dg = Qaintessent.backward(RotationGate(nθ), conj(Δ))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-6)
        end

        begin
            f = (nθ) -> 2*real(sum(Δ .* Qaintessent.sparse_matrix(RotationGate(nθ))))
            # special case: zero vector
            nθ = zeros(3)
            ngrad = ngradient(f, nθ)
            dg = Qaintessent.backward(RotationGate(nθ), conj(Δ))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-6)
        end
    end

    @testset "two qubit gates" begin
        # fictitious gradients of cost function with respect to quantum gate
        Δ = randn(ComplexF64, 4, 4)
        @testset "entanglement gates" begin
            for g in [EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
                f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[]))))
                θ = 2π*rand()
                ngrad = ngradient(f, [θ])
                dg = Qaintessent.backward(g(θ), conj(Δ))
                @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)
            end
        end

        @testset "two qubit gates" begin
            for g in [SwapGate()]
                dg = Qaintessent.backward(g, conj(Δ))
                @test dg ≈ g
            end
        end

    end

    @testset "controlled gates" begin
        # fictitious gradients of cost function with respect to quantum gate
        Δ = randn(ComplexF64, 8, 8)

        for g in [RxGate, RyGate, RzGate]
            f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(ControlledGate{g}(g(θ[]), 2))))
            θ = 2π*rand()
            ngrad = ngradient(f, [θ])
            dg = Qaintessent.backward(ControlledGate{g}(g(θ), 2), conj(Δ))
            @test isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        for g in [EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
            f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(ControlledGate{g}(g(θ[]), 1))))
            θ = 2π*rand()
            ngrad = ngradient(f, [θ])
            dg = Qaintessent.backward(ControlledGate{g}(g(θ), 2), conj(Δ))
            @test isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        for g in [XGate, YGate, ZGate, HadamardGate, SGate, SdagGate, TGate, TdagGate]

            dg = Qaintessent.backward(ControlledGate{g}(g(), 2), conj(Δ))

            @test dg ≈ ControlledGate{g}(g(), 2)
        end
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gradients" begin

    # construct parametrized circuit
    N = 4
    A = rand(ComplexF64, 2 ,2)
    U, R = qr(A)
    U = Array(U)
    i = rand(1:N)

    cgc(θ, ϕ, χ, κ, ωn) = CircuitGate[
        circuit_gate(3, HadamardGate()),
        circuit_gate(2, RzGate(θ), (1, 4)),
        circuit_gate(2, TGate(), 4),
        circuit_gate(2, 3, SwapGate()),
        circuit_gate(2, SdagGate(), 1),
        circuit_gate(3, PhaseShiftGate(ϕ)),
        circuit_gate(1, RyGate(χ)),
        circuit_gate(4, HadamardGate(), 2),
        circuit_gate(3, X, 1),
        circuit_gate(2, Z, 4),
        circuit_gate(4, Y, 1),
        circuit_gate(3, 1, EntanglementYYGate(κ)),
        circuit_gate(3, RotationGate(ωn)),
        circuit_gate(4, SGate()),
        circuit_gate(2, TdagGate(), 3),
        circuit_gate(i, MatrixGate(U)),
        circuit_gate(2, Z),
        circuit_gate(3, X),
        circuit_gate(2, Y, 4),
        circuit_gate(1, SGate()),
        circuit_gate(2, 3, SwapGate(), 4),
        circuit_gate(3, SdagGate()),
        circuit_gate(2, TGate()),
    ]

    # measurement operators
    meas(M) = MeasurementOperator.([Matrix{Float64}(I, 2^N, 2^N), Hermitian(M)], (Tuple(1:N),))

    # parameter values
    θ = 1.5π
    ϕ = 0.3
    χ = √2
    κ = exp(-0.4)
    n = randn(Float64, 3)
    n /= norm(n)
    ωn = 0.2π * n
    M = randn(ComplexF64, 2^N, 2^N)
    M = 0.5*(M + adjoint(M))

    c = Circuit{N}(cgc(θ, ϕ, χ, κ, ωn), meas(M))
    # input quantum state
    ψ = randn(ComplexF64, 2^N)
    ψl = deepcopy(ψ)
    # fictitious gradients of cost function with respect to circuit output
    Δ = [0.3, -1.2]

    dc = Qaintessent.gradients(c, ψ, Δ)[1]

    f(rθ, rϕ, rχ, rκ, ωn, M) = dot(Δ, apply(ψ, Circuit{N}(cgc(rθ[], rϕ[], rχ[], rκ[], ωn), meas(M))))
    # numeric gradients
    ngrad = ngradient(f, [θ], [ϕ], [χ], [κ], ωn, M)
    # symmetrize gradient with respect to measurement operator
    ngrad[end][:] = 0.5*(ngrad[end] + adjoint(ngrad[end]))

    @test all(isapprox.(ngrad,
        (dc.moments[1][2].gate.U.θ,
         dc.moments[4][2].gate.ϕ,
         dc.moments[5][1].gate.θ,
         dc.moments[8][1].gate.θ,
         dc.moments[9][1].gate.nθ,
         sparse_matrix(dc.meas[2])), rtol=1e-5, atol=1e-5))
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gradients with moments" begin

   # construct parametrized circuit
   N = 4
    cgc(θ, ϕ, χ, ωn) = Moment[
        Moment([
            circuit_gate(3, HadamardGate()),
            circuit_gate(2, RzGate(θ), (1, 4)),
        ]),
        Moment(circuit_gate(2, 3, SwapGate())),
        Moment(circuit_gate(3, PhaseShiftGate(ϕ))),
        Moment([
            circuit_gate(3, RotationGate(ωn)),
            circuit_gate(1, RyGate(χ)),
        ]),
    ]

   # measurement operators
   meas(M) = MeasurementOperator.([Matrix{Float64}(I, 2^N, 2^N), Hermitian(M)], (Tuple(1:N),))

   # parameter values
   θ = 1.5π
   ϕ = 0.3
   χ = √2
   n = randn(Float64, 3)
   n /= norm(n)
   ωn = 0.2π * n
   M = randn(ComplexF64, 2^N, 2^N)
   M = 0.5*(M + adjoint(M))

   c = Circuit{N}(cgc(θ, ϕ, χ, ωn), meas(M))

   # input quantum state
   ψ = randn(ComplexF64, 2^N)

   # fictitious gradients of cost function with respect to circuit output
   Δ = [0.3, -1.2]

   dc = Qaintessent.gradients(c, ψ, Δ)[1]

   f(rθ, rϕ, rχ, ωn, M) = dot(Δ, apply(ψ, Circuit{N}(cgc(rθ[], rϕ[], rχ[], ωn), meas(M))))
   # numeric gradients
   ngrad = ngradient(f, [θ], [ϕ], [χ], ωn, M)
   # symmetrize gradient with respect to measurement operator
   ngrad[end][:] = 0.5*(ngrad[end] + adjoint(ngrad[end]))

   @test all(isapprox.(ngrad,
       (dc.moments[1][2].gate.U.θ,
        dc.moments[3][1].gate.ϕ,
        dc.moments[4][2].gate.θ,
        dc.moments[4][1].gate.nθ,
        sparse_matrix(dc.meas[2])), rtol=1e-5, atol=1e-5))
end


text:
using Leadfield
using Winston
using TextPlots

positions, neighbours = readLOC("/Users/rluke/.julia/v0.3/Leadfield/data/FEM_Default50MRI.loc", verbose=true)

println()
println()
println("Back")
plot(vec(positions[1,:]), vec(positions[3,:]), cols=27, rows=15)

println()
println()
println("Side")
plot(vec(positions[2,:]), vec(positions[3,:]), cols=27, rows=15)

println()
println()
println("Top")
plot(vec(positions[1,:]), vec(positions[2,:]), cols=27, rows=15)


text:
##
# This file is a part of Dagger.jl. License is MIT
#
# Based upon https://github.com/google/benchmark, which is licensed under Apache v2:
# https://github.com/google/benchmark/blob/master/LICENSE
#
# In compliance with the Apache v2 license, here are the original copyright notices:
# Copyright 2015 Google Inc. All rights reserved.
##

struct TimeSpec
    tv_sec  :: UInt64 # time_t
    tv_nsec :: UInt64
end

maketime(ts) = ts.tv_sec * UInt(1e9) + ts.tv_nsec

# From bits/times.h on a Linux system
# Check if those are the same on BSD
if Sys.islinux()
    const CLOCK_MONOTONIC          = Cint(1)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(2)
    const CLOCK_THREAD_CPUTIME_ID  = Cint(3)
elseif Sys.isfreebsd() # atleast on FreeBSD 11.1
    const CLOCK_MONOTONIC          = Cint(4)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(14)
elseif Sys.isapple() # Version 10.12 required
    const CLOCK_MONOTONIC          = Cint(6)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(12)
end

@static if Sys.isunix()
@inline function clock_gettime(cid)
    ts = Ref{TimeSpec}()
    ccall(:clock_gettime, Cint, (Cint, Ref{TimeSpec}), cid, ts)
    return ts[]
end

@inline function realtime()
    maketime(clock_gettime(CLOCK_MONOTONIC))
end

@inline function cputime()
    maketime(clock_gettime(CLOCK_PROCESS_CPUTIME_ID))
end
end

@static if Sys.islinux()
@inline function cputhreadtime()
    maketime(clock_gettime(CLOCK_THREAD_CPUTIME_ID))
end
else
@inline cputhreadtime() = time_ns() # HACK
end


text:
using Plots, Test

@testset "Showaxis" begin
    for value in Plots._allShowaxisArgs
        @test plot(1:5, showaxis = value)[1][:yaxis][:showaxis] isa Bool
    end
    @test plot(1:5, showaxis = :y)[1][:yaxis][:showaxis] == true
    @test plot(1:5, showaxis = :y)[1][:xaxis][:showaxis] == false
end

@testset "Magic axis" begin
    @test plot(1, axis=nothing)[1][:xaxis][:ticks] == []
    @test plot(1, axis=nothing)[1][:yaxis][:ticks] == []
end # testset

@testset "Categorical ticks" begin
    p1 = plot('A':'M', 1:13)
    p2 = plot('A':'Z', 1:26)
    p3 = plot('A':'Z', 1:26, ticks = :all)
    @test Plots.get_ticks(p1[1], p1[1][:xaxis])[2] == string.('A':'M')
    @test Plots.get_ticks(p2[1], p2[1][:xaxis])[2] == string.('C':3:'Z')
    @test Plots.get_ticks(p3[1], p3[1][:xaxis])[2] == string.('A':'Z')
end


text:
using NestedTuples: modify, rmap, unwrap, flatten, Leaves, @with, TypelevelExpr
import NestedTuples

export TupleVector
using ElasticArrays, ArraysOfArrays
struct TupleVector{T,X} <: AbstractVector{T}
    data :: X
end

NestedTuples.unwrap(tv::TupleVector) = getfield(tv, :data)

function TupleVector(data::X) where {X <: NamedTuple}
    T = typeof(rmap(first, data))
    return TupleVector{T,X}(data)
end

function TupleVector(a::AbstractVector{T}) where {T}
    a1 = first(a)

    x = TupleVector(undef, a1, length(a))
    x .= a
    return x
end

function Base.similar(a::TupleVector, ::Type{T}, dims::Dims) where {T}
    data = unwrap(a)
    sim = rmap(x->similar(x, dims), data)
    return TupleVector(sim)
end

function TupleVector(::UndefInitializer, x::T, n::Int) where {T<:NamedTuple}

    function initialize(n::Int)
        f(x::T) where {T} = ElasticVector{T}(undef, n)
        f(x::DenseArray{T,N}) where {T,N} = nestedview(ElasticArray{T,N+1}(undef, size(x)..., n), N)
        return f
    end

    data = rmap(initialize(n), x)

    return TupleVector{T, typeof(data)}(data)
end

function TupleVector(; kwargs...)
    return TupleVector(NamedTuple(kwargs))
end

# function TupleVector(x::Union{Tuple, NamedTuple})
#     flattened = flatten(x)
#     @assert allequal(size.(flattened)...)

#     T = typeof(modify(arr -> arr[1], x, Leaves()))
#     N = length(axes(flattened[1]))
#     X = typeof(x)
#     return TupleVector{T,X}(x)
# end

# TupleVector{T}(x...) where {T} = leaf_setter(T)(x...)

import Base

NestedTuples.schema(tv::TupleVector{T}) where {T} = schema(T)

Base.propertynames(tv::TupleVector) = propertynames(unwrap(tv))

function Base.showarg(io::IO, tv::TupleVector{T}, toplevel) where T
    io = IOContext(io, :compact => true)
    print(io, "TupleVector")
    toplevel && println(io, " with schema ", schema(T))
end

Base.show(io::IO, ::MIME"text/plain", tv::TupleVector) = show(io, tv)
Base.show(io::IO, ::MIME"text/html", tv::TupleVector) = show(io, tv)

function Base.show(io::IO, tv::TupleVector)
    summary(io, tv)
    print(io, summarize(tv))
end

function Base.show(io::IO, ::MIME"text/plain", v::Vector{TV}) where {TV <: TupleVector}
    io = IOContext(io, :compact => true)
    n = length(v)
    println(io, n,"-element Vector{$TV}")
    foreach(v) do tv println(io, summarize(tv)) end
end

function Base.getindex(x::TupleVector, j::Int)

    # TODO: Bounds checking doesn't affect performance, am I doing it right?
    function f(arr)
        # @boundscheck all(j .∈ axes(arr))
        return @inbounds arr[j]
    end

    modify(f, unwrap(x), Leaves())
end

function Base.getindex(x::TupleVector, j...)

    # TODO: Bounds checking doesn't affect performance, am I doing it right?
    function f(arr)
        # @boundscheck all(j .∈ axes(arr))
        return @inbounds arr[j...]
    end

    TupleVector(modify(f, unwrap(x), Leaves()))
end

function Base.setindex!(a::TupleVector, x, j::Int)
    a1 = flatten(unwrap(a))
    x1 = flatten(x)

    setindex!.(a1, x1, j)
    return a
end

function Base.length(tv::TupleVector)
    length(flatten(unwrap(tv))[1])
end

function Base.size(tv::TupleVector)
    size(flatten(unwrap(tv))[1])
end

# TODO: Make this pass @code_warntype
Base.getproperty(tv::TupleVector, k::Symbol) = maybewrap(getproperty(unwrap(tv), k))

maybewrap(t::Tuple) = TupleVector(t)
maybewrap(t::NamedTuple) = TupleVector(t)
maybewrap(t) = t

NestedTuples.flatten(tv::TupleVector) = TupleVector(flatten(unwrap(tv)))

# leaf_setter(tv::TupleVector) = TupleVector ∘ leaf_setter(unwrap(tv))

# function TupleVector{T}(::UndefInitializer) where {T}
#     return EmptyTupleVector{T}()
# end

# function Base.push!(::EmptyTupleVector{T}, nt::NamedTuple) where {T}
#     function f(x::t, path) where {t}
#         ea = ElasticArray{t}(undef, 0)
#         push!(ea, x)
#         return ea
#     end

#     function f(x::DenseArray{t}, path) where {t}
#         ea = ElasticArray{t}(undef, size(x)..., 0)
#         nv = nestedview(ea, 1)
#         push!(nv, x)
#         return nv
#     end

#     data = fold(f, nt)
#     X = typeof(data)
#     TupleVector{T,X}(data)
# end

export rmap

function NestedTuples.rmap(f, tv::TupleVector)
    return TupleVector(rmap(f, unwrap(tv)))
end

function Base.resize!(tv::TupleVector, n::Int)
    rmap(x -> resize!(x, n), unwrap(tv))
end


text:
# X and Z components in fT
# conducting 100 Ωm
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 1, 100]

ross_primary = [25.7847     -18.7313]
ross_secondary = [   6.762     -8.0374
  4.7924     -6.6042
  4.0079     -5.9113
  3.4924     -5.3952
  3.1035     -4.9777
  2.7646      -4.608
  2.3721     -4.1681
  1.9079     -3.6087
  1.4085     -2.9386
0.92505     -2.1946
0.52168     -1.4595
0.24675    -0.84605
0.097237    -0.42402
0.032526    -0.18624
0.008896   -0.069836]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
"""
    CircShiftedArray(parent::AbstractArray, shifts)

Custom `AbstractArray` object to store an `AbstractArray` `parent` circularly shifted by `shifts` steps (where `shifts` is
a `Tuple` with one `shift` value per dimension of `parent`).
Use `copy` to collect the values of a `CircShiftedArray` into a normal `Array`.

# Examples

```jldoctest circshiftedarray
julia> v = [1, 3, 5, 4];

julia> s = CircShiftedArray(v, (1,))
4-element CircShiftedArray{Int64,1,Array{Int64,1}}:
 4
 1
 3
 5

julia> copy(s)
4-element Array{Int64,1}:
 4
 1
 3
 5
```
"""
struct CircShiftedArray{T, N, S<:AbstractArray} <: AbstractArray{T, N}
    parent::S
    shifts::NTuple{N, Int}
    function CircShiftedArray(p::AbstractArray{T, N}, n = Tuple(0 for i in 1:N)) where {T, N}
        @assert all(step(x) == 1 for x in axes(p))
        new{T, N, typeof(p)}(p, _padded_tuple(p, n))
    end
end

"""
    CircShiftedVector{T, S<:AbstractArray}

Shorthand for `CircShiftedArray{T, 1, S}`.
"""
const CircShiftedVector{T, S<:AbstractArray} = CircShiftedArray{T, 1, S}

CircShiftedVector(v::AbstractVector, n = (0,)) = CircShiftedArray(v, n)

size(s::CircShiftedArray) = size(parent(s))

@inline function bringwithin(idx::Int, range::AbstractRange)
    a, b = extrema(range)
    n = length(range)
    while idx < a
        idx += n
    end
    while idx > b
        idx -= n
    end
    idx
end

@inline bringwithin(idxs::Tuple, ranges::Tuple) =
    (bringwithin(idxs[1], ranges[1]), bringwithin(Base.tail(idxs), Base.tail(ranges))...)

@inline bringwithin(idxs::Tuple{}, ranges::Tuple{}) = ()

@inline function getindex(s::CircShiftedArray{T, N}, x::Vararg{Int, N}) where {T, N}
    v = parent(s)
    ind = offset(shifts(s), x)
    if checkbounds(Bool, v, ind...)
        @inbounds ret = v[ind...]
    else
        i = bringwithin(ind, axes(v))
        @inbounds ret = v[i...]
    end
    ret
end

@inline function setindex!(s::CircShiftedArray{T, N}, el, x::Vararg{Int, N}) where {T, N}
    v = parent(s)
    ind = offset(shifts(s), x)
    if checkbounds(Bool, v, ind...)
        @inbounds v[ind...] = el
    else
        i = bringwithin(ind, axes(v))
        @inbounds v[i...] = el
    end
end

parent(s::CircShiftedArray) = s.parent

"""
    shifts(s::CircShiftedArray)

Returns amount by which `s` is shifted compared to `parent(s)`.
"""
shifts(s::CircShiftedArray) = s.shifts

checkbounds(::CircShiftedArray, I...) = nothing
checkbounds(::Type{Bool}, ::CircShiftedArray, I...) = true


text:
function summarize_regression_test(field_names, fields, correct_fields)
    for (field_name, φ, φ_c) in zip(field_names, fields, correct_fields)
        Δ = Array(φ) .- φ_c

        Δ_min      = minimum(Δ)
        Δ_max      = maximum(Δ)
        Δ_mean     = mean(Δ)
        Δ_abs_mean = mean(abs, Δ)
        Δ_std      = std(Δ)

        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g",
                       field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std))
    end
end

interior(a, grid) = view(a, grid.Hx+1:grid.Nx+grid.Hx,
                            grid.Hy+1:grid.Ny+grid.Hy,
                            grid.Hz+1:grid.Nz+grid.Hz)

function get_fields_from_checkpoint(filename)
    file = jldopen(filename)

    tracers = keys(file["tracers"])
    tracers = Tuple(Symbol(c) for c in tracers)

    velocity_fields = (
                       u = file["velocities/u/data"],
                       v = file["velocities/v/data"],
                       w = file["velocities/w/data"],
                      )

    tracer_fields = NamedTuple{tracers}(Tuple(file["tracers/$c/data"] for c in tracers))

    current_tendency_velocity_fields = (
                                        u = file["timestepper/Gⁿ/u/data"],
                                        v = file["timestepper/Gⁿ/v/data"],
                                        w = file["timestepper/Gⁿ/w/data"],
                                       )

    current_tendency_tracer_fields = NamedTuple{tracers}(Tuple(file["timestepper/Gⁿ/$c/data"] for c in tracers))

    previous_tendency_velocity_fields = (
                                         u = file["timestepper/G⁻/u/data"],
                                         v = file["timestepper/G⁻/v/data"],
                                         w = file["timestepper/G⁻/w/data"],
                                        )

    previous_tendency_tracer_fields = NamedTuple{tracers}(Tuple(file["timestepper/G⁻/$c/data"] for c in tracers))

    close(file)

    solution = merge(velocity_fields, tracer_fields)
    Gⁿ = merge(current_tendency_velocity_fields, current_tendency_tracer_fields)
    G⁻ = merge(previous_tendency_velocity_fields, previous_tendency_tracer_fields)

    return solution, Gⁿ, G⁻
end

include("regression_tests/thermal_bubble_regression_test.jl")
include("regression_tests/rayleigh_benard_regression_test.jl")
include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")

@testset "Regression" begin
    @info "Running regression tests..."

    for arch in archs
        @testset "Thermal bubble [$(typeof(arch))]" begin
            @info "  Testing thermal bubble regression [$(typeof(arch))]"
            run_thermal_bubble_regression_test(arch)
        end

        @testset "Rayleigh–Bénard tracer [$(typeof(arch))]" begin
            @info "  Testing Rayleigh–Bénard tracer regression [$(typeof(arch))]"
            run_rayleigh_benard_regression_test(arch)
        end

        @testset "Ocean large eddy simulation [$(typeof(arch))]" begin
            for closure in (AnisotropicMinimumDissipation(), ConstantSmagorinsky())
                closurename = string(typeof(closure).name.wrapper)
                @info "  Testing oceanic large eddy simulation regression [$closurename, $(typeof(arch))]"
                run_ocean_large_eddy_simulation_regression_test(arch, closure)
            end
        end
    end
end