text:
# Autogenerated wrapper script for PCL_jll for i686-linux-musl-cxx11
export libpcl_common, libpcl_features, libpcl_filters, libpcl_io, libpcl_io_ply, libpcl_kdtree, libpcl_keypoints, libpcl_ml, libpcl_octree, libpcl_recognition, libpcl_registration, libpcl_sample_consensus, libpcl_search, libpcl_segmentation, libpcl_stereo, libpcl_surface, libpcl_tracking, pcl_add_gaussian_noise, pcl_boundary_estimation, pcl_cluster_extraction, pcl_compute_cloud_error, pcl_compute_hausdorff, pcl_concatenate_points_pcd, pcl_convert_pcd_ascii_binary, pcl_crf_segmentation, pcl_demean_cloud, pcl_elch, pcl_extract_feature, pcl_fast_bilateral_filter, pcl_fpfh_estimation, pcl_generate, pcl_gp3_surface, pcl_grid_min, pcl_hdl_grabber, pcl_icp, pcl_icp2d, pcl_linemod_detection, pcl_local_max, pcl_lum, pcl_marching_cubes_reconstruction, pcl_match_linemod_template, pcl_mls_smoothing, pcl_morph, pcl_ndt2d, pcl_ndt3d, pcl_normal_estimation, pcl_outlier_removal, pcl_passthrough_filter, pcl_pcd2ply, pcl_pcd2vtk, pcl_pcd_change_viewpoint, pcl_pcd_convert_NaN_nan, pcl_pcd_introduce_nan, pcl_pclzf2pcd, pcl_plane_projection, pcl_ply2obj, pcl_ply2pcd, pcl_ply2ply, pcl_ply2raw, pcl_plyheader, pcl_poisson_reconstruction, pcl_progressive_morphological_filter, pcl_radius_filter, pcl_sac_segmentation_plane, pcl_spin_estimation, pcl_train_linemod_template, pcl_train_unary_classifier, pcl_transform_from_viewpoint, pcl_transform_point_cloud, pcl_unary_classifier_segment, pcl_uniform_sampling, pcl_vfh_estimation, pcl_voxel_grid, pcl_xyz2pcd

using CompilerSupportLibraries_jll
using FLANN_jll
using boost_jll
using Eigen_jll
JLLWrappers.@generate_wrapper_header("PCL")
JLLWrappers.@declare_library_product(libpcl_common, "libpcl_common.so.1.12")
JLLWrappers.@declare_library_product(libpcl_features, "libpcl_features.so.1.12")
JLLWrappers.@declare_library_product(libpcl_filters, "libpcl_filters.so.1.12")
JLLWrappers.@declare_library_product(libpcl_io, "libpcl_io.so.1.12")
JLLWrappers.@declare_library_product(libpcl_io_ply, "libpcl_io_ply.so.1.12")
JLLWrappers.@declare_library_product(libpcl_kdtree, "libpcl_kdtree.so.1.12")
JLLWrappers.@declare_library_product(libpcl_keypoints, "libpcl_keypoints.so.1.12")
JLLWrappers.@declare_library_product(libpcl_ml, "libpcl_ml.so.1.12")
JLLWrappers.@declare_library_product(libpcl_octree, "libpcl_octree.so.1.12")
JLLWrappers.@declare_library_product(libpcl_recognition, "libpcl_recognition.so.1.12")
JLLWrappers.@declare_library_product(libpcl_registration, "libpcl_registration.so.1.12")
JLLWrappers.@declare_library_product(libpcl_sample_consensus, "libpcl_sample_consensus.so.1.12")
JLLWrappers.@declare_library_product(libpcl_search, "libpcl_search.so.1.12")
JLLWrappers.@declare_library_product(libpcl_segmentation, "libpcl_segmentation.so.1.12")
JLLWrappers.@declare_library_product(libpcl_stereo, "libpcl_stereo.so.1.12")
JLLWrappers.@declare_library_product(libpcl_surface, "libpcl_surface.so.1.12")
JLLWrappers.@declare_library_product(libpcl_tracking, "libpcl_tracking.so.1.12")
JLLWrappers.@declare_executable_product(pcl_add_gaussian_noise)
JLLWrappers.@declare_executable_product(pcl_boundary_estimation)
JLLWrappers.@declare_executable_product(pcl_cluster_extraction)
JLLWrappers.@declare_executable_product(pcl_compute_cloud_error)
JLLWrappers.@declare_executable_product(pcl_compute_hausdorff)
JLLWrappers.@declare_executable_product(pcl_concatenate_points_pcd)
JLLWrappers.@declare_executable_product(pcl_convert_pcd_ascii_binary)
JLLWrappers.@declare_executable_product(pcl_crf_segmentation)
JLLWrappers.@declare_executable_product(pcl_demean_cloud)
JLLWrappers.@declare_executable_product(pcl_elch)
JLLWrappers.@declare_executable_product(pcl_extract_feature)
JLLWrappers.@declare_executable_product(pcl_fast_bilateral_filter)
JLLWrappers.@declare_executable_product(pcl_fpfh_estimation)
JLLWrappers.@declare_executable_product(pcl_generate)
JLLWrappers.@declare_executable_product(pcl_gp3_surface)
JLLWrappers.@declare_executable_product(pcl_grid_min)
JLLWrappers.@declare_executable_product(pcl_hdl_grabber)
JLLWrappers.@declare_executable_product(pcl_icp)
JLLWrappers.@declare_executable_product(pcl_icp2d)
JLLWrappers.@declare_executable_product(pcl_linemod_detection)
JLLWrappers.@declare_executable_product(pcl_local_max)
JLLWrappers.@declare_executable_product(pcl_lum)
JLLWrappers.@declare_executable_product(pcl_marching_cubes_reconstruction)
JLLWrappers.@declare_executable_product(pcl_match_linemod_template)
JLLWrappers.@declare_executable_product(pcl_mls_smoothing)
JLLWrappers.@declare_executable_product(pcl_morph)
JLLWrappers.@declare_executable_product(pcl_ndt2d)
JLLWrappers.@declare_executable_product(pcl_ndt3d)
JLLWrappers.@declare_executable_product(pcl_normal_estimation)
JLLWrappers.@declare_executable_product(pcl_outlier_removal)
JLLWrappers.@declare_executable_product(pcl_passthrough_filter)
JLLWrappers.@declare_executable_product(pcl_pcd2ply)
JLLWrappers.@declare_executable_product(pcl_pcd2vtk)
JLLWrappers.@declare_executable_product(pcl_pcd_change_viewpoint)
JLLWrappers.@declare_executable_product(pcl_pcd_convert_NaN_nan)
JLLWrappers.@declare_executable_product(pcl_pcd_introduce_nan)
JLLWrappers.@declare_executable_product(pcl_pclzf2pcd)
JLLWrappers.@declare_executable_product(pcl_plane_projection)
JLLWrappers.@declare_executable_product(pcl_ply2obj)
JLLWrappers.@declare_executable_product(pcl_ply2pcd)
JLLWrappers.@declare_executable_product(pcl_ply2ply)
JLLWrappers.@declare_executable_product(pcl_ply2raw)
JLLWrappers.@declare_executable_product(pcl_plyheader)
JLLWrappers.@declare_executable_product(pcl_poisson_reconstruction)
JLLWrappers.@declare_executable_product(pcl_progressive_morphological_filter)
JLLWrappers.@declare_executable_product(pcl_radius_filter)
JLLWrappers.@declare_executable_product(pcl_sac_segmentation_plane)
JLLWrappers.@declare_executable_product(pcl_spin_estimation)
JLLWrappers.@declare_executable_product(pcl_train_linemod_template)
JLLWrappers.@declare_executable_product(pcl_train_unary_classifier)
JLLWrappers.@declare_executable_product(pcl_transform_from_viewpoint)
JLLWrappers.@declare_executable_product(pcl_transform_point_cloud)
JLLWrappers.@declare_executable_product(pcl_unary_classifier_segment)
JLLWrappers.@declare_executable_product(pcl_uniform_sampling)
JLLWrappers.@declare_executable_product(pcl_vfh_estimation)
JLLWrappers.@declare_executable_product(pcl_voxel_grid)
JLLWrappers.@declare_executable_product(pcl_xyz2pcd)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FLANN_jll, boost_jll, Eigen_jll)
    JLLWrappers.@init_library_product(
        libpcl_common,
        "lib/libpcl_common.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_features,
        "lib/libpcl_features.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_filters,
        "lib/libpcl_filters.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_io,
        "lib/libpcl_io.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_io_ply,
        "lib/libpcl_io_ply.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_kdtree,
        "lib/libpcl_kdtree.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_keypoints,
        "lib/libpcl_keypoints.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_ml,
        "lib/libpcl_ml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_octree,
        "lib/libpcl_octree.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_recognition,
        "lib/libpcl_recognition.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_registration,
        "lib/libpcl_registration.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_sample_consensus,
        "lib/libpcl_sample_consensus.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_search,
        "lib/libpcl_search.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_segmentation,
        "lib/libpcl_segmentation.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_stereo,
        "lib/libpcl_stereo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_surface,
        "lib/libpcl_surface.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcl_tracking,
        "lib/libpcl_tracking.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        pcl_add_gaussian_noise,
        "bin/pcl_add_gaussian_noise",
    )

    JLLWrappers.@init_executable_product(
        pcl_boundary_estimation,
        "bin/pcl_boundary_estimation",
    )

    JLLWrappers.@init_executable_product(
        pcl_cluster_extraction,
        "bin/pcl_cluster_extraction",
    )

    JLLWrappers.@init_executable_product(
        pcl_compute_cloud_error,
        "bin/pcl_compute_cloud_error",
    )

    JLLWrappers.@init_executable_product(
        pcl_compute_hausdorff,
        "bin/pcl_compute_hausdorff",
    )

    JLLWrappers.@init_executable_product(
        pcl_concatenate_points_pcd,
        "bin/pcl_concatenate_points_pcd",
    )

    JLLWrappers.@init_executable_product(
        pcl_convert_pcd_ascii_binary,
        "bin/pcl_convert_pcd_ascii_binary",
    )

    JLLWrappers.@init_executable_product(
        pcl_crf_segmentation,
        "bin/pcl_crf_segmentation",
    )

    JLLWrappers.@init_executable_product(
        pcl_demean_cloud,
        "bin/pcl_demean_cloud",
    )

    JLLWrappers.@init_executable_product(
        pcl_elch,
        "bin/pcl_elch",
    )

    JLLWrappers.@init_executable_product(
        pcl_extract_feature,
        "bin/pcl_extract_feature",
    )

    JLLWrappers.@init_executable_product(
        pcl_fast_bilateral_filter,
        "bin/pcl_fast_bilateral_filter",
    )

    JLLWrappers.@init_executable_product(
        pcl_fpfh_estimation,
        "bin/pcl_fpfh_estimation",
    )

    JLLWrappers.@init_executable_product(
        pcl_generate,
        "bin/pcl_generate",
    )

    JLLWrappers.@init_executable_product(
        pcl_gp3_surface,
        "bin/pcl_gp3_surface",
    )

    JLLWrappers.@init_executable_product(
        pcl_grid_min,
        "bin/pcl_grid_min",
    )

    JLLWrappers.@init_executable_product(
        pcl_hdl_grabber,
        "bin/pcl_hdl_grabber",
    )

    JLLWrappers.@init_executable_product(
        pcl_icp,
        "bin/pcl_icp",
    )

    JLLWrappers.@init_executable_product(
        pcl_icp2d,
        "bin/pcl_icp2d",
    )

    JLLWrappers.@init_executable_product(
        pcl_linemod_detection,
        "bin/pcl_linemod_detection",
    )

    JLLWrappers.@init_executable_product(
        pcl_local_max,
        "bin/pcl_local_max",
    )

    JLLWrappers.@init_executable_product(
        pcl_lum,
        "bin/pcl_lum",
    )

    JLLWrappers.@init_executable_product(
        pcl_marching_cubes_reconstruction,
        "bin/pcl_marching_cubes_reconstruction",
    )

    JLLWrappers.@init_executable_product(
        pcl_match_linemod_template,
        "bin/pcl_match_linemod_template",
    )

    JLLWrappers.@init_executable_product(
        pcl_mls_smoothing,
        "bin/pcl_mls_smoothing",
    )

    JLLWrappers.@init_executable_product(
        pcl_morph,
        "bin/pcl_morph",
    )

    JLLWrappers.@init_executable_product(
        pcl_ndt2d,
        "bin/pcl_ndt2d",
    )

    JLLWrappers.@init_executable_product(
        pcl_ndt3d,
        "bin/pcl_ndt3d",
    )

    JLLWrappers.@init_executable_product(
        pcl_normal_estimation,
        "bin/pcl_normal_estimation",
    )

    JLLWrappers.@init_executable_product(
        pcl_outlier_removal,
        "bin/pcl_outlier_removal",
    )

    JLLWrappers.@init_executable_product(
        pcl_passthrough_filter,
        "bin/pcl_passthrough_filter",
    )

    JLLWrappers.@init_executable_product(
        pcl_pcd2ply,
        "bin/pcl_pcd2ply",
    )

    JLLWrappers.@init_executable_product(
        pcl_pcd2vtk,
        "bin/pcl_pcd2vtk",
    )

    JLLWrappers.@init_executable_product(
        pcl_pcd_change_viewpoint,
        "bin/pcl_pcd_change_viewpoint",
    )

    JLLWrappers.@init_executable_product(
        pcl_pcd_convert_NaN_nan,
        "bin/pcl_pcd_convert_NaN_nan",
    )

    JLLWrappers.@init_executable_product(
        pcl_pcd_introduce_nan,
        "bin/pcl_pcd_introduce_nan",
    )

    JLLWrappers.@init_executable_product(
        pcl_pclzf2pcd,
        "bin/pcl_pclzf2pcd",
    )

    JLLWrappers.@init_executable_product(
        pcl_plane_projection,
        "bin/pcl_plane_projection",
    )

    JLLWrappers.@init_executable_product(
        pcl_ply2obj,
        "bin/pcl_ply2obj",
    )

    JLLWrappers.@init_executable_product(
        pcl_ply2pcd,
        "bin/pcl_ply2pcd",
    )

    JLLWrappers.@init_executable_product(
        pcl_ply2ply,
        "bin/pcl_ply2ply",
    )

    JLLWrappers.@init_executable_product(
        pcl_ply2raw,
        "bin/pcl_ply2raw",
    )

    JLLWrappers.@init_executable_product(
        pcl_plyheader,
        "bin/pcl_plyheader",
    )

    JLLWrappers.@init_executable_product(
        pcl_poisson_reconstruction,
        "bin/pcl_poisson_reconstruction",
    )

    JLLWrappers.@init_executable_product(
        pcl_progressive_morphological_filter,
        "bin/pcl_progressive_morphological_filter",
    )

    JLLWrappers.@init_executable_product(
        pcl_radius_filter,
        "bin/pcl_radius_filter",
    )

    JLLWrappers.@init_executable_product(
        pcl_sac_segmentation_plane,
        "bin/pcl_sac_segmentation_plane",
    )

    JLLWrappers.@init_executable_product(
        pcl_spin_estimation,
        "bin/pcl_spin_estimation",
    )

    JLLWrappers.@init_executable_product(
        pcl_train_linemod_template,
        "bin/pcl_train_linemod_template",
    )

    JLLWrappers.@init_executable_product(
        pcl_train_unary_classifier,
        "bin/pcl_train_unary_classifier",
    )

    JLLWrappers.@init_executable_product(
        pcl_transform_from_viewpoint,
        "bin/pcl_transform_from_viewpoint",
    )

    JLLWrappers.@init_executable_product(
        pcl_transform_point_cloud,
        "bin/pcl_transform_point_cloud",
    )

    JLLWrappers.@init_executable_product(
        pcl_unary_classifier_segment,
        "bin/pcl_unary_classifier_segment",
    )

    JLLWrappers.@init_executable_product(
        pcl_uniform_sampling,
        "bin/pcl_uniform_sampling",
    )

    JLLWrappers.@init_executable_product(
        pcl_vfh_estimation,
        "bin/pcl_vfh_estimation",
    )

    JLLWrappers.@init_executable_product(
        pcl_voxel_grid,
        "bin/pcl_voxel_grid",
    )

    JLLWrappers.@init_executable_product(
        pcl_xyz2pcd,
        "bin/pcl_xyz2pcd",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("src/init.jl")# Include all relevant functions

"""
This file contains the script to execute learning containing
the Value and Critic function in seperate Networks
    => Extra for the NIPS challenge
    https://www.crowdai.org/challenges/nips-2018-ai-for-prosthetics-challenge
"""

θa = init_θ_single_norm(params.structure_net_pl)
θv = init_θ_single_norm(params.structure_net_vf)
best_net = init_θ_single_norm(params.structure_net_pl)
res = zeros(params.n_iter)
prec = zeros(params.n_iter)
all_traj  = init_all_traj(params)

#Global Variable of Environment!
@everywhere @pyimport osim.env as osim
@everywhere env = osim.ProstheticsEnv(visualize=true)
#@everywhere env = osim.L2RunEnv(visualize=true)

sleep(1)

println("-----------------------")
println("         start         ")
println("-----------------------")

function actor_critic_single(params::Params, θa::Any, θv::Any, res::Array{Float64,1}, prec::Array{Float64,1}, all_traj::Traj_Data, best_net::Any)
    #Init Part:
    best_net = copy(θa)
    min_cost = Inf
    #------------------------
    #Initialize the Policy
    ĝ_a = copy(θa)
    opt_para_a = optimizers(θa,   Adam; lr = params.lr_pl)

    #Initialize the Value Function
    ĝ_v = copy(θv)
    opt_para_v = optimizers(θv,   Adam; lr = params.lr_vf)

    #Init Progress bar
    progress = Progress(params.n_iter,1)

    # Init rms
    rms_all = init_rms(params.s_init)
    rms_loc = init_rms(params.s_init)

    #saving old data to see differences
    a_t_old = all_traj.traj[1].a_t
    weights_old = zeros(length(net_to_vec(θa)),2)

    for i_iter = 1:params.n_iter
        ĝ_a *= 0
        ĝ_v *= 0

        #execute the parallel_episodes
        (ĝ_a, ĝ_v, res[i_iter], prec[i_iter]) = parallel_episodes_skeleton(params,θa,θv,all_traj,rms_all)

        #update RMS:
        #=
        if i_iter == 1
            two_pass_variance!(all_traj, rms_all, params)
        else
            two_pass_variance!(all_traj, rms_loc, params)
            update_variance!(rms_all,rms_loc)
        end
        =#

        #Reorganice Data in minibatches
        mini_batch = build_minibatches(all_traj,params)

        #First refit the value function for new data batch:
        θv = refit_VF(params, θv, mini_batch, opt_para_v)

        #Now apply the gae for the new data based on updated value function
        gae_all!(all_traj,θv,params)

        #Reorganice Data in minibatches as advantages have changed...
        mini_batch = build_minibatches(all_traj,params)

        #Do the Policy / and Value updates
        (dump, θv) = parallel_PPO(params, θa, θv, mini_batch, opt_para_a, opt_para_v)

        #ĝ_a = parallel_gradient(params, θa, mini_batch)
        #θa  = do_TRPO_update(ĝ_a, mini_batch, θa, params)

        if min_cost > res[i_iter]
            min_cost = res[i_iter]
            best_net = copy(θa)
            println("new_best: ",min_cost)
        end

        #Visualization_stuff
        #-----------------------------------------------------------------------
        ProgressMeter.update!(progress,i_iter) #update progress bar
        print(" ",res[i_iter])
        a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
        weights_old, i_iter,θa, θv, res, prec/params.n_traj, ĝ_a, ĝ_v)
        #-----------------------------------------------------------------------

        @everywhere GC.gc(true)
    end
end

actor_critic_single(params, θa, θv, res, prec, all_traj, best_net)
save_data(res, prec, θa, θv, all_traj)

ext: jl

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# implementation of factorised samplers

using DataStructures
using Statistics
using SparseArrays
using LinearAlgebra

"""
    neighbours(G::Vector{<:Pair}, i) = G[i].second

Return extended neighbourhood of `i` including `i`.
`G`: graphs of neightbourhoods
"""
neighbours(G::Vector{<:Pair}, i) = G[i].second
#need refreshments
hasrefresh(::FactBoomerang) = true
hasrefresh(Z::ZigZag) = Z.λref > 0

"""
    λ(∇ϕ, i, x, θ, Z::ZigZag)
`i`th Poisson rate of the `ZigZag` sampler
"""
function λ(∇ϕi, i, x, θ, Z::ZigZag)
    pos(∇ϕi'*θ[i])
end

"""
    λ(∇ϕ, i, x, θ, Z::FactBoomerang)
`i`th Poisson rate of the `FactBoomerang` sampler
"""
function λ(∇ϕi, i, x, θ, B::FactBoomerang)
    pos((∇ϕi - (x[i] - B.μ[i])*B.Γ[i,i])*θ[i])
end

loosen(c, x) = c + x #+ log(c+1)*abs(x)/100
"""
    ab(G, i, x, θ, c, Flow)

Returns the constant term `a` and linear term `b` when computing the Poisson times
from the upper upper bounding rates λᵢ(t) = max(a + b*t)^2. The factors `a` and `b`
can be function of the current position `x`, velocity `θ`, tuning parameter `c` and
the Graph `G`
"""
function ab(G, i, x, θ, c, Z::ZigZag, args...)
    a = loosen(c[i], (idot(Z.Γ, i, x)  - idot(Z.Γ, i, Z.μ))'*θ[i])
    b = loosen(c[i]/100, θ[i]'*idot(Z.Γ, i, θ))
    a, b
end

function ab(G, i, x, θ, c, Z::FactBoomerang)
    nhd = neighbours(G, i)
    z = sqrt(sum((x[j] - Z.μ[j])^2 + θ[j]^2 for j in nhd))
    z2 = (x[i]^2 + θ[i]^2)
    a = c[i]*sqrt(z2)*z + z2*Z.Γ[i,i]
    b = 0.0
    a, b
end

function adapt!(c, i, factor)
    c[i] *= factor
    c
end

"""
    λ_bar(G, i, x, θ, c, Z)

Computes the bounding rate `λ_bar` at position `x` and velocity `θ`.
"""
#λ_bar(G, i, x, θ, c, Z::ZigZag) = pos(ab(G, i, x, θ, c, Z)[1])
#λ_bar(G, i, x, θ, c, Z::FactBoomerang) = pos(ab(G, i, x, θ, c, Z)[1])

function event(i, t, x, θ, Z::Union{ZigZag,FactBoomerang,JointFlow})
    t, i, x[i], θ[i]
end


text:
__precompile__(false)
module SymbolicControlSystems
using LinearAlgebra
using ControlSystems, SymPy, Latexify
import Symbolics
using InteractiveUtils

export sp,
    sym2num,
    latextf,
    Sym,
    @vars,
    simplify,
    tustin,
    doubleeuler,
    s,
    z,
    print_c_array,
    show_construction

# export symbolics2sym

const sp = SymPy.PyCall.PyNULL()
const s = SymPy.Sym("s")
const z = SymPy.Sym("z")

function __init__()
    # global const s = Sym("s")
    # global const z = Sym("z")
    copy!(sp, SymPy.sympy)
end

function SymPy.Sym(sys::StateSpace{<:Any,Sym})
    A, B, C, D = ControlSystems.ssdata(sys)
    if isdiscrete(sys)
        (C*inv(z * I(size(A, 1)) - A)*B+D)[1]
    else
        (C*inv(s * I(size(A, 1)) - A)*B+D)[1]
    end
end

function SymPy.Sym(sys::TransferFunction)
    if isdiscrete(sys)
        ssys = sp.Poly(numvec(sys)[], z) / sp.Poly(denvec(sys)[], z)
    else
        ssys = sp.Poly(numvec(sys)[], s) / sp.Poly(denvec(sys)[], s)
    end
end

"""
    ControlSystems.tf(sys::Sym, h = nothing)

Convert a symbolic, rational expression into a transfer function. `h` denotes the sample time if the system is discrete.
"""
function ControlSystems.tf(sys::Sym, h = nothing)
    n, d = sp.fraction(sys)
    d = sp.Poly(d)
    # d = d.monic() # Don't do this here
    n = n isa Number ? n : sp.Poly(n)
    if h === nothing
        tf(expand_coeffs(n, s), expand_coeffs(d, s))
    else
        tf(expand_coeffs(n, z), expand_coeffs(d, z), h)
    end
end

function expand_coeffs(n, var; numeric = false)
    n = sp.Poly(n, var)
    deg = n.degree() |> Float64 |> Int
    c = n.all_coeffs() # This fails if the coeffs are symbolic
    numeric && (c = Float64.(c))
    [c; zeros(eltype(c), deg - length(c) + 1)]
end
expand_coeffs(n::Real, args...; numeric = false) = n

function ControlSystems.minreal(sys::TransferFunction{<:Any,<:ControlSystems.SisoTf{<:Sym}})
    Sym(sys) |> simplify |> tf
end

function ControlSystems.tf(sys::StateSpace{<:Any,Sym})
    n, p = simplify.(sp.Poly.(simplify.(sp.fraction(simplify(Sym(sys)))), s))
    tf(simplify(n / p))
end

function Symbolics.Num(sys::StateSpace{<:Any,Symbolics.Num})
    A,B,C,D = ssdata(sys)
    if isdiscrete(sys)
        sz = Symbolics.@variables z
    else
        sz = Symbolics.@variables s
    end
    sz = only(sz)
    ex = (C*inv(sz * I(size(A, 1)) - A)*B + D)
    # if length(ex) == 1
    #     return ex[1] # don't do this as this complicates later function building
    # else
    #     return ex
    # end
    # ex = Symbolics.expand(ex)
    # ds = Symbolics.degree(ex, sz)
    # Symbolics.get_variables(ex)
    # Symbolics.polynormalize(ex)
end

function ControlSystems.tf(sys::StateSpace{<:Any,Symbolics.Num})
    n, p = simplify.(sp.Poly.(simplify.(sp.fraction(simplify(Sym(sys)))), s))
    tf(simplify(n / p))
end

function ControlSystems.minreal(sys::StateSpace{<:Any,Symbolics.Num})
    # sys |> Symbolics.Num .|> Symbolics.symbolics_to_sympy .|> sp.simplify
    nsys = Symbolics.Num(sys)
    nsys = Symbolics.simplify.(nsys)
    nsys = Symbolics.simplify_fractions.(nsys)
end

sym2num(P::TransferFunction, args...) = sym2num(Sym(P), args...)

"""
    sym2num(G, h::Real, pairs::Pair...)

Replace symbols by numbers. `h` indicates the sample rate if the system is discrete.
"""
function sym2num(P::Sym, h::Real, pairs::Pair...)
    P.has(s) && error(
        "Found `z` in symbolic expression, provide sample time as second argument to `sym2num`",
    )
    for (sym, val) in pairs
        P = subs(P, (sym, val))
    end
    n, d = sp.Poly.(sp.fraction(P), z)
    nn = float.(expand_coeffs(n, z))
    nd = float.(expand_coeffs(d, z))
    if nd[1] != 1 # make monic
        nn, nd = monic(nn, nd)
    end
    G = tf(nn, nd isa Number ? [nd] : nd, h)
    G
end

"""
    sym2num(G, pairs::Pair...)

Replace symbols by numbers.
"""
function sym2num(P::Sym, pairs::Pair...)
    # for i = 1:2
    P.has(z) && error(
        "Found `z` in symbolic expression, provide sample time as second argument to `sym2num`",
    )
    P.has(s) || error("Found no `s` in symbolic expression, provide")
    for (sym, val) in pairs
        P = subs(P, (sym, val))
    end
    n, d = sp.Poly.(sp.fraction(P), s)
    nn = float.(expand_coeffs(n, s))
    nd = float.(expand_coeffs(d, s))
    if nd[1] != 1 # make monic
        nn, nd = monic(nn, nd)
    end
    G = tf(nn, nd isa Number ? [nd] : nd)
    G
end

"""
    tustin(G::LTISystem, h)

Discretize `G` using the Tustin (bi-linear) method for sample time `h`.
"""
function tustin(G::LTISystem, h)
    iscontinuous(G) || throw(DomainError("Cannot discretize a discrete-domain system."))
    tustin(Sym(G), h)
end

function tustin(Gs::Sym, h)
    hh = Sym("hh")
    Gt = Gs.subs(s, hh * (z - 1) / (z + 1)) |> simplify
    @info "Perforing pole-zero cancellation (this might take a while)"
    @time Gt = sp.cancel(Gt) # Simplify takes forever with floating point values
    @time Gt = sp.collect(Gt, z) # Simplify takes forever with floating point values
    Gt = Gt.subs(hh, (2 / h))
    Gtf = tf(Gt, h)
    n, d = numvec(Gtf)[], denvec(Gtf)[]
    d[1] == 1 && (return Gtf)
    n, d = monic(n, d)
    tf(n, d, h)
end

"""
    doubleeuler(sys::AbstractStateSpace{<:Continuous}, Ts)

Discretize `sys` with a second-order approximation to the exponential map (Forward Euler is a first-order approximation). This is useful if the symbolic expressions for the true ZoH-discretization becomes too complicated. A second-order approximation is in many cases indistinguishable from the true ZoH-discretization.

# Example
```julia
w = 2pi .* exp10.(LinRange(-1, log10(25), 400))
sys = ssrand(1,1,4)
bodeplot(sys, w, lab="cont")
bodeplot!(c2d(sys, 0.001, :fwdeuler), w, label="fwdeuler")
bodeplot!(doubleeuler(sys, 0.001), w, label="double euler")
```
"""
function doubleeuler(sys::AbstractStateSpace{<:Continuous}, Ts)
    A, B, C, D = ssdata(sys)
    T = promote_type(eltype.((A, B, C, D))...)
    ny, nu = size(sys)
    nx = sys.nx
    Ad, Bd, Cd, Dd = (I + Ts * A + Ts^2 / 2 * A^2), Ts * B, C, D
    ss(Ad, Bd, Cd, Dd, Ts)
end

"""
    ccode(G; simplify = identity, cse = true)

Return a string with C-code for filtering a signal `u` through `G`.

# Arguments:
- `G`: A linear system
- `simplify`: A function for symbolic simplification. You may try `Sympy.simplify`, but for large systems, this will take a long time to compute.
- `cse`: Perform common subexpression elimination. This generally improvems the performance of the generated code.
"""
function ccode(G::TransferFunction; simplify = identity, cse = true)
    P = Sym(G)
    P.has(z) || error("Did not find `z` in symbolic expression")
    P.has(s) && error("Found `s` in symbolic expression, provide expression in `z`")
    n, d = numvec(G)[], denvec(G)[]
    if d[1] != 1
        @info "Denominator polynomial not monic, dividing by the leading coefficient."
        n = simplify.(n ./ d[1])
        d = simplify.(d ./ d[1])
    end
    @info "Calling free_symbols"
    @info vars = P.free_symbols
    vars.remove(z)
    vars = collect(vars)
    vars = sort(vars, by = string)
    var_str = ""
    for var in vars
        var_str *= ", double $(var)"
    end

    nu, ny = length(n), length(d)
    code = """
#include <stdio.h>\n
#include <math.h>\n
double transfer_function(double ui$(var_str)) {
    static double u[$(nu)] = {0};
    static double y[$(ny)] = {0};
    int i;
    for (i=$(nu-1); i > 0; --i) {
        u[i] = u[i-1];
    }
    u[0] = ui;
    for (i=$(ny-1); i > 0; --i) {
        y[i] = y[i-1];
    }
    y[0] = 0;
"""
    if cse
        @info "Finding common subexpressions"
        subex, final = sp.cse([n; d])
        if length(final) == 1
            final = final[]
        end
        n = final[1:length(n)]
        d = final[length(n)+1:end]
        for se in subex
            code *= "    double $(se[1]) = $(sp.ccode(se[2]));\n"
        end
    end
    for (i, n) in enumerate(n)
        @info "Writing numerator term $i/$(length(n))"
        code *= "    y[0] += ($(sp.ccode(n)))*u[$(i-1)];\n"
    end
    for (i, n) in enumerate(d[2:end])
        @info "Writing denominator term $i/$(length(d)-1)"
        code *= "    y[0] += ($(sp.ccode(-n)))*y[$(i)];\n"
    end
    code *= "    return y[0];\n}"
    println(code)
    try
        clipboard(code)
    catch
    end
    code
end

function ccode(sys::StateSpace{<:Discrete}; cse = true, function_name = "transfer_function")
    sys.nu == 1 || throw(ArgumentError("Multiple input not yet supported"))
    nx = sys.nx
    ny = sys.ny
    u = Sym("u")
    x = [Sym("x[$(i-1)]") for i = 1:nx]
    # @show P
    if ControlSystems.numeric_type(sys) <: SymPy.Sym
        P = Sym(sys)
        vars = P.free_symbols
        vars.remove(z)
        vars = collect(vars)
        vars = sort(vars, by = string)
        var_str = ""
        for var in vars
            var_str *= ", double $(var)"
        end
    else
        var_str = ""
    end
    x1 = zeros(Sym, nx) # workaround for strange bug with undefinied referece appearing in Pluto only
    y = zeros(Sym, ny)
    @show x1 = mul!(x1, sys.A, x) + sys.B * u
    @show y = mul!(y, sys.C, x) + sys.D * u
    # @show y = sp.collect.(y, x)

    code = """
#include <stdio.h>\n
#include <math.h>\n
void $(function_name)(double *y, double u$(var_str)) {
    static double x[$(nx)] = {0};  // Current state
    double xp[$(nx)] = {0};        // Next state
    int i;
"""
    if cse
        @info "Finding common subexpressions"
        subex, final = sp.cse([x1; y])
        x1 = final[][1:length(x1)]
        y = final[][length(x1)+1:end]
        for se in subex
            code *= "    double $(se[1]) = $(sp.ccode(se[2]));\n"
        end
    end
    code *= "\n    // Advance the state xp = Ax + Bu\n"
    for (i, n) in enumerate(x1)
        code *= "    xp[$(i-1)] = ($(sp.ccode(n)));\n"
    end
    code *= """

        // Accumulate the output y = C*x + D*u
    """

    for (i, n) in enumerate(y)
        code *= "    y[$(i-1)] = ($(sp.ccode(n)));\n"
    end
    code *= """

        // Make the predicted state the current state
        for (i=0; i < $(nx); ++i) {
            x[i] = xp[i];
        }
    """
    code *= "\n}"

    println(code)
    try
        clipboard(code)
    catch
    end
    code
end

function structured_text(code)

    code = replace(code, "pow" => "EXPT")
    code = replace(code, "static " => "")

    code = replace(code, r"double (\w+) =" => s"\1 : LREAL :=")
    code = replace(code, r"double (\w+)" => s"\1 : LREAL")
    code = replace(code, r"(.{1,20}?) \+= " => s"\1 := \1 + ")
    code = replace(code, " = " => " := ")
    println(code)
    try
        clipboard(code)
    catch
    end
    code
end

latextf(x::TransferFunction, args...) = latextf(Sym(x), args...)

Latexify.Latexify(x::TransferFunction) = latextf(x, args...)

function monic(n, d)
    n = n ./ d[1]
    d = d ./ d[1]
    if n isa Number
        n = [n]
    end
    n, d
end

"""
    latextf(x, mon = true)

Return a latex string representing `x`. If `mon`, then the denominator polynomial will be made monic (leading coefficient = 1).
"""
function latextf(x::Sym, mon = true)
    var = x.has(z) ? z : s
    n, d = sp.fraction(x)
    if mon
        n, d = monic(n, d)
    end
    n, d = sp.collect.((n, d), var)
    str = "\$\\dfrac{$((n))}{$((d))}\$"
    str = replace(str, '*' => "")
    try
        clipboard(code)
    catch
    end
    str
end

show_construction(sys::LTISystem; kwargs...) = show_construction(stdout, sys; kwargs...)
function show_construction(io::IO, sys::LTISystem; letb = true)
    # sys = StateSpace(sys)
    letb && println(io, "sys = let")
    prestr = letb ? "    " : ""
    println(io, prestr*"tempA = ", sys.A)
    println(io, prestr*"tempB = ", sys.B)
    println(io, prestr*"tempC = ", sys.C)
    println(io, prestr*"tempD = ", sys.D)
    letb || print(io, "sys = ")
    if isdiscrete(sys)
        println(io, prestr*"ss(tempA, tempB, tempC, tempD, $(sys.Ts))")
    else
        println(io, prestr*"ss(tempA, tempB, tempC, tempD)")
    end
    letb && println(io, "end")
    nothing
end

"""
sys2vec = @(sys) [
        size(sys.A,1);
        size(sys.B,2);
        size(sys.C,1);
        sys.A(:);
        sys.B(:);
        sys.C(:);
        sys.D(:)
    ]
"""
function vec2sys(v, ts = nothing)
    nx = Int(v[1])
    nu = Int(v[2])
    ny = Int(v[3])
    ai = (1:nx^2) .+ 3
    bi = (1:nx*nu) .+ ai[end]
    ci = (1:nx*ny) .+ bi[end]
    di = (1:nu*ny) .+ ci[end]
    A = reshape(v[ai], nx, nx)
    B = reshape(v[bi], nx, nu)
    C = reshape(v[ci], ny, nx)
    D = reshape(v[di], ny, nu)
    ts === nothing ? ss(A, B, C, D) : ss(A, B, C, D, ts)
end

print_c_array(a::AbstractArray, args...; kwargs...) =
    print_c_array(stdout, a, args...; kwargs...)
function print_c_array(io, a::AbstractVector, name = "vec"; cse = false, s = "", struct_name=nothing)
    l = length(a)
    a = float.(a)
    if cse
        a = write_array_cse(io, a, name, s)
    end
    if struct_name === nothing
        s == "" && println(io, "double $name[$l];")
        struct_name = ""
    else
        struct_name = struct_name*"->"
    end
    for i = 1:l
        println(io, "$(struct_name)" * s * "$name[$(i-1)] = $(a[i]);")
    end
    println(io)
end
function print_c_array(io, a::AbstractMatrix, name = "mat"; cse = false, s = "", struct_name=nothing)
    r, c = size(a)
    a = float.(a)
    if cse
        a = write_array_cse(io, a, name, s)
    end
    if struct_name === nothing
        s == "" && println(io, "double $name[$r][$c];")
        struct_name = ""
    else
        struct_name = struct_name*"->"
    end
    for i = 1:r, j = 1:c
        println(io, "$(struct_name)" * s * "$name[$(i-1)][$(j-1)] = $(a[i,j]);")
    end
    println(io)
end

function print_c_array(io, a::AbstractArray{<:Any,3}, name = "array"; cse = false, s = "", struct_name=nothing)
    r, c, d = size(a)
    a = float.(a)
    if cse
        a = write_array_cse(io, a, name, s)
    end
    if struct_name === nothing
        s == "" && println(io, "double $name[$r][$c][$d];")
        struct_name = ""
    else
        struct_name = struct_name*"->"
    end
    for i = 1:r, j = 1:c, k = 1:d
        println(io, "$(struct_name)" * s * "$name[$(i-1)][$(j-1)][$(k-1)] = $(a[i,j,k]);")
    end
    println(io)
end

# Interpolation
function print_c_array(
    io,
    a::Vector{<:AbstractArray},
    t::AbstractVector,
    name          = "mat";
    cse           = false,
    s             = "",
    print_vector  = true,
    print_logic   = true,
    struct_name::Union{Nothing, String} = nothing,
    struct_type   = nothing,
    ivecname = name * "_interp_vect",
)
    length(a) == length(t) || throw(ArgumentError("length mismatch between a and t"))
    print_vector && print_c_array(io, t, ivecname; cse, s, struct_name)
    # for i in eachindex(t)
    #     iname = name*"_$(i-1)"
    #     print_c_array(io, a[i], iname; cse, s)
    # end
    a = cat(a..., dims = 3)

    if struct_name !== nothing
        struct_type isa String || throw(ArgumentError("Expected a struct type since you provided struct name"))
        print_c_array(io, a, name * "_interp"; cse, s, struct_name)
        println(io, "void interpolate_$(name)($(struct_type) *$(struct_name), double t) {")
        ivecname = struct_name*"->"*ivecname
    else
        print_c_array(io, a, name * "_interp"; cse, s)
        println(
        io,
        "void interpolate_$(name)(double *$(name), double *$(name)_interp, double *$(ivecname), double t) {"
        )
    end

    println(
        io,
        """
    int k = 0;
    for (k = 0; k < $(length(t)-2); ++k) { // Loop to find correct interval
        if (t < $(ivecname)[k+1]) {
            break;
        }
    } // t is now between indices k and k+1 modolu edge cases below
    double t0 = $(ivecname)[k];
    double t1 = $(ivecname)[k+1];
    double l = t1 - t0;      // Length of interval between points
    double alpha = (t-t0)/l; // Between 0 and 1, proportion of k+1 point
    if (t < $(ivecname)[0]) { // edge cases
        alpha = 0;
    }else if (t > $(ivecname)[$(length(t)-1)]) {
        alpha = 1;
    }
$(interpolator_string(a, name, struct_name, struct_type))}
""",
    )
end

function interpolator_string(a::AbstractArray{<:Any,3}, name, struct_name, struct_type)
    if struct_name === nothing
        """
            for (int i = 0; i < $(size(a, 1)); ++i) {     // $name has dimensions $((size(a)[1:2]))
                for (int j = 0; j < $(size(a, 2)); ++j) { // $(name)_interp has dimensions $((size(a)))
                    $(name)[i][j] = (1-alpha)*$(name)_interp[i][j][k] + alpha*$(name)_interp[i][j][k+1];
                }
            }
        """
    else
        name = struct_name*"->"*name
        """
            for (int i = 0; i < $(size(a, 1)); ++i) {     // $name has dimensions $((size(a)[1:2]))
                for (int j = 0; j < $(size(a, 2)); ++j) { // $(name)_interp has dimensions $((size(a)))
                    $(name)[i][j] = (1-alpha)*$(name)_interp[i][j][k] + alpha*$(name)_interp[i][j][k+1];
                }
            }
        """
    end
end

# StateSpace
function print_c_array(io, sys::AbstractStateSpace, args...; s = "", en = "", struct_name=nothing)
    print_c_array(io, sys.A, "A" * en; s, struct_name)
    print_c_array(io, sys.B, "B" * en; s, struct_name)
    print_c_array(io, sys.C, "C" * en; s, struct_name)
    print_c_array(io, sys.D, "D" * en; s, struct_name)
end

# Interpolated StateSpace
function print_c_array(
    io,
    sys::Vector{<:AbstractStateSpace},
    t::AbstractVector,
    name = "sys";
    cse = false,
    s = "",
    en = "",
    struct_name::Union{Nothing, String} = nothing,
    struct_type   = nothing,
)
    length(sys) == length(t) || throw(ArgumentError("length mismatch between a and t"))
    ivecname = name * en * "_interp_vect"
    print_c_array(io, t, ivecname; cse, s, struct_name)
    for p in (:A, :B, :C, :D)
        A = getproperty.(sys, p)
        print_c_array(io, A, t, name * "_" * string(p) * en; cse, s, print_vector = false, struct_name, struct_type, ivecname)
    end
end

function write_array_cse(io, a, name = "x", s = "")
    subex, final = sp.cse(a, symbols = [SymPy.Sym(name * string(i)) for i = 1:200])
    new_a = final[]
    for se in subex
        println(io, "double $(se[1]) = $(sp.ccode(se[2]));")
    end
    new_a
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BinomialProcess(n)

A Binomial point process with `n` points.
"""
struct BinomialProcess <: PointProcess
  n::Int
end

ishomogeneous(p::BinomialProcess) = true

struct BinomialSampling end

default_sampling_algorithm(::BinomialProcess) = BinomialSampling()

function rand_single(p::BinomialProcess, b::Box{Dim,T},
                     ::BinomialSampling) where {Dim,T}
  # region configuration
  lo, up = coordinates.(extrema(b))

  # product of uniform distributions
  U = product_distribution([Uniform(lo[i], up[i]) for i in 1:Dim])

  # return point pattern
  PointSet(rand(U, p.n))
end


text:
function search(env::SearchEnv, request; exclude=nothing)

    # Parse the input from the request
    parsedinput = parse_input(env, request)
    issearch = !isempty(parsedinput.search)
    isfilter = !isempty(parsedinput.filter)

    if !isfilter
        # No filter, search always done
        results = search(env.searchers,
                         parsedinput.search;
                         search_method=request.search_method,
                         max_matches=request.max_matches,
                         max_suggestions=request.max_suggestions,
                         custom_weights=request.custom_weights)
    elseif !issearch
        # No search, filter only
        idxs_filt = indexfilter(env.dbdata,
                                parsedinput.filter;
                                id_key=env.id_key,
                                sort_keys=request.sort_fields,
                                sort_reverse=request.sort_reverse,
                                exclude=exclude)
        result = build_result_from_ids(env.dbdata,
                                       idxs_filt,
                                       env.id_key,
                                       make_id();
                                       id_key=env.id_key,
                                       max_matches=request.max_matches)
        results = [result]
    elseif issearch
        # Search and filter search results
        results = search(env.searchers,
                         parsedinput.search;
                         search_method=request.search_method,
                         max_matches=request.max_matches,
                         max_suggestions=request.max_suggestions,
                         custom_weights=request.custom_weights)
        #TODO(Corneliu) Decide whether to do multiple searches with
        # higher max_matches (needs InternalRequest specification) to
        # match filter
        idxs_filt = indexfilter(env.dbdata,
                                parsedinput.filter;
                                id_key=env.id_key,
                                sort_keys=request.sort_fields,
                                sort_reverse=request.sort_reverse,
                                exclude=exclude)

        # Filter out ids that are not present in the the filtered ids
        for result in results
            filter!(score_idx -> in(score_idx[2], idxs_filt), result.query_matches)
        end
    end

    return results
end


text:
@time using CSV, XLSX, DataFrames, StatsBase, Statistics
@time using Plots, Plots.PlotMeasures, LaTeXStrings, StatsPlots
default()

over50(ndwi) = vec(maximum(ndwi, dims = 1) .> 50)

schedules = DataFrame(XLSX.readtable("C:/Users/rmsms/OneDrive/lab/aespa//schedule.xlsx", "schedule")...)

ntwk = CSV.read("D:/simulated/ntwk.csv", DataFrame)

import_dir = "D:/simulated/"
export_dir = "C:/Users/rmsms/OneDrive/lab/aespa/png/"

color_ = Dict("00" => :black, "0V" => colorant"#C00000", "M0" => colorant"#0070C0", "MV" => colorant"#7030A0")
shape_ = Dict("00" => :+, "0V" => :x, "M0" => :circle, "MV" => :star5)
label_ = Dict("00" => "No control", "0V" => "Vaccination", "M0" => "Restriction", "MV" => "Both control")

te_color_ = Dict(
    "S" => colorant"#FFC000",
    "E" => colorant"#ED7D31",
    "I" => colorant"#C00000",
    "R" => colorant"#A5A5A5",
    "V" => colorant"#70AD47")

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
struct TikzAxis <: AbstractTikzAxis
	plots::Vector{AbstractTikzPlot}
end

function TikzAxis()
	return TikzAxis([])
end

function EmptyTikzAxis()
	return TikzAxis()
end


text:
using Distributions, DiffEqBayes

println("Starting the test")
@test DiffEqBayes.stan_string(Bernoulli(1)) == "bernoulli(1.0)"
@test DiffEqBayes.stan_string(Bernoulli) == "bernoulli"

@test DiffEqBayes.stan_string(Binomial(5,0.3)) == "binomial(5, 0.3)"
@test DiffEqBayes.stan_string(Binomial) == "binomial"

@test DiffEqBayes.stan_string(BetaBinomial(5,1,1)) == "beta_binomial(5, 1.0, 1.0)"
@test DiffEqBayes.stan_string(BetaBinomial) == "beta_binomial"

@test DiffEqBayes.stan_string(Hypergeometric(5,5,3)) == "hypergeometric(5, 5, 3)"
@test DiffEqBayes.stan_string(Hypergeometric) == "hypergeometric"

@test DiffEqBayes.stan_string(NegativeBinomial(5,0.3)) == "neg_binomial(5.0, 0.3)"
@test DiffEqBayes.stan_string(NegativeBinomial) == "neg_binomial"

@test DiffEqBayes.stan_string(Poisson(5)) == "poisson(5.0)"
@test DiffEqBayes.stan_string(Poisson) == "poisson"

@test DiffEqBayes.stan_string(Normal(0,1)) == "normal(0.0, 1.0)"
@test DiffEqBayes.stan_string(Normal) == "normal"

@test DiffEqBayes.stan_string(TDist(5)) == "student_t(5.0,0,1)"
@test DiffEqBayes.stan_string(TDist) == "student_t"

@test DiffEqBayes.stan_string(Cauchy(0,1)) == "cauchy(0.0, 1.0)"
@test DiffEqBayes.stan_string(Cauchy) == "cauchy"

@test DiffEqBayes.stan_string(Laplace(0,1)) == "double_exponential(0.0, 1.0)"
@test DiffEqBayes.stan_string(Laplace) == "double_exponential"

@test DiffEqBayes.stan_string(Distributions.Logistic(0,1)) == "logistic(0.0, 1.0)"
@test DiffEqBayes.stan_string(Distributions.Logistic) == "logistic"

@test DiffEqBayes.stan_string(Gumbel(0,1)) == "gumbel(0.0, 1.0)"
@test DiffEqBayes.stan_string(Gumbel) == "gumbel"

@test DiffEqBayes.stan_string(LogNormal(0,1)) == "lognormal(0.0, 1.0)"
@test DiffEqBayes.stan_string(LogNormal) == "lognormal"

@test DiffEqBayes.stan_string(Chisq(5)) == "chi_square(5.0)"
@test DiffEqBayes.stan_string(Chisq) == "chi_square"

@test DiffEqBayes.stan_string(Exponential(5)) == "exponential(5.0)"
@test DiffEqBayes.stan_string(Exponential) == "exponential"

@test DiffEqBayes.stan_string(Gamma(2,3)) == "gamma(2.0, 3.0)"
@test DiffEqBayes.stan_string(Gamma) == "gamma"

@test DiffEqBayes.stan_string(InverseGamma(2,3)) == "inv_gamma(2.0, 3.0)"
@test DiffEqBayes.stan_string(InverseGamma) == "inv_gamma"

@test DiffEqBayes.stan_string(Weibull(1,1)) == "weibull(1.0, 1.0)"
@test DiffEqBayes.stan_string(Weibull) == "weibull"

@test DiffEqBayes.stan_string(Frechet(1,1)) =="frechet(1.0, 1.0)"
@test DiffEqBayes.stan_string(Frechet) == "frechet"

@test DiffEqBayes.stan_string(Rayleigh(5)) == "rayleigh(5.0)"
@test DiffEqBayes.stan_string(Rayleigh) == "rayleigh"

@test DiffEqBayes.stan_string(Pareto(2,3)) == "pareto(2.0, 3.0)"
@test DiffEqBayes.stan_string(Pareto) == "pareto"

@test DiffEqBayes.stan_string(GeneralizedPareto(0,1,2)) == "pareto_type_2(0.0, 1.0, 2.0)"
@test DiffEqBayes.stan_string(GeneralizedPareto) == "pareto_type_2"

@test DiffEqBayes.stan_string(Beta(3,3)) == "beta(3.0, 3.0)"
@test DiffEqBayes.stan_string(Beta) == "beta"

@test DiffEqBayes.stan_string(Uniform(1,4)) =="uniform(1.0, 4.0)"
@test DiffEqBayes.stan_string(Uniform) == "uniform"

@test DiffEqBayes.stan_string(VonMises(0,2)) == "von_mises(0.0, 2.0)"
@test DiffEqBayes.stan_string(VonMises) == "von_mises"


text:
const AGGR2STR = Dict{Symbol,String}(:add => "∑", :sub => "-∑", :mul => "∏", :div => "1/∏",
                                     :max => "max", :min => "min", :mean => "𝔼[]")

"""
    GCNConv([graph, ]in=>out)
    GCNConv([graph, ]in=>out, σ)

Graph convolutional layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs)
or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs). Is optionnal so you can give a `FeaturedGraph` to
the layer instead of only the features.
- `in`: the dimension of input features.
- `out`: the dimension of output features.
- `bias::Bool=true`: keyword argument, whether to learn the additive bias.

Data should be stored in (# features, # nodes) order.
For example, a 1000-node graph each node of which poses 100 features is constructed.
The input data would be a `1000×100` array.
"""
struct GCNConv{T,F,S<:AbstractFeaturedGraph}
    weight::AbstractMatrix{T}
    bias::AbstractVector{T}
    σ::F
    fg::S
end

function GCNConv(ch::Pair{<:Integer,<:Integer}, σ = identity;
                 init=glorot_uniform, T::DataType=Float32, bias::Bool=true)
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    fg = NullGraph()
    GCNConv(T.(init(ch[2], ch[1])), b, σ, fg)
end

function GCNConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, σ = identity;
                 init=glorot_uniform, T::DataType=Float32, bias::Bool=true)
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    fg = FeaturedGraph(adj)
    GCNConv(T.(init(ch[2], ch[1])), b, σ, fg)
end

function GCNConv(fg::FeaturedGraph, ch::Pair{<:Integer,<:Integer}, σ = identity;
                 init=glorot_uniform, T::DataType=Float32, bias::Bool=true)
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    GCNConv(T.(init(ch[2], ch[1])), b, σ, fg)
end

@functor GCNConv

function (g::GCNConv)(A::AbstractMatrix, X::AbstractMatrix)
    L̃ = normalized_laplacian(A, eltype(X); selfloop=true)
    L̃ = convert(typeof(X), L̃)  # ensure L has the same type as X, especially X::CuArray

    g.σ.(g.weight * X * L̃ .+ g.bias)
end

(g::GCNConv)(A::AbstractMatrix, X::Transpose{S,R}) where {S,R<:AbstractMatrix} = g(A, R(X))

function (g::GCNConv)(X::AbstractMatrix{T}) where {T}
    @assert has_graph(g.fg) "A GCNConv created without a graph must be given a FeaturedGraph as an input."
    A = adjacency_matrix(g.fg)
    g(A, X)
end

function (g::GCNConv)(fg::FeaturedGraph)
    X = node_feature(fg)
    A = adjacency_matrix(fg) # TODO: choose graph from g or fg
    Zygote.ignore() do
        g.fg isa NullGraph || (g.fg.graph = A)
    end
    X_ = g(A, X)
    FeaturedGraph(A, X_)
end

function Base.show(io::IO, l::GCNConv)
    in_channel = size(l.weight, ndims(l.weight))
    out_channel = size(l.weight, ndims(l.weight)-1)
    print(io, "GCNConv(G(V=", nv(l.fg))
    print(io, ", E), ", in_channel, "=>", out_channel)
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ")")
end

"""
    ChebConv([graph, ]in=>out, k)

Chebyshev spectral graph convolutional layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`,
`SimpleWeightedDiGraph` (from SimpleWeightedGraphs). Is optionnal so you can give a `FeaturedGraph` to
the layer instead of only the features.
- `in`: the dimension of input features.
- `out`: the dimension of output features.
- `k`: the order of Chebyshev polynomial.
- `bias::Bool=true`: keyword argument, whether to learn the additive bias.
"""
struct ChebConv{T,S<:AbstractFeaturedGraph}
    weight::AbstractArray{T,3}
    bias::AbstractVector{T}
    fg::S
    k::Integer
    in_channel::Integer
    out_channel::Integer
end

function ChebConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, k::Integer;
                  init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    b = bias ? init(ch[2]) : zeros(T, ch[2])
    fg = FeaturedGraph(adj)
    ChebConv(init(ch[2], ch[1], k), b, fg, k, ch[1], ch[2])
end

function ChebConv(ch::Pair{<:Integer,<:Integer}, k::Integer;
                  init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    b = bias ? init(ch[2]) : zeros(T, ch[2])
    fg = NullGraph()
    ChebConv(init(ch[2], ch[1], k), b, fg, k, ch[1], ch[2])
end

@functor ChebConv

function (c::ChebConv)(L̃::AbstractMatrix{S}, X::AbstractMatrix{T}) where {S<:Real, T<:Real}
    @assert size(X, 1) == c.in_channel "Input feature size must match input channel size."
    @assert size(X, 2) == size(L̃, 1) "Input vertex number must match Laplacian matrix size."

    Z_prev = X
    Z = X * L̃
    Y = view(c.weight,:,:,1) * Z_prev
    Y += view(c.weight,:,:,2) * Z
    for k = 3:c.k
        Z, Z_prev = 2*Z*L̃ - Z_prev, Z
        Y += view(c.weight,:,:,k) * Z
    end
    return Y .+ c.bias
end

function (c::ChebConv)(X::AbstractMatrix{T}) where {T<:Real}
    @assert has_graph(c.fg) "A ChebConv created without a graph must be given a FeaturedGraph as an input."
    g = graph(c.fg)
    L̃ = scaled_laplacian(g, T)
    L̃ = convert(typeof(X), L̃)  # ensure L̃ has the same type as X, especially X::CuArray

    c(L̃, X)
end

(c::ChebConv)(X::Transpose{T,R}) where {T<:Real,R<:AbstractMatrix} = c(R(X))

function (c::ChebConv)(fg::FeaturedGraph)
    @assert has_graph(fg) "A given FeaturedGraph must contain a graph."
    g = graph(fg)
    Zygote.ignore() do
        c.fg isa NullGraph || (c.fg.graph = g)
    end
    X = node_feature(fg)
    L̃ = scaled_laplacian(adjacency_matrix(fg))
    L̃ = convert(typeof(X), L̃)  # ensure L has the same type as X, especially X::CuArray

    X_ = c(L̃, X)
    FeaturedGraph(g, X_)
end

function Base.show(io::IO, l::ChebConv)
    print(io, "ChebConv(G(V=", nv(l.fg))
    print(io, ", E), ", l.in_channel, "=>", l.out_channel)
    print(io, ", k=", l.k)
    print(io, ")")
end

"""
    GraphConv([graph, ]in=>out)
    GraphConv([graph, ]in=>out, σ)
    GraphConv([graph, ]in=>out, σ, aggr)

Graph neural network layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`,
`SimpleWeightedDiGraph` (from SimpleWeightedGraphs). Is optionnal so you can give a `FeaturedGraph` to
the layer instead of only the features.
- `in`: the dimension of input features.
- `out`: the dimension of output features.
- `bias::Bool=true`: keyword argument, whether to learn the additive bias.
- `σ=identity`: activation function.
- `aggr::Symbol=:add`: an aggregate function applied to the result of message function. `:add`, `:max` and `:mean` are available.
"""
struct GraphConv{V<:AbstractFeaturedGraph,T} <: MessagePassing
    fg::V
    weight1::AbstractMatrix{T}
    weight2::AbstractMatrix{T}
    bias::AbstractVector{T}
    σ
    aggr::Symbol
end

function GraphConv(el::AbstractVector{<:AbstractVector{<:Integer}},
                   ch::Pair{<:Integer,<:Integer}, σ=identity, aggr=:add;
                   init = glorot_uniform, bias::Bool=true, T::DataType=Float32)
    w1 = T.(init(ch[2], ch[1]))
    w2 = T.(init(ch[2], ch[1]))
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    fg = FeaturedGraph(el)
    GraphConv(fg, w1, w2, b, σ, aggr)
end

function GraphConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, σ=identity, aggr=:add;
                   init = glorot_uniform, bias::Bool=true, T::DataType=Float32)
    w1 = T.(init(ch[2], ch[1]))
    w2 = T.(init(ch[2], ch[1]))
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    fg = FeaturedGraph(adjacency_list(adj))
    GraphConv(fg, w1, w2, b, σ, aggr)
end

function GraphConv(ch::Pair{<:Integer,<:Integer}, σ=identity, aggr=:add;
                   init = glorot_uniform, bias::Bool=true, T::DataType=Float32)
    w1 = T.(init(ch[2], ch[1]))
    w2 = T.(init(ch[2], ch[1]))
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    GraphConv(NullGraph(), w1, w2, b, σ, aggr)
end

@functor GraphConv

message(g::GraphConv, x_i, x_j::AbstractVector, e_ij) = g.weight2 * x_j
update(g::GraphConv, m::AbstractVector, x::AbstractVector) = g.σ.(g.weight1*x .+ m .+ g.bias)
function (gc::GraphConv)(X::AbstractMatrix)
    @assert has_graph(gc.fg) "A GraphConv created without a graph must be given a FeaturedGraph as an input."
    g = graph(gc.fg)
    _, X = propagate(gc, adjacency_list(g), Fill(0.f0, 0, ne(g)), X, :add)
    X
end
(g::GraphConv)(fg::FeaturedGraph) = propagate(g, fg, :add)

function Base.show(io::IO, l::GraphConv)
    in_channel = size(l.weight1, ndims(l.weight1))
    out_channel = size(l.weight1, ndims(l.weight1)-1)
    print(io, "GraphConv(G(V=", nv(l.fg), ", E=", ne(l.fg))
    print(io, "), ", in_channel, "=>", out_channel)
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ", aggr=", AGGR2STR[l.aggr])
    print(io, ")")
end

"""
    GATConv([graph, ]in=>out)

Graph attentional layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`,
`SimpleWeightedDiGraph` (from SimpleWeightedGraphs). Is optionnal so you can give a `FeaturedGraph` to
the layer instead of only the features.
- `in`: the dimension of input features.
- `out`: the dimension of output features.
- `bias::Bool=true`: keyword argument, whether to learn the additive bias.
- `negative_slope::Real=0.2`: keyword argument, the parameter of LeakyReLU.
"""
struct GATConv{V<:AbstractFeaturedGraph, T <: Real} <: MessagePassing
    fg::V
    weight::AbstractMatrix{T}
    bias::AbstractVector{T}
    a::AbstractMatrix{T}
    negative_slope::Real
    channel::Pair{<:Integer,<:Integer}
    heads::Integer
    concat::Bool
end

function GATConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}; heads::Integer=1,
                 concat::Bool=true, negative_slope::Real=0.2, init=glorot_uniform,
                 bias::Bool=true, T::DataType=Float32)
    w = T.(init(ch[2]*heads, ch[1]))
    b = bias ? T.(init(ch[2]*heads)) : zeros(T, ch[2]*heads)
    a = T.(init(2*ch[2], heads))
    fg = FeaturedGraph(adjacency_list(adj))
    GATConv(fg, w, b, a, negative_slope, ch, heads, concat)
end

function GATConv(ch::Pair{<:Integer,<:Integer}; heads::Integer=1,
                 concat::Bool=true, negative_slope::Real=0.2, init=glorot_uniform,
                 bias::Bool=true, T::DataType=Float32)
    w = T.(init(ch[2]*heads, ch[1]))
    b = bias ? T.(init(ch[2]*heads)) : zeros(T, ch[2]*heads)
    a = T.(init(2*ch[2], heads))
    GATConv(NullGraph(), w, b, a, negative_slope, ch, heads, concat)
end

@functor GATConv

# Here the α that has not been softmaxed is the first number of the output message
function message(g::GATConv, x_i::AbstractVector, x_j::AbstractVector, e_ij)
    x_i = reshape(g.weight*x_i, :, g.heads)
    x_j = reshape(g.weight*x_j, :, g.heads)
    n = size(x_i, 1)
    e = vcat(x_i, x_j+zero(x_j))
    e = sum(e .* g.a, dims=1)  # inner product for each head, output shape: (1, g.heads)
    e = leakyrelu.(e, g.negative_slope)
    vcat(e, x_j)  # shape: (n+1, g.heads)
end

# After some reshaping due to the multihead, we get the α from each message,
# then get the softmax over every α, and eventually multiply the message by α
function apply_batch_message(g::GATConv, i, js, edge_idx, E::AbstractMatrix, X::AbstractMatrix, u)
    e_ij = hcat([message(g, get_feature(X, i), get_feature(X, j), get_feature(E, edge_idx[(i,j)])) for j = js]...)
    n = size(e_ij, 1)
    alphas = Flux.softmax(reshape(view(e_ij, 1, :), g.heads, :), dims=2)
    msgs = view(e_ij, 2:n, :) .* reshape(alphas, 1, :)
    reshape(msgs, (n-1)*g.heads, :)
end

function update_batch_edge(g::GATConv, adj, E::AbstractMatrix, X::AbstractMatrix, u)
    n = size(adj, 1)
    # In GATConv, a vertex must always receive a message from itself
    Zygote.ignore() do
        add_self_loop!(adj, n)
    end

    edge_idx = edge_index_table(adj)
    hcat([apply_batch_message(g, i, adj[i], edge_idx, E, X, u) for i in 1:n]...)
end

# The same as update function in batch manner
function update_batch_vertex(g::GATConv, M::AbstractMatrix, X::AbstractMatrix, u)
    M = M .+ g.bias
    if !g.concat
        N = size(M, 2)
        M = reshape(mean(reshape(M, :, g.heads, N), dims=2), :, N)
    end
    return M
end

function (gat::GATConv)(X::AbstractMatrix)
    @assert has_graph(gat.fg) "A GATConv created without a graph must be given a FeaturedGraph as an input."
    g = graph(gat.fg)
    _, X = propagate(gat, adjacency_list(g), Fill(0.f0, 0, ne(g)), X, :add)
    X
end
(g::GATConv)(fg::FeaturedGraph) = propagate(g, fg, :add)

function Base.show(io::IO, l::GATConv)
    in_channel = size(l.weight, ndims(l.weight))
    out_channel = size(l.weight, ndims(l.weight)-1)
    print(io, "GATConv(G(V=", nv(l.fg), ", E=", ne(l.fg))
    print(io, "), ", in_channel, "=>", out_channel)
    print(io, ", LeakyReLU(λ=", l.negative_slope)
    print(io, "))")
end

"""
    GatedGraphConv([graph, ]out, num_layers)

Gated graph convolution layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`,
`SimpleWeightedDiGraph` (from SimpleWeightedGraphs). Is optionnal so you can give a `FeaturedGraph` to
the layer instead of only the features.
- `out`: the dimension of output features.
- `num_layers` specifies the number of gated recurrent unit.
- `aggr::Symbol=:add`: an aggregate function applied to the result of message function. `:add`, `:max` and `:mean` are available.
"""
struct GatedGraphConv{V<:AbstractFeaturedGraph, T <: Real, R} <: MessagePassing
    fg::V
    weight::AbstractArray{T}
    gru::R
    out_ch::Integer
    num_layers::Integer
    aggr::Symbol
end

function GatedGraphConv(adj::AbstractMatrix, out_ch::Integer, num_layers::Integer;
                        aggr=:add, init=glorot_uniform, T::DataType=Float32)
    w = T.(init(out_ch, out_ch, num_layers))
    gru = GRUCell(out_ch, out_ch)
    fg = FeaturedGraph(adjacency_list(adj))
    GatedGraphConv(fg, w, gru, out_ch, num_layers, aggr)
end

function GatedGraphConv(out_ch::Integer, num_layers::Integer;
                        aggr=:add, init=glorot_uniform, T::DataType=Float32)
    w = T.(init(out_ch, out_ch, num_layers))
    gru = GRUCell(out_ch, out_ch)
    GatedGraphConv(NullGraph(), w, gru, out_ch, num_layers, aggr)
end

@functor GatedGraphConv

message(g::GatedGraphConv, x_i, x_j::AbstractVector, e_ij) = x_j
update(g::GatedGraphConv, m::AbstractVector, x) = m

function (ggc::GatedGraphConv)(X::AbstractMatrix{T}) where {T<:Real}
    @assert has_graph(ggc.fg) "A GraphConv created without a graph must be given a FeaturedGraph as an input."
    ggc(adjacency_list(ggc.fg), X)
end

function (ggc::GatedGraphConv{V,T})(fg::FeaturedGraph) where {V,T<:Real}
    g = graph(fg)
    H = ggc(adjacency_list(g), node_feature(fg))
    FeaturedGraph(g, H)
end

function (ggc::GatedGraphConv)(adj::AbstractVector{T}, X::AbstractMatrix{S}) where {T<:AbstractVector,S<:Real}
    H = X
    m, n = size(H)
    @assert (m <= ggc.out_ch) "number of input features must less or equals to output features."
    (m < ggc.out_ch) && (H = vcat(H, zeros(S, ggc.out_ch - m, n)))

    for i = 1:ggc.num_layers
        M = view(ggc.weight, :, :, i) * H
        _, M = propagate(ggc, adj, Fill(0.f0, 0, ne(adj)), M, :add)
        H, _ = ggc.gru(H, M)  # BUG: FluxML/Flux.jl#1381
    end
    H
end

function Base.show(io::IO, l::GatedGraphConv)
    print(io, "GatedGraphConv(G(V=", nv(l.fg), ", E=", ne(l.fg))
    print(io, "), (=>", l.out_ch)
    print(io, ")^", l.num_layers)
    print(io, ", aggr=", AGGR2STR[l.aggr])
    print(io, ")")
end

"""
    EdgeConv(graph, nn)
    EdgeConv(graph, nn, aggr)

Edge convolutional layer.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `nn`: a neural network
- `aggr::Symbol=:max`: an aggregate function applied to the result of message function. `:add`, `:max` and `:mean` are available.
"""
struct EdgeConv{V<:AbstractFeaturedGraph} <: MessagePassing
    fg::V
    nn
    aggr::Symbol
end

function EdgeConv(adj::AbstractMatrix, nn; aggr::Symbol=:max)
    fg = FeaturedGraph(adjacency_list(adj))
    EdgeConv(fg, nn, aggr)
end

function EdgeConv(nn; aggr::Symbol=:max)
    EdgeConv(NullGraph(), nn, aggr)
end

@functor EdgeConv

message(e::EdgeConv, x_i::AbstractVector, x_j::AbstractVector, e_ij) = e.nn(vcat(x_i, x_j .- x_i))
update(e::EdgeConv, m::AbstractVector, x) = m

function (e::EdgeConv)(X::AbstractMatrix)
    @assert has_graph(e.fg) "A EdgeConv created without a graph must be given a FeaturedGraph as an input."
    g = graph(e.fg)
    _, X = propagate(e, adjacency_list(g), Fill(0.f0, 0, ne(g)), X, e.aggr)
    X
end

(e::EdgeConv)(fg::FeaturedGraph) = propagate(e, fg, e.aggr)

function Base.show(io::IO, l::EdgeConv)
    print(io, "EdgeConv(G(V=", nv(l.fg), ", E=", ne(l.fg))
    print(io, "), ", l.nn)
    print(io, ", aggr=", AGGR2STR[l.aggr])
    print(io, ")")
end


text:
export MultiArmBanditsEnv

mutable struct MultiArmBanditsEnv <: AbstractEnv
    true_reward::Float64
    true_values::Vector{Float64}
    rng::AbstractRNG
    # cache
    reward::Float64
    is_terminated::Bool
end

"""
    MultiArmBanditsEnv(;true_reward=0., k = 10,rng=Random.GLOBAL_RNG)

`true_reward` is the expected reward. `k` is the number of arms. See
[multi-armed bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit) for more
detailed explanation.

This is a **one-shot** game. The environment terminates immediately after taking
in an action. Here we use it to demonstrate how to write a customized
environment with only minimal interfaces defined.
"""
function MultiArmBanditsEnv(; true_reward = 0.0, k = 10, rng = Random.GLOBAL_RNG)
    true_values = true_reward .+ randn(rng, k)
    MultiArmBanditsEnv(true_reward, true_values, rng, 0.0, false)
end

"""
First we need to define the action space. In the [`MultiArmBanditsEnv`](@ref)
environment, the possible actions are `1` to `k` (which equals to
`length(env.true_values)`).

!!! note
    Although we decide to return an action space of `Base.OneTo`  here, it is
    not a hard requirement. You can return anything else (`Tuple`,
    `Distribution`, etc) that is more suitable to describe your problem and
    handle it correctly in the `you_env(action)` function. Some algorithms may
    require that the action space must be of `Base.OneTo`. However, it's the
    algorithm designer's job to do the checking and conversion.
"""
RLBase.action_space(env::MultiArmBanditsEnv) = Space(OneTo(length(env.true_values)))

"""
In our design, the return of taking an action in `env` is **undefined**. This is
the main difference compared to those interfaces defined in
[OpenAI/Gym](https://github.com/openai/gym). We find that the async manner is
more suitable to describe many complicated environments. However, one of the
inconveniences is that we have to cache some intermediate data for future
queries. Here we have to store `reward` and `is_terminated` in the instance of
`env` for future queries.
"""
function (env::MultiArmBanditsEnv)(action)
    env.reward = randn(env.rng) + env.true_values[action]
    env.is_terminated = true
end

RLBase.is_terminated(env::MultiArmBanditsEnv) = env.is_terminated

"""
!!! warn
    If the `env` is not started yet, the returned value is meaningless. The
    reason why we don't throw an exception here is to simplify the code logic to
    keep type consistency when storing the value in buffers.
"""
RLBase.reward(env::MultiArmBanditsEnv) = env.reward

"""
Since `MultiArmBanditsEnv` is just a one-shot game, it doesn't matter what the
state is after each action. So here we can simply set it to a constant `1`.
"""
RLBase.state(env::MultiArmBanditsEnv) = 1

RLBase.state_space(env::MultiArmBanditsEnv) = Space(OneTo(1))

function RLBase.reset!(env::MultiArmBanditsEnv)
    env.is_terminated = false
    # since the reward is meaningless if the game is not started yet,
    # we don't need to reset the reward here.
end

"""
The multi-arm bandits environment is a stochastic environment. The resulted
reward may be different even after taking the same actions each time. So for
this kind of environments, the `Random.seed!(env)` must be implemented to help
increase reproducibility without creating a new instance of the same `rng`.
"""
Random.seed!(env::MultiArmBanditsEnv, x) = seed!(env.rng, x)

# For this simple one-shot environment, the default definitions are enough.
# Here we redefined them to help you compare the traits across different
# environments to gain a better understanding.

RLBase.NumAgentStyle(::MultiArmBanditsEnv) = SINGLE_AGENT
RLBase.DynamicStyle(::MultiArmBanditsEnv) = SEQUENTIAL
RLBase.ActionStyle(::MultiArmBanditsEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::MultiArmBanditsEnv) = IMPERFECT_INFORMATION  # the distribution of noise and original reward is unknown to the agent
RLBase.StateStyle(::MultiArmBanditsEnv) = Observation{Int}()
RLBase.RewardStyle(::MultiArmBanditsEnv) = TERMINAL_REWARD
RLBase.UtilityStyle(::MultiArmBanditsEnv) = GENERAL_SUM
RLBase.ChanceStyle(::MultiArmBanditsEnv) = STOCHASTIC  # the same action lead to different reward each time.


text:
"""
Check if MKL library `libmkl_rt.so` is available in `Libdl.dllist()`
(as is the case when loading [MKL.jl](https://github.com/JuliaLinearAlgebra/MKL.jl) or [MKL_jll.jl](https://github.com/JuliaBinaryWrappers/MKL_jll.jl)).
"""
mkl_is_loaded() = any(endswith(lib, "libmkl_rt.$(Libdl.dlext)") for lib in Libdl.dllist())

"Try to find `libmkl_rt.so` in `Libdl.dllist()`. Returns `nothing` if it can't be found."
function find_mkl()
    for lib in Libdl.dllist()
        if endswith(lib, "libmkl_rt.$(Libdl.dlext)")
            return lib
        end
    end
    return nothing
end

"""
    mkl_get_dynamic()
Wrapper around the MKL function [`mkl_get_dynamic`](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/support-functions/threading-control/mkl-get-dynamic.html).
"""
mkl_get_dynamic() = @ccall find_mkl().mkl_get_dynamic()::Cint

"""
    mkl_set_dynamic(flag::Integer)
Wrapper around the MKL function [`mkl_set_dynamic`](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/support-functions/threading-control/mkl-set-dynamic.html).
"""
mkl_set_dynamic(flag::Integer) = @ccall find_mkl().MKL_Set_Dynamic(flag::Cint)::Cvoid

"Returns the name of the loaded BLAS library (the first, if multiple are loaded)"
BLAS_lib() = basename(first(BLAS.get_config().loaded_libs).libname)


text:
module Registry

import ..Pkg, ..Types, ..API
using ..Pkg: depots1
using ..Types: RegistrySpec, Context

"""
    Pkg.Registry.add(url::String)
    Pkg.Registry.add(registry::RegistrySpec)

Add new package registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.add("General")
Pkg.Registry.add(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
Pkg.Registry.add(RegistrySpec(url = "https://github.com/JuliaRegistries/General.git"))
```
"""
function add end
add(reg::Union{String,RegistrySpec}) = add([reg])
add(regs::Vector{String}) = add([RegistrySpec(name = name) for name in regs])
add(regs::Vector{RegistrySpec}) = add(Context(), regs)
add(ctx::Context, regs::Vector{RegistrySpec}) =
    Types.clone_or_cp_registries(ctx, regs)

"""
    Pkg.Registry.rm(registry::String)
    Pkg.Registry.rm(registry::RegistrySpec)

Remove registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.rm("General")
Pkg.Registry.rm(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
function rm end
rm(reg::Union{String,RegistrySpec}) = rm([reg])
rm(regs::Vector{String}) = rm([RegistrySpec(name = name) for name in regs])
rm(regs::Vector{RegistrySpec}) = rm(Context(), regs)
rm(ctx::Context, regs::Vector{RegistrySpec}) = Types.remove_registries(ctx, regs)

"""
    Pkg.Registry.update()
    Pkg.Registry.update(registry::RegistrySpec)
    Pkg.Registry.update(registry::Vector{RegistrySpec})

Update registries. If no registries are given, update
all available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.update()
Pkg.Registry.update("General")
Pkg.Registry.update(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
function update end
update(reg::Union{String,RegistrySpec}) = update([reg])
update(regs::Vector{String}) = update([RegistrySpec(name = name) for name in regs])
update(regs::Vector{RegistrySpec} = Types.collect_registries(depots1())) =
    update(Context(), regs)
update(ctx::Context, regs::Vector{RegistrySpec} = Types.collect_registries(depots1())) =
    Types.update_registries(ctx, regs; force=true)

"""
    Pkg.Registry.status()

Display information about available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.status()
```
"""
function status()
    regs = Types.collect_registries()
    regs = unique(r -> r.uuid, regs) # Maybe not?
    Types.printpkgstyle(stdout, Symbol("Registry Status"), "")
    if isempty(regs)
        println("  (no registries found)")
    else
        for reg in regs
            printstyled(" [$(string(reg.uuid)[1:8])]"; color = :light_black)
            print(" $(reg.name)")
            reg.url === nothing || print(" ($(reg.url))")
            println()
        end
    end
end

end # module


text:
using Test
using FlightMechanicsSimulator
using FlightMechanicsUtils

#  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
#  and simulation: dynamics, controls design, and autonomous systems. John Wiley
#  & Sons. (page 193 table 3.6-2)
trim_test_data = [
#   TAS   thtl    AOA     DE       thtl_tol    AOA_tol    DE_tol
#   ft/s  unit    deg     deg      unit        deg        deg
    130   0.816   45.6    20.1     0.0005      0.05       0.15
    140   0.736   40.3   -1.36     0.001       0.05       0.05
    150   0.619   34.6    0.173    0.0005      0.05       0.05
    170   0.464   27.2    0.621    0.001       0.05       0.05
    640   0.23    0.742  -0.871    0.0005      0.015      0.0005
    800   0.378   -0.045 -0.943    0.0005      0.001      0.001
    200   0.287   19.7    0.723    0.0005      0.05       0.05
    260   0.148   11.6   -0.09     0.0005      0.05       0.05
    300   0.122   8.49   -0.591    0.0005      0.01       0.005
    350   0.107   5.87   -0.539    0.001       0.005      0.005
    400   0.108   4.16   -0.591    0.0005      0.005      0.005
    440   0.113   3.19   -0.671    0.0005      0.005      0.005
    500   0.137   2.14   -0.756    0.001       0.01       0.005
    540   0.16    1.63   -0.798    0.0005      0.005      0.005
    600   0.2     1.04   -0.846    0.0005      0.01       0.005
    700   0.282   0.382  -0.9      0.0005      0.001      0.0005
]

xcg = 0.35

for case in eachrow(trim_test_data)
    local x, controls, x_trim, controls_trim, x_dot_trim, outputs_trim, cost
    x = [
        case[1]*FT2M,  #-> vt (m/s)
        deg2rad(10.),  # -> alpha (rad)
        0.0,  # -> beta (rad)
        0.0,  # -> phi (rad)
        deg2rad(10.),  # -> theta (rad)
        0.0,  #  -> psi (rad)
        0.0,  # -> P (rad/s)
        0.0,  # -> Q (rad/s)
        0.0,  # -> R (rad/s)
        0.0,  # -> North (m)
        0.0,  # -> East (m)
        0.0,  # -> Altitude (m)
        50.0,  # -> Pow
    ]

    controls = [
        0.5,  # thtl
        0.0,  # elev
        0.0,  # ail
        0.0,  # rudder
    ]

    # TRIM
    dssd, controls_trim, outputs_trim, cost = trim(
        SixDOFAeroEuler(x),
        controls,
        F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
        F16StevensAtmosphere(x[12]),
        LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
        0.0,
        0.0,
    )

    x_trim = get_x(dssd)
    x_dot_trim = get_xdot(dssd)

    @test isapprox(cost, zeros(6), atol=1e-12)
    @test isapprox(controls_trim[1], case[2], atol=case[5])  # THTL
    @test isapprox(rad2deg(x_trim[2]), case[3], atol=case[6])  # AOA
    @test isapprox(controls_trim[2], case[4], atol=case[7])  # DE
end

#  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
#  and simulation: dynamics, controls design, and autonomous systems. John Wiley
#  & Sons. (page 195 table 3.6-3)

# NOMINAL (first column)
xcg = 0.35

x = [
    502*FT2M,  #-> vt (m/s)
    deg2rad(10.),  # -> alpha (rad)
    0.0,  # -> beta (rad)
    0.0,  # -> phi (rad)
    deg2rad(10.),  # -> theta (rad)
    0.0,  #  -> psi (rad)
    0.0,  # -> P (rad/s)
    0.0,  # -> Q (rad/s)
    0.0,  # -> R (rad/s)
    0.0,  # -> North (m)
    0.0,  # -> East (m)
    0.0,  # -> Altitude (m)
    50.0,  # -> Pow
]

controls = [
    0.5,  # thtl
    0.0,  # elev
    0.0,  # ail
    0.0,  # rudder
]

dssd, controls_trim, outputs_trim, cost = trim(
        SixDOFAeroEuler(x),
        controls,
        F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
        F16StevensAtmosphere(x[12]),
        LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
        0.0,
        0.0,
    )

x_trim = get_x(dssd)
x_dot_trim = get_xdot(dssd)

@test isapprox(cost, zeros(6), atol=1e-12)
@test isapprox(x_trim[2], 0.03691, atol=0.00005)  # AOA
@test isapprox(x_trim[3], -4e-9, atol=1e-8)  # AOS
@test isapprox(x_trim[4], 0)  # PHI
@test isapprox(x_trim[5], 0.03691, atol=0.00005)  # THETA
@test isapprox(x_trim[7], 0)  # P
@test isapprox(x_trim[8], 0)  # Q
@test isapprox(x_trim[9], 0)  # R
@test isapprox(controls_trim[1], 0.1385, atol=0.0001)  # THTL
@test isapprox(controls_trim[2], -0.7588, atol=0.0002)  # DE
@test isapprox(controls_trim[3], -1.2e-7, atol=1e-6)  # DA
@test isapprox(controls_trim[4], 6.2e-7, atol=1e-6)  # DR

# XCG = 0.3 (second column)
xcg = 0.3

x = [
    502*FT2M,  #-> vt (m/s)
    deg2rad(10.),  # -> alpha (rad)
    0.0,  # -> beta (rad)
    0.0,  # -> phi (rad)
    deg2rad(10.),  # -> theta (rad)
    0.0,  #  -> psi (rad)
    0.0,  # -> P (rad/s)
    0.0,  # -> Q (rad/s)
    0.0,  # -> R (rad/s)
    0.0,  # -> North (m)
    0.0,  # -> East (m)
    0.0,  # -> Altitude (m)
    50.0,  # -> Pow
]

controls = [
    0.5,  # thtl
    0.0,  # elev
    0.0,  # ail
    0.0,  # rudder
]

dssd, controls_trim, outputs_trim, cost = trim(
        SixDOFAeroEuler(x),
        controls,
        F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
        F16StevensAtmosphere(x[12]),
        LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
        0.0,
        0.0,
    )

x_trim = get_x(dssd)
x_dot_trim = get_xdot(dssd)

@test isapprox(cost, zeros(6), atol=1e-12)
@test isapprox(x_trim[2], 0.03936, atol=0.00005)  # AOA
@test isapprox(x_trim[3], 4.1e-9, atol=1e-8)  # AOS
@test isapprox(x_trim[4], 0)  # PHI
@test isapprox(x_trim[5], 0.03936, atol=0.00005)  # THETA
@test isapprox(x_trim[7], 0)  # P
@test isapprox(x_trim[8], 0)  # Q
@test isapprox(x_trim[9], 0)  # R
@test isapprox(controls_trim[1], 0.1485, atol=0.00005)  # THTL
@test isapprox(controls_trim[2], -1.931, atol=0.0001)  # DE
@test isapprox(controls_trim[3], -7e-8, atol=1e-6)  # DA
@test isapprox(controls_trim[4], 8.3e-7, atol=1e-6)  # DR

# XCG = 0.38 (third column)
xcg = 0.38

x = [
    502*FT2M,  #-> vt (m/s)
    deg2rad(10.),  # -> alpha (rad)
    0.0,  # -> beta (rad)
    0.0,  # -> phi (rad)
    deg2rad(10.),  # -> theta (rad)
    0.0,  #  -> psi (rad)
    0.0,  # -> P (rad/s)
    0.0,  # -> Q (rad/s)
    0.0,  # -> R (rad/s)
    0.0,  # -> North (m)
    0.0,  # -> East (m)
    0.0,  # -> Altitude (m)
    50.0,  # -> Pow
]

controls = [
    0.5,  # thtl
    0.0,  # elev
    0.0,  # ail
    0.0,  # rudder
]

dssd, controls_trim, outputs_trim, cost = trim(
        SixDOFAeroEuler(x),
        controls,
        F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
        F16StevensAtmosphere(x[12]),
        LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
        0.0,
        0.0,
    )

x_trim = get_x(dssd)
x_dot_trim = get_xdot(dssd)

@test isapprox(cost, zeros(6), atol=1e-12)
@test isapprox(x_trim[2], 0.03544, atol=0.00005)  # AOA
@test isapprox(x_trim[3], 3.1e-8, atol=1e-7)  # AOS
@test isapprox(x_trim[4], 0)  # PHI
@test isapprox(x_trim[5], 0.03544, atol=0.00005)  # THETA
@test isapprox(x_trim[7], 0)  # P
@test isapprox(x_trim[8], 0)  # Q
@test isapprox(x_trim[9], 0)  # R
@test isapprox(controls_trim[1], 0.1325, atol=0.0001)  # THTL
@test isapprox(controls_trim[2], -0.05590, atol=0.0005)  # DE
@test isapprox(controls_trim[3], -5.1e-8, atol=1e-6)  # DA
@test isapprox(controls_trim[4], 4.3e-6, atol=1e-5)  # DR

# Coordinated turn (fourth column)
xcg = 0.3

x = [
    502*FT2M,  #-> vt (m/s)
    deg2rad(10.),  # -> alpha (rad)
    0.0,  # -> beta (rad)
    0.0,  # -> phi (rad)
    deg2rad(10.),  # -> theta (rad)
    0.0,  #  -> psi (rad)
    0.0,  # -> P (rad/s)
    0.0,  # -> Q (rad/s)
    0.0,  # -> R (rad/s)
    0.0,  # -> North (m)
    0.0,  # -> East (m)
    0.0,  # -> Altitude (m)
    50.0,  # -> Pow
]

controls = [
    0.5,  # thtl
    0.0,  # elev
    0.0,  # ail
    0.0,  # rudder
]

dssd, controls_trim, outputs_trim, cost = trim(
        SixDOFAeroEuler(x),
        controls,
        F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
        F16StevensAtmosphere(x[12]),
        LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
        0.0,
        0.3,  # rad/s
    )

x_trim = get_x(dssd)
x_dot_trim = get_xdot(dssd)

@test isapprox(cost, zeros(6), atol=1e-12)
@test isapprox(x_trim[2], 0.2485, atol=0.0005)  # AOA
@test isapprox(x_trim[3], 4.8e-4, atol=0.00005)  # AOS
@test isapprox(x_trim[4], 1.367, atol=0.0005)  # PHI
@test isapprox(x_trim[5], 0.05185, atol=0.00005)  # THETA
@test isapprox(x_trim[7], -0.01555, atol=0.00001)  # P
@test isapprox(x_trim[8], 0.2934, atol=0.00005)  # Q
@test isapprox(x_trim[9], 0.06071, atol=0.000005)  # R
@test isapprox(controls_trim[1], 0.8499, atol=0.0005)  # THTL
@test isapprox(controls_trim[2], -6.256, atol=0.001)  # DE
@test isapprox(controls_trim[3], 0.09891, atol=0.00005)  # DA
@test isapprox(controls_trim[4], -0.4218, atol=0.0005)  # DR


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #
function _compute_crr_index_array(numberOfLevels::Int64)

    # ok, so lets build an index array -
    number_items_per_level = [(i+1) for i=0:numberOfLevels]
    tmp_array = Array{Int64,1}()
    theta = 0
    for value in number_items_per_level
        for index = 1:value
            push!(tmp_array,theta)
        end
        theta = theta + 1
    end

    N = sum(number_items_per_level[1:(numberOfLevels-1)])
    index_array = Array{Int64,2}(undef,N,3)
    for row_index = 1:N

        #index_array[row_index,1] = tmp_array[row_index]
        index_array[row_index,1] = row_index
        index_array[row_index,2] = row_index + 1 + tmp_array[row_index]
        index_array[row_index,3] = row_index + 2 + tmp_array[row_index]
    end

    return index_array
end

function _compute_full_binomial_index_array(numberOfLevels::Int64)

    number_of_elements = (2^(Int(numberOfLevels))) - 1
    index_array = Array{Int64,2}(undef,number_of_elements,3)

    # populate the index array -
    for index = 1:number_of_elements

        # compute the child indexs -
        left_child_index = 2*(index - 1) + 2
        right_child_index = 2*(index - 1) + 3

        # grab and go -
        index_array[index,1] = index
        index_array[index,2] = left_child_index
        index_array[index,3] = right_child_index
    end

    # return -
    return index_array
end

# array based -
function _build_binary_lattice_intrinsic_value_array(contractSet::Set{PSAbstractAsset}, binaryPriceArray::Array{Float64,1})::PSResult

    # initialize -
    intrinsic_value_array = similar(binaryPriceArray)

    # ok, so lets go through each price, compute the intrinsic value - and then store
    for (index,price) in enumerate(binaryPriceArray)

        # ok, get the underlying price -
        underlying_price_value = binaryPriceArray[index]

        # compute the intrinsic value -
        result = intrinsic_value(contractSet,underlying_price_value)
        if (isa(result.value,Exception) == true)
            return result
        end
        iv_value = result.value

        # add the iv to the array - index should work out
        intrinsic_value_array[index] = iv_value
    end

    # return -
    return PSResult{Array{Float64,1}}(intrinsic_value_array)
end

function _build_binary_lattice_underlying_price_array(basePrice::Float64, volatility::Float64, timeToExercise::Float64;
    numberOfLevels::Int64 = 100)::PSResult

    # compute up and down perturbations -
    Δt = (timeToExercise/numberOfLevels)
    U = exp(volatility * √Δt)
    D = 1 / U

    # compute the index array =
    index_array = _compute_crr_index_array(numberOfLevels);
    max_element = index_array[end,end]
    N = index_array[end,1]

    # compute price array -
    priceArray = zeros(max_element)
    priceArray[1] = basePrice
    for row_index = 1:N

        # parent index -
        parent_index = index_array[row_index,1]
        left_child_index = index_array[row_index,2]
        right_child_index = index_array[row_index,3]

        # get the basePrice -
        basePrice = priceArray[parent_index]

        # compute the prices -
        down_price = basePrice*D
        priceArray[right_child_index] = down_price

        up_price = basePrice*U
        priceArray[left_child_index] = up_price
    end

    # return the price array -
    return PSResult{Array{Float64,1}}(priceArray)
end

function _build_binary_lattice_option_value_array(intrinsicValueArray::Array{Float64,1}, latticeModel::PSBinaryLatticeModel;
    earlyExcercise::Bool = true, numberOfLevels::Int64 = 100)::PSResult

    # initialize -
    contract_price_array = copy(intrinsicValueArray)

    # get stuff from the lattice model -
    volatility = latticeModel.volatility
    timeToExercise = latticeModel.timeToExercise
    riskFreeRate = latticeModel.riskFreeRate
    dividendRate = latticeModel.dividendRate

    # compute up and down perturbations -
    Δt = (timeToExercise/numberOfLevels)
    U = exp(volatility * √Δt)
    D = 1 / U
    p = (exp((riskFreeRate - dividendRate)*Δt) - D)/(U - D)
    DF = exp(-riskFreeRate*Δt)

    # create a index table -
    index_table = _compute_crr_index_array(numberOfLevels)
    N = index_table[end,1]

    # ok, so now lets compute the value for the nodes -
    for compute_index = 1:N

        # get the indexs -
        parent_node_index = index_table[compute_index,1]
        child_left_index = index_table[compute_index,2]
        child_right_index = index_table[compute_index,3]

        # compute -
        contract_price = DF*(p*contract_price_array[child_left_index]+(1-p)*contract_price_array[child_right_index])
        if (earlyExcercise == false)
            contract_price_array[parent_node_index] = contract_price
        else
            excercise_value = contract_price_array[parent_node_index]
            contract_price_array[parent_node_index] = max(excercise_value,contract_price)
        end
    end

    # setup the results tuple -
    results_tuple = (option_index_table=index_table, option_contract_price_array=contract_price_array, U=U, D=D, PUP=p, PDOWN=(1-p), DF=DF)

    # return -
    return PSResult(results_tuple)
end
# ----------------------------------------------------------------------------------------------------------- #

# --- PUBLIC METHODS ---------------------------------------------------------------------------------------- #

function compute_underlying_price_distribution(timeStepIndex::Int64,latticeModel::PSBinaryLatticeModel, baseUnderlyingPrice::Float64,
    movementFunction::Function)::PooksoftBase.PSResult

    # initialize -
    priceDistributionArray = Array{Float64,2}(undef,(timeStepIndex+1), 4)
    riskFreeRate = latticeModel.riskFreeRate
    kIndexVector = range(0,stop=timeStepIndex,step=1) |> collect

    # compute the up and downmove values -
    (u,d) = movementFunction(latticeModel)

    # compute the probability p of an *up* move -
    p = ((1 + riskFreeRate) - d)/(u - d)

    # main loop -
    for (index,k) in enumerate(kIndexVector)

        # compute the price -
        price = baseUnderlyingPrice*(u^k)*(d^(timeStepIndex - k))

        # compute the probability -
        prob = binomial(timeStepIndex,k)*(p^k)*((1 - p)^(timeStepIndex-k))

        # package -
        priceDistributionArray[index,1] = timeStepIndex
        priceDistributionArray[index,2] = k
        priceDistributionArray[index,3] = prob
        priceDistributionArray[index,4] = price
    end

    # return -
    return PSResult{Array{Float64,2}}(priceDistributionArray)
end

"""
    option_contract_price(contractSet::Set{PSAbstractAsset}, latticeModel::PSBinaryLatticeModel, baseUnderlyingPrice::Float64;
        earlyExercise::Bool = false)::PooksoftBase.PSResult

Estimate the price of a contract using a binary lattice pricing model.
"""
function option_contract_price(contractSet::Set{PSAbstractAsset}, latticeModel::PSBinaryLatticeModel, baseUnderlyingPrice::Float64;
    earlyExercise::Bool = true)::PooksoftBase.PSResult

    # initialize -
    option_contract_price = 0.0

    # we need to get some stuff from the lattice model -
    volatility = latticeModel.volatility
    timeToExercise = latticeModel.timeToExercise
    numberOfLevels = latticeModel.numberOfLevels

    # compute the price array -
    result = _build_binary_lattice_underlying_price_array(baseUnderlyingPrice, volatility, timeToExercise;
        numberOfLevels = numberOfLevels)
    if (isa(result.value,Exception) == true)
        return result
    end
    lattice_price_array = result.value

    # compute the intrinsic value array -
    result = _build_binary_lattice_intrinsic_value_array(contractSet,lattice_price_array)
    if (isa(result.value,Exception) == true)
        return result
    end
    iv_array = result.value

    # ok, let's build the option value array -
    result = _build_binary_lattice_option_value_array(iv_array, latticeModel;
        earlyExcercise = earlyExercise, numberOfLevels = numberOfLevels)
    if (isa(result.value,Exception) == true)
        return result
    end
    cost_calc_tuple = result.value

    # create a named tuple and return the results -
    results_tuple = (lattice_price_array=lattice_price_array, intrinsic_value_array = iv_array, cost_calculation_result=cost_calc_tuple)

    # return -
    return PSResult(results_tuple)
end

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


text:
"""
    AmplitudeConstrainedSet

Abstract set of complex-values `x` with a given absolute value `|x| = amp`. The amplitude can be extracted by function `amp`.

"""
abstract type  AmplitudeConstrainedSet <: FeasibleSet end
export AmplitudeConstrainedSet

amp(s::AmplitudeConstrainedSet) = s.amp #|| error("amplitude is not defined for $(typeof(s))")
getelement(s::AmplitudeConstrainedSet) = complex(float(amp(s))) # element of the set is a complex array

"""
    ConstrainedByAmplitude{T,N}

Set of abstract arryas with element type `T` and dimensions `N` defined by the amplitude constraint `|x| = amp`.

"""
struct ConstrainedByAmplitude{T,N} <: AmplitudeConstrainedSet where {T <: Real, N}
    amp::Array{T,N}  #
end

amp(s::ConstrainedByAmplitude) = s.amp

"""
ConstrainedByAmplitude(a::AbstractArray{T,N})

Construct set defined by the amplitude constraint `|x| = a`. Type and dimension of the set are inhereited from the array.
"""
function ConstrainedByAmplitude(a::AbstractArray{T,N}) where {T,N}
    ConstrainedByAmplitude{T,N}(a)
end

function ConstrainedByAmplitude(a::AbstractArray{Union{T,Nothing},N}) where {T,N}
    ConstrainedByAmplitude{Union{T,Nothing},N}(a)
end

# function ConstrainedByAmplitude(a::AbstractArray{T,N})
#     ConstrainedByAmplitude{T,N}(a)
# end
export ConstrainedByAmplitude

"""
    ConstrainedByAmplitudeMasked(a, mask::Vector{CartesianIndex{N}})
    ConstrainedByAmplitudeMasked(a, AbstractArray{Bool})

The amplitude constrained set only in the indexes given by mask:  `|xᵢ| = aᵢ` for `i ∈ mask`.
"""
struct ConstrainedByAmplitudeMasked{T,N} <: AmplitudeConstrainedSet where {T <: Real, N}
    amp::Array{T,N}  #
    mask::Vector{CartesianIndex{N}}  #
end

ConstrainedByAmplitudeMasked(a, m::AbstractArray{Bool}) = ConstrainedByAmplitudeMasked(a, findall(m))

export ConstrainedByAmplitudeMasked
amp(s::ConstrainedByAmplitudeMasked) = s.amp

# here the backward plan is in place
FourierTransformedSet(s::AmplitudeConstrainedSet) =
    FourierTransformedSet(s, FFTW.plan_fft(getelement(s)), FFTW.plan_ifft!(getelement(s)))

function project!(xp, x, feasset::ConstrainedByAmplitude)
    # @inbounds for i in eachindex(xp)
    #     xp[i] = feasset.amp[i] * exp( im * angle(x[i]))
    # end
    # xp .= feasset.amp .* exp.(1.0im .* angle.(x)) # this is 1000 times slower
    # xp .= _replace_amp.(feasset.amp, x)

    # xp .= feasset.amp .* _unit_amp.(x)
    # xp .= update_amplitude.(feasset.amp,x)

    @inbounds for i in eachindex(xp)
        xp[i] = update_amplitude(feasset.amp[i], x[i])
    end

    return xp
end

function project!(x, feasset::ConstrainedByAmplitude)
    for i in eachindex(x)
        x[i] = update_amplitude(feasset.amp[i], x[i])
    end

    return x
end

function project!(xp, x, feasset::ConstrainedByAmplitudeMasked)
    @inbounds for i in feasset.mask
        xp[i] = update_amplitude(feasset.amp[i], x[i])
    end
    return xp
end

function project!(x, feasset::ConstrainedByAmplitudeMasked)
    @inbounds for i in feasset.mask
        x[i] = update_amplitude(feasset.amp[i], x[i])
    end
    return x
end

update_amplitude(amp, x) = isnothing(amp) ? x : amp * _unit_amp(x)

update_amplitude!(xp, amp, x) = isnothing(amp) ? x : amp * _unit_amp(x)

@inline function _unit_amp(z)
    abs(z)≈0 ? one(z) : z/abs(z)
end

@inline function _replace_amp(amp,z)
    abs(z)≈0 ? zero(z) : z* (amp/abs(z))
end

size(feasset::ConstrainedByAmplitude) = size(feasset.amp)

include("ShapeConstraint.jl")


text:
using PowerModels
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)
    model = PP2PM.get_model(pm["pm_model"])

    solver = PP2PM.get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    report_duals = get(pm, "report_duals", false)  # true or false
    branch_limits = get(pm, "branch_limits", "hard")  # "hard", "soft", "none"
    objective = get(pm, "objective", "default")  # TODO: check best options

    settings = Dict("output" => Dict("branch_flows" => true,
                                     "duals" => report_duals))

    result = PowerModels.run_opf(pm, model, solver, setting=settings)
    return result
end

max_stars_repo_name: venidera/pandapower

max_issues_repo_name: venidera/pandapower

max_forks_repo_name: venidera/pandapower