text:
#RSDeltaSigma: realize continuous-time NTF
#-------------------------------------------------------------------------------

#==Helper functions
===============================================================================#

#==realizeNTF_ct
===============================================================================#
"""`(ABCDc, tdac2) = realizeNTF_ct(NTF, form='FB', tdac, ordering=[1:n], bp=zeros(...), ABCDc)`

Realize a NTF with a continuous-time loop filter.

# Output
 - `ABCDc`: A state-space description of the CT loop filter
 - `tdac2`: A matrix with the DAC timings, including ones that were
   automatically added.

# Input Arguments
 - `NTF`: A noise transfer function in pole-zero form.
 - `form`: `= {:FB, :FF}` Specifies the topology of the loop filter.
 - `-->` For the `:FB` structure, the elements of `Bc` are calculated
   so that the sampled pulse response matches the `L1` impulse respnse.
 - `-->` For the `:FF` structure, `Cc` is calculated.
 - `tdac`: The timing for the feedback DAC(s).
 - `-->` If `tdac[1]>=1`, direct feedback terms are added to the quantizer.
 - `-->` Multiple timings (1 or more per integrator) for the `:FB`
   topology can be specified by making `tdac` a cell array, e.g.
   `tdac = { [1,2]; [1 2]; {[0.5 1],[1 1.5]}; []}`.
 - `-->` In above example, the first two integrators have dacs with `[1,2]`
   timing, the third has a pair of dacs, one with `[0.5 1]` timing and the
   other with `[1 1.5]` timing, and there is no direct feedback DAC to the
   quantizer
 - `ordering`: A vector specifying which NTF zero-pair to use in each
   resonator. Default is for the zero-pairs to be used in the order specified
   in the NTF.
 - `bp`: A vector specifying which resonator sections are bandpass.
   The default (zeros(...)) is for all sections to be lowpass.
 - `ABCDc`: The loop filter structure, in state-space form.
   If this argument is omitted, ABCDc is constructed according to `form`.
"""
function realizeNTF_ct(ntf, form=:FB, tdac=[0 1], ordering=1:0, bp=Int[], ABCDc=Float64[])
	(ntf_z, ntf_p) = _zpkdata(ntf, 1)
	ntf_z = deepcopy(ntf_z); ntf_p = deepcopy(ntf_p) #Don't modify original
	ntf_z = realfirst(ntf_z); ntf_p = realfirst(ntf_p) #VERIFYME: Expected by realizeNTF... here too?
	ntf_z_orig = deepcopy(ntf_z) #VERIFYME: Doesn't get modified to match legacy code. Intentional or omission??
	order = length(ntf_p)
	(order2, isodd, Δodd) = orderinfo(order)
	if isempty(ordering); ordering = 1:order2; end
	if isempty(bp); bp = zeros(Int, order2); end

	#compensate for limited accuracy of zero calculation
	ntf_z[findall( abs.(ntf_z .- 1) .< eps(1.0)^(1/(1+order)) )] .= 1.0
#ntf_z_orig=ntf_z #???

	multi_tmg = (eltype(tdac)<:Array)
	if multi_tmg
		if size(tdac) != (order+1,)
			error("For multi-array tdac, tdac be a vector of length `order+1`")
		end
		if form != :FB
			error("Currently only supporting form=:FB for multi-array tdac")
		end
	else
		if size(tdac) != (1, 2)
			#TODO: Update message
			error("For non cell array tdac, size(tdac) must be (1, 2)")
		end
	end

	n_direct = 0
	n_extra = 0
	tdac2 = Float64[-1 -1]
	if !multi_tmg
		#Need direct terms for every interval of memory in the DAC
		n_direct = ceil(Int, tdac[2]) - 1
		if ceil(Int, tdac[2]) - floor(Int, tdac[1]) > 1
			n_extra = n_direct-1 #tdac pulse spans a sample point
		else
			n_extra = n_direct
		end
		tdac2 = Float64[ -1 -1;
			tdac;
			0.5*ones(n_extra,1)*[-1 1] .+ cumsum(ones(n_extra,2), dims=1)
		]
	end

	if isempty(ABCDc)
		ABCDc = zeros(order+1, order+2);
		#Stuff the A portion
		if isodd
			ABCDc[1,1] = real( log( ntf_z[1] ) )
			ABCDc[2,1] = 1
		end
		for i = 1:order2
			n = bp[i]
			i1 = 2*i + Δodd - 1
			zi = 2*ordering[i] + Δodd - 1
			w = abs( angle( ntf_z[zi] ) )
			ABCDc[i1 .+ [0 1 2], i1 .+ [0 1]] = [
				0  -w^2;
				1   0;
				n  1-n;
			]
		end
		ABCDc[1,order+1] = 1
		ABCDc[1,order+2] = -1 #2006.10.02 Changed to -1 to make FF STF have +ve gain at DC
	end

	Ac = ABCDc[1:order, 1:order]
	local Bc, Cc, Dc, tp

	if :FB == form
		Cc = ABCDc[[order+1],1:order]
		if !multi_tmg
			Bc = [eye(order) zeros(order,1)]
			Dc = [zeros(1,order) 1]
			tp = repeat(tdac,order+1,1)
		else	#Assemble tdac2, Bc and Dc
			tdac2 = Float64[-1 -1]
			Bc = ones(order, 0) #0-column
			Bci = [eye(order) zeros(order,1)]
			Dc = ones(1, 0)
			Dci = [zeros(1,order) 1]
			for i in 1:length(tdac)
				tdi = tdac[i]
				if (eltype(tdi)<:Array)
					for j in 1:length(tdi)
						tdj = tdi[j]
						tdac2 = [tdac2; tdj]
						Bc = [Bc Bci[:,i]]
						Dc = [Dc Dci[:,i]]
					end
				elseif !isempty(tdi)
					tdac2 = [tdac2; tdi]
					Bc = [Bc Bci[:,i]]
					Dc = [Dc Dci[:,i]]
				end
			end
			tp = tdac2[2:end,:]
		end
	elseif :FF == form
		Cc = [eye(order); zeros(1,order)]
		Bc = [-1; zeros(order-1,1)]
		Dc = [zeros(order,1); 1]
		tp = tdac #2008-03-24 fix from Ayman Shabra
	else
		error("Sorry, no code for form ", form)
	end

	#Sample the L1 impulse response
	n_imp = ceil(Int, 2*order + maximum(tdac2[:,2]) + 1)
	y = impL1(ntf,n_imp)

	sys_c = _ss( Ac, Bc, Cc, Dc )
	yy = pulse(sys_c, tp, 1.0, 1.0*n_imp, true)
	yy = squeeze(yy)
	#Endow yy with n_extra extra impulses.
	#These will need to be implemented with n_extra extra DACs.
	#!! Note: if t1=int, matlab says pulse(sys) @t1 ~=0
	#!! This code corrects this problem.
	if n_extra>0
		y_right = padb([zeros(1,n_extra+1); eye(n_extra+1)], n_imp+1)
		#Replace the last column in yy with an ordered set of impulses
		yy = [yy[:,1:end-1] y_right[:,end:-1:1]]
	end

	#Solve for the coefficients:
	x = yy\y
	if norm(yy*x - y) > 1e-4
		@warn("Pulse response fit is poor.")
	end
	local Bc2, Dc2
	if :FB == form
		if !multi_tmg
			Bc2 = [ x[1:order] zeros(order,n_extra) ]
			Dc2 = x[order+1:end]'
		else
			BcDc = [Bc; Dc]
			i = findall(BcDc .!= 0)
			BcDc[i] = x
			Bc2 = BcDc[1:end-1,:]
			Dc2 = BcDc[[end],:]
		end
	elseif :FF == form
		Bc2 = [Bc zeros(order,n_extra)]
		Cc = x[1:order]'
		Dc2 = x[order+1:end]'
	else
		error("No code for form \"$form\"")
	end

	Dc1 = 0
	Dc = [Dc1 Dc2]
	Bc1 = [1; zeros(order-1,1)]
	Bc = [Bc1 Bc2]

	#Scale Bc1 for unity STF magnitude at f0
	fz = angle.(ntf_z_orig)/(2π)
	f1 = fz[1]
	ibz = findall(abs.(fz .- f1) .<= abs.(fz .+ f1))
	fz = fz[ibz]
	f0 = mean(fz)
	if minimum(abs.(fz)) < 3*minimum(abs.(fz .- f0))
		f0 = 0
	end

	#MALaforge: FIXME: ss->zpk conversion appears not to be numerically stable
	L0c = _zpk(_ss(Ac,Bc1,Cc,Dc1))
	#Hard-coded approximation of L0c's poly-ratio representation:
	msg  = "FIXME: Using hard-coded system to get realizeNTF_ct to work."
	msg *= "\nThis is NOT ok; results are probably not useful."
	@warn(msg)
	L0c = ControlSystems.tf([1], [1, 0, 0.0057829713287632966, 0])
	(_z, _p, _k) = ControlSystems.zpkdata(L0c)
	L0c = _zpk(_z, _p, _k)
	G0 = evalTFP(L0c, ntf, f0)
	Bc[:,1] = Bc[:,1]/abs.(G0)

	ABCDc = [Ac Bc; Cc Dc]

	return (ABCDc, tdac2)
end

#Last line


text:
module Ranks
    @enum Rank begin
        NULL=0
        TWO=2
        THREE=3
        FOUR=4
        FIVE=5
        SIX=6
        SEVEN=7
        EIGHT=8
        NINE=9
        TEN=10
        JACK=11
        QUEEN=12
        KING=13
        ACE=14
    end
end

module Suits
    @enum Suit begin
        NULL=0
        CLUBS=1
        DIAMONDS=2
        HEARTS=3
        SPADES=4
    end
end

module Rankings
    @enum Ranking begin
        HIGH=1
        PAIR=2
        DUBS=3
        TRIP=4
        STRAIGHT=5
        FLUSH=6
        FULL=7
        QUAD=8
        SFLUSH=9
        ROYAL=10
    end
end

struct Card
    rank::Ranks.Rank
    suit::Suits.Suit
end

function compare(c1::Card, c2::Card; tie_break=false)
    if c1.rank > c2.rank
        1
    elseif c1.rank < c2.rank
        -1
    elseif ! tie_break
        0
    elseif c1.suit > c2.suit
        1
    elseif c1.suit < c2.suit
        -1
    else
        0
    end
end

mutable struct Deck
    cards::Vector{Card}
    top::Int
end

mutable struct Hole
    playerID::Int
    card1::Card
    card2::Card
end

mutable struct Player
    playerID::Int
    hole::Hole
end

mutable struct Board
    flop1::Card
    flop2::Card
    flop3::Card
    turn::Card
    river::Card
end

mutable struct Hand
    playerID::Int
    ranking::Rankings.Ranking
    card1::Card
    card2::Card
    card3::Card
    card4::Card
    card5::Card
end

function compare(h1::Hand, h2::Hand; tie_break=false)
    if h1.ranking > h2.ranking
        1
    elseif h1.ranking < h2.ranking
        -1
    elseif compare(h1.card1, h2.card1, tie_break=tie_break) != 0
        compare(h1.card1, h2.card1, tie_break=tie_break)
    elseif compare(h1.card2, h2.card2, tie_break=tie_break) != 0
        compare(h1.card2, h2.card2, tie_break=tie_break)
    elseif compare(h1.card3, h2.card3, tie_break=tie_break) != 0
        compare(h1.card3, h2.card3, tie_break=tie_break)
    elseif compare(h1.card4, h2.card4, tie_break=tie_break) != 0
        compare(h1.card4, h2.card4, tie_break=tie_break)
    elseif compare(h1.card5, h2.card5, tie_break=tie_break) != 0
        compare(h1.card5, h2.card5, tie_break=tie_break)
    else
        0
    end
end

ext: jl

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
# SPDX-License-Identifier: MIT

using CUDA
using CUDA.CUSPARSE
using MarkovModels
using SparseArrays
using LinearAlgebra

function makefsm(T, S)
    fsm = FSM{LogSemifield{T}}()

    prev = addstate!(fsm, pdfindex = 1)
    setinit!(prev)
    link!(fsm, prev, prev)
    for s in 2:S
        state = addstate!(fsm, pdfindex = s)
        link!(fsm, prev, state)
        link!(fsm, state, state)
        prev = state
    end
    setfinal!(prev)
    renormalize!(fsm)
    fsm
end

function main(T, B, S, N; warmup)
    # Generate some pseudo-(log)likelihood
    lhs = ones(T, S, N, B)

    # Build the Finite State Machine (FSM).
    fsm = makefsm(T, S)

    precision = T == Float64 ? "double" : "single"

    cfsm = compile(fsm, allocator = spzeros)
    if warmup stateposteriors(cfsm, lhs) end
    etime = time_ns()
    stateposteriors(cfsm, lhs)
    etime = (time_ns() - etime) / 1e9
    println("julia\t$precision\t$B\t$S\t$N\tsemifield\tsparse\tcpu\t$etime")

    cfsm = compile(fsm, allocator = spzeros) |> gpu
    lhs = CuArray(lhs)
    if warmup stateposteriors(cfsm, lhs) end
    etime = time_ns()
    stateposteriors(cfsm, lhs)
    etime = (time_ns() - etime) / 1e9
    println("julia\t$precision\t$B\t$S\t$N\tsemifield\tsparse\tgpu\t$etime")
end

const T = Float32
const Bs = 10:10
const Ss = 300:300:1500
const Ns = 500:500:10000

warmup = true
for B in Bs
    for S in Ss
        for N in Ns
            if N < S continue end
            main(T, B, S, N; warmup)
            global warmup = false
        end
    end
end

max_stars_repo_name: MartinKocour/MarkovModels

max_issues_repo_name: MartinKocour/MarkovModels

max_forks_repo_name: MartinKocour/MarkovModels


text:
using MeshCat
using GeometryBasics

using Dojo
using DojoEnvironments

vis = Visualizer()
render(vis)

include(joinpath(@__DIR__, "..", "..", "..", "DojoEnvironments/src/utilities.jl"))
include(joinpath(@__DIR__, "..", "..", "..", "DojoEnvironments/src/panda/methods/initialize.jl"))

mech = get_panda(
    timestep=0.01,
    gravity=0.0 * 9.81,
    spring=1.0,
    damper=100.0,
    # contact=true,
    model_type=:no_end_effector,
    limits=false,
    object_type=:box,
)

nu = input_dimension(mech)
function ctrl!(m, k; kp=10.0, kd=5.0)
    pref = [0,-0.8,0,1.6,0,-2.4]#,0.5 * π]
    y = get_minimal_state(m)
    p = y[12 .+ (1:2:12)]
    v = y[12 .+ (2:2:12)]
    u = [
            5.0 * kp * [1,1.0e-1,1.0e-1,1.0e-1,1.0e-1,1.0e-1] .* (pref - p) - kd * v;
            [0.0; 0.0; -1.0 * get_body(m, :box).mass * 9.81];
            zeros(3);
    ]
    set_input!(m, u*m.timestep)
    return nothing
end

joint_angles = [1.3,-0.85,0.3,1.3,0.3,-2.0]#,0.5 * π]
joint_velocities = [0,0,0,0,0,0]#,0]

y = zeros(length(get_minimal_state(mech)))
y[1:3] = [0.5; 0.5; 0.125]
for i = 1:6
    y[12 + (i - 1) * 2 .+ (1:2)] = [joint_angles[i]; joint_velocities[i]]
end
set_minimal_state!(mech, y)
# set_minimal_state!(mech, copy(y_start[12 .+ (1:18)]))
# set_minimal_state!(mech, copy(y_start))

box_body = get_body(mech, :box)
box_body.state.q2 =  RotZ(0.25 * π)

storage = simulate!(mech, 5.0,
    ctrl!,
    verbose=true,
    opts=SolverOptions(verbose=true, rtol=1.0e-6, btol=1.0e-6));

vis, anim = visualize(mech, storage, vis=vis, show_contact=false)

setobject!(vis["ee"], GeometryBasics.Sphere(Point3f0(0),
        convert(Float32, mech.contacts[1].model.collision.radius_sphere)),
        MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))

for t in 1:length(storage)
    MeshCat.atframe(anim, t) do
        MeshCat.settransform!(vis["ee"],
        MeshCat.Translation(storage.x[5][t] + vector_rotate(mech.contacts[1].model.collision.origin_sphere, storage.q[5][t])))
    end
end
MeshCat.setanimation!(vis, anim)


text:
## vertices properties ##

"""
    degree

Degree of a given vertex in a given hypergraph. Also referred to as valency.

"""
function degree end

"""
    degrees

Degrees of several vertices in a given hypergraph.

"""
function degrees end

"""
    indegree

"""
function indegree end

"""
    indegrees

"""
function indegrees end

"""
    outdegree

"""
function outdegree end

"""
    outdegrees

"""
function outdegrees end

# internal function to deal with all cases
function _degree(x::T, v, f::Function) where {T<:AbstractHyperGraph}
    has_vertex(x, v) ? sum([num_has_vertex(e, v, f=f) for e in hyperedges(x)]) : error("vertex $v not found in hypergraph vertices")
end

# the degree of a vertex in a hypergraph
degree(x::T, v) where {T<:AbstractHyperGraph} = _degree(x, v, vertices)

# the in- and outdegree of a vertex in an oriented hypergraph
@traitfn indegree(x::T::IsOriented, v) where {T<:AbstractHyperGraph} = _degree(x, v, tgt)
@traitfn outdegree(x::T::IsOriented, v) where {T<:AbstractHyperGraph} = _degree(x, v, src)

# degree functions for vectors of vertices
degrees(x::T, vs::AbstractVector) where {T<:AbstractHyperGraph} = [degree(x, v) for v in vs]
@traitfn indegrees(x::T::IsOriented, vs::AbstractVector) where {T<:AbstractHyperGraph} = [indegree(x, v) for v in vs]
@traitfn outdegrees(x::T::IsOriented, vs::AbstractVector) where {T<:AbstractHyperGraph} = [outdegree(x, v) for v in vs]

# degree functions for hypergraphs
degrees(x::T) where {T<:AbstractHyperGraph} = degrees(x, vertices(x))
@traitfn indegrees(x::T::IsOriented) where {T<:AbstractHyperGraph} = indegrees(x, vertices(x))
@traitfn outdegrees(x::T::IsOriented) where {T<:AbstractHyperGraph} = outdegrees(x, vertices(x))

## hyperedges properties ##

# length, cardinality
Base.length(e::T) where {T<:AbstractHyperEdge} = length(vertices(e))
"""
    cardinality(e)

Cardinality of edge `e`, defined as its number of endpoints.
"""
cardinality(e::T) where {T<:AbstractHyperEdge} = length(e)
cardinalities(es::AbstractVector{T}) where {T<:AbstractHyperEdge} = cardinality.(es)
cardinalities(x::T) where {T<:AbstractHyperGraph} = cardinalities(hyperedges(x))

# number of source and of target objects
@traitfn nsrcs(e::T::IsOriented) where {T<:AbstractHyperEdge} = length(objects(e.src))
@traitfn ntgts(e::T::IsOriented) where {T<:AbstractHyperEdge} = length(objects(e.tgt))

## hypergraphs properties ##

# number of vertices, of hyperedges
nv(x::T) where {T<:AbstractHyperGraph} = length(vertices(x))
nhe(x::T) where {T<:AbstractHyperGraph} = length(hyperedges(x))

# extending Base.size
Base.size(x::T) where {T<:AbstractHyperGraph} = (nv(x), nhe(x))

# maximum degree of a hypergraphs
"""
    delta(x)

The maximum degree delta or Δ of a hypergraph `x` is the maximum degree over its vertices [Berge1989, Hellmuth2012].
"""
delta(x::T) where {T<:AbstractHyperGraph} = maximum(degrees(x))
Δ(x::T) where {T<:AbstractHyperGraph} = delta(x)

# rank, co-rank, and order of hypergraph
"""
    rank(x)

Rank of hypergraph `x`, defined as the maximum cardinality over its hyperedges [Berge1989, Hellmuth2012, Bretto2013, Burgio2020].
"""
rank(x::T) where {T<:AbstractHyperGraph} = maximum(cardinalities(x))
"""
    co_rank(x)

Co-rank of hypergraph `x`, defined as the minimum cardinality over its hyperedges [Hellmuth2012, Bretto2013, Burgio2020]. Sometimes referred to as anti-rank, e.g. in [Berge1989].
"""
co_rank(x::T) where {T<:AbstractHyperGraph} = minimum(cardinalities(x))
order(x::T) where {T<:AbstractHyperGraph} = nv(x)

# size and volume of hypergraph
hypergraph_size(x::T) where {T<:AbstractHyperGraph} = sum(cardinalities(hyperedges(x)))
volume(x::T) where {T<:AbstractHyperGraph} = sum(degrees(x))


text:
import DataStructures: BinaryMinHeap
using LinearAlgebra: norm

struct LSHBase
    alphas::Matrix{Float64} # Columns are linear equations.
    betas::Array{Float64, 2} # One per hash / one per alpha column.
    w::Array{Float64, 2} # Normalizing constant per hash.
end

"""
A gaussian LSH initialization method.

args:
    - n_dims: the number of dimensions in the input data.
    - w: the denominator and the max adjust value. May be an array of numbers.
    - n_hashes: number of hashes.
"""
function LSHBase(n_dims::Int, w=1.0, n_hashes=0)
    if n_hashes == 0
        n_hashes = Int(max(round(sqrt(n_dims)), 2))
    end
    alphas = randn(n_dims, n_hashes)

    # Normalize the alphas.
    alpha_norms = reshape([norm(col, 2) for col in eachcol(alphas)], 1, size(alphas, 2))
    unit_alphas = alphas ./ alpha_norms
    alphas = unit_alphas

    # Make sure that w is a column vector or an int.
    if isa(w, Float64) || isa(w, Int)
        w = reshape([w for _ in 1:n_hashes], 1, n_hashes)
    else
        w = reshape(w, 1, n_hashes)
    end
    betas = rand(1, n_hashes) .* w

    LSHBase(alphas, betas, w)
end

"""
The hash function.
"""
function hash(lsh::LSHBase, data)
    if length(size(data)) == 1
        data = reshape(data, 1, maximum(size(data)))
        hash_address = map(floor, ((data * lsh.alphas) .+ lsh.betas) ./ lsh.w)
        return vec(map(Int, hash_address))
    else
        hash_address = map(floor, ((data * lsh.alphas) .+ lsh.betas) ./ lsh.w)
        try
            return map(Int, hash_address)
        catch
            throw(error(join(["EXCEPTION", hash_address, data, lsh.alphas, lsh.betas, lsh.w], ", ")))
        end
    end
end

"""
Function to aggregate raw data into hash buckets.
"""
function reservoir_function(lshb::LSHBase, data::Array{Float64, 2})
    # Update vals
    vals = Dict{Array{Int, 1}, Any}()
    hashed = hash(lshb, data)
    for i in 1:size(hashed, 1)
        if !haskey(vals, hashed[i, :])
            vals[hashed[i, :]] = Array{Any, 1}()
        end
        push!(vals[hashed[i, :]], data[i, :])
    end
    return vals
end

"""
Function to aggregate raw data into hash buckets as average and count lists.
"""
function mean_count_function(lshb::LSHBase, data::Array{Float64, 2})
    vals = Dict{Array{Int, 1}, Any}()
    hashed = hash(lshb, data)
    for i in 1:size(hashed, 1)
        if !haskey(vals, hashed[i, :])
            vals[hashed[i, :]] = Array{Union{Int, Array{Float64, 1}}, 1}([[0. for _ in data[1, :]], 0])
        end
        vals[hashed[i, :]][1] += data[i, :]
        vals[hashed[i, :]][2] += 1
    end
    return vals
end

mutable struct LSH
    lshb::LSHBase
    aggregate_function::Function
    # Hash address to value.
    vals::Dict{Array{Int, 1}, Any}
end

"""
Create a new randomized LSH
args:
    - data: raw input data.
    - aggregate_function: a function to aggregate raw data points that fall
      into the same hash bucket.
    - n_hashes: the number of hashes.

"""
function LSH(data::Array{Float64, 2}, aggregate_function::Function=reservoir_function, n_hashes=0, w=.1)
    best = LSHBase(size(data, 2), w, n_hashes)

    # Update vals
    vals = aggregate_function(best, data)
    return LSH(best, aggregate_function, vals)
end

"""
Create a new LSH using the values in the provided LSHBase.
"""
function LSH(lshb::LSHBase, data::Array{Float64, 2}, aggregate_function::Function=reservoir_function)
    vals = aggregate_function(lshb, data)
    return LSH(lshb, aggregate_function, vals)
end

function hash(lsh::LSH, v)
    return hash(lsh.lshb, v)
end

"""
Calculate the score of a probe as the sum over the distance of the point from
the hash boundary + delta.

args:
    - query: a hash address.
    - delta: array of ints to apply to the hash address to get the probe
             location.
"""
function probe_score(lsh::LSHBase, query, delta)
    hashed_float = (query' * lsh.alphas .+ lsh.betas) ./ lsh.w
    if delta < 0
        return (map(floor, hashed_float) - hashed_float .+ [delta for _ in lsh.w]).^2
    else
        return (map(ceil, hashed_float) - hashed_float .+ [delta for _ in lsh.w]).^2
    end
end

"""
Generate query deltas in order of their probe_score and evaluate the probes for
neighbors. Return the neighbor hash addresses.

args
    lsh: mplsh object
    query: vector of data in original space
    target_neighbors: how many neighbors to return
    n_probes: the max number of probes to try
    probe_distance: positive int how many buckets away from query to include
                    in search.
"""
function get_neighbor_addresses(lsh::LSH, query, target_neighbors=1, probe_distance=1, n_probes=typemax(Int))
    neighbor_addresses = []
    query_address = hash(lsh, query)

    # Compute the probe_scores for all nonzero integer deltas in
    # [-probe_distance, probe_distance] (for each pair (i, delta_i))
    # perturbation_set will be (score, variable_index, probe_distance)
    perturbation_set = []
    for pd in 1:probe_distance
        perturbation_set = vcat(perturbation_set, vec([(score, hash_index, pd) for (hash_index, score) in enumerate(probe_score(lsh.lshb, query, pd))]))
        perturbation_set = vcat(perturbation_set, vec([(score, hash_index, -pd) for (hash_index, score) in enumerate(probe_score(lsh.lshb, query, -pd))]))
    end

    # Sort the pairs by score.
    scores = [pair[1] for pair in perturbation_set]
    perturbation_set = perturbation_set[sortperm(scores)]

    # Initialize the heap. The value [1] indicates that we'll start by trying
    # the first element in perturbation_set, which is sorted by cost.
    h = BinaryMinHeap{Tuple{Float64, Array{Int, 1}}}()
    push!(h, (perturbation_set[1][1], [1]))

    deltas = []
    count = 0
    n_probes_so_far = 0
    while length(deltas) < n_probes && length(h) > 0
        n_probes_so_far += 1
        count += 1
        # Get the next best probe.
        ai = pop!(h)

        # Create the hash address for this probe.
        probe = copy(query_address)
        for ind in ai[2]
            probe[perturbation_set[ind][2]] += perturbation_set[ind][3]
        end

        # If the bucket this delta indicates is not empty, add it to the neighbors set.
        if haskey(lsh.vals, probe)
            push!(neighbor_addresses, probe)
        end

        # If we've got the desired number of actual nonempty neighbors, return them.
        if length(neighbor_addresses) >= target_neighbors
            return neighbor_addresses
        end

        push!(deltas, ai[2])

        # Shift
        set = copy(ai[2])
        score = ai[1]
        # Remove the score for the current end perturbation.
        score -= perturbation_set[ai[2][end]][1]
        # Add the score for the new end perturbation.
        next = ai[2][end] + 1
        if next <= length(perturbation_set)
            score += perturbation_set[next][1]
            set[end] += 1
            push!(h, (score, set))
        end

        # Expand
        set = copy(ai[2])
        set = vcat(set, set[end] + 1)

        # Already updated the set at this point, offsets reflect.
        if length(set) <= length(lsh.lshb.w) - 1 && set[end] <= length(perturbation_set)
            score = ai[1] + perturbation_set[set[end]][1]
            push!(h, (score, set))
        end
    end

    return neighbor_addresses
end

"""
Find the approximate nearest neighbors with some constraints.

use:
Retrieve the first non-empty neighboring bucket's contents that can be found
adjacent to the query location.
    nearest_neighbors(lsh, [.23, -1, 3])
Retrieve the first non-empty neighboring bucket's contents within +-100 buckets
of the query location.
    nearest_neighbors(lsh, [.23, -1, 3], 1, probe_distance=100)
Retrieve the first 10 non-empty neighboring bucket's contents within +-100 buckets
of the query location.
    nearest_neighbors(lsh, [.23, -1, 3], 10, probe_distance=100)
And so on....

args
    lsh: mplsh object
    query: vector of data in original space
    n_neighbors: how many neighbors to return
    n_probes: the max number of probes to try
    probe_distance: positive int how many buckets away from query to include
                    in search.
"""
function nearest_neighbors(lsh::LSH, query, n_neighbors=1, probe_distance=3, n_probes=typemax(Int))
    nns = get_neighbor_addresses(lsh, query, n_neighbors, probe_distance, n_probes)
    neighbors = []
    for nn in nns
        neighbors = vcat(neighbors, lsh.vals[nn])
    end
    return Array{Float64, 2}(hcat(neighbors...)')
end

max_stars_repo_name: arider/MPLSH

max_issues_repo_name: arider/MPLSH

max_forks_repo_name: arider/MPLSH


text:
function p2fheader(n, es)
  ftype = c_ftype(float_for_posit_size(n))

  pterm = "const $(p(n,es)) pval"

  "  $(ftype) $(to_f(n,es,pterm))"
end

function f2pheader(n, es)
  ftype = c_ftype(float_for_posit_size(n))
  fterm = "const $ftype fval"
  "  $(p(n,es)) $(to_p(n,es,fterm))"
end

function generate_posit_conv_h(io, posit_defs)
  write(io,"""
  #ifndef __POSIT_CONV_H
  #define __POSIT_CONV_H

  #include \"posit.h\"

  #ifdef __cplusplus
  extern "C" {
  #endif

  """)
  for n in sort(collect(keys(posit_defs)))
    for es in posit_defs[n]
      write(io, p2fheader(n, es))
      write(io, ";\n")
      write(io, f2pheader(n, es))
      write(io, ";\n")
    end

    write(io, "\n")
  end
  write(io,"""
  #ifdef __cplusplus
  }
  #endif
  #endif
  """)
end


text:
#=

A call to app/healSelection using tiny image.

>julia
julia> include("test/apps/healSelection.jl")
julia> @time test(image)
=#

# For "load" local  files
using FileIO
# The tested item
include("../../../src/apps/healSelection.jl")

# mock test data
# global so it persists and out of the test
image = load("/home/bootch/git/resynthesizerJulia/test/data/in/tinylighthouse.png")

# selection mask for corpus
mask = falses(size(image))
mask[3:4, 2:3] .= true

function test(image, mask)
    # frisketDepth 1
    healSelection(image, mask, 1)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
@info "Caesar is loading tools related to AprilTags.jl."

using .AprilTags
# using Rotations
# using .Manifolds

include("Pose2AprilTag4Corners.jl")

#


text:
"""
    WeightedPoint3

An object of the type [`WeightedPoint3`](@ref) is a tuple of a three-dimensional
point and a scalar weight.

See also: [`Point3`](@ref)
"""
WeightedPoint3

"""
    WeightedPoint3(::Origin = ORIGIN)

Introduces a weighted point with Cartesian coordinates `(0, 0, 0)` and weight
`0`.
"""
WeightedPoint3(::Origin = ORIGIN)

"""
    WeightedPoint3(p::Point3)

Introduces a weighted point from point `p` and weight `0`.
"""
WeightedPoint3(p::Point3)

"""
    WeightedPoint3(p::Point3, w::Real)

Introduces a weighted point from point `p` and weight `w`.
"""
WeightedPoint3(p::Point3, w::Real)

WeightedPoint3(p, w::Real) = WeightedPoint3(p, convert(FT, w))

"""
    WeightedPoint3(x::Real, y::Real)

Introduces a weighted point with coordinates `x`, `y`, and weight `0`.
"""
WeightedPoint3(x::Real, y::Real) =
    WeightedPoint3(convert(FT, x), convert(FT, y))

"""
    point(p::WeightedPoint3)

Returns the point of the weighted point.
"""
point(p::WeightedPoint3)

"""
    weight(p::WeightedPoint3)

Returns the weight of the weighted point.
"""
weight(p::WeightedPoint3)

"""
    ==(p::WeightedPoint3, q::WeightedPoint3)

Test for equality.

Two points are equal, iff their ``x``, ``y`` and ``z`` coordinates are equal.
The point can be compared with [`ORIGIN`](@ref).

!!! warn "Warning"

    Comparison and equality (==, <, etc.) currently operate directly on the
    underlying bare point. Consequently:

    ```julia-repl
    julia> p = Point3(1.0, 2.0)
    PointC2(1, 2)

    julia> wp, wq = WeightedPoint3.((p, p), (1.0, 2.0))    # same bare point here,
    (Weighted_pointC2(1, 2, 1), Weighted_pointC2(1, 2, 2)) # but different weights

    julia> wp == wq
    true
    ```
"""
==(p::WeightedPoint3, q::WeightedPoint3)

"""
    hx(p::WeightedPoint3)

Returns the homogeneous ``x`` coordinate.
"""
hx(p::WeightedPoint3)

"""
    hy(p::WeightedPoint3)

Returns the homogeneous ``y`` coordinate.
"""
hy(p::WeightedPoint3)

"""
    hz(p::WeightedPoint3)

Returns the homogeneous ``z`` coordinate.
"""
hz(p::WeightedPoint3)

"""
    hw(p::WeightedPoint3)

Returns the homogenizing coordinate.
"""
hw(p::WeightedPoint3)

"""
    x(p::WeightedPoint3)

Returns the Cartesian ``x`` coordinate, that is [`hx()`](@ref)/[`hw()`](@ref).
"""
x(p::WeightedPoint3)

"""
    y(p::WeightedPoint3)

Returns the Cartesian ``y`` coordinate, that is [`hy()`](@ref)/[`hw()`](@ref).
"""
y(p::WeightedPoint3)

"""
    z(p::WeightedPoint3)

Returns the Cartesian ``z`` coordinate, that is [`hz()`](@ref)/[`hw()`](@ref).
"""
z(p::WeightedPoint3)

@doc raw"""
    homogeneous(p::WeightedPoint3, i::Integer)

Returns the `i`ᵗʰ homogeneous coordinate of `p`.

!!! info "Precondition"

    ``0 ≤ i ≤ 2``.
"""
homogeneous(p::WeightedPoint3, i::Integer)

@doc raw"""
    cartesian(p::WeightedPoint3, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of `p`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``.
"""
cartesian(p::WeightedPoint3, i::Integer)

"""
    dimension(p::WeightedPoint3)

Returns the dimension (the constant 2).
"""
dimension(p::WeightedPoint3)

"""
    bbox(p::WeightedPoint3)

Returns a bounding box containing `p`.
"""
bbox(p::WeightedPoint3)

"""
    transform(p::WeightedPoint3, t::AffTransformation3)

Returns the weighted point obtained by applying `t` on `p`.
"""
transform(p::WeightedPoint3, t::AffTransformation3)


text:
using DiffEqGPU, StochasticDiffEq, Test, DiffEqNoiseProcess
using CUDA
using Random

function f!(du,u,p,t)
 @inbounds begin
     du[1] = p[1]*u[1]
     du[2] = -p[1]*u[2]
 end
 nothing
end

function scalar_noise!(du,u,p,t)
 @inbounds begin
     du[1] = p[2]
     du[2] = -p[2]
 end
 nothing
end

function prob_func(prob, i, repeat)
    Random.seed!(seeds[i])
    W = WienerProcess(0.0,0.0,0.0)
    remake(prob,noise=W)
end

# fix seeds
seed = 100
Random.seed!(seed)
numtraj= 100
seeds = rand(UInt, numtraj)
W = WienerProcess(0.0,0.0,0.0)

CUDA.allowscalar(false)
u0 = [-0.1f0,0.0f0]
tspan = (0.0f0,100.0f0)
p = [1.9f0,0.1f0]

prob = SDEProblem(f!,scalar_noise!,u0,tspan,p, noise=W)
ensembleprob = EnsembleProblem(prob, prob_func = prob_func)

@info "scalar SDE"

#Performance check with nvvp
# CUDAnative.CUDAdrv.@profile
# check either on CPU with EnsembleCPUArray() or on GPU with EnsembleGPUArray()
#@test_nowarn sol = @time solve(ensembleprob,DRI1(),EnsembleCPUArray(),trajectories=numtraj)
sol = @time solve(ensembleprob,DRI1NM(),EnsembleGPUArray(),trajectories=numtraj)

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
println("Reloading Modules...")
include("sem.jl")

println("Reading mesh...")
msh = Nektar.nektar("meshes/circle.rea")

println("Creating Basis...")
Q = 15
b = sem.Basis2d(Q)

println("Creating Elements...")
elems = [sem.Element2d(b, el, msh.curves) for el in msh.elems];

Nel = length(elems)
x = sem.new_var(Q, Nel)
y = sem.new_var(Q, Nel)

for e = 1:Nel
    x[:,:,e] = elems[e].x
    y[:,:,e] = elems[e].y
end

fun(x,y) = 1 .+ x.*x .+ y.*y/2 .+ x.*y

dfdx(x,y) = 2x .+ y
dfdy(x,y) = y .+ x

f = fun(x, y)
dx = dfdx(x, y)
dy = dfdy(x,y)
ws = zeros(Q,Q)
dxc = similar(f)
dyc = similar(f)
sem.deriv_x!(elems, f, dxc);
sem.deriv_y!(elems, f, dyc);


text:
using Test
using ReversibleSeismic, NiLang.AD, NiLang
using CUDA

@testset "sum instr" begin
    out = 0.4
    x = randn(5)
    CUDA.allowscalar(true)

    @test PlusEq(sum)(out, abs2, CuArray(x))[1] ≈ sum(abs2, x) + 0.4
    @test check_grad(PlusEq(sum), (out, abs2, x); iloss=1)
    @test check_grad(MinusEq(sum), (out, abs2, x); iloss=1)
    @test check_grad(PlusEq(sum), (out, abs2, x |> CuArray); iloss=1)
    @test check_grad(MinusEq(sum), (out, abs2, x |> CuArray); iloss=1)
    CUDA.allowscalar(false)
end

@testset "PlusEq(identity) instr - CUDA" begin
    x = randn(5) |> CuArray
    y = randn(5) |> CuArray
    @i function f(out, x, y)
        x .+= y
        out += sum((@skip! abs2), x)
    end
    out_, y_, x_ = (~f)(f(0.4, copy(y), copy(x))...)
    @test out_ ≈ 0.4
    @test x_ ≈ x
    @test y_ ≈ y
    CUDA.allowscalar(true)
    @test check_grad(f, (0.4, y, x); iloss=1)
    @test check_grad(~f, (0.4, y, x); iloss=1)

    @i function g(out, x, y)
        x[1:3] .+= y[1:3]
        out += sum((@skip! abs2), x)
    end
    out_, y_, x_ = (~g)(g(0.4, copy(y), copy(x))...)
    @test out_ ≈ 0.4
    @test x_ ≈ x
    @test y_ ≈ y
    @test check_grad(g, (0.4, y, x); iloss=1)
    CUDA.allowscalar(false)
end


text:
mutable struct Domain
    box_exist::Int32                          # 0 = not yet created, 1 = exists
    dimension::Int32                          # 2 = 2d, 3 = 3d
    nonperiodic::Int32                        # 0 = periodic in all 3 dims
                                            # 1 = periodic or fixed in all 6
                                            # 2 = shrink-wrap in any of 6
    triclinic::Int32    # 0 = orthog box, 1 = triclinic
    tiltsmall::Int32    # 1 if limit tilt, else 0
    box_change::Int32           # 1 if any of next 3 flags are set, else 0
    box_change_size::Int32      # 1 if box size changes, 0 if not
    box_change_shape::Int32     # 1 if box shape changes, 0 if not
    box_change_domain::Int32    # 1 if proc sub-domains change, 0 if not
    vdeform::Int32
    deform_flag::Int32        # 1 if fix deform exist, else 0
    deform_vremap::Int32      # 1 if fix deform remaps v, else 0
    deform_groupbit::Int32    # atom group to perform v remap for
    boundary::Vector{Int32} #[3*2];    # settings for 6 boundaries
                         # 0 = periodic
                         # 1 = fixed non-periodic
                         # 2 = shrink-wrap non-periodic
                         # 3 = shrink-wrap non-per w/ min
    pbc::Vector{Int32} #[3];          # xyz periodicity as array

    face::Array{Float64,2}#[6*3];   # unit normals of 6 box faces
    corners::Array{Float64,2}#[8*3];                     # 8 corner points
    faces::Array{Float64,2};#[6*4*3];   # 4 corner pts of 6 prism faces
    boxhi::Vector{Float64}#[3];   # orthogonal box global bounds
    boxlo::Vector{Float64}#[3];   # orthogonal box global bounds
    boxtilt::Vector{Float64}#[3]; # 3 tilt factors
    boxhi_lamda::Vector{Float64}#[3]; # lamda box = (0,1)
    boxlo_lamda::Vector{Float64}#[3]; # lamda box = (0,1)
    boxhi_bound::Vector{Float64}#[3];  # bounding box of tilted domain
    boxlo_bound::Vector{Float64}#[3];  # bounding box of tilted domain
    subhi::Vector{Float64}#[3]; # sub-box bounds on this proc
    sublo::Vector{Float64}#[3]; # sub-box bounds on this proc
    ssubhi::Vector{Float64}#[3]; # shifted sub-box on this proc
    ssublo::Vector{Float64}#[3]; # shifted sub-box on this proc
    bsubhi::Vector{Float64}#[3]; # bounding for sub-box on this proc
    bsublo::Vector{Float64}#[3]; # bounding for sub-box on this proc
    subhi_lamda::Vector{Float64}#[3];  # bounds of subbox in lamda
    sublo_lamda::Vector{Float64}#[3];  # bounds of subbox in lamda
    h::Vector{Float64}#[6]               # shape matrix in Voigt ordering
                                      # Voigt = xx,yy,zz,yz,xz,xy
    h_inv::Vector{Float64}#[6];  # inverse of h
    h_rate::Vector{Float64}#[6]; # rate of box size/shape change

    Domain() = new();
end

function printout(dom::Domain)
    display("boxlo: "); display(dom.boxlo);
    display("boxhi: "); display(dom.boxhi);
    display("boxtilt: "); display(dom.boxtilt);
    display("sublo: "); display(dom.sublo);
    display("subhi: "); display(dom.subhi);
    display("sublo_lamda: "); display(dom.sublo_lamda);
    display("subhi_lamda: "); display(dom.subhi_lamda);
    display("bsublo: "); display(dom.bsublo);
    display("bsubhi: "); display(dom.bsubhi);
    display("ssublo: "); display(dom.ssublo);
    display("ssubhi: "); display(dom.ssubhi);
    display("h: "); display(dom.h);
    display("h_inv: "); display(dom.h_inv);
end

function initdomain(boxhi, boxtilt, pbc=[1; 1; 1], bcs=[1; 1; 1; 1; 1; 1])

    dom = Domain();
    dim = 3
    dom.boxlo_lamda = zeros(3)
    dom.boxhi_lamda = zeros(3)
    dom.sublo_lamda = zeros(3)
    dom.subhi_lamda = zeros(3)
    dom.h = zeros(6)
    dom.h_inv = zeros(6)
    dom.h_rate = zeros(6)

    dom.box_change = 0
    dom.deform_flag = 0
    dom.vdeform = 0

    dom.dimension = dim
    dom.boxlo = [0.0; 0.0; 0.0];
    dom.boxhi = boxhi;
    dom.boxtilt = boxtilt;
    dom.pbc = pbc;
    dom.bcs = bcs;

    if (dom.boxtilt[1] == 0 && dom.boxtilt[2] == 0 && dom.boxtilt[3] == 0)
        dom.triclinic = 0;
    else
        dom.triclinic = 1;
    end

    for i=1:3
        dom.boxlo_lamda[i] = 0.0;
        dom.boxhi_lamda[i] = 1.0;
        dom.sublo_lamda[i] = 0.0;
        dom.subhi_lamda[i] = 1.0;
    end

    return dom
end

function setdomainstruct(pbc, a, b, c)

    dom = Domain();
    dim = 3
    dom.boxlo = zeros(dim)
    dom.boxhi = zeros(dim)
    dom.sublo = zeros(dim)
    dom.subhi = zeros(dim)
    dom.bsublo = zeros(dim)
    dom.bsubhi = zeros(dim)
    dom.ssublo = zeros(dim)
    dom.ssubhi = zeros(dim)
    dom.boxtilt = zeros(dim)
    dom.boxlo_bound = zeros(dim)
    dom.boxhi_bound = zeros(dim)
    dom.boxlo_lamda = zeros(dim)
    dom.boxhi_lamda = zeros(dim)
    dom.sublo_lamda = zeros(dim)
    dom.subhi_lamda = zeros(dim)
    dom.h = zeros(6)
    dom.h_inv = zeros(6)
    dom.h_rate = zeros(6)

    dom.box_change = 0
    dom.deform_flag = 0
    dom.vdeform = 0

    dom.dimension = dim
    dom.pbc = pbc
    dom.boxlo[1] = 0.0;
    dom.boxlo[2] = 0.0;
    dom.boxlo[3] = 0.0;
    dom.boxhi[1] = a[1];
    dom.boxhi[2] = b[2];
    dom.boxhi[3] = c[3];
    dom.boxtilt[1] = b[1];
    dom.boxtilt[2] = c[1];
    dom.boxtilt[3] = c[2];
    for i=1:3
        dom.boxlo_lamda[i] = 0.0;
        dom.boxhi_lamda[i] = 1.0;
        dom.sublo_lamda[i] = 0.0;
        dom.subhi_lamda[i] = 1.0;
    end

    if ((abs(dom.boxtilt[1]) > 1e-10) || (abs(dom.boxtilt[2]) > 1e-10) || (abs(dom.boxtilt[3]) > 1e-10))
        dom.triclinic = 1;
    else
        dom.triclinic = 0;
    end

    #dom.h, dom.h_inv, dom.boxlo_bound, dom.boxhi_bound, dom.boxhi, dom.boxlo, dom.boxtilt, dom.triclinic)
    dom = SetGlobalBox(dom);
    dom = SetLocalOrthBox(dom);
    dom = ShiftedSubbox(dom);
    dom = BoundingSubbox(dom);
    return dom
end

function setlatticedomain(dom::Domain, lat::Lattice, reg::Region)

    dim = 3
    dom.dimension = dim
    lat = SetupLattice(lat, unitstyle, dim);

    dom.triclinic = reg.triclinic;
    dom.boxtilt[0] = lat.spacing[0]*reg.boxtilt[0];
    dom.boxtilt[1] = lat.spacing[0]*reg.boxtilt[1];
    dom.boxtilt[2] = lat.spacing[1]*reg.boxtilt[2];
    for i = 1:3
        dom.boxhi[i] = lat.spacing[i]*reg.boxhi[i];
        dom.boxlo[i] = lat.spacing[i]*reg.boxlo[i];
        dom.boxlo_lamda[i] = 0.0;
        dom.boxhi_lamda[i] = 1.0;
        dom.sublo_lamda[i] = 0.0;
        dom.subhi_lamda[i] = 1.0;
    end

    dom = SetGlobalBox(dom);
    dom = SetLocalOrthBox(dom);
    dom = ShiftedSubbox(dom);
    dom = BoundingSubbox(dom);
    lat = LatticeBoundingBox(lat, dom)
    lat = latticebounds(lat)

    # cpuSetGlobalBox(dom.h, dom.h_inv, dom.boxlo_bound, dom.boxhi_bound,
    #         dom.boxhi, dom.boxlo, dom.boxtilt, dom.triclinic);
    # cpuSetLocalOrthBox(dom.subhi, dom.sublo, dom.boxhi, dom.boxlo,
    #         dom.subhi_lamda, dom.sublo_lamda, 3);
    # epsilon = zeros(3)
    # for i=1:3
    #     epsilon[i] = 1e-6*(dom.boxhi[i]-dom.boxlo[i]);
    # end
    # cpuShiftedSubbox(dom.ssublo, dom.ssubhi, dom.boxlo, dom.boxhi,
    #         dom.boxlo_lamda, dom.boxhi_lamda, dom.sublo, dom.subhi,
    #         dom.sublo_lamda, dom.subhi_lamda, epsilon, pbc, dom.triclinic);
    # cpuBoundingSubbox(dom.bsublo, dom.bsubhi, dom.sublo, dom.subhi,
    #         dom.sublo_lamda, dom.subhi_lamda, dom.boxlo, dom.h, dom.triclinic);
    # cpuLatticeBoundingBox(lat.sublo, lat.subhi, dom.bsublo, dom.bsubhi,
    #         lat.primitive, lat.rotaterow, lat.primitinv, lat.rotatecol,
    #         lat.origin, lat.spacing, lat.scale);
    # lat.latticebounds();

    return dom, lat
end

function CreateLatticeAtom(lat)
    dim = 3
    natom = lat.natom
    x = zeros(dim*natom)
    t = zeros(natom)
    x, t = AtomLattice(x, t, lat, 0, dim)
    x = reshape(x, (dim, natom))
    return x, t
end

function SetAtomTypeFraction(x, t, atomtypefraction, seed, save)
    dim, nlocal = size(x)
    newtype = Int32(atomtypefraction[1]);
    fraction = atomtypefraction[2];
    seed0 = Int32(atomtypefraction[3]);
    t, seed, save, count = SetAtomType(x, fraction, t, seed, save, seed0, newtype, dim, nlocal)
    return t, seed, save, count
end

function CreateInitialVelocity(x, t, createvelocity, atommass, second, seed, save)

    dim, nlocal = size(x)
    t_desired = createvelocity[1];
    seed0 = Int32(createvelocity[2]);
    dist_flag = Int32(createvelocity[3]);
    sum_flag = Int32(createvelocity[4]);
    loop_flag = Int32(createvelocity[5]);
    #momentum_flag = Int32(createvelocity[6]);
    #rotation_flag = Int32(createvelocity[7]);
    mvv2e = createvelocity[8]
    boltz = createvelocity[9]

    mpiRank = 0
    natom=nlocal;
    tdof = (nlocal-1)*dim;
    tfactor = mvv2e / (tdof * boltz);

    ilist = zeros(nlocal)
    amass = zeros(nlocal)
    vcm = zeros(3)

    for i=0:nlocal
         ilist[i+1] = i;
    end
    v = 0.0*x

    VelocityCreate(x, v, atommass, second, seed, save, ilist, t, seed0,
            sum_flag, dist_flag, loop_flag, dim, mpiRank, nlocal, natom);
    masstotal = ComputeMass(amass, atommass, t, ilist, nlocal);
    ComputeVCM(vcm, v, atommass, masstotal, ilist, t, dim, nlocal);
    VelocityZeroMomentum(v, vcm, dim, nlocal);
    temp = ComputeTempScalar(v, atommass, tfactor, t, ilist, dim, nlocal);
    cpuArrayMultiplyScalar(v, sqrt(t_desired/temp), nlocal*dim);

    return v, vcm, masstotal, temp
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using YaoBase

export PhaseGate, phase

"""
    PhiGate

Global phase gate.
"""
mutable struct PhaseGate{T <: Real} <: PrimitiveBlock{1}
    theta::T
end

"""
    phase(theta)

Returns a global phase gate. Defined with following matrix form:

```math
exp(iθ) \\mathbf{I}
```

# Example

You can create a global phase gate with a phase (a real number).

```jldoctest
julia> phase(0.1)
phase(0.1)
```
"""
phase(θ::AbstractFloat) = PhaseGate(θ)
phase(θ::Real) = phase(Float64(θ))

mat(::Type{T}, gate::PhaseGate) where T = exp(T(im * gate.theta)) * IMatrix{2, T}()

# parametric interface
niparams(::Type{<:PhaseGate}) = 1
getiparams(x::PhaseGate) = x.theta
setiparams!(r::PhaseGate, param::Real) = (r.theta = param; r)

YaoBase.isunitary(r::PhaseGate) = true
Base.adjoint(blk::PhaseGate) = PhaseGate(-blk.theta)
Base.copy(block::PhaseGate{T}) where T = PhaseGate{T}(block.theta)
Base.:(==)(lhs::PhaseGate, rhs::PhaseGate) = lhs.theta == rhs.theta

cache_key(gate::PhaseGate) = gate.theta


text:
using CSV
using DataFrames

function createUSDF()

    csvfile = CSV.File("states_daily_4pm_et.csv")
    df = DataFrame(csvfile)
    idx = typeof.(df[:,:recovered]) .== Int
    dfrecovered = df[idx,:]
    states = unique(dfrecovered,:state)[:,:state]

    for i = 1:length(states)
        dfstatetmp = df[df[:,:state] .== states[i],:]
        dfstate = select(dfstatetmp,[:date,:positive,:death,:recovered])
        for i = 1:size(dfstate,2)
            idx = typeof.(dfstate[:,i]) .== Int
            dfstate[.!idx,i] .= 0
        end
        rename!(dfstate,[Symbol("date"), Symbol("cum_case"), Symbol("cum_death"), Symbol("cum_recover")])
        dfstate[:,:active] = dfstate[:,:cum_case] - dfstate[:,:cum_death] - dfstate[:,:cum_recover]

        sort!(dfstate,:cum_case)

        CSV.write("covid_timeseries_US_" * states[i] * ".csv",dfstate)
    end

end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
ElementCardType = CardType{:element}()
const ElementCard = NastranCard{ElementCardType}

type CBEAM <: ElementCard
    id::Int64
    prop_id::Int64
end
function convert(::Type{CBEAM},card)
    id = card[2]::Int64
    prop_id = card[3]::Int64
    CBEAM(id,prop_id)
end

type CTRIAR <: ElementCard
    id::Int64
    prop_id::Int64
    thickness1::Nullable{Float64}
    thickness2::Nullable{Float64}
    thickness3::Nullable{Float64}
end
function convert(::Type{CTRIAR},card)
    id = card[2]::Int64
    prop_id = card[3] == "" ? id : card[3]::Int64
    tflag = 0
    thickness1 = nothing
    thickness2 = nothing
    thickness3 = nothing
    if length(card) > 10
        tflag = card[11] == "" ? 0 : card[11]::Int64
        thickness1 = card[12] == "" ? nothing : card[12]::Float64
        thickness2 = card[13] == "" ? nothing : card[13]::Float64
        thickness3 = card[14] == "" ? nothing : card[14]::Float64
    else
    end
    CTRIAR(id,prop_id,thickness1,thickness2,thickness3)
end

type CQUADR <: ElementCard
    id::Int64
    prop_id::Int64
end
function convert(::Type{CQUADR},card)
    id = card[2]::Int64
    prop_id = card[3]::Int64
    CQUADR(id,prop_id)
end

type CROD <: ElementCard
    id::Int64
    prop_id::Int64
end
function convert(::Type{CROD},card)
    id = card[2]::Int64
    prop_id = card[3] == "" ? id : card[3]::Int64
    CROD(id,prop_id)
end

type CBUSH <: ElementCard
    id::Int64
end
function convert(::Type{CBUSH},card)
    id = card[2]::Int64
    CBUSH(id)
end

type RBE2 <: ElementCard
    id::Int64
end
function convert(::Type{RBE2},card)
    id = card[2]::Int64
    RBE2(id)
end

type RBE3 <: ElementCard
    id::Int64
end
function convert(::Type{RBE3},card)
    id = card[2]::Int64
    RBE3(id)
end

type CELAS1 <: ElementCard
    id::Int64
end
function convert(::Type{CELAS1},card)
    id = card[2]::Int64
    CELAS1(id)
end

type CELAS2 <: ElementCard
    id::Int64
end
function convert(::Type{CELAS2},card)
    id = card[2]::Int64
    CELAS2(id)
end

type CONM2 <: ElementCard
    id::Int64
    grid_id::Int64
    csys_id::Int64
    mass::Float64
    x::Float64
    y::Float64
    z::Float64
    I11::Float64
    I21::Float64
    I22::Float64
    I31::Float64
    I32::Float64
    I33::Float64
end
function convert(::Type{CONM2},card)
    id = card[2]::Int64
    grid_id = card[3]::Int64
    csys_id = card[4] == "" ? 0 : card[4]::Int64
    mass = card[5]::Float64
    x = card[6] == "" ? 0.0 : card[6]::Float64
    y = card[7] == "" ? 0.0 : card[7]::Float64
    z = card[8] == "" ? 0.0 : card[8]::Float64
    I11 = I21 = I22 = I31 = I32 = I33 = 0.0
    if length(card) > 9
        I11 = card[10] == "" ? 0.0 : card[10]::Float64
        I21 = card[11] == "" ? 0.0 : card[11]::Float64
        I22 = card[12] == "" ? 0.0 : card[12]::Float64
        I31 = card[13] == "" ? 0.0 : card[13]::Float64
        I32 = card[14] == "" ? 0.0 : card[14]::Float64
        I33 = card[15] == "" ? 0.0 : card[15]::Float64
    end
    CONM2(id,grid_id,csys_id,mass,x,y,z,I11,I21,I22,I31,I32,I33)
end

type CAERO1 <: ElementCard
    id::Int64
end
function convert(::Type{CAERO1},card)
    id = card[2]::Int64
    CAERO1(id)
end

type SPLINE1 <: ElementCard
    id::Int64
end
function convert(::Type{SPLINE1},card)
    id = card[2]::Int64
    SPLINE1(id)
end

type SPLINE2 <: ElementCard
    id::Int64
end
function convert(::Type{SPLINE2},card)
    id = card[2]::Int64
    SPLINE2(id)
end

type CHEXA <: ElementCard
    id::Int64
    prop_id::Int64
    point1::Int64
    point2::Int64
    point3::Int64
    point4::Int64
    point5::Int64
    point6::Int64
    point7::Int64
    point8::Int64
    point9::Int64
    point10::Int64
    point11::Int64
    point12::Int64
    point13::Int64
    point14::Int64
    point15::Int64
    point16::Int64
    point17::Int64
    point18::Int64
    point19::Int64
    point20::Int64
end
function convert(::Type{CHEXA},card)
    CHEXA(card[2]::Int64,
          card[3]::Int64,
          card[4]::Int64,
          card[5]::Int64,
          card[6]::Int64,
          card[7]::Int64,
          card[8]::Int64,
          card[9]::Int64,
          card[10]::Int64,
          card[11]::Int64,
          card[12] == "" ? 0 : card[12]::Int64,
          card[13] == "" ? 0 : card[13]::Int64,
          card[14] == "" ? 0 : card[14]::Int64,
          card[15] == "" ? 0 : card[15]::Int64,
          card[16] == "" ? 0 : card[16]::Int64,
          card[17] == "" ? 0 : card[17]::Int64,
          length(card) < 18 || card[18] == "" ? 0 : card[18]::Int64,
          length(card) < 18 || card[19] == "" ? 0 : card[19]::Int64,
          length(card) < 18 || card[20] == "" ? 0 : card[20]::Int64,
          length(card) < 18 || card[21] == "" ? 0 : card[21]::Int64,
          length(card) < 18 || card[22] == "" ? 0 : card[22]::Int64,
          length(card) < 18 || card[23] == "" ? 0 : card[23]::Int64)
end


text:
#' @@name AbstractDataFrame
#'
#' @@description
#'
#' An AbstractDataFrame is a Julia abstract type for which all concrete
#' types expose an database-like interface.
abstract AbstractDataFrame

#' @@name colmissing
#'
#' @@description
#'
#' Count the number of missing values in every column of an AbstractDataFrame.
#'
#' @@arg adf::AbstractDataFrame An AbstractDataFrame.
#'
#' @@return missing::Vector{Int} The number of missing values in each column.
#'
#' @@examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' colmissing(df)
function colmissing(adf::AbstractDataFrame) # -> Vector{Int}
    nrows, ncols = size(adf)
    missing = zeros(Int, ncols)
    for j in 1:ncols
        missing[j] = countna(adf[j])
    end
    return missing
end


text:
# %% imports
cd("Julia/CartPole/CMAES_colab")

using PyCall
gym = pyimport("gym")

include("StateSpikeGaussGrid.jl")
include("new_net_arch.jl")
include("CartpoleGridGenerator.jl")

# %% sim functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl > sr) ? action = 0 : action = 1
    if sl == sr
        action = rand(0:1)
    end
    return action, sl, sr
end

# %% fitness function play 100 episodes and check les contraintes a chaque épisode

function play_episode(net::Network, n_steps::Int64 = 200)
    env = gym.make("CartPole-v1")
    env.seed(0)
    obs = env.reset()
    grid_cells = grid_gen()
    win_size = 10 # ms
    reward_tot = 0
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        spikes = latency_spikes(obs, grid_cells, win_size)
        influence = 2
        @inbounds for msec in 1:win_size
            spiked = step!(net, spikes[msec], false)
            if msec >= length(net.arch)-1
                count_tot .+= influence .* count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
            end
        end
        n_left = div(net.arch[end],2)
        n_right = net.arch[end] - n_left
        action, sl, sr = left_right(count_tot, n_left, n_right)
        obs_new, reward, done, _ = env.step(action)
        reward_tot += reward
        obs = obs_new
    end
    env.close()
    env = nothing
    Base.GC.gc()
    return reward_tot
end

# %% objective function ( fitness + contraintes ) ( create new network and play 100 episodes with fitness function )

function rescale(genes::Array{Float64}, tunekeys::Array{String})
    cfg = YAML.load(open("cfg_cmaes.yml"))
    N = length(genes)
    scaled_genes = zeros(N)
    n_param = length(tunekeys)
    for idx in 1:N
        if idx <= n_param
            down = cfg[tunekeys[idx]*"_d"]
            up = cfg[tunekeys[idx]*"_u"]
            scaled_genes[idx] = ((genes[idx]+1)*0.5)*(up-down) + down
        elseif idx > n_param
            scaled_genes[idx] = (genes[idx]+1)*0.5
        end
    end
    return scaled_genes
end

function generate(genes_p::Array{Float64}, tunekeys::Array{String})
    params = Dict(tunekeys .=> rescale(genes_p,tunekeys))
    keys = ["tda", "wmax", "theta", "ttheta", "std", "m"]
    values = [0.995, 1, 0.0, 0.998, 0.1, 0.5]
    for idx in 1:length(keys)
        params[keys[idx]] = values[idx]
    end
    return params
end

function set_weight!(net::Network, genes_w::Array{Float64})
    w = net.s.e
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]])
            w[mat_idx][idx] = genes_w[idx+offset]
        end
    end
end

function objective(genes::Array{Float64})
    tunekeys = ["vthresh", "ae", "b", "c", "de",
                "wei", "wie", "ai", "di"]
    n_params = length(tunekeys)
    genes_p = genes[1:n_params]
    genes_w = genes[n_params + 1 : end]
    params = generate(genes_p, tunekeys)
    n = Network([120,60,80],params)
    set_weight!(n,genes_w)
    fit = play_episode(n)
    @show fit
    -fit
end

# %% tests

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: July 2020
## Update by Ziyi Yin

using JUDI, Images

export setup_model, parse_commandline, setup_geom

"""
Simple 2D model setup used for the tests.
"""

function smooth(v, sigma=3)
    return Float32.(imfilter(v,  Kernel.gaussian(sigma)))
end

"""
Sets up a simple 2D layered model for the wave equation operators tests
"""
function setup_model(tti=false, nlayer=2; n=(301, 151), d=(10., 10.))
    ## Set up model structure
    o = (0., 0.)

    v = ones(Float32,n) .* 1.5f0
    vp_i = range(1.5f0, 3.5f0, length=nlayer)
    for i in range(2, nlayer, step=1)
        v[:, (i-1)*Int(floor(n[2] / nlayer)) + 1:end] .= vp_i[i]  # Bottom velocity
    end

    v0 = smooth(v, 10)
    rho0 = (v .+ .5f0) ./ 2
    # Slowness squared [s^2/km^2]
    m = (1f0 ./ v).^2
    m0 = (1f0 ./ v0).^2
    dm = vec(m - m0)

    # Setup model structure
    if tti
        println("TTI Model")
        epsilon = smooth((v0[:, :] .- 1.5f0)/12f0, 3)
        delta =  smooth((v0[:, :] .- 1.5f0)/14f0, 3)
        theta =  smooth((v0[:, :] .- 1.5f0)/4, 3)
        model0 = Model(n, d, o, m0; epsilon=epsilon, delta=delta, theta=theta)
        model = Model(n, d, o, m; epsilon=epsilon, delta=delta, theta=theta)
    else
        model = Model(n,d,o,m,rho=rho0)
        model0 = Model(n,d,o,m0,rho=rho0)
    end

    return model, model0, dm
end

"""
Sets up a simple 2D acquisition for the wave equation operators tests
"""

function setup_geom(model; nsrc=1, tn=1500f0)
    ## Set up receiver geometry
    nxrec = model.n[1] - 2
    xrec = collect(range(model.d[1], stop=(model.n[1]-2)*model.d[1], length=nxrec))
    yrec = collect(range(0f0, stop=0f0, length=nxrec))
    zrec = collect(range(50f0, stop=50f0, length=nxrec))

    # receiver sampling and recording time
    T = tn   # receiver recording time [ms]
    dt = .75f0    # receiver sampling interval [ms]

    # Set up receiver structure
    recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=T, nsrc=nsrc)

    ## Set up source geometry (cell array with source locations for each shot)
    ex =  (model.n[1] - 1) * model.d[1]
    nsrc > 1 ? xsrc = range(.25f0 * ex, stop=.75f0 * ex, length=nsrc) : xsrc = .5f0 * ex
    xsrc = convertToCell(xsrc)
    ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
    zsrc = convertToCell(range(2*model.d[2], stop=2*model.d[2], length=nsrc))

    # Set up source structure
    srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=T)

    # setup wavelet
    f0 = 0.015f0     # MHz
    wavelet = ricker_wavelet(T, dt, f0)
    q = judiVector(srcGeometry, wavelet)

    ntComp = get_computational_nt(srcGeometry, recGeometry, model; dt=dt)
    info = Info(prod(model.n), nsrc, ntComp)

    return q, srcGeometry, recGeometry, info
end

### Process command line args
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--tti"
            help = "TTI, default False"
            action = :store_true
        "--fs"
            help = "Free surface, default False"
            action = :store_true
        "--isic"
            help = "Inverse scattering imaging condition, default False"
            action = :store_true
        "--nlayer", "-n"
            help = "Number of layers"
            arg_type = Int
            default = 2
        "--parallel", "-p"
            help = "Number of workers"
            arg_type = Int
            default = 1
    end
    return parse_args(s)
end


text:
using Test
using Forecast

@testset "loess" begin
    @test_throws AssertionError loess(rand(5),rand(5); d=3)

    x = loess(rand(10),rand(10); predict = rand(20))
    @test length(x) == 20

    x = loess(rand(10), rand(10); d=1, predict= rand(20))
    @test length(x) == 20

    x = loess(sin.(collect(1:10.)), sin.(collect(1:10.)); d=2)
    @test x == [0.8414709848078965, 0.9092974268256817, 0.14112000805986716, -0.7568024953079282, -0.9589242746631381, -0.2794154981989261, 0.6569865987187891, 0.9893582466233829, 0.4121184852417571, -0.5440211108893702]

    x = loess(sin.(collect(1:10.)), sin.(collect(1:10.)); d=1)
    @test x == [0.8414709848078965, 0.9092974268256815, 0.14112000805986719, -0.7568024953079282, -0.9589242746631382, -0.2794154981989259, 0.6569865987187891, 0.9893582466233817, 0.41211848524175654, -0.5440211108893699]

    x = loess(Int64.(round.(100*rand(100))),
              Int64.(round.(100*rand(100))); predict = Int64.(round.(100*rand(20))))
    @test length(x) == 20

    x = loess(1:2:200,
              vcat(rand(80),repeat([missing],20));
              predict = 1.:2:40)
    @test length(x) == 20

end


text:
# This script illustrates how to use latin hypercube sampling. In the first example, we produce a sample with a non-diagonal covariance matrix to illustrate that the latin property is lost for all dimensions but the first:
using MonteCarloMeasurements, LatinHypercubeSampling, Test
ndims  = 2
N      = 40   # Number of particles
ngen   = 2000 # How long to run optimization
X, fit = LHCoptim(N,ndims,ngen)
m, Σ   = [1,2], [2 1; 1 4] # Desired mean and covariance
particles = transform_moments(X, m, Σ)
@test mean(particles, dims=1)[:] ≈ m
@test cov(particles) ≈ Σ

p = Particles(particles)
plot(scatter(eachcol(particles)..., title="Sample"), plot(fit, title="Fitness vs. iteration"))
vline!(particles[:,1]) # First dimension is still latin
hline!(particles[:,2]) # Second dimension is not

# If we do the same thing with a diagonal covariance matrix, the latin property is approximately preserved in all dimensions provided that the latin optimizer was run sufficiently long.
m, Σ   = [1,2], [2 0; 0 4] # Desired mean and covariance
particles = transform_moments(X, m, Σ)
p         = Particles(particles)
plot(scatter(eachcol(particles)..., title="Sample"), plot(fit, title="Fitness vs. iteration"))
vline!(particles[:,1]) # First dimension is still latin
hline!(particles[:,2]) # Second dimension is not

# We provide a method which absolutely preserves the latin property in all dimensions, but if you use this, the covariance of the sample will be slighlty wrong
particles = transform_moments(X, m, Σ, preserve_latin=true)
p         = Particles(particles)
plot(scatter(eachcol(particles)..., title="Sample"), plot(fit, title="Fitness vs. iteration"))
vline!(particles[:,1]) # First dimension is still latin
hline!(particles[:,2]) # Second dimension is not
@test mean(particles, dims=1)[:] ≈ m
@test cov(particles) ≉ Σ # OBS not the not approx equal!

# We can also visualize the statistics of the sample
using StatsPlots
corrplot(particles)
#
plot(density(p[1]), density(p[2]))


text:
# Define the new set structure
struct MyNewSet <: InfiniteOpt.InfiniteScalarSet
    attr1::Float64 # REPLACE WITH DESIRED TYPE
    attr2::Float64 # REPLACE WITH DESIRED TYPE
    # ADD MORE ATTRIBUTES AS NEEDED
    # constructor
    function MyNewSet(attr1, attr2)
        # PUT APPROPRIATE CHECKS HERE AND OTHER NEEDED OPERATIONS
        return new(convert(Float64, attr1), convert(Float64, attr2))
    end
end

# Extend supports_in_set
function InfiniteOpt.supports_in_set(supports::Union{Real, Vector{<:Real}},
                                     set::MyNewSet)::Bool
    # DETERMINE IF SUPPORTS ARE IN THE DOMAIN OF `set`
    in_set = all(set.attr1 .<= supports .<= set.attr2) # REPLACE WITH ACTUAL CHECK
    return in_set
end

# Extend generate_support_values if possible
function InfiniteOpt.generate_support_values(set::MyNewSet;
                                             num_supports::Int = 10,
                                             sig_digits::Int = 5)::Tuple{Vector{<:Real}, DataType}
    # REPLACE BELOW WITH METHODS TO GENERATE `num_samples` with `sig_fig`
    supports = collect(range(set.attr1, stop = set.attr2, length = num_supports))
    return round.(supports, sigdigits = sig_digits), UniformGrid
end

# Extend JuMP.has_lower_bound (optional if the answer is always false)
function JuMP.has_lower_bound(set::MyNewSet)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.lower_bound(set::MyNewSet)
    return set.attr1 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_lower_bound if appropriate
function JuMP.set_lower_bound(set::MyNewSet, lower::Real)::MyNewSet
    return MyNewSet(lower, set.attr2) # REPLACE WITH ACTUAL CONSTRUCTOR
end

# Extend JuMP.has_supper_bound (optional if the answer is always false)
function JuMP.has_upper_bound(set::MyNewSet)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.upper_bound(set::MyNewSet)
    return set.attr2 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_upper_bound if appropriate
function JuMP.set_upper_bound(set::MyNewSet, upper::Real)::MyNewSet
    return MyNewSet(set.attr1, upper) # REPLACE WITH ACTUAL CONSTRUCTOR
end


text:
N = 100
h = 1.0/N
y = zeros(N)
y[1] = 0
t = linrange(0,2,N)

function func(t,x)
    return x^(1/3.0)
end

for i=1:N-1
    y[i+1] = y[i] + h*func(t[i], y[i])
end
println("y", "  ", "y(t)", "  ", "Error")
for i=1:N
    println(t[i]^(3/2.0), "  ", y[i], " ", abs(t[i]^1.5-y[i])/t[i]^1.5)
end

max_stars_repo_name: NeveIsa/hatex

max_issues_repo_name: NeveIsa/hatex

max_forks_repo_name: saketkc/hatex


text:
include((@__DIR__)*"/../preamble.jl")

include("absorbing_model/model_def.jl")
include("default_params.jl")

point_mass_initial_condition(n,rng) = (X = fill(eps(),n), φ = fill(1,n))
function exp_initial_condition(n,rng)
    i = rand(rng,1:2,n)
    X = Vector{Float64}(undef,n)
    for k in 1:n
        x = Inf
        while x>=1.0
            x = -log(rand(rng))
        end
        X[k]=x
    end
    return (X=X,φ=i)
end
rng = StableRNGs.StableRNG(16092021)

_first_exit = first_exit_x(0.5*eps(),1.0)
function first_exit_or_t10(model,SFM,SFM0)
    t1 = _first_exit(model,SFM,SFM0)
    i2 = SFM.t>10.0
    return (Ind=(i2||t1.Ind), SFM=t1.SFM)
end

sims = simulate(model,first_exit_or_t10,exp_initial_condition(n_sims,rng),rng)
pwd()
CSV.write("hitting_times_model/hitting_times/data/exp/sims.csv",
    DataFrame(t=sims.t,φ=sims.φ,X=sims.X,n=sims.n)
)


text:
import LinearAlgebra.BLAS: BlasInt, @blasfunc

for (hseqr,elty) in ((:zhseqr_,:ComplexF64),)
    @eval function hesseneigvals(M::Matrix{$elty})
        if isempty(M)
            return $elty[]
        end
       A=vec(M)

        N=size(M,1)
        info  =0

        ilo = 1; ihi = N; ldh=N;ldz=N;lwork = N

        z=zero($elty)
        work  = Array{$elty}(undef, N*N)
        w=Array{$elty}(undef, N)

        Ec='E'
        Nc='N'
        ccall((@blasfunc($hseqr),LAPACK.liblapack),
            Nothing,
            (Ref{UInt8}, Ref{UInt8},
        Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, #A
        Ref{BlasInt}, Ptr{$elty}, Ref{$elty}, #z
        Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}),
        Ec, Nc,
        N , ilo, ihi, A,
        ldh, w, z,
        ldz, work, lwork, info)
        w
    end
end

for (hseqr,elty) in ((:dhseqr_,:Float64),)
    @eval function hesseneigvals(M::Matrix{$elty})
        if isempty(M)
            return $elty[]
        end

        A=vec(M)

        N=size(M,1)
        info  =0

        ilo = 1; ihi = N; ldh=N;ldz=N;

        lwork = -1

        z=zero($elty)
        work  = Array{$elty}(undef,  1)
        wr=Array{$elty}(undef, N)
        wi=Array{$elty}(undef, N)

        Ec='E'
        Nc='N'
        for i=1:2
            ccall((@blasfunc($hseqr),LAPACK.liblapack),
                Nothing,
                (Ref{UInt8},Ref{UInt8},
            Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty}, #A
            Ref{BlasInt},Ptr{$elty},Ptr{$elty},Ref{$elty}, #z
            Ref{BlasInt},Ptr{$elty},Ref{BlasInt},Ref{BlasInt}),
            Ec, Nc,
            N , ilo, ihi, A,
            ldh, wr,wi, z,
            ldz, work, lwork, info)

            if lwork < 0
                lwork=Int(real(work[1]))
                work=Array{$elty}(undef, lwork)
            end
        end

        wr+im*wi
    end
end


text:
# Load the NPCF code
println("Loading NPCF module")
using NPCFs, Test, Distributed
Npart = 50 # number of particles
boxsize = 200 # input size

## Single thread runtest, comparing NPCF simple + pairwise outputs for different N and D values
println("Single-worker runtest")
for N in 2:5
    for D in 2:4
        println("Testing Cartesian N=$N, D=$D")
        npcf = NPCF(N=N,D=D,verb=true,periodic=true,volume=boxsize^D,coords="cartesian",r_min=20,r_max=95,nbins=6,lmax=1);

        pos = hcat(rand(Npart,D)*boxsize,ones(Npart)); # columns: [x, y, z, weight]

        npcf_output1 = compute_npcf_simple(pos, npcf);
        npcf_output2 = compute_npcf_pairwise(pos, npcf);

        # Compute the error
        err = sum(npcf_output1[:]-npcf_output2[:])/sum(npcf_output2[:]);
        println("Error is $err\n")

        @test abs(err)<1e-10
    end
end

## As above, but for spherical geometry in 2 dimensions
println("Spherical geometry runtest")
for N in 2:5
    println("Testing spherical N=$N, D=2")
    npcf = NPCF(N=N,D=2,verb=true,periodic=false,volume=4pi,coords="spherical",r_min=-0.8,r_max=0.8,nbins=6,lmax=1);

    phi_arr = rand(Npart)*2pi
    theta_arr = acos.(rand(Npart).*2 .-1)
    pos = hcat(phi_arr,theta_arr,ones(Npart));

    npcf_output1 = compute_npcf_simple(pos, npcf);
    npcf_output2 = compute_npcf_pairwise(pos, npcf);

    # Compute the error
    err = sum(npcf_output1[:]-npcf_output2[:])/sum(npcf_output2[:]);
    println("Error is $err\n")

    @test abs(err)<1e-10
end

## As above, but for an aperiodic geometry
println("Aperiodic box runtest")
for N in 2:5
    for D in 2:4
        println("Testing Cartesian N=$N, D=$D")
        npcf = NPCF(N=N,D=D,verb=true,periodic=false,volume=boxsize^D,coords="cartesian",r_min=20,r_max=95,nbins=6,lmax=1);

        pos = hcat(rand(Npart,D)*boxsize,ones(Npart)); # columns: [x, y, z, weight]

        npcf_output1 = compute_npcf_simple(pos, npcf);
        npcf_output2 = compute_npcf_pairwise(pos, npcf);

        # Compute the error
        err = sum(npcf_output1[:]-npcf_output2[:])/sum(npcf_output2[:]);
        println("Error is $err\n")

        @test abs(err)<1e-10
    end
end

## Multiple worker runtest
# Here we use N = D = 3, and run on two workers
println("Multi-worker runtest")
addprocs(2) # run on two workers
@everywhere using NPCFs
npcf = NPCF(N=3,D=3,verb=true,periodic=true,volume=boxsize^3,coords="cartesian",r_min=20,r_max=95,nbins=6,lmax=1);

pos = hcat(rand(Npart,3)*boxsize,ones(Npart)); # columns: [x, y, z, weight]

npcf_output1 = compute_npcf_simple(pos, npcf);
npcf_output2 = compute_npcf_pairwise(pos, npcf);

# Compute the error
err = sum(npcf_output1[:]-npcf_output2[:])/sum(npcf_output2[:]);
println("Error is $err\n")

@test abs(err)<1e-10

println("All tests passed!")


text:
using Test
#using DelimitedFiles, LinearAlgebra

#using StableRNGs
#rng = StableRNG(123)

import MLJBase
const Mlj = MLJBase
import Distributions
using BetaML

TESTRNG = FIXEDRNG # This could change...
#TESTRNG = StableRNG(123)

println("*** Testing Clustering...")

# ==================================
# New test
# ==================================
println("Testing initRepreserntative...")

Z₀ = initRepresentatives([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.6 38],2,initStrategy="given",Z₀=[1.7 15; 3.6 40])
@test isapprox(Z₀,[1.7  15.0; 3.6  40.0])

# ==================================
# New test
# ==================================
println("Testing kmeans...")

(clIdx,Z) = kmeans([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.3 38; 5.1 -2.3; 5.2 -2.4],3,initStrategy="grid",rng=copy(TESTRNG))
@test clIdx == [2, 2, 2, 2, 3, 3, 3, 1, 1]
#@test (clIdx,Z) .== ([2, 2, 2, 2, 3, 3, 3, 1, 1], [5.15 -2.3499999999999996; 1.5 11.075; 3.366666666666667 36.666666666666664])

# ==================================
# New test
# ==================================
println("Testing kmedoids...")
(clIdx,Z) = kmedoids([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.3 38; 5.1 -2.3; 5.2 -2.4],3,initStrategy="shuffle",rng=copy(TESTRNG))
@test clIdx == [1, 1, 1, 1, 2, 2, 2, 3, 3]

# ==================================
# New test
# ==================================
println("Testing mixture initialisation and log-pdf...")

m1 = SphericalGaussian()
m2 = SphericalGaussian([1.1,2,3])
m3 = SphericalGaussian(nothing,10.2)
mixtures = [m1,m2,m3]
X = [1 10 20; 1.2 12 missing; 3.1 21 41; 2.9 18 39; 1.5 15 25]
initMixtures!(mixtures,X,minVariance=0.25,rng=copy(TESTRNG))
@test sum([sum(m.μ) for m in mixtures]) ≈ 102.2
@test sum([sum(m.σ²) for m in mixtures]) ≈ 19.651086419753085
mask = [true, true, false]
@test lpdf(m1,X[2,:][mask],mask) ≈ -3.461552516784797

m1 = DiagonalGaussian()
m2 = DiagonalGaussian([1.1,2,3])
m3 = DiagonalGaussian(nothing,[0.1,11,25.0])
mixtures = [m1,m2,m3]
initMixtures!(mixtures,X,minVariance=0.25,rng=copy(TESTRNG))
@test sum([sum(m.σ²) for m in mixtures]) ≈ 291.27933333333334
@test lpdf(m1,X[2,:][mask],mask) ≈ -3.383055441795939

m1 = FullGaussian()
m2 = FullGaussian([1.1,2,3])
m3 = FullGaussian(nothing,[0.1 0.2 0.5; 0 2 0.8; 1 0 5])
mixtures = [m1,m2,m3]
initMixtures!(mixtures,X,minVariance=0.25,rng=copy(TESTRNG))
@test sum([sum(m.σ²) for m in mixtures]) ≈ 264.77933333333334
@test lpdf(m1,X[2,:][mask],mask) ≈ -3.383055441795939

# ==================================
# New test
# ==================================
println("Testing gmm...")
X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
clusters = gmm(X,3,verbosity=NONE, initStrategy="grid",rng=copy(TESTRNG))
@test isapprox(clusters.BIC,114.1492467835965)
#clusters.pₙₖ
#clusters.pₖ
#clusters.mixtures
#clusters.BIC

# ==================================
# New test
# ==================================
println("Testing predictMissing...")
X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
out = predictMissing(X,3,mixtures=[SphericalGaussian() for i in 1:3],verbosity=NONE, initStrategy="grid",rng=copy(TESTRNG))
@test isapprox(out.X̂[2,2],14.155186593170251)

X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
out2 = predictMissing(X,3,mixtures=[DiagonalGaussian() for i in 1:3],verbosity=NONE, initStrategy="grid",rng=copy(TESTRNG))
@test out2.X̂[2,2] ≈ 14.588514438886131

X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
out3 = predictMissing(X,3,mixtures=[FullGaussian() for i in 1:3],verbosity=NONE, initStrategy="grid",rng=copy(TESTRNG))
@test out3.X̂[2,2] ≈ 11.166652292936876

# ==================================
# NEW TEST
println("Testing MLJ interface for Clustering models....")
X, y                           = Mlj.@load_iris

model                          = KMeans(rng=copy(TESTRNG))
modelMachine                   = Mlj.machine(model, X)
(fitResults, cache, report)    = Mlj.fit(model, 0, X)
distances                      = Mlj.transform(model,fitResults,X)
yhat                           = Mlj.predict(model, fitResults, X)
acc = BetaML.accuracy(Mlj.levelcode.(yhat),Mlj.levelcode.(y),ignoreLabels=true)
@test acc > 0.8

model                          = KMedoids(rng=copy(TESTRNG))
modelMachine                   = Mlj.machine(model, X)
(fitResults, cache, report)    = Mlj.fit(model, 0, X)
distances                      = Mlj.transform(model,fitResults,X)
yhat                           = Mlj.predict(model, fitResults, X)
acc = BetaML.accuracy(Mlj.levelcode.(yhat),Mlj.levelcode.(y),ignoreLabels=true)
@test acc > 0.8

model                       =  GMMClusterer(mixtures=:diag_gaussian,rng=copy(TESTRNG))
modelMachine                =  Mlj.machine(model, X) # DimensionMismatch
(fitResults, cache, report) =  Mlj.fit(model, 0, X)
yhat_prob                   =  Mlj.predict(model, fitResults, X)  # Mlj.transform(model,fitResults,X)
# how to get this ??? Mlj.predict_mode(yhat_prob)
@test Distributions.pdf(yhat_prob[end],2) ≈ 0.5937443601647852

X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
X = Mlj.table(X)
model                       =  MissingImputator(rng=copy(TESTRNG))
modelMachine                =  Mlj.machine(model,X)
(fitResults, cache, report) = Mlj.fit(model, 0, X)
XD                          =  Mlj.transform(model,fitResults,X)
XDM                         =  Mlj.matrix(XD)
@test isapprox(XDM[2,2],15.441553354222702)
# Use the previously learned structure to imput missings..
Xnew_withMissing            = Mlj.table([1.5 missing; missing 38; missing -2.3; 5.1 -2.3])
XDNew                       = Mlj.transform(model,fitResults,Xnew_withMissing)
XDMNew                      =  Mlj.matrix(XDNew)
@test isapprox(XDMNew[1,2],13.818691793037452)

#=
# Marginally different
XD  == XDNew
Mlj.matrix(XD)  ≈ Mlj.matrix(XDNew)
XDM  == XDMNew

for r in 1:size(XDM,1)
    for c in 1:size(XDM,2)
        if XDM[r,c] != XDMNew[r,c]
            println("($r,$c): $(XDM[r,c]) - $(XDMNew[r,c])")
        end
    end
end
=#

#=
@test Mlj.mean(Mlj.LogLoss(tol=1e-4)(yhat_prob, y)) < 0.0002
Mlj.predict_mode(yhat_prob)
N = size(Mlj.matrix(X),1)
nCl = size(fitResults[2],1)
yhat_matrix = Array{Float64,2}(undef,N,nCl)
[yhat_matrix[n,c]= yhat_prob[n][c] for n in 1:N for c in 1:nCl]
yhat_prob[2]
 Mlj.matrix(yhat_prob)
acc = accuracy(Mlj.levelcode.(yhat),Mlj.levelcode.(y),ignoreLabels=true)
ynorm = Mlj.levelcode.(y)
accuracy(yhat_prob,ynorm,ignoreLabels=true)
@test acc > 0.8
=#

#=
using MLJBase, BetaML
y, _                        =  make_regression(1000, 3, rng=123);
ym                          =  MLJBase.matrix(y)
model                       =  GMM(rng=copy(BetaML.TESTRNG))
(fitResults, cache, report) =  MLJBase.fit(model, 0, nothing, y)
yhat_prob                   =  MLJBase.transform(model,fitResults,y)
yhat_prob                   =  MLJBase.predict(model, fitResults, y)

modelMachine                =  MLJBase.machine(model, nothing, y)
mach                        =  MLJBase.fit!(modelMachine)
yhat_prob                   =  MLJBase.predict(mach, nothing)
=#


text:
# Equality

# equals() now depends on == instead
# of the other way round.
function ==(x::Decimal, y::Decimal)
    a = norm(x)
    b = norm(y)
    a.c == b.c && a.q == b.q && a.s == b.s
end

function ==(x::Number, y::Decimal)
    decimal(x) == y
end

function ==(x::Decimal, y::Number)
    x == decimal(y)
end


text:
using Chemfiles

file = Trajectory("filename.xyz")
frame = read(file)

println("There are $(size(frame)) atoms in the frame")
pos = positions(frame);

# Do awesome science here with the positions

if has_velocities(frame)
    vel = velocities(frame)

    # If the file contains information about the
    # velocities, you will find them here.
end


text:
# Random access lists (Okasaki, 1995)
# (une structure utilisant de la récursion polymorphe)

struct Pair fst; snd end
struct Zero    s0    end
struct One  e; s1    end

function length(l::Nothing) 0                    end
function length(l::Zero   ) 2 * length(l.s0)     end
function length(l::One    ) 2 * length(l.s1) + 1 end

function get(l::Zero, i)
    x = get(l.s0, div(i, 2))
    if i % 2 == 0 x.fst else x.snd end
end
function get(l::One, i)
    if i == 0 return l.e end
    x = get(l.s1, div(i - 1, 2))
    if i % 2 == 1 x.fst else x.snd end
end

function cons(x, l::Nothing) One(x, nothing)                end
function cons(x, l::Zero   ) One(x, l.s0)                   end
function cons(x, l::One    ) Zero(cons(Pair(x, l.e), l.s1)) end

function sequence(i, j)
    if j < i
        nothing
    else
        cons(i, sequence(i + 1, j))
    end
end

s = sequence(34, 55)
for i = 0: 21
    println(get(s, i))
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ResourcesMoveInfo <: SwaggerModel
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    targetResourceGroup::Any # spec type: Union{ Nothing, String } # spec name: targetResourceGroup

    function ResourcesMoveInfo(;resources=nothing, targetResourceGroup=nothing)
        o = new()
        validate_property(ResourcesMoveInfo, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(ResourcesMoveInfo, Symbol("targetResourceGroup"), targetResourceGroup)
        setfield!(o, Symbol("targetResourceGroup"), targetResourceGroup)
        o
    end
end # type ResourcesMoveInfo

const _property_map_ResourcesMoveInfo = Dict{Symbol,Symbol}(Symbol("resources")=>Symbol("resources"), Symbol("targetResourceGroup")=>Symbol("targetResourceGroup"))
const _property_types_ResourcesMoveInfo = Dict{Symbol,String}(Symbol("resources")=>"Vector{String}", Symbol("targetResourceGroup")=>"String")
Base.propertynames(::Type{ ResourcesMoveInfo }) = collect(keys(_property_map_ResourcesMoveInfo))
Swagger.property_type(::Type{ ResourcesMoveInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ResourcesMoveInfo[name]))}
Swagger.field_name(::Type{ ResourcesMoveInfo }, property_name::Symbol) =  _property_map_ResourcesMoveInfo[property_name]

function check_required(o::ResourcesMoveInfo)
    true
end

function validate_property(::Type{ ResourcesMoveInfo }, name::Symbol, val)
end


text:
"""
$(SIGNATURES)
Get number of (solver) controls, accounting for minimum time and infeasible start
# Output
- m̄:  number of non infeasible controls (ie, system controls + time control if minimum time). System controls augmented by one if time is included as a control for minimum time problems.
- mm: total number of solver controls
"""
function get_num_controls(solver::Solver)
    n,m = get_sizes(solver)
    m̄ = m
    solver.state.minimum_time ? m̄ += 1 : nothing
    solver.state.infeasible ? mm = m̄ + n : mm = m̄
    return m̄, mm
end

"""
$(SIGNATURES)
Get number of (solver) stats, accounting for minimum time
# Output
- n̄:  number of non infeasible states (ie, system states + time state if minimum time). System state augmented by one if time is included as a state for minimum time problems.
- nn: total number of solver states
"""
function get_num_states(solver::Solver)
    n,m = get_sizes(solver)
    n̄ = n
    solver.state.minimum_time ? n̄ += 1 : nothing

    # TODO for now:
    nn = n̄
    return n̄, nn
end

"""
$(SIGNATURES)
    Compute the optimal control problem cost
"""
function cost(solver::Solver,vars::DircolVars)
    cost(solver,vars.X,vars.U)
end

function cost(solver::Solver, X::AbstractMatrix, U::AbstractMatrix)
    cost(solver, to_dvecs(X), to_dvecs(U))
end

function cost(solver::Solver,X::Trajectory,U::Trajectory)
    n,m,N = get_sizes(solver)
    J = 0.0
    costfun = solver.obj.cost
    for k = 1:N-1
        J += stage_cost(costfun,X[k][1:n],U[k][1:m])*solver.dt
    end
    J += stage_cost(costfun, X[N][1:n])
end

"""
$(SIGNATURES)
    Compute the optimal control problem unconstrained cost,
    including minimum time and infeasible controls
"""
function _cost(solver::Solver{M,Obj},res::SolverVectorResults,X=res.X,U=res.U) where {M, Obj<:Objective}
    # pull out solver/objective values
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    costfun = solver.obj.cost
    dt = solver.dt
    xf = solver.obj.xf

    J = 0.0
    for k = 1:N-1
        # Get dt if minimum time
        solver.state.minimum_time ? dt = U[k][m̄]^2 : nothing

        # Stage cost
        J += (stage_cost(costfun,X[k][1:n],U[k][1:m]))*dt

        # Minimum time cost
        solver.state.minimum_time ? J += solver.opts.R_minimum_time*dt : nothing

        # Infeasible control cost
        solver.state.infeasible ? J += 0.5*solver.opts.R_infeasible*U[k][m̄.+(1:n)]'*U[k][m̄.+(1:n)] : nothing
    end

    # Terminal Cost
    J += stage_cost(costfun, X[N][1:n])

    return J
end

"""
$(SIGNATURES)
    Compute the Augmented Lagrangian constraints cost
"""
function cost_constraints(solver::Solver, res::ConstrainedIterResults)
    N = solver.N
    J = 0.0
    if solver.opts.al_type == :default
        for k = 1:N
            J += 0.5*res.C[k]'*res.Iμ[k]*res.C[k] + res.λ[k]'*res.C[k]
        end
    elseif solver.opts.al_type == :algencan
        C,λ = res.C, res.λ
        for k = 1:N
            Iμk = Diagonal(res.μ[k])
            Ia = Diagonal(res.active_set[k])
            J += C[k]'Ia*(λ[k] + 0.5*Iμk*C[k]) - 0.5λ[k]'*(Iμk\(I-Ia))*λ[k]
        end
    end
    return J
end

function cost_constraints(solver::Solver, res::UnconstrainedIterResults)
    return 0.
end

function cost(solver::Solver, res::SolverIterResults, X=res.X, U=res.U)
    _cost(solver,res,X,U) + cost_constraints(solver,res)
end

"""
$(SIGNATURES)
    Calculate dynamics and constraint Jacobians (perform prior to the backwards pass)
"""
function update_jacobians!(res::ConstrainedIterResults, solver::Solver, jacobians::Symbol=:both)::Nothing
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    ∇f = discrete_dynamics_jacobian(solver)
    ∇c = constraint_jacobian(solver)

    dt = solver.dt

    for k = 1:N-1
        # Update discrete dynamics Jacobians
        jacobians in [:dynamics,:both] ? ∇f(res.fdx[k],res.fdu[k], res.X[k][1:n], res.U[k]) : nothing

        # Update constraint Jacobians
        if !solver.state.fixed_constraint_jacobians && jacobians in [:constraints,:both]

            ∇c(res.Cx[k], res.Cu[k], res.X[k],res.U[k])

            # Minimum time special case
            if solver.state.minimum_time
                # No equality constraint at first time step
                if k == 1
                    res.Cu[k][p,m̄] = 0.
                    res.Cx[k][p,n̄] = 0.
                end
                # Jacobian for x[n̄]_k+1 = u[m̄]_k
                res.fdu[k][n̄,m̄] = 1.
            end
        end
    end

    # Update terminal constraint Jacobian
    (!solver.state.fixed_terminal_constraint_jacobian && jacobians in [:constraints,:both]) ? ∇c(view(res.Cx[N], 1:p_N,1:n), res.X[N][1:n]) : nothing

    return nothing
end

function update_jacobians!(res::UnconstrainedIterResults, solver::Solver, jacobians::Symbol=:none)::Nothing
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    for k = 1:N-1
        # Update discrete dynamics Jacobians
        solver.Fd(res.fdx[k],res.fdu[k], res.X[k][1:n], res.U[k])
    end

    return nothing
end

function evaluate_trajectory(solver::Solver, X, U)
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    results = init_results(solver,X,U)
    update_jacobians!(results, solver)
    update_constraints!(results, solver)
    return results
end

"""
$(SIGNATURES)
Interpolate a trajectory using cubic interpolation
"""
function interp_traj(N::Int,tf::Float64,X::AbstractMatrix,U::AbstractMatrix)::Tuple{Matrix,Matrix}
    if isempty(X)
        X2 = X
    else
        X2 = interp_rows(N,tf,X)
    end
    U2 = interp_rows(N-1,tf,U)
    return X2, U2
end

interp_traj(N::Int,tf::Float64,X::Trajectory,U::Trajectory) = interp_traj(N,tf,to_array(X),to_array(U))

"""
$(SIGNATURES)
Interpolate the rows of a matrix using cubic interpolation
"""
function interp_rows(N::Int,tf::Float64,X::AbstractMatrix)::Matrix
    n,N1 = size(X)
    t1 = range(0,stop=tf,length=N1)
    t2 = collect(range(0,stop=tf,length=N))
    X2 = zeros(n,N)
    for i = 1:n
        interp_cubic = CubicSplineInterpolation(t1, X[i,:])
        X2[i,:] = interp_cubic(t2)
    end
    return X2
end

"""
$(SIGNATURES)
Generates the correctly sized input trajectory, tacking on infeasible and minimum
time controls, if required. Will interpolate the initial trajectory as needed.
# Arguments
* X0: Matrix of initial states. May be empty. If empty and the infeasible flag is set in the solver, it will initialize a linear interpolation from start to goal state.
* U0: Matrix of initial controls. May either be only the dynamics controls, or include infeasible and minimum time controls (as necessary).
"""
function get_initial_trajectory(solver::Solver, X0::Matrix{Float64}, U0::Matrix{Float64})
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    if size(U0,1) ∉ [m,mm]
        ArgumentError("Size of U0 must be either: system controls OR all expected solver controls (system + infeasible + minimum time)")
    end

    if N-1 != size(U0,2)
        @info "Interpolating initial guess"
        X0,U0 = interp_traj(N,solver.obj.tf,X0,U0)
    end

    if solver.state.minimum_time
        solver.state.infeasible ? sep = " and " : sep = " with "
        solve_string = sep * "Minimum Time..."

        # Initialize controls with sqrt(dt)
        if size(U0,1) == m
            U_init = [U0; ones(1,N-1)*sqrt(get_initial_dt(solver))]
        else
            U_init = copy(U0)
        end
    else
        solve_string = "..."
        U_init = U0
    end

    if solver.state.infeasible
        solve_string =  "Solving Constrained Problem with Infeasible Start" * solve_string

        # Generate infeasible controls
        if size(U0,1) == m
            ui = infeasible_controls(solver,X0,U_init)  # generates n additional control input sequences that produce the desired infeasible state trajectory
            U_init = [U_init; ui]  # augment control with additional control inputs that produce infeasible state trajectory
        end

        # Assign state trajectory
        if isempty(X0)
            X_init = line_trajectory(solver)
        else
            X_init = X0
        end
    else
        if solver.state.constrained
            solve_string = "Solving Constrained Problem" * solve_string
        else
            solve_string = "Solving Unconstrained Problem" * solve_string
        end
        X_init = zeros(n,N)
    end

    @info solve_string

    if solver.state.minimum_time
        X_init = [X_init; zeros(1,N)]
    end

    return X_init, U_init
end

"""
$(SIGNATURES)
    Regularization update scheme
        - see "Synthesis and Stabilization of Complex Behaviors through Online Trajectory Optimization"
"""
function regularization_update!(results::SolverResults,solver::Solver,status::Symbol=:increase)
    if status == :increase # increase regularization
        # @logmsg InnerLoop "Regularization Increased"
        results.dρ[1] = max(results.dρ[1]*solver.opts.bp_reg_increase_factor, solver.opts.bp_reg_increase_factor)
        results.ρ[1] = max(results.ρ[1]*results.dρ[1], solver.opts.bp_reg_min)
        if results.ρ[1] > solver.opts.bp_reg_max
            @warn "Max regularization exceeded"
        end
    elseif status == :decrease # decrease regularization
        results.dρ[1] = min(results.dρ[1]/solver.opts.bp_reg_increase_factor, 1.0/solver.opts.bp_reg_increase_factor)
        results.ρ[1] = results.ρ[1]*results.dρ[1]*(results.ρ[1]*results.dρ[1]>solver.opts.bp_reg_min)
    end
end

"""
$(SIGNATURES)
    Calculate the problem gradient using heuristic from iLQG (Todorov) solver
"""
function gradient_todorov(res::SolverVectorResults)
    N = length(res.X)
    maxes = zeros(N)
    for k = 1:N-1
        maxes[k] = maximum(abs.(res.d[k])./(abs.(res.U[k]).+1))
    end
    mean(maxes)
end

"""
$(SIGNATURES)
    Calculate the infinity norm of the gradient using feedforward term d (from δu = Kδx + d)
"""
function gradient_feedforward(res::SolverVectorResults)
    norm(res.d,Inf)
end

"""
$(SIGNATURES)
    Gradient of the Augmented Lagrangian: ∂L/∂u = Quuδu + Quxδx + Qu + Cu'λ + Cu'IμC
"""
function gradient_AuLa(results::SolverIterResults,solver::Solver)
    N = solver.N
    X = results.X; X_ = results.X_; U = results.U; U_ = results.U_

    Qx = results.bp.Qx; Qu = results.bp.Qu; Qxx = results.bp.Qxx; Quu = results.bp.Quu; Qux = results.bp.Qux

    gradient_norm = 0. # ℓ2-norm

    for k = 1:N-1
        δx = X_[k] - X[k]
        δu = U_[k] - U[k]
        tmp = Quu[k]*δu + Qux[k]*δx + Qu[k]

        gradient_norm += sum(tmp.^2)
    end

    return sqrt(gradient_norm)
end

function update_gradient(results,solver)
    if solver.opts.gradient_type == :todorov
        gradient = gradient_todorov(results)
    elseif solver.opts.gradient_type == :AuLa
        gradient = gradient_AuLa(results,solver)
    elseif solver.opts.gradient_type == :feedforward
        gradient = gradient_feedforward(results)
    end
    return gradient
end

function max_penalty(res::ConstrainedIterResults)
    N = length(res.X)
    val = 0
    if length(res.μ[1]) > 0
        val = maximum(maximum.(res.μ[1:N-1]))
    end
    if length(res.μ[N]) > 0
        val = max(val,maximum(res.μ[N]))
    end
    return val
end

max_penalty(res::UnconstrainedIterResults) = 0


text:
#=##############################################################################################
Filename: post_plot_lift_moment_distribution.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: `<: PostAction` function plots the lift distribution of all lifting surfaces at the specified steps
=###############################################################################################

"""
    post_plot_lift_moment_distribution(aircraft, parameters, step_range, step_symbol) <: PostAction

# Arguments

* `aircraft::Aircraft`: aircraft system struct
* `parameters <: Parameters`: inherits from the `Parameters` type; object containing data required by `<: Action` functions
* `step_range::AbstractArray`: range of simulation times
* `step_symbol::String`: defines the step, e.g. `alpha` or `time`

# Modifies:

* `aircraft::Aircraft`
* `parameters <: Parameters`

# Returns:

* `flag::Bool`: true if an action experiences errors
* saves plots in `parameters.plot_directory`

`parameters <: Parameters` requires the following elements:

* `surfacenames::Vector{String}`: names of each lifting surface to be shown in the legend
* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface
* `plot_directory::String`: path to the folder where plots will be saved
* `plot_base_name::String`: first part of saved figure file names
* `plot_extension::String`: extension of saved figure file names
* `plotstepi::Vector{Int}`: which steps at which to plot

"""
# * `cl_ylim::Vector{Float64}`: y axis limits for ploting c_l
# * `cd_ylim::Vector{Float64}`: y axis limits for ploting c_d
# * `cy_ylim::Vector{Float64}`: y axis limits for ploting c_y

function post_plot_lift_moment_distribution(aircraft, parameters, step_range, step_symbol)

    flags = Vector{Bool}(undef,length(step_range))
    for stepi = 1:length(step_range) # loop over all steps
        # run existing plot function
        flags[stepi] = plot_lift_moment_distribution(aircraft, parameters, nothing, nothing, step_range, stepi, step_symbol)
    end # loop over all steps

    return prod(flags)
end

"""
    post_plot_lift_moment_distribution(system, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft`: system to be simulated
* `step_range::AbstractArray`: defines each step of the simulation

# Returns:

* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface
* `surfacenames::Vector{String}`: names of each lifting surface to be shown in the legend
* `plot_directory::String`: path to the folder where plots will be saved
* `plot_base_name::String`: first part of saved figure file names
* `plot_extension::String`: extension of saved figure file names
* `plotstepi::Vector{Int}`: which steps at which to plot

"""
function post_plot_lift_moment_distribution(aircraft, step_range)

    # extract info
    nsteps = length(step_range)
    nwings = length(aircraft.wing_system.system.surfaces)
    # initialize parameters
    surfacenames = ["wing$i" for i in 1:nwings]
    cfs = [[zeros(3, length(aircraft.wing_system.system.surfaces[i])) for i in 1:nwings] for j in 1:nsteps]
    cms = deepcopy(cfs)
    plot_directory = ""
    plot_base_name = "default"
    plot_extension = ".pdf"
    plotstepi = deepcopy(step_range)

    return cfs, cms, surfacenames, plot_directory, plot_base_name, plot_extension, plotstepi
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
"""
	grdedit(cmd0::String="", arg1=nothing, kwargs...)

Reads the header information in a binary 2-D grid file and replaces the information with
values provided on the command line.

If single input is a G GMTgrid object, it will update the z_min|max values of the G.range member

Full option list at [`grdedit`]($(GMTdoc)grdedit.html)

Parameters
----------

- **A** | **adjust** :: [Type => Bool]

    If necessary, adjust the file’s x_inc, y_inc to be compatible with its domain.
    ($(GMTdoc)grdedit.html#a)
- **C** | **adjust** :: [Type => Bool]

    Clear the command history from the grid header.
    ($(GMTdoc)grdedit.html#c)
- **D** | **header** :: [Type => Str]    ``Arg = [+xxname][+yyname][+zzname][+sscale][+ooffset][+ninvalid][+ttitle][+rremark``

    Change these header parameters.
    ($(GMTdoc)grdedit.html#d)
- **E** | **header** :: [Type => Str]    ``Arg = [a|h|l|r|t|v]``

    Transform the grid in one of six ways and (for l|r|t) interchange the x and y information
    ($(GMTdoc)grdedit.html#e)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdedit(....) form.
    ($(GMTdoc)grdedit.html#g)
- $(GMT.opt_J)
- **N** | **replace** :: [Type => Str | Mx3 array]      ``Arg = replace=fname | replace=Array``

    Read the ASCII (or binary) file table and replace the corresponding nodal values in the
    grid with these x,y,z values. Alternatively, provide a Mx3 matrix with values to be changed.
    ($(GMTdoc)grdedit.html#n)
- $(GMT.opt_R)
- **S** | **wrap** :: [Type => Bool]

    For global, geographical grids only. Grid values will be shifted longitudinally according to
    the new borders given in ``limits`` (R option).
    ($(GMTdoc)grdedit.html#s)
- **T** | **toggle** :: [Type => Bool]

    Make necessary changes in the header to convert a gridline-registered grid to a pixel-registered
    grid, or vice-versa.
    ($(GMTdoc)grdedit.html#t)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)
"""
function grdedit(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdedit", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	arg2 = nothing
	(isa(arg1, GMTgrid) && length(kwargs) == 0) && (arg1.range[5:6] .= extrema(arg1); return arg1)  # Update the z_min|max

	cmd, = parse_common_opts(d, "", [:G :R :J :V_params :bi :di :e :f :w :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :adjust], [:C :clear_history], [:D :header], [:E :flip], [:S :wrap], [:T :toggle]])
	cmd, args, n, = add_opt(d, cmd, "N", [:N :replace], :data, Array{Any,1}([arg1, arg2]), (x="",))
    if (n > 0)  arg1, arg2 = args[:]  end

	common_grd(d, cmd0, cmd, "grdedit ", arg1, arg2)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdedit(arg1, cmd0::String=""; kw...) = grdedit(cmd0, arg1; kw...)


text:
module LoggingFacilities

using Logging
using LoggingExtras
using Dates
using JSON
using Pipe
using TimeZones

export  migrate, inject, remove, mutate, compose,
        KwargsProperty, MessageProperty, LevelProperty,
        LevelLocation, KwargsLocation, BeginningMessageLocation, EndingMessageLocation, WholeMessageLocation,
        OneLineTransformerLogger, TimestampTransformerLogger, JSONTransformerLogger,
        ColorMessageTransformerLogger, ColorSpec,
        FixedMessageWidthTransformerLogger,
        MessageOnlyLogger

include("types.jl")
include("transform.jl")
include("common.jl")
include("sink.jl")

end #module


text:
"""
condentropy \\

for the condEntropy function specifically, calculating conditional entropy \\
Arguments:	marg_dist, marginal distributions \\
						cond_dists, conditional distributions \\
						base, the base for entropy calculations \\
Returns:		Float64, the conditional entropy \\
"""

function condentropy(marg_dist::Array{Float64,1}, cond_dists::Array{Float64}; base::Int=2)
	length(marg_dist) == size(cond_dists)[2] || error("Number of marginal and conditional distributions does not match!")
	ents = [entropy(cond_dists[:, i], base) for i = 1:size(cond_dists)[2]]
	return sum(marg_dist .* ents)
end