text:
module Econlite

using Parameters
using Markdown
using Unrolled
using LinearAlgebra
using Statistics
import Random.seed!

# Single good preferences
export IsoElastic, LogUtility, CARA, RiskNeutral
export util, invutil, dutil, invdutil
include("utility.jl")

# Aggregators
export aggr, invaggr, partial_aggr, invpartial_aggr, gradient_aggr, gradient_aggr!
export Armington, CobbDouglas, AggLinear
include("aggregators.jl")

# Taxes
export net, gross, tax
export GouveiaStrauss, AdValorem, NullTax
include("taxes.jl")

# Maximum
export generalized_argmax!
export SumNormalize, SumRaw
include("maximum.jl")

# Generate grid
export gridpoints
include("gridpoints.jl")

# Iteration
export SolutionObject, IterationSpec, DisplayMod, NoDisplay
export within_tolerance, contract!
include("iteration.jl")

# Markov chains
export IndexedMarkovChain, MarkovChain, MarkovChainExtended
export expect_markov, transition, markov_invariant, draw_next,
    markov_genpath, markov_genpath_seeded, markov_genpath!, mean,
    indeptensor, makemean1!
include("markovchain.jl")

export rouwenhorst, rouwenhorst_mc
include("rouwenhorst.jl")

end


text:
export  model_problem,
        test_problem,
        accumulate_current_maps,
        calculate_cum_current_map,
        calculate_max_current_map

using Test

 """
 Construct nodemap specific to a connected component
 """
function construct_local_node_map(nodemap, component, polymap)
    local_nodemap = zeros(eltype(nodemap), size(nodemap))
    idx = findall(in(component), nodemap)
    local_nodemap[idx] = nodemap[idx]
    if nodemap == local_nodemap
        return local_nodemap
    end
    _construct_local_nodemap(local_nodemap, polymap, idx)
end

function _construct_local_nodemap(local_nodemap, polymap, idx)
    if isempty(polymap)
        i = findall(x ->x!=0, local_nodemap)
        local_nodemap[i] = 1:length(i)
        return local_nodemap
    else
        local_polymap = zeros(eltype(local_nodemap), size(local_nodemap))
        local_polymap[idx] = polymap[idx]
        return construct_node_map(local_nodemap, local_polymap)
    end
end

"""
Define model circuitscape problem - helps in tests
"""
function model_problem(T, s)

    # Cell map is the uniform 1 across the grid
    cellmap = ones(T, s, s)

    # Nodemap is just 1:endof(cellmap)
    nodemap = reshape(1:lastindex(cellmap), s, s) |> collect

    # Construct graph with 4 neighbors
    G = construct_graph(cellmap, nodemap, true, true)

    # Return laplacian
    laplacian(G)
end
model_problem(s::Integer) = model_problem(Float64, s)

# Testing utilities

function test_problem(str)
    base_path = joinpath(dirname(pathof(Circuitscape)), "..", "test", "input")
    str2 = replace(str, ".ini" => "")
    if occursin("sgVerify", str)
        config_path = joinpath(base_path, "raster", "pairwise", replace(str2, "sgVerify" => ""), str)
    elseif occursin("mgVerify", str)
        config_path = joinpath(base_path, "raster", "advanced", replace(str2, "mgVerify" => ""), str)
    elseif occursin("oneToAll", str)
        config_path = joinpath(base_path, "raster", "one_to_all", replace(str2, "oneToAllVerify" => ""), str)
    elseif occursin("allToOne", str)
        config_path = joinpath(base_path, "raster", "all_to_one", replace(str2, "allToOneVerify" => ""), str)
    elseif occursin("sgNetworkVerify", str)
        config_path = joinpath(base_path, "network", str)
    elseif occursin("mgNetworkVerify", str)
        config_path = joinpath(base_path, "network", str)
    end
end

# Utility funciton for calling Circuitscape with Cholmod
function compute_cholmod(str, batch_size = 5)
    cfg = parse_config(str)
    T = cfg["precision"] in SINGLE ? Float32 : Float64
    V = cfg["use_64bit_indexing"] in TRUELIST ? Int64 : Int32
    if T == Float32
        cswarn("Cholmod supports only double precision. Changing precision to double.")
        T = Float64
    end
    cfg["solver"] = "cholmod"
    cfg["cholmod_batch_size"] = string(batch_size)
    _compute(T, V, cfg)
end

function compute_single(str)
    cfg = parse_config(str)
    cfg["precision"] = "single"
    V = cfg["use_64bit_indexing"] in TRUELIST ? Int64 : Int32
    T = Float32
    if cfg["solver"] == "cholmod"
        cswarn("Cholmod supports only double precision. Changing precision to double.")
        T = Float64
    end
    _compute(T, V, cfg)
end

function compute_parallel(str, n_processes = 2)
    cfg = parse_config(str)
    cfg["parallelize"] = "true"
    cfg["max_parallel"] = "$(n_processes)"
    compute(cfg)
end

function get_output_flags(cfg)

    # Output flags
    write_volt_maps = cfg["write_volt_maps"] in TRUELIST
    write_cur_maps = cfg["write_cur_maps"] in TRUELIST
    write_cum_cur_map_only = cfg["write_cum_cur_map_only"] in TRUELIST
    write_max_cur_maps = cfg["write_max_cur_maps"] in TRUELIST
    set_null_currents_to_nodata = cfg["set_null_currents_to_nodata"] in TRUELIST
    set_null_voltages_to_nodata = cfg["set_null_voltages_to_nodata"] in TRUELIST
    compress_grids = cfg["compress_grids"] in TRUELIST
    log_transform_maps = cfg["log_transform_maps"] in TRUELIST

    o = OutputFlags(write_volt_maps, write_cur_maps,
                    write_cum_cur_map_only, write_max_cur_maps,
                    set_null_currents_to_nodata, set_null_voltages_to_nodata,
                    compress_grids, log_transform_maps)
end

# Helps start new processes from the INI file
function myaddprocs(n)
    addprocs(n)
    @everywhere Core.eval(Main, :(using Circuitscape))
end

# Reads the directory with the current maps
# and accumulates all current maps
function accumulate_current_maps(path, f)
    dir = dirname(path)
    base = basename(path)

    # If base file has a dot
    name = split(base, ".out")[1]

    cmap_list = readdir(dir) |>
                    x -> filter(y -> startswith(y, "$(name)_"), x) |>
                    x -> filter(y -> occursin("_curmap_", y), x)
    isempty(cmap_list) && return

    headers = ""
    first_file = joinpath(dir, cmap_list[1])
    nrow = 0
    ncol = 0

    # Read the headers from the first file
    open(first_file, "r") do f

        # Get num cols
        str = readline(f)
        headers = headers * str * "\n"
        ncol = split(str)[2] |> x -> parse(Int, x)

        # Get num rows
        str = readline(f)
        headers = headers * str * "\n"
        nrow = split(str)[2] |> x -> parse(Int, x)

        # Just append the rest
        for i = 3:6
            headers = headers * readline(f) * "\n"
        end
    end

    accum = zeros(nrow, ncol)
    for file in cmap_list
        csinfo("Accumulating $file", cfg["suppress_messages"] in TRUELIST)
        cmap_path = joinpath(dir, file)
        cmap = readdlm(cmap_path, skipstart = 6)
        f_in_place!(accum, cmap, f)
    end
    for i in eachindex(accum)
        if accum[i] < -9999
            accum[i] = -9999
        end
    end

    name =  if isequal(f, +)
                "cum"
            elseif isequal(f, max)
                "max"
            end

    accum_path = joinpath(dir, name * "_$(name)_curmap.asc")
    csinfo("Writing to $accum_path", cfg["suppress_messages"] in TRUELIST)
    open(accum_path, "w") do f
        write(f, headers)
        writedlm(f, round.(accum, digits=8), ' ')
    end

end

function f_in_place!(accum, cmap, f)
    accum .= f.(accum, cmap)
end

calculate_cum_current_map(path) = accumulate_current_maps(path, +)
calculate_max_current_map(path) = accumulate_current_maps(path, max)

function postprocess_cum_curmap!(accum)
    for i in eachindex(accum)
        if accum[i] < -9999
            accum[i] = -9999
        end
    end
end

mycsid() = myid() - minimum(workers()) + 1

function initialize_cum_maps(cellmap::Matrix{T}, max = false) where T
    cum_curr = Vector{SharedMatrix{T}}(undef,nprocs())
    for i = 1:nprocs()
        cum_curr[i] = SharedArray(zeros(T, size(cellmap)...))
    end
    max_curr = Vector{SharedMatrix{T}}()
    if max
        max_curr = Vector{SharedMatrix{T}}(undef,nprocs())
        for i = 1:nprocs()
            max_curr[i] = SharedArray(fill(T(-9999), size(cellmap)...))
        end
    end
    cum_branch_curr = Vector{SharedVector{T}}()
    cum_node_curr = Vector{SharedVector{T}}()

    Cumulative(cum_curr, max_curr,
        cum_branch_curr, cum_node_curr)
end

function initialize_cum_vectors(v::Vector{T}) where T
    cum_curr = Vector{SharedMatrix{T}}()
    max_curr = Vector{SharedMatrix{T}}()
    cum_branch_curr = Vector{SharedVector{T}}(undef,nprocs())
    cum_node_curr = Vector{SharedVector{T}}(undef,nprocs())
    for i = 1:nprocs()
        cum_branch_curr[i] = SharedArray(zeros(T, size(v)...))
        cum_node_curr[i] = SharedArray(zeros(T, size(v)...))
    end

    Cumulative(cum_curr, max_curr,
        cum_branch_curr, cum_node_curr)
end

function runtests(f = compute)

    str = if f == compute_single
            "Single"
          else
            "Double"
          end

    is_single = false
    tol = 1e-6
    str == "Single" && (is_single = true; tol = 1e-4)

    @testset "$str Precision Tests" begin
        @testset "Network Pairwise" begin
            # Network pairwise tests
            for i = 1:3
                @info("Testing sgNetworkVerify$i")
                r = f("input/network/sgNetworkVerify$(i).ini")
                x = readdlm("output_verify/sgNetworkVerify$(i)_resistances.out")
                valx = x[2:end, 2:end]
                valr = r[2:end, 2:end]
                @test sum(abs2, valx - valr) < tol
                pts_x = x[2:end,1]
                pts_r = r[2:end,1]
                @test pts_x .+ 1 == pts_r
                compare_all_output("sgNetworkVerify$(i)", is_single)
                @info("Test sgNetworkVerify$i passed")
            end
        end

        @testset "Network Advanced" begin
            # Network advanced tests
            for i = 1:3
                @info("Testing mgNetworkVerify$i")
                r = f("input/network/mgNetworkVerify$(i).ini")
                x = readdlm("output_verify/mgNetworkVerify$(i)_voltages.txt")
                @. x[:,1] = x[:,1] + 1
                @test sum(abs2, x - r) < tol
                compare_all_output("mgNetworkVerify$(i)", is_single)
                @info("Test mgNetworkVerify$i passed")
            end
        end

        @testset "Raster Pairwise" begin
            # Raster pairwise tests
            for i = 1:17
                # Weird windows 32 stuff
                if i == 16 && Sys.WORD_SIZE == 32
                    continue
                end

                @info("Testing sgVerify$i")
                r = f("input/raster/pairwise/$i/sgVerify$(i).ini")
                x = readdlm("output_verify/sgVerify$(i)_resistances.out")
                _x = readdlm("output/sgVerify$(i)_resistances.out")
                # x = x[2:end, 2:end]
                @test sum(abs2, _x - r) < tol
                @test sum(abs2, x - r) < tol
                compare_all_output("sgVerify$(i)", is_single)
                @info("Test sgVerify$i passed")
            end
        end

        @testset "Raster Advanced" begin
            # Raster advanced tests
            for i in 1:5
                @info("Testing mgVerify$i")
                r = f("input/raster/advanced/$i/mgVerify$(i).ini")
                x = readdlm("output_verify/mgVerify$(i)_voltmap.asc"; skipstart = 6)
                @test sum(abs2, x - r) < 1e-4
                # compare_all_output("mgVerify$(i)")
                @info("Test mgVerify$i passed")
            end
        end

        @testset "Raster One to All" begin
            # Raster one to all test
            for i in 1:13
                @info("Testing oneToAllVerify$i")
                r = f("input/raster/one_to_all/$i/oneToAllVerify$(i).ini")
                x = readdlm("output_verify/oneToAllVerify$(i)_resistances.out")
                # x = x[:,2]
                @test sum(abs2, x - r) < tol
                compare_all_output("oneToAllVerify$(i)", is_single)
                @info("Test oneToAllVerify$i passed")
            end
        end

        @testset "Raster All to One" begin
            # Raster all to one test
            for i in 1:12
                @info("Testing allToOneVerify$i")
                r = f("input/raster/all_to_one/$i/allToOneVerify$(i).ini")
                x = readdlm("output_verify/allToOneVerify$(i)_resistances.out")
                # x = x[:,2]

                @test sum(abs2, x - r) < tol
                compare_all_output("allToOneVerify$(i)", is_single)
                @info("Test allToOneVerify$i passed")
            end
        end
    end
end

function compare_all_output(str, is_single = false)

    gen_list, list_to_comp = generate_lists(str)
    tol = is_single ?  1e-4 : 1e-6

    for f in gen_list
        !occursin("_", f) && continue
        occursin("resistances", f) && continue

        @info("Testing $f")

        # Raster output files
        if endswith(f, "asc")
            r = read_aagrid("output/$f")
            x = get_comp(list_to_comp, f)
            @test compare_aagrid(r, x, tol)
            @info("Test $f passed")

        # Network output files
        elseif occursin("Network", f)

            # Branch currents
            if occursin("branch", f)
                r = read_branch_currents("output/$f")
                x = !startswith(f, "mg") ? get_network_comp(list_to_comp, f) : readdlm("output_verify/$f")
                @test compare_branch(r, x, tol)
                @info("Test $f passed")

            # Node currents
            else
                r = read_node_currents("output/$f")
                x = !startswith(f, "mg") ? get_network_comp(list_to_comp, f) : readdlm("output_verify/$f")
                @test compare_node(r, x, tol)
                @info("Test $f passed")
            end
        end
    end

end

list_of_files(str, pref) = readdir(pref) |> y -> filter(x -> startswith(x, "$(str)_"), y)
generate_lists(str) = list_of_files(str, "output/"), list_of_files(str, "output_verify/")
read_branch_currents(str) = readdlm(str)
read_node_currents(str) = readdlm(str)

read_aagrid(file) = readdlm(file, skipstart = 6) # Will change to 6

compare_aagrid(r::Matrix{T}, x::Matrix{T}, tol = 1e-6) where T = sum(abs2, x - r) < tol

function get_comp(list_to_comp, f)
    outfile = ""
    if f in list_to_comp
        outfile = "output_verify/$f"
    end
    readdlm(outfile; skipstart = 6)
end

function get_network_comp(list_to_comp, f)
    s = split(f, ['_', '.'])
    for i = 1:size(s, 1)
        if all(isnumeric, s[i])
            f = replace(f, "_$(s[i])"=>"_$(parse(Int, s[i]) - 1 |> string)", count=1)
        end
    end
    @assert isfile("output_verify/$f")
    readdlm("output_verify/$f")
end

function compare_branch(r, x, tol = 1e-6)
    @. x[:,1] = x[:,1] + 1
    @. x[:,2] = x[:,2] + 1
    sum(abs2, sortslices(r, dims=1) - sortslices(x, dims=1)) < tol
end

function compare_node(r, x, tol = 1e-6)
    @. x[:,1] = x[:,1] + 1
    sum(abs2, sortslices(r, dims=1) - sortslices(x, dims=1)) < tol
end

ext: jl


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test
using Compat

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_parametric)

import ParametricTypes.TemplateType, ParametricTypes.NonTypeParam

p1 = TemplateType{ParametricTypes.P1, ParametricTypes.P2}()
p2 = TemplateType{ParametricTypes.P2, ParametricTypes.P1}()

println("Dumping object p1:")
dump(p1)

@test ParametricTypes.get_first(p1) == 1
@test ParametricTypes.get_second(p2) == 1
@test typeof(ParametricTypes.get_first(p1)) == Int32
@test typeof(ParametricTypes.get_second(p2)) == Int32

@test ParametricTypes.get_first(p2) == 10.
@test ParametricTypes.get_second(p1) == 10.
@test typeof(ParametricTypes.get_first(p2)) == Float64
@test typeof(ParametricTypes.get_second(p1)) == Float64

nontype1 = ParametricTypes.NonTypeParam{Int32, Int32(1)}()
@test ParametricTypes.get_nontype(nontype1) == 1

nontype2 = ParametricTypes.NonTypeParam{UInt32, UInt32(2)}()
@test ParametricTypes.get_nontype(nontype2) == UInt32(2)

nontype3 = ParametricTypes.NonTypeParam{Int32, Int32(1)}(3)
@test ParametricTypes.get_nontype(nontype3) == 3

nontype4 = ParametricTypes.NonTypeParam{Int64, Int64(64)}()
@test ParametricTypes.get_nontype(nontype4) == Int64(64)


text:
module SettingsBundles

export
    # Types
    SettingsBundle, SettingsYAMLFile, SettingsJSONFile,
    # Functions
    add!, flatten

abstract SettingsSource <: Associative
include("sources/SettingsYAMLFile.jl")
include("sources/SettingsJSONFile.jl")

"""
The SettingsBundle type; a dictionary on staroids

Also, it only accepts symbols as inputs to make searching
for things more uniform. So I guess I should say sources
must be Associative objects with symbol based keys.
"""
type SettingsBundle <: SettingsSource
    sources::Array{Associative}

    SettingsBundle() = new(Array(Associative, 0))
    SettingsBundle(sources...) = new(reverse([i for i in sources])) # Make sure it's not a Tuple
end

Base.showdict(s::SettingsSource; kw...) = Base.showdict(STDOUT, s; kw...)
function Base.showdict(out::IO, s::SettingsSource; limit::Bool=false, compact::Bool=false)
    k = keys(s)
    print_with_color(:bold, out, "$(typeof(s)) with $(length(k)) sources")
    length(k) == 0 ? print_with_color(:bold, out, "\n") : print_with_color(:bold, out, ":\n")
    for key in k
        print_with_color(:bold, out, "  $(key)\t=> $(s[key])\n")
    end
end

"""
Addes a settings source to a SettingsBundle. The most recently
added source is the highest priority
"""
add!(item::SettingsBundle, source::Associative) = insert!(item.sources, 1, source)

"""
Takes and existing settings objects, and produces a dictionary
"""
function flatten(item::SettingsBundle)
    flat = Dict()
    for key in keys(item)
        flat[key] = item[key]
    end
    return flat
end

"""
getindex for a Settings object involve going through the `sources`
array and going until we find a settings object that has the desired
key. Once found, it returns the objects. If non of the sources has
the key, throws an error
"""
function Base.getindex(item::SettingsBundle, key)
    for source in item.sources
        try
            return source[key]
        catch
            # Do nothing. We still have more sources to check :)
        end
    end
    # Not in any of our sources. Shame.
    throw(KeyError(key))
end

"You connot set items in a SettingsBundle by direct assigment"
function Base.setindex!(item::SettingsBundle, key, value)
    throw(error("Elements in SettingsBundles cannot be modified directly"))
end

"Compiles a list of available keys"
function Base.keys(item::SettingsBundle)
    source_keys = map(keys, item.sources)
    return union(source_keys...)
end

end # module


text:
import StatsFuns: log2π

@node NormalMeanPrecision Stochastic [ out, (μ, aliases = [ mean ]), (τ, aliases = [ invcov, precision ]) ]

conjugate_type(::Type{ <: NormalMeanPrecision }, ::Type{ Val{ :out } }) = NormalMeanPrecision
conjugate_type(::Type{ <: NormalMeanPrecision }, ::Type{ Val{ :μ } })   = NormalMeanPrecision
conjugate_type(::Type{ <: NormalMeanPrecision }, ::Type{ Val{ :τ } })   = Gamma

@average_energy NormalMeanPrecision (q_out::Any, q_μ::Any, q_τ::Any) = begin
    μ_mean, μ_var     = mean_var(q_μ)
    out_mean, out_var = mean_var(q_out)
    return (log2π - mean(log, q_τ) + mean(q_τ) * (μ_var + out_var + abs2(μ_mean - out_mean))) / 2
end

@average_energy NormalMeanPrecision (q_out_μ::MultivariateNormalDistributionsFamily, q_τ::Any) = begin
    out_μ_mean, out_μ_cov = mean_cov(q_out_μ)
    return (log2π - mean(log, q_τ) + mean(q_τ) * (out_μ_cov[1,1] + out_μ_cov[2,2] - out_μ_cov[1,2] - out_μ_cov[2,1] + abs2(out_μ_mean[1] - out_μ_mean[2]))) / 2
end


text:
"""
    cudnnNormalizationForward(x, xmean, xvar, bias, scale; o...)
    cudnnNormalizationForward!(y, x, xmean, xvar, bias, scale; o...)

Return batch normalization applied to `x`:

    y .= ((x .- mean(x; dims)) ./ sqrt.(epsilon .+ var(x; dims))) .* scale .+ bias  # training
    y .= ((x .- xmean) ./ sqrt.(epsilon .+ xvar)) .* scale .+ bias                  # inference

Bias and scale are trainable parameters, xmean and xvar are modified to collect statistics
during training and treated as constants during inference. Note that during inference the
values given by xmean and xvar arguments are used in the formula whereas during training the
actual mean and variance of the minibatch are used in the formula: the xmean/xvar arguments
are only used to collect statistics. In the original paper bias is referred to as beta and
scale as gamma (Batch Normalization: Accelerating Deep Network Training by Reducing Internal
Covariate Shift, S. Ioffe, C. Szegedy, 2015).

Keyword arguments:
* `epsilon = 1e-5`: epsilon value used in the normalization formula
* `exponentialAverageFactor = 0.1`: factor used in running mean/variance calculation: `runningMean = runningMean*(1-factor) + newMean*factor`
* `training = false`: boolean indicating training vs inference mode
* `mode::cudnnNormMode_t = CUDNN_NORM_PER_CHANNEL`: Per-channel layer is based on the paper. In this mode `scale` etc. have dimensions (1,1,C,1). The other alternative is `CUDNN_NORM_PER_ACTIVATION` where `scale` etc. have dimensions `(W,H,C,1)`.
* `algo::cudnnNormAlgo_t = CUDNN_NORM_ALGO_STANDARD`: The other alternative, `CUDNN_NORM_ALGO_PERSIST`, triggers the new semi-persistent NHWC kernel when certain conditions are met (see cudnn docs).
* `normOps::cudnnNormOps_t = CUDNN_NORM_OPS_NORM`: Currently the other alternatives, `CUDNN_NORM_OPS_NORM_ACTIVATION` and `CUDNN_NORM_OPS_NORM_ADD_ACTIVATION` are not supported.
* `z = nothing`: for residual addition to the result of the normalization operation, prior to the activation (will be supported when CUDNN_NORM_OPS_NORM_ADD_ACTIVATION is supported)
* `groupCnt = 1`: Place holder for future work, should be set to 1 now
* `alpha = 1; beta = 0`: scaling parameters: return `alpha * new_y + beta * old_y`

"""
cudnnNormalizationForward, cudnnNormalizationForward!

# Public methods
cudnnNormalizationForward(x, xmean, xvar, bias, scale; o...) = cudnnNormalizationForwardWithDefaults(x, xmean, xvar, bias, scale; o...)
cudnnNormalizationForward!(y, x, xmean, xvar, bias, scale; o...) = cudnnNormalizationForwardWithDefaults(x, xmean, xvar, bias, scale; y, o...)

# Private method
function cudnnNormalizationForwardWithDefaults(
    x, mean, variance, bias, scale;

    # Inference parameters:
    y = similar(x),
    z = nothing, # for residual addition to the result of the normalization operation, prior to the activation
    mode::cudnnNormMode_t = CUDNN_NORM_PER_CHANNEL, # Per-channel layer is based on the paper Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift, S. Ioffe, C. Szegedy, 2015.
    normOps::cudnnNormOps_t = CUDNN_NORM_OPS_NORM,  # Currently CUDNN_NORM_OPS_NORM_ACTIVATION and CUDNN_NORM_OPS_NORM_ADD_ACTIVATION are only supported in the NHWC layout (training,backward), not supported (inference)
    algo::cudnnNormAlgo_t = CUDNN_NORM_ALGO_STANDARD, # trigger the new semi-persistent NHWC kernel when CUDNN_NORM_ALGO_PERSIST
    alpha::Real = 1,
    beta::Real = 0,
    epsilon::Real = Cdouble(1e-5), # Has to be >= 0. Should be the same in forward and backward functions.
    groupCnt::Integer = Cint(1),   # Place hold for future work, should be set to 1 now

    # Training-only parameters:
    training = false,
    exponentialAverageFactor::Real = Cdouble(0.1),
    savedMean = nothing, # Optionally save intermediate results from the forward pass here - can be reused to speed up backward pass. NULL if unused.
    savedInvVariance = nothing,

    # Activation parameters:
    activationMode::cudnnActivationMode_t = CUDNN_ACTIVATION_IDENTITY,
    activationReluNanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
    activationCoef::Real = 1,
    activationDesc::Union{Nothing,cudnnActivationDescriptor} = (normOps == CUDNN_NORM_OPS_NORM ? nothing : cudnnActivationDescriptor(activationMode, activationReluNanOpt, Cdouble(activationCoef))),

    # Tensor descriptors:
    format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x; format),
    yDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(y; format),
    zDesc::Union{Nothing,cudnnTensorDescriptor} = (z === nothing ? nothing : cudnnTensorDescriptor(z; format)),
    normScaleBiasDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(scale; format),
    normMeanVarDesc::Union{Nothing,cudnnTensorDescriptor} = (mean === nothing ? nothing : cudnnTensorDescriptor(mean; format)),

    # Temporary space used in training:
    workspace = nothing,
    reserveSpace = nothing,
    dx = Ref{Any}(),
    dscale = Ref{Any}(),
    dbias = Ref{Any}(),
    dz = Ref{Any}(),
)
    @assert epsilon >= 0 && exponentialAverageFactor >= 0  "epsilon and exponentialAverageFactor should be non-negative."
    @assert groupCnt == 1  "Currently only groupCnt=1 is supported."
    @assert normOps === CUDNN_NORM_OPS_NORM "Currently only normOps=CUDNN_NORM_OPS_NORM is supported."
    alpha, beta = (a->scalingParameter(eltype(x),a)).((alpha, beta))
    # Backward called separately on each variable. We will calculate all gradients on first call. Use `dready` to avoid subsequent calls.
    dready = Ref{Bool}(false)   # this will be turned to `true` by the first backward call.
    cudnnNormalizationForwardAD(x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready)
end

# AD method:
function cudnnNormalizationForwardAD(x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready)
    issimilar(x,y) = (typeof(x) === typeof(y) && (x === nothing || size(x) === size(y)))
    if training
        mean === nothing ? savedMean = nothing : savedMean === nothing ? savedMean = similar(mean) : @assert issimilar(mean, savedMean)
        variance === nothing ? savedInvVariance = nothing : savedInvVariance === nothing ? savedInvVariance = similar(variance) : @assert issimilar(variance, savedInvVariance)
        workspaceSize, reserveSpaceSize = cudnnNormalizationTempSpaceSizes(mode, normOps, algo, xDesc, zDesc, yDesc, normScaleBiasDesc, activationDesc, normMeanVarDesc, groupCnt)
        if reserveSpaceSize > 0 && reserveSpace === nothing; reserveSpace = cudnnTempSpace(reserveSpaceSize); end
        @assert sizeof(reserveSpace) >= reserveSpaceSize  "reserveSpace should be at least $(reserveSpaceSize) bytes"
        if workspaceSize > 0 && workspace === nothing; workspace = cudnnTempSpace(workspaceSize); end
        @assert sizeof(workspace) >= workspaceSize  "workspace should be at least $(workspaceSize) bytes"
        cudnnNormalizationForwardTraining(handle(), mode, normOps, algo, alpha, beta, xDesc, x, normScaleBiasDesc, scale, bias, exponentialAverageFactor, something(normMeanVarDesc,C_NULL), something(mean,CU_NULL), something(variance,CU_NULL), epsilon, something(savedMean,CU_NULL), something(savedInvVariance,CU_NULL), something(activationDesc,C_NULL), something(zDesc,C_NULL), something(z,CU_NULL), yDesc, y, something(workspace,CU_NULL), sizeof(workspace), something(reserveSpace,CU_NULL), sizeof(reserveSpace), groupCnt)
    else
        @assert mean !== nothing && variance !== nothing && normMeanVarDesc !== nothing "normalization mean and variance are required in inference mode."
        cudnnNormalizationForwardInference(handle(), mode, normOps, algo, alpha, beta, xDesc, x, normScaleBiasDesc, scale, bias, normMeanVarDesc, mean, variance, something(zDesc,C_NULL), something(z,CU_NULL), something(activationDesc,C_NULL), yDesc, y, epsilon, groupCnt)
    end
    return y
end

# Helper functions
function cudnnNormalizationTempSpaceSizes(mode, normOps, algo, xDesc, zDesc, yDesc, normScaleBiasDesc, activationDesc, normMeanVarDesc, groupCnt)
    workspaceSize, reserveSpaceSize = Ref{Csize_t}(0), Ref{Csize_t}(0)
    cudnnGetNormalizationForwardTrainingWorkspaceSize(handle(), mode, normOps, algo, xDesc, something(zDesc,C_NULL), yDesc, normScaleBiasDesc, something(activationDesc,C_NULL), something(normMeanVarDesc,C_NULL), workspaceSize, groupCnt)
    cudnnGetNormalizationTrainingReserveSpaceSize(handle(), mode, normOps, algo, something(activationDesc,C_NULL), xDesc, reserveSpaceSize, groupCnt)
    workspaceSize[], reserveSpaceSize[]
end


text:
module TSML

include("system.jl")
using .System

include("types.jl")
using .TSMLTypes

include("utils.jl")
using .Utils

include("dataproc.jl")
using .DataProc

include("transformers.jl")
using .TSMLTransformers

include("baseline.jl")
using .BaselineAlgos

if LIB_SKL_AVAILABLE # from System module
    include("scikitlearn.jl")
    using .SKLearners
end

if LIB_CRT_AVAILABLE # from System module
    include("caret.jl")
    using .CaretLearners
end

include("multilearner.jl")
using .MultiLearners

include("decisiontree.jl")
using .DecisionTreeLearners

include("datareader.jl")
using .DataReaders

include("datawriter.jl")
using .DataWriters

include("statifier.jl")
using .Statifiers

include("monotonicer.jl")
using .Monotonicers

include("cliwrapper.jl")
using .CLIWrappers
export tsmlrun

include("tsclassifier.jl")
using .TSClassifiers

include("outliernicer.jl")
using .Outliernicers

end # module


text:
using Test
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.SurfaceFluxes.UniversalFunctions

@testset "SurfaceFluxes" begin
    struct EarthParameterSet <: AbstractEarthParameterSet end
    ps = EarthParameterSet()

    @testset "UniversalFunctions" begin
        @test Pr_0_Businger(ps) ≈ 0.74
        @test a_m_Businger(ps) ≈ 4.7
        @test a_h_Businger(ps) ≈ 4.7
        @test Pr_0_Gryanik(ps) ≈ 0.98
        @test a_m_Gryanik(ps) ≈ 5.0
        @test a_h_Gryanik(ps) ≈ 5.0
        @test b_m_Gryanik(ps) ≈ 0.3
        @test b_h_Gryanik(ps) ≈ 0.4
        @test Pr_0_Grachev(ps) ≈ 0.98
        @test a_m_Grachev(ps) ≈ 5.0
        @test a_h_Grachev(ps) ≈ 5.0
        @test b_m_Grachev(ps) ≈ a_m_Grachev(ps) / 6.5
        @test b_h_Grachev(ps) ≈ 5.0
        @test c_h_Grachev(ps) ≈ 3.0
    end
end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# mutable struct of FixedJoint is defined in Modia3D/src/Composition/object3d.jl,
# since Object3D references FixedJoint and FixedJoint references Object3D

function computeKinematics!(joint::FixedJoint, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D    = obj.parent
   noTranslation::Bool = obj.r_rel ≡ ModiaMath.ZeroVector3D
   noRotation::Bool    = obj.R_rel ≡ ModiaMath.NullRotation
   acc::Bool           = obj.computeAcceleration

   obj.r_abs = noTranslation ? parent.r_abs : parent.r_abs + parent.R_abs'*obj.r_rel
   obj.R_abs = noRotation    ? parent.R_abs : obj.R_rel*parent.R_abs

   if analysis == ModiaMath.DynamicAnalysis
      dynamics::Object3Ddynamics       = obj.dynamics
      parentDynamics::Object3Ddynamics = parent.dynamics

      if noTranslation
         dynamics.v0 = parentDynamics.v0
         dynamics.a0 = parentDynamics.a0
      else
         dynamics.v0 = parentDynamics.v0 + parent.R_abs'*cross(parentDynamics.w, obj.r_rel)
         #if acc
             dynamics.a0 = parentDynamics.a0 + parent.R_abs'*(cross(parentDynamics.z, obj.r_rel) +
                                                              cross(parentDynamics.w, cross(parentDynamics.w, obj.r_rel)))
         #else
         #    dynamics.a0 = parentDynamics.a0
         #end
      end

      if noRotation
         dynamics.w = parentDynamics.w
         dynamics.z = parentDynamics.z
      else
         dynamics.w = obj.R_rel*parentDynamics.w
         #if acc
             dynamics.z = obj.R_rel*parentDynamics.z
         #else
         #    dynamics.z = parentDynamics.z
         #end
      end
   end
   return nothing
end

function computeForceTorqueAndResidue!(joint::FixedJoint, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D = obj.parent

   if obj.R_rel ≡ ModiaMath.NullRotation
      # Fixed translation
      parent.dynamics.f += obj.dynamics.f
      parent.dynamics.t += obj.dynamics.t + cross(obj.r_rel, obj.dynamics.f)
   else
      # Fixed translation and rotation
      parent.dynamics.f += obj.R_rel'*obj.dynamics.f
      parent.dynamics.t += obj.R_rel'*obj.dynamics.t - cross(obj.r_rel, parent.dynamics.f)
   end
   return nothing
end


text:
"""
    inf(A::IntervalMatrix{T}) where {T}

Return the infimum of an interval matrix `A`, which corresponds to taking the
element-wise infimum of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the infima of each element in `A`.
"""
function inf(A::IntervalMatrix{T}) where {T}
    return map(inf, A.mat)
end

"""
    sup(A::IntervalMatrix{T}) where {T}

Return the supremum of an interval matrix `A`, which corresponds to taking the
element-wise supremum of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the suprema of each element in `A`.
"""
function sup(A::IntervalMatrix{T}) where {T}
    return map(sup, A.mat)
end

"""
    mid(A::IntervalMatrix{T}) where {T}

Return the midpoint of an interval matrix `A`, which corresponds to taking the
element-wise midpoint of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the midpoints of each element in `A`.
"""
function mid(A::IntervalMatrix{T}) where {T}
    return map(mid, A.mat)
end

"""
    radius(A::IntervalMatrix{T}) where {T}

Return the radius of an interval matrix `A`, which corresponds to taking the
element-wise radius of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the radii of each element in `A`.
"""
function radius(A::IntervalMatrix{T}) where {T}
    return map(radius, A.mat)
end

"""
    diam(A::IntervalMatrix{T}) where {T}

Return a matrix whose entries describe the diameters of the intervals.

### Input

- `A` -- interval matrix

### Output

A matrix `B` of the same shape as `A` such that `B[i, j] == diam(A[i, j])` for
each `i` and `j`.
"""
function diam(A::IntervalMatrix{T}) where {T}
    return map(diam, A.mat)
end

"""
    midpoint_radius(A::IntervalMatrix{T}) where {T}

Split an interval matrix ``A`` into two scalar matrices ``C`` and ``S``
such that ``A = C + [-S, S]``.

### Input

- `A` -- interval matrix

### Output

A pair `(C, S)` such that the entries of `C` are the central points and the
entries of `S` are the (nonnegative) radii of the intervals in `A`.
"""
midpoint_radius(A::IntervalMatrix) = (mid(A), radius(A))

real(M::IntervalMatrix) = IntervalMatrix(real(M.mat))
imag(M::IntervalMatrix) = IntervalMatrix(imag(M.mat))


text:
# To run: `julia dfs_test.jl`
include("dfs.jl")

nodes = [
     1 => [2, 7, 8],
     2 => [1, 3, 6],
     3 => [2, 4, 5],
     4 => [3],
     5 => [3],
     6 => [2],
     7 => [1],
     8 => [1, 9, 12],
     9 => [8, 10, 11],
    10 => [9],
    11 => [9],
    12 => [8]
]

visited = Number[]
Dfs.visit(1, nodes, n -> push!(visited, n))

print_joined(STDOUT, visited, ", ")
println()

max_stars_repo_name: adlawson/search-algorithms

max_issues_repo_name: adlawson/search-algorithms

max_forks_repo_name: adlawson/search-algorithms


text:
@test SMeasureTest(τ = 1) isa SMeasureTest
@test SMeasureTest(τ = 1) isa CausalityTest

# Some random time series
npts = 100
xr, yr = rand(npts), rand(npts);

# Initialise test, specifying embedding dimension, emebdding lag
# and number of nearest neighbors
test = SMeasureTest(m = 2, τ = 1, K = 2:5)
@test causality(xr, yr, test) isa Vector{T} where T


text:
#=
    MainResult.jl

    Perform ABC-SMC, then ABC-MCMC, to produce data for the main results.

    WARNING: This script takes ~1 day to run on a 4-core machine. For convenience, the
    results are saved in the `JLD2` file `MainResult.jld2`.

=#

using Inference
using Model

using CSV, DataFrames, DataFramesMeta
using Distributions
using JLD2
using Plots
using StatsBase
using StatsPlots

##############################################################
## Setup problem
##############################################################

include("MainResult_Setup.jl")

##############################################################
## ABC SMC
##############################################################

P = abc_smc(dist,prior,maxsteps=30;n=1000,ptarget=0.005,α=0.75)

##############################################################
## Determine discrepancy function
##############################################################

# The threshold is chosen so the "best fit particle" will be accepted ~50% of the time
D = [dist(minimum(P).θ) for i = 1:100]
ε = quantile(D,0.5)
K = d -> d < ε ? 1.0 : 0.0

##############################################################
## ABC MCMC
##############################################################

# Pilot chains (1,000,000 iterations, skip=100)
@time C_pilot = abc_mcmc(dist,prior,
    minimum(P).θ,
    2.38^2 * cov(P) / ndims(P),
    K,
    n=10000,param_names=param_names,nchains=4,thinning=100
)

# Tuned chains (10,000,000 iterations, skip=100)
@time C = abc_mcmc(dist,prior,
    minimum(P).θ,
    2.38^2 * cov(Matrix(C_pilot)') / ndims(P),
    K,
    n=100000,param_names=param_names,nchains=4,thinning=100
)

##############################################################
## Find best fit
##############################################################

Θ = Matrix(C[1:1000:end])
D = [mean(dist(θ[:]) for i = 1:100) for θ in eachcol(Θ)]
θ = Θ[:,findmin(D)[2]]

##############################################################
## Save results
##############################################################

@save "Results/MainResult.jld2" P ε C_pilot C θ

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using Test, SpikeSynchrony

####### ANALYTICAL TEST ########
y1 = [1]
y2 = [2]
ta = [0, 0, 1, 1, 2, 2, 3, 3]
# I calculated S by hand for these trains:
Sa = [0, 0, 5//9, 3//8, 3//8, 5//9, 0, 0]

t, S = SPIKE_distance_profile(y1, y2; tf = 3)
@test t == ta
@test S ≈ Sa    atol = 1e-15

D_S_a = sum( 1//2 * (Sa[2:end] .+ Sa[1:end-1]) .* (ta[2:end] .- ta[1:end-1]))//(ta[end] - ta[1])
D_S = SpikeSynchrony.trapezoid_integral(t, S)/(t[end] - t[1])
@test D_S ≈ D_S_a    atol = 1e-15

####### Second Analytic test ########
y1 = [2.0, 5.0, 8.0]
y2 = [1, 5, 9]
t, S  = SPIKE_distance_profile(y1, y2; t0 = 0, tf = 10)

analytic_t = [0, 1, 2, 5, 8, 9, 10]
analytic_S = [0, 0, 0.555555555555555580,  0.222222222222222210,
11/36, 25/98, 0, 0, 25/98, 11/36, 0.222222222222222210,  0.555555555555555580, 0, 0]

@test S ≈ analytic_S  atol = 1e-15

####### TEST WITH ARBITRARY REAL DATA ########
t0 = 0; tf = 500
y1 = cumsum(rand(10:10:tf, 100))
y2 = cumsum(rand(10:10:tf, 100))
t, S  = SPIKE_distance_profile(y1, y2)
t, S2 = SPIKE_distance_profile(y2, y1)
t, S0 = SPIKE_distance_profile(y1, y1)
# tr, Sr = SPIKE_distance_profile(abs.(reverse(y1) .- y1[end]), abs.(reverse(y1) .- y1[end]),

# Fundamental tests from the definition of S and D_S
@test issorted(t)
@test (0 .≤ S) == trues(length(S))
@test S == S2
@test S0 == zeros(length(S0))

D_S = SpikeSynchrony.trapezoid_integral(t, S)/(t[end]-t[1]) # == SPIKE_distance(y1, y2)
@test 0 ≤ D_S ≤ 1


text:
promote_rule(::Type{ArbComplex{P}}, ::Type{ArbReal{P}}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{ArbFloat{P}}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbReal{P}}, ::Type{ArbFloat{P}}) where {P} = ArbReal{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{Mag}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbReal{P}}, ::Type{Mag}) where {P} = ArbReal{P}
promote_rule(::Type{ArbFloat{P}}, ::Type{Mag}) where {P} = ArbFloat{P}

promote_type(::Type{ArbFloat{P}}, ::Type{BigInt}) where {P} = ArbFloat{P}
promote_type(::Type{ArbReal{P}}, ::Type{BigInt}) where {P} = ArbReal{P}
promote_type(::Type{ArbComplex{P}}, ::Type{BigInt}) where {P} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{BigInt}) where {P} = promote_type(ArbFloat{P}, BigInt)
promote_rule(::Type{ArbReal{P}}, ::Type{BigInt}) where {P} = promote_type(ArbReal{P}, BigInt)
promote_rule(::Type{ArbComplex{P}}, ::Type{BigInt}) where {P} = promote_type(ArbComplex{P}, BigInt)

promote_type(::Type{BigFloat},::Type{ArbFloat{P}}) where {P} = ArbFloat{P}
promote_type(::Type{BigFloat},::Type{ArbReal{P}}) where {P} = ArbReal{P}
promote_type(::Type{BigFloat},::Type{ArbComplex{P}}) where {P} = ArbComplex{P}
promote_type(::Type{ArbFloat{P}}, ::Type{BigFloat}) where {P} = ArbFloat{P}
promote_type(::Type{ArbReal{P}}, ::Type{BigFloat}) where {P} = ArbReal{P}
promote_type(::Type{ArbComplex{P}}, ::Type{BigFloat}) where {P} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbFloat{P}, BigFloat)
promote_rule(::Type{ArbReal{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbReal{P}, BigFloat)
promote_rule(::Type{ArbComplex{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbComplex{P}, BigFloat)

promote_rule(::Type{ArbFloat{P}}, ::Type{I}) where {P, I<:Integer} = ArbFloat{P}
promote_rule(::Type{ArbReal{P}}, ::Type{I}) where {P, I<:Integer} = ArbReal{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{I}) where {P, I<:Integer} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbFloat{P}
promote_rule(::Type{ArbReal{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbReal{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbComplex{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{C}) where {P, C<:Complex} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, x::ArbFloat{Q}) where {P,Q} =
    P>Q ? ArbFloat{P} : ArbFloat{Q}
promote_rule(::Type{ArbReal{P}}, x::ArbReal{Q}) where {P,Q} =
    P>Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbComplex{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}
promote_rule(::Type{ArbFloat{P}}, x::ArbReal{Q}) where {P,Q} =
    P>Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbFloat{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}
promote_rule(::Type{ArbReal{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}

promote_rule(::Type{ArbFloat{P}},::Type{ArbFloat{Q}}) where {P,Q} = P<Q ? ArbFloat{P} : ArbFloat{Q}
promote_rule(::Type{ArbReal{P}},::Type{ArbReal{Q}}) where {P,Q} = P<Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbComplex{P}},::Type{ArbComplex{Q}}) where {P,Q} = P<Q ? ArbComplex{P} : ArbComplex{Q}

convert(::Type{ArbFloat}, x::ArbFloat{P}) where {P} = x
convert(::Type{ArbReal}, x::ArbReal{P}) where {P} = x
convert(::Type{ArbComplex}, x::ArbComplex{P}) where {P} = x
convert(::Type{ArbFloat}, x::ArbReal{P}) where {P} = ArbFloat{P}(x)
convert(::Type{ArbFloat}, x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x))
convert(::Type{ArbReal}, x::ArbComplex{P}) where {P} = real(x)
convert(::Type{ArbReal}, x::ArbFloat{P}) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex}, x::ArbFloat{P}) where {P} = ArbComplex{P}(x)

convert(::Type{ArbReal{P}}, x::ArbComplex{P}) where {P} = real(x)
convert(::Type{ArbFloat{P}}, x::ArbComplex{P}) where {P} = midpoint(real(x), ArbFloat{P})
convert(::Type{Mag}, x::ArbComplex{P}) where {P} = radius(real(x), Mag)

convert(::Type{ArbComplex{P}}, x::ArbReal{P}) where {P} = ArbComplex{P}(x)
convert(::Type{ArbComplex{P}}, x::ArbFloat{P}) where {P} = ArbComplex{P}(ArbReal{P}(x))
convert(::Type{ArbComplex{P}}, x::Mag) where {P} = ArbComplex{P}(ArbFloat{P}(x))

convert(::Type{ArbFloat{P}}, x::ArbReal{P}) where {P} = midpoint(x, ArbFloat{P})
convert(::Type{Mag}, x::ArbReal{P}) where {P} = radius(x, Mag)

convert(::Type{ArbReal{P}}, x::ArbFloat{P}) where {P} = ArbReal{P}(x)
convert(::Type{ArbReal{P}}, x::Mag) where {P} = ArbReal{P}(ArbFloat{P}(x))

convert(::Type{Mag}, x::ArbFloat{P}) where {P} = Mag(x)
convert(::Type{ArbFloat{P}}, x::Mag) where {P} = ArbFloat{P}(x)

convert(::Type{ArbFloat{P}}, x::BigFloat) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::BigFloat) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::BigFloat) where {P} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::BigInt) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::BigInt) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::BigInt) where {P} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::I) where {P, I<:Integer} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::I) where {P, I<:Integer} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::I) where {P, I<:Integer} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::AbstractFloat) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::AbstractFloat) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::AbstractFloat) where {P} = ArbComplex{P}(x)
convert(::Type{ArbComplex{P}}, x::Complex) where {P} = ArbComplex{P}(x)


text:
using Plots
using Printf
# %%
#f = (a::Float64,b::Float64,n::Int64) -> n==1 ? (b+a/b)/2 : (f(a,b,n-1) + a/f(a,b,n-1))/2

function f(a::Float64, b::Float64, n::Int64)
    if n == 1
        return (b+a/b)/2
    else
        return (f(a, b, n-1) + a/f(a, b, n-1))/2
    end
end

function nohint_f(a, b, n)
    if n == 1
        return (b+a/b)/2
    else
        return (f(a, b, n-1) + a/f(a, b, n-1))/2
    end
end

g = n -> f(3.0,3.0,n)
nohint_g = n -> nohint_f(3.0,3.0,n)

for n in 1:20
    @printf("n=%d g(n)=%f\n",n, g(n))
end

#%%
anim = @animate for n=1:20
    seq = collect(1:n)
    plot(seq,g.(seq))
end
gif(anim, "anim_fps10.gif", fps = 10)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using HPAT

HPAT.DistributedPass.set_debug_level(3)

@acc hpat function index_test(numCenter, file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points", file_name)
    D,N = size(points) # number of features, instances
    centroids :: Array{Float64,2} = rand(D, numCenter)

    dist::Array{Array{Float64,1},1} = [ Float64[sqrt(sum((points[:,i+j]-centroids[:,j]).^2)) for j in 1:numCenter] for i in 1:(N-numCenter)]
    labels::Array{Int,1} = [indmin(dist[i]) for i in 1:(N-numCenter)]
    return sum(labels)
end

function main()

end

using HDF5
using MPI
if MPI.Comm_rank(MPI.COMM_WORLD)==0
    h5write("test1_1.hdf5", "/points", [1 2 3 4; 5 6 7 8])
end
a = index_test(2, "test1_1.hdf5")
println(HPAT.get_saved_array_partitionings())

if MPI.Comm_rank(MPI.COMM_WORLD)==0
        rm("test1_1.hdf5")
end


text:
using PyCall

const bpy = PyNULL()

const context = PyNULL()
const data = PyNULL()

const object_ops = PyNULL()
const object_context = PyNULL()

const render_ops = PyNULL()
const render_context = PyNULL()

function __init__()
    copy!(bpy, pyimport("bpy"))

    copy!(context, bpy.context)
    copy!(data, bpy.data)

    copy!(object_ops, bpy.ops.object)
    copy!(object_context, context.scene.objects)

    copy!(render_ops, bpy.ops.render)
    copy!(render_context, context.scene.render)
end

save(filename::AbstractString) = bpy.ops.wm.save_mainfile(filepath=filename)


text:
# For once you have the sensitivity by two methods (e.g  both finite-differencing and  AD)
# the code here checks it is correct.
# Goal is to only call `@isapprox` on things that render well
# Note that this must work well both on Differential types and Primal types

"""
    check_equal(actual, expected; kwargs...)

`@test`'s  that `actual ≈ expected`, but breaks up data such that human readable results
are shown on failures.
Understands things like `unthunk`ing `ChainRuleCore.Thunk`s, etc.
All keyword arguments are passed to `isapprox`.
"""
function check_equal(
    actual::Union{AbstractArray{<:Number}, Number},
    expected::Union{AbstractArray{<:Number}, Number};
    kwargs...
)
    @test isapprox(actual, expected; kwargs...)
end

for (T1, T2) in ((AbstractThunk, Any), (AbstractThunk, AbstractThunk), (Any, AbstractThunk))
    @eval function check_equal(actual::$T1, expected::$T2; kwargs...)
        check_equal(unthunk(actual), unthunk(expected); kwargs...)
    end
end

check_equal(::Zero, x; kwargs...) = check_equal(zero(x), x; kwargs...)
check_equal(x, ::Zero; kwargs...) = check_equal(x, zero(x); kwargs...)
check_equal(x::Zero, y::Zero; kwargs...) = @test true

"""
    _can_pass_early(actual, expected; kwargs...)
Used to check if `actual` is basically equal to `expected`, so we don't need to check deeper;
and can just report `check_equal` as passing.

If either `==` or `≈` return true then so does this.
The `kwargs` are passed on to `isapprox`
"""
function _can_pass_early(actual, expected; kwargs...)
    actual == expected && return true
    try
        return isapprox(actual, expected; kwargs...)
    catch err
        # Might MethodError, might DimensionMismatch, might fail for some other reason
        # we don't care, whatever errored it means we can't quit early
    end
    return false
end

function check_equal(actual::AbstractArray, expected::AbstractArray; kwargs...)
    if _can_pass_early(actual, expected)
        @test true
    else
        @test eachindex(actual) == eachindex(expected)
        @testset "$(typeof(actual))[$ii]" for ii in eachindex(actual)
            check_equal(actual[ii], expected[ii]; kwargs...)
        end
    end
end

function check_equal(actual::Composite{P}, expected::Composite{P}; kwargs...) where P
    if _can_pass_early(actual, expected)
        @test true
    else
        all_keys = union(keys(actual), keys(expected))
        @testset "$P.$ii" for ii in all_keys
            check_equal(getproperty(actual, ii), getproperty(expected, ii); kwargs...)
        end
    end
end

function check_equal(
    ::Composite{ActualPrimal}, expected::Composite{ExpectedPrimal}; kwargs...
) where {ActualPrimal, ExpectedPrimal}
    # this will certainly fail as we have another dispatch for that, but this will give as
    # good error message
    @test ActualPrimal === ExpectedPrimal
end

# Some structual differential and a natural differential
function check_equal(actual::Composite{P, T}, expected; kwargs...) where {T, P}
    if _can_pass_early(actual, expected)
        @test true
    else
        @assert (T <: NamedTuple)  # it should be a structual differential if we hit this

        # We are only checking the properties that are in the Composite
        # the natural differential is allowed to have other properties that we ignore
        @testset "$P.$ii" for ii in propertynames(actual)
            check_equal(getproperty(actual, ii), getproperty(expected, ii); kwargs...)
        end
    end
end
check_equal(x, y::Composite; kwargs...) = check_equal(y, x; kwargs...)

# This catches comparisons of Composites and Tuples/NamedTuple
# and gives an error message complaining about that
const LegacyZygoteCompTypes = Union{Tuple,NamedTuple}
check_equal(::C, ::T; kwargs...) where {C<:Composite,T<:LegacyZygoteCompTypes} = @test C === T
check_equal(::T, ::C; kwargs...) where {C<:Composite,T<:LegacyZygoteCompTypes} = @test T === C

# Generic fallback, probably a tuple or something
function check_equal(actual::A, expected::E; kwargs...) where {A, E}
    if _can_pass_early(actual, expected)
        @test true
    else
        c_actual = collect(actual)
        c_expected = collect(expected)
        if (c_actual isa A) && (c_expected isa E)  # prevent stack-overflow
            throw(MethodError, check_equal, (actual, expected))
        end
        check_equal(c_actual, c_expected; kwargs...)
    end
end

"""
_check_add!!_behaviour(acc, val)

This checks that `acc + val` is the same as `add!!(acc, val)`.
It matters primarily for types that overload `add!!` such as `InplaceableThunk`s.

`acc` is the value that has been accumulated so far.
`val` is a deriviative, being accumulated into `acc`.

`kwargs` are all passed on to isapprox
"""
function _check_add!!_behaviour(acc, val; kwargs...)
    # Note, we don't test that `acc` is actually mutated because it doesn't have to be
    # e.g. if it is immutable. We do test the `add!!` return value.
    # That is what people should rely on. The mutation is just to save allocations.
    acc_mutated = deepcopy(acc)  # prevent this test changing others
    check_equal(add!!(acc_mutated, val), acc + val; kwargs...)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Core definitions for interacting with the libuv library from Julia

include(string(length(Core.ARGS) >= 2 ? Core.ARGS[2] : "", "uv_constants.jl"))  # include($BUILDROOT/base/uv_constants.jl)

# convert UV handle data to julia object, checking for null
function uv_sizeof_handle(handle)
    if !(UV_UNKNOWN_HANDLE < handle < UV_HANDLE_TYPE_MAX)
        throw(DomainError(handle))
    end
    ccall(:uv_handle_size,Csize_t,(Int32,),handle)
end

function uv_sizeof_req(req)
    if !(UV_UNKNOWN_REQ < req < UV_REQ_TYPE_MAX)
        throw(DomainError(req))
    end
    ccall(:uv_req_size,Csize_t,(Int32,),req)
end

for h in uv_handle_types
@eval const $(Symbol("_sizeof_",lowercase(string(h)))) = uv_sizeof_handle($h)
end
for r in uv_req_types
@eval const $(Symbol("_sizeof_",lowercase(string(r)))) = uv_sizeof_req($r)
end

uv_handle_data(handle) = ccall(:jl_uv_handle_data,Ptr{Cvoid},(Ptr{Cvoid},),handle)
uv_req_data(handle) = ccall(:jl_uv_req_data,Ptr{Cvoid},(Ptr{Cvoid},),handle)
uv_req_set_data(req,data) = ccall(:jl_uv_req_set_data,Cvoid,(Ptr{Cvoid},Any),req,data)
uv_req_set_data(req,data::Ptr{Cvoid}) = ccall(:jl_uv_req_set_data,Cvoid,(Ptr{Cvoid},Ptr{Cvoid}),req,data)

macro handle_as(hand, typ)
    quote
        data = uv_handle_data($(esc(hand)))
        data == C_NULL && return
        unsafe_pointer_to_objref(data)::($(esc(typ)))
    end
end

associate_julia_struct(handle::Ptr{Cvoid}, @nospecialize(jlobj)) =
    ccall(:jl_uv_associate_julia_struct, Cvoid, (Ptr{Cvoid}, Any), handle, jlobj)
disassociate_julia_struct(uv) = disassociate_julia_struct(uv.handle)
disassociate_julia_struct(handle::Ptr{Cvoid}) =
    handle != C_NULL && ccall(:jl_uv_disassociate_julia_struct, Cvoid, (Ptr{Cvoid},), handle)

# A dict of all libuv handles that are being waited on somewhere in the system
# and should thus not be garbage collected
const uvhandles = IdDict()
const preserve_handle_lock = Threads.SpinLock()
function preserve_handle(x)
    lock(preserve_handle_lock)
    v = get(uvhandles, x, 0)::Int
    uvhandles[x] = v + 1
    unlock(preserve_handle_lock)
    nothing
end
function unpreserve_handle(x)
    lock(preserve_handle_lock)
    v = uvhandles[x]::Int
    if v == 1
        pop!(uvhandles, x)
    else
        uvhandles[x] = v - 1
    end
    unlock(preserve_handle_lock)
    nothing
end

## Libuv error handling ##

struct IOError <: Exception
    msg::AbstractString
    code::Int32
    IOError(msg::AbstractString, code::Integer) = new(msg, code)
end

showerror(io::IO, e::IOError) = print(io, "IOError: ", e.msg)

function _UVError(pfx::AbstractString, code::Integer)
    code = Int32(code)
    IOError(string(pfx, ": ", struverror(code), " (", uverrorname(code), ")"), code)
end

struverror(err::Int32) = unsafe_string(ccall(:uv_strerror,Cstring,(Int32,),err))
uverrorname(err::Int32) = unsafe_string(ccall(:uv_err_name,Cstring,(Int32,),err))

uv_error(prefix::Symbol, c::Integer) = uv_error(string(prefix),c)
uv_error(prefix::AbstractString, c::Integer) = c < 0 ? throw(_UVError(prefix,c)) : nothing

## event loop ##

eventloop() = uv_eventloop::Ptr{Cvoid}
#mkNewEventLoop() = ccall(:jl_new_event_loop,Ptr{Cvoid},()) # this would probably be fine, but is nowhere supported

function process_events()
    return ccall(:jl_process_events, Int32, (Ptr{Cvoid},), eventloop())
end

function uv_alloc_buf end
function uv_readcb end
function uv_writecb_task end
function uv_return_spawn end
function uv_asynccb end
function uv_timercb end

function reinit_stdio()
    global uv_jl_alloc_buf     = @cfunction(uv_alloc_buf, Cvoid, (Ptr{Cvoid}, Csize_t, Ptr{Cvoid}))
    global uv_jl_readcb        = @cfunction(uv_readcb, Cvoid, (Ptr{Cvoid}, Cssize_t, Ptr{Cvoid}))
    global uv_jl_writecb_task  = @cfunction(uv_writecb_task, Cvoid, (Ptr{Cvoid}, Cint))
    global uv_jl_return_spawn  = @cfunction(uv_return_spawn, Cvoid, (Ptr{Cvoid}, Int64, Int32))
    global uv_jl_asynccb       = @cfunction(uv_asynccb, Cvoid, (Ptr{Cvoid},))
    global uv_jl_timercb       = @cfunction(uv_timercb, Cvoid, (Ptr{Cvoid},))

    global uv_eventloop = ccall(:jl_global_event_loop, Ptr{Cvoid}, ())
    global stdin = init_stdio(ccall(:jl_stdin_stream, Ptr{Cvoid}, ()))
    global stdout = init_stdio(ccall(:jl_stdout_stream, Ptr{Cvoid}, ()))
    global stderr = init_stdio(ccall(:jl_stderr_stream, Ptr{Cvoid}, ()))
end

"""
    stdin

Global variable referring to the standard input stream.
"""
:stdin

"""
    stdout

Global variable referring to the standard out stream.
"""
:stdout

"""
    stderr

Global variable referring to the standard error stream.
"""
:stderr

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
"""
Implementations of particle filters, also referred to as sequential importance
sampling, and variants thereof.

Algorithms based on Bayesian Filtering & Smoothing (Särkka, 2014)

Wouter Kouw
19-12-2019
"""

using Distributions
using Random

function sequential_importance_sampling(observations,
                                        nominal_state_transition,
                                        nominal_likelihood;
                                        imp_processp=1.0,
                                        imp_measurementp=1.0,
                                        num_particles=100)
    """
    Sequential importance sampling (Alg. 7.2)

    This filter is built for the following model:
        p(x_k | x_k-1) ~ nominal_state_transition(x_k-1)
        p(y_k | x_k) ~ nominal_likelihood(x_k)
        π(x_k | x_k-1) ~ N(x_k-1, imp_processp^-1)
        π(y_k | x_k) ~ N(x_k, imp_measurementp^-1)
    """

    # Time horizon
    time_horizon = length(observations)

    # Preallocate importance sample array
    samples = zeros(time_horizon+1, num_particles)

    # Preallocate importance weight array
    weights = zeros(time_horizon+1, num_particles)

    # Initialize importance weight
    weights[1,:] = ones(1, num_particles) / num_particles

    # Draw from prior
    samples[1,:] = rand(Normal(0, 1), num_particles)

    for t = 2:time_horizon+1

        for i = 1:num_particles

            "Update importance distribution"

            # Update parameters
            impp_precision = imp_processp + imp_measurementp
            impp_mean = inv(impp_precision)*(imp_processp*samples[t-1, i] + imp_measurementp*observations[t-1])

            # Define new importance posterior distribution
            pi_ti = Normal(impp_mean, sqrt(inv(impp_precision)))

            # Draw new importance sample
            samples[t, i] = rand(pi_ti, 1)[1]

            "Update importance weights"

            # Evaluate nominal distributions
            px = pdf(nominal_state_transition(samples[t-1, i]), samples[t, i])
            py = pdf(nominal_likelihood(samples[t, i]), observations[t-1])

            # Evaluate probability of importance sample under importance dist.
            pπ = pdf(pi_ti, samples[t,i])

            # Update weight
            weights[t, i] = weights[t-1,i]*py*px/pπ

        end

        # Normalize weights
        weights[t,:] = weights[t,:] / sum(weights[t,:])
    end
    return (samples, weights)
end

function sequential_importance_resampling(observations,
                                          nominal_state_transition,
                                          nominal_likelihood;
                                          imp_processp=1.0,
                                          imp_measurementp=1.0,
                                          num_particles=100,
                                          resampling_threshold=0.5)
    """
    Sequential importance resampling (Alg. 7.4)

    This filter is built for the following model:
        p(x_k | x_k-1) ~ nominal_state_transition(x_k-1)
        p(y_k | x_k) ~ nominal_likelihood(x_k)
        π(x_k | x_k-1) ~ N(x_k-1, imp_processp^-1)
        π(y_k | x_k) ~ N(x_k, imp_measurementp^-1)
    """

    # Time horizon
    time_horizon = length(observations)

    # Preallocate importance sample array
    samples = zeros(time_horizon+1, num_particles)

    # Preallocate importance weight array
    weights = zeros(time_horizon+1, num_particles)

    # Initialize importance weight
    weights[1,:] = ones(1, num_particles) / num_particles

    # Draw from prior
    samples[1,:] = rand(Normal(0, 1), num_particles)

    for t = 2:time_horizon+1

        for i = 1:num_particles

            "Update importance distribution"

            # Update parameters
            impp_precision = imp_processp + imp_measurementp
            impp_mean = inv(impp_precision)*(imp_processp*samples[t-1, i] + imp_measurementp*observations[t-1])

            # Define new importance posterior distribution
            pi_ti = Normal(impp_mean, sqrt(inv(impp_precision)))

            # Draw new importance sample
            samples[t, i] = rand(pi_ti, 1)[1]

            "Update importance weights"

            # Evaluate nominal distributions
            px = pdf(nominal_state_transition(samples[t-1, i]), samples[t, i])
            py = pdf(nominal_likelihood(samples[t, i]), observations[t-1])

            # Evaluate probability of importance sample under importance dist.
            pπ = pdf(pi_ti, samples[t,i])

            # Update weight
            weights[t, i] = weights[t-1,i]*(py * px / pπ)
        end

        # Normalize weights
        weights[t,:] = weights[t,:] / sum(weights[t,:])

        "Resample according to weights"

        # Compute effective sample size
        num_particles_effective = 1 ./ sum(weights[t, :].^2)

        # Check whether resampling is necessary
        if num_particles_effective < resampling_threshold * num_particles

            # Define categorical distribution based on weights
            dd = Categorical(weights[t,:])

            # Draw from discrete distribution
            index_set = rand(dd, num_particles)

            # Select samples according to drawn index set
            samples[t,:] = samples[t, index_set]

            # Weights are reset to 1/N
            weights[t,:] = ones(1,num_particles) / num_particles

        end
    end
    return samples, weights
end

function bootstrap_filter(observations,
                          nominal_state_transition,
                          nominal_likelihood;
                          num_particles=100)
    """
    Bootstrap filter (Alg. 7.5)

    This filter is built for the following model:
        p(x_k | x_k-1) ~ nominal_state_transition(x_k-1)
        w_k ~ nominal_likelihood(x_k)
    """

    # Time horizon
    time_horizon = length(observations)

    # Preallocate importance sample array
    samples = zeros(time_horizon+1, num_particles)

    # Preallocate importance weight array
    weights = zeros(time_horizon+1, num_particles)

    # Initialize importance weight
    weights[1,:] = ones(1, num_particles) / num_particles

    # Draw from prior
    samples[1,:] = rand(Normal(0, 1), num_particles)

    for t = 2:time_horizon+1

        for i = 1:num_particles

            # Sample from state transition
            samples[t, i] = rand(nominal_state_transition(samples[t-1, i]), 1)[1]

            # Compute weights
            weights[t, i] = pdf(nominal_likelihood(samples[t, i]), observations[t-1])

        end

        # Normalize weights
        weights[t,:] = weights[t,:] / sum(weights[t,:])

        # Define categorical distribution based on weights
        dd = Categorical(weights[t,:])

        # Draw indices from discrete distribution defined by weights
        index_set = rand(dd, num_particles)

        # Select samples according to drawn index set
        samples[t,:] = samples[t, index_set]

        # Weights are reset to 1/N
        weights[t,:] = ones(1,num_particles) / num_particles

    end
    return (samples, weights)
end

function rao_blackwellized_particle_filter(observations,
                                           nominal_dynamics,
                                           importance_dist,
                                           A, Q, H, R;
                                           num_particles=100,
                                           resampling_threshold=0.5)
    """
    Rao-Blackwellization of sequential importance resampling.

    This filter is built for the following model:

        p(x_k | x_k-1, u_k-1) ~ N(A_k-1(u_k-1) x_k-1, Q_k-1(u_k-1))
        p(y_k | x_k, u_k) ~ N(H_k(u_k) x_k, R_k(u_k))
        p(u_k | u_k-1) ~ nominal_dynamics(u_k-1)
        π(u_k | u_k-1, y_k) ~ importance_dist(u_k-1, y_k)
        p(u_0) ~ N(0,1)

    where we assume A, Q, H, and R are functions of u_k. For example:

        A_k = A(u_k) = u_k ⋅ 1
        Q_k = Q(u_k) = u_k^2
        H_k = H(u_k) = 1/2 ⋅ u_k
        R_k = R(u_k) = 1/2 ⋅ u_k^2
    """

    # Time horizon
    time_horizon = length(observations)

    # Preallocate particle arrays
    samples = zeros(time_horizon+1, num_particles)
    weights = zeros(time_horizon+1, num_particles)
    means = zeros(time_horizon+1, num_particles)
    cvars = zeros(time_horizon+1, num_particles)

    # Draw initial particles
    weights[1,:] = ones(1, num_particles) / num_particles
    samples[1,:] = rand(Normal(0, 0.1), num_particles)
    means[1,:] = rand(Normal(0,1), num_particles)
    cvars[1,:] = rand(Gamma(1,1), num_particles)

    for t = 2:time_horizon+1

        # Allocate predictions of Kalman filter parameters
        mean_ = zeros(num_particles,)
        cvar_ = zeros(num_particles,)

        for i = 1:num_particles

            # Convenience variables
            A_tmini = A(samples[t-1,i])
            Q_tmini = Q(samples[t-1,i])

            "1. Kalman filter predictions"

            # Move Kalman filter forward
            mean_[i] = A_tmini*means[t-1,i]
            cvar_[i] = A_tmini*cvars[t-1,i]*A_tmini + Q_tmini

            "2. Draw new latent variables u"

            # Draw from given importance distribution
            samples[t,i] = rand(importance_dist(samples[t-1,i], observations[t-1]), 1)[1]

            "3. Calculate new weights"

            # Convencience variable
            H_ti = H(samples[t,i])
            R_ti = R(samples[t,i])

            # Marginalized measurement likelihood
            pyu = pdf(Normal(H_ti*mean_[i], sqrt(H_ti*cvar_[i]*H_ti + R_ti)), observations[t-1])
            puu = pdf(nominal_dynamics(samples[t-1,i]), samples[t,i])
            πuu = pdf(importance_dist(samples[t-1,i], observations[t-1]), samples[t,i])

            # Update weights
            weights[t,i] = (pyu * puu / πuu) * weights[t-1,i]
        end

        # Normalize weights
        weights[t,:] = weights[t,:] / sum(weights[t,:])

        for i = 1:num_particles

            # Convenience variables
            H_ti = H(samples[t,i])
            R_ti = R(samples[t,i])

            "4. Kalman filter updates"

            # Difference between prediction and observation
            v = observations[t-1] - H_ti*mean_[i]

            # Updated precision
            S = H_ti * cvar_[i]* H_ti + R_ti

            # Kalman gain
            K = cvar_[i]*H_ti * inv(S)

            # Update Kalman filter parameters
            means[t,i] = mean_[i] + K*v
            cvars[t,i] = cvar_[i] - K*S*K'
        end

        # Compute effective sample size
        num_particles_effective = 1 ./ sum(weights[t, :].^2)

        # Check whether resampling is necessary
        if num_particles_effective < resampling_threshold * num_particles

            # Define categorical distribution based on weights
            dd = Categorical(weights[t,:])

            # Draw from discrete distribution
            index_set = rand(dd, num_particles)

            # Select samples according to drawn index set
            samples[t,:] = samples[t, index_set]

            # Weights are reset to 1/N
            weights[t,:] = ones(1,num_particles) / num_particles

        end
    end
    return (samples, weights, means, cvars)
end

max_stars_repo_name: wmkouw/kalman-filters

max_issues_repo_name: wmkouw/kalman-filters

max_forks_repo_name: wmkouw/kalman-filters


text:
module SfAstronomy

using Unitful, UnitfulAstro

import ..SfUnits: Angle, to_angle

export absolute_magnitude, apparent_magnitude, diffuse_disc_q, diffuse_sphere_q

"""
    absolute_magnitude(radius::Unitful.Length, geometric_albedo::Real)

Compute the absolute magnitude of lambertian disc with the given `radius` and `geometric_albedo`.
"""
function absolute_magnitude(radius::Unitful.Length, geometric_albedo::Real)
	d = ustrip(radius *2 |> u"km")
	return 5 * log10(1329 / (d * sqrt(geometric_albedo)))
end

"""
	diffuse_sphere_q(α::Angle)

Approximate phase integral of an ideal diffuse reflecting sphere for observer angle `α`.
"""
diffuse_sphere_q(α) = (2/3) * ((1 - (α / 180u"°")) * cos(α) + (1/π) * sin(α))

"""
	diffuse_disc_q(α::Angle)

Phase integral of a diffuse reflecting disc for observer angle `α`.

Value will be ≥ 0 and ≤ 1
"""
diffuse_disc_q(α) = abs(cos(α))

"""
	apparent_magnitude(h::Real, d_bs::Unitful.Length, d_bo::Unitful.Length, d_os::Unitful.Length, q_α::Real)

Compute the apparent magnitude of an object with absolute magnitude `h`, distance to the sun `d_bs`, distance to the observer `d_bo`, distance from the observer to the sun `d_os` and phase integral `q_α`.
"""
function apparent_magnitude(h::Real, d_bs::Unitful.Length, d_bo::Unitful.Length, d_os::Unitful.Length, q_α::Real)
	return h + 5 * log10(ustrip(d_bs * d_bo / d_os^2 |> u"km/km")) - 2.5 * log10(q_α)
end

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
struct Event
end

function Event_CreateEvent(event::String)
    return event
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
# ClimateDataIO.jl
#
#   Module for loading typical data file formats found in soil and atmospheric
#	sciences
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.7.0
# 09.12.2016
# Last Edit: 07.01.2022

__precompile__(true)

"""# ClimateDataIO

Tools for loading various file types common to soil, atmospheric, and climate sciences. This includes SLT files from EddyMeas, Aerodyne QCL data, CampbellScientific TOA5 DAT files, Los Gatos Research data files, and Licor GHG files.

For more specific information see each functions' help.

---

### Aerodyne QCL Files
`AerodyneStatus`: Convert the StatusW column in STC files into boolean values

`str_load`: Load Aerodyne STR files

`stc_load`: Load Aerodyne STC files

### Campbell Scientific Text DAT Files
`csci_textload`: Load CampbellScientific text DAT files

`csci_textread`: Load a CampbellScientific text DAT files

`csci_times`: Load the minimum and maximum dates and times in the listed files

### EddyMeas SLT Files
`slt_load`: Load SLT files

`slt_header`: Load SLT header info

`slt_read`: Read SLT file converting everything except mV signals

`slt_config`: Load SLT configuration data files (CFG)

`slt_timeshift`: Shift a set of SLT files' time

`slt_write`: Write data to an SLT file

`slt_trim`: Remove columns from SLT files

### LGR Laser Data
`lgr_load`: Load  Los Gatos Research (LGR) text data files

`lgr_read`: Load a Los Gatos Research (LGR) text data file

### Licor GHG Files
`ghg_load`: Load Licor GHG files

`ghg_read`: Load a Licor GHG file

`licor_split`: Split Licor text data files into smaller compressed GHG files

`LI7200Diagnostic`: Parse LI-7200 diagnostic values

---

#### Requirements
* CSV
* Crayons.Box
* DataFrames
* Dates
* DelimitedFiles
* Printf
* Statistics
* StatsBase
* Test"""
module ClimateDataIO

using Crayons.Box
using CSV
using DataFrames
using Dates
using DelimitedFiles
using Printf
using Statistics
using StatsBase
using Test

import Base.show

export AerodyneStatus,
    csci_textload,
    csci_textread,
    csci_times,
    ghg_load,
    ghg_pack,
    ghg_read,
    LI7200Diagnostic,
    licor_split,
    lgr_load,
    lgr_read,
    show,
    slt_config,
    slt_header,
    slt_load,
    slt_read,
    slt_timeshift,
    slt_write,
    slt_trim,
    stc_load,
    str_load

dir = splitdir(@__FILE__)[1]
include(joinpath(dir, "aerodyne_parsetime.jl"))
include(joinpath(dir, "AerodyneStatus.jl"))
include(joinpath(dir, "csci_textload.jl"))
include(joinpath(dir, "csci_textread.jl"))
include(joinpath(dir, "csci_times.jl"))
include(joinpath(dir, "dirlist.jl"))
include(joinpath(dir, "findnewton.jl"))
include(joinpath(dir, "ghg_load.jl"))
include(joinpath(dir, "ghg_pack.jl"))
include(joinpath(dir, "ghg_read.jl"))
include(joinpath(dir, "lgr_load.jl"))
include(joinpath(dir, "lgr_read.jl"))
include(joinpath(dir, "LI7200Diagnostic.jl"))
include(joinpath(dir, "licor_split.jl"))
include(joinpath(dir, "str_load.jl"))
include(joinpath(dir, "stc_load.jl"))
include(joinpath(dir, "slt_read.jl"))
include(joinpath(dir, "slt_load.jl"))
include(joinpath(dir, "slt_header.jl"))
include(joinpath(dir, "slt_config.jl"))
include(joinpath(dir, "slt_configload.jl"))
include(joinpath(dir, "slt_timeshift.jl"))
include(joinpath(dir, "slt_write.jl"))
include(joinpath(dir, "slt_trim.jl"))
include(joinpath(dir, "ziptextfiles.jl"))

# Define 7zip location constant
if Sys.iswindows() & (VERSION >= VersionNumber("1.3.0"))
    const exe7z = joinpath(splitdir(Sys.BINDIR)[1], "libexec", "7z.exe")
elseif Sys.iswindows() & (VERSION < VersionNumber("1.3.0"))
    const exe7z = joinpath(splitdir(Sys.BINDIR)[1], "7z.exe")
elseif Sys.isunix() & (VERSION >= VersionNumber("1.3.0"))
    const exe7z = joinpath(splitdir(Sys.BINDIR)[1], "libexec", "7z")
elseif Sys.isunix() & (VERSION < VersionNumber("1.3.0"))
    const exe7z = joinpath(splitdir(Sys.BINDIR)[1], "7z")
end
end # End of module

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO

max_forks_repo_name: JeremyRueffer/ClimateDataIO


text:
export ConnectableObservable, connectable, connect

import Base: ==
import Base: show

"""
    ConnectableObservable{D}(subject, source)

A connectable observable encapsulates the multicasting infrastructure with provided subject, but does not immediately subscribe to the source.
It subscribes to the source when its `connect` method is called.

See also: [`connect`](@ref), [`Subscribable`](@ref)
"""
@subscribable struct ConnectableObservable{D, J, S} <: Subscribable{D}
    subject :: J
    source  :: S
end

function on_subscribe!(observable::ConnectableObservable, actor)
    return subscribe!(observable.subject, actor)
end

"""
    connectable(subject::J, source::S) where J where S

Creates a `ConnectableObservable` with a given subject object and a source observable.

# Example

```jldoctest
using Rocket

c = connectable(Subject(Int; scheduler = AsapScheduler()), from(1:3))

subscribe!(c, logger());

connect(c);
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed

```

See also: [`ConnectableObservable`](@ref), [`connect`](@ref), [`subscribe!`](@ref)
"""
connectable(subject::J, source::S) where { J, S } = as_connectable(as_subject(J), as_subscribable(S), subject, source)

as_connectable(::InvalidSubjectTrait, ::InvalidSubscribableTrait, subject, source) = throw(InvalidSubjectTraitUsageError(subject))
as_connectable(::InvalidSubjectTrait, _,                          subject, source) = throw(InvalidSubjectTraitUsageError(subject))
as_connectable(_,                     ::InvalidSubscribableTrait, subject, source) = throw(InvalidSubscribableTraitUsageError(source))

as_connectable(::ValidSubjectTrait{D1}, ::SimpleSubscribableTrait{D2},    subject, source) where { D1, D2  } = throw(InconsistentActorWithSubscribableDataTypesError{D2, D1}(source, subject))
as_connectable(::ValidSubjectTrait{D1}, ::ScheduledSubscribableTrait{D2}, subject, source) where { D1, D2  } = throw(InconsistentActorWithSubscribableDataTypesError{D2, D1}(source, subject))

as_connectable(::ValidSubjectTrait{D1},  ::SimpleSubscribableTrait{D2},    subject::J, source::S) where { D1, D2 <: D1, J, S } = ConnectableObservable{D2, J, S}(subject, source)
as_connectable(::ValidSubjectTrait{D1},  ::ScheduledSubscribableTrait{D2}, subject::J, source::S) where { D1, D2 <: D1, J, S } = ConnectableObservable{D2, J, S}(subject, source)

"""
    connect(connectable::ConnectableObservable)

When `connect` is called, the subject passed to the multicast operator is subscribed to the source and the subject’s observers receive the multicast notifications,
which fits our basic mental model of stream multicasting. Returns a subscription.

See also: [`connectable`](@ref), [`ConnectableObservable`](@ref)
"""
connect(connectable::ConnectableObservable) = subscribe!(connectable.source, connectable.subject)

Base.:(==)(c1::ConnectableObservable{D},  c2::ConnectableObservable{D})  where { D      } = c1.subject == c2.subject && c1.source == c2.source
Base.:(==)(c1::ConnectableObservable{D1}, c2::ConnectableObservable{D2}) where { D1, D2 } = false

Base.show(io::IO, ::ConnectableObservable{D}) where D = print(io, "ConnectableObservable($D)")


text:
```
Snapshot information for a currency.
```
function ticker(id::String)
  res = HTTP.get(ENDPOINT * "products/$id/ticker")
  rawData = JSON.parse(String(res.body))
  df = DataFrame(rawData)
  df.ask = parse.(Float64, df.ask)
  df.bid = parse.(Float64, df.bid)
  df.size = parse.(Float64, df.size)
  df.volume = parse.(Float64, df.volume)
  df.price = parse.(Float64, df.price)
  #df.time = DateTime.(chop.(df.time), dateformat"yyyy-mm-ddTHH:MM:SS.ssssss")
  df
end


text:
using Term
import Term.Progress: ProgressBar, start!, update!, stop!, with, @track, addjob!, render
import Term.Consoles: clear

tprint(hLine("progress bars"; style = "blue"))

"""
This example shows how to create progress bars visualizations in Term.

Progress bars in Term can handle multiple simultaneous tasks ('jobs')
whose progress we can monitor (think of nested for loops). The
basic workflow involves creating a progress bars, adding jobs
and updating them to update the visual display

"""

pbar = ProgressBar()
job = addjob!(pbar; N = 100)
start!(pbar)
for i in 1:100
    update!(job)
    sleep(0.01)
    i % 25 == 0 && println("We can print from here too")
    render(pbar)
end
stop!(pbar)

"""
Or with multiple jobs
"""

pbar = ProgressBar()
job = addjob!(pbar; N = 100)
job2 = addjob!(pbar; N = 50)
job3 = addjob!(pbar; N = 25)

start!(pbar)
for i in 1:100
    update!(job)
    i % 2 == 0 && update!(job2)
    i % 3 == 0 && update!(job3)
    sleep(0.01)
    render(pbar)
end
stop!(pbar)

"""
As you can see, the progress bar display stays as the bottom while the text is printed above.
To make that happen Term needs to mess around with your terminal, calling
`stop!` ensures that the normal terminal behavior is restored. If you
don't call `stop!` (e.g. because your code gave an error), then the
terminal is not correctly restored. To avoid that we can use `with(pbar)`, this
will ensure that the progress bar is correctly stopped no matter what happens.

"""

pbar = ProgressBar(; expand = true)
job = addjob!(pbar; N = 100)

with(pbar) do
    for i in 1:100
        update!(job)
        sleep(0.01)
        i % 25 == 0 && println("We can print from here too")
        render(pbar)
    end
end

"""
As you can see, we don't need to call `start!` and `stop!` any more, `with` takes
care of that (even if there's an error in the loop).

Also, you might have noticed that `expand=true` makes the progress bar expand
to fill in all the available space. An alternative is to set the `width` value.

In addition to `width`, you can also use a convenient macro called @track
to wrap any loop with a progress bar display.
"""

@track for i in 1:100
    sleep(0.01)
end

"""
You can customiza what kind of information should show in your progress bar and what it
should look like. Each bit of information (the text, the percentage count...) is a
"column" and you can choose which columns to display and in which order:
"""

for details in (:minimal, :default, :detailed)
    pbar = ProgressBar(; columns = details, width = 140)
    with(pbar) do
        job = addjob!(pbar; N = 100)
        for i in 1:100
            update!(job)
            sleep(0.02)
        end
    end
end

"""
Or you can create your choice of columns
"""

import Term.Progress: CompletedColumn, SeparatorColumn, ProgressColumn, DescriptionColumn

mycols = [DescriptionColumn, CompletedColumn, SeparatorColumn, ProgressColumn]
cols_kwargs = Dict(:DescriptionColumn => Dict(:style => "red bold"))

pbar = ProgressBar(; columns = mycols, columns_kwargs = cols_kwargs, width = 140)
with(pbar) do
    job = addjob!(pbar; N = 100)
    for i in 1:100
        update!(job)
        sleep(0.02)
    end
end

"""
As you can see you can use `cols_kwargs` to pass additional info to each column.

Sometimes you are not sure how many iterations your loop should run for.
In that case you can use spinners!
"""

import Term.Progress: SPINNERS

for spinner in keys(SPINNERS)
    columns_kwargs = Dict(
        :SpinnerColumn => Dict(:spinnertype => spinner, :style => "bold green"),
        :CompletedColumn => Dict(:style => "dim"),
    )

    pbar = ProgressBar(; columns = :spinner, columns_kwargs = columns_kwargs)
    with(pbar) do
        job = addjob!(pbar; description = "{orange1}$spinner...")
        for i in 1:250
            update!(job)
            sleep(0.0025)
        end
    end
end

"""
You can make your own column too! You just need to define a type and an update! method for them

"""

using Random

using Term.Progress
import Term.Progress:
    AbstractColumn, DescriptionColumn, CompletedColumn, SeparatorColumn, ProgressColumn
import Term.Segments: Segment
import Term.Measures: Measure

struct RandomColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    style::String

    function RandomColumn(job::ProgressJob; style = "red")
        txt = Segment(randstring(6), style)
        return new(job, [txt], txt.measure, style)
    end
end

function Progress.update!(col::RandomColumn, color::String, args...)::String
    txt = Segment(randstring(6), col.style)
    return txt.text
end

mycols = [DescriptionColumn, CompletedColumn, SeparatorColumn, ProgressColumn, RandomColumn]
cols_kwargs = Dict(:RandomColumn => Dict(:style => "green bold"))

pbar = ProgressBar(; columns = mycols, columns_kwargs = cols_kwargs, width = 140)
with(pbar) do
    job = addjob!(pbar; N = 100)
    for i in 1:100
        update!(job)
        sleep(0.02)
    end
end

"""
This is an example of how to use a progress bar to show file
downloading/loading progress.
"""

import Term.Progress: SeparatorColumn, ProgressColumn, DescriptionColumn, DownloadedColumn

FILESIZE = 2342341
CHUNK = 2048
nsteps = Int64(ceil(FILESIZE / CHUNK))
@info nsteps

mycols = [DescriptionColumn, SeparatorColumn, ProgressColumn, DownloadedColumn]

pbar = ProgressBar(; columns = mycols, width = 140)
job = addjob!(pbar; N = FILESIZE)

with(pbar) do
    for i in 1:nsteps
        update!(job; i = CHUNK)
        sleep(0.001)
    end
end


text:
using LinearAlgebra: I
import TensorOperations

export CompositeQSystem, add_hamiltonian!, add_lindblad!

# tensor products of quantum systems
mutable struct CompositeQSystem
    subsystems::Vector
    fixed_Hs::Vector{Tuple} # tuple of Matrix and expansion indices
    parametric_Hs::Vector{Tuple} # tuple of Functions and expansion indices
    fixed_Ls::Vector{Tuple} # tuple of Matrix and expansion indices for collapse operators
    parametric_Ls::Vector{Tuple} # tuple of Matrix, expansion indices, and time-dependent function for collapse operators
    dimension::Int
end

CompositeQSystem(qs) = CompositeQSystem(qs, [], [], [], [], prod(dimension(q) for q in qs))

# helper functions for CompositeQSystems
dimension(cqs::CompositeQSystem) = cqs.dimension
dim(cqs::CompositeQSystem) = dimension(cqs)

find_indices(cqs::CompositeQSystem, s::QSystem) = findall([sub == s for sub in cqs.subsystems])
find_indices(cqs::CompositeQSystem, s_label::AbstractString) = findall([label(sub) == s_label for sub in cqs.subsystems])
function find_indices(cqs::CompositeQSystem, s::Union{Vector{<:QSystem}, Vector{<:AbstractString}})
    return vcat([find_indices(cqs, _s) for _s in s]...)
end

""" Add a subsystem static Hamiltonian matrix to a CompositeQSystem """
function add_hamiltonian!(cqs::CompositeQSystem, ham::AbstractMatrix{<:Number}, acting_on::Union{Q, Array{Q}}) where {Q<:QSystem}
    idxs = embed_indices(cqs, acting_on)
    push!(cqs.fixed_Hs, (ham, idxs))
end

""" Add a subystem Hamiltonian to a CompositeQSystem """
add_hamiltonian!(cqs::CompositeQSystem, qs::QSystem) = add_hamiltonian!(cqs, hamiltonian(qs), qs)
add_hamiltonian!(cqs::CompositeQSystem, qs::QSystem, ϕ::Real) = add_hamiltonian!(cqs, hamiltonian(qs, ϕ), qs)

# TODO how to do this dispatch vs adding a fixed Hamiltonian - Jameson says not to do this https://discourse.julialang.org/t/functions-and-callable-methods/2983/3
""" Add a time parameterized subsystem Hamiltonian to a CompositeQSystem """
function add_hamiltonian!(cqs::CompositeQSystem, ham::Function, acting_on::Union{Q, Array{Q}}) where Q<:QSystem
    idxs = embed_indices(find_indices(cqs, acting_on), [dimension(s) for s in cqs.subsystems])
    push!(cqs.parametric_Hs, (ham, idxs))
end

""" Add the parametric hamiltonians from a CompositeQSystem to a Hamiltonian matrix, in place. """
function add_parametric_hamiltonians!(ham::Matrix{<:Number}, cqs::CompositeQSystem, t::Real)
    for (ham_func, idxs) = cqs.parametric_Hs
        embed_add!(ham, ham_func(t), idxs)
    end
end

""" Add a subystem collapse static operator matrix to a CompositeQSystem """
function add_lindblad!(cqs::CompositeQSystem, lind_op::AbstractMatrix{T}, acting_on::Union{Q, Array{Q}}) where {T<:Number, Q<:QSystem}
    idxs = embed_indices(cqs, acting_on)
    push!(cqs.fixed_Ls, (lind_op, idxs))
end

""" Add a subystem time dependent collapse operator matrix to a CompositeQSystem """
function add_lindblad!(cqs::CompositeQSystem, lind_op::Function, acting_on::Union{Q, Array{Q}}) where Q<:QSystem
    idxs = embed_indices(cqs, acting_on)
    push!(cqs.parametric_Ls, (lind_op, idxs))
end

""" In place addition of an operator embedded into a larger Hilbert space given a set of expansion indices"""
function embed_add!(op, added_op, expand_idxs)
    # TODO explore performance of @inbounds
    for ct = 1:length(expand_idxs)
        for idx = expand_idxs[ct]
            op[idx] += added_op[ct]
        end
    end
end

"""
    embed(m::Matrix, acting_on::Vector, dims::Vector)

Embed a subsystem operator `m` acting on subsystem indices `acting_on` into a larger tensor product
space with subsystem dimensions `dims`. This could be implemented by starting with a tensor product
ordering of `m` ⊗ I and then swapping the subsystems of `m` into the desired order. E.g. in a three
qubit system CNOT13 can be created with Id⊗SWAP * CNOT⊗Id * Id⊗SWAP.

In practice we use the fact that a tensor product structure is a statement about the strides bewteen
various subsystem elements and using reshaping the array to match the subsystem dimensions to expose
those strides and make the swaps a `permutedims`.
"""
function embed(m::Matrix, acting_on::Vector, dims::Vector)

    @assert size(m, 1) == prod(dims[acting_on]) "Oops! Dimensions of operator do not match dims argument."

    # create the large matrix by tensoring on an identity term to the operator`m`
    d_rest = prod(dims) ÷ size(m,1)
    M = kron(m, Matrix{eltype(m)}(I, d_rest, d_rest))

    # Reshape into a multi-dimensional array with the dimensions given by swapped (i.e. `m` operator
    # first) subsystem dimensions. Because the of the way tensor product indices work (last
    # subsystem is fastest) we need to reverse the order of the dimensions.
    # e.g
    # julia> reshape(["a"; "b"] ⊗ ["a"; "b"; "c"] ⊗ ["a"; "b"; "c"; "d"], (4,3,2))
    # 4×3×2 Array{String,3}:
    # [:, :, 1] =
    #  "aaa"  "aba"  "aca"
    #  "aab"  "abb"  "acb"
    #  "aac"  "abc"  "acc"
    #  "aad"  "abd"  "acd"
    #
    # [:, :, 2] =
    #  "baa"  "bba"  "bca"
    #  "bab"  "bbb"  "bcb"
    #  "bac"  "bbc"  "bcc"
    #  "bad"  "bbd"  "bcd"

    # For a matrix we repeat for rows then columns
    l = length(dims)
    identity_idxs = filter(x -> !(x in acting_on), 1:l)
    swapped_dims = reverse([dims[acting_on]; dims[identity_idxs]])
    M = reshape(M, tuple([swapped_dims; swapped_dims]...))

    # now we want to permute the dimensions to restore the correct subsystem ordering start from the
    # permutation to get the current ordering with the `m` operator first; then invert and again
    # handle the way tensor product indices work e.g. if we wanted the permutation [2 1 3] (i.e.
    # swap 1st and 2nd subsystem) then we actually need to swap the 2nd and 3rd dimensions
    reverse_perm = (l+1) .- reverse(invperm([acting_on; identity_idxs]))
    IC = [reverse_perm; reverse_perm .+ l]
    M = TensorOperations.tensorcopy(M, 1:length(IC), IC)

    # reshape back
    return reshape(M, prod(dims), prod(dims))
end

"""
    embed_indices(acting_on::Vector, dims::Vector)

Return the linear indices that map a subystem operation acting on subystem indices `acting_on` into
a larger tensor product space with subsystem dimensions `dims`.
"""
function embed_indices(acting_on::Vector, dims::Vector)
    # Strategy:
    # 1. create a subystem matrix with elements equal to the linear index e.g. [1 2; 3 4]
    dim_acting_on = prod(dims[acting_on])
    m = reshape(collect(1:dim_acting_on^2), (dim_acting_on,dim_acting_on))
    # 2. embed the matrix
    M = embed(m, acting_on, dims)
    # 3. find where the linear indices got embedded
    f(A) = (LinearIndices(A))[findall(A)]
    return [f(M .== x) for x in 1:dim_acting_on^2]
end

embed_indices(cqs::CompositeQSystem, acting_on) = embed_indices(find_indices(cqs, acting_on), [dimension(s) for s in cqs.subsystems])

""" Calculate the static part of the Hamiltonian of a CompositeQSystem """
function hamiltonian(cqs::CompositeQSystem)
    ham = zeros(ComplexF64, (dimension(cqs), dimension(cqs)))
    for (new_ham, idxs) in cqs.fixed_Hs
        embed_add!(ham, new_ham, idxs)
    end
    return ham
end


text:
using CircularArrays
using Test

v1 = CircularVector(rand(Int64, 5))

@test IndexStyle(CircularArray) == IndexCartesian()
@test IndexStyle(CircularVector) == IndexLinear()

@test size(v1, 1) == 5
@test typeof(v1) == CircularVector{Int64}
@test isa(v1, CircularVector)
@test isa(v1, AbstractVector{Int})
@test !isa(v1, AbstractVector{String})
@test v1[2] == v1[2 + length(v1)]
v1[2] = 0
v1[3] = 0
@test v1[2] == v1[3]
@test_throws MethodError v1[2] = "Hello"

v2 = CircularVector("abcde", 5)

@test prod(v2) == "abcde"^5

@test_throws MethodError push!(v1, 15)

b_arr = [2 4 6 8; 10 12 14 16; 18 20 22 24]
a1 = CircularArray(b_arr)
@test size(a1) == (3, 4)
@test a1[2, 3] == 14
a1[2, 3] = 17
@test a1[2, 3] == 17
@test !isa(a1, CircularVector)
@test !isa(a1, AbstractVector)
@test isa(a1, AbstractArray)

@test size(reshape(a1, (2, 2, 3))) == (2, 2, 3)

a2 = CircularArray(4, (2, 3))
@test isa(a2, CircularArray{Int, 2})


text:
#####
##### `reshape`
#####

function rrule(::typeof(reshape), A::AbstractArray, dims::Tuple{Vararg{Int}})
    return reshape(A, dims), (Rule(Ȳ->reshape(Ȳ, dims)), DNERule())
end

function rrule(::typeof(reshape), A::AbstractArray, dims::Int...)
    Y, (rule, _) = rrule(reshape, A, dims)
    return Y, (rule, fill(DNERule(), length(dims))...)
end

#####
##### `hcat` (🐈)
#####

function rrule(::typeof(hcat), A::AbstractArray, Bs::AbstractArray...)
    Y = hcat(A, Bs...)
    Xs = (A, Bs...)
    rules = ntuple(length(Bs) + 1) do i
        l = mapreduce(j->size(Xs[j], 2), Base.add_sum, 1:i-1; init=0)
        u = l + size(Xs[i], 2)
        dim = u > l + 1 ? (l+1:u) : u
        # NOTE: The copy here is defensive, since `selectdim` returns a view which we can
        # materialize with `copy`
        Rule(Ȳ->copy(selectdim(Ȳ, 2, dim)))
    end
    return Y, rules
end

#####
##### `vcat`
#####

function rrule(::typeof(vcat), A::AbstractArray, Bs::AbstractArray...)
    Y = vcat(A, Bs...)
    n = size(A, 1)
    ∂A = Rule(Ȳ->copy(selectdim(Ȳ, 1, 1:n)))
    ∂Bs = ntuple(length(Bs)) do i
        l = n + mapreduce(j->size(Bs[j], 1), Base.add_sum, 1:i-1; init=0)
        u = l + size(Bs[i], 1)
        Rule(Ȳ->copy(selectdim(Ȳ, 1, l+1:u)))
    end
    return Y, (∂A, ∂Bs...)
end

#####
##### `fill`
#####

function rrule(::typeof(fill), value::Any, dims::Tuple{Vararg{Int}})
    return fill(value, dims), (Rule(sum), DNERule())
end

function rrule(::typeof(fill), value::Any, dims::Int...)
    return fill(value, dims), (Rule(sum), ntuple(_->DNERule(), length(dims))...)
end


text:
#!/usr/bin/env julia

cmds = ARGS
cmd = ARGS[1]
times = try
    parse(Int64, ARGS[2])
    catch exp
    @show exp
    exit(-1)
end
for i in 1:times
    run(`$cmd`)
    sleep(2)
end


text:
import OpenCL.cl

const CL_float2 = Complex64
const CL_double2 = Complex128

clblasInvalidCommandQueue             = cl.CL_INVALID_COMMAND_QUEUE
clblasInvalidContext                  = cl.CL_INVALID_CONTEXT
clblasInvalidMemObject                = cl.CL_INVALID_MEM_OBJECT
clblasInvalidDevice                   = cl.CL_INVALID_DEVICE
clblasInvalidEventWaitList            = cl.CL_INVALID_EVENT_WAIT_LIST
clblasOutOfResources                  = cl.CL_OUT_OF_RESOURCES
clblasOutOfHostMemory                 = cl.CL_OUT_OF_HOST_MEMORY
clblasInvalidOperation                = cl.CL_INVALID_OPERATION
clblasCompilerNotAvailable            = cl.CL_COMPILER_NOT_AVAILABLE
clblasBuildProgramFailure             = cl.CL_BUILD_PROGRAM_FAILURE
clblasNotImplemented                  = -1024
clblasNotInitialized                  = -1023
clblasInvalidMatA                     = -1022
clblasInvalidMatB                     = -1021
clblasInvalidMatC                     = -1020
clblasInvalidVecX                     = -1019
clblasInvalidVecY                     = -1018
clblasInvalidDim                      = -1017
clblasInvalidLeadDimA                 = -1016
clblasInvalidLeadDimB                 = -1015
clblasInvalidLeadDimC                 = -1014
clblasInvalidIncX                     = -1013
clblasInvalidIncY                     = -1012
clblasInsufficientMemMatA             = -1011
clblasInsufficientMemMatB             = -1010
clblasInsufficientMemMatC             = -1009
clblasInsufficientMemVecX             = -1008
clblasInsufficientMemVecY             = -1007

@enum clblasOrder clblasRowMajor clblasColumnMajor
@enum clblasTranspose clblasNoTrans clblasTrans clblasConjTrans
@enum clblasUplo clblasUpper clblasLower
@enum clblasDiag clblasUnit clblasNonUnit
@enum clblasSide clblasLeft clblasRight