text:
# Note: be sure to add the JSON package:
# julia> ] add JSON

using JSON

mutable struct ModelData
    data::Dict{AbstractString,Any}
    sim::Dict{AbstractString,Any}

    neuron::Dict{String,Any}
    dendrite::Dict{String,Any}
    compartment::Dict{String,Any}

    synapses::Array{Any,1}
    synapse::Dict{String,Any}
    active_synapse::Int64

    loaded::Bool
    app_data_loaded::Bool
    sim_data_loaded::Bool

    # Neuron parameters
    threshold::String

    bug::Bool

    sim_changed::Bool
    app_changed::Bool

    # True = apply synapse property to all synapses
    # not just the active synapse.
    apply_to_all::Bool

    function ModelData()
        o = new()

        o.threshold = "\0"^16
        o.loaded = false
        o.bug = true
        o.sim_changed = false
        o.app_changed = false
        o.sim_data_loaded = false
        o.app_data_loaded = false
        o.apply_to_all = true
        o
    end
end

function load!(model::ModelData, file::String)
    println("Loading: ", file)

    # json = ""
    # lines = readlines(file)
    # for line in lines
    #     # strip any comments out
    #     range = findfirst("//", line)
    #     if range ≠ nothing
    #         line = String(SubString(line, 1:range[1] - 1))
    #     end
    #     json = json * line
    # end

    # json = open(file) do fd
    #     read(fd, String)
    # end

    # data = JSON.parse(json)

    data = JSON.parsefile(file)

    model.data = data
    model.app_data_loaded = true;
end

# Called when the app.jl is closing
function unload(model::ModelData)
    if is_app_changed(model)
        save(model, "../data/app.json")
    end

    if is_sim_changed(model)
        save_sim(model)
    end
end

function save(model::ModelData, file::String)
    println("Saving ", file)
    open(file, "w") do f
        JSON.print(f, model.data, 2)
    end
end

function is_loaded(model::ModelData)
    model.sim_data_loaded
end

function is_sim_loaded(model::ModelData)
    model.loaded = model.app_data_loaded && model.sim_data_loaded
    model.loaded
end

function is_app_changed(model::ModelData)
    model.app_changed
end

function is_sim_changed(model::ModelData)
    model.sim_changed
end

function is_changed(model::ModelData)
    is_app_changed(model) || is_sim_changed(model)
end

function load_sim!(model::ModelData)
    sim = simulation(model)
    path = data_path(model)
    file = path * sim
    println("Loading sim model... ", file)

    # json = open(file) do fd
    #     read(fd, String)
    # end

    # json = ""
    # lines = readlines(file)
    # for line in lines
    #     # strip any comments out
    #     range = findfirst("//", line)
    #     if range ≠ nothing
    #         line = String(SubString(line, 1:range[1] - 1))
    #     end
    #     json = json * line
    # end

    # data = JSON.parse(json)

    data = JSON.parsefile(file)

    # println(data)

    model.sim = data;
    if model.sim == nothing
        println("##### WARNING #####")
        println("sim json node not found")
    end

    model.neuron = model.sim["Neuron"]
    model.active_synapse = model.sim["ActiveSynapse"]

    # TODO change json to be array of dendrites
    dendrites = model.neuron["Dendrites"]
    if dendrites == nothing
        println("##### WARNING #####")
        println("dendrites json node not found")
    end
    model.dendrite = dendrites

    compartments = model.dendrite["Compartments"]
    if compartments == nothing
        println("##### WARNING #####")
        println("compartments json node not found")
    end
    model.compartment = compartments[1]
    model.synapses = model.compartment["Synapses"]
    model.synapse = model.synapses[model.active_synapse]

    println("Sim model loaded")
    model.sim_data_loaded = true;
end

function save_sim(model::ModelData)
    sim = simulation(model)
    path = data_path(model)
    file = path * sim
    println("Saving sim data: ", file)
    open(file, "w") do f
        JSON.print(f, model.sim, 2)
    end
end

function set_changed(model::ModelData, what::Int64)
    if what == 0
        # println("App data changed")
        # print_trace(stacktrace(), "")
        model.app_changed = true
    elseif what == 1
        # println("Sim data changed")
        # print_trace(stacktrace(), "")
        model.sim_changed = true
    end
end
# --------------------------------------------------------------
# Field utils
# --------------------------------------------------------------
function prep_field(v, pad::Int64 = 16)
    string(v) * " "^pad
end

function strip_null(v::String)
    range = findfirst("\0", v)
    if range ≠ nothing
        v = String(SubString(v, 1:(range[1] - 1)))
    end
    v
end

# --------------------------------------------------------------
# Calculated Values
# --------------------------------------------------------------
function total_simulation_time(model::ModelData)
    duration = Model.duration(model)
    # span = Model.span_time(model)
    duration #* span
end

# Span is = Duration
function span_time(model::ModelData)
    duration = Model.duration(model)
    spans = Model.spans(model)
    Int64(duration / spans)
end

# --------------------------------------------------------------
# Application Setter/Getters
# --------------------------------------------------------------
function data_path(model::ModelData)
    model.data["DataPath"]
end

function simulation(model::ModelData)
    model.data["Simulation"]
end
function set_simulation!(model::ModelData, v::String)
    model.data["Simulation"] = v
    set_changed(model, 0)
end

# Duration is how long a Span is in units of time.
# Thus the total simulation time = Duration * Spans
function duration(model::ModelData)
    model.data["Duration"]
end
function set_duration!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["Duration"] = parse(Int64, v)
    set_changed(model, 0)
end

function range_start(model::ModelData)
    model.data["RangeStart"]
end
function set_range_start!(model::ModelData, v::String)
    v = strip_null(v)
    pv = parse(Int64, v)
    if pv != model.data["RangeStart"]
        model.data["RangeStart"] = pv
        set_changed(model, 0)
    end
end
function set_range_start!(model::ModelData, v::Int64)
    if v != model.data["RangeStart"]
        model.data["RangeStart"] = v
        set_changed(model, 0)
    end
end

function range_end(model::ModelData)
    model.data["RangeEnd"]
end
function set_range_end!(model::ModelData, v::String)
    v = strip_null(v)
    pv = parse(Int64, v)
    if pv != model.data["RangeEnd"]
        model.data["RangeEnd"] = pv
        set_changed(model, 0)
    end
end
function set_range_end!(model::ModelData, v::Int64)
    if v != model.data["RangeEnd"]
        model.data["RangeEnd"] = v
        set_changed(model, 0)
    end
end

function scroll(model::ModelData)
    model.data["Scroll"]
end
function set_scroll!(model::ModelData, v::Float64)
    if v != model.data["Scroll"]
        model.data["Scroll"] = v
        set_changed(model, 0)
    end
end

function spans(model::ModelData)
    model.data["Spans"]
end
function set_spans!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["Spans"] = parse(Int64, v)
    set_changed(model, 0)
end

function time_scale(model::ModelData)
    model.data["TimeScale"]
end
function set_time_scale!(model::ModelData, v::String)
    v = strip_null(v)
    set_time_scale!(model, parse(Int64, v))
    set_changed(model, 0)
end
function set_time_scale!(model::ModelData, v::Int64)
    model.data["TimeScale"] = v
    set_changed(model, 0)
end

function data_output_path(model::ModelData)
    model.data["DataOutputPath"]
end
function set_data_output_path!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["DataOutputPath"] = v
    set_changed(model, 0)
end

function poisson_files(model::ModelData)
    model.data["OutputPoissonFiles"]
end
function set_poisson_files!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputPoissonFiles"] = v
    set_changed(model, 0)
end

function output_stimulus_files(model::ModelData)
    model.data["OutputStimulusFiles"]
end
function set_output_stimulus_files!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputStimulusFiles"] = v
    set_changed(model, 0)
end

function source_stimulus(model::ModelData)
    model.data["SourceStimulus"]
end
function set_source_stimulus!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["SourceStimulus"] = v
    set_changed(model, 0)
end

function output_cell_spikes(model::ModelData)
    model.data["OutputSomaSpikeFiles"]
end
function set_output_cell_spikes!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSomaSpikeFiles"] = v
    set_changed(model, 0)
end

function output_soma_apFast(model::ModelData)
    model.data["OutputSomaAPFastFiles"]
end
function set_output_soma_apFast!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSomaAPFastFiles"] = v
    set_changed(model, 0)
end

function output_soma_apSlow(model::ModelData)
    model.data["OutputSomaAPSlowFiles"]
end
function set_output_soma_apSlow!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSomaAPSlowFiles"] = v
    set_changed(model, 0)
end

function output_soma_psp(model::ModelData)
    model.data["OutputSomaPSPFiles"]
end
function set_output_soma_psp!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSomaPSPFiles"] = v
    set_changed(model, 0)
end

function output_synapse_weights(model::ModelData)
    model.data["OutputSynapseWeightFiles"]
end
function set_synapse_weights!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSynapseWeightFiles"] = v
    set_changed(model, 0)
end

function output_synapse_surge(model::ModelData)
    model.data["OutputSynapseSurgeFiles"]
end
function set_synapse_surge!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSynapseSurgeFiles"] = v
    set_changed(model, 0)
end

function output_synapse_psp(model::ModelData)
    model.data["OutputSynapsePspFiles"]
end
function set_synapse_psp!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSynapsePspFiles"] = v
    set_changed(model, 0)
end

function output_synapse_input_spike(model::ModelData)
    model.data["OutputSynapseSpikeFiles"]
end
function set_synapse_input_spike!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputSynapseSpikeFiles"] = v
    set_changed(model, 0)
end

function output_dendrite_avg(model::ModelData)
    model.data["OutputDendriteAvgFiles"]
end
function set_dendrite_avg!(model::ModelData, v::String)
    v = strip_null(v)
    model.data["OutputDendriteAvgFiles"] = v
    set_changed(model, 0)
end

# --------------------------------------------------------------
# Simulation Setter/Getters
# --------------------------------------------------------------
# Streams
# Firing rate = spikes over an interval of time or
# Poisson events per interval of time.
# For example, spikes in a 1 sec span.
# A firing rate in unit/ms, for example, 0.2 in 1ms (0.2/1)
# or 200 in 1sec (200/1000ms)
function firing_rate(model::ModelData)
    model.sim["Firing_Rate"]
end
function set_firing_rate!(model::ModelData, v::String)
    v = strip_null(v)
    set_firing_rate!(model, parse(Float64, v))
end
function set_firing_rate!(model::ModelData, v::Float64)
    model.sim["Firing_Rate"] = v
    set_changed(model, 1)
end

function synapses(model::ModelData)
    model.sim["Synapses"]
end
function synapses!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["Synapses"] = parse(Int64, v)
    set_changed(model, 1)
end

function active_synapse(model::ModelData)
    model.sim["ActiveSynapse"]
end
function set_active_synapse!(model::ModelData, v::String)
    v = strip_null(v)
    set_active_synapse!(model, parse(Int64, v))
end
function set_active_synapse!(model::ModelData, v::Int64)
    if v == active_synapse(model)
        return
    end

    model.sim["ActiveSynapse"] = v
    model.active_synapse = model.sim["ActiveSynapse"]

    # Find the correct synapse node
    for synapse in model.synapses
        if synapse["id"] == model.active_synapse
            model.synapse = model.synapses[model.active_synapse]
            break
        end
    end

    set_changed(model, 1)
end

function percent_excititory_synapses(model::ModelData)
    model.sim["PercentOfExcititorySynapses"]
end
function set_percent_excititory_synapses!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["PercentOfExcititorySynapses"] = parse(Int64, v)
    set_changed(model, 1)
end
function set_percent_excititory_synapses!(model::ModelData, v::Float64)
    model.sim["PercentOfExcititorySynapses"] = v
    set_changed(model, 1)
end

# If Hertz = 0 then stimulus is distributed as poisson.
# Hertz is = cycles per second (or 1000ms per second)
# 10Hz = 10 applied in 1000ms or every 100ms = 1000/10Hz
# This means a stimulus is generated every 100ms which also means the
# Inter-spike-interval (ISI) is fixed at 100ms
function hertz(model::ModelData)
    model.sim["Hertz"]
end
function set_hertz!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["Hertz"] = parse(Float64, v)
    set_changed(model, 1)
end

# Shrinks or expands ISI for stimulus
function stimulus_scaler(model::ModelData)
    model.sim["StimulusScaler"]
end
function set_stimulus_scaler!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["StimulusScaler"] = parse(Int64, v)
    set_changed(model, 1)
end

function poisson_pattern_min(model::ModelData)
    model.sim["Poisson_Pattern_min"]
end
function set_poisson_pattern_min!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["Poisson_Pattern_min"] = parse(Int64, v)
    set_changed(model, 1)
end

function poisson_pattern_max(model::ModelData)
    model.sim["Poisson_Pattern_max"]
end
function set_poisson_pattern_max!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["Poisson_Pattern_max"] = parse(Int64, v)
    set_changed(model, 1)
end

function poisson_pattern_spread(model::ModelData)
    model.sim["Poisson_Pattern_spread"]
end
function set_poisson_pattern_spread!(model::ModelData, v::String)
    v = strip_null(v)
    model.sim["Poisson_Pattern_spread"] = parse(Int64, v)
    set_changed(model, 1)
end

# --------------------------------------------------------------
# Neuron Setter/Getters
# --------------------------------------------------------------
function refractory_period(model::ModelData)
    model.neuron["RefractoryPeriod"]
end
function set_refractory_period!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["RefractoryPeriod"] = parse(Float64, v)
    set_changed(model, 1)
end

function ap_max(model::ModelData)
    model.neuron["APMax"]
end
function set_ap_max!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["APMax"] = parse(Float64, v)
    set_changed(model, 1)
end

function threshold(model::ModelData)
    model.neuron["Threshold"]
end
function set_threshold!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["Threshold"] = parse(Float64, v)
    set_changed(model, 1)
end

function fast_surge(model::ModelData)
    model.neuron["nFastSurge"]
end
function set_fast_surge!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["nFastSurge"] = parse(Float64, v)
    set_changed(model, 1)
end

function slow_surge(model::ModelData)
    model.neuron["nSlowSurge"]
end
function set_slow_surge!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["nSlowSurge"] = parse(Float64, v)
    set_changed(model, 1)
end

function slow_surge(model::ModelData)
    model.neuron["nSlowSurge"]
end
function set_slow_surge!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["nSlowSurge"] = parse(Float64, v)
    set_changed(model, 1)
end

function tao(model::ModelData)
    model.neuron["ntao"]
end
function set_tao!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["ntao"] = parse(Float64, v)
    set_changed(model, 1)
end

function tao_j(model::ModelData)
    model.neuron["ntaoJ"]
end
function set_tao_j!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["ntaoJ"] = parse(Float64, v)
    set_changed(model, 1)
end

function tao_s(model::ModelData)
    model.neuron["ntaoS"]
end
function set_tao_s!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["ntaoS"] = parse(Float64, v)
    set_changed(model, 1)
end

function w_max(model::ModelData)
    model.neuron["wMax"]
end
function set_w_max!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["wMax"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_w_max!(model::ModelData, v::Float64)
    model.neuron["wMax"] = v
    set_changed(model, 1)
end

function w_min(model::ModelData)
    model.neuron["wMin"]
end
function set_w_min!(model::ModelData, v::String)
    v = strip_null(v)
    model.neuron["wMin"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_w_min!(model::ModelData, v::Float64)
    model.neuron["wMin"] = v
    set_changed(model, 1)
end

# --------------------------------------------------------------
# Dendrite Setter/Getters
# --------------------------------------------------------------
function tao_eff(model::ModelData)
    model.dendrite["taoEff"]
end
function set_tao_eff!(model::ModelData, v::String)
    v = strip_null(v)
    model.dendrite["taoEff"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_tao_eff!(model::ModelData, v::Float64)
    model.dendrite["taoEff"] = v
    set_changed(model, 1)
end

function dendrite_length(model::ModelData)
    model.dendrite["length"]
end
function set_dendrite_length!(model::ModelData, v::String)
    v = strip_null(v)
    model.dendrite["length"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_dendrite_length!(model::ModelData, v::Float64)
    model.dendrite["length"] = v
    set_changed(model, 1)
end

function dendrite_min_psp(model::ModelData)
    model.dendrite["MinPSPValue"]
end
function set_dendrite_min_psp!(model::ModelData, v::String)
    v = strip_null(v)
    model.dendrite["MinPSPValue"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_dendrite_min_psp!(model::ModelData, v::Float64)
    model.dendrite["MinPSPValue"] = v
    set_changed(model, 1)
end

# --------------------------------------------------------------
# Compartments Setter/Getters
# --------------------------------------------------------------
function weight_max(model::ModelData)
    model.compartment["WeightMax"]
end
function set_weight_max!(model::ModelData, v::String)
    v = strip_null(v)
    model.compartment["WeightMax"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_weight_max!(model::ModelData, v::Float64)
    model.compartment["WeightMax"] = v
    set_changed(model, 1)
end

function weight_divisor(model::ModelData)
    model.compartment["WeightDivisor"]
end
function set_weight_divisor!(model::ModelData, v::String)
    v = strip_null(v)
    model.compartment["WeightDivisor"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_weight_divisor!(model::ModelData, v::Float64)
    model.compartment["WeightDivisor"] = v
    set_changed(model, 1)
end

# --------------------------------------------------------------
# Synapse Setter/Getters
# --------------------------------------------------------------
function id(model::ModelData)
    model.synapse["id"]
end
function set_id!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["id"] = parse(Float64, v)
    set_changed(model, 1)
end

function alpha(model::ModelData)
    model.synapse["alpha"]
end
# Set alpha only on the active synapse
function set_alpha!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["alpha"] = parse(Float64, v)
    set_changed(model, 1)
end
# Set alpha on all synapses
function set_alphas!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["alpha"] = f
    end
    set_changed(model, 1)
end

function ama(model::ModelData)
    model.synapse["ama"]
end
function set_ama!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["ama"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_amas!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["ama"] = f
    end
    set_changed(model, 1)
end

function amb(model::ModelData)
    model.synapse["amb"]
end
function set_amb!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["amb"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_ambs!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["amb"] = f
    end
    set_changed(model, 1)
end

function lambda(model::ModelData)
    model.synapse["lambda"]
end
function set_lambda!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["lambda"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_lambdas!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["lambda"] = f
    end
    set_changed(model, 1)
end

function learning_rate_fast(model::ModelData)
    model.synapse["learningRateFast"]
end
function set_learning_rate_fast!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["learningRateFast"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_learning_rate_fasts!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["learningRateFast"] = f
    end
    set_changed(model, 1)
end

function learning_rate_slow(model::ModelData)
    model.synapse["learningRateSlow"]
end
function set_learning_rate_slow!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["learningRateSlow"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_learning_rate_slows!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["learningRateSlow"] = f
    end
    set_changed(model, 1)
end

function mu(model::ModelData)
    model.synapse["mu"]
end
function set_mu!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["mu"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_mus!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["mu"] = f
    end
    set_changed(model, 1)
end

function taoI(model::ModelData)
    model.synapse["taoI"]
end
function set_taoI!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["taoI"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_taoIs!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["taoI"] = f
    end
    set_changed(model, 1)
end

function taoN(model::ModelData)
    model.synapse["taoN"]
end
function set_taoN!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["taoN"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_taoNs!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["taoN"] = f
    end
    set_changed(model, 1)
end

function taoP(model::ModelData)
    model.synapse["taoP"]
end
function set_taoP!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["taoP"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_taoPs!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["taoP"] = f
    end
    set_changed(model, 1)
end

# Distance from soma
function distance(model::ModelData)
    model.synapse["distance"]
end
function set_distance!(model::ModelData, v::String)
    v = strip_null(v)
    model.synapse["distance"] = parse(Float64, v)
    set_changed(model, 1)
end
function set_distances!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["distance"] = f
    end
    set_changed(model, 1)
end

function weight(model::ModelData)
    model.synapse["w"]
end
function set_weight!(model::ModelData, v::Float64)
    model.synapse["w"] = v
    set_changed(model, 1)
end
function set_weights!(model::ModelData, v::String)
    v = strip_null(v)
    f = parse(Float64, v)
    for synapse in model.synapses
        synapse["w"] = f
    end
    set_changed(model, 1)
end


text:
# tests for the encounter type and related functions.

# using Base.Test

# push!(LOAD_PATH, ".")
# include("includes.jl")

# encounter tests
function test_encounter_constructor()
    e = Encounter()
end

# update! tests
function test_update_actions()
    e = Encounter()
    uav_actions = UAVActions(2)
    uav_actions.actions = deg2rad([-10., -10.])
    utilities = [[[1.,2.] [3.,4.]] [[1.,2.] [3.,4.]]]
    belief_state = spzeros(10000, 1)
    belief_state[1,1] = 10
    belief_states = [belief_state belief_state]

    update!(e, uav_actions, utilities, belief_states)
    @test e.actions == reshape(uav_actions.actions, (2, 1))
    update!(e, uav_actions, utilities, belief_states)
    @test e.actions == [uav_actions.actions uav_actions.actions]
end

function test_update_state()
    e = Encounter()
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(zeros(Float64, 8)...)
    s = State([uav_state_own; uav_state_int])
    update!(e, s)
    @test e.paths[1].x_coordinates == [1]
    @test e.paths[1].y_coordinates == [1]
    @test e.paths[2].x_coordinates == [0]
    @test e.paths[2].y_coordinates == [0]
    update!(e, s)
    @test e.paths[1].x_coordinates == [1, 1]
    @test e.paths[1].y_coordinates == [1, 1]
    @test e.paths[2].x_coordinates == [0, 0]
    @test e.paths[2].y_coordinates == [0, 0]
end
tic()
# Environment tests
test_encounter_constructor()

# update! tests
test_update_actions()
test_update_state()
toc()

ext: jl


text:
module ARules
using DataFrames

export Node, Rule, apriori, frequent, unique_items, shownodes,
       # these below are for pilotting
       has_children, younger_siblings, growtree!, randstr, occurence, gen_support_dict, frequent_item_tree,
       gen_rules!, gen_rules, gen_node_rules, occurrence

include("frequent_itemset_tree.jl")
include("utils.jl")
include("rule_generation.jl")

# package code goes here

end # module


text:
using GeometricEquations
using Test

include("functions.jl")
include("initial_conditions.jl")

@testset "$(rpad("Stochastic Differential Equations (SDE)",80))" begin

    λ  = 2.
    μ  = 1.

    sde_v_params = (t, q, v) -> sde_v(λ, t, q, v)
    sde_B_params = (t, q, v) -> sde_B(μ, t, q, v)

    sde  = SDE(1, 1, sde_v_params, sde_B_params, t₀, x₀)
    sde1 = SDE(1, 1, sde_v_params, sde_B_params, x₀)
    sde2 = SDE(1, 3, sde_v_params, sde_B_params, x₀)
    sde3 = SDE(1, 1, sde_v_params, sde_B_params, x₁ₛ)

    @test ndims(sde) == 2
    @test periodicity(sde) == zero(x₀)
    @test functions(sde) == NamedTuple{(:v,:B)}((sde_v_params, sde_B_params))

    @test sde == sde1
    @test sde != sde2
    @test sde != sde3

    @test hash(sde) == hash(sde1)
    @test hash(sde) != hash(sde2)
    @test hash(sde) != hash(sde3)

    @test sde1.d == 2
    @test sde2.d == 2
    @test sde3.d == 2

    @test sde1.ns == 1
    @test sde2.ns == 3
    @test sde3.ns == 1

    @test nsamples(sde1) == 1
    @test nsamples(sde2) == 1
    @test nsamples(sde3) == 3

    @test sde == similar(sde, t₀, x₀, 1)
    @test sde == similar(sde, t₀, x₀)
    @test sde == similar(sde, x₀, 1)
    @test sde == similar(sde, x₀)

end

@testset "$(rpad("Partitioned Stochastic Differential Equations (PSDE)",80))" begin

    psde  = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, t₀, q₀, p₀)
    psde1 = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, q₀, p₀)
    psde2 = PSDE(1, 3, psde_v, psde_f, psde_B, psde_G, q₀, p₀)
    psde3 = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, q₁ₛ, p₁ₛ)

    @test ndims(psde) == 1
    @test periodicity(psde) == zero(q₀)
    @test functions(psde) == NamedTuple{(:v,:f,:B,:G)}((psde_v,psde_f,psde_B,psde_G))

    @test psde == psde1
    @test psde != psde2
    @test psde != psde3

    @test hash(psde) == hash(psde1)
    @test hash(psde) != hash(psde2)
    @test hash(psde) != hash(psde3)

    @test psde1.d == 1
    @test psde2.d == 1
    @test psde3.d == 1

    @test psde1.ns == 1
    @test psde2.ns == 3
    @test psde3.ns == 1

    @test nsamples(psde1) == 1
    @test nsamples(psde2) == 1
    @test nsamples(psde3) == 3

    @test psde == similar(psde, t₀, q₀, p₀, 1)
    @test psde == similar(psde, q₀, p₀, 1)

    @test psde != similar(psde, t₀, q₁ₛ, p₁ₛ)
    @test psde != similar(psde, q₁ₛ, p₁ₛ)

end

@testset "$(rpad("Split Partitioned Stochastic Differential Equations (SPSDE)",80))" begin

    spsde  = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, t₀, q₀, p₀)
    spsde1 = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₀, p₀)
    spsde2 = SPSDE(1, 3, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₀, p₀)
    spsde3 = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₁ₛ, p₁ₛ)

    @test ndims(spsde) == 1
    @test periodicity(spsde) == zero(q₀)
    @test functions(spsde) == NamedTuple{(:v,:f1,:f2,:B,:G1,:G2)}((spsde_v,spsde_f1,spsde_f2,spsde_B,spsde_G1,spsde_G2))

    @test spsde == spsde1
    @test spsde != spsde2
    @test spsde != spsde3

    @test hash(spsde) == hash(spsde1)
    @test hash(spsde) != hash(spsde2)
    @test hash(spsde) != hash(spsde3)

    @test spsde1.d == 1
    @test spsde2.d == 1
    @test spsde3.d == 1

    @test spsde1.ns == 1
    @test spsde2.ns == 3
    @test spsde3.ns == 1

    @test nsamples(spsde1) == 1
    @test nsamples(spsde2) == 1
    @test nsamples(spsde3) == 3

    @test spsde == similar(spsde, t₀, q₀, p₀, 1)
    @test spsde == similar(spsde, q₀, p₀, 1)

    @test spsde != similar(spsde, t₀, q₁ₛ, p₁ₛ)
    @test spsde != similar(spsde, q₁ₛ, p₁ₛ)

end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie
GLMakie.activate!() # HIDE
let
    t = 0:0.1:15
    fig, _ = lines(sin.(t), cos.(t), t/4, color = t/4, linewidth = 4,
        colormap = :plasma)
    fig
    save(joinpath(@__DIR__, "output", "line3D.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((c+d*x)^3*csch(a+b*x), x) == :(-6*b^-4*d^3*Polylogarithms.polylog(4,-1*exp(a+b*x))+-2*b^-1*(c+d*x)^3*arctanh(exp(a+b*x))+6*b^-4*d^3*Polylogarithms.polylog(4,exp(a+b*x))+-6*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,exp(a+b*x))+-3*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,-1*exp(a+b*x))+3*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(a+b*x))+6*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,-1*exp(a+b*x)))
@test integrate((c+d*x)^2*csch(a+b*x), x) == :(-2*b^-1*(c+d*x)^2*arctanh(exp(a+b*x))+-2*b^-3*d^2*Polylogarithms.polylog(3,exp(a+b*x))+2*b^-3*d^2*Polylogarithms.polylog(3,-1*exp(a+b*x))+-2*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,-1*exp(a+b*x))+2*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,exp(a+b*x)))
@test integrate((c+d*x)*csch(a+b*x), x) == :(d*b^-2*Polylogarithms.polylog(2,exp(a+b*x))+-1*d*b^-2*Polylogarithms.polylog(2,-1*exp(a+b*x))+-1*b^-1*(2c+2*d*x)*arctanh(exp(a+b*x)))
@test integrate((c+d*x)^3*csch(a+b*x)^2, x) == :(-1*b^-1*(c+d*x)^3+-1*b^-1*(c+d*x)^3*coth(a+b*x)+-3//2*b^-4*d^3*Polylogarithms.polylog(3,exp(2a+2*b*x))+3*d*b^-2*(c+d*x)^2*log(1+-1*exp(2a+2*b*x))+3*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(2,exp(2a+2*b*x)))
@test integrate((c+d*x)^2*csch(a+b*x)^2, x) == :(-1*b^-1*(c+d*x)^2+b^-3*d^2*Polylogarithms.polylog(2,exp(2a+2*b*x))+-1*b^-1*(c+d*x)^2*coth(a+b*x)+2*d*b^-2*(c+d*x)*log(1+-1*exp(2a+2*b*x)))
@test integrate(csch(a+b*x)^2*(c+d*x), x) == :(d*b^-2*log(sinh(a+b*x))+-1*b^-1*(c+d*x)*coth(a+b*x))
@test integrate((c+d*x)^3*csch(a+b*x)^3, x) == :(b^-1*(c+d*x)^3*arctanh(exp(a+b*x))+-3*b^-4*d^3*Polylogarithms.polylog(2,-1*exp(a+b*x))+-3*b^-4*d^3*Polylogarithms.polylog(4,exp(a+b*x))+3*b^-4*d^3*Polylogarithms.polylog(2,exp(a+b*x))+3*b^-4*d^3*Polylogarithms.polylog(4,-1*exp(a+b*x))+-6*b^-3*d^2*(c+d*x)*arctanh(exp(a+b*x))+-3*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,-1*exp(a+b*x))+3*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,exp(a+b*x))+-3//2*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(a+b*x))+-3//2*d*b^-2*(c+d*x)^2*csch(a+b*x)+-1//2*b^-1*(c+d*x)^3*coth(a+b*x)*csch(a+b*x)+3//2*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,-1*exp(a+b*x)))
@test integrate((c+d*x)^2*csch(a+b*x)^3, x) == :(b^-1*(c+d*x)^2*arctanh(exp(a+b*x))+b^-3*d^2*Polylogarithms.polylog(3,exp(a+b*x))+-1*b^-3*d^2*Polylogarithms.polylog(3,-1*exp(a+b*x))+-1*b^-3*d^2*arctanh(cosh(a+b*x))+d*b^-2*(c+d*x)*Polylogarithms.polylog(2,-1*exp(a+b*x))+-1*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,exp(a+b*x))+-1*d*b^-2*(c+d*x)*csch(a+b*x)+-1//2*b^-1*(c+d*x)^2*coth(a+b*x)*csch(a+b*x))
@test integrate(csch(a+b*x)^3*(c+d*x), x) == :(b^-1*(c+d*x)*arctanh(exp(a+b*x))+(1/2)*d*b^-2*Polylogarithms.polylog(2,-1*exp(a+b*x))+-1//2*d*b^-2*Polylogarithms.polylog(2,exp(a+b*x))+-1//2*d*b^-2*csch(a+b*x)+-1//2*b^-1*(c+d*x)*coth(a+b*x)*csch(a+b*x))
@test integrate(x*csch(x)^-3//2+1//3*x*csch(x)^(1/2), x) == :(-4//9*csch(x)^-3//2+2//3*x*csch(x)^-1//2*cosh(x))
@test integrate(x*csch(x)^-5//2+3//5*x*csch(x)^-1//2, x) == :(-4//25*csch(x)^-5//2+2//5*x*csch(x)^-3//2*cosh(x))
@test integrate(x*csch(x)^-7//2+-5//21*x*csch(x)^(1/2), x) == :(-4//49*csch(x)^-7//2+20//63*csch(x)^-3//2+-10//21*x*csch(x)^-1//2*cosh(x)+2//7*x*csch(x)^-5//2*cosh(x))
@test integrate(x^2*csch(x)^-3//2+1//3*x^2*csch(x)^(1/2), x) == :(-8//9*x*csch(x)^-3//2+16//27*csch(x)^-1//2*cosh(x)+2//3*x^2*csch(x)^-1//2*cosh(x)+-16//27*im*(im*sinh(x))^(1/2)*csch(x)^(1/2)*Elliptic.F(1//4*pi+-1//2*im*x,2))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^3*cosh(c+d*x), x) == :(a^-1*d^-1*(e+f*x)^3*sinh(c+d*x)+-6*a^-1*d^-4*f^3*cosh(c+d*x)+1//4*b*a^-2*f^-1*(e+f*x)^4+-1*b*a^-2*d^-1*(e+f*x)^3*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-2*d^-1*(e+f*x)^3*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-2*d^-4*f^3*Polylogarithms.polylog(4,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-2*d^-4*f^3*Polylogarithms.polylog(4,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-3*f*a^-1*d^-2*(e+f*x)^2*cosh(c+d*x)+6*a^-1*d^-3*f^2*(e+f*x)*sinh(c+d*x)+-3*b*f*a^-2*d^-2*(e+f*x)^2*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-3*b*f*a^-2*d^-2*(e+f*x)^2*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+6*b*a^-2*d^-3*f^2*(e+f*x)*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+6*b*a^-2*d^-3*f^2*(e+f*x)*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^2*cosh(c+d*x), x) == :(a^-1*d^-1*(e+f*x)^2*sinh(c+d*x)+2*a^-1*d^-3*f^2*sinh(c+d*x)+1//3*b*a^-2*f^-1*(e+f*x)^3+-1*b*a^-2*d^-1*(e+f*x)^2*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-2*d^-1*(e+f*x)^2*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*f*a^-1*d^-2*(e+f*x)*cosh(c+d*x)+2*b*a^-2*d^-3*f^2*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+2*b*a^-2*d^-3*f^2*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*b*f*a^-2*d^-2*(e+f*x)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*b*f*a^-2*d^-2*(e+f*x)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)*cosh(c+d*x), x) == :(a^-1*d^-1*(e+f*x)*sinh(c+d*x)+(1/2)*b*a^-2*f^-1*(e+f*x)^2+-1*f*a^-1*d^-2*cosh(c+d*x)+-1*b*f*a^-2*d^-2*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*f*a^-2*d^-2*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-2*d^-1*(e+f*x)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-2*d^-1*(e+f*x)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*cosh(c+d*x), x) == :(a^-1*d^-1*sinh(c+d*x)+-1*b*a^-2*d^-1*log(b+a*sinh(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^3*cosh(c+d*x)^2, x) == :(1//8*a^-1*f^-1*(e+f*x)^4+-3//8*a^-1*d^-4*f^3*cosh(c+d*x)^2+1//4*a^-3*b^2*f^-1*(e+f*x)^4+3//8*a^-1*d^-2*f^3*x^2+(1/2)*a^-1*d^-1*(e+f*x)^3*cosh(c+d*x)*sinh(c+d*x)+-1*b*a^-2*d^-1*(e+f*x)^3*cosh(c+d*x)+6*b*a^-2*d^-4*f^3*sinh(c+d*x)+-3//4*f*a^-1*d^-2*(e+f*x)^2*cosh(c+d*x)^2+3//4*e*x*a^-1*d^-2*f^2+b*a^-3*d^-1*(e+f*x)^3*(a^2+b^2)^(1/2)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-3*d^-1*(e+f*x)^3*(a^2+b^2)^(1/2)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-3*d^-4*f^3*(a^2+b^2)^(1/2)*Polylogarithms.polylog(4,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-2*d^-3*f^2*(e+f*x)*cosh(c+d*x)+3*b*f*a^-2*d^-2*(e+f*x)^2*sinh(c+d*x)+6*b*a^-3*d^-4*f^3*(a^2+b^2)^(1/2)*Polylogarithms.polylog(4,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+3//4*a^-1*d^-3*f^2*(e+f*x)*cosh(c+d*x)*sinh(c+d*x)+-6*b*a^-3*d^-3*f^2*(a^2+b^2)^(1/2)*(e+f*x)*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-3*b*f*a^-3*d^-2*(e+f*x)^2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+3*b*f*a^-3*d^-2*(e+f*x)^2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+6*b*a^-3*d^-3*f^2*(a^2+b^2)^(1/2)*(e+f*x)*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^2*cosh(c+d*x)^2, x) == :(1//6*a^-1*f^-1*(e+f*x)^3+1//3*a^-3*b^2*f^-1*(e+f*x)^3+1//4*x*a^-1*d^-2*f^2+(1/2)*a^-1*d^-1*(e+f*x)^2*cosh(c+d*x)*sinh(c+d*x)+-1*b*a^-2*d^-1*(e+f*x)^2*cosh(c+d*x)+-2*b*a^-2*d^-3*f^2*cosh(c+d*x)+-1//2*f*a^-1*d^-2*cosh(c+d*x)^2*(e+f*x)+1//4*a^-1*d^-3*f^2*cosh(c+d*x)*sinh(c+d*x)+b*a^-3*d^-1*(e+f*x)^2*(a^2+b^2)^(1/2)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-3*d^-1*(e+f*x)^2*(a^2+b^2)^(1/2)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*b*a^-3*d^-3*f^2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+2*b*f*a^-2*d^-2*(e+f*x)*sinh(c+d*x)+2*b*a^-3*d^-3*f^2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*b*f*a^-3*d^-2*(a^2+b^2)^(1/2)*(e+f*x)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+2*b*f*a^-3*d^-2*(a^2+b^2)^(1/2)*(e+f*x)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*cosh(c+d*x)^2*(e+f*x), x) == :((1/2)*e*x*a^-1+1//4*f*a^-1*x^2+e*x*a^-3*b^2+(1/2)*f*a^-3*b^2*x^2+-1//4*f*a^-1*d^-2*cosh(c+d*x)^2+b*f*a^-2*d^-2*sinh(c+d*x)+(1/2)*a^-1*d^-1*(e+f*x)*cosh(c+d*x)*sinh(c+d*x)+-1*b*a^-2*d^-1*(e+f*x)*cosh(c+d*x)+b*f*a^-3*d^-2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+b*a^-3*d^-1*(a^2+b^2)^(1/2)*(e+f*x)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*f*a^-3*d^-2*(a^2+b^2)^(1/2)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-3*d^-1*(a^2+b^2)^(1/2)*(e+f*x)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*cosh(c+d*x)^2, x) == :((1/2)*x*a^-3*(a^2+2*b^2)+-1//2*a^-2*d^-1*(2b+-1*a*sinh(c+d*x))*cosh(c+d*x)+2*b*a^-3*d^-1*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^3*cosh(c+d*x)^3, x) == :(-40//9*a^-1*d^-4*f^3*cosh(c+d*x)+-2//27*a^-1*d^-4*f^3*cosh(c+d*x)^3+-1//4*b*a^-2*d^-1*(e+f*x)^3+2//3*a^-1*d^-1*(e+f*x)^3*sinh(c+d*x)+a^-3*b^2*d^-1*(e+f*x)^3*sinh(c+d*x)+-6*a^-3*b^2*d^-4*f^3*cosh(c+d*x)+-2*f*a^-1*d^-2*(e+f*x)^2*cosh(c+d*x)+-3//8*b*x*a^-2*d^-3*f^3+-1//2*b*a^-2*d^-1*(e+f*x)^3*sinh(c+d*x)^2+-1//3*f*a^-1*d^-2*(e+f*x)^2*cosh(c+d*x)^3+1//3*a^-1*d^-1*(e+f*x)^3*cosh(c+d*x)^2*sinh(c+d*x)+1//4*b*a^-4*f^-1*(e+f*x)^4*(a^2+b^2)+40//9*a^-1*d^-3*f^2*(e+f*x)*sinh(c+d*x)+-1*b*a^-4*d^-1*(e+f*x)^3*(a^2+b^2)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-4*d^-1*(e+f*x)^3*(a^2+b^2)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-4*d^-4*f^3*(a^2+b^2)*Polylogarithms.polylog(4,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-6*b*a^-4*d^-4*f^3*(a^2+b^2)*Polylogarithms.polylog(4,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-3*f*a^-3*b^2*d^-2*(e+f*x)^2*cosh(c+d*x)+6*a^-3*b^2*d^-3*f^2*(e+f*x)*sinh(c+d*x)+-3//4*b*a^-2*d^-3*f^2*sinh(c+d*x)^2*(e+f*x)+2//9*a^-1*d^-3*f^2*cosh(c+d*x)^2*(e+f*x)*sinh(c+d*x)+3//8*b*a^-2*d^-4*f^3*cosh(c+d*x)*sinh(c+d*x)+-3*b*f*a^-4*d^-2*(e+f*x)^2*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-3*b*f*a^-4*d^-2*(e+f*x)^2*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+6*b*a^-4*d^-3*f^2*(e+f*x)*(a^2+b^2)*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+6*b*a^-4*d^-3*f^2*(e+f*x)*(a^2+b^2)*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+3//4*b*f*a^-2*d^-2*(e+f*x)^2*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*csch(c+d*x))^-1*(e+f*x)^2*cosh(c+d*x)^3, x) == :(2//3*a^-1*d^-1*(e+f*x)^2*sinh(c+d*x)+2//27*a^-1*d^-3*f^2*sinh(c+d*x)^3+14//9*a^-1*d^-3*f^2*sinh(c+d*x)+a^-3*b^2*d^-1*(e+f*x)^2*sinh(c+d*x)+2*a^-3*b^2*d^-3*f^2*sinh(c+d*x)+-4//3*f*a^-1*d^-2*(e+f*x)*cosh(c+d*x)+-2//9*f*a^-1*d^-2*cosh(c+d*x)^3*(e+f*x)+-1//2*b*a^-2*d^-1*(e+f*x)^2*sinh(c+d*x)^2+-1//4*b*a^-2*d^-1*f^2*x^2+-1//4*b*a^-2*d^-3*f^2*sinh(c+d*x)^2+1//3*b*a^-4*f^-1*(e+f*x)^3*(a^2+b^2)+1//3*a^-1*d^-1*(e+f*x)^2*cosh(c+d*x)^2*sinh(c+d*x)+-1*b*a^-4*d^-1*(e+f*x)^2*(a^2+b^2)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-4*d^-1*(e+f*x)^2*(a^2+b^2)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*f*a^-3*b^2*d^-2*(e+f*x)*cosh(c+d*x)+2*b*a^-4*d^-3*f^2*(a^2+b^2)*Polylogarithms.polylog(3,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+2*b*a^-4*d^-3*f^2*(a^2+b^2)*Polylogarithms.polylog(3,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1//2*b*e*f*x*a^-2*d^-1+(1/2)*b*f*a^-2*d^-2*(e+f*x)*cosh(c+d*x)*sinh(c+d*x)+-2*b*f*a^-4*d^-2*(e+f*x)*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-2*b*f*a^-4*d^-2*(e+f*x)*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1*cosh(c+d*x)^3*(e+f*x), x) == :(-2//3*f*a^-1*d^-2*cosh(c+d*x)+-1//9*f*a^-1*d^-2*cosh(c+d*x)^3+1//3*a^-1*d^-1*(2*e+2*f*x)*sinh(c+d*x)+a^-3*b^2*d^-1*(e+f*x)*sinh(c+d*x)+(1/2)*b*a^-4*f^-1*(e+f*x)^2*(a^2+b^2)+-1*f*a^-3*b^2*d^-2*cosh(c+d*x)+-1//2*b*a^-2*d^-1*sinh(c+d*x)^2*(e+f*x)+-1//4*b*f*x*a^-2*d^-1+1//3*a^-1*d^-1*cosh(c+d*x)^2*(e+f*x)*sinh(c+d*x)+-1*b*f*a^-4*d^-2*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*f*a^-4*d^-2*(a^2+b^2)*Polylogarithms.polylog(2,-1*a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-4*d^-1*(e+f*x)*(a^2+b^2)*log(1+a*(b+(a^2+b^2)^(1/2))^-1*exp(c+d*x))+-1*b*a^-4*d^-1*(e+f*x)*(a^2+b^2)*log(1+a*(b+-1*(a^2+b^2)^(1/2))^-1*exp(c+d*x))+1//4*b*f*a^-2*d^-2*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*csch(c+d*x))^-1*cosh(c+d*x)^3, x) == :(1//3*a^-1*d^-1*sinh(c+d*x)^3+a^-3*d^-1*(a^2+b^2)*sinh(c+d*x)+-1//2*b*a^-2*d^-1*sinh(c+d*x)^2+-1*b*a^-4*d^-1*(a^2+b^2)*log(b+a*sinh(c+d*x)))

max_stars_repo_path: test/6 Hyperbolic functions/6.6 Hyperbolic cosecant/6.6.1 (c+d x)^m (a+b csch)^n.jl


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Shell_Heat_Transfer
	Shell_Heat_Transfer()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Reynolds Number",
				:Default=>100,
				:Lower=>1E-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Prandtl Number",
				:Default=>0.7,
				:Lower=>1e-6
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Shell Side Film Coefficient",
				:Default=>1,
				:Lower=>1e-12,
				:Upper=>1e6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Phi Correction",
				:Default=>1,
				:Lower=>1e-3
			)),
			[:Re,:PR,:hshell,:Phi,]
		)
	end
	Re::positive
	PR::positive
	hshell::heat_trans_coeff
	Phi::positive
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Shell_Heat_Transfer
function atributes(in::Shell_Heat_Transfer,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Shell_Heat_Transfer(_::Dict{Symbol,Any})=begin
	newModel=Shell_Heat_Transfer()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
type Job
    cmd :: AbstractString
    jobid :: Int
    pbs_id :: Int
end

function ex(job)
    output = "$(job.pbs_id)-$(myid()-1).out"
    err = "$(job.pbs_id)-$(myid()-1).err"
    println("start $(job.jobid)-th job on process $(myid()-1).")
    open(io->println(io, "start $(job.jobid)-th job on process $(myid()-1)."), output, "a")
    open(io->println(io, "start $(job.jobid)-th job on process $(myid()-1)."), err, "a")
    tic()
    run(pipeline(`sh -c $(job.cmd)`, stdout = output, stderr = err, append=true))
    t = toq()
    println("finish $(job.jobid)-th job on process $(myid()-1) ($t sec)")
    open(io->println(io, "finish $(job.jobid)-th job on process $(myid()-1) ($t sec)"), output, "a")
    open(io->println(io, "finish $(job.jobid)-th job on process $(myid()-1) ($t sec)"), err, "a")
    return nothing
end

max_stars_repo_name: yomichi/SekireiTools

max_issues_repo_name: yomichi/SekireiTools

max_forks_repo_name: yomichi/SekireiTools


text:
module PlotAxes
using AxisArrays
using DataFrames
using Requires
using Dates
using Statistics

export asplotable, plotaxes, linrange, logrange, rangefn

struct ContinuousPlotAxis{N}
  step::N
  scale::Symbol
end

struct QualitativePlotAxis end

PlotAxis(x::Vector{<:Number}) = ContinuousPlotAxis(cleanup(x[2] - x[1]),:linear)
PlotAxis(x::AbstractRange{<:Number}) = ContinuousPlotAxis(cleanup(step(x)),:linear)
PlotAxis(x::Vector{<:TimeType}) = ContinuousPlotAxis(cleanup(x[2] - x[1]),:linear)
PlotAxis(x::AbstractRange{<:TimeType}) = ContinuousPlotAxis(cleanup(step(x)),:linear)
PlotAxis(x) = QualitativePlotAxis()

const current_backend = Ref{Union{Nothing,Symbol}}(nothing)
const available_backends = Dict{Symbol,Function}()

"""
    plotaxes(data,[axis1,axis2,etc...];quantize=(100,100,10,10,...))

A rudimentary, quantized display of large arrays of medium dimensionality (up
to about 4-6 dimensions, depending on the backend).

## Plot layout

By default all axes are plotted, but you can explicitly specifiy the names of
the axes as symbols to look at the data averaged across the unlisted
dimensions.

A single axis is plotted as a line. Multiple axes are plotted as a heatmap.
The first two axes specified are the x and y axes of this heatmap. The
remaining axes are plotted along rows and columns of a grid of plots; some
backends allow a row or column to represent multiple dimensions by wrapping a
dimension (e.g. ggplot).

## Display backends

You can change how the plot is displayed using `PlotAxes.set_backend!`.

## Axis Transformations

You can transform axes by a given function by explicitly specifying
axes, and using a `name => function`. Example:

    using PlotAxes, VegaLite
    data = AxisArray(rand(10,10),Axis{:a}(range(0,1,length=10)),
        Axis{:b}(exp.(range(0,1,length=10))))
    df, = PlotAxes.asplotable(data,:a,:b => logrange)

Will result in a plot with a y axis named `logb` with axis values ranging
between 0 and 1.

## Data Quantization

The data are quantized by default to maintain reasonable performance. You can
change the amount of quantization, specifying the maximum number of bins per
axis as a tuple. The order of a quantization tuple is the same as the
axis arguments passed, which defaults to the natural order of the dimensinos
(rows, cols, etc...).

"""
function plotaxes(args...;kwds...)
  if current_backend[] isa Nothing
    error("No backend defined for plot axes. Call `PlotAxes.set_backend`")
  else
    fn = available_backends[current_backend[]]
    fn(args...;kwds...)
  end
end

"""
    set_backend!(symbol)

Set the backend used to display plots when calling `plotaxes`. Call
`list_backends()` for a list of available backends.

Note that, when a package with a backend is loaded (e.g. `using Gadfly`) this
method will be called automatically for the new backend.
"""
function set_backend!(x::Symbol)
  if x ∉ list_backends()
    error("The symbol `$x` is not an available backend. "*
      "Select from one of $(list_backends()).")
  else
    current_backend[] = x
  end
end

"""
    list_backends()

List all currently available backends for plotting with `plotaxes`. This will
be populated with available backends as packages that are supported by
`PlotAxes` are loaded (e.g. via `using`)

# Supported backends

- Gadfly
- VegaLite
- RCall (via ggplot2)

"""
list_backends() = keys(available_backends)

asplotable(x::AbstractArray,args...;kwds...) =
  asplotable(AxisArray(x),args...;kwds...)
asplotable(x::AxisArray;kwds...) = asplotable(x,axisnames(x)...;kwds...)
default_quantize(x) = (100,)
default_quantize(x,y) = (100,100,)
default_quantize(x,y,args...) = (100,100,fill(10,length(args))...)
bin(i,step) = floor(Int,(i-1)/step)+1
bin(ii::CartesianIndex,steps) = CartesianIndex(bin.(ii.I,steps))
# unbin(i,step) = (i-1)*step + 1, i*step

cleanup(x::Number) = x
cleanup(x::TimeType) = x
cleanup(x) = string(x)
default_value(::Type{T}) where T <: Number = zero(float(T))
default_value(::Type{T}) where T <: TimeType = T(0)
default_value(::Type) = ""

function quantize(x::AbstractRange{<:TimeType},steps::Number)
  step = steps[1]
  qsize = bin(length(x),step)
  if qsize >= length(x)
    return x
  end

  range(first(x),last(x),step=(last(x) - first(x))/(qsize-1))
end

function quantize(x,steps)
  qsize = bin.(size(x),steps)
  if all(qsize .>= size(x))
    return x
  end
  if default_value(eltype(x)) isa String
    error("Cannot quantize non-numeric value of type $(eltype(x)).")
  end

  values = fill(default_value(eltype(x)),qsize)
  # TODO: computation of n could be optimized
  # we're taking a "dumb" approach that is easy to understand but inefficient
  n = fill(0,qsize)

  for I in CartesianIndices(x)
    values[bin(I,steps)] += cleanup(x[I])
    n[bin(I,steps)] += 1
  end
  values ./= n
  values
end

axarg_name(x::Symbol) = x
axarg_name((name,fn)::Pair) = name
axarg_name(x) =
  throw(ArgumentError("Unexpected argument. Must be a Symbol or Pair."))

struct RangeFn{F}
  fn::F
  tol::Float64
end
(r::RangeFn)(x) = r.fn(x)
(r::RangeFn)(;tol) = RangeFn(r.fn,tol)

rawnum(x) = x
rawnum(x::TimePeriod) = Dates.value(x)
torange(xs) = range(first(xs),last(xs), length=length(xs))
torange(xs::AbstractArray{<:TimeType}) =
  range(first(xs),last(xs), step=(last(xs) - first(xs))/(length(xs)-1))
rangeable(::Type{<:Number}) = true
rangeable(::Type{<:TimeType}) = true
rangeable(x) = false

function Base.Broadcast.broadcasted(r::RangeFn{typeof(identity)},
  vals::AbstractRange)

  vals
end

function Base.Broadcast.broadcasted(r::RangeFn{typeof(identity)}, vals)
  if rangeable(eltype(vals))
    fnvals = cleanup.(vals)
    std(rawnum.(diff(fnvals))) < r.tol ? torange(fnvals) : fnvals
  else
    vals
  end
end

function Base.Broadcast.broadcasted(r::RangeFn,vals)
  fnvals = r.fn.(cleanup.(vals))
  if rangeable(eltype(fnvals))
    std(rawnum.(diff(fnvals))) < r.tol ? torange(fnvals) : fnvals
  else
    fnvals
  end
end

"""
    rangefn(fn;tol=1e-8)

Transforms `fn` to a ranged function. A ranged function acts just like the
wrapped, single-arugment function except that when it is broadcast over a
number of values and those values are very close to evenly spaced, the
returned values are returned as a Range of exactly evenly spaced values.

# Example

    julia> fn = rangefn(identity,tol=0.1)
    julia> fn.([1,2,3,4.1,5,6,7,8])
    1.0:1.0:8.0

# Tolerance

What counts as "very close" can be adjusted using `tol`. If `std(diff(x))` is
smaller than tol, a `Range` is returned. The tolerance can also be modified
by passing a single keyword argument to the ranged function. That is, the
following is true:

    f = rangefn(fn)
    f(tol=1e-12) == rangefn(fn,tol=1e-12)

"""
rangefn(x;tol=1e-8) = RangeFn(x,tol)

"""
    logrange(x;tol=1e-8)

Log transform of an axis. This is similar to `log`, but will translate
broadcasted output to a `Range` object when std(diff(xs)) < tol.

# See also
- `linrange`
- `rangefn`
"""
const logrange = rangefn(log)

"""
    linrange(x;tol=1e-8)

Linear transform of an axis. This is similar to `identity`, but will translate
broadcasted output to a `Range` object when std(diff(xs)) < tol.

# See also
- `logrange`
- `rangefn`
"""
const linrange = rangefn(identity)

const exact = identity

fn_prefix(x) = ""
fn_prefix(x::typeof(identity)) = "exact_"
fn_prefix(x::RangeFn{typeof(identity)}) = ""
fn_prefix(x::RangeFn) = fn_prefix(x.fn)
function fn_prefix(x::Function)
  name = string(Base.typename(typeof(x)))
  pattern =
    r"typeof\([-[:alpha:]_\u2207][[:word:]\u207A-\u209C!\u2032\u2207]*\)"
  if occursin(pattern,name)
    replace(name,r"typeof\((.*)\)" => s"\1")
  else
    ""
  end
end

# apply any specified transform to the axes
function transform_axes(xs,showax)
  allax = collect(Any,axisnames(xs))
  showax = collect(showax)
  showax_i = indexin(axarg_name.(showax),allax)
  if any(isnothing,showax_i)
    ax = axarg_name(showax[findfirst(isnothing,showax_i)])
    error("Could not find the axis $ax.")
  end
  allax[showax_i] = showax
  axs = map(ax -> transform_axis(xs,ax),allax)

  result = AxisArray(xs.data,axs...)
  result, axisnames(result)[showax_i]
end

# default to a linear "transformation"
transform_axis(xs,name::Symbol) = transform_axis(xs,name => linrange)
function transform_axis(xs,(name,fn)::Pair)
  allvals = axisvalues(xs)
  axisnames(xs)
  n = axisdim(xs,Axis{name})
  vals = fn.(allvals[n])
  newname = Symbol(string(fn_prefix(fn),name))
  Axis{newname}(vals)
end

function asplotable(x::AxisArray,ax1,axes...;
                    quantize=default_quantize(ax1,axes...))
  show_axes = (ax1,axes...)
  x,show_axes = transform_axes(x,show_axes)
  qs = map(axisnames(x)) do ax
    if ax ∈ show_axes
      min(size(x,Axis{ax}),quantize[findfirst(isequal(ax),show_axes)])
    else
      1
    end
  end

  steps = size(x) ./ qs
  vals = PlotAxes.quantize(x,steps)
  axqvals = PlotAxes.quantize.(axisvalues(x),steps)

  df = DataFrame(value = convert(Array{eltype(vals)},vec(vals)))
  for ax in show_axes
    axi = findfirst(isequal(axarg_name(ax)),axisnames(x))
    df[!,ax] .= default_value(eltype(axqvals[axi]))
    for (j,jj) in enumerate(CartesianIndices(vals))
      df[j,ax] = cleanup(axqvals[axi][jj.I[axi]])
    end
  end

  df, map(axv -> PlotAxis(axv),
    map(ax -> axqvals[findfirst(isequal(ax),axisnames(x))],show_axes))
end

# using Gadfly
# include("gadfly.jl")

const AxisId = Union{Symbol,Pair}

function __init__()
  @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" begin
    using .RCall
    @info "Loading RCall ggplot2 backend for `PlotAxes`"
    include("ggplot2.jl")
  end
  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    using .Unitful
    cleanup(x::Quantity) = ustrip(x)
    default_value(::Type{<:Quantity{T}}) where T = default_value(T)
  end
  @require VegaLite="112f6efa-9a02-5b7d-90c0-432ed331239a" begin
    using .VegaLite
    @info "Loading VegaLite backend for `PlotAxes`"
    include("vegalite.jl")
  end
  @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" begin
    using .Gadfly
    @info "Loading Gadfly backend for `PlotAxes`"
    include("gadfly.jl")
  end
end

end # module


text: ENV["GITHUB_ACCESS_TOKEN"] = "YOUR GITHUB TOKEN HERE"

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
export MNIST

"""
The MNIST database of handwritten digits
- Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
- Website: http://yann.lecun.com/exdb/mnist/
MNIST is a classic image-classification dataset that is often
used in small-scale machine learning experiments. It contains
70,000 images of handwritten digits. Each observation is a 28x28
pixel gray-scale image that depicts a handwritten version of 1 of
the 10 possible digits (0-9).
# Interface
- [`MNIST.traintensor`](@ref), [`MNIST.trainlabels`](@ref), [`MNIST.traindata`](@ref)
- [`MNIST.testtensor`](@ref), [`MNIST.testlabels`](@ref), [`MNIST.testdata`](@ref)
# Utilities
- [`MNIST.download`](@ref)
- [`MNIST.convert2image`](@ref)
"""
module MNIST
    using DataDeps
    using ColorTypes
    using FixedPointNumbers

    bytes_to_type(::Type{UInt8}, A::Array{UInt8}) = A
    bytes_to_type(::Type{N0f8}, A::Array{UInt8}) = reinterpret(N0f8, A)
    bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:Integer = convert(Array{T}, A)
    bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:AbstractFloat = A ./ T(255)
    bytes_to_type(::Type{T}, A::Array{UInt8}) where T<:Number  = convert(Array{T}, reinterpret(N0f8, A))

    function with_accept(f, manual_overwrite)
        auto_accept = if manual_overwrite == nothing
            get(ENV, "DATADEPS_ALWAYS_ACCEPT", false)
        else
            manual_overwrite
        end
        withenv(f, "DATADEPS_ALWAYS_ACCEPT" => string(auto_accept))
    end

    function datadir(depname, dir = nothing; i_accept_the_terms_of_use = nothing)
        with_accept(i_accept_the_terms_of_use) do
            if dir == nothing
                # use DataDeps defaults
                @datadep_str depname
            else
                # use user-provided dir
                if isdir(dir)
                    dir
                else
                    DataDeps.env_bool("DATADEPS_DISABLE_DOWNLOAD") && error("DATADEPS_DISABLE_DOWNLOAD enviroment variable set. Can not trigger download.")
                    DataDeps.download(DataDeps.registry[depname], dir)
                    dir
                end
            end
        end::String
    end

    function datafile(depname, filename, dir = nothing; recurse = true, kw...)
        path = joinpath(datadir(depname, dir; kw...), filename)
        if !isfile(path)
            @warn "The file \"$path\" does not exist, even though the dataset-specific folder does. This is an unusual situation that may have been caused by a manual creation of an empty folder, or manual deletion of the given file \"$filename\"."
            if dir === nothing
                @info "Retriggering DataDeps.jl for \"$depname\""
                download_dep(depname; kw...)
            else
                @info "Retriggering DataDeps.jl for \"$depname\" to \"$dir\"."
                download_dep(depname, dir; kw...)
            end
            if recurse
                datafile(depname, filename, dir; recurse = false, kw...)
            else
                error("The file \"$path\" still does not exist. One possible explaination could be a spelling error in the name of the requested file.")
            end
        else
            path
        end::String
    end

    function download_dep(depname, dir = DataDeps.determine_save_path(depname); kw...)
        DataDeps.download(DataDeps.registry[depname], dir; kw...)
    end

    function download_docstring(modname, depname)
        """
        The corresponding resource file(s) of the dataset is/are
        expected to be located in the specified directory `dir`. If
        `dir` is omitted the directories in
        `DataDeps.default_loadpath` will be searched for an existing
        `$(depname)` subfolder. In case no such subfolder is found,
        `dir` will default to `~/.julia/datadeps/$(depname)`. In the
        case that `dir` does not yet exist, a download prompt will be
        triggered. You can also use `$(modname).download([dir])`
        explicitly for pre-downloading (or re-downloading) the
        dataset. Please take a look at the documentation of the
        package DataDeps.jl for more detail and configuration
        options.
        """
    end

    function _colorview(::Type{T}, array::AbstractArray{<:Number}) where T <: Color
        __images_supported__ ||
            error("Converting to image requires `using ImageCore`")
        ImageCore.colorview(T, array)
    end

    export

        traintensor,
        testtensor,

        trainlabels,
        testlabels,

        traindata,
        testdata,

        convert2image,

        download

    @deprecate convert2features reshape

    const DEPNAME = "MNIST"
    const TRAINIMAGES = "train-images-idx3-ubyte.gz"
    const TRAINLABELS = "train-labels-idx1-ubyte.gz"
    const TESTIMAGES  = "t10k-images-idx3-ubyte.gz"
    const TESTLABELS  = "t10k-labels-idx1-ubyte.gz"

    """
        download([dir]; [i_accept_the_terms_of_use])
    Trigger the (interactive) download of the full dataset into
    "`dir`". If no `dir` is provided the dataset will be
    downloaded into "~/.julia/datadeps/$DEPNAME".
    This function will display an interactive dialog unless
    either the keyword parameter `i_accept_the_terms_of_use` or
    the environment variable `DATADEPS_ALWAYS_ACCEPT` is set to
    `true`. Note that using the data responsibly and respecting
    copyright/terms-of-use remains your responsibility.
    """
    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    include(joinpath("Reader","Reader.jl"))
    include("interface.jl")
    include("utils.jl")

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: THE MNIST DATABASE of handwritten digits
            Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
            Website: http://yann.lecun.com/exdb/mnist/
            [LeCun et al., 1998a]
                Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
                "Gradient-based learning applied to document recognition."
                Proceedings of the IEEE, 86(11):2278-2324, November 1998
            The files are available for download at the offical
            website linked above. Note that using the data
            responsibly and respecting copyright remains your
            responsibility. The authors of MNIST aren't really
            explicit about any terms of use, so please read the
            website to make sure you want to download the
            dataset.
            """,
            "https://ossci-datasets.s3.amazonaws.com/mnist/" .* [TRAINIMAGES, TRAINLABELS, TESTIMAGES, TESTLABELS],
            "0bb1d5775d852fc5bb32c76ca15a7eb4e9a3b1514a2493f7edfcf49b639d7975",
        ))
    end
end


text:
"Generic state difference represented as additions, deletions, and assignments."
struct GenericDiff <: Diff
    add::Vector{Term} # List of additions
    del::Vector{Term} # List of deletions
    ops::Dict{Term,Term} # Dictionary of assignment operations
end

GenericDiff() = GenericDiff(Term[], Term[], Dict{Term,Any}())

function combine!(d1::GenericDiff, d2::GenericDiff)
    append!(d1.add, d2.add)
    append!(d1.del, d2.del)
    merge!(d1.ops, d2.ops)
    return d1
end

function as_term(diff::GenericDiff)
    add = diff.add
    del = [Compound(:not, [t]) for t in diff.del]
    ops = [Compound(:assign, [t, v]) for (t, v) in diff.ops]
    return Compound(:and, [add; del; ops])
end

is_redundant(diff::GenericDiff) =
    issetequal(diff.add, diff.del) && all(k == v for (k, v) in diff.ops)

Base.empty(diff::GenericDiff) = GenericDiff()

"Conditional state difference, represented as paired conditions and sub-diffs."
struct ConditionalDiff{D <: Diff} <: Diff
    conds::Vector{Vector{Term}}
    diffs::Vector{D}
end

ConditionalDiff{D}() where {D} = ConditionalDiff{D}([], Vector{D}())

function combine!(d1::ConditionalDiff{D}, d2::D) where {D}
    for (i, cs) in enumerate(d1.conds)
        if isempty(cs)
            combine!(d1.diffs[i], d2)
            return d1
        end
    end
    push!(d1.conds, Term[])
    push!(d1.diffs, d2)
    return d1
end

function combine!(d1::ConditionalDiff{D}, d2::ConditionalDiff{D}) where {D}
    append!(d1.conds, d2.conds)
    append!(d1.diffs, d2.diffs)
    return d1
end

function as_term(diff::ConditionalDiff)
    branches = map(zip(diff.conds, diff.diffs)) do (cs, d)
        effect = as_term(d)
        isempty(cs) && return effect
        cond = length(cs) == 1 ? cs[1] : Compound(:and, cs)
        return Compound(:when, [cond, effect])
    end
    return Compound(:and, branches)
end

is_redundant(diff::ConditionalDiff) = all(is_redundant.(diff.diffs))

Base.empty(diff::ConditionalDiff{D}) where {D} = ConditionalDiff{D}()


text:
#=
Coherent exciton scattering (CES) approximation (Roden, Eisfeld, Briggs)

* mean-field type approximation
* xact monomer Green function (electronic and vibrational DOF) is replaced by its average in the vibrational ground state
* reproduces absorption and circ. dichroism of H- and J-aggregates

=#

using PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles, MultidimensionalSpectroscopy

# make sure to set script directory as pwd()
cd(@__DIR__)

# to use PyPlot in GUI mode under Linux
pygui(true)
# run once with calc_2d = false to initialize functions
calc_2d = false   # !! use less values for k if running 2D

Ω(k)      = Emon + 2 * V * cos( pi * k / (N + 1) )
Ω2(k1,k2) = Emon + 2 * V * cos( pi * k1 / (N + 1) ) + Emon + 2 * V * cos( pi * k2 / (N + 1) )
δ(l,m) = l == m
# sum over 2 indices that are not equal
[i+j for i in 1:2 for j in 1:2 if i != j]
# https://aip.scitation.org/doi/pdf/10.1063/1.469393
μμ(k)        = μ_mon * sqrt( 2 / (N+1)) *             (1 - (-1)^k ) / 2 * cot( pi * k  / (2 * (N+1)))
μμ2(k1,k2,k) = μ_mon * sqrt( 2 / (N+1)) * ( δ(k2,k) * (1 + (-1)^k1) / 2 * cot( pi * k1 / (2 * (N+1))) -
                                            δ(k1,k) * (1 + (-1)^k2) / 2 * cot( pi * k2 / (2 * (N+1))) +
                                            1/2 * ( δ(k1-k2-k,0) - δ(k1-k2+k,0) ) *
                                                  ( cot(pi * k1 / (2 * (N+1))) + cot(pi * k2 / (2 * (N+1))) ) +
                                            1/2 * ( δ(k1+k2+k,2*(N+1)) - δ(k1+k2-k,0) ) *
                                                  ( cot(pi * k1 / (2 * (N+1))) - cot(pi * k2 / (2 * (N+1))) ) )

Sp(i) = transition(B_TDBC,1+i,i)
Sm(i) = transition(B_TDBC,i,1+i)

k  = [1, 2, 3, 5, 7, 9]
k  = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
k1 = [1, 2]#, 3, 4, 5, 6, 7]
k2 = [1, 2]#, 3, 4, 5, 6, 7]
K  = [(2,1),(3,2), (4,1), (6,1), (5,2)]

V = -2455 / 2 / 8065
μ_mon = 1
Emon = 2.41532 # 513 nm ... https://aip.scitation.org/doi/pdf/10.1063/1.469393
Emon = 2.41032 # 513 nm ... https://aip.scitation.org/doi/pdf/10.1063/1.469393

B_TDBC = NLevelBasis(1 + length(k) + length(K))

try     # need to delete method first when working on it
    m = @which make_Ops(1)
    Base.delete_method(m)
catch
end
function make_Ops(num)
    global N = num
    ##
    ωs_ge  = [Ω(i) for i in k]
    #ωs_ef  = [Ω2(i,j) for i in k for j in k2 if i > j]
    ωs_ef  = [Ω2(i...) for i in K]
    ωs     = vcat(ωs_ge,ωs_ef)

    μs_ge = [μμ(i) for i in k]
    #μs_ge = sqrt.(complex(μs_ge))
    μs_ef = zeros(length(k),length(K))
    for l in 1:length(K)
        for i in 1:length(k)
            μs_ef[i,l] = μμ2(K[l]...,k[i])
        end
    end

    H   = sum(ωs[i] * transition(B_TDBC,1+i,1) * transition(B_TDBC,1,1+i) for i in 1:length(ωs))
    μ12 = sum(μs_ge[i] * transition(B_TDBC,1+i,1) for i in 1:length(k))
    μ23 = sum(μs_ef[i,j] * transition(B_TDBC,1+length(k)+j,1+i) for j in 1:length(K) for i in 1:length(k))

    ##

    Σ⁺ge = μ12
    Σ⁺ef = μ23
    Σ⁺ = μ12 + μ23
    #Σ⁺ = sum(transition(testb,1+i,1) for i in 1:length(k))

    #Σ⁺ = sum(Sp(i) for i in k)
    Σ⁻ = dagger(Σ⁺);
    μ12 = (μ12) + dagger(μ12);
    #μ12.data = sqrt.(μ12.data)
    μ23 = μ23 + dagger(μ23)
    #μ23.data = sqrt.(μ23.data)

    # make initial density matrix
    Psi0 = nlevelstate(B_TDBC,1)
    rho0 = dm(Psi0)

    # make collapse operator
    #ccc = diagonaloperator(b_TLS,[1, 1/4, 1])

    #     relaxation        dephasing
    #    _____________    _______________
    #L = [j12, j13, j23,  j21*j12, j31*j13]
    #Γ = [.01,.00 ,.05,    .18,     .01]
    #L = Γ .* L
    #...redo for new system
    #L = [.05 * cmds.tri(μ12,"L") * cmds.tri(μ12,"U")]
    L1 = 0.01 .* [transition(B_TDBC,1,1+i) for i in 1:length(k)]
    L2 = 0.09 .* [transition(B_TDBC,2,2+i) for i in 1:length(k)]

    L_dephasing     = [.04 * Sp(i) * Sm(i) for i in 1:length(k)]
    L_dephasingESA  = [.06 * Sp(i) * Sm(i) for i in length(k):length(k)+length(K)]

    L = vcat(L1,L2,L_dephasing)
    L = vcat(L1,L2)

    ##

    rho1 = μ12 * rho0 * μ12
    μ12  = μ12 / sqrt(tr(rho1))
    μ23  = μ23 / sqrt(tr(rho1))
    rho1 = μ12 * rho0 * μ12
    rho2 = μ23 * rho1 * μ23
    #μ23  = μ23 / sqrt(tr(rho2))
    rho2 = μ23 * rho1 * μ23

    return H, μ12, μ23, Σ⁺ge, Σ⁺ef, Σ⁺, Σ⁻, μs_ge, μs_ef, ωs_ge, ωs_ef, rho0, rho1, rho2, L

end

H, μ12, μ23, Σ⁺ge, Σ⁺ef, Σ⁺, Σ⁻, μs_ge, μs_ef, ωs_ge, ωs_ef, rho0, rho1, rho2, L = make_Ops(14)

##
#figure()
#bar(ωs_ge,μs_ge.^2,width=0.003)
##

cmp = create_colormap("bright");

Δ = 0.05

C = 2416 / 8065.5
# https://www.pks.mpg.de/fileadmin/user_upload/MPIPKS/group_pages/QuantumAggregates/PDF/pubilcations/Eisfeld_Briggs_2002_CP_61.pdf
# https://www.pks.mpg.de/~eisfeld/wwww/pdf/Eisfeld_Briggs_2006_CP_376.pdf

Emon = 2.41532 # 513 nm ... https://aip.scitation.org/doi/pdf/10.1063/1.469393
E = 2.08
E = Emon - C # Eisfeld
E2 = 2*E + Δ

# or analytically https://aip.scitation.org/doi/pdf/10.1063/1.469393

#V = -2505 / 2
#Eg_k1 = Emon + 2*V * cos(pi/(20+1)) / 8065
#E = Eg_k1

#Ek1_k1k2 = Emon + 2*V * cos(2*pi/(20+1)) / 8065

#E2 = Eg_k1 + Ek1_k1k2

# get eigenstates and eigenvalues of H
#  using Julia LinearAlgebra functions instead of QuantumOptics requires to use H.data
#states   = eigvecs(dense(H).data)
#energies = eigvals(dense(H).data)

# set up figure
fig  = figure(figsize=(15,9));
ax1  = subplot(241); title("Energy levels")
ax1a = subplot(242); title("Noise-power spectrum S(ω)")
ax2  = subplot(243); title("correlation function LB")
ax2a = subplot(244); title("correlation function RF")
ax3  = subplot(223); title("Kinetics in eigen/coupled basis")
ax4  = subplot(247); title("Abs. spectrum LB")
ax4a = subplot(248); title("Abs. spectrum RF")

sca(ax1)
plot_levels(H,0)
hlines(.5*E2,-.45,.45, color="red",linestyle="dashed")
hlines(Emon,-.45,.45, color="green",linestyle="dashed")
hlines(2*Emon,-.45,.45, color="green",linestyle="dashed")

# time list for evaluation of master equation
Tmax = 400
dt = 1.
tlist = [0:dt:Tmax;]

#M = 1
#corr = sum(timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12) for i in 1:M) ./ M
##
D = .1
tc = .9
#gt = D^2 * tc^2 * (exp.(-collect(-150:0.2:150)./tc) .+ collect(-150:0.2:150)./tc .- 1)
gt = D^2 * tc^2 * (exp.(-tlist./tc) .+ tlist./tc .- 1)

Δ = .01 # root-mean-square amplitude of ⟨δω(t)|δω(0)⟩
g_slow = Δ^2 .* tlist.^2 ./2 # static limit
Λ = 100 # inverse correlation time (tc)⁻¹
g_fast = Δ^2 .* tlist ./ Λ   # Markovian limit

# or ...
#Δ_slow = 10 * 0.00414
#Λ_slow = 1
#Δ_fast = 54 * 0.00414
#Λ_fast = 5
#
#g_slow = Δ_slow^2 / Λ_slow^2 .* (exp.(-Λ_slow .* tlist) .+ Λ_slow .* tlist .-1)
#g_fast = Δ_fast^2 / Λ_fast^2 .* (exp.(-Λ_fast .* tlist) .+ Λ_fast .* tlist .-1)

# from MUKAMEL Chapter 8

#gt = 2 * λ * kb * T / (ħ * Λ^2) * (exp.(-Λ .* tlist) .+ Λ .* tlist .- 1) .
#        - 1im * (λ/Λ) .* (exp.(-Λ .* tlist) .+ Λ .* tlist .- 1)

κ = 3
Λ = .02

λ = 0.2 # half the Stokes shift in eV !?
kB = 8.617333262145 * 10^-5 # in eV * K^-1
T = 300 # Kelvin
ħ = 6.582119569 * 10^-16 # in eV s
#Δ = sqrt(2 * λ * kB * T / ħ)

Δ = Λ / κ # unit 1/s

gt = Δ^2 / Λ^2 * (exp.(-Λ .* tlist) .+ Λ .* tlist .- 1)

#gt = g_slow + g_fast

corr = zeros(size(tlist))
list = [12, 13, 14, 15, 16, 17, 18]
list = collect(14:.5:16)
for s in list
    local H, μ12, μ23, Σ⁺ge, Σ⁺ef, Σ⁺, Σ⁻, μs_ge, μs_ef, ωs_ge, ωs_ef, rho0, rho1, rho2, L = make_Ops(s);
    global corr = exp(-(s-15)^2) .* (timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12) + corr) .* exp.(-gt) ./ 2
end

corr = corr ./ maximum(real(corr))

sca(ax2)
##
corr = corr .* exp.(-gt)

tout, rhot = timeevolution.master(tlist,μ12 * rho0,H,L)

# zeropad corr and extrapolate tlist
zp = 11
corr = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)
ω,spec = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true)

#figure(figsize=(6,3));
# add the results to previous figure
sca(ax2)
plot(tnew,real(corr))
plot(tlist,exp.(-gt))

xlabel("time"); ylabel("Corr. func. ⟨μ(t) μ₀₁ ρ₀⟩")

sca(ax4)
plot(-ω,spec)
xlabel("ω"); ylabel("Absorption")

tout, rhot = timeevolution.master(tlist,μ23 * rho1,H,L)
corr = sum(timecorrelations.correlation(tlist, rho1, H, L, μ23, μ23) for i in 1:N) ./ N
corr = corr .* exp.(-gt)

corr = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)
ω,specESA = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true)
#plot(-[E, E], [0, .5])
plot(-ω,-specESA)
xlabel("ω"); ylabel("Absorption")

tout, rhot = timeevolution.master(tlist,μ12 * rho0 * μ12,H,L)

vecs   = eigvecs(dense(H).data)
eivecs = [Ket(H.basis_l,vecs[:,i]) for i in 1:length(vecs[1,:])]

sca(ax3)
# levels: 1 -> GS, 2 -> LP, 3 - N-1 -> ..., N -> UP
for i in 1:length(eivecs)
    es_n = real(expect(dm(eivecs[i]),rhot))
    plot(tout,es_n,linestyle="dashed", label=string(i)); #append!(leg,[string(i)])
    ylim((-0.1, 1))
end
#legend()

#

try     # need to delete method first when working on it
    m = @which spectral_density(1)
    Base.delete_method(m)
catch
end
function spectral_density(ω) # power spectral density, thermal bath noise spectrum
    #w = [(ω .- ωq_TDBC[i]) ./ (.5 * width) for i in 1:length(ωq_TDBC)]
    b = 1 / 0.0259 # eV k_B*T at T=300 K
    f = .5 # 0 ... 1: relative importance of non-radiation vs. dephasing ... f = 0 > no dephasing
    γ = 0.7 # from absorption linewidth γ = γ_non-rad + γ_dephasing
    η = 1
    η = f * γ / (2 * pi * 0.0259)
    ω_cut = .1
    if ω == 0
        Jω = .02
    #elseif !isless(real(ω),0)
    #else
    #    w = [(ω .- ωq_TDBC[i]) ./ (.5 * width) for i in 1:length(ωq_TDBC)]
    #    Jω = sum([ampl_TDBC[i] ./ (1 .+ w[i].^2) for i in 1:length(ωq_TDBC)])  # probably most realistic ... ???
    #    #Jω = 1
    end
    if  !isless(real(ω),0) && ω != 0
        w = [(ω .- ωq_TDBC[i]) ./ (.5 * width) for i in 1:length(ωq_TDBC)]
        Jω = sum([ampl_TDBC[i] ./ (1 .+ w[i].^2) for i in 1:length(ωq_TDBC)])
        Jω = η * ω * exp(-(ω/ω_cut)^2) + Jω
        Jω = Jω * (1/exp(b *  ω - 1) + 1)
    elseif isless(real(ω),0)
        w = [(ω .- ωq_TDBC[i]) ./ (.5 * width) for i in 1:length(ωq_TDBC)]
        Jωa = sum([ampl_TDBC[i] ./ (1 .+ w[i].^2) for i in 1:length(ωq_TDBC)])
        Jω = η * -ω * exp(-(-ω/ω_cut)^2)
        Jω = Jω * (1/exp(b * -ω - 1)) + Jωa
    end
    if ω == 0
        Jω = 0.00000
    end
    return Jω

end

width = 0.0005
# from ... http://www.rsc.org/suppdata/c8/sc/c8sc00171e/c8sc00171e1.pdf
ωq_TDBC = [40, 80, 120, 150, 185, 197] ./ 8065 #.* 0.02998 # cm-1 in Hz #./ 8065 # in eV instead
ωq_TDBC = [0]
# => [0.00496 0.00992 0.01488 0.01860  0.02294 0.02443]
ampl_TDBC = [14, 18,  25,  43,  42, 67, 60] ./ 1000
ampl_TDBC = [0]

Γ = .2
#a_ops = [Γ * (j21+j31+j32) * (j12+j13+j23), spectral_density]
dat = tri(tri(Σ⁺ * Σ⁻,"U"),"L")
a_ops = [Γ * Σ⁺ * Σ⁻, spectral_density]
a_ops = [Γ * Σ⁺ge * dagger(Σ⁺), spectral_density]#, Γ * Σ⁺ef * dagger(Σ⁺ef), spectral_density]
#a_ops = [dat, spectral_density]

R, ekets = timeevolution.bloch_redfield_tensor(H, [a_ops]; J=L)

# linear absorption
tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0,R,H)
corr = expect(μ12,rhot)
#corr = corr .* exp.(-gt)
#plot(tlist,corr .* exp.(-gt))

# calculate and plot spectrum
corr = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)
sca(ax2a)
plot(tnew,real(corr))
ω, spec = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true);
sca(ax4a)
plot(-ω,spec)

# ESA
tout, rhot = timeevolution.master_bloch_redfield(tlist,μ23*rho1,R,H)
corr = expect(μ23,rhot)
#corr = corr .* exp.(-gt)
#plot(tlist,corr .* exp.(-gt))

# calculate and plot spectrum
corr = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)
sca(ax2a)
plot(tnew,real(corr))
ω, specESA = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true);
sca(ax4a)
plot(-ω,-specESA)

#figure()
sca(ax1a)
plot(collect(-1:.001:1),(spectral_density.(collect(-1:.001:1))))

# load TA
TDBC_TA_t0 = readdlm("TDBC_TA_t0.dat",',') # https://aip.scitation.org/doi/pdf/10.1063/1.469393
TDBC_agg_abs = readdlm("TDBCaggr_abs2.dat",',')
calc_ω = ω[ω.<0]
calc_abs = spec[ω.<0]

#convert to energy axis if it is in wavelength
if TDBC_agg_abs[1,1] > 100
    TDBC_agg_abs[:,1] = 6.582119569 * 10^-16 * 3e8 ./ (TDBC_agg_abs[:,1] * 1e-9) * 2 * pi
end

sca(ax4a)
plot(TDBC_agg_abs[:,1],TDBC_agg_abs[:,2],"k--")
bar(ωs_ge,(μs_ge./maximum(μs_ge)).^2,width=0.01,color="k")
#bar(ωs_ef,-(μs_ef[1]./maximum(μs_ef[1])).^2,width=0.01,color="k")

# scale ESA tdm s by ge tdm
μs_ef = μs_ef
temp = [μs_ef[:,i] ./ maximum(μs_ef[:,i]) .* (μs_ge./maximum(μs_ge))[i] for i in 1:length(ωs_ef)]
ωtemp = [[i .- j for j in ωs_ge] for i in ωs_ef]

for i in 1:length(temp)
    #bar(ωs_ef[i] .- ωs_ef[1] .+ ωs_ge,-temp[i], width=0.01)
    bar(ωtemp[i],-temp[i], width=0.01)
end

sca(ax4)
plot(TDBC_agg_abs[:,1],TDBC_agg_abs[:,2],"k--")
bar(ωs_ge,(μs_ge./maximum(μs_ge)).^2,width=0.01,color="k")
for i in 1:length(temp)
    bar(ωtemp[i],-temp[i], width=0.01)
end

tight_layout()

method = "redfield"
if method == "redfield"
    F       = R
    use_sub = false
elseif method == "lindblad"
    F       = L
    use_sub = false
end

#calc_2d = false
if calc_2d

        ## calculate (complex) 3rd order corr function (with T=0)
        zp = 11 # zeropad up to 2^zp

        ## calculate 2D spectra at
        T = [0, 50, 100, 200, 500, 1000] #0.66 fs

        spectra2d = Array{out2d}(undef, length(T))

        N = 50
        Threads.@threads for i = 1:length(T)
        #for i = 1:length(T)
            spectra2d[i] = make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],
                                                method;debug=false,use_sub=use_sub,t2coh="kin",zp=zp);
        end

        ## crop 2D data and increase dw
        spectra2d = [crop2d(spectra2d[i],1.5;w_max=3,step=1) for i = 1:length(T)]
end

if calc_2d
        ## plot 2D spectra for each(?) T
        # what to plot
        rep="abst"
        scal="lin"

        ## make  subplot layout
        nplots = length(T);
        ncols = Int32(ceil(sqrt(nplots)));
        nrows = Int32(ceil(nplots / ncols));

        # determine maximum value in dataset out2d[:].full2d[:,:]
        maxi = maximum([maximum(real(spectra2d[i].full2d)) for i in 1:length(spectra2d)])

        # plot 2D spectra
        fig, ax = subplots(nrows,ncols,sharex=true,sharey=true,figsize=(ncols*3.5,nrows*3))
        if nplots == 1
                ax = [ax]
        end
        fig.suptitle(rep * " 2D spectrum (" * scal * ". scaling)")
        k = 0
        for i = 1:nrows
            for j = 1:ncols
                global k += 1
                if k > nplots
                    continue
                end
                sca(ax[i,j])
                ax[i,j].set_aspect="equal"
                plot2d(spectra2d[k].ω,round.(spectra2d[k].full2d,digits=1);repr=rep,scaling=scal,norm=maxi)
                title("2D spectrum at $(T[k]) fs")
            end
        end
        tight_layout()
        subplots_adjust(top=0.9)

        ω = spectra2d[1].ω
        ## plot additional things, like energy levels of states
        #plot([E, E], [ω[1], ω[end]],"k--",linewidth=1,alpha=0.25)
        #plot([ω[1], ω[end]],[E, E],"k--",linewidth=1,alpha=0.25)
        #plot([energies[2], energies[2]], [ω[1], ω[end]],"g--",linewidth=1,alpha=0.25)
        #plot([ω[1], ω[end]],[energies[2], energies[2]],"g--",linewidth=1,alpha=0.25)

        ## Save data
        #cmds.save_2d([round.(real(out2d[i].full2d),digits=1) for i = 1:length(T)],T,fn)

        ## plot TA (summed 2D spectrum)
        figure()
        ta = [sum(spectra2d[i].full2d,dims=1) for i in 1:length(spectra2d)]
        ta = vcat(ta...)'
        ta = ta ./ maximum(real(ta))
        plot(ω,ta)
        plot(ω,zeros(size(ω)),"k",linewidth = .5)
        plot(1239.8 ./ TDBC_TA_t0[:,1], -1 * TDBC_TA_t0[:,2]/maximum(-TDBC_TA_t0[:,2]) ,"k--")
        plot(TDBC_agg_abs[:,1],TDBC_agg_abs[:,2] ./ maximum(TDBC_agg_abs[:,2]),"r--")
        plot(-calc_ω,calc_abs./maximum(calc_abs[-calc_ω.>0]),"red")
        bar(ωs_ge,(μs_ge./maximum(μs_ge)).^2,width=0.005)
        #[bar(ωs_ef[i] .- ωs_ge,-μs_ef[:,i].^2 .* (μs_ge./maximum(μs_ge)).^2, width=0.005) for i in 1:length(ωs_ef)]
        for i in 1:length(temp)
            bar(ωtemp[i],-temp[i], width=0.01)
        end
        xlabel("Energy/frequency")
        ylabel("Difference absorption")
        tight_layout()

end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, CUDA, LinearAlgebra, Plots

begin
    x0 = -1.0f0
    x1 = 1.0f0
    ncell = 100
    dx = (x1 - x0) / ncell
    deg = 2
    nsp = deg + 1
    cfl = 0.1
    dt = cfl * dx
    t = 0.0f0
    a = 1.0f0
    tspan = (0.0f0, 1.0f0)
    nt = tspan[2] / dt |> floor |> Int
    bc = :period
end

ps = FRPSpace1D(x0, x1, ncell, deg)

u = zeros(Float32, ncell, nsp)
for i = 1:ncell, ppp1 = 1:deg+1
    u[i, ppp1] = exp(-20.0 * ps.xpg[i, ppp1]^2)
end
u = u |> CuArray

prob = FR.FRAdvectionProblem(u, tspan, ps, a, bc)
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

for iter = 1:nt
    step!(itg)
end

ut = CUDA.@allowscalar itg.u[:, 2] |> Array
plot(ps.xpg[:, 2], ut, label = "t=2")
plot!(ps.xpg[:, 2], exp.(-20 .* ps.xpg[:, 2] .^ 2), label = "t=0", line = :dash)


text:
export
    AbstractGibbsNetwork,
    low_energy_spectrum,
    branch_state,
    bound_solution,
    merge_branches,
    Solution,
    Memoize

struct Solution
    energies::Vector{Float64}
    states::Vector{Vector{Int}}
    probabilities::Vector{Float64}
    degeneracy::Vector{Int}
    largest_discarded_probability::Float64
end

empty_solution() = Solution([0.], [[]], [1.], [1], -Inf)

function branch_state(network, σ)
    node = node_from_index(network, length(σ) + 1)
    basis = collect(1:length(local_energy(network, node)))
    vcat.(Ref(σ), basis)
end

function branch_solution(partial_sol::Solution, network::AbstractGibbsNetwork)
    local_dim = length(local_energy(network, node_from_index(network, length(partial_sol.states[1])+1)))
    new_energies = vcat(
            [
                (en .+ update_energy(network, state))
                for (en, state) ∈ zip(partial_sol.energies, partial_sol.states)
            ]
            ...
        )
    new_states = vcat(branch_state.(Ref(network), partial_sol.states)...)
    new_probabilities = vcat(
            [
                partial_sol.probabilities[i] .+ log.(p)
                for (i,p) ∈ enumerate(conditional_probability.(Ref(network), partial_sol.states))
                    ]
            ...
        )
    degeneracies = repeat(partial_sol.degeneracy, inner=local_dim)
    ldp = partial_sol.largest_discarded_probability
    Solution(
        new_energies, new_states, new_probabilities, degeneracies, ldp
    )
end

function merge_branches(network::AbstractGibbsNetwork{S, T}) where {S, T}
    function _merge(partial_sol::Solution)
        node = node_from_index(network, length(partial_sol.states[1])+1)
        boundaries = hcat(boundary_state.(Ref(network), partial_sol.states, Ref(node))...)'
        _unique_boundaries, indices = SpinGlassNetworks.unique_dims(boundaries, 1)

        sorting_idx = sortperm(indices)
        sorted_indices = indices[sorting_idx]
        start = 1

        energies = partial_sol.energies[sorting_idx]
        states = partial_sol.states[sorting_idx]
        probs = partial_sol.probabilities[sorting_idx]
        degeneracy = partial_sol.degeneracy[sorting_idx]

        new_energies = []
        new_states = []
        new_probs = []
        new_degeneracy = []

        while start <= size(boundaries, 1)
            stop = start
            while stop + 1 <= size(boundaries, 1) && sorted_indices[start] == sorted_indices[stop+1]
                stop = stop + 1
            end

            best_idx = argmin(@view energies[start:stop]) + start-1

            push!(new_energies, energies[best_idx])
            push!(new_states, states[best_idx])
            push!(new_probs, probs[best_idx])
            push!(new_degeneracy, degeneracy[best_idx])

            start = stop + 1
        end
        Solution(new_energies, new_states, new_probs, new_degeneracy, partial_sol.largest_discarded_probability)
    end
    _merge
end

no_merge(partial_sol::Solution) = partial_sol

function bound_solution(partial_sol::Solution, max_states::Int, merge_strategy=no_merge)
    if length(partial_sol.probabilities) <= max_states
        probs = vcat(partial_sol.probabilities, -Inf)
        k = length(probs)
    else
        probs = partial_sol.probabilities
        k = max_states + 1
    end

    indices = partialsortperm(probs, 1:k, rev=true)
    new_max_discarded_prob = max(partial_sol.largest_discarded_probability, probs[indices[end]])

    indices = @view indices[1:k-1]

    merge_strategy(Solution(
        partial_sol.energies[indices],
        partial_sol.states[indices],
        partial_sol.probabilities[indices],
        partial_sol.degeneracy[indices],
        new_max_discarded_prob
    ))
end

#TODO: incorporate "going back" move to improve alghoritm
function low_energy_spectrum(network::AbstractGibbsNetwork, max_states::Int, merge_strategy=no_merge)
    # Build all boundary mps
    @showprogress "Preprocesing: " for i ∈ network.nrows:-1:1 dressed_mps(network, i) end

    # Start branch and bound search
    sol = empty_solution()
    @showprogress "Search: " for _ ∈ 1:nv(network_graph(network))
        sol = branch_solution(sol, network)
        sol = bound_solution(sol, max_states, merge_strategy)
        # _clear_cache(network, sol)
    end

    # Translate variable order (from network to factor graph)
    inner_perm = sortperm([
        factor_graph(network).reverse_label_map[idx]
        for idx ∈ network.vertex_map.(iteration_order(network))
    ])

    # Sort using energies as keys
    outer_perm = sortperm(sol.energies)
    Solution(
        sol.energies[outer_perm],
        [σ[inner_perm] for σ ∈ sol.states[outer_perm]],
        sol.probabilities[outer_perm],
        sol.degeneracy[outer_perm],
        sol.largest_discarded_probability
    )
end


text:
#######################################################################
# Slice
#######################################################################
struct Slice{H<:Union{Hamiltonian,ParametricHamiltonian},A<:AbstractSparseMatrix,A´<:AbstractSparseMatrix,X,X´}
    h::H
    fullmat::A
    fullmatflat::A´
    axes::X
    axesflat::X´
end

(s::Slice{<:ParametricHamiltonian})(; params...) =
    Slice(s.h(; params...), s.fullmat, s.fullmatflat, s.axes, s.axesflat)

function Base.show(io::IO, m::MIME"text/plain", s::Slice)
    i = get(io, :indent, "")
    print(io, i, "Slice of ")
    show(io, m, parent(s))
    print(io, i, "\n",
"$i  Sliced size      : $(size(s))
$i  Flat size        : $(sizeflat(s))")
end

orbitalstructure(s::Slice) = orbitalstructure(parent(s))

matrixtype(s::Slice) = matrixtype(parent(s))

blocktype(s::Slice) = blocktype(parent(s))

blockeltype(s::Slice) = blockeltype(parent(s))

axesflat(s::Slice) = s.axesflat

axesflat(s::Slice, n) = axesflat(s)[n]

Base.parent(s::Slice) = s.h

Base.axes(s::Slice) = s.axes

Base.axes(s::Slice, n) = axes(s)[n]

Base.size(s::Slice, n...) = size(view(s.fullmat, axes(s)...), n...)

sizeflat(s::Slice, n...) = size(view(s.fullmatflat, axesflat(s)...), n...)

# bloch of slices

# This is faster than bloch!(similarmatrix(s), parent(s), args...), because
# copy!(sparse, view(sparse,...)) is slower than copy(view(sparse, ...))
function bloch(s::Slice, args...)
    h = parent(s)
    m = bloch(h, args...)
    T = eltype(m)
    axes = axes_for_eltype(s, T)
    m´ = copy(view(m, axes...))
    return m´
end

function bloch!(matrix::AbstractMatrix{T}, s::Slice, args...) where {T}
    bfull = bloch!(fullmat_for_eltype(s, T), s.h, args...)
    axes = axes_for_eltype(s, T)
    fast_sparse_copy!(matrix, bfull, axes)
    return matrix
end

# fastpath for dense arrays
similarmatrix(s::Slice, ::Type{S}) where {T<:SMatrix,S<:StridedMatrix{T}} = S(undef, size(s)...)
similarmatrix(s::Slice, ::Type{S}) where {T<:Number, S<:StridedMatrix{T}} = S(undef, sizeflat(s)...)
similarmatrix(s::Slice, ::Type{<:StridedMatrix}) = Matrix{blocktype(s)}(undef, size(s)...)

# in any other case (sparse), we need to add structural zeros
function similarmatrix(s::Slice, args...)
    mat = similarmatrix(s.h, args...)
    axes = axes_for_eltype(s, eltype(mat))
    v = view(mat, axes...)
    mat´ = copy(v)
    return mat´
end

axes_for_eltype(s, ::Type{<:SMatrix}) = s.axes
axes_for_eltype(s, ::Type{<:Number}) = s.axesflat

fullmat_for_eltype(s, ::Type{<:SMatrix}) = s.fullmat
fullmat_for_eltype(s, ::Type{<:Number}) = s.fullmatflat

# indexing into a slice
Base.getindex(s::Slice{<:Hamiltonian}, dn...) = parent(s)[dn...][axes(s)...]

# indexing of Hamiltonian and ParametricHamiltonian
function Base.getindex(ph::Union{Hamiltonian,ParametricHamiltonian}, is, js)
    h = parent(ph)
    fullmatflat, fullmat = similarmatrix_flat_unflat(h)
    rows = all_unflat_indices(h, is)
    cols = all_unflat_indices(h, js)
    flatrows = all_flat_indices(h, is)
    flatcols = all_flat_indices(h, js)
    axes = (rows, cols)
    axesflat = (flatrows, flatcols)
    return Slice(ph, fullmat, fullmatflat, axes, axesflat)
end

# axes and axesflat need not be vectors of indices. Can be Colon or other containers.
all_unflat_indices(h, i::Integer) = all_unflat_indices(h, i:i)
all_unflat_indices(h, is::Colon) = is
all_unflat_indices(h, is::AbstractRange) = is
all_unflat_indices(h, is::SiteSelector) = collect(siteindices(h.lattice, is))
all_unflat_indices(h, is) = siteindices(h.lattice, is)

all_flat_indices(::Hamiltonian{<:Any,<:Any, <:Number}, is) = is
all_flat_indices(h::Hamiltonian{<:Any,<:Any, <:Number}, is::SiteSelector) = collect(siteindices(h.lattice, is))
all_flat_indices(h::Hamiltonian{<:Any,<:Any, <:Number}, i::Integer) = all_flat_indices(h, i:i)
all_flat_indices(h::Hamiltonian{<:Any,<:Any, <:SMatrix}, i::Integer) = all_flat_indices(h, i:i)
all_flat_indices(::Hamiltonian{<:Any,<:Any, <:SMatrix}, is::Colon) = is

function all_flat_indices(h::Hamiltonian{<:Any,<:Any, <:SMatrix}, is)
    o = orbitalstructure(h)
    is´ = Int[]
    for i in siteindices(h.lattice, is)
        append!(is´, flatindices(i, o))
    end
    return is´
end

function all_flat_indices(h::Hamiltonian{<:Any,<:Any, <:SMatrix}, is::AbstractRange)
    o = orbitalstructure(h)
    is´ = intersect(is, 1:nsites(h))
    is == is´ || throw(ArgumentError("Indices $is are out of range $(1:nsites(h))"))
    istart = flatoffset_site(first(is´), o) + 1
    oend, N = flatoffsetorbs_site(last(is´), o)
    iend = oend+N
    return istart:iend
end

similarmatrix_flat_unflat(h::Hamiltonian{<:Any,<:Any, <:SMatrix}) =
    similarmatrix(h, flatten), similarmatrix(h)

function similarmatrix_flat_unflat(h::Hamiltonian{<:Any,<:Any, <:Number})
    mat = similarmatrix(h)
    return mat, mat
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
module ThirdPartyMock
    using Pretend, Test
    Pretend.activate()

    @mockable Base.sin(x::Real)

    fakesin(x::Real) = 10

    function test()
        apply(sin => fakesin) do
            @test sin(1.0) == 10
        end
        @test sin(1.0) ≈ 0.8414709848078965
    end
end

using .ThirdPartyMock
ThirdPartyMock.test()


text:
using PowerModels

function parse_prob_file(file::String, data::Dict{AbstractString,Any})

    file = string(split(file, ".m")[1], ".prob")
    data_string = readstring(open(file))
    data_lines = split(data_string, '\n')
    parsed_matrices = []
    last_index = length(data_lines)
    index = 1

    while index <= last_index
        line = strip(data_lines[index])

        if length(line) <= 0 || strip(line)[1] == '%'
            index += 1
            continue
        end

        if contains(line, "[")
            matrix = PowerModels.parse_matrix(data_lines, index)
            push!(parsed_matrices, matrix)
            index += matrix["line_count"] - 1
        end
        index += 1
    end

    parsed_matrix = parsed_matrices[1]
    prob = []
    if parsed_matrix["name"] == "mpc.branchprobabilities"
        for row in parsed_matrix["data"]
            push!(prob, parse(Float64, row[3]))
        end
    end

    for i in 1:length(prob)
        for branch in data["branch"]
            if branch["index"] == i
                branch["prob"] = prob[i]
                break
            end
        end
    end

    return
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
__precompile__()

module AdaptiveDistanceFields

using RegionTrees
import RegionTrees: needs_refinement, refine_data
using Interpolations: interpolate!,
                      extrapolate,
                      AbstractInterpolation,
                      BSpline,
                      OnGrid,
                      Linear,
                      Line
using StaticArrays: SVector

export AdaptiveDistanceField,
       evaluate,
       ConvexMesh

include("interpolation.jl")
include("adaptivesampling.jl")

struct AdaptiveDistanceField{C <: Cell} <: Function
    root::C
end

function AdaptiveDistanceField(signed_distance::Function, origin::AbstractArray,
              widths::AbstractArray,
              rtol=1e-2,
              atol=1e-2)
    refinery = SignedDistanceRefinery(signed_distance, atol, rtol)
    boundary = HyperRectangle(origin, widths)
    root = Cell(boundary, refine_data(refinery, boundary))
    adaptivesampling!(root, refinery)
    AdaptiveDistanceField(root)
end

function (field::AdaptiveDistanceField)(point::AbstractArray)
    evaluate(field.root, point)
end

end