text:
##### Diffusion approximation for a turbid parallelepiped #####

#=
@article{kienle2005light,
  title={Light diffusion through a turbid parallelepiped},
  author={Kienle, Alwin},
  journal={JOSA A},
  volume={22},
  number={9},
  pages={1883--1888},
  year={2005},
  publisher={Optical Society of America}
}
=#

### TD Fluence ###
"""
    fluence_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

Compute the time-domain fluence in a parallelepiped [lx, ly, lz] with Eqn. 6 Kienle 05.

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `rd::Array{Float64,1}`: target location within medium [x,y,z] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs::Array{Float64,1}`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L::Array{Float64,1}`: the dimensions [lx, ly, lz] of the parallelepied

# Examples
julia> fluence_DA_paralpip_TD(0:1:5, [0.1,10.0], 1.0, 1.0, [2.5,5.0,1.0], [5.0,5.0], [10.0,10.0,10.0])
"""
function fluence_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})
	n = nmed/ndet
    μa = β[1]
    μsp = β[2]
    D = 1/3μsp
	ν = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x = rd[1]
	y = rd[2]
	z = rd[3]
	xu = rs[1]
	yu = rs[2]
	lx = L[1]
	ly = L[2]
	lz = L[3]

	Rt1 = Array{Float64}(undef, length(t))
	Rt2 = zeros(Float64, length(t))
	Rt3 = zeros(Float64, length(t))
	Rt4 = zeros(Float64, length(t))

	if n == 1.0
		A = 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zo = 1/μsp
    zb = 2A*D

    Threads.@threads for n in eachindex(t)

        Rt1[n] = @. ν*exp(-μa*ν*t[n])/((4π*D*ν*t[n])^(3/2))

		  for m in xs
			    x1l = 2m.*lx .+ 4m.*zb .+ xu
				x2l = 2m.*lx + (4m-2).*zb .- xu
				y1m = 2m.*ly .+ 4m.*zb .+ yu
				y2m = 2m.*ly .+ (4m.-2).*zb .- yu
				z1n = 2m.*lz .+ 4m.*zb .+ zo
				z2n = 2m.*lz .+ (4m.-2).*zb .- zo

    	        Rt2[n] += exp(-(x-x1l)^2 / (4D*ν*t[n])) - exp(-(x-x2l)^2 / (4D*ν*t[n]))
	            Rt3[n] += exp(-(y-y1m)^2 / (4D*ν*t[n])) - exp(-(y-y2m)^2 / (4D*ν*t[n]))
	            Rt4[n] += exp(-(z-z1n)^2 /(4D*ν*t[n])) - exp(-(z-z2n)^2 /(4D*ν*t[n]))
          end

          Rt1[n] *= (Rt2[n]*Rt3[n]*Rt4[n])

        if isnan(Rt1[n])
            Rt1[n] = 0
        end

    end

	return Rt1
end

### TD Reflectance ###
"""
    refl_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

Compute the time-domain reflectance from a parallelepiped [lx, ly, lz] applying Fick's law to Eqn. 6 Kienle 05.

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `rd::Array{Float64,1}`: the location of the detector [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs::Array{Float64,1}`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L::Array{Float64,1}`: the dimenensions [lx, ly, lz] of the parallelepied

# Examples
julia> refl_DA_paralpip_TD(0:1:5, [0.1,10.0], 1.0, 1.0, [2.5,5.0], [5.0,5.0], [10.0,10.0,10.0])
"""
function refl_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})
	n = nmed/ndet
	μa = β[1]
	μsp = β[2]
	D = 1/3μsp
	ν = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x = rd[1]
	y = rd[2]
	xu = rs[1]
	yu = rs[2]
	lx = L[1]
	ly = L[2]
    lz = L[3]

    Rt1 = Array{Float64}(undef, length(t))
	Rt2 = zeros(Float64, length(t))
	Rt3 = zeros(Float64, length(t))
	Rt4 = zeros(Float64, length(t))

	if n == 1.0
		A= 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zo = 1/μsp
	zb = 2A*D

    Threads.@threads for n in eachindex(t)

        Rt1[n] = @. exp(-μa*ν*t[n])/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

		  for m in xs
			    x1l = 2m.*lx .+ 4m.*zb .+ xu
				x2l = 2m.*lx + (4m-2).*zb .- xu
				y1m = 2m.*ly .+ 4m.*zb .+ yu
				y2m = 2m.*ly .+ (4m.-2).*zb .- yu
				z1n = 2m.*lz .+ 4m.*zb .+ zo
				z2n = 2m.*lz .+ (4m.-2).*zb .- zo

				Rt2[n] += exp(-(x-x1l)^2 / (4D*ν*t[n])) - exp(-(x-x2l)^2 / (4D*ν*t[n]))

				Rt3[n] += exp(-(y-y1m)^2 / (4D*ν*t[n])) - exp(-(y-y2m)^2 / (4D*ν*t[n]))

				Rt4[n] += z1n*exp(-z1n^2 /(4D*ν*t[n])) - z2n'*exp(-z2n^2 /(4D*ν*t[n]))
          end

        Rt1[n] *= (Rt2[n]*Rt3[n]*Rt4[n])

        if isnan(Rt1[n])
            Rt1[n] = 0
        end

    end

	return Rt1
end

### TD Transmittance ###
"""
    trans_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

Compute the time-domain transmittance from a parallelepiped [lx, ly, lz] with Eqn. 8 Kienle 05.

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `rd::Array{Float64,1}`: the location of the detector [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs::Array{Float64,1}`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L::Array{Float64,1}`: the dimenensions [lx, ly, lz] of the parallelepied

# Examples
julia> trans_DA_paralpip_TD(0:1:5, [0.1,10.0], 1.0, 1.0, [2.5,5.0], [5.0,5.0], [10.0,10.0,10.0])
"""
function trans_DA_paralpip_TD(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})
	n = nmed/ndet
    μa = β[1]
    μsp = β[2]
    D = 1/3μsp
	ν = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x = rd[1]
	y = rd[2]
	z = L[3]
	xu = rs[1]
	yu = rs[2]
	lx = L[1]
	ly = L[2]
	lz = L[3]

	Rt1 = Array{Float64}(undef, length(t))
	Rt2 = zeros(Float64, length(t))
	Rt3 = zeros(Float64, length(t))
	Rt4 = zeros(Float64, length(t))

	if n == 1.0
		A = 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zo = 1/μsp
    zb = 2A*D

    Threads.@threads for n in eachindex(t)

        Rt1[n] = @. exp(-μa*ν*t[n])/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

		  for m in xs
			    x1l = 2m.*lx .+ 4m.*zb .+ xu
				x2l = 2m.*lx + (4m-2).*zb .- xu
				y1m = 2m.*ly .+ 4m.*zb .+ yu
				y2m = 2m.*ly .+ (4m.-2).*zb .- yu
				z1n = 2m.*lz .+ 4m.*zb .+ zo
				z2n = 2m.*lz .+ (4m.-2).*zb .- zo

    	        Rt2[n] += exp(-(x-x1l)^2 / (4D*ν*t[n])) - exp(-(x-x2l)^2 / (4D*ν*t[n]))
	            Rt3[n] += exp(-(y-y1m)^2 / (4D*ν*t[n])) - exp(-(y-y2m)^2 / (4D*ν*t[n]))
	            Rt4[n] += (lz - z1n)*exp(-(lz-z1n)^2 /(4D*ν*t[n])) - (lz - z2n)*exp(-(lz-z2n)^2 /(4D*ν*t[n]))
          end

          Rt1[n] *= (Rt2[n]*Rt3[n]*Rt4[n])

        if isnan(Rt1[n])
            Rt1[n] = 0
        end

    end

	return Rt1
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
#=
    Used from https://github.com/lanl-ansi/Juniper.jl
=#

function init_original_problem!(op::OriginalProblem, model::MOI.AbstractOptimizer)
    num_variables = length(model.variable_info)
    num_linear_le_constraints = length(model.linear_le_constraints)
    num_linear_ge_constraints = length(model.linear_ge_constraints)
    num_linear_eq_constraints = length(model.linear_eq_constraints)
    num_quadratic_le_constraints = length(model.quadratic_le_constraints)
    num_quadratic_ge_constraints = length(model.quadratic_ge_constraints)
    num_quadratic_eq_constraints = length(model.quadratic_eq_constraints)

    op.status = MOI.OPTIMIZE_NOT_CALLED
    op.relaxation_status = MOI.OPTIMIZE_NOT_CALLED
    op.has_nl_objective = model.nlp_data.has_objective
    op.nlp_evaluator = model.nlp_data.evaluator
    op.objective = model.objective

    op.objval = NaN
    op.best_bound = NaN
    op.solution = fill(NaN, num_variables)
    op.nsolutions = 0
    op.solutions = []
    op.num_disc_var = 0
    op.nintvars = 0
    op.nbinvars = 0
    op.nnodes = 1 # is set to one for the root node

    op.nlp_constr_type = []

    op.start_time = time()

    op.nl_solver = model.options.nl_solver
    nl_vec_opts = Vector{Pair}()
    if isa(op.nl_solver, MOI.OptimizerWithAttributes)
        for arg in model.options.nl_solver.params
            push!(nl_vec_opts, arg)
        end
    else
        op.nl_solver = JuMP.optimizer_with_attributes(op.nl_solver)
    end
    op.nl_solver_options = nl_vec_opts

    if model.options.mip_solver !== nothing
        op.mip_solver = model.options.mip_solver
        mip_vec_opts = Vector{Pair}()
        if isa(op.mip_solver, MOI.OptimizerWithAttributes)
            for arg in model.options.mip_solver.params
                push!(mip_vec_opts, arg)
            end
        else
            op.mip_solver = JuMP.optimizer_with_attributes(op.mip_solver)
        end
        op.mip_solver_options = mip_vec_opts
    end
    op.options = model.options
    if model.sense == MOI.MIN_SENSE
        op.obj_sense = :Min
    else
        op.obj_sense = :Max
    end
    op.l_var = info_array_of_variables(model.variable_info, :lower_bound)
    op.u_var = info_array_of_variables(model.variable_info, :upper_bound)
    integer_bool_arr = info_array_of_variables(model.variable_info, :is_integer)
    binary_bool_arr = info_array_of_variables(model.variable_info, :is_binary)
    primal_start_arr = info_array_of_variables(model.variable_info, :start)
    op.primal_start = primal_start_arr
    op.nintvars = sum(integer_bool_arr)
    op.nbinvars = sum(binary_bool_arr)
    op.num_disc_var = sum(integer_bool_arr)+sum(binary_bool_arr)
    op.num_var = length(model.variable_info)
    op.var_type = [:Cont for i in 1:op.num_var]
    op.var_type[integer_bool_arr .== true] .= :Int
    op.var_type[binary_bool_arr .== true] .= :Bin
    op.disc2var_idx = zeros(op.num_disc_var)
    op.var2disc_idx = zeros(op.num_var)
    int_i = 1
    for i=1:op.num_var
        if op.var_type[i] != :Cont
            op.disc2var_idx[int_i] = i
            op.var2disc_idx[i] = int_i
            int_i += 1
        end
    end
    op.num_l_constr = num_linear_le_constraints+num_linear_ge_constraints+num_linear_eq_constraints
    op.num_q_constr = num_quadratic_le_constraints+num_quadratic_ge_constraints+num_quadratic_eq_constraints
    op.num_nl_constr = length(model.nlp_data.constraint_bounds)
    op.num_constr = op.num_l_constr+op.num_q_constr+op.num_nl_constr
end

function create_original_model!(optimizer::MOI.AbstractOptimizer, op::OriginalProblem; fix_start=false)
    ps = op.options.log_levels

    op.model = Model()
    lb = op.l_var
    ub = op.u_var
    # all continuous we solve relaxation first
    @variable(op.model, lb[i] <= x[i=1:op.num_var] <= ub[i])

    for i=1:op.num_var
        JuMP.set_start_value(x[i], op.primal_start[i])
    end

    register_functions!(op.model, op.options.registered_functions)

    # TODO check whether it is supported
    if optimizer.nlp_data.has_objective
        obj_expr = MOI.objective_expr(optimizer.nlp_data.evaluator)
        obj_expr = expr_dereferencing(obj_expr, op.model)
        try
            JuMP.set_NL_objective(op.model, optimizer.sense, obj_expr)
        catch
            error("Have you registered a function? Then please register the function also for Juniper see: \n
            https://lanl-ansi.github.io/Juniper.jl/stable/options/#registered_functions%3A%3AUnion%7BNothing%2CVector%7BRegisteredFunction%7D%7D-%5Bnothing%5D-1")
        end
    elseif optimizer.objective !== nothing
        MOI.set(op.model, MOI.ObjectiveFunction{typeof(optimizer.objective)}(), optimizer.objective)
        MOI.set(op.model, MOI.ObjectiveSense(), optimizer.sense)
    end

    backend = JuMP.backend(op.model);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    qlc = optimizer.quadratic_le_constraints
    qgc = optimizer.quadratic_ge_constraints
    qec = optimizer.quadratic_eq_constraints
    for constr_type in [llc, lgc, lec, qlc, qgc, qec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    for i in 1:op.num_nl_constr
        constr_expr = MOI.constraint_expr(optimizer.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, op.model)
        push!(op.nlp_constr_type, constr_expr.args[1])
        try
            JuMP.add_NL_constraint(op.model, constr_expr)
        catch
            error("Have you registered a function? Then please register the function also for Juniper see: \n
            https://lanl-ansi.github.io/Juniper.jl/stable/options/#registered_functions%3A%3AUnion%7BNothing%2CVector%7BRegisteredFunction%7D%7D-%5Bnothing%5D-1")
        end
    end

    op.x = x
end

function fix_primal_start!(op::OriginalProblem)
    lb = op.l_var
    ub = op.u_var
    x = op.x
    for i=1:op.num_var
        if op.var_type[i] != :Cont && lb[i] <= op.primal_start[i] <= ub[i]
            JuMP.set_lower_bound(op.x[i], op.primal_start[i])
            JuMP.set_upper_bound(op.x[i], op.primal_start[i])
        else
            JuMP.set_start_value(x[i], op.primal_start[i])
        end
    end
end

function unfix_primal_start!(op::OriginalProblem)
    lb = op.l_var
    ub = op.u_var
    x = op.x
    for i=1:op.num_var
        JuMP.set_lower_bound(op.x[i], lb[i])
        JuMP.set_upper_bound(op.x[i], ub[i])
        JuMP.set_start_value(x[i], op.primal_start[i])
    end
end

function solve_root_incumbent_model(op::OriginalProblem)
    status, backend = optimize_get_status_backend(op.model; solver=op.nl_solver)
    incumbent = nothing
    ps = op.options.log_levels
    if state_is_optimal(status; allow_almost=op.options.allow_almost_solved)
        # set incumbent
        objval = MOI.get(backend, MOI.ObjectiveValue())
        solution = JuMP.value.(op.x)
        if are_type_correct(solution, op.var_type, op.disc2var_idx, op.options.atol)
            if only_almost_solved(status) && op.options.allow_almost_solved_integral
                @warn "Start value incumbent only almost locally solved. Disable with `allow_almost_solved_integral=false`"
            end
        end
        #incumbent = Incumbent(objval, solution, only_almost_solved(status))
    end
    return incumbent
end

ext: jl

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
include("intcode.jl")

# Part 1: run the intcode program on all points in the 50x50 starting grid and count 1's
prog = open(f -> parse.(Int64, split(readline(f), ",")), "day19_input.txt")
result = map(x -> run_intcode!(copy(prog), [x[1], x[2]])[1], CartesianIndices((0:49, 0:49)))
println("Number of on pixels: $(sum(result))")

# Part 2: find top left coord of part of beam covering 100x100 area. Estimate first and last y as fn of x
approx_start_factor = 38/50
approx_end_factor = 49/50
test_point(prog, x, y) = run_intcode!(copy(prog), [x, y])[1]
function find_exact_start_end(prog, x)
    ys = Int32(floor(x * approx_start_factor)) - 1
    while test_point(prog, x, ys) == 0 ys += 1 end
    ye = Int32(ceil(x * approx_end_factor)) + 1
    while test_point(prog, x, ye) == 0 ye -= 1 end
    return ys, ye
end
# Look for the first coordinate which can fit a 100x100 square in
for x in 700:1200
    ys1, ye1 = find_exact_start_end(prog, x)
    ys2, ye2 = find_exact_start_end(prog, x + 99)
    if ye1 - ys2 >= 99
        println("First coord for square: ($x, $ys2). Code: $(x * 10000 + ys2)")
        break
    end
end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
## Direct calls

function adjoint_sensitivities(sol,args...;
                                  sensealg=InterpolatingAdjoint(),
                                  kwargs...)
  _adjoint_sensitivities(sol,sensealg,args...;kwargs...)
end

function _adjoint_sensitivities(sol,sensealg,alg,g,t=nothing,dg=nothing;
                                   abstol=1e-6,reltol=1e-3,
                                   checkpoints=sol.t,
                                   corfunc_analytical=nothing,
                                   callback = nothing,
                                   kwargs...)

  if sol.prob isa ODEProblem
    adj_prob = ODEAdjointProblem(sol,sensealg,g,t,dg,checkpoints=checkpoints,
                                 callback = callback,
                                 abstol=abstol,reltol=reltol)

  elseif sol.prob isa SDEProblem
    adj_prob = SDEAdjointProblem(sol,sensealg,g,t,dg,checkpoints=checkpoints,
                                 callback = callback,
                                 abstol=abstol,reltol=reltol,
                                 corfunc_analytical=corfunc_analytical)
  elseif sol.prob isa RODEProblem
    adj_prob = RODEAdjointProblem(sol,sensealg,g,t,dg,checkpoints=checkpoints,
                                callback = callback,
                                abstol=abstol,reltol=reltol,
                                corfunc_analytical=corfunc_analytical)
  else
    error("Continuous adjoint sensitivities are only supported for ODE/SDE/RODE problems.")
  end

  tstops = ischeckpointing(sensealg, sol) ? checkpoints : similar(sol.t, 0)
  adj_sol = solve(adj_prob,alg;
                  save_everystep=false,save_start=false,saveat=eltype(sol[1])[],
                  tstops=tstops,abstol=abstol,reltol=reltol,kwargs...)

  p = sol.prob.p
  l = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(sol.prob.p)
  du0 = -adj_sol[end][1:length(sol.prob.u0)]

  if eltype(sol.prob.p) <: real(eltype(adj_sol[end]))
    dp = real.(adj_sol[end][(1:l) .+ length(sol.prob.u0)])'
  elseif p === nothing || p === DiffEqBase.NullParameters()
    dp = nothing
  else
    dp = adj_sol[end][(1:l) .+ length(sol.prob.u0)]'
  end

  du0,dp
end

function _adjoint_sensitivities(sol,sensealg::SteadyStateAdjoint,alg,g,dg=nothing;
                                   abstol=1e-6,reltol=1e-3,
                                   kwargs...)
  SteadyStateAdjointProblem(sol,sensealg,g,dg;kwargs...)
end

function _adjoint_sensitivities(sol,sensealg::SteadyStateAdjoint,alg;
                                   g=nothing,dg=nothing,
                                   abstol=1e-6,reltol=1e-3,
                                   kwargs...)
  SteadyStateAdjointProblem(sol,sensealg,g,dg;kwargs...)
end

function second_order_sensitivities(loss,prob,alg,args...;
                                    sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                                    kwargs...)
  _second_order_sensitivities(loss,prob,alg,sensealg,args...;kwargs...)
end

function second_order_sensitivity_product(loss,v,prob,alg,args...;
                                          sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                                          kwargs...)
  _second_order_sensitivity_product(loss,v,prob,alg,sensealg,args...;kwargs...)
end


text:
function rectfunc(x::T, x_c::T)::T where T
    if abs(x) <= x_c
        return one(T)
    end

    return zero(T)
end

function raisedcosinefunction(x::T, β, 𝑇)::T where T
    return raisedcosinefunction(x,convert(T,β),convert(T,𝑇))
end

"""
    raisedcosinefunc(x::T, β::T, 𝑇::T)::T
Let low := (1-β)/(2𝑇),
    |x| := abs(x),
    high := (1+β)/(2𝑇).

Returns 1 if |x| ≤ low,
Returns 0.5*(1+cos(π*𝑇/β *(|x|-low))) if low ≤ |x| ≤ high.
Returns 0 otherwise.
"""
function raisedcosinefunction(x::T, β::T, 𝑇::T)::T where T
    two = one(T) + one(T)
    low = (one(T)-β)/(two*𝑇)
    high = (one(T)+β)/(two*𝑇)

    if abs(x) <= low
        return one(T)
    elseif low < abs(x) <= high
        return one(T)/two*( one(T) + cos(π*𝑇/β *(abs(x)-low)) )
    end

    return zero(T)
end

# transition band is cosine. This allows satisfies the power complementary condition.
# cos( (yp-y)/(yp-ys) * π/2 ) for a falling edge. yp is pass band, ys is stop band.
function cosinetransitionlowpassfunction(y::T, yp::T, ys::T)::T where T
    abs_y = abs(y)

    if abs_y <= yp
        return one(T)
    elseif yp < abs_y <= ys
        #return cos( (yp-y)/(ys-yp)*π/2 )
        return cos( (yp-abs_y)/(ys-yp)*π/2 )
    end

    return zero(T)
end

function cosinetransitionhighpassfunction(x::T, xp::T, xs::T)::T where T
    abs_x = abs(x)

    if abs_x <= xs
        return zero(T)
    elseif xs < abs_x <= xp
        #return cos( (xp-x)/(xp-xs)*π/2 )
        return cos( (xp-abs_x)/(xp-xs)*π/2 )
    end

    return one(T)
end

# rs is rising edge's stopband.
# rp is rising edge's passband.
# fp is falling edge's passband.
# fs is falling edge's stopband.
function cosinetransitionbandpassfunction(y::T, rs, rp, fp, fs)::T where T
    abs_y = abs(y)

    # rising edge.
    if abs_y <= rs
        return zero(T)
    elseif rs < abs_y <= rp
        #return cos( (rp-y)/(rp-rs)*π/2 )
        return cos( (rp-abs_y)/(rp-rs)*π/2 )
    end

    # falling edge.
    if rp <= abs_y <= fp
        return one(T)
    elseif fp < abs_y <= fs
        #return cos( (fp-y)/(fs-fp)*π/2 )
        return cos( (fp-abs_y)/(fs-fp)*π/2 )
    end

    return zero(T)
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
# void cufd(double *res, double *grad_Cp, double *grad_Cs, double *grad_Den,
#           double *grad_stf, const double *Cp, const double *Cs,
#           const double *Den, const double *stf, int calc_id, const int gpu_id,
#           int group_size, const int *shot_ids, const string para_fname);

function obscalc(cp,cs,den,stf,shot_ids,para_fname)
    m, n = size(cp)
    res = zeros(1)
    grad_Cp = zeros(m, n)
    grad_Cs = zeros(m, n)
    grad_Den = zeros(m, n)
    grad_stf = zeros(size(stf)...)
    calc_id = Int32(2)
    gpu_id = Int32(0)
    group_size = length(shot_ids)
    ccall((:cufd, "./Src/build/libCUFD.so"), Cvoid, (Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
            Ref{Cdouble}, Ref{Cdouble}, Cint, Cint, Cint, Ref{Cint}, Cstring),
        res, grad_Cp, grad_Cs, grad_Den, grad_stf, cp, cs, den, stf, calc_id, gpu_id, group_size, shot_ids,
        para_fname)
end
nz = 134
nx = 384
cp = 2500ones(nz, nx)
cs = zeros(nz, nx)
den = 1000ones(nz, nx)
shot_ids = Int32[0 1]
para_fname = "/home/lidongzh/TwoPhaseFlowFWI/Ops/FWI/Src/params/Par_file_obs_data.json"

src = Matrix{Float64}(undef, 1, 2001)
src[1,:] = Float64.(reinterpret(Float32, read("/home/lidongzh/TwoPhaseFlowFWI/Ops/FWI/Src/params/ricker_10Hz.bin")))
stf = repeat(src, outer=30)

obscalc(cp,cs,den,stf,shot_ids,para_fname)


text:
import FinanceLib
using BenchmarkTools

import Dates

# @code_warntype FinanceLib.tMul(0.06/12, -120.0)
# @code_warntype FinanceLib.tMul(0.06, -10.0, 12.0)
# @code_warntype FinanceLib.rateGwth(7.35, 8.52, 5.0)
# @code_warntype FinanceLib.periodGwth(100.0,50.0,0.07)

# @code_warntype FinanceLib.pv(10_000_000., 0.09, 5.0)
# @code_warntype FinanceLib.pv(12_704_891.6109538, 0.06, 4.0, 12.0)
# @code_warntype FinanceLib.pvc(11_735.108709918102, 0.08, 2.0)

# @code_warntype FinanceLib.fv(6_499_313.862983453, 0.09, 5.0)
# @code_warntype FinanceLib.fv(10_000_000.0, 0.06, 4.0, 12.0)
# @code_warntype FinanceLib.fvc(10_000., 0.08, 2.0)

# @code_warntype FinanceLib.pvAnnuity(1000.0, 0.12, 5.0)
# @code_warntype FinanceLib.pvAnnuity(7.33764573879378, 0.08, 30.0, 12.0)

# @code_warntype FinanceLib.pvAnnuity(100.0, 0.05)

# @code_warntype FinanceLib.fvAnnuity(1000.0, 0.05, 5.0)
# @code_warntype FinanceLib.fvAnnuity(2000.0, 0.24, 5.0, 3.0)

# @code_warntype FinanceLib.pmt(3_604.776202345007, 0.12, 5.0)
# @code_warntype FinanceLib.pmt(1000.0, 0.08, 30.0, 12.0)

# @code_warntype FinanceLib.fmt(5_525.631250000007, 0.05, 5.0)
# @code_warntype FinanceLib.fmt(54_304.2278549568, 0.24, 5.0, 3.0)

# @code_warntype FinanceLib.pv(FinanceLib.pvAnnuity(10.0^6,.05,30.0),0.05,9.0)

# @code_warntype FinanceLib.effR(0.08, 2.0)
# @code_warntype FinanceLib.nomR(0.08, 4)

# @code_warntype FinanceLib.expR(0.08)

# @code_warntype FinanceLib.expR(0.07,4)

# @code_warntype FinanceLib.npv(0.08, [0.25,6.25,3.5,4.5,1.25],
# [-6.25,1.2,1.25,3.6,2.5], 0.45)
# @code_warntype FinanceLib.npv(0.08, zip([0.25,6.25,3.5,4.5,1.25],
# [-6.25,1.2,1.25,3.6,2.5]), 0.45)

# @code_warntype FinanceLib.irr([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25])
# @code_warntype FinanceLib.irr(zip([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25]))

# @code_warntype FinanceLib.yearFrac(Dates.Date(2027,2,12), Dates.Date(2018,2,12))
# @code_warntype FinanceLib.invYearFrac(Dates.Date(2027,2,12), -8.999315537303216)

# @code_warntype FinanceLib.xnpv(0.08, [Dates.Date(2012,2,25), Dates.Date(2012,6,28),
#   Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-15, 5, 25, -10, 50], Dates.Date(2012,1,10) )

# @code_warntype FinanceLib.xnpv(0.08, zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28), Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-15, 5, 25, -10, 50.]), Dates.Date(2012,1,10) )

# @code_warntype FinanceLib.xirr([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
#   Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-115, 5, 25, -10, 200] )

# @code_warntype FinanceLib.xirr(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
#   Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-115, 5, 25, -10, 200]) )

# td = collect(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
#   Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-115, 5, 25, -10, 200]))

# @code_warntype FinanceLib.dateToPeriodSeries(Dates.Date(2010,05,12), td)

# td1 = FinanceLib.dateToPeriodSeries(Dates.Date(2010,05,12), td)

# @code_warntype FinanceLib.irr(td1)
# @code_warntype FinanceLib.xirr(td)

# @code_warntype FinanceLib.xirr(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
#   Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
#   [-115, 5, 25, -10, 200]) )

# @code_warntype FinanceLib.npv(FinanceLib.PeriodSeries([(0.5,0.05), (1.25, 0.0575), (2, 0.0485), (3.5, 0.0625), (4.25, 0.055)]), [-150, 20, 15, 80, 100], 0.3)

# @code_warntype FinanceLib.xnpv(FinanceLib.DateSeries([(Dates.Date(2014,9,20),0.05), (Dates.Date(2015,2,1), 0.0575), (Dates.Date(2016,10,5), 0.0485), (Dates.Date(2017,12,5), 0.0625), (Dates.Date(2019,1,5), 0.055)]), [-150, 20, 15, 80, 100], Dates.Date(2014,2,15))

# rC = FinanceLib.RateCurve{FinanceLib.NomRate}([0.05, 0.06, 0.07, 0.08], 2)

# @code_warntype FinanceLib.estimR(rC, 1.5)
# @code_warntype FinanceLib.estimR(rC, 1.2)

# eT = FinanceLib.RateCurve{FinanceLib.NomRate}([0.0016, 0.0021, 0.0027, 0.0033, 0.0037, 0.0041], 2)

# @code_warntype FinanceLib.nomR(FinanceLib.effR(FinanceLib.expR(eT))) # N - X - E - N
# @code_warntype FinanceLib.nomR(FinanceLib.expR(FinanceLib.effR(eT))) # N - E - X - N
# @code_warntype FinanceLib.nomR(FinanceLib.expR(eT)) # N - X - N

# println("Benchmarks")

# println("\nNPV")
# tS = collect(zip([0.25,6.25,3.5,4.5,1.25], [-6.25,1.2,1.25,3.6,2.5]))
# print("Separate Vec of tim and cf:"); @btime FinanceLib.npv(0.08, [0.25,6.25,3.5,4.5,1.25], [-6.25,1.2,1.25,3.6,2.5], 0.45)
# print("Tuple of (tim, cf):"); @btime FinanceLib.npv(0.08, tS, 0.45)

# print("Tuple of (tim, rate), cf:"); @btime FinanceLib.npv(FinanceLib.PeriodSeries([(0.5,0.05), (1.25, 0.0575), (2, 0.0485), (3.5, 0.0625), (4.25, 0.055)]), [-150, 20, 15, 80, 100], 0.3)

# println("\nXNPV")
# tS = collect(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28), Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)], [-15, 5, 25, -10, 50.]))

# print("Separate Vec of tim and cf:"); @btime FinanceLib.xnpv(0.08, [Dates.Date(2012,2,25), Dates.Date(2012,6,28), Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)], [-15, 5, 25, -10, 50], Dates.Date(2012,1,10) )

# print("Tuple of (tim, cf):"); @btime FinanceLib.xnpv(0.08, tS, Dates.Date(2012,1,10) )

# print("Tuple of (tim, rate), cf:"); @btime FinanceLib.xnpv(FinanceLib.DateSeries([(Dates.Date(2014,9,20),0.05), (Dates.Date(2015,2,1), 0.0575), (Dates.Date(2016,10,5), 0.0485), (Dates.Date(2017,12,5), 0.0625), (Dates.Date(2019,1,5), 0.055)]), [-150, 20, 15, 80, 100], Dates.Date(2014,2,15))

# println("\nIRR")
# tS = collect(zip([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25]))

# print("Separate Vec of tim and cf:"); @btime FinanceLib.irr([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25])
# print("Tuple of (tim, cf):"); @btime FinanceLib.irr(tS)

# println("\nXIRR")
# tS = collect(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28), Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],[-115, 5, 25, -10, 200]))

# print("Separate Vec of tim and cf:"); @btime FinanceLib.xirr([Dates.Date(2012,2,25), Dates.Date(2012,6,28), Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)], [-115, 5, 25, -10, 200] )
# print("Tuple of (tim, cf):"); @btime FinanceLib.xirr(tS)

# @btime FinanceLib.estimR(rC, 1.5)
# @btime FinanceLib.estimR(rC, 1.2)

# @btime FinanceLib.nomR(FinanceLib.effR(FinanceLib.expR(eT))) # N - X - E - N
# @btime FinanceLib.nomR(FinanceLib.expR(FinanceLib.effR(eT))) # N - E - X - N
# @btime FinanceLib.nomR(FinanceLib.expR(eT)) # N - X - N

# error("Hiya")

# include("FixedIncomes/perf.jl")
# include("Derivatives/perf.jl")
include("Statements/perf.jl")


text:
include("metric_provider.jl")
include("instruments.jl")


text:
# common facilities

# scalar functions

xlogx(x::Real) = x > zero(x) ? x * log(x) : zero(log(x))

xlogy{T<:Real}(x::T, y::T) = x > zero(T) ? x * log(y) : zero(log(x))
xlogy(x::Real, y::Real) = xlogy(promote(x, y)...)

# logistic: 1 / (1 + exp(-x))
#
logistic(x::Real) = one(x) / (one(x) + exp(-x))

# logit: log(x / (1 - x))
#
logit(x::Real) = log(x / (one(x) - x))

# log1psq: log(1+x^2)
#
log1psq(x::Real) = log1p(abs2(x))
log1psq(x::Union{Float32,Float64}) = (ax = abs(x); ax < maxintfloat(x) ? log1p(abs2(ax)) : 2 * log(ax))

# log1pexp: log(1+exp(x))
#
log1pexp(x::Real) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : oftype(exp(-x), x)
log1pexp(x::Float32) = x < 9.0f0 ? log1p(exp(x)) : x < 16.0f0 ? x + exp(-x) : oftype(exp(-x), x)

# log1mexp: log(1 - exp(x))
#
# See:
#   Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))"
#   http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
#
# Note: different than Maechler (2012), no negation inside parantheses
log1mexp(x::Real) = x < loghalf ? log1p(-exp(x)) : log(-expm1(x))

# log2mexp: log(2 - exp(x))
#
log2mexp(x::Real) = log1p(-expm1(x))

# logexpm1: log(exp(x) - 1)
#
logexpm1(x::Real) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x)
logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : oftype(exp(-x), x)

Compat.@dep_vectorize_1arg Real xlogx
Compat.@dep_vectorize_2arg Real xlogy
Compat.@dep_vectorize_1arg Real logistic
Compat.@dep_vectorize_1arg Real logit
Compat.@dep_vectorize_1arg Real log1psq
Compat.@dep_vectorize_1arg Real log1pexp
Compat.@dep_vectorize_1arg Real log1mexp
Compat.@dep_vectorize_1arg Real log2mexp
Compat.@dep_vectorize_1arg Real logexpm1

const softplus = log1pexp
const invsoftplus = logexpm1

# log1pmx: log(1 + x) - x
#
# use naive calculation or range reduction outside kernel range.
# accurate ~2ulps for all x
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return _log1pmx_ker(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return _log1pmx_ker(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return _log1pmx_ker(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    end
end

# logmxp1: log(x) - x + 1
#
function logmxp1(x::Float64)
    if x <= 0.3
        return (log(x) + 1.0) - x
    elseif x <= 0.4
        u = (x-0.375)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    elseif x <= 0.6
        u = 2.0*(x-0.5)
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        return log1pmx(x - 1.0)
    end
end

# The kernel of log1pmx
# Accuracy within ~2ulps for -0.227 < x < 0.315
function _log1pmx_ker(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @horner(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end

## logsumexp

logsumexp{T<:Real}(x::T, y::T) = x > y ? x + log1p(exp(y - x)) : y + log1p(exp(x - y))
logsumexp(x::Real, y::Real) = logsumexp(promote(x, y)...)

function logsumexp{T<:Real}(x::AbstractArray{T})
    isempty(x) && return -Inf
    u = maximum(x)
    s = 0.
    for i = 1:length(x)
        @inbounds s += exp(x[i] - u)
    end
    log(s) + u
end

## softmax

function softmax!{R<:AbstractFloat,T<:Real}(r::AbstractArray{R}, x::AbstractArray{T})
    n = length(x)
    length(r) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    u = maximum(x)
    s = 0.
    @inbounds for i = 1:n
        s += (r[i] = exp(x[i] - u))
    end
    invs = convert(R, inv(s))
    @inbounds for i = 1:n
        r[i] *= invs
    end
    r
end

softmax!{T<:AbstractFloat}(x::AbstractArray{T}) = softmax!(x, x)
softmax{T<:Real}(x::AbstractArray{T}) = softmax!(Array(Float64, size(x)), x)


text:
export Day8, Day8_1, Day8_2
abstract type Day8 <: Day end
abstract type Day8_1 <: Day8 end
abstract type Day8_2 <: Day8 end

abstract type ResultsType end
abstract type Acc <: ResultsType end
abstract type Jmp <: ResultsType end
abstract type Nop <: ResultsType end

struct Results{T<:ResultsType}
    line
    total
    type::Type{T}
end

line_parser(::Type{<:Day8}, line) = begin
    cmd, val = split(line)
    if startswith(val,"+")
        val = val[2:end]
    end
    val = parse(Int,val)
    if cmd == "nop"
        nop(val)
    elseif cmd == "jmp"
        jmp(val)
    elseif cmd == "acc"
        acc(val)
    end
end

process_program(lines) = begin
    total = 0
    processed = Set{Int}()
    line = 1
    while (line ∉ processed) & (line ≤ length(lines))
        push!(processed, line)
        res = lines[line]
        if res.type == Nop
            line += 1
        else
            line += res.line
        end
        total += res.total
    end
    total, line
end

solve(::Type{Day8_1}, lines::Vector) = process_program(lines)[1]

solve(::Type{Day8_2}, lines::Vector) = begin
    possible_changes = Dict{Int,Pair{Results,Results}}()

    for (i, l) in enumerate(lines)
        if l.type in Set([Nop, Jmp])
            possible_changes[i] = Pair(Results(l.line, l.total,Jmp), Results(l.line, l.total, Nop))
        end
    end
    for (k, values) in possible_changes
        org = lines[k]
        for v in values
            lines[k] = v
            t,l = process_program(lines)
            if l >length(lines)
                return t
            end
        end
        lines[k] = org
    end

end

nop(v) = Results(v,0, Nop)
acc(v) = Results(1,v, Acc)
jmp(v) = Results(v,0, Jmp)


text:
"""
Plots a desnity of a distribution for a selected metric (e.g. effective sample size)
* `df`: dataframe of results
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* `dir`: directory of saved plot. Default pwd.
"""
function plotdensity(df::DataFrame,metric::Symbol,group=(:sampler,);save=false,
    figfmt="png",dir="",options...)
    plots = Plots.Plot[]
    layout = SetLayout(df,group)
    grouping = map(x->df[!,x],group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        p=density(df[!,c],group=grouping,grid=false,xlabel=string(c),
            ylabel="Density",layout=layout,fill=(0,.5),width=1.5;options...)
        push!(plots,p)
        save ? savefig(p,string(dir,"density_",c,".",figfmt)) : nothing
    end
    return plots
end

"""
Plots a summary of a metric computed according to a function, func. func defaults to the mean.
* `df`: dataframe of results
* `xvar`: variable assigned to x-axis
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotsummary(df::DataFrame,xvar::Symbol,metric::Symbol,group=(:sampler,);save=false,
    figfmt="png",func=mean,dir="",options...)
    plots = Plots.Plot[]
    layout = SetLayout(df,group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        summary,yvar = summarize(df,c,[xvar,group...],func)
        dfse,yse = summarize(df,c,[xvar,group...],se)
        grouping = map(x->summary[!,x],group)
        p=plot(summary[!,xvar],summary[!,yvar],group=grouping,grid=false,xlabel=string(xvar),
            ylabel=string(yvar),layout=layout,width=1.5,yerror=dfse[!,yse];options...)
        push!(plots,p)
        save ? savefig(p,string(dir,"summary_",c,".",figfmt)) : nothing
    end
    return plots
end

"""
Generates a scatter plot to visualize the relationship between two
benchmarking metrics.
* `df`: dataframe of results
* `xvar`: variable assigned to x-axis
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotscatter(df::DataFrame,xvar::Symbol,metric::Symbol,group=(:sampler,);save=false,
    figfmt="png",func=mean,dir="",options...)
    plots = Plots.Plot[]
    layout = SetLayout(df,group)
    grouping = map(x->df[!,x],group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        p=scatter(df[!,xvar],df[!,c],group=grouping,grid=false,xlabel=string(xvar),
            ylabel=string(c),layout=layout,width=1.5;options...)
        push!(plots,p)
        save ? savefig(p,string(dir,"scatter_",xvar,"_",c,".",figfmt)) : nothing
    end
    return plots
end

"""
Generates a plot to compare the estimated parameters to true parameter values used to
generate simulated data.
* `df`: dataframe of results
* `parms`: namedtuple of parameter names and true values
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotrecovery(df::DataFrame,parms,group=(:sampler,);save=false,
    figfmt="png",dir="",options...)
    plots = Plots.Plot[]
    layout = SetLayout(df,group)
    for s in groupby(df,:sampler)
        sampler = s[!,:sampler][1]
        grouping = map(x->s[!,x],group)
        for (parm,v) in pairs(parms)
            μ = s[!,Symbol(string(parm,"_mean"))]
            lb = μ .- s[!,Symbol(string(parm,"_hdp_lb"))]
            ub = s[!,Symbol(string(parm,"_hdp_ub"))] .- μ
            p=scatter(μ,group=grouping,grid=false,ylabel=string(parm),layout=layout,
                leg=false,yerror=(lb,ub),width=1,title=string(sampler);options...)
            cnt = 0
            [hline!(p,[v],subplot=cnt+=1) for i in 1:layout[1] for j in 1:layout[2]]
            push!(plots,p)
            save ? savefig(p,string(dir,sampler,"_recovery_",parm,".",figfmt)) : nothing
        end
    end
    return plots
end

"""
Test whether column is a metric, e.g. mu_ess for ess.
"""
function isin(metric,col)
    occursin(string(metric),string(col))
end

"""
Creates a layout = (n,1) where n is the number of factors
in the last grouping variable.
"""
function SetLayout(df,group)
    isempty(group) ? (return (1,1)) : nothing
    length(group) == 1 ? (return (1,1)) : nothing
    col = group[end]
    n = length(unique(df[!,col]))
    return (n,1)
end

se(x)=std(x)/sqrt(length(x))

function summarize(df,metric,grouping,func)
    newDF= by(df,grouping,metric=>func)
    yvar = names(newDF)[end]
    return newDF,yvar
end


text:
# export JULIA_NUM_THREADS=‘4’

using Base.Threads: @threads, nthreads, threadid

# threadpool (generic function with 2 methods)

threadpool(fn_map, array) =
  begin
    array_return = Array[[] for _ in 1:nthreads()]
    @threads for e in array
        push!(array_return[threadid()], fn_map(e))
    end ; array_return = vcat(array_return...)

  array_return
  end

threadpool(fn_aggregate, fn_map, array) =
  begin
      array_return = Array[[] for _ in 1:nthreads()]
      @threads for e in array
          push!(array_return[threadid()], fn_map(e))
      end ; array_return = vcat(array_return...)
      value_return = fn_aggregate(array_return)

  value_return
  end

### ###

# length(procs()) < 2 ? addprocs(3) : ()

using Distributed: @everywhere, @distributed

using Distributed: procs, addprocs

# procpool (generic function with 2 methods)

procpool(fn_map, array) =

    @distributed (vcat) for e in array
        fn_map(e)
    end

procpool(fn_aggregate, fn_map, array) =

    fn_aggregate(
        begin @distributed (vcat) for e in array
                fn_map(e)
            end
        end
    )

### ###

batchify(list, batch_size) =

    [list[((i-1)*batch_size)+1:i*batch_size] for i in 1:convert(Int16,floor(length(list)/batch_size))]

##

indexof(e, list) =

    for (i,e2) in enumerate(list)

        e2 == e ? (return i) : ()

    end

indicesof(e, list) =

    findall(x->x==e,list)

##

len(arr) = length(arr)

resize(arr, size) = reshape(arr, size)

int(nr::AbstractFloat) = convert(Int32,floor(nr))
int(nr::AbstractString) = parse(Int32,nr)

float(nr::Integer) = convert(Float32,nr)
float(nr::AbstractString) = parse(Float32,nr)

str(nr) = "$(nr)"

##

using Random: shuffle, shuffle!

choice(arr) = arr[rand(1:length(arr))]

using Random: randperm

choices(arr,n;detailed=true) = detailed ? arr[randperm(length(arr))[1:n]] : shuffle(arr)[1:n]

enum(arr) = enumerate(arr)

type(e) = typeof(e)

##

input(prompt="") =
begin

    print(prompt)

chomp(readline())
end

##

@info "Utils read. \nthreads  : $(nthreads()) of $(length(Sys.cpu_info())) \nprocesses: $(length(procs())) of $(length(Sys.cpu_info()))"


text:
# Julia wrapper for header: /mnt/ai/home/dyuret/src/cudnn/4.0.4/include/cudnn.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function cudnnGetVersion()
    ccall((:cudnnGetVersion,libcudnn),Cint,())
end

function cudnnGetErrorString(status)
    ccall((:cudnnGetErrorString,libcudnn),Ptr{UInt8},(cudnnStatus_t,),status)
end

function cudnnCreate(handle)
    ccall((:cudnnCreate,libcudnn),cudnnStatus_t,(Ptr{cudnnHandle_t},),handle)
end

function cudnnDestroy(handle)
    ccall((:cudnnDestroy,libcudnn),cudnnStatus_t,(cudnnHandle_t,),handle)
end

function cudnnSetStream(handle,streamId)
    ccall((:cudnnSetStream,libcudnn),cudnnStatus_t,(cudnnHandle_t,Cint),handle,streamId)
end

function cudnnGetStream(handle,streamId)
    ccall((:cudnnGetStream,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Cint}),handle,streamId)
end

function cudnnCreateTensorDescriptor(tensorDesc)
    ccall((:cudnnCreateTensorDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnTensorDescriptor_t},),tensorDesc)
end

function cudnnSetTensor4dDescriptor(tensorDesc,format,dataType,n,c,h,w)
    ccall((:cudnnSetTensor4dDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,cudnnTensorFormat_t,cudnnDataType_t,Cint,Cint,Cint,Cint),tensorDesc,format,dataType,n,c,h,w)
end

function cudnnSetTensor4dDescriptorEx(tensorDesc,dataType,n,c,h,w,nStride,cStride,hStride,wStride)
    ccall((:cudnnSetTensor4dDescriptorEx,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,cudnnDataType_t,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),tensorDesc,dataType,n,c,h,w,nStride,cStride,hStride,wStride)
end

function cudnnGetTensor4dDescriptor(tensorDesc,dataType,n,c,h,w,nStride,cStride,hStride,wStride)
    ccall((:cudnnGetTensor4dDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,Ptr{cudnnDataType_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),tensorDesc,dataType,n,c,h,w,nStride,cStride,hStride,wStride)
end

function cudnnSetTensorNdDescriptor(tensorDesc,dataType,nbDims,dimA,strideA)
    ccall((:cudnnSetTensorNdDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,cudnnDataType_t,Cint,Ptr{Cint},Ptr{Cint}),tensorDesc,dataType,nbDims,dimA,strideA)
end

function cudnnGetTensorNdDescriptor(tensorDesc,nbDimsRequested,dataType,nbDims,dimA,strideA)
    ccall((:cudnnGetTensorNdDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,Cint,Ptr{cudnnDataType_t},Ptr{Cint},Ptr{Cint},Ptr{Cint}),tensorDesc,nbDimsRequested,dataType,nbDims,dimA,strideA)
end

function cudnnDestroyTensorDescriptor(tensorDesc)
    ccall((:cudnnDestroyTensorDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,),tensorDesc)
end

function cudnnTransformTensor(handle,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnTransformTensor,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnAddTensor(handle,alpha,bDesc,b,beta,yDesc,y)
    ccall((:cudnnAddTensor,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,bDesc,b,beta,yDesc,y)
end

function cudnnAddTensor_v3(handle,alpha,bDesc,b,beta,yDesc,y)
    ccall((:cudnnAddTensor_v3,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,bDesc,b,beta,yDesc,y)
end

function cudnnSetTensor(handle,yDesc,y,valuePtr)
    ccall((:cudnnSetTensor,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void}),handle,yDesc,y,valuePtr)
end

function cudnnScaleTensor(handle,yDesc,y,alpha)
    ccall((:cudnnScaleTensor,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void}),handle,yDesc,y,alpha)
end

function cudnnCreateFilterDescriptor(filterDesc)
    ccall((:cudnnCreateFilterDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnFilterDescriptor_t},),filterDesc)
end

function cudnnSetFilter4dDescriptor(filterDesc,dataType,k,c,h,w)
    ccall((:cudnnSetFilter4dDescriptor,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,cudnnDataType_t,Cint,Cint,Cint,Cint),filterDesc,dataType,k,c,h,w)
end

function cudnnSetFilter4dDescriptor_v4(filterDesc,dataType,format,k,c,h,w)
    ccall((:cudnnSetFilter4dDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,cudnnDataType_t,cudnnTensorFormat_t,Cint,Cint,Cint,Cint),filterDesc,dataType,format,k,c,h,w)
end

function cudnnGetFilter4dDescriptor(filterDesc,dataType,k,c,h,w)
    ccall((:cudnnGetFilter4dDescriptor,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,Ptr{cudnnDataType_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),filterDesc,dataType,k,c,h,w)
end

function cudnnGetFilter4dDescriptor_v4(filterDesc,dataType,format,k,c,h,w)
    ccall((:cudnnGetFilter4dDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,Ptr{cudnnDataType_t},Ptr{cudnnTensorFormat_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),filterDesc,dataType,format,k,c,h,w)
end

function cudnnSetFilterNdDescriptor(filterDesc,dataType,nbDims,filterDimA)
    ccall((:cudnnSetFilterNdDescriptor,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,cudnnDataType_t,Cint,Ptr{Cint}),filterDesc,dataType,nbDims,filterDimA)
end

function cudnnSetFilterNdDescriptor_v4(filterDesc,dataType,format,nbDims,filterDimA)
    ccall((:cudnnSetFilterNdDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,cudnnDataType_t,cudnnTensorFormat_t,Cint,Ptr{Cint}),filterDesc,dataType,format,nbDims,filterDimA)
end

function cudnnGetFilterNdDescriptor(filterDesc,nbDimsRequested,dataType,nbDims,filterDimA)
    ccall((:cudnnGetFilterNdDescriptor,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,Cint,Ptr{cudnnDataType_t},Ptr{Cint},Ptr{Cint}),filterDesc,nbDimsRequested,dataType,nbDims,filterDimA)
end

function cudnnGetFilterNdDescriptor_v4(filterDesc,nbDimsRequested,dataType,format,nbDims,filterDimA)
    ccall((:cudnnGetFilterNdDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,Cint,Ptr{cudnnDataType_t},Ptr{cudnnTensorFormat_t},Ptr{Cint},Ptr{Cint}),filterDesc,nbDimsRequested,dataType,format,nbDims,filterDimA)
end

function cudnnDestroyFilterDescriptor(filterDesc)
    ccall((:cudnnDestroyFilterDescriptor,libcudnn),cudnnStatus_t,(cudnnFilterDescriptor_t,),filterDesc)
end

function cudnnCreateConvolutionDescriptor(convDesc)
    ccall((:cudnnCreateConvolutionDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnConvolutionDescriptor_t},),convDesc)
end

function cudnnSetConvolution2dDescriptor(convDesc,pad_h,pad_w,u,v,upscalex,upscaley,mode)
    ccall((:cudnnSetConvolution2dDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Cint,Cint,Cint,Cint,Cint,cudnnConvolutionMode_t),convDesc,pad_h,pad_w,u,v,upscalex,upscaley,mode)
end

function cudnnGetConvolution2dDescriptor(convDesc,pad_h,pad_w,u,v,upscalex,upscaley,mode)
    ccall((:cudnnGetConvolution2dDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{cudnnConvolutionMode_t}),convDesc,pad_h,pad_w,u,v,upscalex,upscaley,mode)
end

function cudnnGetConvolution2dForwardOutputDim(convDesc,inputTensorDesc,filterDesc,n,c,h,w)
    ccall((:cudnnGetConvolution2dForwardOutputDim,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),convDesc,inputTensorDesc,filterDesc,n,c,h,w)
end

function cudnnSetConvolutionNdDescriptor(convDesc,arrayLength,padA,filterStrideA,upscaleA,mode,dataType)
    ccall((:cudnnSetConvolutionNdDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},cudnnConvolutionMode_t,cudnnDataType_t),convDesc,arrayLength,padA,filterStrideA,upscaleA,mode,dataType)
end

function cudnnGetConvolutionNdDescriptor(convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode,dataType)
    ccall((:cudnnGetConvolutionNdDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{cudnnConvolutionMode_t},Ptr{cudnnDataType_t}),convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode,dataType)
end

function cudnnSetConvolutionNdDescriptor_v3(convDesc,arrayLength,padA,filterStrideA,upscaleA,mode,dataType)
    ccall((:cudnnSetConvolutionNdDescriptor_v3,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},cudnnConvolutionMode_t,cudnnDataType_t),convDesc,arrayLength,padA,filterStrideA,upscaleA,mode,dataType)
end

function cudnnGetConvolutionNdDescriptor_v3(convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode,dataType)
    ccall((:cudnnGetConvolutionNdDescriptor_v3,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{cudnnConvolutionMode_t},Ptr{cudnnDataType_t}),convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode,dataType)
end

function cudnnGetConvolutionNdForwardOutputDim(convDesc,inputTensorDesc,filterDesc,nbDims,tensorOuputDimA)
    ccall((:cudnnGetConvolutionNdForwardOutputDim,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,Cint,Ptr{Cint}),convDesc,inputTensorDesc,filterDesc,nbDims,tensorOuputDimA)
end

function cudnnDestroyConvolutionDescriptor(convDesc)
    ccall((:cudnnDestroyConvolutionDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,),convDesc)
end

function cudnnFindConvolutionForwardAlgorithm(handle,xDesc,wDesc,convDesc,yDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
    ccall((:cudnnFindConvolutionForwardAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,Cint,Ptr{Cint},Ptr{cudnnConvolutionFwdAlgoPerf_t}),handle,xDesc,wDesc,convDesc,yDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
end

function cudnnGetConvolutionForwardAlgorithm(handle,xDesc,wDesc,convDesc,yDesc,preference,memoryLimitInBytes,algo)
    ccall((:cudnnGetConvolutionForwardAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionFwdPreference_t,Cint,Ptr{cudnnConvolutionFwdAlgo_t}),handle,xDesc,wDesc,convDesc,yDesc,preference,memoryLimitInBytes,algo)
end

function cudnnGetConvolutionForwardWorkspaceSize(handle,xDesc,wDesc,convDesc,yDesc,algo,sizeInBytes)
    ccall((:cudnnGetConvolutionForwardWorkspaceSize,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionFwdAlgo_t,Ptr{Cint}),handle,xDesc,wDesc,convDesc,yDesc,algo,sizeInBytes)
end

function cudnnConvolutionForward(handle,alpha,xDesc,x,wDesc,w,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,yDesc,y)
    ccall((:cudnnConvolutionForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionFwdAlgo_t,Ptr{Void},Cint,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,wDesc,w,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,yDesc,y)
end

function cudnnConvolutionBackwardBias(handle,alpha,dyDesc,dy,beta,dbDesc,db)
    ccall((:cudnnConvolutionBackwardBias,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,dyDesc,dy,beta,dbDesc,db)
end

function cudnnFindConvolutionBackwardFilterAlgorithm(handle,xDesc,dyDesc,convDesc,dwDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
    ccall((:cudnnFindConvolutionBackwardFilterAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnFilterDescriptor_t,Cint,Ptr{Cint},Ptr{cudnnConvolutionBwdFilterAlgoPerf_t}),handle,xDesc,dyDesc,convDesc,dwDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
end

function cudnnGetConvolutionBackwardFilterAlgorithm(handle,xDesc,dyDesc,convDesc,dwDesc,preference,memoryLimitInBytes,algo)
    ccall((:cudnnGetConvolutionBackwardFilterAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionBwdFilterPreference_t,Cint,Ptr{cudnnConvolutionBwdFilterAlgo_t}),handle,xDesc,dyDesc,convDesc,dwDesc,preference,memoryLimitInBytes,algo)
end

function cudnnGetConvolutionBackwardFilterWorkspaceSize(handle,xDesc,dyDesc,convDesc,gradDesc,algo,sizeInBytes)
    ccall((:cudnnGetConvolutionBackwardFilterWorkspaceSize,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionBwdFilterAlgo_t,Ptr{Cint}),handle,xDesc,dyDesc,convDesc,gradDesc,algo,sizeInBytes)
end

function cudnnConvolutionBackwardFilter(handle,alpha,xDesc,x,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dwDesc,dw)
    ccall((:cudnnConvolutionBackwardFilter,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionBwdFilterAlgo_t,Ptr{Void},Cint,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dwDesc,dw)
end

function cudnnConvolutionBackwardFilter_v3(handle,alpha,xDesc,x,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dwDesc,dw)
    ccall((:cudnnConvolutionBackwardFilter_v3,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionBwdFilterAlgo_t,Ptr{Void},Cint,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dwDesc,dw)
end

function cudnnFindConvolutionBackwardDataAlgorithm(handle,wDesc,dyDesc,convDesc,dxDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
    ccall((:cudnnFindConvolutionBackwardDataAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnFilterDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,Cint,Ptr{Cint},Ptr{cudnnConvolutionBwdDataAlgoPerf_t}),handle,wDesc,dyDesc,convDesc,dxDesc,requestedAlgoCount,returnedAlgoCount,perfResults)
end

function cudnnGetConvolutionBackwardDataAlgorithm(handle,wDesc,dyDesc,convDesc,dxDesc,preference,memoryLimitInBytes,algo)
    ccall((:cudnnGetConvolutionBackwardDataAlgorithm,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnFilterDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionBwdDataPreference_t,Cint,Ptr{cudnnConvolutionBwdDataAlgo_t}),handle,wDesc,dyDesc,convDesc,dxDesc,preference,memoryLimitInBytes,algo)
end

function cudnnGetConvolutionBackwardDataWorkspaceSize(handle,wDesc,dyDesc,convDesc,dxDesc,algo,sizeInBytes)
    ccall((:cudnnGetConvolutionBackwardDataWorkspaceSize,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnFilterDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionDescriptor_t,cudnnTensorDescriptor_t,cudnnConvolutionBwdDataAlgo_t,Ptr{Cint}),handle,wDesc,dyDesc,convDesc,dxDesc,algo,sizeInBytes)
end

function cudnnConvolutionBackwardData(handle,alpha,wDesc,w,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardData,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionBwdDataAlgo_t,Ptr{Void},Cint,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,wDesc,w,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dxDesc,dx)
end

function cudnnConvolutionBackwardData_v3(handle,alpha,wDesc,w,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardData_v3,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionBwdDataAlgo_t,Ptr{Void},Cint,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,wDesc,w,dyDesc,dy,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,dxDesc,dx)
end

function cudnnIm2Col(handle,xDesc,x,wDesc,convDesc,colBuffer)
    ccall((:cudnnIm2Col,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnTensorDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,cudnnConvolutionDescriptor_t,Ptr{Void}),handle,xDesc,x,wDesc,convDesc,colBuffer)
end

function cudnnSoftmaxForward(handle,algo,mode,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnSoftmaxForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnSoftmaxAlgorithm_t,cudnnSoftmaxMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,algo,mode,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnSoftmaxBackward(handle,algo,mode,alpha,yDesc,y,dyDesc,dy,beta,dxDesc,dx)
    ccall((:cudnnSoftmaxBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnSoftmaxAlgorithm_t,cudnnSoftmaxMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,algo,mode,alpha,yDesc,y,dyDesc,dy,beta,dxDesc,dx)
end

function cudnnCreatePoolingDescriptor(poolingDesc)
    ccall((:cudnnCreatePoolingDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnPoolingDescriptor_t},),poolingDesc)
end

function cudnnSetPooling2dDescriptor(poolingDesc,mode,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
    ccall((:cudnnSetPooling2dDescriptor,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnPoolingMode_t,Cint,Cint,Cint,Cint,Cint,Cint),poolingDesc,mode,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
end

function cudnnSetPooling2dDescriptor_v4(poolingDesc,mode,maxpoolingNanOpt,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
    ccall((:cudnnSetPooling2dDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnPoolingMode_t,cudnnNanPropagation_t,Cint,Cint,Cint,Cint,Cint,Cint),poolingDesc,mode,maxpoolingNanOpt,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
end

function cudnnGetPooling2dDescriptor(poolingDesc,mode,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
    ccall((:cudnnGetPooling2dDescriptor,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,Ptr{cudnnPoolingMode_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,mode,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
end

function cudnnGetPooling2dDescriptor_v4(poolingDesc,mode,maxpoolingNanOpt,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
    ccall((:cudnnGetPooling2dDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,Ptr{cudnnPoolingMode_t},Ptr{cudnnNanPropagation_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,mode,maxpoolingNanOpt,windowHeight,windowWidth,verticalPadding,horizontalPadding,verticalStride,horizontalStride)
end

function cudnnSetPoolingNdDescriptor(poolingDesc,mode,nbDims,windowDimA,paddingA,strideA)
    ccall((:cudnnSetPoolingNdDescriptor,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnPoolingMode_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,mode,nbDims,windowDimA,paddingA,strideA)
end

function cudnnSetPoolingNdDescriptor_v4(poolingDesc,mode,maxpoolingNanOpt,nbDims,windowDimA,paddingA,strideA)
    ccall((:cudnnSetPoolingNdDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnPoolingMode_t,cudnnNanPropagation_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,mode,maxpoolingNanOpt,nbDims,windowDimA,paddingA,strideA)
end

function cudnnGetPoolingNdDescriptor(poolingDesc,nbDimsRequested,mode,nbDims,windowDimA,paddingA,strideA)
    ccall((:cudnnGetPoolingNdDescriptor,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,Cint,Ptr{cudnnPoolingMode_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,nbDimsRequested,mode,nbDims,windowDimA,paddingA,strideA)
end

function cudnnGetPoolingNdDescriptor_v4(poolingDesc,nbDimsRequested,mode,maxpoolingNanOpt,nbDims,windowDimA,paddingA,strideA)
    ccall((:cudnnGetPoolingNdDescriptor_v4,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,Cint,Ptr{cudnnPoolingMode_t},Ptr{cudnnNanPropagation_t},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,nbDimsRequested,mode,maxpoolingNanOpt,nbDims,windowDimA,paddingA,strideA)
end

function cudnnGetPoolingNdForwardOutputDim(poolingDesc,inputTensorDesc,nbDims,outputTensorDimA)
    ccall((:cudnnGetPoolingNdForwardOutputDim,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnTensorDescriptor_t,Cint,Ptr{Cint}),poolingDesc,inputTensorDesc,nbDims,outputTensorDimA)
end

function cudnnGetPooling2dForwardOutputDim(poolingDesc,inputTensorDesc,n,c,h,w)
    ccall((:cudnnGetPooling2dForwardOutputDim,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,cudnnTensorDescriptor_t,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),poolingDesc,inputTensorDesc,n,c,h,w)
end

function cudnnDestroyPoolingDescriptor(poolingDesc)
    ccall((:cudnnDestroyPoolingDescriptor,libcudnn),cudnnStatus_t,(cudnnPoolingDescriptor_t,),poolingDesc)
end

function cudnnPoolingForward(handle,poolingDesc,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnPoolingForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnPoolingDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,poolingDesc,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnPoolingBackward(handle,poolingDesc,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
    ccall((:cudnnPoolingBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnPoolingDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,poolingDesc,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
end

function cudnnCreateActivationDescriptor(activationDesc)
    ccall((:cudnnCreateActivationDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnActivationDescriptor_t},),activationDesc)
end

function cudnnSetActivationDescriptor(activationDesc,mode,reluNanOpt,reluCeiling)
    ccall((:cudnnSetActivationDescriptor,libcudnn),cudnnStatus_t,(cudnnActivationDescriptor_t,cudnnActivationMode_t,cudnnNanPropagation_t,Cdouble),activationDesc,mode,reluNanOpt,reluCeiling)
end

function cudnnGetActivationDescriptor(activationDesc,mode,reluNanOpt,reluCeiling)
    ccall((:cudnnGetActivationDescriptor,libcudnn),cudnnStatus_t,(cudnnActivationDescriptor_t,Ptr{cudnnActivationMode_t},Ptr{cudnnNanPropagation_t},Ptr{Cdouble}),activationDesc,mode,reluNanOpt,reluCeiling)
end

function cudnnDestroyActivationDescriptor(activationDesc)
    ccall((:cudnnDestroyActivationDescriptor,libcudnn),cudnnStatus_t,(cudnnActivationDescriptor_t,),activationDesc)
end

function cudnnActivationForward(handle,mode,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnActivationForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnActivationMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,mode,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnActivationForward_v4(handle,activationDesc,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnActivationForward_v4,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnActivationDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,activationDesc,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnActivationBackward(handle,mode,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
    ccall((:cudnnActivationBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnActivationMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,mode,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
end

function cudnnActivationBackward_v4(handle,activationDesc,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
    ccall((:cudnnActivationBackward_v4,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnActivationDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,activationDesc,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
end

function cudnnCreateLRNDescriptor(normDesc)
    ccall((:cudnnCreateLRNDescriptor,libcudnn),cudnnStatus_t,(Ptr{cudnnLRNDescriptor_t},),normDesc)
end

function cudnnSetLRNDescriptor(normDesc,lrnN,lrnAlpha,lrnBeta,lrnK)
    ccall((:cudnnSetLRNDescriptor,libcudnn),cudnnStatus_t,(cudnnLRNDescriptor_t,UInt32,Cdouble,Cdouble,Cdouble),normDesc,lrnN,lrnAlpha,lrnBeta,lrnK)
end

function cudnnGetLRNDescriptor(normDesc,lrnN,lrnAlpha,lrnBeta,lrnK)
    ccall((:cudnnGetLRNDescriptor,libcudnn),cudnnStatus_t,(cudnnLRNDescriptor_t,Ptr{UInt32},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),normDesc,lrnN,lrnAlpha,lrnBeta,lrnK)
end

function cudnnDestroyLRNDescriptor(lrnDesc)
    ccall((:cudnnDestroyLRNDescriptor,libcudnn),cudnnStatus_t,(cudnnLRNDescriptor_t,),lrnDesc)
end

function cudnnLRNCrossChannelForward(handle,normDesc,lrnMode,alpha,xDesc,x,beta,yDesc,y)
    ccall((:cudnnLRNCrossChannelForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnLRNDescriptor_t,cudnnLRNMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,normDesc,lrnMode,alpha,xDesc,x,beta,yDesc,y)
end

function cudnnLRNCrossChannelBackward(handle,normDesc,lrnMode,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
    ccall((:cudnnLRNCrossChannelBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnLRNDescriptor_t,cudnnLRNMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,normDesc,lrnMode,alpha,yDesc,y,dyDesc,dy,xDesc,x,beta,dxDesc,dx)
end

function cudnnDivisiveNormalizationForward(handle,normDesc,mode,alpha,xDesc,x,means,temp,temp2,beta,yDesc,y)
    ccall((:cudnnDivisiveNormalizationForward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnLRNDescriptor_t,cudnnDivNormMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,normDesc,mode,alpha,xDesc,x,means,temp,temp2,beta,yDesc,y)
end

function cudnnDivisiveNormalizationBackward(handle,normDesc,mode,alpha,xDesc,x,means,dy,temp,temp2,beta,dXdMeansDesc,dx,dMeans)
    ccall((:cudnnDivisiveNormalizationBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnLRNDescriptor_t,cudnnDivNormMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void}),handle,normDesc,mode,alpha,xDesc,x,means,dy,temp,temp2,beta,dXdMeansDesc,dx,dMeans)
end

function cudnnDeriveBNTensorDescriptor(derivedBnDesc,xDesc,mode)
    ccall((:cudnnDeriveBNTensorDescriptor,libcudnn),cudnnStatus_t,(cudnnTensorDescriptor_t,cudnnTensorDescriptor_t,cudnnBatchNormMode_t),derivedBnDesc,xDesc,mode)
end

function cudnnBatchNormalizationForwardTraining(handle,mode,alpha,beta,xDesc,x,yDesc,y,bnScaleBiasMeanVarDesc,bnScale,bnBias,exponentialAverageFactor,resultRunningMean,resultRunningInvVariance,epsilon,resultSaveMean,resultSaveInvVariance)
    ccall((:cudnnBatchNormalizationForwardTraining,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnBatchNormMode_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},Cdouble,Ptr{Void},Ptr{Void},Cdouble,Ptr{Void},Ptr{Void}),handle,mode,alpha,beta,xDesc,x,yDesc,y,bnScaleBiasMeanVarDesc,bnScale,bnBias,exponentialAverageFactor,resultRunningMean,resultRunningInvVariance,epsilon,resultSaveMean,resultSaveInvVariance)
end

function cudnnBatchNormalizationForwardInference(handle,mode,alpha,beta,xDesc,x,yDesc,y,bnScaleBiasMeanVarDesc,bnScale,bnBias,estimatedMean,estimatedInvVariance,epsilon)
    ccall((:cudnnBatchNormalizationForwardInference,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnBatchNormMode_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Cdouble),handle,mode,alpha,beta,xDesc,x,yDesc,y,bnScaleBiasMeanVarDesc,bnScale,bnBias,estimatedMean,estimatedInvVariance,epsilon)
end

function cudnnBatchNormalizationBackward(handle,mode,alpha,beta,xDesc,x,dyDesc,dy,dxDesc,dx,dBnScaleBiasDesc,bnScale,dBnScaleResult,dBnBiasResult,epsilon,savedMean,savedInvVariance)
    ccall((:cudnnBatchNormalizationBackward,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnBatchNormMode_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},Ptr{Void},Cdouble,Ptr{Void},Ptr{Void}),handle,mode,alpha,beta,xDesc,x,dyDesc,dy,dxDesc,dx,dBnScaleBiasDesc,bnScale,dBnScaleResult,dBnBiasResult,epsilon,savedMean,savedInvVariance)
end

function cudnnSetConvolutionNdDescriptor_v2(convDesc,arrayLength,padA,filterStrideA,upscaleA,mode)
    ccall((:cudnnSetConvolutionNdDescriptor_v2,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},cudnnConvolutionMode_t),convDesc,arrayLength,padA,filterStrideA,upscaleA,mode)
end

function cudnnGetConvolutionNdDescriptor_v2(convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode)
    ccall((:cudnnGetConvolutionNdDescriptor_v2,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{cudnnConvolutionMode_t}),convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode)
end

function cudnnAddTensor_v2(handle,mode,alpha,bDesc,b,beta,yDesc,y)
    ccall((:cudnnAddTensor_v2,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnAddMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,mode,alpha,bDesc,b,beta,yDesc,y)
end

function cudnnConvolutionBackwardFilter_v2(handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardFilter_v2,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
end

function cudnnConvolutionBackwardData_v2(handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardData_v2,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
end

### To make things work on v2:

function cudnnSetConvolutionNdDescriptor(convDesc,arrayLength,padA,filterStrideA,upscaleA,mode)
    ccall((:cudnnSetConvolutionNdDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},cudnnConvolutionMode_t),convDesc,arrayLength,padA,filterStrideA,upscaleA,mode)
end

function cudnnGetConvolutionNdDescriptor(convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode)
    ccall((:cudnnGetConvolutionNdDescriptor,libcudnn),cudnnStatus_t,(cudnnConvolutionDescriptor_t,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{cudnnConvolutionMode_t}),convDesc,arrayLengthRequested,arrayLength,padA,strideA,upscaleA,mode)
end

function cudnnAddTensor(handle,mode,alpha,bDesc,b,beta,yDesc,y)
    ccall((:cudnnAddTensor,libcudnn),cudnnStatus_t,(cudnnHandle_t,cudnnAddMode_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,mode,alpha,bDesc,b,beta,yDesc,y)
end

function cudnnConvolutionBackwardFilter(handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardFilter,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
end

function cudnnConvolutionBackwardData(handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
    ccall((:cudnnConvolutionBackwardData,libcudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,xDesc,x,dyDesc,dy,convDesc,beta,dxDesc,dx)
end


text:
#=

  Constraint utils in Julia ConstraintSolver.jl

  Many of these decompositions were inspired/ported from
  my Picat implementations (http://hakank.org/picat/)

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

#
# resize_matrix(grid)
#
# Resize a (column) matrix of the form (3,)
#    [[1, 2,3],
#     [4,5,6],
#     [7,8,9]
#     ]
# To a 3x3 (row) matrix
#    [ 1 2 3;
#      4 5 6;
#      7 8 9
#    ]
#
function resize_matrix(grid)
    rows = size(grid)
    # reduce(hcat,grid)'
    transpose(reduce(hcat,grid)) # Suggestion by @Mason Protter @ Zulip
end

#
# Print a grid (row matrix)
#
function print_grid(grid)
    rows, cols = size(grid)
    for row in 1:rows
        println(grid[row,:])
    end
    println()
end

#
# print_matrix_str(x,zstring)
#
# Print each row and convert each element to the corresponding position
# in zstring. "z" is for zero-based, i.e. first char corresponds to 0, etc
#
function print_matrix_str(x,zstring="")
    for row in eachrow(x)
        if zstring == ""
            println(row)
        else
            println(join(row.|>c->zstring[c+1]))
        end
    end
    println()
end

#
# Scalar product(model, s,x,v)
#
# Ensure that
#    s == x.*v
# Where x or v might be an array of decision variables
#
function scalar_product(model,x,v,s)
    @constraint(model, s == sum(x.*v))
end

#
# scalar_product(model,x,v,cond,x)
#
# Ensures that x.*v op s
#
function scalar_product(model,x,v,cond,s)
    if cond == :(>=)
        @constraint(model, sum(x.*v) >= s)
    elseif cond == :(<=)
        @constraint(model, sum(x.*v) <= s)
    elseif cond == :(==)
        @constraint(model, sum(x.*v) == s)
    else
        error("scalar_product(model,s,x,cond,v): unknown cond: $cond")
    end

end

#
# to_num(model, x, base, num)
#
# Converts a number num to/from a list of integer list (x) given a base base
#
function to_num(model, x, base, num)
   len = length(x)
   @constraint(model, num == sum([x[i]*base^(len-i) for i in 1:len]))
end

#
# increasing(model, x)
#
# Ensure that array x in increasing order
#
function increasing(model, x)
    len = length(x)
    for i in 2:len
        @constraint(model, x[i-1] <= x[i])
    end
end

function increasing_strict(model, x)
    len = length(x)
    for i in 2:len
        @constraint(model, x[i-1] <= x[i])
        @constraint(model, x[i-1] != x[i])
    end
end

#
# decreasing(model, x)
#
# Ensure that array x in decreasing order
#
function decreasing(model, x)
    len = length(x)
    for i in 2:len
        @constraint(model, x[i-1] >= x[i])
    end
end

function decreasing_strict(model, x)
    len = length(x)
    for i in 2:len
        @constraint(model, x[i-1] >= x[i])
        @constraint(model, x[i-1] != x[i])
    end
end

#
# all_different_except_c
#
# Ensure that all values (except c) are distinct
# Thanks to Ole who fixed some initial problems I had.
# (See https://github.com/Wikunia/ConstraintSolver.jl/issues/202 for
# details.)
#
function all_different_except_c(model, x, c=0)
    n = length(x)

    # Define the variables we'll use
    b_len = length([1 for i in 2:n for j in 1:i-1 for k in 1:3])
    bs = @variable(model, [1:b_len], Bin) # "Anonymous" variables
    c = 1
    for i in 2:n, j in 1:i-1
        b1 = bs[c]
        b2 = bs[c+1]
        b3 = bs[c+2]
        @constraint(model, b1 := {x[i] != 0})
        @constraint(model, b2 := {x[j] != 0})
        @constraint(model, b3 := {b1 + b2 == 2})
        @constraint(model, b3 => {x[i] != x[j]})
        c += 3
    end
    # return bs so we can print it in the main function
    return bs
end

#
# count_ctr(model, x, val, op, s)
#
# Ensure that there the exactly s occurrences of
#    x[i] op val
#
# TODO: handle < and >
#       i.e. <= && !=  and >= && !=
#
function count_ctr(model, x, op, val, s)
    xflat = vcat(x...)
    len = length(xflat)
    b = @variable(model, [1:len], Bin)
    for i in 1:len
        if op == :(==)
            @constraint(model, b[i] := { xflat[i] == val})
        elseif op == :(<=)
            @constraint(model, b[i] := { xflat[i] <= val})
        elseif op == :(>=)
            @constraint(model, b[i] := { xflat[i] >= val})
        elseif op == :(!=)
            @constraint(model, b[i] := { xflat[i] != val})
        end
    end
    @constraint(model, s == sum(b))
end

#
# my_abs(model, x, y, d)
#
#  d = abs(x-y)
#
function my_abs(model, x, y, d)
    b = @variable(model, [1:1], Bin)
    @constraint(model, b[1] := {x >= y})
    @constraint(model, b[1] => {d == x - y})
    @constraint(model, !b[1] => {d == y - x})
end

# my_abs(model, x, a)
#
#  a = abs(x)
#
function my_abs(model, x, a)
    my_abs(model,x, 0, a)
    #=
    b = @variable(model, [1:1], Bin)
    @constraint(model, b[1] := {x >= 0})
    @constraint(model, b[1] => {a == x})
    @constraint(model, !b[1] => {a == -x})
    =#
end

#
# my_element(model, a, ix, val)
#
# Ensure that a[ix] = val
#
function my_element(model, ix, a, val)
    n = length(a)
    b = @variable(model, [1:n], Bin)
    for i in 1:n
        @constraint(model,b[i] => {a[i] == val})
        @constraint(model,b[i] := {ix == i})
    end
end

#
# Both a and gcc are lists.
#
# This version is bidirectional but limited:
#
# The list a can contain only values 1..gcc_len (i.e. the length of gcc).
# This means that the caller must know the max values of a.
# Or rather: if a contains another values they will not be counted.
#
function global_cardinality_count(model, a, gcc)
    len = length(a)
    gcc_len = length(gcc)
    for i in 1:gcc_len
       # count_ctr(model, x, op, val, s)
       count_ctr(model, a,:(==), i, gcc[i])
    end
    # This works as well
    # (1:length(gcc)).|>i->count_ctr(model, a,:(==),i,gcc[i])
end

#
# modulo(model, x, y, z)
#
# Ensures that  x mod y == z.
#
# Note: experimental.
# For large domains the table will blow up so handle with care!
#
function modulo(model, x, y, z)
    lbx = round(Int, JuMP.lower_bound(x))
    ubx = round(Int, JuMP.upper_bound(x))
    lby = round(Int, JuMP.lower_bound(y))
    uby = round(Int, JuMP.upper_bound(y))

    table = resize_matrix([ [i,j,i % j] for i in lbx:ubx, j in lby:uby if j != 0])
    @constraint(model, [x, y, z] in CS.TableSet(table))
end

#
# my_min(model, x, m)
#
# Ensures that m is the minimum value in x
#
function my_min(model, x, m)
    n = length(x)
    for i in 1:n
        @constraint(model, m <= x[i] )
    end
end

#
# my_max(model, x, m)
#
# Ensures that m is the maximum value in x
#
function my_max(model, x, m)
    n = length(x)
    for i in 1:n
        @constraint(model, m >= x[i] )
    end
end

#
# my_argmin(model, x, m, ix)
#
# Ensure that m is the minimum element and m = x[ix]
#
function my_argmin(model, x, m, ix)
    n = length(x)
    my_element(model, ix,x, m)
    for i in 1:n
        @constraint(model, m <= x[i] )
    end
end

#
# my_argmax(model, x, m, ix)
#
# Ensure that m is the maximum element and m = x[ix]
#
function my_argmax(model, x, m, ix)
    n = length(x)
    my_element(model, ix,x, m)
    for i in 1:n
        @constraint(model, m >= x[i] )
    end
end

#
# Decompositon of cumulative.
#
# Inspired by the MiniZinc implementation:
# http://www.g12.csse.unimelb.edu.au/wiki/doku.php?id=g12:zinc:lib:minizinc:std:cumulative.mzn&s[]=cumulative
# The MiniZinc decomposition is discussed in the paper:
# A. Schutt, T. Feydy, P.J. Stuckey, and M. G. Wallace.
# 'Why cumulative decomposition is not as bad as it sounds.'
# Download:
# http://www.cs.mu.oz.au/%7Epjs/rcpsp/papers/cp09-cu.pdf
# http://www.cs.mu.oz.au/%7Epjs/rcpsp/cumu_lazyfd.pdf
#
#
# Parameters:
#
# s: start_times    assumption: array of IntVar
# d: durations      assumption: array of int
# r: resources      assumption: array of int
# b: resource limit assumption: IntVar or int
#
function cumulative(model, start, duration, resource, limit)
    tasks = [i for i in 1:length(start) if resource[i] > 0 && duration[i] > 0]
    num_tasks = length(tasks)

    times_min = minimum(round.(Int,[JuMP.lower_bound(start[i]) for i in tasks]))
    times_max = maximum(round.(Int,[JuMP.upper_bound(start[i])+duration[i] for i in tasks]))
    println("times: $(times_min)..$(times_max)")
    for t in times_min:times_max
        bs = @variable(model, [1:num_tasks], Bin)
        bt = @variable(model, [1:num_tasks], Bin)
        b  = @variable(model, [1:num_tasks], Bin)
        for i in tasks
            # The following don't work since ConstraintSolver don't
            # support nonlinear constraints
            # @constraint(model,sum([(start[i] <= t) * (t <= start[i] + duration[i])*resource[i] for i in tasks])  <= b)

            # is this task active during this time t?
            @constraint(model, bs[i] := {start[i] <= t})
            @constraint(model, bt[i] := {t <= start[i]+duration[i]-1}) # should be '<'
            @constraint(model, b[i] := { bs[i] + bt[i] == 2}) # is this task active in time t ?
        end
        # Check that there's no conflicts in time t
        @constraint(model,sum([b[i]*resource[i] for i in tasks]) <= limit)
  end

end

#
# circuit(x) succeeds for the array x if it's a circuit.
#
# This implementation use an extra array (z) for the orbit of x[1].
#
function circuit(model, x)
    n = length(x)
    z = @variable(model, 1 <= z[1:n] <= n, Int)

    #
    # The main constraint is that Z[I] must not be 1
    # until I = N, and for I = N it must be 1.
    #
    @constraint(model, x in CS.AllDifferentSet())
    @constraint(model, z in CS.AllDifferentSet())

    # put the orbit of x[1] in z[1..n]
    @constraint(model, x[1] == z[1])

    # when i = n it must be 1
    @constraint(model, z[n] == 1)

    # Get the orbit for z
    for i in 2:n
       my_element(model, z[i-1],x,z[i])
    end

end

#
# circuit_path(model,x,path) succeeds for the array x if it's a circuit
# and path is the corresponding path (orbit).
#
# This implementation use an extra array (z) for the orbit of x[1].
#
function circuit_path(model, x, path)
    n = length(x)

    #
    # The main constraint is that Z[I] must not be 1
    # until I = N, and for I = N it must be 1.
    #
    # @constraint(model, x in CS.AllDifferentSet())
    # @constraint(model, path in CS.AllDifferentSet())

    # put the orbit of x[1] in z[1..n]
    @constraint(model, x[1] == path[1])

    # when i = n it must be 1
    @constraint(model, path[n] == 1)

    # Get the orbit for path
    for i in 2:n
       my_element(model, path[i-1],x,path[i])
    end

end

#
# inverse(model, x)
#
#  From Global Constraint Catalog
#  http://www.emn.fr/z-info/sdemasse/gccat/Cinverse.html
#  """
#  inverse(NODES)
#
# Enforce each vertex of a digraph to have exactly one predecessor and
# one successor. In addition the following two statements are equivalent:
#    - The successor of the ith node is the jth node.
#    - The predecessor of the jth node is the ith node.
# """
#
#  This means that for each element X[I] either
#     - X[I] = I
#     or
#     - X[I] = J <=> X[J] = I
#
# Cf the 2d version assignment_ctr(model, x, y) below
#
function inverse(model, x)
    len = length(x)
    b = @variable(model, [1:len,1:len], Bin)
    # Picat:
    # foreach(I in 1..Len, J in 1..Len)
    #  J = X[I] <=> I = X[J]
    # end
    for i in 1:len, j in 1:len
       @constraint(model, b[i,j] := { x[i] == j})
       @constraint(model, b[j,i] := { x[j] == i})
       @constraint(model, (b[i,j] + b[j,i]) != 1 )
    end
end

#=
From the Picat documentation
"""
assignment(FDVars1,FDVars2): This constraint ensures that FDVars2 is a dual
assignment of FDVars1, i.e., if the ith element of FDVars1 is j,
then the jth element of FDVars2 is i.
"""

The constraint can be defined as:
  foreach(I in 1..N, J in 1..N)
     X[I] = J <=> Y[J] = I
  end.

Note: We assume that both x and y are distinct arrays (respectively).
Also: Compare with the single array constraint inverse (tested in inverse.jl)

=#
function assignment_ctr(model, x, y)
    n = length(x)
    b1 = @variable(model, [1:n,1:n], Bin)
    b2 = @variable(model, [1:n,1:n], Bin)
    for i in 1:n, j in 1:n
       @constraint(model, b1[i,j] := { x[i] == j})
       @constraint(model, b2[j,i] := { y[j] == i})
       @constraint(model, b1[i,j] + b2[j,i] != 1 )
    end
end

#=
   matrix_element(model, x, i ,j, val)

   Ensures that x[i,j] = val

=#
function matrix_element(model, x,i,j,val)
    n, m = size(x)
    bi = @variable(model,[1:n], Bin) # row
    bj = @variable(model,[1:m], Bin) # column
    bij = @variable(model, [1:n,1:m], Bin) # matrix
    @constraint(model, sum(bi) == 1)
    @constraint(model, sum(bj) == 1)
    @constraint(model, sum(bij[:]) == 1)
    for k in 1:n, l in 1:m
        @constraint(model,bi[k] := { k == i }) # fix row
        @constraint(model,bj[l] := { l == j }) # fix column
        @constraint(model,bij[k,l] := { bi[k] + bj[l] == 2 })
        # Connection constraints
        # These don't work now (ConstraintSolver.jl v0.5.3)
        # since for  b => {....} b must be a boolen variable
        # See https://github.com/Wikunia/ConstraintSolver.jl/issues/224
        # @constraint(model, val == x[k,l]  => {bi[k] + bj[l] == 2} )
        # @constraint(model, (bi[k] + bj[l] == 2) => { val == x[k,l] } )

        @constraint(model, bij[k,l] => { val == x[k,l] } )
    end
end

#=
  regular(model,x,q,s,q0,f)

  This is a translation of MiniZinc's regular constraint (defined in
  lib/zinc/globals.mzn), via the Comet code (via Picat).
  All comments are from the MiniZinc code.
  """
  The sequence of values in array 'x' (which must all be in the range 1..S)
  is accepted by the DFA of 'Q' states with input 1..S and transition
  function 'd' (which maps (1..Q, 1..S) -> 0..Q)) and initial state 'q0'
  (which must be in 1..Q) and accepting states 'F' (which all must be in
  1..Q).  We reserve state 0 to be an always failing state.
  """

  x : IntVar array
  Q : number of states
  S : input_max
  d : transition matrix
  q0: initial state
  F : accepting states
=#
function regular(model,x,q,s,d,q0,f)
    # """
    # If x has index set m..n-1, then a[m] holds the initial state
    # (q0), and a[i+1] holds the state we're in after  processing
    # x[i].  If a[n] is in F, then we succeed (ie. accept the string).
    # """
    m = 1
    n2 = length(x)+1
    # The selected states (which is returned)
    a = @variable(model, 1 <= a[1:n2] <= q, Int)

    @constraint(model, a[m] == q0) # Set a[0], initial state
    for i in 1:length(x)
        @constraint(model, 1 <= x[i])      # Do this in case it's a var.
        @constraint(model, x[i] <= s)      # Do this in case it's a var.
        # Here is MiniZinc's infamous matrix element
        #     a[i+1] = d[a[i], x[i]]
        matrix_element(model, d, a[i], x[i], a[i+1])
    end
    # """Check the final state is in F."""
    f_ix = @variable(model, [1:1],CS.Integers(1:length(f))) # index in accepting states
    my_element(model,f_ix[1],f,a[n2])

    return a
end

#
# atmost(model,x,v,m)
#
# Ensures that the number of elements with value v is <= m
#
function atmost(model,x,v,m)
    n = length(x)
    b = @variable(model, [1:n],Bin)
    for i in 1:n
        @constraint(model,b[i]:={x[i] == v})
    end
    @constraint(model, sum(b) <= m)
end

#
# atleast(model,x,v,m)
#
# Ensures that the number of elements with value v is >= m
#
function atleast(model,x,v,m)
    n = length(x)
    b = @variable(model, [1:n],Bin)
    for i in 1:n
        @constraint(model,b[i]:={x[i] == v})
    end
    @constraint(model, sum(b) >= m)
end

#
# exactly(model,x,v,m)
#
# Ensures that the number of elements with value v is == m
#
function exactly(model,x,v,m)
    n = length(x)
    b = @variable(model, [1:n],Bin)
    for i in 1:n
        @constraint(model,b[i]:={x[i] == v})
    end
    @constraint(model, sum(b) == m)
end

#
# latin_square(model, x)
#
# Ensures that the square matrix is a Latin Square.
#
function latin_square(model, x)
    n,_ = size(x)
    for i in 1:n
        @constraint(model, x[i,:] in CS.AllDifferentSet())
        @constraint(model, x[:,i] in CS.AllDifferentSet())
    end
end

#
# no_overlap(model, begins,durations)
#
# Ensure that there is no overlap between the tasks.
#
function no_overlap(model, begins,durations)
    n = length(begins)
    for i in 1:n, j in i+1:n
        b = @variable(model,[1:2], Bin)
        @constraint(model,b[1] := {begins[i] + durations[i] <= begins[j]})
        @constraint(model,b[2] := {begins[j] + durations[j] <= begins[i]})
        @constraint(model, sum(b) >= 1)
    end
end

#
# mult_table2(model,lb=1,ub=9) # 2D
#
# Generate a multiplication table of i*j and i*j*k
# for i,j,.. in lb:ub
#
# i*j
function mult_table2(lb=1,ub=9)
    resize_matrix([[i,j,i*j] for i in lb:ub, j in lb:ub if i != j])
end

# i*j*k
function mult_table3(lb=1,ub=9)
    resize_matrix([[i,j,k,i*j*k] for i in lb:ub, j in lb:ub, k in lb:ub if i != j && i != k && j != k])
end

# i*j*k*l
function mult_table4(lb=1,ub=9)
    resize_matrix([[i,j,k,l,i*j*k*l] for i in lb:ub, j in lb:ub, k in lb:ub, l in lb:ub
                                     if i != j && i != k && i != l &&
                                        j != k && j != l &&
                                        k != l
                                        ])
end

#
# global_contiguity_regular(model, x)
#
# Ensures that x contains 0 and 1s and that the
# 1s are in a contiguous sequence.
#
function global_contiguity_regular(model, x)
    n = length(x)

    # Transition function (MiniZinc style)
    # This use the regular expression "0*1*0*" to
    # require that all 1's (if any) in an array appear contiguously.
    transition = resize_matrix([
                  [1,2], # state 1 (start) input 0 -> state 1, input 1 -> state 2 i.e. 0*
                  [3,2], # state 2: 1*
                  [3,0]  # state 3: 0*
                 ])
    n_states = 3
    input_max = 2
    initial_state = 1
    accepting_states = [1,2,3]

    reg_input = @variable(model, [1:n], CS.Integers(1:input_max)) # 1..2

    # Translate x's 0..1 to reg_input's 1..2
    # (since regular use 0 as an invalid state)
    for i in 1:n
       @constraint(model,reg_input[i] == x[i]+1)
    end

    return regular(model,reg_input,n_states,input_max,transition,initial_state, accepting_states)

end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module BuildOverlap

export buildelecoverlap

using TypesParticles, TypesBasis
using  UtilityFunctions

@doc raw"""
                      buildelecoverlap(basis::Basis)

returns the building of electron overlap integral matrix of expectation values:

     S_{m,n} = langle m | n rangle = int d tau phi_m phi_n

Description: This function calculates the electron-electron overlap matrix which essential captures
the inner product projection of basis function m on n, this is neccessary to track basis-functions that
are not orthogonal to one another.

NOTES: multiple dispatch only for GaussOrbitals

"""
function buildelecoverlap(natoms::Int,basisfunc::Array{GaussOrbitals},basis::Basis)
    numbasisfunc = basis.nbasisfunc; #number of basis functions
    osize = natoms*basis.nbasisfunc; #overlap size
    overlap = zeros(Float64,osize,osize); #overlap matrix Ne*Gaussians x Ne*Gaussians

    #flatbasisfunc=vcat(basisfunc...); #Need to improve
    flatbasisfunc = flattenbasisfunc(natoms,numbasisfunc,basisfunc);

    for n=1:osize
        for m=1:osize
            nbasis = flatbasisfunc[n];
            mbasis = flatbasisfunc[m];
            nprims = length(nbasis.alphas);
            mprims = length(mbasis.alphas);
            for np=1:nprims
                for mp=1:mprims
                    coefs = nbasis.coefs[np]*mbasis.coefs[mp];
                    gaussoverlap = getgaussoverlap(np,mp,nbasis,mbasis);
                    overlap[n,m] += coefs*gaussoverlap;
                end
            end # mprims
        end
    end #osize

    return overlap
end #buildelecoverlap

end #module