text:
all_lines = readlines("day8.txt")
all_lines_vec = hcat(split.(all_lines," | ")...)

inputs = all_lines_vec[1,:]
outputs = all_lines_vec[2,:]

keep_track_ints = zeros(Int64, size(all_lines,1))

for jk in 1:size(all_lines,1)
    full_line = inputs[jk] * " " * outputs[jk]
    all_vales = hcat(split.(full_line," ")...)

    d = Matrix(undef, 10, 2)
    d[:,1] = collect(0:9)
    d[1+1,2] = join(sort(split(all_vales[length.(all_vales) .== 2][1],"")))
    d[4+1,2] = join(sort(split(all_vales[length.(all_vales) .== 4][1],"")))
    d[7+1,2] = join(sort(split(all_vales[length.(all_vales) .== 3][1],"")))
    d[8+1,2] = join(sort(split(all_vales[length.(all_vales) .== 7][1],"")))

    # Numbers created with 6 letters
    sixlengths = hcat(split.(all_vales[length.(all_vales) .== 6],"")...)
    new_6 = hcat([sort(sixlengths[:,i]) for i in 1:size(sixlengths,2)]...)

    # 4 is in 9
    fours = sort(split(d[4+1,2],""))

    nines_orig = sixlengths[:,[issubset(Set(fours), Set(new_6[:,ii])) for ii in 1:size(new_6,2)]]
    d[9+1,2] = join(sort(nines_orig[:,1]))

    # 1 is not in 6
    ones = sort(split(d[1+1,2],""))

    sixes_orig = sixlengths[:,[!issubset(Set(ones), Set(new_6[:,ii])) for ii in 1:size(new_6,2)]]
    d[6+1,2] = join(sort(sixes_orig[:,1]))

    # 4 is not in 0 and 1 is in 0
    zeros_orig = sixlengths[:,[!issubset(Set(fours), Set(new_6[:,ii])) && issubset(Set(ones), Set(new_6[:,ii])) for ii in 1:size(new_6,2)]]
    d[0+1,2] = join(sort(zeros_orig[:,1]))

    # Numbers created with 5 letters
    fivelengths = hcat(split.(all_vales[length.(all_vales) .== 5],"")...)
    new_5 = hcat([sort(fivelengths[:,i]) for i in 1:size(fivelengths,2)]...)

    # 1 is in 3
    threes_orig = fivelengths[:,[issubset(Set(ones), Set(new_5[:,ii])) for ii in 1:size(new_5,2)]]
    d[3+1,2] = join(sort(threes_orig[:,1]))

    # complement of 7 and 9 is in 5
    comple_sev_nine = setdiff(Set(split(d[9+1,2],"")), Set(split(d[7+1,2],"")))

    fives_orig = fivelengths[:,[issubset(Set(comple_sev_nine), Set(new_5[:,ii])) for ii in 1:size(new_5,2)]]
    d[5+1,2] = join(sort(fives_orig[:,1]))

    # complement of 6 and 4 is in 2
    comple_six_four = setdiff(Set(split(d[6+1,2],"")), Set(split(d[4+1,2],"")))

    twos_orig = fivelengths[:,[issubset(Set(comple_six_four), Set(new_5[:,ii])) for ii in 1:size(new_5,2)]]
    d[2+1,2] = join(sort(twos_orig[:,1]))

    checks = join.(sort.(split.(split(outputs[jk]," "),"")))

    keep_track_ints[jk] =  parse(Int64, join([findall(x->x==checks[ii], d)[1][1] - 1 for ii in 1:size(checks,1)]))
end

sum(keep_track_ints)


text:
# Normal strings can't be type parameters hence this
struct StaticString{N}
    chars::NTuple{N, Char}
end
StaticString(s::AbstractString) = StaticString((s...,))
Base.string(s::StaticString) = *(s.chars...,)

_static(s::Vector) = _static(Tuple(s))
_static(s::Tuple) = _static.(s)
_static(s::AbstractString) = StaticString(s)
_static(s::Any) = s

Base.show(io::IO, s::StaticString) = show(io, string(s))
Base.println(s::StaticString) = println(string(s))
Base.print(s::StaticString) = print(string(s))

struct StaticSymbol{N, C<:NTuple{N, Any}}
    elements::C
end
StaticSymbol(s::Symbol) = StaticSymbol((String(s)...,))
Base.Symbol(s::StaticSymbol) = Symbol(s.elements...,)

_static(s::Symbol) = StaticSymbol(s)

Base.show(io::IO, s::StaticSymbol) = show(io, Symbol(s))
Base.println(s::StaticSymbol) = println(Symbol(s))
Base.print(s::StaticSymbol) = print(Symbol(s))

# Keys are type parameters
struct StaticDict{K, V} <: AbstractDict{K,V}
    values::V
end
StaticDict(ps::Base.Generator) = StaticDict(ps...)
StaticDict(ps::Pair...) = StaticDict(ps)
function StaticDict(ps::Tuple{Vararg{Pair}})
    keys = _static(ntuple(i -> ps[i].first, Val(length(ps))))
    values = ntuple(i -> ps[i].second, Val(length(ps)))
    return StaticDict{keys, typeof(values)}(values)
end
StaticDict(keys::Tuple, values::Tuple) = StaticDict{_static(keys), typeof(values)}(values)

Base.Dict(d::StaticDict) = Dict(k => v for (k, v) in d)
@generated function Base.getindex(d::StaticDict{keys}, k) where {keys}
    expr = Expr(:block)
    for (i, _k) in enumerate(QuoteNode.(keys))
        push!(expr.args, quote
            k === $_k && return d.values[$i]
        end)
    end
    push!(expr.args, :(throw(ArgumentError("Key $k not found."))))
    return expr
end

Base.keys(::StaticDict{keys}) where {keys} = keys
Base.values(d::StaticDict) = d.values

Base.length(d::StaticDict) = length(d.values)
Base.eltype(d::StaticDict) = eltype(d.values)
function Base.iterate(d::StaticDict, state=iterate(keys(d)))
    state !== nothing && return (state[1] => d[state[1]]), iterate(keys(d), state[2])
    return nothing
end

Base.show(io::IO, d::StaticDict) = show(io, Dict(d))

export StaticDict

ext: jl


text:
module env

# import Base:
#     ENV,
#     string,
#     _hasenv,
#     length,
#     show

# export

get(key::AbstractString) = has(key) ? ENV[key] : ""
get(key::Symbol) = has(string(key)) ? ENV[string(key)] : ""
set(key::AbstractString, value::AbstractString) = ENV[key] = value
set(key::Symbol, value::Symbol) = ENV[string(key)] = string(value)
has(key::AbstractString) = _hasenv(key)
has(key::Symbol) = _hasenv(string(key))
length() = length(ENV)
show() = show(ENV)

export get, set, has, length, show

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# to run the script, you should have these packages installed
using Plots
using StatsPlots
using BenchmarkTools
using Statistics
using RollingFunctions
using ImageFiltering
using SortFilters

y = randn(1000000)
windows = [5, 15, 51, 151]
t_srt = zeros(size(windows))
t_img = zeros(size(windows))
t_rol = zeros(size(windows))

for i in eachindex(windows)
    w = windows[i]

    t = @benchmark SortFilters.movsort($y, $w, 0.5) samples=3
    t_srt[i] = median(t).time / 10^6

    t = @benchmark ImageFiltering.mapwindow($median, $y, $w) samples=3
    t_img[i] = median(t).time / 10^6

    t = @benchmark RollingFunctions.rollmedian($y, $w) samples=3
    t_rol[i] = median(t).time / 10^6
end

ctg = repeat(["1 - SortFilters . movsort",
              "2 - ImageFiltering . mapwindow median",
              "3 - RollingFunctions . rollmedian",
              ], inner = length(windows))

#bars = log10.(hcat(t_mfl, t_srt, t_img, t_rol, t_ind)) # log scale
bars = hcat(t_srt, t_img, t_rol)
xnames = repeat(1:length(windows), outer = size(bars,2))

groupedbar(xnames, bars, group = ctg, xlabel = "Window lengths", ylabel = "Time, ms",
        title = "Moving median performance, array length = 1000000.", bar_width = 0.7,
        lw = 0, framestyle = :box, legend=:topleft, ylims = (0, 5200))

xticks!(1:4, string.(windows))
#yticks!(collect(0:1:3), "10^{" .* string.(collect(0:1:3)) .* "}") # log scale

png("plot.png")


text:
@testset "PeakResults" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 4])
    @test pr.indices == [0, 9]
    @test pr.peaks == [1, 2]
    @test pr.locations == [3, 4]
    @test pr.prominences == [5, 6]
    @test pr.width_bounds == [1 2; 3 4]

end

@testset "peakheights" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 8])
    peaks = peakheights(pr)
    @test peaks == [1, 2]

end

@testset "peaklocations" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 8])
    locations = peaklocations(pr)
    @test locations == [3, 4]

end

@testset "peakprominences" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 8])
    prominences = peakprominences(pr)
    @test prominences == [5, 6]

end

@testset "peakwidths" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 8])
    widths = peakwidths(pr)
    @test widths == [1, 5]

end

@testset "isempty" begin

    pr = PeakResults([0, 9], [1, 2], [3, 4], [5, 6], [1 2; 3 4])
    @test !isempty(pr)
    @test length(pr) == 2

    pr = PeakResults([], [], [], [], reshape([], 0, 2))
    @test isempty(pr)
    @test length(pr) == 0

end

@testset "getindex" begin

    indices = [0, 2, 4, 6]
    peaks = [1, 3, 5, 7]
    locations = [10, 11, 12, 13]
    prominences = [21, 22, 23, 24]
    width_bounds = [1 1; 2 2; 3 3; 4 4]
    pr = PeakResults(indices, peaks, locations, prominences, width_bounds)
    pr_select = pr[[2, 4]]
    @test pr_select.indices == [2, 6]
    @test pr_select.peaks == [3, 7]
    @test pr_select.locations == [11, 13]
    @test pr_select.prominences == [22, 24]
    @test pr_select.width_bounds == [2 2; 4 4]

end

@testset "sort" begin

    indices = [2, 0, 6, 4]
    peaks = [5, 3, 2, 7]
    locations = [11, 10, 13, 12]
    prominences = [21, 25, 23, 22]
    width_bounds = [1 1; 2 2; 3 3; 4 4]
    pr = PeakResults(indices, peaks, locations, prominences, width_bounds)

    pr_sort = sort(pr)
    @test pr_sort.indices == [6, 0, 2, 4]
    @test pr_sort.peaks == [2, 3, 5, 7]
    @test pr_sort.locations == [13, 10, 11, 12]
    @test pr_sort.prominences == [23, 25, 21, 22]
    @test pr_sort.width_bounds == [3 3; 2 2; 1 1; 4 4]

    pr_sort_rev = sort(pr; rev=true)
    @test pr_sort_rev.indices == [4, 2, 0, 6]
    @test pr_sort_rev.peaks == [7, 5, 3, 2]
    @test pr_sort_rev.locations == [12, 11, 10, 13]
    @test pr_sort_rev.prominences == [22, 21, 25, 23]
    @test pr_sort_rev.width_bounds == [4 4; 1 1; 2 2; 3 3]

    pr_sort_prom = sort(pr; ref="prominence")
    @test pr_sort_prom.indices == [2, 4, 6, 0]
    @test pr_sort_prom.prominences == [21, 22, 23, 25]

    pr_sort_prom_rev = sort(pr; rev=true, ref="prominence")
    @test pr_sort_prom_rev.indices == [0, 6, 4, 2]
    @test pr_sort_prom_rev.prominences == [25, 23, 22, 21]

    pr_sort_loc = sort(pr; ref="location")
    @test pr_sort_loc.indices == [0, 2, 4, 6]
    @test pr_sort_loc.prominences == [25, 21, 22, 23]

    pr_sort_loc_rev = sort(pr; rev=true, ref="location")
    @test pr_sort_loc_rev.indices == [6, 4, 2, 0]
    @test pr_sort_loc_rev.prominences == [23, 22, 21, 25]

    @test_throws ArgumentError sort(pr; ref="unknown")

end

@testset "plot" begin

    data = [1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1]
    pr = findpeaks(data)
    p = plot(pr, data)
    @test p.n == 8

end

max_stars_repo_name: vkumpost/FindPeaks

max_issues_repo_name: vkumpost/FindPeaks

max_forks_repo_name: vkumpost/FindPeaks


text:
export BrickletThermalImagingHighContrastImageLowLevel
struct BrickletThermalImagingHighContrastImageLowLevel
    image_chunk_offset::Integer
    image_chunk_data::Vector{Integer}
end

export BrickletThermalImagingTemperatureImageLowLevel
struct BrickletThermalImagingTemperatureImageLowLevel
    image_chunk_offset::Integer
    image_chunk_data::Vector{Integer}
end

export BrickletThermalImagingStatistics
struct BrickletThermalImagingStatistics
    spotmeter_statistics::Vector{Integer}
    temperatures::Vector{Integer}
    resolution::Integer
    ffc_status::Integer
    temperature_warning::Vector{Bool}
end

export BrickletThermalImagingHighContrastConfig
struct BrickletThermalImagingHighContrastConfig
    region_of_interest::Vector{Integer}
    dampening_factor::Integer
    clip_limit::Vector{Integer}
    empty_counts::Integer
end

export BrickletThermalImagingFluxLinearParameters
struct BrickletThermalImagingFluxLinearParameters
    scene_emissivity::Integer
    temperature_background::Integer
    tau_window::Integer
    temperatur_window::Integer
    tau_atmosphere::Integer
    temperature_atmosphere::Integer
    reflection_window::Integer
    temperature_reflection::Integer
end

export BrickletThermalImagingSPITFPErrorCount
struct BrickletThermalImagingSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletThermalImagingIdentity
struct BrickletThermalImagingIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletThermalImaging
"""
80x60 pixel thermal imaging camera
"""
mutable struct BrickletThermalImaging <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletThermalImaging(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_thermal_imaging")
        deviceInternal = package.BrickletThermalImaging(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_high_contrast_image_low_level
"""
    $(SIGNATURES)

Returns the current high contrast image. See `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Thermal_Imaging.html#high-contrast-image-vs-temperature-image>`__
for the difference between
High Contrast and Temperature Image. If you don't know what to use
the High Contrast Image is probably right for you.

The data is organized as a 8-bit value 80x60 pixel matrix linearized in
a one-dimensional array. The data is arranged line by line from top left to
bottom right.

Each 8-bit value represents one gray-scale image pixel that can directly be
shown to a user on a display.

Before you can use this function you have to enable it with
:func:`Set Image Transfer Config`.
"""
function get_high_contrast_image_low_level(device::BrickletThermalImaging)
    return device.deviceInternal.get_high_contrast_image_low_level()
end

export get_temperature_image_low_level
"""
    $(SIGNATURES)

Returns the current temperature image. See `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Thermal_Imaging.html#high-contrast-image-vs-temperature-image>`__
for the difference between High Contrast and Temperature Image.
If you don't know what to use the High Contrast Image is probably right for you.

The data is organized as a 16-bit value 80x60 pixel matrix linearized in
a one-dimensional array. The data is arranged line by line from top left to
bottom right.

Each 16-bit value represents one temperature measurement in either
Kelvin/10 or Kelvin/100 (depending on the resolution set with :func:`Set Resolution`).

Before you can use this function you have to enable it with
:func:`Set Image Transfer Config`.
"""
function get_temperature_image_low_level(device::BrickletThermalImaging)
    return device.deviceInternal.get_temperature_image_low_level()
end

export get_statistics
"""
    $(SIGNATURES)

Returns the spotmeter statistics, various temperatures, current resolution and status bits.

The spotmeter statistics are:

* Index 0: Mean Temperature.
* Index 1: Maximum Temperature.
* Index 2: Minimum Temperature.
* Index 3: Pixel Count of spotmeter region of interest.

The temperatures are:

* Index 0: Focal Plain Array temperature.
* Index 1: Focal Plain Array temperature at last FFC (Flat Field Correction).
* Index 2: Housing temperature.
* Index 3: Housing temperature at last FFC.

The resolution is either `0 to 6553 Kelvin` or `0 to 655 Kelvin`. If the resolution is the former,
the temperatures are in Kelvin/10, if it is the latter the temperatures are in Kelvin/100.

FFC (Flat Field Correction) Status:

* FFC Never Commanded: Only seen on startup before first FFC.
* FFC Imminent: This state is entered 2 seconds prior to initiating FFC.
* FFC In Progress: Flat field correction is started (shutter moves in front of lens and back). Takes about 1 second.
* FFC Complete: Shutter is in waiting position again, FFC done.

Temperature warning bits:

* Index 0: Shutter lockout (if true shutter is locked out because temperature is outside -10°C to +65°C)
* Index 1: Overtemperature shut down imminent (goes true 10 seconds before shutdown)
"""
function get_statistics(device::BrickletThermalImaging)
    return device.deviceInternal.get_statistics()
end

export set_resolution
"""
    $(SIGNATURES)

Sets the resolution. The Thermal Imaging Bricklet can either measure

* from 0 to 6553 Kelvin (-273.15°C to +6279.85°C) with 0.1°C resolution or
* from 0 to 655 Kelvin (-273.15°C to +381.85°C) with 0.01°C resolution.

The accuracy is specified for -10°C to 450°C in the
first range and -10°C and 140°C in the second range.
"""
function set_resolution(device::BrickletThermalImaging, resolution)
    device.deviceInternal.set_resolution(resolution)
end

export get_resolution
"""
    $(SIGNATURES)

Returns the resolution as set by :func:`Set Resolution`.
"""
function get_resolution(device::BrickletThermalImaging)
    return device.deviceInternal.get_resolution()
end

export set_spotmeter_config
"""
    $(SIGNATURES)

Sets the spotmeter region of interest. The 4 values are

* Index 0: Column start (has to be smaller then Column end).
* Index 1: Row start (has to be smaller then Row end).
* Index 2: Column end (has to be smaller then 80).
* Index 3: Row end (has to be smaller then 60).

The spotmeter statistics can be read out with :func:`Get Statistics`.
"""
function set_spotmeter_config(device::BrickletThermalImaging, region_of_interest)
    device.deviceInternal.set_spotmeter_config(region_of_interest)
end

export get_spotmeter_config
"""
    $(SIGNATURES)

Returns the spotmeter config as set by :func:`Set Spotmeter Config`.
"""
function get_spotmeter_config(device::BrickletThermalImaging)
    return device.deviceInternal.get_spotmeter_config()
end

export set_high_contrast_config
"""
    $(SIGNATURES)

Sets the high contrast region of interest, dampening factor, clip limit and empty counts.
This config is only used in high contrast mode (see :func:`Set Image Transfer Config`).

The high contrast region of interest consists of four values:

* Index 0: Column start (has to be smaller or equal then Column end).
* Index 1: Row start (has to be smaller then Row end).
* Index 2: Column end (has to be smaller then 80).
* Index 3: Row end (has to be smaller then 60).

The algorithm to generate the high contrast image is applied to this region.

Dampening Factor: This parameter is the amount of temporal dampening applied to the HEQ
(history equalization) transformation function. An IIR filter of the form::

 (N / 256) * previous + ((256 - N) / 256) * current

is applied, and the HEQ dampening factor
represents the value N in the equation, i.e., a value that applies to the amount of
influence the previous HEQ transformation function has on the current function. The
lower the value of N the higher the influence of the current video frame whereas
the higher the value of N the more influence the previous damped transfer function has.

Clip Limit Index 0 (AGC HEQ Clip Limit High): This parameter defines the maximum number of pixels allowed
to accumulate in any given histogram bin. Any additional pixels in a given bin are clipped.
The effect of this parameter is to limit the influence of highly-populated bins on the
resulting HEQ transformation function.

Clip Limit Index 1 (AGC HEQ Clip Limit Low): This parameter defines an artificial population that is added to
every non-empty histogram bin. In other words, if the Clip Limit Low is set to L, a bin
with an actual population of X will have an effective population of L + X. Any empty bin
that is nearby a populated bin will be given an artificial population of L. The effect of
higher values is to provide a more linear transfer function; lower values provide a more
non-linear (equalized) transfer function.

Empty Counts: This parameter specifies the maximum number of pixels in a bin that will be
interpreted as an empty bin. Histogram bins with this number of pixels or less will be
processed as an empty bin.
"""
function set_high_contrast_config(device::BrickletThermalImaging, region_of_interest, dampening_factor, clip_limit, empty_counts)
    device.deviceInternal.set_high_contrast_config(region_of_interest, dampening_factor, clip_limit, empty_counts)
end

export get_high_contrast_config
"""
    $(SIGNATURES)

Returns the high contrast config as set by :func:`Set High Contrast Config`.
"""
function get_high_contrast_config(device::BrickletThermalImaging)
    return device.deviceInternal.get_high_contrast_config()
end

export set_image_transfer_config
"""
    $(SIGNATURES)

The necessary bandwidth of this Bricklet is too high to use getter/callback or
high contrast/temperature image at the same time. You have to configure the one
you want to use, the Bricklet will optimize the internal configuration accordingly.

Corresponding functions:

* Manual High Contrast Image: :func:`Get High Contrast Image`.
* Manual Temperature Image: :func:`Get Temperature Image`.
* Callback High Contrast Image: :cb:`High Contrast Image` callback.
* Callback Temperature Image: :cb:`Temperature Image` callback.
"""
function set_image_transfer_config(device::BrickletThermalImaging, config)
    device.deviceInternal.set_image_transfer_config(config)
end

export get_image_transfer_config
"""
    $(SIGNATURES)

Returns the image transfer config, as set by :func:`Set Image Transfer Config`.
"""
function get_image_transfer_config(device::BrickletThermalImaging)
    return device.deviceInternal.get_image_transfer_config()
end

export set_flux_linear_parameters
"""
    $(SIGNATURES)

Sets the flux linear parameters that can be used for radiometry calibration.

See FLIR document 102-PS245-100-01 for more details.

.. versionadded:: 2.0.5\$nbsp;(Plugin)
"""
function set_flux_linear_parameters(device::BrickletThermalImaging, scene_emissivity, temperature_background, tau_window, temperatur_window, tau_atmosphere, temperature_atmosphere, reflection_window, temperature_reflection)
    device.deviceInternal.set_flux_linear_parameters(scene_emissivity, temperature_background, tau_window, temperatur_window, tau_atmosphere, temperature_atmosphere, reflection_window, temperature_reflection)
end

export get_flux_linear_parameters
"""
    $(SIGNATURES)

Returns the flux linear parameters, as set by :func:`Set Flux Linear Parameters`.

.. versionadded:: 2.0.5\$nbsp;(Plugin)
"""
function get_flux_linear_parameters(device::BrickletThermalImaging)
    return device.deviceInternal.get_flux_linear_parameters()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletThermalImaging)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletThermalImaging, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletThermalImaging)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletThermalImaging, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletThermalImaging, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletThermalImaging, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletThermalImaging)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletThermalImaging)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletThermalImaging)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletThermalImaging, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletThermalImaging)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletThermalImaging)
    return device.deviceInternal.get_identity()
end

export get_high_contrast_image
"""
Returns the current high contrast image. See `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Thermal_Imaging.html#high-contrast-image-vs-temperature-image>`__
for the difference between
High Contrast and Temperature Image. If you don't know what to use
the High Contrast Image is probably right for you.

The data is organized as a 8-bit value 80x60 pixel matrix linearized in
a one-dimensional array. The data is arranged line by line from top left to
bottom right.

Each 8-bit value represents one gray-scale image pixel that can directly be
shown to a user on a display.

Before you can use this function you have to enable it with
:func:`Set Image Transfer Config`.
"""
function get_high_contrast_image(device::BrickletThermalImaging)

    image_length = 4800

    lock(device.stream_lock) do
        ret = get_high_contrast_image_low_level(device, )

        if ret.image_chunk_offset == (1 << 16) - 1 # maximum chunk offset -> stream has no data
            image_length = 0
            image_out_of_sync = false
            image_data = ()
        else
            image_out_of_sync = ret.image_chunk_offset != 0
        end

        image_data = ret.image_chunk_data

        while !image_out_of_sync && length(image_data) < image_length
            ret = get_high_contrast_image_low_level(device, )
            image_out_of_sync = ret.image_chunk_offset != length(image_data)
            image_data += ret.image_chunk_data
        end

        if image_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.image_chunk_offset + 62 < image_length
                ret = get_high_contrast_image_low_level(device, )
            end

            throw(TinkerforgeStreamOutOfSyncError("Image stream is out-of-sync"))
        end
    end

        return image_data[:image_length]
end

export get_temperature_image
"""
Returns the current temperature image. See `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Thermal_Imaging.html#high-contrast-image-vs-temperature-image>`__
for the difference between High Contrast and Temperature Image.
If you don't know what to use the High Contrast Image is probably right for you.

The data is organized as a 16-bit value 80x60 pixel matrix linearized in
a one-dimensional array. The data is arranged line by line from top left to
bottom right.

Each 16-bit value represents one temperature measurement in either
Kelvin/10 or Kelvin/100 (depending on the resolution set with :func:`Set Resolution`).

Before you can use this function you have to enable it with
:func:`Set Image Transfer Config`.
"""
function get_temperature_image(device::BrickletThermalImaging)

    image_length = 4800

    lock(device.stream_lock) do
        ret = get_temperature_image_low_level(device, )

        if ret.image_chunk_offset == (1 << 16) - 1 # maximum chunk offset -> stream has no data
            image_length = 0
            image_out_of_sync = false
            image_data = ()
        else
            image_out_of_sync = ret.image_chunk_offset != 0
        end

        image_data = ret.image_chunk_data

        while !image_out_of_sync && length(image_data) < image_length
            ret = get_temperature_image_low_level(device, )
            image_out_of_sync = ret.image_chunk_offset != length(image_data)
            image_data += ret.image_chunk_data
        end

        if image_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.image_chunk_offset + 31 < image_length
                ret = get_temperature_image_low_level(device, )
            end

            throw(TinkerforgeStreamOutOfSyncError("Image stream is out-of-sync"))
        end
    end

        return image_data[:image_length]
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletThermalImaging, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
# # Phase recovery using MaxCut
#
# In this example, we relax the phase retrieval problem similar to the classical
# [MaxCut](http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf) semidefinite
# program and recover the phase of the signal given the magnitude of the linear
# measurements.
#
# Phase recovery has wide applications such as  in X-ray and crystallography
# imaging, diffraction imaging or microscopy and audio signal processing. In all
# these applications, the detectors cannot measure the phase of the incoming wave
# and only record its amplitude i.e complex measurements of a signal
# $x \in \mathbb{C}^p$ are obtained from a linear injective operator $A$, **but we
# can only measure the magnitude vector $Ax$, not the phase of $Ax$**.
#
# Recovering the phase of $Ax$ from $|Ax|$ is a **nonconvex optimization problem**. Using results from [this paper](https://arxiv.org/abs/1206.0102), the problem can be relaxed to a (complex) semidefinite program (complex SDP).
#
# The original reprsentation of the problem is as follows:
#
# $$
# \begin{array}{ll}
#   \text{find} & x \in \mathbb{C}^p \\
#     \text{subject to} & |Ax| = b
# \end{array}
# $$
#
# where $A \in \mathbb{C}^{n \times p}$ and $b \in \mathbb{R}^n$.

# In this example, **the problem is to find the phase of $Ax$ given the value $|Ax|$**.
# Given a linear operator $A$ and a vector $b= |Ax|$ of measured amplitudes,
# in the noiseless case, we can write $Ax = \text{diag}(b)u$ where
# $u \in \mathbb{C}^n$  is a phase vector, satisfying
# $|\mathbb{u}_i| = 1$ for $i = 1,\ldots, n$.
#
# We relax this problem as Complex Semidefinite Programming.
#
# ### Relaxed Problem similar to [MaxCut](http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf)
#
# Define the positive semidefinite hermitian matrix
# $M = \text{diag}(b) (I - A A^*) \text{diag}(b)$. The problem is:
#
# $$
# \begin{array}{ll}
#   \text{minimize} & \langle U, M \rangle \\
#     \text{subject to} & \text{diag}(U) = 1\\
#     & U \succeq 0
# \end{array}
# $$
#
# Here the variable $U$ must be hermitian ($U \in \mathbb{H}_n $),
# and we have a solution to the phase recovery problem if $U = u u^*$
# has rank one. Otherwise, the leading singular vector of $U$ can be used
# to approximate the solution.

using Convex, SCS, LinearAlgebra

n = 20
p = 2
A = rand(n, p) + im * randn(n, p)
x = rand(p) + im * randn(p)
b = abs.(A * x) + rand(n)

M = diagm(b) * (I(n) - A * A') * diagm(b)
U = ComplexVariable(n, n)
objective = inner_product(U, M)
c1 = diag(U) == 1
c2 = U in :SDP
p = minimize(objective, c1, c2)
solve!(p, SCS.Optimizer; silent_solver = true)
evaluate(U)

#-

# Verify if the rank of $U$ is 1:
B, C = eigen(evaluate(U));
length([e for e in B if (abs(real(e)) > 1e-4)])

#-

# Decompose $U = uu^*$ where $u$ is the phase of $Ax$
u = C[:, 1];
for i in 1:n
    u[i] = u[i] / abs(u[i])
end
u


text:
# https://juliabyexample.helpmanual.io/

# simple operations
number1 = 50
number2 = 80

sum = number1 + number2

println(sum)

max_stars_repo_name: ldaniel/Practicing-Julia

max_issues_repo_name: ldaniel/Practicing-Julia

max_forks_repo_name: ldaniel/Practicing-Julia


text:
export IsothermalEnergy

mutable struct IsothermalEnergy <: EnergyComponent{Nothing}
    specific_energies::Vector{Float}

    """
        IsothermalEnergy(
            positions::AbstractVector{SVectorF{N}},
            mass_densities::AbstractVector{Float},
            energy_distribution::InitialEnergyDistribution{N},
            eos::EquationOfState,
        )
    """
    function IsothermalEnergy(
        positions::AbstractVector{SVectorF{N}},
        mass_densities::AbstractVector{Float},
        energy_distribution::InitialEnergyDistribution{N},
        eos::EquationOfState,
    ) where {N}
        specific_energies = energy_distribution(positions, mass_densities, eos)
        new(specific_energies)
    end
end

get_specific_energies(isothermal_energy::IsothermalEnergy) =
    isothermal_energy.specific_energies
initderivatives(::IsothermalEnergy) = nothing
function initderivatives!(::Nothing, ::IsothermalEnergy) end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using DataFrames, CSV

df = CSV.read("couleur.tsv", DataFrame)

unique(df[!, "genre"])

replacements = Dict([item => counter for (counter, item) in enumerate(unique(df[!, "genre"]))]...)

df[!, "genre"] = map(x -> replacements[x], df[!, "genre"])

df[df[!, "genre"] .===2, :]

replace!(df[!, "genre"], missing => mean(filter(x -> x !== missing, df[!,"genre"])))

max_stars_repo_name: monkeyusage/brownbaglunch

max_issues_repo_name: monkeyusage/brownbaglunch

max_forks_repo_name: monkeyusage/brownbaglunch


text:
# Abstract types
export AbstractNode, Value, ValueType, AbstractVariable,
       AbstractArrayVariable, AbstractMatrixVariable, AbstractVectorVariable
# builtin concrete types
export Variable, Node, CachedNode, forward, gradient, backward, value, args, arg, operator
export register

using LinearAlgebra
# export register!

"""
    Operator

Abstract type for operators in the computation graph.
"""
abstract type Operator end

"""
    Trait

This module contains function traits as a subtype of [`Operator`](@ref).
"""
module Trait
import YAAD: Operator

"""
    Method{FT} <: Operator

This trait wraps a callable object in Julia (usually a `Function`).
"""
struct Method{FT} <: Operator
    f::FT
end

(op::Method)(args...; kwargs...) = op.f(args...; kwargs...)

"""
    Broadcasted{FT} <: Operator

This trait wraps a callable object that being broadcasted. It will help to
dispatch different gradient methods overloaded for broadcasted operation
comparing to [`Method`](@ref).
"""
struct Broadcasted{FT} <: Operator
    f::FT
end

(op::Broadcasted)(args...; kwargs...) = op.f.(args...; kwargs...)

end # Trait

"""
    AbstractNode

Abstract type for nodes in computation graph.
"""
abstract type AbstractNode end

"""
    Value{T} <: AbstractNode

Abstract type for nodes contains a value in a computation graph.
"""
abstract type Value{T} <: AbstractNode end

"""
    ValueType{T}

Value types, it can be an [`AbstractVariable`](@ref) contains `T` or just a
value of type `T`.
"""
const ValueType{T} = Union{Value{T}, T}

"""
    AbstractVariable{T} <: Value{T}

Abstract type for variables, variables are types that contains value and gradients.
"""
abstract type AbstractVariable{T} <: Value{T} end

"""
    AbstractArrayVariable{T, N}

Alias for [`AbstractVariable`](@ref), abstract type for variables contains an array.
"""
const AbstractArrayVariable{T, N} = AbstractVariable{AT} where {T, N, AT <: AbstractArray{T, N}}

"""
    AbstractMatrixVariable{T}

Abstract type for variables contains a matrix. See [`AbstractVariable`](@ref) for more.
"""
const AbstractMatrixVariable{T} = AbstractArrayVariable{T, 2}

"""
    AbstractVectorVariable{T}

Abstract type for variables contains a vector. See [`AbstractVariable`](@ref) for more.
"""
const AbstractVectorVariable{T} = AbstractArrayVariable{T, 1}

"""
    Variable{T} <: Value{T}

A kind of leaf node. A general type for variables in a comput-graph.
Similar to PyTorch's Variable, gradient will be accumulated to `var.grad`.
"""
mutable struct Variable{T} <: AbstractVariable{T}
    value::T
    grad::T

    Variable(val::T) where T = new{T}(val, zero(val))
    Variable(val::T, grad::T) where T = new{T}(val, grad)
end

"""
    Node{FT, ArgsT} <: AbstractNode

General node in a comput-graph. It stores a callable operator `f` of type `FT`
and its arguments `args` in type `ArgsT` which should be a tuple.
"""
struct Node{FT <: Operator, ArgsT <: Tuple, KwargsT <: NamedTuple} <: AbstractNode
    f::FT
    args::ArgsT
    kwargs::KwargsT
end

# wrap function to Method
Node(f::Function, args, kwargs) = Node(Trait.Method(f), args, kwargs)
Node(op, args) = Node(op, args, NamedTuple())

"""
    CachedNode{NT, OutT} <: AbstractNode

Stores the cache of output with type `OutT` from a node of
type `NT` in comput-graph. CachedNode is mutable, its output
can be updated by [`forward`](@ref).
"""
struct CachedNode{NT <: AbstractNode, OutT} <: Value{OutT}
    node::NT
    output::OutT
end

function CachedNode(f, args...; kwargs...)
    node = Node(f, args, kwargs.data)
    output = forward(node)
    CachedNode(node, output)
end

Base.size(x::AbstractNode) = size(value(x))
Base.size(x::AbstractNode, d::Int) = size(value(x), d)
Base.similar(x::AbstractNode) = Variable(similar(value(x)))
Base.similar(x::AbstractNode, dims::Dims) = Variable(similar(value(x), dims))
Base.similar(x::AbstractNode, element_type::Type{S}, dims::Dims) where S = Variable(similar(value(x), element_type, dims))
Base.axes(x::AbstractNode) = axes(value(x))

"""
    arg(node, i) -> ArgumentType

Returns the `i`-th argument of the call in `node`.
"""
function arg end

"""
    args(node) -> Tuple

Returns the arguments of the call in `node`.
"""
function args end

"""
    kwargs(node) -> NamedTuple

Returns the keyword arguements of the call in `node`.
"""
function kwargs end

"""
    operator(node) -> YAAD.Operator

Returns the operator called in this node.
"""
function operator end

arg(x::Node, i::Int) = x.args[i]
args(x::Node) = x.args
kwargs(x::Node) = x.kwargs
operator(x::Node) = x.f

arg(x::CachedNode, i::Int) = x.node.args[i]
args(x::CachedNode) = x.node.args
kwargs(x::CachedNode) = x.node.kwargs
operator(x::CachedNode) = x.node.f

Base.eltype(x::AbstractNode) = eltype(value(x))

"""
    value(node)

Returns the value when forwarding at current node. `value` is different
than [`forward`](@ref) method, `value` only returns what the node contains,
it will throw an error, if this node does not contain anything.
"""
function value end

# forward other values
value(x) = x

function value(x::AbstractNode)
    error(
        "Expected value for node type $(typeof(x)) ",
        "check if you defined a non-cached node",
        " or overload value function for your node."
    )
end

value(x::Variable) = x.value
value(x::CachedNode) = value(x.output)

"""
    forward(node) -> output

Forward evaluation of the comput-graph. This method will call the operator
in the comput-graph and update the cache.

    forward(f, args...) -> output

For function calls.
"""
function forward end

"""
    forward!(output, ...) -> output
"""
function forward! end

forward(x) = x
forward(x::NT) where {NT <: AbstractNode} = error("forward method is not implemented for node type: $NT")
forward(x::Colon) = x
forward(node::Value) = value(node)
forward(node::Node) = forward(node.f, map(forward, node.args)...; map(forward, node.kwargs)...)
forward(op::Operator, args...; kwargs...) = op.f(args...; kwargs...)
forward(op::Trait.Broadcasted, args...) = Broadcast.broadcasted(op.f, args...)

# # This allow the evaluation of a static graph
# forward(node::CachedNode) = forward!(node, node.node) # forward to modify this
# forward!(v::CachedNode, node::Node) = forward!(v, node.f, map(forward, node.args)...; map(forward, node.kwargs)...)
# forward!(v::CachedNode, op::Operator, args...; kwargs...) = (v.output = op.f(args...; kwargs...))
# forward!(v::CachedNode, op::Trait.Broadcasted, args...) = (v.output = Broadcast.broadcasted(op.f, args...))

"""
    backward(node) -> nothing

Backward evaluation of the comput-graph.
"""
function backward end

# return nothing for non-node types
backward(x, grad) = x
backward(x::AbstractNode) = backward(x::AbstractNode, one(eltype(x)))

function backward(x::Variable, grad)
    @. x.grad += grad
    nothing
end

backward(node::CachedNode, grad) = backward(node, node.node.f, grad)
backward(node::CachedNode, op::Operator, grad) = backward(node, op.f, grad)

@generated function rm_constant_args(args::Tuple)
    ex = Expr[]
    for (i, each) in enumerate(args.parameters)
        push!(ex, :(tuple($i, args[$i])))
    end

    quote
        tuple($(ex...))
    end
end

function backward(node::CachedNode, f, grad)
    backward_type_assert(node, grad)
    # TODO: replace with @assert when there is a compiler option for it
    @boundscheck backward_size_assert(node, grad)

    grad_inputs = gradient(node, grad)

    for (each, each_grad) in zip(args(node), grad_inputs)
        backward(each, each_grad)
    end
    nothing
end

"""
    backward_type_assert(node, grad)

throw more readable error msg for backward type check.
"""
function backward_type_assert end

# mute the compiler error msg
backward_type_assert(args...) = true

backward_type_assert(node::CachedNode{<:AbstractNode, T}, grad::T) where T = true
# exclude arrays
backward_type_assert(node::CachedNode{<:AbstractNode, T1}, grad::T2) where
    {T, N, T1 <: AbstractArray{T, N}, T2 <: AbstractArray{T, N}} = true
backward_type_assert(node::CachedNode{<:AbstractNode, T1}, grad::T2) where {T1, T2} =
    error("Gradient is expected to have the same",
          " type with outputs, expected $T1",
          " got $T2")

function backward_size_assert(node::CachedNode, grad)
    size(node.output) == size(grad) ||
        error(
            "gradient should have the same size with output,",
            " expect size $(size(node.output)), got $(size(grad))"
        )
end

backward_size_assert(node::CachedNode{<:AbstractNode, T1}, grad::T2) where {T1 <: UniformScaling, T2} = true
backward_size_assert(node::CachedNode{<:AbstractNode, T1}, grad::T2) where {T1, T2 <: UniformScaling} = true

"""
    gradient(node, grad)

Returns the gradient.
"""
function gradient end

## CachedNode
# 1. general interface
gradient(x::CachedNode, grad) = gradient(x.node.f, grad, x.output, map(value, x.node.args)...; map(value, x.node.kwargs)...)

# NOTE: operators help to define different grads when the fn is the same
# e.g Broadcasted{typeof(sin)} and `sin`

# 2. forward operator to function type
# this simplifies some operator's definition
# unwrap the value here for convenience, but remember to define
# for different constant sometimes, e.g *(::Value, ::ConstantType) only need to calculate first
gradient(x::Operator, grad, output, args...; kwargs...) =
    gradient(x.f, grad, output, args...; kwargs...)

gradient(fn, grad, output, args...; kwargs...) =
    error(
        "gradient of operator $fn is not defined\n",
        "Possible Fix:\n",
        "define one of the following:\n",
        "1. gradient(::typeof($fn), grad, output, args...; kwargs...)\n",
        "2. gradient(op::Trait.Method{typeof($fn)}, grad, output, args...; kwargs...)\n",
        "3. gradient(op::Trait.Broadcasted{typeof($fn)}, grad, output, args...; kwargs...)\n"
    )

"""
    register(f, args...; kwargs...)

This is just a alias for constructing a `CachedNode`. But notice this function
is used for register a node in `tape` in the global tape version implementation:

https://github.com/Roger-luo/YAAD.jl/tree/tape
"""
register(f, args...; kwargs...) = CachedNode(f, args...; kwargs...)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using JSON
using Libdl
using LinearAlgebra

################################################################################
# First try to detect and load existing libmxnet
################################################################################
libmxnet_detected = false
libmxnet_curr_ver = get(ENV, "MXNET_COMMIT", "master")
curr_win = "20190608"  # v1.5.0

if haskey(ENV, "MXNET_HOME")
  MXNET_HOME = ENV["MXNET_HOME"]
  @info("MXNET_HOME environment detected: $MXNET_HOME")
  @info("Trying to load existing libmxnet...")
  # In case of macOS, if user build libmxnet from source and set the MXNET_HOME,
  # the output is still named as `libmxnet.so`.
  lib = Libdl.find_library(["libmxnet.$(Libdl.dlext)", "libmxnet.so"],
                           [joinpath(MXNET_HOME, "lib"), MXNET_HOME])
  if !isempty(lib)
    @info("Existing libmxnet detected at $lib, skip building...")
    libmxnet_detected = true
  else
    @info("Failed to load existing libmxnet, trying to build from source...")
  end
end

# Try to find cuda
CUDAPATHS = String[]
if haskey(ENV, "CUDA_HOME")
  push!(CUDAPATHS, joinpath(ENV["CUDA_HOME"], "lib64"))
elseif Sys.islinux()
  append!(CUDAPATHS, ["/opt/cuda/lib64", "/usr/local/cuda/lib64"])
end

if Sys.isunix()
  nvcc_path = Sys.which("nvcc")
  if nvcc_path ≢ nothing
    @info "Found nvcc: $nvcc_path"
    push!(CUDAPATHS, replace(nvcc_path, "bin/nvcc", "lib64"))
  end
end

HAS_CUDA = false
HAS_CUDNN = false
let cudalib = Libdl.find_library(["libcuda", "nvcuda.dll"], CUDAPATHS)
  global HAS_CUDA = !isempty(cudalib) && Libdl.dlopen_e(cudalib) != C_NULL
end

if !HAS_CUDA && Sys.iswindows()
  # TODO: this needs to be improved.
  try
    run(`nvcc --version`)
    global HAS_CUDA = true
  catch
  end
end

if HAS_CUDA  # then check cudnn
  let cudnnlib = Libdl.find_library("libcudnn", CUDAPATHS)
    global HAS_CUDNN = !isempty(cudnnlib) && Libdl.dlopen_e(cudnnlib) != C_NULL
    if HAS_CUDNN && !haskey(ENV, "CUDA_HOME")  # inference `CUDA_HOME`
      ENV["CUDA_HOME"] = dirname(dirname(cudnnlib))
    end
  end
end

if HAS_CUDA
  @info("Found a CUDA installation.")
  if HAS_CUDNN
    @info("Found a CuDNN installation.")
  end
  @info("CUDA_HOME -> $(get(ENV, "CUDA_HOME", "nothing"))")
else
  @info("Did not find a CUDA installation, using CPU-only version of MXNet.")
end

# propagate more build flags from ENV
const CC  = get(ENV, "CC", nothing)
const CXX = get(ENV, "CXX", nothing)
const ADD_CFLAGS  = get(ENV, "ADD_CFLAGS", nothing)
const ADD_LDFLAGS = get(ENV, "ADD_LDFLAGS", nothing)
const USE_JEMALLOC = get(ENV, "USE_JEMALLOC", nothing)  # "0" or "1"

function get_cpucore()
    if haskey(ENV, "TRAVIS")  # on travis-ci
        2
    else
        min(Sys.CPU_THREADS, 32)
    end
end

using BinDeps
@BinDeps.setup
if !libmxnet_detected
  if Sys.iswindows()
    if Sys.ARCH != :x86_64
      @info("Prebuilt windows binaries are only available on 64bit. You will have to built MXNet yourself.")
      return
    end
    @info("Downloading pre-built packages for Windows.")
    base_url = "https://github.com/yajiedesign/mxnet/releases/download/weekly_binary_build_v2/prebuildbase_win10_x64_vc14_v2.7z"

    if libmxnet_curr_ver == "master"
      _cmd = "{
        [System.Net.ServicePointManager]::SecurityProtocol='tls12';
        Invoke-WebRequest -Uri 'https://api.github.com/repos/yajiedesign/mxnet/releases/latest'
        -OutFile 'mxnet.json'}"
      # download_cmd uses powershell 2, but we need powershell 3 to do this
      run(`powershell -NoProfile -Command $_cmd`)
      curr_win = JSON.parsefile("mxnet.json")["tag_name"]
      @info("Can't use MXNet master on Windows, using latest binaries from $curr_win.")
    end
    # TODO: Get url from JSON.
    # TODO: detect cuda version and select corresponding url.
    name = "mxnet_x64_$(HAS_CUDA ? "vc141_gpu_cu101" : "vc14_cpu").7z"
    package_url = "https://github.com/yajiedesign/mxnet/releases/download/$(curr_win)/$(curr_win)_$(name)"

    exe7z = joinpath(Sys.BINDIR, "7z.exe")

    run(download_cmd(package_url, "mxnet.7z"))
    # this command will create the dir "usr\\lib"
    run(`$exe7z e mxnet.7z *\\build\\* *\\lib\\* -y -ousr\\lib`)

    run(download_cmd(base_url, "mxnet_base.7z"))
    run(`$exe7z x mxnet_base.7z -y -ousr`)
    run(`cmd /c copy "usr\\prebuildbase_win10_x64_vc14_v2\\3rdparty\\bin\\*.dll" "usr\\lib"`)

    # testing
    run(`cmd /c dir "usr\\lib"`)
    return
  end

  ################################################################################
  # If not found, try to build automatically using BinDeps
  ################################################################################

  blas_path = Libdl.dlpath(Libdl.dlopen(Base.libblas_name))
  blas_vendor = LinearAlgebra.BLAS.vendor()

  ilp64 = ""
  if blas_vendor == :openblas64
    ilp64 = "-DINTERFACE64"
  end

  FORCE_LAPACK = false
  if blas_vendor == :unknown
    @info("Julia is built with an unkown blas library ($blas_path).")
    @info("Attempting build without reusing the blas library")
    USE_JULIA_BLAS = false
  elseif !(blas_vendor in (:openblas, :openblas64))
    @info("Unsure if we can build against $blas_vendor.")
    @info("Attempting build anyway.")
    USE_JULIA_BLAS = true
  else
    USE_JULIA_BLAS = true
    FORCE_LAPACK = true
  end
  @info("USE_JULIA_BLAS -> $USE_JULIA_BLAS")

  blas_name = blas_vendor == :openblas64 ? "openblas" : string(blas_vendor)
  MSHADOW_LDFLAGS = "MSHADOW_LDFLAGS=-lm $blas_path"

  #--------------------------------------------------------------------------------
  # Build libmxnet
  mxnet = library_dependency("mxnet", aliases=["mxnet", "libmxnet", "libmxnet.so"])

  _prefix = joinpath(BinDeps.depsdir(mxnet), "usr")
  _srcdir = joinpath(BinDeps.depsdir(mxnet), "src")
  _mxdir  = joinpath(_srcdir, "mxnet")
  _libdir = joinpath(_prefix, "lib")
  # We have do eagerly delete the installed libmxnet.so
  # Otherwise we won't rebuild on an update.
  run(`rm -f $_libdir/libmxnet.$(Libdl.dlext)`)
  provides(BuildProcess,
    (@build_steps begin
      CreateDirectory(_srcdir)
      CreateDirectory(_libdir)
      @build_steps begin
        BinDeps.DirectoryRule(_mxdir, @build_steps begin
          ChangeDirectory(_srcdir)
          `git clone https://github.com/apache/incubator-mxnet mxnet`
        end)
        @build_steps begin
          ChangeDirectory(_mxdir)
          `git fetch`
          if libmxnet_curr_ver != "master"
            `git checkout $libmxnet_curr_ver`
          else
            `git checkout origin/$libmxnet_curr_ver`
          end
          `git submodule update --init --recursive`
          `git -C 3rdparty/mshadow checkout -- make/mshadow.mk`
          `cp -v ../../cblas.h include/cblas.h`
          `sed -i -s "s/MSHADOW_CFLAGS = \(.*\)/MSHADOW_CFLAGS = \1 $ilp64/" 3rdparty/mshadow/make/mshadow.mk`

          # Copy config.mk, always override the file
          if Sys.isapple()
            `cp make/osx.mk config.mk`
          else
            `cp make/config.mk config.mk`
          end

          # Configure OpenCV
          `sed -i -s 's/USE_OPENCV = 1/USE_OPENCV = 0/' config.mk`

          # Configure CUDA
          if HAS_CUDA
            @build_steps begin
              `sed -i -s 's/USE_CUDA = 0/USE_CUDA = 1/' config.mk`
              # address https://github.com/apache/incubator-mxnet/pull/7856
              `sed -i -s "s/ADD_LDFLAGS =\(.*\)/ADD_LDFLAGS =\1 -lcublas -lcusolver -lcurand -lcudart/" config.mk`
              if haskey(ENV, "CUDA_HOME")
                `sed -i -s "s@USE_CUDA_PATH = NONE@USE_CUDA_PATH = $(ENV["CUDA_HOME"])@" config.mk`
              end
              if haskey(ENV, "CUDA_HOME")
                # address https://github.com/apache/incubator-mxnet/pull/7838
                flag = "-L$(ENV["CUDA_HOME"])/lib64 -L$(ENV["CUDA_HOME"])/lib"
                `sed -i -s "s@ADD_LDFLAGS =\(.*\)@ADD_LDFLAGS =\1 $flag@" config.mk`
              end
              if HAS_CUDNN
                `sed -i -s 's/USE_CUDNN = 0/USE_CUDNN = 1/' config.mk`
              end
            end
          end

          # Force enable LAPACK build
          # Julia's OpenBLAS has LAPACK functionality already
          if FORCE_LAPACK
            if Sys.isapple()
              MSHADOW_LDFLAGS *= " -framework Accelerate"
            end
            `sed -i -s 's/ADD_CFLAGS =\(.*\)/ADD_CFLAGS =\1 -DMXNET_USE_LAPACK/' config.mk`
          end

          # propagate more build flags from ENV
          if CC != nothing
            `sed -i -s "s@^export CC =\(.*\)@export CC = $CC@" config.mk`
          end
          if CXX != nothing
            `sed -i -s "s@^export CXX =\(.*\)@export CXX = $CXX@" config.mk`
          end
          if ADD_CFLAGS != nothing
            `sed -i -s "s@ADD_CFLAGS =\(.*\)@ADD_CFLAGS =\1 $ADD_CFLAGS@" config.mk`
          end
          if ADD_LDFLAGS != nothing
            `sed -i -s "s@ADD_LDFLAGS =\(.*\)@ADD_LDFLAGS =\1 $ADD_LDFLAGS@" config.mk`
          end
          if USE_JEMALLOC != nothing
            `sed -i -s "s@USE_JEMALLOC =\(.*\)@USE_JEMALLOC = $USE_JEMALLOC@" config.mk`
          end

          if USE_JULIA_BLAS
            `make -j$(get_cpucore()) USE_BLAS=$blas_name $MSHADOW_LDFLAGS`
          else
            `make -j$(get_cpucore())`
          end
        end
        FileRule(joinpath(_libdir, "libmxnet.$(Libdl.dlext)"), @build_steps begin
          # the output file on macos is still in `.so` suffix
          # so we rename it
          `cp $_mxdir/lib/libmxnet.so $_libdir/libmxnet.$(Libdl.dlext)`
        end)
      end
    end), mxnet, installed_libpath=_libdir)

  @BinDeps.install Dict(:mxnet => :mxnet)
end

max_stars_repo_name: SHEwiner/incubator-mxnet

max_issues_repo_name: chongruo/incubator-mxnet

max_forks_repo_name: ElaineBao/incubator-mxnet


text:
#!/usr/bin/env julia

# Read the data file and plot the location of all of the baselines, in k\lambda
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.gridding

import PyPlot.plt
using LaTeXStrings

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species * transition]
model = config["model"]

pars = convert_dict(config["parameters"], config["model"])

# Load the data file so we can use it to compute lnprob
dvarr = DataVis(config["data_file"])
dv = dvarr[1]

fig, ax = plt[:subplots](figsize = (6, 6))
ax[:plot](dv.uu, dv.vv, ".")
ax[:set_xlabel](L"UU [k$\lambda$]")
ax[:set_ylabel](L"VV [k$\lambda$]")

plt[:savefig]("baselines.png")


text:
using Pkg
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("EXPRES-StellarSignals")   end
cd("submissions/2021feb15")
Pkg.activate(".")
starid = 101501
include("line_by_line_fits_generic.jl")


text:
function system_lhs(boat_angles::Boat_angles)
    inertia_lhs(SA[boat_angles.θ_t, boat_angles.θ_k, boat_angles.γ_OA])
end

function system_rhs(u, w, θ′, boat_angles::Boat_angles)
    inertia_rhs(
        SA[
            u,
            w,
            θ′,
            boat_angles.θ_t,
            boat_angles.θ′_t,
            boat_angles.θ′′_t,
            boat_angles.θ_k,
            boat_angles.θ′_k,
            boat_angles.θ′′_k,
            boat_angles.γ_OA,
            boat_angles.γ′_OA,
            boat_angles.γ′′_OA,
        ],
    )
end


text:
# This example shows how to use MPI.Scatterv! and MPI.Gatherv!
# roughly based on the example from
# https://stackoverflow.com/a/36082684/392585

using MPI

"""
    split_count(N::Integer, n::Integer)

Return a vector of `n` integers which are approximately equally sized and sum to `N`.
"""
function split_count(N::Integer, n::Integer)
    q,r = divrem(N, n)
    return [i <= r ? q+1 : q for i = 1:n]
end

MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

root = 0

if rank == root
    M, N = 4, 7

    test = Float64[i for i = 1:M, j = 1:N]
    output = similar(test)

    # Julia arrays are stored in column-major order, so we need to split along the last dimension
    # dimension
    M_counts = [M for i = 1:comm_size]
    N_counts = split_count(N, comm_size)

    # store sizes in 2 * comm_size Array
    sizes = vcat(M_counts', N_counts')
    size_ubuf = UBuffer(sizes, 2)

    # store number of values to send to each rank in comm_size length Vector
    counts = vec(prod(sizes, dims=1))

    test_vbuf = VBuffer(test, counts) # VBuffer for scatter
    output_vbuf = VBuffer(output, counts) # VBuffer for gather
else
    # these variables can be set to `nothing` on non-root processes
    size_ubuf = UBuffer(nothing)
    output_vbuf = test_vbuf = VBuffer(nothing)
end

if rank == root
    println("Original matrix")
    println("================")
    @show test sizes counts
    println()
    println("Each rank")
    println("================")
end
MPI.Barrier(comm)

local_size = MPI.Scatter(size_ubuf, NTuple{2,Int}, root, comm)
local_test = MPI.Scatterv!(test_vbuf, zeros(Float64, local_size), root, comm)

for i = 0:comm_size-1
    if rank == i
        @show rank local_test
    end
    MPI.Barrier(comm)
end

MPI.Gatherv!(local_test, output_vbuf, root, comm)

if rank == root
    println()
    println("Final matrix")
    println("================")
    @show output
end


text:
module ScalingTests

using Interpolations
using Base.Test

# Model linear interpolation of y = -3 + .5x by interpolating y=x
# and then scaling to the new x range

itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid())

sitp = @inferred(scale(itp, -3:.5:1.5))

for (x,y) in zip(-3:.05:1.5, 1:.1:10)
    @test_approx_eq sitp[x] y
end

# Verify that it works in >1D, with different types of ranges

gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)

xs = -5:.5:5
ys = -4:.2:4
zs = Float64[testfunction(x,y) for x in xs, y in ys]

itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid())
sitp2 = @inferred scale(itp2, xs, ys)

for x in xs, y in ys
    @test_approx_eq testfunction(x,y) sitp2[x,y]
end

# Test gradients of scaled grids
xs = -pi:.1:pi
ys = sin(xs)
itp = interpolate(ys, BSpline(Linear()), OnGrid())
sitp = @inferred scale(itp, xs)

for x in -pi:.1:pi
    g = @inferred(gradient(sitp, x))[1]
    @test_approx_eq_eps cos(x) g .05
end

# Verify that return types are reasonable
@inferred(getindex(sitp2, -3.4, 1.2))
@inferred(getindex(sitp2, -3, 1))
@inferred(getindex(sitp2, -3.4, 1))

sitp32 = @inferred scale(interpolate(Float32[testfunction(x,y) for x in -5:.5:5, y in -4:.2:4], BSpline(Quadratic(Flat())), OnGrid()), -5f0:.5f0:5f0, -4f0:.2f0:4f0)
@test typeof(@inferred(getindex(sitp32, -3.4f0, 1.2f0))) == Float32

# Iteration
itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat())), OnCell())
knots = map(d->1:100:201, 1:3)
sitp = @inferred scale(itp, knots...)
function foo!(dest, sitp)
    i = 0
    for s in eachvalue(sitp)
        dest[i+=1] = s
    end
    dest
end
function bar!(dest, sitp)
    for I in CartesianRange(size(dest))
        dest[I] = sitp[I]
    end
    dest
end
rfoo = Array(Float64, Interpolations.ssize(sitp))
rbar = similar(rfoo)
foo!(rfoo, sitp)
bar!(rbar, sitp)
@test_approx_eq rfoo rbar
@test (@elapsed foo!(rfoo, sitp)) < (@elapsed bar!(rbar, sitp))

end


text:
module Gear

using JuMP, BARON, Test

m = Model(BARON.Optimizer)
@variable(m, 12 ≤ x[i=1:4] ≤ 60, start = 24)

@NLconstraints(m, begin
    x[3] ≤ x[4]
    x[2] ≤ x[1]
end)

@NLobjective(m, Min, 6.931 - x[1]*x[2]/(x[3]*x[4])^2 + 1)

optimize!(m)

@test isapprox(value(x[1]), 60, rtol=1e-6)
@test isapprox(value(x[2]), 60, rtol=1e-6)
@test isapprox(value(x[3]), 12, rtol=1e-6)
@test isapprox(value(x[4]), 12, rtol=1e-6)
@test isapprox(objective_value(m), 7.75738888889, rtol=1e-6)

end # module


text:
module AlignmentAlgorithms

export nwscore,
       nwsetupmatrix,
       nwscorematrix,
       swsetupmatrix,
       swscorematrix
#export nwalign

include("needleman_wunsch.jl")
include("smith_waterman.jl")

end


text:
struct GenVariant
    chrom::String
    varid::String
    rsid::String
    pos::UInt
    alleles::Vector{String}
    probs::Matrix{Float64}
    ploidy::Vector{UInt8}
end

"""
    load_gen_data()
load data from "example.gen"
"""
function load_gen_data()
    variants = GenVariant[]
    open(BGEN.datadir("example.gen")) do io
        for l in readlines(io)
            tokens = split(l)
            chrom, varid, rsid, pos, ref, alt = tokens[1:6]
            pos = parse(UInt, pos)
            probs = reshape(parse.(Float64, tokens[7:end]), 3, :)
            nan_col_idx = transpose(sum(probs; dims=1)) .== 0
            for (i, v) in enumerate(nan_col_idx)
                if v
                    probs[:, i] .= NaN
                end
            end
            push!(variants, GenVariant(chrom, varid, rsid, pos, [ref, alt], probs, []))
        end
    end
    variants
end

"""
    parse_vcf_samples(format, samples)
Parse sample data from VCF.
"""
function parse_vcf_samples(format, samples)
    samples = [Dict(zip(split(format, ":"), split(x, ":"))) for x in samples]
    ks = Dict([("GT", r"[/|]"), ("GP", ","), ("HP", ",")])
    samples2 = [Dict{String, Vector{SubString{String}}}() for x in samples]
    for (x, y) in zip(samples, samples2)
        for k in keys(x)
            y[k] = split(x[k], ks[k])
        end
    end
    probs = [occursin("GP", format) ? x["GP"] : x["HP"] for x in samples2]
    probs = map(x -> parse.(Float64, x), probs)
    max_len = maximum(length(x) for x in probs)
    probs_out = Matrix{Float64}(undef, max_len, length(samples))
    fill!(probs_out, NaN)
    for i in 1:length(probs)
        probs_out[1:length(probs[i]), i] .= probs[i]
    end
    ploidy = [length(y["GT"]) for y in samples2]
    probs_out, ploidy
end

"""
    load_vcf_data()
Load data from "complex.vcf" for comparison
"""
function load_vcf_data()
    variants = GenVariant[]
    open(BGEN.datadir("complex.vcf")) do io
        for l in readlines(io)
            if startswith(l, "#")
                continue
            end
            tokens = split(l)
            chrom, pos, varid, ref, alts = tokens[1:5]
            pos = parse(UInt, pos)
            format = tokens[9]
            samples = tokens[10:end]
            varid = split(varid, ",")
            if length(varid) > 1
                rsid, varid = varid
            else
                rsid, varid = varid[1], ""
            end
            probs, ploidy = parse_vcf_samples(format, samples)
            var = GenVariant(string(chrom), varid, string(rsid), pos, vcat(String[ref],
                                                split(alts, ",")), probs, ploidy)
            push!(variants, var)
        end
    end
    variants
end

"""
    load_haps_data()
Load data from "haplotypes.haps" for comparison.
"""
function load_haps_data()
    variants = GenVariant[]
    open(BGEN.datadir("haplotypes.haps")) do io
        for l in readlines(io)
            tokens = split(l)
            chrom, varid, rsid, pos, ref, alt = tokens[1:6]
            pos = parse(UInt, pos)
            probs = tokens[7:end]
            probs = [x == "0" ? [1.0, 0.0] : [0.0, 1.0] for x in probs]
            probs = [probs[pos:pos+1] for pos in 1:2:length(probs)]
            probs = [hcat(x[1], x[2]) for x in probs]
            probs = hcat(probs...)
            probs = reshape(probs, 4, :)
            var = GenVariant(chrom, varid, rsid, pos, [ref, alt], probs, [])
            push!(variants, var)
        end
    end
    variants
end

"""
    epsilon(bit_depth)
Max difference expected for the bit depth
"""
@inline function epsilon(bit_depth)
    return 2 / (2 ^ (bit_depth - 1))
end

"""
    array_equal(truth, parsed, bit_depth)
check if the two arrays are sufficiently equal
"""
@inline function array_equal(truth, parsed, bit_depth)
    eps_abs = 3.2e-5
    all(isapprox.(truth, parsed; atol=max(eps_abs, epsilon(bit_depth)),
        nans=true))
end


text:
function newtonav!(yout::Vector{Float64},iters::Vector{Int64},Vs::Float64,vs::Float64,
    zs::Float64,ts::Float64,tolz::Float64,state::String,Kvo::Float64,Kvc::Float64,
    Plv::Float64,Pa::Float64,Q::Float64,Vlv::Float64,zeta::Float64,W1root::Float64,W2root::Float64,
    E::Float64,V0::Float64,beta::Float64,rho::Float64,c0::Float64,A0::Float64,Aann::Float64,
    leff::Float64,Ks::Float64,th::Float64,h::Float64,f::Function,J::Function,maxiter::Int16,
    maxval::Float64,epsJ::Float64,epsN::Float64)

    # initial guess
    if state == "opening"
        zdot = (1.-zeta).*Kvo.*(Plv.-Pa.*mmHgToPa);
    else
        zdot = zeta.*Kvc.*(Plv.-Pa.*mmHgToPa);
    end
    if abs(zdot) <= tolz
        println("dζ/dt under tolerance at HR = $(60/th) bpm. tol = $tolz.
            ζ = $zeta. Switching to reduced system.")
        x0 = zeros(2);
        xx = zeros(2);
        if Q == 0
            x0[1] = -W2root + 1e-7;
        else
            x0[1] = W1root;
        end
        if Q != 0
            x0[2] = Vlv - h*Q;
        elseif state == "opening" && Q == 0
            x0[2] = Vlv;
        end

        # non-dimensionalize
        x0[1] /= vs;
        x0[2] /= Vs;
        # println(x0)
    else
        x0 = zeros(3);
        xx = zeros(3);
        if Q == 0
            x0[1] = -W2root + 1e-7;
        else
            x0[1] = W1root;
        end
        if Q != 0
            x0[2] = Vlv - h*Q;
        elseif state == "opening" && Q == 0
            x0[2] = Vlv;
        end
        if zeta == 0
            x0[3] = zs;
        elseif state == "opening"
            x0[3] = zeta + h*zdot;
        elseif state == "closing" && abs(Plv/mmHgToPa - Pa) < 20
            x0[3] = zeta + h*zdot;
        elseif state == "closing"
            x0[3] = zeta;
        end

        # non-dimensionalize
        x0[1] /= vs;
        x0[2] /= Vs;
        x0[3] /= zs;
        # println(x0)
    end

    # assign function/Jacobian handles
    # f = CVModule.fav;
    # J = CVModule.Jav;

    # setup for iterations
    xx .= x0;
    x = zeros(length(x0));
    N = 1;
    Plv = 0.;
    Pao = 0.;

    # max step size for line searches
    stpmax = 100*max(sqrt(norm(xx)),length(x0));

    # println(state)
    # Newton iterations
    while N <= maxiter
        Plv = E.*(xx[2].*Vs .- V0);
        Pao = beta.*((2.*rho./beta).*((xx[1].*vs .- W2root)./8 .+ c0).^2 .- sqrt.(A0));
        if Plv > Pao
            state = "opening";
        else
            state = "closing";
            # println(state)
        end
        # println(xx)
        # determine Jacobian, check invertibility
        JJ = J(xx,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        D = diagm(maximum!(zeros(length(xx)),abs.(JJ)).^-1);
        # println(JJ)
        # println(D)
        # println(D*JJ)
        # println(cond(D*JJ))
        # println(det(D*JJ))
        if abs(det(D*JJ)) < epsJ
            println(JJ)
            println(D*JJ)
            print(det(D*JJ))
            error("Newton Jacobian is singular.");
        end
        # compute gradient of line search objective function
        fvec = D*f(xx,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        # println(f(xx,system,n,state))
        # println(fvec)
        fold = 0.5*dot(fvec,fvec);
        g = (D*JJ)'*fvec;
        # compute newton step
        # println(inv(D*JJ))
        s = -inv(D*JJ)*fvec;
        # line search to update state vector
        fn, xn, check = CVModule.linesearch(xx,fold,g,s,stpmax,state,
            Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,f,J,maxiter);
        # println(xn[1]*vs)
        # check if sufficiently close to root
        JJ = J(xn,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        D = diagm(maximum!(zeros(length(xx)),abs.(JJ)).^-1);
        fvec = D*f(xn,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        # println(fvec)
        if norm(fvec) <= epsN*1000
            x = xn;
            x[1] = x[1]*vs;
            x[2] = x[2]*Vs;
            if length(x) == 3
                x[3] = x[3]*zs;
            end
            # println(x)
            # println(system.branches.W2root)
            # println(inv(D*JJ))
            # println(fvec)
            iters[1] += N;
            break
        end
        if check
            test = 0.
            den = max(fn,0.5*length(xn));
            for i = 1:length(xn)
                temp = abs(g[i])*max(abs(xn[i]),1.)/den;
                if temp > test
                    test = temp;
                end
            end
            if test < 1e-6
                check = true;
                println("Warning: gradient of objective function close to zero.")
            else
                check = false;
            end
        end
        # check for divergence
        if norm(fvec,Inf) >= maxval
            iters[1] += N;
            println(xn)
            println(fvec)
            error("Newton iteration diverged.");
        end
        N+=1;
        xx = xn;
        if N == maxiter
            println(JJ)
            # println(D)
            # println(D*JJ)
            println(xn)
            println(fvec)
            println(norm(fvec))
            error("Newton iteration failed to converge.");
        end
    end

    # update based on converged solution
    yout[1] = ((2.*rho./beta).^2.*((x[1] .- W2root)./8 .+ c0).^4);
    yout[2] = 0.5.*yout[1].*(x[1] .+ W2root);
    yout[3] = x[2];
    if length(x) == 3
        if x[3] < 1e-5 && state == "closing"
            yout[4] = 0;
            println("Aortic valve closed.")
        else
            yout[4] = x[3];
        end
    else
        if zeta < 1e-5 && state == "closing"
            yout[4] = 0;
            println("Aortic valve closed.")
        else
            yout[4] = zeta;
        end
    end

    yout[5] = x[1];
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
"""
    Recycler(restartA, inB)

Set up the recycler.  This requires that the functions, restartA and inB, take,
as arguments, the same state type as used in the underlying sampling regime

### Fields
* `restartA!`         - Resamples the source state A
* `inB`               - Detects if in state B
* `n_recycle_iters`   - Number of iterations between recycling check
"""
struct Recycler{TA, TB} <: AbstractRecycler
    restartA!::TA
    inB::TB
    n_recycle_iters::Int
end

"""
    Recycler(restartA!::TA, inB::TB; n_recycle_iters = 1)

Construct a recycler.

### Fields
* `restartA!`         - Resamples the source state A
* `inB`               - Detects if in state B
### Optional Fields
* `n_recycle_iters = 1`   - Number of iterations between recycling check
"""
function Recycler(restartA!::TA, inB::TB; n_recycle_iters = 1) where {TA, TB}

        return Recycler(restartA!, inB, n_recycle_iters)
end

"""
    sample_trajectory!(x, sampler, recycler; options=MDOptions())

In place applciation of the `sampler` to `x`.  Number of iterations are
set using the `options` argument.

### Fields

* `x`         - Starting position for sampler, modified in place
* `recycler`  - Recycling structure for A→B transitions
* `sampler`   - Desired sampler
* `options`   - Sampling options, including number of iteration

"""
function sample_trajectory!(x::Tx, sampler::S, recycler::R; options=MDOptions()) where {Tx, S<:AbstractSampler, R<:AbstractRecycler}

    state = InitState!(x, sampler);
    for i in 1:options.n_iters
        if(mod(i-1,recycler.n_recycle_iters)==0)
            if (recycler.inB(state))
                recycler.restartA!(state)
            end
        end
        UpdateState!(state, sampler);
    end
    x
end

"""
    sample_trajectory(x₀, sampler, recycler; options=MDOptions())

Run the `sampler` starting at `x₀`.  Number of iterations and interval between
saves are set using the `options` argument.  For Metropolis samplers, the
running acceptance rates are also resturned.

### Fields

* `x`         - Starting position for sampler, modified in place
* `recycler`  - Recycling structure for A→B transitions
* `sampler`   - Desired sampler
* `options`   - Sampling options, including number of iteration

"""
function sample_trajectory(x₀::Tx, sampler::S, recycler::R; options=MDOptions()) where {Tx,  S<:MetropolisSampler, R<:AbstractRecycler}

    n_accept = Int(0);

    state = InitState(x₀, sampler);

    # allocate memory for samples
    samples = Tx[similar(x₀) for i = 1:options.n_save];
    acceptance_rates = zeros(options.n_save);
    save_index = 1;
    for i = 1:options.n_iters
        if(mod(i-1,recycler.n_recycle_iters)==0)
            if (recycler.inB(state))
                recycler.restartA!(state)
            end
        end
        UpdateState!(state, sampler);
        n_accept+=state.accept;
        if(mod(i,options.n_save_iters)==0)
            @. samples[save_index] = deepcopy(state.x);
            acceptance_rates[save_index] = n_accept/i;
            save_index+=1;
        end
    end
    return samples, acceptance_rates
end

function sample_trajectory(x₀::Tx, sampler::S, recycler::R; options=MDOptions()) where {Tx,  S<:NonMetropolisSampler, R<:AbstractRecycler}

    state = InitState(x₀, sampler);
    samples = Tx[similar(x₀) for i = 1:options.n_save];
    save_index = 1;
    for i = 1:options.n_iters
        if(mod(i-1,recycler.n_recycle_iters)==0)
            if (recycler.inB(state))
                recycler.restartA!(state)
            end
        end
        UpdateState!(state, sampler);
        if(mod(i,options.n_save_iters)==0)
            @. samples[save_index] = deepcopy(state.x);
            save_index+=1;
        end
    end
    return samples
end

max_stars_repo_name: gideonsimpson/JuBasicMD

max_issues_repo_name: gideonsimpson/JuBasicMD

max_forks_repo_name: gideonsimpson/JuBasicMD


text:
module MCMCChain

import Showoff: showoff
import StatsBase: autocor, autocov, countmap, counts, describe, predict,
       quantile, sample, sem, summarystats
#import Plots: plot
import LinearAlgebra: diag

using RecipesBase
using StatPlots
using Distributions
using SpecialFunctions

export Chains, getindex, setindex!
export plot, traceplot, meanplot, densityplot, histogramplot, mixeddensityplot, autcorplot
export describe

# export diagnostics functions
export discretediag, gelmandiag, gewekediag, heideldiag, rafterydiag

abstract type AbstractChains end

struct Chains <: AbstractChains
	value::Array{Float64, 3}
	range::AbstractRange{Int}
	names::Vector{AbstractString}
	chains::Vector{Int}
end

# imports
include("utils.jl")

include("chains.jl")
include("chainsummary.jl")
include("discretediag.jl")
include("fileio.jl")
include("gelmandiag.jl")
include("gewekediag.jl")
include("heideldiag.jl")
include("mcse.jl")
#include("modelchains.jl")
#include("modelstats.jl")
include("rafterydiag.jl")
include("stats.jl")
include("plot.jl")
#include("plot2.jl")

end # module


text:
include("branch.jl")
include("bus.jl")
include("generator.jl")
include("load.jl")
include("storage.jl")

"""Validates an array of devices."""
function validate_devices(devices::Array{<: Device, 1})::Bool
    is_valid = true

    for device in devices
        if !validate(device)
            is_valid = false
        end
    end

    return is_valid
end

"""Validates the contents of a System."""
function validate(sys::System)::Bool
    is_valid = true

    for field in (:buses, :loads)
        objs = getfield(sys, field)
        if !validate_devices(objs)
            is_valid = false
        end
    end

    if !validate(sys.generators)
        is_valid = false
    end

    for field in (:branches, :storage)
        objs = getfield(sys, field)
        if !isnothing(objs)
            if !validate_devices(objs)
                is_valid = false
            end
        end
    end

    @debug "System validation" is_valid
    return is_valid
end


text:
#=
CDLP:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using DataStructures

#minMode function: given a collection [1, 5, 5, 3, 3, 7] will return 3 because the highest repetition is 2
#and the lowest most repeated number is 3
function minMode(c)

    minVal = -1
    i = 0

    for row in eachrow(c)
    i = i+1
    row = collect(row)
    temp = filter(x -> x!=0, row)
    count = counter(temp)
    sortedCollection = sort(collect(count), by=x->x[2], rev=true)
    minVal = sortedCollection[1][1]
    repetitions = sortedCollection[1][2]
    for (key, value) in sortedCollection
        if(value == repetitions)
            if(key < minVal)
                minVal = key
            end
        end
    end

    end

    return minVal
end

#can be useful to tweak the algorithm
function diag_conversion(lab, n)
for i = 1:n
        for j = 1:n
            if(i==j)
                lab[i,j]=i
            else
                lab[i,j]=0
            end
        end
    end
    return lab
end

#A is the input matrix, n is the number of nodes, t the number of iterations Output: label propagation matrix
function cdlp(A, n, t)

    lab = GBMatrix{Int64}(n,n)

    for i = 1:n
        for j = 1:n
            if(i==j)
                lab[i,j]=i
            else
                lab[i,j]=0
            end
        end
    end

    print("STARTING\n")
    for k = 1:t
            F = mul(A, lab, Semirings.PLUS_TIMES)
            print("\n\n")
            for i = 1:n
            r = F[i,:]
            r = sort(r, dims=2)
            lab[i, i] = minMode(r)
        end

        print("\n\nITERATION COMPLETED\n\n")
        print(sparse(lab))

    end

    print("\n\n\nFINISHED")
    return lab

end


text:
# f(x)=1./(1+exp(-x))
Fsigmoid(x)=(sigmoid(x),[]);
Fsigmoid_inplace(value,aux,x)=copy!(value,sigmoid(x))

Dsigmoid(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)=axpy!(1.0,grad_c.*f_c.*(1.-f_c),grad_n)

if PROC=="GPU"
    #function Fsigmoid(x::CudaArray)
    #    tmp=copy(x)
    #    sigmoid!(x,tmp)
    #    return (tmp,nothing) ## problem here is that I don't know how to free tmp -- this will cause memory leak I think.
    #end
    Fsigmoid_inplace(value,aux,x::CudaArray)=sigmoid!(x,value)
    #Fsigmoid_inplace(value,aux,x::CudaArray{Float32})=sigmoid32!(x,value)

    Dsigmoid(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)=tx1mx!(grad_c,f_c,grad_n)

end

Derivative[Fsigmoid]=Dsigmoid
Inplace[Fsigmoid]=Fsigmoid_inplace

sigmoid(A::ADnode)=ADnode(Fsigmoid,A)
export sigmoid

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using RecipesBase

@recipe function f(tree::Phylo.AbstractTree; treetype = :dendrogram, marker_group = nothing, line_group = nothing, showtips = true, tipfont = (7,))

    linecolor --> :black
    grid --> false
    framestyle --> :none
    legend --> false
    colorbar --> true
    size --> (1000, 1000)

    d, h, n = _findxy(tree)
    adj = 0.03maximum(values(d))
    tipannotations = map(x->(d[x] + adj, h[x], x), getleafnames(tree))

    lz = get(plotattributes, :line_z, nothing)
    mz = get(plotattributes, :marker_z, nothing)
    isnothing(lz) || (line_z := _handlez(lz, tree, n))
    isnothing(mz) || (marker_z := _handlez(mz, tree, n))
    mg = _handlez(marker_group, tree, n)
    lg = _handlez(line_group, tree, n)

    x, y = Float64[], Float64[]
    for node ∈ n
        if hasinbound(tree, node)
            push!(x, d[getparent(tree, node)], d[getparent(tree, node)], d[node], NaN)
            push!(y, h[getparent(tree, node)], h[node], h[node], NaN)
        end
    end

    marker_x, marker_y = _handlemarkers(plotattributes, mg, tree, d, h, n)

    if treetype == :dendrogram
        Dendrogram(x, y, tipannotations, marker_x, marker_y, showtips, tipfont, mg, lg)
    elseif treetype == :fan
        Fan(x, y, tipannotations, marker_x, marker_y, showtips, tipfont, mg, lg)
    else
        throw(ArgumentError("Unsupported `treetype`; valid values are `:dendrogram` or `:fan`"))
    end
end

struct Dendrogram; x; y; tipannotations; marker_x; marker_y; showtips; tipfont; marker_group; line_group; end
struct Fan; x; y; tipannotations; marker_x; marker_y; showtips; tipfont; marker_group; line_group; end

@recipe function f(dend::Dendrogram)
    ex = extrema(filter(isfinite, dend.x))
    xlims --> (ex[1] - 0.05 * ex[2], ex[2] * 1.15)

    # tip annotations
    dend.showtips && (annotations := map(x -> (x[1], x[2], (x[3], :left, dend.tipfont...)), dend.tipannotations))

    sa = get(plotattributes, :series_annotations, nothing)
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing
        label --> ""
    #    primary := false

        lc = _extend(get(plotattributes, :linecolor, nothing), dend.x)
        lc !== nothing && (linecolor := lc)
        la = _extend(get(plotattributes, :linealpha, nothing), dend.x)
        la !== nothing && (linealpha := la)
        lz = _extend(get(plotattributes, :line_z, nothing), dend.x)
        lz !== nothing && (line_z := lz)

        dend.x, dend.y
    end
    if !isempty(dend.marker_x) || sa !== nothing
        if isnothing(dend.marker_group)
            @series begin
                seriestype := :scatter
                sa !== nothing && (series_annotations := sa)
                label --> ""
                dend.marker_x, dend.marker_y
            end
        else
            groups = sort(unique(dend.marker_group))
            for group in groups
                idxs = findall(==(group), dend.marker_group)
                @series begin
                    seriestype := :scatter
                    sa !== nothing && (series_annotations := sa[idxs])
                    label --> string(group)
                    dend.marker_x[idxs], dend.marker_y[idxs]
                end
            end
        end
    end
    primary := false
    label := ""
    nothing
end

@recipe function f(fan::Fan)
    adjust(y) = 2pi*y / (length(fan.tipannotations) + 1)

    aspect_ratio := 1

    # tip annotations
    mx = maximum(filter(isfinite, fan.x))
    if fan.showtips
        xlims --> (1.5 .* (-mx, mx))
        ylims --> (1.5 .* (-mx, mx))
        annotations := map(x -> (_tocirc(x[1], adjust(x[2]))..., (x[3], :left,
            rad2deg(adjust(x[2])), fan.tipfont...)), fan.tipannotations)
    end

    sa = get(plotattributes, :series_annotations, nothing)
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing
        label := ""

        x, y = _circle_transform_segments(fan.x, adjust(fan.y))
        lc = _extend(get(plotattributes, :linecolor, nothing), x)
        lc !== nothing && (linecolor := lc)
        la = _extend(get(plotattributes, :linealpha, nothing), x)
        la !== nothing && (linealpha := la)
        lz = _extend(get(plotattributes, :line_z, nothing), x)
        lz !== nothing && (line_z := lz)
        x, y
    end
    if !isempty(fan.marker_x) || sa !== nothing
        if isnothing(fan.marker_group)
            @series begin
                seriestype := :scatter
                sa !== nothing && (series_annotations := sa)
                label --> ""
                _xcirc.(adjust(fan.marker_y), fan.marker_x), _ycirc.(adjust(fan.marker_y), fan.marker_x)
            end
        else
            groups = sort(unique(fan.marker_group))
            for group in groups
                idxs = findall(==(group), fan.marker_group)
                @series begin
                    seriestype := :scatter
                    sa !== nothing && (series_annotations := sa[idxs])
                    label --> string(group)
                    _xcirc.(adjust(fan.marker_y[idxs]), fan.marker_x[idxs]), _ycirc.(adjust(fan.marker_y[idxs]), fan.marker_x[idxs])
                end
            end
        end
    end
    primary := false
    label := ""
    nothing
end

_handlez(x, tree, names) = x
_handlez(x::Union{String, Symbol}, tree, names) = [getnodedata(tree, name, x) for name in names]
function _handlez(x::Dict, tree, names) where T
    ret = fill(NaN, length(names))
    for (i, n) in enumerate(names)
        name = getnodename(tree, n)
        if haskey(x, name)
            ret[i] = x[name]
        end
    end
    ret
end

_mylength(x) = 1
_mylength(x::AbstractVector) = length(x)
function _handlemarkers(plotattributes, marker_group, tree, d, h, names)
    isnothing(marker_group) || (plotattributes[:marker_group] = marker_group)
    markerfields = filter(x->occursin(r"marker", String(x)), keys(plotattributes))
    isempty(markerfields) && return (Float64[], Float64[])
    maxlength =  maximum([_mylength(plotattributes[k]) for k in markerfields])
    if maxlength ∈ (1, ninternal(tree))
        names = filter(node -> !isleaf(tree, node), names)
    end
    [d[node] for node in names], [h[node] for node in names]
end

function _extend(tmp, x)
    tmp isa AbstractVector && abs(length(tmp) - count(isnan, x)) < 2 || return nothing
    ret = similar(x, eltype(tmp))
    j = 1 + length(tmp) - count(isnan, x)
    for i in eachindex(x)
        ret[i] = tmp[j]
        isnan(x[i]) && (j += 1)
    end
    return ret
end

"""
    sort!(::AbstractTree; rev = false)

Sorts the branches descending from each node by total number of
descendants. This creates a clearer tree for plotting. The
process is also called "ladderizing" the tree. Use `rev=true` to
reverse the sorting order.
"""
function Base.sort!(tree::AbstractTree; rev = false)
    function loc!(clade::String)
        if isleaf(tree, clade)
            return 1
        end

        sizes = map(loc!, getchildren(tree, clade))
        node = getnode(tree, clade)
        node.other .= node.other[sortperm(sizes, rev = rev)]
        sum(sizes) + 1
    end

    loc!(first(nodenamefilter(isroot, tree)))
    tree
end

"""
    sort(::AbstractTree; rev = false)

Copies a tree and sorts its branches. See `sort!` for further details.
"""
Base.sort(tree::AbstractTree; rev = false) = sort!(copy(tree))

function _nodedepths(tree::Phylo.AbstractTree)
    function finddepths!(clade::String)
        if !in(clade, keys(depth))
            for subclade in getchildren(tree, clade)
                finddepths!(subclade)
            end
        end
        if !isleaf(tree, clade)
            ch_depths = [depth[child] for child in getchildren(tree, clade)]
            depth[clade] = (maximum(ch_depths) + minimum(ch_depths)) / 2.
        end
    end

    depth = Dict(tip => float(i) for (i, tip) in enumerate(getnodename(tree, x) for x in traversal(tree, preorder) if isleaf(tree, x)))
    sizehint!(depth, nnodes(tree))
    root = getnodename(tree, getroot(tree))
    finddepths!(root)
    depth
end

function _findxy(tree::Phylo.AbstractTree)
    height = nodeheights(tree)
    depth = _nodedepths(tree)
    height, depth, AxisArrays.axes(height, 1).val
end

function _find_tips(height, depth, tree)
    x, y, l = Float64[], Float64[], String[]
    for k in keys(height)
        if isleaf(tree, k)
            push!(x, height[k])
            push!(y, depth[k])
            push!(l, k)
        end
    end
    x, y, l
end

function _p_circ(start_θ, end_θ, r=1)
    steps = range(start_θ, stop=end_θ, length = 1+ceil(Int, 60abs(end_θ - start_θ)))
    retx = Array{Float64}(undef, length(steps))
    rety = similar(retx)
    for u in eachindex(steps)
        retx[u] = _xcirc(steps[u], r)
        rety[u] = _ycirc(steps[u], r)
    end
    retx, rety
end

_xcirc(x, r) = r*cos(x)
_ycirc(y, r) = r*sin(y)
_tocirc(x, y) = _xcirc(y, x), _ycirc(y, x)

function _circle_transform_segments(xs, ys)
    retx, rety = Float64[], Float64[]
    function _transform_seg(_x, _y)
        tmpx, tmpy = _p_circ(_y[1], _y[2], _x[1])
        append!(retx, tmpx)
        append!(rety, tmpy)
        push!(retx, _xcirc(_y[3], _x[3]), NaN)
        push!(rety, _ycirc(_y[3], _x[3]), NaN)
    end
    i = 1
    while !(i === nothing) && i < length(xs)
        j = findnext(isnan, xs, i) - 1
        _transform_seg(view(xs,i:j), view(ys, i:j))
        i = j + 2
    end
    retx, rety
end

"""
    map_depthfirst(FUN, start, tree, eltype = nothing)

Apply `FUN` to each node in `tree` in depth-first order, and return the result.
`FUN` must take two arguments, `val` and `node`,  where `val` is the result of
applying `FUN` to the previous node, and `node` is the current node. `start`
specifies the initial value of `val`, and `eltype` specifies the type of the
return value of `FUN`.

### Examples
≡≡≡≡≡≡≡≡≡≡≡
Define a function to evolve a trait on the tree according to Brownian motion

julia> evolve(tree) = map_depthfirst((val, node) -> val + randn(), 0., tree, Float64)
"""
function map_depthfirst(FUN, start, tree, eltype = nothing)
    root = first(nodenamefilter(isroot, tree))
    eltype === nothing && (eltype = typeof(FUN(start, root)))
    ret = Vector{eltype}()
    function local!(val, node)
        push!(ret, val)
        for ch in getchildren(tree, node)
            local!(FUN(val, node), ch)
        end
    end
    local!(start, root)
    ret
end


text:
"""
    abstract type DiffusionParameters

Abstract type for all input structures to diffusion functions.
"""
abstract type DiffusionParameters end

"""
    DiffusionKernelParams

Structure that contains the dependent coefficients used in most of the diffusion theory calculations
"""
struct DiffusionKernelParams{T <: AbstractFloat} <: DiffusionParameters
    D::T                    # Diffusion coefficient
    z0::T                   # isotropic source depth
    ν::T                    # speed of light in medium
    zb::T                   # extrapolation length (EBC)
end
DiffusionKernelParams(μsp::T) where {T<:AbstractFloat} = DiffusionKernelParams{T}(D_coeff(μsp), z0_coeff(μsp), zero(eltype(T)), zero(eltype(T)))
DiffusionKernelParams(μsp::T, n_med::T) where {T<:AbstractFloat} = DiffusionKernelParams{T}(D_coeff(μsp), z0_coeff(μsp), ν_coeff(n_med), zero(eltype(T)))
DiffusionKernelParams(μsp::T, n_med::T, n_ext::T) where {T<:AbstractFloat} = DiffusionKernelParams{T}(D_coeff(μsp), z0_coeff(μsp), ν_coeff(n_med), zb_coeff(A_coeff(n_med / n_ext), D_coeff(μsp)))

"""
    A_coeff(n)

Function to approximate the A coefficient (mismatch between external and internal index of refractions) given in Eq. 25 from Ref. [1].
n is the ratio between the refractive index of the diffusing medium and that of the external medium (n_med / n_ext).
Uses the polynomial fit given in Appendix A [1]. Used to calculate the extrapolation length.

References:
[1] Daniele Contini, Fabrizio Martelli, and Giovanni Zaccanti,
    "Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,"
    Appl. Opt. 36, 4587-4599 (1997)
"""
function A_coeff(n)
    if n > 1.0
        A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
        5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
    elseif n < 1.0
        A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
    else
        A = 1.0
    end
    return A
end

"""
    D_coeff(μsp)

Returns the diffusion coefficient as a function of the reduced scattering coefficient.
Note that this is sometimes defined with the absorption coefficient D = 1 / (3 * (μsp + μa)).
"""
D_coeff(μsp) = 1 / (3 * μsp)

"""
    ν_coeff(n_med)

Computes the speed of light in the medium (cm/ns).
This is a good place to change the units of the returning function along with optical properties if interested in mm or ps.
"""
ν_coeff(n_med) = 29.9792458 / n_med

"""
    z0_coeff(μsp)

The location of the isotropic point source within the medium assuming some perpindicular incidence of light beam.
Note that this is sometimes defined with the absorption coefficient z0 = 1 / (μsp + μa)
"""
z0_coeff(μsp) = 1 / μsp

"""
    zb_coeff(A, D)

The extrapolation length calculated with the A coefficient (accounts for mismatch index of refraction at boundary)
and D (diffusion coefficient).
"""
zb_coeff(A, D) = 2 * A * D

max_issues_repo_name: heltonmc/LightPropagation


text:
"""
    System{N}

    graph structure representing a linear system ("Ax = b")

    matrix_entries: matrix entries
    vector_entries: residual entries
    diagonal_inverses: matrix inverses computed during LU factorization
    acyclic_children: contains direct children that are not part of a cycle
    cyclic_children: contains direct and indirect children that are part of a cycle
    parents: contains direct and cycle-opening parents
    dfs_list: list of graph nodes (depth first)
    graph: graph object
    dfs_graph: direct graph object
    dimrow: dimensions of row
    dimcol: dimensions of columns
"""
struct System{N}
    matrix_entries::SparseMatrixCSC{Entry, Int64}
    vector_entries::Vector{Entry}
    diagonal_inverses::Vector{Entry}
    acyclic_children::Vector{Vector{Int64}} # Contains direct children that are not part of a cycle
    cyclic_children::Vector{Vector{Int64}}  # Contains direct and indirect children that are part of a cycle
    parents::Vector{Vector{Int64}}          # Contains direct and cycle-opening parents
    dfs_list::SVector{N,Int64}
    graph::SimpleGraph{Int64}
    dfs_graph::SimpleDiGraph{Int64}
    dimrow::Vector{Int64}
    dimcol::Vector{Int64}

    function System{T}(A, dimrow, dimcol; force_static = false) where T
        N = length(dimrow)
        @assert N == length(dimcol)

        static = force_static || (all(dimrow.<=10) && all(dimcol.<=10))
        full_graph = Graph(A)
        matrix_entries = spzeros(Entry,N,N)

        for i = 1:N
            for j = 1:N
                if i == j
                    matrix_entries[i,j] = Entry{T}(dimrow[i], dimcol[j], static = static)
                elseif j ∈ all_neighbors(full_graph,i)
                    matrix_entries[i,j] = Entry{T}(dimrow[i], dimcol[j], static = static)
                    matrix_entries[j,i] = Entry{T}(dimrow[j], dimcol[i], static = static)
                end
            end
        end

        vector_entries = [Entry{T}(dim, static = static) for dim in dimrow];
        # this is only well-defined for dimrow == dimcol
        diagonal_inverses = [Entry{T}(dim, dim, static = static) for dim in dimrow];

        graphs, roots = split_adjacency(A)
        dfs_list = Int64[]
        acyclic_children = [Int64[] for i=1:N]
        cycles = [Vector{Int64}[] for i=1:N]
        parents = [Int64[] for i=1:N]
        edgelist = Graphs.SimpleEdge{Int64}[]

        for (i,graph) in enumerate(graphs)
            root = roots[i]
            dfs_graph = dfs_tree(graph, root)
            sub_dfs_list, cycle_closures = dfs(graph, root)
            append!(dfs_list, sub_dfs_list)

            cycle_dfs_graph = copy(dfs_graph)
            cycle_dfs_graph_reverse = copy(dfs_graph)
            for cycle_closure in cycle_closures
                add_edge!(cycle_dfs_graph, cycle_closure...)
                add_edge!(cycle_dfs_graph_reverse, reverse(cycle_closure)...)
            end
            append!(edgelist, collect(edges(cycle_dfs_graph_reverse)))

            for v in sub_dfs_list
                acyclic_children[v] = neighbors(dfs_graph, v)
                parents[v] = neighbors(reverse(dfs_graph), v)
            end

            for cyclic_members in simplecycles(cycle_dfs_graph)
                v, cyclic_children = cycle_parent_children(cyclic_members, parents)
                cyclic_children = reverse(cyclic_children)
                lump_cycles!(cycles[v], cyclic_children)

                acyclic_children[v] = setdiff(acyclic_children[v], cyclic_children)
                for c in cyclic_children
                    matrix_entries[v,c] = Entry{T}(dimrow[v], dimcol[c], static = static);
                    matrix_entries[c,v] = Entry{T}(dimrow[c], dimcol[v], static = static);

                    v ∉ parents[c] && push!(parents[c],v)
                end
            end
        end

        full_dfs_graph = SimpleDiGraph(edgelist)
        cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]

        new{N}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph, dimrow, dimcol);
    end

    System(A, dimrow, dimcol; force_static = false) = System{Float64}(A, dimrow, dimcol; force_static = force_static);
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N}) where {N}
    println(io, "System with "*string(N)*" nodes:")
    # SparseArrays._show_with_braille_patterns(io, system.matrix_entries) # TODO this fails on Julia 1.6
end

children(system, v) = outneighbors(system.dfs_graph, v)
connections(system, v) = neighbors(system.graph, v)
parents(system, v) = inneighbors(system.dfs_graph, v)