text:
# NLEQ-Options for Nonlinear Equations solvers

import Base: show

"""
super-type for all types storing options for NLEQ solvers.

NLEQ solvers often have serveral parameters for fine-tuning them.
In this NLEQInterface this parameters are called 'options' and
they are stored in key/value paris. For the key a
`AbstractString` is used. The value can be `Any`-thing.
The key is often called the option-name.

All types for this purpose have this abstract type as super-type.

Required fields are: `name`, `lastchanged`, `options`
"""
abstract AbstractOptionsNLEQ <: Any

"""
Stores options for NLEQ Solver(s) together with a name.
Additionally the time of the last change is saved.

Options can be set at construction time, e.g.

opt=OptionsNLEQ("test",
"loglevel" => NLEQ.LOG_ALL,
"logio"    => STDERR)

or later. For changing single options

oldValue = setOption!(opt,"myopt","new value")
oldValue = setOption!(opt,"myopt" => "new value")

and for changing many options at once:

oldValues = setOption!(opt, "myopt" => "new value",
    "oldopt" => 56)

see also: `setOption!`, `setOptions!`
"""
type OptionsNLEQ <: AbstractOptionsNLEQ
    name        :: AbstractString
    lastchanged :: DateTime
    options     :: Dict{AbstractString,Any}

    function OptionsNLEQ(name::AbstractString="")
        obj = new(name,now(),Dict{AbstractString,Any}())
        return obj
    end
end

"""
function OptionsNLEQ(name::AbstractString, copyOptionsFrom::AbstractOptionsNLEQ)
Copy constructor for OptionsNLEQ type object
"""
function OptionsNLEQ(name::AbstractString, copyOptionsFrom::AbstractOptionsNLEQ)
    opt = OptionsNLEQ(name)
    copyOptions!(opt,copyOptionsFrom)
    return opt
end

"""
function OptionsNLEQ(name::AbstractString, copyOptionsFrom::AbstractOptionsNLEQ)
Copy constructor for OptionsNLEQ type object
"""
function OptionsNLEQ(copyOptionsFrom::AbstractOptionsNLEQ)
    return OptionsNLEQ("",copyOptionsFrom)
end

"""
function OptionsNLEQ(name::AbstractString, copyOptionsFrom::AbstractOptionsNLEQ)
Constructor for OptionsNLEQ type object
"""
function OptionsNLEQ(name::AbstractString, pairs::Pair...)
    opt = OptionsNLEQ(name)
    setOptions!(opt,pairs...)
    return opt
end

"""
function OptionsNLEQ(name::AbstractString, copyOptionsFrom::AbstractOptionsNLEQ)
Constructor for OptionsNLEQ type object
"""
function OptionsNLEQ(pairs::Pair...)
    return OptionsNLEQ("",pairs...)
end

"""
function getOption(opt::AbstractOptionsNLEQ, name::AbstractString,
    default::Any=nothing)

get saved value of option with given `name` or set option to
`default` if option is unknown.
"""
function getOption(opt::AbstractOptionsNLEQ, name::AbstractString,
    default::Any=nothing)
    return haskey(opt.options,name)?opt.options[name]:default
end

"""
function getOption!(opt::AbstractOptionsNLEQ, name::AbstractString,
    default::Any=nothing)

function to get value and set default value in case no value was present
This avoids the recursive calls to setOption and getOption
"""
function getOption!(opt::AbstractOptionsNLEQ, name::AbstractString,
    default::Any=nothing)
    if haskey(opt.options,name)
        return opt.options[name]
    else
        setOption!(opt,name,default)
        return default
    end
end

"""
function setOption!(opt::AbstractOptionsNLEQ, name::AbstractString, value::Any)

set NLEQ-Option with given `name` and return old value
(or `nothing` if there was no old value).
"""
function setOption!(opt::AbstractOptionsNLEQ, name::AbstractString, value::Any)
    oldValue = getOption(opt,name)
    opt.options[name]=value
    opt.lastchanged=now()
    return oldValue
end

"""
function setOption!(opt::AbstractOptionsNLEQ, pair::Pair)

set NLEQ-Option with given (`name`,`value`) pair and return old value
(or `nothing` if there was no old value).
"""
function setOption!(opt::AbstractOptionsNLEQ, pair::Pair)
    return setOption!(opt,pair.first,pair.second)
end

"""
function setOptions!(opt::AbstractOptionsNLEQ, pairs::Pair...)

set many NLEQ-Options and return an array with the old option values.
"""
function setOptions!(opt::AbstractOptionsNLEQ, pairs::Pair...)
    oldValues=Any[]
    for (name,value) in pairs
        push!(oldValues, setOption!(opt,name,value))
    end
    return oldValues
end

"""
function initOption!(opt::AbstractOptionsNLEQ, name::AbstractString, value::Any)

initialize NLEQ-Option with given `name` only if it does not exist already and return the old value if changed
(or `nothing` if there was an old value).
"""
function initOption!(opt::AbstractOptionsNLEQ, name::AbstractString, value::Any)
    return haskey(opt.options,name)?nothing:setOption!(opt,name,value);
end

"""
function initOption!(opt::AbstractOptionsNLEQ, pair::Pair)

initialize NLEQ-Option with given (`name`,`value`) pair only if it does not exist already and return the old value if changed
(or `nothing` if there was an old value).
"""
function initOption!(opt::AbstractOptionsNLEQ, pair::Pair)
    return initOption!(opt,pair.first,pair.second);
end

"""
function initOptions!(opt::AbstractOptionsNLEQ, pairs::Pair...)

initialize many NLEQ-Options and return an array with the old option values if the options did not exist before
"""
function initOptions!(opt::AbstractOptionsNLEQ, pairs::Pair...)
    oldValues=Any[]
    for (name,value) in pairs
        push!(oldValues, initOption!(opt,name,value))
    end
    return oldValues
end

"""
function copyOptions!(dest::AbstractOptionsNLEQ, source::AbstractOptionsNLEQ)

copy all options from other NLEQ-Option object.
"""
function copyOptions!(dest::AbstractOptionsNLEQ, source::AbstractOptionsNLEQ)
    merge!(dest.options,source.options)
    dest.lastchanged=now()
    return dest
end

function show(io::IO, opt::AbstractOptionsNLEQ)
    print(io,typeof(opt)," ")
    isempty(opt.name) || print(io,"'",opt.name,"' ")
    len=length(opt.options)
    print(io,"with ",len," option",len!=1?"s":"",len>0?":":"."); println(io)
    if len>0
        maxLen=2+max( 0,map(length,keys(opt.options))... )
        for key in sort(collect(keys(opt.options)))
            print(io,lpad(key,maxLen),": ")
            show(io,opt.options[key]); println(io)
        end
    end
    print(io,"lasted changed "); show(io,opt.lastchanged);
end


text:
using AdvancedHMC

export AdHMC

"""
    AdHMC
The interface between ComradeSoss and the AdvancedHMC.jl package. To construct this you
should use the function.

```
AdHMC(cm::Soss.ConditionalModel)
```
"""
struct AdHMC{H,P,A,T} <: ComradeSoss.AbstractMCMC
    hamiltonian::H
    proposal::P
    adaptor::A
    transform::T
end

function AdHMC(cm::Soss.ConditionalModel;
              chunksize=5,
              start=nothing,
              minv = nothing,
              Metric=DiagEuclideanMetric,
              Integrator=Leapfrog,
              Treesample=MultinomialTS,
              Nuts=GeneralisedNoUTurn)

    t = xform(cm)
    dim = t.dimension
    function ℓπ(x)
        p, logjac = TV.transform_and_logjac(t, x)
        return logdensity(cm, p) + logjac
    end

    grado = GradObj(ℓπ, rand(dim); chunksize=chunksize)
    function gradℓπ(x)
        res = DiffResults.GradientResult(zeros(dim))
        ForwardDiff.gradient!(res, ℓπ, x, grado.cfg)
        return (DiffResults.value(res), DiffResults.gradient(res))
    end
    gradℓπ(zeros(dim))
    println("A derivative took ")
    @time gradℓπ(zeros(dim))
    println("Adjust your expectations accordingly")

    if isnothing(minv)
        metric = Metric(dim)
    else
        metric = Metric(minv)
    end
    hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
    if isnothing(start)
        x0 = zeros(dim)
    else
        x0 = TV.inverse(t, start)
    end
    initial_ϵ = find_good_stepsize(hamiltonian, x0)
    integrator = Integrator(initial_ϵ)
    proposal = NUTS{Treesample, Nuts}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))
    return AdHMC(hamiltonian, proposal, adaptor, t)
end

function AdHMCRev(cm::Soss.ConditionalModel;
    start=nothing,
    minv = nothing,
    Metric=DiagEuclideanMetric,
    Integrator=Leapfrog,
    Treesample=MultinomialTS,
    Nuts=GeneralisedNoUTurn)

    t = xform(cm)
    dim = t.dimension
    function ℓπ(x)
        p, logjac = TV.transform_and_logjac(t, x)
        return logdensity(cm, p) + logjac
    end

    grado = ReverseDiff.GradientTape(ℓπ, rand(dim))
    ctape = ReverseDiff.compile(grado)
    function gradℓπ(x)
        res = DiffResults.GradientResult(zeros(dim))
        ReverseDiff.gradient!(res, ctape, x)
        return DiffResults.value(res), DiffResults.gradient(res)
    end
    gradℓπ(zeros(dim))
    println("A derivative took ")
    @time gradℓπ(zeros(dim))
    println("Adjust your expectations accordingly")

    if isnothing(minv)
        metric = Metric(dim)
    else
        metric = Metric(minv)
    end
    hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
    if isnothing(start)
        x0 = zeros(dim)
    else
        x0 = TV.inverse(t, start)
    end
    initial_ϵ = find_good_stepsize(hamiltonian, x0)
    integrator = Integrator(initial_ϵ)
    proposal = NUTS{Treesample, Nuts}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))
    return AdHMC(hamiltonian, proposal, adaptor, t)
end

"""
    sample(s::AdHMC, nsamples, nadapt, start, kwargs...)
Sample the posterior using AdvancedHMC with a total of `nsamples` where the first `nadapt`
samples are adaptation. The starting location is given by `start` and this is assumed to
be in model parameter space not the continous one.
"""
function StatsBase.sample(smplr::AdHMC, nsamples, nadapt, θ0::NamedTuple; kwargs...)
    start = TV.inverse(smplr.transform, θ0)
    samples, stats = AdvancedHMC.sample(smplr.hamiltonian, smplr.proposal, start, nsamples, smplr.adaptor, nadapt; kwargs...)
    return TupleVector(TV.transform.(Ref(smplr.transform), samples)), stats
end

ext: jl


text:
function create_overview_stats(u::String; base_path::String=joinpath(Pinguu.base_dir[2:end], "stats/"))
    eval_date = string(today())
    invest_then, invest_now = Pinguu.get_invest_values(u)
    past_days, profits_daily = Pinguu.get_profit_for_past_days(u, 90)
    past_months, profits_monthly = Pinguu.get_profit_for_past_months(u, 24)
    overview = Dict(
        "bot_profit"  => Pinguu.get_total_profit(u),
        "bot_profit_today" => Pinguu.get_profit_for_single_day(u, eval_date),
        "active_bots" => Pinguu.get_num_active_bots(u),
        "bot_volume"  => Pinguu.get_total_bot_volume(u),
        "invest_then" => invest_then,
        "invest_now"  => invest_now,
        "past_days"   => past_days,
        "profits_daily" => profits_daily,
        "past_months"   => past_months,
        "profits_monthly" => profits_monthly
    )
    if u != "*"
        filename = joinpath(base_path, "users/" * u * "/profits.json")
    else
        filename = joinpath(base_path, "all_profits.json")
    end
    !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
    open(filename, "w") do f
        JSON.print(f, overview, 4)
    end
end
export create_overview_stats


text:
function solve24(nums)
    length(nums) != 4 && error("Input must be a 4-element Array")
    syms = [+,-,*,/]
    for x in syms, y in syms, z in syms
        for i = 1:24
            a,b,c,d = nthperm(nums,i)
            if round(x(y(a,b),z(c,d)),5) == 24
                return "($a$y$b)$x($c$z$d)"
            elseif round(x(a,y(b,z(c,d))),5) == 24
                return "$a$x($b$y($c$z$d))"
            elseif round(x(y(z(c,d),b),a),5) == 24
                return "(($c$z$d)$y$b)$x$a"
            elseif round(x(y(b,z(c,d)),a),5) == 24
                return "($b$y($c$z$d))$x$a"
            end
        end
    end
    return "0"
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using SequentialMonteCarlo
using RNGPool
using CoupledConditionalSMC
import SMCExamples.Particles.Float64Particle

@inline function simplelM(::Int64, particle::Float64Particle,
  newParticle::Float64Particle, ::Nothing)
  x::Float64 = particle.x
  y::Float64 = newParticle.x
  t::Float64 = x - y
  return -0.5 * t * t
end

function setupSimpleModel(n::Int64, boxSize::Float64)
  @inline function lG(::Int64, particle::Float64Particle, ::Nothing)
    return abs(particle.x) < boxSize ? 0.0 : -Inf
  end
  @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
    particle::Float64Particle, ::Nothing)
    if p == 1
      newParticle.x = randn(rng)
    else
      newParticle.x = particle.x + randn(rng)
    end
  end
  return SMCModel(M!, lG, n, Float64Particle, Nothing), simplelM
end

model, lM = setupSimpleModel(10000, 10.0)

import Statistics: mean, std

setRNGs(12345)

m = 1000

using DataFrames
box_df = DataFrame(n=Int[], N=Int[], Type=String[], mean=Union{Float64, Missing}[],
  std=Union{Float64, Missing}[])

for n in [500, 1000, 2000, 4000]
  for k in 4:7
    N = 2^k
    println(n, " ", N)
    maxit = 10*n
    vs = CoupledConditionalSMC.couplingTimes(model, N, n, m, true, true, maxit)
    push!(box_df, [n, N, "AT", mean(vs), std(vs)])
    vs = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m, :AS, true, true, maxit)
    push!(box_df, [n, N, "AS", mean(vs), std(vs)])
    vs = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m, :BS, true, true, maxit)
    push!(box_df, [n, N, "BS", mean(vs), std(vs)])
  end
end

using JLD2
@save "paper/box.jld2" box_df

println(box_df)


text:
# function to calculate the volume of a sphere
function sphere_vol(r)
    # julia allows Unicode names (in UTF-8 encoding)
    # so either "pi" or the symbol π can be used
    return 4/3*pi*r^3
end

max_stars_repo_name: AnjuJoon/sphinxcontrib-jupyter

max_issues_repo_name: AnjuJoon/sphinxcontrib-jupyter

max_forks_repo_name: AnjuJoon/sphinxcontrib-jupyter


text:
#####
##### `Thunk`
#####

macro thunk(body)
    return :(Thunk(() -> $(esc(body))))
end

struct Thunk{F}
    f::F
end

@inline (thunk::Thunk{F})() where {F} = (thunk.f)()

#####
##### `forward_chain`
#####

forward_chain(args...) = materialize(_forward_chain(args...))

@inline _forward_chain(ẋ::Nothing, ∂::Thunk) = false
@inline _forward_chain(ẋ, ∂::Thunk) = broadcasted(*, ẋ, ∂())
_forward_chain(ẋ, ∂::Thunk, args...) = broadcasted(+, _forward_chain(ẋ, ∂), _forward_chain(args...))

#####
##### `reverse_chain!`
#####

@inline reverse_chain!(x̄::Nothing, ∂::Thunk) = false

@inline function reverse_chain!(x̄, ∂::Thunk)
    thunk = ∂()
    x̄_value = adjoint_value(x̄)
    casted = should_increment(x̄) ? broadcasted(+, x̄_value, thunk) : thunk
    if should_materialize_into(x̄)
        return materialize!(x̄_value, casted)
    else
        return materialize(casted)
    end
end

adjoint_value(x̄) = x̄

should_increment(::Any) = true

should_materialize_into(::Any) = false

#####
##### miscellanous defaults
#####

# TODO: More defaults, obviously!

markup(::Any) = Ignore()
markup(::Real) = RealScalar()
markup(::Complex) = ComplexScalar()
markup(x::Tuple{Vararg{<:Real}}) = RealTensor(layout(x))
markup(x::Tuple{Vararg{<:Complex}}) = ComplexTensor(layout(x))
markup(x::AbstractArray{<:Real}) = RealTensor(layout(x))
markup(x::AbstractArray{<:Complex}) = ComplexTensor(layout(x))
markup(x::AbstractArray) = error("Cannot infer domain of array from eltype", x)

layout(x::Tuple) = Layout(length(x), (length(x),), CPUDevice(), false)
layout(x::Array) = Layout(length(x), size(x), CPUDevice(), true)

should_materialize_into(::Array) = true


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

# geometry setup in domain [0,1]^2
solid_left = 0.45
solid_right = 0.55
solid_top = 0.5
solid_bottom = 0.52

chip_left = 0.48
chip_right = 0.52
chip_top = 0.5
chip_bottom = 0.505

k_mold = 0.014531
k_chip = 2.60475
k_air = 0.64357
nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

m = 50
n = 50
h = 1/n
NT = 7    # number of iterations for Newton's method, 8 is good for m=400

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
chip_fem_idx = Array{Int64, 1}([])
chip_fvm_idx = Array{Int64, 1}([])
chip_fem_top_idx = Array{Int64, 1}([])

for i = 1:(m+1)
    for j = 1:(n+1)
        if (i-1)*h >= solid_left-1e-9 && (i-1)*h <= solid_right+1e-9 && (j-1)*h >= solid_top-1e-9 && (j-1)*h <= solid_bottom+1e-9
            # print(i, j)
            global solid_fem_idx = [solid_fem_idx; (j-1)*(m+1)+i]
            if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_bottom+1e-9
                global chip_fem_idx = [chip_fem_idx; (j-1)*(m+1)+i]
            end
            if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_top+1e-9
                global chip_fem_top_idx = [chip_fem_top_idx; (j-1)*(m+1)+i]
            end
        end
    end
end

for i = 1:m
    for j = 1:n
        if (i-1)*h + h/2 >= solid_left-1e-9 && (i-1)*h + h/2 <= solid_right+1e-9 &&
            (j-1)*h + h/2 >= solid_top-1e-9 && (j-1)*h + h/2 <= solid_bottom+1e-9
            global solid_fvm_idx = [solid_fvm_idx; (j-1)*m+i]
            if (i-1)*h + h/2 >= chip_left-1e-9 && (i-1)*h + h/2 <= chip_right+1e-9 && (j-1)*h + h/2 >= chip_top-1e-9 && (j-1)*h + h/2<= chip_bottom+1e-9
                global chip_fvm_idx = [chip_fvm_idx; (j-1)*m+i]
            end
        end
    end
end

k_fem = k_air * constant(ones((m+1)*(n+1)))
k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
k_fem = scatter_update(k_fem, chip_fem_idx, k_chip * ones(length(chip_fem_idx)))
kgauss = fem_to_gauss_points(k_fem, m, n, h)

heat_source_fem = zeros((m+1)*(n+1))
heat_source_fem[chip_fem_idx] .= power_source #/ h^2
heat_source_fem[chip_fem_top_idx] .= 82.46295

heat_source_gauss = fem_to_gauss_points(heat_source_fem, m, n, h)

B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
heat_source = constant(compute_fem_source_term1(heat_source_gauss, m, n, h))

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))
# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = fem_to_gauss_points(k, m, n, h)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, m, n, h)

bd = bcnode("all", m, n, h)

# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

bd = [bd; bd .+ (m+1)*(n+1);
     2*(m+1)*(n+1)+1; 2*(m+1)*(n+1)+m;
    #  (2*(m+1)*(n+1)+m*n )+1:(2*(m+1)*(n+1)+m*n )+m+1]
     bd .+ (2*(m+1)*(n+1)+m*n )]

# add solid region into boundary condition for u, v, p
bd = [bd; solid_fem_idx; solid_fem_idx .+ (m+1)*(n+1); solid_fvm_idx .+ 2(m+1)*(n+1)]

function compute_residual(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = fem_to_gauss_points(T, m, n, h)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, m, n, h)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, m, n, h) * T - heat_source
    R = [F;G;H0;T0]
    return R
end

function compute_jacobian(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]

    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, m, n, h))

    gradT = eval_grad_on_gauss_pts1(T, m, n, h)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = -buoyance_coef * constant(compute_fem_mass_matrix1(m, n, h))
    T_mat = [SparseTensor(spzeros((m+1)*(n+1), (m+1)*(n+1))); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(m*n, (m+1)*(n+1)))]

    N = 2*(m+1)*(n+1) + m*n
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros((m+1)*(n+1), m*n)) M]]
end

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(m*n+3*(m+1)*(n+1)))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n+(m+1)*(n+1)))

sess = Session(); init(sess)
output = run(sess, S)

matwrite("xz_chip_unstructured_data.mat",
    Dict(
        "V"=>output[end, :]
    ))

u0, v0, p0, t0 = zeros((m+1)*(n+1)), zeros((m+1)*(n+1)), zeros(m*n), zeros((m+1)*(n+1))

figure(figsize=(10,10))
subplot(221)
title("u velocity")
visualize_scalar_on_fem_points(output[NT+1, 1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
subplot(222)
title("v velocity")
visualize_scalar_on_fem_points(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
subplot(223)
visualize_scalar_on_fvm_points(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std, m, n, h);gca().invert_yaxis()
title("pressure")
subplot(224)
title("temperature")
visualize_scalar_on_fem_points(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty, m, n, h);gca().invert_yaxis()
tight_layout()

print("Solution range:",
    "\n [u velocity] \t min:", minimum(output[NT+1, 1:(m+1)*(n+1)] .* u_std), ",\t max:", maximum(output[NT+1, 1:(m+1)*(n+1)] .* u_std),
    "\n [v velocity] \t min:", minimum(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std), ",\t max:", maximum(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std),
    "\n [pressure]   \t min:", minimum(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std), ",\t max:", maximum(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std),
    "\n [temperature]\t min:", minimum(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty), ",\t\t\t max:", maximum(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty))

# separate plots
# figure();visualize_scalar_on_fem_points(output[NT+1, 1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis(); savefig("U.png")
# figure();visualize_scalar_on_fem_points(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis();savefig("V.png")
# figure();visualize_scalar_on_fvm_points(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std, m, n, h);savefig("P.png")
# figure();visualize_scalar_on_fem_points(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty, m, n, h);gca().invert_yaxis();savefig("T.png")

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

# prev_data = matread("steady_cavity_data.mat")["V"]

# final_u2=prev_data[NT+1, 1:(1+m)*(1+n)]
# final_v2=prev_data[NT+1, (1+m)*(1+n)+1:2*(m+1)*(n+1)]
# final_p2=prev_data[NT+1, 2*(m+1)*(n+1)+1:end]

# u12 = final_u2[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# u22 = final_u2[Int(n/2)+1:m+1:end]

# v12 = final_v2[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# v22 = final_v2[Int(n/2)+1:m+1:end]
####################################################################################
final_u=output[NT+1, 1:(1+m)*(1+n)] .* u_std
final_v=output[NT+1, (1+m)*(1+n)+1:2*(m+1)*(n+1)] .* u_std
final_p=output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std
final_t=output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty

u1 = final_u[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
u2 = final_u[Int(n/2)+1:m+1:end]

v1 = final_v[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
v2 = final_v[Int(n/2)+1:m+1:end]

t1 = final_t[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
t2 = final_t[Int(n/2)+1:m+1:end]
xx = 0:h:1
xx = xx .* 0.0305

figure();plot(xx, u1);#plot(xx, u12);
savefig("u_horizontal.png")

figure();plot(xx, u2);#plot(xx, u22);
savefig("u_vertical.png")

figure();plot(xx, v1);#plot(xx, v12);
savefig("v_horizontal.png")

figure();plot(xx, v2);#plot(xx, v22);
savefig("v_vertical.png")

figure();plot(xx, t1);#plot(xx, t12);
savefig("t_horizontal.png")

figure();plot(xx, t2);#plot(xx, t22);
savefig("t_vertical.png")

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

# p12 = final_p2[(Int(n/2)-1)*m+1: (Int(n/2)-1)*m+m]
# p22 = final_p2[Int(n/2)*m+1: Int(n/2)*m+m]
# p32 = 0.5 * (p12 .+ p22)

# p42 = final_p2[Int(n/2):m:end]
# p52 = final_p2[Int(n/2)+1:m:end]
# p62 = 0.5 * (p42 .+ p52)

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

p1 = final_p[(Int(n/2)-1)*m+1: (Int(n/2)-1)*m+m]
p2 = final_p[Int(n/2)*m+1: Int(n/2)*m+m]
p3 = 0.5 * (p1 .+ p2)

p4 = final_p[Int(n/2):m:end]
p5 = final_p[Int(n/2)+1:m:end]
p6 = 0.5 * (p4 .+ p5)

xx = h/2 :h:1
xx = xx .* 0.0305

figure();plot(xx, p3);#plot(xx, p32);
savefig("p_horizontal.png")

figure();plot(xx, p6);#plot(xx, p62);
savefig("p_vertical.png")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name     = "CGAL"
rversion = "5.3"
version  = VersionNumber(rversion)

# Collection of sources required to build CGAL
sources = [
    ArchiveSource("https://github.com/CGAL/cgal/releases/download/v$rversion/CGAL-$rversion.tar.xz",
                  "2c242e3f27655bc80b34e2fa5e32187a46003d2d9cd7dbec8fbcbc342cea2fb6"),
]

# Bash recipe for building across all platforms
script = raw"""
## pre-build setup
# exit on error
set -eu

cmake -B build \
  `# cmake specific` \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_FIND_ROOT_PATH=$prefix \
  -DCMAKE_INSTALL_PREFIX=$prefix \
  -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN \
  CGAL-*/

## and away we go..
cmake --build build --config Release --target install -- -j$nproc
install_license CGAL-*/LICENSE*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
# CGAL is, as of 5.0, a header-only library, removing support for lib
# compilation in 5.3
products = Product[]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Essential dependencies
    Dependency("boost_jll"; compat="=1.71.0"),
    Dependency("GMP_jll"; compat="6.1.2"),
    Dependency("MPFR_jll"; compat="4.0.2"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    imfilter_cpu(img, kern)

Perform filtering on `img` with kernel `kern` using the FFT algorithm.
"""
function imfilter_cpu(img::AbstractArray{T,N}, kern::AbstractArray{K,N}) where {T<:Real,K<:Real,N}
  imfilter(img, centered(kern), Inner(), Algorithm.FFT())
end

"""
    convdist(img, kern; [weights])

Compute distance between all subimages of `img` and `kern`.
Optionally, activate/deactivate pixels in the kern using `weights`.
"""
function convdist(img::AbstractArray, kern::AbstractArray; weights=nothing)
  # default to uniform weights
  weights == nothing && (weights = ones(kern))

  wkern = weights.*kern

  A² = imfilter_cpu(img.^2, weights)
  AB = imfilter_cpu(img, wkern)
  B² = sum(abs2, wkern)

  D = abs.(A² .- 2AB .+ B²)

  # always return a plain simple array
  parent(D)
end

"""
    selectpatch(imgs, tilesize, center)

Helper function to extract a patch of size
`tilesize` from `imgs` centered at Cartesian
index `center`.
"""
function selectpatch(imgs, tilesize, center)
  N = length(tilesize)

  # patch corners
  start  = ntuple(i -> center[i] - (tilesize[i]-1) ÷ 2, N)
  finish = ntuple(i -> center[i] +     tilesize[i] ÷ 2, N)

  [view(img, [start[i]:finish[i] for i=1:N]...) for img in imgs]
end


text:
#!/usr/bin/env julia --project

include("../shared_source/boilerplate.jl")
include("ThreeDimensionalCompressibleNavierStokesEquations.jl")

ClimateMachine.init()

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("refvals_buoyancy.jl")

    ########
    # Setup physical and numerical domains
    ########
    Ωˣ = IntervalDomain(-2π, 2π, periodic = true)
    Ωʸ = IntervalDomain(-2π, 2π, periodic = true)
    Ωᶻ = IntervalDomain(0, 4π, periodic = false)

    second_order = DiscretizedDomain(
        Ωˣ × Ωʸ × Ωᶻ;
        elements = 5,
        polynomial_order = 2,
        overintegration_order = 1,
    )

    fourth_order = DiscretizedDomain(
        Ωˣ × Ωʸ × Ωᶻ;
        elements = 3,
        polynomial_order = 4,
        overintegration_order = 1,
    )

    grids = Dict("second_order" => second_order, "fourth_order" => fourth_order)

    ########
    # Define timestepping parameters
    ########
    start_time = 0
    end_time = 0.1
    Δt = 0.001
    method = SSPRK22Heuns

    timestepper = TimeStepper(method = method, timestep = Δt)

    callbacks = (
        Info(),
        # VTKState(; iteration = 1, filepath = "output/buoyancy"),
        StateCheck(10),
    )

    ########
    # Define physical parameters and parameterizations
    ########
    Lˣ, Lʸ, Lᶻ = length(Ωˣ × Ωʸ × Ωᶻ)
    parameters = (
        ρₒ = 1, # reference density
        cₛ = sqrt(10), # sound speed
        α = 1e-4, # thermal expansion coefficient
        g = 10, # gravity
        θₒ = 10, # initial temperature value
        Lˣ = Lˣ,
        Lʸ = Lʸ,
        H = Lᶻ,
    )

    physics = FluidPhysics(;
        advection = NonLinearAdvectionTerm(),
        dissipation = ConstantViscosity{Float64}(μ = 0, ν = 0, κ = 0),
        coriolis = nothing,
        buoyancy = Buoyancy{Float64}(α = parameters.α, g = parameters.g),
    )

    ########
    # Define initial conditions
    ########

    u₀(p, x, y, z) = -0
    v₀(p, x, y, z) = -0
    w₀(p, x, y, z) = -0
    θ₀(p, x, y, z) = -p.θₒ * (1 - z / 4π)

    ρ₀(p, x, y, z) =
        p.ρₒ * (1 - (p.α * p.g / p.cₛ^2) / 2 * (-p.θₒ * (1 - z / 4π))^2)
    ρu₀(p, x...) = ρ₀(p, x...) * u₀(p, x...)
    ρv₀(p, x...) = ρ₀(p, x...) * v₀(p, x...)
    ρw₀(p, x...) = ρ₀(p, x...) * w₀(p, x...)
    ρθ₀(p, x...) = ρ₀(p, x...) * θ₀(p, x...)

    ρu⃗₀(p, x...) = @SVector [ρu₀(p, x...), ρv₀(p, x...), ρw₀(p, x...)]
    initial_conditions = (ρ = ρ₀, ρu = ρu⃗₀, ρθ = ρθ₀)

    orientations = Dict(
        "" => ClimateMachine.Orientations.NoOrientation(),
        "_flat" => ClimateMachine.Orientations.FlatOrientation(),
    )

    for (key1, grid) in grids
        for (key2, orientation) in orientations
            key = key1 * key2

            println("running ", key)

            local_physics = FluidPhysics(;
                orientation = orientation,
                advection = physics.advection,
                dissipation = physics.dissipation,
                coriolis = physics.coriolis,
                buoyancy = physics.buoyancy,
            )

            @testset "$(key)" begin
                model = SpatialModel(
                    balance_law = Fluid3D(),
                    physics = local_physics,
                    numerics = (flux = RoeNumericalFlux(),),
                    grid = grid,
                    boundary_conditions = NamedTuple(),
                    parameters = parameters,
                )

                simulation = Simulation(
                    model = model,
                    initial_conditions = initial_conditions,
                    timestepper = timestepper,
                    callbacks = callbacks,
                    time = (; start = start_time, finish = end_time),
                )

                ########
                # Run the model
                ########
                evolve!(
                    simulation,
                    model,
                    refDat = getproperty(refVals, Symbol(key)),
                )
            end
        end
    end
end


text:
Random.seed!(0)

# facts("qrdelcol") do
    m = 100
    A = randn(m,m)
    Q, R = qr(A)
    for i in 100:-1:1
        k = rand(1:i)
        global A = A[:,1:i .!= k]
        global R = qrdelcol(R, k)
        @assert norm( R'*R - A'*A ) < 1e-5
    end
# end


text:
getsize(width, height, depth) =
    2*(width*height+width*depth+height*depth),
    width*height*depth

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# A gaussian kernel function
@inline function Φ(x::Vector{T},
                   y::Vector{T},
                   r::T=1.0) where T<:AbstractFloat
    n  = 1.0 / sqrt(2π*r)
    s  = 1.0 / (2r^2)
    return n*exp(-s*sum((x.-y).^2))
end

# A kernel matrix
function ΦΦ(X::AbstractArray{T},
            r::T=1.0) where T<:AbstractFloat
    n = size(X,1)
    K = zeros(n,n)
    for i=1:n
        for j=1:i
            K[i, j] = Φ(X[i, :], X[j, :],r)
            K[j, i] = K[i, j]
        end
        K[i, i] = Φ(X[i, :], X[i, :],r)
    end
    K
end

# A kernel matrix for test data
function ΦΦ(X::AbstractArray{T},
            Z::AbstractArray{T},
            r::T=1.0) where T<:AbstractFloat
    (nx,mx)    = size(X)
    (nz,mz)    = size(Z)
    K          = zeros(T,nz, nx)
    for i=1:nz
        for j=1:nx
            K[i, j] = Φ(Z[i, :], X[j, :],r)
        end
    end
    K
end

@inline ∑(λ,y,n,K) = sum(λ .* y .* K)

function trainer(model::KernelPerceptron{T},
	              X::AbstractArray{T},
        			  Y::Vector{T}) where T<:AbstractFloat
   Y[Y .== 0]  .= -1 # fix in the future outside this function
   max_epochs  = model.max_epochs
   λ           = model.λ    # langrange multipliers
   K           = ΦΦ(X,model.width) # computing the kernel gram matrix
   n           = size(X,1)
   history     = []
   nerrors     = Inf
   epochs      = 0
   while  nerrors>0 && epochs < max_epochs
   # stops when error is equal to zero or grater than last_error or reached max iterations
       nerrors = 0
       # weight updates for all samples
       for i=1:n
          yp   = sign(∑(λ,Y,n,K[:,i]))
          if Y[i] != yp
             nerrors +=1
			 λ[i]    += 1    # missclassification counter for sample i
          end
		 end
       epochs+=1
       push!(history,nerrors)
   end
   if nerrors > 0
      warn("[Kernel Perceptron] Train not converged. Max epochs $(max_epochs) reached. Error history: $(history) \n Try to increase max_epochs or change kernel params.")
   end

   # storing only the tough samples ("support vectors")
   sv               = λ .> 0
   model.λ          = λ[sv]
   model.sv_x       = vec(X[sv,:])
   model.sv_y       = Y[sv]
   model.history    = history
   model.last_epoch = epochs

   #println("[Kernel perceptron] #$(length(model.λ)) support vectors out of $(n) samples.")

end

function predictor(model::KernelPerceptron{T},
	                X::AbstractArray{T}) where T<:AbstractFloat

   width        = model.width
   sv_x,sv_y,λ  = model.sv_x,model.sv_y,model.λ

   k   = size(sv_y,1)
   n   = size(X,1)
   y   = zeros(T,n)
   for i=1:n
      s = .0
      for j=1:k # can be vectorized in the future.
          s += λ[j] * sv_y[j] * Φ(X[i,:],sv_x[j,:],width) # this is simply a weighted voting into a kernel space
      end
      y[i] = s
   end
   y   = sign.(y)
   y[y .== -1]  .= 0 # fix in the future outside this function!!
   return y

end


text:
# 2520 is the smallest number that can be divided by each of the numbers from 1 to 10
# without any remainder.
# What is the smallest positive number that is evenly
# divisible by all of the numbers from 1 to 20?

# (1,2) (1,3) Means 1 * 2 * 1 * 3
# (2,2) means 2 * 2
include("./primes.jl")
using Primes

function pop(arr, value)
    inx = findfirst(x -> x == value, arr)
    if inx != 0
        deleteat!(arr, inx)
        return value
    end
    return 0
end

function merge(factor_set_1, factor_set_2)
    results = []
    for factor in factor_set_1
        pop(factor_set_2, factor)
        push!(results, factor)
    end

    for factor in factor_set_2
        push!(results, factor)
    end
    return results
end

@assert Primes.prime_factors(2) == [2]
@assert Primes.prime_factors(4) == [2,2]
@assert Primes.prime_factors(8) == [2,2,2]
@assert Primes.prime_factors(21) == [3,7]

@assert merge([2,2], [2,2,3]) == [2,2,3]
@assert merge([2,7], [2,2,3]) == [2,7,2,3]
@assert merge([7], [2,3,3]) == [7,2,3,3]

function min_common_factor(nb)
    factors = map(x -> Primes.prime_factors(x), 2:nb)
    min_factors = reduce(merge, [], factors)
    return prod(min_factors)
end

@assert min_common_factor(10) == 2520
println(min_common_factor(20))

max_stars_repo_name: ddugue/project-euler

max_issues_repo_name: ddugue/project-euler

max_forks_repo_name: ddugue/project-euler


text:
using DDD
using Test, SparseArrays

cd(@__DIR__)
@testset "BCC mobility" begin
    fileDislocationParameters = "../inputs/simParams/sampleDislocationParameters.json"
    fileMaterialParameters = "../inputs/simParams/sampleMaterialParameters.json"
    fileFEMParameters = "./testData/sampleFEMParameters.json"
    fileIntegrationParameters = "../inputs/simParams/sampleIntegrationParameters.json"
    fileSlipSystem = "../data/slipSystems/BCC.json"
    fileDislocationLoop = "../inputs/dln/samplePrismShear.json"
    dlnParams, matParams, femParams, intParams, slipSystems, missing = loadParameters(
        fileDislocationParameters,
        fileMaterialParameters,
        fileFEMParameters,
        fileIntegrationParameters,
        fileSlipSystem,
        fileDislocationLoop,
    )

    prismPentagon = DislocationLoop(;
        loopType = loopPrism(),
        numSides = 5,
        nodeSide = 1,
        numLoops = 1,
        segLen = 10 * ones(5),
        slipSystemIdx = 4,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )

    shearHexagon = DislocationLoop(;
        loopType = loopShear(),
        numSides = 6,
        nodeSide = 1,
        numLoops = 1,
        segLen = 10 * ones(6),
        slipSystemIdx = 4,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1; 1],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )

    network = DislocationNetwork([shearHexagon, prismPentagon], memBuffer = 1)

    regularCuboidMesh = buildMesh(matParams, femParams)
    f = spzeros(regularCuboidMesh.numNode * 3)
    fHat = spzeros(regularCuboidMesh.numNode * 3)
    u = spzeros(regularCuboidMesh.numNode * 3)
    uHat = spzeros(regularCuboidMesh.numNode * 3)
    forceDisplacement = ForceDisplacement(nothing, uHat, u, nothing, fHat, f)

    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)

    network2 = deepcopy(network)
    dlnMobility!(dlnParams, matParams, network2)

    force, vel = dlnMobility(dlnParams, matParams, network)
    testVel = [
        -0.036175006034684 -0.036175006034656 0.088938030619749
        0.060538472465831 0.060538472465831 0.121076825150312
        0.071350296347816 0.071350296347788 -0.018587478143454
        0.036175184446079 0.036175184446051 -0.088937984731268
        -0.060538377421600 -0.060538377421600 -0.121076879461304
        -0.071350353768921 -0.071350353768893 0.018587492151485
        -0.000039076066951 -0.000039076038383 0.000039076070642
        0.000014926772741 0.000014926761829 -0.000014926738200
        0.000048298670216 0.000048298634903 -0.000048298652559
        0.000014926727728 0.000014926716815 -0.000014926751357
        -0.000039076030376 -0.000039076001807 0.000039075998116
    ]

    testForce = [
        -0.228724254571166 -0.226639000305794 0.900128776494974
        0.686441227083849 0.686441227083850 1.372882454167666
        0.677022353942198 0.674937099676816 -0.003532577752966
        0.228724254571162 0.226639000305798 -0.900128776495018
        -0.686441227083830 -0.686441227083827 -1.372882454167693
        -0.677022353942161 -0.674937099676790 0.003532577752984
        -0.543143464917538 1.326707393542668 0.784736210352455
        0.991807802754387 0.277588328715755 1.268948329254914
        1.156114370762920 -1.155148398080229 -0.000482986341346
        -0.277289793905604 -0.991509267944234 -1.269246864065067
        -1.327488914694218 0.542361943765987 -0.783954689200902
    ]

    @test isapprox(vel', testVel, rtol = 1e-5)
    @test isapprox(force', testForce)
    @test isapprox(vel, network2.nodeVel)
    @test isapprox(force, network2.nodeForce)

    idx = rand(1:(network.numNode[1]))
    forceIdx, velIdx = dlnMobility(dlnParams, matParams, network, idx)
    @test isapprox(force[:, idx], forceIdx)
    @test isapprox(vel[:, idx], velIdx)
    network3 = deepcopy(network)
    dlnMobility!(dlnParams, matParams, network3, idx)
    @test isapprox(force[:, idx], network3.nodeForce[:, idx])
    @test isapprox(vel[:, idx], network3.nodeVel[:, idx])

    idx = rand(1:(network.numNode[1]), 5)
    forceIdx, velIdx = dlnMobility(dlnParams, matParams, network, idx)
    @test isapprox(force[:, idx], forceIdx)
    @test isapprox(vel[:, idx], velIdx)
    network4 = deepcopy(network)
    dlnMobility!(dlnParams, matParams, network4, idx)
    @test isapprox(force[:, idx], network4.nodeForce[:, idx])
    @test isapprox(vel[:, idx], network4.nodeVel[:, idx])
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
"""
    bin2int(bin: String) -> Integer

Return the decimal integer for the given `bin`.

This is used for **@piping.**

**Example**
```julia
bin2int("101")  # returns 5
```
"""
bin2int(str) = parse(Int, str; base=2)


text:
import CLBLAS
import OpenCL.cl

const clblas = CLBLAS

function randomMatrixSize(num::Int)
   # Format is (m,k,n) so matrix mul is (m,n)
   matSizes = Array(Tuple, num)
   for i in 1:num
      matSizes[i] = (rand(20:100) * 32, rand(20:100) * 32, rand(20:100) * 32)
   end
   return matSizes
end

function cpu_gemm(A, B, C, alpha, beta)
    startT = time()
    C = alpha*A*B + beta*C
    runtime = time() - startT
    runtime
end

function gpu_gemm(A, B, C, alpha, beta)
     const alphaF = cl.cl_float(alpha)
     const betaF  = cl.cl_float(beta)
     const u_zero = uint32(0)

     startT = time()
     future = clblas.clblasSgemm(u_zero, u_zero, alphaF, A, B, betaF, C)
     cl.api.clWaitForEvents(1, pointer(future.event))
     runtime = time() - startT
     runtime
end

function gpu_clblas_gemm(ctx, queue, AData, BData, CData, alpha, beta)
     A = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=AData)
     B = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=BData)
     C = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=CData)

     m, k = size(AData)
     k, n = size(BData)
     m = convert(Csize_t, m)
     k = convert(Csize_t, k)
     n = convert(Csize_t, n)

     lda = m
     ldb = k
     ldc = m

     event = cl.UserEvent(ctx)
     ptr_event = [pointer(event)]

     startT = time()
     clblas.clblasSgemm(1, 0, 0, m, n, k,
                             alpha, pointer(A), 0, lda,
                             pointer(B), 0, ldb, beta,
                             pointer(C), 0, ldc,
                             1, [pointer(queue)],
                             0, C_NULL, pointer(ptr_event))

     cl.api.clWaitForEvents(1, pointer(ptr_event))
     runtime = time() - startT
     runtime
end

function gpu_gemm_block(ctx, queue, kernel, AData, BData, CData, alpha, beta)
     A = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=AData)
     B = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=BData)
     C = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=CData)

     m, k = size(AData)
     k, n = size(BData)

     blocksize = 16
     localmem1 = cl.LocalMem(Float32, blocksize^2)
     localmem2 = cl.LocalMem(Float32, blocksize^2)

     evt = kernel[queue, (m, n), (blocksize, blocksize)](convert(Csize_t, m), convert(Csize_t, k), alpha, beta, C, A, B, localmem1, localmem2)
     runtime = evt[:profile_duration] / 1e9
     runtime
end

function many_gc(num)
   for id in 1:num
      gc()
   end
end

function run(device, context, queue, kernel, num)
   cpu_times     =  zeros(Float32, num)
   gpu_times     =  zeros(Float32, num)
   gpu_times_fut = zeros(Float32, num)
   gpu_times_blk = zeros(Float32, num)

   matrix_sizes = randomMatrixSize(num)
   alpha = 5.0
   beta = 2.0

   #gc_disable()
   for id in 1:num
      m, k, n = matrix_sizes[1]
      A = rand(Float32, m, k)
      B = rand(Float32, k, n)
      C = rand(Float32, m, n)

      gpu_times[id]     = gpu_clblas_gemm(context, queue, A, B, C, alpha, beta)
      gpu_times_fut[id] = gpu_gemm(A, B, C, alpha, beta)
      gpu_times_blk[id] = gpu_gemm_block(context, queue, kernel, A, B, C, alpha, beta)
      cpu_times[id]     = cpu_gemm(A, B, C, alpha, beta)
   end
   #gc_enable()

   (cpu_times, gpu_times, gpu_times_fut, gpu_times_blk)
end

function simulate(num, verbose = false)
    clblas.setup(true)
    if(verbose)
        for platform in cl.platforms()
            if platform[:name] == "Portable Computing Language"
                warn("Portable Computing Language platform not yet supported")
                continue
            end

            for device in cl.devices(platform)
                @printf("====================================================\n")
                @printf("Platform name:    %s\n",  platform[:name])
                @printf("Platform profile: %s\n",  platform[:profile])
                @printf("Platform vendor:  %s\n",  platform[:vendor])
                @printf("Platform version: %s\n",  platform[:version])
                @printf("----------------------------------------------------\n")
                @printf("Device name: %s\n", device[:name])
                @printf("Device type: %s\n", device[:device_type])
                @printf("Device mem: %i MB\n",           device[:global_mem_size] / 1024^2)
                @printf("Device max mem alloc: %i MB\n", device[:max_mem_alloc_size] / 1024^2)
                @printf("Device max clock freq: %i MHZ\n",  device[:max_clock_frequency])
                @printf("Device max compute units: %i\n",   device[:max_compute_units])
                @printf("Device max work group size: %i\n", device[:max_work_group_size])
                @printf("Device max work item size: %s\n",  device[:max_work_item_size])
           end
        end
    end

    device, context, queue = clblas.get_next_compute_context()
    kernel_source = open(readall, "GEMM-NVIDIA.cl")
    prg  = cl.Program(context, source=kernel_source) |> cl.build!
    mmul = cl.Kernel(prg, "mmul")

    run(device, context, queue, mmul, 2)

    cpu_times, gpu_times, gpu_times_fut, gpu_times_blk = run(device, context, queue, mmul, num)
    cpu_avg     = sum(cpu_times)/num
    gpu_avg     = sum(gpu_times)/num
    gpu_fut_avg = sum(gpu_times_fut)/num
    gpu_blk_avg = sum(gpu_times_blk)/num
    println("Average CPU    :", cpu_avg)
    println("Average GPU    :", gpu_avg)
    println("Average GPU FUT:", gpu_fut_avg)
    println("Average GPU BLK:", gpu_blk_avg)
    println("GPU BLAS xer CPU    :", cpu_avg / gpu_avg)
    println("GPU FUT  xer CPU    :", cpu_avg / gpu_fut_avg)
    println("GPU BLK  xer CPU    :", cpu_avg / gpu_blk_avg)

    clblas.teardown()
end

num = int(ARGS[1])
simulate(num)

#@pyimport matplotlib.pyplot as plt
#plt.plot([1:num], gpuMmulTimes2, color="green", linewidth=2.0, linestyle="--", label="GPU-Mul2")
#plt.plot([1:num], cpuTimes, color="red", linewidth=2.0, linestyle="--", label="CPU")
#plt.plot([1:num], gpuTimes, color="blue", linewidth=2.0, linestyle="--", label="GPU")
##plt.plot([1:num], gpuMmulTimes, color="yellow", linewidth=2.0, linestyle="--", label="GPU-Mul")

#plt.ylabel("time ns")
#plt.title("CPU vs GPU times")
#plt.show()


text:
using PyPlot
# == Set seed and generate a_t sequence == #
srand(123)
n = 100
a_seq = sin.(linspace(0, 5*pi, n)) + 2 + 0.1 * randn(n)

# == Model parameters == #
γ = 0.8
m = 1
d = γ * [1, -1]
h = 1.0

# == Initial conditions == #
y_m = [2.0]

testlq = LQFilter(d, h, y_m)
y_hist, L, U, y = optimal_y(testlq, a_seq)
y = y[end:-1:1]  # reverse y

# == Plot simulation results == #
fig, ax = subplots(figsize=(10, 6.1))
ax[:set_xlabel]("Time")

# == Some fancy plotting stuff -- simplify if you prefer == #
bbox = (0., 1.01, 1., .101)

time = 1:length(y)
ax[:set_xlim](0, maximum(time))
ax[:plot](time, a_seq / h, "k-o", ms=4, lw=2, alpha=0.6, label=L"$a_t$")
ax[:plot](time, y, "b-o", ms=4, lw=2, alpha=0.6, label=L"$y_t$")
ax[:grid]()
ax[:legend](ncol=2, bbox_to_anchor= bbox, loc = 3, mode = "expand", fontsize= 16)


text:
struct Figure1a
    θ₁::Vector{Float64}
    function Figure1a(x=[0.7,0.8,3.3,5,π,4,2.7,0.55,0.26])
        new(x)
    end
end
struct Figure1b
    θ₂::Vector{Float64}
    function Figure1b(x= vcat([0.7,0.8,3.3,5,π,4,2.7,0.55,0.26], [0.2,0.45,0.05,2π/3,0.2]))
        new(x)
    end
end
struct Figure1c
    θ₃::Vector{Float64}
    function Figure1c(x = vcat([0.7,1.2,3.3,5,-π/4,4,2.7,0.55,0.26], [0.1,0.5,0.05,π/3,0.1], [0.2,0.6,0.05,5π/3,0.1]))
        new(x)
    end
end
struct Figure1
    fig1a::Figure1a
    fig1b::Figure1b
    fig1c::Figure1c
    function Figure1()
        new(Figure1a(),Figure1b(),Figure1c())
    end
end

@recipe function f(f::Figure1a)
    ylims --> (0.2,1.8)
    legend --> true
    colorbar --> :right
    clims --> (0,1.25)
    colorbar_ticks --> (0:0.25:1.25)
    levels --> 50
    seriescolor := :thermal
    FrequencyDirectionSpectra.FDJONSWAP(f.θ₁),true
end

@recipe function f(f::Figure1b)
    ylims --> (0.2,1.8)
    legend --> true
    colorbar --> :right
    clims --> (0,1.25)
    colorbar_ticks --> (0:0.25:1.25)
    levels --> 50
    seriescolor := :thermal
    FrequencyDirectionSpectra.FDJONSWAP(f.θ₂[1:9]) + FrequencyDirectionSpectra.FDSwell(f.θ₂[10:end]),true
end

@recipe function f(f::Figure1c)
    ylims --> (0.2,1.8)
    legend --> true
    colorbar --> :right
    clims --> (0,1.25)
    colorbar_ticks --> (0:0.25:1.25)
    levels --> 50
    seriescolor := :thermal
    FrequencyDirectionSpectra.FDJONSWAP(f.θ₃[1:9]) + FrequencyDirectionSpectra.FDSwell(f.θ₃[10:14]) + FrequencyDirectionSpectra.FDSwell(f.θ₃[15:end]),true
end

@recipe function f(f::Figure1)
    ylims --> (0.2,1.8)
    layout := (1,3)
    clims --> (0,1.25)
    levels --> 50
    seriescolor := :thermal
    @series begin
        subplot := 1
        seriestype := :contour
        FrequencyDirectionSpectra.FDJONSWAP(f.fig1a.θ₁),true
    end
    @series begin
        subplot := 2
        seriestype := :contour
        FrequencyDirectionSpectra.FDJONSWAP(f.fig1b.θ₂[1:9]) + FrequencyDirectionSpectra.FDSwell(f.fig1b.θ₂[10:end]),true
    end
    @series begin
        subplot := 3
        seriestype := :contour
        FrequencyDirectionSpectra.FDJONSWAP(f.fig1c.θ₃[1:9]) + FrequencyDirectionSpectra.FDSwell(f.fig1c.θ₃[10:14]) + FrequencyDirectionSpectra.FDSwell(f.fig1c.θ₃[15:end]),true
    end
end

struct ColorFig1
    fig1::Figure1
end

@recipe function f(fc::ColorFig1)
    ylims --> (0.2,1.8)
    layout := 1
    clims --> (0,1.25)
    levels --> 50
    legend := true
    seriescolor := :thermal
    seriestype := :contour
    FrequencyDirectionSpectra.FDJONSWAP(fc.fig1.fig1a.θ₁),true
end


text:
"""
Configuration/Acquisition/Generation Command Set
"""
# sets and gets
const NUMBER_OF_RECORDS                 = U32(0x10000001)
const PRETRIGGER_AMOUNT                 = U32(0x10000002)
const RECORD_LENGTH                     = U32(0x10000003)
const TRIGGER_ENGINE                    = U32(0x10000004)
const TRIGGER_DELAY                     = U32(0x10000005)
const TRIGGER_TIMEOUT                   = U32(0x10000006)
const SAMPLE_RATE                       = U32(0x10000007)
const CONFIGURATION_MODE                = U32(0x10000008)
const DATA_WIDTH                        = U32(0x10000009)
const SAMPLE_SIZE                       = DATA_WIDTH
const AUTO_CALIBRATE                    = U32(0x1000000A)
const TRIGGER_XXXXX                     = U32(0x1000000B)
const CLOCK_SOURCE                      = U32(0x1000000C)
const CLOCK_SLOPE                       = U32(0x1000000D)
const IMPEDANCE                         = U32(0x1000000E)
const INPUT_RANGE                       = U32(0x1000000F)
const COUPLING                          = U32(0x10000010)
const MAX_TIMEOUTS_ALLOWED              = U32(0x10000011)
const OPERATING_MODE                    = U32(0x10000012)
const CLOCK_DECIMATION_EXTERNAL         = U32(0x10000013)
const LED_CONTROL                       = U32(0x10000014)
const ATTENUATOR_RELAY                  = U32(0x10000018)
const EXT_TRIGGER_COUPLING              = U32(0x1000001A)
const EXT_TRIGGER_ATTENUATOR_RELAY      = U32(0x1000001C)
const TRIGGER_ENGINE_SOURCE             = U32(0x1000001E)
const TRIGGER_ENGINE_SLOPE              = U32(0x10000020)
const SEND_DAC_VALUE                    = U32(0x10000021)
const SLEEP_DEVICE                      = U32(0x10000022)
const GET_DAC_VALUE                     = U32(0x10000023)
const GET_SERIAL_NUMBER                 = U32(0x10000024)
const GET_FIRST_CAL_DATE                = U32(0x10000025)
const GET_LATEST_CAL_DATE               = U32(0x10000026)
const GET_LATEST_TEST_DATE              = U32(0x10000027)
const SEND_RELAY_VALUE                  = U32(0x10000028)
const GET_LATEST_CAL_DATE_MONTH         = U32(0x1000002D)
const GET_LATEST_CAL_DATE_DAY           = U32(0x1000002E)
const GET_LATEST_CAL_DATE_YEAR          = U32(0x1000002F)
const GET_PCIE_LINK_SPEED               = U32(0x10000030)
const GET_PCIE_LINK_WIDTH               = U32(0x10000031)
const SETGET_ASYNC_BUFFSIZE_BYTES       = U32(0x10000039)
const SETGET_ASYNC_BUFFCOUNT            = U32(0x10000040)
const SET_DATA_FORMAT                   = U32(0x10000041)
const GET_DATA_FORMAT                   = U32(0x10000042)
const DATA_FORMAT_UNSIGNED              = U32(0x0)
const DATA_FORMAT_SIGNED                = U32(0x1)
const SET_SINGLE_CHANNEL_MODE           = U32(0x10000043)
const GET_SAMPLES_PER_TIMESTAMP_CLOCK   = U32(0x10000044)
const GET_RECORDS_CAPTURED              = U32(0x10000045)
const GET_MAX_PRETRIGGER_SAMPLES        = U32(0x10000046)
const SET_ADC_MODE                      = U32(0x10000047)
const ADC_MODE_DEFAULT                  = U32(0x0)
const ADC_MODE_DES                      = U32(0x1)
const ADC_MODE_DES_WIDEBAND             = U32(0x2)
const ADC_MODE_RESET_ENABLE             = U32(0x8001)      # not clear if should be U32
const ADC_MODE_RESET_DISABLE            = U32(0x8002)      # not clear if should be U32
const ECC_MODE                          = U32(0x10000048)
const ECC_DISABLE                       = U32(0x0)         # not clear if should be U32
const ECC_ENABLE                        = U32(0x1)         # not clear if should be U32
const GET_AUX_INPUT_LEVEL               = U32(0x10000049)
const AUX_INPUT_LOW                     = U32(0x0)         # not clear if should be U32
const AUX_INPUT_HIGH                    = U32(0x1)         # not clear if should be U32
const EXT_TRIGGER_IMPEDANCE             = U32(0x10000065)
const EXT_TRIG_50_OHMS                  = U32(0x0)         # not clear if should be UL
const EXT_TRIG_300_OHMS                 = U32(0x1)         # not clear if should be UL
const GET_CHANNELS_PER_BOARD            = U32(0x10000070)
const GET_CPF_DEVICE                    = U32(0x10000071)
const PACK_MODE                         = U32(0x10000072)
const PACK_DEFAULT                      = Clong(0x0)
const PACK_8_BITS_PER_SAMPLE            = Clong(0x1)
const PACK_12_BITS_PER_SAMPLE           = Clong(0x2)
const GET_FPGA_TEMPERATURE              = U32(0x10000080)
const API_FLAGS                         = U32(0x10000090)
const API_ENABLE_TRACE                  = U32(0x1)         # not clear if should be UL

# gets board specific parameters
const MEMORY_SIZE                       = U32(0x1000002A)
const MEMORY_SIZE_MSAMPLES              = U32(0x1000004A)
const BOARD_TYPE                        = U32(0x1000002B)
const ASOPC_TYPE                        = U32(0x1000002C)
const GET_BOARD_OPTIONS_LOW             = U32(0x10000037)
const GET_BOARD_OPTIONS_HIGH            = U32(0x10000038)
const OPTION_STREAMING_DMA              = <<(U32(1),0)
const OPTION_EXTERNAL_CLOCK             = <<(U32(1),1)
const OPTION_DUAL_PORT_MEMORY           = <<(U32(1),2)
const OPTION_180MHZ_OSCILLATOR          = <<(U32(1),3)
const OPTION_LVTTL_EXT_CLOCK            = <<(U32(1),4)
const OPTION_SW_SPI                     = <<(U32(1),5)
const OPTION_ALT_INPUT_RANGES           = <<(U32(1),6)
const OPTION_VARIABLE_RATE_10MHZ_PLL    = <<(U32(1),7)
const OPTION_2GHZ_ADC                   = <<(U32(1),8)
const OPTION_DUAL_EDGE_SAMPLING         = <<(U32(1),9)
const OPTION_OEM_FPGA                   = <<(U32(1),47)

# sets and gets
# The transfer offset is defined as the place to start
# the transfer relative to trigger. The value is signed.
# -------TO>>>T>>>>>>>>>TE------------
const TRANSFER_OFFET                    = U32(0x10000030)  # typo in C API?
const TRANSFER_LENGTH                   = U32(0x10000031)

# Transfer related constants
const TRANSFER_RECORD_OFFSET            = U32(0x10000032)
const TRANSFER_NUM_OF_RECORDS           = U32(0x10000033)
const TRANSFER_MAPPING_RATIO            = U32(0x10000034)

# only gets
const TRIGGER_ADDRESS_AND_TIMESTAMP     = U32(0x10000035)

# MASTER/SLAVE CONTROL sets/gets
const MASTER_SLAVE_INDEPENDENT          = U32(0x10000036)

# boolean gets
const TRIGGERED                         = U32(0x10000040)
const BUSY                              = U32(0x10000041)
const WHO_TRIGGERED                     = U32(0x10000042)
const GET_ASYNC_BUFFERS_PENDING         = U32(0x10000050)
const GET_ASYNC_BUFFERS_PENDING_FULL    = U32(0x10000051)
const GET_ASYNC_BUFFERS_PENDING_EMPTY   = U32(0x10000052)
const ACF_SAMPLES_PER_RECORD            = U32(0x10000060)
const ACF_RECORDS_TO_AVERAGE            = U32(0x10000061)
const ACF_MODE                          = U32(0x10000062)

"""
Sample rates that the internal clock of a board can generate.

Note: Available sample rates for a given board form a subset of
this class' members. Please see your board's specification as well
as the ATS-SDK manual for more information.
"""
# Sample rate values
const SAMPLE_RATE_1KSPS     = U32(0x1)
const SAMPLE_RATE_2KSPS     = U32(0x2)
const SAMPLE_RATE_5KSPS     = U32(0x5)
const SAMPLE_RATE_10KSPS    = U32(0x8)
const SAMPLE_RATE_20KSPS    = U32(0xA)
const SAMPLE_RATE_50KSPS    = U32(0xC)
const SAMPLE_RATE_100KSPS   = U32(0xE)
const SAMPLE_RATE_200KSPS   = U32(0x10)
const SAMPLE_RATE_500KSPS   = U32(0x12)
const SAMPLE_RATE_1MSPS     = U32(0x14)
const SAMPLE_RATE_2MSPS     = U32(0x18)
const SAMPLE_RATE_5MSPS     = U32(0x1A)
const SAMPLE_RATE_10MSPS    = U32(0x1C)
const SAMPLE_RATE_20MSPS    = U32(0x1E)
const SAMPLE_RATE_25MSPS    = U32(0x21)
const SAMPLE_RATE_50MSPS    = U32(0x22)
const SAMPLE_RATE_100MSPS   = U32(0x24)
const SAMPLE_RATE_125MSPS   = U32(0x25)
const SAMPLE_RATE_160MSPS   = U32(0x26)
const SAMPLE_RATE_180MSPS   = U32(0x27)
const SAMPLE_RATE_200MSPS   = U32(0x28)
const SAMPLE_RATE_250MSPS   = U32(0x2B)
const SAMPLE_RATE_400MSPS   = U32(0x2D)
const SAMPLE_RATE_500MSPS   = U32(0x30)
const SAMPLE_RATE_800MSPS   = U32(0x32)
const SAMPLE_RATE_1000MSPS  = U32(0x35)
const SAMPLE_RATE_1GSPS     = SAMPLE_RATE_1000MSPS
const SAMPLE_RATE_1200MSPS  = U32(0x37)
const SAMPLE_RATE_1500MSPS  = U32(0x3A)
const SAMPLE_RATE_1600MSPS  = U32(0x3B)
const SAMPLE_RATE_1800MSPS  = U32(0x3D)
const SAMPLE_RATE_2000MSPS  = U32(0x3F)
const SAMPLE_RATE_2GSPS     = SAMPLE_RATE_2000MSPS
const SAMPLE_RATE_2400MSPS  = U32(0x6A)
const SAMPLE_RATE_3000MSPS  = U32(0x75)
const SAMPLE_RATE_3GSPS     = SAMPLE_RATE_3000MSPS
const SAMPLE_RATE_3600MSPS  = U32(0x7B)
const SAMPLE_RATE_4000MSPS  = U32(0x80)
const SAMPLE_RATE_4GSPS     = SAMPLE_RATE_4000MSPS

# User-defined sample rate; used with External Clock
const SAMPLE_RATE_USER_DEF  = U32(0x40)

"""
ATS665 Specific Setting for using the PLL.
The base value can be used to create a PLL frequency
in a simple manner.

Ex.
      105 MHz = PLL_10MHZ_REF_100MSPS_BASE + 5000000
      120 MHz = PLL_10MHZ_REF_100MSPS_BASE + 20000000
"""
const PLL_10MHZ_REF_100MSPS_BASE = U32(0x05F5E100)

# ATS665 Specific Decimation constants
const DECIMATE_BY_8              = U32(0x00000008)
const DECIMATE_BY_64             = U32(0x00000040)

#impedancevalues #suchresistance
const IMPEDANCE_1M_OHM           = U32(0x00000001)
const IMPEDANCE_50_OHM           = U32(0x00000002)
const IMPEDANCE_75_OHM           = U32(0x00000004)
const IMPEDANCE_300_OHM          = U32(0x00000008)

"""
Types of clocks that a board can use for acquiring data.
Note: Available sources for a given board form a subset of this
class' members. Please see your board's specification as well as
the ATS-SDK manual for more information.
"""
const INTERNAL_CLOCK            = U32(0x1)
const EXTERNAL_CLOCK            = U32(0x2)
const FAST_EXTERNAL_CLOCK       = U32(0x2)
const MEDIUM_EXTERNAL_CLOCK     = U32(0x3)
const SLOW_EXTERNAL_CLOCK       = U32(0x4)
const EXTERNAL_CLOCK_AC         = U32(0x5)
const EXTERNAL_CLOCK_DC         = U32(0x6)
const EXTERNAL_CLOCK_10MHz_REF  = U32(0x7)
const INTERNAL_CLOCK_10MHz_REF  = U32(0x8)    # 0x8 to 0x12 not ULd in C API
const EXTERNAL_CLOCK_10MHz_PXI  = U32(0xA)
const INTERNAL_CLOCK_DIV_4      = U32(0xF)
const INTERNAL_CLOCK_DIV_5      = U32(0x10)
const MASTER_CLOCK              = U32(0x11)
const INTERNAL_CLOCK_SET_VCO    = U32(0x12)

"""
Direction of the edge from the external clock signal that the board
syncrhonises with.
"""
const CLOCK_EDGE_RISING     = U32(0x0)
const CLOCK_EDGE_FALLING    = U32(0x1)

"""
Board input ranges (amplitudes) identifiers. PM stands for
plus/minus.

Note: Available input ranges for a given board *and* a given
configuration form a subset of this class' members. Please see
your board's specification as well as the ATS-SDK manual for more
information.
"""
const INPUT_RANGE_PM_20_MV	= U32(0x00000001)     # was not in Python API
const INPUT_RANGE_PM_40_MV  = U32(0x2)
const INPUT_RANGE_PM_50_MV  = U32(0x3)
const INPUT_RANGE_PM_80_MV  = U32(0x4)
const INPUT_RANGE_PM_100_MV = U32(0x5)
const INPUT_RANGE_PM_200_MV = U32(0x6)
const INPUT_RANGE_PM_400_MV = U32(0x7)
const INPUT_RANGE_PM_500_MV = U32(0x8)
const INPUT_RANGE_PM_800_MV = U32(0x9)
const INPUT_RANGE_PM_1_V    = U32(0xA)
const INPUT_RANGE_PM_2_V    = U32(0xB)
const INPUT_RANGE_PM_4_V    = U32(0xC)
const INPUT_RANGE_PM_5_V    = U32(0xD)
const INPUT_RANGE_PM_8_V    = U32(0xE)
const INPUT_RANGE_PM_10_V   = U32(0xF)
const INPUT_RANGE_PM_20_V   = U32(0x10)
const INPUT_RANGE_PM_40_V   = U32(0x11)
const INPUT_RANGE_PM_16_V   = U32(0x12)
const INPUT_RANGE_HIFI      = U32(0x20)
const INPUT_RANGE_PM_1_V_25 = U32(0x21)
const INPUT_RANGE_PM_2_V_5  = U32(0x25)
const INPUT_RANGE_PM_125_MV = U32(0x28)
const INPUT_RANGE_PM_250_MV = U32(0x30)

"""Coupling types identifiers for all boards input"""
const AC_COUPLING           = U32(1)
const DC_COUPLING           = U32(2)

"""Trigger engine identifiers."""
const TRIG_ENGINE_J         = U32(0)
const TRIG_ENGINE_K         = U32(1)

"""Trigger engine operation identifiers."""
const TRIG_ENGINE_OP_J            = U32(0)
const TRIG_ENGINE_OP_K            = U32(1)
const TRIG_ENGINE_OP_J_OR_K       = U32(2)
const TRIG_ENGINE_OP_J_AND_K      = U32(3)
const TRIG_ENGINE_OP_J_XOR_K      = U32(4)
const TRIG_ENGINE_OP_J_AND_NOT_K  = U32(5)
const TRIG_ENGINE_OP_NOT_J_AND_K  = U32(6)

"""Types of input that the board can trig on."""
const TRIG_CHAN_A   = U32(0x00000000)
const TRIG_CHAN_B   = U32(0x00000001)
const TRIG_CHAN_C   = U32(0x00000004)
const TRIG_CHAN_D   = U32(0x00000005)
const TRIG_CHAN_E   = U32(0x00000006)
const TRIG_CHAN_F   = U32(0x00000007)
const TRIG_CHAN_G   = U32(0x00000008)
const TRIG_CHAN_H   = U32(0x00000009)
const TRIG_CHAN_I   = U32(0x0000000A)
const TRIG_CHAN_J   = U32(0x0000000B)
const TRIG_CHAN_K   = U32(0x0000000C)
const TRIG_CHAN_L   = U32(0x0000000D)
const TRIG_CHAN_M   = U32(0x0000000E)
const TRIG_CHAN_N   = U32(0x0000000F)
const TRIG_CHAN_O   = U32(0x00000010)
const TRIG_CHAN_P   = U32(0x00000011)
const TRIG_PXI_STAR = U32(0x00000100)
const TRIG_EXTERNAL = U32(0x00000002)
const TRIG_DISABLE  = U32(0x00000003)

"""Edge of the external trigger signal that the board syncrhonises with."""
const TRIGGER_SLOPE_POSITIVE = U32(1)
const TRIGGER_SLOPE_NEGATIVE = U32(2)

"""
Board input channel identifiers

Note: The channels available for a given board form a subset of this
class' members. Please see your board's specification as well as
the ATS-SDK manual for more information.
"""
const CHANNEL_ALL   = U8(0x0000)
const CHANNEL_A     = U8(0x0001)
const CHANNEL_B     = U8(0x0002)
const CHANNEL_C     = U8(0x0004)
const CHANNEL_D     = U8(0x0008)
const CHANNEL_E     = U8(0x0010)
const CHANNEL_F     = U8(0x0020)
const CHANNEL_G     = U8(0x0040)
const CHANNEL_H     = U8(0x0080)

# Master/slave configuration
const BOARD_IS_INDEPENDENT  = U32(0x00000000)
const BOARD_IS_MASTER       = U32(0x00000001)
const BOARD_IS_SLAVE        = U32(0x00000002)
const BOARD_IS_LAST_SLAVE   = U32(0x00000003)

# LED control
const LED_OFF               = U32(0x00000000)
const LED_ON                = U32(0x00000001)

# Attenuator relay
const AR_X1                 = U32(0x00000000)
const AR_DIV40              = U32(0x00000001)

# External Trigger Attenuator Relay
const ETR_DIV5              = U32(0x00000000)
const ETR_X1                = U32(0x00000001)
const ETR_5V                = U32(0x00000000)
const ETR_1V                = U32(0x00000001)
const ETR_TTL               = U32(0x00000002)
const ETR_2V5               = U32(0x00000003)

# Device sleep state
const POWER_OFF             = U32(0x00000000)
const POWER_ON              = U32(0x00000001)

# Software events control
const SW_EVENTS_OFF         = U32(0x00000000)
const SW_EVENTS_ON          = U32(0x00000001)

# Timestamp value reset control
const TIMESTAMP_RESET_FIRSTTIME_ONLY  = U32(0x00000000)
const TIMESTAMP_RESET_ALWAYS          = U32(0x00000001)

"""
DAC Names used by API AlazarDACSettingAdjust
"""
# DAC Names Specific to the ATS460
const ATS460_DAC_A_GAIN           = U32(0x00000001)
const ATS460_DAC_A_OFFSET         = U32(0x00000002)
const ATS460_DAC_A_POSITION       = U32(0x00000003)
const ATS460_DAC_B_GAIN           = U32(0x00000009)
const ATS460_DAC_B_OFFSET         = U32(0x0000000A)
const ATS460_DAC_B_POSITION       = U32(0x0000000B)
const ATS460_DAC_EXTERNAL_CLK_REF = U32(0x00000007)

# DAC Names Specific to the ATS660
const ATS660_DAC_A_GAIN           = U32(0x00000001)
const ATS660_DAC_A_OFFSET         = U32(0x00000002)
const ATS660_DAC_A_POSITION       = U32(0x00000003)
const ATS660_DAC_B_GAIN           = U32(0x00000009)
const ATS660_DAC_B_OFFSET         = U32(0x0000000A)
const ATS660_DAC_B_POSITION       = U32(0x0000000B)
const ATS660_DAC_EXTERNAL_CLK_REF = U32(0x00000007)

# DAC Names Specific to the ATS665
const ATS665_DAC_A_GAIN           = U32(0x00000001)
const ATS665_DAC_A_OFFSET         = U32(0x00000002)
const ATS665_DAC_A_POSITION       = U32(0x00000003)
const ATS665_DAC_B_GAIN           = U32(0x00000009)
const ATS665_DAC_B_OFFSET         = U32(0x0000000A)
const ATS665_DAC_B_POSITION       = U32(0x0000000B)
const ATS665_DAC_EXTERNAL_CLK_REF = U32(0x00000007)

# Error return values
const SETDAC_INVALID_SETGET	   =  660
const SETDAC_INVALID_CHANNEL   =  661
const SETDAC_INVALID_DACNAME   =  662
const SETDAC_INVALID_COUPLING  =  663
const SETDAC_INVALID_RANGE	   =  664
const SETDAC_INVALID_IMPEDANCE =  665
const SETDAC_BAD_GET_PTR       =  667
const SETDAC_INVALID_BOARDTYPE =  668

const CSO_DUMMY_CLOCK_DISABLE	           = 0
const CSO_DUMMY_CLOCK_TIMER		           = 1
const CSO_DUMMY_CLOCK_EXT_TRIGGER	       = 2
const CSO_DUMMY_CLOCK_TIMER_ON_TIMER_OFF = 3

"""
Auxiliary IO; this section had unsigned declaration in C API
"""
const AUX_OUT_TRIGGER               = U32(0)
const AUX_OUT_PACER                 = U32(2)
const AUX_OUT_BUSY                  = U32(4)
const AUX_OUT_CLOCK                 = U32(6)
const AUX_OUT_RESERVED              = U32(8)
const AUX_OUT_CAPTURE_ALMOST_DONE	= U32(10)
const AUX_OUT_AUXILIARY             = U32(12)
const AUX_OUT_SERIAL_DATA           = U32(14)
const AUX_OUT_TRIGGER_ENABLE        = U32(16)

const AUX_IN_TRIGGER_ENABLE         = U32(1)
const AUX_IN_DIGITAL_TRIGGER        = U32(3)
const AUX_IN_GATE					= U32(5)
const AUX_IN_CAPTURE_ON_DEMAND		= U32(7)
const AUX_IN_RESET_TIMESTAMP		= U32(9)
const AUX_IN_SLOW_EXTERNAL_CLOCK	= U32(11)
const AUX_IN_AUXILIARY				= U32(13)
const AUX_IN_SERIAL_DATA			= U32(15)

const AUX_INPUT_AUXILIARY			= AUX_IN_AUXILIARY
const AUX_INPUT_SERIAL_DATA			= AUX_IN_SERIAL_DATA

# AlazarSetExternalTriggerOperationForScanning
const STOS_OPTION_DEFER_START_CAPTURE	= U32(1)

# Data skipping
const SSM_DISABLE             =	U32(0)
const SSM_ENABLE              =	U32(1)

# Coprocessor
const CPF_REG_SIGNATURE       = U32(0)
const CPF_REG_REVISION        = U32(1)
const CPF_REG_VERSION         = U32(2)
const CPF_REG_STATUS          = U32(3)

const CPF_OPTION_DMA_DOWNLOAD = U32(1)

const CPF_DEVICE_UNKNOWN      = U32(0)
const CPF_DEVICE_EP3SL50      = U32(1)
const CPF_DEVICE_EP3SE260     = U32(2)

const LSB_DEFAULT             = U32(0)
const LSB_EXT_TRIG	          = U32(1)
const LSB_AUX_IN_0	          = U32(2)
const LSB_AUX_IN_1	          = U32(3)

"""
AutoDMA acquisitions flags
Control Flags for AutoDMA used in AlazarStartAutoDMA

Note: Not all AlazarTech devices are capable of dual-ported
acquisitions. Please see your board's specification for more
information.
"""
const ADMA_EXTERNAL_STARTCAPTURE  = U32(0x00000001)
const ADMA_ENABLE_RECORD_HEADERS  = U32(0x00000008)
const ADMA_SINGLE_DMA_CHANNEL     = U32(0x00000010)
const ADMA_ALLOC_BUFFERS          = U32(0x00000020)
const ADMA_TRADITIONAL_MODE       = U32(0x00000000)
const ADMA_CONTINUOUS_MODE        = U32(0x00000100)
const ADMA_NPT                    = U32(0x00000200)
const ADMA_TRIGGERED_STREAMING    = U32(0x00000400)
const ADMA_FIFO_ONLY_STREAMING    = U32(0x00000800)
const ADMA_INTERLEAVE_SAMPLES     = U32(0x00001000)
const ADMA_GET_PROCESSED_DATA     = U32(0x00002000)
const ADMA_DSP                    = U32(0x00004000)
const ADMA_ENABLE_RECORD_FOOTERS  = U32(0x00010000)

const ADMA_CLOCKSOURCE            = U32(0x00000001)
const ADMA_CLOCKEDGE              = U32(0x00000002)
const ADMA_SAMPLERATE             = U32(0x00000003)
const ADMA_INPUTRANGE             = U32(0x00000004)
const ADMA_INPUTCOUPLING          = U32(0x00000005)
const ADMA_IMPUTIMPEDENCE         = U32(0x00000006)       # typooo
const ADMA_EXTTRIGGERED           = U32(0x00000007)
const ADMA_CHA_TRIGGERED          = U32(0x00000008)
const ADMA_CHB_TRIGGERED          = U32(0x00000009)
const ADMA_TIMEOUT                = U32(0x0000000A)
const ADMA_THISCHANTRIGGERED      = U32(0x0000000B)
const ADMA_SERIALNUMBER           = U32(0x0000000C)
const ADMA_SYSTEMNUMBER           = U32(0x0000000D)
const ADMA_BOARDNUMBER            = U32(0x0000000E)
const ADMA_WHICHCHANNEL           = U32(0x0000000F)
const ADMA_SAMPLERESOLUTION       = U32(0x00000010)
const ADMA_DATAFORMAT             = U32(0x00000011)

"""Boards"""
const ATS850  = 1
const ATS310  = 2
const ATS330  = 3
const ATS855  = 4
const ATS315  = 5
const ATS335  = 6
const ATS460  = 7
const ATS860  = 8
const ATS660  = 9
const ATS665  = 10
const ATS9462 = 11
const ATS9434 = 12
const ATS9870 = 13
const ATS9350 = 14
const ATS9325 = 15
const ATS9440 = 16
const ATS9410 = 17
const ATS9351 = 18
const ATS9310 = 19
const ATS9461 = 20
const ATS9850 = 21
const ATS9625 = 22
const ATG6500 = 23
const ATS9626 = 24
const ATS9360 = 25
const AXI9870 = 26
const ATS9370 = 27
const ATU7825 = 28
const ATS9373 = 29
const ATS9416 = 30

"""Record average options"""
const CRA_MODE_DISABLE         = U32(0)
const CRA_MODE_ENABLE_FPGA_AVE = U32(1)
const CRA_OPTION_UNSIGNED      = U32(0)
const CRA_OPTION_SIGNED        = U32(1)

"""AlazarDSP constants"""
const DSP_WINDOW_NONE                = U32(0)
const DSP_WINDOW_HANNING             = U32(1)
const DSP_WINDOW_HAMMING             = U32(2)
const DSP_WINDOW_BLACKMAN            = U32(3)
const DSP_WINDOW_BLACKMAN_HARRIS     = U32(4)
const DSP_WINDOW_BARTLETT            = U32(5)
const NUM_DSP_WINDOW_ITEMS           = U32(6)

const DSP_MODULE_NONE                = U32(0x0FFFF)
const DSP_MODULE_FFT                 = U32(0x10000)
const DSP_MODULE_PCD                 = U32(0x10001)

const FFT_OUTPUT_FORMAT_U32          = U32(0x0000)
const FFT_OUTPUT_FORMAT_U16_LOG      = U32(0x0001)
const FFT_OUTPUT_FORMAT_U16_AMP2     = U32(0x0101)
const FFT_OUTPUT_FORMAT_U8_LOG       = U32(0x0002)
const FFT_OUTPUT_FORMAT_U8_AMP2      = U32(0x0102)
const FFT_OUTPUT_FORMAT_REAL_S32     = U32(0x0003)
const FFT_OUTPUT_FORMAT_IMAG_S32     = U32(0x0004)
const FFT_OUTPUT_FORMAT_FLOAT_AMP2   = U32(0x000A)
const FFT_OUTPUT_FORMAT_FLOAT_LOG    = U32(0x000B)
const FFT_OUTPUT_FORMAT_RAW_PLUS_FFT = U32(0x1000)

const FFT_FOOTER_NONE                = U32(0x0)
const FFT_FOOTER_NPT                 = U32(0x1)


text:
export DNN_CTC
export DNN_Batch_CTC
export RNN_Batch_CTC
export CRNN_Batch_CTC

"""
    DNN_CTC(p::Variable{Array{T}}, seq; blank=1, weight=1.0)

for case batchsize==1 for test case

# Inputs
`p`      : 2-D Variable, probability or weighted probability\n
`seq`    : 1-D Array, input sequence's label.\n
`weight` : weight for CTC loss\n

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function DNN_CTC(p::Variable{Array{T}}, seq; blank=1, weight=1.0) where T
    L = length(seq) * 2 + 1
    r, loglikely = CTC(ᵛ(p), seq, blank=blank)
    if p.backprop
        function DNN_CTC_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, DNN_CTC_Backward)
    end
    return loglikely / L
end

"""
    DNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T

a batch of concatenated input sequence is processed by neural networks into `p`

# Inputs
`p`         : 2-D Variable, probability or weighted probability\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : records each input sequence's length, like [20,17,...]\n
`weight`    : weight for CTC loss\n

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function DNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T
    batchsize = length(inputLengths)
    loglikely = zeros(T, batchsize)
    I, F = indexbounds(inputlens)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        span = I[b]:F[b]
        r[:,span], loglikely[b] = CTC(p.value[:,span], seqlabels[b], blank=blank)
        loglikely[b] /= length(seqlabels[b]) * 2 + 1
    end

    if p.backprop
        function DNN_Batch_CTC_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, DNN_Batch_CTC_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    RNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T

a batch of padded input sequence is processed by neural networks into `p`

# Inputs
`p`         : 3-D Variable with shape (featdims,timesteps,batchsize), probability or weighted probability\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : each input's length, like [19,97,...]\n
`weight`    : weight for CTC loss\n

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function RNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T
    batchsize = length(inputlens)
    loglikely = zeros(T, batchsize)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        Tᵇ = inputlens[b]
        Lᵇ = length(seqlabels[b])
        r[:,1:Tᵇ,b], loglikely[b] = CTC(p.value[:,1:Tᵇ,b], seqlabels[b], blank=blank)
        loglikely[b] /= Lᵇ * 2 + 1
    end

    if p.backprop
        function RNN_Batch_CTC_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, RNN_Batch_CTC_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    CRNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector) where T -> LogLikely

a batch of padded input sequence is processed by neural networks into `p`

# Inputs
`p`         : 3-D Variable with shape (featdims,timesteps,batchsize), probability or weighted probability.\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`weight`    : weight for CTC loss\n

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function CRNN_Batch_CTC(p::Variable{Array{T}}, seqlabels::Vector; blank=1, weight=1.0) where T
    featdims, timesteps, batchsize = size(p)
    loglikely = zeros(T, batchsize)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        r[:,:,b], loglikely[b] = CTC(p.value[:,:,b], seqlabels[b], blank=blank)
        loglikely[b] /= length(seqlabels[b]) * 2 + 1
    end

    if p.backprop
        function CRNN_Batch_CTC_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, CRNN_Batch_CTC_Backward)
    end
    return sum(loglikely)/batchsize
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Very simple implementation of some adaptive MCMC algorithms.
module AdaptiveMCMC

import Base.get

export adaptive_rwm, AdaptiveMetropolis, RobustAdaptiveMetropolis,
AdaptiveScalingMetropolis, AdaptiveScalingWithinAdaptiveMetropolis,
AdaptState, RWMState,
draw!, accept!, adapt!, adapt_rb!,
StepSize, PolynomialStepSize, RAMStepSize, get

using Random, LinearAlgebra, StaticArrays

# Abstract type of the state of an adaptive random walk Metropolis MCMC.
abstract type AdaptState end

include("linearalgebra_helpers.jl")
include("rwm.jl")       # Random walk Metropolis functionality
include("stepsize.jl")  # Step size sequences
include("ram.jl")       # RobustAdaptiveMetropolis
include("am.jl")        # AdaptiveMetropolis
include("asm.jl")       # AdaptiveScalingMetropolis and AdaptiveScalingWithinAdaptiveMetropolis
include("adaptive_rwm.jl")      # Adaptive random walk Metropolis interface `arwm`

end # module


text:
module TestForeach

using Test
using ThreadPools

include("util.jl")

@testset "foreaches" begin

    @testset "tforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += 1
        end
        tforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test primary
        @inferred tforeach(fn!, objs)
    end

    @testset "bforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += 1
        end
        bforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @inferred bforeach(fn!, objs)
    end

    @testset "qforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += 1
        end
        qforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test primary
        @inferred qforeach(fn!, objs)
    end

    @testset "qbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += 1
        end
        qbforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @inferred qbforeach(fn!, objs)
    end

    @testset "logtforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += 1
        end
        p = logtforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test length(p.recs) == N*2
        @test primary
        @inferred logtforeach(fn!, objs)
    end

    @testset "logbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += 1
        end
        p = logbforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test length(p.recs) == N*2
        @inferred logbforeach(fn!, objs)
    end

    @testset "logqforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += 1
        end
        p = logqforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test length(p.recs) == N*2
        @test primary
        @inferred logqforeach(fn!, objs)
    end

    @testset "logqbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += 1
        end
        p = logqbforeach(fn!, objs)
        @test [x.data for x in objs] == collect(2:N+1)
        @test length(p.recs) == N*2
        @inferred logqbforeach(fn!, objs)
    end

end

end # module


text:
using Turing
using Distributions

import Turing: ParticleContainer, weights, resample!, effectiveSampleSize, TraceC, TraceR, Trace, current_trace

function f()
  t = TArray(Float64, 1);
  t[1] = 0;
  while true
    rand(current_trace(), Normal(0,1))
    produce(0)
    rand(current_trace(), Normal(0,1))
    t[1] = 1 + t[1]
  end
end

pc = ParticleContainer{TraceC}(f)

push!(pc, TraceC(pc.model))
push!(pc, TraceC(pc.model))
push!(pc, TraceC(pc.model))

Base.@assert weights(pc)[1] == [1/3, 1/3, 1/3]
Base.@assert weights(pc)[2] ≈ log(3)
Base.@assert pc.logE ≈ log(1)

Base.@assert consume(pc) == log(1)

resample!(pc)
Base.@assert pc.num_particles == length(pc)
Base.@assert weights(pc)[1] == [1/3, 1/3, 1/3]
Base.@assert weights(pc)[2] ≈ log(3)
Base.@assert pc.logE ≈ log(1)
Base.@assert effectiveSampleSize(pc) == 3

Base.@assert consume(pc) ≈ log(1)
resample!(pc)
Base.@assert consume(pc) ≈ log(1)


text:
"""
$(SIGNATURES)

Default postprocess perfoms duality-based bound tightening on the `y`.
"""
function postprocess!(t::ExtensionType, m::GlobalOptimizer)
    if m._parameters.dbbt_depth > m._iteration_count
        variable_dbbt!(m._current_node, m._lower_lvd, m._lower_uvd,
                       m._lower_objective_value, m._global_upper_bound,
                       m._branch_variable_count)
    end
    return
end
postprocess!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = postprocess!(_ext_typ(m), m)

"""
$(SIGNATURES)

Checks to see if current node should be reprocessed.
"""
repeat_check(t::ExtensionType, m::GlobalOptimizer) = false
repeat_check(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = repeat_check(_ext_typ(m), m)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# All constants defined in the file are exported

for name in [
    ### Matrix types from cvdef.h ###
    :CV_CN_MAX,
    :CV_CN_SHIFT,
    :CV_DEPTH_MAX,

    :CV_8U,
    :CV_8S,
    :CV_16U,
    :CV_16S,
    :CV_32S,
    :CV_32F,
    :CV_64F,
    :CV_USRTYPE1,

    :CV_MAT_DEPTH_MASK,

    :CV_8UC1,
    :CV_8UC2,
    :CV_8UC3,
    :CV_8UC4,

    :CV_16UC1,
    :CV_16UC2,
    :CV_16UC3,
    :CV_16UC4,

    :CV_16SC1,
    :CV_16SC2,
    :CV_16SC3,
    :CV_16SC4,

    :CV_32SC1,
    :CV_32SC2,
    :CV_32SC3,
    :CV_32SC4,

    :CV_32FC1,
    :CV_32FC2,
    :CV_32FC3,
    :CV_32FC4,

    :CV_64FC1,
    :CV_64FC2,
    :CV_64FC3,
    :CV_64FC4,

    :CV_MAT_CN_MASK,
    :CV_MAT_TYPE_MASK,
    :CV_MAT_CONT_FLAG_SHIFT,
    :CV_MAT_CONT_FLAG,
    :CV_SUBMAT_FLAG_SHIFT,
    :CV_SUBMAT_FLAG,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), string(name, ";"))
    @eval begin
        global const $name = $ex
        export $name
    end
end

const BorderTypes = Cxx.CppEnum{Symbol("cv::BorderTypes"),UInt32}
const CmpTypes = Cxx.CppEnum{Symbol("cv::CmpTypes"),UInt32}
const DecompTypes = Cxx.CppEnum{Symbol("cv::DecompTypes"),UInt32}
const DftFlags = Cxx.CppEnum{Symbol("cv::DftFlags"),UInt32}
const GemmFlags = Cxx.CppEnum{Symbol("cv::GemmFlags"),UInt32}

for name in [
    :ACCESS_READ,
    :ACCESS_WRITE,
    :ACCESS_RW,
    :ACCESS_MASK,
    :ACCESS_FAST,

    :BORDER_CONSTANT,
    :BORDER_REPLICATE,
    :BORDER_REFLECT,
    :BORDER_WRAP,
    :BORDER_REFLECT_101,
    :BORDER_TRANSPARENT,
    :BORDER_REFLECT101,
    :BORDER_DEFAULT,
    :BORDER_ISOLATED,

    :CMP_EQ,
    :CMP_GT,
    :CMP_GE,
    :CMP_LT,
    :CMP_LE,
    :CMP_NE,

    :DECOMP_LU,
    :DECOMP_SVD,
    :DECOMP_EIG,
    :DECOMP_CHOLESKY,
    :DECOMP_QR,
    :DECOMP_NORMAL,

    :DFT_INVERSE,
    :DFT_SCALE,
    :DFT_ROWS,
    :DFT_COMPLEX_OUTPUT,
    :DFT_REAL_OUTPUT,
    :DCT_INVERSE,
    :DCT_ROWS,

    :GEMM_1_T,
    :GEMM_2_T,
    :GEMM_3_T,

    :NORM_INF,
    :NORM_L1,
    :NORM_L2,
    :NORM_L2SQR,
    :NORM_HAMMING,
    :NORM_HAMMING2,
    :NORM_TYPE_MASK,
    :NORM_RELATIVE,
    :NORM_MINMAX,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), string("cv::", name, ";"))
    @eval global const $name = $ex
end


text:
using IntervalContractors
import Documenter

Documenter.makedocs(
    modules = [IntervalContractors],
    format = :html,
    sitename = "IntervalContractors.jl",
    root = joinpath(dirname(dirname(@__FILE__)), "docs"),
    pages = Any["Home" => "index.md"],
    strict = true,
    linkcheck = true,
    checkdocs = :exports,
    authors = "David Sanders"
)

Documenter.deploydocs(
    julia = "0.6",
    repo = "github.com/JuliaIntervals/IntervalContractors.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
# Dual infeasible test problems
# Problems with m1 equality constraints, m3-dimensional second order constraint and m4-dimensional psd constraint
# x1 is unbounded below --> hence problem is dual infeasible

nn = 1
rng = Random.MersenneTwister(1313)

sum_detected = 0
@testset "Dual infeasible QP problems - Testset 2" begin
  for iii = 1:nn

   # choose size of problem
    n = rand(rng,10:50)
    m1 = rand(rng,2:10) #equality constraints
    m2 = 1 #inequality constraints
    m3 = rand(rng,3:10) # second order cone constraint
    r = rand(rng,4:10)
    m4 = r^2 # one sdp cone constraint of size r

    m = m1+m2+m3+m4

    # make problem primal feasible
    A = sprandn(rng,m,n,0.8)*50
    xtrue = rand(rng,n,1)*50
    s1 = zeros(m1,1)
    s2 = rand(rng)
    s3 = randn(rng,m3-1,1)
    s3 = [norm(s3,2)+1;s3]
    s4 = vec(generate_pos_def_matrix(rng, r))
    s = [s1;s2;s3;s4]

    # make problem unbounded in x1
    A[:,1] .= 0
    A[m1+1,:] = [-1;zeros(n-1)]
    b = A*xtrue+s
    b[m1+1] = 0
    b = vec(b)

    P = spzeros(n,n)
    q = vec([-1;randn(rng,n-1)])

    A1 = -A[1:m1,:]
    A2 = -A[m1+1:m1+m2,:]
    A3 = -A[m1+m2+1:m1+m2+m3,:]
    A4 = -A[m1+m2+m3+1:end,:]
    b1 = b[1:m1]
    b2 = b[m1+1:m1+m2]
    b3 = b[m1+m2+1:m1+m2+m3]
    b4 = b[m1+m2+m3+1:end]

    cs1 = COSMO.Constraint(A1,b1,COSMO.ZeroSet)
    cs2 = COSMO.Constraint(A2,b2,COSMO.Nonnegatives)
    cs3 = COSMO.Constraint(A3,b3,COSMO.SecondOrderCone)
    cs4 = COSMO.Constraint(A4,b4,COSMO.PsdCone)

    settings = COSMO.Settings(max_iter=10000,eps_abs = 1e-5,eps_rel=1e-5)
    model = COSMO.Model()
    assemble!(model,P,q,[cs1;cs2;cs3;cs4], settings = settings)
    res = COSMO.optimize!(model);

    @test res.status == :Dual_infeasible
  end
end
nothing


text:
export ModelKitSystem

"""
    ModelKitSystem(F:System)

Construct a system from the given [`System`](@ref) `F`.
The difference to `F` is that this compiles a straight line programm for the fast
evaluation of `F` and that `ModelKitSystem <: AbstractSystem`.
"""
struct ModelKitSystem{S} <: AbstractSystem
    system::ModelKit.CompiledSystem{S}
end
ModelKitSystem(F::ModelKit.System; optimizations::Bool = true) =
    ModelKitSystem(ModelKit.compile(F; optimizations = optimizations))

Base.size(F::ModelKitSystem) = size(F.system)
ModelKit.variables(F::ModelKitSystem) = variables(F.system.system)
ModelKit.parameters(F::ModelKitSystem) = parameters(F.system.system)
ModelKit.variable_groups(F::ModelKitSystem) = variable_groups(F.system.system)

function Base.show(io::IO, F::ModelKitSystem)
    println(io, typeof(F), ":")
    println(io, F.system)
end

(F::ModelKitSystem)(x, p = nothing) = F.system.system(x, p)
evaluate!(u, F::ModelKitSystem, x, p = nothing) = ModelKit.evaluate!(u, F.system, x, p)
evaluate_and_jacobian!(u, U, F::ModelKitSystem, x, p = nothing) =
    ModelKit.evaluate_and_jacobian!(u, U, F.system, x, p)
taylor!(u, v::Val, F::ModelKitSystem, tx, p = nothing) =
    ModelKit.taylor!(u, v, F.system, tx, p)


text:
bin_edges_v = 0:1:20000
h = Histogram(bin_edges_v)

@test all(Bool[isapprox(edges(h)[i], bin_edges_v[i]) for i = 1:length(bin_edges_v)])
bin_centers_v = (bin_edges_v[2:end]+bin_edges_v[1:end-1])*0.5
@test all(Bool[isapprox(midpoints(h)[i], bin_centers_v[i]) for i = 1:length(bin_edges_v)-1])
@test isapprox(binsize(h),1.0)

addcounts!(h, ones(100)*10)
@test counted(h)==100
@test donethru(h) == 100
@test sum(counts(h))==100

addcounts!(h, ones(100)*-10)
@test counted(h)==100
@test misses(h) == 100
@test donethru(h) == 200

s = Bool[1,1,1,1,0,0,0,0]
events = [-1,2,3,4,5,6,7,8]
update_histogram!(h, s, events)
@test counted(h) == 103
@test donethru(h) == 208
@test misses(h) == 105

@test counts(h)[10]==100
@test counts(h)[2]==1

# test interface to KernelDensity
uvh = uvhist(1:100, collect(1:100))
@test edges(uvh)==1:100
@test density(uvh)==ones(99)
conv(h, Normal(0,1))


text:
using BinDeps
@BinDeps.setup

version = "4.25.1"
url = "ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib-$version.tar.bz2"

wcs = library_dependency("libwcs")
provides(Sources, URI(url), wcs, unpacked_dir="wcslib-$version")
depsdir = BinDeps.depsdir(wcs)
builddir = joinpath(depsdir, "src/wcslib-$version")
prefix = joinpath(depsdir, "usr")
configopts = ["--disable-fortran", "--without-cfitsio", "--without-pgplot"]
@unix_only libfilename = "libwcs.so.$version"
@osx_only libfilename = "libwcs.dylib"
provides(BuildProcess,
         (@build_steps begin
            GetSources(wcs)
            @build_steps begin
              ChangeDirectory(builddir)
              FileRule(joinpath(prefix,"lib",libfilename),
                       @build_steps begin
                         `./configure --prefix=$prefix $configopts`
                         `make install`
                       end)
            end
          end),
         wcs)

@BinDeps.install [:libwcs => :libwcs]