text:
# -*- coding: utf-8 -*-

#
# Knowledge Tracing
# ベースとなるKnowledge Tracingのプログラムです．
# データ(カラムは['result', 'user_id', 'item_id', 'skill', 'date'])を読み込み，
# KTのパラメータを学習，スキル習得確率と設問に対する応答結果を予測します．
# 読み込む前に日付順にソートする必要があります．これも別スクリプトで可能です．
# AUCや正答率は別ファイルのevalution.pyで求めてください．
#

using DataFrames
include("./baum_welch_lib.jl")
using Hmm

## パラメータ学習の初期値
p_init = 0.1
ρ = [p_init 1-p_init]
p_learn = 0.1
p_forget = 0
A = [1-p_learn p_learn;  p_forget 1-p_forget]
p_slip = 0.1
p_guess = 0.1
B = [1-p_guess p_guess; p_slip 1-p_slip]

# データファイルの読み込み
function read_data(path)

  # データの読み込み
  df = readtable(path)

  # このカラムは必要ない
  delete!(df, :date)
  delete!(df, :item_id)

  # 正解:1、不正解:0　→ 正解:0、不正解:1　にする
  df2 = df[:result]
  df2[df2[:] .== 0] = 2
  df2[df2[:] .== 1] = 0
  df2[df2[:] .== 2] = 1
  df[:result] = df2

  # ユーザー一覧
  # 先行研究だとmost_commonでやってるから回答数の多いユーザーから順になってる
  # リストになってるけど俺のはなってない
  user_list = unique(df[:user_id])

  # スキル一覧
  skill_list = sort(unique(df[:skill]))

  return df, user_list, skill_list
end

# データの分割
function parse_data(df, user, skill)

  println("Parsing data ...")
  # 同じスキルを使う問題を抽出
  s_data = df[df[:skill] .== skill, :]

  # ユーザー数
  n = length(user)

  # 結果格納用
  test_data_df = df[1,:]
  deleterows!(test_data_df, 1)
  train_data_df = df[1,:]
  deleterows!(train_data_df, 1)
  train_data_p_u = []

  for i in 1:n
    # 同じユーザーが回答したものを抽出
    u_data = s_data[s_data[:user_id] .== user[i],:]
    # それをカウント
    u_data_count = length(u_data[:user_id])

    # テストデータを作る場合
    # 回答数が2回以下なら使えないデータ
    if (u_data_count <= 2)
      continue
    end

    # トレーニングデータ用の回答結果をリストに格納
    push!(train_data_p_u, u_data[1:end-1, :result])
    # トレーニングデータを抽出
    append!(train_data_df, u_data[1:end-1, :])
    # テストデータを抽出
    append!(test_data_df, u_data[end, :])
  end

  # テストデータが作れなかった場合は，仕方ないのでスキップする
  if length(test_data_df) == 0
    print("テストデータを作るためのデータが足りません...")
  end

  return train_data_p_u, test_data_df, train_data_df
end

# 成績予測
function predicter(train_data, test_data, hmm_ρ, hmm_A, hmm_B)

  # パラメータをセット
  pinit = hmm_ρ[1]
  plearn = hmm_A[1, 2]
  pguess = hmm_B[1, 2]
  pslip = hmm_B[2, 1]

  # 予測するユーザー数をカウント(テストデータの行数と一致)
  s_num = length(test_data[:result])

  # 正解:0、不正解:1　→ 正解:1、不正解:0　に戻す
  # トレーニングデータについて
  df2 = train_data[:result]
  df2[df2[:] .== 0] = 2
  df2[df2[:] .== 1] = 0
  df2[df2[:] .== 2] = 1
  train_data[:result] = df2
  # # テストデータについて
  df2 = test_data[:result]
  df2[df2[:] .== 0] = 2
  df2[df2[:] .== 1] = 0
  df2[df2[:] .== 2] = 1
  test_data[:result] = df2

  # 予測結果格納用
  # カラム名をつけておく
  p_te_all = DataFrame(result=[], user_id=[], skill=[], predict=[], probability=[])
  p_tr_all =  DataFrame(result=[], user_id=[], skill=[], predict=[], probability=[])

  # 各ユーザーについて求める
  for s in 1:s_num
    # ユーザー名
    user_name = test_data[s, :user_id]

    # そのユーザーの回答結果(トレーニングデータの中身)
    user_data = train_data[train_data[:user_id] .== user_name, :]
    cases = user_data[:result]
    # 回答回数
    n = length(cases)

    # スキル習得確率を計算
    pL = [pinit]
    for i in 1:n
      if cases[i] == 1.0
        num = (pL[i] * (1 - pslip)) / ((pL[i] * (1-pslip)) + (1-pL[i]) * pguess)
      else
        num = (pL[i] * pslip) / ((pL[i] * pslip) + (1-pL[i]) * (1-pguess))
      end
      numpL = num + (1-num) * plearn
      pL = vcat(pL, numpL)
    end

    # 設問に正解する確率を計算
    pC = []  # 確率
    CoI = [] # Correct or Incorrect
    for i in 1:n+1
      numC = (pL[i] * (1-pslip)) + ((1-pL[i]) * pguess)
      append!(pC, numC)
      append!(CoI, ((pC[i] > 0.5)?1:0))
    end

    # 各ユーザーに対する予測結果(トレーニングデータ)
    p_tr = hcat(user_data, CoI[1:end-1], pC[1:end-1])
    rename!(p_tr, :x1, :predict)
    rename!(p_tr, :x1_1, :probability)
    # 各ユーザーに対する予測結果(テストデータ)
    p_te = hcat(test_data[s, :], CoI[end], pC[end])
    rename!(p_te, :x1, :predict)
    rename!(p_te, :x1_1, :probability)

    # すべてのユーザーの予測結果に追加(トレーニングデータ)
    p_tr_all = vcat(p_tr_all, p_tr)
    # すべてのユーザーの予測結果に追加(テストデータ)
    p_te_all = vcat(p_te_all, p_te)
  end

  return p_tr_all, p_te_all
end

### MAIN ###
function main()

  # コマンドライン引数を読み取る
  # argc = length(ARGS)
  #
  # # 引数の個数チェック
  # if (argc-1 < 1)
  #   print("Usage: \$ julia", ARGS[1], "datafile.csv")
  #   exit()
  # end

  # データファイルの読み込み
  filepath = ARGS[1]
  df, user_list, skill_list = read_data(filepath)

  # HMMパラメータの読み込み
  HMM_model_flag = false
  # if argc-1 == 2
  #   HMM_model_flag = True
  #   HMM_model = pd.read_csv(argvs[2])
  # end

  # スキル数
  s_num = length(skill_list)

  # モデルのパラメータを格納する変数
  model = DataFrame(skill = [], init = [], learn = [], guess = [], slip = [])
  global predict_test = DataFrame([])
  global predict_train = DataFrame([])
  for i in 1:s_num
    println("+--------------------------------------------+**...")
    println("|  TRAINING SKILL ", skill_list[i], "")
    println("+--------------------------------------------+**...")

    # データの分割
    global train = 0
    train, test_df, train_df = parse_data(df, user_list, skill_list[i])
    if length(test_df) == 0
      println("skipping...")
      continue
    end

    # HMMのパラメータ学習
    if HMM_model_flag == true
      println("Loading HMM Parameters from Model File...")
      println("ここはまだやってない")
      # hmm = BaumWelch(A, B, pi)
      # params = HMM_model[HMM_model['skill'] == skill_list[i]]
      # # init
      # hmm.pi[0] = float(params['init'].values)
      # # learn
      # hmm.A[1, 0] = float(params['learn'].values)
      # # slip
      # hmm.B[0, 1] = float(params['slip'].values)
      # # guess
      # hmm.B[1, 0] = float(params['guess'].values)
    else
      # HMMのパラメータ学習
      println("Training HMM Parameters ...")
      hmm = Hmm.hmm_initialization(A, B, ρ)
      Hmm.train(hmm,train,1e-6,1000)
    end

    println("***  Learned Params  ***")
    println("init:", hmm.ρ[1], " learn:", hmm.A[1, 2], " guess:", hmm.B[1, 2], " slip:", hmm.B[2, 1])

    # 予測
    println("Predicting Future Student Performance ...")
    tr_res, te_res = predicter(train_df, test_df, hmm.ρ, hmm.A, hmm.B)

    # 結果を保持
    predict_test = vcat([predict_test, te_res])
    predict_train = vcat([predict_train, tr_res])

    # パラメータの記録
    tmp =  [skill_list[i], hmm.ρ[1], hmm.A[1, 2], hmm.B[1, 2], hmm.B[2, 1]]
    push!(model, tmp)
  end

  # # パラメータをファイルに保存
  # time = datetime.now().strftime("%Y%m%d-%H%M%S")
  # if  HMM_model_flag == False:
  #     model.to_csv("model_" + time + ".csv", index = False)
	writetable("model_params.csv", model)
  #
  # # 予測結果の出力
  # predict_test.to_csv("predict_test_" + time + ".csv", index = False)
  # predict_train.to_csv("predict_train_" + time + ".csv", index = False)
  println(predict_train)
  println(predict_test)
  println(model)
end

main()


text:
function system(components::Array{String,1}, method::String; kwargs...)
    # possible kwargs... are filepaths for
    # customdatabase_like, customdatabase_unlike, customdatabase_assoc
    set_components = [Set([components[i]]) for i in 1:length(components)]
    if method == "PCSAFT"
        raw_params = retrieveparams(components, method; kwargs...)
        model = PCSAFT(set_components, create_PCSAFTParams(raw_params))
    elseif method == "SAFTVRMie"
        raw_params = retrieveparams(components, method; kwargs...)
        model = SAFTVRMie(set_components, create_SAFTVRMieParams(raw_params))
    elseif method == "ogSAFT"
        raw_params = retrieveparams(components, method; kwargs...)
        model = ogSAFT(set_components, create_ogSAFTParams(raw_params))
    elseif method == "sPCSAFT"
        raw_params = retrieveparams(components, method; kwargs...)
        model = sPCSAFT(set_components, create_sPCSAFTParams(raw_params))
    end
    return model
end

ext: jl


text:
module GPbox

const Float = Float64

using Gadfly: layer, Geom
using Distributions
using Distances

export PlotGPsamples1D, estGP, GPkernel, genCovs

####Kernel functions####
struct GPkernel
    kernel_type::String
    param::Dict
end

####Covariance Matrix Generation####
function genCovs(x::Matrix{T},y::Matrix{T},gpk::GPkernel) where T<:Float
    if(gpk.kernel_type=="squared exp")
        nx = size(x,2);
        ny = size(y,2);
        lensq::T = (get(gpk.param,"len",1))^2;
        len_fact::T = 1/(2*lensq);
        covs = pairwise(SqEuclidean(),x,y);
        covs = exp.((-covs).*len_fact);
    end
    return covs
end

####Estimating posterior GP (mean and covarainces) given training point, test points and GP kernel####
function estGP(test_pts::AbstractArray{T},train_pts::AbstractArray{T},train_vals::AbstractArray{T},gpk::GPkernel) where T<:Float
    cv = genCovs(test_pts,test_pts,gpk);
    num_testpts = size(test_pts,2);
    num_fdim = size(train_vals,1);
    mn = zeros(num_fdim,num_testpts);
    traincv = genCovs(train_pts,train_pts,gpk);
    crosscv = genCovs(test_pts,train_pts,gpk);
    itraincv = inv(traincv);
    gp_mn = crosscv*itraincv*train_vals';
    temp = crosscv*itraincv*crosscv';
    gp_cv = cv - (temp + temp')/2;  #hack for positive definite CV
    return gp_mn, gp_cv
end

####Plotting Functions####

#Draw samples from 1D GP and output as array of Gadfly plotting layers
function PlotGPsamples1D(test_pts,nsam,gp_mn,gp_cov)
    dist = MultivariateNormal(gp_mn[:],gp_cov)
    f = rand(dist,nsam);
    layers = [layer(
                    x=test_pts,
                    y=f[:,i],
                    Geom.line) for i in 1:nsam];
    return layers
end

end


text:
@testset "Transitivity" begin
    completedg = CompleteDiGraph(4)
    circledg = PathDiGraph(4)
    add_edge!(circledg, 4, 1)
    for circle in testdigraphs(circledg)
        T = eltype(circle)
        complete = DiGraph{T}(completedg)
        newcircle = @inferred(transitiveclosure(circle))
        @test newcircle == complete
        @test ne(circle) == 4
        @test newcircle == @inferred(transitiveclosure!(circle))
        @test ne(circle) == 12
    end

    loopedcompletedg = copy(completedg)
    for i in vertices(loopedcompletedg)
        add_edge!(loopedcompletedg, i, i)
    end
    circledg = PathDiGraph(4)
    add_edge!(circledg, 4, 1)
    for circle in testdigraphs(circledg)
        T = eltype(circle)
        loopedcomplete = DiGraph{T}(loopedcompletedg)
        newcircle = @inferred(transitiveclosure(circle, true))
        @test newcircle == loopedcomplete
        @test ne(circle) == 4
        @test newcircle == @inferred(transitiveclosure!(circle, true))
        @test ne(circle) == 16
    end

    # transitivereduction
    let
        # transitive reduction of the nullgraph is again a nullgraph
        nullgraph = SimpleDiGraph(0)
        for g in testdigraphs(nullgraph)
            @test g == @inferred(transitivereduction(g))
        end

        # transitive reduction of a path is a path again
        pathgraph = PathDiGraph(10)
        for g in testdigraphs(pathgraph)
            @test g == @inferred(transitivereduction(g))
        end

        # transitive reduction of the transitive closure of a path is a path again
        for g in testdigraphs(pathgraph)
            gclosure = transitiveclosure(g)
            @test g == @inferred(transitivereduction(gclosure))
        end

        # Transitive reduction of a complete graph should be s simple cycle
        completegraph = CompleteDiGraph(9)
        for g in testdigraphs(completegraph)
            greduced = @inferred(transitivereduction(g))
            @test length(strongly_connected_components(greduced)) == 1 &&
                ne(greduced) == nv(g)
        end

        # transitive reduction a graph with no edges is the same graph again
        noedgegraph = SimpleDiGraph(7)
        for g in testdigraphs(noedgegraph)
            @test g == @inferred(transitivereduction(g))
        end

        # transitve reduction should maintain a selfloop only when selflooped==true
        selfloopgraph = SimpleDiGraph(1)
        add_edge!(selfloopgraph, 1, 1)
        for g in testdigraphs(selfloopgraph)
            @test g == @inferred(transitivereduction(g; selflooped=true));
            @test g != @inferred(transitivereduction(g; selflooped=false));
        end

        # transitive should not maintain selfloops for strongly connected components
        # of size > 1
        selfloopgraph2 = SimpleDiGraph(2)
        add_edge!(selfloopgraph2, 1, 1)
        add_edge!(selfloopgraph2, 1, 2)
        add_edge!(selfloopgraph2, 2, 1)
        add_edge!(selfloopgraph2, 2, 2)
        for g in testdigraphs(selfloopgraph2)
            @test g != @inferred(transitivereduction(g; selflooped=true));
        end

        # directed barbell graph should result in two cycles connected by a single edge
        barbellgraph = SimpleDiGraph(9)
        for i in 1:4, j in 1:4
            i == j && continue
            add_edge!(barbellgraph, i, j)
            add_edge!(barbellgraph, j, i)
        end
        for i in 5:9, j in 5:9
            i == j && continue
            add_edge!(barbellgraph, i, j)
            add_edge!(barbellgraph, j, i)
        end
        add_edge!(barbellgraph, 1, 5);
        for g in testdigraphs(barbellgraph)
            greduced = @inferred(transitivereduction(g))
            scc = strongly_connected_components(greduced)
            @test Set(scc) == Set([[1:4;], [5:9;]])
            @test ne(greduced) == 10
            @test length(weakly_connected_components(greduced)) == 1
        end
    end
end


text:
@testset "Ticker" begin
    res = ticker("BTC-USD")
    @test isa(res, DataFrame)

end


text:
using Clang
using Clang.Generators
using CSFML.LibCSFML.CSFML_jll

cd(@__DIR__)

include_dir = joinpath(CSFML_jll.artifact_dir, "include") |> normpath
sfml_dir = joinpath(include_dir, "SFML")
audio_h = joinpath(sfml_dir, "Audio.h")
graphics_h = joinpath(sfml_dir, "Graphics.h")
network_h = joinpath(sfml_dir, "Network.h")
system_h = joinpath(sfml_dir, "System.h")
window_h = joinpath(sfml_dir, "Window.h")
headers = [audio_h, graphics_h, network_h, system_h, window_h]

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

ctx = create_context(headers, args, options)

build!(ctx)


text:
"""
    to_vec(x)

Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
transformation.
"""
function to_vec(x::Real)
    function Real_from_vec(x_vec)
        return first(x_vec)
    end
    return [x], Real_from_vec
end

function to_vec(z::Complex)
    function Complex_from_vec(z_vec)
        return Complex(z_vec[1], z_vec[2])
    end
    return [real(z), imag(z)], Complex_from_vec
end

# Base case -- if x is already a Vector{<:Real} there's no conversion necessary.
to_vec(x::Vector{<:Real}) = (x, identity)

function to_vec(x::AbstractVector)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = cumsum(map(length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return oftype(x, x_Vec)
    end
    return vcat(x_vecs...), Vector_from_vec
end

function to_vec(x::AbstractArray)
    x_vec, from_vec = to_vec(vec(x))

    function Array_from_vec(x_vec)
        return oftype(x, reshape(from_vec(x_vec), size(x)))
    end

    return x_vec, Array_from_vec
end

# Some specific subtypes of AbstractArray.
function to_vec(x::Base.ReshapedArray{<:Any, 1})
    x_vec, from_vec = to_vec(parent(x))
    function ReshapedArray_from_vec(x_vec)
        p = from_vec(x_vec)
        return Base.ReshapedArray(p, x.dims, x.mi)
    end

    return x_vec, ReshapedArray_from_vec
end

# To return a SubArray we would endup needing to copy the `parent` of `x` in `from_vec`
# which doesn't seem particularly useful. So we just convert the view into a copy.
# we might be able to do something more performant but this seems good for now.
to_vec(x::Base.SubArray) = to_vec(copy(x))

function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    function AbstractTriangular_from_vec(x_vec)
        return T(reshape(back(x_vec), size(x)))
    end
    return x_vec, AbstractTriangular_from_vec
end

function to_vec(x::T) where {T<:LinearAlgebra.HermOrSym}
    x_vec, back = to_vec(Matrix(x))
    function HermOrSym_from_vec(x_vec)
        return T(back(x_vec), x.uplo)
    end
    return x_vec, HermOrSym_from_vec
end

function to_vec(X::Diagonal)
    x_vec, back = to_vec(Matrix(X))
    function Diagonal_from_vec(x_vec)
        return Diagonal(back(x_vec))
    end
    return x_vec, Diagonal_from_vec
end

function to_vec(X::Transpose)
    x_vec, back = to_vec(Matrix(X))
    function Transpose_from_vec(x_vec)
        return Transpose(permutedims(back(x_vec)))
    end
    return x_vec, Transpose_from_vec
end

function to_vec(X::Adjoint)
    x_vec, back = to_vec(Matrix(X))
    function Adjoint_from_vec(x_vec)
        return Adjoint(conj!(permutedims(back(x_vec))))
    end
    return x_vec, Adjoint_from_vec
end

function to_vec(X::PermutedDimsArray{T, D, P, Pinv} where {T, D}) where {P, Pinv}
    x_vec, back = to_vec(collect(X))
    function PermutedDimsArray_from_vec(x_vec)
        X_parent = collect(PermutedDimsArray(back(x_vec), Pinv))
        return PermutedDimsArray(X_parent, P)
    end
    return x_vec, PermutedDimsArray_from_vec
end

# Non-array data structures

function to_vec(x::Tuple)
    x_vecs, x_backs = zip(map(to_vec, x)...)
    sz = cumsum(collect(map(length, x_vecs)))
    function Tuple_from_vec(v)
        return ntuple(length(x)) do n
            return x_backs[n](v[sz[n] - length(x_vecs[n]) + 1:sz[n]])
        end
    end
    return reduce(vcat, x_vecs), Tuple_from_vec
end

function to_vec(x::NamedTuple)
    x_vec, back = to_vec(values(x))
    function NamedTuple_from_vec(v)
        v_vec_vec = back(v)
        return typeof(x)(v_vec_vec)
    end
    return x_vec, NamedTuple_from_vec
end

# Convert to a vector-of-vectors to make use of existing functionality.
function to_vec(d::Dict)
    d_vec, back = to_vec(collect(values(d)))
    function Dict_from_vec(v)
        v_vec_vec = back(v)
        return Dict(key => v_vec_vec[n] for (n, key) in enumerate(keys(d)))
    end
    return d_vec, Dict_from_vec
end

# ChainRulesCore Differentials
function FiniteDifferences.to_vec(x::Composite{P}) where{P}
    x_canon = canonicalize(x)  # to be safe, fill in every field and put in primal order.
    x_inner = ChainRulesCore.backing(x_canon)
    x_vec, back_inner = FiniteDifferences.to_vec(x_inner)
    function Composite_from_vec(y_vec)
        y_back = back_inner(y_vec)
        return Composite{P, typeof(y_back)}(y_back)
    end
    return x_vec, Composite_from_vec
end

function FiniteDifferences.to_vec(x::AbstractZero)
    function AbstractZero_from_vec(x_vec::Vector)
        return x
    end
    return Bool[], AbstractZero_from_vec
end


text:
@doc raw"""
    HopfieldRBM(g, θ, γ, w)
    HopfieldRBM(g, w)

Construct an RBM with spin visible units and Gaussian hidden units.
If not given, `θ = 0` and `γ = 1` by default.

```math
E(v, h) = -g'v - θ'h + \sum_\mu \frac{γ_\mu}{2} h_\mu^2 - v'wh
```
"""
function HopfieldRBM(g::AbstractArray, θ::AbstractArray, γ::AbstractArray, w::AbstractArray)
    @assert size(θ) == size(γ)
    @assert size(w) == (size(g)..., size(θ)...)
    return RBM(Spin(g), Gaussian(θ, γ), w)
end

function HopfieldRBM(g::AbstractArray{T,N}, w::AbstractArray) where {T,N}
    @assert size(g) == size(w)[1:N]
    m = size(w)[(N + 1):end]
    θ = zeros(eltype(g), m)
    γ = ones(eltype(g), m)
    return HopfieldRBM(g, θ, γ, w)
end

function HopfieldRBM(::Type{T}, n::Union{Int,TupleN{Int}}, m::Union{Int,TupleN{Int}}) where {T}
    g = zeros(T, n)
    θ = zeros(T, m)
    γ = ones(T, m)
    w = zeros(T, n..., m...)
    return HopfieldRBM(g, θ, γ, w)
end

HopfieldRBM(n::Union{Int,TupleN{Int}}, m::Union{Int,TupleN{Int}}) = HopfieldRBM(Float64, n, m)


text:
module BrmlFunctions
# General helper functions:

#VERSION < v"0.4-" && using Docile
#@docstrings

import Base: vec, findmax

# Type Unions:
export IntOrIntArray
IntOrIntArray=Union{Integer,Array{Integer},Array{}} # include also the ranges and vectors

export vec
@doc """
 Creates a vector with 1 element of value v
""" ->
function vec(v::Integer)
    return vec([v])
end

export maxarray
@doc """
 Finds the largest non-zero value of matrix **`A`**
""" ->
function maxarray(A)
    p,q,r=findnz(A)
    return maximum(r)
end

export SparseIntMatrix
@doc """
 Initialises a new matrix of size **M x N**

* `M`:     number of rows
* `N`:     number of columns
""" ->
function SparseIntMatrix(M,N)
    #return Int(spzeros(M,N))
    return round.(Int64,spzeros(M,N))
end

export fillmatrix!
@doc """
 Initialises a subarray of **`A`** with value **`val`**

* `A`:         array to fill
* `indsi`:     vector/range of row indeces
* `indsj`:     vector/range of column indeces
* `val`:       value to fill with the elements
""" ->
function  fillmatrix!(A,indsi,indsj,val)
    for i=indsi
        for j=indsj
            A[i,j]=val
        end
    end
    return A
end

export states
@doc """
 Enumerates all states given the number of states of each variable from a given set similar to MATLAB _ind2subv_

* `ns`:    vector with the total number of states for each variable
""" ->
function states(ns)
    # enumerate all states eg states([2 2 3]) (like ind2subv.m)
    # first index changes the most (like matlab)
    #ns=fliplr(ns)
    ns=flipdim(ns,2)
    n=length(ns)
    nstates=prod(vec(ns))
    c=cumprod(vec(ns[end:-1:1]))
    c=vcat(1,c[1:end-1])
    s=Array(Int64,nstates,n)
    for i=1:nstates
        tmp=i-1
        for j=n:-1:1
            s[i,n+1-j]=floor(tmp/c[j])
            tmp=tmp-s[i,n+1-j]*c[j]
        end
    end
    #return fliplr(s.+1)
    return flipdim(s.+1,2)
end

export standardise
@doc """
 Transforms an one-dimensional array **`A`** of numerical values into a column vector,
  otherwise it leaves the input unchanged
""" ->
function standardise{T<:Number}(A::Array{T,})
    # If the array is a vector, make this a column vector, otherwise leave unchanged
    if isavector(A)
        A=vec(A)
    end
    return A
end

export mysize
@doc """
 Returns the size for each dimension of a vector or matrix

* `a`:     vector or matrix to analyse
""" ->
function mysize(a)
    a=standardise(a)
    s=size(a)
    if isavector(a)
        sz=s[1]
    else
        sz=s[1:end]
    end
    return sz
end

export memberinds
@doc """
 Returns a vector mapping the indices of **`x`** and **`y`** for which the corresponding elements have the same value

* `x`:     item to compare
* `y`:     item to compare
""" ->
function memberinds(x,y)
    ind=zeros(Int64,length(x),1)
    for i=1:length(x)
        for j=1:length(y)
            if y[j]==x[i]
                ind[i]=j;
            end
        end
    end
    return ind
end

export myvcat
@doc """
 Wrapper for standard **`vcat`** used for
 backward compatibility from porting the MATLAB version of BRML code
""" ->
function myvcat(x,y)
    if length(x)==0
        return y
    end

    if length(y)==0
        return x
    end

    return vcat(x,y)
end

export isavector
@doc """
 Verifies if matrix **`A`** has one of the dimensions of size **1**
""" ->
function isavector(A)

    if length(size(A))>2
        return false
    elseif  (size(A,1)==1 && size(A,2)>1) || (size(A,2)==1 && size(A,1)>1)
        return true
    else
        return false
    end

end

export numstates
@doc """
 Returns the number of elements of matrix **`A`**
""" ->
function numstates(A)
    if isavector(A)
        n=prod(size(A));
    else
        n=size(A);
    end
    return n
end

export replace!
@doc """
 Finds the elements having value **`f`** in matrix **`A`** and replaces their value with **`r`**
""" ->
function replace!(A,f,r)
    A[find(A.==f)]=r
end

export findone
@doc """
 Returns the indices of the first element in matrix **`A`** that has the value **`c`**
""" ->
function findone(A,c)
    A, B=findn(A.==c)
    A=A[1]; B=B[1]
    return A,B
end

# Basic graph routines

export neighboursize
@doc """
 Returns the number of neighbours in an graph

     nsize = neighboursize(A,<node>)

* `A`: adjacency matrix of the graph
* `node`: (_optional_) the node for which to find the number of neighbours

If **`node`** is missing, return the neighbour sizes (including self) of each node.
If **`A`** is directed, returns the number of parents of the specified node.
""" ->
function neighboursize(A,nodes=[])
    if isempty(nodes)
        nsize=sum(A,1);
    else
        nsize=sum(A[:,vec(nodes)],1);
    end
    return nsize
end

export istree
@doc """
 Checks if the input graph is singly-connected (a polytree)

    tree, elimseq, schedule = istree(A, <root>=[]; <ReturnElimSeq>=false)

##### Input :
* `A`: graph's adjacency matrix (zeros on diagonal)
* `root`: (_optional_) root node of the graph

##### Outputs:
* `tree`: _true_ if graph is singly connected, otherwise _false_
* `elimseq`: a variable elimination sequence in which simplical nodes of the tree are listed,
 as each simplical node is removed from the tree.
* `schedule`: the sequence of messages from node to node corresponding to elimseq

If **`A`** is a directed graph, the elimination schedule begins with the nodes with no children.
 If root is specified, the last node eliminated is root.

If the graph is connected and the number of edges is less than the number of nodes,
 it must be a tree. However, to deal with the general case in which it is unknown if the graph
 is connected w check using elimination.

A tree/singly-connected graph must admit a recursive simplical node elimination. That is at
 any stage in the elimination there must be a node with either 0 or 1 neighbours in the remaining graph.
""" ->
function istree(A,root=[];ReturnElimSeq=false)
    C = size(A,1); # number of nodes in the graph
    schedule=zeros(Int,C,2)
    tree=true; # assume A is singly connected
    AA=copy(A); # adjacency of the eliminated graph
    elimseq=[]; # set of variables eliminated (in sequence)

    for node=1:C
        # now find the number of neighbours:
        nn=(C+1)*ones(1,C)  # ensures that we don't pick eliminated nodes
        s=1:C; r=zeros(1,C)
        r[elimseq]=1; # need to check this
        s=s[find(r.==0)];
        nn[s]=neighboursize(AA',s)
        if !isempty(root)
            nn[root]=C+1 # ensures we don't pick root
        end

        val, elim=findmin(nn) # find node with least number of neighbours
        neigh = find(AA[:,vec(elim)]) # find the non zero elements

        if length(neigh)>1; # if least has more than 1 neighbour, cannot be a tree
            tree=false
            break;
        end

        AA[vec(elim),:]=0;  AA[:,vec(elim)]=0; # eliminate node from graph
        elimseq=[elimseq... elim]; # add eliminated node to elimination set

        if isempty(neigh);  schedule[node,:]=[elim elim];
        else
            schedule[node,:]=[elim neigh];
        end

    end
    if !tree;
        if ReturnElimSeq
            return false, [],[]
        else
            return false
        end
    end

    c=[]
    for i=1:size(schedule,1)
        if schedule[i,1]!=schedule[i,2] # remove self elimination
            c=[c... i]
        end
    end
     if ReturnElimSeq
         return tree, elimseq, schedule[vec(c),:]
     else
         return tree
     end
end

export findmax
@doc """
 Maximises a multi-dimensional array over a set of dimensions

    [maxval maxstate] = findmax(A, variables)

* `A`:         array to fill
* `variables`: <???>

Finds the values and states that maximize the multi-dimensional
 array **`A`** over the dimensions in <???>maxover
""" ->
function findmax(A,variables;Ind2Sub=false)

    maxval,maxind=findmax(A,variables)
    if !Ind2Sub
        return maxval, maxind
    end

    s=ind2sub(size(A),maxind[:]) # compatibilty with matlab BRML
    maxstate=zeros(Int64,length(s[1]),length(s))
    for i=1:length(s)
        maxstate[:,i]=s[i]
    end
    return maxval, maxstate

end

export randgen
@doc """
 Returns a random value from distribution **`p`**
""" ->
function randgen(p)
    p=p./sum(p)
    f=find(rand().<cumsum(vec(p)))
    return(Int(minimum(f)))
end

export myipermutedims
@doc """
 Wrapper for **`ipermutedims`** used for backward compatibility
 from porting the MATLAB version of BRML code
""" ->
function myipermutedims{T<:AbstractArray}(A::T,d)
    if length(d)==1 && d[1]==1
        return A
    else
        #ipermutedims(A,d)
        permutedims(A,invperm(d[:]))
    end
end

export DictToArray
@doc """
 Copies the dictionary **D** values into a new vector
 with the size equal to the number of keys in the dictionary
""" ->
function DictToArray(D)
    if isa(D,Dict)
        L=length(collect(keys(D)))
        pot=Array(Any,L)
        ky=collect(keys(D))
        for k=1:length(keys(D))
            pot[k]=D[ky[k]]
        end
        return pot
    else
        return D
    end
end

#export iskey
#function iskey(k::Any,D::Dict)
#    allkeys=collect(keys(D))
#    return any(k.==allkeys)
#end

export condp
@doc """
 Creates a conditional distribution from an array

    pnew = condp(pin, varargin)

##### Input:
* `pin`:      a positive matrix
* `varargin`: optional input specifying which indices form the distribution variables

##### Output:
* `pnew`:    a new matrix with **`sum(pnew, 1) = ones(1, size(p, 2))`**

##### Example:
    r = rand(4, 2, 3);
    p = condp(r, DistributionIndices=[3, 1]);

**`p`** is now an array of the same size as **`r`**, but with **`sum(sum(p,3),1) = 1`** for each of the dimensions of the 2nd index.
""" ->
function condp(p; DistributionIndices::IntOrIntArray=[]) ## FIXME! This doesn't work when p is more than an 2D array
    p=p+realmin(); # in case all unnormalised probabilities are zero

    if isempty(DistributionIndices)
        p=p./repmat(sum(p,1),size(p,1),1)
        return p
    else
    if DistributionIndices==0
        #p=p./sum(p[:])
        p=p./sum(p)
        return p
    end
        allvars=1:length(size(p))
        sizevars=size(p)
        condvars=setdiff(allvars,DistributionIndices)
        newp=deepcopy(permutedims(p,vcat(DistributionIndices,condvars)))
        newp=reshape(newp,prod(sizevars[DistributionIndices]),prod(sizevars[condvars]))
        newp=newp./repmat(sum(newp,1),size(newp,1),1)
        pnew=reshape(newp,sizevars[vcat(DistributionIndices,condvars)])
        #pnew=ipermutedims(pnew,vcat(DistributionIndices,condvars))
        pnew=ipermutedims(pnew,invperm(vcat(DistributionIndices,condvars)[:]))
        return pnew
    end
end

export normp
@doc """
 Creates a normalised array from an array

    pnew = normp(pin)

##### Input:
* `pin`:      an array

##### Output:
* `pnew`:    a new array with **`sum(pnew) =1`**

##### Example:
    r = rand(4, 2, 3);
    p = normp(r);

**`p`** is now an array of the same size as **`r`**, but with **`sum(p) = 1`**
""" ->
function normp(p) ## FIXME! This doesn't work when p is more than an 2D array
    p=p+realmin(); # in case all unnormalised probabilities are zero
    p=p./sum(p)
    return p
end

#function condp(p, DistributionIndices=[]) ## FIXME! This doesn't work when p is more than an 2D array

#    if isempty(DistributionIndices)
#        p=p./repmat(sum(p,1),size(p,1),1)
#        return p
#    else
#    if DistributionIndices==0
#        p=p./sum(p[:])
#        return p
#    end
#        p=p+realmin(); # in case all unnormalised probabilities are zero
#        allvars=1:length(size(p))
#        sizevars=size(p)
#        condvars=setdiff(allvars,DistributionIndices)
#        newp=deepcopy(permutedims(p,vcat(DistributionIndices,condvars)))
#        newp=reshape(newp,prod(sizevars[DistributionIndices]),prod(sizevars[condvars]))
#        newp=newp./repmat(sum(newp,1),size(newp,1),1)
#        pnew=reshape(newp,sizevars[vcat(DistributionIndices,condvars)])
#        pnew=ipermutedims(pnew,vcat(DistributionIndices,condvars))
#        return pnew
#    end
#end

export condexp
@doc """
 Computes **`p`** proportional to **`exp(log p)`**
""" ->
function condexp(logp)
    pmax=maximum(logp,1)
    P =size(logp,1)
    return condp(exp.(logp-repmat(pmax,P,1)))
end

export validgridposition
@doc """
 Returns `true` if point (x, y) is on a defined grid (1:Gx, 1:Gy)

    v = validgridposition(x, y, Gx, Gy)
""" ->
function validgridposition(x, y, Gx, Gy)
    if x > Gx || x < 1
        return false
    end
    if y > Gy || y < 1
        return false
    end
    return true
end

end #module


text:
function makegrid(layers=3,startlocs=[(0,0,2)],groundlevel=false)
	grid=Set{Tuple}()
	push!(grid,startlocs...)
	connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0), (0,0,1),(1,0,1),(0,1,1),(0,0,-1),(1,0,-1),(1,-1,-1)]
	if groundlevel
		connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0)]
	end
	for layer in 1:layers
		tgrid=Array{Tuple,1}()
		for loc in grid
			if loc[3]==2
				for c in connections
					x,y,z=loc
					x+=c[1];y+=c[2];z+=c[3]
					push!(tgrid,(x,y,z))
				end
			end
		end
		for t in tgrid
			push!(grid,t)
		end
	end
	return grid
end
function adjacent(hex,spacing=1,layer=false)
	connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0), (0,0,1),(1,0,1),(0,1,1),(0,0,-1),(1,0,-1),(1,-1,-1)]
	if layer
		connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0)]
	end
	if hex[3]==1
		if layer
			connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0)]
		else
			connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0), (0,0,1),(-1,0,1),(-1,1,1)]
		end
	elseif hex[3]==3
		if layer
			connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0)]
		else
			connections=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(1,-1,0),(-1,1,0),(0,0,-1),(-1,0,-1),(0,-1,-1)]
		end
	end
	adj=Array{Tuple,1}()
	for c in connections
		x,y,z=hex
		x+=spacing*c[1];y+=spacing*c[2];z+=spacing*c[3]
		push!(adj,(x,y,z))
	end
	return adj
end

function freelocs(game,layer=2)
	free=0
	tot=0
	for (loc,col) in game.map
		if loc[3]==layer
			tot+=1
			if col==0
				free+=1
			end
		end
	end
	return (free,tot)
end

function numcolors(rgb)
	nc=0
	for c in rgb
		if c>0
			nc+=1
		end
	end
	return nc
end

function irlocs(unit::Unit)
	return makegrid(unit.ir,[unit.loc],unit.groundlevel)
end

function hex_to_pixel(q,r,size)
    x = size * sqrt(3) * (q + r/2)
    y = size * 3/2 * r
    return x, y
end
function pixel_to_hex(x,y,size)
    q = (x * sqrt(3)/3 - y / 3) / size
    r = y * 2/3 / size
    return (q, r)
end

function triangle(ctx,x,y,size,up=-1)
	polygon(ctx, [Point(x,y),Point(x+size,y),Point(x+size/2,y+up*size)])
	fill(ctx)
end
function hexlines(ctx,x,y,size)
	size*=2
	move_to(ctx,x-size/4,y-size*sin(pi/3)/2)
	rel_line_to(ctx,size/2,size*sin(pi/3))
	move_to(ctx,x-size/2,y)
	rel_line_to(ctx,size,0)
	move_to(ctx,x+size/4,y-size*sin(pi/3)/2)
	rel_line_to(ctx,-size/2,size*sin(pi/3))
	stroke(ctx)
end
function samegroup(group1::Group,group2::Group)
	for unit in group1.body
		if !in(unit,group2.body)
			return false
		end
	end
	return true
end
function hasgroup(grouparray,group::Group)
	for g in grouparray
		if samegroup(group,g)
			return true
		end
	end
	false
end


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ dbb88400-065d-11eb-3fe8-9b5b7f873d3e
begin
	using MortalityTables
	using PlutoUI
	using Plots
end

# ╔═╡ de790392-065d-11eb-047d-a10af7f2fdc7
md"
# Mortality Table Comparison Tool

*For full features, use MortalityTables.jl version `2.1.2` or higher.*
"

# ╔═╡ 5322a010-065f-11eb-0864-3d69c92cff6c
md"**Issue Age:** $(@bind issue_age Slider(18:100,show_value=true,default=50))"

# ╔═╡ 968d79c0-0695-11eb-2586-0565817b37a8
md"
## The rest of the code

The non-plotting stuff.
"

# ╔═╡ 4694f4aa-47dc-4148-99e0-62ad3cde88a3
function life_expectancies(table,select,age_range)
	map(age_range) do age
		if checkbounds(Bool, select ? table.select : table.ultimate, age)
			age + life_expectancy(select ? table.select[age] : table.ultimate,age)
		else
			missing
		end
	end
end

# ╔═╡ ad6c9830-068e-11eb-2934-6f460c421895
function rate_ratio(a, a_select,b, b_select, issue_age,duration)
		att_age = issue_age + duration - 1
		a = a_select ? a.select[issue_age][att_age] : a.ultimate[att_age]
		b = b_select ? b.select[issue_age][att_age] : b.ultimate[att_age]
		return b / a
end

# ╔═╡ df760e20-0692-11eb-16bb-e5b41f8a5255
available_names = sort!([x.name for x in MortalityTables.table_source_map])

# ╔═╡ 26f84770-065e-11eb-03ca-435dfb6ed5da
md"
**Table A:** $(@bind name_a Select(
available_names,
default=\"1975-80 Modified Basic Table with Milliman Extension - Male, ALB\"
))

*or* Mort.SOA.org table ID (leave blank to use dropdown): $(@bind lookup_a TextField())

Use select rates (if available): $(@bind sel_a CheckBox(default=true))
"

# ╔═╡ b6fbd7f0-065f-11eb-1af7-5141be0a585f
if lookup_a == ""

	table_a = MortalityTables.table(name_a)
else
	table_a = MortalityTables.table(parse(Int,lookup_a))
end

# ╔═╡ 7ae73f7e-065e-11eb-1472-b7bc0e5cff29
md"
**Table B:** $(@bind name_b Select(
available_names,
default=\"2015 VBT Male Non-Smoker RR100 ALB\"
))

*or* Mort.SOA.org table ID (leave blank to use dropdown): $(@bind lookup_b TextField())

Use select rates (if available): $(@bind sel_b CheckBox(default=true))
"

# ╔═╡ bb15f410-065f-11eb-2342-9319b0ff6237
if lookup_b == ""
	table_b = MortalityTables.table(name_b)
else
	table_b = MortalityTables.table(parse(Int,lookup_b))
end

# ╔═╡ d60bc480-065e-11eb-2072-b3ac00e39598
let
	if sel_a
		rates_a = table_a.select[issue_age][issue_age:end]
	else
		rates_a = table_a.ultimate[issue_age:end]
	end
	if sel_b
		rates_b = table_b.select[issue_age][issue_age:end]
	else
		rates_b = table_b.ultimate[issue_age:end]
	end
	# left plot
	p1 = plot(
		rates_a,
		label="$name_a",
		title="unscaled",
		xaxis="duration",
		legend=:topleft
	)
	plot!(p1, rates_b,label="$name_b")

	# right plot
	p2 = plot(
		rates_a,
		label="$name_a",
		yaxis=:log,
		title="log scaled",
		xaxis="duration",
		legend=false
	)
	plot!(p2, rates_b,label="$name_b")

	plot(p1,p2,size=(980,400),titlefontsize=9)

end

# ╔═╡ 54c92a0c-c8fe-40b5-bcf7-a791c8999b23
let
	p = plot(ylim=(65,100),legend=:bottomright,title="Life Expectancy",palette=:PuOr_4,size=(600,400),xlabel="issue age")
	age_range = 20:100
	le_a = life_expectancies(table_a,sel_a,age_range)
	le_b = life_expectancies(table_b,sel_b,age_range)
	scatter!(age_range,le_a, m = (:white,stroke(:dodgerblue3)),label="$name_a")
	scatter!(age_range,le_b, m = (:white,stroke(:salmon1)),label="$name_b")
	vline!([issue_age],alpha=0.5,color=:grey,linewidth=3,label="selected issue age",titlefontsize=9)
end

# ╔═╡ 378a33f2-0695-11eb-3dee-af6f6298edad
# the matrix of ratios for the contour plots
begin
	issue_ages = 25:80
	ω = min(omega(table_b.ultimate),omega(table_a.ultimate))
	durations = 1:(ω-80)
	# durations = 1:40

	ratio = [rate_ratio(table_a,sel_a,table_b,sel_b,ia,dur) for ia in issue_ages, dur in durations]
end;

# ╔═╡ e38d2460-068f-11eb-3d19-e1902fc7b28c
begin

	c1 = contour(durations,
	        issue_ages,
	        ratio,
			title="Ratio",
	        xlabel="duration",ylabel="select issue age",
			titlefontsize=8,
	        fill=true,
		size=(675,400),
		contour_labels=true
	        )
	hline!(c1,[issue_age],alpha=0.5,color=:grey,linewidth=3,legend=false)

	c2 = contour(durations,
	        issue_ages,
	        log.(ratio),
	        xlabel="duration",ylabel="issue age",
	        title="Ratio (LOG SCALED)",
			titlefontsize=8,
	        fill=true,
		size=(675,400),
		color=:turbo,
		clim=(-3,3),
		contour_labels=true,
		levels=50
	        )
	hline!(c2,[issue_age],alpha=0.5,color=:grey,linewidth=3,legend=false)

	 l = @layout [
    		title{0.05h}
			grid(1,2)
            ]
	titleplot=plot(title="$name_b over \n $name_a\n(Grey bar is issue age selected above)", leg = false, ticks = nothing, border = :none,titlefontsize=9)

	plot(titleplot,c1,c2,size=(980,400),layout= l)

end

# ╔═╡ 03dabe48-7a4a-11eb-252f-5bd03b2e0201
# make the notebook just a little wider
html"""<style>
pluto-notebook {
    width: 1000px;
}
"""

# ╔═╡ Cell order:
# ╟─de790392-065d-11eb-047d-a10af7f2fdc7
# ╟─26f84770-065e-11eb-03ca-435dfb6ed5da
# ╟─7ae73f7e-065e-11eb-1472-b7bc0e5cff29
# ╟─5322a010-065f-11eb-0864-3d69c92cff6c
# ╟─d60bc480-065e-11eb-2072-b3ac00e39598
# ╟─e38d2460-068f-11eb-3d19-e1902fc7b28c
# ╟─54c92a0c-c8fe-40b5-bcf7-a791c8999b23
# ╟─968d79c0-0695-11eb-2586-0565817b37a8
# ╠═dbb88400-065d-11eb-3fe8-9b5b7f873d3e
# ╠═b6fbd7f0-065f-11eb-1af7-5141be0a585f
# ╠═4694f4aa-47dc-4148-99e0-62ad3cde88a3
# ╠═bb15f410-065f-11eb-2342-9319b0ff6237
# ╠═ad6c9830-068e-11eb-2934-6f460c421895
# ╠═378a33f2-0695-11eb-3dee-af6f6298edad
# ╟─df760e20-0692-11eb-16bb-e5b41f8a5255
# ╟─03dabe48-7a4a-11eb-252f-5bd03b2e0201

max_stars_repo_name: EarthGoddessDude/Learn

max_issues_repo_name: EarthGoddessDude/Learn

max_forks_repo_name: EarthGoddessDude/Learn


text:
function create_sysimage(sysout::AbstractString, add_dlmreader::Bool = true; add_package::Bool = true)
    juliaCMD = Base.julia_cmd()
    tmpIMD = mktemp()
    close(tmpIMD[2])
    if add_package
        run(`$juliaCMD -e "using Pkg; Pkg.add(\"PackageCompiler\")"`)
        if add_dlmreader
            run(`$juliaCMD -e "using Pkg; Pkg.add(\"DLMReader\")"`)
        end
    end
    warmupcode = "using InMemoryDatasets; ds = IMD.warmup();println(ds)"
    if add_dlmreader
        tmpDLM = mktemp()
        close(tmpDLM[2])
        warmupcode *= "; using DLMReader; DLMReader.warmup(); repeat!(ds, 10000); ds.x1 = rand(10000); ds.x2 = rand(Int, 10000); filewriter(\"$(escape_string(tmpDLM[1]))\", ds); filereader(\"$(escape_string(tmpDLM[1]))\");"
    end

    run(`$juliaCMD --trace-compile=$(tmpIMD[1]) -e "$warmupcode"`)
    tmpIMD_out = mktemp()
    close(tmpIMD_out[2])
    f = open(tmpIMD[1])
    fout = open(tmpIMD_out[1], "w")
    # TODO should make this part less hard coding
    write(fout, "using REPL, Pkg, InMemoryDatasets, InMemoryDatasets.PooledArrays, InMemoryDatasets.DataAPI, InMemoryDatasets.PrettyTables, InMemoryDatasets.Tables, PackageCompiler, TOML, Logging, SuiteSparse
    ")
    if add_dlmreader
        write(fout, "; using DLMReader, DLMReader.InlineStrings, DLMReader.InlineStrings.Parsers
        ")
    end
    for line in eachline(f)
        # not sure why this happens
        write(fout, replace(line, "000000"=>""))
    end
    close(fout)
    close(f)
    run(`$juliaCMD -e "using PackageCompiler;
    create_sysimage(:InMemoryDatasets, sysimage_path=\"$(escape_string(sysout))\", precompile_execution_file=\"$(escape_string(tmpIMD_out[1]))\")"`)
    println("Now exit julia and re-run it again using $(sysout) as `--sysimage`, e.g. in unix type OS you can use the following command:
    julia --sysimage $sysout")
end


text:
using Compat

runtest(filename) = (println(filename); include(filename))

runtest("api.jl")


text:
function bob(stimulus::AbstractString)
    stimulus = strip(stimulus)

    if issilence(stimulus)
        return "Fine. Be that way!"
    elseif isquestion(stimulus) && isshouting(stimulus)
        return "Calm down, I know what I'm doing!"
    elseif isshouting(stimulus)
        return "Whoa, chill out!"
    elseif isquestion(stimulus)
        return "Sure."
    else
        return "Whatever."
    end
end

issilence(stimulus::AbstractString) = isempty(stimulus)
isquestion(stimulus::AbstractString) = endswith(stimulus, '?')

function isshouting(stimulus::AbstractString)
    all(isuppercase, stimulus) && return true
    !any(isletter, stimulus) && return false

    for c in stimulus
        # ignore all non-letter chars
        if isletter(c) && !isuppercase(c)
            return false
        end
    end

    return true
end


text:
using Test
using PyCall
using BenchmarkTools
using Fermi
using Fermi.Wavefunction
using Fermi.CoupledCluster: RCCSD,RCCD,DFRCCD

psi4.core.be_quiet() #turn off output
psi4.set_num_threads(6)
using LinearAlgebra
BLAS.set_num_threads(6)
# > setup
tol = 1E-14
psi4.set_options(Dict("D_CONVERGENCE" => 14,
                      "E_CONVERGENCE" => 12,
					  "scf_type" => "pk"))
mol2 = psi4.geometry("""
                     O
                     H 1 1.1
                     H 1 1.1 2 104.0
					 symmetry c1
					 """)
e2,wfn2 = psi4.energy("hf/cc-pvtz",mol=mol2,return_wfn=true)
psi4.set_options(Dict("D_CONVERGENCE" => 14,
                      "E_CONVERGENCE" => 12,
                      "scf_type"      => "df"))
e3,wfn3 = psi4.energy("hf/cc-pvtz",mol=mol2,return_wfn=true)
JuWfn3 = Wfn(wfn3)
JuWfn2 = Wfn(wfn2)
printdo=false
#println(@btime RCCD.do_rccd(JuWfn2; doprint=printdo))
println(@btime RCCSD.do_rccsd(JuWfn2; doprint=printdo))
#println(@btime DFRCCD.do_df_rccd(JuWfn3; doprint=printdo))
#println(psi4.energy("ccsd/sto-3g",mol=mol2) - e2)


text:
function trs(P, q::AbstractVector{T}, r::T; kwargs...) where {T}
	output = trs_boundary(P, q, r; kwargs...)
	return check_interior!(output..., P, q)
end

function trs(P, q::AbstractVector{T}, r::T, C::AbstractMatrix{T}; kwargs...) where {T}
	output = trs_boundary(P, q, r, C; kwargs...)
	return check_interior!(output..., P, q)
end

function check_interior!(X::AbstractMatrix, info::TRSinfo, P, q::AbstractVector; direct=false)
	if info.λ[1] < -1e-9 # Global solution is in the interior
		x1 = X[:, 1]
		if !direct
			cg!(x1, P, -q, reltol=(eps(real(eltype(q)))/2)^(2/3))
		else
			x1 = -(P\q)
		end
		X[:, 1] = x1
		info.λ[1] = 0
	end
	if size(X, 2) > 1 && info.λ[2] < -1e-9
		# No local-no-global minimiser can exist in the interior; discard it.
		X = X[:, 1:1]; info.λ = info.λ[1:1]
	end
	return X, info
end


text:
# Define a Half Arrow Matrix Type
mutable struct HalfArrow{T} <: AbstractMatrix{T}
    D::Vector{T} # diagonal
    z::Vector{T} # last column
end

# Define its size

size(A::HalfArrow{T}, dim::Integer) where T = dim==1 ? max(length(A.D),length(A.z)) :
length(A.D)+1
size(A::HalfArrow{T}) where T = size(A,1), size(A,2)

# Index into a SymArrow
function getindex(A::HalfArrow{T},i::Integer,j::Integer) where T
    n=length(A.D)+1
    if i==j<n return A.D[i]
    elseif j==n; return A.z[i]
    else return zero(typeof(A.D[1]))
    end
end

# Dense version of HalfArrow
Matrix(A::HalfArrow{T}) where T =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]


text:
module Jobs
# This module handles all the Job related functionality
using Parameters, StructTypes, CodeTracking, JLD2, Dates, JSON3
using ..Calculations
using ..Structures
using ..Utils
using ..DFControl

include("environment.jl")
include("job.jl")
include("versioning.jl")
include("registry.jl")
export Job, Environment
end


text:
using ACE, JuLIP, ACEatoms, ACEbase, Test, LinearAlgebra
using ACE: evaluate, evaluate_d, SymmetricBasis, NaiveTotalDegree, PIBasis
using ACEbase.Testing: fdtest

# I2LModel = ACEatoms.Models.I2LModel

##

# construct the 1p-basis
@info("Construcing L1")
D = NaiveTotalDegree()
maxdeg = 5
ord = 3
species = [:Ti, :Al]
B1p = ACEatoms.ZμRnYlm_1pbasis(; species = species, maxdeg=maxdeg, D = D,
                                 rin = 1.2, rcut = 5.0)
ACE.init1pspec!(B1p, maxdeg = maxdeg, Deg = ACE.NaiveTotalDegree())
φ = ACE.Invariant()
pibasis = PIBasis(B1p, ord, maxdeg; property = φ)
basis = SymmetricBasis(pibasis, φ)
cTi = randn(length(basis))
cAl = randn(length(basis))
L1 = Dict(:Ti => ACE.LinearACEModel(basis, cTi; evaluator = :standard),
          :Al => ACE.LinearACEModel(basis, cAl; evaluator = :standard) )

##
@info("Construcing L2")
Pk_in = ACE.transformed_jacobi(maxdeg, ACE.Transforms.IdTransform(),
                               1.0, 0.0; pin = 0, pcut = 0)
Pk = ACE.Scal1pBasis(Pk_in, :u, :k)
Bnlmk = B1p * Pk
ACE.init1pspec!(Bnlmk, maxdeg = maxdeg, Deg = ACE.NaiveTotalDegree())
length(Bnlmk)
φ = ACE.Invariant()
pibasis2 = PIBasis(Bnlmk, ord, maxdeg; property = φ)
basis2 = SymmetricBasis(pibasis2, φ)
cTi = randn(length(basis2))
cAl = randn(length(basis2))
L2 = Dict(:Ti => ACE.LinearACEModel(basis2, cTi; evaluator = :standard),
          :Al => ACE.LinearACEModel(basis2, cAl; evaluator = :standard) )

##

@info("Construct I2LModel")
model = ACEatoms.Models.I2LModel(L1, L2)

@info("Basic parameter tests")
println(@test( ACE.nparams(model) == 2 * length(basis) + 2 * length(basis2) ))
c2 = randn(ACE.nparams(model))
ACE.set_params!(model, c2)
println(@test( ACE.params(model) == c2 ))

##

@info("construct random atomic structure")

zTi = AtomicNumber(:Ti)
zAl = AtomicNumber(:Al)
at = bulk(:Ti, cubic=true) * 3
at.Z[2:3:end] .= zAl
at = rattle!(at, 0.1)

##

Us = rand(length(at))
nlist = neighbourlist(at, 5.0)
env1, _ = ACEatoms.environment(model.L1, at, nlist, 1)
env2 = ACEatoms.Models.AtEnvU1(at, nlist, 1, Us)
ACE.evaluate(model.L1.models[at.Z[1]], env1)
ACE.evaluate(model.L2.models[at.Z[1]], env2)

##

energy(model, at)

## forces

## derivatives w.r.t. parameters


text:
#
###
# Lagrange polynomial function spaces
###

""" Lagrange polynomial spectral space """
struct LagrangePolynomialSpace{T,D,Tpts,
                               Tdom<:AbstractDomain{T,D},
                               Tquad,Tgrid,
                               Tmass<:AbstractField{T,D},
                               Tderiv, #Tl2g,
                              } <: AbstractSpectralSpace{T,D}
    domain::Tdom
    npoints::Tpts
    quadratures::Tquad
    grid::Tgrid
    mass_mat::Tmass
    deriv_mats::Tderiv
#   local2global::Tl2g # TODO
end

function LagrangePolynomialSpace(domain::AbstractDomain{<:Number,1}, n;
                                 quadrature = gausslobatto,
                                 T = Float64,
                                )

    """ reset deformation to map from [-1,1]^D """
    ref_domain = unit_square(2)
    domain = ref_domain # map_from_ref(domain, ref_domain) # TODO

    z, w = quadrature(n)

    D = lagrange_deriv_mat(z)

    z = T.(z)
    w = T.(w)

    npoints = (n,)
    quadratures = ((z, w),)
    grid = (Field(z),)
    mass_mat = Field(w)
    deriv_mats = (D,)

    space = LagrangePolynomialSpace(
                                    domain, npoints, quadratures, grid,
                                    mass_mat, deriv_mats,
                                   )
    domain isa DeformedDomain ? deform(space, mapping) : space
end

GaussLobattoLegendre1D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausslobatto, kwargs...)
GaussLegendre1D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausslegendre, kwargs...)
GaussChebychev1D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausschebyshev, kwargs...)

function LagrangePolynomialSpace(domain::AbstractDomain{<:Number,2}, nr, ns;
                                 quadrature = gausslobatto,
                                 T = Float64,
                                )
    msg = "spectral polynomials work with logically rectangular domains"
    @assert domain isa BoxDomain msg

    """ reset deformation to map from [-1,1]^D """
    ref_domain = unit_square(2)
    domain = ref_domain # map_from_ref(domain, ref_domain) # TODO

    zr, wr = quadrature(nr)
    zs, ws = quadrature(ns)

    zr, wr = T.(zr), T.(wr)
    zs, ws = T.(zs), T.(ws)

    r, s = ndgrid(zr,zs)

    Dr = lagrange_deriv_mat(zr)
    Ds = lagrange_deriv_mat(zs)

    npoints = (nr, ns,)
    quadratures = ((zr, wr), (zs, ws),)
    grid = (r, s)
    mass_mat = Field(wr * ws')
    deriv_mats = (Dr, Ds,)

    space = LagrangePolynomialSpace(
                                    domain, npoints, quadratures, grid,
                                    mass_mat, deriv_mats,
                                   )
    domain isa DeformedDomain ? deform(space, mapping) : space
end

GaussLobattoLegendre2D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausslobatto, kwargs...)
GaussLegendre2D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausslegendre, kwargs...)
GaussChebychev2D(args...; kwargs...) = LagrangePolynomialSpace(args...; quadrature=gausschebyshev, kwargs...)

get_grid(space::LagrangePolynomialSpace) = space.grid
get_domain(space::LagrangePolynomialSpace) = space.domain
get_numpoints(space::LagrangePolynomialSpace) = space.points

function massOp(space::LagrangePolynomialSpace)
    @unpack mass_mat = space

    DiagonalOp(mass_mat)
end

function gradOp(space::LagrangePolynomialSpace{<:Number,1})
    (Dr,) = space.deriv_mats

    Dx = MatrixOp(Dr)

    [Dx]
end

function gradOp(space::LagrangePolynomialSpace{<:Number,2})
    (nr, ns) = space.npoints
    (Dr, Ds) = space.deriv_mats

    Ir = sparse(I, nr, nr)
    Is = sparse(I, ns, ns)

    Dx = TensorProductOp2D(Dr, Is)
    Dy = TensorProductOp2D(Ir, Ds)

    [Dx
     Dy]
end

function gradOp(space::LagrangePolynomialSpace{<:Number,3})
    (Dr, Ds, Dt) = space.deriv_mats
    (nr, ns, nt) = space.npoints

    Ir = sparse(I, nr, nr)
    Is = sparse(I, ns, ns)
    It = sparse(I, nt, nt)

    Dx = TensorProductOp3D(Dr, Is, It)
    Dy = TensorProductOp3D(Ir, Ds, It)
    Dz = TensorProductOp3D(Ir, Is, Dt)

    [Dx
     Dy
     Dz]
end

function interpOp(space1::LagrangePolynomialSpace{<:Number,1},
                  space2::LagrangePolynomialSpace{<:Number,1},
                 )
    r1, _ = space1.quadratures[1]
    r2, _ = space2.quadratures[1]

    J = lagrange_interp_mat(r2, r1) # from 1 to 2

    MatrixOp(J)
end

function interpOp(space1::LagrangePolynomialSpace{<:Number,2},
                  space2::LagrangePolynomialSpace{<:Number,2},
                 )
    r1, _ = space1.quadratures[1]
    r2, _ = space2.quadratures[1]

    s1, _ = space1.quadratures[2]
    s2, _ = space2.quadratures[2]

    Jr = lagrange_interp_mat(r2, r1) # from 1 to 2
    Js = lagrange_interp_mat(s2, s1)

    TensorProductOp3D(Jr, Js)
end

function interpOp(space1::LagrangePolynomialSpace{<:Number,3},
                  space2::LagrangePolynomialSpace{<:Number,3},
                 )
    r1, _ = space1.quadratures[1]
    r2, _ = space2.quadratures[1]

    s1, _ = space1.quadratures[2]
    s2, _ = space2.quadratures[2]

    t1, _ = space1.quadratures[3]
    t2, _ = space2.quadratures[3]

    Jr = lagrange_interp_mat(r2, r1) # from 1 to 2
    Js = lagrange_interp_mat(s2, s1)
    Jt = lagrange_interp_mat(t2, t1)

    TensorProductOp3D(Jr, Js, Jt)
end
#


text:
@reexport module CirculationGeneration

using LinearAlgebra, SparseArrays
using Unitful

#=============================================
Tools for making simple circulations
=============================================#

"""
    T_advection(ϕ, orig, dest, v3D, nb)

Returns the sparse matrix transport operator, `T`.

`T` is such that it gives the flux divergence due to
the volumetric flow rate `ϕ` (in m³ s⁻¹) from box `orig` to box `dest`.
"""
function T_advection(ϕ::Tϕ, orig, dest, v3D::Array{Tv,3}, nb) where {Tϕ, Tv}
    unit(Tϕ) / unit(Tv) ≠ unit(u"1/s") && error("Unit problem")
    return sparse([orig, dest], [orig, orig], ustrip.([ϕ / v3D[orig], -ϕ / v3D[dest]]), nb, nb)
end
"""
    T_advection(ϕ, sequence, v3D, nb)

Returns the sparse matrix transport operator, `T`.

`T` is such that it gives the flux divergence due to
the volumetric flow rate `ϕ` (in m³ s⁻¹) going through all the boxes in `sequence`.
"""
function T_advection(ϕ, sequence, v3D, nb)
    origs = sequence[1:end-1]
    dests = sequence[2:end]
    sum(T_advection(ϕ, orig, dest, v3D, nb) for (orig,dest) in zip(origs, dests))
end
"""
    T_diffusion(ν, i, j, v3D, nb)

Returns the sparse matrix transport operator, `T`.

`T` is such that it gives the flux divergence due to
the volumetric mixing rate `ν` (in m³ s⁻¹) between boxes `i` and `j`.
"""
function T_diffusion(ν, i, j, v3D, nb)
    T_advection(ν, i, j, v3D, nb) + T_advection(ν, j, i, v3D, nb)
end

end


text:
module StringSearch

include("sais.jl")
include("lcp.jl")
include("align.jl")

end


text:
"""
    vm_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart)

construct varimax NGWP and GP_star.tag in place.

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `GP_star::GraphPart`: GraphPart object of the dual graph

# Output Argument
- `wavelet_packet::Array{Float64,3}`: the varimax NGWP. The first index is for
    selecting wavelets at a fixed level; the second index is for selecting the
    level `j`; the third index is for selecting elements in the wavelet vector.

"""
function vm_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart)
    rs = GP_star.rs
    inds = GP_star.inds
    (N, jmax) = Base.size(inds)

    GP_star.tag = zeros(Int, N, jmax)
    GP_star.tag[:, 1] = Vector{Int}(0:(N - 1))

    wavelet_packet = zeros(N, jmax, N)
    wavelet_packet[:, 1, :] = Matrix{Float64}(I, N, N)

    for j = 2:jmax
        regioncount = count(!iszero, rs[:, j]) - 1
        for r = 1:regioncount
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            GP_star.tag[indr, j] = Vector{Int}(0:(length(indr) - 1))
            wavelet_packet[indr, j, :] = varimax(𝚽[:, inds[indr, j]])'
        end
    end

    return wavelet_packet
end


text:
# Load X and y variable
using JLD
data = load("quantum.jld")
(X,y) = (data["X"],data["y"])

# Add bias variable, initialize w, set regularization and optimization parameters
(n,d) = size(X)
lambda = 1

# Initialize
maxPasses = 10
progTol = 1e-4
verbose = true
w = zeros(d,1)
lambda_i = lambda/n # Regularization for individual example in expectation

# Start running stochastic gradient
w_old = copy(w);
for k in 1:maxPasses*n

    # Choose example to update 'i'
    i = rand(1:n)

    # Compute gradient for example 'i'
    r_i = -y[i]/(1+exp(y[i]*dot(w,X[i,:])))
    g_i = r_i*X[i,:] + (lambda_i)*w

    # Choose the step-size
    alpha = 1/(lambda_i*k)

    # Take thes stochastic gradient step
    w -= alpha*g_i

    # Check for lack of progress after each "pass"
    if mod(k,n) == 0
        yXw = y.*(X*w)
        f = sum(log.(1 + exp.(-y.*(X*w)))) + (lambda/2)norm(w)^2
        delta = norm(w-w_old,Inf);
        if verbose
            @printf("Passes = %d, function = %.4e, change = %.4f\n",k/n,f,delta);
        end
        if delta < progTol
            @printf("Parameters changed by less than progTol on pass\n");
            break;
        end
        w_old = copy(w);
    end
end


text:
"""
Ipopt wrapper
"""

"""
    minimize_solver_specific(solver::IPOPT, cache, x0, lx, ux, lg, ug, rows, cols, outputfile)

Minimize function interfacing to IPOPT

# Arguments
    - `cache::GradOrJacCache`: cache created using _create_cache()

Nanami's Function = {red wine + accounting + driving + Bay Area + singing + eating out + coffee + pesto genovese + strawberry + mango + cinnamon roll + french toast + Korean food^(Faily and Friednds)}
The most difficult function on the planet
"""
function minimize_ipopt(options::Dict, cache, x0, lx, ux, lg, ug, rows, cols, outputfile::Bool)

    # initialize
    nx = length(x0)
    ng = length(lg)

    # initialize data for caching results since ipopt separates functions out
    xlast = 2*x0
    f = 0.0
    g = zeros(ng)
    df = zeros(nx)
    dg = zeros(length(rows))

    function ipcon(x, g)  # ipopt calls this function first
        f = evaluate!(g, df, dg, x, cache)
        xlast = x
        return nothing
    end

    function ipobj(x)
        if !isequal(x, xlast)
            f = func!(g, x)
        end
        return f
    end

    function ipgrad_obj(x, grad)
        if !isequal(x, xlast)
            f = evaluate!(g, df, dg, x, cache)
            xlast = x
        end
        grad[:] = df
        return nothing
    end

    function ipjac_con(x, mode, r, c, values)
        if mode == :Structure
            r[:] = rows
            c[:] = cols
        else
            if !isequal(x, xlast)
                f = evaluate!(g, df, dg, x, cache)
                xlast = x
            end
            values[:] = dg
        end
        return nothing
    end

    nzJ = length(rows)
    nH = 1  # irrelevant for quasi-newton
    prob = Ipopt.createProblem(nx, lx, ux, ng, lg, ug, nzJ, nH,
        ipobj, ipcon, ipgrad_obj, ipjac_con)

    # set Ipopt options
    Ipopt.addOption(prob, "hessian_approximation", "limited-memory")
    for (key, value) in options
        Ipopt.addOption(prob, key, value)
    end

    # open output file
    if outputfile==true
        # initialize options
        filename = "ipopt.out"
        print_level = 5
        if haskey(options, "output_file")
            filename = options["output_file"]
        end
        if haskey(options, "print_level")
            print_level = options["print_level"]
        end
        # create Ipopt output file
        Ipopt.openOutputFile(prob, filename, print_level)
    end

    # solve problem
    prob.x = x0
    status = solveProblem(prob)

    # return xstar, fstar, info
    xstar = prob.x
    fstar = prob.obj_val
    info  = Ipopt.ApplicationReturnStatus[status]
    return xstar, fstar, info
end

"""
    process_ipopt_out(filename::String)

Process ipopt output file
"""
function process_ipopt_out(filename::String)
    # get lines
    lines = []
    open(filename) do f
        # read till end of file
        while ! eof(f)
            # read a new / next line for every iteration
            line = readline(f)
            push!(lines, line)
        end
    end
    # process
    nlines = length(lines)
    res = Dict(
        "n_iter" => parse(Int64, lines[nlines-20][26:end]),
        "n_obj_func_eval" =>  parse(Int64, lines[nlines-10][56:end]),
        "n_obj_grad_eval" =>  parse(Int64, lines[nlines-9][56:end]),
        "n_neqc_eval" =>      parse(Int64, lines[nlines-8][56:end]),
        "n_nieqc_eval" =>     parse(Int64, lines[nlines-7][56:end]),
        "n_neqc_jac_eval" =>  parse(Int64, lines[nlines-6][56:end]),
        "n_nieqc_jac_eval" => parse(Int64, lines[nlines-5][56:end]),
        "n_lagrange_hessian" => parse(Int64, lines[nlines-4][56:end]),
        "t_cpu_sec_ipopt" => parse(Float64, lines[nlines-3][56:end]),
        "t_cpu_sec_nlp"   => parse(Float64, lines[nlines-2][56:end]),
        "EXIT"   => lines[nlines][7:end],
        "t_cpu" => parse(Float64, lines[nlines-2][56:end])+parse(Float64, lines[nlines-3][56:end]),
        "objective" => parse(Float64, lines[end-17][53:end]),
    )
    return lines, res
end

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
const TimeNS = typeof(time_ns())

struct RecordLocation
    pkgid::Union{UUID,Nothing}
    fullname::Tuple{Vararg{Symbol}}
    line::Int
    file::Union{Symbol,Nothing}
    tag::Union{Symbol,Nothing}
end

mutable struct Record
    location::RecordLocation
    total_times::Vector{typeof(Ref{TimeNS}())}
    total_calls::Vector{typeof(Ref{UInt}())}
end
# TODO: use cacheline-aligned state

function Record(location::RecordLocation)
    total_times = [Ref{TimeNS}(0) for _ in 1:Threads.nthreads()]
    total_calls = [Ref{UInt}(0) for _ in 1:Threads.nthreads()]
    return Record(location, total_times, total_calls)
end

function Base.empty!(record::Record)
    empty!(record.total_times)
    empty!(record.total_calls)
    # sizehint!(record.total_times, 0)
    # sizehint!(record.total_calls, 0)
    return record
end

const __RECORDS = Vector{Record}[]

"""
    StopWatches.clear()

Clear all records. The caller must ensure that no recordings are happening
concurrently.
"""
function clear()
    foreach(empty!, Iterators.flatten(__RECORDS))
    return
end

function init_records()
    empty!(__RECORDS)
    append!(__RECORDS, (Record[] for _ in 1:Threads.nthreads()))
end

@generated function _allocate_record(
    ::Val{pkgid},
    ::Val{fullname},
    ::Val{line},
    ::Val{file},
    ::Val{tag},
) where {pkgid,fullname,line,file,tag}
    location = RecordLocation(pkgid, fullname, line, file, tag)
    record = Record(location)
    # I am so sorry :)
    push!(__RECORDS[Threads.threadid()], record)
    quote
        $(QuoteNode(record))
    end
end
# TODO: Use StaticStorages.jl instead of the generated function hack.

function allocate_record_expr(__source__, __module__, tag::Union{Symbol,Nothing})
    pkgid = Base.PkgId(__module__).uuid
    fn = fullname(__module__)
    quote
        _allocate_record(
            $(Val(pkgid)),
            $(Val(fn)),
            $(Val(__source__.line)),
            $(Val(__source__.file)),
            $(Val(tag)),
        )
    end
end

macro allocate_record(tag::QuoteNode = QuoteNode(nothing))
    allocate_record_expr(__source__, __module__, tag.value)
end

# Used in test:
_sample_record1() = @allocate_record
_sample_record2() = @allocate_record :sample_record2
_sample_record3() = @allocate_record

"""
    @stopwatch(tag::Symbol, code)
    @stopwatch(code)

Record the start and stop times for executing `code`.

`tag` must be a literal symbol; i.e., `@stopwatch :mytag ...` instead of
`@stopwatch mytag ...`.

# Examples

```jldoctest
julia> using StopWatchesBase

julia> @stopwatch begin
           sleep(0.01)
           1
       end
1
```

The timing can now be printed with `using StopWatches` then
`StopWatches.summary()`.
"""
macro stopwatch(tag::QuoteNode, code)
    tag.value isa Symbol || error("`tag` must be a symbol; got: $(tag.value)")
    return stopwatch_impl(__source__, __module__, tag.value, code)
end

macro stopwatch(code)
    return stopwatch_impl(__source__, __module__, nothing, code)
end

macro stopwatch(tag::Symbol, code)
    error("use `@stopwatch :$tag ...`")
end

function stopwatch_impl(__source__, __module__, tag, code)
    record = allocate_record_expr(__source__, __module__, tag)
    return quote
        local t0 = time_ns()
        local ans = $(esc(Expr(:block, __source__, code)))
        local t1 = time_ns()
        local record = $record
        $push_record!(record, t1 - t0)
        ans
    end
end
# TODO: Add an option to use the geometric distribution trick to only measure
# some random subset?

function push_record!(record::Record, duration)
    # hope three is no yield :)
    tid = Threads.threadid()
    record.total_times[tid][] += duration
    record.total_calls[tid][] += 1
end
# TODO: Support other statistics? Use OnlineStats.jl?


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: marketplace_catalog
using AWS.Compat
using AWS.UUIDs
"""
    CancelChangeSet()

Used to cancel an open change request. Must be sent before the status of the request changes to APPLYING, the final stage of completing your change request. You can describe a change during the 60-day request history retention period for API calls.

# Required Parameters
- `catalog`: Required. The catalog related to the request. Fixed value: AWSMarketplace.
- `changeSetId`: Required. The unique identifier of the StartChangeSet request that you want to cancel.

"""

cancel_change_set(catalog, changeSetId; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("PATCH", "/CancelChangeSet", Dict{String, Any}("catalog"=>catalog, "changeSetId"=>changeSetId); aws_config=aws_config)
cancel_change_set(catalog, changeSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("PATCH", "/CancelChangeSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("catalog"=>catalog, "changeSetId"=>changeSetId), args)); aws_config=aws_config)

"""
    DescribeChangeSet()

Provides information about a given change set.

# Required Parameters
- `catalog`: Required. The catalog related to the request. Fixed value: AWSMarketplace
- `changeSetId`: Required. The unique identifier for the StartChangeSet request that you want to describe the details for.

"""

describe_change_set(catalog, changeSetId; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("GET", "/DescribeChangeSet", Dict{String, Any}("catalog"=>catalog, "changeSetId"=>changeSetId); aws_config=aws_config)
describe_change_set(catalog, changeSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("GET", "/DescribeChangeSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("catalog"=>catalog, "changeSetId"=>changeSetId), args)); aws_config=aws_config)

"""
    DescribeEntity()

Returns the metadata and content of the entity.

# Required Parameters
- `catalog`: Required. The catalog related to the request. Fixed value: AWSMarketplace
- `entityId`: Required. The unique ID of the entity to describe.

"""

describe_entity(catalog, entityId; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("GET", "/DescribeEntity", Dict{String, Any}("catalog"=>catalog, "entityId"=>entityId); aws_config=aws_config)
describe_entity(catalog, entityId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("GET", "/DescribeEntity", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("catalog"=>catalog, "entityId"=>entityId), args)); aws_config=aws_config)

"""
    ListChangeSets()

Returns the list of change sets owned by the account being used to make the call. You can filter this list by providing any combination of entityId, ChangeSetName, and status. If you provide more than one filter, the API operation applies a logical AND between the filters. You can describe a change during the 60-day request history retention period for API calls.

# Required Parameters
- `Catalog`: The catalog related to the request. Fixed value: AWSMarketplace

# Optional Parameters
- `FilterList`: An array of filter objects.
- `MaxResults`: The maximum number of results returned by a single call. This value must be provided in the next call to retrieve the next set of results. By default, this value is 20.
- `NextToken`: The token value retrieved from a previous call to access the next page of results.
- `Sort`: An object that contains two attributes, SortBy and SortOrder.
"""

list_change_sets(Catalog; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/ListChangeSets", Dict{String, Any}("Catalog"=>Catalog); aws_config=aws_config)
list_change_sets(Catalog, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/ListChangeSets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Catalog"=>Catalog), args)); aws_config=aws_config)

"""
    ListEntities()

Provides the list of entities of a given type.

# Required Parameters
- `Catalog`: The catalog related to the request. Fixed value: AWSMarketplace
- `EntityType`: The type of entities to retrieve.

# Optional Parameters
- `FilterList`: An array of filter objects. Each filter object contains two attributes, filterName and filterValues.
- `MaxResults`: Specifies the upper limit of the elements on a single page. If a value isn't provided, the default value is 20.
- `NextToken`: The value of the next token, if it exists. Null if there are no more results.
- `Sort`: An object that contains two attributes, SortBy and SortOrder.
"""

list_entities(Catalog, EntityType; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/ListEntities", Dict{String, Any}("Catalog"=>Catalog, "EntityType"=>EntityType); aws_config=aws_config)
list_entities(Catalog, EntityType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/ListEntities", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Catalog"=>Catalog, "EntityType"=>EntityType), args)); aws_config=aws_config)

"""
    StartChangeSet()

This operation allows you to request changes for your entities. Within a single ChangeSet, you cannot start the same change type against the same entity multiple times. Additionally, when a ChangeSet is running, all the entities targeted by the different changes are locked until the ChangeSet has completed (either succeeded, cancelled, or failed). If you try to start a ChangeSet containing a change against an entity that is already locked, you will receive a ResourceInUseException. For example, you cannot start the ChangeSet described in the example below because it contains two changes to execute the same change type (AddRevisions) against the same entity (entity-id@1).

# Required Parameters
- `Catalog`: The catalog related to the request. Fixed value: AWSMarketplace
- `ChangeSet`: Array of change object.

# Optional Parameters
- `ChangeSetName`: Optional case sensitive string of up to 100 ASCII characters. The change set name can be used to filter the list of change sets.
- `ClientRequestToken`: A unique token to identify the request to ensure idempotency.
"""

start_change_set(Catalog, ChangeSet; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/StartChangeSet", Dict{String, Any}("Catalog"=>Catalog, "ChangeSet"=>ChangeSet); aws_config=aws_config)
start_change_set(Catalog, ChangeSet, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_catalog("POST", "/StartChangeSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Catalog"=>Catalog, "ChangeSet"=>ChangeSet), args)); aws_config=aws_config)


text:
import Base: ==, convert, vec

import JuLIP
using JuLIP: Atoms, JVec, JMat, AbstractCalculator, mat, vecs,
             numbers, positions, cell, pbc

using JuLIP.MLIPs: IPBasis

import JuLIP: read_dict, write_dict

export LsqDB, Dat, basis, configs, config

"""
`Dat`: store one atomistic configuration and the associated
observations (e.g. E, F, V, ...)
"""
mutable struct Dat
   at::Atoms                         # configuration
   configtype::String                # group identifier
   D::Dict{String, Vector{Float64}}  # list of observations
   rows::Dict{String, Vector{Int}}   # row indices for LSQ system
   info::Dict{String, Any}           # anything else...
end

# should D.info also be compared?
==(d1::Dat, d2::Dat) = (
      (d1.configtype == d2.configtype) && (d1.D == d2.D) &&
      (d1.info == d2.info) &&
      all( f(d1.at) == f(d2.at) for f in (positions, numbers, cell, pbc) )
   )

function Dat(at::Atoms, config_type::AbstractString; kwargs...)
   dat = Dat(at, config_type, Dict{String, Vector{Float64}}(),
             Dict{String, Vector{Int}}(), Dict{String, Any}())
   for (key, val) in kwargs
      str_key = string(key)
      if !(ismissing(val) || (val == nothing))
         # convert the observation to a basic vector and store it
         dat.D[str_key] = vec_obs(str_key, val)
      end
   end
   return dat
end

write_dict(d::Dat) =
   Dict("__id__" => "IPFitting.Dat",
         "at" => write_dict(d.at),
         "configtype" => d.configtype,
         "D" => d.D,
         "rows" => d.rows,
         "info" => d.info)

function Dat(D::Dict)
   at = Atoms(D["at"])
   return Dat(at,
              D["configtype"],
              Dict{String, Any}(D["D"]),
              Dict{String, Any}(D["rows"]),
              Dict{String, Any}(D["info"]))
end

read_dict(::Val{Symbol("IPFitting.Dat")}, D::Dict) = Dat(D)

observation(d::Dat, key::String) = d.D[key]
hasobservation(d::Dat, key::String) = haskey(d.D, key)
observation(key::String, d::Dat) = d.D[key]

# -----------------------------------------------------------------

"""
`mutable struct LsqDB{TD, TB}`

A representation of a least-squares system stored on disk, which can
be extended by adding data, or basis functions
"""

mutable struct LsqDB
   basis::IPBasis
   configs::Vector{Dat}
   Ψ::Matrix{Float64}
   dbpath::String
end

basis(db::LsqDB) = db.basis
basis(db::LsqDB, i::Integer) = db.basis[i]
configs(db::LsqDB) = db.configs
config(db::LsqDB, i::Integer) = db.configs[i]

# prototypes

"""
convert some real data, in some generic format, into a vector to be stored
in a `Dat` or Lsq system. E.g.,
```
F = forces(...)::Vector{JVecF}
vec_obs(::Val{:F}, F) = mat(F)[:]
```
or equivalently, `vec_obs("F", F)`
"""
function vec_obs end

"""
convert a Vector{T} to some real (atomistic) data, e.g.,
```
x::Vector{Float64}
devec_obs(::Val{:F}, x) = vecs( resize(x, 3, length(x) ÷ 3) )
```
"""
function devec_obs end

"""
evaluate a specific observation type
"""
function eval_obs end

"""
create special weights for different observations
"""
weighthook(::Val, d::Dat) = 1.0
weighthook(s::String, d::Dat) = weighthook(Val(Symbol(s)), d)

"""
create special weights for different observations, similar to
`weighthooks` but this is used for computing errors rather than
for fitting
"""
err_weighthook(::Val, d::Dat) = 1.0
err_weighthook(s::String, d::Dat) = err_weighthook(Val(Symbol(s)), d)


text:
using SampleChains
using Test

@testset "SampleChains.jl" begin
    # Write your tests here.
end


text:
######################################################################
# plots.jl
######################################################################
using PyPlot
import PyPlot.plot

# ensures the fonts will be Type 1 if we save as pdf
# if not, fonts might be Type 3, which most conferences don't allow
rc("text", usetex=true)

# specifying font family and size desires
rc("font", family="Times New Roman", size=16)

# must be called before plot
function axis_font(af::Real)
    rc("font", size=af)
end

# include helper files
include("plot_trajectory.jl")

# export my functions
export plot, plot_trajectory, axis_font

# export PyPlot functions user might want
export figure, savefig, xlabel, ylabel, zlabel, title, hold
export xlim, ylim, zlim

"""
`plot(em::ErgodicManager, xd::VVF; alpha=1.0, cmap="Greys", show_score=true, lw=1.0, ms=9)`

`plot(em::ErgodicManager; alpha=1.0, cmap="Greys")`

The "Greys" cmap is dark where there is most density.
The "gray" cmap option is light where there is most density.

An `alpha` value closest to 1.0 is darker; less is more transparent.
"""
function plot(em::ErgodicManager, xd::VVF;
              alpha=1.0,
              cmap="Greys",
              show_score::Bool=true,
              ls::String="-",
              lw::Real=1.5,
              mew::Real=1,
              mfc::String="w",
              ms::Real=10,
              no_domain::Bool=false
             )

    # If it is in R3, let the trajectory know
    dims = 2
    if typeof(em) == ErgodicManagerR3
        dims = 3
    end

    # plot the trajectory
    plot_trajectory(xd, ls="-", lw=lw, mew=mew, mfc=mfc, ms=ms, dims=dims)

    # plot domain
    plot(em, alpha=alpha, cmap=cmap, no_domain=no_domain)

    # determines if ergodic score should be shown
    if show_score
        es = ergodic_score(em, xd)
        title_string = "es = $(round(es,digits=5))"
        title(title_string)
    end
end

# multi-agent version
# vtm is a vector of trajectory managers
function plot(em::ErgodicManager, xd::VVF, vtm::VTM;
              alpha=1.0,
              cmap="Greys",
              show_score::Bool=true,
              c::Vector{String}=ones(String, length(vtm)),
              lw::Float64=1.0,
              ms::Float64=8.0,
              no_domain::Bool=false
             )

    xds = vvf2vvvf(xd, vtm)
    num_agents = length(vtm)

    # if no colors are defined, just make them all blue
    if c[1] == ""
        for j = 1:num_agents
            c[j] = "b"
        end
    end

    # plot the trajectory
    for j = 1:num_agents
        xd = xds[j]
        plot_trajectory(xd, c=c[j], lw=lw, ms=ms)
    end

    plot(em, alpha=alpha, cmap=cmap, no_domain=no_domain)

    # determines if ergodic score should be shown
    #if show_score
    #	es = ergodic_score(em, xd)
    #	title_string = "es = $(round(es,5))"
    #	title(title_string)
    #end
end

# special plotting function for R2T
# input n is a Julia index; it goes from 1 to N+1 (instead of 0 to N)
function plot(em::ErgodicManager, xd::VVF, n::Int;
              alpha=1.0,
              cmap="Greys",
              show_score::Bool=true,
              lw::Float64=1.0,
              ms::Float64=6.0,
              no_domain::Bool=false
             )

    # plotting the trajectory
    dims = 2
    plot_trajectory(xd[1:n], lw=lw, ms=ms, dims=dims)

    # plotting the distribution
    a = [x_min(em), x_max(em), y_min(em), y_max(em)]
    imshow(em.phi', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
end

function plot(em::ErgodicManagerR2T, xd::VVF, n::Int;
              alpha=1.0,
              cmap="Greys",
              show_score::Bool=true,
              lw::Float64=1.0,
              mfc::String="w",
              ms::Float64=6.0,
              no_domain::Bool=false
             )

    # plotting the trajectory
    dims = 2
    plot_trajectory(xd[1:n], lw=lw, mfc=mfc, ms=ms, dims=dims)

    # plotting the distribution
    a = [x_min(em), x_max(em), y_min(em), y_max(em)]
    imshow(em.phi[:,:,n]', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
end
function plot(em::ErgodicManagerR2T, xd::VVF, n::Int, vtm::VTM;
              alpha=1.0,
              cmap="Greys",
              show_score::Bool=true,
              lw::Float64=1.0,
              ms::Float64=6.0,
              no_domain::Bool=false
             )

    # plotting the trajectory
    dims = 2
    plot_trajectory(xd[1:n], vtm, lw=lw, ms=ms, onlyMarks=onlyMarks, dims=dims)

    # plotting the distribution
    a = [x_min(em), x_max(em), y_min(em), y_max(em)]
    imshow(em.phi[:,:,n]', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
end

function plot(em::ErgodicManagerR2; alpha=1.0, cmap="Greys",no_domain=false)
    a = [x_min(em), x_max(em), y_min(em), y_max(em)]
    if no_domain
        rar = zeros(size(em.phi'))
        imshow(rar, interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    else
        imshow(em.phi', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    end
    axis(a)
    tick_params(direction="in")
end

# TODO: this is not ready to go yet
function plot(em::ErgodicManagerR3; alpha=1.0, cmap="Greys",no_domain=false)
    #if no_domain
    #	rar = zeros(size(em.phi'))
    #	imshow(rar, interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    #else
    #	imshow(em.phi', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    #end

    # this is super sketch but let's try it
    x_temp, y_temp, z_temp = get_pts2(em)
    scatter3D(x_temp,y_temp,z_temp)
    xlim(x_min(em), x_max(em))
    ylim(y_min(em), y_max(em))
    zlim(z_min(em), z_max(em))
end

function get_pts2(em::ErgodicManagerR3)
    w = Weights(vec(em.phi))
    nx,ny,nz = x_cells(em), y_cells(em), z_cells(em)
    x_temp, y_temp, z_temp = Float64[], Float64[], Float64[]
    xs,ys,zs = x_size(em), y_size(em), z_size(em)

    for i = 1:1000
        xi, yi, zi = ind2sub(size(em.phi), sample(w))
        push!(x_temp, (xi-.5) * xs)
        push!(y_temp, (yi-.5) * ys)
        push!(z_temp, (zi-.5) * zs)
    end
    return x_temp, y_temp, z_temp
end

function get_pts(em::ErgodicManagerR3)
    x_temp = Float64[]
    y_temp = Float64[]
    z_temp = Float64[]
    for xi = 1:x_cells(em)
        for yi = 1:y_cells(em)
            for zi = 1:z_cells(em)
                if em.phi[xi,yi,zi] > 0.5
                    push!(x_temp, (xi - .5)*x_size(em))
                    push!(y_temp, (yi - .5)*y_size(em))
                    push!(z_temp, (zi - .5)*z_size(em))
                end
            end
        end
    end
    return x_temp, y_temp, z_temp
end

function plot(em::ErgodicManagerSE2;
              alpha=1.0,
              cmap="Greys",
              no_domain=false
             )

    a = [x_min(em), x_max(em), y_min(em), y_max(em)]
    rar = size(em.phi)
    temp_phi = reshape(sum(em.phi,3), x_cells(em), y_cells(em))
    imshow(temp_phi', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    axis(a)
end

function plot(d::Domain, phi::Matrix{Float64}; alpha=1.0, cmap="Greys")
    a = [x_min(d), x_max(d), y_min(d), y_max(d)]
    imshow(phi', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
    axis(a)
    tick_params(direction="in")
end


text:
md"""
# Describing the Cuts using SVG

We use (locally modified, pull request pending) `NativeSVG.jl` to
generate SVG code.
"""

"""
The CSS stylesheet we use for SVG rendering in reports.
"""
const STYLESHEET = """
  g.everything {
  }
  g.everything * {
    vector-effect: non-scaling-stroke;
  }
  .cut {
    stroke-width: 1px;
    stroke: rgb(0%,50%, 50%);
    stroke-dasharray: 4 4;
  }
  .factory-edge {
    stroke-width: 1px;
    stroke: blue;
    fill: none;
  }
  .finished {
    stroke: none;
    fill: rgba(0%, 50%, 0%, 50%);
  }
  text.finished {
    color: white;
    text-anchor: middle;
    font-family: sans-serif;
    font-size: 2px;
    vector-effect: non-scaling-stroke;
  }
"""

#= this macro does not play well with Pluto notebooks
# Abstract measurements in svg user space from the dimensions we use
# for measuring panels:
@Unitful.unit svgd "svgd" SVGDistance 0.01u"inch" false
Unitful.register(@thismodule)
=#

"""
Turn a Unitful length quantity to a floating point number we can use in SVG.
"""
function svgdistance(d)::Real
    ustrip(Real, u"inch", d)
end

"""
Draw an SVG `rect` representing a panel.
"""
function panelrect(io::IO, panel::AbstractPanel, cssclass::String)
    # It's confusing that panel.width corresponds to SVG length
    # and panel.length corresponds to SVG width.  Sorry.
    # This is a consequence of the x and y coordinaces of a panel
    # corresponding with the panel's length and width respectively.
    g(io) do
	write(io, string("<!-- $(panel.label): ",
			 "$(panel.width) by $(panel.length), ",
			 "at $(panel.x), $(panel.y) -->\n"))
	rect(io; class=cssclass,
	     x=svgdistance(panel.x),
	     y=svgdistance(panel.y),
	     width=svgdistance(panel.length),
	     height=svgdistance(panel.width)) do
		 if panel isa FinishedPanel
		     title(io) do
			 print(io, "$(panel.length) × $(panel.width)")
		     end
		 end
	     end
	if panel isa FinishedPanel
	    NativeSVG.text(io;
			   class = cssclass,
			   x = svgdistance(panel.x + panel.length / 2),
			   y = svgdistance(panel.y + panel.width / 2)) do
			       write(io, "$(panel.length) × $(panel.width)")
		           end
	end

    end
end

export panelrect, panelrect

"""
The space between panels in an SVG drawing, and space between panels
and SVG edge.
"""
const SVG_PANEL_MARGIN = 2u"inch"

function toSVG(state::SearchState)
    buf = IOBuffer()
    toSVG(buf, state)
    return take!(buf)
end

export SVG_PANEL_MARGIN, toSVG

function toSVG(io::IO, state::SearchState)::Nothing
    rpg = makePanelGraph(state)
    write(io, """<?xml version="1.0" encoding="UTF-8"?>\n""")
    write(io, """<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"\n""")
    write(io, """          "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n""")
    # Outermost SVG:
    vpwidth = svgdistance(maximum(major.(keys(rpg))) + 2 * SVG_PANEL_MARGIN)
    vpheight = svgdistance(sum(minor.(filter(p -> p isa BoughtPanel, keys(rpg)))) +
	2 * SVG_PANEL_MARGIN)
    svg(io; xmlns="http://www.w3.org/2000/svg",
	width="90%",
	viewBox="0 0 $(vpwidth) $(vpheight)",
	style="background-color: pink") do
	    style(io; type="text/css") do
		write(io, STYLESHEET)
	    end
	    g(io; class="everything") do
		y = SVG_PANEL_MARGIN
		for stock in filter((p) -> p isa BoughtPanel, keys(rpg))
		    # We want to have the longer dimension of panel run
		    # horizontally.  If so, we can apply a 90 degree rotation.
		    # Here we just translate successive stock panels (BoughtPanel)
		    # by its minor dimension and margins to space them out.
		    # The toSVG method of BoughtPanel will deal with rotation.
		    tx = svgdistance(SVG_PANEL_MARGIN)
		    ty = svgdistance(y)
		    g(io; transform="translate($(tx), $(ty))") do
			toSVG(io, stock, rpg)
			y += minor(stock) + SVG_PANEL_MARGIN
		    end
		end
	    end
	end
end

function toSVG(io::IO, panel::AbstractPanel, rpg::PanelGraph)::Nothing
    # Do nothing, We only draw certain types of panel.
end

function toSVG(io::IO, panel::BoughtPanel, rpg::PanelGraph)::Nothing
    # We want to have the longer dimension of panel run horizontally.
    # This is already anticipated above wnere we calculate the SVG viewBox.
    transform = ""
    if panel.length != major(panel)
	tx = svgdistance(0u"inch")
	ty = svgdistance(panel.width)
	transform = "rotate(90) translate($tx $ty)"
    end
    g(io; class="BoughtPanel",
      transform=transform) do
	  panelrect(io, panel, "factory-edge")
	  for p in rpg[panel]
	      toSVG(io, p, rpg)
	  end
      end
end

function toSVG(io::IO, panel::Panel, rpg::PanelGraph)::Nothing
    g(io; class="Panel") do
	for p in rpg[panel]
	    toSVG(io, p, rpg)
	end
	endX = panel.x + panel.length
	endY = panel.y + panel.width
	if panel.cut_axis isa LengthAxis
	    startX = endX
	    startY = panel.y
	else
	    startX = panel.x
	    startY = endY
	end
	startX, startY, endX, endY = svgdistance.((startX, startY, endX, endY))
	d = "M $startX $startY L $endX $endY"
	NativeSVG.path(io; class="cut", d=d)
    end
end

function toSVG(io::IO, panel::FinishedPanel, rpg::PanelGraph)::Nothing
    g(io; class="FinishedPanel") do
	panelrect(io, panel, "finished")
    end
end


text:
"""
    NodalFieldModule

Module for nodal fields.
"""
module NodalFieldModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FieldModule: Field, nents
import FinEtools.FieldModule.@add_Field_fields

"""
    NodalField{T<:Number}

Nodal field.
"""
mutable struct NodalField{T<:Number} <: Field
    @add_Field_fields()
end

"""
    NodalField(data::FMat{T}=[]) where {T<:Number}

Constructor of nodal field.
"""
function NodalField(data::FMat{T}=[]) where {T<:Number}
    values = deepcopy(data)
    dofnums = 0*similar(values,FInt)
    is_fixed = similar(values,Bool)
    fill!(is_fixed, 0)
    fixed_values = zeros(T,size(values))
    nfreedofs = 0
    return NodalField(values, dofnums, is_fixed, fixed_values, nfreedofs)
end

function NodalField(data::FVec{T}) where {T<:Number}
    return NodalField(reshape(data, length(data), 1))
end

"""
    nnodes(self::NodalField)::FInt = nents(self)

Provide the number of nodes  in the nodal field.
"""
nnodes(self::NodalField)::FInt = nents(self)

end


text:
# tensortoolbox.jl
#
# Main file for module TensorToolbox, a Julia package for working with
# with tensors and tensor operations

module TensorToolbox

# Exports
#---------
# Types:
export VectorSpace, ElementarySpace, ElementaryHilbertSpace, EuclideanSpace
export ComplexSpace, CartesianSpace, GeneralSpace
export UnitaryRepresentationSpace, AbelianSpace
export U1Charge, ZNCharge, Parity
export CompositeSpace, ProductSpace, InvariantSpace
export IndexSpace, TensorSpace, AbstractTensor, Tensor, InvariantTensor
export TruncationScheme
export AbstractTensorMap
export SpaceError, IndexError

# general vector space methods
export space, issubspace, dim, dual, cnumber, iscnumber, directsum, fuse, basis
export ⊗, ℂ, ℝ # some unicode

# vector spaces with symmetries
export sectors, invariant

# tensor characteristics
export spacetype, tensortype, order
# tensor constructors
export tensor, tensorcat, to_range
# index manipulations
export insertind, deleteind, fuseind, splitind
# tensor operations
#export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct
#export tensorcopy!, tensoradd!, tensortrace!, tensorcontract!, tensorproduct!
# tensor factorizations
export leftorth, rightorth
export leftorth!, rightorth!, svd!

# truncation schemes
export notrunc, truncerr, truncdim, truncspace

# tensor maps
export domain, codomain, hermitian, posdef, tensormap

# tensor networks
export AbstractTensorNetwork, TensorNetwork
export network, optimizecontract

# General imports:
#------------------
importall Base
using Base.Cartesian

using TensorOperations
import TensorOperations.numind
import TensorOperations.scalar

# Exception types:
#------------------
abstract TensorException <: Exception

type SpaceError <: TensorException
    message::AbstractString
end
SpaceError()=SpaceError("Space mismatch")
# Exception type for all errors related to vector space mismatch

type IndexError <: TensorException
    message::AbstractString
end
IndexError()=IndexError("Invalid index specification")
# Exception type for all errors related to invalid tensor index specification.

# Types and methods for vector spaces
#-------------------------------------
abstract VectorSpace
space(V::VectorSpace) = V # just returns the space
issubspace(V1::VectorSpace,V2::VectorSpace) = V1==V2 # default, only identical spaces are subspaces
⊗(a,b,c...)=⊗(a,⊗(b,c...)) # introduce ⊗ as operator
*(V1::VectorSpace,V2::VectorSpace) = ⊗(V1,V2) # for convenience, product of vector spaces is tensor product

#include("basis.jl") # basis: iterator over a set of basis vectors spanning the space

include("elementaryspace.jl") # elementary finite-dimensional vector spaces
include("compositespace.jl") # composing elementary vector spaces

# Types and methods for tensors
#-------------------------------
import TensorOperations
#import TensorOperations: TCBuffer, defaultcontractbuffer
# intentionally shadow original TensorOperation methods for StridedArray objects

# define truncation schemes for tensors
include("tensors/truncation.jl")

# general definitions
include("tensors/abstracttensor.jl")

# Implementations:
include("tensors/tensor.jl") # generic tensor living in a ProductSpace without special properties
include("tensors/invarianttensor.jl") # generic tensor living in a ProductSpace without special properties

# Tensor networks: contract a network of tensors
#------------------------------------------------
include("tensornetworks/abstracttensornetwork.jl")
include("tensornetworks/tensornetwork.jl")
include("tensornetworks/optimize.jl")

# Tensor maps: linear maps acting on tensors
#--------------------------------------------
include("tensormaps/abstracttensormap.jl")
include("tensormaps/linearcombination.jl")
include("tensormaps/composition.jl")
include("tensormaps/shifted.jl")
include("tensormaps/dual.jl")
#include("tensormaps/adjoint.jl")
include("tensormaps/hermitian.jl") # guaranteed self adjoint operator
include("tensormaps/posdef.jl") # guaranteed positive definite operator
include("tensormaps/tensormap.jl") # a dense tensor map implemented using a tensor

end


text:
function subvar(sym, shift::Int)
    if shift == 0
        return sym
    elseif shift >= 1
        return Symbol(string(sym, "__f_", shift))
    else
        return Symbol(string(sym, "__b_", -shift))
    end
end

function walktree(tree, subvar, table)
#     if typeof(tree)==Symbol
#         return tree
    if typeof(tree)==Expr
        if (tree.head == :call) && (length(tree.args)==2)
            sym = tree.args[1]
            shift = tree.args[2]
            if (sym,shift) in keys(table)
                return subvar(sym,shift)
            else
                return Expr(tree.head, [walktree(e, subvar, table) for e in tree.args]...)
            end
        else
            return Expr(tree.head, [walktree(e, subvar, table) for e in tree.args]...)
        end
    else
        return tree
    end
end

function construct_table(symbols, api)
    d = Dict()
    for (symbol_group,shift) in api
        for ind in 1:length(symbols[symbol_group])
            variable = symbols[symbol_group][ind]
            if shift == 0
                t_symbol_group = symbol_group
            elseif shift == 1
                t_symbol_group = symbol(string(symbol_group,"_f"))
            elseif shift == -1
                t_symbol_group = symbol(string(symbol_group,"_b"))
            end
            d[(variable, shift)] = (t_symbol_group, ind)
        end
    end
    return d
end

function rewrite_function(funcode, symbols, api)

    new_args = Any[]
    for e in api
        if e[2] == 1
            push!(new_args, Symbol(string(e[1],:_f)))
        elseif e[2] == -1
            push!(new_args, Symbol(string(e[1],:_b)))
        else
            push!(new_args, e[1])
        end
    end

    # new_args = [e[1] for e in api]
    if length(funcode.args[1].args)>=2
        # println("Inserting")
        splice!(funcode.args[1].args,2:1, new_args)
    else
        append!(funcode.args[1].args,new_args)
    end

    funbody = funcode.args[2]
    table = construct_table(symbols, api)

    # # substitute_variable(expr) = walktree(expr, )

    for i = 1:length(funbody.args)
        funbody.args[i] = walktree(funbody.args[i], subvar, table)
    end

    new_definitions = Expr[]
    for ((symbol,tshift),(t_symbol_group,ind)) in construct_table(symbols, api)
        new_sym = (subvar(symbol, tshift))  # e__f_1
        new_ins = :($new_sym = $t_symbol_group[$ind])
        push!(new_definitions, new_ins)
    end

    prepend!(funbody.args, new_definitions)

    return funcode

end

macro quickdef(symbols, api, content)
    funcode = copy((content).args[2])
    new_funcode = rewrite_function(funcode, Core.eval(symbols), Core.eval(api))
    # if debug
        # println(new_funcode)
    # end
    Core.eval(new_funcode)
end

symbols = Dict(
    :states => [:a, :b],
    :controls => [:x, :y],
    :shocks => [:e],
    :parameters => [:p1]
)

recipes = Dict(
    :transition => [
        (:markov_states, -1),
        (:states, -1),
        (:controls, -1),
        (:markov_states, 0),
        (:parameters, 0),
#         (:states, 1),
    ],
    :lower_bound => [
        (:markov_states, 0),
        (:states, 0),
        (:parameters, 0),
    ],
    :upper_bound => [
        (:markov_states, -1),
        (:states, 0),
        (:parameters, 0),
    ],
    :arbitrage => [
        (:markov_states,0),
        (:states,0),
        (:controls, 0),
        (:markov_states,1),
        (:states,1),
        (:controls, 1),
        (:parameters, 0),
    ],
    :auxiliary => [
        (:markov_states,0),
        (:states,0),
        (:controls, 0),
        (:parameters, 0),
    ],
    :felicity => [
        (:markov_states, 0),
        (:states, 0),
        (:controls, 0),
        (:parameters, 0)
    ]
)

#api = [
#    (:states, -1),
#    (:controls, -1),
#    (:shocks, 0),
#    (:parameters, 0),
#    (:states, 1),
#]

#@quickdef symbols api begin
#function transition()
#        a(1) = a(-1) + b(-1)
#         return out
#end
#end

max_issues_repo_name: EconForge/Dolang