text:
#-------------------------------------------------------------------------------
#             GNSSTools.jl Signal Simulation and Data Processing Demo
#-------------------------------------------------------------------------------
"""
    demo(;sigtype="l1ca", include_carrier=true, include_adc=true,
          include_noise=true, include_databits=true, simulatedata=true,
          saveto=missing, T="short", showplot=true, f_d=800.,
          fd_rate=0., prn=26, n0=1000., t_length=1e-3, M=4000,
          fd_range=5000., dll_b=2, state_num=3, dynamickf=true,
          cov_mult=1., q_a=1, figsize=missing, CN0=45., plot3d=true,
          show_acq_plot=true, doppler_curve=missing, doppler_t=missing,
          fd_center=missing, sig_freq=missing, signal=missing, q_mult=1,
          print_steps=true, σω=1000)

Runs a demo of `GNSSTools` showing major capabilities such as course/fine acquisition
and code/carrier phase and Doppler frequency tracking on simulated and real
data.
"""
function demo(;sigtype="l1ca", include_carrier=true, include_adc=true,
               include_databits=true, simulatedata=true,
               saveto=missing, T="short", showplot=true, f_d=800.,
               fd_rate=0., prn=26, n0=1000., t_length=30, M=1,
               fd_range=5000., dll_b=5, state_num=3, dynamickf=true,
               cov_mult=1., q_a=1, figsize=missing, CN0=45., plot3d=true,
               show_acq_plot=false, doppler_curve=missing, doppler_t=missing,
               fd_center=0., sig_freq=missing, signal=missing, q_mult=1,
               print_steps=true, σω=1000, channel="I", file_name=missing,
               include_thermal_noise=true, include_phase_noise=true,
               fine_acq_method=:fft, skip_to=0.1, tracking_T=missing,
               save_interval=1, file_location="", file_prefix="SimulatedSignal",
               h_parms=h_parms_tcxo[1])
    if print_steps
        println("Running GNSSTools Signal Simulation and Data Processing Demo")
    end
    # L5Q parameters
    if sigtype == "l5"
        f_s = 25e6  # Hz
        f_if = 0.  # Hz
        Tsys = 535.  # K
        phi = π/4  # rad
        nADC = 8  # bits
        B = 2*L5_chipping_rate  # Hz
        if channel == "I"
            fine_acq_T = 10e-3
        elseif channel == "Q"
            fine_acq_T = 20e-3
        else
            error("Invalid channel specified")
        end
        if ismissing(file_name)
            file_name = "hi_e06_20190411_092347_004814_1176.45M_25.0M_USRP8_X300_LB-SJ-10100-SF_Dish-LinZ.sc4"  # L5
        end
        f_center = 1176.45e6
        if ismissing(sig_freq)
            sig_freq = L5_freq
        end
        if include_databits
            signal_type = define_l5_code_type(t_length; prns=prn,
                                              sig_freq=sig_freq)
        else
            signal_type = define_l5_code_type(;prns=prn, sig_freq=sig_freq)
        end
    elseif sigtype == "l1ca"
        f_s = 5e6  # Hz
        f_if = 1.25e6  # Hz
        # f_if = 0.  # Hz
        Tsys = 535.  # K
        phi = π/4  # rad
        nADC = 8  # bits
        B = 2*l1ca_chipping_rate  # Hz
        fine_acq_T = 10e-3
        if ismissing(sig_freq)
            sig_freq = L1_freq
        end
        if ismissing(file_name)
            file_name = "hi_e06_20190411_092347_004814_1575.42M_5.0M_USRP4_X300_LB-SJ-10100-SF_Dish-LinZ.sc4"  # L1
        end
        f_center = 1575.42e6
        if include_databits
            signal_type = define_l1ca_code_type(t_length; prns=prn,
                                                sig_freq=sig_freq)
        else
            signal_type = define_l1ca_code_type(;prns=prn, sig_freq=sig_freq)
        end
    else
        error("Invalid signal type specified.")
    end
    if simulatedata
        # Simulate data
        if ~ismissing(doppler_curve) && ~ismissing(doppler_t)
            zero_idx = findfirst(t -> t == 0, doppler_t)
            f_d = doppler_curve[zero_idx]
            fd_rate = (doppler_curve[zero_idx+1]-doppler_curve[zero_idx]) /
                      (doppler_t[zero_idx+1]-doppler_t[zero_idx])
        end
        if ismissing(signal)
            data = definesignal(signal_type, f_s, save_interval; prn=prn,
                                f_if=f_if, f_d=f_d, fd_rate=fd_rate, Tsys=Tsys,
                                CN0=CN0, phi=phi, nADC=nADC,
                                include_carrier=include_carrier,
                                include_adc=include_adc,
                                include_thermal_noise=include_thermal_noise,
                                include_phase_noise=include_phase_noise,
                                code_start_idx=n0,
                                receiver_h_parms=h_parms)
            if print_steps
                p = Progress(floor(Int, t_length/save_interval), 1, "Generating PRN $(prn) $(sigtype) signal...")
            end
            file_name = ""
            operation = "replace"
            if include_thermal_noise
                new_thermal_noise = true
            else
                new_thermal_noise = false
            end
            if include_phase_noise
                new_phase_noise = true
            else
                new_phase_noise = false
            end
            for i in 1:floor(Int, t_length/save_interval)
                generatesignal!(data; doppler_curve=doppler_curve, doppler_t=doppler_t)
                file_name = signal_to_sc_data_file(data, sig_freq, file_location,
                                                file_prefix;
                                                operation=operation)
                operation = "append"
                phase_noise_init_phase = mean(real.(real(data.phase_noise[end-100:end])))
                definesignal!(data; new_thermal_noise=new_thermal_noise,
                             new_phase_noise=new_phase_noise,
                             phase_noise_init_phase=phase_noise_init_phase,
                             start_t=i*save_interval)
                if print_steps
                    next!(p)
                end
            end
            data = missing
            GC.gc()
            data = loaddata(:sc8, file_name, f_s, sig_freq-f_if,
                            sig_freq, t_length)
        else
            data = signal
        end
    else
        # Load data
        if print_steps
            print("Loading $(sigtype) data...")
        end
        file_dir = "/media/Srv3Pool2/by-location/hi/"
        file_path = string(file_dir, file_name)
        data_type = :sc4
        start_t = 1e-3
        data = loaddata(data_type, file_path, f_s, f_center, sig_freq,
                        t_length; skip_to=skip_to)
        if print_steps
            println("Done")
        end
    end
    if T == "long"
        # Use 20ms coherent integration for L5Q signal and 10ms
        # coherent integration for L5I signal
        if sigtype == "l5"
            if channel == "I"
                acquisition_T = 10e-3
                if ismissing(tracking_T)
                    tracking_T = 10e-3
                end
            elseif channel == "Q"
                acquisition_T = 20e-3
                if ismissing(tracking_T)
                    tracking_T = 20e-3
                end
            else
                error("Invalid signal type specified.")
            end
        else
            acquisition_T = 1e-3
            if ismissing(tracking_T)
                tracking_T = 1e-3
            end
        end
    else
        acquisition_T = 1e-3
        if ismissing(tracking_T)
            tracking_T = 1e-3
        end
    end
    # Process signal
    if print_steps
        print("Processing signal...")
    end
    results = process(data, signal_type, prn, channel;
                      σω=σω, fd_center=fd_center, fd_range=fd_range, M=M,
                      acquisition_T=acquisition_T, fine_acq_T=fine_acq_T,
                      tracking_T=tracking_T, cov_mult=cov_mult,
                      q_a=q_a, q_mult=q_mult, dynamickf=dynamickf, dll_b=dll_b,
                      state_num=state_num, fd_rate=fd_rate, show_plot=false,
                      fine_acq_method=fine_acq_method, return_corrresult=true)
    acqresults, trackresults, corr_result, SNR_est = results
    # Plot results and save if `saveto` is a string
    n0_est = acqresults.n0_idx_course
    fd_est = acqresults.fd_course
    println("Done")
    if showplot
        if print_steps
            print("Generating figures...")
        end
        # Course Acquisition Results
        if (sigtype == "l1ca") && show_acq_plot
            if ~plot3d
                fig, ax = subplots(1, 1)
                img = ax.imshow(corr_result, interpolation="none", aspect="auto",
                                extent=[0, size(corr_result)[2], fd_range, -fd_range],
                                vmin=minimum(corr_result), vmax=maximum(corr_result)/4)
                scatter(n0_est, fd_est, color="r", s=400, marker="o",
                        facecolor="none")
            else
                fig = figure()
                ax = Axes3D(fig)
                fd_bins = Array(-fd_range+fd_center:1/replica.t_length:fd_range+fd_center)
                x, y = meshgrid(Array(1:replica.sample_num), fd_bins./1000.)
                surf(x, y, corr_result)
                max_val_idx = argmax(collect(Iterators.flatten(corr_result)))
                scatter3D(collect(Iterators.flatten(x))[max_val_idx],
                          collect(Iterators.flatten(y))[max_val_idx],
                          collect(Iterators.flatten(corr_result))[max_val_idx],
                          c="r")
            end
            xlabel("Samples")
            ylabel("Doppler (Hz)")
            suptitle("PRN $(prn)\nfd course: $(fd_est)Hz; n₀ course: $(n0_est)")
        end
        # Tracking results
        plotresults(trackresults; saveto=saveto, figsize=figsize)
        if print_steps
            println("Done")
        end
    end
    return (trackresults, data)
end

"""
    demo(a, plane_num, satellite_per_plane, user_lla=(40.01, -105.2437, 1655);
         sigtype="l1ca", include_carrier=true, include_adc=true, include_noise=true,
         include_databits=true, T="short", showplot=true, prn=26, n0=1000.,
         t_length=1e-3, M=4000, fd_range=5000., dll_b=1., state_num=3,
         dynamickf=true, cov_mult=1., q_a=2., figsize=missing, CN0=45.,
         plot3d=true, show_acq_plot=true, saveto=missing, incl=56.,
         sig_freq=missing, t_start=3/60/24, ΔΩ=360/plane_num, q_mult=1,
         print_steps=true, eop=get_eop(), σω=1000.)

Runs a demo of `GNSSTools` showing major capabilities such as course/fine acquisition
and code/carrier phase and Doppler frequency tracking on simulated data based
off a user defined constellation.

    - `a`: semi-major axis in meters
    - `plane_num`: number of orbital planes in constellation
    - `satellite_per_plane`: number of satellites in each orbital plane
    - `user_lla`: user [latitude, longitude, altitude] in [deg, deg, meters]
                  defaults to (40.01, -105.2437, 1655), or Boulder, CO, USA
"""
function demo(a, plane_num, satellite_per_plane, user_lla=(40.01, -105.2437, 1655);
              sigtype="l1ca", include_carrier=true, include_adc=true,
              include_noise=true,
              include_databits=true, T="short", showplot=true, prn=26, n0=1000.,
              t_length=60, M=1, fd_range=5000., dll_b=1., state_num=3,
              dynamickf=true, cov_mult=1., q_a=1., figsize=missing, CN0=45.,
              plot3d=true, show_acq_plot=false, saveto=missing, incl=56.,
              sig_freq=missing, t_start=3/60/24, ΔΩ=360/plane_num, q_mult=1,
              print_steps=true, eop=get_eop(), σω=1000., channel="I",
              include_thermal_noise=true, include_phase_noise=true,
              fine_acq_method=:fft, tracking_T=missing,
              Δf_per_plane=360/satellite_per_plane/2, save_interval=1,
              file_location="", file_prefix="SignalFromLEO",
              h_parms=h_parms_tcxo[1])
    if print_steps
        println("Running GNSSTools Constellation Demo")
    end

    # L5Q parameters
    if sigtype == "l5"
        if ismissing(sig_freq)
            sig_freq = L5_freq
        end
        f_s = 25e6  # Hz
        f_if = 0.  # Hz
        Tsys = 535.  # K
        phi = π/4  # rad
        nADC = 8  # bits
        B = 2.046e7  # Hz
        if channel == "I"
            fine_acq_T = 10e-3
        elseif channel == "Q"
            fine_acq_T = 20e-3
        else
            error("Invalid channel specified")
        end
        if include_databits
            signal_type = define_l5_code_type(t_length; prns=prn,
                                              sig_freq=sig_freq)
        else
            signal_type = define_l5_code_type(;prns=prn, sig_freq=sig_freq)
        end
    elseif sigtype == "l1ca"
        if ismissing(sig_freq)
            sig_freq = L1_freq
        end
        # f_s = 25e6  #
        f_s = 5e6  # Hz
        f_if = 1.25e6  # Hz
        # f_if = 0.  # Hz
        Tsys = 535.  # K
        phi = π/4  # rad
        nADC = 8  # bits
        B = 2.046e6  # Hz
        fine_acq_T = 10e-3
        if include_databits
            signal_type = define_l1ca_code_type(t_length; prns=prn,
                                                sig_freq=sig_freq)
        else
            signal_type = define_l1ca_code_type(;prns=prn, sig_freq=sig_freq)
        end
    else
        error("Invalid signal type specified.")
    end

    # Simulate data
    if print_steps
        print("Setting signal parameters and generating constellation...")
    end
    data = definesignal(signal_type, f_s, save_interval; prn=prn,
                        f_if=f_if, f_d=0., fd_rate=0., Tsys=Tsys,
                        CN0=CN0, phi=phi, nADC=nADC,
                        include_carrier=include_carrier,
                        include_adc=include_adc,
                        include_thermal_noise=include_thermal_noise,
                        include_phase_noise=include_phase_noise,
                        code_start_idx=n0,
                        receiver_h_parms=h_parms)
    constellation_t = Array(data.start_t:1:t_length)
    if t_length >= 1
        t_step = 0.25
    else
        t_step = 0.001
    end
    doppler_t = Array(data.start_t:t_step:t_length+t_step)
    constellation = define_constellation(a, plane_num, satellite_per_plane,
                                         incl, constellation_t;
                                         show_plot=false, obs_lla=user_lla,
                                         eop=eop, t_start=t_start, ΔΩ=ΔΩ,
                                         Δf_per_plane=Δf_per_plane)
    elevations = []
    max_elevations = []
    for satellite in constellation.satellites
        elevation = calcelevation(satellite, user_lla;
                                  name=string(satellite.id))
        push!(max_elevations, maximum(elevation.sat_elevation))
        push!(elevations, elevation)
    end
    max_idx = argmax(max_elevations)
    user_ecef = GeodetictoECEF(deg2rad(user_lla[1]), deg2rad(user_lla[2]),
                               user_lla[3])
    doppler_curve = zeros(length(doppler_t))
    for i in 1:length(doppler_curve)
        julian_date = doppler_t[i]/60/60/24 .+ t_start
        doppler_curve[i] = calcdoppler(constellation.satellites[max_idx].init_orbit,
                                       julian_date, eop, user_ecef, sig_freq)
    end
    sat_rae = calcelevation(constellation.satellites[max_idx], user_lla)
    sat_elevations = sat_rae.sat_elevation
    sat_azimuth = sat_rae.sat_azimuth
    fd_rate = (doppler_curve[2]-doppler_curve[1])/0.001
    definesignal!(data; f_d=doppler_curve[1], fd_rate=fd_rate)
    if print_steps
        println("Done")
    end
    if print_steps
        p = Progress(floor(Int, t_length/save_interval), 1, "Generating PRN $(prn) $(sigtype) signal...")
    end
    file_name = ""
    operation = "replace"
    if include_thermal_noise
        new_thermal_noise = true
    else
        new_thermal_noise = false
    end
    if include_phase_noise
        new_phase_noise = true
    else
        new_phase_noise = false
    end
    for i in 1:floor(Int, t_length/save_interval)
        generatesignal!(data; doppler_curve=doppler_curve, doppler_t=doppler_t)
        file_name = signal_to_sc_data_file(data, sig_freq, file_location,
                                           file_prefix;
                                           operation=operation)
        operation = "append"
        phase_noise_init_phase = mean(real.(real(data.phase_noise[end-100:end])))
        definesignal!(data; new_thermal_noise=new_thermal_noise,
                      new_phase_noise=new_phase_noise,
                      phase_noise_init_phase=phase_noise_init_phase,
                      start_t=i*save_interval)
        if print_steps
            next!(p)
        end
    end
    data = missing
    GC.gc()
    data = loaddata(:sc8, file_name, f_s, sig_freq-f_if, sig_freq, t_length)
    if T == "long"
        # Use 20ms coherent integration for L5Q signal and 10ms
        # coherent integration for L5I signal
        if sigtype == "l5"
            if channel == "I"
                acquisition_T = 10e-3
                if ismissing(tracking_T)
                    tracking_T = 10e-3
                end
            elseif channel == "Q"
                acquisition_T = 20e-3
                if ismissing(tracking_T)
                    tracking_T = 20e-3
                end
            else
                error("Invalid signal type specified.")
            end
        else
            acquisition_T = 1e-3
            if ismissing(tracking_T)
                tracking_T = 1e-3
            end
        end
    else
        acquisition_T = 1e-3
        if ismissing(tracking_T)
            tracking_T = 1e-3
        end
    end
    # Process signal
    if print_steps
        print("Processing signal...")
    end
    Δfd = 1/acquisition_T  # Hz
    fd_center = round(doppler_curve[1]/Δfd)*Δfd  # Hz
    results = process(data, signal_type, prn, channel;
                      σω=σω, fd_center=fd_center,
                      fd_range=fd_range, M=M,
                      acquisition_T=acquisition_T,
                      fine_acq_T=fine_acq_T,
                      tracking_T=tracking_T,
                      cov_mult=cov_mult,
                      q_a=q_a, q_mult=q_mult, dynamickf=dynamickf, dll_b=dll_b,
                      state_num=state_num, fd_rate=fd_rate, show_plot=false,
                      fine_acq_method=fine_acq_method, return_corrresult=true)
    acqresults, trackresults, corr_result, SNR_est = results
    # Plot results and save if `saveto` is a string
    n0_est = acqresults.n0_idx_course
    fd_est = acqresults.fd_course
    println("Done")
    # Plot results and save if `saveto` is a string
    if showplot
        if print_steps
            print("Generating figures...")
        end
        # Course Acquisition Results
        if (sigtype == "l1ca") && show_acq_plot
            if ~plot3d
                fig, ax = subplots(1, 1)
                img = ax.imshow(corr_result, interpolation="none", aspect="auto",
                                extent=[0, size(corr_result)[2], fd_range, -fd_range],
                                vmin=minimum(corr_result), vmax=maximum(corr_result)/4)
                scatter(n0_est, fd_est, color="r", s=400, marker="o",
                        facecolor="none")
            else
                fig = figure()
                ax = Axes3D(fig)
                fd_bins = Array(-fd_range+fd_center:1/replica.t_length:fd_range+fd_center)
                x, y = meshgrid(Array(1:replica.sample_num), fd_bins./1000.)
                surf(x, y, corr_result)
                max_val_idx = argmax(collect(Iterators.flatten(corr_result)))
                scatter3D(collect(Iterators.flatten(x))[max_val_idx],
                          collect(Iterators.flatten(y))[max_val_idx],
                          collect(Iterators.flatten(corr_result))[max_val_idx],
                          c="r")
            end
            xlabel("Samples")
            ylabel("Doppler (kHz)")
            suptitle("PRN $(prn)\nfd course: $(fd_est)Hz; n₀ course: $(n0_est)")
        end
        # Tracking results
        A = sqrt(2*GNSSTools.k*Tsys)*10^(CN0/20)
        σ = sqrt(GNSSTools.k*B*Tsys)
        truth_SNR = 10log10(f_s * tracking_T * A^2 / σ^2)
        plotresults(trackresults; saveto=saveto, figsize=figsize,
                    doppler_curve=doppler_curve, doppler_t=doppler_t,
                    truth_SNR=truth_SNR)
        if print_steps
            println("Done")
        end
    end
    return (trackresults, data, doppler_curve, doppler_t, sat_elevations,
            sat_azimuth)
end


text:
## Following Qt

abstract Object

abstract Widget <: Object

abstract Observable <: Object
abstract AbstractModel <: Observable

abstract Container <: Widget
abstract Layout <: Object

abstract Control <: Widget
abstract WidgetModel <: Control
abstract WidgetVectorModel <: WidgetModel
abstract Style <: Widget

## imlement these?
abstract AbstractButton <: Widget
abstract AbstractSlider <: Widget
abstract AbstractSpinBox <: Widget

ext: jl


text:
## Contexts ##

# the low-level context constructor
function _ctx_new()
    p = ccall((:zmq_ctx_new, libzmq), Ptr{Cvoid},  ())
    if p == C_NULL
        throw(StateError(jl_zmq_error_str()))
    end
    return p
end

mutable struct Context
    data::Ptr{Cvoid}

    # need to keep a list of weakrefs to sockets for this Context in order to
    # close them before finalizing (otherwise zmq_term will hang)
    sockets::Vector{WeakRef}

    function Context()
        zctx = new(_ctx_new(), WeakRef[])
        finalizer(close, zctx)
        return zctx
    end
    function Context(::UndefInitializer)
        zctx = new(C_NULL, WeakRef[])
        finalizer(close, zctx)
        return zctx
    end
end

function Context(f::Function, args...)
    ctx = Context(args...)
    try
        f(ctx)
    finally
        close(ctx)
    end
end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, c::Context) = getfield(c, :data)

# define a global context that is initialized lazily
# and is used by default in Socket constructors, to
# save 99% of users from the low-level need to set up
# a context
const _context = Context(undef)

"""
    context()

Return the default ZMQ context (of type `Context`), initializing
it if this has not been done already.  (This context is automatically
closed when Julia exits.)
"""
function context()
    if !isopen(_context)
        setfield!(_context, :data, _ctx_new())
    end
    return _context
end

@deprecate Context(n::Integer) Context()

Base.isopen(ctx::Context) = getfield(ctx, :data) != C_NULL
function Base.close(ctx::Context)
    if isopen(ctx) # don't close twice!
        for w in getfield(ctx, :sockets)
            s = w.value
            if s isa Socket && isopen(s)
                s.linger = 0 # allow socket to shut down immediately
                close(s)
            end
        end
        empty!(getfield(ctx, :sockets))
        rc = ccall((:zmq_ctx_destroy, libzmq), Cint,  (Ptr{Cvoid},), ctx)
        setfield!(ctx, :data, C_NULL)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
    end
end
@deprecate term(ctx::Context) close(ctx)

function _get(ctx::Context, option::Integer)
    val = ccall((:zmq_ctx_get, libzmq), Cint, (Ptr{Cvoid}, Cint), ctx, option)
    if val < 0
        throw(StateError(jl_zmq_error_str()))
    end
    return val
end
function _set(ctx::Context, option::Integer, value::Integer)
    rc = ccall((:zmq_ctx_set, libzmq), Cint, (Ptr{Cvoid}, Cint, Cint), ctx, option, value)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

const ctxopts = (:io_threads, :max_sockets, :ipv6)
Base.propertynames(::Context) = ctxopts
@eval function Base.getproperty(value::Context, name::Symbol)
    $(propexpression(ctxopts) do p
        :(_get(value, $(Symbol(uppercase(String(p))))))
    end)
end
@eval function Base.setproperty!(value::Context, name::Symbol, x::Integer)
    $(propexpression(ctxopts) do p
        :(_set(value, $(Symbol(uppercase(String(p)))), x))
    end)
    return x
end

function Base.get(ctx::Context, option::Integer)
    Base.depwarn("get(ctx, option) is deprecated; use ctx.option instead", :get)
    return _get(ctx, option)
end
function set(ctx::Context, option::Integer, value::Integer)
    Base.depwarn("set(ctx, option, val) is deprecated; use ctx.option = val instead", :set)
    return _set(ctx, option, value)
end


text:
using Documenter
using FeatureDescriptors

makedocs(;
    modules=[FeatureDescriptors],
    authors="Invenia Technical Computing Corporation",
    repo="https://github.com/invenia/FeatureDescriptors.jl/blob/{commit}{path}#L{line}",
    sitename="FeatureDescriptors.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://invenia.github.io/FeatureDescriptors.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    checkdocs=:exports,
    strict=true,
)

deploydocs(;
    repo="github.com/invenia/FeatureDescriptors.jl",
    devbranch = "main",
    push_preview = true,
)

max_stars_repo_name: invenia/FeatureDescriptors

max_issues_repo_name: invenia/FeatureDescriptors

max_forks_repo_name: invenia/FeatureDescriptors


text:
"""
    Default behavior of a component. If there is no base_power field, assume is in the system's base power.
"""
function get_base_power(c::Component)
    return get_internal(c).units_info.base_value
end

function _get_multiplier(c::T) where {T <: Component}
    setting = get_internal(c).units_info
    if isnothing(setting)
        return 1.0
    elseif setting.unit_system == IS.UnitSystem.DEVICE_BASE
        return 1.0
    elseif setting.unit_system == IS.UnitSystem.SYSTEM_BASE
        numerator = get_base_power(c)
        denominator = setting.base_value
    elseif setting.unit_system == IS.UnitSystem.NATURAL_UNITS
        numerator = get_base_power(c)
        denominator = 1.0
    else
        error("Undefined Conditional")
    end
    return numerator / denominator
end

function get_value(c::Component, value::Float64)
    return _get_multiplier(c) * value
end

function get_value(c::Component, value::Min_Max)
    m = _get_multiplier(c)
    return (min = value.min * m, max = value.max * m)
end

function get_value(c::Component, value::StartUp_ShutDown)
    m = _get_multiplier(c)
    return (startup = value.startup * m, shutdown = value.shutdown * m)
end

function get_value(c::Component, value::Up_Down)
    m = _get_multiplier(c)
    return (up = value.up * m, down = value.down * m)
end

function get_value(c::Component, value::FromTo_ToFrom_Float)
    m = _get_multiplier(c)
    return (from_to = value.from_to * m, to_from = value.to_from * m)
end

function get_value(c::Component, value::Nothing)
    return value
end

function get_value(c::T, value::V) where {T <: Component, V}
    @warn("conversion not implemented for $(V) in component $(T)")
    return value::V
end

function get_value(::Nothing, _)
    return
end

function set_value(c::Component, value::Float64)
    return (1 / _get_multiplier(c)) * value
end

function set_value(c::Component, value::Min_Max)
    m = 1 / _get_multiplier(c)
    return (min = value.min * m, max = value.max * m)
end

function set_value(c::Component, value::StartUp_ShutDown)
    m = 1 / _get_multiplier(c)
    return (startup = value.startup * m, shutdown = value.shutdown * m)
end

function set_value(c::Component, value::Up_Down)
    m = 1 / _get_multiplier(c)
    return (up = value.up * m, down = value.down * m)
end

function set_value(c::Component, value::FromTo_ToFrom_Float)
    m = 1 / _get_multiplier(c)
    return (from_to = value.from_to * m, to_from = value.to_from * m)
end

function set_value(c::Component, value::Nothing)
    return value
end

function set_value(c::T, value::V) where {T <: Component, V}
    @warn("conversion not implemented for $(V) in component $(T)")
    return value::V
end

function set_value(::Nothing, _)
    return
end


text:
function make_fqparser( filename; forcegzip=false )
   fopen = open( filename, "r" )
   if isgzipped( filename ) || forcegzip
      to_open = ZlibInflateInputStream( fopen, reset_on_end=true )
   else
      to_open = BufferedInputStream( fopen )
   end
   FASTQ.Reader( to_open, fill_ambiguous=DNA_A )
end

# modified for Bio v0.2 with tryread_bool!
@inline function read_chunk!( chunk, parser )
   i = 1
   while i <= length(chunk) && !eof(parser)
      read!( parser, chunk[i] )
      i += 1
   end
   while i <= length(chunk)
      pop!(chunk) # clean up if we are at the end
   end
   parser
end

function allocate_chunk( parser; size=10000 )
  chunk = Vector{eltype(parser)}( size )
  for i in 1:length(chunk)
     chunk[i] = eltype(parser)()
  end
  chunk
end

function allocate_fastq_records( size::Int=10000 )
   chunk = Vector{FASTQRecord}( undef, size )
   for i in 1:length(chunk)
      chunk[i] = FASTQRecord()
   end
   chunk
end

function process_reads!( parser, param::AlignParam, lib::GraphLib, quant::GraphLibQuant,
                         multi::MultiMapping{SGAlignSingle}, mod::B;
                         bufsize=150, sam=false, qualoffset=33 ) where B <: BiasModel

   reads  = allocate_fastq_records( bufsize )
   mean_readlen = 0.0
   total        = 0
   mapped       = 0
   if sam
      stdbuf = BufferedOutputStream( stdout )
      write_sam_header( stdbuf, lib )
   end
   while length(reads) > 0
      read_chunk!( reads, parser )
      total += length(reads)
      @inbounds for i in 1:length(reads)
         fill!( reads[i], qualoffset )
         align = ungapped_align( param, lib, reads[i] )
         if !isnull( align )
            biasval = count!( mod, reads[i].sequence )
            if length( align.value ) > 1
               push!( multi, align.value, biasval, quant, lib )
               sam && write_sam( stdbuf, reads[i], align.value, lib, qualoffset=qualoffset )
            else
               count!( quant, align.value[1], biasval )
               sam && write_sam( stdbuf, reads[i], align.value[1], lib, qualoffset=qualoffset )
            end
            mapped += 1
            @fastmath mean_readlen += (length(reads[i].sequence) - mean_readlen) / mapped
         end
      end
      if total % 100000 == 0
         GC.gc()
      end
   end # end while
   if sam
      close(stdbuf)
   end
   mapped,total,mean_readlen
end

function process_paired_reads!( fwd_parser, rev_parser, param::AlignParam,
                                lib::GraphLib, quant::GraphLibQuant,
                                multi::MultiMapping{SGAlignPaired}, mod::B;
                                bufsize=50, sam=false, qualoffset=33 ) where B <: BiasModel

   fwd_reads  = allocate_fastq_records( bufsize )
   rev_reads  = allocate_fastq_records( bufsize )
   mean_readlen = 0.0
   total        = 0
   mapped       = 0
   if sam
      stdbuf = BufferedOutputStream( stdout )
      write_sam_header( stdbuf, lib )
   end
   while length(fwd_reads) > 0 && length(rev_reads) > 0
      read_chunk!( fwd_reads, fwd_parser )
      read_chunk!( rev_reads, rev_parser )
      total += length(fwd_reads)
      @inbounds for i in 1:length(fwd_reads)
         fill!( fwd_reads[i], qualoffset )
         fill!( rev_reads[i], qualoffset )
         fwd_aln,rev_aln = ungapped_align( param, lib, fwd_reads[i], rev_reads[i] )
         if !isnull( fwd_aln ) && !isnull( rev_aln )
            biasval = count!( mod, fwd_reads[i].sequence, rev_reads[i].sequence )
            if length( fwd_aln.value ) > 1
               push!( multi, fwd_aln.value, rev_aln.value, biasval, quant, lib )
               sam && write_sam( stdbuf, fwd_reads[i], rev_reads[i], fwd_aln.value, rev_aln.value, lib,
                                 paired=true, is_pair_rc=param.is_pair_rc, qualoffset=qualoffset )
            else
               count!( quant, fwd_aln.value[1], rev_aln.value[1], biasval )
               sam && write_sam( stdbuf, fwd_reads[i], fwd_aln.value[1], lib,
                                 paired=true, fwd_mate=true, is_pair_rc=param.is_pair_rc,
                                 qualoffset=qualoffset )
               sam && write_sam( stdbuf, rev_reads[i], rev_aln.value[1], lib,
                                 paired=true, fwd_mate=false, is_pair_rc=param.is_pair_rc,
                                 qualoffset=qualoffset )
            end
            mapped += 1
            @fastmath mean_readlen += (length(fwd_reads[i].sequence) - mean_readlen) / mapped
         end
      end
   end # end while
   if sam
      close(stdbuf)
   end
   mapped,total,mean_readlen
end


text:
using Makie, GeometryTypes, GLAbstraction
using Makie: to_signal, add_pan, add_zoom, add_mousebuttons, add_mousedrag, selection_rect
using Base: RefValue

scene = Scene()

add_mousebuttons(scene)
add_mousedrag(scene)
scene[:keyboardbuttons] = lift_node(scene[:buttons_pressed]) do x
    map(Keyboard.Button, x)
end

cam = Scene(
    :area => lift_node(FRect, scene[:window_area]),
    :projection => eye(Mat4f0),
    :view => eye(Mat4f0)
)
Makie.update_cam!(cam, to_value(cam, :area))
add_zoom(cam, scene)
add_pan(cam, scene)

projview = lift_node(*, cam[:projection], cam[:view])

function addcam(scat, cam)
    for robj in extract_renderable(scat.visual)
        robj[:view] = to_signal(cam[:view])
        robj[:projection] = to_signal(cam[:projection])
        robj[:projectionview] = to_signal(projview)
    end
end

scat = scatter(rand(100) .* 500f0, rand(100) .* 500f0, markersize = 20)
addcam(scat, cam)
scat = scatter(rand(100) .* 500f0, rand(100) .* 500f0, markersize = 20)
addcam(scat, cam)
scat = scatter(rand(100) .* 500f0, rand(100) .* 500f0, markersize = 20)
addcam(scat, cam)
a = axis(linspace(0, 500, 5), linspace(0, 500, 5))
addcam(a, cam)
rectviz, rect = selection_rect(scene, cam)

addcam(rectviz, cam)
extract_renderable(a.visual)

using GLVisualize, TextParse
using TextParse
using TextParse: Record, Field, Numeric, tryparsenext

x = assetpath("cat.obj")

obj = readstring(x)

function nv_parse(str, pos, len, opts = nothing)
    str[pos] == 'v' && str[pos + 1] == ' ' || return Nullable{Point3f0}()
    p, pos = tryparsenext(PointSpace, str, pos + 2, len)
    isnull(p) && error("Invalid object")
    Point3f0(get(p)), pos
end
CustomParser(nv_parse, Point)
findnext(obj, '\n', 7)

obj[pos - 1]
p, pos = nv_parse(obj, pos - 1, length(obj))
nv_parse(obj, pos, length(obj))

field = Field(Numeric(Float32), ignore_init_whitespace = false, ignore_end_whitespace = false, delim = ' ')
PointSpace = Record((
        field,
        field,
        Field(Numeric(Float32), eoldelim = true)
))
tryparsenext(PointSpace, str, pos + 1, length(str))

sign, i = TextParse.tryparsenext_sign(str, pos+2, length(str))

c, ii = next(str, i)
TextParse.tryparsenext_base10(Int, str, i, length(str))
str[(pos+1):pos+27]

str = obj
pos = 23
str[pos+1]
str[pos] == 'n' && str[pos + 1] == 'v' && str[pos + 2] == ' '
nv_parse()


text:
using Test
using UnROOT, LorentzVectors
using StaticArrays
using InteractiveUtils
using MD5

@static if VERSION > v"1.5.0"
    import Pkg
    Pkg.add("Polyester")
    using ThreadsX, Polyester
end

const SAMPLES_DIR = joinpath(@__DIR__, "samples")

# @stack
struct A
    n::Int32
    m::Int64
end
struct B
    o::Float16
    p
    q::Bool
end
struct C
    r
    s
    t
    u
end
@UnROOT.stack D A B C
expected_fieldnames = (:n, :m, :o, :p, :q, :r, :s, :t, :u)
expected_fieldtypes = [Int32, Int64, Float16, Any, Bool, Any, Any, Any, Any]

@testset "utils" begin
    @test fieldnames(D) == expected_fieldnames
    @test [fieldtype(D, f) for f in fieldnames(D)] == expected_fieldtypes
end

# io.jl
UnROOT.@io struct Foo
    a::Int32
    b::Int64
    c::Float32
    d::SVector{5, UInt8}
end

@testset "@io" begin

    d = SA{UInt8}[1, 2, 3, 4, 5]

    foo = Foo(1, 2, 3, d)

    @test foo.a == 1
    @test foo.b == 2
    @test foo.c ≈ 3
    @test d == foo.d

    @test 32 == sizeof(Foo)
    @test 21 == UnROOT.packedsizeof(Foo)

    buf = IOBuffer(Vector{UInt8}(1:sizeof(Foo)))
    foo = UnROOT.unpack(buf, Foo)

    @test foo.a == 16909060
    @test foo.b == 361984551142689548
    @test foo.c ≈ 4.377526f-31
    @test foo.d == UInt8[0x11, 0x12, 0x13, 0x14, 0x15]
end

struct Bar
    x::Int8
    y::UInt16
end

@testset "io functions" begin
    @test 21 == UnROOT.packedsizeof(Foo)
    @test 3 == UnROOT.packedsizeof(Bar)
end

@testset "Header and Preamble" begin
    fobj = open(joinpath(SAMPLES_DIR, "km3net_online.root"))
    file_preamble = UnROOT.unpack(fobj, UnROOT.FilePreamble)
    @test "root" == String(file_preamble.identifier)

    file_header = UnROOT.unpack(fobj, UnROOT.FileHeader32)
    @test 100 == file_header.fBEGIN
end

@testset "ROOTFile" begin
    ROOTFile(joinpath(SAMPLES_DIR, "tree_with_histos.root")) do rootfile
        @test 100 == rootfile.header.fBEGIN
        @test 1 == length(rootfile.directory.keys)
        @test "t1" ∈ keys(rootfile)
        for key in keys(rootfile)
            rootfile[key]
        end
    end

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_custom_struct.root"))
    @test 100 == rootfile.header.fBEGIN
    @test 1 == length(rootfile.directory.keys)
    @test "T" ∈ keys(rootfile)
    for key in keys(rootfile)
        rootfile[key]
    end
    close(rootfile)

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "histograms.root"))
    for branch in ["one", "two", "three"]
        @test branch in keys(rootfile)
    end
    close(rootfile)

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    @test 100 == rootfile.header.fBEGIN
    @test 10 == length(rootfile.directory.keys)
    @test "E" ∈ keys(rootfile)
    @test "META" ∈ keys(rootfile)
    @test "JTRIGGER::JTriggerParameters" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L0" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L1" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L2" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_SN" ∈ keys(rootfile)
    @test "KM3NET_EVENT" ∈ keys(rootfile)
    @test "KM3NET_SUMMARYSLICE" ∈ keys(rootfile)
    # for key in keys(rootfile)
    #     @test !ismissing(rootfile[key])
    # end
    close(rootfile)
end

@testset "readbasketsraw()" begin
    array_md5 = [0xb4, 0xe9, 0x32, 0xe8, 0xfb, 0xff, 0xcf, 0xa0, 0xda, 0x75, 0xe0, 0x25, 0x34, 0x9b, 0xcd, 0xdf]
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    data, offsets = UnROOT.array(rootfile, "KM3NET_EVENT/KM3NET_EVENT/snapshotHits"; raw=true)
    @test array_md5 == md5(data)
    close(rootfile)

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_jagged_array.root"))
    data, offsets = UnROOT.array(rootfile, "t1/int32_array"; raw=true)

    @test data isa Vector{UInt8}
    @test offsets isa Vector{Int32}
    @test data[1:3] == UInt8[0x40, 0x00, 0x00]

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_vector_multiple_baskets.root"))
    data, offsets = UnROOT.array(rootfile, "t1/b1"; raw=true)
    @test unique(diff(offsets)) == [18]
    close(rootfile)
end

@testset "No (basket) compression" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "uncomressed_lz4_int32.root"))
    arr = UnROOT.array(rootfile, "t1/int32_array")
    @test length(arr) == 3
    @test all(arr .== [[1,2], [], [3]])
    close(rootfile)
end

@testset "Uncompressed trees" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue87_uncompressed_a.root"))
    @test LazyTree(rootfile,"Events").Jet_pt ≈ [[27.3245, 24.8896, 20.8534],Float32[],[5.3306]]

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue87_uncompressed_b.root"))
    @test LazyTree(rootfile,"Events").myval[2:5] ≈ [[0.1], [0.2, 1.2], [0.3, 1.3, 2.3], [0.4, 1.4, 2.4, 3.4]]
end

@testset "Compressions" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array_lzma.root"))
    @test rootfile isa ROOTFile
    arr = UnROOT.array(rootfile, "t1/float_array")
    @test 100000 == length(arr)
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ arr[1:5] atol=1e-7
    close(rootfile)

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array_lz4.root"))
    arr = collect(rootfile["t1/float_array"])
    @test 100000 == length(arr)
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ arr[1:5] atol=1e-7
    close(rootfile)

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_int_array_zstd.root"))
    arr = collect(rootfile["t1/a"])
    @test arr == 0:99
    close(rootfile)
end

@testset "ROOTDirectoryHeader" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_histos.root"))
    header = rootfile.directory.header
    @test 5 == header.fVersion
    @test 1697049339 == header.fDatimeC
    @test 1697049339  == header.fDatimeM
    @test 111 == header.fNbytesKeys
    @test 78 == header.fNbytesName
    @test 100 == header.fSeekDir
    @test 0 == header.fSeekParent
    @test 1398 == header.fSeekKeys
    close(rootfile)

    # rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    # header = rootfile.directory.header
    # @test 5 == header.fVersion
    # @test 1658540644 == header.fDatimeC
    # @test 1658540645 == header.fDatimeM
    # @test 629 == header.fNbytesKeys
    # @test 68 == header.fNbytesName
    # @test 100 == header.fSeekDir
    # @test 0 == header.fSeekParent
    # @test 1619244 == header.fSeekKeys
    # close(rootfile)
end

@testset "LazyBranch and LazyTree" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array.root"))
    branch = rootfile["t1"]["int32_array"]
    arr = UnROOT.array(rootfile, branch)
    arr2 = UnROOT.arrays(rootfile, "t1")[1]

    @test hash(branch) == hash(rootfile["t1"]["int32_array"])
    @test hash(branch) != hash(rootfile["t1"]["float_array"])
    @test arr == arr2

    table = LazyTree(rootfile, "t1")
    BA = LazyBranch(rootfile, branch)
    @test length(arr) == length(BA)
    @test BA[1] == arr[1]
    @test BA[end] == arr[end]
    @test BA[20:30] == arr[20:30]
    @test BA[1:end] == arr
    @test table.int32_array[20:30] == BA[20:30]
    @test table[:, :int32_array][20:30] == BA[20:30]
    @test table[23, :int32_array] == BA[23]
    @test table[20:30, :int32_array] == BA[20:30]
    @test table[:].int32_array[20:30] == BA[20:30]
    @test [row.int32_array for row in table[20:30]] == BA[20:30]
    @test sum(table.int32_array) == sum(row.int32_array for row in table)
    @test [row.int32_array for row in table] == BA
    close(rootfile)
end

@testset "TLorentzVector" begin
    # 64bits T
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "TLorentzVector.root"))
    branch = rootfile["t1/LV"]
    tree = LazyTree(rootfile, "t1")

    @test branch[1].x == 1.0
    @test branch[1].t == 4.0
    @test eltype(branch) === LorentzVectors.LorentzVector{Float64}
    @test tree[1].LV.x == 1.0
    @test tree[1].LV.t == 4.0
    close(rootfile)

    # jagged LVs
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "Jagged_TLorentzVector.root"))
    branch = rootfile["t1/LVs"]
    tree = LazyTree(rootfile, "t1")

    @test eltype(branch) <: AbstractVector{LorentzVectors.LorentzVector{Float64}}
    @test eltype(branch) <: SubArray
    @test length.(branch[1:10]) == 0:9
    close(rootfile)
end

@testset "TNtuple" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "TNtuple.root"))
    arrs = [collect(rootfile["n1/$c"]) for c in "xyz"]
    @test length.(arrs) == fill(100, 3)
    @test arrs[1] ≈ 0:99
    @test arrs[2] ≈ arrs[1] .+ arrs[1] ./ 13
    @test arrs[3] ≈ arrs[1] .+ arrs[1] ./ 17
    close(rootfile)
end

@testset "Singly jagged branches" begin
    # 32bits T
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_jagged_array.root"))
    data = rootfile["t1/int32_array"]
    @test data[1] == Int32[]
    @test data[1:2] == [Int32[], Int32[0]]
    @test data[end] == Int32[90, 91, 92, 93, 94, 95, 96, 97, 98]
    close(rootfile)

    # 64bits T
    T = Float64
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_jagged_array_double.root"))
    data = rootfile["t1/double_array"]
    @test data isa AbstractVector
    @test eltype(data) <: AbstractVector{T}
    @test data[1] == T[]
    @test data[1:2] == [T[], T[0]]
    @test data[end] == T[90, 91, 92, 93, 94, 95, 96, 97, 98]
    close(rootfile)
end

@testset "Doubly jagged branches" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_doubly_jagged.root"))
    vvi = [[[2], [3, 5]], [[7, 9, 11], [13]], [[17], [19], []], [], [[]]]
    vvf = [[[2.5], [3.5, 5.5]], [[7.5, 9.5, 11.5], [13.5]], [[17.5], [19.5], []], [], [[]]]
    @test UnROOT.array(rootfile, "t1/bi") == vvi
    @test rootfile["t1/bi"] == vvi
    @test eltype(eltype(eltype(rootfile["t1/bi"]))) === Int32
    @test UnROOT.array(rootfile, "t1/bf") == vvf
    @test rootfile["t1/bf"] == vvf
    @test eltype(eltype(eltype(rootfile["t1/bf"]))) === Float32
    close(rootfile)
end

@testset "NanoAOD" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root"))
    event = UnROOT.array(rootfile, "Events/event")
    @test event[1:3] == UInt64[12423832, 12423821, 12423834]
    Electron_dxy = rootfile["Events/Electron_dxy"]
    @test eltype(Electron_dxy) == SubArray{Float32, 1, Vector{Float32}, Tuple{UnitRange{Int64}}, true}
    @test Electron_dxy[1:3] ≈ [Float32[0.0003705], Float32[-0.00981903], Float32[]]
    HLT_Mu3_PFJet40 = UnROOT.array(rootfile, "Events/HLT_Mu3_PFJet40")
    @test eltype(HLT_Mu3_PFJet40) == Bool
    @test HLT_Mu3_PFJet40[1:3] == [false, true, false]
    tree = LazyTree(rootfile, "Events", [r"Muon_(pt|eta|phi)$", "Muon_charge", "Muon_pt"])
    @test sort(propertynames(tree) |> collect) == sort([:Muon_pt, :Muon_eta, :Muon_phi, :Muon_charge])
    tree = LazyTree(rootfile, "Events", r"Muon_(pt|eta)$")
    @test sort(propertynames(tree) |> collect) == sort([:Muon_pt, :Muon_eta])
    @test occursin("LazyEvent", repr(first(iterate(tree))))
    @test sum(rootfile["Events/HLT_Mu3_PFJet40"]) == 443
    close(rootfile)
end

@testset "Branch filtering" begin
    # Branch selection behavior: if not regex, require exact name match
    treebranches = ["Muon_pt", "Muon_eta", "Muon_phi", "Muon_charge", "Muon_ptErr",
                   "Muon_", "_pt", "Muon.pt"]
    _m(s::AbstractString) = isequal(s)
    _m(r::Regex) = Base.Fix1(occursin, r)
    filter_branches(selected) = Set(mapreduce(b->filter(_m(b), treebranches), ∪, selected))
    @test (filter_branches([r"Muon_(pt|eta|phi)$", "Muon_charge", "Muon_pt"]) ==
           Set(["Muon_pt", "Muon_eta", "Muon_phi", "Muon_charge"]))
    @test filter_branches(["Muon_pt"]) == Set(["Muon_pt"])
    @test filter_branches(["Muon.pt"]) == Set(["Muon.pt"])
end

@testset "Displaying files" begin
    files = filter(x->endswith(x, ".root"), readdir(SAMPLES_DIR))
    _io = IOBuffer()
    for f in files
        r = ROOTFile(joinpath(SAMPLES_DIR, f))
        show(_io, r)
        close(r)
    end

    # test that duplicate trees (but different cycle numbers)
    # are only displayed once, and that histograms show up
    f = UnROOT.samplefile("tree_cycles_hist.root")
    @test length(collect(eachmatch(r"Events", repr(f)))) == 1
    @test length(collect(eachmatch(r"myTH2F", repr(f)))) == 1
    close(f)
end

@testset "Displaying trees" begin
    f = UnROOT.samplefile("NanoAODv5_sample.root")
    t = LazyTree(f, "Events", ["nMuon","MET_pt","Muon_pt"])
    _io = IOBuffer()
    show(_io, t)
    show(_io, t[1:10])
    show(_io, t.Muon_pt)
    show(_io, t.Muon_pt[1:10])
    s = repr(t[1:10])
    @test length(collect(eachmatch(r"Float32\[", s))) == 0
    _io = IOBuffer()
    show(_io, t)
    _iostring = String(take!(_io))
    @test length(split(_iostring,'\n')) > length(t)
    @test occursin("───────", _iostring)
    @test !occursin("NamedTuple", _iostring)
    show(_io, t; crop=:both)
    @test length(split(String(take!(_io)),'\n')) <= Base.displaysize()[1]
    _io = IOBuffer()
    show(_io, MIME"text/html"(), t)
    _iostring = String(take!(_io))
    @test occursin("</table>", _iostring)
    close(f)
end

# Custom bootstrap things

@testset "custom boostrapping" begin
    # manual interpretation (splitting)
    f_manual = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))

    data, offsets = UnROOT.array(f_manual, "KM3NET_EVENT/KM3NET_EVENT/KM3NETDAQ::JDAQEventHeader"; raw=true)
    headers_manual = UnROOT.splitup(data, offsets, UnROOT._KM3NETDAQEventHeader; jagged=false)

    data, offsets = UnROOT.array(f_manual, "KM3NET_EVENT/KM3NET_EVENT/snapshotHits"; raw=true)
    event_hits_manual = UnROOT.splitup(data, offsets, UnROOT._KM3NETDAQHit; skipbytes=10)

    data, offsets = UnROOT.array(f_manual, "KM3NET_EVENT/KM3NET_EVENT/triggeredHits"; raw=true)
    event_thits_manual = UnROOT.splitup(data, offsets, UnROOT._KM3NETDAQTriggeredHit; skipbytes=10)

    close(f_manual)  # we can close, everything is in memory

    # automatic interpretation
    customstructs = Dict(
            "KM3NETDAQ::JDAQEvent.snapshotHits" => Vector{UnROOT._KM3NETDAQHit},
            "KM3NETDAQ::JDAQEvent.triggeredHits" => Vector{UnROOT._KM3NETDAQTriggeredHit},
            "KM3NETDAQ::JDAQEvent.KM3NETDAQ::JDAQEventHeader" => UnROOT._KM3NETDAQEventHeader
    )
    f_auto = UnROOT.ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"), customstructs=customstructs)
    headers_auto = f_auto["KM3NET_EVENT/KM3NET_EVENT/KM3NETDAQ::JDAQEventHeader"]
    event_hits_auto = f_auto["KM3NET_EVENT/KM3NET_EVENT/snapshotHits"]
    event_thits_auto = f_auto["KM3NET_EVENT/KM3NET_EVENT/triggeredHits"]

    for event_hits ∈ [event_hits_manual, event_hits_auto]
        @test length(event_hits) == 3
        @test length(event_hits[1]) == 96
        @test length(event_hits[2]) == 124
        @test length(event_hits[3]) == 78
        @test event_hits[1][1].dom_id == 806451572
        @test event_hits[1][1].tdc == 30733918
        @test event_hits[1][end].dom_id == 809544061
        @test event_hits[1][end].tdc == 30735112
        @test event_hits[3][1].dom_id == 806451572
        @test event_hits[3][1].tdc == 63512204
        @test event_hits[3][end].dom_id == 809544061
        @test event_hits[3][end].tdc == 63512892
    end
    for event_thits ∈ [event_thits_manual, event_thits_auto]
        @test length(event_thits) == 3
        @test length(event_thits[1]) == 18
        @test length(event_thits[2]) == 53
        @test length(event_thits[3]) == 9
        @test event_thits[1][1].dom_id == 806451572
        @test event_thits[1][1].tdc == 30733918
        @test event_thits[1][end].dom_id == 808972598
        @test event_thits[1][end].tdc == 30733192
        @test event_thits[3][1].dom_id == 808447186
        @test event_thits[3][1].tdc == 63511558
        @test event_thits[3][end].dom_id == 809526097
        @test event_thits[3][end].tdc == 63511708
    end

    for headers ∈ [headers_manual, headers_auto]
        @test length(headers) == 3
        for header in headers
            @test header.run == 6633
            @test header.detector_id == 44
            @test header.UTC_seconds == 0x5dc6018c
        end
        @test headers[1].frame_index == 127
        @test headers[2].frame_index == 127
        @test headers[3].frame_index == 129
        @test headers[1].UTC_16nanosecondcycles == 0x029b9270
        @test headers[2].UTC_16nanosecondcycles == 0x029b9270
        @test headers[3].UTC_16nanosecondcycles == 0x035a4e90
        @test headers[1].trigger_counter == 0
        @test headers[2].trigger_counter == 1
        @test headers[3].trigger_counter == 0
        @test headers[1].trigger_mask == 22
        @test headers[2].trigger_mask == 22
        @test headers[3].trigger_mask == 4
        @test headers[1].overlays == 6
        @test headers[2].overlays == 21
        @test headers[3].overlays == 0
    end

    close(f_auto)
end

# Histograms
@testset "histograms" begin
    f = ROOTFile(joinpath(SAMPLES_DIR, "histograms1d2d.root"))
    for k in ["myTH1F", "myTH1D"]
        @test f[k][:fName] == k
        @test f[k][:fEntries] == 4.0
        @test f[k][:fSumw2] == [0.0, 800.0, 2.0, 0.0]
        @test f[k][:fXaxis_fXmin] == -2.0
        @test f[k][:fXaxis_fXmax] == 2.0
        @test f[k][:fXaxis_fXbins] == []
        @test f[k][:fXaxis_fNbins] == 2
        @test f[k][:fN] == [0.0, 40.0, 2.0, 0.0]
    end

    k = "myTH1D_nonuniform"
    @test f[k][:fName] == k
    @test f[k][:fEntries] == 4.0
    @test f[k][:fSumw2] == [0.0, 800.0, 2.0, 0.0]
    @test f[k][:fXaxis_fXmin] == -2.0
    @test f[k][:fXaxis_fXmax] == 2.0
    @test f[k][:fXaxis_fXbins] == [-2, 1, 2]
    @test f[k][:fXaxis_fNbins] == 2
    @test f[k][:fN] == [0.0, 40.0, 2.0, 0.0]

    for k in ["myTH2F", "myTH2D"]
        @test f[k][:fName] == k
        @test f[k][:fEntries] == 4.0
        @test f[k][:fSumw2] == [0.0, 0.0, 0.0, 0.0, 0.0, 400.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 400.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        @test f[k][:fXaxis_fXmin] == -2.0
        @test f[k][:fXaxis_fXmax] == 2.0
        @test f[k][:fXaxis_fXbins] == []
        @test f[k][:fXaxis_fNbins] == 2
        @test f[k][:fYaxis_fXmin] == -2.0
        @test f[k][:fYaxis_fXmax] == 2.0
        @test f[k][:fYaxis_fXbins] == []
        @test f[k][:fYaxis_fNbins] == 4
        @test f[k][:fN] == [0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    end

    @test UnROOT.parseTH(f["myTH1F"]) == ([40.0, 2.0], (-2.0:2.0:2.0,), [800.0, 2.0])
    @test UnROOT.parseTH(f["myTH2D"]) == ([20.0 0.0 0.0 20.0; 1.0 0.0 0.0 1.0], (-2.0:2.0:2.0, -2.0:1.0:2.0), [400.0 0.0 0.0 400.0; 1.0 0.0 0.0 1.0])
    @test UnROOT.parseTH(f["myTH1D_nonuniform"]) == ([40.0, 2.0], ([-2.0, 1.0, 2.0],), [800.0, 2.0])

    close(f)

    f = ROOTFile(joinpath(SAMPLES_DIR, "cms_ntuple_wjet.root"))
    binlabels = ["Root", "Weight", "Preselection", "SelectGenPart", "GoodRunsList", "EventFilters", "SelectLeptons", "SelectJets", "Trigger", "ObjectsSelection", "SSPreselection", "NjetGeq4", "AK4CategTagHiggsJets", "AK4CategTagVBSJets", "AK4CategChannels", "AK4CategPresel"]
    @test f["AK4CategPresel_cutflow"][:fXaxis_fModLabs].objects == binlabels
    close(f)
end

# Issues

@testset "issues" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue7.root"))
    @test 2 == length(keys(rootfile))
    @test [1.0, 2.0, 3.0] == UnROOT.array(rootfile, "TreeD/nums")
    @test [1.0, 2.0, 3.0] == UnROOT.array(rootfile, "TreeF/nums")
    close(rootfile)

    # issue 55
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "cms_ntuple_wjet.root"))
    pts1 = UnROOT.array(rootfile, "variable/met_p4/fCoordinates/fCoordinates.fPt"; raw=false)
    pts2 = LazyTree(rootfile, "variable", [r"met_p4/fCoordinates/.*", "mll"])[!, Symbol("met_p4/fCoordinates/fCoordinates.fPt")]
    pts3 = rootfile["variable/good_jets_p4/good_jets_p4.fCoordinates.fPt"]
    @test 24 == length(pts1)
    @test Float32[69.96958, 25.149912, 131.66693, 150.56802] == pts1[1:4]
    @test pts1 == pts2
    @test pts3[1:2] == [[454.0, 217.5, 89.5, 30.640625], [184.375, 33.28125, 32.28125, 28.46875]]
    close(rootfile)

    # issue 61
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue61.root"))
    @test rootfile["Events/Jet_pt"][:] == Vector{Float32}[[], [27.324587, 24.889547, 20.853024], [], [20.33066], [], []]
    close(rootfile)

    # issue 78
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue61.root"))
    arr = LazyTree(rootfile,"Events").Jet_pt;
    _ = length.(arr);
    @test length.(arr.buffer) == length.(arr.buffer_range)
    close(rootfile)

    # issue 108
    # unsigned short -> Int16, ulong64 -> UInt64
    # file minified with `rooteventselector --recreate -l 2 "trackntuple.root:trackingNtuple/tree" issue108_small.root`
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue108_small.root"))
    @test rootfile["tree/trk_algoMask"][2] == [0x0000000000004000, 0x0000000000004000, 0x0000000000004000, 0x0000000000004000]
    @test rootfile["tree/pix_ladder"][3][1:5] == UInt16[0x0001, 0x0001, 0x0001, 0x0001, 0x0003]
    close(rootfile)

    # issue 116
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue116.root"))
    @test length(rootfile["fTree"].fBranches.elements) == 112
    close(rootfile)
end

@testset "jagged subbranch type by leaf" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_offline.root"))
    times = UnROOT.array(rootfile, "E/Evt/trks/trks.t")
    @test times[1][1] ≈ 7.0311446e7
    @test times[10][11] ≈ 5.4956456e7

    ids_jagged = UnROOT.array(rootfile, "E/Evt/trks/trks.id")
    @test all(ids_jagged[1] .== collect(1:56))
    @test all(ids_jagged[9] .== collect(1:54))

    close(rootfile)
end

@testset "Type stability" begin
    function isfullystable(func)
        io = IOBuffer()
        print(io, (@code_typed func()).first);
        typed = String(take!(io))
        return !occursin("::Any", typed)
    end

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root"))
    t = LazyTree(rootfile, "Events", ["MET_pt"])[1:10]

    function f1()
        s = 0.0f0
        for evt in t
            s += evt.MET_pt
        end
        s
    end
    f2() = sum(t.MET_pt)

    @test isfullystable(f1)
    @test isfullystable(f2)

    close(rootfile)
end

@testset "Parallel and enumerate interface" begin
    t = LazyTree(ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root")), "Events", ["Muon_pt"])
    @test eachindex(enumerate(t)) == eachindex(t)
    nmu = 0
    for evt in t
        nmu += length(evt.Muon_pt)
    end
    @test nmu == 878

    nmu = 0
    for (i,evt) in enumerate(t)
        nmu += length(evt.Muon_pt)
    end
    @test nmu == 878

    if get(ENV, "CI", "false") == "true"
        @test Threads.nthreads() > 1
    end
    nmus = zeros(Int, Threads.nthreads())
    Threads.@threads for i in 1:length(t)
        nmus[Threads.threadid()] += length(t.Muon_pt[i])
    end
    @test sum(nmus) == 878

    et = enumerate(t)
    @test firstindex(et) == firstindex(t)
    @test lastindex(et) == lastindex(t)
    test_i, test_evt = et[2]
    @test test_i == 2
    @test test_evt isa UnROOT.LazyEvent
    @test !isempty(hash(t.Muon_pt.b))
end

@static if VERSION > v"1.5.1"
    t = LazyTree(ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root")), "Events", ["Muon_pt"])
    @testset "Multi threading" begin
        nmus = zeros(Int, Threads.nthreads())
        Threads.@threads for (i, evt) in enumerate(t)
            nmus[Threads.threadid()] += length(t.Muon_pt[i])
        end
        @test sum(nmus) == 878

        nmus .= 0
        Threads.@threads for evt in t
            nmus[Threads.threadid()] += length(evt.Muon_pt)
        end
        @test count(>(0), nmus) > 1 # test @threads is actually threading
        @test sum(nmus) == 878

        nmus .= 0
        @batch for evt in t
            nmus[Threads.threadid()] += length(evt.Muon_pt)
        end
        @test count(>(0), nmus) > 1 # test @batch is actually threading
        @test sum(nmus) == 878

        nmus .= 0
        t_dummy = LazyTree(ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root")), "Events", ["Muon_pt"])
        @batch for evt in vcat(t,t_dummy) # avoid using the same underlying file handler
            nmus[Threads.threadid()] += length(evt.Muon_pt)
        end
        @test sum(nmus) == 2*878

        for j in 1:3
            inds = [Vector{Int}() for _ in 1:Threads.nthreads()]
            Threads.@threads for (i, evt) in enumerate(t)
                push!(inds[Threads.threadid()], i)
            end
            @test sum([length(inds[i] ∩ inds[j]) for i=1:length(inds), j=1:length(inds) if j>i]) == 0
        end
    end
end

@testset "TDirectory" begin
    f = UnROOT.samplefile("tdir_complicated.root")
    @test length(keys(f["mydir"])) == 4
    @test sort(keys(f["mydir"])) == ["Events", "c", "d", "mysubdir"]
    @test sort(keys(f["mydir/mysubdir"])) == ["e", "f"]
    @test sum(length.(LazyTree(f, "mydir/Events").Jet_pt)) == 4
    @test sum(length.(f["mydir/Events/Jet_pt"])) == 4

    f = UnROOT.samplefile("issue11_tdirectory.root")
    @test sum(f["Data/mytree/Particle0_E"]) ≈ 1012.0
end

@testset "Basic C++ types" begin
    f = UnROOT.samplefile("tree_basictypes.root")
    onesrow = LazyTree(f,"t")[2] |> collect |> values .|> first .|> Int
    @test all(onesrow .== 1)
end

@testset "basketarray_iter()" begin
    f = UnROOT.samplefile("tree_with_vector_multiple_baskets.root")
    t = LazyTree(f,"t1")
    @test (UnROOT.basketarray_iter(f, f["t1"]["b1"]) .|> length) == [1228, 1228, 44]
    @test (UnROOT.basketarray_iter(t.b1) .|> length) == [1228, 1228, 44]
    @test length(UnROOT.basketarray(t.b1, 1)) == 1228
end

@testset "Cluster ranges" begin
    t = LazyTree(UnROOT.samplefile("tree_with_clusters.root"),"t1");
    @test all(UnROOT._clusterbytes(t; compressed=true) .< 10000)
    @test all(UnROOT._clusterbytes(t; compressed=false) .< 10000)
    @test UnROOT._clusterbytes([t.b1,t.b2]) == UnROOT._clusterbytes(t)
    @test length(UnROOT._clusterranges([t.b1])) == 157
    @test length(UnROOT._clusterranges([t.b2])) == 70
    @test length(UnROOT._clusterranges(t)) == 18 # same as uproot4
    @test sum(UnROOT._clusterbytes([t.b1]; compressed=true)) == 33493.0 # same as uproot4
    @test sum(UnROOT._clusterbytes([t.b2]; compressed=true)) == 23710.0 # same as uproot4
end

@testset "Vcat/chaining" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root"))
    t = LazyTree(rootfile, "Events", ["nMuon", "Muon_pt"])
    tt = vcat(t,t)
    @test (@allocated vcat(t,t)) < 1000
    @test length(tt) == 2*length(t)
    s1 = sum(t.nMuon)
    s2 = sum(tt.nMuon)
    @test s2 == 2*s1
    alloc1 = @allocated sum(length, t.Muon_pt)
    alloc2 = @allocated sum(evt->length(evt.nMuon), tt)
    @test alloc2 < 2.1 * alloc1
    close(rootfile)
end

@static if VERSION > v"1.5.0"
    @testset "Broadcast fusion" begin
        rootfile = ROOTFile(joinpath(SAMPLES_DIR, "NanoAODv5_sample.root"))
        t = LazyTree(rootfile, "Events", "nMuon")
        testf(evt) = evt.nMuon == 4
        testf2(evt) = evt.nMuon == 4
        alloc1 = @allocated a1 = testf.(t)
        alloc1 += @allocated a2 = testf2.(t)
        alloc1 += @allocated idx1 = findall(a1 .& a2)
        alloc2 = @allocated idx2 = findall(@. testf(t) & testf2(t))
        @assert !isempty(idx1)
        @test idx1 == idx2
        # compiler optimization is good on 1.8
        @test alloc1 > 1.4*alloc2
    end
end


text:
using Spring

sp=SpringParam(ls=0.1,ms=1,ks=1,ns=400)
ip=InitParam(λ=1.0,shift=1.5,pos=sp.ls*sp.ns/2)
ap=AnimParam(δt=1.0,nδt=1000,nδtperframe=10)

animate_spring(sp,ip,ap)


text:
"""
    FRAPMesh{T} <: Mesh{T}

A QBD-RAP (FRAM) discretisation scheme for a DiscretisedFluidQueue.

# Arguments
- `nodes`: the cell edges
- `me`: the MatrixExponential used to approximate model the fluid queue on each cell.
"""
struct FRAPMesh{T} <: Mesh{T}
    nodes::T
    me::AbstractMatrixExponential
end
FRAPMesh{T}(nodes::T,n_bases::Int) where T = FRAPMesh{T}(nodes,build_me(cme_params[n_bases]))
FRAPMesh(nodes::AbstractVector{Float64},n_bases::Int) = FRAPMesh(nodes,build_me(cme_params[n_bases]))
function FRAPMesh()
    FRAPMesh(Array{Float64,1}(undef,0),0)
end

"""

    n_bases_per_cell(mesh::FRAPMesh)

Number of bases in a cell
"""
n_bases_per_cell(mesh::FRAPMesh) = _order(mesh.me)

"""

    cell_nodes(mesh::FRAPMesh)

The cell centre
"""
cell_nodes(mesh::FRAPMesh) = Array(((mesh.nodes[1:end-1] + mesh.nodes[2:end]) / 2 )')

"""

    basis(mesh::FRAPMesh)

Constant ""
"""
basis(mesh::FRAPMesh) = ""

function build_lazy_generator(
    dq::DiscretisedFluidQueue{FRAPMesh{T}};
    v::Bool=false,
) where T
    me = dq.mesh.me
    blocks = (me.s*me.a, me.S, me.s*me.a)
    boundary_flux = OneBoundaryFlux(me.s[:],me.a[:])
    D = me.D
    out = LazyGenerator(dq,blocks,boundary_flux,D)
    v && println("UPDATE: LazyGenerator object created with keys ", keys(out))
    return out
end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
using LinearAlgebra
using SparseArrays

export Laplacian2D,sparse_corr

struct Laplacian2D
    lpx::Laplacian1D
    lpy::Laplacian1D
    nx::Int
    ny::Int
    Lx::Float64
    Ly::Float64
    bc::NTuple{2,NTuple{2,BoundaryCondition}}
    allneumann::Bool
    function Laplacian2D(nx,ny,Lx,Ly,bcleft,bcright,bcbottom,bctop)
        bc=((bcleft,bcright),(bcbottom,bctop))
        allneumann = bc[1][1]==bc[1][1]==bc[1][1]==bc[1][1]==neumann
        allneumann && println("allneumann !!")
        new(Laplacian1D(nx,Lx,bcleft,bcright),Laplacian1D(ny,Ly,bcbottom,bctop),nx,ny,Lx,Ly,bc,allneumann)
    end

end
function push_element!(Is,Js,Vs,I,J,v)
    push!(Is,I)
    push!(Js,J)
    push!(Vs,v)
end

function SparseArrays.sparse(a::Laplacian2D)
    Is=Int[]
    Js=Int[]
    Vs=Float64[]
    lpx,lpy=a.lpx,a.lpy
    nx,ny=lpx.n_,lpy.n_

    Is=Int[]
    Js=Int[]
    Vs=Float64[]

    for i=1:nx
        for j=1:ny
            I=i+nx*(j-1)
            push_element!(Is,Js,Vs,I,I,lpx[i,i]+lpy[j,j])# sp[I,I]=Lx[i,i]

            i>1  && push_element!(Is,Js,Vs,I,I-1,lpx[i,i-1])#(sp[I,I-1]=Lx[i,i-1])
            i<nx && push_element!(Is,Js,Vs,I,I+1,lpx[i,i+1])#(sp[I,I+1]=Lx[i,i+1])
            j>1 &&  push_element!(Is,Js,Vs,I,I-ny,lpy[j,j-1])#(sp[I,I-ny]+=Ly[j,j-1])
            j<ny && push_element!(Is,Js,Vs,I,I+ny,lpy[j,j+1])#(sp[I,I+ny]+=Ly[j,j+1])

        end
    end
    sp=sparse(Is,Js,Vs)
    sp
end

function sparse_corr(a::Laplacian2D)
    sp=SparseArrays.sparse(a::Laplacian2D)
    if a.allneumann
        sp[1,1]=(3/2)sp[1][1]
    end
    sp
end


text:
__precompile__(true)

module OrthoMatchingPursuit

using LinearAlgebra
using ElasticArrays

include("qr.jl")
include("elasticqr.jl")
include("lsomp.jl")
include("pivotedqr.jl")
include("qromp.jl")
include("batchqromp.jl")
include("wrapper.jl")

end # module

max_forks_repo_name: mleprovost/QROMP


text:
using Solaris, Test, Flux, DiffEqFlux, Optim, Plots, Tracker

include("test_io.jl")
include("test_layer.jl")
include("test_train.jl")
include("test_widget.jl")


text:
export linear_constitutive_law
"""
    nn_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1}, grad_input::Int64, grad_θ::Int64)

Computes the stress with the given neural network weights θ.
# Inputs
`input` is a ``n\\times 9`` array, θ is a `p` dimensional array
`grad_input = 0`: do not compute gradients with respect to input; otherwise compute the gradients, in this case, `g` must have the same size as σ
`grad_θ = 0`    : do not compute gradients with respect to θ; otherwise compute the gradients

output:
- σ
- ``\\frac{\\partial \\sigma}{\\partial ipt}``
- ``\\frac{\\partial \\sigma\\cdot g}{\\partial \\theta}``
"""
#=
function linear_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1},
    g::Union{Array{Float64,2}, Nothing}=nothing, grad_input::Bool=false, grad_θ::Bool=false)

    @assert size(θ,1)==6

    H0 = [θ[1] θ[2] θ[3];
          θ[2] θ[4] θ[5];
          θ[3] θ[5] θ[6]]

    n = size(input,1)
    @assert size(input,2)==9

    ε = input[:, 1:3]

    σ = zeros(Float64, n, 3)
    for i = 1:n
        σ[i,:] = H0 * ε[i, :]
    end

    dinput = zeros()
    if(grad_input)
        dinput = zeros(Float64, n, 9, 3)
        for i = 1:n
            dinput[i,1:3,:] = H0'
        end

        grad_input = 1
    else
        grad_input = 0
    end

    dθ = zeros()
    if(grad_θ)
        @assert size(g,1)==n && size(g,2)==3

        dθ = zeros(length(θ))

        for i = 1:n

            dθ[1] += ε[i,1] * g[i,1]
            dθ[2] += ε[i,1] * g[i,2] + ε[i,2] * g[i,1]
            dθ[3] += ε[i,1] * g[i,3] + ε[i,3] * g[i,1]
            dθ[4] += ε[i,2] * g[i,2]
            dθ[5] += ε[i,2] * g[i,3] + ε[i,3] * g[i,2]
            dθ[6] += ε[i,3] * g[i,3]

        end

        grad_θ = 1
    else
        grad_θ = 0
    end

    return σ, dinput, dθ
end
=#

function linear_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1}, config::Array{Int64,1},
    g::Union{Array{Float64,2}, Nothing}=nothing, grad_input::Bool=false, grad_θ::Bool=false)

    @assert size(θ,1)==6

    H0 = [θ[1] θ[2] θ[3];
          θ[2] θ[4] θ[5];
          θ[3] θ[5] θ[6]]

    n = size(input,1)
    @assert size(input,2)==9

    ε ,ε0, σ0 = input[:, 1:3],input[:, 4:6],input[:, 7:9]

    σ = zeros(Float64, n, 3)
    for i = 1:n
        σ[i,:] = H0 * (ε[i, :] -  ε0[i, :]) +  σ0[i, :]
        #σ[i,:] = H0 * ε[i, :]  +  σ0[i, :]
    end

    I3 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    dinput = zeros()
    if(grad_input)
        dinput = zeros(Float64, n, 9, 3)
        for i = 1:n
            dinput[i,1:3,:] =  H0'
            dinput[i,4:6,:] = -H0'
            dinput[i,7:9,:] =  I3'
        end

        grad_input = 1
    else
        grad_input = 0
    end

    dθ = zeros()
    if(grad_θ)
        @assert size(g,1)==n && size(g,2)==3

        dθ = zeros(length(θ))

        for i = 1:n

            # dθ[1] += ε[i,1] * g[i,1]
            # dθ[2] += ε[i,1] * g[i,2] + ε[i,2] * g[i,1]
            # dθ[3] += ε[i,1] * g[i,3] + ε[i,3] * g[i,1]
            # dθ[4] += ε[i,2] * g[i,2]
            # dθ[5] += ε[i,2] * g[i,3] + ε[i,3] * g[i,2]
            # dθ[6] += ε[i,3] * g[i,3]

            dθ[1] += (ε[i,1]-ε0[i,1]) * g[i,1]
            dθ[2] += (ε[i,1]-ε0[i,1]) * g[i,2] + (ε[i,2]-ε0[i,2]) * g[i,1]
            dθ[3] += (ε[i,1]-ε0[i,1]) * g[i,3] + (ε[i,3]-ε0[i,3]) * g[i,1]
            dθ[4] += (ε[i,2]-ε0[i,2]) * g[i,2]
            dθ[5] += (ε[i,2]-ε0[i,2]) * g[i,3] + (ε[i,3]-ε0[i,3]) * g[i,2]
            dθ[6] += (ε[i,3]-ε0[i,3]) * g[i,3]

        end

        grad_θ = 1
    else
        grad_θ = 0
    end

    return σ, dinput, dθ
end

# function linear_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1},
#     g::Union{Array{Float64,2}, Nothing}=nothing, grad_input::Bool=false, grad_θ::Bool=false)

#     @assert size(θ,1)==9

#     H0 = [θ[1] θ[2] θ[3];
#           θ[4] θ[5] θ[6];
#           θ[7] θ[8] θ[9]]

#     n = size(input,1)
#     @assert size(input,2)==9

#     ε ,ε0, σ0 = input[:, 1:3],input[:, 4:6],input[:, 7:9]

#     σ = zeros(Float64, n, 3)
#     for i = 1:n
#         σ[i,:] = H0 * (ε[i, :] -  ε0[i, :]) +  σ0[i, :]
#         #σ[i,:] = H0 * ε[i, :]  +  σ0[i, :]
#     end

#     I3 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
#     dinput = zeros()
#     if(grad_input)
#         dinput = zeros(Float64, n, 9, 3)
#         for i = 1:n
#             dinput[i,1:3,:] =  H0'
#             dinput[i,4:6,:] = -H0'
#             dinput[i,7:9,:] =  I3'
#         end

#         grad_input = 1
#     else
#         grad_input = 0
#     end

#     dθ = zeros()
#     if(grad_θ)
#         @assert size(g,1)==n && size(g,2)==3

#         dθ = zeros(length(θ))

#         for i = 1:n

#             # dθ[1] += ε[i,1] * g[i,1]
#             # dθ[2] += ε[i,1] * g[i,2] + ε[i,2] * g[i,1]
#             # dθ[3] += ε[i,1] * g[i,3] + ε[i,3] * g[i,1]
#             # dθ[4] += ε[i,2] * g[i,2]
#             # dθ[5] += ε[i,2] * g[i,3] + ε[i,3] * g[i,2]
#             # dθ[6] += ε[i,3] * g[i,3]
#             for j = 1:3
#                 for k = 1:3
#                 #dθ[j, k]
#                 dθ[(j-1)*3 + k] += g[i,j] * ε[i,k]
#                 end
#             end

#         end

#         grad_θ = 1
#     else
#         grad_θ = 0
#     end

#     return σ, dinput, dθ
# end


text:
export AbstractSpace, israndomized

"""
    AbstractSpace{T, S}

General abstract type for space with `T` as its content
in state `S`.

In fact, an instance of a space should be a type with a
temporary memory to store an element of the space in type `T`.
"""
abstract type AbstractSpace{T, S<:SpaceState} end

israndomized(x::AbstractSpace{T, Randomized}) where {T} = true
israndomized(x::AbstractSpace{T, S}) where {T, S} = false

# Space Interface
export data, reset!, copy, acquire, shake!, randomize!

"""
    data(space::AbstractSpace{T, S}) -> T

get the current data of this space
"""
function data(space::AbstractSpace) end

"""
    reset!(space) -> space # with initilized state

reset the temporary data of this space to initial position
be careful that the state of the space will also be changed
assign it to your variable manually.

```julia
space = reset!(space)
```
"""
function reset!(space::AbstractSpace) end

"""
    copy(space) -> space

copy the space instance
"""
function copy(space::AbstractSpace) end

"""
    acquire(space) -> element

get a copy of current element of the space
"""
function acquire(space::AbstractSpace) end

"""
    shake!(space) -> space

shake the space and move current element to a random position.
the space has to be a space with `Randomized` state. or it will
an `UnRandomizedError`.
"""
function shake!(space::AbstractSpace) end

"""
    randomize!(space) -> space

randomized the space, move current element to a random position
and change the state of this space to `Initialized`
"""
function randomize!(space::AbstractSpace) end

######################################################
export AbstractSpaceTraverser, traverse

abstract type AbstractSpaceTraverser end

"""
    traverse(space) -> iterator

get an iterator that traverses the whole space.
"""
function traverse(space::AbstractSpace) end

################################
# Some default definitions
################################

reset!(space::AbstractSpace{T, UnRandomized}) where T = space
randomize!(space::AbstractSpace{T, Randomized}) where T = space

function shake!(space::AbstractSpace{T, UnRandomized}) where T
    throw(UnRandomizedError())
end


text:
struct Instance <: LavaAbstraction
    handle::Vk.Instance
    layers::Vector{String}
    extensions::Vector{String}
    messenger::Optional{Vk.DebugUtilsMessengerEXT}
end

vk_handle_type(::Type{Instance}) = Vk.Instance

Instance(handle::Vk.Instance) = Instance(handle, [], [], nothing)

function Instance(layers, extensions, messenger_info::Optional{Vk.DebugUtilsMessengerCreateInfoEXT} = nothing; application_info = C_NULL)
    handle = unwrap(create(Instance, Vk.InstanceCreateInfo(layers, extensions; application_info, next = something(messenger_info, C_NULL))))
    if !isnothing(messenger_info)
        messenger = debug_messenger(handle, messenger_info)
    else
        messenger = nothing
    end
    Instance(handle, layers, extensions, messenger)
end

const REQUIRED_DEVICE_EXTENSIONS = [
    "VK_KHR_synchronization2",
]

function init(;
    instance_layers = String[],
    instance_extensions = String[],
    application_info = Vk.ApplicationInfo(v"1", v"1", v"1.2"),
    device_extensions = String[],
    device_specific_features::AbstractVector{Symbol} = Symbol[],
    device_vulkan_features::AbstractVector{Symbol} = Symbol[],
    queue_config = dictionary([
        Vk.QUEUE_GRAPHICS_BIT | Vk.QUEUE_COMPUTE_BIT => 1,
        Vk.QUEUE_TRANSFER_BIT => 1,
    ]),
    with_validation = true,
    debug = true,
    message_types = Vk.DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT | Vk.DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT
)

    if with_validation && "VK_LAYER_KHRONOS_validation" ∉ instance_layers
        push!(instance_layers, "VK_LAYER_KHRONOS_validation")
    end
    if debug && "VK_EXT_debug_utils" ∉ instance_extensions
        push!(instance_extensions, "VK_EXT_debug_utils")
    end

    available_layers = unwrap(Vk.enumerate_instance_layer_properties())
    unsupported_layers = filter(!in(getproperty.(available_layers, :layer_name)), instance_layers)
    if !isempty(unsupported_layers)
        error("Requesting unsupported instance layers: $unsupported_layers")
    end

    available_extensions = unwrap(Vk.enumerate_instance_extension_properties())
    unsupported_extensions = filter(!in(getproperty.(available_extensions, :extension_name)), instance_extensions)
    if !isempty(unsupported_extensions)
        error("Requesting unsupported instance extensions: $unsupported_extensions")
    end

    dbg_info = if debug
        Vk.DebugUtilsMessengerCreateInfoEXT(
            |(
                Vk.DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
                Vk.DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
                Vk.DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
                Vk.DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
            ),
            message_types,
            debug_callback_c[],
        )
    else
        nothing
    end

    instance = Instance(instance_layers, instance_extensions, dbg_info; application_info)

    union!(device_vulkan_features, [:buffer_device_address, :descriptor_indexing, :descriptor_binding_partially_bound])
    synchronization_features = physical_device_features(Vk.PhysicalDeviceSynchronization2FeaturesKHR, [:synchronization2])
    vulkan_features = physical_device_features(Vk.PhysicalDeviceVulkan12Features, device_vulkan_features; next = synchronization_features)
    device_features = physical_device_features(Vk.PhysicalDeviceFeatures, device_specific_features)
    enabled_features = Vk.PhysicalDeviceFeatures2(device_features; next = vulkan_features)

    physical_device = pick_supported_device(unwrap(Vk.enumerate_physical_devices(instance)), enabled_features)

    available_extensions = unwrap(Vk.enumerate_device_extension_properties(physical_device))
    union!(device_extensions, REQUIRED_DEVICE_EXTENSIONS)
    unsupported_extensions = filter(!in(getproperty.(available_extensions, :extension_name)), device_extensions)
    if !isempty(unsupported_extensions)
        error("Requesting unsupported device extensions: $unsupported_extensions")
    end

    device = Device(physical_device, device_extensions, queue_config, enabled_features; next = enabled_features)

    instance, device
end

function physical_device_features(@nospecialize(T), features; next = C_NULL)
    names = fieldnames(T)
    unknown = filter(!in(names), features)
    if !isempty(unknown)
        error("Trying to set unknown features: $unknown")
    end
    fields = map(in(features), filter(≠(:next), names))
    if :next in names
        T(fields...; next)
    else
        T(fields...)
    end
end

function pick_supported_device(physical_devices, features)
    unsupported = nothing
    for pdevice in physical_devices
        # TODO: fix Vk.get_physical_device_features_2 in Vulkan.jl

        # initialize structure
        original = Vk._PhysicalDeviceFeatures2(
            Vk._PhysicalDeviceFeatures(ntuple(Returns(false), fieldcount(Vk.PhysicalDeviceFeatures))...);
            next = Vk._PhysicalDeviceVulkan12Features(ntuple(Returns(false), fieldcount(Vk.PhysicalDeviceVulkan12Features) - 1)...)
        )

        # build reference
        original_vk = Base.unsafe_convert(Vk.core.VkPhysicalDeviceFeatures2, original)
        ref = Ref(original_vk)

        # fill it with data
        GC.@preserve original Vk.core.vkGetPhysicalDeviceFeatures2(pdevice, ref)
        filled = ref[]

        # load next chain
        filled_next = unsafe_load(convert(Ptr{Vk.core.VkPhysicalDeviceVulkan12Features}, filled.pNext))

        # reconstruct wrapper structs (equivalent to `Vk.from_vk` that is defined for other structs, but by hand)
        vulkan_features = Vk.PhysicalDeviceVulkan12Features((getproperty(filled_next, name) for name in fieldnames(Vk.core.VkPhysicalDeviceVulkan12Features)[2:end])...)
        device_specific_features = Vk.PhysicalDeviceFeatures((getproperty(filled.features, name) for name in fieldnames(Vk.core.VkPhysicalDeviceFeatures))...)

        pdevice_features = Vk.PhysicalDeviceFeatures2(device_specific_features; next = vulkan_features)

        unsupported = unsupported_features(features, pdevice_features)
        isempty(unsupported) && return pdevice
    end
    throw("Physical device features $unsupported are required but not available on any device.")
end

function unsupported_features(requested::Vk.PhysicalDeviceFeatures2, available::Vk.PhysicalDeviceFeatures2)
    d = Dictionary{Symbol,Vector{Symbol}}()
    unsupported_vulkan = unsupported_features(requested.next, available.next)
    isempty(unsupported_vulkan) || insert!(d, :vulkan, unsupported_vulkan)
    unsupported_device = unsupported_features(requested.features, available.features)
    isempty(unsupported_device) || insert!(d, :device, unsupported_device)
    d
end

function unsupported_features(requested::T, available::T) where {T}
    filter(collect(fieldnames(T))) do name
        name == :next && return false
        getproperty(requested, name) && !getproperty(available, name)
    end
end

function debug_messenger(instance, info::Vk.DebugUtilsMessengerCreateInfoEXT)
    unwrap(create(Vk.DebugUtilsMessengerEXT, instance, info))
end


text:
mutable struct Ast2CnfResult
  """
  `identifier_dict` is a dictionary which maps original variable identifiers to variables in the resulting CNF formula.
  """
  identifier_dict::Dict{String, UInt}
  """
  Number of variables in the resulting formula.
  """
  variable_count::UInt
  only_left_to_right::Bool
  clauses::Vector{Clause}
end

function Ast2CnfResult(only_left_to_right::Bool)
  Ast2CnfResult(
    Dict{String, UInt}(),
    0,
    only_left_to_right,
    Clause[]
  )
end

"""
    ast2cnf(root::AstNode; <keyword arguments>)::Ast2CnfResult

Converts a parsed formula in NNF into a formula in CNF.
Returns `Ast2CnfResult`.

# Arguments
- `only_left_to_right::Bool=false`: only generate left to right implications in Tseytin encoding.

See also [`Ast2CnfResult`](@ref).
"""
function ast2cnf(root::AstNode; only_left_to_right=false)::Ast2CnfResult
  res = Ast2CnfResult(only_left_to_right)
  single_lit = ast2cnf!(res, root)

  if 0 == length(res.clauses)
    res.clauses = [Clause(single_lit)]
  end

  res
end

function identifier_to_variable!(res::Ast2CnfResult, identifier::String)::UInt
  get!(res.identifier_dict, identifier) do
    res.variable_count += 1
  end
end

function ast2cnf!(res::Ast2CnfResult, node::VariableNode)::Int
  identifier_to_variable!(res, node.identifier)
end

function ast2cnf!(res::Ast2CnfResult, node::NotNode)::Int
  -signed(identifier_to_variable!(res, node.identifier))
end

"""
    extract_binary_conjunctive(res::Ast2CnfResult, node::AstNode)::Tuple{Int, Int, Int}

Returns literals for LHS and RHS of `node` and the literal of `node` (in this order).
"""
function extract_binary_conjunctive!(res::Ast2CnfResult,
                                     node::AstNode)::Tuple{Int, Int, Int}
  lhs::Int = ast2cnf!(res, node.lhs)
  rhs::Int = ast2cnf!(res, node.rhs)
  clause_identifier = res.variable_count += 1

  lhs, rhs, clause_identifier
end

function ast2cnf!(res::Ast2CnfResult, node::AndNode)::UInt
  lhs::Int, rhs::Int, cl_id::Int = extract_binary_conjunctive!(res, node)

  push!(res.clauses, Clause(-cl_id, lhs))
  push!(res.clauses, Clause(-cl_id, rhs))
  if !res.only_left_to_right
    push!(res.clauses, Clause(cl_id, -lhs, -rhs))
  end

  cl_id
end

function ast2cnf!(res::Ast2CnfResult, node::OrNode)::UInt
  lhs::Int, rhs::Int, cl_id::Int = extract_binary_conjunctive!(res, node)

  push!(res.clauses, Clause(-cl_id, lhs, rhs))
  if !res.only_left_to_right
    push!(res.clauses, Clause(cl_id, -lhs))
    push!(res.clauses, Clause(cl_id, -rhs))
  end

  cl_id
end


text:
using CUDA #hide

# Assuming that the ProbabilisticCircuits Julia package has been installed with `julia -e 'using Pkg; Pkg.add("ProbabilisticCircuits")'`, we can start using it as follows.

using ProbabilisticCircuits

# ### Reasoning with manually constructed circuits

# We begin by creating three positive literals (boolean variables) and manually construct a probabilistic circuit that encodes a Naive Bayes (NB) distribution with the following form: `Pr(rain, rainbow, wet) = Pr(rain) * Pr(rainbow|rain) * Pr(wet|rain)`.
rain, rainbow, wet = pos_literals(ProbCircuit, 3)
rain_pos = (0.7 * rainbow + 0.3 * (-rainbow)) * (0.9 * wet + 0.1 * (-wet)) # Pr(rainbow|rain=1) * Pr(wet|rain=1)
rain_neg = (0.2 * rainbow + 0.8 * (-rainbow)) * (0.3 * wet + 0.7 * (-wet)) # Pr(rainbow|rain=0) * Pr(wet|rain=0)
circuit = 0.4 * (rain * rain_pos) + 0.6 * ((-rain) * rain_neg); # Pr(rain, rainbow, wet)

# Just like any probability distribution, we can evaluate the probabilistic circuit on various inputs. Note that since log probabilities are used in probabilistic circuits for numerical stability, we need to take exponent of the evaluation output to get the probabilities.
exp(circuit(true, true, true)) # Pr(rain=1, rainbow=1, wet=1)
#-
exp(circuit(true, false, false)) # Pr(rain=1, rainbow=0, wet=0)

# From the above examples, we see that it is less likely to rain if we do not see rainbows and the streets are not wet.

# The purpose of this package is to offer a unified tool for efficient learning and inference (i.e., answering probabilistic queries such as marginals and MAP) over probabilistic circuits, which subsume a large class of tractable probabilistic models. We first use the above manually constructed circuit to demonstrate several queries that can be answered efficiently. Similar to [logic circuits](https://github.com/Juice-jl/LogicCircuits.jl), answering the following queries require *decomposability* and *determinism*, which is already satisfied by construction:
isdecomposable(circuit) && isdeterministic(circuit)

# Decomposability allows us to compute marginal probabilities given partial evidence efficiently (linear time w.r.t. the circuit size). For example, we want to ask the probability of observing rainbows. That is, we want to marginalize out the variables rain and wet. This can be done by evaluating the circuit with partial evidence:
exp(circuit(missing, true, missing)) # Pr(rainbow=1)

# Being able to compute marginals immediately offers the ability to compute conditional probabilities. For example, to compute the probability of raining given rainbow=1 and wet=1, we simply take the quotient of Pr(rain=1, rainbow=1, wet=1) and Pr(rainbow=1, wet=1):
exp(circuit(true, true, true) - circuit(missing, true, true)) # Pr(rain=1|rainbow=1, wet=1)

# If we are additionally supplied with the structural property *determinism*, we can answer some more advanced queries. For example, we can to compute the maximum a posteriori (MAP) query of the distribution:
assignments, log_prob = MAP(circuit, [missing, missing, missing])
print("The MAP assignment of the circuit is (rain=$(assignments[1]), rainbow=$(assignments[2]), wet=$(assignments[3])), with probability $(exp(log_prob)).")

# Besides the above examples, ProbabilisticCircuits.jl provides functionalities for a wide variety of queries, which are detailed in [this manual](https://juice-jl.github.io/ProbabilisticCircuits.jl/stable/manual/queries/).

# ### Building complex circuit structures

# ProbabilisticCircuits.jl provides tools to compile classic Probabilistic Graphical Models (PGMs) and Tractable Probabilistic Models (TPMs) into probabilistic circuits efficiently. For example, we can compile a factor graph (FG) into a probabilistic circuit with one line of code:
fg = fromUAI(zoo_fg_file("asia.uai")) # Load example factor graph
fg_circuit = ProbCircuit(compile_factor_graph(fg)[1]) # Compile the FG to a PC
print("`fg_circuit` contains $(num_edges(fg_circuit)) edges and $(num_parameters(fg_circuit)) parameters.")

# ### Learning probabilistic circuits from data

# ProbabilisticCircuits.jl offers various parameter learning and structure learning algorithms. It further support mini-batch learning on both CPUs and GPUs, which makes learning large models from large datasets very efficient.

# We use the binarized MNIST dataset to demonstrate example probabilistic circuit learning functionalities.
train_data, valid_data, test_data = twenty_datasets("binarized_mnist");

# We start with learning the parameters of a *decomposable* and *deterministic* probabilistic circuit. We first load the structure of the circuit from file:
circuit = zoo_psdd("mnist.psdd")
print("The loaded circuit contains $(num_edges(circuit)) edges and $(num_parameters(circuit)) parameters.")
#-
print("Structural properties of the circuit: decomposability: $(isdecomposable(circuit)), determinism: $(isdeterministic(circuit)).")

# Given that the circuit is decomposable and deterministic, the maximum likelihood estimation (MLE) of its parameters is in closed-form. That is, we can learn the MLE parameters deterministically:
estimate_parameters(circuit, train_data; pseudocount = 0.1) #hide
t = @elapsed estimate_parameters(circuit, train_data; pseudocount = 0.1)
print("Learning the parameters on a CPU took $(t) seconds.")

# Optionally, we can use GPUs to speedup the learning process:
if CUDA.functional() #hide
estimate_parameters(circuit, train_data; pseudocount = 0.1, use_gpu = true) #hide
t = @elapsed estimate_parameters(circuit, train_data; pseudocount = 0.1, use_gpu = true)
print("Learning the parameters on a GPU took $(t) seconds.")
end #hide

# Note that the insignificant speedup is due to the fact that the circuit is too small to make full use of the GPU. For large circuits the speedup could be at least ~10x.

# After the learning process, we can evaluate the model on the validation/test dataset. Here we use average log-likelihood per sample as the metric (we again utilize GPUs for efficiency):
avg_ll = log_likelihood_avg(circuit, test_data)
print("The average test data log-likelihood is $(avg_ll).")

# Besides `estimate_parameters`, ProbabilisticCircuits.jl offers iterative parameter learning algorithms such as Expectation-Maximization (EM) (i.e., `estimate_parameters_em`) and Stochastic Gradient Descent (SGD) (i.e., `sgd_parameter_learning`).

# ProbabilisticCircuits.jl also offers functionalities for learning the circuit structure and parameters simultaneously. For example, the Strudel structure learning algorithm is implemented natively in the package, and can be used with a few lines of code:
circuit_strudel = learn_circuit(train_data; maxiter = 100, verbose = false)
avg_ll = log_likelihood_avg(circuit_strudel, test_data)
print("The learned circuit contains $(num_edges(circuit)) edges and $(num_parameters(circuit)) parameters.\n")
print("The average test data log-likelihood is $(avg_ll).")


text:
# ---
# title: 149. Max Points on a Line
# id: problem149
# author: Indigo
# date: 2021-01-17
# difficulty: Hard
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/max-points-on-a-line/description/>
# hidden: true
# ---
#
# Given _n_ points on a 2D plane, find the maximum number of points that lie on
# the same straight line.
#
# **Example 1:**
#
#
#
#     Input: [[1,1],[2,2],[3,3]]
#     Output: 3
#     Explanation:
#     ^
#     |
#     |        o
#     |     o
#     |  o
#     +------------->
#     0  1  2  3  4
#
#
# **Example 2:**
#
#
#
#     Input: [[1,1],[3,2],[5,3],[4,1],[2,3],[1,4]]
#     Output: 4
#     Explanation:
#     ^
#     |
#     |  o
#     |     o        o
#     |        o
#     |  o        o
#     +------------------->
#     0  1  2  3  4  5  6
#
#
# **NOTE:**  input types have been changed on April 15, 2019. Please reset to
# default code definition to get new method signature.
#
#
## @lc code=start
using LeetCode
function max_points(points::Vector{Vector{Int}})
    res = 0
    for i in 1:length(points)
        if res > length(points) - i
            return res
        end
        cnt = DefaultDict{Rational, Int}(1)
        for j in i+1:length(points)
            cnt[(points[i][2] - points[j][2]) // (points[i][1] - points[j][1])] += 1
        end
        res = max(res, maximum(p.second for p in cnt))
    end
    res
end
## @lc code=end


text:
using Lasso, FactCheck
DATADIR = joinpath(dirname(@__FILE__), "data")

# Test against fused lasso results from glmgen
lakehuron = readcsv(joinpath(DATADIR, "LakeHuron.csv"); header=true)[1][:, 3]
facts("FusedLasso") do
	for lambda in (10, 1, 0.1)
		context("λ = $(lambda)") do
			@fact coef(fit(FusedLasso, lakehuron, lambda)) => roughly(vec(readcsv(joinpath(DATADIR, "LakeHuron_lambda_$lambda.csv"))))
		end
	end
end


text:
"""
    AIMStep!(::Threaded, p::AIMProblem{N,T}, c::AIMCache{N,U}) where {N <: Unsigned, T <: Number, U <: Any}

Performs a single step of the AIM algorithm in parallel using threads:
1. The initial data arrays are not altered.
2. The previous arrays receive the values of the current arrays.
3. The results of the next step computed using the initial and current arrays. Results are stored in the buffer arrays.
4. The current arrays receive the contents of the buffer arrays.

# Input
- `::Threaded`: Instance of the `Threaded` object.
- `p::AIMProblem`: The problem data to use in the computation.
- `c::AIMCache`: The cache of arrays that corresponds to the problem p.

# Output
    nothing
"""
function AIMStep!(::Threaded, p::AIMProblem{N,T}, c::AIMCache{N,U}) where {N <: Unsigned, T <: Number, U <: Any}

    zeroN = zero(N)
    oneN = one(N)
    twoN = oneN + oneN

    # The currently computed coefficients will be the previous iteration coefficients
    # after the iteration is done, so we store the current values to the previous array
    @inbounds copy!(c.pcda, c.ccda)
    @inbounds copy!(c.pdda, c.cdda)

    # Apply the AIM formula to compute the coeficients writing the results to the
    # buffer arrays and reading data from all the other arrays.
    # The use of a buffer array as temporary storage to the coefficients allow for the
    # computation of the coefficients to be parallelized
    Threads.@threads for i in zeroN:(get_niter(p) - oneN)

        c_sum = zero(T)
        d_sum = zero(T)

        for k in zeroN:i
            @inbounds c_sum += c.icda[k + oneN] * c.ccda[i - k + oneN]
            @inbounds d_sum += c.idda[k + oneN] * c.ccda[i - k + oneN]
        end

        @inbounds c.bcda[i + oneN] = (i + oneN) * c.ccda[i + twoN] + c.cdda[i + oneN] + c_sum
        @inbounds c.bdda[i + oneN] = (i + oneN) * c.cdda[i + twoN] + d_sum
    end

    # Copy the data in the buffer arrays to the current arrays
    @inbounds copy!(c.ccda, c.bcda)
    @inbounds copy!(c.cdda, c.bdda)

    return nothing
end

"""
    AIMStep!(::Serial, p::AIMProblem{N,T}, c::AIMCache{N,U}) where {N <: Unsigned, T <: Number, U <: Any}

Performs a single step of the AIM algorithm sequentially:
1. The initial data arrays are not altered.
2. The previous arrays receive the values of the current arrays.
3. The results of the next step computed using the initial and current arrays. Results are stored in the buffer arrays.
4. The current arrays receive the contents of the buffer arrays.

# Input
- `::Serial`: Instance of `Serial` object.
- `p::AIMProblem`: The problem data to use in the computation.
- `c::AIMCache`: The cache of arrays that corresponds to the problem p.

# Output
    nothing
"""
function AIMStep!(::Serial, p::AIMProblem{N,T}, c::AIMCache{N,U}) where {N <: Unsigned, T <: Number, U <: Any}

    zeroN = zero(N)
    oneN = one(N)
    twoN = oneN + oneN

    # The currently computed coefficients will be the previous iteration coefficients
    # after the iteration is done, so we store the current values to the previous array
    @inbounds copy!(c.pcda, c.ccda)
    @inbounds copy!(c.pdda, c.cdda)

    # Apply the AIM formula to compute the coeficients writing the results to the
    # buffer arrays and reading data from all the other arrays.
    # The use of a buffer array as temporary storage to the coefficients allow for the
    # computation of the coefficients to be parallelized
    for i in zeroN:(get_niter(p) - oneN)

        c_sum = zero(T)
        d_sum = zero(T)

        for k in zeroN:i
            @inbounds c_sum += c.icda[k + oneN] * c.ccda[i - k + oneN]
            @inbounds d_sum += c.idda[k + oneN] * c.ccda[i - k + oneN]
        end

        @inbounds c.bcda[i + oneN] = (i + oneN) * c.ccda[i + twoN] + c.cdda[i + oneN] + c_sum
        @inbounds c.bdda[i + oneN] = (i + oneN) * c.cdda[i + twoN] + d_sum
    end

    # Copy the data in the buffer arrays to the current arrays
    @inbounds copy!(c.ccda, c.bcda)
    @inbounds copy!(c.cdda, c.bdda)

    return nothing
end


text:
using Logging, Test, HSM

struct HandledEvent <: HSM.AbstractHsmEvent
end

struct UnhandledEvent <: HSM.AbstractHsmEvent
end

"""
Handle Event tests state machine.
"""
mutable struct TestUnhandledEventMachine <: HSM.AbstractHsmState
    state_info::HSM.HsmStateInfo
    handled_event::Bool

    TestUnhandledEventMachine(parent_state) = new(HSM.HsmStateInfo(parent_state), false)
end

# TestUnhandledEventMachineState event handlers.
# This test deliberately does not handle all events in the root machine state.
# In practice, your root machine state must have `on_event()` event handlers
# for all possible state machine events, otherwise the machine may raise an
# exception, and end up in an unexpected state.
function HSM.on_event!(state::TestUnhandledEventMachine, event::HandledEvent)
    @debug "on_event!(TestUnhandledEventMachine, HandledEvent)"
    state.handled_event = true
    return true
end

@testset "handle_event!() -- HSM guards against unhandled events." begin
    machine = TestUnhandledEventMachine(nothing)

    # Unhandled event should throw an error and not set `handled_event`.
    @test_throws HSM.HsmUnhandledEventError HSM.handle_event!(machine, UnhandledEvent())
    @test !machine.handled_event

    # Handled event should not throw an error and set `handled_event`.
    HSM.handle_event!(machine, HandledEvent())
    @test machine.handled_event
end