text:
module CBSOData3

import HTTP
import JSON
import Tables

export get_tables, get_meta, ODataTable

BASE_URL = "http://opendata.cbs.nl"
API = "ODataFeed/odata"
CATALOG = "ODataCatalog"

function get_odata(url)
    result = HTTP.get(url)
    data = JSON.parse(String(result.body))
end

"""
    get_tables(base, catalog)

Get all the tables in the catalog as a list of dicts. The `Identifier` can be
used to get the actual data. If no parameters are given, the data is retrieved
from the (CBS)[https://www.cbs.nl] OData3 portal.

The optional parameters are:

  - base::String, basename of the OData3 server
  - catalog::String, path part of the OData3 server for catalog information

"""
function get_tables(base=BASE_URL, catalog=CATALOG)
    url = base * "/" * catalog * "/Tables?\$format=json"
    return get_odata(url)["value"]
end

"""
    get_meta(table; base, api)

Get the metadata for a given `table`. The result is a dict with the keys `"TableInfos"`, `"DataProperties"` and all the classifications.

Optional parameters are:

  - base::String, basename of the OData3 server
  - api::String, path part of the OData3 server for regular data
"""
function get_meta(table; base=BASE_URL, api=API)
    url = base * "/" * api * "/" * table * "?\$format=json"
    metalinks = get_odata(url)["value"]
    meta = Dict{String,Any}()
    for link in metalinks
        name = link["name"]
        if name in ("UntypedDataSet", "TypedDataSet")
            continue
        end
        linkurl = link["url"] * "?\$format=json"
        meta[name] = get_odata(linkurl)["value"]
    end
    return meta
end

struct ODataTable
    table::String
    block::Vector{Any}
    nextlink::Union{Nothing, String}
    names::Vector{Symbol}
    types::Vector{Type}
    length::Int
end

function firstline(t::ODataTable)
    datablock = getfield(t, :block)
    nextlink = getfield(t, :nextlink)
    row = 1
    return (OdataRow(datablock[row], t), (datablock, nextlink, row))
end

function nextline(t::ODataTable, status)
    datablock, nextlink, row = status
    if row == length(datablock) && !isnothing(nextlink)
        datablock, nextlink = get_block(nextlink)
        row = 0
    end
    row += 1
    if row > length(datablock)
        return nothing
    else
        return (OdataRow(datablock[row], t), (datablock, nextlink, row))
    end
end

struct OdataRow <: Tables.AbstractRow
    row::Dict{String, Any}
    table::ODataTable
end

Tables.istable(::ODataTable) = true
names(t::ODataTable) = getfield(t, :names)
Tables.columnnames(t::ODataTable) = getfield(t, :names)
types(t::ODataTable) = getfield(t, :types)
Tables.schema(t::ODataTable) = Tables.Schema(names(t), types(t))

Tables.rowaccess(::ODataTable) = true
Tables.rows(t::ODataTable) = t
Base.eltype(t::ODataTable) = OdataRow
Base.length(t::ODataTable) = getfield(t, :length)
Base.iterate(t::ODataTable) = firstline(t)
Base.iterate(t::ODataTable, st) = nextline(t, st)

Tables.getcolumn(r::OdataRow, s::String) = something(getfield(r, :row)[s], missing)
Tables.getcolumn(r::OdataRow, i::Int) = something(collect(values(getfield(r, :row)))[getfield(getfield(r, :table), :names)[i]], missing)
Tables.getcolumn(r::OdataRow, s::Symbol) = something(getfield(r, :row)[String(s)], missing)

Tables.columnnames(r::OdataRow) = collect(keys(getfield(r, :row)))
Base.NamedTuple(r::OdataRow) = NamedTuple{Tuple(Symbol.(keys(getfield(r, :row))))}(values(getfield(r, :row)))

function get_block(url)
    data = get_odata(url)
    nextlink = "odata.nextLink" in keys(data) ? data["odata.nextLink"] : nothing
    return (data["value"], nextlink)
end

function get_table(table; base=BASE_URL, api=API, kind="UntypedDataSet", columns=[], filter="")
    if length(columns) == 0
        select = ""
    else
        select = "&\$select=" * join(columns, ",")
    end
    if length(filter) > 0
        filter = "&\$filter=" * HTTP.escape(filter)
    end
    url = base * "/" * api * "/" * table * "/" * kind * "/?\$format=json" * select * filter
    data = get_odata(url)
    nextlink = "odata.nextLink" in keys(data) ? data["odata.nextLink"] : nothing
    return (data["value"], nextlink)
end

"""
    ODataTable(table; typed, columns, filters, base, api, catalog)

Get a `Tables` object for the dataset defined by `table`.
Optional parameters are:

  - typed::Bool, do we want a TypedDataSet (numbers as ints or floats) or an
                 UntypedDataSet (everything string)
  - columns::Vector{String}, exact columnnames to select, empty for all columns
  - filters::String, OData3 specification for row filter
  - base::String, basename of the OData3 server
  - api::String, path part of the OData3 server for regular data
  - catalog::String, path part of the OData3 server for catalog information

# Examples
```julia
tbl = CBSOData3.ODataTable("82811NED", columns=["Perioden", "Onderzoekspopulatie_1", "Innovatoren_2"], filter="Perioden eq '2010X000'");
```
"""
function ODataTable(table; typed=true, columns=[], filter="", base=BASE_URL, api=API, catalog=CATALOG)
    kind = typed ? "TypedDataSet" : "UntypedDataSet"
    if length(filter) == 0
        tabellen = get_tables(base, catalog)
        tabelinfo = [x for x in tabellen if x["Identifier"] == table]
        if length(tabelinfo) == 1
            recordcount = tabelinfo[1]["RecordCount"]
        else
            recordcount = None
        end
        nextblock, nextlink = get_table(table, base=base, api=api, kind=kind, columns=columns)
    else
        # bij opgegeven filter alle blokken inlezen
        nextblock, nextlink = get_table(table, base=base, api=api, kind=kind, columns=columns, filter=filter)
        while !isnothing(nextlink)
            nxtblk, nextlink = get_block(nextlink)
            nextblock = vcat(nextblock, nxtblk)
        end
        recordcount = length(nextblock)
    end
    if length(columns) > 1
        names = Symbol.(columns)
    else
        meta = get_meta(table)
        names = Symbol.(collect(keys(nextblock[1])))
        dd = [Symbol(x["Key"])=>x["Position"] for x in meta["DataProperties"] if Symbol(x["Key"]) in names]
        sort!(dd, by=x->x.second)
        dd = [x.first=>nr for (nr, x) in enumerate(dd)]
        names = [x.first for x in dd]
    end
    types = [typeof(nextblock[1][String(name)]) for name in names]
    if typed
        types = [Union{Missing, x} for x in types]
    end
    return ODataTable(table, nextblock, nextlink, names, types, recordcount)
end

end # module


text:
using Plots, LaTeXStrings, VortexDistributions
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

using FourierGPE

## Initialize simulation
L = (200.,200.)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

## set simulation parameters
μ = 1.0
g = 0.01
γ = 0.5
ti = 0.0
tf = .3pi
Nt = 150
t = LinRange(ti,tf,Nt)

## useful state functions
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))

## make initial state
x,y = X
kx,ky = K
ψi = ψ0.(x,y',μ,g)

## ground state
ψg = copy(ψi)
showpsi(x,y,ψg)

## initial dipole
ψv = copy(ψg)
d = 30
ξv = healinglength(0.,0.,μ,g)
pv = PointVortex(0.,d/2,1)
nv = PointVortex(0,-d/2,-1)

dipole = [nv;pv]
psi = Torus(ψv,x,y) # set field topology for VortexDistributions
vortex!(psi,dipole) # make dipole
ψv = psi.ψ # pull out the field
showpsi(x,y,ψv)

## set simulation parameters
c = sqrt(μ)
tf = L[1]/c/4
γ = 0.1

t = LinRange(ti,tf,Nt)
ϕi = kspace(ψv,sim)
reltol = 1e-7 #uniform requires slightly smaller tolerance
# alg = Vern7()
simd = Sim(L,N)
@pack_Sim! simd

## remove edge artifacts with small γ
sold = runsim(simd)

## plot
ψd = xspace(sold[end],sim)
showpsi(x,y,ψd)

anim = psimovie(sold,simd)
gif(anim,"./examples/dipole_damping_small.gif",fps=25)

K2 = k2(K)

## energy totals

# all kinetic
heatmap(x,y,log10.(abs2.(fftshift(ϕ))),aspectratio=1)

## single vortex
ψv = copy(ψg)
ξv = healinglength(0.,0.,μ,g)
pv = PointVortex(0.,0.,1)

psi = Torus(ψv,x,y)
vortex!(psi,pv)
ψv = psi.ψ
showpsi(x,y,ψv)

psi = XField(copy(ψv),X,K,K2)
et,ei,ec = ikespectrum(psi,X,K,Kw)
Ei = fftshift(fft(ei))*prod(DX)

heatmap(x,y,log10.(abs2.(fftshift(ϕ))),aspectratio=1)

showenergiesk(psi)

ext: jl


text:
# =====  Configurations =====  #
"""
    ComplexGradientConfig(f,z)

Create a configuration for calculation of the gradient of a function
`f` with respect to complex array `z` by automatic differentiation.
"""
struct ComplexGradientConfig{T,V,N,D,M,C} <: ForwardDiff.AbstractConfig{N}
    rseeds::NTuple{N,Partials{M,V}}
    iseeds::NTuple{N,Partials{M,V}}
    duals::D
    cache::C
end

function ComplexGradientConfig(f::F,
                        x::AbstractArray{Complex{V}},
                        ::Chunk{N} = Chunk(x,length(x)),
                        ::T = Tag(f, Complex{V}),
                        cache = nothing) where {F,V,N,T}
    rseeds, iseeds = construct_complex_seeds(Partials{2N,V})
    duals = similar(x, ComplexDual{T,V,2N})
    return ComplexGradientConfig{T,V,N,typeof(duals),2N,typeof(cache)}(rseeds, iseeds, duals, cache)
end


text:
function tcm_estimateForest_RTS(N, assortments, transaction_counts, K, depth, perTol_EM, time_limit_EM)

	forest_left = Array{Int64,1}[]
	forest_right = Array{Int64,1}[]
	forest_product = Array{Int64,1}[]
	forest_isLeaf = Array{Bool,1}[]

	for k in 1:K
		tree_left, tree_right, tree_product, tree_isLeaf = tcm_randomCompleteTree( N, depth-1 )
		push!(forest_left, tree_left)
		push!(forest_right, tree_right)
		push!(forest_product, tree_product)
		push!(forest_isLeaf, tree_isLeaf)
	end

	initial_lambda = 1/K * ones(K)

	# time_limit_EM =

	println("Starting EM...")
	# perTol_EM = 1e-5

    lambda, v, loglik, elapsed_time = tcm_forest_EM_RTS(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, perTol_EM, initial_lambda)

    return lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time
end

function tcm_estimateForest_RTS_withInitial(N, assortments, transaction_counts, K_sample, depth, perTol_EM, time_limit_EM,
											initial_left, initial_right, initial_product, initial_isLeaf)

	forest_left = deepcopy(initial_left) #Array{Int64,1}[]
	forest_right = deepcopy(initial_right) #Array{Int64,1}[]
	forest_product = deepcopy(initial_product) #Array{Int64,1}[]
	forest_isLeaf = deepcopy(initial_isLeaf) #Array{Bool,1}[]

	for k in 1:K_sample
		tree_left, tree_right, tree_product, tree_isLeaf = tcm_randomCompleteTree( N, depth-1 )
		push!(forest_left, tree_left)
		push!(forest_right, tree_right)
		push!(forest_product, tree_product)
		push!(forest_isLeaf, tree_isLeaf)
	end

	K = length(initial_left) + K_sample;

	initial_lambda = 1/K * ones(K)

	# time_limit_EM =

	println("Starting EM...")
	# perTol_EM = 1e-5

    lambda, v, loglik, elapsed_time = tcm_forest_EM_RTS(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, perTol_EM, initial_lambda)

    return lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time
end

function tcm_forest_EM_RTS(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, perTol_EM, initial_lambda)
	lambda = copy(initial_lambda); #1/K * ones(K)
    v, v_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
    logv = log.(v)
    logv[ isinf.(logv)] = 0.0;
    loglik = dot( convert(Array{Float64,}, transaction_counts), logv)

    # @show v
    # @show loglik

	v = zeros(size(transaction_counts))

	# perTol_EM = 1e-6
    iter = 0;

	start_time = time()

	while (time() - start_time < time_limit_EM)
		iter += 1;
		h, v, v_by_class = tcm_forest_E(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
		lambda = tcm_forest_M_update_lambda(K, h, transaction_counts)

		v, v_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
		logv = log.(v)
		logv[ isinf.(logv)] = 0.0;
		new_loglik = dot( convert(Array{Float64,}, transaction_counts), logv)

		# println("\t Ranking EM: iter ", iter, " -- time elapsed = ", time() - start_time, " log likelihood = ", new_loglik)

	    if ( (new_loglik - loglik) / abs(loglik) < perTol_EM )
	    	temp = (new_loglik - loglik) / abs(loglik)
	    	# @show temp
	    	# @show new_loglik
	    	# @show loglik

	    	loglik = new_loglik
			break
		end

        # if (iter > 100)
        #     error("Stop here")
        # end

		loglik = new_loglik
	end

	elapsed_time = time() - start_time

	return lambda, v, loglik, elapsed_time

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
export convertCsvToBasketsExample, dppLearningExample,
  dppLearningBayesianExample, predictionExample,
  predictionForMCMCSamplesExample, predictionMetricsExample

dataDir = "../data/Amazon-baby-registry"

# An example of converting a file in sparse CSV basket format to the data
# format expected by the DPP modules.
function convertCsvToBasketsExample()
  csvBasketDataFileName = "$dataDir/1_100_100_100_apparel_regs.csv"
  trainingBasketsDictFileName = "$dataDir/apparel-regs-training-basketsDict.jls"
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jls"
  trainingSetSizePercent = 0.8

  convertSparseCsvToBaskets(csvBasketDataFileName, trainingBasketsDictFileName,
    testBasketsDictFileName, trainingSetSizePercent)
end

# An example of using stochastic gradient ascent for DPP learning
function dppLearningExample()
  trainingBasketsDictFileName = "$dataDir/apparel-regs-training-basketsDict.jls"
  trainingBasketsDictObjectName = "trainingBasketsDict"
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jls"
  testBasketsDictObjectName = "testBasketsDict"
  learnedModelOutputDirName = "example-SGA-model"
  numItemTraits = 30
  alpha = 1.0

  doDPPLearningSparseVectorData(trainingBasketsDictFileName,
    trainingBasketsDictObjectName, testBasketsDictFileName,
    testBasketsDictObjectName, learnedModelOutputDirName, numItemTraits, alpha)

end

# An example of using stochastic gradient HMC for DPP learning
function dppLearningBayesianExample()
  trainingBasketsDictFileName = "$dataDir/apparel-regs-training-basketsDict.jls"
  trainingBasketsDictObjectName = "trainingBasketsDict"
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jls"
  testBasketsDictObjectName = "testBasketsDict"
  learnedModelOutputDirName = "example-SGHMC-model"
  numItemTraits = 30

  doDPPBayesianLearningSparseVectorData(trainingBasketsDictFileName,
    trainingBasketsDictObjectName, testBasketsDictFileName,
    testBasketsDictObjectName, learnedModelOutputDirName, numItemTraits,
    runStochasticGradientHamiltonianMonteCarloSampler)
end

# An example of computing predictions using a DPP model learned by stochastic
# gradient ascent
function predictionExample()
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jls"
  testBasketsDictObjectName = "testBasketsDict"
  numItemTraits = 30
  alpha = 1.0
  learnedDPPParamsFileName = "example-SGA-model/learnedDPPParamsMatrix-k$numItemTraits-lambdaPop$alpha.jls"
  resultsForTestInstancesDictFileName =
    "example-SGA-model/resultsForTestInstancesDict-k$numItemTraits-lambdaPop$alpha-Amazon-apparel-regs.jls"

  computePredictionsSparseVectorData(testBasketsDictFileName,
    testBasketsDictObjectName, learnedDPPParamsFileName,
    resultsForTestInstancesDictFileName)
end

# An example of computing predictions using a DPP model learned by stochastic
# gradient HMC
function predictionForMCMCSamplesExample()
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jld"
  testBasketsDictObjectName = "testBasketsDict"
  numItemTraits = 30
  sampleLag = 1
  collectedMCMCSamplesDirPathName =
    "example-SGHMC-model/learnedDPPParams-k$numItemTraits-SGHMC-collectedSamples"
  resultsForTestInstancesDictFileName =
    "example-SGHMC-model/resultsForTestInstancesDict-k$numItemTraits-collectedSamples-SGHMC-Amazon-apparel-regs.jld"

  computePredictionsForMCMCSamples(testBasketsDictFileName,
    testBasketsDictObjectName, collectedMCMCSamplesDirPathName, sampleLag,
    resultsForTestInstancesDictFileName)
end

# An example of computing prediction metrics using DPP model prediction results
function predictionMetricsExample()
  trainingBasketsDictFileName = "$dataDir/apparel-regs-training-basketsDict.jls"
  trainingBasketsDictObjectName = "trainingBasketsDict"
  testBasketsDictFileName = "$dataDir/apparel-regs-test-basketsDict.jls"
  testBasketsDictObjectName = "testBasketsDict"
  numItemTraits = 30
  alpha = 1.0
  resultsForTestInstancesDictFileName =
    "example-SGA-model/resultsForTestInstancesDict-k$numItemTraits-lambdaPop$alpha-Amazon-apparel-regs.jls"
  # resultsForTestInstancesDictFileName =
  #   "example-SGHMC-model/resultsForTestInstancesDict-k$numItemTraits-collectedSamples-SGHMC-Amazon-apparel-regs.jls"

  computePredictionMetricsSparseVectorData(trainingBasketsDictFileName,
    trainingBasketsDictObjectName, testBasketsDictFileName,
    testBasketsDictObjectName, resultsForTestInstancesDictFileName)
end

# Optimization-based learning pipeline:
# convertCsvToBasketsExample()
# dppLearningExample()
# predictionExample()
# predictionMetricsExample()

# Bayesian learning pipeline:
# dppLearningBayesianExample()
# predictionForMCMCSamplesExample()
# predictionMetricsExample()


text:
# Is this overhead material or not?
# Maybe the benefit of using a constant is not so significant
# in real use cases?

variable = 10
function sum_variable_many_times(N)
    total = rand(variable)
    for i in 1:N
        total .+= rand(variable)
    end
    return total
end
@btime sum_variable_many_times(100);
#=
julia> @btime sum_variable_many_times(100);
  54.986 μs (301 allocations: 20.47 KiB)
=#

const constant = 10
function sum_constant_many_times(N)
    total = rand(constant)
    for i in 1:N
        total .+= rand(constant)
    end
    return total
end
@btime sum_constant_many_times(100);
#=
julia> @btime sum_constant_many_times(100);
  8.282 μs (101 allocations: 15.78 KiB)
=#

#= How bad?
julia> 54.986 / 8.282
6.6392175802946145
=#


text:
#unit-ubound.jl
#test that ubound constructors work

#test that the resolve_as_utype! directive works.

x = Unum{2,2}(0x000000000000000C, 0x0000000000000000, 0x0001, 0x0003, 0x0003)
y = Unum{2,2}(0x000000000000000C, 0x1000000000000000, 0x0001, 0x0003, 0x0003)
z = Unum{2,2}(0x000000000000000C, 0x0000000000000000, 0x0001, 0x0003, 0x0002)
@test Unums.resolve_as_utype!(x, y) == z

x = Unum{2,2}(0x000000000000000C, 0x2000000000000000, 0x0001, 0x0003, 0x0003)
y = Unum{2,2}(0x000000000000000C, 0x3000000000000000, 0x0001, 0x0003, 0x0003)
z = Unum{2,2}(0x000000000000000C, 0x2000000000000000, 0x0001, 0x0003, 0x0002)
@test Unums.resolve_as_utype!(x, y) == z

#disable this test for now.
#=
x = sss(Unum{3,5})
y = Unum{3,5}(0x0000000000000000, 0xFFFFFFFF00000000, 0x0001, 0x0007, 0x001F)
q = Unums.resolve_as_utype!(x, y)
@test q.lower == Unum{3,5}(0x0000000000000000, 0x0000000000000000, 0x0001, 0x0007, 0x0000)
@test q.upper == Unum{3,5}(0x0000000000000000, 0x8000000000000000, 0x0001, 0x0007, 0x0000)
=#

#x = Unum{3,5}(0x0000000000000000, 0x8000000000000000, 0x0001, 0x0007, 0x001F)
#y = Unum{3,5}(0x0000000000000000, 0xF000000000000000, 0x0001, 0x0007, 0x0004)
#q = Unums.resolve_as_utype!(x, y)
#@test q.lower == Unum{3,5}(0x0000000000000000, 0x8000000000000000, 0x0001, 0x0007, 0x0001)
#@test q.upper == Unum{3,5}(0x0000000000000000, 0xF000000000000000, 0x0001, 0x0007, 0x0004)

#test comparisons against similar unums

x = Unum{2,2}(0x000000000000000C, 0x2000000000000000, 0x0001, 0x0003, 0x0002)
b = Ubound(Unum{2,2}(0x000000000000000C, 0x2000000000000000, 0x0001, 0x0003, 0x0003),Unum{2,2}(0x000000000000000C, 0x3000000000000000, 0x0001, 0x0003, 0x0003))

@test x == b

################################################################################
## open-ubound-helper:  A short function which returns the open interval version
## of a bounding unum for a ubound.
#=
# first test warlpiri one as a lower bound should yield walpiri some.
@test Unums.__open_ubound_helper(one(Unum{0,0}), true) == Unum{0,0}(z16, z16, o16, t64, z64)
# as an upper bound should yield walpiri few
@test Unums.__open_ubound_helper(one(Unum{0,0}), false) == Unum{0,0}(z16, z16, o16, z64, z64)
# check to see that warlpiri few is unchanged.
@test Unums.__open_ubound_helper(pos_sss(Unum{0,0}), true) == pos_sss(Unum{0,0})
@test Unums.__open_ubound_helper(pos_sss(Unum{0,0}), false) == pos_sss(Unum{0,0})
# let's do some awkward zero tests.
@test Unums.__open_ubound_helper(zero(Unum{0,0}), true) == pos_sss(Unum{0,0})
@test Unums.__open_ubound_helper(zero(Unum{0,0}), false) == neg_sss(Unum{0,0})
@test Unums.__open_ubound_helper(-zero(Unum{0,0}), true) == pos_sss(Unum{0,0})
@test Unums.__open_ubound_helper(-zero(Unum{0,0}), false) == neg_sss(Unum{0,0})
#and finally, test wierd infinities
#@test_throws ArgumentError Unums.__open_ubound_helper(pos_inf(Unum{0,0}), true)
#@test_throws ArgumentError Unums.__open_ubound_helper(neg_inf(Unum{0,0}), false)
@test Unums.__open_ubound_helper(pos_inf(Unum{0,0}), false) == pos_mmr(Unum{0,0})
@test Unums.__open_ubound_helper(neg_inf(Unum{0,0}), true) == neg_mmr(Unum{0,0})

################################################################################
## open-ubound:  Do we generate open ubounds correctly?

#first on warlpiri ubounds, are we correctly generating open ubounds?  Start
#doubly-closed end input.
uresult = Ubound(neg_sss(Unum{0,0}), pos_sss(Unum{0,0}))
@test open_ubound(zero(Unum{0,0}), one(Unum{0,0})) == pos_sss(Unum{0,0})
@test open_ubound(neg_one(Unum{0,0}), one(Unum{0,0})) == uresult
#half-open inputs
@test open_ubound(neg_sss(Unum{0,0}), one(Unum{0,0})) == uresult
@test open_ubound(neg_one(Unum{0,0}), pos_sss(Unum{0,0})) == uresult
#open inputs
@test open_ubound(neg_sss(Unum{0,0}), pos_sss(Unum{0,0})) == uresult
#in unum{1,1} test oddly overlapping ubounds. This is a ubound pair that looks
#as follows:
# (a)
# ( b )
# and should appropriately generate b as the ubound.
olu_b = Unum{1,1}(z16, z16, o16, z64, o64)
olu_a = Unum{1,1}(o16, z16, o16, z64, o64)
println(olu_b)
println(olu_a)
#println(open_ubound(olu_a, olu_b))
println(Unums.ubound_resolve(Ubound(olu_a, olu_b)))
@test open_ubound(olu_a, olu_b) == olu_b
=#

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using HttpServer
using WebSockets

#global Dict to store open connections in
global connections = Dict{Int,WebSocket}()
global usernames   = Dict{Int,String}()

function decodeMessage( msg )
    bytestring(msg)
end

wsh = WebSocketHandler() do req, client
    global connections
    @show connections[client.id] = client
    while true
        msg = read(client)
        msg = decodeMessage(msg)
        if beginswith(msg, "setusername:")
            println("SETTING USERNAME: $msg")
            usernames[client.id] = msg[13:]
        end
        if beginswith(msg, "say:")
            println("EMITTING MESSAGE: $msg")
            for (k,v) in connections
                if k != client.id
                    write(v, (usernames[client.id] * ": " * msg[5:]))
                end
            end
        end
    end
end

onepage = readall("./examples/chat-client.html")
httph = HttpHandler() do req::Request, res::Response
  Response(onepage)
end

server = Server(httph, wsh)
println("Chat server listening on 8000...")
run(server,8000)


text:
# Julia wrapper for header: /opt/local/include/gdal_version.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0


text:
# volatility functions
@testset "Volatility" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = cumsum(randn(N))
        X = cumsum(randn(N, 3), dims=1)
        tmp = bbands(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = tr(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = atr(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = keltner(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
    end
    @testset "Temporal" begin
        x = TS(cumsum(randn(N)))
        X = TS(cumsum(randn(N, 3), dims=1))
        tmp = bbands(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = tr(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = atr(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = keltner(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
    end
end


text:
"""
`module  JAC.RadiativeOpacity`
    ... a submodel of JAC that contains all methods for computing radiative -- spectral & mean -- opacities.
"""
module RadiativeOpacity

    using Printf, JAC, ..Basics, ..Defaults, ..ManyElectron, ..Radial, ..TableStrings

    """
    `struct  RadiativeOpacity.Settings  <:  AbstractPropertySettings`
        ... defines a type for the details and parameters of computing radiative opacities.

        + approach           ::AbstractOpacityApproach
            ... Determines the approach for calculating radiative opacities: {SincleConf(), ...}.
        + NoExcitations      ::Int64
            ... Number of excitations with regard to the given (reference) configurations.
        + NoLevels           ::Int64
            ... maximum number of levels (to terminate long computations in a regular fashion).
        + maxExcitation      ::Float64
            ... maximum excitation energy w.r.t. the ground level for which levels are considered.
        + deltaLambda        ::Float64          ... Delta Lambda in the expansion opacity.
        + time               ::Float64          ... Time after merger in the expansion opacity.
        + densities          ::Array{Float64,1} ... densities [g/cm^3] in the expansion opacity.
        + temperatures       ::Array{Float64,1} ... temperatures [K] in the expansion opacity.
    """
    struct Settings  <:  AbstractPropertySettings
        approach             ::Int64
        NoExcitations        ::Int64
        NoLevels             ::Int64
        maxExcitation        ::Float64
        DeltaLambda          ::Float64
        time                 ::Float64
        densities            ::Array{Float64,1}
        temperatures         ::Array{Float64,1}
    end

    """
    `RadiativeOpacity.Settings()`
        ... constructor for an `empty` instance of RadiativeOpacity.Settings for the computation of radiative opacities.
    """
    function Settings()
        Settings( )
    end

    # `Base.show(io::IO, settings::RadiativeOpacity.Settings)`  ... prepares a proper printout of settings::RadiativeOpacity.Settings.
    function Base.show(io::IO, settings::RadiativeOpacity.Settings)
        println(io, "approach:              $(settings.approach)  ")
        println(io, "NoExcitations:         $(settings.NoExcitations)  ")
        println(io, "NoLevels:              $(settings.NoLevels)  ")
        println(io, "maxExcitation:         $(settings.maxExcitation)  ")
        println(io, "DeltaLambda:           $(settings.DeltaLambda)  ")
        println(io, "time:                  $(settings.time)  ")
        println(io, "densities:             $(settings.densities)  ")
        println(io, "temperatures:          $(settings.temperatures)  ")
    end

    """
    `struct  RadiativeOpacity.Outcome`
        ... defines a type to keep the outcome of a opacity computation as well other results.

        + spectralOpacity       ::Float64            ... Spectral opacity.
        + expansionOpacity      ::Float64            ... Expansion opacity.
        + rosselandOpacity      ::Float64            ... Mean Rosseland opacity.
    """
    struct Outcome
        spectralOpacity         ::Float64
        expansionOpacity        ::Float64
        rosselandOpacity        ::Float64
    end

    """
    `RadiativeOpacity.Outcome()`
        ... constructor for an `empty` instance of RadiativeOpacity.Outcome for the computation of opacities.
    """
    function Outcome()
        Outcome(0., 0., 0.)
    end

    # `Base.show(io::IO, outcome::RadiativeOpacity.Outcome)`  ... prepares a proper printout of outcome::RadiativeOpacity.Outcome.
    function Base.show(io::IO, outcome::RadiativeOpacity.Outcome)
        println(io, "spectralOpacity:         $(outcome.spectralOpacity)  ")
        println(io, "expansionOpacity:        $(outcome.expansionOpacity)  ")
        println(io, "rosselandOpacity:        $(outcome.rosselandOpacity)  ")
    end

    ## Introduce some sort Opacity.Line to keep the data E_i, E_f, level_i, level_f, J_i, J_f, gf ... for all calculated
    ## lines ... but not the full representation of amplitudes and levels
    ## Perhaps, levels::Array{Level,1}  ... to keep a number of levels
    ## calculate multiplet of configurations
    ## generate a list of configurations for which levels are taken into account.

end # module


text:
#' # Weaved Julia Script

#' This page was generated from a **julia script**, automatically converted to markdown with [Weave.jl](https://github.com/JunoLab/Weave.jl). The conversion happens automatically when the menu is constructed.

#' ## Julia Dots

#' PNG with path relative to site, i.e. `/assets/images/julia-logo-dots-small.png`, but this only works in the Franklin-rendered site:

#' ![Julia dots](/assets/images/julia-logo-dots-small.png)

#' ## Julia Speeder

#' GIF with path relative to file but within `_assets/`, i.e. `../../_assets/images/juliaspeeder32x32.gif`, and postprocessed accordingly, so it works both in Franklin and in Weave:

#' ![Julia speeder](../../_assets/images/juliaspeeder32x32.gif)

#' ## Math

#' Inline, as in $e^{i\pi} = -1$, or in a block:
#' $$
#' \exp(i\pi) + 1 = 0
#' $$

#' ## Code chunks

#' Here is a julia chunk:

x = 1

#' and another:

println("hello world!")

#' ## Load file

#' Keep in mind that Weave renders this from the directory where the file will be weaved to, which is different from the way Literate handles it (*c.f.* Section [literatejl.jl](/pages/literated/literatejl)):

pwd()

#' With that settled, we can use the proper path:

open("../../_assets/data/wow.txt") do io
    read(io, String)
end

#' ## Plot

using Plots

x = 0.0:0.01:2π
y = sin.(2x) + sin.(5x)

plot(x, y, title = "A nice wave", titlefont=12, label = false)

max_stars_repo_name: rmsrosa/booksjl-franklin-template

max_issues_repo_name: rmsrosa/booksjl-franklin-template

max_forks_repo_name: rmsrosa/booksjl-franklin-template


text:
################################################################################
#
#   DEFINITION OF ABSTRACT HAMILTONIAN
#   - label type L (of bond labels)
#   - unitcell type U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
#   - bond hamiltonian type HB<:AbstractBondHamiltonian{L,NS} where {NS}
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractHamiltonian{
        L,
        UC<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
        HB<:AbstractBondHamiltonian{L,NS} where {NS}
    } end

# export the type
export AbstractHamiltonian

# DOCSTRING
"""
    abstract type AbstractHamiltonian{L,UC,HB}

The abstract Hamiltonian type that describes all quadratic hamiltonian implementations in reciprocal space. It is parametric in three types,
- `L` the label type of the bonds it covers
- `UC` the unitcell that it describes
- `HB` the bond hamiltonian that describes interactions on the bonds

# Examples

... not implemented yet ...

"""
AbstractHamiltonian

################################################################################
#
#	INTERFACING / ACCESSING HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# Function to obtain the unitcell
import LatPhysBase.unitcell
function unitcell(
            h :: H
        ) :: UC where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'unitcell' for Hamiltonian type " * string(typeof(h)))
end

# export the function
export unitcell

# Function to obtain the bond hamiltonian
function bondHamiltonian(
            h :: H
        ) :: HB where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'bondHamiltonian' for Hamiltonian type " * string(typeof(h)))
end

# export the function
export bondHamiltonian

function saveHamiltonian(
        h  :: H,
        fn :: AbstractString,
        group :: AbstractString = "hamiltonian"
        ;
        append :: Bool = false
    ) where {L,UC<:AbstractUnitcell,N,HB<:AbstractBondHamiltonian{L,N}, H<:AbstractHamiltonian{L,UC,HB}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'saveHamiltonian' for Hamiltonian type " * string(H))
end

function loadHamiltonian(
        ::Type{H},
        fn :: AbstractString,
        group :: AbstractString = "hamiltonian"
    ) where {L,UC<:AbstractUnitcell,N,HB<:AbstractBondHamiltonian{L,N}, H<:AbstractHamiltonian{L,UC,HB}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'loadHamiltonian' for Hamiltonian type " * string(H))
end

export loadHamiltonian, saveHamiltonian

################################################################################
#
#	MATRIX CALCULATION
#
################################################################################

# function to compute the matrix for any bond term dimension
function matrixAtK(
            h :: H,
            k :: Vector{<:Real}
        ) :: Matrix{Complex{Float64}} where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # new matrix
    h_matrix = zeros(Complex{Float64}, dim(h),dim(h))

    # add all bonds to the matrix
    @inbounds @simd for b in bonds(unitcell(h))
        # get the bond indices
        i_from  = from(b) :: Int64
        i_to    = to(b)   :: Int64
        delta_r = vector(b, unitcell(h))
        # get the interaction matrix and add it to the general matrix
        h_matrix[(i_from-1)*(NS)+1:(i_from)*(NS), (i_to-1)*(NS)+1:(i_to)*(NS)] .+= 0.5 .* bondterm(bondHamiltonian(h), b) .* exp(im*dot(k,delta_r))
    end

    # return the matrix
    return h_matrix
end

# function to compute the matrix for bond term dimension 1
function matrixAtK(
            h :: H,
            k :: Vector{<:Real}
        ) :: Matrix{Complex{Float64}} where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,1},H<:AbstractHamiltonian{L,UC,HB}}

    # new matrix
    h_matrix = zeros(Complex{Float64}, dim(h),dim(h))

    # scalar products
    k_dot_dx = zeros(Float64, length(bonds(unitcell(h))))
    # fill in the dot products
    @inbounds @simd for bi in 1:length(bonds(unitcell(h)))
        k_dot_dx[bi] = vector_dot_k(bonds(unitcell(h))[bi], unitcell(h), k)
    end

    # add all bonds to the matrix
    @inbounds @simd for bi in 1:length(bonds(unitcell(h)))
        # get the interaction matrix and add it to the general matrix
        h_matrix[from(bonds(unitcell(h))[bi]), to(bonds(unitcell(h))[bi])] += exp(im*k_dot_dx[bi]) * 0.5 * bondterm(bondHamiltonian(h), bonds(unitcell(h))[bi])
    end

    # return the matrix
    return h_matrix
end

# function to compute the matrix for any bond term dimension into an already allocated matrix
function matrixAtK!(
            h :: H,
            k :: Vector{<:Real},
            h_matrix :: Array{Complex{Float64}, 2}
        ) :: Matrix{Complex{Float64}} where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # check matrix
    @assert size(h_matrix) == (dim(h),dim(h)) "Wrong matrix size passend: size(h_matrix)=$(size(h_matrix)), required=$((dim(h),dim(h)))"

    # erase elements in matrix
    h_matrix .= 0.0 + 0.0im

    # add all bonds to the matrix
    @inbounds @simd for b in bonds(unitcell(h))
        # get the bond indices
        i_from  = from(b) :: Int64
        i_to    = to(b)   :: Int64
        delta_r = vector(b, unitcell(h))
        # get the interaction matrix and add it to the general matrix
        h_matrix[(i_from-1)*(NS)+1:(i_from)*(NS), (i_to-1)*(NS)+1:(i_to)*(NS)] .+= 0.5 .* bondterm(bondHamiltonian(h), b) .* exp(im*dot(k,delta_r))
    end

    # return the matrix
    return h_matrix
end

# function to compute the matrix for bond term dimension 1 into an already allocated matrix
function matrixAtK!(
            h :: H,
            k :: Vector{<:Real},
            h_matrix :: Array{Complex{Float64}, 2}
        ) :: Matrix{Complex{Float64}} where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,1},H<:AbstractHamiltonian{L,UC,HB}}

    # check matrix
    @assert size(h_matrix) == (dim(h),dim(h)) "Wrong matrix size passend: size(h_matrix)=$(size(h_matrix)), required=$((dim(h),dim(h)))"

    # erase elements in matrix
    h_matrix .= 0.0 + 0.0im

    # add all bonds to the matrix
    @inbounds @simd for b in bonds(unitcell(h))
        # get the interaction matrix and add it to the general matrix
        h_matrix[from(b), to(b)] += 0.5 * bondterm(bondHamiltonian(h), b) * exp(im*dot(k, vector(b, unitcell(h))))
    end

    # return the matrix
    return h_matrix
end

# export matrix function
export matrixAtK, matrixAtK!

# function to find a k where a certain energy is hit
function findKAtEnergy(
            h :: H,
            energy_cut :: Real,
            k_start :: Vector{<:Real}
            ;
            max_variation ::Real = 20*pi,
            max_errors :: Integer = 10,
            max_steps  :: Integer = 100,
            precision_energy :: Real = 1e-6,
            diff_step_size_k :: Real = 1e-6,
        ) :: Vector{Float64} where {L,LS,D,S<:AbstractSite{LS,D},N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # search until there are enough points
    for __try__ in 1:max_errors
        # assume k to be k_start
        k = k_start .+ Float64[2*(rand()-0.5)*max_variation*(__try__-1)/(max_errors-1) for i in 1:length(k_start)]
        H_eps = zeros(D)
        # start with the initial energy
        evs = eigvals(matrixAtK(h,k))
        e0  = minimum(abs.(evs.-energy_cut))
        # iterate i over 100 newton steps (maximum)
        for i in 1:max_steps
            # check if the energy is already converged
            if e0 < precision_energy
                #println("$(__try__-1) failed attempts, $(i) steps needed in attempt $(__try__)")
                # just return the k vector
                return k
            end
            # the current energy
            H_0 = e0
            # the energy when going along primary directions
            @simd for j in 1:D
                k_grad = deepcopy(k)
                @inbounds k_grad[j] += diff_step_size_k
                evs = eigvals(matrixAtK(h,k_grad))
                @inbounds H_eps[j] = minimum(abs.(evs.-energy_cut))
            end
            # the gradient of the energy
            dH = (H_eps .- H_0) ./ diff_step_size_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH .* (H_0 / dHdH)
            k .-= dk
            # calculate a new energy
            evs = eigvals(matrixAtK(h,k))
            e0  = minimum(abs.(evs.-energy_cut))
        end
    end

    # raise error
    error("Too many errors occured consecutively without finding points --> maybe no points with energy = $(energy_cut)?")
end

# export finding function
export findKAtEnergy

# dimension calculation
function dim(
            h :: H
        ) :: Int64 where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS},H<:AbstractHamiltonian{L,UC,HB}}

    # return the size of the hamiltonian
    return NS*numSites(unitcell(h))
end

# export dimension function
export dim


text:
import Base.+, Base.-, Base./, Base.*
import SchumakerSpline.evaluate
function evaluate(f::PE_Function, x::Real)
    diff = x - f.base_
    return f.a_ * exp(f.b_ * diff) * (diff)^f.d_
end
function (s::PE_Function)(x::Union{Real,Date,DateTime,DatePeriod})
    return evaluate(s, x)
end

function derivative(f::PE_Function)
    if f.d_ == 0
        return PE_Function(f.a_ * f.b_, f.b_, f.base_, f.d_)
    else
        return PE_Function(f.a_ * f.d_, f.b_, f.base_, f.d_ - 1) +  PE_Function(f.a_ * f.b_, f.b_, f.base_, f.d_)
    end
end

function indefinite_integral(f::PE_Function)
    if (abs(f.b_) < tol) & (f.d_ == 0)
        return PE_Function(f.a_, 0.0, 0.0, 1)
    elseif f.d_ == 0
        return PE_Function(f.a_/f.b_, f.b_, f.base_, 0)
    elseif abs(f.b_) < tol
        return PE_Function(f.a_/(1.0+f.d_), 0.0, f.base_, f.d_ + 1)
    end
    # Integration by parts:      $\int u v' = uv - \int u'v$
    # I use u = (x-c)^d and v' = ae^(b(x-c))
    v = PE_Function(f.a_/f.b_, f.b_, f.base_, 0)
    u = PE_Function(1.0,0.0,f.base_,f.d_)
    return u * v - indefinite_integral( derivative(u) *v    )
end

function +(f::PE_Function,number::Real)
    constant_function = PE_Function(number, 0.0,0.0,0)
    return Sum_Of_Functions([f, constant_function])
end
function -(f::PE_Function, number::Real)
    return +(f, -number)
end
function *(f::PE_Function, number::Real)
    return PE_Function(f.a_*number, f.b_, f.base_, f.d_)
end
function /(f::PE_Function, number::Real)
    return *(f, 1/number )
end

function +(f1::PE_Function, f2::PE_Function)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::PE_Function, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::PE_Function, f2::Piecewise_Function)
    return Piecewise_Function(f2.starts_, f1 .+ f2.functions_)
end

function +(f1::Sum_Of_Functions, f2::PE_Function)
    return +(f2,f1)
end
function +(f1::Piecewise_Function, f2::PE_Function)
    return +(f2,f1)
end

function -(f1::PE_Function, f2::PE_Function)
    return Sum_Of_Functions([f1,-1*f2])
end
function -(f1::PE_Function, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,-1*f2])
end
function -(f1::PE_Function, f2::Piecewise_Function)
    return Piecewise_Function(f2.starts_, f1 .- f2.functions_)
end

function -(f1::Sum_Of_Functions, f2::PE_Function)
    return +(f1,1*f2)
end
function -(f1::Piecewise_Function, f2::PE_Function)
    return +(f1,-1*f2)
end

function *(f1::PE_Function,f2::PE_Function)
    base_func = f1
    other_func = f2
    diff = f2.base_ - f1.base_
    if abs(diff) < tol
        return PE_Function(f1.a_ * f2.a_, f1.b_ + f2.b_, f1.base_, f1.d_ + f2.d_ )
    elseif (f1.d_ < tol) & (f1.b_ < tol)
        return f1.a_ * f2
    elseif (f2.d_ < tol) & (f2.b_ < tol)
        return f2.a_ * f1
    end
    if diff < 0
        base_func = f2
        other_func = f1
    end
    converted_other = change_base_of_PE_Function(other_func, base_func.base_)
    return base_func * converted_other
end
function *(f1::PE_Function, f2::Sum_Of_Functions)
    multiplied_functions = f1 .* f2.functions_
    return Sum_Of_Functions(multiplied_functions)
end
function *(f1::PE_Function, f2::Piecewise_Function)
    return Piecewise_Function(f2.starts_, f1 .* f2.functions_)
end

function *(f1::Sum_Of_Functions, f2::PE_Function)
    return *(f2,f1)
end
function *(f1::Piecewise_Function, f2::PE_Function)
    return *(f2,f1)
end


text:
"""
    GPAR(fs)

"""
struct GPAR{Tfs} <: AbstractGPs.AbstractGP # this is a bit of a hack, as a GPAR isn't a GP.
    fs::Tfs
end

dim_out(f::GPAR) = length(f.fs)

(f::GPAR)(x::ColVecs, Σs) = AbstractGPs.FiniteGP(f, x, Σs)

const FiniteGPAR = AbstractGPs.FiniteGP{<:GPAR}

extract_data(fx::FiniteGPAR) = fx.f, fx.x.X, fx.Σy

function AbstractGPs.rand(rng::AbstractRNG, fx::FiniteGPAR)
    f, X, Σs = extract_data(fx)
    Y = Matrix{Float64}(undef, 0, length(fx.x))

    for p in 1:dim_out(f)
        x_p = ColVecs(vcat(X, Y))
        y_p = rand(rng, f.fs[p](x_p, Σs[p]))
        Y = vcat(Y, y_p')
    end
    return ColVecs(Y)
end

function AbstractGPs.logpdf(fx::FiniteGPAR, y::ColVecs)
    f, X, Σs = extract_data(fx)
    Y = y.X

    l = logpdf(f.fs[1](ColVecs(X), Σs[1]), Y[1, :])
    for p in 2:dim_out(f)
        x_p = ColVecs(vcat(X, Y[1:p-1, :]))
        l += logpdf(f.fs[p](x_p, Σs[p]), Y[p, :])
    end
    return l
end

function AbstractGPs.posterior(fx::FiniteGPAR, y::ColVecs)
    f, X, Σs = extract_data(fx)
    Y = y.X
    fs_post = map(enumerate(f.fs)) do (p, f_p)
        x_p = ColVecs(vcat(X, Y[1:p-1, :]))
        return posterior(f_p(x_p, Σs[p]), Y[p, :])
    end
    return GPAR(fs_post)
end


text:
#####
##### Generate a flame graph for self-play using ProfileView
#####

import AlphaZero

ENV["TRAINING_MODE"] = "debug"

include("../games.jl")
GAME = get(ENV, "GAME", "connect-four")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

network = Training.Network{Game}(Training.netparams)
env = AlphaZero.Env{Game}(Training.params, network)

using Profile
using ProfileView

# The neural network should not be used at iteration 0
AlphaZero.self_play!(env, nothing) # To compile every function
Profile.clear()
@profile AlphaZero.self_play!(env, nothing)
ProfileView.svgwrite("self-play-profile.svg")


text:
FitFunction(s::Symbol) = FitFunction(Float64, Val(s))
FitFunction(T::Type{<:AbstractFloat}, s::Symbol) = FitFunction(T, Val(s))

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:Gauss})
    ff = FitFunction{T}( Gauss, 1, 3 )
    set_parameter_names!(ff, ["A", "σ", "μ"])
    set_initial_parameters!(ff, [1, 1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end
function FitFunction(T::Type{<:AbstractFloat}, s::Val{:PDF_Gauss})
    ff = FitFunction{T}( PDF_Gauss, 1, 2 )
    set_parameter_names!(ff, ["σ", "μ"])
    set_initial_parameters!(ff, [1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:GaussPlusConstantBackground})
    ff = FitFunction{T}( Gauss_plus_const_background, 1, 4 )
    set_parameter_names!(ff, ["A", "σ", "μ", "offset"])
    set_initial_parameters!(ff, [1, 1, 0, 1])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end
function FitFunction(T::Type{<:AbstractFloat}, s::Val{:GaussPlusLinearBackground})
    ff = FitFunction{T}( Gauss_plus_linear_background, 1, 5 )
    set_parameter_names!(ff, ["A", "σ", "μ", "offset", "linear_slope"])
    set_initial_parameters!(ff, [1, 1, 0, 1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end
function FitFunction(T::Type{<:AbstractFloat}, s::Val{:GaussPlusLinearBackgroundPlusConstBackground})
    ff = FitFunction{T}( Gauss_plus_linear_background_plus_const_background, 1, 6 )
    set_parameter_names!(ff, ["A", "σ", "μ", "offset", "linear_slope", "const_background"])
    set_initial_parameters!(ff, Float64[1, 1, 0, 1, 0, nextfloat(0.0)])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:Cauchy})
    ff = FitFunction{T}( Cauchy, 1, 3 )
    set_parameter_names!(ff, ["scale", "σ", "μ"])
    set_initial_parameters!(ff, [1, 1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:Linear})
    ff = FitFunction{T}( Linear, 1, 2 )
    set_parameter_names!(ff, ["offset", "slope"])
    set_initial_parameters!(ff, [0, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function min_0(x::T) where {T <: Real}
    return x < 0 ? zero(typeof(x)) : x
end
function min_0(x)
    min_0.(x)
end

function get_inf(x::T) where {T <: Real}
    return T(Inf)
end
function get_inf(x)
    return get_inf.(x)
end
function Cauchy(x, p)
    scale = p[1]
    σ = p[2]
    μ  = p[3]
    return @fastmath @. scale * σ / (π * (σ^2 + (x - μ)^2))
end

function Linear(x, p)
    offset  = p[1]
    slope   = p[2]
    return @fastmath @. offset + slope * x
end

function Gauss(x, p)
    scale = p[1]
    σ     = p[2]
    μ     = p[3]
    return @fastmath @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2))
end

function PDF_Gauss(x, p)
    σ = p[1]
    μ = p[2]
    return min_0(@fastmath @. exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)))
end

function Gauss_plus_linear_background(x, p)
    A = p[1]
    σ = p[2]
    μ = p[3]
    offset = p[4]
    lin_slope = p[5]
    if σ <= 0 || A < 0 return get_inf(x) end
    return @fastmath @. A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2) + min_0(lin_slope * (x - μ)) + offset
end

function Gauss_plus_linear_background_plus_const_background(x, p)
    A = p[1]
    σ = p[2]
    μ = p[3]
    offset = p[4]
    lin_slope = p[5]
    const_background = p[6]
    if σ <= 0 || A < 0 return get_inf(x) end
    return @. min_0(@fastmath @. A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2) + lin_slope * (x - μ) + offset) + const_background
end

function Gauss_plus_const_background(x, p)
    A = p[1]
    σ = p[2]
    μ = p[3]
    offset = p[4]
    if σ <= 0 || A < 0 return get_inf(x) end
    return min_0(@fastmath @. A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2) + offset)
end


text:
"""
    GPIO(i::Int32)
Lowest form of communication with the GPIO pins. The available pins are
listed in the "channel" parameter, and appear as directories in /sys/class/gpio
after being exported. For instance, to setup a GPIO on "gpio112", configure it as an output pin and set
it to high, the following code would be used.

    `gpio = GPIO(1)`
    `write!(gpio, (2,"out"))`
    `write!(gpio, (1, "1"))`

The operation of reading the current output value of the GPIO is done by

    `read(gpio, 1)`

See the test/BeagleBone/GPIO_test.jl for more examples.
"""
type GPIO <: IO_Object
  i::Int32
  basedir::String
  filestreams::Array{IOStream,1}
  function GPIO(i::Int32)
    (i <= 0 || i > length(gpio_channels)) && error("Invalid GPIO index: $i")

    # If the tests re being run, a dummy filesystem is created
    basedir = export_gpio(i)

    # setup IOstreams
    value_filestream = open("$(basedir)/$(gpio_channels[i])/value","r+")
    direction_filestream = open("$(basedir)/$(gpio_channels[i])/direction","r+")
    edge_filestream = open("$(basedir)/$(gpio_channels[i])/edge","r")

    # Initialize object
    return new(i, basedir, [value_filestream, direction_filestream, edge_filestream])
  end
end

"""
    write!(gpio::GPIO, args::Tuple{Int32,String}, debug::Bool=false)
Writes an entry to an operation on a GPIO, of the form args = (operation, entry).
"""
function write!(gpio::GPIO, args::Tuple{Int32,String}, debug::Bool=false)
  debug && return
  operation, entry = args[1], args[2]
  # Only filestreams 1 and 2 are writable
  operation ∉ [1,2] && error("Invalid GPIO operation $operation for writing")
  if entry in gpio_operations[operation]
    seekstart(gpio.filestreams[operation])
    write(gpio.filestreams[operation], "$entry\n")
    flush(gpio.filestreams[operation])
  else
    error("Invalid entry for GPIO operation $(operation): $(entry)")
  end
end

"""
    l = read(gpio::GPIO, operation::Int32, debug::Bool=false)
Reads the current value from an operation on a GPIO.
"""
function read(gpio::GPIO, operation::Int32, debug::Bool=false)
  debug && return
  # Filestreams 1, 2 and 3 are readable
  operation ∉ [1,2,3] && error("Invalid GPIO operation: $operation for reading")
  seekstart(gpio.filestreams[operation])
  l = readline(gpio.filestreams[operation])
  return l
end

"""
    teardown(gpio::GPIO, debug::Bool=false)
Closes all open streams on the GPIO, and unexports it from the file system.
"""
function teardown(gpio::GPIO, debug::Bool=false)
  debug && return

  #Close all IOStreams
  for stream in gpio.filestreams
    close(stream)
  end

  #Unexport filestructure
  if isdefined(:RUNNING_TESTS)
    # Remove the dummy file system for testing
    basedir = "$(pwd())/testfilesystem/gpio"
    try
      rm("$(gpio.basedir)/$(gpio_channels[gpio.i])"; recursive=true)
    catch
      error("Could not remove the requested GPIO testfiles for channel $(gpio_channels[i]).")
    end
  else
    # Remove the file system
    filename = "/sys/class/gpio/unexport"
    #TODO Verify if this is the correct command to send to unexport
    write(filename, gpio_channels[gpio.i])
  end
end

"""
    export_gpio(i::Int32, debug::Bool=false)
Export the GPIO file system, either for real-time or testing usecases.
"""
function export_gpio(i::Int32)
  if isdefined(:RUNNING_TESTS)
    # Export a dummy file system for testing
    basedir = "$(pwd())/testfilesystem/gpio"
    try
      #println("$(basedir)/$(gpio_channels[i])")
      mkpath("$(basedir)/$(gpio_channels[i])")
    catch
      error("Could not export the GPIO device for channel $(gpio_channels[i]) for testing as the directory $(basedir)/$(gpio_channels[i]) already exists.")
    end
    try
      f = open("$(basedir)/$(gpio_channels[i])/value", "w"); write(f,"0"); close(f);
      f = open("$(basedir)/$(gpio_channels[i])/direction", "w"); write(f,"0"); close(f);
      f = open("$(basedir)/$(gpio_channels[i])/edge", "w"); write(f,"0"); close(f);
    catch
      error("Could not open the requested GPIO testfiles for channel $(gpio_channels[i]).")
    end
  else
    basedir = "/sys/class/gpio"
    # TODO Export for real-time use
  end
  return basedir
end

"""
    to_string(gpio::GPIO, debug::Bool=false)

Generates a string representation of the GPIO device.
"""
function to_string(gpio::GPIO, debug::Bool=false)
  debug && return
  message = "\nID: $(gpio.i)\n\nAvailable filestreams:\n"
  for ii = 1:length(gpio.filestreams)
    message = string(message, "  index=$(ii) - name=$(gpio.filestreams[ii].name) - write/read=$(iswritable(gpio.filestreams[ii]))/$(isreadable(gpio.filestreams[ii]))\n")
  end
  return message
end


text:
function multitime(callable; numruns=20, seconds=2)
    benchmark_f = () -> @elapsed callable()
    runtime = benchmark_f()
    if runtime < seconds
        runtime = mapreduce(j -> benchmark_f(), min, 2:numruns; init = runtime)
    end
    return runtime
end

function chi2(gaussian_estimate, actual_value)
    diff = gaussian_estimate.μ - actual_value
    chi2_vals = diff' * (gaussian_estimate.Σ \ diff)
    return chi2_vals
end

function make_savable(sol; appxsol=nothing)
    # sol = solve(args...; kwargs...)
    @assert DiffEqBase.has_analytic(sol.prob.f) || !isnothing(appxsol)

    y_ref(t) =
        DiffEqBase.has_analytic(sol.prob.f) ?
        sol.prob.f.analytic(sol.prob.u0, sol.prob.p, t) :
        appxsol(t)
    densetimes = collect(range(sol.t[1], stop=sol.t[end], length=1000))
    return (
        x=sol.x,
        t=sol.t,
        pu=sol.pu,
        u_ref=y_ref.(sol.t),
        t_dense=densetimes,
        pu_dense=sol.p(densetimes),
        u_ref_dense=y_ref.(densetimes),
    )
end

function compute_errors(savable_sol)
    # savable_sol should be an output of `make_savable`
    ssol = savable_sol
    diffs = ssol.pu.μ - ssol.u_ref
    dense_diffs = ssol.pu_dense.μ - ssol.u_ref_dense

    errors = Dict{Symbol, Float64}()

    errors[:final] = norm(mean(abs.(diffs[end])))
    errors[:final_est] = norm(mean(abs.(sqrt.(diag(ssol.pu.Σ[end])))))

    errors[:l2] = norm(sqrt.(mean.([float.(d) .^ 2 for d in diffs])))
    errors[:l2_est] = norm(sqrt.(mean.(diag.(ssol.pu.Σ))))
    errors[:l∞] = norm(maximum.([float.(d) .^ 2 for d in diffs]))
    errors[:l∞_est] = norm(maximum.(diag.(ssol.pu.Σ)))

    errors[:L2] = norm(sqrt.(mean.([float.(d) .^ 2 for d in dense_diffs])))
    errors[:L2_est] = norm(sqrt.(mean.(diag.(ssol.pu_dense.Σ))))
    errors[:L∞] = norm(maximum.([float.(d) .^ 2 for d in dense_diffs]))
    errors[:L∞_est] = norm(maximum.(diag.(ssol.pu_dense.Σ)))

    errors[:final_χ²] = chi2(ssol.pu[end], ssol.u_ref[end])
    errors[:χ²] = mean(remove_initial_nan(chi2.(ssol.pu, ssol.u_ref)))
    errors[:Χ²] = mean(remove_initial_nan(chi2.(ssol.pu_dense, ssol.u_ref_dense)))

    return errors
end

function remove_initial_nan(arr)
    @assert issorted(.!isnan.(arr))
    return arr[.!isnan.(arr)]
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
using Pkg
Pkg.activate("test")

using Clustering
using DataFrames
using Gurobi
using JuMP
using CSV
using MathOptInterface
using MLDatasets
using Test
using Random

Random.seed!(1);
include("../src/OptimalTrees.jl")
using .OptimalTrees
global OT = OptimalTrees
SOLVER_SILENT = with_optimizer(Gurobi.Optimizer, OutputFlag = 0, Gurobi.Env())

include("utilities.jl")


text:
using Base.Test

using TALib
using TALib: TA_RetCode, TA_SUCCESS, _TA_COS

@testset "level 0" begin
    epsilon = 1e-10
    eps_price = 1e-6

    N_INDICATORS = 158
    @test length(D_INDICATORS) == N_INDICATORS
    @test length(INDICATORS) == N_INDICATORS
    @test length(D_GROUPS) == 10

    inReal = [0, pi/2, pi, 3pi/2, 0, pi/2, pi, 3pi/2]
    N = length(inReal)
    start_idx = 0
    end_idx = N - 1
    outBegIdx = Ref{Cint}(0)
    outNbElement = Ref{Cint}(0)
    outReal = fill(NaN, N)
    ret_code = _TA_COS(start_idx, end_idx, inReal, outBegIdx, outNbElement, outReal)
    errorCode = TA_RetCode(ret_code)
    println(outBegIdx[])
    println(outReal)
    outReal = circshift(outReal, outBegIdx[])
    println(outReal)
    @test errorCode == TA_SUCCESS::TA_RetCode
    println("angles: $inReal")
    println("COS(angles): $outReal")
    @test sum(outReal) ≈ 0 atol=epsilon
    @test outReal == cos.(inReal)

    #=
    filename = joinpath(basepath(), "ford_2012.csv")
    data = readdlm(filename, ',', skipstart=1)

    typealias Price Float64
    typealias Volume Int

    dt = Date[]
    opn = Price[]
    hig = Price[]
    low = Price[]
    cls = Price[]
    price = Price[]
    vol = Volume[]

    for i in 1:size(data, 1)
        push!(dt, Date(data[i, 1]))
        push!(opn, data[i, 2])
        push!(hig, data[i, 3])
        push!(low, data[i, 4])
        push!(cls, data[i, 5])
    end

    println(dt)
    println(opn)
    println(hig)
    println(low)
    println(cls)
    =#

end


text:
using TableView
using Test, WebIO, Tables, JSON

@testset "installation" begin
    version = readchomp(joinpath(@__DIR__, "..", "ag-grid.version"))
    @test isfile(joinpath(@__DIR__, "..", "deps", "ag-grid-$(version)", "ag-grid.js"))
end
@testset "named tuple table" begin
    nttable = [
        (a = 2.0, b = 3),
        (a = 3.0, b = 12)
    ]
    @test showtable(nttable) isa WebIO.Scope
end
@testset "named tuple table with missing and nothing" begin
    nttable = [
        (a = 2.0, b = 3, c = missing),
        (a = 3.0, b = 12, c = nothing)
    ]
    @test showtable(nttable) isa WebIO.Scope
end
@testset "inf, nan, and missing serializing" begin
    rows = Tables.table([NaN Inf -Inf 0 missing])
    names = [:a, :b, :c, :d, :e]
    types = vcat([Float64 for _ in 1:4], [Missing])
    Base.show(io::IO, x::Missing) = print(io, "test_missing")
    json = TableView.table2json(Tables.Schema(names, types), rows, types)
    firstrow = JSON.parse(json)[1]
    @test firstrow["a"] == "NaN"
    @test firstrow["b"] == "Inf"
    @test firstrow["c"] == "-Inf"
    @test firstrow["d"] == 0
    @test firstrow["e"] == "test_missing"
end
@testset "normal array" begin
    array = rand(10, 10)
    @test showtable(array) isa WebIO.Scope
end


text:
__precompile__()

"""
Merge functions from different modules into a single function.

This module provides two exported macros that can be used to merge different functions.

- [`@merge`](@ref)
- [`@kwmerge`](@ref)

"""
module MergedMethods

using Base.Meta, Compat

export @merge, @kwmerge

"""
Merge functions from different modules.

**Signature**

```julia
@merge(func, modules...)
```

**Examples**

```julia
module A

type T end

f(::T) = T

end

module B

type T end

f(::T) = T

end

@merge f A B

f(A.T()) == A.T
f(B.T()) == B.T
```

If several modules provide methods suitable for the specified arguments then the first
listed, in the above example it would be `A`, is selected.

**Generated Code**

`@merge` uses `@generated` to build specialised code for each tuple of arguments passed to
the merged function. The resulting code will typically have little to no overhead due to
inlining of the dispatch function.

```julia
module A

f(x::Int, y::Float64) = x + 2y

end

module B

f(x::Float64, y::Int) = x - 3y

end

@merge f A B

test(x, y) = f(x, y) + f(y, x)
@code_llvm test(1, 1.0)
```

```llvm
define double @julia_test_24571(i64, double) {
top:
  %2 = call double @"julia_*.1245"(i64 2, double %1)
  %3 = call double @"julia_+_24572"(i64 %0, double %2)
  %4 = mul i64 %0, 3
  %5 = call double @julia_-.1288(double %1, i64 %4)
  %6 = fadd double %3, %5
  ret double %6
}
```

*See also:* [`@kwmerge`](@ref).
"""
macro merge(func, modules...) buildmerge(func, false, modules) end

"""
Variant of `@merge` with support for passing keyword arguments.

**Signature**

```julia
@kwmerge(func, modules...)
```

**Examples**

```julia
@kwmerge f A B
f(1, a = 2, b =  3)
```

**Note:**

This macro should only be used when keywords are actually needed since the generated code
will probably not be as efficient as that of `@merge`.

*See also:* [`@merge`](@ref).
"""
macro kwmerge(func, modules...) buildmerge(func, true, modules) end

function buildmerge(f, keywords, mods)
    mapping = [(m, getfield(m, f)) for m in [getfield(current_module(), m) for m in mods]]
    m = :(($pickmodule)($mapping, args))
    q = quot(f)
    expr = keywords ?
        :($(f)(args...; kwargs...) = :($($m).$($q)(args...; kwargs...))) :
        :($(f)(args...)            = :($($m).$($q)(args...)))
    esc(:(@generated $expr))
end

function pickmodule(choices, args)
    for (m, fn) in choices
        method_exists(fn, args) && return module_name(m)
    end
    error("No suitable method found with arguments '$args'.")
end

end # module


text:
function get_inside_point(obj::InsideFinder)
    return jcall(obj, "getInsidePoint", S2Point, ())
end

function visit_internal_node(obj::InsideFinder, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::InsideFinder, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::InsideFinder, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_sharedataset

"""
    vtk_sharedataset(;kwargs...)
    vtk_sharedataset(children::Any;kwargs...)
    vtk_sharedataset(children_maker::Function;kwargs...)

A ShareDataSet component.
ShareDataSet capture a dataset or a source and allow it to use it in another
pipeline or representation.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `name` (String; optional): Unique dataset name to cross reference
- `port` (Real; optional): downstream connection port
"""
function vtk_sharedataset(; kwargs...)
        available_props = Symbol[:children, :id, :name, :port]
        wild_props = Symbol[]
        return Component("vtk_sharedataset", "ShareDataSet", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_sharedataset(children::Any; kwargs...) = vtk_sharedataset(;kwargs..., children = children)
vtk_sharedataset(children_maker::Function; kwargs...) = vtk_sharedataset(children_maker(); kwargs...)

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
use_mpi() = true
using .MPI


text:
using QXTools
using QXTools.Circuits
using ArgParse
using Logging

function parse_commandline(ARGS)
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--rows", "-r"
            help = "Number of rows"
            default = 3
            arg_type = Int64
        "--cols", "-c"
            help = "Number of columns"
            default = 3
            arg_type = Int64
        "--depth", "-d"
            help = "Number layers"
            default = 8
            arg_type = Int64
        "--seed"
            help = "Seed to use for circuit generation, contraction planning and selecting amplitudes to compute."
            default = nothing
        "--decompose"
            help = "Set if two qubit gates should be decomposed."
            default = true
            arg_type = Bool
        "--hypergraph"
            help = "Set if the hypergraph structure of the tensor network should be used."
            default = true
            arg_type = Bool
        "--time"
            help = "The number of seconds to run flow cutter for."
            default = 30
            arg_type = Int64
        "--sliced_bonds", "-s"
            help = "Number bonds to slice"
            default = 2
            arg_type = Int64
        "--score_function"
            help = "Function to maximise when selecting vertices to remove."
            default = :tree_trimming
            arg_type = Symbol
        "--prefix", "-p"
            help = "Prefix to use for output files"
            required = true
            arg_type = String
        "--verbose", "-v"
            help = "Verbose output"
            action = :store_true
    end
    import_settings!(s, create_output_parser())

    return parse_args(ARGS, s)
end

function main(ARGS)
    parsed_args = parse_commandline(ARGS)

    # Citcuit parameters.
    rows = parsed_args["rows"]
    cols = parsed_args["cols"]
    depth = parsed_args["depth"]

    # Tensor network parameters.
    decompose = parsed_args["decompose"]
    hypergraph = parsed_args["hypergraph"]

    # Contraction and slicing parameters.
    time = parsed_args["time"]
    score_function = parsed_args["score_function"]
    number_bonds_to_slice = parsed_args["sliced_bonds"]

    # Other parameters.
    seed = parsed_args["seed"]
    if seed !== nothing
        seed = parse(Int64, seed)
    end
    output_prefix = parsed_args["prefix"]
    verbose = parsed_args["verbose"]

    output_args = process_output_args(rows*cols, parsed_args)

    @info("Create circuit with $(rows * cols) qubits")
    circ = create_rqc_circuit(rows, cols, depth, seed)
    @info("Circuit created with $(circ.circ_ops.len) gates")

    generate_simulation_files(circ,
                              output_prefix,
                              number_bonds_to_slice;
                              decompose=decompose,
                              seed=seed,
                              output_args=output_args,
                              hypergraph=hypergraph,
                              time=time,
                              score_function=score_function)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text:
#! julia

using Pkg
Pkg.activate(@__DIR__)
Pkg.develop(PackageSpec(path=dirname(@__DIR__))) # adds the package this script is called from
Pkg.instantiate()
Pkg.update()

run = true
while run
    include("make.jl")

    using LiveServer; serve(dir=joinpath(@__DIR__, "build"))

    println("Run again? Enter! Exit witn 'q'.")
    if readline() == "q"
        global run = false
    end
end


text:
# using PyPlot, AcousticWave
# # =============dot product test for joint inversion================================
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
# v  = 3000. * ones(nz, nx);
#
# # initialize source term
# isx = []; isz = []; pic = zeros(nz, nx)
# for ix = 1 : nx
#     for iz = 1 : nz
#         if rand() < 0.001
#            push!(isx, ix); push!(isz, iz)
#            pic[iz, ix] = 1.0
#         end
#     end
# end
# pos_src = convert(Array{Int64,2}, hcat(isz, isx)); f0 = 30.0; ot = zeros(size(pos_src,1));
# src = InitSources(pos_src, f0, ot, dt);
# (tl, tu) = SrcRange(src)
# ntw = round(Int64, (tu-tl)/dt+1)
# path = "/Users/wenlei/Desktop/srcSpt.bin"
# src2srcSpts(path, src, nz, nx)
#
# # finite difference coefficients
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
#
# # receiver positions
# irx = []; irz = []; pic_rec = zeros(nz, nx)
# for ix = 1 : nx
#     for iz = 1 : nz
#         if rand() < 0.01
#            push!(irx, ix); push!(irz, iz)
#            pic_rec[iz, ix] = 1.0
#         end
#     end
# end
# pos_rec = convert(Array{Int64,2}, hcat(irz, irx));
#
# # modeling
# shot = MultiStepForward(pos_rec, path, fidMtx, tmax=0.5);
# shot1= MultiStepForward(pos_rec, src, fidMtx, tmax=0.5);
#
# nt = shot.nt
# shot2 = InitShot(0, 0, pos_rec, 0.0, nt, dt)
# nrec  = size(pos_rec, 1)
# w = Ricker(40.0, dt)
# ls = floor(Int64, length(w)/2)
# for irec = 1 : nrec
#     tmp  = conv(rand(nt), w)
#     shot2.d[:, irec] = tmp[ls+1:end-ls]
# end
# path_adj = "/Users/wenlei/Desktop/adjoint.bin"
# MultiStepAdjoint(ntw, path_adj, shot2, fidMtx)
#
# tmp = dot(vec(shot.d), vec(shot2.d))
# tmp1 = 0.0
# for it = 1 : ntw
#     p = readSrcSpt(path, it)
#     p1 = readSrcSpt(path_adj, it)
#     tmp1 = tmp1 + dot(vec(p), vec(p1))
# end

using PyPlot, AcousticWave
# =============the result of adjoint================================
nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
# initialize source term
pos_src = [101 75; 101 125]; f0 = 30.0; ot = zeros(size(pos_src,1));
src = InitSources(pos_src, f0, ot, dt);
(tl, tu) = SrcRange(src)
ntw = round(Int64, (tu-tl)/dt+1)

vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
pos_rec = hcat(iz, ix);
shot = MultiStepForward(pos_rec, src, fidMtx, tmax=0.5);
# path_adj = "/Users/wenlei/Desktop/adj.bin"
# MultiStepAdjoint(ntw, path_adj, shot, fidMtx)

path_x = "/Users/wenlei/Desktop/x.bin"
path_s = "/Users/wenlei/Desktop/s.bin"
path_p = "/Users/wenlei/Desktop/p.bin"
CGjoint(ntw, path_x, path_s, path_p, shot, fidMtx, pos_rec, tmax=0.5, mu=0.1, niter=30)
SeisPlot(readSrcSpt(path_x, 10), dx=10, dy=10, xlabel="Length", xunits="(m)", ylabel="Depth", yunits="(m)")

shot_pre = MultiStepForward(pos_rec, path_x, fidMtx, tmax=0.5)
it = 40
plot(shot.d[:,it])
plot(shot_pre.d[:,it])

nz = 100; nx = 77; ot = 0.0; dt = 0.001; nt = 33;
path1 = "/Users/wenlei/Desktop/tmp1.bin"
path2 = "/Users/wenlei/Desktop/tmp2.bin"
fid1 = open(path1, "w");
fid2 = open(path2, "w");
write(fid1, nz, nx, ot, dt, nt)
write(fid2, nz, nx, ot, dt, nt)
for it = 1 : nt
    write(fid1,     ones(nz*nx))
    write(fid2, 2.0*ones(nz*nx))
end
close(fid1)
close(fid2)

W = ones(nz*nx)
mu = -0.5
obtainS(path1, mu, W, path2)

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
abstract type SimulationState end

system(::SimulationState) = error("Unimplemented")
Base.time(::SimulationState) = error("Unimplemented")
positions(state::SimulationState) = positions(system(state))
velocities(state::SimulationState) = velocities(system(state))
periods(state::SimulationState) = periods(system(state))

abstract type SimulationTape{TM <: MosiModel} end

struct SimpleTape{T, TM <: MosiModel{T}} <: SimulationTape{TM}
    times::Vector{Float64}
    positions::Vector{Vector{T}}
    velocities::Vector{Vector{T}}
    periods::Vector{Vector{T}}
    model::TM
end
SimulationTape(ts, rs) = SimpleTape(ts, rs, [], [], UnknownModel())
SimulationTape(ts, rs, vs, ps) = SimpleTape(ts, rs, vs, ps, UnknownModel())
mosi_model(tape::SimulationTape) = tape.model
natoms(tape::SimulationTape{UnknownModel}) = length(tape.positions[1])
natoms(tape::SimulationTape) = natoms(tape.model)
Base.length(tape::SimulationTape) = length(tape.times)
times(tape::SimulationTape) = tape.times
positions(tape::SimulationTape) = tape.positions
velocities(tape::SimulationTape) = tape.velocities
periods(tape::SimulationTape) = tape.periods

times(tape::SimulationTape, i) = tape.times[i]
positions(tape::SimulationTape, i) = tape.positions[i]
velocities(tape::SimulationTape, i) = tape.velocities[i]
periods(tape::SimulationTape{T}, i) where T = length(tape.periods) ≡ 0 ? T[] : tape.periods[i]

function get_configuration_func(
    tape::SimulationTape;
    atom_range = 1:natoms(tape)
)
    rss = positions(tape)
    pss = periods(tape)
    model = tape.model
    box = pbc_box(model)
    (i::Int) -> ConfigurationSystem(
        view(rss[i], atom_range),
        view(pss[i], atom_range),
        box
    )
end

abstract type SimulationSetup end
mosi_model(::SimulationSetup) = UnknownModel()
init_state(::SimulationSetup) = error("Unimplemented")

abstract type SimulationResult end
tape(result::SimulationResult) = result.tape
observables(result::SimulationResult) = result.observables

struct SimulationError <: Exception
    msg::String
end

Base.showerror(io::IO, e::SimulationError) = print(io, e.msg)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Plots
using Plots.PlotMeasures
pyplot(fmt = :svg)
bar(rand(10), c=:white, right_margin=1cm)
plot!(twinx(), rand(10), label="v1")
plot!(twinx(), rand(10) * 10, label="v2")

# %%
using Plots
using Plots.PlotMeasures
pyplot(fmt = :svg)
bar(rand(10), c=:white, right_margin=1cm)
plot!(twinx(), [rand(10) rand(10) * 10], label=["v1" "v2"])

# %%
using Plots
using Plots.PlotMeasures
pyplot(fmt = :svg)
plot(rand(10), label="v1", right_margin=1cm)
plot!(rand(10) * 10, label="v2")
bar!(twinx(), rand(10), label="", fillalpha=0)

# %%
using Plots
using Plots.PlotMeasures
pyplot(fmt = :svg)
bar(rand(10), c=:white, right_margin=1cm)
ax2 = twinx()
plot!(ax2, rand(10), label="v1")
plot!(ax2, rand(10) * 10, label="v2")

# %%

max_stars_repo_path: 0010/How to use Plots.twinx.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0010/How to use Plots.twinx.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0010/How to use Plots.twinx.jl

max_forks_repo_name: genkuroki/public


text:
__precompile__()
module HilbertSpaceFillingCurve

    export hilbert

    const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    if isfile(depsfile)
        include(depsfile)
    else
        error("HilbertSpaceFillingCurve not properly installed. Please run Pkg.build(\"HilbertSpaceFillingCurve\") then restart Julia.")
    end

    const bitmask_t = Culonglong
    global const bits_per_byte = 8

    function hilbert(d::T, ndims, nbits = 32) where T <: Integer
        @assert ndims*nbits <= sizeof(bitmask_t) * bits_per_byte

        p = bitmask_t.(zeros(ndims))
        ccall((:hilbert_i2c, libhilbert), Nothing, (Int,Int,bitmask_t,Ptr{bitmask_t}),ndims,nbits,d,p)
        Int.(p)
    end

    function hilbert(p::Vector{T}, ndims, nbits = 32) where T <: Integer
        @assert ndims*nbits <= sizeof(bitmask_t) * bits_per_byte
        @assert length(p) == ndims

        Int(ccall((:hilbert_c2i, libhilbert), bitmask_t, (Int,Int,Ptr{bitmask_t}),ndims,nbits, bitmask_t.(p)))
    end

end # module


text:
struct PDM <: MethodType
    c::Float64
    o::Float64
    norms::Vector{Float64}
end

@inline function PDM(c::Float64, o::Float64)
    norms = Vector{Float64}(undef,12000)
    previous = 0.0

    for cur_pos in 0:11999
        previous = cur_pos == 0 ? 0.0 : norms[cur_pos]
        norms[cur_pos + 1] = previous + exp(-c * cur_pos)
    end

    return PDM(c, o, norms)
end

computeDocFreq(method::PDM, docFreq) = docFreq.doc_freq

gap_score(method::PDM, pos, ndf) = 0.0
match_score(method::PDM, q_pos, c_pos, ndf, q_gap_score, c_gap_score) =  exp(-method.c * (min(q_pos, c_pos) - 1.0)) * exp(-method.o * abs(q_pos - c_pos))

# @inline function compute_final_similarity!(similarities, method::PDM, query, candidates, df_vec)
#     query_len = query.stats.doc_len

#     @simd for i = eachindex(similarities)
#         norm = method.norms[min(query_len, candidates[i].stats.doc_len)]
#         @inbounds similarities[i] /= norm
#     end

#     return similarities
# end

function similarity(method::PDM, query::Vector{UInt32}, candidate::Vector{UInt32}, doc_freq::Vector{Float64})
    q_len = length(query)
    c_len = length(candidate)

    row = zeros(c_len)

    @inbounds for (q_pos, q_func) in enumerate(query)
        diagonal = 0.0

        for (c_pos, c_func) in enumerate(candidate)
            # Align gap to query position
            above = row[c_pos]

            # Align gap to candidate position
            left = c_pos == 1 ? 0.0 : row[c_pos - 1]

            if q_func == c_func
                # Matching
                diagonal += match_score(method, q_pos, c_pos, 0.0, 0.0, 0.0)
            end

            diagonal_tmp = row[c_pos]

            # Replace max with IF ELSE considerably improves the performance =O
            # row[c_pos] = max(left, above, diagonal)
            m = left

            if above > m
                m = above
            end

            if diagonal > m
                m = diagonal
            end

            row[c_pos] = m

            diagonal = diagonal_tmp
        end
    end

    @inbounds sim = row[c_len]
    norm = method.norms[min(q_len, c_len)]

    return sim / norm
end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
## This file contains im2col-backed implementations of convolution for 2d and 3d
## convolutions.  Expect to see a lot of indexing.

# Helper functions for flipkernel-induced dyslexia
@inline function kernel_index(w, h, d, cdims::ConvDims{N, S, P, D, false}) where {N, S, P, D}
    kernel_w, kernel_h, kernel_d = kernel_size(cdims)
    return (kernel_w - w + 1, kernel_h - h + 1, kernel_d - d + 1)
end
@inline function kernel_index(w, h, d, cdim::ConvDims{N, S, P, D, true}) where {N, S, P, D}
    return (w, h, d)
end

"""
    conv_im2col!(y, x, w, cdims, col=similar(x); alpha=1, beta=0)

Perform a convolution using im2col and GEMM, store the result in `y`.  The  kwargs
`alpha` and `beta` control accumulation behavior; internally this operation is
implemented as a matrix multiply that boils down to `y = alpha * x * w + beta * y`, thus
by setting `beta` to a nonzero value, multiple results can be accumulated into `y`, or
by setting `alpha` to a nonunitary value, various gain factors can be applied.

Note for the particularly performance-minded, you can provide a pre-allocated `col`,
which should eliminate any need for large allocations within this method.
"""
@timeit_debug to function conv_im2col!(
                y::AbstractArray{T,5}, x::AbstractArray{T,5},
                w::AbstractArray{T,5}, cdims::DenseConvDims;
                col::AbstractArray{T,2}=similar(x, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where {T}
    check_dims(size(x), size(w), size(y), cdims)

    #   COL   *    W    ->    Y
    # [M x K] * [K x N] -> [M x N]
    #
    #  M: output spatial resolution
    #  N: output channels
    #  K: size of input "patch" (kernel size and input channels combined)
    #
    # In english, we're grabbing each input patch and laying them out along
    # the M dimension in `col`, so that the GEMM call below multiplies each
    # kernel (which is kernel_h * kernel_w * channels_in elments long) is
    # dotproducted with that input patch, effectively computing a convolution
    # in a somewhat memory-wasteful but easily-computed way (since we already
    # have an extremely highly-optimized GEMM call available in BLAS).
    M = prod(output_size(cdims))
    N = channels_out(cdims)
    K = prod(kernel_size(cdims))*channels_in(cdims)

    @inbounds for batch_idx in 1:size(x,5)
        # We invoke `@timeit_debug` on the outside of `im2col!()` because inference
        # doesn't like us putting it on the inside.
        @timeit_debug to "im2col!" im2col!(col, view(x, :, :, :, :, batch_idx), cdims)
        GC.@preserve col, w, y, begin
            col_ptr = pointer(col)
            w_ptr = pointer(w)
            y_ptr = pointer(y, (batch_idx - 1)*M*N + 1)
            @timeit_debug to "gemm!" gemm!(Val(false), Val(false), M, N, K, alpha, col_ptr, w_ptr, beta, y_ptr)
        end
    end
    return y
end

"""
    ∇conv_filter_im2col!(dw, x, dy, cdims, col=similar(dw); alpha=1, beta=0)

Conv backward pass onto the weights using im2col and GEMM; stores the result in `dw`.
See the documentation for `conv_im2col!()` for explanation of optional parameters.
"""
@timeit_debug to function ∇conv_filter_im2col!(
                dw::AbstractArray{T,5}, x::AbstractArray{T,5},
                dy::AbstractArray{T,5}, cdims::DenseConvDims;
                col::AbstractArray{T,2} = similar(dw, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where {T}
    check_dims(size(x), size(dw), size(dy), cdims)

    #   COL'   *   dY   ->    dW
    # [M x K] * [K x N] -> [M x N]
    #
    #  M: size of input "patch" (kernel size and input channels combined)
    #  N: output channels
    #  K: output spatial resolution
    #
    # In english, we're grabbing each input patch and laying them out along
    # the K dimension in `col`, then multiplying in `dY` to compute a dot
    # product between all pixels in the input that were multiplied by a single
    # position in the W kernel, and all output pixels of the same location,
    # across output channels.  This slice of `col` therefore constitutes every
    # input pixel that touched a particular element of the kernel.
    #
    # This is identical to a convolution between x and a dimension-permuted dY,
    # where we

    M = prod(kernel_size(cdims))*channels_in(cdims)
    N = channels_out(cdims)
    K = prod(output_size(cdims))

    @inbounds for batch_idx in 1:size(x,5)
        # We invoke `@timeit_debug` on the outside of `im2col!()` because inference
        # doesn't like us putting it on the inside.
        @timeit_debug to "im2col!" im2col!(col, view(x, :, :, :, :, batch_idx), cdims)
        GC.@preserve col, dw, dy, begin
            col_ptr = pointer(col)
            dy_ptr = pointer(dy,(batch_idx - 1)*K*N + 1)
            dw_ptr = pointer(dw)
            @timeit_debug to "gemm!" gemm!(Val(true), Val(false), M, N, K, alpha, col_ptr, dy_ptr, beta, dw_ptr)
        end

        # Because we accumulate over batches in this loop, we must set `beta` equal
        # to `1.0` from this point on.
        beta = T(1)
    end
    return dw
end

"""
    ∇conv_data_im2col!(dx, w, dy, cdims, col=similar(dx); alpha=1, beta=0)

Conv2d backward pass onto the input using im2col and GEMM; stores the result in `dx`.
See the documentation for `conv_im2col!()` for explanation of other parameters.
"""
@timeit_debug to function ∇conv_data_im2col!(
                dx::AbstractArray{T,5}, dy::AbstractArray{T,5},
                w::AbstractArray{T,5}, cdims::DenseConvDims;
                col::AbstractArray{T,2} = similar(dx, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where {T}
    check_dims(size(dx), size(w), size(dy), cdims)

    #    dY        W'   ->    dX
    # [M x K] * [K x N] -> [M x N]
    #
    #  M: output spatial resolution
    #  N: size of input "patch" (kernel size and input channels combined)
    #  K: output channels
    #
    # In english, we're taking the output image and laying it out by pixel,
    # with channels lying along the `K` dimension in `col`.  We then multiply
    # in `W'` to compute a dot product between each pixel location and the
    # entire kernel.  This dot product therefore constitutes every output pixel
    # that was a function of a particular input pixel.
    #
    # This is identical to a transposed convolution between dY and W

    M = prod(output_size(cdims))
    N = prod(kernel_size(cdims))*channels_in(cdims)
    K = channels_out(cdims)

    @inbounds for batch_idx in 1:size(dx, 5)
        GC.@preserve col, w, dy, begin
            dy_ptr = pointer(dy, (batch_idx - 1)*M*K + 1)
            w_ptr = pointer(w)
            col_ptr = pointer(col)
            @timeit_debug to "gemm!" gemm!(Val(false), Val(true), M, N, K, alpha, dy_ptr, w_ptr, T(0), col_ptr)
        end
        @timeit_debug to "col2im!" col2im!(view(dx, :, :, :, :, batch_idx), col, cdims)
    end
    return dx
end

"""
    im2col!(col, x, cdims)

Converts a 3d image `x` into a matrix `col` for usage with GEMM-calculated convolution.
Patches of `x` of size (kernel_w, kernel_h, kernel_d, C_in) will be extracted and laid
out along the rows of `col`, one for each output pixel.  This routine is used by all
im2col-based convolutions, just with extra singleton dimensions added in the case of `2d`
or `1d` images.
"""
function im2col!(col::AbstractArray{T,2}, x::AbstractArray{T,4},
                                          cdims::ConvDims) where {T}
    if spatial_dims(cdims) != 3
        throw(DimensionMismatch("im2col!() only accepts 3d convoluitional inputs"))
    end

    # Extract those nice, compile-time constant type parameters from `cdims`.
    width, height, depth = input_size(cdims)
    kernel_w, kernel_h, kernel_d = kernel_size(cdims)
    C_in = channels_in(cdims)
    pad_w_lo, pad_w_hi, pad_h_lo, pad_h_hi, pad_d_lo, pad_d_hi = padding(cdims)
    dil_w, dil_h, dil_d = dilation(cdims)
    stride_w, stride_h, stride_d = stride(cdims)
    out_width, out_height, out_depth = output_size(cdims)

    # Reshape col for easy access.
    col_reshaped = reshape(col, (
        # Output resolution
        out_width,
        out_height,
        out_depth,

        # By input patch size
        kernel_w,
        kernel_h,
        kernel_d,
        C_in,
    ))

    padded_regions, central_region = calc_padding_regions(cdims)

    # A helper function to project from output (w, h) to input (input_w, input_h)
    @inline project(idx, stride, pad) = (idx - 1)*stride - pad + 1

    # We begin by copying the central region of the image which requires no padding at all.
    # Eliminating the branches of the fully generalized version below gives us a nice
    # speedup on the majority of the data.
    @timeit_debug to "im2col!() - central region" begin
        @inbounds for c in 1:C_in
            # Unpack "central region"
            w_region, h_region, d_region = central_region

            for kd in 1:kernel_d,
                kh in 1:kernel_h,
                kw in 1:kernel_w,
                d in d_region,
                h in h_region,
                w in w_region

                input_kd = project(d, stride_d, pad_d_lo) + (kd - 1)*dil_d
                input_kh = project(h, stride_h, pad_h_lo) + (kh - 1)*dil_h
                input_kw = project(w, stride_w, pad_w_lo) + (kw - 1)*dil_w
                kidxs = kernel_index(kw, kh, kd, cdims)

                xval::T = x[input_kw, input_kh, input_kd, c]
                col_reshaped[w, h, d, kidxs..., c] = xval
            end
        end
    end

    # For each "padded region", we run the fully general version
    @timeit_debug to "im2col!() - padded region" begin
        @inbounds for (w_region, h_region, d_region) in padded_regions
            for c in 1:C_in,
                d in d_region,
                h in h_region,
                w in w_region,
                kd in 1:kernel_d,
                kh in 1:kernel_h,
                kw in 1:kernel_w

                input_kd = project(d, stride_d, pad_d_lo) + (kd - 1)*dil_d
                input_kh = project(h, stride_h, pad_h_lo) + (kh - 1)*dil_h
                input_kw = project(w, stride_w, pad_w_lo) + (kw - 1)*dil_w

                kidxs = kernel_index(kw, kh, kd, cdims)

                # If this d is off the edge, then deal with the entire plane
                # in one fell swoop, like a ravenous flock of crows.  CAW CAW.
                if input_kd <= 0 || input_kd > depth
                    for kh in 1:kernel_h,
                        kw in 1:kernel_w
                        col_reshaped[w, h, d, kidxs..., c] = T(0)
                    end
                    continue
                end

                # Same for `h`, but in this case it's only a line, not a plane.
                # This results in slightly less caw'ing.
                if input_kh <= 0 || input_kh > height
                    for kw in 1:kernel_w
                        col_reshaped[w, h, d, kidxs..., c] = T(0)
                    end
                    continue
                end

                # If this `w` is off the edge it and only it gets cleared out
                if input_kw <= 0 || input_kw > width
                    col_reshaped[w, h, d, kidxs..., c] = T(0)
                    continue
                end

                # Copy the data over
                xval::T = x[input_kw, input_kh, input_kd, c]
                col_reshaped[w, h, d, kidxs..., c] = xval
            end
        end
    end
end

"""
    col2im!(x, col, cdims)

Does the inverse of `im2col!()`, converting `col` back into a 3d image, used for backward
passes, transposed convolutions, etc...

Note that this method has not been optimized in the same way as `im2col()` has, because
it is slightly more complicated due to the more chaotic data access patterns, and I'm not
desperate enough yet.
"""
col2im!

function col2im!(x::AbstractArray{T,4}, col::AbstractArray{T,2},
                                  cdims::ConvDims) where T
    if spatial_dims(cdims) != 3
        throw(DimensionMismatch("col2im!() only accepts 3d convoluitional inputs"))
    end

    # Extract those nice, compile-time constant type parameters from `cdims`.
    width, height, depth = input_size(cdims)
    kernel_w, kernel_h, kernel_d = kernel_size(cdims)
    C_in = channels_in(cdims)
    pad_w_lo, pad_w_hi, pad_h_lo, pad_h_hi, pad_d_lo, pad_d_hi = padding(cdims)
    dil_w, dil_h, dil_d = dilation(cdims)
    stride_w, stride_h, stride_d = stride(cdims)
    out_width, out_height, out_depth = output_size(cdims)

    # TODO: Rewrite this method so we don't have this fill!() at the beginning!
    # Calculate each output pixel once rather than accumulating into it?
    fill!(x, T(0))

    # Reshape col for easy access.
    col_reshaped = reshape(col, (
        # Output resolution
        out_width,
        out_height,
        out_depth,

        # By input patch size
        kernel_w,
        kernel_h,
        kernel_d,
        C_in,
    ))

    # A helper function to project from output (w, h) to input (input_w, input_h)
    @inline project(idx, stride, pad) = (idx - 1)*stride - pad + 1

    @inbounds for c in 1:C_in
        for kd in 1:kernel_d,
            kh in 1:kernel_h,
            kw in 1:kernel_w

            for d in 1:out_depth
                input_kd = project(d, stride_d, pad_d_lo) + (kd - 1)*dil_d

                # If this d is off the edge, then deal with the entire plane
                # in one fell swoop, like a ravenous flock of crows.  CAW CAW.
                if input_kd <= 0 || input_kd > depth
                    continue
                end

                for h in 1:out_height
                    input_kh = project(h, stride_h, pad_h_lo) + (kh - 1)*dil_h

                    # Same for `h`, but in this case it's only a line, not a plane.
                    # This results in slightly less caw'ing.
                    if input_kh <= 0 || input_kh > height
                        continue
                    end

                    for w in 1:out_width
                        input_kw = project(w, stride_w, pad_w_lo) + (kw - 1)*dil_w

                        # If this `w` is off the edge, only it gets cleared out.
                        if input_kw <= 0 || input_kw > width
                            continue
                        end

                        # Copy the data over
                        kidxs = kernel_index(kw, kh, kd, cdims)
                        cval::T = col_reshaped[w, h, d, kidxs..., c]
                        x[input_kw, input_kh, input_kd, c] += cval
                    end
                end
            end
        end
    end
end


text:
include("dependencies_for_runtests.jl")
include("dependencies_for_poisson_solvers.jl")
using GPUArrays

#####
##### Run some basic AMD tests
##### - this test file is meant to be used for early development for
##### - quick CI using amdgpu buildkite queue. Using standard tests produces
##### - so much error chatter it is hard to unravel what to fix. Standard tests
##### - are also  possibly overly exhaustive which is a bottleneck in turn around.
#####

@testset "Quick AMD" begin
    @info "Testing AMD basics..."

    for arch in archs
        @testset "Quick AMD set field [$(typeof(arch))]" begin

            @info "  Quick AMD set field from func [$(typeof(arch))]..."
            N = 3
            topo = (Bounded, Bounded, Bounded)
            FT=Float64
            f(x, y, z) = 1 + exp(x) * sin(y) * tanh(z)
            grid = RectilinearGrid(arch, FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            u = XFaceField(grid)
            set!(u, f )
            ans=@allowscalar u .- 0
            # Check against CPU
            grid = RectilinearGrid(CPU(), FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            u = XFaceField(grid)
            set!(u, f )
            ans_c=@allowscalar u .- 0
            @test ans ≈ ans_c

            @info "  Quick AMD set field from GPU array [$(typeof(arch))]..."
            N = 3
            topo = (Bounded, Bounded, Bounded)
            FT=Float64
            grid = RectilinearGrid(arch, FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            Nx, Ny, Nz = size(grid)
            f = rand(Nx, Ny, Nz)
            f_d = ROCArray(f)
            default_bcs = FieldBoundaryConditions()
            u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
            u = CenterField(grid,boundary_conditions=u_bcs)
            set!(u, f_d )
            ans=@allowscalar u .- 0
            # Check against CPU
            grid = RectilinearGrid(CPU(), FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            Nx, Ny, Nz = size(grid)
            u = CenterField(grid,boundary_conditions=u_bcs)
            set!(u, f )
            ans_c=@allowscalar u .- 0
            @test ans ≈ ans_c

            @info "  Quick AMD set field from CPU array [$(typeof(arch))]..."
            N = 3
            topo = (Bounded, Bounded, Bounded)
            FT=Float64
            grid = RectilinearGrid(arch, FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            Nx, Ny, Nz = size(grid)
            f = rand(Nx, Ny, Nz)
            default_bcs = FieldBoundaryConditions()
            u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
            u = CenterField(grid,boundary_conditions=u_bcs)
            set!(u, f )
            ans=@allowscalar u .- 0
            # Check against CPU
            grid = RectilinearGrid(CPU(), FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))
            default_bcs = FieldBoundaryConditions()
            u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
            u = CenterField(grid,boundary_conditions=u_bcs)
            set!(u, f )
            ans_c=@allowscalar u .- 0

            @test ans ≈ ans_c

        end
    end
end


text:
using JSCall, WebIO, JSExpr
using Test
using JSExpr: jsexpr
using JSCall
using WebIO
THREE, document = JSModule(
    :THREE,
    "https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js",
)

scene = THREE.new.Scene()
width, height = 500, 500
# Create a basic perspective camera
camera = THREE.new.PerspectiveCamera(75, width / height, 0.1, 1000)
camera.position.z = 4
renderer = THREE.new.WebGLRenderer(Dict(:antialias => true))
renderer.setClearColor("#ffffff")
geometry = THREE.new.BoxGeometry(1, 1, 1)
material = THREE.new.MeshBasicMaterial(color = "#433F81")
cube = THREE.new.Mesh(geometry, material);
scene.add(cube)
container = document.querySelector("#container")
container.appendChild(renderer.domElement)
display((scope(THREE))(dom"div#container"()))

@jsfun function render(unused)
    @var cube = $cube
    @var renderer = $renderer
    @var scene = $scene
    @var camera = $camera
    cube.rotation.x += 0.01
    cube.rotation.y += 0.01
    renderer.render(scene, camera)
    requestAnimationFrame(render)
end

render(scene)