text:
using Revise
using AlphaZero
using ProgressLogging
using TensorBoardLogger
using Logging
using ParameterSchedulers: Scheduler, Cos
import Flux
import FileIO
import Random
import CUDA

import LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

include("../mu_game_wrapper.jl")
# include("../architectures/simplemlp.jl")
# include("../architectures/resnet.jl")
include("../network.jl")

# include("../alphazerolike.jl")
include("../trace.jl")
include("../play.jl")
include("../simulations.jl")
include("../training.jl")
include("../learning.jl")
include("../benchmark.jl")
include("../probe-games.jl")

# 3x Faster than julia default implementation #! may break something
Base.hash(x::Array{Float32,3}) = hash(vec(x))
Base.hash(x::Tuple{Array{Float32,3},Int}) = hash(x[1]) >> x[2]

device = Flux.cpu
# device = Flux.gpu
# CUDA.allowscalar(false)

gspec = Examples.games["connect-four"]

##

n=1
self_play = (;
  sim=SimParams(
    num_games=100,
    num_workers=1,
    # TODO make batch sizes for each network
    batch_size=1,
    use_gpu=true,
    # todevice=Flux.cpu, #?
    reset_every=4, #not used, mcts resets everytime
    flip_probability=0.,
    alternate_colors=false),
  mcts = MctsParams(
    num_iters_per_turn=128, #1000 benchmark
    cpuct=1.25,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.25,
    dirichlet_noise_α=0.1),
  device=device)

arena = (;
  sim=SimParams(
    num_games=3,
    num_workers=3,
    batch_size=2,
    use_gpu=true,
    reset_every=1,
    flip_probability=0., #0.5
    alternate_colors=true),
  mcts = MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.3),
    dirichlet_noise_ϵ=0.1),
  update_threshold=0.4,
  device=device)

learning_params = (;
  num_unroll_steps=10, #if =0, g is not learning, muzero-general=20
  td_steps=50, # with max length=9, always go till the end of the game, rootvalues don't count
  discount=0.997,
  #// value_loss_weight = 0.25, #TODO
  l2_regularization=1f-4, #Float32
  #// l2_regularization=0f0, #Float32
  loss_computation_batch_size=1024,
  batches_per_checkpoint=10,
  num_checkpoints=1,
  opt=Scheduler(
    Cos(λ0=5e-3, λ1=1e-5, period=10000), # cosine annealing, google 2e4, generat doesn't use any
    Flux.ADAM()
  ),
  model_type=:mlp,
  device=device
)

benchmark_sim = SimParams(
    num_games=400,
    num_workers=64,
    batch_size=32,
    use_gpu=true,
    reset_every=1,
    flip_probability=0.0, #0.5
    alternate_colors=true)

bench_mcts = MctsParams(
  num_iters_per_turn=500, #1000 benchmark
  cpuct=1.25,
  temperature=ConstSchedule(0.3),
  dirichlet_noise_ϵ=0.25,
  dirichlet_noise_α=0.1)

benchmark = (;
  vanilla_mcts=Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MctsRollouts(bench_mcts),
    benchmark_sim),
  # minmax_d5=Benchmark.Duel(
  #   Mu(arena.mcts),
  #   Benchmark.MinMaxTS(depth=5, amplify_rewards=true, τ=1.),
  #   benchmark_sim)
)

  μparams = MuParams(self_play, learning_params, arena, 1_000, 10_000)

  hs_shape = 128
# muzero-general params
simμNetworkHP = MuNetworkHP(
  gspec,
  PredictionSimpleHP(hiddenstate_shape=hs_shape, width=256, depth_common=0,
    depth_vectorhead=1, depth_scalarhead=0, use_batch_norm=true, batch_norm_momentum=0.8f0),
  DynamicsSimpleHP(hiddenstate_shape=hs_shape, width=256, depth_common=-1,
   depth_vectorhead=1, depth_scalarhead=0,  use_batch_norm=true, batch_norm_momentum=0.8f0),
  RepresentationSimpleHP(hiddenstate_shape=hs_shape, width=0, depth=-1))

  num_filters=32
resμNetworkHP = MuNetworkHP(
  gspec,
  PredictionResnetHP(
    num_blocks=3,
    num_filters=num_filters,
    conv_kernel_size=(3,3),
    num_policy_head1_filters=2,
    num_value_head2_filters=2,
    batch_norm_momentum=0.6f0
  ),
  DynamicsResnetHP(
    num_blocks=3,
    num_filters=num_filters,
    conv_kernel_size=(3,3),
    num_reward_head2_filters=32,
    batch_norm_momentum=0.6f0
  ),
  RepresentationResnetHP(
    num_blocks=2,
    num_filters=num_filters,
    conv_kernel_size=(3,3),
    batch_norm_momentum=0.6f0
  )
)

  envsim = MuEnv(gspec, μparams, MuNetwork(simμNetworkHP))
  envres = MuEnv(gspec, μparams, MuNetwork(resμNetworkHP))
  # env = MuEnv(gspec, μparams, deepcopy(env.bestnns))
  # env = MuEnv(gspec, μparams, MuNetwork(μNetworkHP), experience=FileIO.load("results/c4_5x4/memory_mctsrollout500.jld2", "mem"))
  # env = MuEnv(gspec, μparams, μNetwork, experience=FileIO.load("results/c4_5x4/memory_mctsrollout500.jld2", "mem"))
##
|
# cpu
@timed self_play_step!(envsim)
@timed learning_step!(envsim)

path = "results/c4/" * format(now(),"yyyy-mm-ddTHHMM") * "_custominit_batch256/"
tblogger=TBLogger(path, min_level=Logging.Info)

with_logger(tblogger) do
  @info "params" params=μparams
  # @info "alphalike params" αnns.f.net.hyper
  # @timed learning_step!(env)
  @info "network params" resμNetworkHP
  # @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
  @info "Benchmark" run_duel(envres, benchmark)...
  mutrain!(envsim, benchmark=benchmark, path=path)
end

# # simple test
# game = GI.init(gspec)

# GI.play!(game, 3)
# # hiddenstate = h(game)
# onehotaction = GI.encode_action(gspec,a)
# cat(hiddenstate, onehotaction, dims=3)

# state1 = GI.current_state(game)
# GI.play!(game, 3)
# state2 = GI.current_state(game)
# GI.play!(game, 3)
# state3 = GI.current_state(game)

# states = [state1, state2, state3]

# vstates = Flux.batch(GI.vectorize_state(gspec, s) for s in states)
# hiddenstates = forward(h, vstates)

# actions = [3,3,1]
# onehotactions = Flux.batch(GI.encode_action(gspec, a) for a in actions)

# g((hiddenstates[:,:,:,1],3))
# hiddenstates[:,:,:,1]
# stateactions = cat(hiddenstates, onehotactions, dims=3)
# forward(g, stateactions)

# forward(f, hiddenstates)

batch_r = [([9.103096f37, 6.4342253f37, 1.2172573f37, -5.9345403f37, 4.335317f37, -1.922554f37, -1.1086562f37, 5.4448747f37, -4.133477f37, -4.6036476f37, -5.2983734f37, -1.0046158f37, 2.0917472f37, 5.018153f37, 6.926142f36, 5.2010117f37, -4.739836f37, -6.837085f35, 2.9408774f37, 8.1636866f37, 5.5545777f37, 5.1190206f37, -3.3175674f37, 2.8772867f36, 6.9162795f36, -3.0321363f37, -1.3474731f38, -2.6346488f37, -1.7418001f37, 3.6669965f37, 4.0769748f37, 5.642504f37, -4.3100587f37, -4.536272f37, 9.1013085f36, 7.2361563f37, 6.699831f37, 5.213427f37, 6.788594f36, 2.5127895f37, -3.7990249f37, 1.6626149f37, -8.751232f36, 6.1124434f37, 2.8890792f37, -5.3267003f37, 6.6228264f37, -4.0651613f37, -8.750854f37, 1.0138958f38, -5.0922444f36, 5.532515f37, 2.0873239f37, -6.4766546f37, 2.3834298f37, 1.0271321f37, 2.3720243f36, 4.5693703f37, -3.903427f37, -5.037562f37, -7.147416f37, 3.9698593f37, -1.0508346f37, -3.3813056f37, 3.5642704f37, -7.0731227f37, 1.851113f37, 1.7753027f37, 1.5221927f37, -1.3454395f37, 1.0799189f37, -5.2645626f37, 9.382386f37, -6.22857f37, -4.6733255f36, -9.5401476f36, 1.2582538f37, -5.741112f37, -1.6655605f37, -1.3851418f37, 1.2454471f37, 3.844346f37, -1.1275677f37, 3.9049994f37, -8.1779766f37, 3.785576f37, -2.2327318f37, -8.806831f37, 7.617063f37, -1.7628736f37, -2.319465f37, -1.1858449f38, -3.098632f37, 2.773229f37, 5.9916525f37, 3.9901686f37, 2.7722208f37, 3.9903658f37, 2.4513655f37, 6.464902f36, 2.2294734f37, 6.5114775f37, 6.6826453f37, -7.40496f37, -2.9538444f37, -1.0950436f37, -1.7714793f37, 1.3216431f37, -1.1667053f38, -7.4938174f37, 4.5241126f37, 1.0923443f37, -2.8192844f35, 7.2007087f37, 1.03249977f37, 3.3486109f37, 1.427612f37, -1.7262653f37, 5.0123397f37, 9.7588985f36, -2.2738784f37, -8.310829f37, 8.0603386f37, 6.633395f37, 5.777899f37, 8.0670303f37, -7.0415035f37, 5.5774223f37], 1)] #([5.848612f35, 4.1338997f35, 7.820708f34, -3.8128573f35, 2.7853797f35, -1.2352136f35, -7.122959f34, 3.498256f35, -2.655701f35, -2.9577786f35, -3.4041307f35, -6.454512f34, 1.343918f35, 3.2240936f35, 4.449947f34, 3.341577f35, -3.0452786f35, -4.3927318f33, 1.8894722f35, 5.2450545f35, 3.568738f35, 3.2888985f35, -2.1314901f35, 1.8486182f34, 4.443612f34, -1.9481046f35, -8.657326f35, -1.6927247f35, -1.119082f35, 2.3559934f35, 2.6193991f35, 3.6252284f35, -2.7691518f35, -2.914491f35, 5.847462f34, 4.649129f35, 4.304549f35, 3.349554f35, 4.361579f34, 1.6144316f35, -2.4408196f35, 1.0682068f35, -5.6225434f34, 3.9271586f35, 1.8561924f35, -3.4223306f35, 4.255073f35, -2.611809f35, -5.6223014f35, 6.514136f35, -3.2716955f34, 3.5545624f35, 1.3410764f35, -4.161159f35, 1.5313202f35, 6.5991766f34, 1.523997f34, 2.9357568f35, -2.5078966f35, -3.236563f35, -4.592114f35, 2.5505786f35, -6.751465f34, -2.172442f35, 2.2899934f35, -4.5443816f35, 1.18931425f35, 1.1406067f35, 9.779875f34, -8.644255f34, 6.9383266f34, -3.3824075f35, 6.0280498f35, -4.0017685f35, -3.0025483f34, -6.12941f34, 8.084105f34, -3.6885824f35, -1.0700993f35, -8.899341f34, 8.001822f34, 2.4699382f35, -7.244465f34, 2.5089071f35, -5.2542347f35, 2.43218f35, -1.4344985f35, -5.658263f35, 4.8938567f35, -1.1326215f35, -1.4902236f35, -7.618887f35, -1.9908276f35, 1.7817602f35, 3.8495534f35, 2.5636274f35, 1.7811127f35, 2.563754f35, 1.5749673f35, 4.153608f34, 1.4324054f35, 4.1835322f35, 4.2935052f35, -4.7575842f35, -1.8978035f35, -7.0354965f34, -1.1381512f35, 8.491371f34, -7.495921f35, -4.814673f35, 2.9066789f35, 7.0181594f34, -1.8113296f33, 4.6263535f35, 6.6336646f34, 2.1514356f35, 9.172204f34, -1.1091012f35, 3.220357f35, 6.269955f34, -1.4609349f35, -5.3395903f35, 5.178655f35, 4.2618623f35, 3.7122193f35, 5.1829523f35, -4.5240683f35, 3.583416f35], 1)]

recur = RecurrentOracle(envsim.curnns)
@enter evaluate_batch(recur, batch_r)
evaluate_batch(recur, batch_r)

S_A = Float32[9.103096f37, 6.4342253f37, 1.2172573f37, -5.9345403f37, 4.335317f37, -1.922554f37, -1.1086562f37, 5.4448747f37, -4.133477f37, -4.6036476f37, -5.2983734f37, -1.0046158f37, 2.0917472f37, 5.018153f37, 6.926142f36, 5.2010117f37, -4.739836f37, -6.837085f35, 2.9408774f37, 8.1636866f37, 5.5545777f37, 5.1190206f37, -3.3175674f37, 2.8772867f36, 6.9162795f36, -3.0321363f37, -1.3474731f38, -2.6346488f37, -1.7418001f37, 3.6669965f37, 4.0769748f37, 5.642504f37, -4.3100587f37, -4.536272f37, 9.1013085f36, 7.2361563f37, 6.699831f37, 5.213427f37, 6.788594f36, 2.5127895f37, -3.7990249f37, 1.6626149f37, -8.751232f36, 6.1124434f37, 2.8890792f37, -5.3267003f37, 6.6228264f37, -4.0651613f37, -8.750854f37, 1.0138958f38, -5.0922444f36, 5.532515f37, 2.0873239f37, -6.4766546f37, 2.3834298f37, 1.0271321f37, 2.3720243f36, 4.5693703f37, -3.903427f37, -5.037562f37, -7.147416f37, 3.9698593f37, -1.0508346f37, -3.3813056f37, 3.5642704f37, -7.0731227f37, 1.851113f37, 1.7753027f37, 1.5221927f37, -1.3454395f37, 1.0799189f37, -5.2645626f37, 9.382386f37, -6.22857f37, -4.6733255f36, -9.5401476f36, 1.2582538f37, -5.741112f37, -1.6655605f37, -1.3851418f37, 1.2454471f37, 3.844346f37, -1.1275677f37, 3.9049994f37, -8.1779766f37, 3.785576f37, -2.2327318f37, -8.806831f37, 7.617063f37, -1.7628736f37, -2.319465f37, -1.1858449f38, -3.098632f37, 2.773229f37, 5.9916525f37, 3.9901686f37, 2.7722208f37, 3.9903658f37, 2.4513655f37, 6.464902f36, 2.2294734f37, 6.5114775f37, 6.6826453f37, -7.40496f37, -2.9538444f37, -1.0950436f37, -1.7714793f37, 1.3216431f37, -1.1667053f38, -7.4938174f37, 4.5241126f37, 1.0923443f37, -2.8192844f35, 7.2007087f37, 1.03249977f37, 3.3486109f37, 1.427612f37, -1.7262653f37, 5.0123397f37, 9.7588985f36, -2.2738784f37, -8.310829f37, 8.0603386f37, 6.633395f37, 5.777899f37, 8.0670303f37, -7.0415035f37, 5.5774223f37, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

@enter forward(envsim.curnns.g, S_A)

c = Float32[9.103096f37, 6.4342253f37, 1.2172573f37, -5.9345403f37, 4.335317f37, -1.922554f37, -1.1086562f37, 5.4448747f37, -4.133477f37, -4.6036476f37, -5.2983734f37, -1.0046158f37, 2.0917472f37, 5.018153f37, 6.926142f36, 5.2010117f37, -4.739836f37, -6.837085f35, 2.9408774f37, 8.1636866f37, 5.5545777f37, 5.1190206f37, -3.3175674f37, 2.8772867f36, 6.9162795f36, -3.0321363f37, -1.3474731f38, -2.6346488f37, -1.7418001f37, 3.6669965f37, 4.0769748f37, 5.642504f37, -4.3100587f37, -4.536272f37, 9.1013085f36, 7.2361563f37, 6.699831f37, 5.213427f37, 6.788594f36, 2.5127895f37, -3.7990249f37, 1.6626149f37, -8.751232f36, 6.1124434f37, 2.8890792f37, -5.3267003f37, 6.6228264f37, -4.0651613f37, -8.750854f37, 1.0138958f38, -5.0922444f36, 5.532515f37, 2.0873239f37, -6.4766546f37, 2.3834298f37, 1.0271321f37, 2.3720243f36, 4.5693703f37, -3.903427f37, -5.037562f37, -7.147416f37, 3.9698593f37, -1.0508346f37, -3.3813056f37, 3.5642704f37, -7.0731227f37, 1.851113f37, 1.7753027f37, 1.5221927f37, -1.3454395f37, 1.0799189f37, -5.2645626f37, 9.382386f37, -6.22857f37, -4.6733255f36, -9.5401476f36, 1.2582538f37, -5.741112f37, -1.6655605f37, -1.3851418f37, 1.2454471f37, 3.844346f37, -1.1275677f37, 3.9049994f37, -8.1779766f37, 3.785576f37, -2.2327318f37, -8.806831f37, 7.617063f37, -1.7628736f37, -2.319465f37, -1.1858449f38, -3.098632f37, 2.773229f37, 5.9916525f37, 3.9901686f37, 2.7722208f37, 3.9903658f37, 2.4513655f37, 6.464902f36, 2.2294734f37, 6.5114775f37, 6.6826453f37, -7.40496f37, -2.9538444f37, -1.0950436f37, -1.7714793f37, 1.3216431f37, -1.1667053f38, -7.4938174f37, 4.5241126f37, 1.0923443f37, -2.8192844f35, 7.2007087f37, 1.03249977f37, 3.3486109f37, 1.427612f37, -1.7262653f37, 5.0123397f37, 9.7588985f36, -2.2738784f37, -8.310829f37, 8.0603386f37, 6.633395f37, 5.777899f37, 8.0670303f37, -7.0415035f37, 5.5774223f37, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

s1 = envsim.curnns.g.statehead[1](c)
s2 = envsim.curnns.g.statehead[2](s1)
dense2 = envsim.curnns.g.statehead[2]

dense2.weight[:,1]
s1


text:
@with_kw mutable struct History
    # Number of completed training epochs
    epochs::Int = 0
    # total number of completed steps
    steps::Int = 0
    # number of completed steps in current epoch
    stepsepoch::Int = 0
end

"""
    Recorder()

Maintains a [`History`](#). It's stored in `learner.cbstate.history`.
"""
struct Recorder <: Callback end

stateaccess(::Recorder) = (
    cbstate = (history = Write(),),
    step = Read(),
    )

function init!(::Recorder, learner)
    if !haskey(learner.cbstate, :history)
        learner.cbstate.history = DefaultDict{Phase, History}(() -> History())
    end
end

function on(::EpochBegin, phase::Phase, recorder::Recorder, learner)
    learner.cbstate.history[phase].stepsepoch = 0
end

function on(::StepEnd, phase::Phase, recorder::Recorder, learner)
    history = learner.cbstate.history[phase]
    history.steps += 1
    history.stepsepoch += 1
end

function on(::EpochEnd, phase::Phase, recorder::Recorder, learner)
    history = learner.cbstate.history[phase]
    history.epochs += 1
end

ext: jl


text:
"""
    point_source(medium::Acoustic, source_position, amplitude=1)::RegularSource

Create 2D [`Acoustic`](@ref) point [`RegularSource`](@ref) (zeroth Hankel function of first type)
"""
function point_source(medium::Acoustic{T,2}, source_position::AbstractVector, amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,2}} where T <: AbstractFloat

    # Convert to SVector for efficiency and consistency
    source_position = SVector{2,T}(source_position)

    if typeof(amplitude) <: Number
        amp(ω) = amplitude
    else
        amp = amplitude
    end
    source_field(x,ω) = (amp(ω)*im)/4 * hankelh1(0,ω/medium.c * norm(x-source_position))

    function source_coef(order,centre,ω)
        k = ω/medium.c
        r, θ = cartesian_to_radial_coordinates(centre - source_position)

        # using Graf's addition theorem
        return (amp(ω)*im)/4 * [hankelh1(-n,k*r) * exp(-im*n*θ) for n = -order:order]
    end

    return RegularSource{T,Acoustic{T,2},WithoutSymmetry{2}}(medium, source_field, source_coef)
end

# If we replaced 3 with Dim below this could should work for all dimensions! Test carefully after changing.
function point_source(medium::Acoustic{T,3}, source_position, amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,3}} where T <: AbstractFloat

    # Convert to SVector for efficiency and consistency
    source_position = SVector{3,T}(source_position)

    if typeof(amplitude) <: Number
        amp(ω) = amplitude
    else
        amp = amplitude
    end

    # source_field(x,ω) = amp(ω) / (T(4π) * norm(x-source_position)) * exp(im * ω/medium.c * norm(x-source_position))
    # source_field(x,ω) = amp(ω)/sqrt(4π) * shankelh1(0, ω/medium.c * norm(x-source_position))
    source_field(x,ω) = amp(ω) * outgoing_basis_function(medium, ω)(0,x-source_position)[1]

    # centre = 12.0 .* rand(3)
    # x = centre + rand(3)

    # order = 5
    # U = outgoing_translation_matrix(medium, order, ω,  centre);
    # vs = regular_basis_function(medium, ω)(order, x - centre);
    # us = outgoing_basis_function(medium, ω)(0, x);
    # sum(U[1,:] .* vs) - us[1]

    function source_coef(order,centre,ω)
        U = outgoing_translation_matrix(medium, order, ω,  centre);
        return amp(ω) * U[1,:]
    end

    return RegularSource{T,Acoustic{T,3},WithoutSymmetry{3}}(medium, source_field, source_coef)
end

function plane_source(medium::Acoustic{T,Dim}; position::AbstractArray{T} = SVector(zeros(T,Dim)...),
        direction = SVector(one(T), zeros(T,Dim-1)...),
        amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,Dim}} where {T, Dim}

    plane_source(medium, position, direction, amplitude)
end

"""
    plane_source(medium::Acoustic, source_position, source_direction=[1,0], amplitude=1)::RegularSource

Create an [`Acoustic`](@ref) planar wave [`RegularSource`](@ref)
"""
function plane_source(medium::Acoustic{T,2}, position::AbstractArray{T}, direction::AbstractArray{T} = SVector(one(T),zero(T)), amplitude::Union{T,Complex{T}} = one(T))::RegularSource{T,Acoustic{T,2}} where {T}

    # Convert to SVector for efficiency and consistency
    position = SVector(position...)
    direction = SVector((direction ./ norm(direction))...) # unit direction

    S = (abs(dot(direction,azimuthalnormal(2))) == one(T)) ? PlanarAzimuthalSymmetry{2} : PlanarSymmetry{2}

    # This pseudo-constructor is rarely called, so do some checks and conversions
    if iszero(norm(direction))
        throw(DomainError("RegularSource direction must not have zero magnitude."))
    end

    if typeof(amplitude) <: Number
        amp(ω) = amplitude
    else
        amp = amplitude
    end

    source_field(x,ω) = amp(ω)*exp(im*ω/medium.c*dot(x-position, direction))

    function source_coef(order,centre,ω)
        # Jacobi-Anger expansion
        θ = atan(direction[2],direction[1])
        source_field(centre,ω) * [exp(im * n *(T(pi)/2 -  θ)) for n = -order:order]
    end

    return RegularSource{T,Acoustic{T,2},S}(medium, source_field, source_coef)
end

function plane_source(medium::Acoustic{T,3}, position::AbstractArray{T}, direction::AbstractArray{T} = SVector(zero(T),zero(T),one(T)), amplitude::Union{T,Complex{T}} = one(T)) where {T}

    # Convert to SVector for efficiency and consistency
    position = SVector(position...)
    direction = SVector( (direction ./ norm(direction))...) # unit direction

    S = (abs(dot(direction,azimuthalnormal(3))) == one(T)) ? PlanarAzimuthalSymmetry{3} : PlanarSymmetry{3}

    # This pseudo-constructor is rarely called, so do some checks and conversions
    if iszero(norm(direction))
        throw(DomainError("RegularSource direction must not have zero magnitude."))
    end

    if typeof(amplitude) <: Number
        amp(ω) = amplitude
    else
        amp = amplitude
    end

    source_field(x,ω) = amp(ω) * exp(im * ω/medium.c * dot(x-position, direction))

    function source_coef(order,centre,ω)
        # plane-wave expansion for complex vectors
        r, θ, φ  = cartesian_to_radial_coordinates(direction)
        Ys = spherical_harmonics(order, θ, φ)
        lm_to_n = lm_to_spherical_harmonic_index

        return T(4pi) * source_field(centre,ω) .*
        [
            Complex{T}(im)^l * (-one(T))^m * Ys[lm_to_n(l,-m)]
        for l = 0:order for m = -l:l]
    end

    return RegularSource{T,Acoustic{T,3},S}(medium, source_field, source_coef)
end

function regular_spherical_coefficients(psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}) where {T,Dim}

    source = plane_source(psource.medium;
        amplitude = psource.amplitude[1],
        position = psource.position,
        direction = psource.direction
    )

    source.coefficients
end


text:
function _evolution(mdl, body)
    nt = Base.return_types(mdl, ())[1]
    mdl_kwargs = fieldnames(nt)
    mdl_type = typeof(mdl)
    callex = Expr(:call, :(AnimalBehavior.evol!), Expr(:parameters, mdl_kwargs...), :(M::T), :s, :a, :r)
    whereex = Expr(:where, callex, :(T<:$mdl_type))
    ex = Expr(:(=), whereex, body)
    return ex
end

macro evolution(mdl, expr)
    body = QuoteNode(expr)
    return esc(quote
        eval(AnimalBehavior._evolution($mdl, $body))
    end)
end

function evol! end

function _observation(mdl, body)
    nt = Base.return_types(mdl, ())[1]
    mdl_kwargs = fieldnames(nt)
    mdl_type = typeof(mdl)
    callex = Expr(:call, :(AnimalBehavior.observ), Expr(:parameters, mdl_kwargs...), :(M::T), :s)
    whereex = Expr(:where, callex, :(T<:$mdl_type))
    ex = Expr(:(=), whereex, body)
    return ex
end

macro observation(mdl, expr)
    body = QuoteNode(expr)
    return esc(quote
        eval(AnimalBehavior._observation($mdl, $body))
    end)
end

function observ end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
@test isapprox(AutoCore._mod2pi2(0.0), 0.0)
@test isapprox(AutoCore._mod2pi2(0.5), 0.5)
@test isapprox(AutoCore._mod2pi2(2pi + 1), 1.0)
@test isapprox(AutoCore._mod2pi2(1 - 2pi), 1.0)


text:
include("Hopfield-model.jl")
include("Analysis_tools.jl")

function equiliabration(J,xi,β,steps)  # equilibration of a single pattern in a given network J
    for ii in 1:steps

        switchup=rand(1:length(xi)) # choose a given spinflip

        if (exp(-β*(2*(xi[ switchup])*(transpose(J[:,switchup])*xi)/length(xi)))>rand()) # test the energy diffrence from this spinflip

            xi[switchup]*=-1
        end
    end
        return xi
    end

    function Vec_equiliabration_compartment(J,xi,β,βeq,steps)   # equilibration of a set of patterns for multiple compartments
         xi_final=copy(xi)
        for ii in 1:length(xi_final[1,:])

            xi_final[:,ii]=equiliabration(J[:,:,givecompartment(J,xi[:,ii],β)],xi[:,ii],βeq,steps)
        end
        return xi_final
    end

    function Vec_equiliabration(J,xi,β,steps)  # equilibration of a set of patterns for a single network
         xi_final=copy(xi)
        for ii in 1:length(xi_final[1,:])
            xi_final[:,ii]=equiliabration(J,xi[:,ii],β,steps)
        end
        return xi_final
    end

function energies_to_compartments(J,xi)  # energies for the C compartments (compartments)
    es=Vector{Float64}()
    for ii in 1:length(J[1,1,:])
        append!(es,energie(J[:,:,ii],xi))
    end
    return es
end

function givecompartment(J,xi,β) # choose a compartment to update/ find stored state
    probs=exp.(β .*energies_to_compartments(J,xi))
    d = Categorical(probs./sum(probs))
    return rand(d)
end

function givecompartment_stats(J,xi,β) # choose a compartment to update/ find stored state and return the energy
    es=energies_to_compartments(J,xi)
    probs=exp.(β .* es)
    d = Categorical(probs./sum(probs))
    pos=rand(d)
    return pos,es[pos]
end

function update_compartments(J,xi,γ,β) # update a system with many compartments
    compartment=givecompartment(J,xi,β)    # find the compartment
    J[:,:,compartment]=new_J(J[:,:,compartment],xi,γ) # update the compartment
    return J
end

function update_compartments_stats(J,xi,γ,β) # update a system with many compartments and give the energy

    compartment,e_take=givecompartment_stats(J,xi,β)

    J[:,:,compartment]=new_J(J[:,:,compartment],xi,γ)
    return J,compartment,e_take
end

function start_compartment(xi,γ,compartments) # initiate system with multiple compartmetns only intitial state J(0)
    l=length(xi[:,1])
    J=zeros(l,l,compartments)
    for ii in 1:compartments
        J[:,:,ii]=start_J(xi[:,ii],1- 0.95*compartments/length(xi[1,:]) )
    end
    for ii in compartments:length(xi[1,:])
        b=convert(Int64,((ii+(compartments-1)) % compartments) +1 )
        J[:,:,b]=new_J(J[:,:,b],xi[:,ii],1- 0.95*compartments/length(xi[1,:]) )
    end
    return J
end

function build_compartment(l,n,x,γ,β,runin,compartments) # give states and system with multiple compartments after inital steps
        if (compartments > n)
            println(" to many compartments the maximum is the number of patters here $(n) you have selected $(compartments) ")
            return -1
        end

        xi= get_patterns(l,n,0.5)
            xi,J = let_run_in_compartment(xi,γ,x,β,runin,n,compartments)
    return J,xi
end

function let_run_in_compartment(xi,γ,x,β,runin,n,compartments) # initian of system and patterns and making it independed of starting conditions

        J= start_compartment(xi,γ,compartments)

    for ii in 1:runin
        activestate=updateactivestate(n)
        xi[:]=new_XI_rand(xi[:], x)
        J=update_compartments(J,xi[:,activestate],γ,β)
    end
    return xi,J
end

function get_false_compartment_withe_compartment(l,n,x,γ,β,runin,compartments,steps,reps)

    states=Vector{Int16}()
    compartment=Vector{Int16}()
    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)

            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)
                append!(states,activestate)
                append!(compartment,buc )
                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

    end
    return states,compartment,energies_take ,energies_rand
end

function get_false_compartment_stats_compartment(l,n,x,γ,β,runin,compartments,steps,reps)

    states=Vector{Int16}()
    compartment=Vector{Int16}()
    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)

            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)
                append!(states,activestate)
                append!(compartment,buc )
                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

    end
    return sum(sign.(abs.((convert.(Int64,(states.+(compartments-1)) .% compartments .+1 )) .-compartment)))/length(compartment),mean(1 ./ (1 .+exp.(β .*(energies_take.-energies_rand))))
end

function get_false_compartment_stats_compartment_alongx(l,n,xs,γ,β,runin,compartments,steps,reps)
    d1=get_false_compartment_stats_compartment.(l,n,xs,γ,β,runin,compartments,steps,reps)
    e1s=zeros(length(xs))
e2s=zeros(length(xs))
for ii in 1 :(length(xs))
    e1s[ii] = d1[ii][1]
    e2s[ii] = d1[ii][2]
end
return e1s,e2s
end

function get_false_compartment_stats_compartment_alonggamma(l,n,x,γs,β,runin,compartments,steps,reps)
    d1=get_false_compartment_stats_compartment.(l,n,x,γs,β,runin,compartments,steps,reps)
    e1s=zeros(length(γs))
e2s=zeros(length(γs))
for ii in 1 :(length(γs))
    e1s[ii] = d1[ii][1]
    e2s[ii] = d1[ii][2]
end
return e1s,e2s
end

function get_statistics_withe_compartment(l,n,x,γ,β,runin,compartments,steps,reps,eqsteps,βeq)

    states=Vector{Int16}()
    compartment=Vector{Int16}()
    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)
            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)
                append!(states,activestate)
                append!(compartment,buc )
                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return states,compartment,energies_take ,energies_rand,overlaps
end

function get_statistics_compartment(l,n,x,γ,β,runin,compartments,steps,reps,eqsteps,βeq)

    states=Vector{Int16}()
    compartment=Vector{Int16}()
    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)
            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)
                append!(states,activestate)
                append!(compartment,buc )
                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return states,compartment,sum(sign.(energies_take .-maximum(energies_rand)).-1)./(length(energies_take)*(-2)),sum(sign.(energies_rand .-minimum(energies_take)).+1)/(2*length(energies_rand)),overlaps
end

function get_statistics_compartment_no_use_estats(l,n,x,γ,β,runin,compartments,steps,reps,eqsteps,βeq)

    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)
            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)

                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return sum(sign.(energies_take .-maximum(energies_rand)).-1)./(length(energies_take)*(-2)),sum(sign.(energies_rand .-minimum(energies_take)).+1)/(2*length(energies_rand)),mean(energies_take),std(energies_take),lower_std(energies_take),upper_std(energies_take),overlaps
end

function get_statistics_compartment_no_use(l,n,x,γ,β,runin,compartments,steps,reps,eqsteps,βeq)

    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_compartment(l,n,x,γ,β,runin,compartments)
            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J,buc,e_take=update_compartments_stats(J,xi[:,activestate],γ,β)

                append!(energies_take,e_take)
                append!(energies_rand,energie(J[:,:,buc],start_pattern(l,0.5)))
            end

            finalXI_vec=Vec_equiliabration_compartment(J,xi,β,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return sum(sign.(energies_take .-maximum(energies_rand)).-1)./(length(energies_take)*(-2)),sum(sign.(energies_rand .-minimum(energies_take)).+1)/(2*length(energies_rand)),overlaps
end

function get_statistics_One(l,n,x,γ,runin,steps,reps,eqsteps,βeq)

    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_matrix_xi_x_runin(l,n,x,γ,runin)
            finalXI_vec=Vec_equiliabration(J,xi,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                append!(energies_take,energie(J,xi[:,activestate]))
                append!(energies_rand,energie(J,start_pattern(l,0.5)))
                J=new_J(J,xi[:,activestate],γ)

            end

            finalXI_vec=Vec_equiliabration(J,xi,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return sum(sign.(energies_take .-maximum(energies_rand)).-1)./(length(energies_take)*(-2)),sum(sign.(energies_rand .-minimum(energies_take)).+1)/(2*length(energies_rand)),overlaps
end

function get_statistics_One_estats(l,n,x,γ,runin,steps,reps,eqsteps,βeq)

    energies_take=Vector{Float64}()
    energies_rand=Vector{Float64}()
    overlaps=Vector{Int16}()

    for run in 1:reps
            J,xi=build_matrix_xi_x_runin(l,n,x,γ,runin)
            finalXI_vec=Vec_equiliabration(J,xi,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                append!(energies_take,energie(J,xi[:,activestate]))
                append!(energies_rand,energie(J,start_pattern(l,0.5)))
                J=new_J(J,xi[:,activestate],γ)

            end

            finalXI_vec=Vec_equiliabration(J,xi,βeq,eqsteps)
            data=abs.(transpose(xi)*(finalXI_vec))
            for ii in 1:length(data[:,1])
                append!(overlaps, data[ii,ii])
            end
    end
    return sum(sign.(energies_take .-maximum(energies_rand)).-1)./(length(energies_take)*(-2)),sum(sign.(energies_rand .-minimum(energies_take)).+1)/(2*length(energies_rand)),mean(energies_take),std(energies_take),lower_std(energies_take),upper_std(energies_take),overlaps
end

function get_mean_overlaps_One(l,n,x,γ,runin,steps,reps,eqsteps,βeq)

     dumi1,dumi2,dumi3,dumi4,dumi5,dumi6,overlaps=get_statistics_One_estats(l,n,x,γ,runin,steps,reps,eqsteps,βeq)
     return mean(overlaps)/l
end

function get_statistics_One_allstats(l,n,x,γ,runin,steps,reps,eqsteps,βeq)   # produce all statistics for one compartment at a given point

    dist_start=Vector{Int16}()
    dist_next=Vector{Int16}()
    open_paths=Vector{Int16}()

    dist_start_rand=Vector{Int16}()
    dist_next_rand=Vector{Int16}()
    open_paths_rand=Vector{Int16}()

    partision_mean_vecs=Vector{Float64}()
    partision_std_vecs=Vector{Float64}()
    partision_first_vecs=Vector{Float64}()
    partision_last_vecs=Vector{Float64}()

    partision_mean_patts=Vector{Float64}()
    partision_std_patts=Vector{Float64}()
    partision_first_patts=Vector{Float64}()
    partision_last_patts=Vector{Float64}()

    best_overlaps=Vector{Array{Float64,1}}()
    mean_current_direction=Vector{Array{Float64,1}}()
    best_next_direciton=Vector{Array{Float64,1}}()

    for run in 1:reps
            J,xi=build_matrix_xi_x_runin(l,n,x,γ,runin)

            xi[:]=new_XI_rand(xi[:], x)
            finalXI_vec=Vec_equiliabration(J,copy(xi),βeq,eqsteps)

            F=eigen(J)
                 evecs=zeros(l,n)
                 evlues=zeros(n)
                 for jj in 0:n-1
                    evecs[:,jj+1] = F.vectors[:,l-jj]
                    evlues[jj+1]=F.values[l-jj]
                  end

            dist=abs.(transpose(finalXI_vec)*xi)
            for ii in 1:length(dist[:,1])
                append!(dist_start, convert(Int16,dist[ii,ii]))
                 dist[ii,ii]=0
                 append!(dist_next,convert(Int16,maximum(abs.(dist[ii,:]))))
                 append!(open_paths,convert(Int16,sum(-sign.(xi[ :,ii].*(transpose(J[:,:])*xi[:,ii])).+1)/2))
                 d1,d2,d3=produce_overlap_stats(xi[:,ii],evecs,evlues)
                 append!(best_overlaps,[d1])
                 append!(mean_current_direction,[d2])
                 append!(best_next_direciton,[d3])
            end
            energies=Vector{Float64}()
            for ii in 1:n
               append!(energies,energie(J,xi[:,ii]))
            end

            p1,p2,p3,p4,p5,p6,p7,p8=Support_partisipation(xi,evecs,energies)

            append!(partision_mean_vecs,p1)
            append!(partision_std_vecs,p2)
            append!(partision_first_vecs,p3)
            append!(partision_last_vecs,p4)

            append!(partision_mean_patts,p5)
            append!(partision_std_patts,p6)
            append!(partision_first_patts,p7)
            append!(partision_last_patts,p8)

            xi_rand= get_patterns(l,n,0.5)
            finalXI_rand=Vec_equiliabration(J,copy(xi_rand),βeq,eqsteps)

            dist_rand=abs.(transpose(finalXI_rand)*xi_rand)

            dist_rand_pats=abs.(transpose(finalXI_rand)*xi)
            for ii in 1:length(dist[:,1])
                append!(dist_start_rand, convert(Int16,dist_rand[ii,ii]))

                append!(dist_next_rand,convert(Int16,maximum(abs.(dist_rand_pats[ii,:]))))
                append!(open_paths_rand,convert(Int16,sum(-sign.(xi_rand[ :,ii].*(transpose(J[:,:])*xi_rand[:,ii])).+1)/2))

            end
            for ii in 1:steps

                activestate=updateactivestate(n)
                xi[:]=new_XI_rand(xi[:], x)
                J=new_J(J,xi[:,activestate],γ)

            end
            xi[:]=new_XI_rand(xi[:], x)
            finalXI_vec=Vec_equiliabration(J,xi,βeq,eqsteps)

            F=eigen(J)
                 evecs=zeros(l,n)
                 evlues=zeros(n)
                 for jj in 0:n-1
                    evecs[:,jj+1] = F.vectors[:,l-jj]
                    evlues[jj+1]=F.values[l-jj]
                  end

            dist=abs.(transpose(finalXI_vec)*xi)
            for ii in 1:length(dist[:,1])
                append!(dist_start, convert(Int16,dist[ii,ii]))
                 dist[ii,ii]=0
                 append!(dist_next,convert(Int16,maximum(abs.(dist[ii,:]))))
                 append!(open_paths,convert(Int16,sum(-sign.(xi[ :,ii].*(transpose(J[:,:])*xi[:,ii])).+1)/2))
                 d1,d2,d3=produce_overlap_stats(xi[:,ii],evecs,evlues)
                 append!(best_overlaps,[d1])
                 append!(mean_current_direction,[d2])
                 append!(best_next_direciton,[d3])
            end
            energies=Vector{Float64}()
           for ii in 1:n
               append!(energies,energie(J,xi[:,ii]))
           end
            p1,p2,p3,p4,p5,p6,p7,p8=Support_partisipation(xi,evecs,energies)

            append!(partision_mean_vecs,p1)
            append!(partision_std_vecs,p2)
            append!(partision_first_vecs,p3)
            append!(partision_last_vecs,p4)

            append!(partision_mean_patts,p5)
            append!(partision_std_patts,p6)
            append!(partision_first_patts,p7)
            append!(partision_last_patts,p8)

    end
    return dist_start,dist_next,open_paths,     dist_start_rand,dist_next_rand,open_paths_rand,   mean(partision_mean_vecs),sqrt(mean(partision_std_vecs.^2)),mean(partision_first_vecs),mean(partision_last_vecs),     mean(partision_mean_patts),sqrt(mean(partision_std_patts.^2)),mean(partision_first_patts),mean(partision_last_patts),best_overlaps,mean_current_direction,best_next_direciton

end

function produce_overlap_stats(xi,evecs,evlues)
    n=length(evlues)
    l=length(xi)
    α=n/l
    states=convert.(Int64,range(1,stop=n,length=n))
    best_overlaps=Vector{Float64}()
    best_next_direciton=Vector{Float64}()
    mean_current_direction=Vector{Float64}()
    o_squared_curent=0
    lam_current=0

    overlaps=(transpose(xi)*evecs)./(sqrt(l))

    for ii in 1:min(10,n-4)

        positionnow=states[argmax(abs.(overlaps[states]))]

        lam_current=copy(lam_current)+overlaps[positionnow]^2*evlues[positionnow]
        o_squared_curent=copy(o_squared_curent)+copy(overlaps[positionnow])^2

        append!(best_overlaps,o_squared_curent)
        append!(mean_current_direction,lam_current*(o_squared_curent))
        states=setdiff(states,positionnow)
        append!(best_next_direciton,maximum(abs.(transpose(evlues[states]).*evecs[:,states]))*sqrt(abs.(1-o_squared_curent))/sqrt(α)*(sqrt(n/(n-ii))))
    end
    return best_overlaps,mean_current_direction,best_next_direciton
end

function build_rec_for_gamma_One(l,n,x,runin,steps,βeq,equil_steps,γ_steps,γ_min,reps)

    pointsgamtake=convert.(Int16,range(1,stop=20,length=20))
    pointsgamtake2=convert.(Int16,range(1,stop=10,length=10))
        γ=Vector{Float64}()

        if(γ_min< 1 - 0.02)
            append!(γ,range(0.9*γ_min,stop=γ_min,length=21)[pointsgamtake])
            append!(γ,range(γ_min,stop=0.99999,length=convert(Int,γ_steps-20)))

        else

            append!(γ,range(0.9*γ_min,stop=1-0.02,length=21)[pointsgamtake])

            append!(γ,range(1-0.02,stop=min(0.999,1 - (1-γ_min)/4),length=convert(Int,γ_steps-30)))
            append!(γ,range(min(0.999,1 - (1-γ_min)/4),stop=0.99999,length=convert(Int,11))[pointsgamtake2])
        end

    e_rac=zeros(γ_steps)
    e_rand=zeros(γ_steps)
    runin_e=max.(10*n,min.(floor.(Int64,log(0.00001)./ log.(γ) .+1),runin))
    steps_e=min.(runin_e,steps)
    overlaps=zeros(Int16,n*reps*2,γ_steps)
    for ii in 1:γ_steps
        e_rac[ii],e_rand[ii],overlaps[:,ii]=get_statistics_One(l,n,x,γ[ii],runin_e[ii],steps_e[ii],reps,equil_steps,βeq)
    end
    return e_rac,e_rand,overlaps,γ
end

function build_rec_for_gamma_compartment(l,n,x,runin,steps,β,βeq,equil_steps,γ_steps,γ_min,reps,compartments)
    pointsgamtake=convert.(Int16,range(1,stop=20,length=20))
    pointsgamtake2=convert.(Int16,range(1,stop=10,length=5))
        γ=Vector{Float64}()

        if(γ_min< 1 - 0.02)
            append!(γ,range(0.9*γ_min,stop=γ_min,length=21)[pointsgamtake])
            append!(γ,range(γ_min,stop=0.99999,length=convert(Int,γ_steps-20)))

        else

            append!(γ,range(0.9*γ_min,stop=1-0.02,length=21)[pointsgamtake])

            append!(γ,range(1-0.02,stop=min(0.999,1 - (1-γ_min)/4),length=convert(Int,γ_steps-25)))
            append!(γ,range(min(0.999,1 - (1-γ_min)/4),stop=0.99999,length=convert(Int,6))[pointsgamtake2])
        end

    pattern_used=zeros(Int16,steps*reps,γ_steps)
    compartment_used=zeros(Int16,steps*reps,γ_steps)
    runin_e=max.(10*n,min.(floor.(Int64,2*compartments*(log(0.00001)./ log.(γ) .+1)),runin))
    e_rac=zeros(γ_steps)
    e_rand=zeros(γ_steps)
    overlaps=zeros(Int16,n*reps*2,γ_steps)
    for ii in 1:γ_steps
        pattern_used[:,ii],compartment_used[:,ii],e_rac[ii],e_rand[ii],overlaps[:,ii]=get_statistics_compartment(l,n,x,γ[ii],β,runin_e[ii],compartments,steps,reps,equil_steps,βeq)
    end
    return pattern_used,compartment_used,e_rac,e_rand,overlaps,γ
end

function build_rec_for_gamma_compartment_no_use(l,n,x,runin,steps,β,βeq,equil_steps,γ_steps,γ_min,reps,compartments)

    pointsgamtake=convert.(Int16,range(1,stop=20,length=20))
    pointsgamtake2=convert.(Int16,range(2,stop=6,length=5))
        γ=Vector{Float64}()

        if(γ_min< 1 - 0.02)
            append!(γ,range(0.001,stop=γ_min,length=21)[pointsgamtake])
            append!(γ,range(γ_min,stop=0.99999,length=convert(Int,γ_steps-20)))

        else

            append!(γ,range(0.01,stop=1-0.02,length=21)[pointsgamtake])

            append!(γ,range(1-0.02,stop=min(0.999,1 - (1-γ_min)/4),length=convert(Int,γ_steps-25)))
            append!(γ,range(min(0.999,1 - (1-γ_min)/4),stop=0.99999,length=convert(Int,6))[pointsgamtake2])
        end

    runin_e=max.(10*n,min.(floor.(Int64,2*compartments*(log(0.00001)./ log.(γ) .+1)),runin))
    steps_e=min.(runin_e,steps)
    e_rac=zeros(γ_steps)
    e_rand=zeros(γ_steps)
    overlaps=zeros(Int16,n*reps*2,γ_steps)
    for ii in 1:γ_steps
        e_rac[ii],e_rand[ii],overlaps[:,ii]=get_statistics_compartment_no_use(l,n,x,γ[ii],β,runin_e[ii],compartments,steps_e[ii],reps,equil_steps,βeq)
    end
    return e_rac,e_rand,overlaps,γ
end

function build_rec_for_gamma_x_One(l,n,runin,steps,βeq,equil_steps,γ_steps,x_max,xsteps,reps)
    x=Vector{Float64}()
    append!(x,0.0)
    append!(x,exp.(range(log(0.0001),stop=log(x_max),length=xsteps-1)))
    x=x./n

    γ_min=0.01

        γ_min= max.(1-1.999 /n,1 .-3 .*sqrt.(8 .*x./(n -1)))

    γ_min[1]=0.95  # reset for x=0

    e_rac=zeros(γ_steps,xsteps)
    e_rand=zeros(γ_steps,xsteps)
    overlaps=zeros(Int16,n*reps*2,γ_steps,xsteps)
    γ=zeros(γ_steps,xsteps)
    for ii in 1:xsteps
        e_rac[:,ii],e_rand[:,ii],overlaps[:,:,ii],γ[:,ii]=build_rec_for_gamma_One(l,n,x[ii],runin,steps,βeq,equil_steps,γ_steps,γ_min[ii],reps)
    end

    return e_rac,e_rand,overlaps,γ,x
end

function build_rec_for_gamma_x_compartment(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
    x=Vector{Float64}()
    append!(x,0.0)
    append!(x,exp.(range(log(0.0001),stop=log(x_max),length=xsteps-1)))   # set the smallest number of mutation rates
    x=x./n

    γ_min=0.01 .*ones(xsteps)

    if(compartments<n)
        γ_min= max.(1-1.999 * compartments/n,1 .-3 .*sqrt.(8 .*x./(n/compartments -1)))
    end

    γ_min[1]=0.95  # reset for x=0

    pattern_used=zeros(Int16,steps*reps,γ_steps,xsteps)
    compartment_used=zeros(Int16,steps*reps,γ_steps,xsteps)
    e_rac=zeros(γ_steps,xsteps)
    e_rand=zeros(γ_steps,xsteps)
    overlaps=zeros(Int16,n*reps*2,γ_steps,xsteps)
        γ=zeros(γ_steps,xsteps)
    for ii in 1:xsteps
        e_rac[:,ii],e_rand[:,ii],overlaps[:,:,ii],γ[:,ii]=build_rec_for_gamma_compartment_no_use(l,n,x[ii],runin,steps,β,βeq,equil_steps,γ_steps,γ_min[ii],reps,compartments)
    end

    return e_rac,e_rand,overlaps,γ,x
end

function build_rec_for_opt_gamma_One(l,n,x,γ,runin,steps,βeq,equil_steps,reps)  # give the results and statistics along an (optimal) learning rate
    γ_steps=length(γ)
    e_rac=zeros(γ_steps)
    e_rand=zeros(γ_steps)
    e_mean=zeros(γ_steps)
    e_std_low=zeros(γ_steps)
    e_std_high=zeros(γ_steps)
    overlaps=zeros(Int16,n*reps*2,γ_steps)
    runin_e=max.(10*n,min.(floor.(Int64,log(0.00001)./ log.(γ) .+1),runin))
    steps_e=min.(runin_e,steps)
    for ii in 1:γ_steps
        e_rac[ii],e_rand[ii],e_mean[ii],stdnormal,e_std_low[ii],e_std_high[ii],overlaps[:,ii]=get_statistics_One_estats(l,n,x[ii],γ[ii],runin_e[ii],steps_e[ii],reps,equil_steps,βeq)
    end
    return e_rac,e_rand,e_mean,e_std_low,e_std_high,overlaps
end

function build_rec_for_opt_gamma_compartment(l,n,x,γ,runin,steps,β,βeq,equil_steps,reps,compartments)  # give the results and statistics along an (optimal) learning rate  for multiple compartments
    γ_steps=length(γ)
    pattern_used=zeros(Int16,steps*reps,γ_steps)
    compartment_used=zeros(Int16,steps*reps,γ_steps)

    e_rac=zeros(γ_steps)
    e_rand=zeros(γ_steps)
    overlaps=zeros(Int16,n*reps*2,γ_steps)
    for ii in 1:γ_steps
        pattern_used[:,ii],compartment_used[:,ii],e_rac[ii],e_rand[ii],overlaps[:,ii]=get_statistics_compartment(l,n,x[ii],γ[ii],β,runin,compartments,steps,reps,equil_steps,βeq)
    end
    return pattern_used,compartment_used,e_rac,e_rand,overlaps
end

function distribution_only_allEV(l,n,x,γ,runin,reps)

    evalue_stats=Vector{Array{Float64,1}}()
        for run in 1:reps

                J,xi=build_matrix_xi_x_runin(l,n,x,γ,runin)

                 F=eigen(J)
                 evlues=zeros(l)
                 for jj in 0:l-1
                    evlues[jj+1]=F.values[l-jj]
                  end
                append!(evalue_stats,[evlues])
            end
    evalue_means=zeros(l)

    for ii in 1:l
        evslocal=Vector{Float64}()
        for run in 1:reps
            append!(evslocal,evalue_stats[run][ii])
        end
        evalue_means[ii]=mean(evslocal)
    end

    return evalue_stats,evalue_means
end

function distribution_EV(l,n,x,γ,runin,reps)

    evalue_stats=Vector{Array{Float64,1}}()
    entropies=Vector{Tuple{Array{Float64,1},Float64,Float64,Array{Float64,1},Float64,Float64,Array{Float64,1},Float64,Float64}}()
    for run in 1:reps

                J,xi=build_matrix_xi_x_runin(l,n,x,γ,runin)

                 F=eigen(J)
                 evecs=zeros(l,n)
                 evlues=zeros(n)
                 for jj in 0:n-1
                    evecs[:,jj+1] = F.vectors[:,l-jj]
                    evlues[jj+1]=F.values[l-jj]
                  end
                append!(evalue_stats,[evlues])

                energies=Vector{Float64}()
               for ii in 1:n
                   append!(energies,energie(J,xi[:,ii]))
               end
                local_entropy=Support_entropy(xi,evecs,energies)

                append!(entropies,[local_entropy])

    end

    return evalue_stats,entropies
end

function distribution_EV_for_x(l,n,xsteps,xmin,xmax,runin_max,reps)
    x=exp.(range(log(xmin),stop=log(0.01),length=4))
    if (xsteps>5)
        append!(x,range(0.015,stop=xmax,length=xsteps-4))
    end
    γ=1 .- sqrt.(8 .*x./(n.-1))
    runin=min.(runin_max,floor.(Int64,log(0.00001)./ log.(γ) .+1))
    data=distribution_EV.(l,n,0,γ,runin,reps)
    return data
end


text:
using Makie: GLNormalMesh
using GeometryTypes: HyperRectangle, HyperSphere, HomogenousMesh
using GLAbstraction: rotationmatrix_x, rotationmatrix_y, rotationmatrix_z,
      translationmatrix, Point3f0, Vec3f0
using Colors: Colorant, RGBA, red, green, blue

export mesh, mesharray

struct MeshColorScheme
  segment   ::Colorant
  node      ::Colorant
  plane     ::Colorant
  planenode ::Colorant
end
const defaultmeshcolorscheme = MeshColorScheme(
  RGBA(0.2f0, 0.2f0, 1.0f0, 0.5f0),
  RGBA(1.0f0, 0.0f0, 0.0f0, 0.5f0),
  RGBA(0.0f0, 1.0f0, 0.0f0, 0.2f0),
  RGBA(0.0f0, 1.0f0, 0.0f0, 0.5f0)
)

struct VisualizationParameters
  meshcolorscheme :: MeshColorScheme
  nodesize :: Float32
end

function nodesize(element::Element, context::Context)
  cumulative_mindim = Inf32
  for e in element
    canidate_mindim = mindim(e, context)
    if canidate_mindim < cumulative_mindim
      cumulative_mindim = canidate_mindim
    end
  end
  return cumulative_mindim/6.0f0
end
mindim(::Element, ::Context) = Inf32
function mindim(s::Segment, context::Context)
  (w,h) = width_height(s,context)
  if h>w
    return w
  else
    return h
  end
end

function rxyz(lengthvector::Array{Float64,1}, widthvector::Array{Float64,1})
  zangle = atan2(lengthvector[2],lengthvector[1])
  rot_rz = rotationmatrix_z(zangle)
  nrot_rz = rotationmatrix_z(-zangle)
  yangle = π/2+acos(lengthvector[3]/norm(lengthvector))
  rot_ry = rotationmatrix_y(yangle)
  nrot_ry = rotationmatrix_y(-yangle)
  rot_wxyz = nrot_ry * nrot_rz * vcat(widthvector,1.0)
  xangle = atan2(rot_wxyz[3],rot_wxyz[2])
  rot_rx = rotationmatrix_x(xangle)
  return rot_rz * rot_ry * rot_rx
end

function nodemesh(xyz::Array{Float64,1}, size::Float32, color::Colorant)
  GLNormalMesh((HyperSphere(Point3f0(xyz[1],xyz[2],xyz[3]), size), color))
end

appendmesh!(::Array{HomogenousMesh,1},::Element,::Context,
            ::VisualizationParameters) = nothing
function appendmesh!(allmesh::Array{HomogenousMesh,1},
                     node::Node,
                     context::Context,
                     vp::VisualizationParameters)
  push!(allmesh, nodemesh(xyz1(node,context)[1:3], vp.nodesize, vp.meshcolorscheme.node))
  return nothing
end

centeroftwopoints(a::Array{Float64,1}, b::Array{Float64,1}) =
  a[1:3] + (b[1:3]-a[1:3])./2

function appendmesh!(allmesh::Array{HomogenousMesh,1},
                     segment::Segment,
                     context::Context,
                     vp::VisualizationParameters)
   (n1_xyz1,n2_xyz1) = nodes_xyz1(segment, context)
   n1 = n1_xyz1[1:3]
   n2 = n2_xyz1[1:3]
   (width,height) = width_height(segment, context)
   lengthvector= (n2-n1)
   length = norm(lengthvector)
   mesh = GLNormalMesh((
            HyperRectangle(
              Vec3f0(-0.5f0*length,-0.5f0*width,-0.5f0*height),
              Vec3f0(length,width,height)
            ),
          vp.meshcolorscheme.segment))
   c = n1 + lengthvector./2
   correction = translationmatrix(Vec3f0(c...))*
                rxyz(lengthvector, wxyz(segment,context))
   push!(allmesh, correction*mesh)
   return nothing
end

function appendmesh!(allmesh::Array{HomogenousMesh,1},
                     plane::UniformPlane,
                     context::Context,
                     vp::VisualizationParameters)
  (c1,c2,c3,thick) = corners_xyz1_thick(plane, context)
  lxyz = c1[1:3]-c2[1:3]
  length = norm(lxyz)
  wxyz = c3[1:3]-c2[1:3]
  width = norm(wxyz)
  height = thick
  center = centeroftwopoints(c1,c3)
  uncorrectedplanemesh = GLNormalMesh((
    HyperRectangle(
      Vec3f0(-0.5f0*length,-0.5f0*width,-0.5f0*height),
      Vec3f0(length,width,height)),vp.meshcolorscheme.plane
  ))
  correction =
    translationmatrix(Vec3f0(center...)) *
    rxyz(lxyz,wxyz)
  push!(allmesh, correction * uncorrectedplanemesh)
  for xyz1 in nodes_xyz1(plane, context)
    push!(allmesh, nodemesh(xyz1[1:3],vp.nodesize, vp.meshcolorscheme.planenode))
  end
  return nothing
end

"""
    mesharray(element)

Returns an array of `HomogenousMesh` objects for use with `GLVisualize`.

`mesharray` will accept any `Element` as it's argument, but only makes sense
 for `Group`'s.  Planes will not show their holes.
"""
function mesharray(element::Element,mcs::MeshColorScheme=defaultmeshcolorscheme)
  context = Context(element)
  vp = VisualizationParameters(mcs, nodesize(element,context))
  allmesh = Array{HomogenousMesh}(0)
  for e in element
    appendmesh!(allmesh,e,context,vp)
  end
  return allmesh
end

"""
    mesh(element)

Equivalent to `merge(mesharray(element))`.

Use the following to visualize `element`.
```
m = mesh(element)

using GLVisualize
window = glscreen()
_view(visualize(m), window)
renderloop(window)
```
"""
mesh(element::Element) = merge(mesharray(element))


text:
"""
    @expr <expression>

Return the original expression object.

# Example

```julia
julia> ex = @expr x + 1
:(x + 1)
```
"""
macro expr(ex)
    return QuoteNode(ex)
end

"""
    @test_expr <type> <ex>

Test if the syntax type generates the same expression `ex`. Returns the
corresponding syntax type instance. Requires `using Test` before using
this macro.

# Example

```julia
def = @test_expr JLFunction function (x, y)
    return 2
end
@test is_kw_fn(def) == false
```
"""
macro test_expr(type, ex)
    @gensym def generated_expr original_expr
    quote
        $def = Expronicon.@expr $type $ex
        println($def)
        $generated_expr = $prettify($codegen_ast($def))
        $original_expr = $prettify($(Expr(:quote, ex)))
        @test $compare_expr($generated_expr, $original_expr)
        $def
    end |> esc
end

"""
    @test_expr <expr> == <expr>

Test if two expression is equivalent semantically, this uses `compare_expr`
to decide if they are equivalent, ignores things such as `LineNumberNode`
generated `Symbol` in `Expr(:curly, ...)` or `Expr(:where, ...)`.
"""
macro test_expr(ex::Expr)
    ex.head === :call && ex.args[1] === :(==) || error("expect <expr> == <expr>, got $ex")
    lhs, rhs = ex.args[2], ex.args[3]
    quote
        $__source__
        @test $compare_expr($prettify($lhs), $prettify($rhs))
    end |> esc
end

"""
    @expr <type> <expression>

Return the expression in given type.

# Example

```julia
julia> ex = @expr JLKwStruct struct Foo{N, T}
           x::T = 1
       end
#= kw =# struct Foo{N, T}
    #= /home/roger/code/julia/Expronicon/test/analysis.jl:5 =#
    x::T = 1
end
```
"""
macro expr(type, ex)
    quote
        $type($(Expr(:quote, ex)))
    end |> esc
end

struct AnalysisError <: Exception
    expect::String
    got
end

anlys_error(expect, got) = throw(AnalysisError(expect, got))

function Base.show(io::IO, e::AnalysisError)
    print(io, "expect ", e.expect, " expression, got ", e.got, ".")
end

"""
    compare_expr(lhs, rhs)

Compare two expression of type `Expr` or `Symbol` semantically, which:

1. ignore the detail value `LineNumberNode` in comparision
2. ignore the detailed name of typevars in `Expr(:curly, ...)` or `Expr(:where, ...)`

This gives a way to compare two Julia expression semantically which means
although some details of the expression is different but they should
produce the same lowered code.
"""
function compare_expr(lhs, rhs)
    @switch (lhs, rhs) begin
        @case (::Symbol, ::Symbol)
            lhs === rhs
        @case (Expr(:curly, name, lhs_vars...), Expr(:curly, &name, rhs_vars...))
            all(map(compare_vars, lhs_vars, rhs_vars))
        @case (Expr(:where, lbody, lparams...), Expr(:where, rbody, rparams...))
            compare_expr(lbody, rbody) &&
                all(map(compare_vars, lparams, rparams))
        @case (Expr(head, largs...), Expr(&head, rargs...))
                isempty(largs) && isempty(rargs) ||
            (length(largs) == length(rargs) && all(map(compare_expr, largs, rargs)))
        # ignore LineNumberNode
        @case (::LineNumberNode, ::LineNumberNode)
            true
        @case _
            lhs == rhs
    end
end

"""
    compare_vars(lhs, rhs)

Compare two expression by assuming all `Symbol`s are variables,
thus their value doesn't matter, only where they are matters under
this assumption. See also [`compare_expr`](@ref).
"""
function compare_vars(lhs, rhs)
    @switch (lhs, rhs) begin
        @case (::Symbol, ::Symbol)
            true
        @case (Expr(head, largs...), Expr(&head, rargs...))
            all(map(compare_vars, largs, rargs))
        # ignore LineNumberNode
        @case (::LineNumberNode, ::LineNumberNode)
            true
        @case _
            lhs == rhs
    end
end

"""
    is_literal(x)

Check if `x` is a literal value.
"""
function is_literal(x)
    !(x isa Expr || x isa Symbol || x isa GlobalRef)
end

"""
    is_gensym(s)

Check if `s` is generated by `gensym`.

!!! note
    Borrowed from [MacroTools](https://github.com/FluxML/MacroTools.jl).
"""
is_gensym(s::Symbol) = occursin("#", string(s))
is_gensym(s) = false

"""
    gensym_name(x::Symbol)

Return the gensym name.

!!! note
    Borrowed from [MacroTools](https://github.com/FluxML/MacroTools.jl).
"""
function gensym_name(x::Symbol)
    m = Base.match(r"##(.+)#\d+", String(x))
    m === nothing || return m.captures[1]
    m = Base.match(r"#\d+#(.+)", String(x))
    m === nothing || return m.captures[1]
    return "x"
end

"""
    support_default(f)

Check if field type `f` supports default value.
"""
support_default(f) = false
support_default(f::JLKwField) = true

function has_symbol(@nospecialize(ex), name::Symbol)
    ex isa Symbol && return ex === name
    ex isa Expr || return false
    return any(x->has_symbol(x, name), ex.args)
end

"""
    has_kwfn_constructor(def[, name = struct_name_plain(def)])

Check if the struct definition contains keyword function constructor of `name`.
The constructor name to check by default is the plain constructor which does
not infer any type variables and requires user to input all type variables.
See also [`struct_name_plain`](@ref).
"""
function has_kwfn_constructor(def, name = struct_name_plain(def))
    any(def.constructors) do fn::JLFunction
        isempty(fn.args) && fn.name == name
    end
end

"""
    has_plain_constructor(def, name = struct_name_plain(def))

Check if the struct definition contains the plain constructor of `name`.
By default the name is the inferable name [`struct_name_plain`](@ref).

# Example

```julia
def = @expr JLKwStruct struct Foo{T, N}
    x::Int
    y::N

    Foo{T, N}(x, y) where {T, N} = new{T, N}(x, y)
end

has_plain_constructor(def) # true

def = @expr JLKwStruct struct Foo{T, N}
    x::T
    y::N

    Foo(x, y) = new{typeof(x), typeof(y)}(x, y)
end

has_plain_constructor(def) # false
```

the arguments must have no type annotations.

```julia
def = @expr JLKwStruct struct Foo{T, N}
    x::T
    y::N

    Foo{T, N}(x::T, y::N) where {T, N} = new{T, N}(x, y)
end

has_plain_constructor(def) # false
```
"""
function has_plain_constructor(def, name = struct_name_plain(def))
    any(def.constructors) do fn::JLFunction
        fn.name == name || return false
        fn.kwargs === nothing || return false
        length(def.fields) == length(fn.args) || return false
        for (f, x) in zip(def.fields, fn.args)
            f.name === x || return false
        end
        return true
    end
end

"""
    is_function(def)

Check if given object is a function expression.
"""
function is_function(@nospecialize(def))
    @match def begin
        ::JLFunction => true
        Expr(:function, _, _) => true
        Expr(:(=), _, _) => true
        Expr(:(->), _, _) => true
        _ => false
    end
end

"""
    is_kw_function(def)

Check if a given function definition supports keyword arguments.
"""
function is_kw_function(@nospecialize(def))
    is_function(def) || return false

    if def isa JLFunction
        return def.kwargs !== nothing
    end

    _, call, _ = split_function(def)
    @match call begin
        Expr(:tuple, Expr(:parameters, _...), _...) => true
        Expr(:call, _, Expr(:parameters, _...), _...) => true
        Expr(:block, _, ::LineNumberNode, _) => true
        _ => false
    end
end

@deprecate is_kw_fn(def) is_kw_function(def)
@deprecate is_fn(def) is_function(def)

"""
    is_struct(ex)

Check if `ex` is a struct expression.
"""
function is_struct(@nospecialize(ex))
    ex isa Expr || return false
    return ex.head === :struct
end

"""
    is_struct_not_kw_struct(ex)

Check if `ex` is a struct expression excluding keyword struct syntax.
"""
function is_struct_not_kw_struct(ex)
    is_struct(ex) || return false
    body = ex.args[3]
    body isa Expr && body.head === :block || return false
    any(is_field_default, body.args) && return false
    return true
end

"""
    is_ifelse(ex)

Check if `ex` is an `if ... elseif ... else ... end` expression.
"""
function is_ifelse(@nospecialize(ex))
    ex isa Expr || return false
    return ex.head === :if
end

"""
    is_for(ex)

Check if `ex` is a `for` loop expression.
"""
function is_for(@nospecialize(ex))
    ex isa Expr || return false
    return ex.head === :for
end

"""
    is_field(ex)

Check if `ex` is a valid field expression.
"""
function is_field(@nospecialize(ex))
    @match ex begin
        :($name::$type = $default) => false
        :($(name::Symbol) = $default) => false
        name::Symbol => true
        :($name::$type) => true
        _ => false
    end
end

"""
    is_field_default(ex)

Check if `ex` is a `<field expr> = <default expr>` expression.
"""
function is_field_default(@nospecialize(ex))
    @match ex begin
        :($name::$type = $default) => true
        :($(name::Symbol) = $default) => true
        _ => false
    end
end

"""
    split_doc(ex::Expr) -> line, doc, expr

Split doc string from given expression.
"""
function split_doc(ex::Expr)
    if ex.head === :macrocall && ex.args[1] == GlobalRef(Core, Symbol("@doc"))
        return ex.args[2], ex.args[3], ex.args[4]
    else
        return nothing, nothing, ex
    end
end

"""
    split_function(ex::Expr) -> head, call, body

Split function head declaration with function body.
"""
function split_function(ex::Expr)
    @match ex begin
        Expr(:function, call, body) => (:function, call, body)
        Expr(:(=), call, body) => (:(=), call, body)
        Expr(:(->), call, body) => (:(->), call, body)
        _ => anlys_error("function", ex)
    end
end

"""
    split_function_head(ex::Expr) -> name, args, kw, whereparams, rettype

Split function head to name, arguments, keyword arguments and where parameters.
"""
function split_function_head(ex::Expr)
    @match ex begin
        Expr(:tuple, Expr(:parameters, kw...), args...) => (nothing, args, kw, nothing, nothing)
        Expr(:tuple, args...) => (nothing, args, nothing, nothing, nothing)
        Expr(:call, name, Expr(:parameters, kw...), args...) => (name, args, kw, nothing, nothing)
        Expr(:call, name, args...) => (name, args, nothing, nothing, nothing)
        Expr(:block, x, ::LineNumberNode, Expr(:(=), kw, value)) => (nothing, Any[x], Any[Expr(:kw, kw, value)], nothing, nothing)
        Expr(:block, x, ::LineNumberNode, kw) => (nothing, Any[x], Any[kw], nothing, nothing)
        Expr(:(::), call, rettype) => begin
            name, args, kw, whereparams, _ = split_function_head(call)
            (name, args, kw, whereparams, rettype)
        end
        Expr(:where, call, whereparams...) => begin
            name, args, kw, _, rettype = split_function_head(call)
            (name, args, kw, whereparams, rettype)
        end
        _ => anlys_error("function head expr", ex)
    end
end

"""
    split_struct_name(ex::Expr) -> name, typevars, supertype

Split the name, type parameters and supertype definition from `struct`
declaration head.
"""
function split_struct_name(@nospecialize(ex))
    return @match ex begin
        :($name{$(typevars...)}) => (name, typevars, nothing)
        :($name{$(typevars...)} <: $type) => (name, typevars, type)
        ::Symbol => (ex, [], nothing)
        :($name <: $type) => (name, [], type)
        _ => anlys_error("struct", ex)
    end
end

"""
    split_struct(ex::Expr) -> ismutable, name, typevars, supertype, body

Split struct definition head and body.
"""
function split_struct(ex::Expr)
    ex.head === :struct || error("expect a struct expr, got $ex")
    name, typevars, supertype = split_struct_name(ex.args[2])
    body = ex.args[3]
    return ex.args[1], name, typevars, supertype, body
end

function split_ifelse(ex::Expr)
    conds, stmts = [], []
    otherwise = split_ifelse!((conds, stmts), ex)
    return conds, stmts, otherwise
end

function split_ifelse!((conds, stmts), ex::Expr)
    ex.head in [:if, :elseif] || return ex
    push!(conds, ex.args[1])
    push!(stmts, ex.args[2])

    if length(ex.args) == 3
        return split_ifelse!((conds, stmts), ex.args[3])
    end
    return
end

function split_forloop(ex::Expr)
    ex.head === :for || error("expect a for loop expr, got $ex")
    lhead = ex.args[1]
    lbody = ex.args[2]
    return split_for_head(lhead)..., lbody
end

function split_for_head(ex::Expr)
    if ex.head === :block
        vars, itrs = [], []
        for each in ex.args
            each isa Expr || continue # skip other things
            var, itr = split_single_for_head(each)
            push!(vars, var)
            push!(itrs, itr)
        end
        return vars, itrs
    else
        var, itr = split_single_for_head(ex)
        return Any[var], Any[itr]
    end
end

function split_single_for_head(ex::Expr)
    ex.head === :(=) || error("expect a single loop head, got $ex")
    return ex.args[1], ex.args[2]
end

function uninferrable_typevars(def::Union{JLStruct, JLKwStruct}; leading_inferable::Bool=true)
    typevars = name_only.(def.typevars)
    field_types = [field.type for field in def.fields]

    if leading_inferable
        idx = findfirst(typevars) do t
            !any(map(f->has_symbol(f, t), field_types))
        end
        idx === nothing && return []
    else
        idx = 0
    end
    uninferrable = typevars[1:idx]

    for T in typevars[idx+1:end]
        any(map(f->has_symbol(f, T), field_types)) || push!(uninferrable, T)
    end
    return uninferrable
end

"""
    JLFunction(ex::Expr)

Create a `JLFunction` object from a Julia function `Expr`.

# Example

```julia
julia> JLFunction(:(f(x) = 2))
f(x) = begin
    #= REPL[37]:1 =#
    2
end
```
"""
function JLFunction(ex::Expr)
    line, doc, expr = split_doc(ex)
    head, call, body = split_function(expr)
    name, args, kw, whereparams, rettype = split_function_head(call)
    JLFunction(head, name, args, kw, rettype, whereparams, body, line, doc)
end

"""
    JLStruct(ex::Expr)

Create a `JLStruct` object from a Julia struct `Expr`.

# Example

```julia
julia> JLStruct(:(struct Foo
           x::Int
       end))
struct Foo
    #= REPL[38]:2 =#
    x::Int
end
```
"""
function JLStruct(ex::Expr)
    line, doc, expr = split_doc(ex)
    ismutable, typename, typevars, supertype, body = split_struct(expr)

    fields, constructors, misc = JLField[], JLFunction[], []
    field_doc, field_line = nothing, nothing

    body = flatten_blocks(body)

    for each in body.args
        @switch each begin
            @case :($name::$type)
                push!(fields, JLField(name, type, field_doc, field_line))
            @case name::Symbol
                push!(fields, JLField(name, Any, field_doc, field_line))
            @case ::String
                field_doc = each
            @case ::LineNumberNode
                field_line = each
            @case GuardBy(is_function)
                if name_only(each) === typename
                    push!(constructors, JLFunction(each))
                else
                    push!(misc, each)
                end
            @case _
                push!(misc, each)
        end
    end
    JLStruct(typename, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLKwStruct(ex::Expr, typealias=nothing)

Create a `JLKwStruct` from given Julia struct `Expr`, with an option to attach
an alias to this type name.

# Example

```julia
julia> JLKwStruct(:(struct Foo
           x::Int = 1
       end))
#= kw =# struct Foo
    #= REPL[39]:2 =#
    x::Int = 1
end
```
"""
function JLKwStruct(ex::Expr, typealias=nothing)
    line, doc, expr = split_doc(ex)
    ismutable, typename, typevars, supertype, body = split_struct(expr)

    fields, constructors, misc = JLKwField[], JLFunction[], []
    field_doc, field_line = nothing, nothing
    body = flatten_blocks(body)
    for each in body.args
        @switch each begin
            @case :($name::$type = $default)
                push!(fields, JLKwField(name, type, field_doc, field_line, default))
            @case :($(name::Symbol) = $default)
                push!(fields, JLKwField(name, Any, field_doc, field_line, default))
            @case name::Symbol
                push!(fields, JLKwField(name, Any, field_doc, field_line, no_default))
            @case :($name::$type)
                push!(fields, JLKwField(name, type, field_doc, field_line, no_default))
            @case ::String
                field_doc = each
            @case ::LineNumberNode
                field_line = each
            @case GuardBy(is_function)
                if name_only(each) === typename
                    push!(constructors, JLFunction(each))
                else
                    push!(misc, each)
                end
            @case _
                push!(misc, each)
        end
    end
    JLKwStruct(typename, typealias, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLIfElse(ex::Expr)

Create a `JLIfElse` from given Julia ifelse `Expr`.

# Example

```julia
julia> ex = :(if foo(x)
             x = 1 + 1
         elseif goo(x)
             y = 1 + 2
         else
             error("abc")
         end)
:(if foo(x)
      #= REPL[41]:2 =#
      x = 1 + 1
  elseif #= REPL[41]:3 =# goo(x)
      #= REPL[41]:4 =#
      y = 1 + 2
  else
      #= REPL[41]:6 =#
      error("abc")
  end)

julia> JLIfElse(ex)
if foo(x)
    begin
        #= REPL[41]:2 =#
        x = 1 + 1
    end
elseif begin
    #= REPL[41]:3 =#
    goo(x)
end
    begin
        #= REPL[41]:4 =#
        y = 1 + 2
    end
else
    begin
        #= REPL[41]:6 =#
        error("abc")
    end
end
```
"""
function JLIfElse(ex::Expr)
    ex.head === :if || error("expect an if ... elseif ... else ... end expression")
    conds, stmts, otherwise = split_ifelse(ex)
    return JLIfElse(conds, stmts, otherwise)
end

"""
    JLFor(ex::Expr)

Create a `JLFor` from given Julia for loop expression.

# Example

```julia
julia> ex = @expr for i in 1:10, j in 1:j
           M[i, j] += 1
       end
:(for i = 1:10, j = 1:j
      #= REPL[3]:2 =#
      M[i, j] += 1
  end)

julia> jl = JLFor(ex)
for i in 1 : 10,
    j in 1 : j
    #= loop body =#
    begin
        #= REPL[3]:2 =#
        M[i, j] += 1
    end
end

julia> jl.vars
2-element Vector{Any}:
 :i
 :j

julia> jl.iterators
2-element Vector{Any}:
 :(1:10)
 :(1:j)
```
"""
function JLFor(ex::Expr)
    vars, itrs, body = split_forloop(ex)
    return JLFor(vars, itrs, body)
end


text:
"""
$(SIGNATURES)

Calls clip
"""
function clip_in_tet_coordinates(p::poly_eight{4,T}) where {T}
    if length(p) == 3
        return clip_in_tet_coordinates(p.v[1], p.v[2], p.v[3])
    elseif length(p) == 4
        return clip_in_tet_coordinates(p.v[1], p.v[2], p.v[3], p.v[4])
    else
        error("something is wrong")
    end
end

function clip_in_tet_coordinates(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}) where {T}
    return clip(z1, z2, z3, 1)
end

function clip_in_tet_coordinates(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}) where {T}
    return clip(z1, z2, z3, z4, 1)
end

"""
$(SIGNATURES)

Together with cut_clip, implements the clipping of polygons by tetrahedra using the Sutherland–Hodgman algorithm.
The Sutherland-Hodgman algorithm essentially clips a polygon by one plane at a time.
You should not call this function directly.
Clipped polygons are assumed to have at most 8 sides.
This implementation should correctly handle the case where 2 of fewer vertices lie exactly on a plane.
"""
function clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, i::Int64) where {T}
    # Polygons are represented as vertices in tetrahedral coordinates.
    # If all components of a vertex are positive this means that the vertex lies within the tetrahedron
    # Negative components mean that a vertix lies outside a face.
    # This function permutes the vertices until the first is non-positive and the second is positive for cut_clip.

    (i == 5) && (return poly_eight(3, (z1, z2, z3, z1, z1, z1, z1, z1)))  # there is no 5th plane return polygon

    s = SVector{3,T}(z1[i], z2[i], z3[i])
    is_non_pos = (s .<= 0.0)
    all(is_non_pos) && (return poly_eight{4,T}())  # entire polygon clipped away
    if all(0.0 .<= s)  # non_negative
        return clip(z1, z2, z3, i+1)  # all points inside this plane go to next plane
    else
        is_non_pos[1] && !is_non_pos[2] && (return cut_clip(z1, z2, z3, i))
        is_non_pos[2] && !is_non_pos[3] && (return cut_clip(z2, z3, z1, i))
        is_non_pos[3] && !is_non_pos[1] && (return cut_clip(z3, z1, z2, i))
    end
    error("Non-finite vertex likely")
end

function clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, i::Int64) where {T}
    (i == 5) && (return poly_eight(4, (z1, z2, z3, z4, z1, z1, z1, z1)))

    s = SVector{4,T}(z1[i], z2[i], z3[i], z4[i])
    is_non_pos = (s .<= 0.0)
    all(is_non_pos) && (return poly_eight{4,T}())
    if all(0.0 .<= s)  # non_negative
        return clip(z1, z2, z3, z4, i + 1)
    else
        is_non_pos[1] && !is_non_pos[2] && (return cut_clip(z1, z2, z3, z4, i))
        is_non_pos[2] && !is_non_pos[3] && (return cut_clip(z2, z3, z4, z1, i))
        is_non_pos[3] && !is_non_pos[4] && (return cut_clip(z3, z4, z1, z2, i))
        is_non_pos[4] && !is_non_pos[1] && (return cut_clip(z4, z1, z2, z3, i))
    end
    error("Non-finite vertex likely")
end

function clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, i::Int64) where {T}
    (i == 5) && (return poly_eight(5, (z1, z2, z3, z4, z5, z1, z1, z1)))

    s = SVector{5,T}(z1[i], z2[i], z3[i], z4[i], z5[i])
    is_non_pos = (s .<= 0.0)
    all(is_non_pos) && (return poly_eight{4,T}())
    if all(0.0 .<= s)  # non_negative
        return clip(z1, z2, z3, z4, z5, i + 1)
    else
        is_non_pos[1] && !is_non_pos[2] && (return cut_clip(z1, z2, z3, z4, z5, i))
        is_non_pos[2] && !is_non_pos[3] && (return cut_clip(z2, z3, z4, z5, z1, i))
        is_non_pos[3] && !is_non_pos[4] && (return cut_clip(z3, z4, z5, z1, z2, i))
        is_non_pos[4] && !is_non_pos[5] && (return cut_clip(z4, z5, z1, z2, z3, i))
        is_non_pos[5] && !is_non_pos[1] && (return cut_clip(z5, z1, z2, z3, z4, i))
    end
    error("Non-finite vertex likely")
end

function clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, z6::SVector{4,T}, i::Int64) where {T}
    (i == 5) && (return poly_eight(6, (z1, z2, z3, z4, z5, z6, z1, z1)))

    s = SVector{6,T}(z1[i], z2[i], z3[i], z4[i], z5[i], z6[i])
    is_non_pos = (s .<= 0.0)
    all(is_non_pos) && (return poly_eight{4,T}())
    if all(0.0 .<= s)  # non_negative
        return clip(z1, z2, z3, z4, z5, z6, i + 1)
    else
        is_non_pos[1] && !is_non_pos[2] && (return cut_clip(z1, z2, z3, z4, z5, z6, i))
        is_non_pos[2] && !is_non_pos[3] && (return cut_clip(z2, z3, z4, z5, z6, z1, i))
        is_non_pos[3] && !is_non_pos[4] && (return cut_clip(z3, z4, z5, z6, z1, z2, i))
        is_non_pos[4] && !is_non_pos[5] && (return cut_clip(z4, z5, z6, z1, z2, z3, i))
        is_non_pos[5] && !is_non_pos[6] && (return cut_clip(z5, z6, z1, z2, z3, z4, i))
        is_non_pos[6] && !is_non_pos[1] && (return cut_clip(z6, z1, z2, z3, z4, z5, i))
    end
    error("Non-finite vertex likely")
end

function clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, z6::SVector{4,T}, z7::SVector{4,T}, i::Int64) where {T}
    (i == 5) && (return poly_eight(7, (z1, z2, z3, z4, z5, z6, z7, z1)))

    s = SVector{7,T}(z1[i], z2[i], z3[i], z4[i], z5[i], z6[i], z7[i])
    is_non_pos = (s .<= 0.0)
    all(is_non_pos) && (return poly_eight{4,T}())
    if all(0.0 .<= s)  # non_negative
        return clip(z1, z2, z3, z4, z5, z6, z7, i + 1)
    else
        is_non_pos[1] && !is_non_pos[2] && (return cut_clip(z1, z2, z3, z4, z5, z6, z7, i))
        is_non_pos[2] && !is_non_pos[3] && (return cut_clip(z2, z3, z4, z5, z6, z7, z1, i))
        is_non_pos[3] && !is_non_pos[4] && (return cut_clip(z3, z4, z5, z6, z7, z1, z2, i))
        is_non_pos[4] && !is_non_pos[5] && (return cut_clip(z4, z5, z6, z7, z1, z2, z3, i))
        is_non_pos[5] && !is_non_pos[6] && (return cut_clip(z5, z6, z7, z1, z2, z3, z4, i))
        is_non_pos[6] && !is_non_pos[7] && (return cut_clip(z6, z7, z1, z2, z3, z4, z5, i))
        is_non_pos[7] && !is_non_pos[1] && (return cut_clip(z7, z1, z2, z3, z4, z5, z6, i))
    end
    error("Non-finite vertex likely")
end

"""
$(SIGNATURES)

Removes a non-positive vertex from a polygon.
"""
function cut_clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, i::Int64) where {T}
    # z1 is gauranteed to be non-positive
    # z2 is gauranteed to be positive

    z_start = clip_node(z1, z2, i)
    if 0.0 < z3[i]
        z_end = clip_node(z1, z3, i)
        return clip(z_start, z2, z3, z_end, i + 1)
    else
        z_end = clip_node(z3, z2, i)
        return clip(z_start, z2, z_end, i + 1)
    end
end

function cut_clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, i::Int64) where {T}
    (z3[i] <= 0.0) && (return cut_clip(z1, z2, z3, i))
    z_start = clip_node(z1, z2, i)
    if 0.0 < z4[i]
        z_end = clip_node(z1, z4, i)
        return clip(z_start, z2, z3, z4, z_end, i + 1)
    else
        z_end = clip_node(z4, z3, i)
        return clip(z_start, z2, z3, z_end, i + 1)
    end
end

function cut_clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, i::Int64) where {T}
    (z4[i] <= 0.0) && (return cut_clip(z1, z2, z3, z4, i))
    z_start = clip_node(z1, z2, i)
    if 0.0 < z5[i]
        z_end = clip_node(z1, z5, i)
        return clip(z_start, z2, z3, z4, z5, z_end, i + 1)
    else
        z_end = clip_node(z5, z4, i)
        return clip(z_start, z2, z3, z4, z_end, i + 1)
    end
end

function cut_clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, z6::SVector{4,T}, i::Int64) where {T}
    (z5[i] <= 0.0) && (return cut_clip(z1, z2, z3, z4, z5, i))
    z_start = clip_node(z1, z2, i)
    if 0.0 <= z6[i]
        z_end = clip_node(z1, z6, i)
        return clip(z_start, z2, z3, z4, z5, z6, z_end, i + 1)
    else
        z_end = clip_node(z6, z5, i)
        return clip(z_start, z2, z3, z4, z5, z_end, i + 1)
    end
end

function cut_clip(z1::SVector{4,T}, z2::SVector{4,T}, z3::SVector{4,T}, z4::SVector{4,T}, z5::SVector{4,T}, z6::SVector{4,T}, z7::SVector{4,T}, i::Int64) where {T}
    (z6[i] <= 0.0) && (return cut_clip(z1, z2, z3, z4, z5, z6, i))
    z_start = clip_node(z1, z2, i)
    if 0.0 <= z7[i]
        z_end = clip_node(z1, z7, i)
        return poly_eight(8, (z_start, z2, z3, z4, z5, z6, z7, z_end))
    else
        z_end = clip_node(z7, z6, i)
        return poly_eight(7, (z_start, z2, z3, z4, z5, z6, z_end, z_start))
    end
end

function clip_node(z_non::SVector{4,T}, z_pos::SVector{4,T}, i::Int64) where {T}
    v_non = z_non[i]
    v_pos = z_pos[i]
    return weightPoly(z_non, z_pos, v_non, v_pos)
end


text:
#1-dim CAs

#Elementary & totalistic CA with r-nearest neighbor and k states
struct CellularAutomaton

    #user given values
    N::Int #code
    ruleset::Array{Int8,1} #set of rules
    k::Int #number of states
    r::Int #r-nearest neighbor

    #internal variables
    cells::Array{Int8,2}

    function CellularAutomaton(N::Int, init::Array{Int,1}, gen::Int; k::Int=2, r::Int=1)
        ruleset = rule(N, k, r)

        w = length(init)
        cells = zeros(Int8, gen, w)
        cells[1,:] = Array{Int8}(init[:])'

        if k == 2
            #Elementary CA
            mp = reverse(Int8[k^i for i = 0:2r+1])
        else
            #Totalistic CA
            mp = Array{Int8}(ones(2r+1))
        end

        for i = 2:gen
            for j = 1:w
                ind = 0 ; s = 1
                for q = (j-r):(j+r)

                    #Cyclic boundary conditions
                    if q < 1
                        q = w-q
                    end
                    if q > w
                        q = q-w
                    end

                    ind += mp[s]*cells[i-1, q]
                    s += 1
                end
                cells[i,j] = ruleset[ind+1]
            end
        end

        new(N, ruleset, k, r, cells)
    end

    #Quick setup for one seed cell in the center
    function CellularAutomaton(N::Int, gen::Int; kvs...)
        init = Array{Int}(zeros(2*gen))
        init[gen] = 1

        CellularAutomaton(N, init, gen; kvs...)
    end
end

#Parse rule numbering according to the Wolfram code
function rule(n::Int, k=2, r=1)
    if k < 2
        error("number of states (k) must be larger than 1; you gave $k")
    elseif k == 2
        #Elementary CA
        return digits(n, base=k, pad=k^(2r+1))
    else
        #Totalistic CA
        return digits(n, base=k, pad=(2r+1)k-2)
    end
end


text:
using NLsolve
using Dolang
using DataFrames

function newton(f, init; tol=1e-8, maxit=500, backsteps=5, verbose=false)

    line = 0.5.^(0:(backsteps-1))
    x0 = init[:]
    x = copy(x0)
    N = length(x0)
    it = 0
    err = 1.0
    while (it<maxit) & (err>tol)
        it += 1
        ress, jacc = f(reshape(x0,size(init)...))
        res = reshape(ress,N)
        jac = reshape(jacc,N,N)
        err = maximum(abs(ress))
        ul = 0.1
        if err>tol
            dx = -jac\res
            nerr = 2.0
            for l in line
                ul = l
                try
                    x = x0+l*dx
                    nr, junk = f(reshape(x,size(init)...))
                    nerr = maximum(abs(nr))
                catch
                    println("skip")
                end
                if nerr<err
                    break
                end
            end
            if nerr>err
                println("Backtracking didn't work")
            end
        end
        x0 = x
        if verbose
            println(it, ": ", err, ":", ul)
        end
    end
    return reshape(x0,size(init)...)
end

function deterministic(model, exogenous_series; verbose=false, T=100)

    N=T

    symbols = model.symbols
    endogenous = symbols[:endogenous]
    exogenous = symbols[:exogenous]
    parameters = symbols[:parameters]

    it = Dolang.IncidenceTable(model.equations)

    fut_syms = it.by_date[1]
    cur_syms = it.by_date[0]
    past_syms = it.by_date[-1]

    fut_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[1]]
    cur_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[0]]
    past_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[-1]]
    shocks_inds = [i for (i,v) in enumerate(exogenous) if v in it.by_date[0]]
    parms_inds = [i for (i,v) in enumerate(parameters) if v in it.by_date[0]]

    P = length(symbols[:endogenous])
    indices = fut_inds
    indices = cat(1,indices, P+cur_inds)

    indices = cat(1,indices, 2*P+past_inds)
    findices = cat(1,indices, 3*P+shocks_inds)

    ne = length(model.symbols[:exogenous])

    T = size(exogenous_series,1)
    exogenous = zeros(N, ne)
    exogenous[1:T,:] = exogenous_series
    for t=T+1:N
        exogenous[t,:] = exogenous[T,:]
    end

    f_s = model.functions.f_static
    f_d = model.functions.f_dynamic

    y0 = [model.calibration[e] for e in model.symbols[:endogenous]]
    e0 = [model.calibration[e] for e in model.symbols[:exogenous]]
    p = [model.calibration[e] for e in model.symbols[:parameters]]

    e_start = exogenous[1,:]
    e_final = exogenous[end,:]

    function fobj_start(x, out)
        out[:] = f_s(cat(1, x, e_start) , p)
    end

    function fobj_final(x, out)
        out[:] = f_s(cat(1, x, e_final) , p)
    end

    fdiff_start = NLsolve.DifferentiableMultivariateFunction(fobj_start)
    sol_start = nlsolve(fdiff_start,y0)
    y_start = sol_start.zero
    fdiff_final = NLsolve.DifferentiableMultivariateFunction(fobj_final)
    sol_final = nlsolve(fdiff_final,y0)
    y_final = sol_final.zero

    y_init_0 = [(1-l)*y_start + l*y_final for l in linspace(0,1,N)]
    y_init = cat(1, [e' for e in y_init_0]...) :: Array{Float64, 2}

    n_y = length(model.symbols[:endogenous])
    n_e = length(model.symbols[:exogenous])

    function residual(y)

        a = y[1:end-2,:]   # y{t-1}
        b = y[2:end-1,:]   # y_t
        c = y[3:end,:]     # y_{t+1}
        d = exogenous[1:end-2,:]

        x = cat(2,c,b,a,d)
        xx = x[:,findices]

        N,P = size(y)

        res = zeros(N,P)
        jac0 = zeros(N,P,3*P)
        for n = 2:N-1
            res[n,:] = f_d(Der{0},xx[n-1,:],p)
            jac0[n,:,indices] = f_d(Der{1},xx[n-1,:],p)[:,1:length(indices)]
        end
        jac = zeros(N,P,N,P)
        for n in 2:N-1
            jac[n,:,n,:] += jac0[n,:,(P+1):2*P]  # y_t
            jac[n,:,n+1,:] += jac0[n,:,1:P] # y_{t+1}
            jac[n,:,n-1,:] +=  jac0[n,:,2*P+1:3*P]  # y_{t-1}
        end
        jac[1,:,1,:] = eye(P)
        jac[end,:,end,:] = eye(P)
        res[1,:] = y[1,:]-y_start
        res[end,:] = y[end,:]-y_final
        return res, jac
    end
    a,b = residual(y_init)
    if verbose
        println("Initial: ", maximum(abs(a)))
    end
    sol = newton(residual, y_init, verbose=verbose)

    tvec = collect(0:(size(sol,1)-1))
    sol = cat(2,tvec,sol,exogenous)

    colnames = cat(1, [:t], model.symbols[:endogenous], model.symbols[:exogenous])
    columns = Dict(colnames[i] => sol[:,i] for i=1:size(sol,2))
    return DataFrame(columns)
end


text:
using ElasticPDMats, PDMats, LinearAlgebra
using Test

@testset "PDMat Tests" begin
    a = rand(10, 10); m = a*a';
    epdmat = ElasticPDMat(m, capacity = 2000);
    test_pdmat(epdmat, m, verbose = 0)
end

@testset "append!" begin
    a = rand(10, 10); m = a*a';
    epdmat = ElasticPDMat(m[1:9, 1:9])
    append!(epdmat, m[:, 10])
    @test cholesky(m).U ≈ view(epdmat.chol).U

    epdmat = ElasticPDMat(m[1:6, 1:6])
    append!(epdmat, m[:, 7:10])
    @test cholesky(m).U ≈ view(epdmat.chol).U

    epdmat = ElasticPDMat(m[1:6, 1:6], capacity = 6, stepsize = 6)
    append!(epdmat, m[:, 7:10])
    @test epdmat.mat.m.capacity[1] == epdmat.chol.capacity ==
        size(epdmat.mat.m.data, 1) == size(epdmat.chol.c.factors, 1) == 12
    @test cholesky(m).U ≈ view(epdmat.chol).U
end

@testset "deleteat!" begin
    a = rand(10, 10); m = a*a';
    epdmat = ElasticPDMat(m)
    deleteat!(epdmat, 3)
    m2 = m[[1:2; 4:10], [1:2; 4:10]]
    @test cholesky(m2).U ≈ view(epdmat.chol).U

    a = rand(10, 10); m = a*a';
    epdmat = ElasticPDMat(m)
    deleteat!(epdmat, [3, 8, 7])
    m2 = m[[1:2; 4:6; 9:10], [1:2; 4:6; 9:10]]
    @test cholesky(m2).U ≈ view(epdmat.chol).U
end

@testset "util" begin
    a = rand(10, 10); m = a*a';
    epdmat = ElasticPDMat(m, capacity = 100, stepsize = 50)
    epdmat2 = ElasticPDMat(m, cholesky(m), capacity = 100, stepsize = 50)
    epdmat3 = ElasticPDMat()
    @test epdmat3.mat.m.capacity[1] == epdmat3.chol.capacity == epdmat3.mat.m.stepsize[1] == epdmat3.chol.stepsize == 10^3
    @test epdmat.mat.m.data == epdmat2.mat.m.data
    @test epdmat.chol.c.factors == epdmat2.chol.c.factors
    @test size(epdmat.mat) == (10, 10)
    setstepsize!(epdmat, 100)
    @test epdmat.mat.m.stepsize[1] == epdmat.chol.stepsize == 100
    setcapacity!(epdmat, 20)
    @test epdmat.mat.m.capacity[1] == epdmat.chol.capacity == 20
    @test size(epdmat.mat.m.data) == size(epdmat.chol.c.factors) == (20, 20)
    @test view(epdmat.mat)[:] == m[:]
    @test view(epdmat.chol).U == cholesky(m).U
end


text:
using NumericalMethodsforEngineers, Test

sv = [16., 4., 5., 8., -4., 22.]
kdiag = [1, 3, 6]

amat = fromskyline(kdiag, sv)
amat |> display
println()

println("\n'Julia' version:\n")
a = amat
b = [4., 2., 5.]

a[2, 2] += 1.0e20
b[2] = a[2, 2] * 5.0

a = sparse(a)

f = cholesky(a)
y = f.L \ b
c = f.U \ y

println("Solution Vector: \n $c")

@test c == a \ b

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

println("\n'Fortran' version:\n")
b = [4., 2., 5.]

nfix = 1
no = [2;]
val = [5.0;]
penalty = 1.0e20

sv[kdiag[no]] .+= penalty
b[no] = sv[kdiag[no]] .* val

sparin!(sv, kdiag)
spabac!(sv, b, kdiag)
println("Solution Vector: \n $b")

sv = [16., 4., 5., 8., -4., 22.]
b = [4., 2., 5.]
sv[kdiag[no]] .+= penalty
b[no] = sv[kdiag[no]] .* val

svt = sparin(sv, kdiag)
bt = spabac(svt, b, kdiag)
println("Solution Vector 2: \n $bt")


text:
struct Molecules{T<:AbstractSequencingView}
    genome::Vector{LongSequence{DNAAlphabet{2}}}
    views::Vector{T}
end

"""
    Molecules(gen::Vector{LongSequence{DNAAlphabet{2}}}, ng::Int = 1)

Create a pool of `ng` copies of a genome defined by the `gen` vector of sequences.
"""
function Molecules(gen::Vector{LongSequence{DNAAlphabet{2}}}, ng::Int = 1)
    vs = Vector{BasicSequencingView}()
    for (i, seq) in enumerate(gen)
        push!(vs, BasicSequencingView(i, 1, length(seq)))
    end
    if ng > 1
        vs = amplify(vs, ng)
    end
    return Molecules{BasicSequencingView}(gen, vs)
end

"""
    Molecules(rdr::FASTA.Reader, ng::Int = 1)

Create a pool of `ng` copies of the genome read in from the `FASTA.Reader`.
"""
function Molecules(rdr::FASTA.Reader, ng::Int = 1)
    gen = Vector{LongSequence{DNAAlphabet{2}}}()
    rec = eltype(rdr)()
    while !eof(rdr)
        read!(rdr, rec)
        push!(gen, FASTA.sequence(LongSequence{DNAAlphabet{2}}, rec))
    end
    return Molecules(gen, ng)
end

"""
    Molecules(file::String, ng::Int)

Create a pool of `ng` copies of the genome in the fasta formatted `file`.

!!! note
    The argument `iscircular` is currently not used.
"""
function Molecules(file::String, ng::Int = 1)
    open(FASTA.Reader, file) do rdr
        Molecules(rdr, ng)
    end
end

@inline views(u::Molecules) = u.views
@inline genome(u::Molecules) = u.genome
@inline nmolecules(m::Molecules) = length(views(m))
@inline Base.firstindex(m::Molecules) = 1
@inline Base.lastindex(m::Molecules) = nmolecules(m)
@inline Base.eachindex(m::Molecules) = Base.OneTo(lastindex(m))
@inline Base.eltype(m::Molecules{T}) where {T<:AbstractSequencingView} = Molecule{T}
@inline Base.length(m::Molecules) = nmolecules(m)

@inline function Base.getindex(m::Molecules{T}, i::Int) where {T<:AbstractSequencingView}
    @boundscheck begin
        if i ∉ eachindex(m)
            throw(BoundsError())
        end
    end
    return Molecule{T}(m, convert(Int, i))
end

@inline function Base.iterate(m::Molecules, state = firstindex(m))
    state ∈ eachindex(m) ? (@inbounds m[state], state + 1) : nothing
end

function Base.show(io::IO, p::Molecules{T}) where {T<:AbstractSequencingView}
    println(io, "Pool of $(nmolecules(p)) DNA molecules:")
    mx, av, mn = summarize_lengths(views(p))
    if mx === mn
        println(io, " All molecules are of the same size: $mx")
    else
        println(io, " Maximum molecule size: $mx")
        println(io, " Average molecule size: $av")
        println(io, " Minimum molecule size: $mn")
    end
    if T === TaggedSequencingView
        tags = summarize_tags(views(p))
        if length(tags) > 0
            println(io, " Number of distinct tags: $(length(keys(tags)))")
        end
    end
end

struct Molecule{T<:AbstractSequencingView}
    parent::Molecules{T}
    idx::Int
end

@inline genome(m::Molecule) = genome(m.parent)
@inline Base.length(m::Molecule) = @inbounds length(views(m.parent)[m.idx])

@inline function view(m::Molecule)
    v = views(m.parent)
    @inbounds return v[m.idx] # Assume Molecules can't be made with bad idx values.
end

function flip!(m::Molecule)
    v = views(m.parent)
    @inbounds begin
        v[m.idx] = reverse(v[m.idx])
    end
    return m
end

function amplify!(m::Molecule, n::Int)
    M = m.parent
    for _ in 1:n
        push!(M, m)
    end
    return m
end

function sequence(m::Molecule)
    return extract_sequence(genome(m.parent), view(m))
end

function Base.show(io::IO, m::Molecule)
    v = view(m)
    println(io, "A ", length(v), "bp DNA molecule")
    show(io, sequence(m))
    println(io, "\nOrigin: Haplotype ", seqid(v), ", First bp: ", first(v), ", Last bp: ", last(v))
end

# Transformations
# ---------------

function Base.push!(M::Molecules{T}, m::Molecule{T}) where {T<:AbstractSequencingView}
    if genome(M) ≢ genome(m)
        throw(ArgumentError("Can't add molecule: non-equivalent origin genomes"))
    end
    push!(views(M), view(m))
    return M
end

function Base.deleteat!(M::Molecules, i::Int...)
    deleteat!(views(M), i)
    return M
end

"""
    flip(p::Molecules)

Create a copy of flipped DNA molecules from some input pool.

The flip function lets you randomly flip some of the molecules in a pool to the
opposite orientation they are in.
"""
function flip(pred::F, M::Molecules) where {F<:Function}
    M′ = typeof(M)(genome(M), copy(views(M)))
    return flip!(pred, M′)
end

function flip!(pred::Function, M::Molecules)
    @inbounds for i ∈ eachindex(M)
        m = M[i]
        if pred(m)
            flip!(m)
        end
    end
end

function amplify!(pred::Function, M::Molecules, n::Int)
    for i ∈ eachindex(M)
        @inbounds m = M[i]
        if pred(m)
            amplify!(m, n)
        end
    end
    return M
end

function amplify(pred::F, M::Molecules, n::Int) where {F<:Function}
    M′ = typeof(M)(genome(M), copy(views(M)))
    return amplify!(pred, M′, n)
end

"""
    fragment(p::Molecules, meansize::Int)

Create a new pool of DNA molecules, by breaking up the DNA molecules in some input pool.

This method breaks up a DNA molecule in a pool `p`, such that the average
length of the fragments is approximately `meansize`.

It fragments a molecule by scattering an appropriate number of breakpoints
across the molecule, before cutting the molecule at those breakpoints.

!!! note
    Breakpoints are scattered entirely at random across a molecule.
    No two or more breakpoints can fall in exactly the same place, as those
    positions are sampled without replacement.

!!! note
    The appropriate number of breakpoints to scatter across a molecule is
    calculated as:

    ```math
    \\frac{L}{S} - 1
    ```

    Where ``L`` is the length of the molecule being fragmented, and ``S`` is
    the desired expected fragment size.
    This calculation assumes breakpoints fall randomly across the molecule
    (see above note).

!!! note
    If a DNA molecule being fragmented is smaller than the desired `meansize`,
    then it will not be broken, it will simply be included in the new pool.
"""
function fragment(p::Molecules, meansize::Int)
    np = typeof(p)(genome(p), fragment(views(p), meansize::Int))
    return np
end

"""
    subsample(p::Molecules, n::Int)

Create a new pool of DNA molecules by sampling an input pool.

!!! note
    DNA molecules in the input pool `p` are selected according to the
    uniform distribution; no one molecule is more or less likely to be selected
    than another.

!!! note
    Sampling is done without replacement, so it is impossible for the new
    pool that is created to recieve one molecule the input pool twice.
"""
function subsample(p::Molecules, n::Int)
    np = typeof(p)(genome(p), subsample(views(p), n))
    return np
end

function subsample(p::Molecules, cov::ExpectedCoverage, readlen)
    n = needed_sample_size(cov, sum(length.(genome(p)))bp, readlen)
    return subsample(p, n)
end

"""
    tag(u::Molecules{BasicSequencingView}, ntags::Int)

Create a pool of tagged DNA molecules from some input pool.

The new pool has the same DNA molecules as the input pool.
However, each DNA molecule in the new tagged pool will be assigned a tag
in the range of `1:ntags`.

For any tagged molecules in a pool, any other molecules that
are derived from that tagged molecule will inherit the same tag.
For example, if a DNA fragment in a pool is tagged, and then it is
subsequently fragmented during a [`fragment`](@ref) transform, then all the
smaller fragments derived from that long fragment will inherit that long
fragment's tag.

!!! note
    Which fragment gets a certain tag is completely random.
    It is possible for two distinct DNA molecules in a pool to be assigned
    the same tag. The likelihood of that happening depends on the size of the
    tag pool (`ntags`), and the number of fragments in the pool.
"""
function tag(p::Molecules{BasicSequencingView}, ntags::Int)
    @info "Attaching $(ntags) tags randomly to each molecule in pool..."
    np = Molecules{TaggedSequencingView}(genome(p), tag(views(p), ntags))
    @info "Done"
    return np
end

function select(fn::Function, p::Molecules)
    np = typeof(p)(genome(p), select(fn, views(p)))
end