text:
using LinearAlgebra, StaticArrays, Attitude, SatelliteDynamics
using Infiltrator, ProgressMeter
using MATLAB

function initial_conditions()

    # Declare initial state in terms of osculating orbital elements
    rp = R_EARTH + 200e3;ra = R_EARTH + 35000e3
    a = (ra+rp)/2
    e = (ra-rp)/(ra+rp)
    i = deg2rad(27)
    # i = 0
    ω = deg2rad(34)
    Ω = deg2rad(133)
    M = deg2rad(180)
    oe0  = [a,e,i,Ω,ω,M]

    # Convert osculating elements to Cartesean state
    eci0 = sOSCtoCART(oe0, use_degrees=false)

    # modified equinoctial elements
    p = a*(1-e^2)
    f = e*cos(ω + Ω)
    g = e*sin(ω + Ω)
    h = tan(i/2)*cos(Ω)
    k = tan(i/2)*sin(Ω)
    L = wrap_to_2pi(ω + Ω + M)

    eoe0 = [p;f;g;h;k;L]

    r_eci = eci0[1:3]
    v_eci = eci0[4:6]

    p_0 = u_from_x(r_eci)                          # u
    p_prime_0 = uprime_from_xdot(v_eci,p_0)     # du/ds
    h_0 = GM_EARTH/norm(r_eci) - .5*norm(v_eci)^2         # m^2/s^2
    t_0 = 0

    ks0 = [p_0;p_prime_0;h_0;t_0]

    return eci0, eoe0, ks0
end
function rk4(f::Function, x_n, dt)
    k1 = dt*f(x_n)
    k2 = dt*f(x_n+k1/2)
    k3 = dt*f(x_n+k2/2)
    k4 = dt*f(x_n+k3)
    return (x_n + (1/6)*(k1+ 2*k2 + 2*k3 + k4))
end
function fode_j2(r_eci)
    r = norm(r_eci)
    x,y,z = r_eci
    J2 = J2_EARTH
    μ = GM_EARTH

    # precompute repeated stuff
    Re_r_sqr = 1.5*J2*(R_EARTH/r)^2
    five_z_sqr = 5*z^2/r^2

    return  SA[(μ/r^3)*x*(Re_r_sqr*(five_z_sqr - 1)),
               (μ/r^3)*y*(Re_r_sqr*(five_z_sqr - 1)),
               (μ/r^3)*z*(Re_r_sqr*(five_z_sqr - 3))]
end
function cartesian_dynamics(x)
    r_eci = SA[x[1],x[2],x[3]]
    v_eci = SA[x[4],x[5],x[6]]

    # relavent parameters
    J2 = J2_EARTH
    μ = GM_EARTH

    # eci position stuff
    r = norm(r_eci)
    x,y,z = r_eci

    # precompute repeated stuff
    Re_r_sqr = 1.5*J2*(R_EARTH/r)^2
    five_z_sqr = 5*z^2/r^2

    return  SA[v_eci[1],v_eci[2],v_eci[3],
               (-μ/r^3)*x*(1 - Re_r_sqr*(five_z_sqr - 1)),
               (-μ/r^3)*y*(1 - Re_r_sqr*(five_z_sqr - 1)),
               (-μ/r^3)*z*(1 - Re_r_sqr*(five_z_sqr - 3))]
end
function J2_from_EOE(x)
    """J2 acceleration in RTN from equinoctial orbital elements.
    Args:
        p: a*(1-e^2) semi-latus rectum       (m)
        f: e*cos(ω + Ω)                      (non-dimensional)
        g: e*sin(ω + Ω)                      (non-dimensional)
        h: tan(i/2)*cos(Ω)                   (non-dimensional)
        k: tan(i/2)*sin(Ω)                   (non-dimensional)
        L: ω + Ω + θ argument of latitude    (rad)
    Returns:
        a_rtn: acceleration in the RTN frame  (m/s)
    """

    # unpack state
    p,f,g,h,k,L = x
    J2 =  J2_EARTH

    # precompute trig functions
    sinL = sin(L)
    cosL = cos(L)

    # constants for EOE
    mu = GM_EARTH
    w = 1+f*cosL + g*sinL
    r = p/w
    Re = R_EARTH

    # accelerations
    a_r = -((3*mu*J2*Re^2)/(2*r^4))*(1 -   (  (12*(h*sinL-k*cosL)^2 )
                                                    /   ((1+h^2 + k^2)^2))  )
    a_t = -((12*mu*J2*Re^2)/(r^4))*((  ((h*sinL - k*cosL)*(h*cosL+k*sinL) )
                                                    /  ((1+h^2 + k^2)^2))  )
    a_n = -((6*mu*J2*Re^2)/(r^4))*((  ((h*sinL - k*cosL)*(1-h^2-k^2) )
                                                    /  ((1+h^2 + k^2)^2))  )
    return a_r, a_t, a_n
end
function eoe_dynamics(x)
    """Equations of motion for equinoctial elements in 2-body problem.
    Args:
        epc: time(::epoch)
        x: state
            - p: a*(1-e^2) semi-latus rectum       (m)
            - f: e*cos(ω + Ω)                      (non-dimensional)
            - g: e*sin(ω + Ω)                      (non-dimensional)
            - h: tan(i/2)*cos(Ω)                   (non-dimensional)
            - k: tan(i/2)*sin(Ω)                   (non-dimensional)
            - L: ω + Ω + θ argument of latitude    (rad)
            - mass: spacecraft mass                (kg)
        u: thruster acceleration                   (m/s^2)
        thruster_on: bool for thruster             ()
    Returns:
        ̇x: state derivative
    """

    # unpack state
    p,f,g,h,k,L = x

    # equinoctial parameters for GVE's
    sinL = sin(L)
    cosL = cos(L)
    alpha_2 = h^2 - k^2
    s2 = 1+h^2 + k^2
    w = 1+f*cosL + g*sinL
    r = p/w

    # precompute to save calculations
    sqrtpu = sqrt(p/GM_EARTH)

    # J2
    u_r, u_t, u_n = J2_from_EOE(x)

    # GVE's for equinoctial orbital elements
    p_dot = (2*p/w)*sqrtpu*u_t
    f_dot = sqrtpu*( u_r*sinL    +    ((w+1)*cosL + f)*u_t/w    -
                                            (h*sinL - k*cosL)*g*u_n/w  )
    g_dot = sqrtpu*(-u_r*cosL   +    ((w+1)*sinL + g)*u_t/w    +
                                            (h*sinL - k*cosL)*f*u_n/w  )
    h_dot = sqrtpu*(s2*u_n/(2*w))*cosL
    k_dot = sqrtpu*(s2*u_n/(2*w))*sinL
    L_dot = sqrt(GM_EARTH*p)*(w/p)^2 + (1/w)*sqrtpu*(h*sinL - k*cosL)*u_n

    # return state derivative
    return SA[p_dot,f_dot,g_dot,h_dot,k_dot,L_dot]
end
function rv_from_eoe(eoe_vec)
    """Position and Velocity resolved in ECI
    Args:
        p: a*(1-e^2) semi-latus rectum     (m)
        f: e*cos(ω + Ω)                    (non-dimensional)
        g: e*sin(ω + Ω)                    (non-dimensional)
        h: tan(i/2)*cos(Ω)                 (non-dimensional)
        k: tan(i/2)*sin(Ω)                 (non-dimensional)
        L: ω + Ω + θ argument of latitude  (rad)
    Returns:
        r: position vector from eci origin to spacecraft, resolved in eci. (m)
        v: velocity vector of spacecraft in eci, resolved in eci.          (m/s)
    """

    # unpack eoe
    p,f,g,h,k,L = eoe_vec

    # precompute sin and cosine
    sinL = sin(L)
    cosL = cos(L)

    # various eoe specific parameters
    alf2 = h^2 - k^2
    s2 = 1 + h^2 + k^2
    w = 1 + f*cosL + g*sinL
    r = p/w

    # eci position vector
    r = SA[(r/s2) * (  cosL + alf2*cosL + 2*h*k*sinL),
         (r/s2) * (  sinL - alf2*sinL + 2*h*k*cosL),
         (2*r/s2)*(h*sinL -   k *cosL)]

    # ECI velocity vector
    v = sqrt(GM_EARTH/p)*SA[-(1/(s2)) * (sinL + alf2*sinL - 2*h*k*cosL + g - 2*f*h*k + alf2*g),
         -(1/(s2))* (-cosL + alf2*cosL + 2*h*k*sinL - f + 2*g*h*k + alf2*f),
          (2/s2) * (h*cosL +    k*sinL + f*h + g*k)]

    # @infiltrate
    return SA[r[1],r[2],r[3],v[1],v[2],v[3]]
end

function ks_dynamics(x)

    # unpack state
    p = SVector(x[1],x[2],x[3],x[4])
    p_prime = SVector(x[5],x[6],x[7],x[8])
    h = x[9]
    t = x[10]

    # position in eci
    r_eci = x_from_u(p)

    # j2 accel
    a_j2 = fode_j2(r_eci)

    L = L_fx(p)

    # append 0 to end of P_j2 acceleration vector
    P_j2 = SA[a_j2[1],a_j2[2],a_j2[3],0]

    p_dp = -(h/2)*p + ((dot(p,p))/(2))*L'*P_j2

    return SVector(p_prime[1],p_prime[2],p_prime[3],p_prime[4],
                   p_dp[1],p_dp[2],p_dp[3],p_dp[4],
                   -2*dot(p_prime,L'*P_j2),
                   dot(p,p))
end
function rv_from_ks(p,p_prime)
    x = x_from_u(p)
    xdot = xdot_from_u(p,p_prime)
    return SA[x[1],x[2],x[3],xdot[1],xdot[2],xdot[3]]
end
function ϵ(rv)
    r = SA[rv[1],rv[2],rv[3]]
    v = SA[rv[4],rv[5],rv[6]]
    return dot(v,v)/2 - GM_EARTH/norm(r)
end
function P2(x)
    return .5*(3*x^2 -1)
end
function true_anom_from_rv(r,v)
    h = cross(r,v)
    e = (v × h)/GM_EARTH - r/norm(r)

    ν = acos(dot(e,r)/(norm(e)*norm(r)))
    if dot(r,v)<0
        ν = 2*pi - ν
    end
    return ν
end
function sma_from_eci(rv)
    rv = convert(Array{Float64},rv)

    try
        r = rv[1:3]
        v = rv[4:6]
        oe = sCARTtoOSC(rv,use_degrees = false)
        a,e,i,Ω,ω,M = oe
        θ = true_anom_from_rv(r,v)
        cosϕ = sin(i)*sin(θ + ω)
        U_j2 = - (GM_EARTH/norm(r))*J2_EARTH*(R_EARTH/norm(r))^2*P2(cosϕ)
        ϵ = dot(v,v)/2 - GM_EARTH/norm(r) - U_j2
        # return oe[1]
        return ϵ
    catch err
        return NaN
    end

    return oe[1]
end
function driver()
    eci0, eoe0, ks0 = initial_conditions()

    num_orbits = 6
    tf = num_orbits*36951.67397573928
    # N = 1000
    N = 600
    # N = Int(round((tf/(24*3600))*80))
    @show N
    dt = tf/N

    t_vec = [(i-1)*dt for i = 1:N]

    X_eci = [@SVector zeros(6) for i = 1:N]
    X_eci[1] = SVector{6}(eci0)
    X_eoe = [@SVector zeros(6) for i = 1:N]
    X_eoe[1] = eoe0
    X_ks = [@SVector zeros(10) for i = 1:N]
    X_ks[1] = ks0

    # one day in fictitious time (this was guess and check)
    sf = num_orbits*0.001563905
    ds = sf/N

    @inbounds for i = 1:(N-1)
        X_eci[i+1] = rk4(cartesian_dynamics,X_eci[i],dt)
        X_eoe[i+1] = rk4(eoe_dynamics,X_eoe[i],dt)
        X_ks[i+1] = rk4(ks_dynamics,X_ks[i],ds)
    end

    ksm = mat_from_vec(X_ks)
    goal_tf = 36951.67397573928

    @show goal_tf - X_ks[end][10]

    ks_t_vec = vec(ksm[10,:])

    X_eoe_eci = [rv_from_eoe(X_eoe[i]) for i = 1:length(X_eoe)]
    X_ks_eci = [rv_from_ks(X_ks[i][1:4],X_ks[i][5:8]) for i = 1:length(X_ks)]

    eoe_m = mat_from_vec(X_eoe_eci)
    eci_m = mat_from_vec(X_eci)
    ks_m = mat_from_vec(X_ks_eci)

    # SMA from all
    eoe_a= [sma_from_eci(X_eoe_eci[i]) for i = 1:length(X_eoe_eci)]
    eci_a = [sma_from_eci(X_eci[i]) for i = 1:length(X_eci)]
    ks_a = [sma_from_eci(X_ks_eci[i]) for i = 1:length(X_ks_eci)]

    T_period = 36951.67397573928
    t_vec /= T_period
    ks_t_vec /= T_period
    eps_0 = eci_a[1]
    pl1 = 1000*100*(1 .- eoe_a /eps_0)
    pl2 = 1000*100*(1 .- eci_a /eps_0)
    pl3 = 1000*100*(1 .- ks_a /eps_0)
    mat"
    figure
    hold on
    plot($t_vec , $pl1,'linewidth',3)
    plot($t_vec , $pl2)
    plot($ks_t_vec, $pl3)
    legend('EOE','Cartesian','KS')
    %xlim([0 5])
    %ylim([-1.5 2.8])
    %xticks([0 1 2 3 4 5])
    %xlabel('Number of Orbits')
    %ylabel('Specific Orbital Energy (kJ/kg)')
    %set(gca, 'XScale', 'log')
    addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
    matlab2tikz('energy_test3.tikz')
    hold off
    "
end

driver()

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
using Flux
using Flux: mse, back!
using Flux.Optimise: _update_params!
using Statistics
using Formatting
using Random
using StatsBase

"""
VanillaPG

    # Notes
    - VPG chooses a random action according to its policy function, π.
        - Can learn a (nearly) deterministic policy if that is optimal.
"""
mutable struct VanillaPG
    model  # model
    opt  # optimizer ("watching" model)
    batch_size  # episodes per batch
    γ  # discount rate
end

action(agent::VanillaPG, env::AbstractEnvironment) = sample(env.action_space, weights(agent.model(env.state)))

function compute_loss(agent::VanillaPG, states, actions, rewards)
    states = hcat(states...)
    index = CartesianIndex.(actions, axes(states, 2))
    return sum(cumsum(rewards) .* log.(agent.model(states)[index]))
end

function update!(agent::VanillaPG, loss)
    loss = -mean(loss)
    back!(loss)
    _update_params!(agent.opt, Flux.params(agent.model))
end

function run_episode(agent::VanillaPG, env::AbstractEnvironment; rendered = false, verbose = false)
    total_reward = 0.0
    steps = 0
    s′, done = reset!(env)
    rendered && render(env)
    states = []
    actions = []
    rewards = []
    while !done
        s = deepcopy(s′)
        a = action(agent, env)
        push!(states, s)
        push!(actions, a)
        s′, r, done, info = step!(env, a)
        steps += 1
        total_reward += r
        push!(rewards, agent.γ^steps * r)  # NOTE: **discounted** rewards
        verbose && @info "s = $s, a = $a, r = $r, s′ = $(s′), done = $done"
        steps == env.maxsteps && break
    end
    verbose && @info "episode finished in $steps timesteps"
    loss = compute_loss(agent, states, actions, rewards)
    return total_reward, steps, loss
end

function evaluate(agent::VanillaPG, env::AbstractEnvironment, n = 100; verbose = false, rendered = false)
    env.rendered = rendered
    history = []
    @progress for i = 1:n
        total_reward, steps, _ = run_episode(agent, env, verbose = verbose)
        push!(history, (total_reward = total_reward, steps = steps))
    end
    total_reward = mean([e.total_reward for e in history])
    steps = mean([e.steps for e in history])
    env.rendered = false
    return total_reward, steps
end

function train!(agent::VanillaPG, env::AbstractEnvironment; max_episodes = Inf, eval_freq = 1, eval_episodes = 100, verbose = false, patience = 5)

    @info "training agent..."
    recorder = ExperimentRecorder("./log-$(randstring(5))")
    current_best_score = -Inf
    n_evals_since_best_score = 0
    improving = true
    batch_size = agent.batch_size
    batches = 0
    episodes = 0

    while improving && episodes < max_episodes

        loss = Array{Flux.Tracker.TrackedReal{Float64}, 1}(undef, batch_size)
        for i = 1:batch_size
            # training episode
            total_reward, steps, loss[i] = run_episode(agent, env)
            episodes += 1
            verbose && @info "<training status message>"
            message = (episode=episodes, score=total_reward, steps=steps, epsilon=NaN)
            stamp!(recorder, message, :train)
        end

        # network update
        update!(agent, loss)
        batches += 1

        # agent.opt.eta = agent.opt.eta * 0.99

        # fitness evaluation
        if batches % eval_freq == 0
            total_reward, steps = evaluate(agent, env, eval_episodes)
            @info "episode: $episodes, score: $total_reward, η: $(agent.opt.eta)"
            message = (episode=episodes, score=total_reward, steps=steps, epsilon=NaN)
            stamp!(recorder, message, :test)

            # snapshots
            # if score > current_best_score
            #     current_best_score = score
            #     # save a snapshot of agent...
            # else
            #     n_evals_since_best_score += 1
            #     n_evals_since_best_score == patience && break
            # end
        end
    end
    verbose && @info "agent finished training in $episodes episodes"
    return recorder
end

ext: jl


text:
mutable struct DynLine <: DynBranch
    name::Symbol
    available::Bool
    arc::PSY.Arc
    r::Float64  # System per-unit value
    x::Float64  # System per-unit value
    b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}  # System per-unit value
    n_states::Int64
    states::Vector{Symbol}

                function DynLine(name::String,
                                 available::Bool,
                                 arc::PSY.Arc,
                                 r::Float64, # System per-unit value
                                 x::Float64, # System per-unit value
                                 b::NamedTuple{(:from, :to), Tuple{Float64, Float64}})

                    n_states = 2
                    states = Vector{Symbol}()
                    states = [:Il_R
                              :Il_I]

                    new(Symbol(name),
                        available,
                        arc,
                        r,
                        x,
                        b,
                        n_states,
                        states)

                end

end


text:
#Data generation Cubic
batchsize = Integer(1e4)
generateHistogram(intensityCube; qcut=intensityCube.rmax, K2=KMAX, K3=K, N=N, max_pictures=batchsize, batchsize = batchsize, number_incident_photons=calculate_incident_photons(10), numprocesses=1, file="/tmp/correlations.dat", noise=GaussianNoise(0.1, 0.5, 14), lambda=2.0)
loadHistograms(8,8, "/tmp/correlations.dat")

# Data generation spherical harmonics
# volume = getSurfaceVolume(intensity)
# volume.radial_interp = false
# try rm("/tmp/correlations.dat") end
# @profile generateHistogram(volume; qcut=intensity.rmax, K2=KMAX, K3=K, N=N, max_pictures=batchsize, batchsize = batchsize, number_incident_photons=calculate_incident_photons(10), numprocesses=nworkers(), file="/tmp/correlations.dat", noise=GaussianNoise(0.1, intensity.rmax/4.0, 10))
# loadHistograms(8,8, "/tmp/correlations.dat")


text:
using Test
using Clustering

@testset "V-measure" begin
    @testset "reproducing fig.2" begin
        # Tests are taken from the fig. 2 of the referenced paper:
        # V-Measure: A conditional entropy-based external cluster evaluation measure,
        # Andrew Rosenberg and Julia Hirschberg

        clus = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
        v = vmeasure(clus, clus)
        @test v == 1.0

        clas = [1, 1, 1, 2, 3, 3, 3, 3, 1, 2, 2, 2, 2, 1, 3]
        v = vmeasure(clas, clus)
        @test v ≈ 0.14 atol=1e-2

        clas = [1, 1, 1, 2, 2, 3, 3, 3, 1, 1, 2, 2, 2, 3, 3]
        v = vmeasure(clas, clus)
        @test v ≈ 0.39 atol=1e-2

        clus = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 5, 6, 6]
        clas = [1, 1, 1, 2, 2, 3, 3, 3, 1, 1, 2, 2, 2, 3, 3, 1, 2, 3, 1, 2, 3]
        v = vmeasure(clas, clus)
        @test v ≈ 0.30 atol=1e-2

        clus = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 6, 7, 8, 9]
        v = vmeasure(clas, clus)
        @test v ≈ 0.41 atol=1e-2

        @test_throws ArgumentError vmeasure(clas, clus, β = -1.0)
    end

    @testset "comparing 2 k-means clusterings" begin
        Random.seed!(34568) # set random seed for RNG used by kmeans()
        rng = StableRNG(34568)
        m = 3
        n = 1000
        k = 10

        # non-weighted
        v = mean([begin
            x = rand(rng, m, n)
            vmeasure(kmeans(x, k; maxiter=50),
                     kmeans(x, k; maxiter=50))
        end for _ in 1:200])
        @test 0.5 < v < 1.0
        @test v ≈ 0.75 atol=1e-2 # FIXME why 0.75?
    end

    @testset "comparing 2 random label assignments" begin
        rng = StableRNG(34568)
        k = 10
        n = 10000

        a1 = rand(rng, 1:k, n)
        a2 = rand(rng, 1:k, n)
        v = vmeasure(a1, a2)
        @test v ≈ 0.0 atol=1e-2 # should be close to zero
    end
end


text:
#=
hdf5utils.jl
Utilities for reading HDF5 data files for the 2019 fastMRI challenge.
Code by Steven Whitaker
Documentation by Jeff Fessler
=#

export h5_get_keys, h5_get_attributes, h5_get_ismrmrd
export h5_get_ESC, h5_get_RSS, h5_get_kspace

using HDF5

"""
`h5_get_keys(filename::String)`
Get `names` from file.  Returns a `?`
"""
function h5_get_keys(filename::String)

	return h5open(filename, "r") do file
		names(file)
	end

end

"""
`h5_get_attributes(filename::String)`
Get `attrs` from file.  Returns a `Dict`
"""
function h5_get_attributes(filename::String)

	return h5open(filename, "r") do file
		a = attrs(file)
		attr = Dict{String,Any}()
		for s in names(a)
			attr[s] = read(a[s])
		end
		attr
	end

end

"""
`h5_get_ismrmrd(filename::String)`
Get ISMRM header data from file.  Returns a `?`
"""
function h5_get_ismrmrd(filename::String)

	return h5read(filename, "ismrmrd_header")

end

"""
`h5_get_ESC(filename::String; T::DataType = ComplexF32)`
Return `Array` of ESC (emulated single coil) data from file.

HDF5.jl reads the data differently than Python's h5py.
This is significant because the fastMRI paper
says that the datasets have certain dimensionality,
but the dimensions are permuted in Julia compared to Python,
hence the calls to `permutedims` in the following functions.
"""
function h5_get_ESC(filename::String; T::DataType = ComplexF32)

	data = T.(h5read(filename, "reconstruction_esc"))
	return permutedims(data, ndims(data):-1:1)

end

"""
`h5_get_RSS(filename::String; T::DataType = ComplexF32)`
Return `Array` of RSS (root sum of squares) data from file.
"""
function h5_get_RSS(filename::String; T::DataType = ComplexF32)

	data = T.(h5read(filename, "reconstruction_rss"))
	return permutedims(data, ndims(data):-1:1)

end

"""
`h5_get_kspace(filename::String; T::DataType = ComplexF32))`
Return `Array` of kspace data from file.
"""
function h5_get_kspace(filename::String; T::DataType = ComplexF32)

	data = h5open(filename, "r") do file
		readmmap(file["kspace"], Array{getcomplextype(file["kspace"])})
	end
	return permutedims(T.(data), ndims(data):-1:1)

end

#=
Copied (basically) from
https://github.com/MagneticParticleImaging/MPIFiles.jl/blob/79711bf7af389f9e2dd4b0370e64040e5da1e193/src/Utils.jl#L33
=#
function getcomplextype(dataset)
	T = HDF5.hdf5_to_julia_eltype(
		HDF5Datatype(HDF5.h5t_get_member_type(datatype(dataset).id, 0)))
	return Complex{T}
end


text:
#

using Flux
using RoME

# loading data
fg = initfg()
loadDFG("/tmp/caesar/2020-05-11T02:25:53.702/fg_204_resolve.tar.gz", Main, fg)

# first two poses well separated from the others

#


text:
using Unmarshal
using Test
using JSON
using LazyJSON

# write your own tests here

# Unmarshal tests

println("Start testing with library JSON")
include("unmarshal.jl")
println("Done testing with library JSON")

println("Start testing with library LazyJSON")
include("unmarshalLazyJSON.jl")
println("Done testing with library LazyJSON")


text:
# host array

export VEArray

@enum ArrayState begin
  ARRAY_UNMANAGED
  ARRAY_MANAGED
  ARRAY_FREED
end

mutable struct VEArray{T,N} <: AbstractGPUArray{T,N}
  buf::Union{Nothing,Mem.Device}
  dims::Dims{N}

  state::ArrayState

  #ctx::VEContext
  #dev::VEDevice

  function VEArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
    Base.isbitsunion(T) && error("VEArray does not yet support union bits types")
    Base.isbitstype(T)  || error("VEArray only supports bits types") # allocatedinline on 1.3+
    #ctx = context()
    #dev = device()
    buf = Mem.Device(prod(dims) * sizeof(T)) # not needed # , Base.datatype_alignment(T))
    obj = new{T,N}(buf, dims, ARRAY_MANAGED) # add these later # , ctx, dev)
    finalizer(unsafe_free!, obj)
    return obj
  end
end

function unsafe_free!(xs::VEArray)
  # this call should only have an effect once, because both the user and the GC can call it
  if xs.state == ARRAY_FREED
    return
  elseif xs.state == ARRAY_UNMANAGED
    throw(ArgumentError("Cannot free an unmanaged buffer."))
  end

  finalize(xs.buf)
  xs.state = ARRAY_FREED

  # the object is dead, so we can also wipe the pointer
  xs.buf = nothing

  return
end

#device(a::VEArray) = a.dev
#context(a::VEArray) = a.ctx

## alias detection

Base.dataids(A::VEArray) = (UInt(pointer(A)),)

Base.unaliascopy(A::VEArray) = copy(A)

## convenience constructors

VEVector{T} = VEArray{T,1}
VEMatrix{T} = VEArray{T,2}
VEVecOrMat{T} = Union{VEVector{T},VEMatrix{T}}

# type and dimensionality specified, accepting dims as series of Ints
VEArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = VEArray{T,N}(undef, dims)

# type but not dimensionality specified
VEArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = VEArray{T,N}(undef, dims)
VEArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    VEArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
VEArray{T,1}() where {T} = VEArray{T,1}(undef, 0)

# do-block constructors
for (ctor, tvars) in (:VEArray => (), :(VEArray{T}) => (:T,), :(VEArray{T,N}) => (:T, :N))
  @eval begin
    function $ctor(f::Function, args...) where {$(tvars...)}
      xs = $ctor(args...)
      try
        f(xs)
      finally
        unsafe_free!(xs)
      end
    end
  end
end

Base.similar(a::VEArray{T,N}) where {T,N} = VEArray{T,N}(undef, size(a))
Base.similar(a::VEArray{T}, dims::Base.Dims{N}) where {T,N} = VEArray{T,N}(undef, dims)
Base.similar(a::VEArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = VEArray{T,N}(undef, dims)

function Base.copy(a::VEArray{T,N}) where {T,N}
  b = similar(a)
  @inbounds copyto!(b, a)
end

## array interface

Base.elsize(::Type{<:VEArray{T}}) where {T} = sizeof(T)

Base.size(x::VEArray) = x.dims
Base.sizeof(x::VEArray) = Base.elsize(x) * length(x)

## derived types

export VEDenseArray, VEDenseVector, VEDenseMatrix, VEDenseVecOrMat,
       VEDenseArray, VEDenseVector, VEDenseMatrix, VEDenseVecOrMat,
       VEWrappedArray, VEWrappedVector, VEWrappedMatrix, VEWrappedVecOrMat

VEContiguousSubArray{T,N,A<:VEArray} = Base.FastContiguousSubArray{T,N,A}

# dense arrays: stored contiguously in memory
VEDenseReinterpretArray{T,N,A<:Union{VEArray,VEContiguousSubArray}} = Base.ReinterpretArray{T,N,S,A} where S
VEDenseReshapedArray{T,N,A<:Union{VEArray,VEContiguousSubArray,VEDenseReinterpretArray}} = Base.ReshapedArray{T,N,A}
DenseSubVEArray{T,N,A<:Union{VEArray,VEDenseReshapedArray,VEDenseReinterpretArray}} = Base.FastContiguousSubArray{T,N,A}
VEDenseArray{T,N} = Union{VEArray{T,N}, DenseSubVEArray{T,N}, VEDenseReshapedArray{T,N}, VEDenseReinterpretArray{T,N}}
VEDenseVector{T} = VEDenseArray{T,1}
VEDenseMatrix{T} = VEDenseArray{T,2}
VEDenseVecOrMat{T} = Union{VEDenseVector{T}, VEDenseMatrix{T}}

# strided arrays
VEStridedSubArray{T,N,A<:Union{VEArray,VEDenseReshapedArray,VEDenseReinterpretArray},
                  I<:Tuple{Vararg{Union{Base.RangeIndex, Base.ReshapedUnitRange,
                                        Base.AbstractCartesianIndex}}}} = SubArray{T,N,A,I}
VEStridedArray{T,N} = Union{VEArray{T,N}, VEStridedSubArray{T,N}, VEDenseReshapedArray{T,N}, VEDenseReinterpretArray{T,N}}
VEStridedVector{T} = VEStridedArray{T,1}
VEStridedMatrix{T} = VEStridedArray{T,2}
VEStridedVecOrMat{T} = Union{VEStridedVector{T}, VEStridedMatrix{T}}

Base.pointer(x::VEStridedArray{T}) where {T} = Base.unsafe_convert(VEPtr{T}, x)
@inline function Base.pointer(x::VEStridedArray{T}, i::Integer) where T
    Base.unsafe_convert(VEPtr{T}, x) + Base._memory_offset(x, i)
end

# wrapped arrays: can be used in kernels
VEWrappedArray{T,N} = Union{VEArray{T,N}, WrappedArray{T,N,VEArray,VEArray{T,N}}}
VEWrappedVector{T} = VEWrappedArray{T,1}
VEWrappedMatrix{T} = VEWrappedArray{T,2}
VEWrappedVecOrMat{T} = Union{VEWrappedVector{T}, VEWrappedMatrix{T}}

## interop with other arrays

@inline function VEArray{T,N}(xs::AbstractArray{<:Any,N}) where {T,N}
  A = VEArray{T,N}(undef, size(xs))
  copyto!(A, convert(Array{T}, xs))
  return A
end

# underspecified constructors
VEArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = VEArray{T,N}(xs)
(::Type{VEArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = VEArray{S,N}(x)
VEArray(A::AbstractArray{T,N}) where {T,N} = VEArray{T,N}(A)

# idempotency
VEArray{T,N}(xs::VEArray{T,N}) where {T,N} = xs

# Level Zero references
VERef(x::Any) = VERefArray(VEArray([x]))
VERef{T}(x) where {T} = VERefArray{T}(VEArray(T[x]))
VERef{T}() where {T} = VERefArray(VEArray{T}(undef, 1))

## conversions

Base.convert(::Type{T}, x::T) where T <: VEArray = x

## interop with C libraries

Base.unsafe_convert(::Type{Ptr{T}}, x::VEArray{T}) where {T} =
  throw(ArgumentError("cannot take the host address of a $(typeof(x))"))
Base.unsafe_convert(::Type{VEPtr{T}}, x::VEArray{T}) where {T} = reinterpret(VEPtr{T}, pointer(x.buf))

## interop with GPU arrays

function Base.unsafe_convert(::Type{VEDeviceArray{T,N,AS.Global}}, a::VEDenseArray{T,N}) where {T,N}
    VEDeviceArray{T,N,AS.Global}(size(a), reinterpret(LLVMPtr{T,AS.Global}, a.buf.ptr))
end

Adapt.adapt_storage(::Adaptor, xs::VEArray{T,N}) where {T,N} =
  Base.unsafe_convert(VEDeviceArray{T,N,AS.Global}, xs)

# we materialize ReshapedArray/ReinterpretArray/SubArray/... directly as a device array
Adapt.adapt_structure(::Adaptor, xs::VEDenseArray{T,N}) where {T,N} =
  Base.unsafe_convert(VEDeviceArray{T,N,AS.Global}, xs)

## interop with CPU arrays

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{VEArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(VEArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:VEArray{T}}, xs::AbstractArray) where {T} =
  isbits(xs) ? xs : convert(VEArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::VEArray) = convert(Array, xs)

Base.collect(x::VEArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.copyto!(dest::VEArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::VEDenseArray{T}, src::Array{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::VEDenseArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::Array{T}, src::VEDenseArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::VEDenseArray{T}, doffs::Integer, src::VEDenseArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  #@assert device(dest) == device(src) && context(dest) == context(src)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::VEDenseArray{T}, src::VEDenseArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.unsafe_copyto!(      # not needed # ctx::ZeContext, dev::ZeDevice,
                             dest::VEDenseArray{T}, doffs, src::Array{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(      # not needed # ctx::ZeContext, dev::ZeDevice,
                             dest::Array{T}, doffs, src::VEDenseArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end

  # copies to the host are synchronizing
  #synchronize(global_queue(context(src), device(src)))
  synchronize()

  return dest
end

function Base.unsafe_copyto!(         # not needed # ctx::ZeContext, dev::ZeDevice,
                             dest::VEDenseArray{T}, doffs, src::VEDenseArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

## utilities

zeros(T::Type, dims...) = fill!(VEArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(VEArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float64, dims...)
ones(dims...) = ones(Float64, dims...)
fill(v, dims...) = fill!(VEArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(VEArray{typeof(v)}(undef, dims...), v)

function Base.fill!(A::VEDenseArray{T}, val) where T
  #B = [convert(T, val)]
  #unsafe_fill!(pointer(A), pointer(B), length(A))
  Mem.set!(pointer(A), convert(T, val), length(A))
  A
end

## views

@inline function Base.view(A::VEArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(Array, j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(VEArray, j), J)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(J_gpu...)), J_gpu...)
end

#device(a::SubArray) = device(parent(a))
#context(a::SubArray) = context(parent(a))

# contiguous subarrays
function Base.unsafe_convert(::Type{VEPtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Base.RangeIndex}}}) where {T,N,P}
    return Base.unsafe_convert(VEPtr{T}, parent(V)) +
           Base._memory_offset(V.parent, map(first, V.indices)...)
end

# reshaped subarrays
function Base.unsafe_convert(::Type{VEPtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Union{Base.RangeIndex,Base.ReshapedUnitRange}}}}) where {T,N,P}
   return Base.unsafe_convert(VEPtr{T}, parent(V)) +
          (Base.first_index(V)-1)*sizeof(T)
end

## reshape

#device(a::Base.ReshapedArray) = device(parent(a))
#context(a::Base.ReshapedArray) = context(parent(a))

Base.unsafe_convert(::Type{VEPtr{T}}, a::Base.ReshapedArray{T}) where {T} =
  Base.unsafe_convert(VEPtr{T}, parent(a))

## reinterpret

#device(a::Base.ReinterpretArray) = device(parent(a))
#context(a::Base.ReinterpretArray) = context(parent(a))

Base.unsafe_convert(::Type{VEPtr{T}}, a::Base.ReinterpretArray{T,N,S} where N) where {T,S} =
  VPEtr{T}(Base.unsafe_convert(VEPtr{S}, parent(a)))


text:
function deparse(ex::Expr)
    if ex.head != :call
        return string(ex)
    else
        if contains([:+, :-, :*, :/, :^], ex.args[1])
            if length(ex.args) == 2
                return strcat(ex.args[1], deparse(ex.args[2]))
            else
                return join(map(x -> deparse(x), ex.args[2:end]), strcat(" ", string(ex.args[1]), " "))
            end
        else
            return strcat(ex.args[1], "(", join(map(x -> deparse(x), ex.args[2:end]), ", "), ")")
        end
    end
end
deparse(other::Any) = string(other)

# TODO: Examine string contents of inputs, insert parentheses if added:
# + (CONTAINS * OR /)
# - (CONTAINS * OR /)


text:
printstyled("  time\n", color=:light_green)

# Timestamp
t0 = time()
ts0 = String(split(string(u2d(t0)), '.')[1])
ts1 = String(split(timestamp(t0), '.')[1])
ts2 = String(split(timestamp(u2d(t0)), '.')[1])
ts3 = String(split(timestamp(string(u2d(t0))), '.')[1])
ts4 = String(split(timestamp(DateTime(string(u2d(t0)))), '.')[1])
@test ts1 == ts2 == ts3 == ts4 == ts0

# j2md
@test j2md(1, 1) == (1,1)
@test j2md(2000, 1) == (1,1)
@test j2md(2000, 60) == (2,29)
@test j2md(2000, 366) == (12,31)
@test_throws BoundsError j2md(2000, 367)
@test j2md(2012, 366) == (12,31)
@test j2md(2013, 60) == (3,1)
@test j2md(2015, 60) == (3,1)
@test j2md(2015, 365) == (12, 31)
@test j2md(2016, 61) == (3,1)
@test_throws BoundsError j2md(2013, 366)
@test_throws BoundsError j2md(2015, 366)
@test j2md(2100, 60) == (3,1)

# md2j
@test md2j(2000,1,1) == 1
@test md2j(2000,2,29) == 60
@test md2j(2013,3,1) == 60
@test md2j(2100,3,1) == 60
@test md2j(2001,1,1) == 1
@test 60 == md2j(2015,3,1)
@test 61 == md2j(2016,3,1)
@test 365 == md2j(2015,12,31)
@test 365 == md2j(1900,12,31)

# parsetimewin
t1 = t0 + 86400.0
dt0 = u2d(t0)
dt1 = u2d(t1)
st0 = string(dt0)
st1 = string(dt1)

d0, d1 = parsetimewin(dt0, dt1);  @test d1 > d0                      # dt, dt
d0, d1 = parsetimewin(dt0, t1);   @test d1 > d0                      # dt, r
d0, d1 = parsetimewin(dt0, st1);  @test d1 > d0                      # dt, s
d0, d1 = parsetimewin(t0, dt1);   @test d1 > d0                      # r, dt
d0, d1 = parsetimewin(t0, t1);    @test d1 > d0                      # r, r
d0, d1 = parsetimewin(t0, st1);   @test d1 > d0                      # r, s
d0, d1 = parsetimewin(st0, dt1);  @test d1 > d0                      # s, dt
d0, d1 = parsetimewin(st0, t1);   @test d1 > d0                      # s, r
d0, d1 = parsetimewin(st0, st1);  @test d1 > d0                      # s, s

# Checking that parsetimewin sorts correctly; type mismatch intentional
d0, d1 = parsetimewin(0, -600);  @test d1 > d0
d0, d1 = parsetimewin(-600.0, 0);  @test d1 > d0
d0, d1 = parsetimewin(600, 0.0);  @test d1 > d0
d0, d1 = parsetimewin(0.0, 600.0);  @test d1 > d0

@test ≈(600000, (DateTime(d1)-DateTime(d0)).value)
d0, d1 = parsetimewin("2016-02-29T23:30:00", "2016-03-01T00:30:00")
@test ≈(3600000, (DateTime(d1)-DateTime(d0)).value)
t = DateTime(now())
s = t-Hour(2)
d0, d1 = parsetimewin(s, t)
@test ≈(7200000, (DateTime(d1)-DateTime(d0)).value)

# t_collapse, t_expand
T = Int64[1 1451606400000000; 100001 30000000; 250001 12330000; 352303 99000000; 360001 0]
fs = 100.0
Δ = round(Int64, sμ/fs)
t_long = t_expand(T, fs)
@test ≈(T, t_collapse(t_long, fs))

T1 = hcat(cumsum(ones(Int64,size(T,1))), cumsum(T[:,2]))
fs1 = 0.0
@test ≈(T1, t_collapse(t_expand(T1, fs1), fs1))

# endtime
printstyled(stdout, "    endtime\n", color=:light_green)
@test endtime(T, fs) == last(t_long)
@test endtime(T, Δ) == endtime(T, fs)
@test endtime(Array{Int64,2}(undef, 0, 0), Δ) == 0
@test endtime(Array{Int64,2}(undef, 0, 0), 100Δ) == 0
@test endtime(Array{Int64,2}(undef, 0, 0), fs) == 0

printstyled(stdout, "    t_win, w_time\n", color=:light_green)
printstyled(stdout, "      faithful representation of gaps\n", color=:light_green)
fs = 100.0
Δ = round(Int64, sμ/fs)
t = [1 0; 6 980000; 8 100000; 10 0]
# t1 = t_expand(t, fs)
 #       0
 #   10000
 #   20000
 #   30000
 #   40000
 # 1030000
 # 1040000
 # 1150000
 # 1160000
 # 1170000

t2 = t_win(t, Δ)
 # first      last
 # sample     sample
 # time       time
 #       0    40000 ===> so that collect(w[i,1]:Δ:w[i,2]) for each window == t_expand
 # 1030000  1040000
 # 1150000  1170000
 #
@test t2[1,2] == 40000
@test t2[2,:] == [1030000, 1040000]
@test t2[3,:] == [1150000, 1170000]

printstyled(stdout, "      arbitrary windows containing gaps\n", color=:light_green)
t = [1 999999998990000; 101 10000; 297 0]
@test w_time(t_win(t, Δ), Δ) == t

t = [1 999999998990000; 101 10000; 297 1000000; 303 40000; 500 1000000; 10000 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      length-0 gap\n", color=:light_green)
t = [1 0; 6 0; 8 0; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      negative gap\n", color=:light_green)
t = [1 0; 6 -2Δ; 8 -10Δ; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      single-point gap\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      non-null gap at end\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 Δ]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      negative gap at end\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 -5Δ]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "    mktime\n", color=:light_green)
fv = 0.0005
iv = Array{Int32,1}([1980, 082, 10, 35, 39, 890])
(m,d) = j2md(iv[1], iv[2])

ts_0 = round(Int64, d2u(DateTime(iv[1], m, d, iv[3], iv[4], iv[5]))*sμ) +
       iv[6]*1000 +
       round(Int64,fv*1000.0)
ts_1 = Date(iv[1], m, d).instant.periods.value * 86400000000 +
       div(Time(iv[3], iv[4], iv[5]).instant.value, 1000) +
       iv[6]*1000 +
       round(Int64,fv*1000.0) -
       SeisIO.dtconst
ts_2 = mktime(iv[1], iv[2], iv[3], iv[4], iv[5], iv[6]*Int32(1000)) +
       round(Int64, fv*1000.0)
iv[6] *= Int32(1000)
ts_3 = mktime(iv) + round(Int64, fv*1000.0)
@test ts_0 == ts_1 == ts_2 == ts_3
# timespec()

printstyled(stdout, "    int2tstr, tstr2int\n", color=:light_green)
s = "2018-01-01T00:00:00.000001"
t = "2018-01-04T00:00:00.003900"
si = tstr2int(s)
ti = tstr2int(t)
j = loop_time(si, ti)

t = "2018-01-04T00:00:00.39"
for (n,s) in enumerate(["2018-01-01T00:00:00.000001",
                        "2018-01-01T00:00:00",
                        "2018-01-01T00:00:00.035",
                        "2016-02-29T00:00:00.02",
                        "2018-02-28T00:00:00.33"])

  s_str = identity(s)
  if length(s) == 19
    s_str *= "."
  end
  s_str = rpad(s_str, 26, '0')
  @test (int2tstr(tstr2int(s))) == s_str
  if n < 4
    @test loop_time(tstr2int(s), tstr2int(t)) == 4
    @test loop_time(tstr2int(s), tstr2int(t), ti=43200000000) == 7
    @test loop_time(tstr2int(s), tstr2int(t), ti=3600000000) == 73
  elseif n == 4
    @test loop_time(tstr2int(s), tstr2int("2016-03-02T00:00:02")) == 3
  elseif n == 5
    @test loop_time(tstr2int(s), tstr2int("2018-03-02T00:00:02")) == 3
  end
end


text:
using EzXML

# struct for storing xml data
mutable struct Voc_xml
	folder::Array{String}
	filename::Array{String}
	width::Array{Int64}
	height::Array{Int64}
	depth::Array{Int64}
	segmented::Array{Int64}
	name::Array{String}
	pose::Array{String}
	truncated::Array{Int64}
	difficult::Array{Int64}
	xmin::Array{Int64}
	ymin::Array{Int64}
	xmax::Array{Int64}
	ymax::Array{Int64}
	Voc_xml() = new([],[],[],[],[],[],[],[],[],[],[],[],[],[])
end

function get_xml_meta(children)
	folder 	  = children[1].content
	filename  = children[2].content
	width 	  = elements(children[5])[1].content
	height 	  = elements(children[5])[1].content
	depth 	  = elements(children[5])[3].content
	segmented = children[6].content
	return (folder, filename, width, height, depth, segmented)
end

function upd_xml_data!(voc_xml_data, voc_xml_meta, children)
	for i = 7:length(children) # first 6 fields are meta
		push!(voc_xml_data.folder, voc_xml_meta[1])
		push!(voc_xml_data.filename, voc_xml_meta[2])
		push!(voc_xml_data.width, parse(Int64,voc_xml_meta[3]))
		push!(voc_xml_data.height, parse(Int64,voc_xml_meta[4]))
		push!(voc_xml_data.depth, parse(Int64,voc_xml_meta[5]))
		push!(voc_xml_data.segmented, parse(Int64,voc_xml_meta[6]))
		push!(voc_xml_data.name, elements(children[i])[1].content)
		push!(voc_xml_data.pose, elements(children[i])[2].content)
		push!(voc_xml_data.truncated, parse(Int64,elements(children[i])[3].content))
		push!(voc_xml_data.difficult, parse(Int64,elements(children[i])[4].content))
		if elements(children[i])[5].name == "bndbox"# skip all the person: head, hand etc....
			push!(voc_xml_data.xmin, parse(Int64,elements(elements(children[i])[5])[1].content))
			push!(voc_xml_data.ymin, parse(Int64,elements(elements(children[i])[5])[2].content))
			push!(voc_xml_data.xmax, parse(Int64,elements(elements(children[i])[5])[3].content))
			push!(voc_xml_data.ymax, parse(Int64,elements(elements(children[i])[5])[4].content))
		end
	end
	return nothing
end

function upd_xml_data!(voc_xml_data, voc_xml_path)
	for (root, dirs, files) in walkdir(voc_xml_path)
		for file in files
			xml_raw = String(read(voc_xml_path*"\\"*file))
			# Parse an XML string
			doc = parsexml(xml_raw)
			# Get the root element from `doc`.
			rt = doc.root
			children  = elements(rt)
			voc_xml_meta = get_xml_meta(children)
			upd_xml_data!(voc_xml_data, voc_xml_meta, children)
		end
	end
return nothing
end

voc_xml_path = raw"C:\Projects\TrafficSign\DataSet\VOCdevkit\VOC2007\Annotations";
# init empty struct
voc_xml_data = Voc_xml()
@time upd_xml_data!(voc_xml_data, voc_xml_path)


text:
# Undergraduate Student: Arturo Burgos
# Professor: João Rodrigo Andrade
# Federal University of Uberlândia - UFU, Fluid Mechanics Laboratory - MFLab, Block 5P, Uberlândia, MG, Brazil

using Printf # Just the same as import a module
using Statistics
using CPUTime

println("Choose a N value: ")
N = chomp(readline()) # Here I get the char from the user
n = parse(Int32,N)  # Here I convert the char into an integer(32bits)

A = rand(1:200,n,n);
B = rand(1:200,n,n);

D = A*B; # This is the vectorized way for this problem, like MATLAB
C = A.*B; # This multiplies point by point

R = zeros(n,n) # Here I create a zero matrix just to get the multiply of A and B

# println("\n\nThe D Matrix is: \n\n",D) # Compare with low n values

# i represents lines and j represents columns. k represents
# the line or column 11element11 indices

function multiloop(mtx1,mtx2,mtx3)

    for i = 1:n
        for j = 1:n
            soma = 0;
            for k = 1:n
                soma = soma + mtx1[i,k]*mtx2[k,j];
                mtx3[i,j] = soma;
            end
        end
    end
    return mtx3 # Returns a value to the variable
end

R = multiloop(A,B,R);
# println("\n\nThe R Matrix is: \n\n",R) # Compare with low n values

println("\n")

@time @CPUtime multiloop(A,B,R)
# Test done!

max_stars_repo_name: arturofburgos/Comparing-Languages

max_issues_repo_name: arturofburgos/Comparing-Languages

max_forks_repo_name: arturofburgos/Comparing-Languages


text:
function SimpleSDMPredictor(::Type{WorldClim}, ::Type{Elevation}, layer::Integer=1; resolution::Float64=10.0, kwargs...)
    @assert resolution in [0.5, 2.5, 5.0, 10.0]
    file = _get_raster(WorldClim, Elevation, layer, resolution)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end


text:
### Turing interface

import .Turing
import .Turing.DynamicPPL as DynPPL
import ComponentArrays: ComponentVector

export TuringMuseProblem

struct TuringMuseProblem{A<:AD.AbstractBackend, M<:Turing.Model} <: AbstractMuseProblem

    autodiff :: A
    model :: M
    model_for_prior
    vi_z′_θ
    vi_z′_θ′
    vi_θ
    vi_θ′
    x
    observed_vars
    latent_vars
    hyper_vars

end

@doc doc"""

    TuringMuseProblem(model; params, autodiff = Turing.ADBACKEND)

Specify a MUSE problem with a [Turing](https://turing.ml) model.

The Turing model should be conditioned on the variables which comprise
the "data", and all other variables should be unconditioned. By
default, any parameter which doesn't depend on another parameter will
be estimated by MUSE, but this can be overridden by passing `params`
as a list of symbols. All other non-conditioned and non-`params`
variables will be considered the latent space.

The `autodiff` parameter should be either
`MuseInference.ForwardDiffBackend()` or
`MuseInference.ZygoteBackend()`, specifying which library to use for
automatic differenation. The default uses whatever the global
`Turing.ADBACKEND` is currently set to.

## Example

```julia
# toy hierarchical model
Turing.@model function toy()
    σ ~ Uniform(0, 1)
    θ ~ Normal(0, σ)
    z ~ MvNormal(zeros(512), exp(σ/2)*I)
    w ~ MvNormal(z, I)
    x ~ MvNormal(w, I)
    y ~ MvNormal(x, I)
    (;σ, θ, z, w, x, y)
end
sim = toy()()
model = toy() | (;sim.x, sim.y)
prob = TuringMuseProblem(model, params=(:σ, :θ))

# get solution
result = muse(prob, (σ=0.5, θ=0))
```

Here we have chosen `(σ, θ)` to be the parameters which MuseInference
will estimate (note that the default would have only chosen `σ`). The
observed data are `(x,y)` and the remaining `(z,w)` are the latent
space.

!!! note

    You can also call [`muse`](@ref), etc... directly on the model, e.g.
    `muse(model, (σ=0.5, θ=0))`, in which case the parameter names `params`
    will be read from the keys of provided the starting point.

!!! note

    The model function cannot have any of the random variables as
    arguments, although it can have other parameters as arguments. E.g.,

    ```julia
    # not OK
    @model function toy(x)
        x ~ Normal()
        ...
    end

    # OK
    @model function toy(σ)
        x ~ Normal(σ)
        ...
    end
    ```

"""
function TuringMuseProblem(
    model;
    params = (:θ,),
    autodiff = nothing,
)

    # set backend based on Turing's by default
    if autodiff == nothing
        if Turing.ADBACKEND[] == :zygote
            autodiff = ZygoteBackend()
        elseif Turing.ADBACKEND[] == :forwarddiff
            autodiff = ForwardDiffBackend()
        else
            error("Unsupposed backend from Turing: $(Turing.ADBACKEND)")
        end
    end
    # prevent this constructor from advancing the default RNG for more clear reproducibility
    rng = copy(Random.default_rng())
    # model is expected to be passed in conditioned on x
    x = ComponentVector(model.context.values)
    # figure out variable names
    observed = keys(x)
    latent = keys(delete(_namedtuple(DynPPL.VarInfo(rng, model)), (observed..., params...)))
    # VarInfo for (z,θ) with both transformed
    vi_z′_θ′ = DynPPL.VarInfo(rng, model)
    DynPPL.settrans!.((vi_z′_θ′,), true, _VarName.((latent..., params...)))
    # VarInfo for (z,θ) with only z transformed
    vi_z′_θ = DynPPL.VarInfo(rng, model)
    DynPPL.settrans!.((vi_z′_θ,), true, _VarName.(latent))
    # model with all vars free
    model = DynPPL.decondition(model)
    # model for computing prior, just need any values for (x,z) to condition on here
    vars = _namedtuple(DynPPL.evaluate!!(model, rng)[2])
    model_for_prior = model | select(vars, (observed..., latent...))
    # VarInfo for θ
    vi_θ = DynPPL.VarInfo(rng, model_for_prior)
    # VarInfo for transformed θ
    vi_θ′ = deepcopy(vi_θ)
    DynPPL.settrans!.((vi_θ′,), true, _VarName.(params))

    TuringMuseProblem(
        autodiff,
        model,
        model_for_prior,
        vi_z′_θ,
        vi_z′_θ′,
        vi_θ,
        vi_θ′,
        x,
        observed,
        latent,
        params
    )

end

function transform_θ(prob::TuringMuseProblem, θ)
    vi = deepcopy(prob.vi_θ)
    DynPPL.setval!(vi, θ)
    DynPPL.link!(vi, DynPPL.SampleFromPrior())
    ComponentVector(vi)
end

function inv_transform_θ(prob::TuringMuseProblem, θ)
    vi = deepcopy(prob.vi_θ)
    DynPPL.setval!(vi, θ)
    for k in keys(θ)
        DynPPL.settrans!(vi, true, _VarName(k))
    end
    DynPPL.invlink!(vi, DynPPL.SampleFromPrior())
    ComponentVector(vi)
end

standardizeθ(prob::TuringMuseProblem, θ::NamedTuple) = 1f0 * ComponentVector(θ) # ensure at least Float32
standardizeθ(prob::TuringMuseProblem, θ::Number) =
    length(prob.hyper_vars) == 1 ? standardizeθ(prob, (;θ)) : error("Invalid θ type for this problem.")

function logPriorθ(prob::TuringMuseProblem, θ, θ_space)
    vi = is_transformed(θ_space) ? prob.vi_θ′ : prob.vi_θ
    DynPPL.logprior(prob.model_for_prior, DynPPL.VarInfo(vi, θ))
end

function ∇θ_logLike(prob::TuringMuseProblem, x, z, θ, θ_space)
    model = DynPPL.condition(prob.model, x)
    vi = is_transformed(θ_space) ? prob.vi_z′_θ′ : prob.vi_z′_θ
    first(AD.gradient(prob.autodiff, θ -> DynPPL.logjoint(model, DynPPL.VarInfo(vi, z, θ)), θ))
end

function ẑ_at_θ(prob::TuringMuseProblem, x, z₀, θ; ∇z_logLike_atol)
    model = DynPPL.condition(prob.model, x)
    neglogp(z) = -DynPPL.logjoint(model, DynPPL.VarInfo(prob.vi_z′_θ, z, θ))
    soln = Optim.optimize(optim_only_fg!(neglogp, prob.autodiff), z₀, Optim.LBFGS(), Optim.Options(g_tol=∇z_logLike_atol))
    _check_optim_soln(soln)
    soln.minimizer, soln
end

function sample_x_z(prob::TuringMuseProblem, rng::AbstractRNG, θ)
    model = DynPPL.condition(prob.model, θ)
    vi = DynPPL.VarInfo(rng, model)
    vars_untransformed = map(copy, _namedtuple(vi))
    DynPPL.link!(vi, DynPPL.SampleFromPrior())
    vars_transformed = map(copy, _namedtuple(vi))
    (;
        x = ComponentVector(select(vars_untransformed, prob.observed_vars)),
        z = ComponentVector(select(vars_transformed,   prob.latent_vars))
    )
end

# helped to extract parameters from a sampled model. feels like there
# should be a less hacky way to do this...
function _namedtuple(vi::DynPPL.VarInfo)
    map(DynPPL.TypedVarInfo(vi).metadata) do m
        if m.vns[1] isa DynPPL.VarName{<:Any,Setfield.IdentityLens} && length(m.vals)==1
            m.vals[1]
        else
            m.vals
        end
    end
end

ComponentVector(vi::DynPPL.VarInfo) = ComponentVector(_namedtuple(vi))

function DynPPL.VarInfo(vi::DynPPL.TypedVarInfo, x::Union{NamedTuple,ComponentVector}, xs::Union{NamedTuple,ComponentVector}...)
    DynPPL.VarInfo(vi, merge(map(_namedtuple, (x, xs...))...))
end

function DynPPL.VarInfo(vi::DynPPL.TypedVarInfo, x::NamedTuple)
    T = promote_type(map(eltype, values(x))..., map(eltype, _namedtuple(values(vi)))...) # if x is ForwardDiff Duals
    DynPPL.VarInfo(
        NamedTuple{keys(vi.metadata)}(map(keys(vi.metadata),values(vi.metadata)) do k,v
            DynPPL.Metadata(
                v.idcs,
                v.vns,
                v.ranges,
                atleast1d(getfield(x,k)),
                v.dists,
                v.gids,
                v.orders,
                v.flags,
            )
        end),
        Base.RefValue{T}(DynPPL.getlogp(vi)),
        vi.num_produce
    )
end

DynPPL.condition(model::DynPPL.Model, x::ComponentVector) = DynPPL.condition(model, _namedtuple(x))
_VarName(x::Symbol) = DynPPL.VarName{x}()

atleast1d(x::Number) = [x]
atleast1d(x::AbstractVector) = x

_params_from_θ₀(θ₀::Number) = (:θ,)
_params_from_θ₀(θ₀) = keys(θ₀)

function muse!(result::MuseResult, model::DynPPL.Model, θ₀ = result.θ; kwargs...)
    muse!(result, TuringMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end
function get_J!(result::MuseResult, model::DynPPL.Model, θ₀ = result.θ; kwargs...)
    get_J!(result, TuringMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end
function get_H!(result::MuseResult, model::DynPPL.Model, θ₀ = result.θ; kwargs...)
    get_H!(result, TuringMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end


text:
module ModelParameters

import Base.show
import Base.getproperty

include("Types.jl")
include("Set.jl")
include("GetBounds.jl")
include("Show.jl")

export AbstractParameter,
    AbstractParameterSet,
    ParameterGroup,
    ParameterSpace,
    Parameter,
    BoundedParameter,
    getproperty,
    description,
    get_bounds,
    get_initial_guess,
    get_lower_bounds,
    get_upper_bounds,
    set,
    ←,
    show

end


text:
@testset "Bridge" begin
    gint = SimpleGraph(13)
    add_edge!(gint, 1, 7)
    add_edge!(gint, 1, 2)
    add_edge!(gint, 1, 3)
    add_edge!(gint, 12, 13)
    add_edge!(gint, 10, 13)
    add_edge!(gint, 10, 12)
    add_edge!(gint, 12, 11)
    add_edge!(gint, 5, 4)
    add_edge!(gint, 6, 4)
    add_edge!(gint, 8, 9)
    add_edge!(gint, 6, 5)
    add_edge!(gint, 1, 6)
    add_edge!(gint, 7, 5)
    add_edge!(gint, 7, 3)
    add_edge!(gint, 7, 8)
    add_edge!(gint, 7, 10)
    add_edge!(gint, 7, 12)

    for g in testgraphs(gint)
        brd = @inferred(bridges(g))
        ans = [
            Edge(1, 2),
            Edge(8, 9),
            Edge(7, 8),
            Edge(11, 12),
        ]
        @test brd == ans
    end
    for level in 1:6
        btree = LightGraphs.binary_tree(level)
        for tree in [btree, Graph{UInt8}(btree), Graph{Int16}(btree)]
            brd = @inferred(bridges(tree))
            ans = collect(edges(tree))
            @test Set(brd) == Set(ans)
        end
    end

    hint = blockdiag(wheel_graph(5), wheel_graph(5))
    add_edge!(hint, 5, 6)
    for h in (hint, Graph{UInt8}(hint), Graph{Int16}(hint))
        @test @inferred(bridges(h)) == [
            Edge(5, 6),
        ]
    end

    dir = SimpleDiGraph(10, 10)
    @test_throws MethodError bridges(dir)
end


text:
using Random

if scene_mode == "synthetic"
    sim_rng_seed = prediction_rng_seed;
    sim_rng = MersenneTwister(sim_rng_seed);
    scene_param = SyntheticSceneParameter(sim_rng);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = BICControlParameter(u_norm_max, tcalc, dtr, min_dist);

elseif scene_mode == "data"
    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = BICControlParameter(u_norm_max, tcalc, dtr, min_dist);

else
    @error @error "scene_mode: $(scene_mode) is not supported!"
end


text:
# run tests related to parallel functionality that are not actually parallel
"""
  Run serial rk4 and newton cases, and generate input files for the
  equivalent parallel cases
"""
function test_parallel()
  start_dir = pwd()
  fname = "input_vals_vortex3.jl"

  # rk4
  cd("./rk4/serial")
  mesh, sbp, eqn, opts = solvePDE(fname)
  cd("../parallel")

  opts["smb_name"] = "SRCMESHES/psquare2.smb"
  make_input(opts, "input_vals_parallel")

  # staggered grid
  cd("../staggered_serial")
  mesh, sbp, eqn, opts = solvePDE(fname)

  cd("../staggered_parallel")
  opts["smb_name"] = "SRCMESHES/psquare2.smb"
  make_input(opts, "input_vals_parallel")

  cd("../../")
  # lserk
  if !isdir("./lserk")
    mkdir("./lserk")
  end
  if !isdir("./lserk/serial")
    mkdir("./lserk/serial")
  end
  if !isdir("lserk/parallel")
    mkdir("./lserk/parallel")
  end

  opts = Input.read_input_file("./rk4/serial/input_vals_vortex3.jl")
  cd("./lserk/serial")
  opts["run_type"] = 30
  make_input(opts, "input_vals_vortex3")
  mesh, sbp, eqn,  opts = solvePDE(fname)

  cd("../parallel")
  opts["smb_name"] = "SRCMESHES/psquare2.smb"
  make_input(opts, "input_vals_parallel")

  # newton
  cd("../../newton/serial")
  fname = "input_vals_vortex3.jl"
  mesh, sbp, eqn, opts = solvePDE(fname)

  cd("../parallel")
  opts["smb_name"] = "SRCMESHES/psquare2.smb"
  make_input(opts, "input_vals_parallel")

  cd(start_dir)

  return nothing
end

#test_parallel()
add_func1!(EulerTests, test_parallel, [TAG_PARALLEL, TAG_SHORTTEST])


text:
"""
    VariationalConv(filter, in => out, σ = identity;
                    stride = 1, pad = 0, dilation = 1,
                    groups = 1,
                    weight_dist = TrainableMvNormal,
                    bias_dist = TrainableMvNormal,
                    [bias, weight, init])
    VariationalConvBE(σ, weight_sampler, bias_sampler,
                    stride, pad, dilation, groups)

Creates a variational conv layer. Computes variational bayesian
approximation to the distribution over the parameters of the conv layer.
The stochasticity is during the forward pass, instead of using point
estimates for weights and biases, we sample from the distribution over
weights and biases. Gradients of the distribution's learnable parameters
are trained using the reparameterization trick.

Reference - https://arxiv.org/abs/1505.05424
We use DistributionsAD - https://github.com/TuringLang/DistributionsAD.jl
to help us with backprop.

# Fields
- `σ`: Activation function, applies to logits after layer transformation
- `weight_sampler`: A trainable distribution from which weights are sampled
                    in every forward pass
- `bias_sampler`: A trainable distribution from which biases are sampled in
                    every forward pass
- `stride`: Convolution stride
- `pad`
- `dilation`
- `groups`

# Arguments
- `filter::NTuple{N,Integer}`: Kernel dimensions, eg, (5, 5)
- `ch::Pair{<:Integer,<:Integer}`: Input channels => output channels
- `σ::F=identity`: Activation of the dense layer, defaults to identity
- `weight_dist=TrainableMvNormal`: Initialization function for weights.
- `bias_dist=TrainableMvNormal`: Initialization function for biases.

"""
struct VariationalConv{N,M,F,A,V}
    σ::F
    weight_sampler::A
    bias_sampler::V
    stride::NTuple{N,Int}
    pad::NTuple{M,Int}
    dilation::NTuple{N,Int}
    groups::Int
end

function VariationalConv(
    k::NTuple{N,Integer},
    ch::Pair{<:Integer,<:Integer},
    σ = identity;
    stride = 1,
    pad = 0,
    dilation = 1,
    groups = 1,
    init = glorot_normal,
    weight_dist = TrainableMvNormal,
    bias_dist = TrainableMvNormal,
    weight = Flux.convfilter(k, (ch[1] ÷ groups => ch[2])),
    bias = true,
    device = cpu,
) where {N}

    stride = Flux.expand(Val(N), stride)
    dilation = Flux.expand(Val(N), dilation)
    pad = Flux.calc_padding(VariationalConv, pad, size(weight)[1:N], dilation, stride)
    # Distribution from which weights are sampled
    weight_sampler = weight_dist((k..., ch...), init = init, device = device)
    # Distribution from which biases are sampled
    if bias
        bias = create_bias(weight, bias, size(weight, N + 2))
        bias_sampler = bias_dist(size(bias), init = init, device = device)
    else
        bias_sampler = () -> Flux.Zeros()
    end
    return VariationalConv(σ, weight_sampler, bias_sampler, stride, pad, dilation, groups)
end

@functor VariationalConv (weight_sampler, bias_sampler)

function (c::VariationalConv)(x)
    # Sample weights and biases
    weight = c.weight_sampler()
    bias = c.bias_sampler()
    σ, b = c.σ, reshape(bias, ntuple(_ -> 1, length(c.stride))..., :, 1)
    output = Flux.NNlib.conv(
        x,
        weight;
        stride = c.stride,
        pad = c.pad,
        dilation = c.dilation,
        groups = c.groups,
    )
    σ.(output .+ b)
end


text:
# Autogenerated wrapper script for OpenCV_jll for i686-w64-mingw32-cxx03
export libopencv_calib3d, libopencv_core, libopencv_dnn, libopencv_features2d, libopencv_flann, libopencv_gapi, libopencv_highgui, libopencv_imgcodecs, libopencv_imgproc, libopencv_julia, libopencv_objdetect, libopencv_stitching, libopencv_video, libopencv_videoio

using Qt5Base_jll
using Libglvnd_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("OpenCV")
JLLWrappers.@declare_library_product(libopencv_calib3d, "libopencv_calib3d452.dll")
JLLWrappers.@declare_library_product(libopencv_core, "libopencv_core452.dll")
JLLWrappers.@declare_library_product(libopencv_dnn, "libopencv_dnn452.dll")
JLLWrappers.@declare_library_product(libopencv_features2d, "libopencv_features2d452.dll")
JLLWrappers.@declare_library_product(libopencv_flann, "libopencv_flann452.dll")
JLLWrappers.@declare_library_product(libopencv_gapi, "libopencv_gapi452.dll")
JLLWrappers.@declare_library_product(libopencv_highgui, "libopencv_highgui452.dll")
JLLWrappers.@declare_library_product(libopencv_imgcodecs, "libopencv_imgcodecs452.dll")
JLLWrappers.@declare_library_product(libopencv_imgproc, "libopencv_imgproc452.dll")
JLLWrappers.@declare_library_product(libopencv_julia, "libopencv_julia.dll")
JLLWrappers.@declare_library_product(libopencv_objdetect, "libopencv_objdetect452.dll")
JLLWrappers.@declare_library_product(libopencv_stitching, "libopencv_stitching452.dll")
JLLWrappers.@declare_library_product(libopencv_video, "libopencv_video452.dll")
JLLWrappers.@declare_library_product(libopencv_videoio, "libopencv_videoio452.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll, Libglvnd_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libopencv_calib3d,
        "bin\\libopencv_calib3d452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_core,
        "bin\\libopencv_core452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_dnn,
        "bin\\libopencv_dnn452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_features2d,
        "bin\\libopencv_features2d452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_flann,
        "bin\\libopencv_flann452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_gapi,
        "bin\\libopencv_gapi452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_highgui,
        "bin\\libopencv_highgui452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgcodecs,
        "bin\\libopencv_imgcodecs452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgproc,
        "bin\\libopencv_imgproc452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_julia,
        "bin\\libopencv_julia.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_objdetect,
        "bin\\libopencv_objdetect452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_stitching,
        "bin\\libopencv_stitching452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_video,
        "bin\\libopencv_video452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_videoio,
        "bin\\libopencv_videoio452.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function fuel(mass)::Int
    fuelreq = round(mass/3, RoundDown) - 2
    if fuelreq <= 0
        return 0
    end
    return fuelreq + fuel(fuelreq)
end

masses = open("mass.csv") do f
    parse.(Int, readlines(f))
end
println(sum(round.(masses/3, RoundDown) .- 2))
println(sum(fuel.(masses)))


text:
import RowEchelon

#=
Homework 6
Section 2.5
Problems 2
=#

# Create the matrix and the LU matrices
A = [2 -6 4 ; -4 8 0 ; 0 -4 6]
Ab = [2 -6 4 -7; -4 8 0 5; 0 -4 6 2]
L = [1 0 0 ; -2 1 0; 0 1 1]
U = [2 -6 4 ; 0 -4 8 ; 0 0 -2]

# Row Reduce A
Arref = RowEchelon.rref(Ab)
print("ROW REDUCED SOLUTION: \t")
println(Arref)

# Find Ly = b
# Find Ux = y
# Append b to L
Lb = [1 0 0 -7; -2 1 0 5; 0 1 1 2]
Lbrref = RowEchelon.rref(Lb)

# Append y to U
Uy = [2 -6 4 Lbrref[1,4]; 0 -4 8 Lbrref[2,4]; 0 0 -2 Lbrref[3,4]]
Uyrref = RowEchelon.rref(Uy)
print("LU SOLUTION: \t \t")
println(Uyrref)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using Test, AsterReader
using AsterReader: aster_parse_nodes, aster_read_mesh, MEDFile, get_element_sets

datadir = first(splitext(basename(@__FILE__)))

@testset "test med file exported from gmsh, issue #16" begin
    mesh_file = joinpath(datadir, "t1.med")
    mesh = aster_read_mesh(mesh_file)
end


text:
function murmur64(h::UInt64)
	h ⊻= h >> 33
	h *= 0xff51afd7ed558ccd
	h ⊻= h >> 33
	h *= 0xc4ceb9fe1a85ec53
	h ⊻= h >> 33
	return h
end

# returns random number, modifies the seed
function splitmix64(seed::UInt64)
	seed += 0x9e3779b97f4a7c15
	z = (seed ⊻ (seed >> 30)) * 0xbf58476d1ce4E5b9
	z = (z ⊻ (z >> 27)) * 0x94d049bb133111eb
	return seed, z ⊻ (z >> 31)
end

function mixsplit(key::UInt64, seed::UInt64)
	return murmur64(key + seed)
end

function rotl64(n::UInt64, c::Int)
	return (n << UInt64(c & 63)) | (n >> UInt64((-c) & 63))
end

function reduce(hash, n::UInt32)
	# http://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-modulo-reduction/
	return UInt32((UInt64(hash) * UInt64(n)) >> 32)
end

# xor fold
@inline function fingerprint(::Type{UInt8}, hash::UInt64)
	hash = (hash >> 32) ⊻ ((hash << 32) >> 32) # fold to 32 bits
	hash = (hash >> 16) ⊻ ((hash << 48) >> 48) # fold to 16 bits
	hash = (hash >> 8) ⊻ ((hash << 56) >> 56) # fold to 8 bits
	return UInt8(hash)
end

@inline function fingerprint(::Type{UInt16}, hash::UInt64)
	hash = (hash >> 32) ⊻ ((hash << 32) >> 32) # fold to 32 bits
	hash = (hash >> 16) ⊻ ((hash << 48) >> 48) # fold to 16 bits
	return UInt16(hash)
end

@inline function fingerprint(::Type{UInt32}, hash::UInt64)
	hash = (hash >> 32) ⊻ ((hash << 32) >> 32) # fold to 32 bits
	return UInt32(hash)
end

function mod3(x)
	if (x > 2)
	    x -= 3
    end
	return x
end

function binary_fuse_mulhi(a::UInt64, b::UInt64)
    return UInt64((UInt128(a) * UInt128(b)) >> 64)
end


text:
function nlfeast!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0, 0.0), r=1.0, debug=false, ϵ=10e-12, store=true, spurious=1e-5,
    factorizer=lu, left_divider=ldiv!)

    N, m₀ = size(X)
    Λ, res = zeros(ComplexF64, m₀), Array{Float64}(undef, m₀)
    θ = LinRange(π / nodes, 2 * π - π / nodes, nodes)
    Q₀, Q₁, R = similar(X, ComplexF64), similar(X, ComplexF64), similar(X, ComplexF64)
    A, B = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)

    qt, rt = qr!(X)
    X .= Matrix(qt)

    l = ReentrantLock()

    if store
        facts = Array{Factorization}(undef, nodes)
        Threads.@threads for i = 1:nodes
            z = (r * exp(θ[i] * im) + c)
            tempfact = factorizer(T(z))
            lock(l) do
                facts[i] = tempfact
            end
            if debug print("*") end
        end
        if debug println() end
    end

    for nit = 0:iter

        Q₀ .= 0
        Q₁ .= 0
        # l = ReentrantLock()

        Threads.@threads for i = 1:nodes
            z = (r * exp(θ[i] * im) + c)
            Tinv = similar(X, ComplexF64)
            if nit == 0
                if store
                    left_divider(Tinv, facts[i], X)
                    Tinv .*= (r * exp(θ[i] * im) / nodes)
                else
                    Tinv .= (T(z) \ X) .* (r * exp(θ[i] * im) / nodes)
                end
            else
                resolvent = (1 ./ (z .- Λ)) .* (r * exp(θ[i] * im) / nodes)
                if store
                    left_divider(Tinv, facts[i], R)
                    Tinv .= X - Tinv
                else
                    Tinv .= (X - T(z) \ R)
                end
                rmul!(Tinv,  Diagonal(resolvent))
            end
            lock(l) do
                Q₀ .+= Tinv
                Q₁ .+= Tinv .* z
            end
    		if debug print(".") end
        end
		if debug println() end

		beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)

        update_R!(X, R, Λ, T)
        res .= residuals(R, Λ, T)

        if debug
            iter_debug_print(nit, Λ, res, c, r, spurious)
        end

		res_inside = res[in_contour.(Λ, c, r)]
        if size(res_inside, 1) > 0 && maximum(res_inside) < ϵ
            break
        end
		if nit > 1 && sum(res_inside .< spurious) > 0 && maximum(res_inside[res_inside .< spurious]) < ϵ
			break
		end
    end

    normalize!(X)
    Λ, X, res
end

function nlfeast_it!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0, 0.0), r=1.0, debug=false, ϵ=0.05)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, m₀)
    θ = LinRange(π / nodes, 2 * π - π / nodes, nodes)
    Q₀, Q₁ = zeros(ComplexF64, N, m₀), zeros(ComplexF64, N, m₀)
    R = similar(X, ComplexF64)
    Temp = similar(X, ComplexF64)
    resolvent = similar(Λ)
    Tinv = rand(ComplexF64, nodes, N, m₀)
    A, B = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
    res = Array{Float64}(undef, size(R, 2))

    for i = 1:nodes
        z = (r * exp(θ[i] * im) + c)
        # Temp .= (T(z)\X) .* (r*exp(θ[i]*im)/nodes)
        for j = 1:m₀
            # Tinv[i,:,j], log = gmres(T(z), X[:,j]; tol=10e-8, log=true)
            Tinv[i,:,j] .= bicgstabl(T(z), X[:,j], 2; tol=1e-3)
        end
        Temp .= Tinv[i,:,:] .* (r * exp(θ[i] * im) / nodes)
        Q₀ .+= Temp
        Q₁ .+= Temp .* z
		if debug print(".") end
    end
	if debug println() end

    beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
    # beyn_qr_step!(Q₀, Q₁, X, Λ)

    update_R!(X, R, Λ, T)

    if (iter == 0)
        res .= residuals(R, Λ, T)
    end
    if debug
        iter_debug_print(0, Λ, res, c, r)
    end

    for nit = 1:iter

        Q₀ .= 0
        Q₁ .= 0

        for i = 1:nodes
            z = (r * exp(θ[i] * im) + c)
            resolvent .= (1 ./ (z .- Λ))
            # Tinv[i,:,:] .= T(z)\R
            for j = 1:m₀
                # Tinv[i,:,j], log = gmres!(Tinv[i,:,j], T(z), R[:,j]; log=true, maxiter=1000, tol=1e-3)
                # Tinv[i,:,j], log = gmres(T(z), R[:,j]; log=true, maxiter=1000, tol=1e-3)
                Tinv[i,:,j], log = bicgstabl!(Tinv[i,:,j], T(z), R[ :,j], 1; log=true, tol=1e-8)
                # Tinv[i,:,j], log = bicgstabl(T(z), R[:,j], 2; log=true, max_mv_products=1000, tol=1e-3)
                # if !log.isconverged println(log) end
                # print(log.iters)
                # print(":\t")
                # println(maximum(log[:resnorm]))
            end
            Temp .= X - Tinv[i,:,:]
            rmul!(Temp, (r * exp(θ[i] * im) / nodes) .* Diagonal(resolvent))
            Q₀ .+= Temp
            Q₁ .+= Temp .* z
    		if debug print(".") end
        end
		if debug println() end

        beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_qr_step!(Q₀, Q₁, X, Λ)

        update_R!(X, R, Λ, T)
        res .= residuals(R, Λ, T)

        if debug
            iter_debug_print(nit, Λ, res, c, r)
        end

        if maximum(res[in_contour.(Λ, c, r)]) < ϵ
            break
        end
    end

    normalize!(X)
    Λ, X, residuals(R, Λ, T)
end

function nlfeast_moments!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0, 0.0), r=1.0, debug=false, ϵ=10e-12, moments=2, store=true, spurious=1e-5)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, moments * m₀)
    θ = LinRange(π / nodes, 2 * π - π / nodes, nodes)
    Q = zeros(ComplexF64, 2 * moments, N, m₀)
    Y = zeros(ComplexF64, N, m₀ * moments)
    Q₀, Q₁ = zeros(ComplexF64, moments * N, moments * m₀), zeros(ComplexF64, moments * N, moments * m₀)
    Tinv, Temp, R = similar(X, ComplexF64), similar(X, ComplexF64), similar(Y, ComplexF64)
    resolvent = zeros(ComplexF64, m₀)
    A, B = zeros(ComplexF64, moments * m₀, moments * m₀), zeros(ComplexF64, moments * m₀, moments * m₀)
    res = Array{Float64}(undef, size(R, 2))
    Y = zeros(ComplexF64, N, m₀ * moments)
	if store
	    M = Array{ComplexF64}(undef, (nodes, N, N))
	    facts = Array{LU}(undef, nodes)
	end

    for i = 1:nodes
        z = (r * exp(θ[i] * im) + c)

		if store
			# M[i,:,:] .= ComplexF32.(T(z))
			M[i,:,:] .= T(z)
	        facts[i] = lu!(M[i,:,:])
	        Tinv .= X
	        ldiv!(facts[i], Tinv)
	        Temp .= (r * exp(θ[i] * im) / nodes) .* Tinv
		else
			Temp .= (T(z) \ X) .* (r * exp(θ[i] * im) / nodes)
			# Tinv .= (ComplexF32.(T(z))\ComplexF32.(X)) .* (r*exp(θ[i]*im)/nodes)
		end

        Q[1,:,:] .+= Temp
        for j = 2:2 * moments
            Q[j,:,:] .+= Temp .* z^(j - 1)
        end
    end
    for i in 1:moments, j in 1:moments
        Q₀[(i - 1) * N + 1:i * N, (j - 1) * m₀ + 1:j * m₀] .= Q[i + j - 1,:,:]
        Q₁[(i - 1) * N + 1:i * N, (j - 1) * m₀ + 1:j * m₀] .= Q[i + j,:,:]
    end

    S = svd!(Q₀)
    mul!(A, S.U', Q₁)
    mul!(B, A, S.V)
    mul!(A, B, Diagonal(1 ./ S.S))
    F = eigen!(A)
    mul!(Y, S.U[1:N,:], F.vectors)
    Λ .= F.values

    # beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
    # beyn_qr_step!(Q₀, Q₁, X, Λ)

	update_R_moments!(Y, R, Λ, res, T, c, r)
	X .= Y[:,1:m₀]
    # update_R!(Y, R, Λ, T)
    # res .= residuals(R, Λ, T)
    # p = sortperm(res)
	# p = sortperm(abs.(Λ .- c))
    # res .= res[p]
    # X .= Y[:, p][:,1:m₀]
    # Λ .= Λ[p]
    # R .= R[:, p]

    if debug
        iter_debug_print(0, Λ[1:m₀], res[1:m₀], c, r, spurious)
		# iter_debug_print(0, Λ, res, c, r, spurious)
    end

    for nit = 1:iter

        Q .= 0

        for i = 1:nodes
            z = (r * exp(θ[i] * im) + c)
            resolvent .= (1 ./ (z .- Λ[1:m₀])) .* (r * exp(θ[i] * im) / nodes)

			if store
	            Tinv .= R[:, 1:m₀]
	            ldiv!(facts[i], Tinv)
	            Temp .= X - Tinv
    			else
				Temp .= (X - T(z) \ R[:, 1:m₀])
				# Tinv .= (X - (ComplexF32.(T(z))\ComplexF32.(R[:, 1:m₀])))
			end
            rmul!(Temp, Diagonal(resolvent))

            Q[1,:,:] .+= Temp
            for j = 2:2 * moments
                Q[j,:,:] .+= Temp .* z^(j - 1)
            end
        end

        for i in 1:moments, j in 1:moments
            Q₀[(i - 1) * N + 1:i * N, (j - 1) * m₀ + 1:j * m₀] .= Q[i + j - 1,:,:]
            Q₁[(i - 1) * N + 1:i * N, (j - 1) * m₀ + 1:j * m₀] .= Q[i + j,:,:]
        end

        S = svd!(Q₀)
        mul!(A, S.U', Q₁)
        mul!(B, A, S.V)
        mul!(A, B, Diagonal(1 ./ S.S))
        F = eigen!(A)
        mul!(Y, S.U[1:N,:], F.vectors)
        Λ .= F.values

        # beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_rr_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_qr_step!(Q₀, Q₁, X, Λ)

        update_R_moments!(Y, R, Λ, res, T, c, r)
		X .= Y[:,1:m₀]
        # res .= residuals(R, Λ, T)
        # # p = sortperm(res)
        # # res .= res[p]
        # # X .= Y[:, p][:,1:m₀]
        # # Λ .= Λ[p]
        # # R .= R[:, p]

        if debug
            iter_debug_print(nit, Λ[1:m₀], res[1:m₀], c, r, spurious)
			# iter_debug_print(nit, Λ, res, c, r, spurious)
        end
		res_inside = res[1:m₀][in_contour.(Λ[1:m₀], c, r)]
        if size(res_inside, 1) > 0 && maximum(res_inside) < ϵ
            break
        end
		if nit > 1 && sum(res_inside .< spurious) > 0 && maximum(res_inside[res_inside .< spurious]) < ϵ
			break
		end
    end

	p = sortperm(res)
	res .= res[p]
	X .= Y[:, p][:,1:m₀]
	Λ .= Λ[p]
	R .= R[:, p]

    normalize!(X)
    # Λ[1:m₀], X, res[1:m₀]
	Λ, Y, res
end


text:
module rendahl

using Plots: _show_pdfbackends
include("solver.jl")

export solve, simdata, compute_irfs

end


text:
include("check_inputs.jl")

println("== parsing \"configuration.toml\"...")
configuration = parse_configuration(joinpath(project_root,"configuration.toml"))
check_configuration_file(configuration)

op = configuration["options"]

# prefix input
temp_errors = Dict{String,String}()
temp_prefix = input(op,"prefix", temp_errors, def = "")
check_prefix(temp_prefix)
# cd(temp_prefix)
# const PREFIX = pwd() # making prefix an absolute address
# cd(user_original_directory)
const PREFIX = realpath(temp_prefix) # making prefix an absolute address

# other input
input_errors = Dict{String,String}()
const MY_REGISTRY_UUID          = input(op,"new_registry_uuid",   input_errors, def = "")
const MY_REGISTRY_UUID_OBJECT   = check_uuid!(MY_REGISTRY_UUID,    input_errors)
const DOWNLOAD_ALL_PACKAGES     = input(op,"download_all_packages",input_errors, def = true)
const MY_REGISTRY_NAME          = input(op,"new_registry_name",   input_errors, def = "")
const MY_GIT_USER_NAME          = input(op,"git_user_name",       input_errors, def = "")
const MY_GIT_USER_EMAIL         = input(op,"git_user_email",      input_errors, def = "")
const MY_INFO                   = input(op,"your_info",           input_errors, def = "",    req = false)
const GIT_COMMIT_INFO           = input(op,"git_commit_info",     input_errors, def = "",    req = false)
const LOGGING_TO_FILE           = input(op,"log_to_file",         input_errors, def = true,  req = false)
const SILENT                    = input(op,"silent_installation", input_errors, def = false)
const UPDATE_MODE               = input(op,"update_mode",         input_errors, def = true,  req = false)
const REFRESH_CONF_FILES        = input(op,"refresh_conf_files",  input_errors, def = true,  req = false)
const MAKE_BACKUPS              = input(op,"make_backups",        input_errors, def = true,  req = false)
const CLEAN_INSTALLATION        = input(op,"clean_installation",  input_errors, def = false, req = false)

# network inputs
const INSTALL_ON_NETWORK        = input(op,"install_on_network",  input_errors, def = true,  req = false)
if INSTALL_ON_NETWORK
    const NETWORK_PREFIX = input(op,"network_access_address",input_errors, def = "", req = true)
end
const MAKE_LOCAL_ACCESS_REG = input(op, "make_local_access_reg", input_errors,  def = true, req = false)
const MINIMAL_REG           = input(op, "minimal_reg",           input_errors,  def = false, req = false)

# Internal constants
const HEADER_MASSAGE = """
Mirror registry"""

const SCREIPTS_DIR_NAME = "modules"
const BUILD_INFORMATION_DIR = "build_info"
const BUILD_LOG_FILE_NAME = "BUILD_LOG.log"
const BUILD_TIME_FILE_NAME = "BUILD_TIME.log"
const PACKAGES_DIR = "packages"
const NEW_REG_DIR = "registry"
const SRC_DIR = "op_include"
const REG_DIR = "sp_registries"
const PROJECTS_DIR = "sp_projects"
const MANIFESTS_DIR = "sp_manifests"
const KEEP_ONLINE_DIR = "keep_online"
const REMOVE_FROM_REPO_DIR = "remove_from_repo"
const DATA_DIR = "META_DATA"
const REPO_PATH_EXTENTION = ".git"
const LOCAL_REG_DIR = "local_access"
const PUBLIC_REG_DIR = "public_access"
const ERRORS_DIR = "configuration_errors"
# const REPO = "repo"
# const REG = "registries"
const PACKAGES_PATH          = joinpath(PREFIX,PACKAGES_DIR)
const NEW_REG_PATH           = joinpath(PREFIX,NEW_REG_DIR)
const BUILD_INFORMATION_PATH = joinpath(PREFIX,BUILD_INFORMATION_DIR)
const META_DATA_PATH         = joinpath(PREFIX,BUILD_INFORMATION_DIR,DATA_DIR)
const ERRORS_PATH            = joinpath(PREFIX,BUILD_INFORMATION_DIR,ERRORS_DIR)
const SRC_PATH               = joinpath(PREFIX,SRC_DIR)
const SP_REG_PATH            = joinpath(PREFIX,SRC_DIR,REG_DIR)
const SP_PROJECTS_PATH       = joinpath(PREFIX,SRC_DIR,PROJECTS_DIR)
const SP_MANIFESTS_PATH      = joinpath(PREFIX,SRC_DIR,MANIFESTS_DIR)
const KEEP_ONLINE_PATH       = joinpath(PREFIX,SRC_DIR,KEEP_ONLINE_DIR)
const REMOVE_FROM_REPO_PATH  = joinpath(PREFIX,SRC_DIR,REMOVE_FROM_REPO_DIR)

# reporting the errors.
report_errors(input_errors,path = ERRORS_PATH)


text:
# By default, generic Julia functions are not considered as constraint:
is_constraint(::Function) = false

# Generic inequality constraints
# We add constraint only on vmag_pq
function state_constraint(polar::PolarForm, g, buffer)
    index_pq = polar.indexing.index_pq
    g .= @view buffer.vmag[index_pq]
    return
end
is_constraint(::typeof(state_constraint)) = true
size_constraint(polar::PolarForm{T, IT, VT, AT}, ::typeof(state_constraint)) where {T, IT, VT, AT} = PS.get(polar.network, PS.NumberOfPQBuses())
function bounds(polar::PolarForm, ::typeof(state_constraint))
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    fr_ = npq + npv + 1
    to_ = 2*npq + npv
    return polar.x_min[fr_:to_], polar.x_max[fr_:to_]
end
# State Jacobian: Jx_i = [0, ..., 1, ... 0] where
function jacobian(polar::PolarForm, ::typeof(state_constraint), i_cons, ∂jac, buffer)
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    fr_ = npq + npv

    # Adjoint / State
    fill!(∂jac.∇fₓ, 0)
    ∂jac.∇fₓ[fr_ + i_cons] = 1.0
    # Adjoint / Control
    fill!(∂jac.∇fᵤ, 0)
end
function jtprod(polar::PolarForm, ::typeof(state_constraint), ∂jac, buffer, v)
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    fr_ = npq + npv + 1
    # Adjoint / Control
    fill!(∂jac.∇fᵤ, 0)
    # Adjoint / State
    fill!(∂jac.∇fₓ, 0)
    ∂jac.∇fₓ[fr_:end] .= v
end

# Here, the power constraints are ordered as:
# g = [P_ref; Q_ref; Q_pv]
function power_constraints(polar::PolarForm, g, buffer)
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())
    ref = polar.indexing.index_ref
    pv = polar.indexing.index_pv
    ref_to_gen = polar.indexing.index_ref_to_gen
    pv_to_gen = polar.indexing.index_pv_to_gen
    Vm, Va, pbus, qbus = buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj

    cnt = 1
    # Constraint on P_ref (generator) (P_inj = P_g - P_load)
    # NB: Active power generation has been updated previously inside buffer
    g[cnt] = buffer.pg[ref_to_gen[1]]
    cnt += 1
    # Refresh reactive power generation in buffer
    refresh!(polar, PS.Generator(), PS.ReactivePower(), buffer)
    # Constraint on Q_ref (generator) (Q_inj = Q_g - Q_load)
    # Careful: g could be a view
    if isa(g, SubArray)
        gg = g.parent
        shift = nref + g.indices[1].start - 1
    else
        gg = g
        shift = nref
    end
    if isa(gg, Array)
        kernel! = load_power_constraint_kernel!(CPU(), 4)
    else
        kernel! = load_power_constraint_kernel!(CUDADevice(), 256)
    end
    ev = kernel!(
        gg, buffer.qg, ref_to_gen, pv_to_gen, nref, npv, shift,
        ndrange=nref+npv,
    )
    wait(ev)
    return
end
is_constraint(::typeof(power_constraints)) = true
function size_constraint(polar::PolarForm{T, IT, VT, AT}, ::typeof(power_constraints)) where {T, IT, VT, AT}
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())
    return 2*nref + npv
end
function bounds(polar::PolarForm{T, IT, VT, AT}, ::typeof(power_constraints)) where {T, IT, VT, AT}
    ngen = PS.get(polar.network, PS.NumberOfGenerators())
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())

    # Get all bounds (lengths of p_min, p_max, q_min, q_max equal to ngen)
    p_min, p_max = PS.bounds(polar.network, PS.Generator(), PS.ActivePower())
    q_min, q_max = PS.bounds(polar.network, PS.Generator(), PS.ReactivePower())

    index_ref = polar.indexing.index_ref
    index_pv = polar.indexing.index_pv
    index_gen = polar.indexing.index_generators
    pv_to_gen = polar.indexing.index_pv_to_gen
    ref_to_gen = polar.indexing.index_ref_to_gen

    # Remind that the ordering is
    # g = [P_ref; Q_ref; Q_pv]
    MT = polar.AT
    pq_min = [p_min[ref_to_gen]; q_min[ref_to_gen]; q_min[pv_to_gen]]
    pq_max = [p_max[ref_to_gen]; q_max[ref_to_gen]; q_max[pv_to_gen]]
    return convert(MT, pq_min), convert(MT, pq_max)
end
# Jacobian
function jacobian(
    polar::PolarForm,
    ::typeof(power_constraints),
    i_cons,
    ∂jac,
    buffer
)
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())
    index_pv = polar.indexing.index_pv
    index_pq = polar.indexing.index_pq
    index_ref = polar.indexing.index_ref
    pv_to_gen = polar.indexing.index_pv_to_gen

    vmag = buffer.vmag
    vang = buffer.vang
    adj_x = ∂jac.∇fₓ
    adj_u = ∂jac.∇fᵤ
    adj_vmag = ∂jac.∂vm
    adj_vang = ∂jac.∂va
    adj_pg = ∂jac.∂pg
    fill!(adj_pg, 0.0)
    fill!(adj_vmag, 0.0)
    fill!(adj_vang, 0.0)
    fill!(adj_x, 0.0)
    fill!(adj_u, 0.0)

    adj_inj = 1.0
    if i_cons <= nref
        # Constraint on P_ref (generator) (P_inj = P_g - P_load)
        bus = index_ref[i_cons]
        put_active_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, polar.ybus_re, polar.ybus_im)
    elseif i_cons <= 2*nref
        # Constraint on Q_ref (generator) (Q_inj = Q_g - Q_load)
        bus = index_ref[i_cons - nref]
        put_reactive_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, polar.ybus_re, polar.ybus_im)
    else
        # Constraint on Q_pv (generator) (Q_inj = Q_g - Q_load)
        bus = index_pv[i_cons - 2* nref]
        put_reactive_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, polar.ybus_re, polar.ybus_im)
    end
    if isa(adj_x, Array)
        kernel! = put_adjoint_kernel!(CPU(), 1)
    else
        kernel! = put_adjoint_kernel!(CUDADevice(), 256)
    end
    ev = kernel!(adj_u, adj_x, adj_vmag, adj_vang, adj_pg,
                 index_pv, index_pq, index_ref, pv_to_gen,
                 ndrange=nbus)
    wait(ev)
end
# Jacobian-transpose vector product
function jtprod(
    polar::PolarForm,
    ::typeof(power_constraints),
    ∂jac,
    buffer,
    v::AbstractVector,
)
    m = size_constraint(polar, power_constraints)
    jvx = similar(∂jac.∇fₓ) ; fill!(jvx, 0)
    jvu = similar(∂jac.∇fᵤ) ; fill!(jvu, 0)
    for i_cons in 1:m
        if !iszero(v[i_cons])
            jacobian(polar, power_constraints, i_cons, ∂jac, buffer)
            jx, ju = ∂jac.∇fₓ, ∂jac.∇fᵤ
            jvx .+= jx .* v[i_cons]
            jvu .+= ju .* v[i_cons]
        end
    end
    ∂jac.∇fₓ .= jvx
    ∂jac.∇fᵤ .= jvu
end


text:
"""
	prepdata(x[, y])
Remove NaNs. Will return reduced vector if any NaN in input data
"""
function prepdata(x::Vector{Float64})::Vector{Float64}
	return any(isnan.(x)) ? filter(!isnan,x) : copy(x)
end

function prepdata(x::Vector{Float64},y::Vector{Float64})
	m = isnan.(x.+y);
	if any(m)
		return x[.!m],y[.!m]
	else
		return copy(x),copy(y)
	end
end

"""
	histfit(x,d;bins_method="sqrt",pdf_points,new_fig)
Plot (normed) histogram superimposed by fitted distribution

**Input**
* `x`: input samples
* `d`: distribution (e.g. Distributions.Normal())
* `bins_method`: method for number of bins (see `?histbins`)
* `pdf_points`: number of points used to plot PDF
* `new_fig`: plot to new figure

**Example**
```
d = Distributions.Normal(1,3);
histfit(Distributions.rand(d,1000),d)
```
"""
function histfit(x::Vector{Float64},d;bins_method="Rice",pdf_points::Int=100,new_fig=true)
	nb = histbins(x,method=bins_method);
	new_fig ? PyPlot.figure() : nothing
	PyPlot.plt[:hist](x,nb,density=true);
	PyPlot.plot(range(nb[1],stop=nb[end],length=pdf_points),
		Distributions.pdf.(d,range(nb[1],stop=nb[end],length=pdf_points)),"r-",
		linewidth=2);
end

"""
	histbins(x;method)
Compute number of bins and create bin vector using various methods
See https://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width

**Input**
* datain: input vector
* method = "Scott" | "Rice" | "sqrt" otherwise (range/n) | set to integer with required nr.of bins

**Output**
* vector for bin intervals (start/stop) => length(output)-1 = number of bins

**Example**
```
b = histbins(rand(100),"Scott");
```
"""
function histbins(datain;method="Rice")::Vector{Float64}
	xi = prepdata(datain);
	data_range = extrema(xi);
	n = length(xi);
	if method=="Scott"
		nb = 3.5*Distributions.std(xi)/(n^(1/3));
	elseif method == "Rice"
		nb = 2.0*n^(1/3);
	elseif method == "sqrt"
		nb = sqrt(n);
	elseif eltype(method) == Int
		nb = method;
	else # default
		nb = (data_range[2]-data_range[1])*0.5;
	end
	return collect(range(data_range[1],stop=data_range[2],length=trunc(Int,round(nb))+1));
end

"""
	theobins(d,b)
Count theoretical number of occurances in given bins (`b`) for given
distribution(`d`) using input data (`x`).

**Output**
* bc: expected number of occurences in each bin. Float to avoid summing issues (`sum(bc)==n`)

**Example**
d = Distributions.Normal()
b = collect(-3.:1:3);
n = 100;
bc = HypoTest.theobins(d,b,n);
"""
function theobins(d,b::Vector{Float64},n::Int)::Vector{Float64}
	bc = zeros(Float64,length(b)-1);
	for i in 1:length(bc)
		bc[i] = n*(Distributions.cdf(d,b[i+1])-Distributions.cdf(d,b[i]));
	end
	return bc
end

"""
Count number of occurrences within given bins/intervals (`b`) using input
vector (`x`)

**Output**
* bc: bin counts

**Example**
x = collect(1.:1.:10);
bc = HypoTest.countbins(x,[0.,1.,9.,10.]);
# same as PyPlot.plt[:hist](x,[0.,1.,9.,10.])
"""
function countbins(x::Vector{Float64},b::Vector{Float64})::Vector{Int}
	bc = zeros(Int,length(b)-1);
	for i = 1:length(bc)
		bc[i] = i<length(bc) ? length(findall(y-> y .>= b[i] && y .< b[i+1],x)) :
						   length(findall(y-> y .>= b[i] && y .<= b[i+1],x));
	end
	return bc
end

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
# This code is mostly copied from JSON3: https://github.com/quinnj/JSON3.jl

struct PointerString <: AbstractString
    ptr::Ptr{UInt8}
    len::Int
end

function Base.hash(s::PointerString, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

import Base: ==
function ==(x::String, y::PointerString)
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), pointer(x), y.ptr, y.len) == 0
end
==(y::PointerString, x::String) = x == y

Base.codeunit(s::PointerString) = UInt8
Base.ncodeunits(s::PointerString) = s.len
@inline function Base.codeunit(s::PointerString, i::Integer)
    @boundscheck checkbounds(s, i)
    unsafe_load(s.ptr + i - 1)
end
Base.String(x::PointerString) = unsafe_string(x.ptr, x.len)

function reverseescapechar(b)
    b == UInt8('"')  && return UInt8('"')
    b == UInt8('\\') && return UInt8('\\')
    b == UInt8('/')  && return UInt8('/')
    b == UInt8('b')  && return UInt8('\b')
    b == UInt8('f')  && return UInt8('\f')
    b == UInt8('n')  && return UInt8('\n')
    b == UInt8('r')  && return UInt8('\r')
    b == UInt8('t')  && return UInt8('\t')
    return 0x00
end

utf16_is_surrogate(c::UInt16) = (c & 0xf800) == 0xd800
utf16_get_supplementary(lead::UInt16, trail::UInt16) = Char(UInt32(lead-0xd7f7)<<10 + trail)

charvalue(b) = (UInt8('0') <= b <= UInt8('9')) ? b - UInt8('0') :
               (UInt8('a') <= b <= UInt8('f')) ? b - (UInt8('a') - 0x0a) :
               (UInt8('A') <= b <= UInt8('F')) ? b - (UInt8('A') - 0x0a) :
               throw(ArgumentError("JSON invalid unicode hex value"))

@noinline invalid_escape(str) = throw(ArgumentError("encountered invalid escape character in json string: \"$(String(str))\""))
@noinline unescaped_control(b) = throw(ArgumentError("encountered unescaped control character in json: '$(escape_string(Base.string(Char(b))))'"))

function unescape(s)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == UInt8('\\')
                i += 1
                i > n && invalid_escape(s)
                b = codeunit(s, i)
                if b == UInt8('u')
                    c = 0x0000
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    if utf16_is_surrogate(c)
                        i += 3
                        i > n && invalid_escape(s)
                        c2 = 0x0000
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        ch = utf16_get_supplementary(c, c2)
                    else
                        ch = Char(c)
                    end
                    st = codeunits(Base.string(ch))
                    for j = 1:length(st)-1
                        @inbounds buf[len] = st[j]
                        len += 1
                    end
                    b = st[end]
                else
                    b = reverseescapechar(b)
                    b == 0x00 && invalid_escape(s)
                end
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end

function unescape_blockstr(s)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == UInt8('\\')
                i + 3 > n && invalid_escape(s)
                if codeunit(s, i+1) == UInt8('"') && codeunit(s, i+2) == UInt8('"') && codeunit(s, i+3) == UInt8('"')
                    b = UInt8('"')
                    buf[len] = b
                    len += 1
                    i += 1
                    buf[len] = b
                    len += 1
                    i += 1
                    # Third " added below
                    i += 1 # skip /
                end
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end

"""
    format_block_string(str)

Perform the formatting defined [in the specification](https://spec.graphql.org/October2021/#sec-String-Value.Semantics).
"""
function format_block_string(str)
    lines = split(str, r"\r|\n") # split by line terminators

    common_indent_set = false
    common_indent = 0
    for (i, line) in enumerate(lines)
        i == 1 && continue
        len = length(line)
        indent = count_leading_whitespace(line)
        if indent < len
            if !common_indent_set || indent < common_indent
                common_indent_set = true
                common_indent = indent
            end
        end
    end

    if common_indent_set
        for (i, line) in enumerate(lines)
            i == 1 && continue
            lines[i] = @view(line[common_indent + 1:end])
        end
    end

    i = 1
    lines_to_remove = Int[]
    while all_whitespace(lines[i]) || isempty(lines[i])
        push!(lines_to_remove, i)
        i += 1
    end
    i = length(lines)
    while all_whitespace(lines[i]) || isempty(lines[i])
        push!(lines_to_remove, i)
        i -= 1
    end
    deleteat!(lines, lines_to_remove)
    formatted = ""
    for (i,line) in enumerate(lines)
        if i != 1
            formatted *= '\n'
        end
        formatted *= line
    end
    return formatted
end

function count_leading_whitespace(str)
    count = 0
    for char in str
        if char == '\t' || char == ' '
            count += 1
        else
            break
        end
    end
    return count
end

function all_whitespace(str)
    for char in str
        if !(char == '\t' || char == ' ')
            return false
        end
    end
    return true
end