text:
abstract type CSDistribution{T} end
abstract type Delta{T} <: CSDistribution{T} end

struct HardDelta{T<:Function} <: Delta{T}
    fact::T
end
struct SoftDelta{T<:Function} <: Delta{T}
    fact::T
    negative_evidence::Float64
end

Distributions.rand(d::Delta) = d.fact()
Distributions.logpdf(d::HardDelta{T}, x) where T = d.fact(x) ? 0.0 : -Inf
Distributions.logpdf(d::SoftDelta{T}, x) where T = d.fact(x) ? log(1.0 - exp(d.negative_evidence)) : d.negative_evidence

export HardDelta, SoftDelta

abstract type Program{T} <: CSDistribution{T} end
struct TypedSF{T} <: CSDistribution{T}
    sf::T
    input_type::Union{DataType,Union}
    output_type::Union{DataType,Union}
end
struct CompiledSF{T, S} <: CSDistribution{T}
    p::TypedSF{T}
    input::S
    addresses::Array{String}
    input_is_evidence::Bool
    function CompiledSF{T,S}(p::TypedSF{T}, input::S, addresses::Vector{String}) where {T,S}
        !(input isa p.input_type) && error("Wrong input type")
        input_usage = guess_input_is_evidence(p.sf, input)
        new{T,S}(p, input, addresses, input_usage)
    end
end
compile(p::TypedSF{T}, input::S, addresses::Vector{String}) where {T,S} = CompiledSF{T,S}(p, input, addresses)

(p::TypedSF)(x) = x isa p.input_type ? (y = p.sf(trace(), x); y isa p.output_type ? y : error("Wrong output type")) : error("Wrong input type")
(p::TypedSF)(x...) = x isa p.input_type ? (y = p.sf(trace(), x...); y isa p.output_type ? y : error("Wrong output type")) : error("Wrong input type")

function guess_input_is_evidence(f, i) :: Bool
    t = trace()
    f(t, i)
    for (address, node) in t.trace
        node.observed ? (return true) : nothing
    end
    false
end

function select_csf_sample(csf::CompiledSF, samples)
    if length(csf.addresses) == 1
        return samples[csf.addresses[1]][end]
    else
        return Dict(
            a => samples[a][end] for a in csf.addresses
        )
    end
end

function sample(csf::CompiledSF)
    params = csf.input === nothing ? () : (csf.input,)
    if csf.input_is_evidence
        mh(csf.p.sf; params = params)
    else
        forward_sampling(csf.p.sf; params = params)
    end
end

function Distributions.rand(csf::CompiledSF)
    samples = sample(csf)
    select_csf_sample(csf, samples)
end

function node(value, address :: A, dist :: CompiledSF, is_obs :: Bool, i :: Interpretation) where A
    T = typeof(value)
    samples = sample(dist)
    lp = logprob(samples, dist.addresses)
    ParametricNode{A, CompiledSF, T}(address, dist, value, lp, lp, is_obs, Array{Node, 1}(), Array{Node, 1}(), i, i)
end

export TypedSF, CompiledSF, compile


text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

include("normopt_initialization_template.jl")

const gurobi_env = Gurobi.Env()

const numVariables = 100
const numJCC = 100
const conRHS = (numVariables*1.0)^2

# compute probability that constraints aren't satisfied
function computeRiskLevel(x::Array{Float64})

	x_sqr = Array{Float64}(0)
	for j = 1:numVariables
		if(x[j] > 0)
			push!(x_sqr,(x[j])^2)
		end
	end

	n::Int64 = size(x_sqr,1)
	x_sqr = sort(x_sqr)

	if(x_sqr[n] > x_sqr[1])
		baseFactor1::Float64 = 1.0
		for i = 1:n-1
			baseFactor1 *= sqrt(x_sqr[n]/x_sqr[i])
		end

		k::Int64 = 1
		error_tol::Float64 = 1E-09/numJCC
		error::Float64 = 1.0
		while(error > error_tol)
			k *= 2
			baseFactor2 = exp(lgamma(n/2.0 + k) - lgamma(n/2.0))
			baseFactor3 = exp(k*log(x_sqr[n]/x_sqr[1] - 1.0) - k*log(2) - lfact(k))
			baseFactor4 = chisqcdf(n+2*k,conRHS/x_sqr[n])
			error = baseFactor1*baseFactor2*baseFactor3*baseFactor4
		end

		p::Float64 = 2.0/(1.0/x_sqr[1] + 1.0/x_sqr[n])
		g_base = 1.0 - p./x_sqr
		g = zeros(k-1)
		for i = 1:k-1
			g[i] = sum(g_base.^i)
		end

		c = zeros(Float64,k)
		c[1] = 1.0
		for j = 1:n
			c[1] *= sqrt(p/x_sqr[j])
		end
		for j = 1:k-1
			for r = 1:j
				c[j+1] += g[j+1-r]*c[r]
			end
			c[j+1] /= (2.0*j)
		end

		prob = 0.0
		for j = 0:k-1
			prob += c[j+1]*chisqcdf(n+2*j,conRHS/p)
		end

	else
		prob = chisqcdf(n,conRHS/x_sqr[1])
	end

	riskLevel::Float64 = 1.0 - (prob)^(numJCC)

	return riskLevel
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64},scalings::Array{Float64})

	con = (sum((xi.*x').^2,2) - conRHS*ones(numJCC))./scalings

	return con
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64})

	con = evaluateConstraint(x,xi,ones(numJCC))

	return con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraintGradients(x,xi,ones(numJCC))
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Array{Float64})

	grad_con = zeros(Float64,numJCC,numVariables)

	for i = 1:numJCC
		grad_con[i,:] = (2*(xi[i,:]).^2.*x)/scalings[i]
	end

	return grad_con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	grad_con = (2*(xi.^2).*x)/scalings

	return grad_con
end

# generate samples of the random variables
function generateRandomSamples(numSamples::Int64)

	xi = rand(Normal(0.0,1.0), numJCC, numVariables, numSamples)
	return xi
end

# project the current iterate onto the feasible set
function projectOntoFeasibleSet(x_curr::Array{Float64},objectiveBound::Float64,proj_params::Array{Float64})

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[j=1:numVariables] >= 0)

	@objective(mod, Min, sum((y[i]-x_curr[i])^2 for i = 1:numVariables))
	@constraint(mod, -sum(y[j] for j=1:numVariables) <= objectiveBound)

	solve(mod)
	x_proj = getvalue(y)

	return x_proj, proj_params
end

# get the stochastic gradient of the approximation
function getStochasticGradient(x::Array{Float64},numSamples::Int64,mu::Float64,tau::Float64,scalings::Array{Float64})

	xi = generateRandomSamples(numSamples)

	stochgrad = zeros(Float64,numVariables)

	for samp = 1:numSamples
		con = evaluateConstraint(x,xi[:,:,samp],scalings)

		max_info = findmax(con)
		max_con::Float64 = max_info[1]
		max_index::Int64 = max_info[2]

		baseFactor1::Float64 = tau/(sqrt(mu)*exp(0.5*tau*max_con) + sqrt(1.0/mu)*exp(-0.5*tau*max_con))^2

		grad_max_con = evaluateConstraintGradients(x,xi[max_index,:,samp],scalings[max_index])

		stochgrad += baseFactor1*grad_max_con/numSamples
	end

	return stochgrad

end

# get the stochastic gradient of the approximation for fixed random sample
function getStochasticGradientForFixedSample(x::Array{Float64},xi::SubArray{Float64,3,Array{Float64,3},Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true},mu::Float64,tau::Float64,scalings::Array{Float64})

	numSamples::Int64 = size(xi,3)

	stochgrad = zeros(Float64,numJCC,numVariables)

	for iter = 1:numSamples

		con = evaluateConstraint(x,xi[:,:,iter],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,:,iter],scalings)

		for i = 1:numJCC
			baseFactor1::Float64 = tau/(sqrt(mu)*exp(0.5*tau*con[i]) + sqrt(1.0/mu)*exp(-0.5*tau*con[i]))^2

			stochgrad[i,:] += baseFactor1*grad_con[i,:]/numSamples
		end
	end

	return stochgrad

end

# estimate constraint scalings
function estimateScalings(x::Array{Float64},numSamples::Int64)

	scaling_tol::Float64 = 1E-06

	xi = generateRandomSamples(numSamples)

	constraint_values = zeros(Float64,numJCC,numSamples)

	for samp = 1:numSamples
		constraint_values[:,samp] = evaluateConstraint(x,xi[:,:,samp])
	end

	scalings = zeros(Float64,numJCC)
	for index = 1:numJCC
		scalings[index] = max(median(abs.(constraint_values[index,:])),scaling_tol)
	end

	return scalings
end

# estimate Lipschitz constant of the gradient of the approximation
function estimateCompositeLipschitzConstant(x_ref::Array{Float64},numSamplesForLipVarEst::Int64,numGradientSamples::Int64,batchSize::Int64,mu::Float64,tau::Float64,scalings::Array{Float64},objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	xi = generateRandomSamples(batchSize)

	L_avg::Float64 = 0.0
	for batch = 1:batchSize
		xi_start::Int64 = batch
		xi_end::Int64 = batch

		L_max::Float64 = 0.0
		for samp = 1:numSamplesForLipVarEst
			x_1 = pickPointOnSphere(x_ref,sample_radius)
			x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

			x_2 = pickPointOnSphere(x_ref,sample_radius)
			x_2,~ = projectOntoFeasibleSet(x_2,objectiveBound,proj_params)

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,:,xi_start:xi_end],mu,tau,scalings)
			@views grad_2 = getStochasticGradientForFixedSample(x_2,xi[:,:,xi_start:xi_end],mu,tau,scalings)

			L_est::Float64 = norm(grad_2 - grad_1)/norm(x_2 - x_1)
			if(L_est > L_max)
				L_max = L_est
			end
		end

		L_avg += L_max/batchSize
	end

	return L_avg
end

# estimate step length factor
function estimateStepLength(x_ref::Array{Float64},numGradientSamples::Int64,
							numSamplesForLipVarEst::Int64,batchSize::Int64,maxNumIterations::Int64,
							minNumReplicates::Int64,mu::Float64,tau::Float64,
							scalings::Array{Float64},objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	rho_est::Float64 = estimateCompositeLipschitzConstant(x_ref,numSamplesForLipVarEst,numGradientSamples,batchSize,mu,tau,scalings,objectiveBound)

	variance_est::Float64 = 0.0
	for samp = 1:numSamplesForLipVarEst
		x_1 = pickPointOnSphere(x_ref,sample_radius)
		x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

		variance_1::Float64 = 0.0
		for batch = 1:batchSize
			grad_1 = getStochasticGradient(x_1,numGradientSamples,mu,tau,scalings)

			variance_1 += norm(grad_1)^2/batchSize
		end

		if(variance_1 > variance_est)
			variance_est = variance_1
		end
	end

	gamma::Float64 = 1/sqrt(rho_est*variance_est)

	stepLength::Float64 = gamma/sqrt((maxNumIterations+1.0)*minNumReplicates)

	return stepLength
end

# given a reference point, pick a random point within a given radius
function pickPointOnSphere(x_orig::Array{Float64},radius::Float64)

	if(radius < 1E-09)
		error("Sampling radius in pickPointOnSphere is almost zero! Set manually.")
	end

	n::Int64 = size(x_orig,1)
	unifRand = rand(1)
	factor1::Float64 = radius*(unifRand[1])^(1.0/n)
	normRand = rand(Normal(0,1), n)
	factor2::Float64 = norm(normRand)

	x_samp = x_orig + factor1*normRand/factor2

	return x_samp
end

ext: jl

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using LinearAlgebra
using CUDA

"""
	sim_activity_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32};
		dt=1f-1, t_bound=1f3, s=0f0, savestep=1, minimum_norm=0)

Integrates a tanh rate network with coupling matrix J and gain param g from
initial condition h0 and returns the matrix of single unit activities.
"""
function sim_activity_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32};
		dt=1f-1, t_bound=1f3, s=0f0, savestep=1, minimum_norm=0)

	N = size(J)[1]
	# Check dimensional consistency
	@assert(N==size(J)[2]==length(h0))

	maxstep = round(Int64, t_bound / dt)

	# Allocate memory for output
	out = CuArray{Float32,2}(undef, N, maxstep ÷ savestep)

	h = CuArray{Float32,1}(h0) # Set initial condition
	r = CuArray{Float32,1}(undef, N)

	J = CuArray{Float32,2}(g .* J + s * I)

	step = 0
	while step < maxstep
		step += 1
		map!(tanh, r, h)
		mul!(h, J, r, dt, 1 - dt)
		if step % savestep == 0 # Save every specified time interval
			savept = step ÷ savestep
			out[:,savept] .= h
		end
		if norm(h) < minimum_norm
			D = Matrix(-I + J .* ((sech.(h)).^2)')
			if all(real(eigvals(D)) .< 0)
				break
			end
		end
	end

	return Matrix{Float32}(out)
end

"""
	sim_lyapunov_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32},
		dh::Vector{Float32}, Q0::Matrix{Float32}; dt=1f-1, t_bound=1f3,
		t_ONS=1f0, wait=0f0, return_max_LE=false, s=0f0, printstuff=false)

Integrates a tanh network with a given weight matrix and gain parameter from a
specified initial condition, and returns a Lyapunov spectrum calculated from
Jacobians along the solution trajectory.

Can also return the max Lyapunov exponent calculated from a direct numerical
simulation of a perturbed trajectory.
"""
function sim_lyapunov_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32},
		dh::Vector{Float32}, Q0::Matrix{Float32}; dt=1f-1, t_bound=1f3,
		t_ONS=1f0, wait=0f0, return_max_LE=false, s=0f0, printstuff=false)

	# Check dimensional consistency
	N = size(J, 1)
	@assert size(J) == (length(h0), length(dh)) == size(Q0) == (N,N)

	# Discretization of time
	maxstep = round(Int64, t_bound / dt)
	startstep = round(Int64, wait / dt)
	tONSstep = round(Int64, t_ONS / dt)
	@assert(maxstep % tONSstep == 0 && startstep % tONSstep == 0,
			"The reorthonormalization interval must divide both the total " *
			"simulation time and wait period"
	)

	# Allocate memory for output
	spectrum = CUDA.zeros(N)

	# Allocate vector for activity and firing rate of current state
	h = CuArray{Float32,1}(h0)
	r = CuArray{Float32,1}(undef, N)

	# Load network structure and ICs to GPU
	J = CuArray{Float32,2}(g .* J + s * I)
	Q = CuArray{Float32,2}(Q0)

	# Allocate mem for calculating Jacobian and evolving the orthonormal basis
	sech2(x::Float32)::Float32 = 1f0 / (CUDA.cosh(x))^2
	sech2h = CuArray{Float32,2}(undef, 1, N)
	D = CuArray{Float32,2}(undef, N, N)
	dQ = CuArray{Float32,2}(undef, N, N)

	# Things for performing a direct numerical sim of a perturbation
	if return_max_LE
		max_LE = 0.0f0
		p_factor = 1f-6  # Order of the perturbation handled by the sim
		dh = CuArray{Float32,1}(dh ./ norm(dh))
		h2 = h .+ p_factor .* dh
		r2 = CuArray{Float32,1}(undef, N)
	end

	step = 0
	while step < maxstep
		step += 1

		#Evolve the trajectory and the orthonormal basis in the tangent space
		map!(tanh, r, h)
		mul!(h, J, r, dt, 1 - dt)

		map!(sech2, sech2h, h)
		copyto!(D, J)
		D .*= sech2h
		D[diagind(D)] .-= 1.0

		mul!(dQ, D, Q, dt, 0f0)
        Q .+= dQ

		#Evolve the perturbed trajectory and calculate how much it has grown
		if return_max_LE
			map!(tanh, r2, h2)
			mul!(h2, J, r2, dt, 1 - dt)
			copyto!(dh, h2)
			dh .-= h
			dh ./= p_factor
			scale = norm(dh)

			# If the perturbation scale has grown or shrunk too much, then
			# reset it to 1 and track the cumulative sum of log scales
			if scale > 100 || scale < 0.01 || step == startstep
				copyto!(h2, h)
				dh ./= (scale / p_factor)
				h2 .+= dh
				if step > startstep
					max_LE += log(scale)
				end
			end
		end

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(dQ, Q)
			F = qr!(dQ)
			fill!(Q, 0)
			diags = view(F.factors, diagind(F.factors))
			map!(sign, view(Q, diagind(Q)), diags)
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
            lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	actual_time = dt * (maxstep - startstep)
	if return_max_LE
		max_LE += log(norm((h2 .- h) ./ p_factor))
		spectrum ./= actual_time
		max_LE /= actual_time
		return (Vector{Float32}(spectrum), max_LE)
	else
		spectrum ./= actual_time
		return Vector{Float32}(spectrum)
	end
end


text:
function thermal_electron_electron_collsions()
end

function electron_ion_collisions()
end

"""
# Arguments
 - `μ`:
 - `μ′`:
 - `Z`:
 - `Z′`:
 - `Tᵢ`:
 - `Tᵢ′`:
 - `nᵢ `:
 - `nᵢ′`:
"""
function mixed_ion_collisions(μ, μ′, Z, Z′, Tᵢ, Tᵢ′, nᵢ, nᵢ′)
	T1,	T2 = ustrip.(u"eV", (Tᵢ, Tᵢ′))
	n1,	n2 = ustrip.(u"cm^-3", (nᵢ, nᵢ′))
    λᵢᵢ = 23 - log(
        ((Z * Z′ * (μ + μ′)) / (μ * T2 + μ′ * T1)) *
        √((n1 * Z^2 / T1) + (n2 * Z′^2 / T2)))

    return λᵢᵢ
end

function counter_streaming_ion_collisions(μ, μ′, Z, Z′, Tᵢ, Tᵢ′,nᵢ, nᵢ′)
end


text:
@blasfun clblasSaxpy(n::Csize_t, alpha::cl.CL_float,
                     X::cl.CL_mem, offx::Csize_t, incx::Cint,
                     Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                     n_queues::cl.CL_uint,
                     queues::Ptr{cl.CL_command_queue},
                     n_events_in_wait_list::cl.CL_uint,
                     event_wait_list::Ptr{cl.CL_event},
                     events::Ptr{cl.CL_event})

@blasfun2 clblasSaxpy(n::Csize_t, alpha::cl.CL_float,
                      X::cl.CL_mem, offx::Csize_t, incx::Cint,
                      Y::cl.CL_mem, offy::Csize_t, incy::Cint)

@blasfun clblasDaxpy(n::Csize_t, alpha::cl.CL_double,
                     X::cl.CL_mem, offx::Csize_t, incx::Cint,
                     Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                     n_queues::cl.CL_uint,
                     queues::Ptr{cl.CL_command_queue},
                     n_events_in_wait_list::cl.CL_uint,
                     event_wait_list::Ptr{cl.CL_event},
                     events::Ptr{cl.CL_event})

@blasfun2 clblasDaxpy(n::Csize_t, alpha::cl.CL_double,
                      X::cl.CL_mem, offx::Csize_t, incx::Cint,
                      Y::cl.CL_mem, offy::Csize_t, incy::Cint)

@blasfun clblasCaxpy(n::Csize_t, alpha::CL_float2,
                     X::cl.CL_mem, offx::Csize_t, incx::Cint,
                     Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                     n_queues::cl.CL_uint,
                     queues::Ptr{cl.CL_command_queue},
                     n_events_in_wait_list::cl.CL_uint,
                     event_wait_list::Ptr{cl.CL_event},
                     events::Ptr{cl.CL_event})

@blasfun2 clblasCaxpy(n::Csize_t, alpha::CL_float2,
                      X::cl.CL_mem, offx::Csize_t, incx::Cint,
                      Y::cl.CL_mem, offy::Csize_t, incy::Cint)

@blasfun clblasZaxpy(n::Csize_t, alpha::CL_double2,
                     X::cl.CL_mem, offx::Csize_t, incx::Cint,
                     Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                     n_queues::cl.CL_uint,
                     queues::Ptr{cl.CL_command_queue},
                     n_events_in_wait_list::cl.CL_uint,
                     event_wait_list::Ptr{cl.CL_event},
                     events::Ptr{cl.CL_event})

@blasfun2 clblasZaxpy(n::Csize_t, alpha::CL_double2,
                      X::cl.CL_mem, offx::Csize_t, incx::Cint,
                      Y::cl.CL_mem, offy::Csize_t, incy::Cint)


text:
type MSECriterion{A <: AbstractArray, B <: AbstractArray} <: Criterion{A, B}
    output :: B
    gradInput :: A
end

MSECriterion(T :: Type, batch_size, input_size) =
    MSECriterion(Array(T, (batch_size,)), Array(T, (input_size, batch_size)))

function updateOutput!{A,B}(self :: MSECriterion{A,B}, input :: A, target :: A)
    local input_size = size(self.gradInput)[1]
    local batch_size = size(self.gradInput)[2]
    local sum = zero(eltype(A))
    for b = 1:batch_size
        sum = 0
        for i = 1:input_size
            local r = input[i,b] - target[i,b]
            sum += r * r
        end
        sum /= input_size
        self.output[b] = sum
    end
    return self.output
end

function updateGradInput!{A,B}(self :: MSECriterion{A,B},
                               input :: A,
                               target :: A)
    local input_size = size(self.gradInput)[1]
    local batch_size = size(self.gradInput)[2]
    local sum = zero(eltype(A))
    local d = eltype(A)(2) / input_size
    for b = 1:batch_size
        sum = 0
        for i = 1:input_size
            self.gradInput[i,b] = (input[i,b] - target[i,b]) * d
        end
    end
    return self.gradInput
end

export MSECriterion


text:
import JSON

type NBodySystem
  N
  time
  bodies::Array{Body}

  initial_energy
  soften_len
end

NBodySystem(n, current_time=0) = NBodySystem(n, current_time, Array{Body}(n), 0, 0)

function NBodySystem(d::Dict{String, Any})
  bodies = map(Body, d["bodies"])
  NBodySystem(length(bodies), d["time"], bodies, get(d, "initial_energy", 0), get(d, "soften_len", 0))
end

function set_ids(nb::NBodySystem)
  id = 0
  for b in nb.bodies
    b.id = id
    id += 1
  end
end

function diff(a::NBodySystem, b::NBodySystem)
  res = NBodySystem(a.N, b.time)
  bodies = filter!(x->!isnull(x), map(a.bodies) do body
    other_body_inds = find(x-> x.id == body.id ,b.bodies)
    if length(other_body_inds) == 0
      return Nullable{Body}()
    end

    other_body = b.bodies[other_body_inds[1]]

    Body(
      mass = body.mass,
      pos = other_body.pos - body.pos,
      vel = other_body.vel - body.vel,
      id = body.id
    )
  end)

  res.bodies = bodies
  res
end

function abs(nb::NBodySystem)
  sqrt(
    mapreduce(
      x -> dot(x.pos, x.pos) + dot(x.vel, x.vel),
      +,
      nb.bodies))
end

function abs_pos(nb::NBodySystem)
  sqrt(
    mapreduce(
      x -> dot(x.pos, x.pos),
      +,
      nb.bodies))
end

####### Energy #########

function kin_energy(nb::NBodySystem)
  return mapreduce(kin_energy, +, nb.bodies)
end

function pot_energy(nb::NBodySystem)
  return mapreduce(x -> pot_energy(x, nb), +, nb.bodies)  / 2
end

function total_energy(nb::NBodySystem)
  return kin_energy(nb) + pot_energy(nb)
end

function init_energy!(nb::NBodySystem)
  nb.initial_energy = total_energy(nb)
end

function collision_time_scale(nb::NBodySystem)
  time_scale = Inf
  for b in nb.bodies
    indv_time_scale = collision_time_scale(b, nb.bodies)
    if indv_time_scale < time_scale
      time_scale = indv_time_scale
    end
  end
  return time_scale
end

function accel(body::Body, nb::NBodySystem)
  return accel(body, nb.bodies, nb.soften_len)
end

function jerk(body::Body, nb::NBodySystem)
  return jerk(body, nb.bodies, nb.soften_len)
end

function pot_energy(body::Body, nb::NBodySystem)
  return pot_energy(body, nb.bodies, nb.soften_len)
end

######## Utils #######

function ordered_radii(nb::NBodySystem)
  sort(
    map(x-> dot(x.pos, x.pos), nb.bodies)
  )
end

function quartiles(nb::NBodySystem)
  sorted = ordered_radii(nb)
  r1 = sqrt(sorted[ceil(Int, nb.N/4)])
  r2 = sqrt(sorted[ceil(Int, nb.N/2)])
  r3 = sqrt(sorted[ceil(Int, nb.N*3/4)])
  r4 = sqrt(sorted[end])
  (r1, r2, r3, r4)
end

function shift_to_center_of_mass!(nb::NBodySystem)
  mass_total = mapreduce(b -> b.mass, +, nb.bodies)
  com_pos = mapreduce(b -> b.pos*b.mass, +, nb.bodies) / mass_total
  com_vel = mapreduce(b -> b.vel*b.mass, +, nb.bodies) / mass_total

  for b in nb.bodies
    b.pos -= com_pos
    b.vel -= com_vel
  end
end

function adjust_units!(nb::NBodySystem)
  alpha = -pot_energy(nb) / 0.5
  beta = kin_energy(nb) / 0.25

  for b in nb.bodies
    b.pos *= alpha
    b.vel /= sqrt(beta)
  end
end

######## IO ##########

read_nbody() = read_nbody_json()

function read_nbody_json(stream::IO = STDIN)
  strs = readlines(stream)
  return read_nbody_json(join(strs))
end

function read_nbody_json(str::String)
  js = JSON.parse(str)
  s = NBody.NBodySystem(js)
  return s
end

function read_nbody_dml()
  n = parse(Int, readline())
  current_time = parse(Float64, readline())

  nb = NBodySystem(n, current_time)
  for i in 1:n
    body = read_body()
    nb.bodies[i] = body
  end
  return nb
end

function write_snapshot(nb::NBodySystem)
  write_snapshot(STDOUT, nb)
end

function write_snapshot(io::IO, nb::NBodySystem)
  println(io, JSON.json(nb))
end

function ppx(nb::NBodySystem, io = STDERR)
  # print(io, "N: ", length(nb.bodies), "\n")
  @printf(io, "time: %24.16e\n", nb.time)
  for a in nb.bodies
    show(io, a)
    acc = accel(a, nb)
    print(io, "accel: ", join(acc, ", "), "\n\n")
  end
end

# import Base.show
# function show(io::IO, nb::NBodySystem)
#   print(io, "N: ", length(nb.bodies), "\n")
#   @printf(io, "time: %24.16e\n", nb.time)
#   for a in nb.bodies
#     pp(io, a, nb)
#   end
# end

# import Base.print
# function print(io::IO, nb::NBodySystem)
#   print(io, length(nb.bodies), "\n")
#   @printf(io, "%24.16e\n", nb.time)
#   for a in nb.bodies
#     print(io, a)
#   end
# end

function write_stats(nb::NBodySystem, steps=0, x_info=false)
  tot_energy = total_energy(nb)

  s = """
        Time  $(@sprintf("%.3g", nb.time))
       Steps  $steps
  ========================
       E_Kin  $(@sprintf("%.3g", kin_energy(nb)))
       E_Pot  $(@sprintf("%.3g", pot_energy(nb)))
       E_Tot  $(@sprintf("%.3g", tot_energy))
       E Err  $(@sprintf("%.3g", (tot_energy - nb.initial_energy)))
    Err/Init  $(@sprintf("%.3g", (tot_energy - nb.initial_energy) / nb.initial_energy))

  """
  write(STDERR, s)

  if x_info
    ppx(nb)
  end
end


text:
"""
    getFCProd(agg, FCFact, aggIndices, vibindices)

Get product of Franck-Condon factors. This way the trace over bath will be faster

"""
function getFCProd(agg, FCFact, aggIndices, vibindices)
    elLen = length(agg.molecules)
    aggIndLen = length(aggIndices)
    vibLen = length(vibindices[2])
    elLen += 1
    FCProd = zeros(eltype(FCFact), aggIndLen, aggIndLen)
    # println("FCFact")
    # println(FCFact)

    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)

        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            elOrder2 = OpenQuantumSystems.elIndOrder(elind2)
            #=
            K1 = vibindices[elOrder1][1]
            K2 = vibindices[elOrder1][end]
            L1 = vibindices[elOrder2][1]
            L2 = vibindices[elOrder2][end]
            FCProd[I, J] = sum(FCFact[K1:K2, I] .* FCFact[J, L1:L2])
            =#

            for m = 1:vibLen
                # according to quantarhei, trace_over_vibrations()
                K = vibindices[1][m]
                L = vibindices[1][m]
                FCProd[I, J] += FCFact[K, I] * FCFact[J, L]
            end
        end
    end
    return FCProd
end

"""
    trace_bath(rho, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from `rho`

``\\rho_\\text{tr} = \\operatorname{tr}_B \\{\\rho\\} =
\\sum_{k} \\langle k \\vert \\left( \\sum_{ab} \\rho_{am, bn} \\vert am \\rangle \\langle bn \\vert \\right)\\vert k \\rangle``

"""
function trace_bath(rho::Array, agg, FCProd, aggIndices, vibindices)
    elLen = length(agg.molecules)
    aggIndLen = length(aggIndices)
    vibLen = length(vibindices[2])
    elLen += 1
    rho_traced = zeros(eltype(rho), elLen, elLen)

    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)

        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            elOrder2 = OpenQuantumSystems.elIndOrder(elind2)

            rho_traced[elOrder1, elOrder2] += rho[I, J] * FCProd[I, J]
        end
    end
    return rho_traced
end

function trace_bath(
    rho::T,
    agg,
    FCProd,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    rho_traced =
        trace_bath(rho.data, agg, FCProd, aggIndices, vibindices)
    basisLen = size(rho_traced, 1)
    basis = GenericBasis([basisLen])
    return DenseOperator(basis, basis, rho_traced)
end

"""
    trace_bath_slow(rho, agg, FCFact, aggIndices, vibindices)

Trace out bath degrees of freedom from `rho` without the product of Franck-Condon factors.

"""
function trace_bath_slow(
    rho::Array,
    agg,
    FCFact,
    aggIndices,
    vibindices
)
    elLen = length(agg.molecules)
    aggIndLen = length(aggIndices)
    vibLen = length(vibindices[2])
    elLen += 1
    rho_traced = zeros(eltype(rho), elLen, elLen)

    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)

        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            elOrder2 = OpenQuantumSystems.elIndOrder(elind2)

            for m = 1:vibLen
                # according to quantarhei, trace_over_vibrations()
                K = vibindices[1][m]
                L = vibindices[1][m]
                rho_traced[elOrder1, elOrder2] +=
                    FCFact[K, I] * rho[I, J] * FCFact[J, L]
            end
        end
    end
    return rho_traced
end

function trace_bath_slow(
    rho::T,
    agg,
    FCFact,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    rho_traced = trace_bath_slow(
        rho.data,
        agg,
        FCFact,
        aggIndices,
        vibindices
    )
    basisLen = size(rho_traced, 1)
    basis = GenericBasis([basisLen])
    return DenseOperator(basis, basis, rho_traced)
end

"""
    trace_bath(rho, a, b, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from `rho` without the product of Franck-Condon factors.
The trace will be done only on the Hilber space for electric bra part `a` and ket part `b`.
Input density matrix `rho` is for the whole Hilber space. This method returns number.

"""
function trace_bath(rho::Array, a, b, agg, FCProd, aggIndices, vibindices)
    aggIndLen = length(aggIndices)
    vibLen = length(vibindices[2])
    rho_traced = eltype(rho)(0)

    a_vibindices = vibindices[a]
    b_vibindices = vibindices[b]

    for I in a_vibindices
        for J in b_vibindices
            rho_traced += rho[I, J] * FCProd[I, J]
        end
    end
    return rho_traced
end

function trace_bath(
    rho::T,
    a,
    b,
    agg,
    FCProd,
    aggIndices,
    vibindices,
) where {B<:Basis,T<:Operator{B,B}}
    rho_traced = trace_bath(rho.data, a, b, agg, FCProd, aggIndices, vibindices)
    return rho_traced
end

"""
    trace_bath_part(rho, a, b, agg, FCProd, aggIndices, vibindices)

Trace out bath degrees of freedom from `rho` without the product of Franck-Condon factors.
The trace will be done only on the Hilber space for electric bra part `a` and ket part `b`.
Input density matrix `rho` is only for the subspace. This method returns number.

"""
function trace_bath_part(
    rho::Array,
    a,
    b,
    agg,
    FCProd,
    aggIndices,
    vibindices
)
    vibLen = length(vibindices[2])
    rho_traced = eltype(rho)(0)

    a_vibindices = vibindices[a]
    b_vibindices = vibindices[b]

    for a_vib = 1:vibLen
        I = a_vibindices[a_vib]
        for b_vib = 1:vibLen
            J = b_vibindices[b_vib]
            rho_traced += rho[a_vib, b_vib] * FCProd[I, J]
        end
    end
    return rho_traced
end

function trace_bath_part(
    rho::T,
    a,
    b,
    agg,
    FCProd,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    rho_traced = trace_bath(
        rho.data,
        a,
        b,
        agg,
        FCProd,
        aggIndices,
        vibindices
    )
    return rho_traced
end

"""
    get_rho_bath(rho, agg, FCProd, aggIndices, vibindices; justCopy=false)

This method will return the bath part of `rho` knowing the result of [`trace_bath`](@ref) defined as follows

`` \\rho_\\text{bath} = \\operatorname{tr}_S \\{\\rho\\} ``

`` \\rho_{\\text{bath}, ab} = \\rho_{ab} / \\langle a \\vert \\operatorname{tr}_B \\{ \\rho \\}\\vert b \\rangle``

"""
function get_rho_bath(
    rho::Array,
    agg,
    FCProd,
    aggIndices,
    vibindices;
    justCopy = false,
)
    rho_traced =
        trace_bath(rho, agg, FCProd, aggIndices, vibindices)
    vibLen = length(vibindices[end])
    aggIndLen = length(aggIndices)
    elLen = length(agg.molecules)
    elLen += 1
    rho_bath = zeros(eltype(rho), aggIndLen, aggIndLen)
    rho_bath_ref = zeros(eltype(rho), vibLen, vibLen)

    el1_p = 0
    el2_p = 0
    for el1 = 1:elLen, el2 = 1:elLen
        if abs(rho_traced[el1, el2]) != 0
            el1_p = el1
            el2_p = el2
            break
        end
    end
    vib11 = vibindices[el1_p][1]
    vib12 = vibindices[el1_p][end]
    vib21 = vibindices[el2_p][1]
    vib22 = vibindices[el2_p][end]
    rho_bath_ref[:, :] = rho[vib11:vib12, vib21:vib22] / rho_traced[el1_p, el2_p]
    for el1 = 1:elLen, el2 = 1:elLen
        vib11 = vibindices[el1][1]
        vib12 = vibindices[el1][end]
        vib21 = vibindices[el2][1]
        vib22 = vibindices[el2][end]
        if abs(rho_traced[el1, el2]) != 0
            rho_bath[vib11:vib12, vib21:vib22] =
                rho[vib11:vib12, vib21:vib22] / rho_traced[el1, el2]
        else
            if justCopy
                rho_bath[vib11:vib12, vib21:vib22] = rho[vib11:vib12, vib21:vib22]
            else
                rho_bath[vib11:vib12, vib21:vib22] = rho_bath_ref[:, :]
            end
        end
    end

    return rho_bath
end

function get_rho_bath(
    rho::T,
    agg,
    FCProd,
    aggIndices,
    vibindices;
    justCopy = false,
) where {B<:Basis,T<:Operator{B,B}}
    rho_data = get_rho_bath(
        rho.data,
        agg,
        FCProd,
        aggIndices,
        vibindices;
        justCopy = justCopy,
    )
    return DenseOperator(rho.basis_l, rho.basis_r, rho_data)
end

"""
    ad(rho_traced, rho_bath, agg, FCProd, aggIndices, vibindices)

This is the inverse operation to the trace over bath [`trace_bath`](@ref) and [`get_rho_bath`](@ref)
defined as follows

`` \\rho = \\operatorname{ad}\\{\\rho_\\text{tr}, \\rho_\\text{bath} \\} ``

"""
function ad(
    rho_traced::Array,
    rho_bath::Array,
    agg,
    FCProd,
    aggIndices,
    vibindices;
)
    elLen = length(agg.molecules)
    aggIndLen = length(aggIndices)
    vibLen = length(vibindices[2])
    elLen += 1
    W = zero(rho_bath)

    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)

        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            elOrder2 = OpenQuantumSystems.elIndOrder(elind2)

            W[I, J] = rho_traced[elOrder1, elOrder2] * rho_bath[I, J]
        end
    end
    return W
end

function ad(
    rho_traced::T,
    W_bath::Array,
    agg,
    FCFact,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    W = ad(
        rho_traced.data,
        W_bath,
        agg,
        FCFact,
        aggIndices,
        vibindices
    )
    basisLen = size(W, 1)
    basis = GenericBasis([basisLen])
    return DenseOperator(basis, basis, W)
end

function ad(
    rho_traced::Array,
    W_bath::T,
    agg,
    FCFact,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    W = ad(
        rho_traced,
        W_bath.data,
        agg,
        FCFact,
        aggIndices,
        vibindices
    )
    basisLen = size(W, 1)
    basis = GenericBasis([basisLen])
    return DenseOperator(basis, basis, W)
end

function ad(
    rho_traced::T,
    W_bath::T,
    agg,
    FCFact,
    aggIndices,
    vibindices
) where {B<:Basis,T<:Operator{B,B}}
    W = ad(
        rho_traced.data,
        W_bath.data,
        agg,
        FCFact,
        aggIndices,
        vibindices
    )
    basisLen = size(W, 1)
    basis = GenericBasis([basisLen])
    return DenseOperator(basis, basis, W)
end

"""
    correlation_function(t, rho0_bath, Ham_0, Ham_I, agg, FCProd, aggInds, vibindices)

Get time dependent correlation function for a specified time `t` using following definition

`` C(t) = \\operatorname{tr}_B \\{ \\hat{H}_I^{(I)}(t) \\hat{H}_I \\rho_\\text{bath}(0) \\} ``

"""
function correlation_function(
    t,
    rho0_bath,
    Ham_0,
    Ham_I,
    agg,
    FCProd,
    aggInds,
    vibindices
)
    Ham_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)
    prod = Ham_II_t * Ham_I * rho0_bath
    return trace_bath(prod, agg, FCProd, aggInds, vibindices)
end


text:
# ---
# title: 1360. Number of Days Between Two Dates
# id: problem1360
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/number-of-days-between-two-dates/description/>
# hidden: true
# ---
#
# Write a program to count the number of days between two dates.
#
# The two dates are given as strings, their format is `YYYY-MM-DD` as shown in
# the examples.
#
#
#
# **Example 1:**
#
#
#
#     Input: date1 = "2019-06-29", date2 = "2019-06-30"
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: date1 = "2020-01-15", date2 = "2019-12-31"
#     Output: 15
#
#
#
#
# **Constraints:**
#
#   * The given dates are valid dates between the years `1971` and `2100`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
type Sour{T<:Real}
    pH :: T
end

issour(A) = false
issour(A::Sour) = A.pH < 7


text:
# Binary heap

struct MutableBinaryHeapNode{T}
    value::T
    handle::Int
end

#################################################
#
#   core implementation
#
#################################################

function _heap_bubble_up!(ord::Ordering,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int}, nd_id::Int) where T

    @inbounds nd = nodes[nd_id]
    v::T = nd.value

    swapped = true  # whether swap happens at last step
    i = nd_id

    while swapped && i > 1  # nd is not root
        p = i >> 1
        @inbounds nd_p = nodes[p]

        if Base.lt(ord, v, nd_p.value)
            # move parent downward
            @inbounds nodes[i] = nd_p
            @inbounds nodemap[nd_p.handle] = i
            i = p
        else
            swapped = false
        end
    end

    if i != nd_id
        nodes[i] = nd
        nodemap[nd.handle] = i
    end
end

function _heap_bubble_down!(ord::Ordering,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int}, nd_id::Int) where T

    @inbounds nd = nodes[nd_id]
    v::T = nd.value

    n = length(nodes)
    last_parent = n >> 1

    swapped = true
    i = nd_id

    while swapped && i <= last_parent
        il = i << 1

        if il < n   # contains both left and right children
            ir = il + 1

            # determine the better child
            @inbounds nd_l = nodes[il]
            @inbounds nd_r = nodes[ir]

            if Base.lt(ord, nd_r.value, nd_l.value)
                # consider right child
                if Base.lt(ord, nd_r.value, v)
                    @inbounds nodes[i] = nd_r
                    @inbounds nodemap[nd_r.handle] = i
                    i = ir
                else
                    swapped = false
                end
            else
                # consider left child
                if Base.lt(ord, nd_l.value, v)
                    @inbounds nodes[i] = nd_l
                    @inbounds nodemap[nd_l.handle] = i
                    i = il
                else
                    swapped = false
                end
            end

        else  # contains only left child
            nd_l = nodes[il]
            if Base.lt(ord, nd_l.value, v)
                @inbounds nodes[i] = nd_l
                @inbounds nodemap[nd_l.handle] = i
                i = il
            else
                swapped = false
            end
        end
    end

    if i != nd_id
        @inbounds nodes[i] = nd
        @inbounds nodemap[nd.handle] = i
    end
end

function _binary_heap_pop!(ord::Ordering,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int}, nd_id::Int=1) where T

    # extract node
    rt = nodes[nd_id]
    v = rt.value
    @inbounds nodemap[rt.handle] = 0

    if length(nodes) == 1
        # clear
        empty!(nodes)
    else
        # move the last node to the position of the removed node
        @inbounds nodes[nd_id] = new_rt = nodes[end]
        pop!(nodes)
        @inbounds nodemap[new_rt.handle] = nd_id

        if length(nodes) > 1
            if Base.lt(ord, new_rt.value, v)
                _heap_bubble_up!(ord, nodes, nodemap, nd_id)
            else
                _heap_bubble_down!(ord, nodes, nodemap, nd_id)
            end
        end
    end
    return v
end

function _make_mutable_binary_heap(ord::Ordering, ty::Type{T}, values) where T
    # make a static binary index tree from a list of values

    n = length(values)
    nodes = Vector{MutableBinaryHeapNode{T}}(undef, n)
    nodemap = Vector{Int}(undef, n)

    i::Int = 0
    for v in values
        i += 1
        @inbounds nodes[i] = MutableBinaryHeapNode{T}(v, i)
        @inbounds nodemap[i] = i
    end

    for i = 1 : n
        _heap_bubble_up!(ord, nodes, nodemap, i)
    end
    return nodes, nodemap
end

#################################################
#
#   Binary Heap type and constructors
#
#################################################

mutable struct MutableBinaryHeap{T, O <: Base.Ordering} <: AbstractMutableHeap{T, Int}
    ordering::O
    nodes::Vector{MutableBinaryHeapNode{T}}
    node_map::Vector{Int}

    function MutableBinaryHeap{T, O}() where {T, O}
        ordering = O()
        nodes = Vector{MutableBinaryHeapNode{T}}()
        node_map = Vector{Int}()
        new{T, O}(ordering, nodes, node_map)
    end

    function MutableBinaryHeap{T, O}(xs::AbstractVector{T}) where {T, O}
        ordering = O()
        nodes, node_map = _make_mutable_binary_heap(ordering, T, xs)
        new{T, O}(ordering, nodes, node_map)
    end
end

const MutableBinaryMinHeap{T} = MutableBinaryHeap{T, Base.ForwardOrdering}
const MutableBinaryMaxHeap{T} = MutableBinaryHeap{T, Base.ReverseOrdering}

MutableBinaryMinHeap(xs::AbstractVector{T}) where T = MutableBinaryMinHeap{T}(xs)
MutableBinaryMaxHeap(xs::AbstractVector{T}) where T = MutableBinaryMaxHeap{T}(xs)

function Base.show(io::IO, h::MutableBinaryHeap)
    print(io, "MutableBinaryHeap(")
    nodes = h.nodes
    n = length(nodes)
    if n > 0
        print(io, string(nodes[1].value))
        for i = 2 : n
            print(io, ", $(nodes[i].value)")
        end
    end
    print(io, ")")
end

#################################################
#
#   interfaces
#
#################################################

Base.length(h::MutableBinaryHeap) = length(h.nodes)

Base.isempty(h::MutableBinaryHeap) = isempty(h.nodes)

function Base.push!(h::MutableBinaryHeap{T}, v) where T
    nodes = h.nodes
    nodemap = h.node_map
    i = length(nodemap) + 1
    nd_id = length(nodes) + 1
    push!(nodes, MutableBinaryHeapNode(convert(T, v), i))
    push!(nodemap, nd_id)
    _heap_bubble_up!(h.ordering, nodes, nodemap, nd_id)
    return i
end

function Base.sizehint!(h::MutableBinaryHeap, s::Integer)
    sizehint!(h.nodes, s)
    sizehint!(h.node_map, s)
    return h
end

@inline Base.first(h::MutableBinaryHeap) = h.nodes[1].value

"""
    top_with_handle(h::MutableBinaryHeap)

Returns the element at the top of the heap `h` and its handle.
"""
function top_with_handle(h::MutableBinaryHeap)
    el = h.nodes[1]
    return el.value, el.handle
end

Base.pop!(h::MutableBinaryHeap{T}) where {T} = _binary_heap_pop!(h.ordering, h.nodes, h.node_map)

"""
    update!{T}(h::MutableBinaryHeap{T}, i::Int, v::T)

Replace the element at index `i` in heap `h` with `v`.
This is equivalent to `h[i]=v`.
"""
function update!(h::MutableBinaryHeap{T}, i::Int, v) where T
    nodes = h.nodes
    nodemap = h.node_map
    ordering = h.ordering

    nd_id = nodemap[i]
    v0 = nodes[nd_id].value
    x = convert(T, v)
    nodes[nd_id] = MutableBinaryHeapNode(x, i)
    if Base.lt(ordering, x, v0)
        _heap_bubble_up!(ordering, nodes, nodemap, nd_id)
    else
        _heap_bubble_down!(ordering, nodes, nodemap, nd_id)
    end
end

"""
    delete!{T}(h::MutableBinaryHeap{T}, i::Int)

Deletes the element with handle `i` from heap `h` .
"""
function Base.delete!(h::MutableBinaryHeap{T}, i::Int) where T
     nd_id = h.node_map[i]
    _binary_heap_pop!(h.ordering, h.nodes, h.node_map, nd_id)
    return h
end

Base.setindex!(h::MutableBinaryHeap, v, i::Int) = update!(h, i, v)
Base.getindex(h::MutableBinaryHeap, i::Int) = h.nodes[h.node_map[i]].value


text:
# This file is a part of SimilaritySearch.jl

export compute_distances, select_sss, select_tournament
using Random
"""
    compute_distances(dist::PreMetric, db::AbstractVector{T}, refs::AbstractVector{Int}, q::T) where T

Computes the distances of `q` to the set of references `refs` (each index point to an item in `db`)
It returns an array of tuples `(distance, refID)`
"""
function compute_distances(dist::PreMetric, db::AbstractVector{T}, refs::AbstractVector{Int}, q::T) where T
    [(evaluate(dist, q, db[refs[refID]]), refID) for refID in 1:length(refs)]
end

"""
    compute_distances(dist::PreMetric, refs::AbstractVector{T}, q::T) where T

Computes the distances of `q` to the set of references `refs`
It returns an array of tuples `(distance, refID)`
"""
function compute_distances(dist::PreMetric, refs::AbstractVector{T}, q::T) where T
    [(evaluate(dist, q, ref), refID) for (refID, ref) in enumerate(refs)]
end

"""
    select_tournament(dist::PreMetric, db::AbstractVector{T}, numrefs::Int, tournamentsize::Int) where T

select_tournament selects `numrefs` references from `db` using a tournament criterion; each
individual is selected among`tournamentsize` individuals.

It returns a set of pivots as a list of integers pointing to elements in `db`
"""
function select_tournament(dist::PreMetric, db::AbstractVector{T}, numrefs::Int, tournamentsize::Int) where T
    @info "select_tournament: db=$(typeof(db)), numrefs=$(numrefs), distance=$(dist), tournamentsize=$(tournamentsize)"

    refs = Vector{Int}()
    perm = 1:length(db) |> collect
    shuffle!(perm)
    push!(refs, perm[end])

    for i=1:(numrefs-1)
        sample = perm[(i-1)*tournamentsize+1:i*tournamentsize]
        distant = 0
        distantRef = 0
        for x in sample
            M = compute_distances(dist, db, refs, db[x])
            d, refID = minimum(M)
            if d > distant
                distant = d
                distantRef = x
            end
        end
        # println(sample, distant, distantRef)
        if distantRef > 0
            push!(refs, distantRef)
        end
    end

    return refs
end

"""
    select_sss(dist::PreMetric, db::AbstractVector{T}, alpha::Float64; shuf::Bool=true) where T

select_sss selects the necessary pivots to fulfill the SSS criterion using :param:alpha.
If :param:shuffle_db is true then the database is shuffled before the selection process; in any case,
the estimation of the maximum distance introduces indeterminism, however it could be too small.
If you need better a better random selection set :param:shuf as true

It returns a set of pivots as a list of integers pointing to elements in :param:db
"""
function select_sss(dist::PreMetric, db::AbstractVector{T}, alpha::Float64; shuf::Bool=true) where T
    @info "select_sss: db=$(typeof(db)), alpha=$(alpha), distance=$(dist), shuf=$(shuf)"
    dmax::Float64 = 0.0
    s = Int(round(sqrt(length(db))/2))
    sample1, sample2 = rand(db, s), rand(db, s)
    dmax = maximum([evaluate(dist, sample1[i], sample2[i]) for i=1:s])
    # plot(x=X, Geom.histogram(bincount=100))

    @info "the maximum distance estimated as $(dmax), now selecting pivots"
    xdb = Array(1:length(db))
    if shuf
        shuffle!(xdb)
    end
    #pivots::Array{T,1} = Array(T, 1)
    #pivots[1] = xdb[1]
    pivots = Int[xdb[1]]

    for i=2:length(xdb)
        obj = db[xdb[i]]
        if minimum([evaluate(dist, db[pivID], obj) for pivID in pivots]) / dmax < alpha
            continue
        end

        push!(pivots, i)
    end

    pivots
end


text:
#DS1000Z scope driver

mutable struct DS1000Z
    handle
	initialized
    instr_dict
end
#"MEAN" "PERIod" "PK2pk" "MINImum" "MAXImum"
DS1000Z(handle) = DS1000Z(
                handle,
				false,
                Dict(
                    "CH1"=>"1",
                    "CH2"=>"2",
					"CH3"=>"3",
                    "CH4"=>"4",
					"Source"=>":TRIGger:EDGe:SOURce",
                    "Level"=>":TRIGger:EDGe:LEVel",
                    "Mode"=>":TRIGger:SWEep",
					"AUTO"=>"AUTO",
					"NORMAL"=>"NORMAL",
                    "MEAS1"=>1,
                    "MEAS2"=>2,
                    "MEAS3"=>3,
                    "MEAS4"=>4,
                    "MEAS5"=>5,
		"MEAN"=>"VAVG",
		"PERIod"=>"PERIod",
		"PK2pk"=>"VPP",
		"MINImum"=>"VMIN",
		"MAXImum"=>"VMAX",
		"CRMs" => "VRMS",
		"FREQuency" => "FREQuency",
		"RISe" => "RTIMe",
		"FALL" => "FTIMe",
		"PWIdth" => "PWIDth|",
		"NWIdth" => "NWIDth|"
                    )
                    )

#Misc
function connect!(obj::DS1000Z, handle)
    obj.handle = handle
    obj.connected = 1
    return nothing
    end
function disconnect!(obj::DS1000Z)
    obj.handle = 0
    obj.connected = 0
    return nothing
    end

reset_instr(obj::DS1000Z) = write(obj.handle, "*RST")
clear_instr(obj::DS1000Z) = write(obj.handle, "*CLS")
get_idn(obj::DS1000Z) = query(obj.handle, "*IDN?")

# VERTICAL
function set_ch_coupling(obj::DS1000Z, ch, cpl)
	ch = obj.instr_dict[ch]
	write(obj.handle, ":CHANnel$ch:COUPling $cpl")
end
function get_ch_coupling(obj::DS1000Z, ch)
	ch = obj.instr_dict[ch]
	strip(query(obj.handle, ":CHANnel$ch:COUPling?"))
end

function set_vertical_scale(obj::DS1000Z, ch, value)
    ch = obj.instr_dict[ch]
    cmd = ":CHANnel$ch:SCAle $value"
    write(obj.handle, cmd)
end
function get_vertical_scale(obj::DS1000Z, ch)
    ch = obj.instr_dict[ch]
    cmd = ":CHANnel$ch:SCAle?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

function set_ch_position(obj::DS1000Z, ch, lev)
	ch = obj.instr_dict[ch]
    cmd = ":CHANnel$ch:OFFSet $lev"
    value = write(obj.handle, cmd)
end
function get_ch_position(obj::DS1000Z, ch)
	ch = obj.instr_dict[ch]
    cmd = ":CHANnel$ch:OFFSet?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

# HORIZONTAL
function set_horizontal_scale(obj::DS1000Z, value)
    cmd = ":TIMebase:MAIn:SCAle $value"
    write(obj.handle, cmd)
end
function get_horizontal_scale(obj::DS1000Z, time)
    #time = obj.instr_dict[time]
    cmd = ":TIMebase:MAIn:SCAle?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

# TRIGGER
function set_trig_ch(obj::DS1000Z, ch)
    ch = obj.instr_dict[ch]
    cmd = ":TRIGger:EDGe:SOURce CHANnel$ch"
    write(obj.handle, cmd)
end

function set_trig_level(obj::DS1000Z, level)
    cmd = ":TRIGger:EDGe:LEVel $level"
    write(obj.handle, cmd)
end

function set_trig_mode(obj::DS1000Z, mode)
	mode = obj.instr_dict[mode]
    cmd = ":TRIGger:SWEep $mode"            # AUTO or NORMal or SINGLE
    write(obj.handle, cmd)
end

function get_trig_data(obj::DS1000Z, trig) #???
    trig = obj.instr_dict[trig]
    cmd = "$trig?"
    value = strip(query(obj.handle, cmd))
end

# MEASUREMENT
function set_meas_ch(obj::DS1000Z, Meas_Nr1, ch)
    obj.instr_dict[Meas_Nr1] = obj.instr_dict[ch]
end

function set_meas_type(obj::DS1000Z, Meas_Nr1, Type)
    ch = obj.instr_dict[Meas_Nr1]
	Type = obj.instr_dict[Type]
	obj.instr_dict[Meas_Nr1] = obj.instr_dict[Meas_Nr1] * Type
	cmd = ":MEASure:ITEM $Type,CHANnel$ch"
    write(obj.handle, cmd)
end

function set_meas(obj::DS1000Z, Meas_Nr1, ch, Type)
    set_meas_ch(obj,Meas_Nr1, ch)
    set_meas_type(obj, Meas_Nr1, Type)
end
function get_meas_data(obj::DS1000Z, Meas_Nr1)
    nr = obj.instr_dict[Meas_Nr1]
    #cmd = "*OPC?"
    #value = strip(query(obj.handle, cmd))
	Type = obj.instr_dict[Meas_Nr1][2:end]
	ch = obj.instr_dict[Meas_Nr1][1]
    cmd = ":MEASure:ITEM? $Type,CHANnel$ch"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

# ACQUISITION
# select channel for acq and config
function conf_acq_ch(obj::DS1000Z, ch)
	ch = obj.instr_dict[ch]
    write(obj.handle, ":WAVeform:SOURce CHANnel$ch")
	write(obj.handle, ":WAVeform:MODE NORMal")
	write(obj.handle, ":WAVeform:FORMat BYTE") #WORD|BYTE|ASCii} 225
	write(obj.handle, ":WAVeform:DATA?")
end

function auto_set(obj::DS1000Z)
	write(obj.handle, ":AUToscale")
	@info "Scope auto-set in progress"
	sleep(1)
end

function Trigger_Aquistion(obj::DS1000Z, ch)
	ch = obj.instr_dict[ch]
	TMCdatadescription = zeros(UInt8, 11)
	viRead!(obj.handle, TMCdatadescription)
	effectivebytes = parse(Int32, string(Char.(TMCdatadescription[8:11])...))

	y_buffer = zeros(UInt8, effectivebytes)
	viRead!(obj.handle, y_buffer)

	YINCrement = parse(Float64, query(obj.handle, ":WAVeform:YINCrement?"))
	YORigin = parse(Float64, query(obj.handle, ":WAVeform:YORigin?"))
	YREFerence = parse(Float64, query(obj.handle, ":WAVeform:YREFerence?"))

	dt = parse(Float64, query(obj.handle, ":WAVeform:XINCrement?"))

	y = (y_buffer .- YORigin .- YREFerence) .* YINCrement
	t = collect(0:dt:dt*(length(y)-1))

	return t, y
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
using Test
using BlackBoxOptimizationBenchmarking, Optim
import BlackBoxOptimizationBenchmarking: minimizer, minimum, optimize
const BBOB = BlackBoxOptimizationBenchmarking

import Base.string

#interface for Optim
pinit(D) = 10*rand(D).-5
optimize(opt::Optim.AbstractOptimizer, f, D, run_length) =
    Optim.optimize(f, pinit(D), opt, Optim.Options(f_calls_limit=run_length, g_tol=1e-12))

string(opt::Optim.AbstractOptimizer) = string(typeof(opt).name.name)

test_functions = BBOB.list_functions()

b = BBOB.benchmark(
    NelderMead(), BBOB.sphere, [100, 500, 1000], 20, 2, 1e-6,
)

b = BBOB.benchmark(
    NelderMead(), test_functions[1], [100, 500, 1000], 20, 2, 1e-6,
)

b = BBOB.benchmark(
    GradientDescent(), test_functions[1:2], [100, 500, 1000], 20, 2, 1e-6,
)


text:
# Gradient regression
function get_gain_gpu(::L, ∑::AbstractVector{T}, λ::T) where {L <: GradientRegression,T <: AbstractFloat}
    gain = ∑[1]^2 / (∑[2] + λ * ∑[3]) / 2
    return gain
end

# Gaussian regression
function get_gain_gpu(::L, ∑::AbstractVector{T}, λ::T) where {L <: GaussianRegression,T <: AbstractFloat}
    gain = ∑[1]^2 / (∑[3] + λ * ∑[5]) / 2 + ∑[2]^2 / (∑[4] + λ * ∑[5]) / 2
    return gain
end

#####################
# linear
#####################
function kernel_linear_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(y)
        @inbounds δ[1,i] = 2 * (p[i] - y[i]) * δ[3,i]
        @inbounds δ[2,i] =  2 * δ[3,i]
    end
    return
end

function update_grads_gpu!(loss::Linear, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, (length(y)) / threads)
    @cuda blocks = blocks threads = threads kernel_linear_δ!(δ, p, y)
    CUDA.synchronize()
    return
end

#####################
# Logistic
#####################
function kernel_logistic_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(y)
        @inbounds pred = sigmoid(p[1,i])
        @inbounds δ[1,i] = (pred - y[i]) * δ[3,i]
        @inbounds δ[2,i] = pred * (1 - pred) * δ[3,i]
    end
    return
end

function update_grads_gpu!(loss::Logistic, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, (length(y)) / threads)
    @cuda blocks = blocks threads = threads kernel_logistic_δ!(δ, p, y)
    CUDA.synchronize()
    return
end

################################################################################
# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ)
################################################################################
function kernel_gauss_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    @inbounds if i <= length(y)
        # first order gradients
        δ[1,i] = (p[1,i] - y[i]) / exp(2 * p[2,i]) *  δ[5,i]
        δ[2,i] = (1 - (p[1,i] - y[i])^2 / exp(2 * p[2,i])) *  δ[5,i]
        # second order gradients
        δ[3,i] =  δ[5,i] / exp(2 * p[2,i])
        δ[4,i] = 2 *  δ[5,i] / exp(2 * p[2,i]) * (p[1,i] - y[i])^2
    end
    return
end

function update_grads_gpu!(loss::Gaussian, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, (length(y)) / threads)
    @cuda blocks = blocks threads = threads kernel_gauss_δ!(δ, p, y)
    CUDA.synchronize()
    return
end

function update_childs_∑_gpu!(::L, nodes, n, bin, feat) where {L}
    nodes[n << 1].∑ .= nodes[n].hL[:, bin, feat]
    nodes[n << 1 + 1].∑ .= nodes[n].hR[:, bin, feat]
    return nothing
end


text:
module M0_0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    9.00239411475607e-11*l1.^2 - 1.10330950130372e-10*l1.*l2 - 7.07348217074767e-11*l1.*l3 + 1.54854445732105e-13*l1 + 3.38058845865336e-11*l2.^2 + 4.33428102083718e-11*l2.*l3 - 9.49119158842984e-14*l2 + 1.38960161444059e-11*l3.^2 - 6.08170432755409e-14*l3 + 6.666876091763e-17
end
end


text:
using PlotlyLight
using JSON3
using Test

@testset "sanity check" begin
    @test Plot(Config(x = 1:10)) isa Plot
    @test Plot(Config(x = 1:10), Config(title="Title")) isa Plot
    @test Plot(Config(x = 1:10), Config(title="Title"), Config(displaylogo=true)) isa Plot
end
@testset "src" begin
    p = Plot(Config(y=1:10))

    PlotlyLight.src!(:cdn)
    @test occursin("plotly-latest.min.js", repr("text/html", p))

    PlotlyLight.src!(:none)
    @test !occursin("plotly-latest.min.js", repr("text/html", p))

    PlotlyLight.src!(:standalone)
    @test length(repr("text/html", p)) > 1000

    PlotlyLight.src!(:local)
    @test occursin("deps", repr("text/html", p))
end


text:
export  skipv,
        advance!,
        locate_keyword!,
        chomp_space!,
chomp_eol!,
_peekb

import Base: peek

@inline function chomp_space!(ps::IO)
    n = 0
    while !eof(ps) && ps |> _peekb |> ispdfspace
        skip(ps, 1)
        n += 1
    end
    return n
end

@inline chomp_eol!(ps::IO) =
    while !eof(ps) && ps |> _peekb |> is_crorlf skip(ps,1) end

@inline function skipv(ps::IO, c::UInt8)
    ch = 0xff
    !eof(ps) && ((ch = ps |> _peekb) == c) && return skip(ps, 1)
    error("Found '$(Char(ch))($(UInt(ch)))' Expected '$(Char(c))' here")
end

_peekb(io::IO) = UInt8(peek(io))

@inline skipv(ps::IO, cs::Vector{UInt8}) = for c in cs skipv(ps, c) end

@inline advance!(ps::IO) = read(ps, UInt8)

function kmp_preprocess(P)
    m = length(P)
    pi = zeros(Int, m)
    pi[1] = 0
    k = 0
    for q = 2:m
        while k > 0 && P[k+1] != P[q]
            k = pi[k]
        end
        if P[k+1] == P[q]
            k += 1
        end
        pi[q] = k
    end
    return pi
end

# This cannot be called on a marked stream.
function locate_keyword!(ps::IO, keyword, maxoffset=length(keyword))
    m = length(keyword)
    pi = kmp_preprocess(keyword)

    q = 0
    found=false
    offset = 0
    @assert !ismarked(ps)
    mark(ps)
    while(true)
        c = advance!(ps)
        offset += 1
        while q > 0 && keyword[q+1] != c
            q = pi[q]
        end
        if keyword[q+1] == c
            q +=1
        end
        if q == m
            found = true
            q = pi[q]
            break
        end
        if (offset >= maxoffset)
            break
        end
    end
    if found
        unmark(ps)
        return (offset - length(keyword))
    else
        reset(ps)
        peek(ps)
        return -1
    end
end


text:
@testset "cell linked list range test" begin
  range = 1000
  dataset = range * rand(2, 100) - range / 2
  cll = CellLinkedList(dataset, 3.)

  teststatus = true
  for c in 1:100
    pt = 1.1 * range * rand(2) - range * 1.1 / 2
    r = range / 3
    cllres = cll_inrange(cll, pt, r)
    naiveres = [p for p in 1:size(dataset, 2) if norm(pt - dataset[:, p]) < r + 1e-9]
    teststatus &= sort(cllres) == sort(naiveres)
    if !teststatus
      @show dataset
      @show pt
      @show r
      break
    end
  end
  @test teststatus == true
end


text:
function raindrops(number::Int)

end


text:
#
# Abstract Interface.
#

"""
Abbreviation objects are used to automatically generate context-dependant markdown content
within documentation strings. Objects of this type interpolated into docstrings will be
expanded automatically before parsing the text to markdown.

$(:FIELDS)
"""
abstract Abbreviation

"""
$(:SIGNATURES)

Expand the [`Abbreviation`](@ref) `abbr` in the context of the `DocStr` `doc` and write
the resulting markdown-formatted text to the `IOBuffer` `buf`.
"""
format(abbr, buf, doc) = error("`format` not implemented for `$typeof(abbr)`.")

# Only extend `formatdoc` once with our abstract type. Within the package use a different
# `format` function instead to keep things decoupled from `Base` as much as possible.
Docs.formatdoc(buf::IOBuffer, doc::Docs.DocStr, part::Abbreviation) = format(part, buf, doc)

#
# Implementations.
#

#
# `TypeFields`
#

"""
The singleton type for [`FIELDS`](@ref) abbreviations.

$(:FIELDS)
"""
immutable TypeFields <: Abbreviation end

"""
An [`Abbreviation`](@ref) to include the names of the fields of a type as well as any
documentation that may be attached to the fields.

# Examples

The generated markdown text should look similar to to following example where a
type has three fields (`x`, `y`, and `z`) and the last two have documentation
attached.

```markdown

  - `x`

  - `y`

    Unlike the `x` field this field has been documented.

  - `z`

    Another documented field.
```
"""
const FIELDS = TypeFields()

function format(::TypeFields, buf, doc)
    local docs = get(doc.data, :fields, Dict())
    local binding = doc.data[:binding]
    local object = Docs.resolve(binding)
    local fields = fieldnames(object)
    if !isempty(fields)
        println(buf)
        for field in fields
            print(buf, "  - `", field, "`\n")
            # Print the field docs if they exist and aren't a `doc"..."` docstring.
            if haskey(docs, field) && isa(docs[field], AbstractString)
                println(buf)
                for line in split(docs[field], "\n")
                    println(buf, isempty(line) ? "" : "    ", rstrip(line))
                end
            end
            println(buf)
        end
        println(buf)
    end
    return nothing
end

#
# `ModuleExports`
#

"""
The singleton type for [`EXPORTS`](@ref) abbreviations.

$(:FIELDS)
"""
immutable ModuleExports <: Abbreviation end

"""
An [`Abbreviation`](@ref) to include all the exported names of a module is a sorted list of
`Documenter.jl`-style `@ref` links.

!!! note

    The names are sorted alphabetically and ignore leading `@` characters so that macros are
    *not* sorted before other names.

# Examples

The markdown text generated by the `EXPORTS` abbreviation looks similar to the following:

```markdown

  - [`bar`](@ref)
  - [`@baz`](@ref)
  - [`foo`](@ref)

```
"""
const EXPORTS = ModuleExports()

function format(::ModuleExports, buf, doc)
    local binding = doc.data[:binding]
    local object = Docs.resolve(binding)
    local exports = names(object)
    if !isempty(exports)
        println(buf)
        # Sorting ignores the `@` in macro names and sorts them in with others.
        for name in sort(exports, by = s -> lstrip(string(s), '@'))
            # Skip the module itself, since that's always exported.
            name === module_name(object) && continue
            # We print linked names using Documenter.jl cross-reference syntax
            # for ease of integration with that package.
            println(buf, "  - [`", name, "`](@ref)")
        end
        println(buf)
    end
    return nothing
end

#
# `ModuleImports`
#

"""
The singleton type for [`IMPORTS`](@ref) abbreviations.

$(:FIELDS)
"""
immutable ModuleImports <: Abbreviation end

"""
An [`Abbreviation`](@ref) to include all the imported modules in a sorted list.

# Examples

The markdown text generated by the `IMPORTS` abbreviation looks similar to the following:

```markdown

  - `Foo`
  - `Bar`
  - `Baz`

```
"""
const IMPORTS = ModuleImports()

function format(::ModuleImports, buf, doc)
    local binding = doc.data[:binding]
    local object = Docs.resolve(binding)
    local imports = unique(ccall(:jl_module_usings, Any, (Any,), object))
    if !isempty(imports)
        println(buf)
        for mod in sort(imports, by = string)
            println(buf, "  - `", mod, "`")
        end
        println(buf)
    end
end

#
# `MethodList`
#

"""
The singleton type for [`METHODLIST`](@ref) abbreviations.

$(:FIELDS)
"""
immutable MethodList <: Abbreviation end

"""
An [`Abbreviation`](@ref) for including a list of all the methods that match a documented
`Method`, `Function`, or `DataType` within the current module.

# Examples

The generated markdown text will look similar to the following example where a function
`f` defines two different methods (one that takes a number, and the other a string):

````markdown
```julia
f(num)
```

defined at [`<path>:<line>`](<github-url>).

```julia
f(str)
```

defined at [`<path>:<line>`](<github-url>).
````
"""
const METHODLIST = MethodList()

function format(::MethodList, buf, doc)
    local binding = doc.data[:binding]
    local typesig = doc.data[:typesig]
    local modname = doc.data[:module]
    local func = Docs.resolve(binding)
    local groups = methodgroups(func, typesig, modname; exact = false)
    if !isempty(groups)
        println(buf)
        for group in groups
            println(buf, "```julia")
            for method in group
                printmethod(buf, binding, func, method)
                println(buf)
            end
            println(buf, "```\n")
            if !isempty(group)
                local method = group[1]
                local file = string(method.file)
                local line = method.line
                local path = cleanpath(file)
                local URL = url(method)
                isempty(URL) || println(buf, "defined at [`$path:$line`]($URL).")
            end
            println(buf)
        end
        println(buf)
    end
    return nothing
end

#
# `MethodSignatures`
#

"""
The singleton type for [`SIGNATURES`](@ref) abbreviations.

$(:FIELDS)
"""
immutable MethodSignatures <: Abbreviation end

"""
An [`Abbreviation`](@ref) for including a simplified representation of all the method
signatures that match the given docstring. See [`printmethod`](@ref) for details on
the simplifications that are applied.

# Examples

The generated markdown text will look similar to the following example where a function `f`
defines method taking two positional arguments, `x` and `y`, and two keywords, `a` and the `b`.

````markdown
```julia
f(x, y; a, b...)
```
````
"""
const SIGNATURES = MethodSignatures()

function format(::MethodSignatures, buf, doc)
    local binding = doc.data[:binding]
    local typesig = doc.data[:typesig]
    local modname = doc.data[:module]
    local func = Docs.resolve(binding)
    local groups = methodgroups(func, typesig, modname)
    if !isempty(groups)
        println(buf)
        println(buf, "```julia")
        for group in groups
            for method in group
                printmethod(buf, binding, func, method)
                println(buf)
            end
        end
        println(buf, "\n```\n")
    end
end

#
# `TypeSignature`
#

"""
The singleton type for [`TYPEDEF`](@ref) abbreviations.
"""
immutable TypeDefinition <: Abbreviation end

"""
An [`Abbreviation`](@ref) for including a summary of the signature of a type definition.
Some of the following information may be included in the output:

  * whether the object is an `abstract` type or a `bitstype`;
  * mutability (either `type` or `immutable` is printed);
  * the unqualified name of the type;
  * any type parameters;
  * the supertype of the type if it is not `Any`.

# Examples

The generated output for a type definition such as:

```julia
\"""
\$(TYPEDEF)
\"""
immutable MyType{S, T <: Integer} <: AbstractArray
    # ...
end
```

will look similar to the following:

````markdown
```julia
immutable MyType{S, T<:Integer} <: AbstractArray
```
````

!!! note

    No information about the fields of the type is printed. Use the [`FIELDS`](@ref)
    abbreviation to include information about the fields of a type.
"""
const TYPEDEF = TypeDefinition()

function format(::TypeDefinition, buf, doc)
    local binding = doc.data[:binding]
    local object = Docs.resolve(binding)
    if isa(object, DataType)
        println(buf, "\n```julia")
        if isbitstype(object)
            print(buf, "bitstype ", sizeof(object) * 8, " ")
        elseif isabstracttype(object)
            print(buf, "abstract ")
        else
            print(buf, object.mutable ? "type " : "immutable ")
        end
        print(buf, object.name.name)
        if !isempty(object.parameters)
            print(buf, "{")
            join(buf, object.parameters, ", ")
            print(buf, "}")
        end
        local super = supertype(object)
        super == Any ? println(buf) : println(buf, " <: ", super)
        println(buf, "```\n")
    end
end


text:
# the nord theme
# https://www.nordtheme.com/docs
# note that the colors in the original scheme are 0-indexed 😱

loadcolorscheme(:nord, [
    RGB(0.1796875, 0.203125, 0.25),          # nord0
    RGB(0.23046875, 0.2578125, 0.3203125),   # nord1
    RGB(0.26171875, 0.296875, 0.3671875),    # nord2
    RGB(0.296875, 0.3359375, 0.4140625),     # nord3
    RGB(0.84375, 0.8671875, 0.91015625),     # nord4
    RGB(0.89453125, 0.91015625, 0.9375),     # nord5
    RGB(0.921875, 0.93359375, 0.953125),     # nord6
    RGB(0.55859375, 0.734375, 0.73046875),   # nord7
    RGB(0.53125, 0.75, 0.8125),              # nord8
    RGB(0.50390625, 0.62890625, 0.75390625), # nord9
    RGB(0.3671875, 0.50390625, 0.671875),    # nord10
    RGB(0.74609375, 0.37890625, 0.4140625),  # nord11
    RGB(0.8125, 0.52734375, 0.4375),         # nord12
    RGB(0.91796875, 0.79296875, 0.54296875), # nord13
    RGB(0.63671875, 0.7421875, 0.546875),    # nord14
    RGB(0.703125, 0.5546875, 0.67578125),    # nord15
], "nord", "the Nord theme, Polar Night, Snow Storm, Frost, Aurora")

max_stars_repo_name: cormullion/ColorSchemes

max_issues_repo_name: cormullion/ColorSchemes

max_forks_repo_name: cormullion/ColorSchemes


text:
function splitdata(data::AbstractArray, percentval::Number, rng::StableRNG)
    pivot = floor(Int, length(data) * percentval)
    perm = randperm(rng, length(data))
    valind = sort(perm[1:pivot])
    trainind = sort(perm[(pivot + 1):end])
    (train = data[trainind], val = data[valind])
end

max_stars_repo_name: kpmrozowski/wut-mini-deep-learning

max_issues_repo_name: kpmrozowski/wut-mini-deep-learning

max_forks_repo_name: kpmrozowski/wut-mini-deep-learning


text:
@testset "QuadratureVI" begin
    C = 10.0
    i = QuadratureVI(; clipping=C)
    l = GaussianLikelihood()
    @test AGP.grad_quad(l, 200.0, 0.0, 1.0, i) == (C, -C)
end


text:
using Contour, GLVisualize, GeometryTypes, GLAbstraction, Colors, FileIO

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0.0, 1.0, 360))
end
# create a rotation from the time signal
rotation = map(timesignal) do t
    rotationmatrix_z(Float32(t*2pi)) # -> 4x4 Float32 rotation matrix
end

xrange = -5f0:0.02f0:5f0
yrange = -5f0:0.02f0:5f0

z = Float32[sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x) for x in xrange, y in yrange]
mini = minimum(z)
maxi = maximum(z)
color_ramp = map(x->RGBA{Float32}(x, 1.0), colormap("Blues"))

for h in mini:0.2f0:maxi
    c = contour(xrange, yrange, z, h)
    for elem in c.lines
        points = map(elem.vertices) do p
            Point3f0(p, h)
        end
        line_renderable = visualize(
            points, :lines,
            color=color_lookup(color_ramp, h, mini, maxi),
            model=rotation
        )
        _view(line_renderable, window, camera=:perspective)
    end
end

if !isdefined(:runtests)
    renderloop(window)
end


text:
# This code developed by Connor Schembor

module SimplePosetDrawings

using SimplePosets
using PyPlot

import PyPlot.draw
import Base.show

export generateHasse, draw, SimplePosetDrawing

const DEFAULT_MARKER_SIZE = 10
MARKER_SIZE = DEFAULT_MARKER_SIZE

"""
A `SimplePosetDrawing` is a data type representing a Hasse diagram
of a poset.
"""
struct SimplePosetDrawing
  P::SimplePoset
  directConnections::Array{}
  xy::Dict{Any,Vector{Float64}}

  function SimplePosetDrawing(PP::SimplePoset)
    hasse = generateHasse(PP)
    d = createLocations(hasse)
    connections = getDirectConnections(PP)
    new(PP,connections,d)
  end
end

"""Gets only the direct connections of the poset, and returns the direct connections as a list"""
function getDirectConnections(PP::SimplePoset)
  pairs = relations(PP)
  mins = minimals(PP)
  maxs = maximals(PP)

  relationList = []
  for pair in pairs
    tempList = []
    push!(tempList, pair[1])
    push!(tempList, pair[2])
    push!(relationList, tempList)
    for rel in relationList
      if rel[end] == tempList[1]
        #print("inserting after rel: ", rel, "\n")
        t = deepcopy(rel)
        push!(relationList, t)
        push!(relationList[end], tempList[2])
      end
      if rel[1] == tempList[end]
        #print("inserting after pair: ", rel, "\n")
        t = deepcopy(relationList[end])
        push!(relationList, t)
        push!(relationList[end], rel[end])
      end
      #print(relationList, "\n")
    end
  end

  #print("Relation List: ", relationList, "\n")
  """Remove the shorter distance, or 'unnecessary', relations from the relationList"""
  for relation in relationList
    j = length(pairs)
    i = 1
    while i <= j
      if relation[1] == pairs[i][1] && relation[end] == pairs[i][2]
        if length(relation) > length(pairs[i])
          deleteat!(pairs, i)
        end
      end
      j = length(pairs)
      i += 1
    end
  end
  #print(pairs)
  return pairs
end

"""Creates a structure to represent the Hasse diagram. This structure is a dictionary
   that maps each element to all of the elements that 'cover' it"""
function generateHasse(PP::SimplePoset)
  pairs = relations(PP)
  mins = minimals(PP)
  maxs = maximals(PP)

  """Remove reflexive relation pairs"""
  for pair in pairs
    if pair[1] == pair[2]
      filter!(p->p != pair, pairs)
    end
  end

  pairs = getDirectConnections(PP)
  """Construct data structure for the graph"""
  hasseDiagram = Dict()
  for pair in pairs
    if haskey(hasseDiagram, pair[1])
      push!(hasseDiagram[pair[1]], pair[2])
    else
      hasseDiagram[pair[1]] = [pair[2]]
    end
  end
  # print("\n", "hasse diagram: ", hasseDiagram, "\n")
  return hasseDiagram
end

"""Creates coordinate locations for each element in the poset, and returns a coordinate dictionary
   which maps each element to its proper coordinate location"""
function createLocations(hasseDiagram::Dict)
  currentXcoordinate = 0
  currentYcoordinate = 0
  coordinateDict = Dict{Any,Vector{Float64}}()
  FACTOR = 4     #Factor can be adjusted to the desired amount between nodes
  plotted = []
  """Create dictionary of coordinates for each node based on their relation to other nodes"""
  for k in sort(collect(keys(hasseDiagram)))
    #print(k, " - ", hasseDiagram[k])
    len = length(hasseDiagram[k])
    x_tracker = 0 - ((len - 1) * 2)

    if !(k in plotted) && !(hasseDiagram[k] in plotted)
      newX = currentXcoordinate + FACTOR
      positionTaken = true
      while positionTaken
        """Check to see if the position of the node has not been taken"""
        if !([newX,0] in values(coordinateDict))
          coordinateDict[k] = [newX,0]
          push!(plotted, k)
          positionTaken = false
        else
          """If there is already a node there, move the new position in the positive x-direction until
             there is an empty space for it"""
          newX = newX + FACTOR
        end
      end
    elseif !(k in plotted)
      coordinateDict[k] = [currentXcoordinate, currentYcoordinate]
      push!(plotted, k)
    end

    currentXcoordinate = coordinateDict[k][1]
    currentYcoordinate = coordinateDict[k][2]

    for next in hasseDiagram[k]
      if !(next in plotted)
        x = currentXcoordinate + x_tracker
        y = currentYcoordinate + FACTOR
        positionTaken = false
        while !positionTaken
          if ([x,y] in values(coordinateDict))
            x = x + FACTOR
          else
            coordinateDict[next] = [x,y]
            push!(plotted, next)
            positionTaken = true
          end
        end
      end
      x_tracker = x_tracker + FACTOR
    end

  end
  #print(coordinateDict)
  return coordinateDict
end

function edraw(v::Vector{Float64}, w::Vector{Float64})
    plot([v[1],w[1]], [v[2],w[2]],
         color="black", marker="o", markersize=MARKER_SIZE,
         markerfacecolor="white", label="HI", clip_on=false)
end

"""
`draw(X::SimplePosetDrawing)` draws the Hasse diagram of the poset.
"""
function draw(X::SimplePosetDrawing)
  # code (with calls to PyPlot commands) to draw this poset drawing
  conns = X.directConnections
  #print("Direct Conns: ", conns)
  coords = X.xy
  for e in conns
    if length(e) == 2
      edraw(coords[e[1]], coords[e[2]])
    end
  end
  axis("off")
  axis("equal")
  nothing
end

function show(io::IO, X::SimplePosetDrawing)
    print(io,"Drawing of $(X.P)")
end

end # end of module


text:
if VERSION < v"0.6.5"
    Base.load_juliarc()
else
    Base.load_julia_startup()
    using Pkg
    ## needed by console
    const STDIN = stdin
    ## needed by rmarkdown stdout capture
    const STDOUT = stdout
    ## needed by functions like det
    using LinearAlgebra
    ## needed by functions like mean
    using Statistics
end

module JuliaCall

const julia07 = VERSION > v"0.6.5"
const julia14 = VERSION > v"1.3.2"

if julia07
    using Pkg
    ## needed by eval_string function
    const parse = Meta.parse
    ## needed by system checking
    const is_windows = Sys.iswindows
    ## needed by display system
    const Display = AbstractDisplay
    const readstring(s) = read(s, String)
    ## needed by console
    using REPL
    const REPLCompletions = REPL.REPLCompletions

    ## needed by round and signif
    const round1(x::Number, digits) = round(x; digits=digits)
    const signif(x::Number, digits) = round(x; sigdigits=digits)

    ## needed by IRjulia display
    using Base64
else
    ## in julia06
    ## needed by console
    const REPL = Base.REPL
    const REPLCompletions = Base.REPLCompletions

    ## needed by round
    const round1 = round
end

function installed(name)
    @static if julia14
        for p in values(Pkg.dependencies())
            if p.name == name && p.is_direct_dep
                return p.version
            end
        end
        nothing
    elseif julia07
        get(Pkg.installed(), name, nothing)
    else
        Pkg.installed(name)
    end
end

# gc_enable(false)

# Pkg.update()

# if Pkg.installed("RCall") == nothing Pkg.add("RCall") end

using Suppressor

if !is_windows() && !julia07
    @suppress_err begin
        @eval Base JULIA_HOME = joinpath(dirname(JULIA_HOME), "bin")
        @eval Base julia_cmd() = julia_cmd(joinpath(JULIA_HOME, julia_exename()))
    end
end

using RCall

## Fix https://github.com/JuliaInterop/RCall.jl/issues/289
unsafe_store!(cglobal((:R_CStackLimit,RCall.libR),Csize_t), typemax(Csize_t))

const need_display = length(Base.Multimedia.displays) < 2

import Base.Multimedia.display
if need_display
    include("./display/basic.jl")
end
include("./display/IRjulia.jl")
include("./display/RmdJulia.jl")
include("./display/plotsViewer.jl")
include("REPLhook.jl")
include("incomplete_console.jl")
include("RmdStd.jl")

include("convert.jl")
include("JuliaObject.jl")
include("asR.jl")
include("dispatch.jl")

function error_msg(e, bt)
    m = IOBuffer()
    showerror(m, e, bt)
    seek(m, 0)
    readstring(m)
end

function Rerror(e, bt)
    s1 = join(["Error happens in Julia.\n"])
    s2 = error_msg(e, bt)
    s = join([s1 s2])
    rcall(:simpleError, s)
end

function funcfind(name)
    r = Main
    ns = split(name, ".")
    for n in ns
        r = getfield(r, Symbol(n))
    end
    r
end

# function call_decompose(call1)
#     call = rcopy(RObject(Ptr{RCall.VecSxp}(call1)))
#     (call[:fname], call[:named_args], call[:unamed_args], call[:need_return], call[:show_value])
# end

function call_decompose(call1)
    call = Ptr{RCall.VecSxp}(call1)
    # fname = rcopy(String, call[:fname])
    fname = rcopy(String, call[1]) :: String
    # named_args = Any[(rcopy(Symbol, k), rcopy(i)) for (k, i) in enumerate(call[:named_args])]
    # unamed_args = Any[rcopy(i) for i in call[:unamed_args]]
    args = call[2] :: Ptr{RCall.VecSxp}
    symbols = rcopy(Vector{Symbol}, getnames(args)) :: Vector{Symbol}
    es = rcopy(Vector{Any}, args) :: Vector{Any}
    if length(symbols) == 0
        unamed_args = es
        named_args = Any[]
    else
        unamed_args = Any[es[i] for i in 1:length(es) if symbols[i] == Symbol("")]
        named_args = Any[(symbols[i], es[i]) for i in 1:length(es) if symbols[i] != Symbol("")]
    end

    # named_args = Any[]
    # unamed_args = Any[]
    # # for (k,a) in enumerate(call[:args])
    # for (k,a) in enumerate(args)
    #     if isa(k.p, Ptr{RCall.NilSxp})
    #         push!(unamed_args, rcopy(a))
    #     else
    #         push!(named_args, (rcopy(Symbol,k), rcopy(a)))
    #     end
    # end
    # need_return = rcopy(String, call[:need_return])
    need_return = rcopy(String, call[3]) :: String
    # show_value = rcopy(Bool, call[:show_value])
    show_value = rcopy(Bool, call[4]) :: Bool
    (fname, named_args, unamed_args, need_return, show_value)
end

function docall(call1)
    try
        fname, named_args, unamed_args, need_return, show_value = call_decompose(call1);
        if endswith(fname, ".")
            fname = chop(fname);
            # f = eval(Main, parse(fname));
            f = funcfind(fname);
            r = broadcast(f, unamed_args...);
        else
            # f = eval(Main, parse(fname));
            f = funcfind(fname);
            r = f(unamed_args...; named_args...);
        end
        if show_value && r != nothing
            display(r)
        end
        # @static if need_display
        #     proceed(basic_display_manager)
        # end
        if need_return == "R"
            sexp(r);
        elseif need_return == "Julia"
            sexp(JuliaObject(r));
        else
            sexp(nothing);
        end;
    catch e
        Rerror(e, stacktrace(catch_backtrace())).p;
    end;
end

include("interface1.jl")

function exists(x)
    isdefined(Main, Symbol(x))
end

function eval_string(x)
    Core.eval(Main, parse(strip(x)))
end

function installed_package(pkg_name)
    try
        string(installed(pkg_name))
    catch e
        "nothing";
    end
end

function help(fname)
    string(eval_string(join(["@doc " fname])))
end

function assign(name, x)
    Core.eval(Main, Expr(:(=), Symbol(name), x))
end

function str_typeof(x)
    string(typeof(x))
end

function show_string(x)
    buf = IOBuffer()
    show(IOContext(buf, :limit=>true), x)
    return String(take!(buf))
end

## Needed by julia_source
function include1(fname)
    @static if julia07
        Base.include(Main, fname)
    else
        include(fname)
    end
end

# function apply(f, args...)
#     f(args...)
# end

function apply(f, unamed_args...; named_args...)
    f(unamed_args...; named_args...)
end

end

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: KevCaz/JuliaCall

max_forks_repo_name: KevCaz/JuliaCall


text:
include("src/nfxp.jl")
nfxp.maxlike()
p = nfxp.Param()
nfxp.VFI!(p)
data = nfxp.simdata(50,119,p)
include("test/runtests.jl")
nfxp.max_optim()
nfxp.max_nlopt()

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack