text:
#### used in batch_graphs.jl.

mutable struct GISSMOEntryType{}
    #
    entry::String
    molecule_name::String
    eDATA_URL::String # includes FID data, cs, J-coupling values.
    STAR_URL::String # contains concentration info.

    #κ_range
    # each row corresponds to a spin_group.
    # each consecutive pairs of columns correspond to the start
    #   and end of an interval for an unique κ.
    # a negative number means the entire spin group has one κ.
    # see cs_config.txt to help visualize κ_range.

    #ι_range::Vector{Tuple{Float64,Float64}}
    #N_protons::Int
end
# TODO update graphs and batch_graphs to pull concentration data from STAR_URL.
#  near Sample_component.Sample_ID

function getGISSMOentriesall()::Vector{GISSMOEntryType}

    entries = Vector{GISSMOEntryType}(undef, 0)

    tmp = GISSMOEntryType("bmse000795",
        "DSS",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000795/simulation_1/bmse000795_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000795/simulation_1/bmse000795-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000051",
        "L-Tyrosine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000051/simulation_1/bmse000051_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000051/simulation_1/bmse000051-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000042",
        "L-Leucine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000042/simulation_1/bmse000042_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000042/simulation_1/bmse000042-simulation_1.str")
    push!(entries, tmp)

    # tmp = GISSMOEntryType("bmse000089",
    #     "Glycine",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000089/simulation_1/bmse000089_simulation_1_nmredata.zip",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000089/simulation_1/bmse000089-simulation_1.str")
    # push!(entries, tmp)
    tmp = GISSMOEntryType("bmse000977",
        "Glycine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000977/simulation_1/bmse000977_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000977/simulation_1/bmse000977-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000041",
        "L-Isoleucine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000041/simulation_1/bmse000041_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000041/simulation_1/bmse000041-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000860",
        "L-Valine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000860/simulation_1/bmse000860_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000860/simulation_1/bmse000860-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000050",
        "L-Tryptophan",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000050/simulation_1/bmse000050_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000050/simulation_1/bmse000050-simulation_1.str")
    push!(entries, tmp)

    # tmp = GISSMOEntryType("bmse000049",
    #     "L-Threonine",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000049/simulation_1/bmse000049_simulation_1_nmredata.zip",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000049/simulation_1/bmse000049-simulation_1.str",
    #     [1.0 2.5 3.8 3.9 4.5 ])
    # push!(entries, tmp)
    tmp = GISSMOEntryType("bmse000859",
        "L-Threonine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000859/simulation_1/bmse000859_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000859/simulation_1/bmse000859-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000867",
        "L-Serine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000867/simulation_1/bmse000867_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000867/simulation_1/bmse000867-simulation_1.str")
    push!(entries, tmp)
    # [2/3; 1.0] # κ.

    tmp = GISSMOEntryType("bmse000900",
        "L-Phenylalanine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000900/simulation_1/bmse000900_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000900/simulation_1/bmse000900-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000044",
        "L-Methionine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000044/simulation_1/bmse000044_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000044/simulation_1/bmse000044-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000038",
        "L-Glutamine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000038/simulation_1/bmse000038_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000038/simulation_1/bmse000038-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000208",
        "L-(+) Lactic acid",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000208/simulation_1/bmse000208_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000208/simulation_1/bmse000208-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000015",
        "D-(+)-Glucose",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000015/simulation_1/bmse000015_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000015/simulation_1/bmse000015-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000206",
        "Caffeine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000206/simulation_1/bmse000206_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000206/simulation_1/bmse000206-simulation_1.str")
    push!(entries, tmp)

    # tmp = GISSMOEntryType("bmse000039",
    #     "L-Histidine",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000039/simulation_1/bmse000039_simulation_1_nmredata.zip",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000039/simulation_1/bmse000039-simulation_1.str",
    #     ((-1.2), (-1.2)),
    #     [(3.0, 3.35); (3.8, 4.17); (6.83, 8.2)] )
    # push!(entries, tmp)
    tmp = GISSMOEntryType("bmse000976",
        "L-Histidine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000976/simulation_1/bmse000976_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000976/simulation_1/bmse000976-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000947",
        "L-Proline",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000947/simulation_1/bmse000947_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000947/simulation_1/bmse000947-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000953", # 600 MHz.
        "Choline",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000953/simulation_1/bmse000953_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000953/simulation_1/bmse000953-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000912",
        "L-Asparagine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000912/simulation_1/bmse000912_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000912/simulation_1/bmse000912-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000875",
        "L-Aspartic acid",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000875/simulation_1/bmse000875_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000875/simulation_1/bmse000875-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000914",
        "L-Lysine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000914/simulation_1/bmse000914_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000914/simulation_1/bmse000914-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000028",
        "L-Alanine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000028/simulation_1/bmse000028_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000028/simulation_1/bmse000028-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000899",
        "L-Arginine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000899/simulation_1/bmse000899_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000899/simulation_1/bmse000899-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000299",
        "Folate",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000299/simulation_1/bmse000299_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000299/simulation_1/bmse000299-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000975",
        "L-Cysteine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000975/simulation_1/bmse000975_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000975/simulation_1/bmse000975-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000109",
        "Putrescine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000109/simulation_1/bmse000109_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000109/simulation_1/bmse000109-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000183",
        "Succinic acid",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000183/simulation_1/bmse000183_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000183/simulation_1/bmse000183-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000454",
        "Purine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000454/simulation_1/bmse000454_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000454/simulation_1/bmse000454-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000913",
        "L-Glutamic acid",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000913/simulation_1/bmse000913_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000913/simulation_1/bmse000913-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000976",
        "Ethanol",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000297/simulation_1/bmse000297_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000297/simulation_1/bmse000297-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000316",
        "Epinephrine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000316/simulation_1/bmse000316_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000316/simulation_1/bmse000316-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000404",
        "(-)-Norepinephrine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000404/simulation_1/bmse000404_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000404/simulation_1/bmse000404-simulation_1.str")
    push!(entries, tmp)

    # L-DOPA
    tmp = GISSMOEntryType("bmse000322",
        "3,4-Dihydroxy-L-phenylalanine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000322/simulation_1/bmse000322_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000322/simulation_1/bmse000322-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000933",
        "Dopamine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000933/simulation_1/bmse000933_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000933/simulation_1/bmse000933-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000457",
        "5-Hydroxy-L-tryptophan",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000457/simulation_1/bmse000457_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000457/simulation_1/bmse000457-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000757",
        "Serotonin",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000757/simulation_1/bmse000757_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000757/simulation_1/bmse000757-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000172",
        "L-Kynurenine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000172/simulation_1/bmse000172_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000172/simulation_1/bmse000172-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000207",
        "Tryptamine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000207/simulation_1/bmse000207_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000207/simulation_1/bmse000207-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000950",
        "Creatine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000950/simulation_1/bmse000950_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000950/simulation_1/bmse000950-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000155",
        "Creatinine",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000155/simulation_1/bmse000155_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000155/simulation_1/bmse000155-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000901",
        "Myo-Inositol",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000901/simulation_1/bmse000901_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000901/simulation_1/bmse000901-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000856",
        "Glycerol",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000856/simulation_1/bmse000856_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000856/simulation_1/bmse000856-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000871",
        "Gamma-Aminobutyric acid",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000871/simulation_1/bmse000871_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000871/simulation_1/bmse000871-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000905",
        "3-Hydroxybutyrate",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000905/simulation_1/bmse000905_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000905/simulation_1/bmse000905-simulation_1.str")
    push!(entries, tmp)

    tmp = GISSMOEntryType("bmse000877",
        "HEPES",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000877/simulation_1/bmse000877_simulation_1_nmredata.zip",
        "http://gissmo.nmrfam.wisc.edu/entry/bmse000877/simulation_1/bmse000877-simulation_1.str")
    push!(entries, tmp)

    # # uses TMS!
    # tmp = GISSMOEntryType("bmse000923",
    #     "L-Thyroxine",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000923/simulation_1/bmse000923_simulation_1_nmredata.zip",
    #     "http://gissmo.nmrfam.wisc.edu/entry/bmse000923/simulation_1/bmse000923-simulation_1.str")
    # push!(entries, tmp)

    return entries
end

function getGISSMOentriesDMEM()::Vector{GISSMOEntryType}

    molecule_names = Vector{String}(undef, 0)
    push!(molecule_names, "L-Tyrosine")
    push!(molecule_names, "L-Leucine")
    push!(molecule_names, "Glycine")
    push!(molecule_names, "L-Isoleucine")
    push!(molecule_names, "L-Valine")
    push!(molecule_names, "L-Tryptophan")
    push!(molecule_names, "L-Threonine")
    push!(molecule_names, "L-Serine")
    push!(molecule_names, "L-Phenylalanine")
    push!(molecule_names, "L-Methionine")
    push!(molecule_names, "L-Glutamine")
    push!(molecule_names, "D-(+)-Glucose")
    push!(molecule_names, "L-Histidine")
    push!(molecule_names, "Choline")
    push!(molecule_names, "L-Lysine")
    push!(molecule_names, "L-Arginine")

    return getGISSMOentries(molecule_names)
end

function extractfields(X::Vector{GISSMOEntryType}, z::String)::Vector{String}

    return collect( getproperty(X[i], Symbol(z)) for i = 1:length(X))
end
# example usage:
# DMEM_entries = getGISSMOentriesDMEM()
# A = extractfields(DMEM_entries, "molecule_name")
# B = extractfields(DMEM_entries, "entry")
# display([A B])

function exportallGISSMOnamestofile(save_path)
    entries = getGISSMOentriesall()
    #molecule_names = collect( entries[i].molecule_name for i = 1:length(entries))
    molecule_names = extractfields(entries, "molecule_name")

    sort!(molecule_names)

    DelimitedFiles.writedlm(save_path, molecule_names)
end

function getGISSMOentry(target_name::String)::GISSMOEntryType

    entries = getGISSMOentriesall()
    #molecule_names = collect( entries[i].molecule_name for i = 1:length(entries))
    molecule_names = extractfields(entries, "molecule_name")

    inds, _ = searchmoleculelist(molecule_names, target_name)

    k = 0
    if length(inds) > 1
        k = inds[1]
        message_string = "Warning, $(target_name) found multiple entries. Using the first one: entry $(entries[k].entry)."
        println(message_string)

    elseif isempty(inds)
        message_string = "Warning, $(target_name) did not find any entries."
        println(message_string)

        return GISSMOEntryType("", "", "", "")

    else
        k = inds[1]
    end

    return entries[k]
end

function getGISSMOentries(target_names::Vector{String})::Vector{GISSMOEntryType}

    entries = collect( getGISSMOentry(target_names[i]) for i = 1:length(target_names) )

    return entries
end

##### for batch_graphs.jl.

function downloadGISSMOentries(entries::Vector{GISSMOEntryType},
                            base_dir::String,
                            save_dir::String;
                            δ_lb::T = 0.1,
                            δ_ub::T = 0.1,
                            unique_cs_tol::T = 1e-6,
                            delay_time = 0.2,
                            default_cs_delta::T = 0.1) where T <: Real
    #
    N = length(entries)

    out_string = Vector{String}(undef, 0)

    for i = 1:N

        print_string = "$(i): $(entries[i].molecule_name)"
        println(print_string)

        sleep(delay_time) # replace with retry later.
        css_sys, cs_singlets_compact = downloadGISSMOentry(entries[i],
            base_dir, save_dir;
                δ_lb = δ_lb,
                δ_ub = δ_ub,
                unique_cs_tol = unique_cs_tol)

        println()

        addentrydefaultconfigtobuffer!(out_string,
                entries[i].molecule_name,
                css_sys,
                cs_singlets_compact;
                default_cs_delta = default_cs_delta)
    end

    # special cases: add D2O.
    entry = "0"
    molecule_name = "D2O"
    singlet_ppm = 4.8
    css_sys, cs_singlets_compact = specialsingletentry(entry,
        molecule_name,
        base_dir,
        save_dir,
        singlet_ppm;
            δ_lb = δ_lb,
            δ_ub = δ_ub,
            unique_cs_tol = unique_cs_tol)

    # add to buffer.
    addentrydefaultconfigtobuffer!(out_string,
            molecule_name,
            css_sys,
            cs_singlets_compact;
            default_cs_delta = default_cs_delta)

    # write buffer to file.
    file_path = joinpath(save_dir, "default_cs_config.txt")
    DelimitedFiles.writedlm(file_path, out_string)

    return nothing
end

# zero_tol_sigdigits and unique_cs_tol are used to determine unique cs with tolerance.
# base_dir is where the GISSMO entry will be downloaded to.
# save_dir is where the cs and J-value JLD files will be saved.
function downloadGISSMOentry(X::GISSMOEntryType,
                            base_dir::String,
                            save_dir::String;
                            δ_lb::T = 0.1,
                            δ_ub::T = 0.1,
                            unique_cs_tol::T = 1e-6,
                            zero_tol_sigdigits = 6) where T <: Real

    # parse.
    molecule_name = X.molecule_name
    entry = X.entry
    eDATA_URL = X.eDATA_URL
    STAR_URL = X.STAR_URL

    # storage.
    record_name = "$(molecule_name)_$(entry)"
    data_dir = joinpath(base_dir, record_name)
    t = @task begin; isdir(data_dir) || mkdir(data_dir); end
    schedule(t); wait(t)

    # Download the STAR file.
    file_name = "$(record_name).str"
    dest = joinpath(data_dir, file_name)
    t = @task begin; isfile(dest) || download(STAR_URL, dest); end
    schedule(t); wait(t)

    # Download zip file.
    file_name = "$(record_name).zip"
    dest = joinpath(data_dir, file_name)
    t = @task begin; isfile(dest) || download(eDATA_URL, dest); end
    schedule(t); wait(t)

    #q = filter(x->endswith(x, ".zip"), readdir(data_dir))

    # unzip.
    InfoZIP.unzip(dest, data_dir)

    # remove zip file.
    #rm(dest)

    # get SDF filename.
    tmp = filter(isdir, readdir(data_dir; join=true))
    @assert length(tmp) == 1 # should be only one folder, which is the unzipped one.
    folder_path = tmp[1]

    tmp = filter(x->endswith(x, ".sdf"), readdir(folder_path))
    @assert length(tmp) == 1 # should be exactly 1 SDF file.
    SDF_name = tmp[1]
    SDF_path = joinpath(folder_path, SDF_name)

    file_strings = readlines(SDF_path)

    #### get absolute concentrations.
    solubility_gL_solute = solubilityinwater(molecule_name)
    molar_mass_solute = getmolarmass(base_dir, molecule_name)

    compound_names, concentration_values, concentration_units,
    compound_types, solute_concentration_string,
        solute_concentration_unit = parsecompositionfromNMReDATA(file_strings)

    #
    solute_mM, ref_mM, molar_mass_ref, solubility_gL_ref,
    ref_molecule_name = convertsampleconcentration(compound_names, concentration_values, concentration_units,
                            compound_types,
                            solute_concentration_string,
                            solute_concentration_unit,
                            molar_mass_solute, solubility_gL_solute)

    # assemble output. Units in mM<.
    println("solute_concentration = ", solute_mM)
    println("reference_concentration = ", ref_mM)

    #### find pH info.
    pH, temperature = parseconditionsfromNMReDATA(readlines(SDF_path))

    #### find the section on chemical shifts.
    cs_string = "> <NMREDATA_ASSIGNMENT>"

    H_IDs, H_css = parsecssfromNMReDATA(file_strings)

    #J_string = "> <NMREDATA_J>"
    J_IDs, J_vals = parseJsfromNMReDATA(file_strings)

    ### attempt on auto-grouping.

    J_IDs_sys, J_vals_sys, H_IDs_sys = partitionJcouplings(J_IDs, J_vals)

    css_sys = getcssys(H_IDs_sys, H_IDs, H_css)

    𝐽_IDs_sys = getJsys(J_IDs_sys, J_vals_sys, H_IDs_sys)

    H_singlets, cs_singlets,
        cs_singlets_compact = createsingletsysems(H_IDs, H_IDs_sys,
                                    H_css;
                                    zero_tol = unique_cs_tol)
    #
    ### remove spin groups that only have one unique chem shift, but have J-coupling between the atoms.
    # Each of such groups are are really one singlet group.
    removeisolatedcs!(css_sys,
                    cs_singlets,
                    cs_singlets_compact,
                    H_singlets,
                    H_IDs_sys,
                    J_IDs_sys,
                    J_vals_sys,
                    𝐽_IDs_sys;
                    zero_tol_sigdigits = zero_tol_sigdigits)
    #
    println("Entry: ", entry)
    println("cs_singlets: ", cs_singlets)
    println("cs_singlets_compact: ", cs_singlets_compact)
    println("H_singlets: ", H_singlets)

    println("css_sys: ", css_sys)
    println("J_vals_sys: ", J_vals_sys)
    println("𝐽_IDs_sys: ", 𝐽_IDs_sys)
    println()
    #### end code for moving spin group to singlet group.

    #
    cs_LUT, p_cs_sys = constructLUTcss(css_sys; zero_tol = unique_cs_tol)

    cs_len_sys = getcslengthfromLUT(cs_LUT)

    ### cs and J-values upper and lower bounds.
    perturblbfunc = xx->additiveperturbation(xx, δ_lb)
    perturbubfunc = xx->additiveperturbation(xx, δ_ub)

    #
    J_lb_sys = getJperturbed(J_vals_sys, perturblbfunc)
    J_ub_sys = getJperturbed(J_vals_sys, perturbubfunc)

    cs_lb_sys = getJperturbed(css_sys, perturblbfunc)
    cs_ub_sys = getJperturbed(css_sys, perturbubfunc)

    # store.
    #save_path = "$(save_dir)/$(record_name).jld"
    save_path = "$(save_dir)$(record_name).jld"

    JLD.save(save_path,
        "css_sys", css_sys,
        "J_IDs_sys", 𝐽_IDs_sys,
        "J_vals_sys", J_vals_sys,
        "cs_singlets", cs_singlets,

        "J_lb_sys", J_lb_sys,
        "J_ub_sys", J_ub_sys,
        "cs_lb_sys", cs_lb_sys,
        "cs_ub_sys", cs_ub_sys,
        "cs_LUT", cs_LUT,

        "ref_mM", ref_mM,
        "solute_mM", solute_mM,
        "solubility_gL_solute", solubility_gL_solute,
        "solubility_gL_ref", solubility_gL_ref,
        "molar_mass_ref", molar_mass_ref,
        "molar_mass_solute", molar_mass_solute,
        "pH", pH,
        "temperature", temperature,
        "ref_molecule_name", ref_molecule_name)

    return css_sys, cs_singlets_compact
end

function specialsingletentry(  entry,
                    molecule_name,
                    base_dir::String,
                    save_dir::String,
                    singlet_ppm;
                    unique_cs_tol = 1e-6,
                    δ_lb = 0.1,
                    δ_ub = 0.1)
    #
    H_IDs = Vector{Int}(undef, 1)
    H_IDs[1] = 1

    H_css = Vector{Float64}(undef, 1)
    H_css[1] = singlet_ppm # we're setting 4.8 ppm for D2O.

    J_IDs = Vector{Tuple{Int,Int}}(undef, 0)

    J_vals = Vector{Float64}(undef, 0)

    J_IDs_sys, J_vals_sys, H_IDs_sys = partitionJcouplings(J_IDs, J_vals)

    css_sys = getcssys(H_IDs_sys, H_IDs, H_css)

    𝐽_IDs_sys = getJsys(J_IDs_sys, J_vals_sys, H_IDs_sys)

    H_singlets, cs_singlets,
        cs_singlets_compact = createsingletsysems(H_IDs, H_IDs_sys,
                                    H_css;
                                    zero_tol = unique_cs_tol)

    cs_LUT, p_cs_sys = constructLUTcss(css_sys; zero_tol = unique_cs_tol)

    cs_len_sys = getcslengthfromLUT(cs_LUT)

    perturblbfunc = xx->additiveperturbation(xx, δ_lb)
    perturbubfunc = xx->additiveperturbation(xx, δ_ub)

    J_lb_sys = getJperturbed(J_vals_sys, perturblbfunc)
    J_ub_sys = getJperturbed(J_vals_sys, perturbubfunc)

    cs_lb_sys = getJperturbed(css_sys, perturblbfunc)
    cs_ub_sys = getJperturbed(css_sys, perturbubfunc)

    record_name = "$(molecule_name)_$(entry)"
    save_path = "$(save_dir)$(record_name).jld"

    JLD.save(save_path,
        "css_sys", css_sys,
        "J_IDs_sys", 𝐽_IDs_sys,
        "J_vals_sys", J_vals_sys,
        "cs_singlets", cs_singlets,

        "J_lb_sys", J_lb_sys,
        "J_ub_sys", J_ub_sys,
        "cs_lb_sys", cs_lb_sys,
        "cs_ub_sys", cs_ub_sys,
        "cs_LUT", cs_LUT,

        "ref_mM", NaN,
        "solute_mM", NaN,
        "solubility_gL_solute", NaN,
        "solubility_gL_ref", NaN,
        "molar_mass_ref", NaN,
        "molar_mass_solute", NaN,
        "pH", NaN,
        "temperature", NaN,
        "ref_molecule_name", NaN)

    return css_sys, cs_singlets_compact
end

function addentrydefaultconfigtobuffer!(out_string::Vector{String},
                                        molecule_name::String,
                                        css_sys::Vector{Vector{T}},
                                        cs_singlets_compact::Vector{T};
                                        default_cs_delta = 0.1) where T <: Real
    #
    # save heading to template.
    push!(out_string, molecule_name)

    # save the default cs configuration info.
    N_spin_sys = length(css_sys)
    for i = 1:N_spin_sys
        push!(out_string, "Group $(i)")

        for k = 1:length(css_sys[i])
            line_string = "cs $(css_sys[i][k]): $(default_cs_delta), $(k)"
            push!(out_string, line_string)
        end

        push!(out_string, "end group")
    end

    for i = 1:length(cs_singlets_compact)
        push!(out_string, "Group $(i+N_spin_sys)")
        line_string = "cs $(cs_singlets_compact[i]): $(default_cs_delta), 1"
        push!(out_string, line_string)
        push!(out_string, "end group")
    end

    push!(out_string, "end compound")
    push!(out_string, "")

    return nothing
end

function getGISSMOexperimentpaths(record_names, record_entries, dir_list)

    exp_paths = Vector{String}(undef, length(record_names))
    for k = 1:length(record_names)

        inds = findall(xx->occursin("$(record_names[k])_$(record_entries[k])", xx), dir_list)
        @assert length(inds) > 0

        current_folder_path = dir_list[inds[1]]
        next_folder = filter(isdir, readdir(current_folder_path; join = true))[1]

        while !isempty(next_folder) && !occursin("pdata", next_folder)

            current_folder_path = next_folder
            next_folder = filter(isdir, readdir(current_folder_path; join = true))[1]
        end

        exp_paths[k] = current_folder_path
    end

    return exp_paths
end


text:
# template log-on and fetch graph

using Caesar

include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
session = "SESSHAUVDEV2"
addrdict["session"] = session
cg, addrdict = standardcloudgraphsetup(addrdict=addrdict)

@time all = ls(cg, "SESSTURT21")

@time dd = ls(cg, session, sym=:x1)

dd2 = ls(cg, session, exvid=5)

findAllBinaryFactors(cg, session)

using IncrementalInference

fsym = collect(keys(dd))[2]
cv = getCloudVert(cg, session, sym=fsym)

ty = getfnctype(cv)

fieldnames(ty.name)

@show ty.name.name

#

ext: jl


text:
using PooksoftAssetModelingKit
using CSV
using DataFrames
using Test
using Statistics
using KernelDensity
using Dates

function estimate_random_walk_model_test()

    # test setup -
    ticker_symbol =
    path_test_data_file = "./data/Test.csv"

    # load the data file -
    df = CSV.read(path_test_data_file,DataFrame)

    # what dates are we going to look at?
	date_start = Date("2016-01-04")
    date_stop = Date("2020-10-09")

    # filter the date range we want -
	tmp = filter(r->(r.timestamp >= date_start && r.timestamp <= date_stop),df)

    # compute the difference for each ticker -
	result = compute_linear_return_array(tmp)
	if (isa(result.value, Exception) == true)
		return false
	end
    price_return_array = result.value

    # grab the model - is this the right type?
	model = kde_lscv(price_return_array)
    if (isa(model, UnivariateKDE) == false)
        return false
    end

    # we got here - all is ok
    return true
end

function run_gbm_model_simulation()

    # test setup -
    path_test_data_file = "./data/Test.csv"

    # load the data file -
    data_table = CSV.read(path_test_data_file,DataFrame)

    # estimated the average return
    result = compute_linear_return_array(data_table)
    if (isa(result.value,Exception) == true)
        return false
    end
    μ = mean(result.value)

    # estimate the average volatility -
    result = compute_return_volatility(data_table)
    if (isa(result.value,Exception) == true)
        return false
    end
    σ = result.value

    # build the model object -
    result = build_geometric_brownian_motion_model(μ,σ)
    if (isa(result.value,Exception) == true)
        return false
    end
    model = result.value

    # setup and run the simulation -
    tspan = (0.0,14.0)
    stepSize = 0.0001
    initial_condition = 83.81
    number_of_trials = 100

    # run -
    result = evaluate(model,initial_condition,tspan,stepSize;number_of_trials=number_of_trials)
    if (isa(result.value,Exception) == true)
        return false
    end
    simulation_result = result.value

    # check types on the return -
    if (haskey(simulation_result,:X) == true && haskey(simulation_result,:T) == true)
        return true
    end

    # default is to fail (just like everyhting I do in life ...)
    return false
end

@testset "asset_modeling_test_set" begin
    @test run_gbm_model_simulation() == true
    @test estimate_random_walk_model_test() == true
end


text:
include("../CovarianceFunctions/CovarianceFunctions.jl")
include("ExtendKalmanFilter.jl")
# For linear algebra
using LinearAlgebra
# logging
using Logging
using PyPlot

## PDEs type
abstract type AbstractPDEs end
struct NonlinElliptic2d{Tα,Tm,TΩ} <: AbstractPDEs
    # eqn: -Δu + α*u^m = f in [Ω[1,1],Ω[2,1]]*[Ω[1,2],Ω[2,2]]
    α::Tα
    m::Tm
    Ω::TΩ
    bdy::Function
    rhs::Function
end

## sample points
function sample_points_rdm(eqn::NonlinElliptic2d, N_domain, N_boundary)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]

    X_domain = hcat(rand(Float64,(N_domain, 1))*(x1r-x1l).+x1l,rand(Float64,(N_domain, 1))*(x2r-x2l).+x2l)

    N_bd_each=convert(Int64, N_boundary/4)
    if N_boundary != 4 * N_bd_each
        println("[sample points] N_boundary not divided by 4, replaced by ", 4 * N_bd_each)
        N_boundary = 4 * N_bd_each
    end

    X_boundary = zeros((N_boundary, 2))
    # bottom face
    X_boundary[1:N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2l*ones(N_bd_each))
    # right face
    X_boundary[N_bd_each+1:2*N_bd_each, :] = hcat(x1r*ones(N_bd_each),(x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    # top face
    X_boundary[2*N_bd_each+1:3*N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2r*ones(N_bd_each))

    # left face
    X_boundary[3*N_bd_each+1:N_boundary, :] = hcat(x1l*ones(N_bd_each), (x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    return X_domain', X_boundary'
end

function sample_points_grid(eqn::NonlinElliptic2d, h_in, h_bd)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]
    x = x1l + h_in:h_in:x1r-h_in
    y = x2l + h_in:h_in:x2r-h_in
    X_domain = reduce(hcat,[[x[i], y[j]] for i in 1:length(x) for j in 1:length(x)])

    l = length(x1l:h_bd:x1r-h_bd)
    X_boundary = vcat([x1l:h_bd:x1r-h_bd x2l*ones(l)], [x1r*ones(l) x2l:h_bd:x2r-h_bd], [x1r:-h_bd:x1l+h_bd x2r*ones(l)], [x1l*ones(l) x2r:-h_bd:x1l+h_bd])
    return X_domain, X_boundary'
end

# get prior Gram matrices for Θ
function get_Σ0_init(eqn::NonlinElliptic2d, cov::AbstractCovarianceFunction, X_domain, X_boundary)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    Δδ_coefs = -1.0

    # get linearized PDEs correponding measurements
    meas_bdyδ = [PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_intΔδ = [ΔδPointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), Δδ_coefs, 0.0) for i = 1:N_domain]
    meas_intδ = [PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]

    Σ0_init = zeros(2*N_domain+N_boundary,2*N_domain+N_boundary)

    measurements = Vector{Vector{<:AbstractPointMeasurement}}(undef,3)
    measurements[1] = meas_bdyδ; measurements[2] = meas_intδ
    measurements[3] = meas_intΔδ
    cov(Σ0_init, reduce(vcat,measurements))

    return Σ0_init
end

# get_G
function get_G(z,N_boundary,N_domain,α,m,rhs,bdy)
    y = zeros(N_boundary+N_domain)
    y[1:N_boundary] = z[1:N_boundary]-bdy
    y[N_boundary+1:end] = z[N_boundary+N_domain+1:end] + α*z[N_boundary+1:N_boundary+N_domain].^m-rhs
    return y
end

# get dG
function get_dG(z,N_boundary,N_domain,α,m)
    return [LinearAlgebra.I zeros(N_boundary,2*N_domain); zeros(N_domain,N_boundary)  α*m*Diagonal(z[N_boundary+1:N_boundary+N_domain].^(m-1)) LinearAlgebra.I]
end

function get_Σ_η(eqn::NonlinElliptic2d, cov_rhs::AbstractCovarianceFunction, cov_bdy::AbstractCovarianceFunction, X_domain, X_boundary)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    meas_bdδ = [PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_intδ = [PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]

    Σ_η_rhs = zeros(N_domain, N_domain)
    Σ_η_bd = zeros(N_boundary, N_boundary)

    cov_rhs(Σ_η_rhs, meas_intδ)
    cov_bdy(Σ_η_bd, meas_bdδ)

    Σ_η = zeros(N_boundary+N_domain, N_boundary+N_domain)
    @show size(Σ_η), N_boundary, N_domain
    Σ_η[1:N_boundary,1:N_boundary] = Σ_η_bd
    Σ_η[N_boundary+1:end,N_boundary+1:end] = Σ_η_rhs
    return Σ_η
end

α = 100.0
m = 3
Ω = [[0,1] [0,1]]
h_in = 0.02
h_bd = 0.02
lengthscale_u = 0.3
kernel_u = "Matern7half"
cov_u = MaternCovariance7_2(lengthscale_u)
kernel_bdy = "Matern5half"
lengthscale_bdy = 0.3
cov_bdy = MaternCovariance5_2(lengthscale_bdy)
kernel_rhs = "Matern5half"
lengthscale_rhs = 0.3
cov_rhs = MaternCovariance5_2(lengthscale_rhs)

# ground truth solution
freq = 600
s = 4
function fun_u(x)
    ans = 0
    @inbounds for k = 1:freq
        ans += sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
        # H^t norm squared is sum 1/k^{2s-2t}, so in H^{s-1/2}
    end
    return ans
end

# right hand side
function fun_rhs(x)
    ans = 0
    @inbounds for k = 1:freq
        ans += (2*k^2*pi^2)*sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
    end
    return ans + α*fun_u(x)^m
end

# boundary value
function fun_bdy(x)
    return fun_u(x)
end

@info "[solver started] NonlinElliptic2d"
@info "[equation] -Δu + $α u^$m = f"
eqn = NonlinElliptic2d(α,m,Ω,fun_bdy,fun_rhs)
X_domain, X_boundary = sample_points_grid(eqn, h_in, h_bd)
N_domain = size(X_domain,2)
N_boundary = size(X_boundary,2)
@info "[sample points] grid size $h_in"
@info "[sample points] N_domain is $N_domain, N_boundary is $N_boundary"
@info "[prior on u] choose $kernel_u, lengthscale $lengthscale_u\n"
@info "[prior on rhs] choose $kernel_rhs, lengthscale $lengthscale_rhs\n"
@info "[prior on bdy] choose $kernel_bdy, lengthscale $lengthscale_bdy\n"

Σ0_init = get_Σ0_init(eqn, cov_u, X_domain, X_boundary)
Σ_η_init = get_Σ_η(eqn, cov_rhs, cov_bdy, X_domain, X_boundary)
z_init = zeros(2*N_domain+N_boundary) # initial solution
rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]

truth = [fun_u(X_domain[:,i]) for i in 1:N_domain]

N_iter = 20
N_θ = 2*N_domain+N_boundary

r_0, Σ_0 = zeros(Float64, N_θ), Σ0_init
θ0_mean, θθ0_cov = r_0, Σ_0

N_y = N_boundary+N_domain
Σ_η = zeros(N_boundary+N_domain, N_boundary+N_domain)
Σ_η[1:N_boundary,1:N_boundary] = 1e-2*diagm(diag(Σ_0[1:N_boundary,1:N_boundary]))
Σ_η[N_boundary+1:end,N_boundary+1:end] = 1e-2*diagm(diag(Σ_0[N_boundary+N_domain+1:end,N_boundary+N_domain+1:end]))
y = zeros(Float64, N_y)

struct Setup_Param{IT<:Int, FT<:AbstractFloat}
    N_boundary::IT
    N_domain::IT
    α::FT
    m::IT
    rhs::Array{FT,1}
    bdy::Array{FT,1}
    N_θ::IT
    N_y::IT
end

function forward(s_param, θ)
    N_boundary,N_domain,α,m,rhs,bdy = s_param.N_boundary, s_param.N_domain, s_param.α, s_param.m, s_param.rhs, s_param.bdy
    G = get_G(θ,N_boundary,N_domain,α,m,rhs,bdy)
    dG = get_dG(θ,N_boundary,N_domain,α,m)
    return G, dG
end

s_param = Setup_Param(N_boundary, N_domain, α, m, rhs, bdy, N_θ, N_y)

method = "ExKI"
# UKI-1
exki_obj = ExKI_Run(s_param, forward,
method,
y, Σ_η,
r_0,
Σ_0,
θ0_mean, θθ0_cov,
N_iter,)

# @info "ERRORs are :" , norm(exki_obj.θ_mean[end] - θ_post), norm(exki_obj.θθ_cov[end] - Σ_post)
@info norm(truth - exki_obj.θ_mean[end][N_boundary+1:N_boundary+N_domain])/sqrt(N_domain)
# figure()
# plot(truth, "-o", markersize=2, label="Reference")
sol_mean = exki_obj.θ_mean[end][N_boundary+1:N_boundary+N_domain]
sol_std = sqrt.(diag(exki_obj.θθ_cov[end][N_boundary+1:N_boundary+N_domain, N_boundary+1:N_boundary+N_domain]))
# plot(sol_mean, label="Prediction", color="C1")
# plot(sol_mean - sol_std, "--", color="C1")
# plot(sol_mean + sol_std, "--", color="C1")

# legend()

Nh = convert(Int,sqrt(N_domain))
figure()
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(truth,Nh,Nh), label = "Reference")
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_mean - sol_std,Nh,Nh), color="C1", label = "Lower")
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_mean + sol_std,Nh,Nh), color="C1", label = "Upper")
display(gcf())

figure()
contourf(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(abs.(truth-sol_mean),Nh,Nh))
colorbar()
display(gcf())

figure()
contourf(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_std,Nh,Nh))
colorbar()
display(gcf())


text:
global drag_get=Vector{Function}(),
drag_set=Vector{Function}(),
drag_type=Vector{Symbol}()
#= c for common
r for read-only
i for infinite =#
function drag_work(from::Int,to::Int;shift::Bool=false,ctrl::Bool=false)::Bool
	if (drag_type[to]!=:c)||(from==to)
		return false
	end
	f1=deepcopy(drag_get[from](from))::Pair
	f2=deepcopy(drag_get[to](to))::Pair
	if isa(f2.first,EI) # =>空格
		if drag_type[from]==:i
			drag_set[to](to,Pair(f1.first,ctrl ? stack(f1.first) : 0x1))
		else
			if ctrl # 分半
				mi=f1.second>>1
				drag_set[from](from,Pair(mi==0 ? EI() : f1.first,mi))
				drag_set[to](to,Pair(deepcopy(f1.first),f1.second-mi))
			elseif shift # 分一
				if f1.second==0x1
					return false
				end
				drag_set[from](from,Pair(f1.first,0x1))
				drag_set[to](to,Pair(deepcopy(f1.first),f1.second-0x1))
			else # 全移
				drag_set[from](from,Pair(EI(),0x0))
				drag_set[to](to,f1)
			end
		end
	elseif id(f2.first)==id(f2.first) # 同类
		st=stack(f2.first)
		if st==f2.second
			return false
		end
		if drag_type[from]==:i # 填满
			drag_set[to](to,Pair(f2.first,st))
		else
			sum=f1.second+f2.second
			if sum>st
				drag_set[from](from,Pair(f1.first,sum-st))
				drag_set[to](to,Pair(f2.first,st))
			else
				drag_set[from](from,Pair(EI(),0x0))
				drag_set[to](to,Pair(f2.first,sum))
			end
		end
	elseif drag_type[from]==:c
		drag_set[from](from,f2)
		drag_set[to](to,f1)
	else
		return false
	end
	return true
end


text:
using IterTools

cd(dirname(@__FILE__))
lines = open("nine.txt", "r") do f
    readlines(f)
end
values = map(x -> parse(Int, x), lines)

const N = 25
deq = values[1:N]

pair_sums(v) = (i+j for (i,j) in subsets(v,2))

function check_sum(value, deq)
    for _ in (s for s in pair_sums(deq) if s == value)
        return true
    end

    return false
end

weak = 0
len = length(values)
for n in N+1:len
    value = values[n]
    if !check_sum(value, deq)
        global weak = value
        break
    end
    popfirst!(deq)
    push!(deq, value)
end

println(weak)

function acc(init)
    acc = init
    x -> (acc += x; acc)
end

# cumsum() is eager, but we need a lazy version for takewhile()
lazy_cumsum(values) =
    (acc(0)(x) for x in values)

function sum_range(limit, n, m, values)
    sum = 0
    count = 0
    for x in takewhile(x -> x <= limit, lazy_cumsum(values[n:m]))
        sum = x
        count += 1
    end
    (sum == limit, n, n + count - 1)
end

found, i, j = first(filter(x -> x[1] == true, [sum_range(weak, i, len, values) for i in 1:len-1]))
s = values[i:j]
println(s)
println(sum(s))
println(minimum(s) + maximum(s))


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Does not search. Takes an input sequence and replays it, producing emulated logs
"""
module Ref

export RefESParams, RefESResult, mc_search, exprsearch, SearchParams, SearchResult, get_derivtree

using Reexport
using ExprSearch
using RLESUtils, GitUtils, CPUTimeUtils, Observers
@reexport using LinearDerivTrees  #pretty_string
using GrammaticalEvolution
using Iterators

import LinearDerivTrees: get_derivtree
import ..ExprSearch: exprsearch, get_fitness
import Base: isless, copy!

type RefESParams <: SearchParams
  #tree params
  maxsteps::Int64

  #Ref
  actions::Vector{Int64}
  n_samples::Int #output is constant, only to emulate logs

  observer::Observer
end
RefESParams(maxsteps::Int64, actions::Vector{Int64}) = RefESParams(maxsteps, actions, Observer())

type RefESResult <: SearchResult
  tree::LinearDerivTree
  actions::Vector{Int64}
  fitness::Float64
  expr
  best_at_eval::Int64
  totalevals::Int64
end

exprsearch(p::RefESParams, problem::ExprProblem, userargs...) = ref_search(p, problem, userargs...)

get_derivtree(result::RefESResult) = get_derivtree(result.tree)

function ref_search(p::RefESParams, problem::ExprProblem, userargs...)
  @notify_observer(p.observer, "verbose1", ["Starting Ref search"])
  @notify_observer(p.observer, "computeinfo", ["starttime", string(now())])

  grammar = get_grammar(problem)
  tree_params = DerivTreeParams(grammar, p.maxsteps)
  tree = DerivationTree(tree_params)

  tstart = CPUtime_start()
  fitness, expr = playsequence!(tree, problem, p.actions) #tree is modified to match expr
  result = RefESResult(tree, p.actions, fitness, expr, 1, p.n_samples)
  elapsed_cpu_s1 = CPUtime_elapsed_s(tstart) #elapsed time for 1 eval in seconds

  #emulate logs...
  for i = 1:p.n_samples
    @notify_observer(p.observer, "iteration", [i])

    @notify_observer(p.observer, "elapsed_cpu_s", [i, i * elapsed_cpu_s1])
    @notify_observer(p.observer, "current_best", [i, result.fitness, string(result.expr)])
  end

  @notify_observer(p.observer, "result", [result.fitness, string(result.expr), result.best_at_eval, result.totalevals])

  #meta info
  @notify_observer(p.observer, "computeinfo", ["endtime",  string(now())])
  @notify_observer(p.observer, "computeinfo", ["hostname", gethostname()])
  @notify_observer(p.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
  @notify_observer(p.observer, "parameters", ["maxsteps", p.maxsteps])
  @notify_observer(p.observer, "parameters", ["actions", string(p.actions)])
  @notify_observer(p.observer, "parameters", ["totalevals", p.n_samples])

  return result
end

#tree is modified
function playsequence!(tree::DerivationTree, problem::ExprProblem, sequence::Vector{Int64})
  play!(tree, sequence)
  expr = get_expr(tree)
  fitness = get_fitness(problem, expr)

  return fitness, expr
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module CellFreeModelGenerationKit

# include -
include("Include.jl")

# export types -
export VLJuliaModelObject

# parser methods -
export parse_vff_model_document
export parse_vff_metabolic_section
export parse_vff_sequence_section
export parse_vff_bio_types_section
export parse_vff_species_bounds_section
export parse_vff_grn_section
export read_model_document
export generate_default_project

# export methods -
export build_julia_model_object
export generate

# general methods -
export generate_stoichiometric_matrix

# code strategy methods -
export build_data_dictionary_program_component
export build_control_program_component
export build_kinetics_program_component

# advanced: API call back methods -
export reorder_reaction_symbol_array
export reorder_molecular_symbol_array

# export keys -
export ir_master_reaction_table_key
export ir_list_of_molecular_species_key
export ir_list_of_reaction_tags_key
export ir_master_species_bounds_table_key
export ir_sequence_section_table_key

# export Minerva -
export minerva_parser
export minerva_scanner
export minerva_grn_scan_function
export minerva_biological_type_assignment_scan_function
export minerva_bound_type_assignment_scan_function
export minerva_species_bounds_record_parser
export check

end # module

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
@testset "Voigt form" begin
    for dim in 1:3, T in (Float32, Float64)
        A = rand(Tensor{NTuple{2, dim}, T})
        AA = rand(Tensor{NTuple{4, dim}, T})
        A_sym = rand(Tensor{Tuple{@Symmetry{dim, dim}}, T})
        AA_sym = rand(Tensor{NTuple{2, @Symmetry{dim, dim}}, T})
        # tovoigt
        @test (@inferred tovoigt(AA)) * (@inferred tovoigt(A)) ≈ tovoigt(AA ⊡ A)
        @test (@inferred tovoigt(AA_sym)) * (@inferred tovoigt(A_sym, offdiagscale=T(2))) ≈ tovoigt(AA_sym ⊡ A_sym)
        # fromvoigt
        @test (@inferred fromvoigt(SecondOrderTensor{dim}, tovoigt(A)))::SecondOrderTensor{dim, T} ≈ A
        @test (@inferred fromvoigt(FourthOrderTensor{dim}, tovoigt(AA)))::FourthOrderTensor{dim, T} ≈ AA
        @test (@inferred fromvoigt(SymmetricSecondOrderTensor{dim}, tovoigt(A_sym, offdiagscale=T(2)), offdiagscale=T(2)))::SymmetricSecondOrderTensor{dim, T} ≈ A_sym
        @test (@inferred fromvoigt(SymmetricFourthOrderTensor{dim}, tovoigt(AA_sym, offdiagscale=T(2)), offdiagscale=T(2)))::SymmetricFourthOrderTensor{dim, T} ≈ AA_sym
        @test (@inferred fromvoigt(SecondOrderTensor{dim, Float64}, tovoigt(A)))::SecondOrderTensor{dim, Float64} ≈ A
        @test (@inferred fromvoigt(FourthOrderTensor{dim, Float64}, tovoigt(AA)))::FourthOrderTensor{dim, Float64} ≈ AA
        @test (@inferred fromvoigt(SymmetricSecondOrderTensor{dim, Float64}, tovoigt(A_sym, offdiagscale=T(2)), offdiagscale=T(2)))::SymmetricSecondOrderTensor{dim, Float64} ≈ A_sym
        @test (@inferred fromvoigt(SymmetricFourthOrderTensor{dim, Float64}, tovoigt(AA_sym, offdiagscale=T(2)), offdiagscale=T(2)))::SymmetricFourthOrderTensor{dim, Float64} ≈ AA_sym
        # tomandel
        @test (@inferred tomandel(AA_sym)) * (@inferred tomandel(A_sym)) ≈ tomandel(AA_sym ⊡ A_sym)
        # frommandel
        @test (@inferred frommandel(SymmetricSecondOrderTensor{dim}, tomandel(A_sym)))::SymmetricSecondOrderTensor{dim, T} ≈ A_sym
        @test (@inferred frommandel(SymmetricFourthOrderTensor{dim}, tomandel(AA_sym)))::SymmetricFourthOrderTensor{dim, T} ≈ AA_sym
        @test (@inferred frommandel(SymmetricSecondOrderTensor{dim, Float64}, tomandel(A_sym)))::SymmetricSecondOrderTensor{dim, Float64} ≈ A_sym
        @test (@inferred frommandel(SymmetricFourthOrderTensor{dim, Float64}, tomandel(AA_sym)))::SymmetricFourthOrderTensor{dim, Float64} ≈ AA_sym
        # error
        @test_throws Exception Tensorial.default_voigt_order(Val(4))
        @test_throws Exception tovoigt(rand(Mat{4,4}))
    end
end


text:
using Distributed: @everywhere, pmap

@everywhere begin
    using LinearAlgebra: qr!
    using Dates: @dateformat_str, now, format
    using Printf: @sprintf
    using SharedArrays

    push!(LOAD_PATH, @__DIR__)
    using MultiprecisionFloats
    using RKTK
    using DZMisc: approx_norm, orthonormalize_columns!
    using GoldenSectionSearch: golden_section_search

    const AccurateReal = Float64x2
    setprecision(256)
end

@everywhere function log(args...)
    println(args...)
    flush(stdout)
end

@everywhere function approx_coulomb_force_energy(
        points::Matrix{T}, i::Int) where {T <: Real}
    dim, num_points = size(points)
    displ = Vector{Float64}(undef, dim)
    force = zeros(Float64, dim)
    energy = zero(Float64)
    for j = 1 : num_points
        if i != j
            @simd ivdep for k = 1 : dim
                @inbounds displ[k] = (
                    Float64(points[k,i]) - Float64(points[k,j]))
            end
            inv_dist = inv(approx_norm(displ))
            energy += inv_dist
            inv_dist *= abs2(inv_dist)
            @simd ivdep for k = 1 : dim
                @inbounds force[k] += inv_dist * displ[k]
            end
        end
    end
    force, energy
end

################################################################################

@everywhere const EVALUATOR = RKOCEvaluator{AccurateReal}(10, 16)
@everywhere const NUM_VARS = 136
@everywhere const NUM_CONSTRS = 1205

@everywhere const CONSTR_IDXS = [
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 18, 19, 20, 21, 22, 23,
    24, 25, 26, 27, 28, 31, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 51, 52,
    58, 59, 62, 75, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
    100, 101, 106, 107, 110, 115, 116, 123, 134, 135, 201, 202, 203, 205, 206,
    207, 208, 210, 211, 214, 215, 221, 222, 225, 226, 227, 230, 231, 238, 249,
    250, 269, 270, 316, 404, 406, 487, 488, 489, 490, 491, 492, 493, 496, 497,
    498, 499, 500, 501, 502, 507, 508, 511, 516, 517, 524, 535, 536, 539, 544,
    555, 556, 602, 603, 650, 690, 984]

@everywhere const CONSTR_DIM = length(CONSTR_IDXS)

@everywhere const TEMP_RES = Vector{AccurateReal}(undef, NUM_CONSTRS)
@everywhere const SHORT_RES = Vector{AccurateReal}(undef, CONSTR_DIM)
@everywhere const APX_SHORT_RES = Vector{Float64}(undef, CONSTR_DIM)
@everywhere const TEMP_JAC = Matrix{AccurateReal}(undef, NUM_CONSTRS, NUM_VARS)
@everywhere const SHORT_JAC = Matrix{AccurateReal}(undef, NUM_VARS, CONSTR_DIM)
@everywhere const SHORT_JCT = Matrix{AccurateReal}(undef, CONSTR_DIM, NUM_VARS)
@everywhere const APX_SHORT_JAC = Matrix{Float64}(undef, NUM_VARS, CONSTR_DIM)

################################################################################

@everywhere function approx_objective(x)
    evaluate_residual!(TEMP_RES, x, EVALUATOR)
    result = zero(Float64)
    for i = 1 : NUM_CONSTRS
        result += abs2(Float64(TEMP_RES[i]))
    end
    result
end

@everywhere function compute_residual(x)
    evaluate_residual!(TEMP_RES, x, EVALUATOR)
    @simd ivdep for j = 1 : length(CONSTR_IDXS)
        @inbounds SHORT_RES[j] = TEMP_RES[CONSTR_IDXS[j]]
    end
end

@everywhere function compute_jacobian(x)
    evaluate_jacobian!(TEMP_JAC, x, EVALUATOR)
    for j = 1 : length(CONSTR_IDXS)
        k = CONSTR_IDXS[j]
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds SHORT_JCT[j, i] = TEMP_JAC[k, i]
        end
    end
end

@everywhere function compute_orthonormalized_jacobian(x)
    evaluate_jacobian!(TEMP_JAC, x, EVALUATOR)
    for j = 1 : length(CONSTR_IDXS)
        k = CONSTR_IDXS[j]
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds SHORT_JAC[i, j] = TEMP_JAC[k, i]
        end
    end
    orthonormalize_columns!(SHORT_JAC)
    for j = 1 : length(CONSTR_IDXS)
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds APX_SHORT_JAC[i, j] = Float64(SHORT_JAC[i, j])
        end
    end
end

@everywhere function approx_coulomb_force_energy(
        points::SharedArray{T,2}, i::Int) where {T <: Real}
    force, energy = approx_coulomb_force_energy(sdata(points), i)
    compute_orthonormalized_jacobian(points[:,i])
    force - APX_SHORT_JAC * (transpose(APX_SHORT_JAC) * force), energy
end

################################################################################

@everywhere function constrain(x)
    x_old, obj_old = x, approx_objective(x)
    while true
        # log("        Computing Jacobian and residual...")
        compute_residual(x_old)
        compute_jacobian(x_old)
        # log("        Computing step direction...")
        direction = qr!(SHORT_JCT) \ SHORT_RES
        x_new = x_old - direction
        obj_new = approx_objective(x_new)
        if obj_new < obj_old
            # log("        Accepted step: ", @sprintf("%g", obj_new))
            x_old, obj_old = x_new, obj_new
        else
            # log("        Rejected step: ", @sprintf("%g", obj_new))
            step_size, obj_new = golden_section_search(
                h -> approx_objective(x_new - h * direction), 0, 1, 10)
            x_new = x_old - step_size * direction
            if obj_new < obj_old
                # log("        Accepted GSS step: ", @sprintf("%g", obj_new),
                #     " (step size ", @sprintf("%g", step_size), ")")
                x_old, obj_old = x_new, obj_new
            else
                # log("        Rejected GSS step: ", @sprintf("%g", obj_new))
                break
            end
        end
    end
    # log("        Final objective value: ", @sprintf("%g", obj_old))
    x_old, obj_old
end

@everywhere const EPS_THREE_HALVES = Float64(BigFloat(2)^-130)

@everywhere function perturb(x, direction, multiplier, i)
    x_new, obj = constrain(x + multiplier * direction)
    if obj < EPS_THREE_HALVES
        # log("Successfully moved point ", i, " (", obj,
        #     ") by ", multiplier * approx_norm(direction), ".")
        x_new
    else
        multiplier /= 16
        if multiplier < 0.0625
            multiplier = 0.0
            log("WARNING: Failed to move point ", i, " (", obj,
                ") by ", approx_norm(direction), ".")
        end
        perturb(x, direction, multiplier, i)
    end
end

@everywhere function perturb_wrapper(pts, direction, i)
    x = pts[:,i]
    # log("    Moving point ", i, " by distance ",
    #     @sprintf("%g", approx_norm(direction)), ".")
    x_new = perturb(x, direction, 1.0, i)
    # log("    Successfully moved point ", i, " by distance ",
    #     @sprintf("%g", approx_norm(x_new - x)), ".")
    x_new
end

################################################################################

const INPUT_FILENAME = maximum(filename
    for filename in readdir()
    if isfile(filename) && startswith(filename, "RKTK-POINTS-")
                        && endswith(filename, ".txt"))

log("Reading initial points from data file: ", INPUT_FILENAME)
const INPUT_POINTS = [AccurateReal.(BigFloat.(point))
    for point in split.(split(read(INPUT_FILENAME, String), "\n\n"))]
@assert all(length(p) == NUM_VARS for p in INPUT_POINTS)
const NUM_POINTS = length(INPUT_POINTS)

const POINTS = SharedArray{AccurateReal}((NUM_VARS, NUM_POINTS))
for j = 1 : NUM_POINTS
    pt = INPUT_POINTS[j]
    @simd ivdep for i = 1 : NUM_VARS
        @inbounds POINTS[i, j] = pt[i]
    end
end
log("Successfully read ", NUM_POINTS, " initial points.")

################################################################################

while true
    # log("Computing forces...")
    forces_energies = pmap(approx_coulomb_force_energy,
        [POINTS for _ = 1 : NUM_POINTS], 1 : NUM_POINTS)
    log("Total energy: ", sum(p[2] for p in forces_energies) / 2)
    forces = [p[1] for p in forces_energies]
    # log("Moving points...")
    forces *= sqrt(NUM_POINTS) / approx_norm(vcat(forces...)) / 5000
    new_points = pmap(perturb_wrapper,
        [POINTS for _ = 1 : NUM_POINTS], forces, 1 : NUM_POINTS)
    for i = 1 : NUM_POINTS
        pt = new_points[i]
        @simd ivdep for j = 1 : NUM_VARS
            @inbounds POINTS[j, i] = pt[j]
        end
    end
    filename = format(now(), dateformat"RKTK-POINT\S-yyyymmdd-HHMMSS-sss.txt")
    log("Writing points to file: ", filename)
    write(filename, join([join(string.(BigFloat.(POINTS[:,i])), "\n")
        for i = 1 : NUM_POINTS], "\n\n") * "\n")
end


text:
using StatsFuns, Test
using ChainRulesCore
using ChainRulesTestUtils
using Random

@testset "chainrules" begin
    x = exp(randn())
    y = exp(randn())
    z = logistic(randn())
    test_frule(betalogpdf, x, y, z)
    test_rrule(betalogpdf, x, y, z)

    x = exp(randn())
    y = exp(randn())
    z = exp(randn())
    test_frule(gammalogpdf, x, y, z)
    test_rrule(gammalogpdf, x, y, z)

    x = exp(randn())
    y = exp(randn())
    test_frule(chisqlogpdf, x, y)
    test_rrule(chisqlogpdf, x, y)

    x = exp(randn())
    y = exp(randn())
    z = exp(randn())
    test_frule(fdistlogpdf, x, y, z)
    test_rrule(fdistlogpdf, x, y, z)

    x = exp(randn())
    y = randn()
    test_frule(tdistlogpdf, x, y)
    test_rrule(tdistlogpdf, x, y)

    x = rand(1:100)
    y = logistic(randn())
    z = rand(1:x)
    test_frule(binomlogpdf, x, y, z)
    test_rrule(binomlogpdf, x, y, z)

    x = exp(randn())
    y = rand(1:100)
    test_frule(poislogpdf, x, y)
    test_rrule(poislogpdf, x, y)

    # test special case λ = 0
    _, pb = rrule(poislogpdf, 0.0, 0)
    _, x̄1, _ = pb(1)
    @test x̄1 == -1
    _, pb = rrule(poislogpdf, 0.0, 1)
    _, x̄1, _ = pb(1)
    @test x̄1 == Inf
end


text:
module SpikingNeuralNetworks

export SNN
const SNN = SpikingNeuralNetworks

using LinearAlgebra
using SparseArrays
using Requires
using UnPack

include("unit.jl")
include("main.jl")
include("util.jl")

include("neuron/if.jl")
include("neuron/if2.jl")
include("neuron/noisy_if.jl")
include("neuron/poisson.jl")
include("neuron/iz.jl")
include("neuron/adexp.jl")
include("neuron/hh.jl")
include("neuron/rate.jl")

include("synapse/spiking_synapse.jl")
include("synapse/rate_synapse.jl")
include("synapse/fl_synapse.jl")
include("synapse/fl_sparse_synapse.jl")
include("synapse/pinning_synapse.jl")
include("synapse/pinning_sparse_synapse.jl")

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
end

end

max_stars_repo_name: russelljjarvis/backupSNN

max_issues_repo_name: russelljjarvis/backupSNN

max_forks_repo_name: russelljjarvis/backupSNN


text:
module SparseND

import BEAST

struct Banded3D{T} <: AbstractArray{T,3}
    k0::Array{Int,2}
    data::Array{T,3}
    maxk0::Int
    maxk1::Int
    function Banded3D{T}(k0,data,maxk1) where T
        @assert size(k0) == size(data)[2:3]
        return new(k0,data,maximum(k0),maxk1)
    end
end

Banded3D(k0,data::Array{T},maxk1) where {T} = Banded3D{T}(k0,data,maxk1)

bandwidth(A::Banded3D) = size(A.data,1)

import Base: size, getindex, setindex!

# size(A::Banded3D) = size(A.data)
size(A::Banded3D) = tuple(size(A.k0)..., A.maxk1)

function getindex(A::Banded3D, m::Int, n::Int, k::Int)
    k0 = A.k0[m,n]
    k0 == 0 && return zero(eltype(A))
    l = k - k0 + 1
    l < 1 && return zero(eltype(A))
    l > bandwidth(A) && return zero(eltype(A))
    A.data[l,m,n]
end

function setindex!(A::Banded3D, v, m, n, k)
    k0 = A.k0[m,n]
    @assert k0 != 0
    @assert A.k0[m,n] <= k <= A.k0[m,n] + bandwidth(A) - 1
    A.data[k-A.k0[m,n]+1,m,n] = v
end

# """
#     convolve{T}(A::Banded3D, x::Array{T,2})
#
# Compute the *space-time* convolution of A and x.
#
# Returns an array `y` of size `(size(A,1),size(x,2))` such that
# `y[m,i] = sum([A[]])``
# """
# function convolve(A::Banded3D, x::Array{T,2}) where T
#     V = promote_type(eltype(A),eltype(x))
#     y = zeros(V,size(A,1),size(x,2))
#     bw = bandwidth(A)
#     for i in 1:size(A,1)
#         for l in 1 : size(y,2)
#             for j in 1:size(A,2)
#                 k0 = A.k0[i,j]
#                 for k in k0 : k0 + bw - 1
#                     l-k+1 < 1 && continue
#                     y[i,l] += A.data[k,i,j] * x[j,l-k+1]
#                 end
#             end
#         end
#     end
#     return y
# end

function Base.:+(A::Banded3D{T}, B::Banded3D{T}) where {T}

    M = size(A,1)
    N = size(A,2)

    @assert M == size(B,1)
    @assert N == size(B,2)

    Abw = size(A.data,1)
    Bbw = size(B.data,1)

    # keep track of empty columns
    Az = findall(A.k0 .== 0)
    Bz = findall(B.k0 .== 0)

    K0 = min.(A.k0, B.k0)
    K0[Az] = B.k0[Az]
    K0[Bz] = A.k0[Bz]

    AK1 = A.k0 .+ (Abw - 1); AK1[Az] .= -1
    BK1 = B.k0 .+ (Bbw - 1); BK1[Bz] .= -1
    K1 = max.(AK1, BK1)

    bw = maximum(K1 - K0 .+ 1)
    @assert bw > 0
    @assert bw >= Abw
    @assert bw >= Bbw
    data = zeros(T, bw, M, N)
    for m in axes(A.data,2)
        for n in axes(A.data,3)
            k0 = K0[m,n]
            @assert k0 != 0

            Ak0 = A.k0[m,n]
            if Ak0 != 0
                for Al in 1:Abw
                    k = Ak0 + Al - 1
                    l = k - k0 + 1
                    data[l,m,n] += A.data[Al,m,n]
                end
            end

            Bk0 = B.k0[m,n]
            if Bk0 != 0
                for Bl in 1:Bbw
                    k = Bk0 + Bl - 1
                    l = k - k0 + 1
                    l < 1 && @show m, n, k0, k, Bk0, Bl
                    d = data[l,m,n]
                    data[l,m,n] = d + B.data[Bl,m,n]
                end
            end
        end
    end

    Banded3D(K0, data, max(A.maxk1, B.maxk1))
end

function BEAST.convolve(Z::SparseND.Banded3D,x,j,k_start)
    T = promote_type(eltype(Z), eltype(x))
    M,N,L = size(Z)
    K = size(Z.data,1)
    @assert M == size(x,1)
    y = zeros(T,M)
    for n in 1:N
        for m in 1:M
            k0 = Z.k0[m,n] # k0 is 1-based
            l0 = max(1, k_start - k0 + 1)
            l1 = min(K, j - k0 + 1)
            for l in l0 : l1
                k = k0 + l - 1
                # j - k + 1 < 1 && break
                y[m] += Z.data[l,m,n] * x[n,j - k + 1]
            end
        end
    end
    return y
end

struct MatrixOfConvolutions{T} <: AbstractArray{Vector{T},2}
    banded::Banded3D{T}
end

function Base.eltype(x::MatrixOfConvolutions{T}) where {T}
    Vector{T}
end
Base.size(x::MatrixOfConvolutions) = size(x.banded)[1:2]
function Base.getindex(x::MatrixOfConvolutions, m, n)
    return x.banded[m,n,:]
end

struct SpaceTimeData{T} <: AbstractArray{Vector{T},1}
    data::Array{T,2}
end

Base.eltype(x::SpaceTimeData{T}) where {T} = Vector{T}
Base.size(x::SpaceTimeData) = size(x.data)[1]
Base.getindex(x::SpaceTimeData, i::Int) = x.data[i,:]

end # module


text:
"""
    get_SOA_table(id)
    get_SOA_table(table_name)

Given the id or name of a `mort.SOA.org` table, grab it and return it as a `MortalityTable`.

!!! Remember that not all tables have been tested to work.
"""
function get_SOA_table(id::Int)
    readXTbML(joinpath(artifact"mort.soa.org","t$id.xml"))
end

function get_SOA_table(table_name::String;source_map=table_source_map)
    entry = findfirst(x-> x.name == table_name, source_map)
    if entry === nothing
        search_method = StringDistances.Partial(StringDistances.Levenshtein())
        suggestion,_ = StringDistances.findnearest(table_name,[x.name for x in table_source_map],search_method)
        throw(ArgumentError("table name \"table_name\" not found in table set; " *
            "most similar available name is: \"$suggestion\""))
    end
    readXTbML(joinpath(artifact"mort.soa.org","t$(source_map[entry].id).xml"))
end

"""
    table(id)
    table(name)

Given the id or name of a `mort.SOA.org` table, grab it and return it as a `MortalityTable`.

!!! Remember that not all tables have been tested to work.
"""
table = get_SOA_table


text:
export OffHilbert,OffSingularIntegral,Stieltjes,Cauchy

#############
# OffHilbert implements:
#
#       OH f(z) := 1/π\int_Γ f(t)/(t-z) dt,  z ∉ Γ
#
# OffSingularIntegral implements:
#
#       OSI f(z) := 1/π\int_Γ f(t)/(t-z) ds(t),  z ∉ Γ
#
#############

for Op in (:OffHilbert,:OffSingularIntegral)
    @eval begin
        struct $Op{D<:Space,R<:Space,T} <: Operator{T}
            data::BandedMatrix{T,Matrix{T},Base.OneTo{Int}}
            domainspace::D
            rangespace::R
            order::Int
        end

        getindex(C::$Op,k::Integer,j::Integer) =
            k ≤ size(C.data,1) && j ≤ size(C.data,2) ? C.data[k,j] : zero(eltype(C))

        convert(::Type{Operator{T}},OH::$Op) where {T} =
            $Op{typeof(OH.domainspace),
                typeof(OH.rangespace),
                T}(OH.data,OH.domainspace,OH.rangespace,OH.order)

        $Op(ds::Space,rs::Space) = $Op(ds,rs,1)
        $Op(data::BandedMatrix,ds::Space,rs::Space) = $Op(data,ds,rs,1)

        $Op(ds::PeriodicDomain,rs::PeriodicDomain,order) = $Op(Laurent(ds),Laurent(rs),order)
        $Op(ds::PeriodicDomain,rs::PeriodicDomain) = $Op(Laurent(ds),Laurent(rs))

        $Op(ds::Space,rs::Number,order) = $Op(ds,Space(rs),order)
        $Op(ds::Space,rs::Number) = $Op(ds,ConstantSpace(Point(rs)))

        domainspace(C::$Op) = C.domainspace
        rangespace(C::$Op) = C.rangespace
        bandwidths(C::$Op) = bandwidths(C.data)
        # divide size by blocklengths
        blockbandwidths(C::$Op) = ((size(C.data,1)-1)÷maximum(blocklengths(rangespace(C))),
                                 (size(C.data,2)-1)÷maximum(blocklengths(domainspace(C))))
    end
end

function default_OffHilbert(ds::Space,rs::Space,order::Int)
    @assert order==1
    tol=1E-13

    vv=Vector{Vector{ComplexF64}}()
    m=min(100,dimension(rs))
    for k=1:2:1000
        b=Fun(ds,[zeros(k-1);1.])
        v1=Fun(x->-stieltjes(b,x)/π,rs,m)
        b=Fun(ds,[zeros(k);1.])
        v2=Fun(x->-stieltjes(b,x)/π,rs,m)
        if m ≥ 2 && (abs(v1.coefficients[end-1])>100tol || abs(v1.coefficients[end])>100tol ||
                        abs(v2.coefficients[end-1])>100tol || abs(v2.coefficients[end])>100tol)
            warn("OffHilbert not resolved with $m rows")
        end
        if norm(v1.coefficients,Inf)<tol &&
            norm(v2.coefficients,Inf)<tol
            C=zeros(ComplexF64,mapreduce(length,max,vv),length(vv))
            for j=1:length(vv)
                @inbounds C[1:length(vv[j]),j] = vv[j]
            end
            return OffHilbert(convert(BandedMatrix,C),ds,rs,order)
        end

        push!(vv,v1.coefficients)
        push!(vv,v2.coefficients)
    end

    warn("Max Iteration Reached for OffHilbert from "*string(ds)*" to "*string(rs))
    OffHilbert(convert(BandedMatrix,C),ds,rs,order)
end

function default_OffSingularIntegral(::Type{T}, ds::Space, rs::Space, order::Int) where T
    tol=1E-13

    vv=Vector{Vector{T}}()
    m=100
    for k=1:2:1000
        b=Fun(ds,[zeros(k-1);1.])
        v1=Fun(x->singularintegral(order,b,x),rs,m)
        b=Fun(ds,[zeros(k);1.])
        v2=Fun(x->singularintegral(order,b,x),rs,m)
        if isinf(dimension(rs)) && (abs(v1.coefficients[end-1])>100tol || abs(v1.coefficients[end])>100tol ||
            abs(v2.coefficients[end-1])>100tol || abs(v2.coefficients[end])>100tol)
            warn("OffSingularIntegral not resolved with $m rows")
        end
        if norm(v1.coefficients,Inf)<tol &&
            norm(v2.coefficients,Inf)<tol
            C=zeros(T,mapreduce(length,max,vv),length(vv))
            for j=1:length(vv)
                @inbounds C[1:length(vv[j]),j]=vv[j]
            end
            return OffSingularIntegral(convert(BandedMatrix,C),ds,rs,order)
        end

        push!(vv,v1.coefficients)
        push!(vv,v2.coefficients)
    end

    warn("Max Iteration Reached for OffHilbert from "*string(ds)*" to "*string(rs))
    C=zeros(T,mapreduce(length,max,vv),length(vv))
    for j=1:length(vv)
        @inbounds C[1:length(vv[j]),j]=vv[j]
    end
    OffSingularIntegral(convert(BandedMatrix,C),ds,rs,order)
end

# TODO: SingularInteral( , 0 ) is really different
default_OffSingularIntegral(ds::Space, rs::Space, order::Int) =
    default_OffSingularIntegral(order == 0 ? Float64 : ComplexF64, ds, rs, order)

OffHilbert(ds::Space,rs::Space,order::Int) = default_OffHilbert(ds,rs,order)
OffSingularIntegral(ds::Space,rs::Space,order::Int) = default_OffSingularIntegral(ds,rs,order)

## JacobiWeight

for (Op,Len) in ((:OffHilbert,:complexlength),(:OffSingularIntegral,:arclength))
    @eval begin
        function $Op(ds::JacobiWeight{Ultraspherical{Int,DD,RR},DD},rs::Space,ord::Int) where {DD<:IntervalOrSegment,RR}
            @assert order(ds.space) == 1
            @assert ds.α==ds.β==0.5
            d = domain(ds)
            C = (.5*$Len(d))^(1-ord) # probably this is right for all ords ≥ 2. Certainly so for 0,1.

            if ord == 0
                z=Fun(identity,rs)
                x=mobius(ds,z)
                y=joukowskyinverse(Val{true},x)
                yk,ykp1=y,y*y
                ret=Array{typeof(y)}(undef,300)
                ret[1]=-.5logabs(2y)+.25real(ykp1)
                n,l,u = 1,ncoefficients(ret[1])-1,0
                while norm(ret[n].coefficients)>100eps()
                    n+=1
                    if n > length(ret) resize!(ret,2length(ret)) end  # double preallocated ret
                    ykp1*=y
                    ret[n]=chop!(.5*(real(ykp1)/(n+1)-real(yk)/(n-1)) ,100eps())  #will be length 2n-1
                    yk*=y
                    u+=1   # upper bandwidth
                    l=max(l,ncoefficients(ret[n])-n)
                end
            elseif ord == 1
                y=Fun(z->joukowskyinverse(Val{true},mobius(ds,z)),rs)
                ret=Array{typeof(y)}(undef,300)
                ret[1]=-y
                n,l,u = 1,ncoefficients(ret[1])-1,0
                while norm(ret[n].coefficients)>100eps()
                    n+=1
                    if n > length(ret) resize!(ret,2length(ret)) end  # double preallocated ret
                    ret[n]=chop!(y*ret[n-1],100eps())  #will be length 2n-1
                    u+=1   # upper bandwidth
                    l=max(l,ncoefficients(ret[n])-n)
                end
            else
                error("Not implemented for order=$ord")
            end

            M=BandedMatrix{promote_type(typeof(C),cfstype(y))}(Zeros(l+1,n), (l,u))
            for k=1:n,j=1:min(l+1,ncoefficients(ret[k]))
                M[j,k]=C*ret[k].coefficients[j]
            end
            $Op(M,ds,rs,ord)
        end

        function $Op(ds::JacobiWeight{ChebyshevDirichlet{1,1,DD,RR},DD},rs::PolynomialSpace,ord::Int) where {DD<:IntervalOrSegment,RR}
            @assert ds.α==ds.β==-0.5
            d = domain(ds)
            C = (.5*$Len(d))^(1-ord) # probably this is right for all orders ≥ 2. Certainly so for 0,1.

            if ord == 0
                z=Fun(identity,rs)
                x=mobius(ds,z)
                y=joukowskyinverse(Val{true},x)
                yk,ykp1=y,y*y
                ret=Array{typeof(y)}(undef,300)
                ret[1]=-logabs(2y/C)
                ret[2]=-real(yk)
                ret[3]=chop!(-ret[1]-.5real(ykp1),100eps())
                n,l,u = 3,max(ncoefficients(ret[1])-1,ncoefficients(ret[2])-2,ncoefficients(ret[3])-3),2
                while norm(ret[n].coefficients)>100eps()
                    n+=1
                    if n > length(ret) resize!(ret,2length(ret)) end  # double preallocated ret
                    ykp1*=y
                    ret[n]=chop!(real(yk)/(n-3)-real(ykp1)/(n-1),100eps())  #will be length 2n-1
                    yk*=y
                    u+=1   # upper bandwidth
                    l=max(l,ncoefficients(ret[n])-n)
                end
            elseif ord == 1
                z=Fun(identity,rs)
                x=mobius(ds,z)
                y=joukowskyinverse(Val{true},x)
                ret=Array{typeof(y)}(undef,300)
                ret[1]=-1/sqrtx2(x)
                ret[2]=x*ret[1]+1
                ret[3]=2y
                n,l,u = 3,max(ncoefficients(ret[1])-1,ncoefficients(ret[2])-2,ncoefficients(ret[3])-3),2
                while norm(ret[n].coefficients)>100eps()
                    n+=1
                    if n > length(ret) resize!(ret,2length(ret)) end  # double preallocated ret
                    ret[n]=chop!(y*ret[n-1],100eps())  #will be length 2n-1
                    u+=1   # upper bandwidth
                    l=max(l,ncoefficients(ret[n])-n)
                end
            end

            M=BandedMatrix{promote_type(typeof(C),cfstype(y))}(Zeros(l+3,n),(l,u))
            for k=1:n,j=1:min(l+3,ncoefficients(ret[k]))
                M[j,k]=C*ret[k].coefficients[j]
            end
            $Op(M,ds,rs,ord)
        end

    end
end

function OffHilbert(DS::Laurent{D1,R1},RS::Laurent{D2,R2},ord::Int) where {D1<:Circle,D2<:Circle,R1,R2}
    ds=domain(DS);rs=domain(RS)
    @assert ord==1

    # Correct for orrientation
    if !ds.orientation
        rDS=reverseorientation(DS)
        return (-OffHilbert(rDS,RS))*Conversion(DS,rDS)
    end

    if !rs.orientation
        rRS=reverseorientation(RS)
        return Conversion(rRS,RS)*(-OffHilbert(DS,rRS))
    end

    c2=rs.center;c1=ds.center
    r2=rs.radius;r1=ds.radius

    # make positive orientation
    if r1>r2&&abs(c1-c2)<r1  # we are inside the circle, use Taylor series
        M=interior_cauchy(Circle(c1,r1),Circle(c2,r2))
    elseif r1<r2&&abs(c1-c2)<r2 # we surround the domain, use Hardy{False} series
        M=exterior_cauchy(Circle(c1,r1),Circle(c2,r2))
    else
        M=disjoint_cauchy(Circle(c1,r1),Circle(c2,r2))
    end

    OffHilbert(2im*M,DS,RS)
end

function OffHilbert(DS::Fourier{D1,R1},RS::Fourier{D2,R2},ord::Int) where {D1<:Circle,D2<:Circle,R1,R2}
    LD=Laurent(domain(DS))
    LR=Laurent(domain(RS))
    Conversion(LR,RS)*OffHilbert(LD,LR,ord)*Conversion(DS,LD)
end

## Special cases

function exterior_cauchy(b::Circle,a::Circle)
    c=b.center
    r=b.radius

    S=Fun(a,[0.0,0,1])  # Shift to use bandedness
    ret=Array{Fun{Laurent{typeof(a),ComplexF64},ComplexF64}}(undef,300)
    ret[1]=Fun(z->(r/(z-c)),a)
    n=1
    m=ncoefficients(ret[1])-2
    f1=ret[1]*S
    while norm(ret[n].coefficients)>100eps()
        n+=1
        if n > length(ret)
            resize!(ret,2length(ret))
        end
        ret[n]=chop!(f1*ret[n-1],100eps())
        m=max(m,ncoefficients(ret[n])-2)
    end

    M=BandedMatrix{Complex{Float64}}(Zeros(2n,2n),(m,0))
    #j+2k-2≤2n
    #j≤2(n-k)+2
    for k=1:n,j=2:2:min(length(ret[k].coefficients),2(n-k)+2)
        M[j+2k-2,2k]=-ret[k].coefficients[j]
    end
    M
end

function interior_cauchy(a::Circle,b::Circle)
    z=mappoint(a,Circle(),Fun(b))

    ret=Array{Fun{Laurent{typeof(b),ComplexF64},ComplexF64}}(undef,300)
    ret[1]=ones(b)
    n=1
    m=0

    while norm(ret[n].coefficients)>100eps()
        n+=1
        if n > length(ret)
            # double preallocated ret
            resize!(ret,2length(ret))
        end
        ret[n]=z*ret[n-1]  #will be length 2n-1

        # find bandwidth by checking how many coefficients are zero
        # we jump over negative coefficients
        for j=1:2:2n-1
            if norm(ret[n].coefficients[j])>100eps()
                m=max(m,2n-j-1)
                break
            end
        end
    end

    M=BandedMatrix{Complex{Float64}}(Zeros(2n-1,2n-1),(0,m))
    for k=1:n,j=max(1,2k-1-m):2:2k-1
        M[j,2k-1]=ret[k].coefficients[j]
    end

    M
end

function disjoint_cauchy(a::Circle, b::Circle)
    c=a.center
    r=a.radius

    f=Fun(z->r/(z-c),b)

    ret=Array{Fun{Laurent{typeof(b),ComplexF64},ComplexF64}}(undef,300)
    ret[1]=f
    n=1

    l=ncoefficients(f)-2   #lower bandwidth
    u=1             #upper bandwidth

    while norm(ret[n].coefficients)>100eps()
        n+=1
        if n > length(ret)
            # double preallocated ret
            resize!(ret,2length(ret))
        end
    ret[n]=chop!(f*ret[n-1],100eps())  #will be length 2n-1
    u=max(u,ncoefficients(ret[n])-2n)   # upper bandwidth

        # find bandwidth by checking how many coefficients are zero
        # we jump over negative coefficients
        for j=1:2:ncoefficients(ret[n])
            if norm(ret[n].coefficients[j])>100eps()
                l=max(l,2n-j)
                break
            end
        end
    end

    M=BandedMatrix{Complex{Float64}}(Zeros(2n-1,2n),(l,u))
    for k=1:n,j=max(1,2k-u):2:min(ncoefficients(ret[k]),2n-1)
            M[j,2k]=-ret[k].coefficients[j]
    end
    M
end

#############
# Cauchy implements the Cauchy operator corresponding to evaluating the Cauchy transform
#
#       C f(z) := 1/(2πi)\int_Γ f(t)/(t-z) dt
#
# It is given in terms of the Stieltjes operator
#
#       S f(z) := \int_Γ f(t)/(z-t) dt = -2πi*C f(z)
#
# note that the domain of domainspace must be different than the domain of rangespace
#
# The notion of C^± for the left/right limits of the Cauchy operator
# with the domains matching is represented
# using the Hilbert operator and the formulae
#
#    C^+  -  C^- = I
#    C^+  +  C^- = -im*H
#
#   Or for the Stieltjes operator
#
#    S^+ - S^- = -2πi*I
#    S^+ + S^- = -2π*H
#
############
Stieltjes(d,r,ord) = (ord==0 ? π : -π)*OffHilbert(d,r,ord)
Stieltjes(d,r) = (-π)*OffHilbert(d,r)

Cauchy(s::Bool,d)=(s ? 0.5 : -0.5)*I + (-0.5im)*Hilbert(d)
Cauchy(s::Int,d)=Cauchy(s==1,d)
Cauchy(s::Union{Int,Bool})=Cauchy(s,UnsetSpace())
Cauchy(ds,rs,ord)=(1/(2*im))*OffHilbert(ds,rs,ord)
Cauchy(ds,rs)=Cauchy(ds,rs,1)

## Stiejles Functional

function hornervector(y0)
    r=Array{typeof(y0)}(undef,200)
    r[1]=y0
    k=1
    tol=eps()
    while(abs(r[k])>tol)
        k+=1
        if k>length(r)
            resize!(r,2length(r))
        end
        r[k]=r[k-1]*y0
    end

    r[1:k]
end

## OffHilbert Functional

function HornerFunctional(y0,sp,cs)
    v = hornervector(y0)
    FiniteOperator(transpose(v),sp,cs)
end

function OffHilbert(sp::JacobiWeight{Ultraspherical{Int,DD,RR},DD},cs::ConstantSpace{<:Point},k::Int) where {DD<:IntervalOrSegment,RR}
    z = convert(Number, domain(cs))
    @assert k == 1
    @assert order(sp.space) == 1
    if sp.α == sp.β == 0.5
        # this translates the following cauchy to a functional
        #    0.5im*hornersum(cfs,joukowskyinverse(Val{true},mobius(u,z)))
        # which consists of multiplying by 2*im
        -HornerFunctional(joukowskyinverse(Val{true},mobius(sp,z)),sp,cs)
    else
        # calculate directly
        r=Vector{eltype(z)}()
        for k=1:10000
            push!(r,-stieltjes(Fun(sp,[zeros(k-1);1.]),z)/π)
            if abs(last(r)) < eps()
                break
            end
        end
        FiniteOperator(transpose(r),sp,cs)
    end
end

for Op in (:OffHilbert, :OffSingularIntegral)
    defOp = Meta.parse("default_"*string(Op))
    @eval begin
        function $Op(sp::JacobiWeight{Chebyshev{DD,RR},DD},z::Space,k::Int) where {DD,RR}
            if sp.β == sp.α == 0.5
                #try converting to Ultraspherical(1)
                us=JacobiWeight(sp.β,sp.α,Ultraspherical(1,domain(sp)))
                $Op(us,z,k)*Conversion(sp,us)
            else
                $defOp(sp,z,k)
            end
        end
    end
end

function OffHilbert(sp::JacobiWeight{ChebyshevDirichlet{1,1,DD,RR},DD},cs::ConstantSpace{<:Point},k::Int) where {DD<:IntervalOrSegment,RR}
    z = Number(domain(cs))
    @assert k == 1
    if sp.α == sp.β == -0.5
        z=mobius(sp,z)

        sx2z=sqrtx2(z)
        sx2zi=1/sx2z

        FiniteOperator(transpose([-sx2zi;1-sx2zi;2*hornervector(z-sx2z)]),sp,cs)
    else
        # try converting to Canonical
        us=JacobiWeight(sp.β,sp.α,Chebyshev(domain(sp)))
        OffHilbert(us,cs)*Conversion(sp,us)
    end
end


text:
d, batch_size = 5, 100
limits = repeat([-1.0 1.0], d, 1)
mf = (x -> 100, x -> -100)
csampler = UniformContexts(limits)
policy = RandomPolicy(length(mf))
rsampler = GaussianRewards(mf)
metrics = (FunctionalRegret(mf),)
driver = StandardDriver(csampler, policy, rsampler, metrics)

X, a, r = driver(batch_size)
@test size(X) == (d, batch_size)
@test all(r[1, a .== 1] .> 90)
@test all(r[1, a .== 2] .< -90)
@test all(metrics[1].regret[a .== 1] .== 0.0)
@test all(metrics[1].regret[a .== 2] .== 200.0)

driver = StandardDriver(csampler, policy, rsampler)
X, a, r = driver(batch_size)
@test size(X) == (d, batch_size)
@test all(r[1, a .== 1] .> 90)
@test all(r[1, a .== 2] .< -90)

@test_throws(
    ArgumentError("num_batches and batch_size must be positive"), run!(0, 0, driver)
)

# Test latent driver
tform = z -> 100 + z[1]
driver = LatentDriver(csampler, policy, rsampler, tform)
X, a, r = driver(1)
@test 99 <= X[1, 1] <= 101


text:
const Maybe{T} = Union{Nothing, T}

"""
    NoDefault

Type describes a field should have no default value.
"""
struct NoDefault end

"""
    const no_default = NoDefault()

Constant instance for [`NoDefault`](@ref) that
describes a field should have no default value.
"""
const no_default = NoDefault()

"""
    abstract type JLExpr end

Abstract type for Julia syntax type.
"""
abstract type JLExpr end

"""
    JLFunction <: JLExpr

Type describes a Julia function declaration expression.

# Example

Construct a function expression

```julia
julia> JLFunction(;name=:foo, args=[:(x::T)], body= quote 1+1 end, head=:function, whereparams=[:T])
function foo(x::T) where {T}
    #= REPL[25]:1 =#
    1 + 1
end
```

Decompose a function expression

```julia
julia> ex = :(function foo(x::T) where {T}
           #= REPL[25]:1 =#
           1 + 1
       end)
:(function foo(x::T) where T
      #= REPL[26]:1 =#
      #= REPL[26]:3 =#
      1 + 1
  end)

julia> jl = JLFunction(ex)
function foo(x::T) where {T}
    #= REPL[26]:1 =#
    #= REPL[26]:3 =#
    1 + 1
end
```

Generate `Expr` from `JLFunction`

```julia
julia> codegen_ast(jl)
:(function foo(x::T) where T
      #= REPL[26]:1 =#
      #= REPL[26]:3 =#
      1 + 1
  end)
```
"""
mutable struct JLFunction <: JLExpr
    head::Symbol  # function def must have a head
    name::Any  # name can be nothing, Symbol, Expr
    args::Vector{Any}
    kwargs::Maybe{Vector{Any}}
    rettype::Any
    whereparams::Maybe{Vector{Any}}
    body::Any
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
end

function JLFunction(;
        head=:function, name=nothing,
        args=[], kwargs=nothing,
        rettype=nothing,
        whereparams=nothing, body=Expr(:block),
        line=nothing, doc=nothing
    )
    JLFunction(head, name, args, kwargs, rettype, whereparams, body, line, doc)
end

"""
    JLField <: JLExpr

Type describes a Julia field in a Julia struct.
"""
mutable struct JLField <: JLExpr
    name::Symbol
    type::Any
    doc::Maybe{String}
    line::Maybe{LineNumberNode}
end

"""
    JLField(;name, type=Any, doc=nothing, line=nothing)

Create a `JLField` instance.
"""
function JLField(;name, type=Any, doc=nothing, line=nothing)
    JLField(name, type, doc, line)
end

"""
    JLKwField <: JLExpr

Type describes a Julia field that can have a default value in a Julia struct.
"""
mutable struct JLKwField <: JLExpr
    name::Symbol
    type::Any
    doc::Maybe{String}
    line::Maybe{LineNumberNode}
    default::Any
end

"""
    JLKwField(;kw...)

Create a `JLKwField` instance.
"""
function JLKwField(;name, type=Any, doc=nothing, line=nothing, default=no_default)
    JLKwField(name, type, doc, line, default)
end

"""
    JLStruct <: JLExpr

Type describes a Julia struct.

# Example

Construct a Julia struct.

```julia
julia> JLStruct(;name=:Foo, typevars=[:T], fields=[JLField(;name=:x, type=Int)])
struct Foo{T}
    x::Int64
end
```

Decompose a Julia struct expression

```julia
julia> ex = :(struct Foo{T}
           x::Int64
       end)
:(struct Foo{T}
      #= REPL[31]:2 =#
      x::Int64
  end)

julia> jl = JLStruct(ex)
struct Foo{T}
    #= REPL[31]:2 =#
    x::Int64
end
```

Generate a Julia struct expression

```julia
julia> codegen_ast(jl)
:(struct Foo{T}
      #= REPL[31]:2 =#
      x::Int64
  end)
```
"""
mutable struct JLStruct <: JLExpr
    name::Symbol
    ismutable::Bool
    typevars::Vector{Any}
    supertype::Any
    fields::Vector{JLField}
    constructors::Vector{JLFunction}
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
    misc::Any
end

"""
    JLStruct(;kw...)

Create a `JLStruct` instance.
"""
function JLStruct(;
    name, ismutable=false,
    typevars=[], supertype=nothing,
    fields=JLField[], constructors=JLFunction[],
    line=nothing, doc=nothing, misc=nothing)
    JLStruct(name, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLKwStruct <: JLExpr

Type describes a Julia struct that allows keyword definition of defaults.
"""
mutable struct JLKwStruct <: JLExpr
    name::Symbol
    typealias::Maybe{String}
    ismutable::Bool
    typevars::Vector{Any}
    supertype::Any
    fields::Vector{JLKwField}
    constructors::Vector{JLFunction}
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
    misc::Any
end

"""
    JLKwStruct(;kw...)

Create a `JLKwStruct` instance.
"""
function JLKwStruct(;name, typealias=nothing,
    ismutable=false, typevars=[], supertype=nothing,
    fields=JLField[], constructors=JLFunction[],
    line=nothing, doc=nothing, misc=nothing)
    JLKwStruct(name, typealias, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLIfElse <: JLExpr

`JLIfElse` describes a Julia `if ... elseif ... else ... end` expression. It allows one to easily construct
such expression by inserting condition and code block via a map.

# Example

### Construct JLIfElse object

One can construct an `ifelse` as following

```julia
julia> jl = JLIfElse()
nothing

julia> jl[:(foo(x))] = :(x = 1 + 1)
:(x = 1 + 1)

julia> jl[:(goo(x))] = :(y = 1 + 2)
:(y = 1 + 2)

julia> jl.otherwise = :(error("abc"))
:(error("abc"))

julia> jl
if foo(x)
    x = 1 + 1
elseif goo(x)
    y = 1 + 2
else
    error("abc")
end
```

### Generate the Julia `Expr` object

to generate the corresponding `Expr` object, one can call [`codegen_ast`](@ref).

```julia
julia> codegen_ast(jl)
:(if foo(x)
      x = 1 + 1
  elseif goo(x)
      y = 1 + 2
  else
      error("abc")
  end)
```
"""
mutable struct JLIfElse <: JLExpr
    conds::Vector{Any}
    stmts::Vector{Any}
    otherwise::Any
end

"""
    JLIfElse()

Create an emptry `ifelse` syntax type instance.
"""
JLIfElse() = JLIfElse([], [], nothing)

function Base.getindex(jl::JLIfElse, cond)
    idx = findfirst(jl.conds) do x
        cond == x
    end
    idx === nothing && error("cannot find condition: $cond")
    return jl.stmts[idx]
end

function Base.setindex!(jl::JLIfElse, stmt, cond)
    idx = findfirst(jl.conds) do x
        x == cond
    end
    if idx === nothing
        push!(jl.conds, cond)
        push!(jl.stmts, stmt)
    else
        jl.stmts[idx] = stmt
    end
    return stmt
end

Base.length(jl::JLIfElse) = length(jl.conds)
function Base.iterate(jl::JLIfElse, st=1)
    st > length(jl) && return
    jl.conds[st] => jl.stmts[st], st + 1
end

"""
    JLFor <: JLExpr

Syntax type for Julia for loop.
"""
struct JLFor <: JLExpr
    vars::Vector{Any}
    iterators::Vector{Any}
    kernel::Any
end

"""
    JLFor(;vars=[], iterators=[], kernel=nothing)

Generate a `JLFor` object.

# Kwargs

- `vars`: loop variables.
- `iterators`: loop iterators.
- `kernel`: loop kernel.
"""
JLFor(;vars=[], iterators=[], kernel=nothing) = JLFor(vars, iterators, kernel)


text:
using Glob, ParquetFiles, DataFrames, Statistics, Queryverse

function nunique(x)
    size(unique(x), 1)
end

function run()
    ;files = glob("*.snappy.parquet", "/data/julia/")
    files = glob("*.feather", "/data/julia/")
    dfs = DataFrame.(load.(files))
    df = reduce(vcat, dfs)
    df.sales = df.price .* df.quantity
    gdf = groupby(df, :member_id)
    summary = combine(gdf, :sales => sum => :total_spend,
                           :sales => mean => :avg_basket_size,
                           :price => mean => :avg_price,
                           nrow => :n_transactions,
                           :date => nunique => :n_visits,
                           :brand_id => nunique => :n_brands,
                           :style_id => nunique => :n_styles)
    summary |> save("julia.feather")
end

@time run()

; Parquet
; 1 Part: 86.532695 seconds (502.42 M allocations: 24.133 GiB, 6.02% gc time)
; 12 Parts: run: 868.278368 seconds (5.25 G allocations: 252.436 GiB, 19.38% gc time)

; Feather
; 1 Part: 16.211124 seconds (55.16 M allocations: 3.250 GiB, 6.64% gc time)
; 12 Parts: 41.322144 seconds (109.18 M allocations: 11.552 GiB, 25.66% gc time)

; function save_feather(file)
;     df = DataFrame(load(file))
;     new_file = replace(file, "snappy.parquet" => "feather")
;     df |> save(new_file)
; end
;
; files = glob("*.snappy.parquet", "/data/julia/")
; for file in files
;     save_feather(file)
; end

df = DataFrame(load("julia.feather"))
size(df, 1)
size(df, 2)


text:
using FeedbackParticleFilters, Random, Statistics, LinearAlgebra, Distributions, StatsBase

println("Testing BPF.jl:")

@testset "BPF.jl" begin

    print("  inner constructor for BPFState")
    ens   = WeightedParticleEnsemble(randn(3,10))
    state = BPFState(ens)
    print(".")
    @test state.ensemble == ens
    println("DONE.")

    print("  inner constructor for BPF")
    f(x)   = -x
    g(x)   = hcat([1.; -2.; 1.])
    h(x)   = [x[1] * x[2], x[2] * x[3]]
    init   = [0., 0., 0.]
    st_mod = DiffusionStateModel(f, g, init)
    ob_mod = DiffusionObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    N      = 10
    alpha  = 0.5
    filter = BPF(st_mod, ob_mod, N, alpha)
    print(".")
    @test filter.state_model == st_mod
    print(".")
    @test filter.obs_model == ob_mod
    print(".")
    @test filter.N == N
    print(".")
    @test filter.alpha == alpha

    print("  method initial_condition")
    print(".")
    @test initial_condition(filter) ≈ [0., 0., 0.]
    println("DONE.")

    print("  method no_of_particles")
    print(".")
    @test no_of_particles(filter) == N
    print(".")
    @test no_of_particles(state) == N
    println("DONE.")

    print("  method state_model")
    print(".")
    @test state_model(filter) == st_mod
    println("DONE.")

    print("  method obs_model")
    print(".")
    @test obs_model(filter) == ob_mod
    println("DONE.")

    print("  method for Statistics.mean")
    print(".")
    @test Statistics.mean(state) ≈ mean(ens)
    println("DONE.")

    print("  method for Statistics.cov")
    print(".")
    @test Statistics.cov(state) ≈ cov(ens)
    println("DONE.")

    print("  method for Statistics.var")
    print(".")
    @test Statistics.var(state) ≈ var(ens)
    println("DONE.")

    print("  outer constructors for BPFState")
    state2 = BPFState(filter)
    print(".")
    @test no_of_particles(state2) == N
    print(".")
    @test state2.ensemble.positions == zeros(3,N)
    filt_prob = FilteringProblem(st_mod, ob_mod)
    state3 = BPFState(filt_prob, 2*N)
    print(".")
    @test no_of_particles(state3) == 2*N
    println("DONE.")

    print("  outer constructors for BPF")
    filter2 = BPF(filt_prob, N, alpha)
    print(".")
    @test no_of_particles(filter2) == N
    print(".")
    @test filter2.alpha == alpha
    println("DONE.")

    print("  method initialize")
    state4 = initialize(filter)
    print(".")
    @test state4 isa BPFState
    print(".")
    @test state4.ensemble.positions == zeros(3,10)
    print(".")
    @test state4.ensemble.weights == StatsBase.ProbabilityWeights(fill(1/N, N))
    println("DONE.")

    print("  method propagate!")
    f(x)     = -x
    g(x)     = hcat([1.; -2.; 1.])
    h(x)     = [x[1] * x[2], x[2] * x[3]]
    init     = [0., 0., 0.]
    st_mod   = DiffusionStateModel(f, g, init)
    ob_mod   = DiffusionObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    N        = 2
    alpha    = 0.5
    filter   = BPF(st_mod, ob_mod, N, alpha)
    state    = initialize(filter)
    Random.seed!(0)
    propagate!(state, filter, 0.01)
    print(".")
    @test state.ensemble.positions ≈ [0.06791074260357777 -0.013485387193052173; -0.1656826965800072 -0.11732341492662196; -0.03530074003005963 0.029733585084941616]
    println("DONE.")

    print("  method assimilate!")
    ens      = WeightedParticleEnsemble(randn(3,10))
    state    = BPFState(ens)
    oldstate = deepcopy(state)
    f(x)     = -x
    g(x)     = hcat([1.; -2.; 1.])
    h(x)     = [x[1] * x[2], x[2] * x[3]]
    init     = [0., 0., 0.]
    st_mod   = DiffusionStateModel(f, g, init)
    ob_mod   = DiffusionObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    N        = 10
    alpha    = 0.
    filter   = BPF(st_mod, ob_mod, N, alpha)
    dY       = [0.1, 0.2]
    dt       = 0.01
    assimilate!(dY, state, filter, dt)
    print(".")
    @test state.ensemble.positions == oldstate.ensemble.positions
    print(".")
    @test state.ensemble.weights != oldstate.ensemble.weights

    h(x)     = [x[1]^2, x[2]^2]
    ob_mod   = CountingObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    filter   = BPF(st_mod, ob_mod, N, alpha)
    dN       = [1, 0]
    dt       = 0.01
    state    = BPFState(ens)
    oldstate = deepcopy(state)
    assimilate!(dN, state, filter, dt)
    print(".")
    @test state.ensemble.positions == oldstate.ensemble.positions
    print(".")
    @test state.ensemble.weights != oldstate.ensemble.weights
    println("DONE.")

    print("  method update!")
    state  = initialize(filter)
    state2 = initialize(filter)
    Random.seed!(0)
    update!(state, filter, [0.01, -0.02], 0.01)
    Random.seed!(0)
    propagate!(state2, filter, 0.01)
    assimilate!([0.01, -0.02], state2, filter, 0.01)
    print(".")
    @test state.ensemble.positions ≈ state2.ensemble.positions
    print(".")
    @test state.ensemble.weights ≈ state2.ensemble.weights
    println("DONE.")

end; #BPF.jl