text:
using Rayuela
using HDF5
using Printf

# === Saving functions ===
function save_results_pq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, train_error, recall)
  for i = 1:length(C)
    h5write(bpath, "$(trial)/C_$i", C[i])
  end
  h5write(bpath, "$(trial)/B", convert(Matrix{UInt8}, B.-1))
  h5write(bpath, "$(trial)/train_error", train_error)
  h5write(bpath, "$(trial)/recall", recall)
end

function save_results_pq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, train_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_pq_query_base(bpath, trial, C, B, train_error, recall)
end

function save_results_opq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, recall)
  h5write(bpath, "$(trial)/R", R)
  save_results_pq_query_base(bpath, trial, C, B, train_error, recall)
end

function save_results_opq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_opq_query_base(bpath, trial, C, B, R, train_error, recall)
end

function save_results_lsq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, opq_error, recall)
  h5write(bpath, "$(trial)/opq_base", opq_error)
  save_results_opq_query_base(bpath, trial, C, B, R, train_error, recall)
end

function save_results_lsq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, opq_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_lsq_query_base(bpath, trial, C, B, R, train_error, opq_error, recall)
end

function load_chainq(fname::String, m::Integer, trial::Integer)
  B = h5read(fname, "$trial/B"); B = convert(Matrix{Int16}, B); B.+=1
  R = h5read(fname, "$trial/R")
  chainq_error = h5read(fname, "$trial/train_error")
  C = Vector{Matrix{Float32}}(undef, m)
  for i=1:(m); C[i] = h5read(fname, "$trial/C_$i"); end
  return C, B, R, chainq_error
end

function load_rvq(fname::String, m::Integer, trial::Integer)
  B = h5read(fname, "$trial/B"); B = convert(Matrix{Int16}, B); B.+=1
  rvq_error = h5read(fname, "$trial/train_error")
  C = Vector{Matrix{Float32}}(m)
  for i=1:(m); C[i] = h5read(fname, "$trial/C_$i"); end
  return C, B, rvq_error
end

# === Data loading function ===
function load_experiment_data(
  dataset_name::String,
  ntrain::Integer, nbase::Integer, nquery::Integer, V::Bool=false)

  Xt = Rayuela.read_dataset(dataset_name, ntrain, V)
  Xb = Rayuela.read_dataset(dataset_name * "_base", nbase, V)
  Xq = Rayuela.read_dataset(dataset_name * "_query", nquery, V)[:,1:nquery]
  gt = Rayuela.read_dataset(dataset_name * "_groundtruth", nquery, V)
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt .+ 1
  end

  if dataset_name != "Deep1M"
    gt = convert( Vector{UInt32}, gt[1,1:nquery] )
  else
    gt = convert( Vector{UInt32}, gt[1:nquery] )
  end

  return Xt, Xb, Xq, gt
end

# === experiment functions ===
function run_demos(
  dataset_name="SIFT1M",
  ntrain::Integer=Int(1e5),
  m::Integer=8, h::Integer=256, niter::Integer=25)

  nquery, nbase, knn = 0, 0, 0
  if dataset_name == "SIFT1M" || dataset_name == "Deep1M" || dataset_name == "Convnet1M"
    nquery, nbase, knn = Int(1e4), Int(1e6), Int(1e3)
  else
    error("dataset unknown")
  end

  verbose = true

  # Load data
  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  d, _    = size( Xt )

  ntrials = 1 # Increase this to 10 to reproduce ECCV'18 results
  # (Semi-)orthogonal methods: PQ, OPQ, ChainQ
  for trial = 1:ntrials
    C, B, train_error, B_base, recall = Rayuela.experiment_pq(Xt, Xb, Xq, gt, m, h, niter, knn, verbose)
    save_results_pq("./results/$(lowercase(dataset_name))/pq_m$(m)_it$(niter).h5", trial, C, B, train_error, B_base, recall)
  end
  for trial = 1:ntrials
    init = "natural"
    C, B, R, train_error, B_base, recall = Rayuela.experiment_opq(Xt, Xb, Xq, gt, m, h, init, niter, knn, verbose)
    save_results_opq("./results/$(lowercase(dataset_name))/opq_m$(m)_it$(niter).h5", trial, C, B, R, train_error, B_base, recall)
  end

  # Cheap non-orthogonal methods: RVQ, ERVQ
  for trial = 1:ntrials
    C, B, train_error, B_base, recall = Rayuela.experiment_rvq( Xt,       Xb, Xq, gt, m-1, h, niter, knn, verbose)
    save_results_pq("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5",  trial, C, B, train_error, B_base, recall)
  end
  for trial = 1:ntrials
    C, B, train_error = load_rvq("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, train_error, B_base, recall = Rayuela.experiment_ervq(Xt, B, C, Xb, Xq, gt, m-1, h, niter, knn, verbose)
    save_results_pq("./results/$(lowercase(dataset_name))/ervq_m$(m-1)_it$(niter).h5", trial, C, B, train_error, B_base, recall)
  end

  # Precompute init for LSQ/SR
  for trial = 1:ntrials
    C, B, R, train_error = Rayuela.train_opq(Xt, m-1, h, niter, "natural", verbose)
    save_results_opq("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, ones(UInt16,1,1), [0f0])
  end
  for trial = 1:ntrials
    C, B, R, _ = load_chainq("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, chainq_error = train_chainq(    Xt, m-1, h, R, B, C, niter, verbose)
    save_results_opq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", trial, C, B, R, chainq_error, ones(UInt16,1,1), [0f0])
  end

  nsplits_train =  m == 8 ? 1 : 1
  nsplits_base  =  m == 8 ? 2 : 4
  @show nsplits_train, nsplits_base

  ilsiter = 8
  icmiter = 4
  randord = true
  npert   = 4

  for trial = 1:ntrials

    # C, B, R, _ = load_chainq("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(25).h5", m-1, trial)
    # B_base = h5read("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(25).h5", "$trial/B_base")
    # B_base = convert(Matrix{Int16}, B_base) .+ 1
    # @show qerror(Xt, B, C), qerror(Xb, B_base, C)

    C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
    # C, B, R, train_error, B_base, recall = Rayuela.experiment_lsq(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, knn, verbose)
    C, B, R, train_error, B_base, recall = Rayuela.experiment_lsq_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h,
      niter, ilsiter, icmiter, randord, npert, knn, nsplits_train, nsplits_base, verbose)
    save_results_lsq("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)
  end

  # for trial = 1:ntrials
  #   sr_method = "SR_D"
  #   C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   C, B, R, train_error, B_base, recall = Rayuela.experiment_sr_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, knn, nsplits_train, nsplits_base, sr_method, verbose)
  #   save_results_lsq("./results/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)
  # end
  #
  # for trial = 1:ntrials
  #   sr_method = "SR_C"
  #   C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   C, B, R, train_error, B_base, recall = Rayuela.experiment_sr_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, knn, nsplits_train, nsplits_base, sr_method, verbose)
  #   save_results_lsq("./results/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)
  # end

  # GPU methods
  # nsplits_train = 1
  # nsplits_base  = 2
  # Rayuela.experiment_lsq_cuda(Xt, Xb, Xq, gt, m, h, niter, knn, nsplits_train, nsplits_base, verbose)
  # Rayuela.experiment_sr_cuda(Xt, Xb, Xq, gt, m, h, niter, knn, nsplits_train, nsplits_base, verbose)

  # GPU methods with random inputs
  # B = convert(Matrix{Int16}, rand(1:h, m, size(Xt,2)))
  # C = Vector{Matrix{Float32}}(m); for i=1:m; C[i]=zeros(Float32,d,h); end
  # R = eye(Float32, d)
  # Rayuela.experiment_lsq_cuda(Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, verbose)
  # Rayuela.experiment_sr_cuda( Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, verbose)

  # return C, B, train_error, B_base, recall
  # return C, B, R, train_error, B_base, recall
end

# === demos query base ===
function run_demos_query_base(
  dataset_name="labelme",
  ntrain::Integer=Int(2e3),
  m::Integer=8, h::Integer=256, niter::Integer=25, ntrials::Integer=10)

  nquery, nbase, knn = 0, 0, 0
  if dataset_name == "MNIST"
    nquery, nbase, knn = Int(10e3), Int(60e3), Int(1e3)
  elseif dataset_name == "labelme"
    nquery, nbase, knn = Int(2e3), Int(20019), Int(1e3)
  else
    error("dataset unknown")
  end

  verbose = true

  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  d, _ = size( Xt )

  # # (Semi-)orthogonal methods: PQ, OPQ, ChainQ
  # for trial = 1:ntrials
  #   C, B, train_error, recall = Rayuela.experiment_pq_query_base(Xt, Xq, gt, m, h, niter, knn, verbose)
  #   save_results_pq_query_base("./results/$(lowercase(dataset_name))/pq_m$(m)_it$(niter).h5", trial, C, B, train_error, recall)
  # end
  # for trial = 1:ntrials
  #   C, B, R, train_error, recall = Rayuela.experiment_opq_query_base(Xt, Xq, gt, m, h, niter, knn, verbose)
  #   save_results_opq_query_base("./results/$(lowercase(dataset_name))/opq_m$(m)_it$(niter).h5", trial, C, B, R, train_error, recall)
  # end

  # # Cheap non-orthogonal methods: RVQ, ERVQ
  # for trial = 1:ntrials
  #   C, B, train_error, recall = Rayuela.experiment_rvq_query_base(Xt,        Xq, gt, m-1, h, niter, knn, verbose)
  #   save_results_pq_query_base("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5",  trial, C, B, train_error, recall)
  # end
  # for trial = 1:ntrials
  #   C, B, _ = load_rvq("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   C, B, train_error, recall = Rayuela.experiment_ervq_query_base(Xt, B, C, Xq, gt, m-1, h, niter, knn, verbose)
  #   save_results_pq_query_base("./results/$(lowercase(dataset_name))/ervq_m$(m-1)_it$(niter).h5", trial, C, B, train_error, recall)
  # end

    # More expensive non-orthogonal methods
    # Rayuela.experiment_chainq(Xt, Xb, Xq, gt, m, h, niter, knn, verbose)
    # Rayuela.experiment_lsq(Xt, Xb, Xq, gt, m, h, niter, knn, verbose)
    # Rayuela.experiment_sr(Xt, Xb, Xq, gt, m, h, niter, knn, verbose)

  # # Precompute init for LSQ/SR
  # for trial = 1:ntrials
  #   C, B, R, train_error = Rayuela.train_opq(Xt, m-1, h, niter, "natural", verbose)
  #   save_results_opq_query_base("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, [0f0])
  # end
  for trial = 1:ntrials
    C, B, R, _ = load_chainq("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, chainq_error = Rayuela.train_chainq(Xt, m-1, h, R, B, C, niter, verbose)
    save_results_opq_query_base("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", trial, C, B, R, chainq_error, [0f0])
  end

  # nsplits_train = 1
  # for trial = 1:ntrials
  #   C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   # C, B, R, train_error, recall = Rayuela.experiment_lsq_cuda_query_base(Xt, B, C, R, Xq, gt, m-1, h, niter, knn, nsplits_train, verbose)
  #   C, B, R, train_error, recall = Rayuela.experiment_lsq_query_base(Xt, B, C, R, Xq, gt, m-1, h, niter, knn, verbose)
  #   save_results_lsq_query_base("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  # end

  # for trial = 1:ntrials
  #   sr_method = "SR_D"
  #   C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base( Xt, B, C, R, Xq, gt, m-1, h, niter, knn, nsplits_train, sr_method, verbose)
  #   save_results_lsq_query_base("./results/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  # end
  #
  # for trial = 1:ntrials
  #   sr_method = "SR_C"
  #   C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
  #   C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base( Xt, B, C, R, Xq, gt, m-1, h, niter, knn, nsplits_train, sr_method, verbose)
  #   save_results_lsq_query_base("./results/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  # end

  # GPU methods
  # nsplits_train = 1
  # (C, B, R, train_error, recall), opq_error = Rayuela.experiment_lsq_cuda_query_base(Xt, Xq, gt, m-1, h, niter, knn, nsplits_train, verbose)
  # save_results_lsq_query_base("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, opq_error, recall)

  # sr_method = "SR_D"
  # (C, B, R, train_error, recall), opq_error = Rayuela.experiment_sr_cuda_query_base( Xt, Xq, gt, m-1, h, niter, knn, nsplits_train, sr_method, verbose)
  # save_results_lsq_query_base("./results/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, opq_error, recall)
  #
  # sr_method = "SR_C"
  # (C, B, R, train_error, recall), opq_error = Rayuela.experiment_sr_cuda_query_base( Xt, Xq, gt, m-1, h, niter, knn, nsplits_train, sr_method, verbose)
  # save_results_lsq_query_base("./results/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, opq_error, recall)

  # GPU methods with random inputs
  # B = convert(Matrix{Int16}, rand(1:h, m, size(Xt,2)))
  # C = Vector{Matrix{Float32}}(m); for i=1:m; C[i]=zeros(Float32,d,h); end
  # R = eye(Float32, d)
  # Rayuela.experiment_lsq_cuda(Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, verbose)
  # Rayuela.experiment_sr_cuda( Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, verbose)

end

function high_recall_experiments(
  dataset_name="SIFT1M",
  method::String="lsq",
  ntrain::Integer=Int(1e5),
  m::Integer=8,
  h::Integer=256,
  niter::Integer=25)

  if !(method in ["lsq", "src", "srd"])
    error("Method $method unkonwn")
  end

  nquery, nbase, knn = Int(1e4), Int(1e6), Int(1e3)
  verbose = V = true
  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  d, _    = size( Xt )

  ilsiters = [1, 2, 4, 8, 16, 32, 64, 128, 256]
  icmiter = 4
  randord = true
  npert   = 4
  p       = 0.5f0
  nsplits_base = m <= 8 ? 2 : 4

  ntrials = 10
  for trial = 1:ntrials

    @show trial

    fname = "./results/$(lowercase(dataset_name))/$(method)_m$(m)_it$(niter).h5"
    C, B, R, chainq_error = load_chainq(fname, m, trial)
    norms_B, norms_C = get_norms_codebook(B, C)

    # === Encode the base set ===
    B_base = convert(Matrix{Int16}, rand(1:h, m, size(Xb,2)))

    # recalls = zeros(Float32, length(ilsiters), knn)
    Bs_base, _ = encode_icm_cuda(Xb, B_base, C, ilsiters, icmiter, npert, randord, nsplits_base, V)

    for (idx, ilsiter) in enumerate(ilsiters)
      B_base = Bs_base[idx]
      base_error = qerror(Xb, B_base, C)
      if V; @printf("Error in base is %e\n", base_error); end

      # Compute and quantize the database norms
      B_base_norms, db_norms_X = quantize_norms( B_base, C, norms_C )
      db_norms = vec( norms_C[ B_base_norms ] )

      if V; print("Querying m=$m, ilsiter=$ilsiter ... "); end
      @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, eye(Float32, d), knn)
      if V; println("done"); end

      recall = eval_recall(gt, idx, knn)

      println("Saving to $fname...")
      h5write(fname, "$trial/recall_$(ilsiter)", recall)
      # h5write(bpath, "$trial/B_Base_$(ilsiter)", convert(Matrix{UInt8}, B_base.-1))
    end
  end

end

# high_recall_experiments("Convnet1M", "lsq", Int(1e5), 15, 256, 25)
# for method in ["lsq", "src", "srd"], niter = [50, 100], m = [15]
# for method in ["srd"], niter = [25], m = [15]
  # high_recall_experiments("SIFT1M", method, Int(1e5), m, 256, niter)
  # high_recall_experiments("Deep1M", method, Int(1e5), m, 256, niter)
  # high_recall_experiments("Convnet1M", method, Int(1e5), m, 256, niter)
# end

for niter = [25]
  # Increase the second parameter to Int(1e5) to reproduce ECCV'18 results

  run_demos("SIFT1M", Int(1e4),  8, 256, niter)
  # run_demos("SIFT1M", Int(1e5), 16, 256, niter)
  # run_demos("Convnet1M", Int(1e5),   8, 256, niter)
  # run_demos("Convnet1M", Int(1e5),  16, 256, niter)
  # run_demos("Deep1M", Int(1e5),  8, 256, niter)
  # run_demos("Deep1M", Int(1e5), 16, 256, niter)
end

# run_demos_query_base("MNIST",   Int(60e3), 8,  256, 5, 1)
# run_demos_query_base("MNIST",   Int(60e3), 16, 256, 5, 1)
# run_demos_query_base("labelme", Int(20e3), 8,  256, 5, 1)
# run_demos_query_base("labelme", Int(20e3), 16, 256, 5, 1)

# run_demos
# for niter = [100]
#   # run_demos_query_base("labelme", Int(20e3), 8, 256, niter)
#   run_demos_query_base("labelme", Int(20e3), 16, 256, niter)
#   # run_demos_query_base("MNIST",   Int(60e3), 8, 256, niter)
#   run_demos_query_base("MNIST",   Int(60e3), 16, 256, niter)
# end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

pathToFMU = "https://github.com/ThummeTo/FMI.jl/raw/main/model/" * ENV["EXPORTINGTOOL"] * "/IO.fmu"

myFMU = fmi2Load(pathToFMU)
comp = fmi2Instantiate!(myFMU; loggingOn=false)
@test comp != 0

@test fmi2SetupExperiment(comp, 0.0) == 0

@test fmi2EnterInitializationMode(comp) == 0

realValueReferences = ["p_real", "u_real"]
integerValueReferences = ["p_integer", "u_integer"]
booleanValueReferences = ["p_boolean", "u_boolean"]
stringValueReferences = ["p_string", "p_string"]

#########################
# Testing Single Values #
#########################

rndReal = 100 * rand()
rndInteger = round(Integer, 100 * rand())
rndBoolean = rand() > 0.5
rndString = Random.randstring(12)

cacheReal = 0.0
cacheInteger = 0
cacheBoolean = false
cacheString = ""

@test fmi2SetReal(comp, realValueReferences[1], rndReal) == 0
@test fmi2GetReal(comp, realValueReferences[1]) == rndReal
@test fmi2SetReal(comp, realValueReferences[1], -rndReal) == 0
@test fmi2GetReal(comp, realValueReferences[1]) == -rndReal

@test fmi2SetInteger(comp, integerValueReferences[1], rndInteger) == 0
@test fmi2GetInteger(comp, integerValueReferences[1]) == rndInteger
@test fmi2SetInteger(comp, integerValueReferences[1], -rndInteger) == 0
@test fmi2GetInteger(comp, integerValueReferences[1]) == -rndInteger

@test fmi2SetBoolean(comp, booleanValueReferences[1], rndBoolean) == 0
@test fmi2GetBoolean(comp, booleanValueReferences[1]) == rndBoolean
@test fmi2SetBoolean(comp, booleanValueReferences[1], !rndBoolean) == 0
@test fmi2GetBoolean(comp, booleanValueReferences[1]) == !rndBoolean

@test fmi2SetString(comp, stringValueReferences[1], rndString) == 0
@test fmi2GetString(comp, stringValueReferences[1]) == rndString

if ENV["EXPORTINGTOOL"] != "OpenModelica/v1.17.0"
    fmi2Set(comp,
            [realValueReferences[1], integerValueReferences[1], booleanValueReferences[1], stringValueReferences[1]],
            [rndReal,                rndInteger,                rndBoolean,                rndString])
    @test fmi2Get(comp,
                  [realValueReferences[1], integerValueReferences[1], booleanValueReferences[1], stringValueReferences[1]]) ==
                  [rndReal,                rndInteger,                rndBoolean,                rndString]
end

##################
# Testing Arrays #
##################

rndReal = [100 * rand(), 100 * rand()]
rndInteger = [round(Integer, 100 * rand()), round(Integer, 100 * rand())]
rndBoolean = [(rand() > 0.5), (rand() > 0.5)]
tmp = Random.randstring(8)
rndString = [tmp, tmp]

cacheReal = [0.0, 0.0]
cacheInteger =  [fmi2Integer(0), fmi2Integer(0)]
cacheBoolean = [fmi2Boolean(false), fmi2Boolean(false)]
cacheString = [pointer(""), pointer("")]

@test fmi2SetReal(comp, realValueReferences, rndReal) == 0
@test fmi2GetReal(comp, realValueReferences) == rndReal
fmi2GetReal!(comp, realValueReferences, cacheReal)
@test cacheReal == rndReal
@test fmi2SetReal(comp, realValueReferences, -rndReal) == 0
@test fmi2GetReal(comp, realValueReferences) == -rndReal
fmi2GetReal!(comp, realValueReferences, cacheReal)
@test cacheReal == -rndReal

@test fmi2SetInteger(comp, integerValueReferences, rndInteger) == 0
@test fmi2GetInteger(comp, integerValueReferences) == rndInteger
fmi2GetInteger!(comp, integerValueReferences, cacheInteger)
@test cacheInteger == rndInteger
@test fmi2SetInteger(comp, integerValueReferences, -rndInteger) == 0
@test fmi2GetInteger(comp, integerValueReferences) == -rndInteger
fmi2GetInteger!(comp, integerValueReferences, cacheInteger)
@test cacheInteger == -rndInteger

@test fmi2SetBoolean(comp, booleanValueReferences, rndBoolean) == 0
@test fmi2GetBoolean(comp, booleanValueReferences) == rndBoolean
fmi2GetBoolean!(comp, booleanValueReferences, cacheBoolean)
@test cacheBoolean == rndBoolean
not_rndBoolean = collect(!b for b in rndBoolean)
@test fmi2SetBoolean(comp, booleanValueReferences, not_rndBoolean) == 0
@test fmi2GetBoolean(comp, booleanValueReferences) == not_rndBoolean
fmi2GetBoolean!(comp, booleanValueReferences, cacheBoolean)
@test cacheBoolean == not_rndBoolean

@test fmi2SetString(comp, stringValueReferences, rndString) == 0
@test fmi2GetString(comp, stringValueReferences) == rndString
fmi2GetString!(comp, stringValueReferences, cacheString)
@test unsafe_string.(cacheString) == rndString

# this is not suppoerted by OMEdit-FMUs in the repository
if ENV["EXPORTINGTOOL"] != "OpenModelica/v1.17.0"
    # Testing input/output derivatives
    dirs = fmi2GetRealOutputDerivatives(comp, ["y_real"], ones(Int, 1))
    @test dirs == -Inf # at this point, derivative is undefined
    @test fmi2SetRealInputDerivatives(comp, ["u_real"], ones(Int, 1), zeros(1)) == 0

    @test fmi2ExitInitializationMode(comp) == 0
    @test fmi2DoStep(comp, 0.1) == 0

    dirs = fmi2GetRealOutputDerivatives(comp, ["y_real"], ones(Int, 1))
    @test dirs == 0.0
else
    @test fmi2ExitInitializationMode(comp) == 0
    @test fmi2DoStep(comp, 0.1) == 0
end

############
# Clean up #
############

fmi2Unload(myFMU)

ext: jl


text:
using Extremes, DataFrames, Distributions, Gadfly, Dates. Mamba

data = load("wooster")

x = collect(Date(1983,1,1):Day(1):Date(1987,12,31))

data[!,:Date] = x
select!(data, [:Date, :Temperature])

data

plot(data, x=:Date, y=:Temperature, Geom.point)

df = copy(data)

df[!,:Temperature] = -data[:,:Temperature]

filter!(row -> month(row.Date) ∈ (1,2,11,12), df)

plot(df, x=:Date, y=:Temperature, Geom.point)

y = df[:,:Temperature]

threshold =  -10

cluster = getcluster(y, threshold, runlength=4)

z = maximum.(cluster) .- threshold

fm = gpfit(z)

r = returnlevel(fm, threshold, size(df,1), 120, 100, .95)

r.value[]
r.cint[]

data = load("dowjones")

X = data[:,:Index]

X̃ = 100*(log.(X[2:end]) - log.(X[1:end-1]))

df = DataFrame(Date = data[2:end, :Date], Variation = X̃)

plot(df, x=:Date, y=:Variation, Geom.line)

threshold = 2.0

cluster = getcluster(df[:,:Variation], threshold, runlength=4)

z = maximum.(cluster) .- threshold

fm = gpfit(z)

data = load("fremantle")

fm = gevfitbayes(data, :SeaLevel, locationcovid = [:Year, :SOI])

p = plot(fm.sim[:,1,1])

p[2]

p = plot(fm.sim)

hstack([p[1], p[2]])

# Lis all files in a directory
filenames = String[]

for (root, dirs, files) in walkdir("src/validationplots/")
    for file in files
        push!(filenames, joinpath(root, file))
        # println(joinpath(root, file)) # path to files
    end
end


text: struct println a; div; b end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test
using LavoisierCore

model = GERG2008(:C2,:C3)
x0 = [0.5,0.5]
T0 = 300
@testset "initial" begin
    @test length(model) == 2
end


text:
"""
    const Point{D,T} = VectorValue{D,T}

Type representing a point of D dimensions with coordinates of type T.
Fields are evaluated at vectors of `Point` objects.
"""
const Point{D,T} = VectorValue{D,T}

"""
    abstract type Field <: Map

Abstract type representing a physical (scalar, vector, or tensor) field. The
domain is a `Point` and the range a scalar (i.e., a sub-type of Julia `Number`),
a `VectorValue`, or a `TensorValue`.

These different cases are distinguished by the return value obtained when evaluating them. E.g.,
a physical field returns a vector of values when evaluated at a vector of points, and a basis of `nf` fields
returns a 2d matrix (`np` x `nf`) when evaluated at a vector of `np` points.

The following functions (i.e., the `Map` API) need to be overloaded:

- [`evaluate!(cache,f,x)`](@ref)
- [`return_cache(f,x)`](@ref)

and optionally

- [`return_type(f,x)`](@ref)

A `Field` can also provide its gradient if the following function is implemented
- [`gradient(f)`](@ref)

Higher derivatives can be obtained if the resulting object also implements this method.

The next paragraph is out-of-date:

Moreover, if the [`gradient(f)`](@ref) is not provided, a default implementation that uses the
following functions will be used.

- [`evaluate_gradient!(cache,f,x)`](@ref)
- [`return_gradient_cache(f,x)`](@ref)

Higher order derivatives require the implementation of

- [`evaluate_hessian!(cache,f,x)`](@ref)
- [`return_hessian_cache(f,x)`](@ref)

These four methods are only designed to be called by the default implementation of [`field_gradient(f)`](@ref) and thus
cannot be assumed that they are available for an arbitrary field. For this reason, these functions are not
exported. The general way of evaluating a gradient of a field is to
build the gradient with [`gradient(f)`](@ref) and evaluating the resulting object. For evaluating
the hessian, use two times `gradient`.

The interface can be tested with

- [`test_field`](@ref)

For performance, the user can also consider a _vectorised_ version of the
`Field` API that evaluates the field in a vector of points (instead of only one
point). E.g., the `evaluate!` function for a vector of points returns a vector
of scalar, vector or tensor values.

"""
abstract type Field <: Map end

evaluate!(c,f::Field,x::Point) = @abstractmethod

# Differentiation

function gradient end
const ∇ = gradient
∇∇(f) = gradient(gradient(f))

gradient(f,::Val{1}) = ∇(f)
gradient(f,::Val{2}) = ∇∇(f)

evaluate!(cache,::Broadcasting{typeof(∇)},a::Field) = ∇(a)
evaluate!(cache,::Broadcasting{typeof(∇∇)},a::Field) = ∇∇(a)
lazy_map(::Broadcasting{typeof(∇)},a::AbstractArray{<:Field}) = lazy_map(∇,a)
lazy_map(::Broadcasting{typeof(∇∇)},a::AbstractArray{<:Field}) = lazy_map(∇∇,a)

push_∇(∇a::Field,ϕ::Field) = pinvJt(∇(ϕ))⋅∇a

function pinvJt(Jt::MultiValue{Tuple{D,D}}) where D
  inv(Jt)
end

function pinvJt(Jt::MultiValue{Tuple{D1,D2}}) where {D1,D2}
  @check D1 < D2
  J = transpose(Jt)
  transpose(inv(Jt⋅J)⋅Jt)
end

function push_∇∇(∇∇a::Field,ϕ::Field)
  @notimplemented """\n
  Second order derivatives of quantities defined in the reference domain not implemented yet.

  This is a feature that we want to have at some point in Gridap.
  If you are ready to help with this implementation, please contact the
  Gridap administrators.
  """
end

"""
    gradient_type(::Type{T},x::Point) where T
"""
function gradient_type(::Type{T},x::Point) where T
  typeof(outer(zero(x),zero(T)))
end

"""
Type that represents the gradient of a field. The wrapped field must
implement `evaluate_gradient!` and `return_gradient_cache` for this gradient
to work.

N is how many times the gradient is applied
"""
struct FieldGradient{N,F} <: Field
  object::F
  FieldGradient{N}(object::F) where {N,F} = new{N,F}(object)
end

gradient(f::Field) = FieldGradient{1}(f)
gradient(f::FieldGradient{N}) where N = FieldGradient{N+1}(f.object)

testargs(f::FieldGradient,x::Point) = testargs(f.object,x)
return_value(f::FieldGradient,x::Point) = evaluate(f,testargs(f,x)...)
return_cache(f::FieldGradient,x::Point) = nothing
evaluate!(cache,f::FieldGradient,x::Point) = @abstractmethod
testvalue(::Type{FieldGradient{N,T}}) where {N,T} = FieldGradient{N}(testvalue(T))

# Default methods for arrays of points

function testargs(f::Field,x::AbstractArray{<:Point})
  y = copy(x)
  broadcast!(xi->first(testargs(f,xi)),y,x)
  (y,)
end

function return_cache(f::Field,x::AbstractArray{<:Point})
  T = return_type(f,testitem(x))
  s = size(x)
  ab = zeros(T,s)
  cb = CachedArray(ab)
  cf = return_cache(f,testitem(x))
  cb, cf
end

function evaluate!(c,f::Field,x::AbstractArray{<:Point})
  cb, cf = c
  sx = size(x)
  setsize!(cb,sx)
  r = cb.array
  for i in eachindex(x)
    @inbounds r[i] = evaluate!(cf,f,x[i])
  end
  r
end

# GenericField

"""
A wrapper for objects that can act as fields, e.g., functions which implement the `Field` API.
"""
struct GenericField{T} <: Field
  object::T
end

#Field(f) = GenericField(f)
GenericField(f::Field) = f

testargs(a::GenericField,x::Point) = testargs(a.object,x)
return_value(a::GenericField,x::Point) = return_value(a.object,x)
return_cache(a::GenericField,x::Point) = return_cache(a.object,x)
evaluate!(cache,a::GenericField,x::Point) = evaluate!(cache,a.object,x)

function return_cache(f::FieldGradient{N,<:GenericField},x::Point) where N
  return_cache(FieldGradient{N}(f.object.object),x)
end

function evaluate!(c,f::FieldGradient{N,<:GenericField},x::Point) where N
  evaluate!(c,FieldGradient{N}(f.object.object),x)
end

# Make Field behave like a collection

Base.length(::Field) = 1
Base.size(::Field) = ()
Base.axes(::Field) = ()
Base.IteratorSize(::Type{<:Field}) = Base.HasShape{0}()
Base.eltype(::Type{T}) where T<:Field = T
Base.iterate(a::Field) = (a,nothing)
Base.iterate(a::Field,::Nothing) = nothing
Base.getindex(a::Field,i::Integer) =  (@check i == 1; a)
testitem(a::Field) = a

# Zero field

Base.zero(a::Field) = ZeroField(a)

"""
It represents `0.0*f` for a field `f`.
"""
struct ZeroField{F} <: Field
  field::F
end

return_cache(z::ZeroField,x::Point) = zero(return_type(z.field,x))
evaluate!(cache,z::ZeroField,x::Point) = cache
testvalue(::Type{ZeroField{F}}) where F = ZeroField(testvalue(F))

function return_cache(z::ZeroField,x::AbstractArray{<:Point})
  E = return_type(z.field,testitem(x))
  c = zeros(E,size(x))
  CachedArray(c)
end

function evaluate!(c,f::ZeroField,x::AbstractArray{<:Point})
  nx = size(x)
  if size(c) != nx
    setsize!(c,nx)
    fill!(c.array,zero(eltype(c)))
  end
  c.array
end

gradient(z::ZeroField) = ZeroField(gradient(z.field))

# Make Number behave like Field
#

# Number itself does not implement the Field interface since Number objects are not callable in Julia.
# Thus, wrapping a Number in a GenericField does not makes sense since the wrapped object
# is assumed to implement the Field interface.
# I think it is conceptually better to have ConstantField as struct otherwise we break the invariant
# "for any object wrapped in a GenericField we can assume that it implements the Field interface"
struct ConstantField{T<:Number} <: Field
  value::T
end

constant_field(a) = ConstantField(a)

Base.zero(::Type{ConstantField{T}}) where T = ConstantField(zero(T))

function evaluate!(c,f::ConstantField,x::Point)
  f.value
end

function return_cache(f::ConstantField,x::AbstractArray{<:Point})
  nx = size(x)
  c = fill(f.value,nx)
  CachedArray(c)
end

function evaluate!(c,f::ConstantField,x::AbstractArray{<:Point})
  nx = size(x)
  # This optimization is a bug if we include several ConstantField with different states
  # in the same array and we try to reuse cache between them.
  #if size(c) != nx
    setsize!(c,nx)
    fill!(c.array,f.value)
  #end
  c.array
end

function return_cache(f::FieldGradient{N,<:ConstantField},x::Point) where N
  gradient(f.object.value,Val(N))(x)
end

evaluate!(c,f::FieldGradient{N,<:ConstantField},x::Point) where N = c

function return_cache(f::FieldGradient{N,<:ConstantField},x::AbstractArray{<:Point}) where N
  CachedArray(gradient(f.object.value,Val(N)).(x))
end

function evaluate!(c,f::FieldGradient{N,<:ConstantField},x::AbstractArray{<:Point}) where N
  nx = size(x)
  if size(c) != nx
    setsize!(c,nx)
    fill!(c.array,zero(eltype(c)))
  end
  c.array
end

function lazy_map(::Operation{typeof(inv)},a::LazyArray{<:Fill{typeof(constant_field)}})
  v = a.args[1]
  vinv = lazy_map(inv,v)
  lazy_map(constant_field,vinv)
end

## Make Function behave like Field

return_cache(f::FieldGradient{N,<:Function},x::Point) where N = gradient(f.object,Val(N))
evaluate!(c,f::FieldGradient{N,<:Function},x::Point) where N = c(x)

# Operations

"""
A `Field` that is obtained as a given operation over a tuple of fields.
"""
struct OperationField{O,F} <: Field
  op::O
  fields::F
end

function return_value(c::OperationField,x::Point)
  fx = map(f -> return_value(f,x),c.fields)
  return_value(c.op,fx...)
end

function return_cache(c::OperationField,x::Point)
  cl = map(fi -> return_cache(fi,x),c.fields)
  lx = map(fi -> return_value(fi,x),c.fields)
  ck = return_cache(c.op,lx...)
  ck, cl
end

function evaluate!(cache,c::OperationField,x::Point)
  ck, cf = cache
  lx = map((ci,fi) -> evaluate!(ci,fi,x),cf,c.fields)
  evaluate!(ck,c.op,lx...)
end

function return_cache(c::OperationField,x::AbstractArray{<:Point})
  cf = map(fi -> return_cache(fi,x),c.fields)
  lx = map((ci,fi) -> evaluate!(ci,fi,x),cf,c.fields)
  ck = return_cache(c.op,map(testitem,lx)...)
  r = c.op.(lx...)
  ca = CachedArray(r)
  ca, ck, cf
end

function evaluate!(cache,c::OperationField,x::AbstractArray{<:Point})
  ca, ck, cf = cache
  sx = size(x)
  setsize!(ca,sx)
  lx = map((ci,fi) -> evaluate!(ci,fi,x),cf,c.fields)
  r = ca.array
  for i in eachindex(x)
    @inbounds r[i] = evaluate!(ck,c.op,map(lxi -> lxi[i], lx)...)
  end
  r
end

evaluate!(cache,op::Operation,x::Field...) = OperationField(op.op,x)
return_value(op::Broadcasting{<:Operation},x::Field...) = OperationField(op.f.op,x)
evaluate!(cache,op::Broadcasting{<:Operation},x::Field...) = OperationField(op.f.op,x)

# Define some well known operations

for op in (:+,:-,:*,:/,:⋅,:⊙,:⊗,:inv,:det,:meas,:pinvJt,:tr,:grad2curl,:symmetric_part,:transpose)
  @eval ($op)(a::Field...) = Operation($op)(a...)
end

transpose(f::Field) = f

for op in (:+,:-,:*,:/,:⋅,:⊙,:⊗)
  @eval ($op)(a::Field,b::Number) = Operation($op)(a,ConstantField(b))
  @eval ($op)(a::Number,b::Field) = Operation($op)(ConstantField(a),b)
end

#*(A::Number, B::Field) = ConstantField(A)*B
#*(A::Field, B::Number) = A*ConstantField(B)
#⋅(A::Number, B::Field) = ConstantField(A)⋅B
#⋅(A::Field, B::Number) = A⋅ConstantField(B)

#*(A::Function, B::Field) = GenericField(A)*B
#*(A::Field, B::Function) = GenericField(B)*A

# Gradient of the sum
for op in (:+,:-)
  @eval begin
    function gradient(a::OperationField{typeof($op)})
      f = a.fields
      g = map( gradient, f)
      $op(g...)
    end
  end
end

# Gradient of the product

function product_rule(fun,f1,f2,∇f1,∇f2)
  msg = "Product rule not implemented for product $fun between types $(typeof(f1)) and $(typeof(f2))"
  @notimplemented msg
end

function product_rule(::typeof(*),f1::Real,f2::Real,∇f1,∇f2)
  ∇f1*f2 + f1*∇f2
end

for op in (:*,:⋅)
  @eval begin
     function product_rule(::typeof($op),f1::Real,f2::VectorValue,∇f1,∇f2)
       ∇f1⊗f2 + ∇f2*f1
     end

     function product_rule(::typeof($op),f1::VectorValue,f2::Real,∇f1,∇f2)
       product_rule(*,f2,f1,∇f2,∇f1)
     end
  end
end

function product_rule(::typeof(⋅),f1::VectorValue,f2::VectorValue,∇f1,∇f2)
  ∇f1⋅f2 + ∇f2⋅f1
end

function product_rule(::typeof(⋅),f1::TensorValue,f2::VectorValue,∇f1,∇f2)
  ∇f1⋅f2 + ∇f2⋅transpose(f1)
end

for op in (:*,:⋅,:⊙,:⊗)
  @eval begin
    function gradient(a::OperationField{typeof($op)})
      f = a.fields
      @notimplementedif length(f) != 2
      f1, f2 = f
      g1, g2 = map(gradient, f)
      k(F1,F2,G1,G2) = product_rule($op,F1,F2,G1,G2)
      Operation(k)(f1,f2,g1,g2)
    end
  end
end

# Chain rule
function gradient(f::OperationField{<:Field})
  a = f.op
  @notimplementedif length(f.fields) != 1
  b, = f.fields
  x = ∇(a)∘b
  y = ∇(b)
  y⋅x
end

# Composition

"""
    f∘g

It returns the composition of two fields, which is just `Operation(f)(g)`
"""
Base.:∘(f::Field,g::Field) = Operation(f)(g)
evaluate!(cache,::Broadcasting{typeof(∘)},f::Field,g::Field) = f∘g

# Integration

"""
Integration of a given field in the "physical" space
"""
function integrate(a::Field,x::AbstractVector{<:Point},w::AbstractVector{<:Real})
  cache = return_cache(integrate,a,x,w)
  evaluate!(cache,integrate,a,x,w)
end

"""
Integration of a given field in the "reference" space
"""
function integrate(a::Field,q::AbstractVector{<:Point},w::AbstractVector{<:Real},j::Field)
  cache = return_cache(integrate,a,q,w,j)
  evaluate!(cache,integrate,a,q,w,j)
end

function return_cache(::typeof(integrate),a,x,w)
  ca = return_cache(a,x)
  ax = return_value(a,x)
  ck = return_cache(IntegrationMap(),ax,w)
  ca, ck
end

function evaluate!(cache,::typeof(integrate),a,x,w)
  ca, ck = cache
  ax = evaluate!(ca,a,x)
  evaluate!(ck,IntegrationMap(),ax,w)
end

function return_cache(::typeof(integrate),a,q,w,j)
  ca = return_cache(a,q)
  cj = return_cache(j,q)
  aq = return_value(a,q)
  jq = return_value(j,q)
  ck = return_cache(IntegrationMap(),aq,w,jq)
  ca, cj, ck
end

function evaluate!(cache,::typeof(integrate),a,q,w,j)
  ca, cj, ck = cache
  aq = evaluate!(ca,a,q)
  jq = evaluate!(cj,j,q)
  evaluate!(ck,IntegrationMap(),aq,w,jq)
end

struct IntegrationMap <: Map end

function evaluate!(cache,k::IntegrationMap,ax::AbstractVector,w)
  T = typeof( testitem(ax)*testitem(w) + testitem(ax)*testitem(w) )
  z = zero(T)
  r = z
  @check length(ax) == length(w)
  @inbounds for i in eachindex(ax)
    r += ax[i]*w[i]
  end
  r
end

function evaluate!(cache,k::IntegrationMap,aq::AbstractVector,w,jq::AbstractVector)
  T = typeof( testitem(aq)*testitem(w)*meas(testitem(jq)) + testitem(aq)*testitem(w)*meas(testitem(jq)) )
  z = zero(T)
  @check length(aq) == length(w)
  @check length(aq) == length(jq)
  @inbounds for i in eachindex(aq)
    z += aq[i]*w[i]*meas(jq[i])
  end
  z
end

function return_cache(k::IntegrationMap,ax::AbstractArray,w)
  T = typeof( testitem(ax)*testitem(w) + testitem(ax)*testitem(w) )
  r = zeros(T,size(ax)[2:end])
  CachedArray(r)
end

function evaluate!(cache,k::IntegrationMap,ax::AbstractArray,w)
  setsize!(cache,size(ax)[2:end])
  r = cache.array
  @check size(ax,1) == length(w)
  @inbounds for j in CartesianIndices(r)
    rj = zero(eltype(r))
    for p in 1:length(w)
      rj += ax[p,j]*w[p]
    end
    r[j] = rj
  end
  r
end

function return_value(k::IntegrationMap,aq::AbstractArray,w,jq::AbstractVector)
  if size(aq,1) == length(w) && size(aq,1) == length(jq)
    evaluate(k,aq,w,jq)
  else
    c = return_cache(k,aq,w,jq)
    c.array
  end
end

function return_cache(k::IntegrationMap,aq::AbstractArray,w,jq::AbstractVector)
  T = typeof( testitem(aq)*testitem(w)*meas(testitem(jq)) + testitem(aq)*testitem(w)*meas(testitem(jq)) )
  r = zeros(T,size(aq)[2:end])
  CachedArray(r)
end

function evaluate!(cache,k::IntegrationMap,aq::AbstractArray,w,jq::AbstractVector)
  setsize!(cache,size(aq)[2:end])
  r = cache.array
  @check size(aq,1) == length(w) || size(aq,1) == 0
  @check size(aq,1) == length(jq) || size(aq,1) == 0
  fill!(r,zero(eltype(r)))
  cis = CartesianIndices(r)
  @inbounds for p in 1:length(w)
    dV = meas(jq[p])*w[p]
    for j in cis
      r[j] += aq[p,j]*dV
    end
  end
  r
end

function return_cache(k::IntegrationMap,aq::AbstractArray{S,3} where S,w,jq::AbstractVector)
  T = typeof( testitem(aq)*testitem(w)*meas(testitem(jq)) + testitem(aq)*testitem(w)*meas(testitem(jq)) )
  r = zeros(T,size(aq)[2:end])
  s = zeros(typeof(meas(testitem(jq))),length(jq))
  CachedArray(r), CachedArray(s)
end

function evaluate!(cache,k::IntegrationMap,aq::AbstractArray{S,3} where S, w,jq::AbstractVector)
  cache_r, cache_s = cache
  np, ni, nj = size(aq)
  setsize!(cache_r,(ni,nj))
  setsize!(cache_s,(np,))
  r = cache_r.array
  dV = cache_s.array
  @check np == length(w) || np == 0
  @check np == length(jq) || np == 0
  @inbounds for p in 1:np
    dV[p] = meas(jq[p])*w[p]
  end
  #fill!(r,zero(eltype(r)))
  @inbounds for j in 1:nj
    for i in 1:ni
      rij = zero(eltype(aq))
      for p in 1:np
        rij += aq[p,i,j]*dV[p]
      end
      r[i,j] = rij
    end
  end
  r
end

function evaluate!(cache,k::IntegrationMap,aq::AbstractMatrix, w,jq::AbstractVector)
  np, ni = size(aq)
  setsize!(cache,(ni,))
  r = cache.array
  @check np == length(w) || np == 0
  @check np == length(jq) || np == 0
  fill!(r,zero(eltype(r)))
  @inbounds for p in 1:np
    dV = meas(jq[p])*w[p]
    for i in 1:ni
      r[i] += aq[p,i]*dV
    end
  end
  r
end

# Testers

"""
    test_field(
      f::Union{Field,AbstractArray{<:Field}},
      x,
      v,
      cmp=(==);
      grad=nothing,
      gradgrad=nothing)

Function used to test the field interface. `v` is an array containing the expected
result of evaluating the field `f` at the point or vector of points `x`. The comparison is performed using
the `cmp` function. For fields objects that support the `gradient` function, the keyword
argument `grad` can be used. It should contain the result of evaluating `gradient(f)` at x.
Idem for `gradgrad`. The checks are performed with the `@test` macro.
"""
function test_field(f::Field, x, v, cmp=(==); grad=nothing, gradgrad=nothing)
  test_map(v,f,x;cmp=cmp)
  if grad != nothing
    test_map(grad,∇(f),x;cmp=cmp)
  end
  if gradgrad != nothing
    test_map(gradgrad,∇∇(f),x;cmp=cmp)
  end
  true
end

struct VoidFieldMap <: Map
  isvoid::Bool
end

Arrays.evaluate!(cache,k::VoidFieldMap,b) = VoidField(b,k.isvoid)

struct VoidField{F} <: Field
  field::F
  isvoid::Bool
end

function return_cache(f::VoidField,x::Point)
  c = return_cache(f.field,x)
  fx = evaluate!(c,f.field,x)
  c, zero(fx)
end

function evaluate!(cache,f::VoidField,x::Point)
  c,z = cache
  if f.isvoid
    z
  else
    evaluate!(c,f.field,x)
  end
end

function return_cache(f::VoidField,x::AbstractVector{<:Point})
  c = return_cache(f.field,x)
  fx = evaluate!(c,f.field,x)
  z = similar(fx)
  c, CachedArray(z)
end

function evaluate!(cache,f::VoidField,x::AbstractVector{<:Point})
  c,z = cache
  if f.isvoid
    setsize!(z,size(x))
    fill!(z,zero(eltype(z)))
    z.array
  else
    evaluate!(c,f.field,x)
  end
end

testvalue(::Type{VoidField{F}}) where F = VoidField(testvalue(F),false)

gradient(z::VoidField) = VoidField(gradient(z.field),z.isvoid)

function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{VoidFieldMap}},x::AbstractArray)
  p = a.maps.value
  @notimplementedif p.isvoid
  lazy_map(evaluate,a.args[1],x)
end

struct VoidBasisMap <: Map
  isvoid::Bool
end

Arrays.evaluate!(cache,k::VoidBasisMap,b) = VoidBasis(b,k.isvoid)

struct VoidBasis{T,N,A} <: AbstractArray{T,N}
  basis::A
  isvoid::Bool
  function VoidBasis(basis::AbstractArray{T,N},isvoid::Bool) where {T,N}
    new{T,N,typeof(basis)}(basis,isvoid)
  end
end

function Base.size(a::VoidBasis)
  if a.isvoid
    0 .* size(a.basis)
  else
    size(a.basis)
  end
end

Base.IndexStyle(::Type{<:VoidBasis}) = IndexLinear()

function Base.getindex(a::VoidBasis,i::Integer)
  if a.isvoid
    @unreachable "Unable to access 0-length array"
  else
    a.basis[i]
  end
end

Arrays.testitem(a::VoidBasis) = testitem(a.basis)

function _zero_size(a::VoidBasis{T,1} where T)
  (0,)
end

function _zero_size(a::VoidBasis{T,2} where T)
  @check size(a,1) == 1
  (1,0)
end

function Fields.return_cache(a::VoidBasis,x::Point)
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  zs = _zero_size(a)
  r = similar(bx,zs)
  cb,r
end

function Fields.return_cache(a::VoidBasis,x::Field)
  @notimplementedif ndims(a) != 1
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  r = similar(bx,(0,))
  cb,r
end

function Fields.return_cache(a::VoidBasis,x::AbstractVector{<:Point})
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  zs = _zero_size(a)
  r = similar(bx,(length(x),zs...))
  cb,r
end

function Fields.return_cache(a::VoidBasis,v::AbstractVector{<:Field})
  @notimplementedif ndims(a) != 1
  cb = return_cache(a.basis,v)
  bx = return_value(a.basis,v)
  r = similar(bx,(0,length(v)))
  cb,r
end

for T in (:Point,:Field,:(AbstractVector{<:Point}),:(AbstractVector{<:Field}))
  @eval begin

    function Fields.evaluate!(cache,a::VoidBasis,x::$T)
      cb, r = cache
      if a.isvoid
        r
      else
        evaluate!(cb,a.basis,x)
      end
    end

  end
end

function Fields.evaluate!(cache,k::Broadcasting{typeof(∇)},a::VoidBasis)
  VoidBasis(k(a.basis),a.isvoid)
end

function Fields.evaluate!(cache,k::Broadcasting{typeof(∇∇)},a::VoidBasis)
  VoidBasis(k(a.basis),a.isvoid)
end

function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{VoidBasisMap}},x::AbstractArray)
  p = a.maps.value
  @notimplementedif p.isvoid
  lazy_map(evaluate,a.args[1],x)
end


text:
using ODE

T = 1.0
Q = 0.1
dtsampling = 0.0125
dt=0.0025
nsteps = 50000 #If lyapunov_exponent function is chosen, total time = nsteps*dt
nsimulations = 50 #If trajectory function is chosen, total time = nsteps*dt*nsimulations.
thermo = Thermostat("logistic", z-> exp(z/Q)/(Q*( 1 +exp(z/Q))^2.))
potential = Potential("quartic", x->x^4./4.)
integrator = Integrator("RK4", ode4)

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
struct PairFeatureDescriptor <: AbstractPairFeatureDescriptor
    name::String
    encode_f::Any
    decode_f::Any
    length::Int # maybe, maybe not (does constrain/assume vector encoding)
    # probably needs some other stuff...
end


text:
import SimpleMolecule.AbstractMolecularContainer
import SimpleMolecule.clone

function test_clone{T<:AbstractMolecularContainer}(
    x::T,
    new_props::Dict,
    solution::T,
    msg::AbstractString=""
    )
    println("----------------------------------------")
    print_with_color(:blue, "Test clone()\n\n")

    if msg != ""
        print_with_color(:blue, msg * "\n")
    end

    answer = clone(x, new_props)

    @test answer == solution

    print_with_color(:green, "\nVERIFIED! clone()\n")
    println("----------------------------------------\n\n")
end


text:
using FermiCG, NPZ, LinearAlgebra, Plots
using Clustering, SpectralClustering
#load integrals from disk
ints = InCoreInts(
    npzread("integrals_h0.npy"),
    npzread("integrals_h1.npy"),
    npzread("integrals_h2.npy")
);
C = npzread("mo_coeffs_act.npy")
S = npzread("overlap_mat.npy");
D = npzread("density_mat.npy");

function cluster_orbitals(A, n_clusters)
    """
    A is adjacency matrix
    """

    L = Diagonal([sum(A[i,:]) for i in 1:size(A,1)]) - A;
    F = eigen(L);
    perm = sortperm(F.values, by=abs)
    F.values .= F.values[perm]
    F.vectors .= F.vectors[:,perm]

    print(" Laplacian Eigenvalues\n")
    for (fi_idx, fi) in enumerate(F.values)
        @printf(" %4i %12.8f\n", fi_idx, fi)
    end

    p1 = plot(F.values)

    A = abs.(F.vectors[:,1:n_clusters]*F.vectors[:,1:n_clusters]')
    p2 = heatmap(abs.(A), aspect_ratio=:equal, title="Projector", yflip = true)

    clustering = clusterize(KMeansClusterizer(n_clusters), A).assignments;

    perm = sortperm(clustering)
    p3 = heatmap(abs.(A[perm,perm]), aspect_ratio=:equal, title="Projector", yflip = true)

    clusters = [ [] for i in 1:n_clusters]
    for (idx,i) in enumerate(clustering)
        push!(clusters[i],idx)
    end

    fig = plot(p1, p2, p3, layout = @layout [a; b c])
    savefig(fig, "clustering_projector.pdf") # save the most recent fig as fn

    return perm, clusters
end

n_clusters = 2
adjacency = abs.(C'*D*C)
perm, cluster_list = cluster_orbitals(adjacency, n_clusters)

p1 = heatmap(adjacency, aspect_ratio=:equal, title="Adjacency", yflip = true)
p2 = heatmap(adjacency[perm, perm], aspect_ratio=:equal, title="Adjacency (sorted)", yflip = true)

fig = plot(p1, p2, layout = @layout [a b])
savefig(fig, "clustering_adjacency.pdf") # save the most recent fig as fn

npzwrite("clustering_perm.npy", perm)

max_stars_repo_name: nmayhall-vt/FermiCG-data

max_issues_repo_name: nmayhall-vt/FermiCG-data

max_forks_repo_name: nmayhall-vt/FermiCG-data


text:
#!/usr/bin/env julia

#
# Hello World client in Julia
# Connects REQ socket to tcp://localhost:5555
# Sends "Hello" to server, expects "World" back
#

using ZMQ

context = Context()

# Socket to talk to server
println("Connecting to hello world server...")
socket = Socket(context, REQ)
ZMQ.connect(socket, "tcp://localhost:5555")

for request in [1:10]
    println("Sending request $request ...")
    ZMQ.send(socket, "Hello")

    # Get the reply.
    message = bytestring(ZMQ.recv(socket))
    println("Received reply $request [ $message ]")
end

# Making a clean exit.
ZMQ.close(socket)
ZMQ.close(context)

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
# Create example grid data
x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
# Wireframe plot
wireframe(x, y, z)


text:
function v_analytic(m::Model,p::Param,id,it)
    vf = m.v[id,it]
    c = m.c[id,it]
    if it == p.nT
        cons = m.avec
        cash = m.avec
        convert(Point,cash,vfun(id,it,cons,cash,vf,p))
    else
        cons = scaleGrid(p.cfloor_plot,gety(c.env)[2],p.k,logorder = 1)
        cash = scaleGrid(p.a_low,getx(vf.env)[2],p.k,logorder = 1)
        # deleteat!(cons,length(cons))
        # deleteat!(cash,length(cash))
        pts = convert(Point,cash,vfun(id,it,cons,cash,vf,p))
        vcat(pts, vf.env.v[2:end])  # connect at second point
    end
end

function v_analytic(m::BModel,p::Param,id,iy,it)
    vf = m.v[id,iy,it]
    c = m.c[id,iy,it]
    cons = scaleGrid(p.cfloor_plot,gety(c.env)[2],p.k,logorder = 1)
    cash = scaleGrid(p.nT == it ? p.a_lowT : p.a_low,getx(vf.env)[2],p.k,logorder = 1)
    # deleteat!(cons,length(cons))
    # deleteat!(cash,length(cash))
    pts = convert(Point,cash,vfun(x->u(x,id==2,p),it,cons,cash,vf,p))
    vcat(pts, vf.env.v[2:end])  # connect at second point
    # vf.env.v[2:end]  # connect at second point
end

function v_analytic(m::Model,p::Param,id,iy,it)
    vf = m.v[id,iy,it]
    c = m.c[id,iy,it]
    cons = scaleGrid(p.cfloor_plot,gety(c.env)[2],p.k,logorder = 1)
    cash = scaleGrid(p.a_low,getx(vf.env)[2],p.k,logorder = 1)
    # deleteat!(cons,length(cons))
    # deleteat!(cash,length(cash))
    pts = convert(Point,cash,vfun(id,it,cons,cash,vf,p))
    vcat(pts, vf.env.v[2:end])  # connect at second point
    # vf.env.v[2:end]  # connect at second point
end

function plot_s(s::Simulation)

    # inc, cons, w
    py = plot(s.inc',leg = false, title = "income",yticks = [0,5,10,15],ylims = (-1,19))
    if s.p.ν > 0
        pc = plot(s.cons[:,1:(s.p.nT-1)]',leg = false, title = "consumption",ylims = (0,7))
        pw1 = plot(s.w1[:,1:(s.p.nT-1)]',leg = false, title = "w1",ylims = (0,maximum(s.p.a_high) + 0.1*maximum(s.p.a_high)))

    else
        pc = plot(s.cons',leg = false, title = "consumption",ylims = (0,7))
        pw1 = plot(s.w1',leg = false, title = "w1",ylims = (0,maximum(s.p.a_high) + 0.1*maximum(s.p.a_high)))
    end
    i_retires = findall(s.ret_age .> 0)
    scatter!(pc, s.ret_age[i_retires], [s.cons[i,s.ret_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))
    pw0 = plot(s.w0',leg = false, title = "w0",ylims = (0,maximum(s.p.a_high) + 0.1*maximum(s.p.a_high)))
    scatter!(pw0, s.ret_age[i_retires], [s.w0[i,s.ret_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))

    ppr = plot(s.prob_work',leg = false, title = "p(work)")
    scatter!(ppr, s.ret_age[i_retires], [s.prob_work[i,s.ret_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))

    pbeq = if s.p.ν > 0
        histogram(s.w0[:,s.p.nT], title = "Bequests", bins = range(0,stop = 30.0, length = 50), normalize = false, legend = false)
    else
        histogram(s.w1[:,s.p.nT], title = "Bequests", bins = range(0,stop = 30.0, length = 50), normalize =false, legend = false)
    end
    pretage = histogram(s.ret_age, title = "Ret Age", bins = range(0,stop = s.p.nT, length = s.p.nT), normalize = false, legend = false)
    mutil = map(x -> mean(s.util[:,x]), 1:s.p.nT)
    pmutil = plot(1:s.p.nT, mutil, title ="mean utility",ylims = (0,1),yticks = range(0,1,length = 6),legend = false)

    pempty = plot(legend=false,grid=false,foreground_color_subplot=:white)
    asize = 10
    # annotate!(pempty, [(0.3,1,Plots.text("gamma = $(s.p.gamma)",     :left, asize)),
    #                    (0.3,0.9,Plots.text("R = $(s.p.R)",     :left, asize)),
    #                    (0.3,0.8,Plots.text("beta = $(round(s.p.beta,digits=2))", :left, asize)),
    #                    (0.3,0.7,Plots.text("alpha = $(s.p.alpha)",   :left, asize)),
    #                    (0.3,0.6,Plots.text("sigma = $(s.p.sigma)",   :left, asize)),
    #                    (0.3,0.5,Plots.text("lambda = $(s.p.lambda)", :left, asize))#,
    #                    #(0.3,0.4,Plots.text("rho = $(s.p.ρ)",         :left, asize))
    #                    ])
    plot(py,pc,pw0,pretage,
        ppr,pw1,pbeq , pmutil,
          layout = (2,4), size = (800,400))
end

function plot_s(s::BSimulation)

    # inc, cons, w
    py = plot(s.inc',leg = false, title = "income")
    pc = plot(s.cons',leg = false, title = "consumption")
    i_retires = findall(s.file_age .> 0)
    scatter!(pc, s.file_age[i_retires], [s.cons[i,s.file_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))
    pw0 = plot(s.w0',leg = false, title = "w0")
    scatter!(pw0, s.file_age[i_retires], [s.w0[i,s.file_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))

    ppr = plot(s.prob_nofile',leg = false, title = "p(nofile)",ylims = (-0.05,1.05))
    scatter!(ppr, s.file_age[i_retires], [s.prob_nofile[i,s.file_age[i]] for i in i_retires], m = (:rect, 2, 0.6, :white))

    pw1 = plot(s.w1',leg = false, title = "w1")
    pempty = plot(legend=false,grid=false,foreground_color_subplot=:white)
    asize = 10
    annotate!(pempty, [(0.1,1,Plots.text("gamma = $(s.p.gamma)",     :left, asize)),
                       (0.1,0.9,Plots.text("R = $(round(s.p.R,digits=2))",     :left, asize)),
                       (0.1,0.8,Plots.text("beta = $(round(s.p.beta,digits=2))", :left, asize)),
                       (0.1,0.7,Plots.text("alpha = $(s.p.alpha)",   :left, asize)),
                       (0.1,0.6,Plots.text("sigma = $(s.p.sigma)",   :left, asize)),
                       (0.1,0.5,Plots.text("lambda = $(s.p.lambda)", :left, asize)),
                       (0.1,0.4,Plots.text("rho = $(s.p.ρ)",         :left, asize))])
    plot(py,pc,pw0,ppr,pw1,
          pempty ,
          layout = (2,3))
end

@recipe function f(m::FModel,p::Param;id=1)
    grid --> true
    xticks := true
    legend --> true
    cg = cgrad(:plasma)
    # c1 = colorant"red"
    # c2 = colorant"blue"
    # alow,ahi = extrema(m.avec)
    # aspect = (ahi-alow)/(ahi - 0.0)

    nT = size(m.v)[2]
    # cols = range(c1,stop=c2,length=nT)
    cols = cg[range(0.0,stop=1.0,length=nT)]
    xrange = zeros(2)
    xa = extrema(m.avec)
    xrange[1] = xa[1] .- diff(vcat(xa...))[1] .* 0.01
    xrange[2] = xa[2]

    layout := grid(1,2)
    for i in 1:nT
        vt = v_analytic(m,p,id,i)
        lab = ((i==1)|(i==nT)) ? "$i" : ""
        @series begin
            seriestype --> :path
            linewidth --> 1
            legend --> :bottomright
            seriescolor --> cols[i]
            label := lab
            subplot := 1  # value function
            yguide := "value"
            xguide := "Cash on Hand M"
            xlims := xrange
            ylims := (-15,15)
            getx(vt),gety(vt)
        end
        @series begin
            seriestype --> :path
            linewidth --> 1
            legend --> :bottomright
            seriescolor --> cols[i]
            # subplot := 2  #
            label := lab
            # xlims := (alow,ahi)
            # ylims := (0,ahi)
            ylims := extrema(m.avec)
            xlims := extrema(m.avec)
            yguide := "consumption"
            xguide := "Cash on Hand M"
            # aspect_ratio := aspect
            getx(m.c[id,i].env),gety(m.c[id,i].env)
        end
    end
end

@recipe function f(m::GModel,p::Param;id=1,iy=nothing,it=nothing)
    grid --> true
    xticks := true
    legend --> false
    cg = cgrad(:plasma)
    # c1 = colorant"red"
    # c2 = colorant"blue"

    xrange = zeros(2)
    xa = extrema(m.avec)
    xrange[1] = xa[1] .- diff(vcat(xa...))[1] .* 0.01
    xrange[2] = xa[2]

    nT = p.nT
    cols = cg[range(0.0,stop=1.0,length=nT)]

    layout := grid(1,2)

    if isnothing(iy) & !isnothing(it)
        legend := true
        title --> ["value period $it" "consumption period $it"]
        for jy in 1:p.ny
            vt = v_analytic(m,p,id,jy,it)
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend := :bottomright
                label := "y$jy"
                # seriescolor --> cols[i]
                subplot := 1  # value function
                yguide := "value"
                xguide := "Cash on Hand M"
                xlims := xrange
                # ylims --> (-15,15)
                getx(vt),gety(vt)
                # getx(m.v[id,iy,it].env),gety(m.v[id,iy,it].env)
            end
        end
        for jy in 1:p.ny
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :bottomright
                # seriescolor --> cols[i]
                subplot := 2  #
                label := "y$jy"
                xlims := xa
                ylims := xa
                yguide := "consumption"
                xguide := "Cash on Hand M"
                # aspect_ratio := aspect
                getx(m.c[id,jy,it].env),gety(m.c[id,jy,it].env)
            end
        end
    elseif isnothing(it) & !isnothing(iy)
        for i in 1:nT
            vt = v_analytic(m,p,id,iy,i)
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :bottomright
                seriescolor --> cols[i]
                label := ((i == 1) || (i == nT)) ? "$it" : ""
                # label := lab
                xlims := xrange
                # ylims --> (-15,15)
                subplot := 1  # value function
                yguide := "value"
                xguide := "Cash on Hand M"
                getx(vt),gety(vt)
            end
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :bottomright
                seriescolor --> cols[i]
                subplot := 2  #
                label := ((i == 1) || (i == nT)) ? "$it" : ""
                # label := lab
                xlims := xa
                ylims := xa
                yguide := "consumption"
                xguide := "Cash on Hand M"
                # aspect_ratio := aspect
                getx(m.c[id,iy,i].env),gety(m.c[id,iy,i].env)
            end
        end
    else
        println("you need to either give it or iy. not both. not none.")
    end

end

@recipe function f(m::BModel,p::Param;id=1,iy=nothing,it=nothing)
    grid --> true
    xticks := true
    # legend --> false
    cg = cgrad(:plasma)
    cols = cg[range(0.0,stop=1.0,length=p.nT)]

    # cg = cgrad(:inferno)
    # c1 = colorant"red"
    # c2 = colorant"blue"
    xrange = zeros(2)
    xa = extrema(m.avec)
    xrange[1] = xa[1] .- diff(vcat(xa...))[1] .* 0.01
    xrange[2] = xa[2]
    nT = p.nT
    # cols = range(c1,stop=c2,length=nT)

    layout := grid(1,3)

    rates = map(z->mean(z),m.rgrid[1,:,1:(nT-1)])

    if isnothing(iy) & !isnothing(it)
        legend := true
        title --> ["value period $it" "consumption period $it"]
        for jy in 1:p.ny
            vt = v_analytic(m,p,id,jy,it)
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :bottomright
                label := "y$jy"
                # seriescolor --> cols[i]
                subplot := 1  # value function
                yguide := "value"
                xguide := "Cash on Hand M"
                xlims := xrange
                # ylims --> (-15,15)
                getx(vt),gety(vt)
                # getx(m.v[id,iy,it].env),gety(m.v[id,iy,it].env)
            end
        end
        for jy in 1:p.ny
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :bottomright
                # seriescolor --> cols[i]
                subplot := 2  #
                label := "y$jy"
                xlims := xa
                ylims := xa
                yguide := "consumption"
                xguide := "Cash on Hand M"
                # aspect_ratio := aspect
                getx(m.c[id,jy,it].env),gety(m.c[id,jy,it].env)
            end
        end
        for jy in 1:p.ny
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> :topleft
                label := "y$jy"

                # seriescolor --> cols[i]
                subplot := 3  #

                yguide := "interest rate"
                xguide := "age"
                ylims := (0,maximum(rates))

                # aspect_ratio := aspect
                # println(rates[jy,:])
                collect(1:(p.nT-1)),rates[jy,:]
            end
        end
    elseif isnothing(it) & !isnothing(iy)
        for i in 1:nT
            vt = v_analytic(m,p,id,iy,i)
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> false
                seriescolor --> cols[i]
                # label := lab
                label := ((i == 1) || (i == nT)) ? "$it" : ""

                xlims := xrange
                # ylims --> (-15,15)
                subplot := 1  # value function
                yguide := "value"
                xguide := "Cash on Hand M"
                getx(vt),gety(vt)
            end
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> false

                seriescolor --> cols[i]
                subplot := 2  #
                # label := lab
                label := ((i == 1) || (i == nT)) ? "$it" : ""

                xlims := xa
                ylims := xa
                yguide := "consumption"
                xguide := "Cash on Hand M"
                # aspect_ratio := aspect
                getx(m.c[id,iy,i].env),gety(m.c[id,iy,i].env)
            end
        end

        for jy in 1:p.ny
            @series begin
                seriestype --> :path
                linewidth --> 1
                legend --> true
                label := "y$jy"
                # seriescolor --> cols[i]
                subplot := 3  #
                ylims := (0,maximum(rates))

                yguide := "interest rate"
                xguide := "age"
                # aspect_ratio := aspect
                # println(rates[jy,:])
                collect(1:(p.nT-1)),rates[jy,:]
            end
        end
    else
        println("you need to either give it or iy. not both. not none.")
    end
end

struct tester
    m :: Matrix
end

@recipe function f(tt::tester)
    n,m = size(tt.m)
    for i in 1:n
        @series begin
            seriestype --> :path
            series_annotations := ["$i" for i in 1:m]
            (1:m,tt.m[i,:])
        end
    end
end

@recipe function f(l::Vector{Point{T}};numerate=false,marker=false) where T
    # defaults
    grid --> true
    xticks := true
    @series begin
        seriestype --> :path
        linecolor --> :black
        linewidth --> 1
        if marker
            markershape --> :circle
            markerstrokecolor --> :black
            markercolor --> :white
            markersize --> 1
        end
        if numerate
            series_annotations := ["$i" for i in 1:length(l.v)]
        end
        (getx(l),gety(l))
    end
end

@recipe function f(l::MLine;numerate=false,marker=false)
    # defaults
    grid --> true
    xticks := true
    legend := false
    @series begin
        seriestype --> :path
        linecolor --> :black
        linewidth --> 1
        if marker
            markershape --> :circle
            markerstrokecolor --> :black
            markercolor --> :white
            markersize --> 1
        end
        if numerate
            series_annotations := ["$i" for i in 1:length(l.v)]
        end
        (getx(l),gety(l))
    end
end

@recipe function f(L::Vector{MLine{T}};numerate=false,mrk=true) where T
    for l in L
        @series begin
            # subplot := 1
            seriestype := :path
            linewidth := 1
            if mrk
                # println("illegal")
                markershape --> :circle
                markerstrokecolor --> :black
                markercolor --> :white
                markersize := 3
            end
            if numerate
                series_annotations := ["$i" for i in sortperm(getx(l))]
            end
            (getx(l),gety(l))
        end
    end
end
#
@recipe function f(x::Envelope;numerate=false,mrk=true)

    # defaults
    grid --> true
    xticks := true

    any_isec = length(x.isects) > 0

    # println("marker = $mrk, title = $(get!(plotattributes,:title,""))")
    # println("marker = $mrk, numerate = $numerate, title = $(get!(plotattributes,:title,""))")

    # actual envelope
    @series begin
        # subplot := 1
        seriestype := :path
        linecolor --> :red
        linewidth --> 2
        if mrk
            markershape := :circle
            markercolor := :white
            # markeralpha := 0.5
            markerstrokecolor := :black
            markersize := 3
        end
        if numerate
            series_annotations := ["$i" for i in 1:length(getx(x.env))]
        end
        (getx(x.env),gety(x.env))
    end
    if any_isec && mrk
        @series begin
            seriestype := :scatter
            markershape := :star6
            markerstrokecolor := :black
            markercolor := :yellow
            markersize := 5
            (getx(x.isects),gety(x.isects))
        end
    end
end
#
@recipe function f(x::Envelope{T},L::Vector{MLine{T}};numerate=false,removed=false,mrk=true) where T

    # defaults
    grid --> true
    xticks := true
    legend --> false

    any_isec = length(x.isects) > 0

    # println("marker = $mrk, title = $(get!(plotattributes,:title,""))")
    # println("marker = $mrk, numerate = $numerate, title = $(get!(plotattributes,:title,""))")

    # lien arrya
    for l in L
        @series begin
            # subplot := 1
            seriestype := :path
            linewidth := 1
            if mrk
                # println("illegal")
                markershape --> :circle
                markerstrokecolor --> :black
                markercolor --> :white
                markersize := 3
            end
            if numerate
                series_annotations := ["$i" for i in sortperm(getx(l))]
            end
            (getx(l),gety(l))
        end
    end

    # actual envelope
    @series begin
        # subplot := 1
        seriestype := :path
        linecolor --> :red
        linewidth --> 2
        if mrk
            markershape := :circle
            markercolor := :white
            # markeralpha := 0.5
            markerstrokecolor := :black
            markersize := 3
        end
        if numerate
            series_annotations := ["$i" for i in 1:length(getx(x.env))]
        end
        (getx(x.env),gety(x.env))
    end
    if any_isec && mrk
        @series begin
            seriestype := :scatter
            markershape := :star6
            markerstrokecolor := :black
            markercolor := :yellow
            markersize := 5
            (getx(x.isects),gety(x.isects))
        end
    end
end

# plot env and initial mline next to each other, showing new points and removed ones
@recipe function f(x::Envelope{T},L::MLine{T};numerate=false) where T

    # defaults
    grid --> true
    xticks := true
    legend --> false
    layout := (1,2)
    title --> ["Backwards-Bending" "Envelope"]

    any_isec = length(x.isects) > 0
    any_rmv = length(x.removed) > 0

    # println("marker = $mrk, title = $(get!(plotattributes,:title,""))")
    # println("marker = $mrk, numerate = $numerate, title = $(get!(plotattributes,:title,""))")

    # dirty line
    @series begin
        subplot := 1
        seriestype := :path
        linewidth := 1
        markershape --> :circle
        markerstrokecolor --> :black
        markercolor --> :white
        markersize := 3
        if numerate
            series_annotations := ["$i" for i in 1:length(getx(L))]
        end
        (getx(L),gety(L))
    end
    if any_rmv
        @series begin
            subplot := 1
            seriestype = :scatter
            markershape := :rect
            markersize := 3
            markerstrokecolor := :black
            markercolor := :white
            markeralpha := 0.5
            (getx(L)[x.removed],gety(L)[x.removed])
        end
    end

    # actual envelope
    @series begin
        subplot := 2
        seriestype := :path
        linecolor --> :red
        linewidth --> 1
        markershape --> :circle
        markerstrokecolor --> :black
        markercolor --> :white
        markersize := 3
        (getx(x.env),gety(x.env))
    end
    if any_isec
        @series begin
            subplot := 2
            seriestype := :scatter
            markershape := :star6
            markerstrokecolor := :black
            markercolor := :yellow
            markersize := 5
            (getx(x.isects),gety(x.isects))
        end
    end

end

#
# # @recipe function f(x::Envelope; removed=false,num=false,marker=false)
# @recipe function f(x::Envelope;numerate=false,removed=false,mrk=true)
#
#     # defaults
#     grid --> true
#     xticks := true
#
#     if !x.env_clean
#         legend --> :bottomright
#     else
#         legend --> false
#     end
#     any_isec = length(x.isects) > 0
#
#     # println("marker = $mrk, title = $(get!(plotattributes,:title,""))")
#     # println("marker = $mrk, numerate = $numerate, title = $(get!(plotattributes,:title,""))")
#
#     # if line array exists, plot
#     if length(x.L) > 0
#         for l in x.L
#             @series begin
#                 # subplot := 1
#                 seriestype := :path
#                 linewidth := 1
#                 if mrk
#                     # println("illegal")
#                     markershape --> :circle
#                     markerstrokecolor --> :black
#                     markercolor --> :white
#                     markersize := 3
#                 end
#                 if numerate
#                     series_annotations := ["$i" for i in sortperm(getx(l))]
#                 end
#                 (getx(l),gety(l))
#             end
#         end
#     end
#     # plot envelope, if exists
#     if x.env_clean
#         if removed
#             for l in 1:length(x.L)
#                 ir = x.removed[l]
#                 if length(ir) > 0
#                     @series begin
#                         seriestype = :scatter
#                         markershape := :rect
#                         markersize := 3
#                         markerstrokecolor := :black
#                         markercolor := :white
#                         markeralpha := 0.5
#                         (getx(x.L[l])[ir],gety(x.L[l])[ir])
#                     end
#                 end
#             end
#         end
#         @series begin
#             # subplot := 1
#             seriestype := :path
#             linecolor --> :red
#             linewidth --> 2
#             if mrk
#                 markershape := :circle
#                 markercolor := :white
#                 # markeralpha := 0.5
#                 markerstrokecolor := :black
#                 markersize := 3
#             end
#             if numerate
#                 series_annotations := ["$i" for i in 1:length(getx(x.env))]
#             end
#             (getx(x.env),gety(x.env))
#         end
#         if any_isec && mrk
#             @series begin
#                 seriestype := :scatter
#                 markershape := :star6
#                 markerstrokecolor := :black
#                 markercolor := :yellow
#                 markersize := 5
#                 (getx(x.isects),gety(x.isects))
#             end
#         end
#     end
# end

function f3a()

    fs = Function[x->ones(length(x));x->0.5x;x->x .- 2;x->2x .- 8]
    xs = [
         collect(range(-1 , stop = 0.9 ,length= 6))   ,
         collect(range(1  , stop = 7   ,length= 19))  ,
         collect(range(2  , stop = 7   ,length= 15))  ,
         collect(range(4  , stop = 8   ,length= 25))]
    ls = [MLine(i[1],i[2](i[1])) for i in zip(xs,fs)]
    return ls
end
function f3b()

    fs = Function[x->ones(length(x));x->0.5x;x->x .- 2;x->2x .- 8]
    xs = [
         collect(range(-1 , stop = 3.1 , length = 6))   ,
         collect(range(1  , stop = 7   , length = 19))  ,
         collect(range(2  , stop = 7   , length = 15))  ,
         collect(range(4  , stop = 8   , length = 25))]
    ls = [MLine(i[1],i[2](i[1])) for i in zip(xs,fs)]
    # create an envelope
    return ls
end
function f3c()

    fs = Function[x->ones(length(x));x -> 0.5x;x->x .- 2;x -> 2x .- 8.5]
    xs = [
         collect(range(-1 , stop = 3.1 , length= 6))   ,
         collect(range(1  , stop = 7   , length= 4))  ,
         collect(range(2  , stop = 7   , length= 15))  ,
         collect(range(4  , stop = 8   , length= 25))]
    ls = [MLine(i[1],i[2](i[1])) for i in zip(xs,fs)]
    # create an envelope
    return ls
end

function tplot1()
    n = 15
    x1 = collect(range(0 , stop = 10,length = n))
    x2 = collect(range(-1, stop = 9 ,length = n))
    L1 = MLine(x1,x1)
    L2 = MLine(x2,ones(n)*5)
    e = [L1,L2]
    x1,x2,e #,plot(e)
end

function tplot2()
    n = 15
    x1 = collect(range(0  , stop=10 ,length= n))
    x2 = collect(range(-1 , stop=9  ,length= n))
    L1 = MLine(x1,x1)
    L2 = MLine(x2,ones(n)*5)
    LL = [L1,L2]

    # a = splitLine(LL)
    e = upper_env(LL)
    p2 = plot(e,LL)
    savefig(p2,joinpath(@__DIR__,"..","images","tplot2.png"))
    p2
end
function tplot3a()

    en = f3a()

    p1 = plot(en,title = "non-overlapping grids")

    e = upper_env(en)
    p2 = plot(e,en,title = "Envelope")
    plot(p1,p2)

end

function tplot3b()

    en = f3b()

    p1 = plot(en,title = "overlapping grids")

    e = upper_env(en)
    p2 = plot(en)

    plot(p1,p2)

end
function tplot3c()

    en = f3c()

    p1 = plot(en)

    e = upper_env(en)
    removed!(en)
    p2 = plot(en,removed=true)

    plot(p1,p2)
    gui()
    return en

end

function tplot4()

    fs = Function[x->ones(length(x));x->0.5x;x->x .- 2;x->2x .- 8]
    xs = [
         collect(range(0.1 , stop =  1.5 ,length = 5))   ,
         collect(range(1   , stop =  7   ,length = 19))  ,
         collect(range(2   , stop =  7   ,length = 15))  ,
         collect(range(4   , stop =  8   ,length = 25))]
    ls = [MLine(i[1],i[2](i[1])) for i in zip(xs,fs)]
    e = DCEGM.upper_env(ls)

    p1 = plot(ls)
    p2 = plot(e,ls,title="envelope")
    p = plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","tplot4.png"))
    p
end

function tplot5()

    x1 = collect(-0.9:0.3:2.7)
    L1 = MLine(x1, x1)
    x2 = collect(0.0:0.1:1)
    L2 = MLine(x2, 2 .* x2)
    x3 = collect(1.0:0.45:2.9)
    L3 = MLine(x3, (0.1 .* x3) .+ 1.9 )
    ls = [L1,L2,L3]
    e = DCEGM.upper_env(ls)
    p1 = plot(ls)
    p2 = plot(e,ls,title="envelope")

    p = plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","tplot5.png"))
    p

end

function splitf()
    x = [1,2,3,1.5,2.1,2.9]
    y = [1,1.5,1.7,1.2,1.8,2.1]
    L = MLine(x,y)
    p1 = plot(L,title="original",numerate=true)
    e = splitLine(L)
    p2 = plot(e,title="split MLine",numerate=true)
    p = plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","split.png"))
    p
end

function splitf2()
    x = [1,2,3,2.9,2.5,1.9,1.8,1.5,2.1,2.9]
    y = [1,1.5,1.7,1.6,1.55,1.4,1.3,1.2,1.8,2.1]
    L = MLine(x,y)
    p1 = plot(L,title="original",numerate=true)
    e = splitLine(L)
    p2 = plot(e,title="split MLine",numerate=true)
    p = plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","split2.png"))
    p
end

function splitf3()
    x = [1,2,3,2.9,2.5,1.9,1.8,1.5,2.1,2.9]
    y = [1,1.5,1.7,1.6,1.55,1.4,1.3,1.2,1.8,2.1]
    L = MLine(x,y)
    e = DCEGM.secondary_envelope(L)
    p2 = plot(e,L,numerate=true)
    savefig(p2,joinpath(@__DIR__,"..","images","split3.png"))
    p2
end

function splitf3a()
    x = [1,2,3,2.9,2.5,1.9,1.8,1.5,2.1,2.9]
    y = [1,1.5,1.7,1.6,1.55,1.4,1.3,1.2,1.8,2.1]
    L = MLine(x,y)
    e = secondary_envelope(L)
    plot(e,L)

end

function test_upper_env_dec()
    n = 15
    x1 = collect(range(0,stop = 10,length = n))
    x2 = collect(range(-1,stop = 9,length = n))
    x = vcat(x1,x2)
    y = vcat(x1[end:-1:1],ones(n)*5)
    L = MLine(x,y)
    e = splitLine(L)
    p1 = plot(L,numerate=true,title = "initial line")
    p2 = plot(e,marker=true,title = "split line")
    e = upper_env(e)
    p3 = plot(e,marker=true,title = "Envelope w new pts")
    p4 = plot(p1,p2,p3,layout = (1,3),size=(1200,500))
    savefig(p4,joinpath(@__DIR__,"..","images","descending_upper.png"))
    L,e,x1,x2,p4
end

function split_test()
    f1(x) = ones(length(x))
    f2(x) = 0.5x
    f3(x) = x .-2
    f4(x) = 2x .-8
    x1 = collect(range(0.9,stop = 2.1,length = 14))
    x2 = collect(range(1,stop = 7,length = 19))
    x3 = collect(range(1,stop = 7,length = 15))
    x4 = collect(range(1,stop = 8,length = 25))
    X = [x1...,x2...,x3...,x4...]
    L = MLine([x1...,x2...,x3...,x4...],vcat(f1(x1),f2(x2),f3(x3),f4(x4)))
    en = splitLine(L)
    en,L
end

function demo(;n = 10,k = 10)
    a = Array{MLine{Float64}}(undef,10)
    for i in 1:2:n
        a[i] = MLine(sort([0.0; 10*rand(k); 10]), [0.0; rand(k+1).*collect(range(0.5,stop=10,length=k+1))])
        a[i+1] = MLine([0.0; 10],[10.0-i ; -(10-i)*1.5/i+10-i] )
    end
    p1 = plot(a,legend=false,title = "$n lines")
    e = upper_env(a)
    p2 = plot(e,a,title = "envelope")
    p =plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","demo.png"))
    p
end

function demo2(;n = 10,k = 10)
    a = Array{MLine{Float64}}(undef,10)
    slope = 1.5
    for i in 1:2:n
        a[i] = MLine(sort([0.0; 10*rand(k); 10]), [0.0; rand(k+1).*collect(range(0.5,stop=10,length=k+1))])
        a[i+1] = MLine([0.0; 10],[10.0-i ; -(10-i)*slope/i+10-i] )
    end
    # add a particularly tricky line
    tricky = [a[2].v[1]]
    push!(tricky,DCEGM.Point(a[2].v[2].x,a[2].v[2].y + slope))
    push!(a, MLine(tricky))

    p1 = plot(a,legend=false,title = "$n lines")
    e = upper_env(a)
    p2 = plot(e,a,title = "first intersection on grid")
    p = plot(p1,p2)
    savefig(p,joinpath(@__DIR__,"..","images","demo2.png"))
    p,e
end

function allplots()
    tplot1()
    tplot2()
    # savefig(joinpath(p,"tplot3a.png"))
    # tplot3b()
    # savefig(joinpath(p,"tplot3b.png"))
    # tplot3c()
    # savefig(joinpath(p,"tplot3c.png"))
    tplot4()
    tplot5()
    splitf()
    splitf2()
    splitf3()
    demo()
    demo2()
    test_upper_env_dec()
end


text:
module test_millboard_grid

using Millboard
using Test

Millboard.set_table_mode(:grid)

Millboard.display_style[:prepend_newline] = true

board = [11 12 13; 21 22 23]
@test sprint(show, table(board)) == """

+---+----+----+----+
|   | 1  | 2  | 3  |
+===+====+====+====+
| 1 | 11 | 12 | 13 |
+---+----+----+----+
| 2 | 21 | 22 | 23 |
+---+----+----+----+"""

Millboard.display_style[:prepend_newline] = false

end # module test_millboard_grid


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random

# test the basic floating point functions

@testset "flipsign" begin
    for elty in (Float32,Float64)
        x = convert(elty,-2.0)
        x = flipsign(x,-1.0)
        @test flipsign(x,big(-1.0)) == convert(elty,-2.0)
    end
end

@testset "maxintfloat" begin
    @test maxintfloat(Float16) === Float16(2048f0)
    for elty in (Float16,Float32,Float64)
        @test maxintfloat(rand(elty)) === maxintfloat(elty)
    end
    @test maxintfloat() === maxintfloat(Float64)
    @test maxintfloat(Float64, Int32) === 2147483647.0
    @test maxintfloat(Float32, Int32) === maxintfloat(Float32)
    @test maxintfloat(Float64, Int16) === 32767.0
    @test maxintfloat(Float64, Int64) === maxintfloat(Float64)
end

@testset "isinteger" begin
    for elty in (Float16, Float32, Float64)
        @test !isinteger(elty(1.2))
        @test isinteger(elty(12))
        @test isinteger(zero(elty))
        @test isinteger(-zero(elty))
        @test !isinteger(nextfloat(zero(elty)))
        @test !isinteger(prevfloat(zero(elty)))
        @test isinteger(maxintfloat(elty))
        @test isinteger(-maxintfloat(elty))
        @test !isinteger(elty(Inf))
        @test !isinteger(-elty(Inf))
        @test !isinteger(elty(NaN))
    end
end

@testset "ispow2" begin
    for T in (Float16,Float32,Float64,BigFloat)
        for x in (0.25, 1.0, 4.0, exp2(T(exponent(floatmax(T)))), exp2(T(exponent(floatmin(T)))))
            @test ispow2(T(x))
        end
        for x in (1.5, 0.0, 7.0, NaN, Inf)
            @test !ispow2(T(x))
        end
    end
end

@testset "round" begin
    for elty in (Float32, Float64)
        x = rand(elty)
        A = fill(x,(10,10))
        @test round.(A,RoundToZero) == fill(trunc(x),(10,10))
        @test round.(A,RoundUp) == fill(ceil(x),(10,10))
        @test round.(A,RoundDown) == fill(floor(x),(10,10))
        A = fill(x,(10,10,10))
        @test round.(A,RoundToZero) == fill(trunc(x),(10,10,10))
        @test round.(A,RoundUp) == fill(ceil(x),(10,10,10))
        @test round.(A,RoundDown) == fill(floor(x),(10,10,10))
        for elty2 in (Int32,Int64)
            A = fill(x,(10,))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,))
            A = fill(x,(10,10))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10))
            A = fill(x,(10,10,10))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10,10))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10,10))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10,10))
            @test round.(elty2,A) == fill(round(elty2,x),(10,10,10))
        end
    end
end

@testset "Types" begin
    for x in (Int16(0), 1, 2f0, pi, 3//4, big(5//6), 7.8, big(9), big(ℯ))
        @test float(typeof(x)) == typeof(float(x))
        @test float(typeof(complex(x, x))) == typeof(float(complex(x, x)))
    end
end

@testset "significant digits" begin
    # (would be nice to have a smart vectorized
    # version of signif)
    @test round(123.456, sigdigits=1) ≈ 100.
    @test round(123.456, sigdigits=3) ≈ 123.
    @test round(123.456, sigdigits=5) ≈ 123.46
    @test round(123.456, sigdigits=8, base = 2) ≈ 123.5
    @test round(123.456, sigdigits=2, base = 4) ≈ 128.0
    @test round(0.0, sigdigits=1) === 0.0
    @test round(-0.0, sigdigits=1) === -0.0
    @test round(1.2, sigdigits=2) === 1.2
    @test round(1.0, sigdigits=6) === 1.0
    @test round(0.6, sigdigits=1) === 0.6
    @test round(7.262839104539736, sigdigits=2) === 7.3
    @test isinf(round(Inf, sigdigits=3))
    @test isnan(round(NaN, sigdigits=3))
    @test round(1.12312, sigdigits=1000) === 1.12312
    @test round(Float32(7.262839104539736), sigdigits=3) === Float32(7.26)
    @test round(Float32(7.262839104539736), sigdigits=4) === Float32(7.263)
    @test round(Float32(1.2), sigdigits=3) === Float32(1.2)
    @test round(Float32(1.2), sigdigits=5) === Float32(1.2)
    @test round(Float16(0.6), sigdigits=2) === Float16(0.6)
    @test round(Float16(1.1), sigdigits=70) === Float16(1.1)
end

@testset "literal pow matches runtime pow matches optimized pow" begin
    two = 2
    @test 1.0000000105367122^2 == 1.0000000105367122^two
    @test 1.0041504f0^2 == 1.0041504f0^two

    function g2(start, two, N)
        x = start
        n = 0
        for _ in 1:N
           n += (x^2 !== x^two)
           x = nextfloat(x)
        end
        return n
    end
    @test g2(1.0, 2, 100_000_000) == 0
    @test g2(1.0f0, 2, 100_000_000) == 0
    g2′(start, N) = g2(start, 2, N)
    @test g2′(1.0, 100_000_000) == 0
    @test g2′(1.0f0, 100_000_000) == 0

    function g3(start, three, N)
        x = start
        n = 0
        for _ in 1:N
           n += (x^3 !== x^three)
           x = nextfloat(x)
        end
        return n
    end
    @test g3(1.0, 3, 100_000_000) == 0
    @test g3(1.0f0, 3, 100_000_000) == 0
    g3′(start, N) = g3(start, 3, N)
    @test g3′(1.0, 100_000_000) == 0
    @test g3′(1.0f0, 100_000_000) == 0

    function ginv(start, inv, N)
        x = start
        n = 0
        for _ in 1:N
           n += (x^-1 !== x^inv)
           x = nextfloat(x)
        end
        return n
    end
    @test ginv(1.0, -1, 100_000_000) == 0
    @test ginv(1.0f0, -1, 100_000_000) == 0
    ginv′(start, N) = ginv(start, -1, N)
    @test ginv′(1.0, 100_000_000) == 0
    @test ginv′(1.0f0, 100_000_000) == 0

    f(x, p) = x^p
    finv(x) = f(x, -1)
    f2(x) = f(x, 2)
    f3(x) = f(x, 3)
    x = 1.0000000105367122
    @test x^2 == f(x, 2) == f2(x) == x*x == Float64(big(x)*big(x))
    @test x^3 == f(x, 3) == f3(x) == x*x*x == Float64(big(x)*big(x)*big(x))
    x = 1.000000007393669
    @test x^-1 == f(x, -1) == finv(x) == 1/x == inv(x) == Float64(1/big(x)) == Float64(inv(big(x)))
end

@testset "curried approximation" begin

    @test ≈(1.0; atol=1).(1.0:3.0) == [true, true, false]

end


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

#Compiler: Julia 1.5

# Short description of this file: SIMPLE PRELIMINARY EXAMPLE OF NEURAL ODE using
# Julia's library DiffEqFlux.

Pkg.add("DiffEqSensitivity")
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, DiffEqSensitivity
using Plots

u0 = Float32[2.0; 0.0]
datasize = 30
tspan = (0.0f0, 1.5f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))

dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, tanh),
                  FastDense(50, 2))
neural_ode_f(u,p,t) = dudt2(u,p)
pinit = initial_params(dudt2)
prob = ODEProblem(neural_ode_f, u0, tspan, pinit)

function predict_neuralode(p)
  tmp_prob = remake(prob,p=p)
  Array(solve(tmp_prob,Tsit5(),saveat=tsteps,sensealg=BacksolveAdjoint(autojacvec=ReverseDiffVJP(true))))
end

function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end

callback = function (p, l, pred; doplot = true) #callback function to observe
  #training
  display(l)
  # plot current prediction against data
  plt =Plots.scatter(tsteps, ode_data[1,:], label = "data")
  scatter!(plt, tsteps, pred[1,:], label = "prediction")
  if doplot
    display(Plots.plot(plt))
  end
  return false
end

## Solve ODE and train

@time result_neuralode = DiffEqFlux.sciml_train(loss_neuralode, pinit,
                                          ADAM(0.05), cb = callback,
                                          maxiters = 500)

##Predict solution
result_neuralode

@show result_neuralode.minimizer
#=
2.687161 seconds (17.79 M allocations: 1002.418 MiB, 7.41% gc time)

* Status: success

* Candidate solution
  Final objective value:     2.761669e-02

* Found with
  Algorithm:     ADAM

* Convergence measures
  |x - x'|               = NaN ≰ 0.0e+00
  |x - x'|/|x'|          = NaN ≰ 0.0e+00
  |f(x) - f(x')|         = NaN ≰ 0.0e+00
  |f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
  |g(x)|                 = NaN ≰ 0.0e+00

* Work counters
  Seconds run:   3  (vs limit Inf)
  Iterations:    500
  f(x) calls:    500
  ∇f(x) calls:   500
=#


text:
function varargout = optimisation(varargin)
#  OPTIMISATION MATLAB code for optimisation.fig
#       OPTIMISATION, by itself, creates a new OPTIMISATION or raises the existing
#       singleton*.
#
#       H = OPTIMISATION returns the handle to a new OPTIMISATION or the handle to
#       the existing singleton*.
#
#       OPTIMISATION('CALLBACK',hObject,eventData,handles,...) calls the local
#       function named CALLBACK in OPTIMISATION.M with the given input arguments.
#
#       OPTIMISATION('Property','Value',...) creates a new OPTIMISATION or raises the
#       existing singleton*.  Starting from the left, property value pairs are
#       applied to the GUI before optimisation_OpeningFcn gets called.  An
#       unrecognized property name or invalid value makes property application
#       stop.  All inputs are passed to optimisation_OpeningFcn via varargin.
#
#       *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
#       instance to run (singleton)".
#
#  See also: GUIDE, GUIDATA, GUIHANDLES

#  Edit the above text to modify the response to help optimisation

#  Last Modified by GUIDE v2.5 13-Nov-2013 10:44:25

#  Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @optimisation_OpeningFcn, ...
                   'gui_OutputFcn',  @optimisation_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
#  End initialization code - DO NOT EDIT

#  --- Executes just before optimisation is made visible.
function optimisation_OpeningFcn(hObject, ~, handles, varargin)
#  This function has no output args, see OutputFcn.
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)
#  varargin   command line arguments to optimisation (see VARARGIN)

#  Choose default command line output for optimisation
handles.output = hObject;

#  Update handles structure
guidata(hObject, handles);

#  UIWAIT makes optimisation wait for user response (see UIRESUME)
#  uiwait(handles.figure1);

# File: optimisation.jl
# Author:
# Date: 29-Jan-2014

# import XL data
[names,~,returns] = getXLData('hfma_matlab_data.xlsx','10 Hedge Funds');

# parameters
N = size(returns,2); # size of portfolio (N)
R = fMean(returns,12)'; # annual returns (R)
W0 = repmat(1/N,N,1); # initial weights (W0)

covar = cov(returns); # covariance matrix
# annualise variance for monthly returns
for i=1:N
    for j=1:N
        if i == j
            covar(i,j)=covar(i,j)*12;
        end
    end
end
# VCV matrix
VCV = covar;

# annual portfolio statistics
pR = W0'*R; # return
pVar  = W0'*VCV*W0; # variance

# target return, linear contstraints and bounds
rstar = 9.00; # target return
Aeq = [R';ones(1,N)];
beq = [rstar;1];
lb = zeros(N,1);
ub = ones(N,1)*0.5;

# minimum variance optimisation
options = optimset('LargeScale', 'off', 'Algorithm','sqp');
[W,fval,ExitFlag] = fmincon(@(W) minVar(W,VCV),W0,[],[],Aeq,beq,lb,ub,[],options); # calls minVar.jl
    disp('Minimum Variance');
    if (ExitFlag ~= 1)
        disp('Optimization did not converge!');
    else
        disp('Optimization converged');
    end

# PLOT
# weights bar chart

# OUTPUT
# format data
R = num2cell(R);
    for i = 1:numel(R)
        R{i} = sprintf('# 12.3f',R{i});
    end
# set table data
set(handles.table_mean_returns,'ColumnName','Annual|Returns|(# )','RowName',names,'Data',R);
# format data
W0 = num2cell(W0);
    for i = 1:numel(W0)
        W0{i} = sprintf('# 12.2f',W0{i});
    end
set(handles.table_initial_weights,'ColumnName','Initial|Weights|(W0)','RowName',names,'Data',W0);

# format data
pStats = [pR;pVar];
pStats = num2cell(pStats);
    for i = 1:numel(pStats)
        pStats{i} = sprintf('# 10.3f',pStats{i});
    end
# set table data
set(handles.table_portfolio_statistics,'ColumnName','Annual|(# )','RowName',{'Port. Mean','Port. Variance'},'Data',pStats);
# format data
VCV = num2cell(VCV);
    for i = 1:numel(VCV)
        VCV{i} = sprintf('# 8.3f',VCV{i});
    end
# set size column widths
cwidth = num2cell(repmat(50,1,10));
# set table output
set(handles.table_VCV_matrix,'ColumnName',names,'ColumnWidth',cwidth,'RowName',names,'Data',VCV);

# format data
W = num2cell(W);
    for i = 1:numel(W)
        W{i} = sprintf('# 12.3f',W{i});
    end
set(handles.table_optimised_weights,'ColumnName','Optimised|Weights|(W)','RowName',names,'Data',W);

# format data
fval = num2cell(fval);
    for i = 1:numel(fval)
        fval{i} = sprintf('# 10.3f',fval{i});
    end
# set table data
set(handles.table_optimised_results,'ColumnName','Minimised|(# )','RowName','Port. Variance','Data',fval);

#  --- Outputs from this function are returned to the command line.
function varargout = optimisation_OutputFcn(~, ~, handles)
#  varargout  cell array for returning output args (see VARARGOUT);
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)

#  Get default command line output from handles structure
varargout{1} = handles.output;