text:
using Test
using MPI
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using LinearAlgebra
using Printf
using Random

include("advection_diffusion_model.jl")

struct Box{dim} end
struct Sphere end

vertical_unit_vector(::Box{2}, ::SVector{3}) = SVector(0, 1, 0)
vertical_unit_vector(::Box{3}, ::SVector{3}) = SVector(0, 0, 1)
vertical_unit_vector(::Sphere, coord::SVector{3}) = coord / norm(coord)

projection(::EveryDirection, ::SVector{3}) = I
projection(::VerticalDirection, k::SVector{3}) = k * k'
projection(::HorizontalDirection, k::SVector{3}) = I - k * k'

struct TestProblem{adv, diff, dir, topo} <: AdvectionDiffusionProblem end

initial_ρ(::Box, x) = prod(sin.(π * x))
function initial_ρ(::Sphere, x)
    r = norm(x)
    φ = atan(x[2], x[1])
    θ = atan(sqrt(x[1]^2 + x[2]^2), x[3])
    return sin(π * (r - 1)) * sin(φ) * sin(θ)
end

velocity(::Box, x) = sin.(π * x)
function velocity(::Sphere, x)
    r = norm(x)
    φ = atan(x[2], x[1])
    θ = atan(sqrt(x[1]^2 + x[2]^2), x[3])
    return sin(π * (r - 1)) .*
           SVector(cos(φ) * cos(θ), sin(φ) * cos(θ), cos(φ) * sin(θ))
end

function init_velocity_diffusion!(
    ::TestProblem{adv, diff, dir, topo},
    aux::Vars,
    geom::LocalGeometry,
) where {adv, diff, dir, topo}
    k = vertical_unit_vector(topo, geom.coord)
    P = projection(dir, k)
    aux.advection.u =
        !isnothing(adv) ? P * velocity(topo, geom.coord) : zeros(SVector{3})
    aux.diffusion.D =
        !isnothing(diff) ? SMatrix{3, 3}(P) / 200 : zeros(SMatrix{3, 3})
end

function initial_condition!(
    ::TestProblem{adv, diff, dir, topo},
    state,
    aux,
    localgeo,
    t,
) where {adv, diff, dir, topo}
    state.ρ = initial_ρ(topo, localgeo.coord)
end

function create_topology(::Box{dim}, mpicomm, Ne, FT) where {dim}
    brickrange = ntuple(j -> range(FT(0); length = Ne + 1, stop = 1), dim)
    periodicity = ntuple(j -> false, dim)
    bc = ntuple(j -> (3, 3), dim)
    StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = bc,
    )
end

function create_topology(::Sphere, mpicomm, Ne, FT)
    vert_range = grid1d(FT(1), FT(2), nelem = Ne)
    StackedCubedSphereTopology(mpicomm, Ne, vert_range, boundary = (3, 3))
end

create_dg(model, grid, direction) = DGModel(
    model,
    grid,
    RusanovNumericalFlux(),
    CentralNumericalFluxSecondOrder(),
    CentralNumericalFluxGradient(),
    direction = direction,
)

function test_run(
    adv,
    diff,
    topo,
    mpicomm,
    ArrayType,
    FT,
    polynomialorder,
    Ne,
    level,
)
    topology = create_topology(topo(), mpicomm, Ne, FT)
    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
        meshwarp = topo == Sphere ? cubedshellwarp : (x...) -> identity(x),
    )

    problems = (
        p = TestProblem{adv, diff, EveryDirection(), topo()}(),
        vp = TestProblem{adv, diff, VerticalDirection(), topo()}(),
        hp = TestProblem{adv, diff, HorizontalDirection(), topo()}(),
    )

    models = map(AdvectionDiffusion{3}, problems)

    dgmodels = map(models) do m
        (
            dg = create_dg(m, grid, EveryDirection()),
            vdg = create_dg(m, grid, VerticalDirection()),
            hdg = create_dg(m, grid, HorizontalDirection()),
        )
    end

    Q = init_ode_state(dgmodels.p.dg, FT(0), init_on_cpu = true)
    # do one Euler step to trigger numerical fluxes in subsequent evaluations
    let
        dt = 1e-3
        dQ = similar(Q)
        dgmodels.p.dg(dQ, Q, nothing, FT(0))
        Q .+= dt .* dQ
    end

    # evaluate all combinations
    dQ = map(
        x -> map(dg -> (dQ = similar(Q); dg(dQ, Q, nothing, FT(0)); dQ), x),
        dgmodels,
    )

    # set up tolerances
    atolm = 6e-13
    atolv = 3e-4 / 5^(level - 1)
    atolh = 0.0003 / 5^(level - 1)

    @testset "total" begin
        atol = topo <: Box || diff == nothing ? atolm : atolh
        @test isapprox(norm(dQ.p.dg .- dQ.p.vdg .- dQ.p.hdg), 0, atol = atol)
        @test isapprox(norm(dQ.vp.dg .- dQ.vp.vdg .- dQ.vp.hdg), 0, atol = atol)
        @test isapprox(norm(dQ.hp.dg .- dQ.hp.vdg .- dQ.hp.hdg), 0, atol = atol)
    end

    @testset "vertical" begin
        atol = topo <: Box ? atolm : atolv
        @test isapprox(norm(dQ.vp.dg .- dQ.vp.vdg), 0, atol = atol)
        @test isapprox(norm(dQ.vp.hdg), 0, atol = atol)
        @test isapprox(norm(dQ.vp.dg .- dQ.p.vdg), 0, atol = atol)
    end

    @testset "horizontal" begin
        atol = topo <: Box ? atolm : atolh
        @test isapprox(norm(dQ.hp.dg .- dQ.hp.hdg), 0, atol = atol)
        @test isapprox(norm(dQ.hp.vdg), 0, atol = atol)
        @test isapprox(norm(dQ.hp.dg - dQ.p.hdg), 0, atol = atol)
    end
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD
    FT = Float64
    numlevels = 2
    base_num_elem = 4

    # This test doesn't do any heavy computational work,
    # but compiles a lot of model/discretization combinations.
    # Compilation times on the GPU are longer, so we only run one
    # variable-degree case on the GPU
    if ArrayType == Array
        polynomialorders = ((4, 4), (4, 2))
    else
        polynomialorders = ((4, 2),)
    end

    @info @sprintf """Test parameters:
    ArrayType                   = %s
    FloatType                   = %s
    Polynomial orders           = %s
      """ ArrayType FT polynomialorders

    @testset "$(@__FILE__)" begin
        @testset for polyorders in polynomialorders
            @testset for topo in (Box{2}, Box{3}, Sphere)
                @testset for (adv, diff) in (
                    (Advection, nothing),
                    (nothing, Diffusion),
                    (Advection, Diffusion),
                )
                    @testset for level in 1:numlevels
                        Ne = 2^(level - 1) * base_num_elem
                        test_run(
                            adv,
                            diff,
                            topo,
                            mpicomm,
                            ArrayType,
                            FT,
                            polyorders,
                            Ne,
                            level,
                        )
                    end
                end
            end
        end
    end
end
nothing

max_forks_repo_name: jm-c/CLIMA


text:
include("affine.jl")
include("quadratic.jl")

ext: jl


text:
"""
    stress_centrality(g[, vs])
    stress_centrality(g, k)

Calculate the [stress centrality](http://med.bioinf.mpi-inf.mpg.de/netanalyzer/help/2.7/#stressDist)
of a graph `g` across all vertices, a specified subset of vertices `vs`, or a random subset of `k`
vertices. Return a vector representing the centrality calculated for each node in `g`.

The stress centrality of a vertex ``n`` is defined as the number of shortest paths passing through ``n``.

### References
- Barabási, A.L., Oltvai, Z.N.: Network biology: understanding the cell's functional organization. Nat Rev Genet 5 (2004) 101-113
- Shimbel, A.: Structural parameters of communication networks. Bull Math Biophys 15 (1953) 501-507.

# Examples
```jldoctest
julia> using LightGraphs

julia> stress_centrality(star_graph(3))
3-element Array{Int64,1}:
 2
 0
 0

julia> stress_centrality(cycle_graph(4))
4-element Array{Int64,1}:
 2
 2
 2
 2
```
"""
function stress_centrality(g::AbstractGraph, vs::AbstractVector=vertices(g))
    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    stress = zeros(Int, n_v)
    for s in vs
        if degree(g, s) > 0
            state = dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            _stress_accumulate_basic!(stress, state, g, s)
        end
    end
    return stress
end

stress_centrality(g::AbstractGraph, k::Integer) =
    stress_centrality(g, sample(vertices(g), k))

function _stress_accumulate_basic!(stress::Vector{<:Integer},
    state::DijkstraState,
    g::AbstractGraph,
    si::Integer)

    n_v = length(state.parents) # this is the ttl number of vertices
    δ = zeros(Int, n_v)
    P = state.predecessors

    laststress = copy(stress)
    # make sure the source index has no parents.
    P[si] = []
    # we need to order the source vertices by decreasing distance for this to work.
    S = reverse(state.closest_vertices) #Replaced sortperm with this
    for w in S  # w is the farthest vertex from si
        for v in P[w]  # get the predecessors of w
            if v > 0
                δ[v] +=  δ[w] + 1 # increment sp of pred
            end
        end
        δ[w] *= length(P[w]) # adjust the # of sps of vertex
        if w != si
            stress[w] += δ[w]
        end
    end
    return nothing
end


text:
using Images
using Distributions
function ReLU(x)
    if(x>0)
        return x
        else
            return 0.0
    end
end
function softmax(u::Array,output_vector_size)
    z=similar(u)
    u_sum=sum(exp.(u-maximum(u)))
    for k in 1:output_vector_size
        z[k]=exp(u[k]-maximum(u))/u_sum
    end
    return z
end

#畳み込み層
function convolution(l)
    x=Int(0)
    y=Int(0)
    for si in 0:S:datay-H
        y+=1
        x=0
        for sj in 0:S:datax-H
            x+=1
            for k in 1:K
                for i in 1:H
                    for j in 1:H
                        u[l,y,x]+=imgdata[k,i+si,j+sj]*h[i,j]+b[l,y,x]
                    end
                end
            end
            z[l,y,x]=ReLU(u[l,y,x])
        end
    end
    return z
end
#チャネルなし　2回め以降
function convolution(l,u::Array, z::Array)
    x=Int(0)
    y=Int(0)
    for si in 0:S:datay-H
        y+=1
        x=0
        for sj in 0:S:datax-H
            x+=1
                for i in 1:H
                    for j in 1:H
                        u[l,y,x]+=z[l-1,i+si,j+sj]*h[i,j]+b[l,y,x]
                    end
                end
            z[l,y,x]=ReLU(u[l,y,x])
        end
    end
    return z
end

function nuwral_net(l,z::Array)
    u=similar(z[l,:])
    u=w[l,:,:]*z[l-1,:]
    z[l,:]=ReLU.(u)
    return z
end

function nuwral_net(l,z::Array,output_vector_size)
    u=similar(z[l,1:output_vector_size])
    u=w[l,1:output_vector_size,:]*z[l-1,:]
    z[l,1:output_vector_size]=softmax(u,output_vector_size)
    return z
end

function maxpooling(l,z::Array)
    x=Int(0)
    y=Int(0)
    for si in 0:S:datay-P
        y+=1
        x=0
        for sj in 0:S:datax-P
            x+=1
            z[l,y,x]=maximum(z[l-1,1+si:P+si,1+sj:P+sj])
        end
    end
    return z

end

function dReLU(x)
    if(x>0)
        return 1.0
        else
        return 0.0
    end
end

function cnn(imgdata,datay,datax)

end

function main()
	m=Normal(0,1)
	rand(m,4,4)
	#チャネル数
	K=3 #RGB
	train=[0.0 for i in 1:40]
	train[38]=1.0
	data=load("phot/1berserk.jpg")
	datax=3402 #画像サイズ
	datay=3402
	const Layernum=21
	imgdata=channelview(data) #[RGB,4608,3456]
end


text:
"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{T, ThermalStandardUCOutages},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(T, sys)

    if !PSI.validate_available_devices(T, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReactivePowerVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.OnVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.StartVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.StopVariable,
        devices,
        ThermalStandardUCOutages(),
    )

    # Initial Conditions
    PSI.initial_conditions!(optimization_container, devices, ThermalStandardUCOutages())

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.commitment_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.ramp_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.time_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    return
end

"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{T, ThermalStandardUCOutages},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    if !PSI.validate_available_devices(T, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.OnVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.StartVariable,
        devices,
        ThermalStandardUCOutages(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.StopVariable,
        devices,
        ThermalStandardUCOutages(),
    )

    # Initial Conditions
    PSI.initial_conditions!(optimization_container, devices, ThermalStandardUCOutages())

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.commitment_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.ramp_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.time_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalDispatchOutages},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerVariable, devices, ThermalDispatchOutages())
    add_variables!(optimization_container, ReactivePowerVariable, devices, ThermalDispatchOutages())

    # Initial Conditions

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Cost Function
    cost_function!(optimization_container, devices, model, S, get_feedforward(model))

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalDispatchOutages},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerVariable, devices, ThermalDispatchOutages())

    # Initial Conditions

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Cost Function
    cost_function!(optimization_container, devices, model, S, get_feedforward(model))

    return
end


text:
using QuizQuestions
using Test

@testset "QuizQuestions.jl" begin
    @test numericq(1).tol == 0
    @test numericq(1, 1e-3; hint="hint").hint == "hint"

    @test stringq(r"abc").re == r"abc"
    @test_throws MethodError stringq("abc") # need regular expression

    @test radioq(["one","two", "three"], 1, keep_order=true).answer == 1
    @test radioq(1:3, 1, keep_order=true).answer == 1

    @test multiq(("one","two","three"), (2, 3), label="Primes?").label == "Primes?"
end


text:
######### Stan program example  ###########

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  bernoullidata = [
    Dict("N" => 1, "y" => [0]),
    Dict("N" => 1, "y" => [0]),
    Dict("N" => 1, "y" => [0]),
    Dict("N" => 1, "y" => [0]),
  ]

  monitor = ["theta", "lp__", "accept_stat__"]

  global stanmodel, rc, sim1, sim
  stanmodel = Stanmodel(num_samples=1200, thin=2, name="bernoulli", model=bernoullimodel);

  println("\nStanmodel that will be used:")
  stanmodel |> display
  println()
  println("Input observed data dictionary:")
  bernoullidata |> display
  println()

  rc, sim1 = stan(stanmodel, bernoullidata, ProjDir, diagnostics=false, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    ## Subset Sampler Output to variables suitable for describe().
    sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
    describe(sim)
    println()

    ## Brooks, Gelman and Rubin Convergence Diagnostic
    try
      gelmandiag(sim, mpsrf=true, transform=true) |> display
    catch e
      #println(e)
      gelmandiag(sim, mpsrf=false, transform=true) |> display
    end
    println()

    ## Geweke Convergence Diagnostic
    gewekediag(sim) |> display

    ## Highest Posterior Density Intervals
    hpd(sim) |> display

    ## Cross-Correlations
    cor(sim) |> display

    ## Lag-Autocorrelations
    autocor(sim) |> display

    ## Deviance Information Criterion
    #dic(sim) |> display

    ## Plotting
    p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
    draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
    draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

  	# Pairwise contour plots
  	simmod = ModelChains(sim,Mamba.Model())
  	p = plot(simmod, :contour)
  	draw(p, nrow=2, ncol=2, filename="$(stanmodel.name)-contourplot.svg")
  	draw(p, nrow=2, ncol=2, filename="$(stanmodel.name)-contourplot.pdf", fmt=:pdf)

  end # rc == 0
end # cd


text:
allbindings(pat, bs) =
  isbinding(pat) || (isslurp(pat) && pat ≠ :__) ? push!(bs, bname(pat)) :
  isa(pat, TypeBind) ? push!(bs, pat.name) :
  isa(pat, OrBind) ? (allbindings(pat.pat1, bs); allbindings(pat.pat2, bs)) :
  istb(pat) ? push!(bs, tbname(pat)) :
  isexpr(pat, :$) ? bs :
  isa(pat, Expr) ? map(pat -> allbindings(pat, bs), [pat.head, pat.args...]) :
  bs

allbindings(pat) = (bs = Any[]; allbindings(pat, bs); bs)

function bindinglet(bs, body)
  ex = :(let $(esc(:env)) = env
           $body
         end)
  for b in bs
    push!(ex.args, :($(esc(b)) = get(env, $(Expr(:quote, b)), nothing)))
  end
  return ex
end

function makeclause(line, els = nothing)
  env = trymatch(:(pat_ => yes_), line)
  env == nothing && error("Invalid match clause $line")
  pat, yes = env[:pat], env[:yes]
  bs = allbindings(pat)
  pat = subtb(subor(pat))
  quote
    env = trymatch($(Expr(:quote, pat)), ex)
    if env != nothing
      $(bindinglet(bs, esc(yes)))
    else
      $els
    end
  end
end

macro match(ex, lines)
  @assert isexpr(lines, :block)
  result = quote
    ex = $(esc(ex))
  end
  body = nothing
  for line in reverse(rmlines(lines).args)
    isline(result) && push!(result, line)
    body = makeclause(line, body)
  end
  push!(result.args, body)
  return result
end

macro capture(ex, pat)
  bs = allbindings(pat)
  pat = subtb(subor(pat))
  quote
    $([:($(esc(b)) = nothing) for b in bs]...)
    env = trymatch($(Expr(:quote, pat)), $(esc(ex)))
    if env == nothing
      false
    else
      $([:($(esc(b)) = get(env, $(Expr(:quote, b)), nothing)) for b in bs]...)
      true
    end
  end
end


text:
export BrickletCurrent25CurrentCallbackThreshold
struct BrickletCurrent25CurrentCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletCurrent25AnalogValueCallbackThreshold
struct BrickletCurrent25AnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletCurrent25Identity
struct BrickletCurrent25Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletCurrent25
"""
Measures AC/DC current between -25A and +25A
"""
mutable struct BrickletCurrent25 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletCurrent25(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_current25")
        deviceInternal = package.BrickletCurrent25(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_current
"""
    $(SIGNATURES)

Returns the current of the sensor.

If you want to get the current periodically, it is recommended to use the
:cb:`Current` callback and set the period with
:func:`Set Current Callback Period`.
"""
function get_current(device::BrickletCurrent25)
    return device.deviceInternal.get_current()
end

export calibrate
"""
    $(SIGNATURES)

Calibrates the 0 value of the sensor. You have to call this function
when there is no current present.

The zero point of the current sensor
is depending on the exact properties of the analog-to-digital converter,
the length of the Bricklet cable and the temperature. Thus, if you change
the Brick or the environment in which the Bricklet is used, you might
have to recalibrate.

The resulting calibration will be saved on the EEPROM of the Current
Bricklet.
"""
function calibrate(device::BrickletCurrent25)
    device.deviceInternal.calibrate()
end

export is_over_current
"""
    $(SIGNATURES)

Returns *true* if more than 25A were measured.

.. note::
 To reset this value you have to power cycle the Bricklet.
"""
function is_over_current(device::BrickletCurrent25)
    return device.deviceInternal.is_over_current()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

.. note::
 The value returned by :func:`Get Current` is averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletCurrent25)
    return device.deviceInternal.get_analog_value()
end

export set_current_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Current` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Current` callback is only triggered if the current has changed since
the last triggering.
"""
function set_current_callback_period(device::BrickletCurrent25, period)
    device.deviceInternal.set_current_callback_period(period)
end

export get_current_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Current Callback Period`.
"""
function get_current_callback_period(device::BrickletCurrent25)
    return device.deviceInternal.get_current_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletCurrent25, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletCurrent25)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_current_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Current Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the current is *outside* the min and max values"
 "'i'",    "Callback is triggered when the current is *inside* the min and max values"
 "'<'",    "Callback is triggered when the current is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the current is greater than the min value (max is ignored)"
"""
function set_current_callback_threshold(device::BrickletCurrent25, option, min, max)
    device.deviceInternal.set_current_callback_threshold(option, min, max)
end

export get_current_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Current Callback Threshold`.
"""
function get_current_callback_threshold(device::BrickletCurrent25)
    return device.deviceInternal.get_current_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletCurrent25, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletCurrent25)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Current Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Current Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletCurrent25, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletCurrent25)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletCurrent25)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletCurrent25, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
################
# FieldStorage #
################

abstract type FieldStorage{D} end

# Mutable #
#---------#

mutable struct Mutable{D} <: FieldStorage{D}
    data::D
    Mutable{D}(data) where D = new{D}(data)
end

@inline Mutable(data::D) where {D} = Mutable{D}(data)

@inline Base.getindex(x::Mutable) = x.data

@inline Base.setindex!(x::Mutable, y) = (x.data = y)

# Immutable #
#-----------#

struct Immutable{D} <: FieldStorage{D}
    data::D
end

@inline Base.getindex(x::Immutable) = x.data

@inline Base.setindex!(x::Immutable, y) = error("cannot mutate immutable field")

###########################
# metadata type hierarchy #
###########################
#=
TODO: clean up/formalize these notes

The optimal dualtype specification for a given non-leaftype is the union of the dualtypes of
its subtypes, and thus can be implemented as:

```
function dualtype(::Type{C}, ::Type{T}) where {C,T}
    if isempty(subtypes(T))
        return metatype(C, T)
    else
        return Union{map(X -> dualtype(C, X), subtypes(T))...}
    end
end
```

However, this is an expensive implementation of this operation. Thus, our actual `dualtype`
implementation returns a still-correct, but an extremely pessimistic dualtype with the
benefit that dualtype computation is very fast. If, in the future, `dualtype` is
parameterized on world age, then we can call `subtypes` at compile time, and compute a more
optimally bounded dualtype.
=#

const MetaTree = Union{NamedTuple,Array,Unused}

struct MetaNode{D,T<:MetaTree}
    data::Union{D,Unused}
    tree::T
    MetaNode{D}(data, tree::T) where {D,T} = new{D,T}(data, tree)
end

@inline function metanodetype(::Type{C}, ::Type{T}) where {C<:Context,T}
    if isconcretetype(T)
        return Union{MetaNode{metatype(C, T),subtreetype(C, T)},Unused}
    end
    return Union{MetaNode,Unused}
end

@generated function subtreetype(::Type{C}, ::Type{T}) where {C,T}
    if !(isconcretetype(T))
        return :(error("cannot call subtreetype on non-concrete type ", $T))
    end
    if T <: Array
        body = :(Array{metanodetype(C, $(eltype(T))),$(ndims(T))})
    elseif fieldcount(T) == 0
        body = :(Unused)
    else
        S = isimmutable(T) ? :Immutable : :Mutable
        fnames = Expr(:tuple, map(Meta.quot, fieldnames(T))...)
        ftypes = map(F -> :($S{metanodetype(C, $F)}), T.types)
        body = :(NamedTuple{$fnames,Tuple{$(ftypes...)}})
    end
    return quote
        $(Expr(:meta, :inline))
        $body
    end
end

@generated function initmetanode(::C, ::V, metadata::D) where {C<:Context,V,D}
    if V <: Array
        return quote
            $(Expr(:meta, :inline))
            MetaNode{metatype(C, V)}(metadata, similar(value, metanodetype(C, eltype(V))))
        end
    elseif fieldcount(V) == 0
        return quote
            $(Expr(:meta, :inline))
            MetaNode{metatype(C, V)}(metadata, unused)
        end
    else
        S = isimmutable(V) ? :Immutable : :Mutable
        fnames = fieldnames(V)
        fdatas = map(F -> :($S{metanodetype(C, $F)}(unused)), V.types)
        fields = Expr(:tuple)
        for i in 1:fieldcount(V)
            push!(fields.args, :($(fnames[i]) = $(fdatas[i])))
        end
        return quote
            $(Expr(:meta, :inline))
            MetaNode{metatype(C, V)}(metadata, $fields)
        end
    end
end

#######
# Box #
#######

#=
Here, `U` is the innermost, "underlying" type of the value being wrapped. This parameter is
precomputed so that Cassette can directly dispatch on it in the signatures generated for
contextual primitives.
=#
struct Box{C,U,V,D,T}
    context::C
    value::V
    meta::MetaNode{D,T}
    function Box(context::C, value::V, meta::MetaNode{D,T}) where {C<:Context,V,D,T}
        U = _underlying_type(V)
        return new{C,U,V,D,T}(context, value, meta)
    end
end

function Box(context::C, value, metadata = unused) where {C<:Context}
    return Box(context, value, initmetanode(context, value, metadata))
end

_underlying_type(::Type{V}) where {V} = V
_underlying_type(::Type{<:Box{<:Context,U}}) where {U} = U

#######
# API #
#######

@inline isboxed(::Context, ::Any) = false
@inline isboxed(::C,       ::Box{C}) where {C<:Context} = true
@inline isboxed(::Type{C}, ::DataType) where {C<:Context} = false
@inline isboxed(::Type{C}, ::Type{<:Box{C}}) where {C<:Context} = true

@inline unbox(ctx::Context, x) = isboxed(ctx, x) ? x.value : x
@inline unbox(::Type{C},    T::DataType) where {C<:Context,X} = T
@inline unbox(::Type{C},    ::Type{<:Box{C,U,V}}) where {C<:Context,U,V} = V

@inline meta(::C,       x::Box{C}) where {C<:Context} = x.meta.data
@inline meta(::Type{C}, ::Type{<:Box{C,U,V,D}}) where {C<:Context,U,V,D} = D

@inline hasmeta(::Context, ::Any) = false
@inline hasmeta(ctx::C,    x::Box{C}) where {C<:Context} = isa(meta(ctx, x), Unused)

@inline metatype(::Type{<:Context}, ::DataType) = Unused

@generated function unboxcall(ctx::Context, f, args...)
    N = fieldcount(typeof(args))
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, :f, [:(unbox(ctx, args[$i])) for i in 1:N]...))
    end
end

#########
# `new` #
#########

@generated function _newbox(ctx::C, ::Type{T}, args...) where {C<:Context,T}
    unboxed_args = [:(unbox(ctx, args[$i])) for i in 1:length(args)]
    if length(args) == fieldcount(T)
        fields = Expr(:tuple)
        fnames = fieldnames(T)
        boxcount = 0
        S = isimmutable(T) ? :Immutable : :Mutable
        for i in 1:fieldcount(T)
            arg = args[i]
            fname = fnames[i]
            ftype = :(metanodetype(C, $(T.types[i])))
            if arg <: Box{C}
                boxcount += 1
                argmetanode = :(args[$i].meta)
            else
                argmetanode = :unused
            end
            push!(fields.args, :($fname = $S{$ftype}($argmetanode)))
        end
        if !(boxcount == 0 && isimmutable(T))
            return quote
                $(Expr(:meta, :inline))
                Box(ctx, _new(T, $(unboxed_args...)), MetaNode{metatype(C, T)}(unused, $fields))
            end
        end
    end
    return quote
        $(Expr(:meta, :inline))
        _new(T, $(unboxed_args...))
    end
end

@generated function _new(::Type{T}, args...) where {T}
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:new, T, [:(args[$i]) for i in 1:nfields(args)]...))
    end
end

##########################
# `getfield`/`setfield!` #
##########################

@inline _getfield(x, name) = getfield(x, name)
@inline _getfield(x::Box, name) = Box(x.context, getfield(x.value, name), getfield(x.meta.tree, name)[])

@inline _setfield!(x, name, y) = setfield!(x, name, y)
@inline _setfield!(x::Box, name, y) = _setfield!(x, name, y, y, unused)
@inline _setfield!(x::Box{C}, name, y::Box{C}) where {C} = _setfield!(x, name, y, y.value, y.meta)

@inline function _setfield!(x::Box, name, y, y_value, y_meta)
    setfield!(x.value, name, y_value)
    setindex!(getfield(x.meta.tree, name), y_meta)
    return y
end

#########################
# `arrayref`/`arrayset` #
#########################

@inline _arrayref(check, x, i) = arrayref(check, x, i)
@inline _arrayref(check, x::Box, i) = Box(x.context, arrayref(check, x.value, i), arrayref(check, x.meta.tree, i))

@inline _arrayset(check, x, y, i) = arrayset(check, x, y, i)
@inline _arrayset(check, x::Box, y, i) = _arrayset(check, x, y, unused, i)
@inline _arrayset(check, x::Box{C}, y::Box{C}, i) where {C} = _arrayset(check, x, y.value, y.meta, i)

@inline function _arrayset(check, x::Box, y_value, y_meta, i)
    arrayset(check, x.value, y_value, i)
    arrayset(check, x.meta, y_meta, i)
    return x
end


text:
# function ⋃(Sets...)
#     return reduce(∪, Sets)
# end
const ⋃ = union

function ∂(Simplex)
    # if !(Simplex |> issorted) @warn "Simplex $Simplex is not sorted" end
    faces = [setdiff(Simplex, v) for v in reverse(Simplex)]
    return faces
end

function faces(Simplex::Set)
    return Set([setdiff(Simplex, v) for v in Simplex])
end
faces(Set([3,5,3,2]))
∂([1,2,3])
∂([1,3,2])


text:
#thermo module, meta/package to reference all other subpackages
module Thermo
include("UniversalConstants.jl") #provides constants
include("Compounds.jl") #provides an structure to store and retrieve compound properties
end


text:
##Recall that rho(p,q,z) = exp(-beta*H)*f(z)

"""
      Random initial vector with components (q,p,zeta)
                """
function initialcondition(beta::Float64, Q::Float64)
    sigma = 1./sqrt(beta)  #Standard deviation
    q,p = rand(Normal(0.0, sigma), 2)
    z = rand(Logistic(0.0,Q))
    return [q,p,z]
end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg
exampledir = joinpath(@__DIR__, "..", "examples")
Pkg.activate(exampledir)
Pkg.add("Literate")
include(joinpath(exampledir, "generate_notebooks.jl"))


text:
export AbstractApproximator, AbstractQApproximator, AbstractVApproximator,
       update!

"""
    AbstractApproximator

`AbstractApproximator` is a supertype of a collection of approximators(including Tabular, Linear, DNN etc.)
"""
abstract type AbstractApproximator end

"""
    AbstractVApproximator{Ts} <: AbstractApproximator

A supertype of state value approximators with the state of type `Ts`.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `V(s)` | `V`, an instance of `AbstractVApproximator`, must be a functional object which takes in state `s` and return the estimated value of that state |
| `update!(V, s, e)` | Update the value approximator of state `s` by a difference of `e` |
"""
abstract type AbstractVApproximator{Ts} <: AbstractApproximator end

"""
    AbstractQApproximator{Ts, Ta} <: AbstractApproximator

A supertype of action value approximators with the state of type `Ts` and action of type `Ta`.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `Q(s, a)` | `Q`, an instance of `AbstractQApproximator` must be a functional object which takes in state `s` and action `a`. The estimated value of that (state, action) will be returned. |
| `Q(s)` | Return the estimations of all actions at state `s` |
| `update!(Q, s, a, e)` | Update the approximator `Q` of state `s` and action `a` with difference `e` |
| **Optional Methods** | |
| `Q(s, Val(:max))` | Return the maximum value of `Q(s)` |
| `Q(s, Val(:argmax))` | Return the argmax of `Q(s)` (randomly break ties) |
| `update!(Q, s, errors)` | Update the approximator `Q` of state `s` with the difference `errors` of each action |
"""
abstract type AbstractQApproximator{Ts, Ta} <: AbstractApproximator end


text:
# main file for computing dRdq transposed products

import PDESolver.evalResidual_revq

function evalResidual_revq(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                     opts::Dict, t::Number=0.0)

  time = eqn.params.time
  eqn.params.t = t  # record t to params
  myrank = mesh.myrank

  # parallel data exchange must be started before this function is called
  # (both q and res_bar)

  # Forward sweep
  time.t_dataprep += @elapsed dataPrep_for_revq(mesh, sbp, eqn, opts)

  time.t_volume += @elapsed if opts["addVolumeIntegrals"]
    evalVolumeIntegrals_revq(mesh, sbp, eqn, opts)
  end

  if opts["use_GLS"]
    error("GLS not supported for revq product")
  end

  time.t_bndry += @elapsed if opts["addBoundaryIntegrals"]
    evalBoundaryIntegrals_revq(mesh, sbp, eqn, opts)
  end

  time.t_stab += @elapsed if opts["addStabilization"]
     addStabilization_revq(mesh, sbp, eqn, opts)
  end

  time.t_face += @elapsed if mesh.isDG && opts["addFaceIntegrals"]
    evalFaceIntegrals_revq(mesh, sbp, eqn, opts)
  end

  time.t_source += @elapsed evalSourceTerm_revq(mesh, sbp, eqn, opts)

  time.t_sharedface += @elapsed if mesh.commsize > 1
    evalSharedFaceIntegrals_revq(mesh, sbp, eqn, opts)
  end

  time.t_shock += @elapsed if opts["addShockCapturing"]
    evalShockCapturing_revq(mesh, sbp, eqn, opts)
  end

  # apply inverse mass matrix to eqn.res, necessary for CN
  if opts["use_Minv"]
    error("use_Minv not supported for revq product")
    applyMassMatrixInverse3D(mesh, sbp, eqn, opts, eqn.res)
  end

  time.t_dataprep += @elapsed dataPrep_revq(mesh, sbp, eqn, opts)

  return nothing

end  # end evalResidual

"""
  Dataprep-type function that should be called at the beginning of a reverse
  mode product, rather than the regular `dataPrep` function.
"""
function dataPrep_for_revq(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  # apply filtering to input
  if eqn.params.use_filter
    applyFilter(mesh, sbp, eqn, eqn.q, opts)
  end

  getAuxVars(mesh, eqn)

  if opts["use_staggered_grid"]
    error("staggered grid not supported for revq product")
    aux_vars = zeros(Tres, 1, mesh.mesh2.numNodesPerElement)
    for i=1:mesh.numEl
      qs = ro_sview(eqn.q, :, :, i)
      qf = sview(eqn.q_flux, :, :, i)
      interpolateElementStaggered(eqn.params, mesh, qs, aux_vars, qf)
    end
  end

  # only eqn.q_bndry is required, all other methods are no-precompute
  if mesh.isDG
    interpolateBoundary(mesh, sbp, eqn, opts, eqn.q, eqn.q_bndry, eqn.aux_vars_bndry)
  end

  if eqn.params.use_edgestab
    stabscale(mesh, sbp, eqn)
  end

  return nothing
end

"""

Reverse mode of dataPrep w.r.t mesh metrics

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""
function dataPrep_revq(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  # nothing to do here

  return nothing
end

"""

Reverse mode of evalVolumeIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""
function evalVolumeIntegrals_revq(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh,  Tsol, Tres, Tdim}

  integral_type = opts["volume_integral_type"]
  if integral_type == 1
    if opts["Q_transpose"] == true
      calcVolumeIntegrals_nopre_revq(mesh, sbp, eqn, opts)
    else
      error("Q_transpose = false not supported for revq product")
    end  # end if
  elseif integral_type == 2
    calcVolumeIntegralsSplitForm_revq(mesh, sbp, eqn, opts, eqn.volume_flux_func,
                                 eqn.volume_flux_func_revq)
  else
    throw(ErrorException("Unsupported volume integral type = $integral_type"))
  end

  return nothing
end  # end evalVolumeIntegrals

"""
### EulerEquationMod.evalBoundaryIntegrals_revq

Reverse mode of evalBoundaryIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""
function evalBoundaryIntegrals_revq(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}

  getBCFluxes_revq(mesh, sbp, eqn, opts)

  return nothing

end  # end evalBoundaryIntegrals

function addStabilization_revq(mesh::AbstractMesh, sbp::AbstractOperator,
                               eqn::EulerData, opts)

  if opts["use_lps"]
    applyLPStab_revq(mesh, sbp, eqn, opts)
  end

  return nothing
end

"""

Reverse mode of evalFaceIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""
function evalFaceIntegrals_revq(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol}, opts) where {Tmsh, Tsol}

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1
    calcFaceIntegral_nopre_revq(mesh, sbp, eqn, opts, eqn.flux_func_revq,
                                mesh.interfaces)

  elseif face_integral_type == 2

    getFaceElementIntegral_revq(mesh, sbp, eqn, eqn.face_element_integral_func,
                           eqn.flux_func_revq, mesh.sbpface, mesh.interfaces)

  else
    throw(ErrorException("Unsupported face integral type = $face_integral_type"))
  end

  # do some output here?
  return nothing
end

"""

  Reverse mode evalSharedFaceIntegrals with respect to q.  Note that this
  uses the parallel approach of [`finishExchangeData_rev2`](@ref), so the
  functions called here are not precisely the reverse mode of the primal
  functions.

"""
function evalSharedFaceIntegrals_revq(mesh::AbstractDGMesh, sbp, eqn, opts)

  if mesh.npeers == 0
    return nothing
  end

  if getParallelData(eqn.shared_data) != PARALLEL_DATA_ELEMENT ||
     getParallelData(eqn.shared_data_bar) != PARALLEL_DATA_ELEMENT

    error("""parallel data setting must be PARALLEL_DATA_ELEMENT """)
  end

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1
    finishExchangeData_rev2(mesh, sbp, eqn, opts, eqn.shared_data,
                      eqn.shared_data_bar, calcSharedFaceIntegrals_element_revq)
  elseif face_integral_type == 2
    finishExchangeData_rev2(mesh, sbp, eqn, opts, eqn.shared_data,
                      eqn.shared_data_bar, calcSharedFaceElementIntegrals_element_revq)
  else
    throw(ErrorException("unsupported face integral type = $face_integral_type"))
  end

  return nothing
end

"""

Reverse mode of evalSourceTerm w.r.t mesh metrics.

*This function has not been written properly. The developer must update this
documentation as and when this code is developed*

"""
function evalSourceTerm_revq(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # placeholder for multiple source term functionality (similar to how
  # boundary conditions are done)
  if opts["use_src_term"]
    error("source terms not supported for revq product")
  end

  return nothing
end  # end function

"""

  Reverse mode of `evalShockCapturing`

  **Inputs**:

   * mesh : Abstract mesh type
   * sbp  : Summation-by-parts operator
   * eqn  : Euler equation object
   * opts : options dictonary
"""
function evalShockCapturing_revq(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # currently there is only one choice for the shock capturing scheme.
  # If there are more, need something more sophisiticated to choose.
  # Perhaps add an abstract typed field to eqn containing the structs

  capture = eqn.shock_capturing
  applyShockCapturing_revq(mesh, sbp, eqn, opts, capture)

  return nothing
end


text: calc_n_20(reactor, R_0, B_0) = K_n(reactor) / ( B_0 * sqrt(R_0^5) )


text:
"""
    semi_elementary_monomial(x,y)

First, given ``n`` variables, we say ``e_i`` is the elementary symmetric polynomial or order i with coefficient 1.\n
For example, when ``n=3``,

``e_1 = x+y+z``

``e_2 = xy+xz+yz``

``e_3 = xyz``

When ``n=4``,

``e_1 = w+x+y+z``

``e_2 = wx+wy+wz+xy+xz+yz``

``e_3 = wxy+wxz+wyz+xyz``

``e_4 = wxyz``

and so on

From the users perspective, usually the ``n`` is obvious, so it is omitted.

Then, ``S(a_1,a_2...a_n)`` is the symmetric polynomial in ``n`` variables with powers ``0 \\leq a_1 \\leq a_2...\\leq a_n`` and coefficient 1\n
For example,

``S(0,1,2) = x^2y + x^2z + y^2z``

``S(0,0,1,2) = w^2x + w^2y + w^2z + x^2y + x^2z + y^2z``

``S(0,0,0,0,5) = v^5 + w^5 + x^5 + y^5 + z^5``

and so on

A semi elementary monomial is the product of\n
    1. A non elementary symmetric polynomial with power 1
    2. An elementary monomial

For example, the follow are semi elementary

``S(0,0,3)e_1^3e_2^1e_3^0``

``S(0,0,1,2)e_1^3e_2^0e_3^1e_4^2``

These would be represented as two tuples,

``S(0,0,3)e_1^3e_2^1e_3^0           => (0,0,3),(3,1,0)``

``S(0,0,1,2)e_1^3e_2^0e_3^1e_4^2  => (0,0,1,2),(3,0,1,2)``

This is where the name semi comes from, in that nearly all terms in the monomial are elementary except 1.

When the user wants to create a semi elementary monomial, he or she needs to provide two tuples with non-negative integers.
The first tuple will automatically be sorted by the constructor, the second will be provided as it is.

    semi_elementary_monomial((0,2,1),(3,4,5))

returns the same object as

    semi_elementary_monomial((0,1,2),(3,4,5))

Besides sorting, some trivial simplifications will also be made.
First, ``e_n`` terms will be automatically pulled out, for example

    semi_elementary_monomial((2,3,4),(0,0,0))

will return

    semi_elementary_monomial((0,1,2),(0,0,2))

Which is basically saying

``(x^2y^3z^4 + x^2z^3y^4 + y^2x^3z^4 + y^2z^3x^4 + z^2x^3y^4 + z^2x^4y^3) = (x^2y + x^2z + y^2z) * (xyz)^2``

Second, if the first term is already elementary, it will be merged into the second term, for example

    semi_elementary_monomial((0,1,1),(0,0,0))

will return

    semi_elementary_monomial((0,0,0),(0,1,0))
"""
struct semi_elementary_monomial{N}
    sp_term::NTuple{N,Int64}
    powers::NTuple{N,Int64}

    function semi_elementary_monomial(x::NTuple{N,Int64},y::NTuple{N,Int64}) where {N}
        t = TupleTools.sort(x)
        if t[end] - t[1] == 1
            c = sum(t) - N * t[1]
            tmp2 = ntuple(N) do i
                    if i==c
                        return y[i] + 1
                    elseif i==N
                        return y[i] + t[1]
                    else
                        return y[i]
                    end
                end
            tmp1 = ntuple(i->0,N)
        elseif t[1] > 0
            tmp2 = ntuple(i->i==N ? y[i]+t[1] : y[i],N)
            tmp1 = t .- t[1]
        else
            tmp1,tmp2 = t,y
        end
        return new{N}(tmp1,tmp2)
    end
end

dim(x::semi_elementary_monomial{N}) where {N} = N
is_elementary(x::semi_elementary_monomial) = x.sp_term[end]==0
Base.:(==)(x::semi_elementary_monomial,y::semi_elementary_monomial) = (x.sp_term==y.sp_term) && (x.powers==y.powers)
function Base.isless(x::semi_elementary_monomial{N},y::semi_elementary_monomial{N}) where {N}
    for i=N:-1:1
        x.sp_term[i] != y.sp_term[i] && return x.sp_term[i] < y.sp_term[i]
    end
    s1,s2 = sum(x.powers),sum(y.powers)
    return s1!=s2 ? s1 < s2 : x.powers < y.powers
    #return x.powers < y.powers
end

function to_string(x::semi_elementary_monomial)
    if !is_elementary(x)
        head = "S"*string(x.sp_term)
    else
        head = ""
    end
    for (i,n) in enumerate(x.powers)
        if n == 0
            continue
        elseif n == 1
            head *= "e$i "
        else
            head *= "e$i^$n "
        end
    end
    return head
end
Base.show(io::IO, ::MIME{Symbol("text/plain")}, x::semi_elementary_monomial) = join(io,to_string(x))

max_stars_repo_path: src/semi elementary monomial.jl

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_path: src/semi elementary monomial.jl

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_path: src/semi elementary monomial.jl

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
####################################
### Data preparation's functions ###
####################################

"""
The `colnames` dictionary is used to create the `SimData` object. The role of this dictionary is to ensure that the included files contains proper data.

* `:DAs` : a list of columns  in the csv file describing DAs:
    * `:DA_ID` : ID of a DA
    * `:LONGITUDE` : longitude of a DA's centroid
    * `:LATITUDE` : latitude of a DA's centroid

* `demo_stats` : a list of columns which must be included in the csv file with the demographic data used in simulation. **Warning:** Columns must be described in a separate

* `:flows`  : a list of columns  in the csv file with the flows between DAs.
    * `:DA_home` : ID of start DA
    * `:DA_work` : ID of finish DA
    * `:Flow_Volume` : flow volume
"""
 const colnames = Dict(
 :DAs => [:DA_ID, :LONGITUDE, :LATITUDE],
 :flows  => [:DA_home, :DA_work, :FlowVolume]
 )

#get features
"""
Find the nearest node for every feature

**Arguments**
* `frame` : data frame with all features
* `nodes` : a list of nodes from .osm file
* `bounds` : map bounds
"""
function features_to_nodes(frame::DataFrames.DataFrame,
                            nodes::Dict{Int,OpenStreetMapX.ENU},
                            bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA})
    coordinates = values(nodes)
    features = Dict{Int,Tuple{String,String}}()
    sizehint!(features,DataFrames.nrow(frame))
    for i = 1:DataFrames.nrow(frame)
        loc = OpenStreetMapX.ENU(OpenStreetMapX.LLA(frame[:LATITUDE][i], frame[:LONGITUDE][i]), OpenStreetMapX.center(bounds))
        if !OpenStreetMapX.inbounds(loc,OpenStreetMapX.ENU(bounds, OpenStreetMapX.center(bounds)))
            loc = OpenStreetMapX.boundary_point(loc, OpenStreetMapX.ENU(center(bounds), center(bounds)), OpenStreetMapX.ENU(bounds, center(bounds)))
        end
        node = OpenStreetMapX.add_new_node!(nodes,loc)
        features[node] = (frame[:CATEGORY][i], frame[:NAME][i])
    end
    return features
end

"""
Prepare features data objects.

**Arguments**
* `datapath` : path with osm file
* `mapfile` : a `map file name
"""
function get_features_data(datapath::String, mapfile::String, map_data::OpenStreetMapX.MapData ;use_cache::Bool = true)
	cachefile = joinpath(datapath,"features.cache")
	if use_cache && isfile(cachefile)
		f=open(cachefile,"r");
		res=Serialization.deserialize(f);
		close(f);
		@info "Read features data from cache $cachefile"
	else
		mapdata = OpenStreetMapX.parseOSM(joinpath(datapath,mapfile))
		OpenStreetMapX.crop!(mapdata,crop_relations = false)
		feature_classes = Dict("amenity" => 3, "leisure" => 13,"shop" => 23, "tourism" => 25)
		features = OpenStreetMapX.filter_features(mapdata.features, levels = Set(values(feature_classes)))
		feature_to_intersections = OpenStreetMapX.features_to_graph(merge(OpenStreetMapX.ENU(mapdata.nodes,map_data.bounds),map_data.nodes), features, map_data.network)
		res = (features, feature_classes, feature_to_intersections)
		if use_cache
				f=open(cachefile,"w");
				Serialization.serialize(f,res);
				@info "Saved features data to cache $cachefile"
				close(f);
		end
	end
    return res
end

#finding closest network node for each DA
"""
Read a csv file with informations about DAs and then for each DA find the nearest route network node.

**Arguments**
* `datapath` : path with csv files
* `filename` : name of csv file
* `colnames` : an array of columns which must be included in each file
* `map_data` : a `OpenStreetMapX.MapData` object
"""
function DAs_to_nodes(datapath::String, filename::String,
                    colnames::Array{Symbol,1},
                    map_data::OpenStreetMapX.MapData)
    DAframe = Nanocsv.read_csv(joinpath(datapath,filename))
	  #DataFrame(CSVFiles.load(joinpath(datapath,filename)))
    if !all(in(col, DataFrames.names(DAframe)) for col in colnames)
        error("Wrong column names! Data Frame should contain $(String.(colnames).*" "... ) columns!")
    end
    DAs_to_nodes = Dict{Int,Int}()
    sizehint!(DAs_to_nodes,DataFrames.nrow(DAframe))
    for i = 1:DataFrames.nrow(DAframe)
        coords = OpenStreetMapX.ENU(OpenStreetMapX.LLA(DAframe[:LATITUDE][i], DAframe[:LONGITUDE][i]), OpenStreetMapX.center(map_data.bounds))
        DAs_to_nodes[DAframe[:DA_ID][i]] = OpenStreetMapX.nearest_node(map_data.nodes,coords,map_data.network)
    end
    return DAs_to_nodes
end

"""
Include csv file with flow data and create a flow matrix

**Arguments**
* `datapath` : path with csv files
* `filename` : name of csv file
* `colnames` : an array of columns which must be included in each file
"""
function get_flow_data(datapath::String, filename::String, colnames::Array{Symbol,1})
    flows = Nanocsv.read_csv(joinpath(datapath,filename))

    if !all(in(col, DataFrames.names(flows)) for col in colnames)
        error("Wrong column names! Data Frame should contain $(String.(colnames).*" "... ) columns!")
    end
    flows = flows[colnames]
    vals = unique(vcat(flows[:DA_home],flows[:DA_work]))
    flow_dictionary =  Dict(vals[i] => i for i = 1:length(vals))
    flow_matrix = sparse([flow_dictionary[val] for  val in flows[:DA_home]],[flow_dictionary[val] for  val in flows[:DA_work]],flows[:FlowVolume])
    return flow_dictionary, flow_matrix
end

function elapsed(startt::Dates.DateTime)::Int
    Int(round((Dates.now()-startt).value/1000))
end

"""
Read data files and create a `SimData` object

**Arguments**
* `datapath` : path with csv files
* `filename` : name of csv file
* `colnames` : an array of columns which must be included in each file
* `road_levels` : a set with the road categories (see: OpenStreetMapX.ROAD_CLASSES for more informations)
* `google` : boolean variable; indicate whether google Distances API key will be included
"""
function get_sim_data(datapath::String;
                    road_levels::Set{Int} = Set(1:length(OpenStreetMapX.ROAD_CLASSES)),
					google::Bool = false)::OpenStreetMapXSim.SimData
    startt = Dates.now()
    files = collect(values(filenames))
    files = vcat(files...)
	files_in_dir = Set(readdir(datapath))
	found_error = false
	for filename in files
		check = filename != filenames[:googleapi_key] || (filename == filenames[:googleapi_key] && google)
		if !in(filename, files_in_dir) && check
			@error "The file $filename is missing in the directory $datapath"
			found_error = true
		end
	end
	found_error && error("Some file(s) not found in $datapath")
    @info "All config files found [$(elapsed(startt))s]"
    mapfile = filenames[:osm]
    map_data = OpenStreetMapX.get_map_data(datapath, mapfile; road_levels = road_levels)
    @info "Got map_data( data [$(elapsed(startt))s]"
    features, feature_classes, feature_to_intersections = get_features_data(datapath, mapfile, map_data)
    @info "Got feature_to_intersections data [$(elapsed(startt))s]"
    DAs_data = filenames[:DAs]
    DAs_to_intersection = DAs_to_nodes(datapath, DAs_data, colnames[:DAs], map_data)
	demographic_data = Dict(key => nothing for key in keys(DAs_to_intersection))
    @info "Got DAs_to_nodes data [$(elapsed(startt))s]"
	flow_stats = filenames[:flows]
	flow_dictionary, flow_matrix = get_flow_data(datapath,flow_stats, colnames[:flows])
	@info "Got flow_matrix data [$(elapsed(startt))s]"
	googleapi_key = nothing
	if google
		if !haskey(filenames, :googleapi_key)
			error("Google API key not defined! Declare one or run simulation based on OSM Data only!")
		end
		apikey = filenames[:googleapi_key]
		googleapi_key = open(joinpath(datapath,apikey)) do file
			read(file, String)
		end
	end
	@info "All data have been read with total of $(length(map_data.nodes)) map nodes  [$(elapsed(startt))s]"
    return OpenStreetMapXSim.SimData(map_data,
					features, feature_classes,
                    feature_to_intersections,
                    DAs_to_intersection,
					demographic_data,
					nothing,
					flow_dictionary,
					flow_matrix,
					googleapi_key)

end


text:
type Fourier2d <: SuperRes #1D Static model
    freqs :: Vector{Tuple{Int, Int}}
    filter :: Vector{Float64}
    x_max :: Float64
    z_max :: Float64
    approx_grid_x :: Vector{Float64}
    approx_grid_z :: Vector{Float64}
    bounds :: Tuple{Array{Float64}, Array{Float64}}
    function Fourier2d(x_max, z_max, filter, n_approx_x, n_approx_z)
        ind = find(filter)
        f_c_x = div(size(filter, 1) - 1, 2)
        f_c_z = div(size(filter, 2) - 1, 2)
        freqs_x = -f_c_x:f_c_x
        freqs_z = -f_c_z:f_c_z
        freqs2d = [(freq1, freq2) for freq1 in freqs_x, freq2 in freqs_z]
        freqs2d = freqs2d[ind]
        filter2d = filter[ind]
        grid_x = linspace(0, x_max, n_approx_x)
        grid_z = linspace(0, z_max, n_approx_z)
        return new(freqs2d,
                   filter2d,
                   x_max,
                   z_max,
                   grid_x,
                   grid_z,
                   ([0.0; 0.0], [x_max; z_max]))
    end
    function Fourier2d(x_max, filter, n_approx_x, n_approx_z)
        return Fourier2d(x_max, x_max, filter, n_approx_x, n_approx_z)
    end
end
function setBounds(model :: Fourier2d, bounds)
    model.bounds = bounds
    model.approx_grid_x = linspace(bounds[1][1], bounds[2][1], length(model.approx_grid_x))
    model.approx_grid_z = linspace(bounds[1][2], bounds[2][2], length(model.approx_grid_z))
end

function psi(model :: Fourier2d, theta :: Vector{Float64})
  # This function computes the direct problem for a single point theta
  psi_cpx = vec(Complex128[model.filter[k] *
    exp(-2im * pi * (model.freqs[k][1] / model.x_max * theta[1] + model.freqs[k][2] / model.z_max * theta[2] )) for k in eachindex(model.freqs)])
  return [real(psi_cpx); imag(psi_cpx)]
end

function dpsi(model :: Fourier2d, theta :: Vector{Float64})
  # This function computes the gradient of psi
  dx1 = vec(Complex128[-model.filter[k] *
    model.freqs[k][1] / model.x_max * 2im * pi *
    exp(-2im * pi * (model.freqs[k][1]  * theta[1]/model.x_max +  model.freqs[k][2] * theta[2]/model.z_max )) for k in eachindex(model.freqs)])
  dx2 = vec(Complex128[-model.filter[k] *
    model.freqs[k][2] / model.z_max * 2im * pi *
    exp(-2im * pi * (model.freqs[k][1] * theta[1]/model.x_max + model.freqs[k][2] * theta[2]/model.z_max )) for k in eachindex(model.freqs)])
  return [[real(dx1); imag(dx1)] [real(dx2); imag(dx2)]]
end

function getStartingPoint(model :: Fourier2d, v :: Vector{Float64})
  # This function gets an initial value for the location
  values = Float64[dot(v,psi(model, [x; y])) for x = model.approx_grid_x, y = model.approx_grid_z]
  ind = indmin(values);
  i,j = ind2sub(values, ind)
  return vec([model.approx_grid_x[i]; model.approx_grid_z[j]])
end

parameterBounds(model :: Fourier2d) =
    model.bounds

dim(model :: Fourier2d) = 2

type DynamicFourier2d <: DynamicSuperRes #1D Space-time model
    static :: Fourier2d
    times :: Vector{Float64}
    v_max :: Float64
    grid_v :: Vector{Float64}
    bounds :: Tuple{Array{Float64}, Array{Float64}}
    DynamicFourier2d(static, v_max, tau, K, num_v = 20) = new(static,
                                                              linspace(-tau*K, tau*K, 2*K+1),
                                                              v_max,
                                                              linspace(-v_max, v_max, num_v),
                                                              ([0; 0; -v_max;v_max], [static.x_max; static.z_max; v_max; v_max]))

end

function getStartingPoint(model :: DynamicFourier2d, v :: Vector{Float64})
    # This function gets an initial value for the location
    vals = Float64[dot(v,psi(model, [x; y; vel_x; vel_z])) for x = model.static.approx_grid_x, y = model.static.approx_grid_z, vel_x = model.grid_v, vel_z = model.grid_v]
    ind = indmin(vals);
    i,j,k,l = ind2sub(vals, ind)
    return [model.static.approx_grid_x[i]; model.static.approx_grid_z[j]; model.grid_v[k]; model.grid_v[l]]
end

parameterBounds(model :: DynamicFourier2d) =
    model.bounds


text:
using Test
using AllanDeviations
using Random

arrInt = zeros(Int, 5)
arr32 = zeros(Float32, 5)
arr64 = zeros(Float64, 5)

resInt = allandev(arrInt, 1.0, taus = AllTaus)
res32 = allandev(arr32, 1.0, taus = AllTaus)
res64 = allandev(arr64, 1.0, taus = AllTaus)

#Basic result type test
@test isa(resInt.deviation, Array{Int})
@test isa(res32.deviation, Array{Float32})
@test isa(res64.deviation, Array{Float64})
@test isa(res32.tau, Array{Float64})
@test isa(res64.tau, Array{Float64})
@test isa(res32.error, Array{Float32})
@test isa(res64.error, Array{Float64})
@test isa(res32.count, Array{Int})
@test isa(res64.count, Array{Int})

#Allandev from zero arrays is zero
@test sum(res32.deviation) == 0
@test sum(res64.deviation) == 0
#which also holds then for their errors
@test sum(res32.error) == 0
@test sum(res64.error) == 0

#and the result should be the same lenth for both types
@test length(res32.deviation) == length(res64.deviation)
@test length(res32.error) == length(res64.error)
@test length(res32.tau) == length(res64.tau)
@test length(res32.count) == length(res64.count)

#rate tests
arr64 = [1.0, 2.0, 1.0, 2.0, 1.5, 1.5, 2.0, 1.75]
res64 = allandev(arr64, 1.0)
res64r = allandev(arr64, 0.5)

@test sum(abs.(res64.deviation .- (2.0 .* res64r.deviation))) <= 2e-16 #half the rate means half the allan deviation
@test sum(abs.(res64.tau .- (0.5 .* res64r.tau))) <= 2e-16 #and double the tau
@test sum(abs.(res64.count .- res64r.count)) == 0 #but the count stays

#too few data points
@test_throws ErrorException allandev([1.0], 1.0)
@test_throws ErrorException allandev(zeros(Float64, 0), 1.0)
@test_throws ErrorException mallandev([1.0], 1.0)
@test_throws ErrorException hadamarddev([1.0], 1.0)
@test_throws ErrorException timedev([1.0], 1.0)
@test_throws ErrorException totaldev([1.0], 1.0)
@test_throws ErrorException mtie([1.0], 1.0)

#tau errors
@test_throws ErrorException allandev(arr64, 1.0, taus = 1.0)
@test_throws ErrorException allandev(arr64, 1.0, taus = 0.5)
@test_throws ErrorException allandev(arr64, 1.0, taus = -2.0)

#result and comparison tests
resallan = allandev(arr64, 1.0, taus = AllTaus)
resmallan = mallandev(arr64, 1.0, taus = AllTaus)
reshadamard = hadamarddev(arr64, 1.0, taus = AllTaus)
restime = timedev(arr64, 1.0, taus = AllTaus)
restotal = totaldev(arr64, 1.0, taus = AllTaus)
resmtie1 = mtie(arr64, 1.0, taus = AllTaus)
resmtie2 = mtie([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], 1.0, taus = AllTaus)
resmtie3 = mtie([8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0], 1.0, taus = AllTaus)

@test abs(resallan.deviation[1] - 0.97093168314425360) < 2e-16
@test abs(resallan.deviation[2] - 0.18221724671391565) < 2e-16
@test abs(resallan.deviation[3] - 0.20833333333333334) < 2e-16

@test abs(resmallan.deviation[1] - 0.9709316831442536) < 2e-16
@test abs(resmallan.deviation[2] - 0.0919975090242484) < 2e-16

@test abs(reshadamard.deviation[1] - 1.0616418102794056) < 2e-16
@test abs(reshadamard.deviation[2] - 0.1943203969393503) < 2e-16

@test abs(restime.deviation[1] - 0.56056766862807130) < 2e-16
@test abs(restime.deviation[2] - 0.10622957319984969) < 2e-16

@test abs(restotal.deviation[1] - 0.97093168314425360) < 2e-16
@test abs(restotal.deviation[2] - 0.35539004394233290) < 2e-16
@test abs(restotal.deviation[3] - 0.19982631347136330) < 2e-16
@test abs(restotal.deviation[4] - 0.27540216805004764) < 2e-16
@test abs(restotal.deviation[5] - 0.13919410907075053) < 2e-16
@test abs(restotal.deviation[6] - 0.14877975892797604) < 2e-16

@test length(resmtie1.deviation) == length(arr64) - 2
for i = 1:length(resmtie1.deviation)
	@test abs(resmtie1.deviation[i] - 1.0) < 2e-16
	@test abs(resmtie2.deviation[i] - i) < 2e-16
	@test abs(resmtie3.deviation[i] - i) < 2e-16
end

@test resallan.deviation[1] == resmallan.deviation[1] #first element of allan deviation and modified allan deviation is the same
@test abs(resallan.deviation[1] / sqrt(3) - restime.deviation[1]) < 2e-16 #first element of allan deviation and modified allan deviation is the same
@test reshadamard.count[1] < resallan.count[1] #the hadamarddeviation iterates over four terms
@test reshadamard.count[1] < resmallan.count[1] #the hadamarddeviation iterates over four terms
@test abs(resmallan.deviation[2] / sqrt(3) * 2 - restime.deviation[2]) < 2e-16 #test time deviation calculation

#overlapping tests
Random.seed!(0xA11E4DE71A7104_00)
arr64 = rand(512)
resallan = allandev(arr64, 1.0, taus = AllTaus).count
resallan_o = allandev(arr64, 1.0, overlapping = true, taus = AllTaus).count
resallan_c = allandev(arr64, 1.0, overlapping = false, taus = AllTaus).count
resmallan_o = mallandev(arr64, 1.0, overlapping = true, taus = AllTaus).count
resmallan_c = mallandev(arr64, 1.0, overlapping = false, taus = AllTaus).count
reshadamard_o = hadamarddev(arr64, 1.0, overlapping = true, taus = AllTaus).count
reshadamard_c = hadamarddev(arr64, 1.0, overlapping = false, taus = AllTaus).count
restotal_o = totaldev(arr64, 1.0, overlapping = true, taus = AllTaus).count
resmtie_o = mtie(arr64, 1.0, overlapping = true, taus = AllTaus).count

@test length(resallan) == length(resallan_o) #overlapping is standard
@test length(resallan_o) == 255
@test length(resallan_c) == 170
@test sum(resallan_o) == 65280
@test sum(resallan_c) == 2674

@test length(resmallan_o) == 170
@test length(resmallan_c) == 170
@test sum(resmallan_o) == 43605
@test sum(resmallan_c) == 2674

@test length(reshadamard_o) == 170
@test length(reshadamard_c) == 127
@test sum(reshadamard_o) == 43435
@test sum(reshadamard_c) == 2461

@test length(restotal_o) == 510
@test sum(restotal_o) == 260100

@test length(resmtie_o) == 510
@test sum(resmtie_o) == 130815

#frequency conversion tests
resallan = allandev(arr64, 1.0, taus = AllTaus)
resallan_p = allandev(arr64, 1.0, frequency = false, taus = AllTaus)
resallan_f = allandev(arr64, 1.0, frequency = true, taus = AllTaus)
resallan_f_r = allandev(arr64, 0.5, frequency = true, taus = AllTaus)
arr64_f = zeros(Float64, length(arr64) + 1)
arr64_f[2:length(arr64_f)] = cumsum(arr64)
arr64_f_r = zeros(Float64, length(arr64) + 1)
arr64_f_r[2:length(arr64_f)] = cumsum(arr64 * 2.0)
resallan_f_a = allandev(arr64_f, 1.0, frequency = false, taus = AllTaus)
resallan_f_a_r = allandev(arr64_f_r, 1.0, frequency = false, taus = AllTaus)

resmallan = mallandev(arr64, 1.0, frequency = true, taus = AllTaus)
reshadamard = hadamarddev(arr64, 1.0, frequency = true, taus = AllTaus)
restime = timedev(arr64, 1.0, frequency = true, taus = AllTaus)

@test sum(abs.(resallan.deviation .- resallan_p.deviation)) < 2e-16 #phase is standard
@test sum(abs.(resallan_p.deviation .- resallan_f.deviation)) > 2e-16 #frequency is not the same
@test sum(abs.(resallan_f.deviation .- resallan_f_a.deviation)) < 2e-13 #frequency conversion
@test sum(abs.(resallan_f_r.deviation .- (resallan_f_a_r.deviation .* 0.5))) < 2e-13 #frequency conversion with different rate

@test sum(abs(resallan_f.deviation[1] - resmallan.deviation[1])) < 2e-15 #also test for mallandev
@test abs(reshadamard.deviation[1] - 0.29893879132526296) < 2e-16 #hadamarddev
@test abs(restime.deviation[1] - 0.17291213240777910) < 2e-16 #timedev

#taus tests
resallan_a = allandev(arr64, 1.0, taus = AllTaus).count
resallan_qo = allandev(arr64, 1.0, taus = QuarterOctave).count
resallan_ho = allandev(arr64, 1.0, taus = HalfOctave).count
resallan_o = allandev(arr64, 1.0, taus = Octave).count
resallan_hd = allandev(arr64, 1.0, taus = HalfDecade).count
resallan_d = allandev(arr64, 1.0, taus = Decade).count

#all tau arguments
taus_o_a = 2.0 .^(0:floor(log2(length(arr64))))
resallan_o_c = allandev(arr64, 1.0, taus = 2.0).count
resallan_o_a = allandev(arr64, 1.0, taus = taus_o_a).count
#all tau arguments mallandev
resmallan = mallandev(arr64, 1.0, taus = Octave).count
resmallan_o_c = mallandev(arr64, 1.0, taus = 2.0).count
resmallan_o_a = mallandev(arr64, 1.0, taus = taus_o_a).count
#all tau arguments hadamarddev
reshadamard = hadamarddev(arr64, 1.0, taus = Octave).count
reshadamard_o_c = hadamarddev(arr64, 1.0, taus = 2.0).count
reshadamard_o_a = hadamarddev(arr64, 1.0, taus = taus_o_a).count

@test length(resallan_a) == 255
@test sum(resallan_a) == 65280

@test length(resallan_qo) == 21
@test sum(resallan_qo) == 8680

@test length(resallan_ho) == 13
@test sum(resallan_ho) == 5506

@test length(resallan_o) == 8
@test sum(resallan_o) == 3586
@test length(resallan_o_c) == 8 #float tau argument
@test sum(resallan_o_c) == 3586 #float tau argument
@test length(resallan_o_a) == 8 #array tau argument
@test sum(resallan_o_a) == 3586 #array tau argument

@test length(resallan_hd) == 4
@test sum(resallan_hd) == 1736

@test length(resallan_d) == 3
@test sum(resallan_d) == 1314

#mallandev tau tests
@test length(resmallan) == length(resmallan_o_c)
@test length(resmallan_o_a) == length(resmallan_o_c)
@test sum(resmallan) == sum(resmallan_o_c)
@test sum(resmallan_o_a) == sum(resmallan_o_c)

#hadamarddev tau tests
@test length(reshadamard) == length(reshadamard_o_c)
@test length(reshadamard_o_a) == length(reshadamard_o_c)
@test sum(reshadamard) == sum(reshadamard_o_c)
@test sum(reshadamard_o_a) == sum(reshadamard_o_c)


text:
# ------------ Trace typing system ------------ #

# Utilities.
keys(::Type{NamedTuple{K, V}}) where {K, V} = K
keys(::NamedTuple{K, V}) where {K, V} = K
value_types(::Type{NamedTuple{K, V}}) where {K, V} = V
NamedTuple{K, V}() where {K, V} = NamedTuple{K}(map(V.parameters) do p
                                                    p()
                                                end)

function Base.:(<<)(nt1::Type{NamedTuple{K1, V1}}, nt2::Type{NamedTuple{K2, V2}}) where {K1, K2, V1, V2}
    for k in K1
        k in K2 || return false
    end
    for (v, k) in zip(V1.parameters, V2.parameters)
        v << k || return false
    end
    true
end

# ------------ Support types ------------ #

abstract type SupportType end

abstract type BaseLebesgue <: SupportType end
pretty(::Type{BaseLebesgue}) = :lebesgue

abstract type BaseCounting <: SupportType end
pretty(::Type{BaseCounting}) = :counting

# Base types.
struct Reals{N} <: BaseLebesgue end
Base.:(<<)(::Type{Reals{N}}, ::Type{Reals{N}}) where N = true
Base.:(<<)(::Type{Reals}, a) = false

struct RealInterval <: BaseLebesgue
    a
    b
end

struct DiscreteInterval{A, B} <: BaseCounting end

struct PositiveReals{N} <: BaseLebesgue end
Base.:(<<)(::Type{PositiveReals{N}}, ::Type{PositiveReals{N}}) where N = true
Base.:(<<)(::Type{PositiveReals}, a) = false

struct Integers{N} <: BaseCounting end
Base.:(<<)(::Type{Integers{N}}, ::Type{Reals{N}}) where N = true
Base.:(<<)(::Type{Integers{N}}, ::Type{Integers{N}}) where N = true
Base.:(<<)(::Type{Integers}, a) = false

struct Discrete{N} <: BaseCounting end

# Derived types.
struct List{N} <: SupportType
    tt::N
end
List{N}() where N = List{N}(N())

struct Sum{T <: Tuple} <: SupportType
    tt::T
end
Sum{T}() where T <: Tuple = Sum{T}(map(T) do a
                                       a()
                                   end)

# Inferred record type is just a NamedTuple.
const TraceType = NamedTuple

# ------------ Tracer ------------ #

# Type inference.
function infer_support_types(ctx, fn, arg_types...)
    ir = lower_to_ir(fn, arg_types...)
    dh = detect_dynamic_addresses(ir)
    !isempty(dh) && begin
        display(dh)
        println()
        return missing
    end
    try
        tr = ir |> prepare_ir! |> ir -> infer!(ctx, ir)
        tr
    catch e
        @info "Failed to trace $(fn).\nCause: $e\n\nReverting to untyped IR."
        return missing
    end
end

# Extract trace type from inferred IR.
function trace_type(tr)
    ks = Set{Symbol}([])
    types = []
    for (v, st) in tr
        st.expr isa Expr || continue
        st.expr.head == :call || continue
        st.expr.args[1] == trace || continue
        st.expr.args[2] isa QuoteNode || continue
        push!(ks, st.expr.args[2].value)
        if !(st.type isa Type)
            push!(types, st.type)
        elseif st.type isa Type && (st.type <: SupportType || st.type <: Missing)
            push!(types, st.type())
        elseif _lift(st.type) <: NamedTuple
            new = NamedTuple{keys(st.type)}(map(value_types(st.type).parameters) do p
                                                p()
                                            end)
            push!(types, new)
        else
            push!(types, missing)
        end
    end
    NamedTuple{tuple(ks...)}(types)
end


text:
## Return the linear predictor vector
linpred(p::LinPred, f::Real=1.) = p.X.m * (f == 0. ? p.beta0 : fma(p.beta0, p.delbeta, f))

## Install beta0 + f*delbeta as beta0 and zero out delbeta
function installbeta!(p::LinPred, f::Real=1.)
    fma!(p.beta0, p.delbeta, f)
    fill!(p.delbeta, 0.)
    p.beta0
end

typealias BlasReal Union(Float32,Float64)

type DensePredQR{T<:BlasReal} <: DensePred
    X::ModelMatrix{T}             # model matrix
    beta0::Vector{T}              # base coefficient vector
    delbeta::Vector{T}            # coefficient increment
    qr::QR{T}
    function DensePredQR(X::ModelMatrix{T}, beta0::Vector{T})
        n,p = size(X); length(beta0) == p || error("dimension mismatch")
        new(X, beta0, zeros(T,p), qrfact(X.m))
    end
end
DensePredQR{T<:BlasReal}(X::ModelMatrix{T}) = DensePredQR{T}(X, zeros(T,size(X,2)))

delbeta!{T<:BlasReal}(p::DensePredQR{T}, r::Vector{T}) = (p.delbeta = p.qr\r; p)

type DensePredChol{T<:BlasReal} <: DensePred
    X::ModelMatrix{T}                   # model matrix
    beta0::Vector{T}                    # base vector for coefficients
    delbeta::Vector{T}                  # coefficient increment
    chol::Cholesky{T}
    function DensePredChol(X::ModelMatrix{T}, beta0::Vector{T})
        n,p = size(X); length(beta0) == p || error("dimension mismatch")
        new(X, beta0, zeros(T,p), cholfact(X.m'X.m))
    end
end
DensePredChol{T<:BlasReal}(X::ModelMatrix{T}) = DensePredChol{T}(X, zeros(T,size(X,2)))

solve!{T<:BlasReal}(C::Cholesky{T}, B::StridedVecOrMat{T}) = potrs!(C.uplo, C.UL, B)

function delbeta!{T<:BlasReal}(p::DensePredChol{T}, r::Vector{T})
    solve!(p.chol, gemv!('T', 1.0, p.X.m, r, 0.0, p.delbeta))
    p
end

function delbeta!{T<:BlasReal}(p::DensePredChol{T}, r::Vector{T}, wt::Vector{T}, scr::Matrix{T})
    vbroadcast!(Multiply(), scr, p.X.m, wt, 1)
    fac, info = potrf!('U', gemm!('T', 'N', 1.0, scr, p.X.m, 0.0, p.chol.UL))
    info == 0 || error("Singularity detected at column $info of weighted model matrix")
    solve!(p.chol, gemv!('T', 1.0, scr, r, 0.0, p.delbeta))
    p
end

coef(x::LinPred) = x.beta0
coef(x::LinPredModel) = coef(x.pp)

df_residual(x::LinPredModel) = df_residual(x.pp)
df_residual(x::DensePred) = size(x.X, 1) - length(x.beta0)

vcov(x::LinPredModel) = scale(x,true) * vcov(x.pp)
vcov(x::DensePredChol) = inv(x.chol)
vcov(x::DensePredQR) = symmetrize!(potri!('U', x.qr[:R]))

stderr(x::LinPredModel) = sqrt(diag(vcov(x)))

function show(io::IO, obj::LinPredModel)
    println(io, "\n$(obj.ff)\n\nCoefficients:\n")
    println(io, coeftable(obj))
end

## function show(io::IO, obj::GlmMod)
##     cc = coef(obj)
##     se = stderr(obj)
##     zz = cc ./ se
##     pp = 2.0 * ccdf(Normal(), abs(zz))
##     @printf("\n%s\n\nCoefficients:\n", obj.fr.formula)
##     @printf("         Term    Estimate  Std. Error     t value    Pr(>|t|)\n")
##     N = length(cc)
##     for i = 1:N
##         @printf(" %12s%12.5f%12.5f%12.3f%12.3f %-3s\n",
##                 obj.mm.model_colnames[i],
##                 cc[i],
##                 se[i],
##                 zz[i],
##                 pp[i],
##                 p_value_stars(pp[i]))
##     end
##     println("---\nSignif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n")
##     @printf("R-squared: %0.4f\n", 0.0) # TODO: obj.r_squared)
## end

## function p_value_stars(p_value::Float64)
##     if p_value < 0.001
##         return "***"
##     elseif p_value < 0.01
##         return "**"
##     elseif p_value < 0.05
##         return "*"
##     elseif p_value < 0.1
##         return "."
##     else
##         return " "
##     end
## end

ModelFrame(obj::LinPredModel) = obj.fr
ModelMatrix(obj::LinPredModel) = obj.pp.X
model_response(obj::LinPredModel) = obj.rr.y

fitted(m::LinPredModel) = m.rr.mu
predict(mm::LinPredModel) = fitted(mm)
formula(obj::LinPredModel) = ModelFrame(obj).formula
nobs(obj::LinPredModel) = length(model_response(obj))
residuals(obj::LinPredModel) = residuals(obj.rr)


text:
using Bombe
using Documenter

DocMeta.setdocmeta!(Bombe, :DocTestSetup, :(using Bombe); recursive=true)

makedocs(;
    modules=[Bombe],
    authors="Carlos Parada <paradac@carleton.edu>",
    repo="https://github.com/ParadaCarleton/Bombe.jl/blob/{commit}{path}#{line}",
    sitename="Bombe.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ParadaCarleton.github.io/Bombe.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ParadaCarleton/Bombe.jl",
)

max_stars_repo_name: ParadaCarleton/Bombe

max_issues_repo_name: ParadaCarleton/Bombe

max_forks_repo_name: ParadaCarleton/Bombe


text:
module GP

using Distances, LinearAlgebra, Optim

# structs
struct hyperparameters
    σ::AbstractVector
    ℓ::AbstractVector
    σ_n::Float64
    ℓx::Vector{Int64}
end

struct hyperparameters_multi
    σ::AbstractVector
    ℓ::AbstractVector
    σ_n::Float64
    ℓx::Vector{Int64}
end

struct GP_obj
    X::AbstractArray
    Lxx::AbstractArray
    α::AbstractArray
    Θ::hyperparameters
end

struct GP_obj_multi
    X::AbstractArray
    Lxx::AbstractArray
    α::AbstractArray
    Θ::hyperparameters_multi
end

struct SS
    A::AbstractArray
    B::AbstractArray
    C::AbstractArray
    D::AbstractArray
end

# squared exponential kernel
function k_SE(R, σ)
    K = (σ^2) * exp.(-(1/2) * R)
    return K
end

# calculate distance matrix
function r(x₁, x₂)
    r = pairwise(SqEuclidean(), x₁, x₂, dims=2)
    return r
end

# multiply kernels
function r_SE_mul(X₁, X₂, ℓ)
    R = zeros(size(X₁,1), size(X₂,1))
    for i in 1:size(X₁,2)
        R += r(reshape(X₁[:, i], 1, size(X₁,1)), reshape(X₂[:, i], 1, size(X₂,1))) ./ (ℓ[i]^2)
    end
    return R
end

# analytical derivative of posterior mean for dimension d
function δμδxsf(GP_obj,Xs,Σs,k,d)
    return -((GP_obj.Θ.ℓ[GP_obj.Θ.ℓx[d]].^(-2))*(GP_obj.X[:,d] .- Xs[:,d]') .* k(GP_obj.X,Xs,GP_obj.Θ,d))'*GP_obj.α
end

function εf(GP_obj,Xs,Σs,k)
    δμδxs = zeros(size(Xs))
    for d in 1:size(Xs,2)
        δμδxs[:,d] = δμδxsf(GP_obj,Xs,Σs,k,d)
    end
    return [(δμδxs[i,:]' * + diagm(Σs[i,:]) * δμδxs[i,:]) for i in 1:size(Xs,1)]
end

function μf(X, xs, k, Θ, α)
    return (k(X,xs,Θ)'*α)[]
end

# calculate posterior distribution
function posterior(GP_obj, Xs, Σs, k, calc_σ)
    Kxs = k(GP_obj.X,Xs,GP_obj.Θ)
    Kss = k(Xs,Xs,GP_obj.Θ)
    μ = (Kxs'*GP_obj.α)
    # uncertainty propagation
    if calc_σ==true
        σ = sqrt.(diag(Kss) - sum((GP_obj.Lxx \ Kxs).^2, dims=1)' + εf(GP_obj,Xs,Σs,k))
    else
        σ = zeros(size(Xs,1))
    end
    return μ, σ
end

# training kernel cholesky factorization
function calcLxx(X, y, k, Θ::hyperparameters)
    # kernel
    Kxx = k(X,X,Θ) + Θ.σ_n*I
    Cxx = cholesky(Kxx)
    Lxx = Cxx.L
    Lxx_y = (Lxx \ y)
    α = Lxx' \ Lxx_y
    return Lxx, α
end

# negative log marginal likelihood (for optimization of hyperparameters)
function nlml(X, y, k, Θₒₚₜ, Θ₀::hyperparameters)
    nσ, nℓ = size(Θ₀.σ,1), size(Θ₀.ℓ,1)
    Θ = hyperparameters(Θₒₚₜ[1:nσ], Θₒₚₜ[nσ+1:nσ+nℓ], Θ₀.σ_n, Θ₀.ℓx)
    # cholesky factorization of k(X,X')
    Lxx, α = calcLxx(X, y, k, Θ)
    # log marginal likelihood
    n = size(X,1)
    logml = -(1/2) * y' * α  - sum(log.(diag(Lxx))) - (n/2) * log(2*pi)
    print(".")
    return -logml[1]
end

# greedy learning of sparse approximation
function sparse_subset(X, y, ix, k, Θ, m)
    # dimensions of input
    d = size(X,2)
    # induced points
    Xm = X[1,:]'
    ym = y[1]
    Θm = Θ
    ixm = ix[1]
    # not-yet-chosen points
    Xn = X[2:end, :]
    yn = y[2:end]

    print("[")
    for j = 1:(m-1)
        Hmax=-1e5; xmax=zeros(size(X,2)); ymax=0; imax=0; xs=zeros(size(X,2)); ys=0
        print(".")
        for i in 1:size(yn,1)
            # new test point
            xs = Xn[i, :]'
            ys = yn[i]
            Xms = vcat(Xm, xs)
            yms = vcat(ym, ys)
            # calculate the entropy
            Lxx, α = calcLxx(Xms, yms, k, Θ)
            Hs = (1/2)*sum(log.(((2*pi*exp(1))^d)*(diag(Lxx))))
            #Kxx = k(Xms,Xms,Θ) + 0.001*I
            #Hs = (1/2)*sum(log.(((2*pi*exp(1))^d)*(det(Kxx))))
            # use point that maximizes differential entropy
            if Hs>Hmax
                #println("$Hs > $Hmax")
                Hmax = Hs
                xmax = xs
                ymax = ys
                imax = i
            end
        end
        if imax!=0
            Xm = vcat(Xm, xmax)
            ym = vcat(ym, ymax)
            ixm = vcat(ixm,ix[imax])
            Xn = Xn[[i for i in 1:size(Xn,1) if i!=imax],:]
            yn = yn[[i for i in 1:size(Xn,1) if i!=imax]]
        end
    end
    println("]")
    return Xm, ym, ixm
end

# optimize hyperparameters
function optimize_hyperparameters(X, y, k, Θ₀::hyperparameters)
    Θ₀_opt = vcat(Θ₀.σ, Θ₀.ℓ)
    nσ, nℓ = length(Θ₀.σ), length(Θ₀.ℓ)
    lower = vcat(zeros(nσ), zeros(nℓ))
    upper = [Inf for i in 1:(nσ+nℓ)]

    print("[")
    results = optimize(Θ -> nlml(X, y, k, Θ, Θ₀), lower, upper, Θ₀_opt, Fminbox(LBFGS())) #LBFGS()
    println("]")
    print(results)

    Θ₁ = results.minimizer
    Θp = hyperparameters(Θ₁[1:nσ], Θ₁[nσ+1:nσ+nℓ], Θ₀.σ_n, Θ₀.ℓx)

    return Θp
end

end

max_stars_repo_name: smnssn/GPDH

max_issues_repo_name: smnssn/GPDH

max_forks_repo_name: smnssn/GPDH


text:
ENV["PYTHON"] = ""
import Pkg
Pkg.add("PyPlot")
Pkg.build("PyPlot")


text:
using Clang.wrap_c
using Clang.cindex

llvmcfg = "llvm-config"

@static if Sys.isapple()
    using Homebrew
    !Homebrew.installed("llvm") && Homebrew.add("llvm")
    llvmcfg = joinpath(Homebrew.prefix(), "opt/llvm/bin/llvm-config")
end

const LLVM_VERSION = readchomp(`$llvmcfg --version`)
const LLVM_LIBDIR  = readchomp(`$llvmcfg --libdir`)
const LLVM_INCLUDE = joinpath(LLVM_LIBDIR, "clang", LLVM_VERSION, "include")

const BK_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include") |> normpath
const BK_HEADERS = [joinpath(root, file) for (root, dirs, files) in walkdir(BK_INCLUDE) for file in files]

function rewriter(ex::Expr)
    if Meta.isexpr(ex, :type)
        block = ex.args[3]
        isempty(block.args) && (typename = ex.args[2]; return :($typename = Void);)
    end
    Meta.isexpr(ex, :function) || return ex
    signature = ex.args[1]
    for i = 2:length(signature.args)
        func_arg = signature.args[i]
        if !(func_arg isa Symbol) && Meta.isexpr(func_arg, :(::))
            signature.args[i] = func_arg.args[1]
        end
    end
    return ex
end
rewriter(A::Array) = [rewriter(a) for a in A]
rewriter(arg) = arg

function wrap_header(top_hdr::String, cursor_header::String)
    startswith(dirname(cursor_header), BK_INCLUDE) && (top_hdr == cursor_header)
end

wc = wrap_c.init(;
                headers = BK_HEADERS,
                output_file = joinpath(@__DIR__, "api", "bk_api.jl"),
                common_file = joinpath(@__DIR__, "api", "bk_common.jl"),
                clang_includes = vcat(LLVM_INCLUDE, BK_INCLUDE),
                header_wrapped = wrap_header,
                header_library = x->"libbkmaxflow",
                clang_diagnostics = true,
                rewriter=rewriter)

# wrap_structs, immutable_structs
wc.options = wrap_c.InternalOptions(true, true)
run(wc)

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow