text:
# Put in a separate page so it can be used by SciMLDocs.jl

pages=[
    "Home" => "index.md",
    "Example: Nice DiffEq Syntax Without A DSL" => "Example_dsl.md",
    "SLArrays" => "SLArrays.md",
    "LArrays" => "LArrays.md",
    "Relation to NamedTuples" => "NamedTuples_relation.md",
    "Note_labelled_slices.md"
]


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example

import ..InParentNamespace
import FlatBuffers

FlatBuffers.@with_kw mutable struct Monster{T}
#=
#  an example documentation comment: monster object
=#
    pos::Union{Vec3, Nothing} = nothing
    mana::Int16 = 150
    hp::Int16 = 100
    name::String = ""
    inventory::Vector{UInt8} = []
    color::Color = 8
    test_type::UInt8 = 0
    test::T = nothing
    test4::Vector{Test} = []
    testarrayofstring::Vector{String} = []
#=
#  an example documentation comment: this will end up in the generated code
#  multiline too
=#
    testarrayoftables::Vector{Monster{T}} = []
    enemy::Union{Monster{T}, Nothing} = nothing
    testnestedflatbuffer::Vector{UInt8} = []
    testempty::Union{Stat, Nothing} = nothing
    testbool::Bool = false
    testhashs32_fnv1::Int32 = 0
    testhashu32_fnv1::UInt32 = 0
    testhashs64_fnv1::Int64 = 0
    testhashu64_fnv1::UInt64 = 0
    testhashs32_fnv1a::Int32 = 0
    testhashu32_fnv1a::UInt32 = 0
    testhashs64_fnv1a::Int64 = 0
    testhashu64_fnv1a::UInt64 = 0
    testarrayofbools::Vector{Bool} = []
    testf::Float32 = 3.14159
    testf2::Float32 = 3.0
    testf3::Float32 = 0.0
    testarrayofstring2::Vector{String} = []
    testarrayofsortedstruct::Vector{Ability} = []
    flex::Vector{UInt8} = []
    test5::Vector{Test} = []
    vector_of_longs::Vector{Int64} = []
    vector_of_doubles::Vector{Float64} = []
    parent_namespace_test::Union{InParentNamespace, Nothing} = nothing
    vector_of_referrables::Vector{Referrable} = []
    single_weak_reference::UInt64 = 0
    vector_of_weak_references::Vector{UInt64} = []
    vector_of_strong_referrables::Vector{Referrable} = []
    co_owning_reference::UInt64 = 0
    vector_of_co_owning_references::Vector{UInt64} = []
    non_owning_reference::UInt64 = 0
    vector_of_non_owning_references::Vector{UInt64} = []
    #= any_unique_type::UInt8 = 0 =#
    #= any_unique::AnyUniqueAliases = nothing =#
    #= any_ambiguous_type::UInt8 = 0 =#
    #= any_ambiguous::AnyAmbiguousAliases = nothing =#
    vector_of_enums::Vector{Color} = []
end
FlatBuffers.@ALIGN(Monster, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:Monster} = [
    0x00000004, 0x00000006, 0x00000008, 0x0000000A,
    0x0000000E, 0x00000010, 0x00000012, 0x00000014,
    0x00000016, 0x00000018, 0x0000001A, 0x0000001C,
    0x0000001E, 0x00000020, 0x00000022, 0x00000024,
    0x00000026, 0x00000028, 0x0000002A, 0x0000002C,
    0x0000002E, 0x00000030, 0x00000032, 0x00000034,
    0x00000036, 0x00000038, 0x0000003A, 0x0000003C,
    0x0000003E, 0x00000040, 0x00000042, 0x00000044,
    0x00000046, 0x00000048, 0x0000004A, 0x0000004C,
    0x0000004E, 0x00000050, 0x00000052, 0x00000054,
    0x00000056, 0x00000058, 0x0000005A, 0x0000005C,
    0x0000005E, 0x00000060, 0x00000062
]
FlatBuffers.root_type(::Type{T}) where {T<:Monster} = true
FlatBuffers.file_identifier(::Type{T}) where {T<:Monster} = "MONS"
FlatBuffers.file_extension(::Type{T}) where {T<:Monster} = "mon"

Monster{T}(buf::AbstractVector{UInt8}) where {T} = FlatBuffers.read(Monster{T}, buf)
Monster{T}(io::IO) where {T} = FlatBuffers.deserialize(io, Monster{T})

ext: jl


text:
"""
    function read_pdb(file_pdb)

This will read a .pdb file, getting all the coordinate information.
- `file_pdb`: path of target .pdb file.

# Output
- `id_vec`
- ``

# Example
```julia-repl
file_pdb = "gb1.pdb"
coord = read_pdb(file_pdb)
```
"""
function read_pdb(file_pdb)

    # Setting variabls
    col_range_id = 7:11
    col_range_atom_name = 13:16
    col_range_mol = 23:26
    col_range_x = 31:38
    col_range_y = 39:46
    col_range_z = 47:54

    # Reading file
    io = open(file_pdb)
    str = read(io, String)
    reg = r"ATOM..*[A-Z]"
    res = findall(reg, str)
    res = [str[n] for n in res]

    # Ansys
    num_atoms = length(res)
    atom_vec = [String(strip(res[atom][col_range_atom_name], ' ')) for atom=1:num_atoms]
    mol_vec = zeros(Int, num_atoms, 1)
    id_vec = zeros(Int, num_atoms, 1)
    coord_mat = zeros(num_atoms, 3)

    for atom = 1:num_atoms
        id_vec[atom, 1] = parse(Int, res[atom][col_range_id])
        mol_vec[atom, 1] = parse(Int, res[atom][col_range_mol])
        coord_mat[atom, 1] = parse(Float64, res[atom][col_range_x])
        coord_mat[atom, 2] = parse(Float64, res[atom][col_range_y])
        coord_mat[atom, 3] = parse(Float64, res[atom][col_range_z])
    end

    return id_vec, atom_vec, mol_vec, coord_mat
end

"""
    function read_pdb()

This will read a .pdb file, getting all the coordinate information.
- `file_pdb`: path of target .pdb file.

# Example
```julia-repl
file_pdb = "gb1.pdb"
coord = read_pdb(file_pdb)
```
"""
function read_psf(file_psf)
    # Read psf file
    io = open(file_psf)
    raw_info = split(read(io, String), "\n")
    raw_start = findall(x->occursin("!N", x), raw_info)[2:7] .+ 1 # Ignore !NTITLE
    raw_end = raw_start .- 3

    # Atom info
    col_type = 6
    col_charge = 7
    num_atoms = raw_end[2] - raw_start[1] + 1
    type_vec = [String(split(raw_info[i])[col_type]) for i=raw_start[1]:raw_end[2]]
    charge_vec = [parse(Float64, String(split(raw_info[i])[col_charge])) for i=raw_start[1]:raw_end[2]]

    # Bond info
    bond_topo = [strvec2intvec(split(raw_info[i])) for i=raw_start[2]:raw_end[3]]
    bond_topo = Array(reshape(vcat(bond_topo...), (2, :))')

    # Angle info
    angle_topo = [strvec2intvec(split(raw_info[i])) for i=raw_start[3]:raw_end[4]]
    angle_topo = Array(reshape(vcat(angle_topo...), (3, :))')

    # Dihedral info
    dihedral_topo = [strvec2intvec(split(raw_info[i])) for i=raw_start[4]:raw_end[5]]
    dihedral_topo = Array(reshape(vcat(dihedral_topo...), (4, :))')

    # Improper info
    improper_topo = [strvec2intvec(split(raw_info[i])) for i=raw_start[5]:raw_end[6]]
    improper_topo = Array(reshape(vcat(improper_topo...), (4, :))')

    return type_vec, charge_vec, bond_topo, angle_topo, dihedral_topo, improper_topo
end

function assign_para_atom(id_vec, mol_vec, type_vec, charge_vec, coord_mat, potential::T) where T <:Potential
    num_atoms = size(id_vec, 1)
    type_list = unique(type_vec)
    mol_vec = mol_vec .- minimum(mol_vec) .+ 1  # Mol id start from 1
    num_atom_types = size(type_list, 1)

    vec_atom = Vector{Atom}(undef, num_atoms)

    for atom = 1 : num_atoms
        vec_atom[atom] = Atom(id_vec[atom], mol_vec[atom], findfirst(x->x==type_vec[atom], type_list),
                        charge_vec[atom], coord_mat[atom, :])
    end

    para_pair = [findfirst(x->judgeequal(x, [type_list[i] for i in id_vec[atom]]),
                             potential.para_pair) for atom = 1:num_atom_types]
    para_pair = vcat([pair.para for pair in potential.para_pair[para_pair]]...)

    para_pair = reshape(para_pair, (4, :))'

    para_mass = [findfirst(x->judgeequal(x, [type_list[i] for i in id_vec[atom]]),
    potential.para_mass) for atom = 1:num_atom_types]
    para_mass = vcat([pair.para for pair in potential.para_mass[para_mass]]...)

    return num_atoms, num_atom_types, vec_atom, para_pair, para_mass, type_list
end

function assign_para_bond(type_vec, bond_topo, potential)
    num_bonds = size(bond_topo, 1)
    num_bond_types = 0
    vec_bond = Vector{Bond}(undef, num_bonds)
    type_list = [0]
    for id = 1:num_bonds
        para_type_id = findfirst(x->judgeequal(x, [type_vec[i] for i in bond_topo[id, :]]), potential.para_bond)
        if para_type_id == nothing
            print([type_vec[i] for i in bond_topo[id, :]])
            error("No matching bond type is found in Charmm36 force filed.")
        end
        if !in(para_type_id, type_list)
            num_bond_types += 1
            append!(type_list, para_type_id)
        end
        vec_bond[id] = Bond(id, findfirst(x->x==para_type_id, type_list) - 1, bond_topo[id, :])
    end
    para_bond = [bond.para for bond in potential.para_bond[type_list[2:end]]]
    para_bond = Array(hcat(para_bond...)')
    return num_bonds, num_bond_types, vec_bond, para_bond
end

function assign_para_angle(type_vec, angle_topo, potential)
    num_angles = size(angle_topo, 1)
    num_angle_types = 0
    vec_angle = Vector{Angle}(undef, num_angles)
    type_list = [0]
    for id = 1:num_angles
        para_type_id = findfirst(x->judgeequal(x, [type_vec[i] for i in angle_topo[id, :]]), potential.para_angle)
        if para_type_id == nothing
            print([type_vec[i] for i in angle_topo[id, :]])
            error("No matching angle type is found in Charmm36 force filed.")
        end
        if !in(para_type_id, type_list)
            num_angle_types += 1
            append!(type_list, para_type_id)
        end
        vec_angle[id] = Angle(id, findfirst(x->x==para_type_id, type_list)-1, angle_topo[id, :])
    end
    para_angle = [angle.para for angle in potential.para_angle[type_list[2:end]]]
    para_angle = Array(hcat(para_angle...)')
    return num_angles, num_angle_types, vec_angle, para_angle
end

function assign_para_dihedral(type_vec, dihedral_topo, potential)
    num_dihedrals = size(dihedral_topo, 1)
    num_dihedral_types = 0
    vec_dihedral = Vector{Dihedral}(undef, num_dihedrals)
    type_list = [0]
    for id = 1:num_dihedrals
        para_type_id = findfirst(x->judgeequal(x, [type_vec[i] for i in dihedral_topo[id, :]]), potential.para_dihedral)
        if para_type_id == nothing
            print([type_vec[i] for i in dihedral_topo[id, :]])
            error("No matching dihedral type is found in Charmm36 force filed.")
        end
        if !in(para_type_id, type_list)
            num_dihedral_types += 1
            append!(type_list, para_type_id)
        end
        vec_dihedral[id] = Dihedral(id, findfirst(x->x==para_type_id, type_list)-1, dihedral_topo[id, :])
    end
    para_dihedral = [dihedral.para for dihedral in potential.para_dihedral[type_list[2:end]]]
    para_dihedral = Array(hcat(para_dihedral...)')
    return num_dihedrals, num_dihedral_types, vec_dihedral, para_dihedral
end

function assign_para_improper(type_vec, improper_topo, potential)
    num_impropers = size(improper_topo, 1)
    num_improper_types = 0
    vec_improper = Vector{Improper}(undef, num_impropers)
    type_list = [0]
    for id = 1:num_impropers
        para_type_id = findfirst(x->judgeequal(x, [type_vec[i] for i in improper_topo[id, :]]), potential.para_improper)
        if para_type_id == nothing
            print([type_vec[i] for i in improper_topo[id, :]])
            error("No matching improper type shwon above is found in Charmm36 force filed.")
        end
        if !in(para_type_id, type_list)
            num_improper_types += 1
            append!(type_list, para_type_id)
        end
        vec_improper[id] = Improper(id, findfirst(x->x==para_type_id, type_list)-1, improper_topo[id, :])
    end
    para_improper = [improper.para for improper in potential.para_improper[type_list[2:end]]]
    para_improper = Array(hcat(para_improper...)')
    return num_impropers, num_improper_types, vec_improper, para_improper
end

"""
    function convert_vmd(file_pdb::String, file_psf::String, potential::Potential)
This will convert a set of files which created by VMD (.pdb and .psf) file into a `Data` variable that can be used to write a .data file.
- `file_pdb`: String of path of .pdb file.
- `file_psf`: String of path of .psf file.
- `poetntial`: Subclass of `Potential` which contains all information of forcefiled to automatically generate parameters, e.g. pair_coeff, bond_coeff etc.
"""
function converter_vmd(file_pdb::String, file_psf::String, potential::Potential, box_padding=2)
    # Reading Input
    id_vec, atom_vec, mol_vec, coord_mat = read_pdb(file_pdb)
    type_vec, charge_vec, bond_topo, angle_topo, dihedral_topo, improper_topo = read_psf(file_psf)

    box_size = Array(cat(minimum(coord_mat, dims=1), maximum(coord_mat, dims=1), dims=1)')
    #box_size[:, 1] .-= box_padding / 2
    #box_size[:, 2] .+= box_padding / 2
    box_tilt = [0, 0, 0]

    # Assign parameters
    num_atoms, num_atom_types, vec_atom, para_pair, para_mass, type_list = assign_para_atom(id_vec,
            mol_vec, type_vec, charge_vec, coord_mat, potential_charmm36)
    num_bonds, num_bond_types, vec_bond, para_bond = assign_para_bond(type_vec, bond_topo, potential)
    num_angles, num_angle_types, vec_angle, para_angle = assign_para_angle(type_vec, angle_topo, potential)
    num_dihedrals, num_dihedral_types, vec_dihedral, para_dihedral = assign_para_dihedral(type_vec, dihedral_topo, potential)
    num_impropers, num_improper_types, vec_improper, para_improper = assign_para_improper(type_vec, improper_topo, potential)

    # Judgement
    if num_bonds == 0
        vec_bond = 0
    end
    if num_angles == 0
        vec_angle = 0
    end
    if num_dihedrals == 0
        vec_dihedral = 0
    end
    if num_impropers == 0
        vec_improper = 0
    end

    # Creat Data instance
    data_basic = Data_Basic(num_atoms, num_bonds, num_angles, num_dihedrals, num_impropers,
    num_atom_types, num_bond_types, num_angle_types, num_dihedral_types, num_improper_types, box_size, box_tilt)
    data_str = Structure_VMD(type_list, num_atom_types, num_bond_types, num_angle_types, num_dihedral_types, num_improper_types,
        para_mass, para_pair, para_bond, para_angle, para_dihedral, para_improper)
    data = Data_Unit(data_basic, data_str,
        vec_atom, vec_bond, vec_angle, vec_dihedral, vec_improper)
    return data
end


text:
# Autogenerated wrapper script for github_release_jll for i686-linux-musl
export github_release

JLLWrappers.@generate_wrapper_header("github_release")
JLLWrappers.@declare_executable_product(github_release)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        github_release,
        "bin/github-release",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, AverageShiftedHistograms

makedocs(
    format = Documenter.HTML(),
    sitename = "AverageShiftedHistograms.jl",
    authors = "Josh Day",
    clean = true,
    pages = [
        "index.md",
        "api.md"
    ]
)

deploydocs(
    repo   = "github.com/joshday/AverageShiftedHistograms.jl.git",
)


text:
using COSMO
using SketchyCGAL
using LinearAlgebra, SparseArrays
using BSON

include("../utils.jl")

blas_threads = BLAS.get_num_threads()
# Run experiments single threaded
blas_threads = BLAS.set_num_threads(1)

## Data Load
G = graph_from_file(joinpath(dirname(@__DIR__), "data/gset/G72"))
# n = size(G, 1)
# C = -0.25*(Diagonal(G*ones(n)) - G)

## COSMO (try to get true opt)
cosmo_tt = @timed solve_with_COSMO(C; settings = COSMO.Settings(
    verbose = true,
    eps_abs = 1e-5,
    eps_rel = 1e-5,
    decompose = true,
    merge_strategy = COSMO.CliqueGraphMerge,
    max_iter = 5_000,
    # rho = 1e-5,
    # alpha = 1.6,
    kkt_solver = COSMO.QdldlKKTSolver,
    # kkt_solver = with_options(COSMO.MKLPardisoKKTSolver, msg_level_on = true),
    verbose_timing = true,
))
cosmo_results = tt.value
cosmo_time = tt.time
cosmo_alloc = tt.bytes/1e6 #(MB)

## SCGAL (run for awhile to get "true" opt)
## ---------- Problem setup ----------
function setup_maxcut_scgal(G::SparseMatrixCSC)
    # Parameters
    n = size(G, 1)
    d = n

    # Data
    # objective = min -1/4( ∑_ij w_ij*(1-yi*yj) ) = -1/4⟨diag(∑ᵢ w_ij) - W), Y⟩
    # note that this reformulation uses the fact that Tr(Y) = 1
    C = -0.25*(Diagonal(G*ones(n)) - G)
    b = ones(n)

    # Scaling variables -- so trace is bounded by 1
    scale_C = 1 / norm(C)
    scale_X = 1 / n
    # b .= b .* scale_X
    # const C_const = C * scale_C

    # Linear map
    # AX = diag(X)
    function A!(y, X)
        n = size(X, 1)
        for i in 1:n
            y[i] = X[i,i]
        end
        return nothing
    end

    # Adjoint
    # A*z = Diagonal(z)
    function A_adj!(S::SparseMatrixCSC, z)
        for (i, j, v) ∈ zip(findnz(S)...)
            S[i,j] = 0.0
        end
        S[diagind(S)] .= z
        return nothing
     end

    # primative 1: u -> C*u
    function C_u!(v, u)
        mul!(v, C_const, u)
    end

    #primative 2: (u, z) -> (A*z)u
    function Aadj_zu!(v, u, z)
        v .= u .* z
    end

    #primative 3: u -> A(uu^T)
    function A_uu!(z, u)
        z .= u.*u
    end

    return C, b, A!, A_adj!, n, d, scale_X, scale_C
end

## Long Run
function maxcut_scgal_long(G)
    C, b, A!, A_adj!, n, d, scale_X, scale_C = setup_maxcut_scgal(G)
    R = 10
    ηt(t) = 2.0/(t + 1.0)
    δt(t) = 1.0
    tt = @timed scgal_full(
        C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
        max_iters=100_000,
        print_iter=1_000,
        R=R,
        # logging=true,
        # logging_primal=true,
        ηt=ηt,
        δt=δt
    )

    # soln = tt.value
    # trial_time = tt.time
    # trial_alloc = tt.bytes/1e6 #(MB)

    return tt
end

BLAS.set_num_threads(1)
long_run = maxcut_scgal_long(G)
BSON.bson(joinpath(@__DIR__, "output/long_run.bson"), Dict("data" => long_run))
# 100000  -7.743749e+03   5.253496e-01   7.271706e-05      1067.732

## SCGAL Trials for R, Fig 7.1 Recreation
Rs = [10; 25; 50; 100]
function run_trials_R(G, Rs)
    trial_data = Dict("time" => Dict(), "iter" => Dict())
    C, b, A!, A_adj!, n, d, scale_X, scale_C = setup_maxcut_scgal(G)
    ηt(t) = 2.0/(t + 1.0)
    δt(t) = 1.0

    for R in Rs
        tt = @timed scgal_full(
            C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
            max_iters=10_000,
            print_iter=1_000,
            R=R,
            logging=true,
            # logging_primal=true,
            ηt=ηt,
            δt=δt
        )
        trial_data["time"][R] = tt

        tt2 = @timed scgal_full(
            C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
            max_iters=10_000,
            print_iter=1_000,
            R=R,
            logging=true,
            logging_primal=true,
            ηt=ηt,
            δt=δt
        )
        trial_data["iter"][R] = tt2

    end
    return trial_data
end

trial_data = run_trials_R(G, Rs)
BSON.bson(joinpath(@__DIR__, "output/trials_R_G72.bson"), trial_data)

## SCGAL trials for weight (fix R = 10)
function run_trials_weights(G; R=10, max_iters=10_000, print_iter=1000, rseed=0, δs=[])
    trial_data = Dict()
    C, b, A!, A_adj!, n, d, scale_X, scale_C = setup_maxcut_scgal(G)

    tt = @timed scgal_full(
        C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
        max_iters=max_iters,
        print_iter=print_iter,
        R=R,
        logging=true,
        # logging_primal=true,
        ηt=t->2.0/(t + 1.0),
        δt=t->1.0,
        rseed=rseed
    )
    trial_data["std"] = tt

    tt_sa = @timed scgal_full(
        C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
        max_iters=max_iters,
        print_iter=print_iter,
        R=R,
        logging=true,
        # logging_primal=true,
        ηt=t->2.0/(t + 2.0),
        δt=t->1.0t/(1.0t + 1),
        rseed=rseed
    )
    trial_data["sa"] = tt_sa

    for δ in δs
        tt_ea = @timed scgal_full(
            C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
            max_iters=max_iters,
            print_iter=print_iter,
            R=R,
            logging=true,
            # logging_primal=true,
            ηt=t->2.0/(t + 3.0),
            δt=t->δ,
            rseed=rseed
        )
        trial_data["ea_$δ"] = tt_ea
    end

    return trial_data
end

δs = [0.01, 0.05, 0.1, 0.3, 0.5, 0.8, 0.9, 0.95, 0.99]
trial_data_weights = run_trials_weights(G; max_iters=10_000, print_iter=1_000, rseed=0, δs=δs)
BSON.bson(joinpath(@__DIR__, "output/trials_weights_G72.bson"), trial_data_weights)

# ##
# function condition_check(δ, c0)
#     return (1 - (1-δ)*(2+1)^2/2^2)*δ*c0^2
# end
#
# condition_check.(δs, 4)


text:
#!/usr/bin/env julia --project=.

using JuliaEchoServer

port = get(ENV, "PORT", "8000")

server = JuliaEchoServer.EchoServer(parse(Int, port))
JuliaEchoServer.listen(server)

max_stars_repo_name: nmatsui/JuliaEchoServer

max_issues_repo_name: nmatsui/JuliaEchoServer

max_forks_repo_name: nmatsui/JuliaEchoServer


text:
"""
    MH(n_iters::Int)

Metropolis-Hasting sampler.

Usage:

```julia
MH(100, (:m, (x) -> Normal(x, 0.1)))
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

sample(gdemo([1.5, 2]), MH(1000, (:m, (x) -> Normal(x, 0.1)), :s)))
```
"""
mutable struct MH{T} <: InferenceAlgorithm
  n_iters   ::  Int       # number of iterations
  proposals ::  Dict{Symbol,Any}  # Proposals for paramters
  space     ::  Set{T}    # sampling space, emtpy means all
  gid       ::  Int       # group ID
end
function MH(n_iters::Int, space...)
  new_space = Set()
  proposals = Dict{Symbol,Any}()

  # parse random variables with their hypothetical proposal
  for element in space
      if isa(element, Symbol)
        push!(new_space, element)
      else
        @assert isa(element[1], Symbol) "[MH] ($element[1]) should be a Symbol. For proposal, use the syntax MH(N, (:m, (x) -> Normal(x, 0.1)))"
        push!(new_space, element[1])
        proposals[element[1]] = element[2]
      end
  end
  set = Set(new_space)
  MH{eltype(set)}(n_iters, proposals, set, 0)
end
MH{T}(alg::MH, new_gid::Int) where T = MH{T}(alg.n_iters, alg.proposals, alg.space, new_gid)

Sampler(alg::MH) = begin
  alg_str = "MH"

  # Sanity check for space
  if alg.gid == 0 && !isempty(alg.space)
    @assert issubset(Turing._compiler_[:pvars], alg.space) "[$alg_str] symbols specified to samplers ($alg.space) doesn't cover the model parameters ($(Turing._compiler_[:pvars]))"
    if Turing._compiler_[:pvars] != alg.space
      warn("[$alg_str] extra parameters specified by samplers don't exist in model: $(setdiff(alg.space, Turing._compiler_[:pvars]))")
    end
  end

  info = Dict{Symbol, Any}()
  info[:accept_his] = []
  info[:total_eval_num] = 0
  info[:proposal_ratio] = 0.0
  info[:prior_prob] = 0.0
  info[:violating_support] = false

  Sampler(alg, info)
end

propose(model::Function, spl::Sampler{<:MH}, vi::VarInfo) = begin
  spl.info[:proposal_ratio] = 0.0
  spl.info[:prior_prob] = 0.0
  spl.info[:violating_support] = false
  runmodel!(model, vi ,spl)
end

step(model::Function, spl::Sampler{<:MH}, vi::VarInfo, is_first::Bool) = begin
  if is_first
    push!(spl.info[:accept_his], true)
    vi

  else
    if spl.alg.gid != 0 # Recompute joint in logp
      runmodel!(model, vi, nothing)
    end
    old_θ = copy(vi[spl])
    old_logp = getlogp(vi)

    @debug "Propose new parameters from proposals..."
    propose(model, spl, vi)

    @debug "computing accept rate α..."
    is_accept, logα = mh_accept(-old_logp, -getlogp(vi), spl.info[:proposal_ratio])

    @debug "decide wether to accept..."
    if is_accept && !spl.info[:violating_support]  # accepted
      push!(spl.info[:accept_his], true)
    else                      # rejected
      push!(spl.info[:accept_his], false)
      vi[spl] = old_θ         # reset Θ
      setlogp!(vi, old_logp)  # reset logp
    end

    vi
  end
end

function sample(model::Function, alg::MH;
                save_state=false,         # flag for state saving
                resume_from=nothing,      # chain to continue
                reuse_spl_n=0,            # flag for spl re-using
                )

  spl = reuse_spl_n > 0 ?
        resume_from.info[:spl] :
        Sampler(alg)
  alg_str = "MH"

  # Initialization
  time_total = 0.0
  n = reuse_spl_n > 0 ?
      reuse_spl_n :
      alg.n_iters
  samples = Array{Sample}(undef, n)
  weight = 1 / n
  for i = 1:n
    samples[i] = Sample(weight, Dict{Symbol, Any}())
  end

    vi = if resume_from == nothing
        vi_ = VarInfo()
        Base.invokelatest(model, vi_, HamiltonianRobustInit())
        vi_
    else
        resume_from.info[:vi]
    end

  if spl.alg.gid == 0
    runmodel!(model, vi, spl)
  end

  # MH steps
  PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0))
  for i = 1:n
    @debug "$alg_str stepping..."

    time_elapsed = @elapsed vi = step(model, spl, vi, i == 1)
    time_total += time_elapsed

    if spl.info[:accept_his][end]     # accepted => store the new predcits
      samples[i].value = Sample(vi, spl).value
    else                              # rejected => store the previous predcits
      samples[i] = samples[i - 1]
    end
    samples[i].value[:elapsed] = time_elapsed

    PROGRESS[] && (ProgressMeter.next!(spl.info[:progress]))
  end

  println("[$alg_str] Finished with")
  println("  Running time        = $time_total;")
  accept_rate = sum(spl.info[:accept_his]) / n  # calculate the accept rate
  println("  Accept rate         = $accept_rate;")

  if resume_from != nothing   # concat samples
    pushfirst!(samples, resume_from.value2...)
  end
  c = Chain(0, samples)       # wrap the result by Chain
  if save_state               # save state
    save!(c, spl, model, vi)
  end

  c
end

assume(spl::Sampler{<:MH}, dist::Distribution, vn::VarName, vi::VarInfo) = begin
    if isempty(spl.alg.space) || vn.sym in spl.alg.space
      if ~haskey(vi, vn) error("[MH] does not handle stochastic existence yet") end
      old_val = vi[vn]

      if vn.sym in keys(spl.alg.proposals) # Custom proposal for this parameter
        proposal = spl.alg.proposals[vn.sym](old_val)

        if typeof(proposal) == Distributions.Normal{Float64} # If Gaussian proposal
          σ = std(proposal)
          lb = support(dist).lb
          ub = support(dist).ub
          stdG = Normal()
          r = rand(TruncatedNormal(proposal.μ, proposal.σ, lb, ub))
          # cf http://fsaad.scripts.mit.edu/randomseed/metropolis-hastings-sampling-with-gaussian-drift-proposal-on-bounded-support/
          spl.info[:proposal_ratio] += log(cdf(stdG, (ub-old_val)/σ) - cdf(stdG,(lb-old_val)/σ))
          spl.info[:proposal_ratio] -= log(cdf(stdG, (ub-r)/σ) - cdf(stdG,(lb-r)/σ))

        else # Other than Gaussian proposal
          r = rand(proposal)
          if (r < support(dist).lb) | (r > support(dist).ub) # check if value lies in support
            spl.info[:violating_support] = true
            r = old_val
          end
          spl.info[:proposal_ratio] -= logpdf(proposal, r) # accumulate pdf of proposal
          reverse_proposal = spl.alg.proposals[vn.sym](r)
          spl.info[:proposal_ratio] += logpdf(reverse_proposal, old_val)
        end

      else # Prior as proposal
        r = rand(dist)
        spl.info[:proposal_ratio] += (logpdf(dist, old_val) - logpdf(dist, r))
      end

      spl.info[:prior_prob] += logpdf(dist, r) # accumulate prior for PMMH
      vi[vn] = vectorize(dist, r)
      setgid!(vi, spl.alg.gid, vn)
    else
      r = vi[vn]
    end

    # acclogp!(vi, logpdf(dist, r)) # accumulate pdf of prior
    r, logpdf(dist, r)
end

assume(spl::Sampler{<:MH}, dists::Vector{D}, vn::VarName, var::Any, vi::VarInfo) where D<:Distribution =
  error("[Turing] MH doesn't support vectorizing assume statement")

observe(spl::Sampler{<:MH}, d::Distribution, value::Any, vi::VarInfo) =
  observe(nothing, d, value, vi)  # accumulate pdf of likelihood

observe(spl::Sampler{<:MH}, ds::Vector{D}, value::Any, vi::VarInfo)  where D<:Distribution =
  observe(nothing, ds, value, vi) # accumulate pdf of likelihood


text:
using AirSeaFluxes
using Test

@testset "AirSeaFluxes.jl" begin
    mld=10 #mixed layer depth (m)
    tim=86400*30 #relaxation time scale (s)
    pisvel=mld/tim #piston velocity (m/s)
    Co=0.0 #ocean value (e.g. concentration of some compound)
    Ca=1.0 #atmospeheric value (e.g. equivalent compound concentration)
    flx=simpleflux(Ca,Co,pisvel)
    @test isapprox(flx,3.858024691358025e-6; atol=1e-2)

    atemp=300.;aqh=0.001;speed=1.; sst=10.
    all=bulkformulae(atemp,aqh,speed,sst)
    @test isapprox(all["hl"],-3.0606099804357885; rtol=1e-2)
    @test isapprox(all["hs"],2.0282408526727473; rtol=1e-2)
    @test isapprox(all["evap"],1.2244888899523058e-9; rtol=1e-2)
    @test isapprox(all["tau"],0.00915548218468587; rtol=1e-2)
end


text:
# The code to simulate 200 artificial participants given the SurNoR parameters
# fitted to behavior
using SurNoR_2020
using MAT
using PyPlot
using CSV
using JLD

i = 3
# --------------------------------------------------------------------------
# Algorithm
update_form = "VarSMiLe"

if update_form == "SMiLe"
    surprise_form="CC"
else
    surprise_form="BF"
end

novelty_form="log-count"
update_form_model_free = "Q_learning"
gamma_form="satur"

# --------------------------------------------------------------------------
# Defining Hyper Parameters
HyperParam = SurNoR_2020.Str_HyperParam(; gamma_form=gamma_form,
                                          surprise_form=surprise_form,
                                          novelty_form=novelty_form,
                                          update_form=update_form,
                                          update_form_model_free=update_form_model_free)

# --------------------------------------------------------------------------
# Fitting the models
Block_Set = 1:2
Epi_Set = 1:5

# --------------------------------------------------------------------------
# Defining Parameters
Path_Params = "data/Fitted_parameters/Overall/OverAll_SurNoR_BestParams_diff5.csv"
Param_Data = CSV.read(Path_Params)

η0 = Array(Param_Data[i,2:end])
σ0 = [0.01,1e-5,2,0.1,0.1,0.01,0.1,0.01,0.1,0.1,0.1,1,5,0.1,0.1,0.1,1,2]
noise_free = true

# ------------------------------------------------------------------------------
# Loading the environment
Path = string("data/HardEnvironment11s.mat")
EnvGraph_MAT = matread(Path)
EnvGraph1 = SurNoR_2020.Str_EnvGraph(; TranMat = EnvGraph_MAT["Hard_NoSwitch"]["tm"].-1,
                                            InitObs = EnvGraph_MAT["Hard_NoSwitch"]["initS"].-1,
                                            Goal=0)

EnvGraph2 = SurNoR_2020.Str_EnvGraph(; TranMat = EnvGraph_MAT["Hard_Switch411"]["tm"].-1,
                                            InitObs = EnvGraph_MAT["Hard_Switch411"]["initS"].-1,
                                            Goal=0)

InitObs = [9,6,8,4,2]
for i = 1:5
  EnvGraph1.InitObs[i] = InitObs[i]
  EnvGraph2.InitObs[i] = InitObs[i]
end

EnvGraph_Array = [EnvGraph1, EnvGraph2]

# ------------------------------------------------------------------------------
# Simulate
N = 200
Seed_start = 2021
N_lost = [0]
SimAgent = []
for n = 1:N
  println("----------------------------------------")
  @show n
  @show Seed_start + n
  Param = SurNoR_2020.Func_param_generator(η0;
                                           σ = σ0./10,
                                           noise_free=noise_free,
                                           Seed= Seed_start + n)
  Temp_Age = SurNoR_2020.Func_simulate_Agent_Blocks(EnvGraph_Array;
                                                  Epi_Set=1:5, Block_Set=1:2,
                                                  Param=Param,
                                                  HyperParam=HyperParam,
                                                  Greedy = false,
                                                  Seed= Seed_start + n,
                                                  T_stop = 500)
  if size(Temp_Age) == (2,5)
    @show length(Temp_Age[1,1].Input.act')
    @show length(Temp_Age[2,1].Input.act')
    push!(SimAgent,Temp_Age)
  else
    N_lost[1] = N_lost[1]+1
  end
end

N = N - N_lost[1]
Inputs = Array{SurNoR_2020.Str_Input,3}(undef,N,2,5)
for n = 1:N
  for Block = 1:2
    for Epi = 1:5
      Inputs[n,Block,Epi] = SimAgent[n][Block,Epi].Input
    end
  end
end
save("data/Simulated_data/Simulated_data.jld", "Inputs", Inputs, "N_lost", N_lost)


text:
export Material

# --------- #
# Materials #
# --------- #
"""
    Material

This Type stores information about the material of an [`Object`](@ref).
We compute the color of a point by dispatching on the type of the material.
That is why we store `nothing` when a particular field is absent.

!!! note
    `texture_*` and `uv_coordinates` fields are applicable only if the
    object used is a [`Triangle`](@ref). If any of the texture fields are
    present, `uv_coordinates` must be present.

### Fields:

* `color_ambient`     - The Ambient Color of the Material.
* `color_diffuse`     - The Diffuse Color of the Material.
* `color_specular`    - The Specular Color of the Material.
* `specular_exponent` - Specular Exponent of the Material.
* `reflection`        - The reflection coefficient of the material.
* `texture_ambient`   - Stores the ambient texture image (if present) in Vec3 format.
* `texture_diffuse`   - Stores the diffuse texture image (if present) in Vec3 format.
* `texture_specular`  - Stores the specular texture image (if present) in VEc3 format.
* `uv_coordinates`    - The uv coordinates are stored as a vector (of length 3) of tuples.

### Available Constructors

```
Material(;color_ambient = Vec3(1.0f0), color_diffuse = Vec3(1.0f0),
          color_specular = Vec3(1.0f0), specular_exponent::Real = 50.0f0,
          reflection::Real = 0.5f0, texture_ambient = nothing,
          texture_diffuse = nothing, texture_specular = nothing,
          uv_coordinates = nothing)
```
"""
struct Material{T<:AbstractArray, R<:AbstractVector, U<:Union{Vec3, Nothing},
                V<:Union{Vec3, Nothing}, W<:Union{Vec3, Nothing},
                S<:Union{Vector, Nothing}}
    # Color Information
    color_ambient::Vec3{T}
    color_diffuse::Vec3{T}
    color_specular::Vec3{T}
    # Surface Properties
    specular_exponent::R
    reflection::R
    # Texture Information
    texture_ambient::U
    texture_diffuse::V
    texture_specular::W
    # UV coordinates (relevant only for triangles)
    uv_coordinates::S
end

Material(;color_ambient = Vec3(1.0f0), color_diffuse = Vec3(1.0f0),
          color_specular = Vec3(1.0f0), specular_exponent::Real = 50.0f0,
          reflection::Real = 0.5f0, texture_ambient = nothing,
          texture_diffuse = nothing, texture_specular = nothing,
          uv_coordinates = nothing) =
    Material(color_ambient, color_diffuse, color_specular, [specular_exponent],
             [reflection], texture_ambient, texture_diffuse, texture_specular,
             uv_coordinates)

@diffops Material

function Base.zero(m::Material)
    texture_ambient = isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient)
    texture_diffuse = isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse)
    texture_specular = isnothing(m.texture_specular) ? nothing : zero(m.texture_specular)
    uv_coordinates = isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)
    return Material(zero(m.color_ambient), zero(m.color_diffuse), zero(m.color_specular),
                    zero(m.specular_exponent), zero(m.reflection), texture_ambient,
                    texture_diffuse, texture_specular, uv_coordinates)
end

"""
    get_color(m::Material, pt::Vec3, ::Val{T}, obj::Object)

Returns the color at the point `pt`. We use `T` and the type of the Material `m`
for efficiently dispatching to the right function. The right function is determined
by the presence/absence of the texture field. The possible values of `T` are
`:ambient`, `:diffuse`, and `:specular`.
"""
get_color(m::Material{T, R, Nothing, V, W, S}, pt::Vec3,
          ::Val{:ambient}, obj) where {T<:AbstractArray, R<:AbstractVector,
                                       V<:Union{Vec3, Nothing}, W<:Union{Vec3, Nothing},
                                       S<:Union{Vector, Nothing}} =
    m.color_ambient

get_color(m::Material{T, R, U, Nothing, W, S}, pt::Vec3,
          ::Val{:diffuse}, obj) where {T<:AbstractArray, R<:AbstractVector,
                                       U<:Union{Vec3, Nothing}, W<:Union{Vec3, Nothing},
                                       S<:Union{Vector, Nothing}} =
    m.color_diffuse

get_color(m::Material{T, R, U, V, Nothing, S}, pt::Vec3,
          ::Val{:specular}, obj) where {T<:AbstractArray, R<:AbstractVector,
                                        U<:Union{Vec3, Nothing}, V<:Union{Vec3, Nothing},
                                        S<:Union{Vector, Nothing}} =
    m.color_specular

# This function is only available for triangles. Maybe I should put a
# type constraint
# The three functions are present only for type inference. Ideally Julia
# should figure it out not sure why it is not the case.
function get_color(m::Material, pt::Vec3, ::Val{:ambient}, obj)
    v1v2 = obj.v2 - obj.v1
    v1v3 = obj.v3 - obj.v1
    normal = cross(v1v2, v1v3)
    denom = l2norm(normal)

    edge2 = obj.v3 - obj.v2
    vp2 = pt - obj.v2
    u_bary = dot(normal, cross(edge2, vp2)) ./ denom

    edge3 = obj.v1 - obj.v3
    vp3 = pt - obj.v3
    v_bary = dot(normal, cross(edge3, vp3)) ./ denom

    w_bary = 1 .- u_bary .- v_bary

    u = u_bary .* m.uv_coordinates[1][1] .+
        v_bary .* m.uv_coordinates[2][1] .+
        w_bary .* m.uv_coordinates[3][1]
    v = u_bary .* m.uv_coordinates[1][2] .+
        v_bary .* m.uv_coordinates[2][2] .+
        w_bary .* m.uv_coordinates[3][2]

    # which_texture, which_color = infer_index(f)
    image = Zygote.literal_getproperty(m, Val(:texture_ambient))
    width, height = size(image)
    x_val = mod1.((Int.(ceil.(u .* width)) .- 1), width)
    y_val = mod1.((Int.(ceil.(v .* height)) .- 1), height)
    # Zygote doesn't like me using a splat operation here :'(
    img = image[x_val .+ (y_val .- 1) .* stride(image.x, 2)]
    return Zygote.literal_getproperty(m, Val(:color_ambient)) * Vec3(img.x, img.y, img.z)
end

function get_color(m::Material, pt::Vec3, ::Val{:diffuse}, obj)
    v1v2 = obj.v2 - obj.v1
    v1v3 = obj.v3 - obj.v1
    normal = cross(v1v2, v1v3)
    denom = l2norm(normal)

    edge2 = obj.v3 - obj.v2
    vp2 = pt - obj.v2
    u_bary = dot(normal, cross(edge2, vp2)) ./ denom

    edge3 = obj.v1 - obj.v3
    vp3 = pt - obj.v3
    v_bary = dot(normal, cross(edge3, vp3)) ./ denom

    w_bary = 1 .- u_bary .- v_bary

    u = u_bary .* m.uv_coordinates[1][1] .+
        v_bary .* m.uv_coordinates[2][1] .+
        w_bary .* m.uv_coordinates[3][1]
    v = u_bary .* m.uv_coordinates[1][2] .+
        v_bary .* m.uv_coordinates[2][2] .+
        w_bary .* m.uv_coordinates[3][2]

    # which_texture, which_color = infer_index(f)
    image = Zygote.literal_getproperty(m, Val(:texture_diffuse))
    width, height = size(image)
    x_val = mod1.((Int.(ceil.(u .* width)) .- 1), width)
    y_val = mod1.((Int.(ceil.(v .* height)) .- 1), height)
    # Zygote doesn't like me using a splat operation here :'(
    img = image[x_val .+ (y_val .- 1) .* stride(image.x, 2)]
    return Zygote.literal_getproperty(m, Val(:color_diffuse)) * Vec3(img.x, img.y, img.z)
end

function get_color(m::Material, pt::Vec3, ::Val{:specular}, obj)
    v1v2 = obj.v2 - obj.v1
    v1v3 = obj.v3 - obj.v1
    normal = cross(v1v2, v1v3)
    denom = l2norm(normal)

    edge2 = obj.v3 - obj.v2
    vp2 = pt - obj.v2
    u_bary = dot(normal, cross(edge2, vp2)) ./ denom

    edge3 = obj.v1 - obj.v3
    vp3 = pt - obj.v3
    v_bary = dot(normal, cross(edge3, vp3)) ./ denom

    w_bary = 1 .- u_bary .- v_bary

    u = u_bary .* m.uv_coordinates[1][1] .+
        v_bary .* m.uv_coordinates[2][1] .+
        w_bary .* m.uv_coordinates[3][1]
    v = u_bary .* m.uv_coordinates[1][2] .+
        v_bary .* m.uv_coordinates[2][2] .+
        w_bary .* m.uv_coordinates[3][2]

    # which_texture, which_color = infer_index(f)
    image = Zygote.literal_getproperty(m, Val(:texture_specular))
    width, height = size(image)
    x_val = mod1.((Int.(ceil.(u .* width)) .- 1), width)
    y_val = mod1.((Int.(ceil.(v .* height)) .- 1), height)
    # Zygote doesn't like me using a splat operation here :'(
    img = image[x_val .+ (y_val .- 1) .* stride(image.x, 2)]
    return Zygote.literal_getproperty(m, Val(:color_specular)) * Vec3(img.x, img.y, img.z)
end

"""
    specular_exponent(m::Material)

Returns the `specular_exponent` of the Material `m`.
"""
specular_exponent(m::Material) = m.specular_exponent[]

"""
    reflection(m::Material)

Returns the `reflection coefficient` of the Material `m`.
"""
reflection(m::Material) = m.reflection[]


text:
#jl #! format: off
# # Bar and stack plots with [PowerGraphics.jl](github.com/nrel-siip/PowerGraphics.jl)
# PowerGraphics also provides some basic specifications for plotting `SimulationResults`.
# This example demonstrates some simple plotting capabilities using different Plots.julia
# backends.
#
# The plotting capabilities use the Julia Plots package which can generate plots using
# several different graphics packages. We'll use GR.jl and PlotlyJS.jl.
#
# ## Dependencies
using SIIPExamples #for path locations
pkgpath = pkgdir(SIIPExamples)
using PowerSystems #to load results
using PowerSimulations #to load results
using PowerGraphics
using PowerSystemCaseBuilder

# ### Results file
# If you have already run some of the other examples, you should have generated some results
# (If you haven't run some of the other simulations, you can run
include(
    joinpath(pkgpath, "test", "3_PowerSimulations_examples", "03_5_bus_mkt_simulation.jl"),
)

# Alternatively, you can load the results into memory with:
# ```julia
# simulation_folder = joinpath(pkgdir(SIIPExamples), "rts-test")
# simulation_folder =
#     joinpath(simulation_folder, "$(maximum(parse.(Int64,readdir(simulation_folder))))")
#
# results = SimulationResults(simulation_folder);
# uc_results = get_decision_problem_results(results, "UC")
# ```

# Since some of the plotting capabilities rely on input data as well as output data (e.g. fuel plots)
# but the result deserialization doesn't load the `System`, we can add the `System` to the `results`
# so that the plotting routines can find the requisite data.
set_system!(uc_results, sys_DA)

# ## Plots
# By default, PowerGraphics uses the GR graphics package as the backend for Plots.jl to
# generate figures. This creates static plots and should execute without any extra steps.
# For example, we can create a plot of a particular variable in the `uc_results` object:
gr() # loads the GR backend
timestamps = get_realized_timestamps(uc_results)
variable = read_realized_variable(uc_results, "ActivePowerVariable__ThermalStandard")

plot_dataframe(variable, timestamps)

# However, interactive plotting can generate much more insightful figures, especially when
# creating somewhat complex stacked figures. So, we can use the PlotlyJS backend for Plots,
# but it requires that PlotlyJS.jl is installed in your Project.toml (if in a notebook,
# WebIO.jl is required too). To startup the PlotlyJS backend, run:
plotlyjs()

# PowerGraphics creates an un-stacked line plot by default, but supports kwargs to
# create a variety of different figure styles. For example, a stacked area figure can be
# created with the `stack = true` kwarg:

plot_dataframe(variable, timestamps; stack = true)

# Or a bar chart can be created with `bar = true`:
plot_dataframe(variable, timestamps; bar = true)

# Or a stacked bar chart...
plot_dataframe(variable, timestamps; bar = true, stack = true)

# PowerGraphics also supports some basic aggregation to create cleaner plots. For example,
# we can create a plot of the different variables:
generation = get_generation_data(uc_results)
plot_pgdata(generation, stack = true)

reserves = get_service_data(uc_results)
plot_pgdata(reserves)

# Another standard aggregation is available to plot demand values:
plot_demand(uc_results)

# The `plot_demand` function can also be called with the `System` rather than the `StageResults`
# to inspect the input data. This method can also display demands aggregated by a specified
# `<:Topology`:
plot_demand(uc_results.system, aggregation = Area)

# Another standard aggregation exists based on the fuel categories of the generators in the
# `System`
plot_fuel(uc_results)

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
#=
test_parser3:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-03-05
=#

using Test
using HyperCollate, LightGraphs
include("util.jl")

@testset "hypercollate" begin

    @testset "test 1" begin
        xml = "<text><s><subst>Dit kwam van een<add>De</add></subst> te streng doorgedreven rationalisatie</s></text>"
        @show(xml)
        tokens = tokenize(xml)
        @test map(string_value,tokens) == ["<text>", "<s>", "<subst>", "", "Dit kwam van een", "", "<add>", "De", "</add>", "</subst>", " te streng doorgedreven rationalisatie", "</s>", "</text>"]

        g = to_graph(xml)
        @show(g.graph)
        ts = topological_sort_by_dfs(g)
        @show(ts)
        dot = to_dot(g)
        expected="""
        digraph VariantGraph {
            rankdir=LR
            labelloc=b
            color=white
            edge [arrowsize=0.5]
            v1[shape=circle;width=0.05;label=""]
            v2[shape=box;label="Dit kwam van een"]
            v3[shape=box;label="De"]
            v4[shape=circle;width=0.05;label=""]
            v5[shape=box;label=" te streng doorgedreven rationalisatie"]
            v1 -> v2
            v1 -> v3
            v2 -> v4
            v3 -> v4
            v4 -> v5
        }
        """
        _test_normalized_strings_are_equal(dot,expected)
        _print_dot(dot)
    end

    @testset "test 2" begin
        expected = """
        digraph VariantGraph {
            rankdir=LR
            labelloc=b
            color=white
            edge [arrowsize=0.5]
            v1[shape=box;label="To be, or "]
            v2[shape=circle;width=0.05;label=""]
            v3[shape=box;label="whatever"]
            v4[shape=box;label="not to "]
            v5[shape=circle;width=0.05;label=""]
            v6[shape=box;label="butterfly"]
            v7[shape=box;label="be"]
            v8[shape=circle;width=0.05;label=""]
            v9[shape=circle;width=0.05;label=""]
            v1 -> v2
            v2 -> v3
            v2 -> v4
            v3 -> v9
            v4 -> v5
            v5 -> v6
            v5 -> v7
            v6 -> v8
            v7 -> v8
            v8 -> v9
        }
        """
        dot = (to_dot(to_graph("<xml>To be, or <subst>whatever<add>not to <subst>butterfly<add>be</add></subst></add></subst></xml>")))
        _test_normalized_strings_are_equal(dot,expected)
        _print_dot(dot)
    end

    @testset "test 3" begin
        expected = """
        digraph VariantGraph {
            rankdir=LR
            labelloc=b
            color=white
            edge [arrowsize=0.5]
            v1[shape=box;label="pre "]
            v2[shape=circle;width=0.05;label=""]
            v3[shape=box;label="one "]
            v4[shape=box;label="golden goose"]
            v5[shape=box;label=" waddling"]
            v6[shape=box;label="two "]
            v7[shape=box;label="cooked hens"]
            v8[shape=box;label=" smelling"]
            v9[shape=box;label="three "]
            v10[shape=box;label="roasted ducks"]
            v11[shape=box;label=" cooling"]
            v12[shape=circle;width=0.05;label=""]
            v13[shape=box;label=" post"]
            v1 -> v2
            v2 -> v3
            v2 -> v6
            v2 -> v9
            v3 -> v4
            v4 -> v5
            v5 -> v12
            v6 -> v7
            v7 -> v8
            v8 -> v12
            v9 -> v10
            v10 -> v11
            v11 -> v12
            v12 -> v13
        }
        """
        dot = (to_dot(to_graph("<x>pre <app><rdg>one golden goose waddling</rdg><rdg>two cooked hens smelling</rdg><rdg>three roasted ducks cooling</rdg></app> post</x>")))
        _test_normalized_strings_are_equal(dot,expected)
        _print_dot(dot)
    end

    @testset "tokenize" begin
        xml = """
        <xml>

                newlines are fun!

        </xml>
        """
        tokens = tokenize(xml)
        @show(tokens)
        @test length(tokens) == 7
    end

    @testset "solitary add" begin
        xml = "<x>bla1 <add>bla2</add> bla3</x>"
        expected = """
        digraph VariantGraph {
            rankdir=LR
            labelloc=b
            color=white
            edge [arrowsize=0.5]
            v1[shape=box;label="bla1 "]
            v2[shape=circle;width=0.05;label=""]
            v3[shape=box;label="bla2"]
            v4[shape=circle;width=0.05;label=""]
            v5[shape=box;label=" bla3"]
            v1 -> v2
            v2 -> v3
            v2 -> v4
            v3 -> v4
            v4 -> v5
        }
        """
        dot = (to_dot(to_graph(xml)))
        _test_normalized_strings_are_equal(dot,expected)
        _print_dot(dot)
    end

    @testset "solitary del" begin
        xml = "<x>bla1 bla2 bla3</x>"
        expected = """
        digraph VariantGraph {
            rankdir=LR
            labelloc=b
            color=white
            edge [arrowsize=0.5]
            v1[shape=box;label="bla1 "]
            v2[shape=circle;width=0.05;label=""]
            v3[shape=box;label="bla2"]
            v4[shape=circle;width=0.05;label=""]
            v5[shape=box;label=" bla3"]
            v1 -> v2
            v2 -> v3
            v2 -> v4
            v3 -> v4
            v4 -> v5
        }
        """
        dot = (to_dot(to_graph(xml)))
        _test_normalized_strings_are_equal(dot,expected)
        _print_dot(dot)
    end

#     println(to_dot(to_graph("<x>De kat krabt de krullen van de trap</x>")))
#     println(to_dot(to_graph("<x>Donald smacked <choice><option>Huey</option><option>Dewey</option><option>Louie</option></choice> on his beak.</x>")))
end


text:
@testset "episode_turn_buffer" begin
    b = episode_RTSA_buffer(; state_eltype = Int)
    @test length(b) == 0
    @test isfull(b) == false
    @test isempty(b) == true

    push!(b; reward = 0.0, terminal = true, state = 1, action = 1)

    @test length(b) == 0
    @test isfull(b) == false
    @test isempty(b) == false

    push!(b; reward = 1.0, terminal = false, state = 2, action = 2)

    @test length(b) == 1
    @test isfull(b) == false
    @test isempty(b) == false

    push!(b; reward = 2.0, terminal = true, state = 3, action = 3)

    @test length(b) == 2
    @test isfull(b) == true
    @test isempty(b) == false
    @test b[end] == (
        state = 2,
        action = 2,
        reward = 2.0,
        terminal = true,
        next_state = 3,
        next_action = 3,
    )

    push!(b; reward = 0.0, terminal = false, state = 3, action = 3)
    push!(b; reward = 3.0, terminal = false, state = 4, action = 4)
    @test b[end] == (
        state = 3,
        action = 3,
        reward = 3.0,
        terminal = false,
        next_state = 4,
        next_action = 4,
    )
    @test length(b) == 2
end


text:
# call this function to complete all the iteration update steps
function updateParameters!(sys::System, ens::Ensemble)
	# calculate useful vectors using CE-update formula
    rxnActual, rxnExpect = computeRxnStats(sys, ens)

    # update hyper-parameters
    meanHyper, varHyper = hyperUpdate(sys, ens)

    # update kinetic rate parameters
    θstep = computeRateMeans!(sys, rxnActual, rxnExpect, meanHyper)

    # update variances
    computeRateVar!(sys, ens, rxnActual, rxnExpect, θstep, varHyper)
    append!(sys.state.θs, sys.state.θ)
end

function computeRxnStats(sys::System, ens::Ensemble)

    len = length(ens)

    pl = length(sys.model.ps)
    # ---- pre-allocate
    rxnActual = zeros(Float64, (pl, len))
    rxnExpect = zeros(Float64, (pl, len))
    a  = zeros(Float64, sys.model.nr)

    for k in eachindex(ens)
        sys.model.F(ens[k])
        for j in sys.model.ps
            @inbounds for j2 in sys.model.mapping[j]
                rxnActual[j, k] += ens[k].ra[j2]
                rxnExpect[j, k] += ens[k].fly[j2]
                rxnExpect[j, k] += (sys.times[end] - ens[k].t) .*  ens[k].a[j2]
            end
        end
        rxnActual[:, k] ./= getCostVal(ens[k])
        rxnExpect[:, k] ./= getCostVal(ens[k])
    end
    return rxnActual, rxnExpect
end

function hyperUpdate(sys::System, ens::Ensemble)
    n = length(sys.model.bounds[1]) - length(sys.model.ps)
    meanHyper = Vector{Float64}(undef, n)
    varHyper = Vector{Float64}(undef, n)

    for j in eachindex(meanHyper)
        v = Float64[]
        for k in eachindex(ens)
            push!(v,log(ens[k].θ[j+length(sys.model.ps)]))
        end
        meanHyper[j] = mean(v)
        varHyper[j] = var(v)
    end
    return meanHyper, varHyper
end

function computeRateMeans!(sys::System, rxnActual::Matrix{Float64}, rxnExpect::Matrix{Float64}, arm)
    θstep = vec(sum(rxnActual, dims = 2) ./ sum(rxnExpect, dims = 2))
    if sys.routine.sampling == :log
        θstep = log.(θstep)
    end
    append!(θstep, arm)
    if sys.state.i == 1
        sys.state.θ = θstep
        return θstep
    else
        sys.state.θ .= sys.routine.smoothingParameter .* θstep .+ (1.0 - sys.routine.smoothingParameter) .* sys.state.θ
        return θstep
    end
end

function computeRateVar!(sys::System, ens::Ensemble, rxnActual::Matrix{Float64}, rxnExpect::Matrix{Float64}, θstep::Vector{Float64}, ars)

    if sys.routine.sampling == :log
        tmp = log.((rxnActual.+0.01) ./ rxnExpect)
        σ2step = vec(var(tmp, dims = 2))
    else
        σ2step = vec(var(tmp, dims = 2))
    end
    append!(σ2step, ars)

    if sys.state.i == 1
        sys.state.σ2 = diagm(σ2step)
    else
        β = 0.8 - 0.8 * (1.0 - 1.0 / sys.state.i)^5
        sys.state.σ2 .= β .* diagm(σ2step) .+ (1.0 - β) .* sys.state.σ2
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
module BDPlusRest
    using HTTP
    using JSON

    include("BDPlusExceptions.jl")

    # Current URL to discovery the bigquery REST API documents (JSON)
    const DISCOVERY_DOCUMENT = "https://bigquery.googleapis.com/discovery/v1/apis/bigquery/v2/rest"

    try
        bigquery_docs = HTTP.request("GET", DISCOVERY_DOCUMENT)
    catch
        raise_discoverydoc_err()
        exit()
    end

    struct rest_resources

    end
end
const

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
include("./Legendre_polynomials.jl")

function Spherical_harmonics_a!(x::Float64, lmax::Int64;
                                Yₗ = nothing, Pₗ = nothing)
    #### Returns an array of size lmax_half. Each cell contains one (even axisymmetric) spherical harmonic at point cos(θ) = x ####
    #### The "a" in the name is here to outline that this function only gives the "a"xisymmetric harmonics (i.e. m = 0).         ####
    # getting the constants we will need
    lmax_half = div(lmax, 2)

    # let's preallocate the arrays that haven't been given
    if isnothing(Yₗ)
        Yₗ = zeros(lmax_half) # WIll hold the spherical harmonics
    end
    if isnothing(Pₗ)
        Pₗ = zeros(lmax_half) # WIll hold the Legendre polynomials
    end

    # let's comptute the (even) Legendre polynomials at point x
    Legendre_polynomials!(x, lmax, Pₗ = Pₗ)

    # let's compute the spherical harmonics
    for l_half = 1:lmax_half
        l = 2 * l_half
        Yₗ[l_half] = Pₗ[l_half] * sqrt((2 * l + 1) / (4 * π)) # Yₗ₀ = Pₗ₀ * sqrt((2 l + 1) / 4 π)
    end

    return(Yₗ)
end

function Spherical_harmonics!(θ::Float64, ϕ::Float64, lmax::Int64;
                                Yₗₘ = nothing, Pₗₘ = nothing)
    #### Returns an array of size lmax_half * (2 lmax + 1). Each cell contains one (even) spherical harmonic at point (θ, ϕ) ####
    #### This function has no "a" in its title because it returns ALL the spherical harmonics (i.e. m = 0 and m != 0)        ####
    # getting the constants we will need
    lmax_half = div(lmax, 2)

    # let's preallocate the arrays that haven't been given
    if isnothing(Yₗₘ)
        Yₗₘ = zeros(2 * lmax + 1, lmax_half) # Will hold the spherical harmonics
    end
    if isnothing(Pₗₘ)
        Pₗₘ = zeros(lmax + 1, lmax_half) # Will hold the normalized associated Legendre polynomials
    end

    # Let's compute the normalzed associated Legendre polynomials
    Normalized_associated_Legendre_polynomials!(cos(θ), lmax, P = Pₗₘ)

    # Let's compute the spherical harmonics
    for l_half = 1:lmax_half
        l = 2 * l_half

        for m = 0:l
            # The passage from Pₗₘ to Yₗₘ is not the same if m is positive, null or negative.
            # Cf Wikipedia : https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form
            if (m != 0)
                Yₗₘ[lmax + 1 - m, l_half] = (-1)^(m) * sqrt(2) * Pₗₘ[m + 1, l_half] * sin(m * ϕ)
                Yₗₘ[lmax + 1 + m, l_half] = (-1)^(m) * sqrt(2) * Pₗₘ[m + 1, l_half] * cos(m * ϕ)
            else
                Yₗₘ[lmax + 1, l_half] = Pₗₘ[1, l_half]
            end
        end
    end

    return(Yₗₘ)
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
using ModelObjectsLH, Test

mo = ModelObjectsLH;

function single_id_test()
    @testset "SingleId" begin
        id1 = SingleId(:id1, [1, 2])
        println(id1);
        @test mo.index(id1) == [1,2]
        @test has_index(id1)

        id11 = SingleId(:id1, [1, 2]);
        @test isequal(id1, id11)

        id2 = SingleId(:id2, 3)
        @test mo.index(id2) == [3]

        id3 = SingleId(:id3);
        @test isempty(mo.index(id3))
        @test !has_index(id3)
        @test isequal(id3, SingleId(:id3))

        @test isequal([id1, id2], [id1, id2])
        @test !isequal([id1, id1], [id2, id1])
        @test !isequal([id1, id2], [id1])

        id4 = SingleId(:id4);
        s4 = make_string(id4);
        @test isequal(s4, "id4")
        id4a = make_single_id(s4);
        @test isequal(id4, id4a)

        id5 = SingleId(:id5, [4, 2]);
        s5 = make_string(id5);
        id5a = make_single_id(s5);
        @test isequal(id5, id5a)

        id6 = SingleId(:id6, 4)
        s6 = make_string(id6);
        id6a = make_single_id(s6);
        @test isequal(id6, id6a)
    end
end

function dict_single_id_test()
    @testset "Dict Single Id" begin
        d = Dict{SingleId, Any}([SingleId(:id1) => 1, SingleId(:id2, 1) => 2]);
        for (id1, val) in d
            @test id1 isa SingleId;
            @test val isa Integer;
        end
    end
end

function broadcast_test()
    @testset "Broadcasting" begin
        sV = [SingleId(:x), SingleId(:y)];
        eqV = SingleId(:x) .== sV;
        @test eqV[1]
        @test !eqV[2]
    end
end

@testset "SingleId" begin
    single_id_test();
    dict_single_id_test();
    broadcast_test();
end

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


text:
module SAOImageDS9Tests
using Test
using SAOImageDS9, XPA
using XPA: TupleOf

proc = Base.Process[]

function ds9start(timeout::Real = 30.0)
    global proc
    elapsed = 0.0
    seconds = 0.5
    while true
        try
            return SAOImageDS9.connect()
        catch
            if length(proc) < 1
                push!(proc, run(`/usr/bin/ds9`; wait=false))
            end
            if elapsed > timeout
                error("cannot connect to SAOImage/DS9")
            end
            seconds = min(timeout - elapsed, 2*seconds)
            sleep(seconds)
        end
    end
end

function ds9kill()
    global proc
    if length(proc) ≥ 1
        kill(pop!(proc))
    end
    nothing
end

ds9start()

@testset "Get requests" begin
    @test typeof(SAOImageDS9.get(VersionNumber)) === VersionNumber
    @test typeof(SAOImageDS9.get(Int, "width")) === Int
    @test typeof(SAOImageDS9.get(Int, "height")) === Int
end

ds9kill()

end # module


text:
using GLPlot; GLPlot.init()
using Plots; glvisualize()
using GLVisualize, FileIO, Colors, Images

imfolder = filter(readdir(GLVisualize.assetpath())) do path
    endswith(path, "jpg") || endswith(path, "png")
end
images = map(imfolder) do impath
    convert(Matrix{RGBA{N0f8}}, restrict(loadasset(impath)))
end;

mean_colors = map(mean, images)
x = map(comp1, mean_colors); y = map(comp2, mean_colors); z = map(comp3, mean_colors);
p1 = scatter(
    x,y,z, markercolor=mean_colors,
    shape = :circle, markerstrokewidth = 1, markerstrokecolor = "white",
    hover = images, ms = 12
)
p2 = scatter(x,y,z, markerstrokecolor = mean_colors, shape = images, ms = 15)
plot(p1, p2)
gui()


text:
using VrpAnt
using Test

@testset "VrpAnt.jl" begin
    m = Map(12)

    @testset "Creation of the Map" begin
        @test length(m.cities) == 12

        @testset "Solution" begin
            @test m.solution.length == Inf
            @test m.solution.path == Vector{City}[]
            @test m.solution.state == InProgress()
            @test m.solution.usedSpace == 0
            @test m.solution.usedTime == 0
        end
        @test m.ways == Dict{VrpAnt.CityIndex,Dict{VrpAnt.CityIndex,VrpAnt.Way}}()
        @test length(m.listModification) == 12
        @test length(m.transitions) == 12

    end
    optimize!(m)

    @testset "After the optimize function" begin

    end
end


text:
using ModernGL, GLWindow, GLAbstraction, GLFW, ColorTypes, Reactive, GeometryTypes

const window = create_glcontext("Example")

const vert = vert"""
{{GLSL_VERSION}}

in vec2 vertices;
in vec2 texturecoordinates; // must be this name, because collect_for_gl assumes them

out vec2 f_uv;
void main() {
	f_uv = texturecoordinates;
	gl_Position = vec4(vertices, 0.0, 1.0);
}
"""

# you can also load the shader from a file, which you can then edit in any editor and the changes will show up in your opengl program.
#using FileIO; prrogram = TemplateProgram(load("path_to_frag.frag"), load("path_to_vert.vert"))
const frag = frag"""
{{GLSL_VERSION}}

out vec4 outColor;
uniform sampler2D image;
in vec2 f_uv;

void main() {

	outColor = texture(image, f_uv);
}
"""
program = LazyShader(vert, frag) #= {{GLSL_VERSION}} is a template, you could add your own with the kwarg view=Dict{Compat.UTF8String, Compat.UTF8String}(key->replacement) =#

tex = Texture([RGBA{U8}(x,y,sin(x*pi), 1.0) for x=0:0.1:1., y=0:0.1:1.]) #automatically creates the correct texture
data = merge(Dict(
	:image => tex,
    :primitive => GLUVMesh2D(SimpleRectangle{Float32}(-1,-1,2,2))
)) # Transforms the rectangle into a 2D mesh with uv coordinates and then extracts the buffers for the shader

robj = std_renderobject(data, program) # creates a renderable object from the shader and the data.

glClearColor(0, 0, 0, 1)

while isopen(window)
  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    render(robj)

  	swapbuffers(window)
  	pollevents()
end


text:
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
ex, ey, ez = (3,5,7)[1:ndims(Ω)]
(npx, npy, npz) = (1,1,1)[1:ndims(Ω)]
ClimateMachine.gpu_allowscalar(true)
grid = DiscontinuousSpectralElementGrid(Ω, elements = (ex, ey, ez), polynomialorder = (npx, npy, npz), array = ArrayType)
gridhelper = GridHelper(grid)
x, y, z = coordinates(grid)
ϕ =  ScalarField(copy(x), gridhelper)
@testset "Check Interpolation 3D" begin
    ϕ .= x
    @test ϕ(0.1, 0.2, 0.3) ≈ 0.1
    ϕ .= y
    @test ϕ(0.1, 0.2, 0.3) ≈ 0.2
    ϕ .= z
    @test ϕ(0.1, 0.2, 0.3) ≈ 0.3
end
##
Ω = Circle(0,2) × Circle(-4,2)
elements =  (3,5,7)[1:ndims(Ω)]
polynomialorder = (1,1,1)[1:ndims(Ω)]
ClimateMachine.gpu_allowscalar(true)
grid = DiscontinuousSpectralElementGrid(Ω, elements = elements, polynomialorder = polynomialorder, array = ArrayType)
gridhelper = GridHelper(grid)
x, y, z = coordinates(grid)
ϕ =  ScalarField(copy(x), gridhelper)
@testset "Check Interpolation 2D" begin
    ϕ .= x
    @test ϕ(0.1, 0.2) ≈ 0.1
    ϕ .= y
    @test ϕ(0.1, 0.2) ≈ 0.2
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using Bolt
using Test
using DelimitedFiles
using LinearAlgebra

@testset "FFTLog" begin
    N = 64
    μ = 0
    q = 0.0
    r₀ = 1.0
    L = 8.0
    Nhalf = N ÷ 2
    n = range(-Nhalf,Nhalf,length=N)
    r = r₀ .* 10 .^ (n .* L ./ N )
    pl = Bolt.plan_fftlog(r, μ, q, 1.0; kropt=true)
    aₙ = r .^ (μ + 1) .* exp.(-r.^2 / 2)
    y = similar(r, ComplexF64)
    fftdata = readdlm("data/fftlog_example.txt", ' ', Float64, '\n')

    # test forward
    mul!(y, pl, aₙ)
    f_ref = fftdata[:,2]
    @test all(abs.(y .- f_ref) .< 1e-15)
    @test isapprox(y, f_ref)

    # test backward
    y2 = similar(r, ComplexF64)
    ldiv!(y2, pl, y)
    @test all(abs.(y2 .- aₙ) .< 1e-15)

end

# test U with more naive result
# U_μ(μ, x) = 2^x * gamma((μ + 1 + x)/2) / gamma((μ + 1 - x)/2)

## hand-written ℋ derivative for testing
# function ℋ′(x, par::AbstractCosmoParams)
#     a = x2a(x)
#     return -H₀(par) * (2par.Ω_r + (par.Ω_b + par.Ω_m) * a - 2Ω_Λ(par) * a^4) /
#         (2 * a * √(par.Ω_r + (par.Ω_b + par.Ω_m) * a + Ω_Λ(par) * a^4))
# end


text:
export PlusOperator, TimesOperator, mul_coefficients

struct PlusOperator{T,BI} <: Operator{T}
    ops::Vector{Operator{T}}
    bandwidths::BI
    function PlusOperator{T,BI}(opsin::Vector{Operator{T}},bi::BI) where {T,BI}
        n,m=size(first(opsin))
        for k=2:length(opsin)
            @assert size(opsin[k],1)==n && size(opsin[k],2)==m
        end
        new{T,BI}(opsin,bi)
    end
end

Base.size(P::PlusOperator,k::Integer) = size(first(P.ops),k)

PlusOperator(opsin::Vector{Operator{T}},bi::Tuple{UT,VT}) where {T,UT,VT} =
    PlusOperator{T,typeof(bi)}(opsin,bi)

bandwidths(P::PlusOperator) = P.bandwidths

israggedbelow(P::PlusOperator) = isbandedbelow(P) || all(israggedbelow,P.ops)

for (OP,mn) in ((:colstart,:min),(:colstop,:max),(:rowstart,:min),(:rowstop,:max))
    defOP = Meta.parse("default_"*string(OP))
    @eval function $OP(P::PlusOperator,k::Integer)
        if isbanded(P)
            $defOP(P,k)
        else
            mapreduce(op->$OP(op,k),$mn,P.ops)
        end
    end
end

function PlusOperator(ops::Vector)
    # calculate bandwidths
    b1,b2=-720,-720  # approximates ∞,-∞
    for op in ops
        br=bandwidths(op)
        b1=max(br[1],b1)
        b2=max(br[end],b2)
    end
    PlusOperator(ops,(b1,b2))
end

function convert(::Type{Operator{T}},P::PlusOperator) where T
    if T==eltype(P)
        P
    else
        PlusOperator{T,typeof(P.bandwidths)}(Vector{Operator{T}}(P.ops),P.bandwidths)
    end
end

function promoteplus(opsin::Vector{Operator{T}}) where T
    ops=Vector{Operator{T}}()
    # prune zero ops
    for op in opsin
        if !iszeroop(op)
            push!(ops,op)
        end
    end
    PlusOperator(promotespaces(ops))
end

for OP in (:domainspace,:rangespace)
    @eval $OP(P::PlusOperator) = $OP(first(P.ops))
end

domain(P::PlusOperator) = commondomain(P.ops)

+(A::PlusOperator,B::PlusOperator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B.ops...])
+(A::PlusOperator,B::PlusOperator,C::PlusOperator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B),eltype(C))}[A.ops...,B.ops...,C.ops...])
+(A::PlusOperator,B::Operator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B])
+(A::PlusOperator,B::ZeroOperator) = A
+(A::PlusOperator,B::Operator,C::Operator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B),eltype(C))}[A.ops...,B,C])
+(A::Operator,B::PlusOperator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B))}[A,B.ops...])
+(A::ZeroOperator,B::PlusOperator) = B
+(A::Operator,B::Operator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B))}[A,B])
+(A::Operator,B::Operator,C::Operator) =
    promoteplus(Operator{promote_type(eltype(A),eltype(B),eltype(C))}[A,B,C])

Base.stride(P::PlusOperator)=mapreduce(stride,gcd,P.ops)

function getindex(P::PlusOperator{T},k::Integer...) where T
    ret=P.ops[1][k...]::T
    for j=2:length(P.ops)
        ret+=P.ops[j][k...]::T
    end
    ret
end

for TYP in (:RaggedMatrix,:Matrix,:BandedMatrix,
            :BlockBandedMatrix,:BandedBlockBandedMatrix)
    @eval begin
        $TYP(P::SubOperator{T,PP,NTuple{2,BlockRange1}}) where {T,PP<:PlusOperator} =
            convert_axpy!($TYP,P)   # use axpy! to copy
        $TYP(P::SubOperator{T,PP}) where {T,PP<:PlusOperator} =
            convert_axpy!($TYP,P)   # use axpy! to copy
        $TYP(P::SubOperator{T,PP,NTuple{2,UnitRange{Int}}}) where {T,PP<:PlusOperator} =
            convert_axpy!($TYP,P)   # use axpy! to copy
    end
end

function BLAS.axpy!(α,P::SubOperator{T,PP},A::AbstractMatrix) where {T,PP<:PlusOperator}
    for op in parent(P).ops
        BLAS.axpy!(α, view(op,P.indexes[1],P.indexes[2]), A)
    end

    A
end

+(A::Operator,f::Fun) = A+Multiplication(f,domainspace(A))
+(f::Fun,A::Operator) = Multiplication(f,domainspace(A))+A
-(A::Operator,f::Fun) = A+Multiplication(-f,domainspace(A))
-(f::Fun,A::Operator) = Multiplication(f,domainspace(A))-A

for TYP in (:ZeroOperator,:Operator)
    @eval function +(A::$TYP,B::ZeroOperator)
        if spacescompatible(A,B)
            A
        else
            promotespaces(A,B)[1]
        end
    end
end
+(A::ZeroOperator,B::Operator) = B+A

# We need to support A+1 in addition to A+I primarily for matrix case: A+Matrix(I,2,2)
for OP in (:+,:-)
    @eval begin
        $OP(c::Union{UniformScaling,Number},A::Operator) =
            $OP(convert(Operator{promote_type(eltype(A),eltype(c))},c),A)
        $OP(A::Operator,c::Union{UniformScaling,Number}) =
            $OP(A,convert(Operator{promote_type(eltype(A),eltype(c))},c))
    end
end

## Times Operator

struct ConstantTimesOperator{B,T} <: Operator{T}
    λ::T
    op::B
    ConstantTimesOperator{B,T}(c,op) where {B,T} = new{B,T}(c,op)
end
function ConstantTimesOperator(c::Number,op::Operator{TT}) where TT<:Number
    T=promote_type(typeof(c),eltype(op))
    B=convert(Operator{T},op)
    ConstantTimesOperator{typeof(B),T}(T(c),B)
end

ConstantTimesOperator(c::Number,op::ConstantTimesOperator) =
    ConstantTimesOperator(c*op.λ,op.op)

@wrapperstructure ConstantTimesOperator
@wrapperspaces ConstantTimesOperator

convert(::Type{T},C::ConstantTimesOperator) where {T<:Number} = T(C.λ)*convert(T,C.op)

choosedomainspace(C::ConstantTimesOperator,sp::Space) = choosedomainspace(C.op,sp)

for OP in (:promotedomainspace,:promoterangespace),SP in (:UnsetSpace,:Space)
    @eval $OP(C::ConstantTimesOperator,k::$SP) = ConstantTimesOperator(C.λ,$OP(C.op,k))
end

function convert(::Type{Operator{T}},C::ConstantTimesOperator) where T
    if T==eltype(C)
        C
    else
        op=convert(Operator{T},C.op)
        ConstantTimesOperator{typeof(op),T}(T(C.λ),op)
    end
end

getindex(P::ConstantTimesOperator,k::Integer...) =
    P.λ*P.op[k...]

for TYP in (:RaggedMatrix,:Matrix,:BandedMatrix,
            :BlockBandedMatrix,:BandedBlockBandedMatrix)
    @eval begin
        $TYP(S::SubOperator{T,OP,NTuple{2,BlockRange1}}) where {T,OP<:ConstantTimesOperator} =
            convert_axpy!($TYP, S)
        $TYP(S::SubOperator{T,OP,NTuple{2,UnitRange{Int}}}) where {T,OP<:ConstantTimesOperator} =
            convert_axpy!($TYP, S)
        $TYP(S::SubOperator{T,OP}) where {T,OP<:ConstantTimesOperator} =
            convert_axpy!($TYP, S)
    end
end

BLAS.axpy!(α,S::SubOperator{T,OP},A::AbstractMatrix) where {T,OP<:ConstantTimesOperator} =
    unwrap_axpy!(α*parent(S).λ,S,A)

struct TimesOperator{T,BI} <: Operator{T}
    ops::Vector{Operator{T}}
    bandwidths::BI

    function TimesOperator{T,BI}(ops::Vector{Operator{T}},bi::BI) where {T,BI}
        # check compatible
        for k=1:length(ops)-1
            @assert size(ops[k],2) == size(ops[k+1],1)
        end

        # remove TimesOperators buried inside ops
        hastimes = false
        for k=1:length(ops)-1
            @assert spacescompatible(domainspace(ops[k]),rangespace(ops[k+1]))
            hastimes = hastimes || isa(ops[k],TimesOperator)
        end

        if hastimes
            newops=Vector{Operator{T}}(0)
            for op in ops
               if isa(op,TimesOperator)
                    for op2 in op.ops
                        push!(newops,op2)
                    end
                else
                    push!(newops,op)
                end
            end
            ops=newops
        end

        new{T,BI}(ops,bi)
    end
end

function bandwidthsum(P,k)
    ret=0
    for op in P
        ret+=bandwidths(op)[k]
    end
    ret
end

bandwidthssum(P) = (bandwidthsum(P,1),bandwidthsum(P,2))

TimesOperator(ops::Vector{Operator{T}},bi::Tuple{N1,N2}) where {T,N1,N2} =
    TimesOperator{T,typeof(bi)}(ops,bi)

TimesOperator(ops::Vector{Operator{T}}) where {T} = TimesOperator(ops,bandwidthssum(ops))
TimesOperator(ops::Vector{OT}) where {OT<:Operator} =
    TimesOperator(convert(Vector{Operator{eltype(OT)}},ops),bandwidthssum(ops))

TimesOperator(A::TimesOperator,B::TimesOperator) =
    TimesOperator(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B.ops...])
TimesOperator(A::TimesOperator,B::Operator) =
    TimesOperator(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B])
TimesOperator(A::Operator,B::TimesOperator) =
    TimesOperator(Operator{promote_type(eltype(A),eltype(B))}[A,B.ops...])
TimesOperator(A::Operator,B::Operator) =
    TimesOperator(Operator{promote_type(eltype(A),eltype(B))}[A,B])

==(A::TimesOperator,B::TimesOperator)=A.ops==B.ops

function convert(::Type{Operator{T}},P::TimesOperator) where T
    if T==eltype(P)
        P
    else
        TimesOperator(Operator{T}[P.ops...])
    end
end

function promotetimes(opsin::Vector{B},dsp) where B<:Operator
    ops=Vector{Operator{mapreduce(eltype,promote_type,opsin)}}(undef,0)

    for k=length(opsin):-1:1
        if !isa(opsin[k],Conversion)
            op=promotedomainspace(opsin[k],dsp)
            if op==()
                # do nothing
            elseif isa(op,TimesOperator)
                for j=length(op.ops):-1:1
                    push!(ops,op.ops[j])
                end
                dsp=rangespace(op)
            else
                push!(ops,op)
                dsp=rangespace(op)
            end
        end
    end
    if isempty(ops)
        ConstantOperator(1.0,dsp)
    elseif length(ops)==1
        first(ops)
    else
        TimesOperator(reverse!(ops))  # change order in TImesOperator if this is slow
    end
end

promotetimes(opsin::Vector{B}) where {B<:Operator}=promotetimes(opsin,domainspace(last(opsin)))

domainspace(P::TimesOperator)=domainspace(last(P.ops))
rangespace(P::TimesOperator)=rangespace(first(P.ops))

domain(P::TimesOperator)=commondomain(P.ops)

bandwidths(P::TimesOperator) = P.bandwidths

israggedbelow(P::TimesOperator) = isbandedbelow(P) || all(israggedbelow,P.ops)

Base.stride(P::TimesOperator) = mapreduce(stride,gcd,P.ops)

for OP in (:rowstart,:rowstop)
    defOP=Meta.parse("default_"*string(OP))
    @eval function $OP(P::TimesOperator,k::Integer)
        if isbanded(P)
            return $defOP(P,k)
        end
        for j=eachindex(P.ops)
            k=$OP(P.ops[j],k)
        end
        k
    end
end

for OP in (:colstart,:colstop)
    defOP=Meta.parse("default_"*string(OP))
    @eval function $OP(P::TimesOperator, k::Integer)
        if isbanded(P)
            return $defOP(P, k)
        end
        for j=reverse(eachindex(P.ops))
            k=$OP(P.ops[j],k)
        end
        k
    end
end

getindex(P::TimesOperator,k::Integer,j::Integer) = P[k:k,j:j][1,1]
function getindex(P::TimesOperator,k::Integer)
    @assert isafunctional(P)
    P[1:1,k:k][1,1]
end

function getindex(P::TimesOperator,k::AbstractVector)
    @assert isafunctional(P)
    vec(Matrix(P[1:1,k]))
end

for TYP in (:Matrix, :BandedMatrix, :RaggedMatrix)
    @eval function $TYP(V::SubOperator{T,TO,Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,TO<:TimesOperator}
        P = parent(V)

        if isbanded(P)
            if $TYP ≠ BandedMatrix
                return $TYP(BandedMatrix(V))
            end
        elseif isbandedblockbanded(P)
            N = block(rangespace(P), last(parentindices(V)[1]))
            M = block(domainspace(P), last(parentindices(V)[2]))
            B = P[Block(1):N, Block(1):M]
            return $TYP(view(B, parentindices(V)...), _colstops(V))
        end

        kr,jr = parentindices(V)

        (isempty(kr) || isempty(jr)) && return $TYP(Zeros, V)

        if maximum(kr) > size(P,1) || maximum(jr) > size(P,2) ||
            minimum(kr) < 1 || minimum(jr) < 1
            throw(BoundsError())
        end

        @assert length(P.ops) ≥ 2
        if size(V,1)==0
            return $TYP(Zeros, V)
        end

        # find optimal truncations for each operator
        # by finding the non-zero entries
        krlin = Matrix{Union{Int,Infinity}}(undef,length(P.ops),2)

        krlin[1,1],krlin[1,2]=kr[1],kr[end]
        for m=1:length(P.ops)-1
            krlin[m+1,1]=rowstart(P.ops[m],krlin[m,1])
            krlin[m+1,2]=rowstop(P.ops[m],krlin[m,2])
        end
        krlin[end,1]=max(krlin[end,1],colstart(P.ops[end],jr[1]))
        krlin[end,2]=min(krlin[end,2],colstop(P.ops[end],jr[end]))
        for m=length(P.ops)-1:-1:2
            krlin[m,1]=max(krlin[m,1],colstart(P.ops[m],krlin[m+1,1]))
            krlin[m,2]=min(krlin[m,2],colstop(P.ops[m],krlin[m+1,2]))
        end

        krl = Matrix{Int}(krlin)

        # Check if any range is invalid, in which case return zero
        for m=1:length(P.ops)
            if krl[m,1]>krl[m,2]
                return $TYP(Zeros, V)
            end
        end

        # The following returns a banded Matrix with all rows
        # for large k its upper triangular
        BA = convert($TYP{T}, P.ops[end][krl[end,1]:krl[end,2],jr])
        for m = (length(P.ops)-1):-1:1
            BA = convert($TYP{T}, P.ops[m][krl[m,1]:krl[m,2],krl[m+1,1]:krl[m+1,2]])*BA
        end

        $TYP{T}(BA)
    end
end

for TYP in (:BlockBandedMatrix, :BandedBlockBandedMatrix)
    @eval function $TYP(V::SubOperator{T,TO,Tuple{BlockRange1,BlockRange1}}) where {T,TO<:TimesOperator}
        P = parent(V)
        KR,JR = parentindices(V)

        @assert length(P.ops) ≥ 2
        if size(V,1)==0 || isempty(KR) || isempty(JR)
            return $TYP(Zeros, V)
        end

        if Int(maximum(KR)) > nblocks(P,1) || Int(maximum(JR)) > nblocks(P,2) ||
            Int(minimum(KR)) < 1 || Int(minimum(JR)) < 1
            throw(BoundsError())
        end

        # find optimal truncations for each operator
        # by finding the non-zero entries
        KRlin = Matrix{Union{Block,Infinity}}(undef,length(P.ops),2)

        KRlin[1,1],KRlin[1,2] = first(KR),last(KR)
        for m=1:length(P.ops)-1
            KRlin[m+1,1]=blockrowstart(P.ops[m],KRlin[m,1])
            KRlin[m+1,2]=blockrowstop(P.ops[m],KRlin[m,2])
        end
        KRlin[end,1]=max(KRlin[end,1],blockcolstart(P.ops[end],first(JR)))
        KRlin[end,2]=min(KRlin[end,2],blockcolstop(P.ops[end],last(JR)))
        for m=length(P.ops)-1:-1:2
            KRlin[m,1]=max(KRlin[m,1],blockcolstart(P.ops[m],KRlin[m+1,1]))
            KRlin[m,2]=min(KRlin[m,2],blockcolstop(P.ops[m],KRlin[m+1,2]))
        end

        KRl = Matrix{Block{1}}(KRlin)

        # Check if any range is invalid, in which case return zero
        for m=1:length(P.ops)
            if KRl[m,1]>KRl[m,2]
                return $TYP(Zeros, V)
            end
        end

        # The following returns a banded Matrix with all rows
        # for large k its upper triangular
        BA = convert($TYP, view(P.ops[end],KRl[end,1]:KRl[end,2],JR))
        for m = (length(P.ops)-1):-1:1
            BA = convert($TYP, view(P.ops[m],KRl[m,1]:KRl[m,2],KRl[m+1,1]:KRl[m+1,2]))*BA
        end

        convert($TYP, BA)
    end
end

## Algebra: assume we promote

for OP in (:(adjoint),:(transpose))
    @eval $OP(A::TimesOperator)=TimesOperator(reverse!(map($OP,A.ops)))
end

*(A::TimesOperator,B::TimesOperator) =
    promotetimes(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B.ops...])
function *(A::TimesOperator,B::Operator)
    if isconstop(B)
        promotedomainspace(convert(Number,B)*A,domainspace(B))
    else
        promotetimes(Operator{promote_type(eltype(A),eltype(B))}[A.ops...,B])
    end
end
function *(A::Operator,B::TimesOperator)
    if isconstop(A)
        promoterangespace(convert(Number,A)*B,rangespace(A))
    else
        promotetimes(Operator{promote_type(eltype(A),eltype(B))}[A,B.ops...])
    end
end
function *(A::Operator,B::Operator)
    if isconstop(A)
        promoterangespace(convert(Number,A)*B,rangespace(A))
    elseif isconstop(B)
        promotedomainspace(convert(Number,B)*A,domainspace(B))
    else
        promotetimes(Operator{promote_type(eltype(A),eltype(B))}[A,B])
    end
end

# Conversions we always assume are intentional: no need to promote

*(A::ConversionWrapper{TO1},B::ConversionWrapper{TO}) where {TO1<:TimesOperator,TO<:TimesOperator} =
    ConversionWrapper(TimesOperator(A.op,B.op))
*(A::ConversionWrapper{TO},B::Conversion) where {TO<:TimesOperator} =
    ConversionWrapper(TimesOperator(A.op,B))
*(A::Conversion,B::ConversionWrapper{TO}) where {TO<:TimesOperator} =
    ConversionWrapper(TimesOperator(A,B.op))

*(A::Conversion,B::Conversion) = ConversionWrapper(TimesOperator(A,B))
*(A::Conversion,B::TimesOperator) = TimesOperator(A,B)
*(A::TimesOperator,B::Conversion) = TimesOperator(A,B)
*(A::Operator,B::Conversion) =
    isconstop(A) ? promoterangespace(convert(Number,A)*B,rangespace(A)) : TimesOperator(A,B)
*(A::Conversion,B::Operator) =
    isconstop(B) ? promotedomainspace(convert(Number,B)*A,domainspace(B)) : TimesOperator(A,B)

^(A::Operator, p::Integer) = Base.power_by_squaring(A, p)

+(A::Operator) = A
-(A::Operator) = ConstantTimesOperator(-1,A)
-(A::Operator,B::Operator) = A+(-B)

function *(f::Fun, A::Operator)
    if isafunctional(A) && (isinf(bandwidth(A,1)) || isinf(bandwidth(A,2)))
        LowRankOperator(f,A)
    else
        TimesOperator(Multiplication(f,rangespace(A)),A)
    end
end

function *(c::Number,A::Operator)
    if c==1
        A
    elseif c==0
        ZeroOperator(domainspace(A),rangespace(A))
    else
        ConstantTimesOperator(c,A)
    end
end
*(A::Operator,c::Number) = A*(c*one(domainspace(A)))

/(B::Operator,c::Number) = (1.0/c)*B
/(B::Operator,c::Fun) = (1.0/c)*B

## Operations
function mul_coefficients(A::Operator,b)
    n=size(b,1)
    ret = n>0 ? mul_coefficients(view(A,FiniteRange,1:n),b) : b
end

function mul_coefficients(A::TimesOperator,b)
    ret = b
    for k=length(A.ops):-1:1
        ret = mul_coefficients(A.ops[k],ret)
    end

    ret
end

function *(A::Operator, b)
    ds = domainspace(A)
    rs = rangespace(A)
    if isambiguous(ds)
        promotedomainspace(A,space(b))*b
    elseif isambiguous(rs)
        error("Assign spaces to $A before multiplying.")
    else
        Fun(rs,
            mul_coefficients(A,coefficients(b,ds)))
    end
end

mul_coefficients(A::PlusOperator,b::Fun) =
    mapreduce(x->mul_coefficients(x,b),+,A.ops)

*(A::Operator, b::AbstractMatrix{<:Fun}) = A*Fun(b)
*(A::Vector{<:Operator}, b::Fun) = map(a->a*b,convert(Array{Any,1},A))

## promotedomain

function promotedomainspace(P::PlusOperator{T},sp::Space,cursp::Space) where T
    if sp==cursp
        P
    else
        ops = [promotedomainspace(op,sp) for op in P.ops]
        promoteplus(Vector{Operator{mapreduce(eltype,promote_type,ops)}}(ops))
    end
end

function choosedomainspace(P::PlusOperator,sp::Space)
    ret=UnsetSpace()
    for op in P.ops
        sp2=choosedomainspace(op,sp)
        if !isa(sp2,AmbiguousSpace)  # we will ignore this result in hopes another opand
                                     # tells us a good space
            ret=union(ret,sp2)
        end
    end
    ret
end

function promotedomainspace(P::TimesOperator,sp::Space,cursp::Space)
    if sp==cursp
        P
    elseif length(P.ops)==2
        P.ops[1]*promotedomainspace(P.ops[end],sp)
    else
        promotetimes([P.ops[1:end-1];promotedomainspace(P.ops[end],sp)])
    end
end

function choosedomainspace(P::TimesOperator,sp::Space)
    for op in P.ops
        sp=choosedomainspace(op,sp)
    end
    sp
end


text:
export prior

"""
    prior(m, xs...)

Returns the minimal model required to sample random variables `xs...`. Useful for extracting a prior distribution from a joint model `m` by designating `xs...` and the variables they depend on as the prior and hyperpriors.

# Example

```jldoctest
m = @model n begin
    α ~ Gamma()
    β ~ Gamma()
    θ ~ Beta(α,β)
    x ~ Binomial(n, θ)
end;
Soss.prior(m, :x)

# output
@model begin
        β ~ Gamma()
        α ~ Gamma()
        θ ~ Beta(α, β)
    end
```
"""
prior(m::Model, xs...) = before(m, xs..., inclusive = false, strict = true)

export likelihood
"""
    likelihood(m, xs...)

Return a model with only the specified variables in the body. Required dependencies will be included as arguments.

"""
function likelihood(m::Model, xs...)
    M = getmodule(m)
    result = foldl(merge, (Model(M, findStatement(m,x)) for x in xs))

    for x in xs
        for pa in parents(digraph(m), x)
            pa ∈ xs && continue
            result = merge(result, Model(M, Arg(pa)))
        end
    end

    return result
end

export prune

"""
    prune(m, xs...)

Returns a model transformed by removing `xs...` and all variables that depend on `xs...`. Unneeded arguments are also removed.

# Examples

```jldoctest
m = @model n begin
    α ~ Gamma()
    β ~ Gamma()
    θ ~ Beta(α,β)
    x ~ Binomial(n, θ)
end;
prune(m, :θ)

# output
@model begin
        β ~ Gamma()
        α ~ Gamma()
    end
```

```jldoctest
m = @model n begin
    α ~ Gamma()
    β ~ Gamma()
    θ ~ Beta(α,β)
    x ~ Binomial(n, θ)
end;
prune(m, :n)

# output
@model begin
        β ~ Gamma()
        α ~ Gamma()
        θ ~ Beta(α, β)
    end
```
"""
prune(m::Model, xs...) = before(m, xs..., inclusive = false, strict = false)

export predictive
"""
    predictive(m, xs...)

Returns a model transformed by adding `xs...` to arguments with a body containing only statements that depend on `xs`, or statements that are depended upon by children of `xs` through an open path. Unneeded arguments are trimmed.

# Examples
```jldoctest
m = @model (n, k) begin
    β ~ Gamma()
    α ~ Gamma()
    θ ~ Beta(α, β)
    x ~ Binomial(n, θ)
    z ~ Binomial(k, α / (α + β))
end;
predictive(m, :θ)

# output
@model (n, θ) begin
        x ~ Binomial(n, θ)
    end
```

"""
predictive(m::Model, xs...) = _predictive(m, namedtuple(xs)(xs))
# predictive(m::Model, xs...) = after(m, xs..., strict = true)

@generated function _predictive(m::Model, xs)
    return after(type2model(m), getntkeys(xs)...; strict=true)
end

export Do
"""
    Do(m, xs...)

Returns a model transformed by adding `xs...` to arguments. The remainder of the body remains the same, consistent with Judea Pearl's "Do" operator. Unneeded arguments are trimmed.

# Examples
```jldoctest
m = @model (n, k) begin
    β ~ Gamma()
    α ~ Gamma()
    θ ~ Beta(α, β)
    x ~ Binomial(n, θ)
    z ~ Binomial(k, α / (α + β))
end;
Do(m, :θ)

# output
@model (n, k, θ) begin
        β ~ Gamma()
        α ~ Gamma()
        x ~ Binomial(n, θ)
        z ~ Binomial(k, α / (α + β))
    end
```

"""
Do(m::Model, xs...) = after(m, xs..., strict = false)


text:
using CFTime
using Dates

# only for daily scale
function dates_miss(dates)
    date_begin = first(dates)
    date_end = last(dates)
    dates_full = date_begin:Dates.Day(1):date_end
    setdiff(dates_full, dates)
end

# only for daily scale
function dates_nmiss(dates)
    date_begin = first(dates)
    date_end = last(dates)

    n_full = (date_end - date_begin) / convert(Dates.Millisecond, Dates.Day(1)) + 1 |> Int
    n_full - length(dates) # n_miss
end

year = Dates.year
month = Dates.month
day = Dates.day

Year = Dates.Year
Month = Dates.Month
Day = Dates.Day

make_datetime = DateTime
make_date = DateTime

export dates_miss, dates_nmiss,
    year, month, day,
    Year, Month, Day,
    make_datetime, make_date


text:
import FractionalCalculus.FracDiffAlg

"""

Grünwald–Letnikov sense fractional derivative algorithms, please refer to [Grünwald–Letnikov derivative](https://en.wikipedia.org/wiki/Gr%C3%BCnwald%E2%80%93Letnikov_derivative) for more details
"""
abstract type GL <: FracDiffAlg end

"""
# Grünwald–Letnikov sense fractional dervivative.

    fracdiff(f, α, start_point, end_point, GL_Direct())

### Example:

```julia-repl
julia> fracdiff(x->x^5, 0, 0.5, GL_Direct())
```
!!! info "Scope"
    Please note Grunwald-Letnikov sense fracdiff only support ``0 < \\alpha < 1``.

Please refer to [Grünwald–Letnikov derivative](https://en.wikipedia.org/wiki/Gr%C3%BCnwald%E2%80%93Letnikov_derivative) for more details.

Grunwald Letnikov direct compute method to obtain fractional derivative, precision are guaranteed but cause more memory allocation and compilation time.
"""
struct GL_Direct <: GL end

"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f, α, end_point, h, GL_Multiplicative_Additive())

Grünwald–Letnikov multiplication-addition-multiplication-addition··· method to approximate fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, GL_Multiplicative_Additive())
1.127403405642918
```

!!! danger "Inaccurate"
    The `GL_Multiplicative_Additive` method is not accruate, please use it at your own risk.

```tex
@book{oldham_spanier_1984,
title={The fractional calculus: Theory and applications of differentiation and integration to arbitrary order},
author={Oldham, Keith B. and Spanier, Jerome},
year={1984}}
```
"""
struct GL_Multiplicative_Additive <: GL end

"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f, α, end_point, h, GL_Lagrange_Three_Point_Interp())

Using Lagrange three poitns interpolation to approximate the fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, GL_Lagrange_Three_Point_Interp())
1.1283963376811044
```

!!! danger "Inaccurate"
    The `GL_Lagrange_Three_Point_Interp` method is not accruate, please use it at your own risk.
"""
struct GL_Lagrange_Three_Point_Interp <: GL end

"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f::Union{Function, Number}, α::AbstractArray, end_point, h, ::GL_Finite_Difference)::Vector

Use finite difference method to obtain Grünwald Letnikov sense fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.01, GL_Finite_Difference())
1.1269695801851276
```

!!! danger "Inaccurate"
    `GL_Finite_Difference` method is not accruate, please use it at your own risk.
"""
struct GL_Finite_Difference <: GL end

"""
fracdiff(f, α, point, h, GL_High_Precision())

Use the high precision algorithms to compute the Grunwald letnikov fractional derivative.

!!! note
    The value interval passing in the function should be a array!
"""
struct GL_High_Precision <: GL end

################################################################
###                    Type defination done                  ###
################################################################

#=
Grunwald Letnikov direct method
=#
function fracdiff(f::Union{Function, Number}, α::Float64, start_point, end_point, ::GL_Direct)
    #checks(f, α, start_point, end_point)

    g(τ) = (f(end_point)-f(τ)) ./ (end_point - τ) .^ (1+α)
    result = f(end_point)/(gamma(1-α) .* (end_point - start_point) .^ α) .+ quadgk(g, start_point, end_point) .* α ./ gamma(1-α)
    return result
end

#TODO: Use the improved alg!! This algorithm is not accurate
#This algorithm is not good, still more to do
function fracdiff(f::Union{Function, Number}, α, end_point, h, ::GL_Multiplicative_Additive)::Float64

    summation = zero(Float64)
    n = Int64(floor(end_point/h))

    for i ∈ 0:n-1
        summation += gamma(i-α)/gamma(i+1)*f(end_point-i*h)
    end
    result = summation*end_point^(-α)*n^α/gamma(-α)

    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::GL_Multiplicative_Additive)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, GL_Multiplicative_Additive()))
    end
    return ResultArray
end

#TODO: This algorithm is same with the above one, not accurate!!!
#This algorithm is not good, still more to do
function fracdiff(f::Union{Function, Number}, α::Float64, end_point, h, ::GL_Lagrange_Three_Point_Interp)::Float64
    #checks(f, α, 0, end_point)

    n = Int64(floor(end_point/h))
    summation = zero(Float64)

    for i ∈ 0:n-1
        summation += gamma(i-α)/gamma(i+1)*(f(end_point-i*h)+1/4*α*(f(end_point-(i-1)*h)-f(end_point-(i+1)*h))+1/8*α^2*(f(end_point-(i-1)*h)-2*f(end_point-i*h)+f(end_point-(i+1)*h)))
    end

    result = summation*end_point^(-α)*n^α/gamma(-α)

    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::GL_Lagrange_Three_Point_Interp)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, GL_Lagrange_Three_Point_Interp()))
    end
    return ResultArray
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::Real, h, ::GL_Finite_Difference)::Float64

    n = Int64(floor(end_point/h))
    result = zero(Float64)

    @fastmath @simd for i ∈ 0:n
        result += (-1)^i/(gamma(i+1)*gamma(α-i+1))*f(end_point-i*h)
    end

    result1=result*h^(-α)*gamma(α+1)
    return result1
end

function fracdiff(f::Union{Function, Number}, α::AbstractArray, end_point, h, ::GL_Finite_Difference)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, GL_Finite_Difference()))
    end
    return ResultArray
end

function fracdiff(f, α, t, p, ::GL_High_Precision)
    if isa(f, Function)
        y=f.(t)
    elseif isa(f, Vector)
        y=f[:]
    end
    h = t[2] - t[1]

    t = t[:]
    n = length(t)
    u = 0
    du = 0
    r = collect(0:p)*h
    R = reverse(Vandermonde(r), dims=2)
    c = inv(R)*y[1:p+1]
    for i = 1:p+1
        u = u.+c[i]*t.^(i-1)
        du = du.+c[i]*t.^(i-1-α)*gamma(i)/gamma(i-α)
    end

    v = y.-u
    g = genfun(p)
    w = getvec(α, n, g)
    dv = zeros(n)
    for i=1:n
        dv[i] = w[1:i]'*v[i:-1:1]/h^α
    end

    dy = dv .+ du
    if abs(y[1])<1e-10
        dy[1]=0
    end
    return dy
end

"""
P-th precision polynomial generate function

```math
g_p(z)=\\sum_{k=1}^p \\frac{1}{k}(1-z)^k
```
"""
function genfun(p)
    a=collect(1:p+1)
    A=Vandermonde(a)'
    return (1 .-a')*inv(A')
end

function getvec(α, n, g)
    p = length(g)-1
    b = 1 + α
    g0 = g[1]
    w = Float64[]
    push!(w, g[1]^α)

    for m = 2:p
        M = m-1
        dA = b/M
        temp = (-(g[2:m] .*collect((1-dA):-dA:(1-b))))' *w[M:-1:1]/g0
        push!(w, temp)
    end

    for k = p+1:n
        M = k-1
        dA = b/M
        temp = (-(g[2:(p+1)] .*collect((1-dA):-dA:(1-p*dA))))' *w[M:-1:(k-p)]./g0 #FIXME: When p is greater than 2 threw DimensionMismatch
        push!(w, temp)
    end

    return w
end


text:
## Display code

"""
Access SymPy's docstrings

There are several calling styles, as finding the underlying SymPy object from
a Julia object is a bit tricky.

Examples
```
@vars x
import Base.Docs.doc
doc(sin(x))        #
doc(sympy[:sin])   # explicit module lookup
doc(SymPy.mpmath[:hypercomb]) # explicit module lookup
doc(Poly(x^2,x), :coeffs) # coeffs is an object method of the poly instance
doc([x 1;1 x], :LUsolve)  # LUsolve is a matrix method
```
"""
Base.Docs.doc(x::SymbolicObject) = Base.Docs.doc(PyObject(x))
Base.Docs.doc(x::SymbolicObject, s::Symbol) = Base.Docs.doc(PyObject(x)[s])
Base.Docs.doc(x::Array{T,N}, s::Symbol) where {T <: SymbolicObject, N} = Base.Docs.doc(PyObject(x)[s])

## Add some of SymPy's displays
## Some pretty printing

#doc(x::SymbolicObject) = print(x[:__doc__])

"Map a symbolic object to a string"
_str(s::SymbolicObject) = s[:__str__]()

"Map an array of symbolic objects to a string"
_str(a::AbstractArray{SymbolicObject}) = map(_str, a)

"call SymPy's pretty print"
pprint(s::SymbolicObject, args...; kwargs...) = sympy_meth(:pprint, s, args...; kwargs...)

"Call SymPy's `latex` function. Not exported. "
latex(s::SymbolicObject, args...; kwargs...)  = sympy_meth(:latex, s, args...; kwargs...)

"create basic printed output"
function jprint(x::SymbolicObject)
    out = PyCall.pycall(pybuiltin("str"), String, PyObject(x))
    if ismatch(r"\*\*", out)
        out = replace(out, "**", "^")
    end
    out
end
jprint(x::AbstractArray) = map(jprint, x)

## show is called in printing tuples, ...
## we would like to use pprint here, but it does a poor job on complicated multi-line expressions
Base.show(io::IO, s::Sym) = print(io, jprint(s))

## We add show methods for the REPL (text/plain) and IJulia (text/latex)

## text/plain
show(io::IO, ::MIME"text/plain", s::SymbolicObject) =  print(io, sympy["pretty"](s))

show(io::IO, ::MIME"text/latex", x::Sym) = print(io, latex(x, mode="equation*", itex=true))
function  show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
    function toeqnarray(x::Vector{Sym})
        a = join([latex(x[i]) for i in 1:length(x)], "\\\\")
        "\\begin{bmatrix}$a\\end{bmatrix}"
    end
    function toeqnarray(x::AbstractArray{Sym,2})
        sz = size(x)
        a = join([join(map(latex, x[i,:]), "&") for i in 1:sz[1]], "\\\\")
        "\\begin{bmatrix}$a\\end{bmatrix}"
    end
    print(io, toeqnarray(x))
end

## Pretty print dicts
function show(io::IO, ::MIME"text/latex", d::Dict{T,S}) where {T<:Sym, S<:Any}
    Latex(x::Sym) = latex(x)
    Latex(x) = sprint(Base.showcompact, x)

    out = "\\begin{equation*}\\begin{cases}"
    for (k,v) in d
        out = out * Latex(k) * " & \\text{=>} &" * Latex(v) * "\\\\"
    end
    out = out * "\\end{cases}\\end{equation*}"
    print(io, out)
end

## Convert SymPy symbol to Julia expression
convert(::Type{Expr}, x::SymbolicObject) = parse(jprint(x))