text:
using Knet
using ArgParse
using JLD
using HDF5
using JSON

include("lib/vocab.jl")
include("lib/base.jl")
include("lib/convnet.jl")
include("lib/init.jl")
include("lib/model.jl")
include("lib/data.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "Caption generation script for the model."

    @add_arg_table s begin
        ("--images"; help="data file contains images/features")
        ("--captions"; help="zip file contains caption data (karpathy)")
        ("--modelfile"; help="trained model file")
        ("--vocabfile"; help="vocabulary file")
        ("--cnnfile"; help="convnet file for non-finetuned model")
        ("--savefile"; help="save generations in JSON format")
        ("--beamsize"; arg_type=Int; default=1)
        ("--datasplit"; default="test"; help="data split is going to be used")
        ("--extradata"; action=:store_true)
        ("--maxlen"; arg_type=Int; default=25; help="max sentence length")
        ("--nogpu"; action=:store_true)
        ("--testing"; action=:store_true)
        ("--shuffle"; action=:store_true)
        ("--amount"; arg_type=Int; default=20;
         help="generation amount in case of testing")
        ("--feedback"; arg_type=Int; default=100;
         help="feedback in every N generation")
    end

    # parse args
    println("Datetime: ", now()); flush(STDOUT)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)

    savefile = abspath(o[:savefile])

    # load data
    vocab = load(o[:vocabfile], "vocab")
    entries = get_entries(o[:captions], [o[:datasplit]])[1]
    @printf("Data loaded [%s]\n", now()); flush(STDOUT)

    # load weights
    atype = o[:nogpu] ? Array{Float32} : KnetArray{Float32}
    w = load(o[:modelfile], "w")
    w = map(i->convert(atype, i), w)
    score = load(o[:modelfile], "bestscore")
    s = initstate(atype, size(w[3], 1), 1)

    wcnn = load(o[:modelfile], "wcnn")
    if wcnn != nothing
        wcnn = map(i->convert(atype, i), wcnn)
    end

    @printf("Model loaded [%s]\n", now()); flush(STDOUT)

    # generate captions
    ti = now()
    captions = []
    @printf("Generation started (score=%g,date=%s)\n", score, now())
    flush(STDOUT)
    for i = 1:length(entries)
        o[:testing] && counter >= o[:amount] && break

        entry = entries[i]
        filename = entry["filename"]
        image = h5open(o[:images], "r") do f
            read(f, filename)
        end
        caption = generate(w, wcnn, copy(s), image, vocab;
                           maxlen=o[:maxlen], beamsize=o[:beamsize])
        references = map(s->s["raw"], entry["sentences"])
        new_entry = Dict(
            "filename" => filename,
            "hypothesis" => caption,
            "references" => references)
        push!(captions, new_entry)

        if o[:feedback] > 0 && i % o[:feedback] == 0
            @printf("\n%d captions generated so far [%s]\n", i, now())
            flush(STDOUT)
        end
    end

    meta = Dict(
        "model" => o[:modelfile],
        "data" => o[:images],
        "split" => o[:datasplit],
        "score" => score,
        "beamsize" => o[:beamsize],
        "date" => now())

    open(savefile, "w") do f
        write(f, json(Dict("meta"=>meta,"captions"=>captions)))
    end

    tf = now()
    @printf("\nTime elapsed: %s [%s]\n", tf-ti, tf)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
#!/usr/bin/env julia
using Random
function getMass(nBodies)
    mList=zeros(nBodies)
    N=(1.5^(-1.3)-150^(-1.3))/1.3 #crude approximation of IMF integral assuming alpha = 2.3, stellar mass range of 0.5:150 solar masses
    rescale=1e6
    max=floor(Int,N*rescale)
    for i=1:nBodies
        intTarget=rand(0:max,1)[1]/rescale
        m=(1.5^(-1.3)-intTarget*1.3)^(-1/1.3) #just algebra from above
        mList[i]=round(m,digits=2)
    end
    return mList
end

ext: jl


text:
"""
    Samplers

A module in which we subtype 'Sampleable' to facilitate random draws.
"""
module Samplers
using Distributions: AbstractRNG, Discrete, Multivariate, rand, Sampleable, Univariate
using LightGraphs: AbstractGraph, neighbors, nv, vertices
using StochasticProcesses

export

# concrete types
ElementVectorSampler, EVecSampler, MultiSampler, NeighborSampler,

# methods
nsamplers, rand

# include scripts
include("elementvectorsamplers.jl")
include("neighborsamplers.jl")
include("multisamplers.jl")  # requires neighborsamplers.jl
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
using Oceananigans.Simulations:
    stop, iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded, TimeStepWizard

@testset "Simulations" begin
    @info "Testing simulations..."

    for arch in archs

        grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
        model = IncompressibleModel(architecture=arch, grid=grid)

        for Δt in (3, TimeStepWizard(Δt=5.0))
            model.clock.time = 0.0
            model.clock.iteration = 0

            simulation = Simulation(model, Δt=Δt, stop_iteration=1)

            # Just make sure we can construct a simulation without any errors.
            @test simulation isa Simulation

            @test iteration_limit_exceeded(simulation) == false
            @test stop(simulation) == false

            run!(simulation)

            # Just make sure run! executes without any errors.
            @test simulation isa Simulation

            # Some basic tests
            @test iteration_limit_exceeded(simulation) == true
            @test stop(simulation) == true

            t = Δt isa Number ? 3 : 5
            @test model.clock.time ≈ t
            @test model.clock.iteration == 1
            @test simulation.run_time > 0

            @test stop_time_exceeded(simulation) == false
            simulation.stop_time = 1e-12
            @test stop_time_exceeded(simulation) == true

            @test wall_time_limit_exceeded(simulation) == false
            simulation.wall_time_limit = 1e-12
            @test wall_time_limit_exceeded(simulation) == true
        end
    end
end


text:
# This script contains the plotting of timer output
using Plots

function plotTimerOutput(archs, dims, to::TimerOutput, modelName)

    x = fill(0, length(dims))
    y = fill(0.0, length(dims))

    # Build a list of architectures in a single string, to include in the title
    archList = ""
    for arch in archs
        if arch != archs[end]
            archList *= arch*" vs. "
        else
            archList *= arch*" "
        end
    end

    # Make the plot title
    plotTitle = "Performance comparison for " * archList * "for "* modelName * " model"

    # Start the empty plot with title
    p = plot(title=plotTitle)

    # Create placeholder for the reference time (used for computing the speed-up)
    refTime = 0.0

    # For now allow up to 5 different architectures.
    # When more are needed, expand this list with more markers.
    mark = [(:circle), (:square), (:cross), (:pentagon), (:hexagon)]

    # Now loop over all architecture and grid resolutions to build up the plot
    archIdx = 0
    for arch in archs
        archIdx += 1
        dimIdx = 0
        for dim in dims
            dimIdx += 1
            dimstr = string(dim)
            time = TimerOutputs.time(to[arch][dimstr])

            # Determine total number of cells for this resolution
            numcells = 1
            for idim in dim
                numcells *=idim
            end
            x[dimIdx] = numcells
            # Time was given in nano-seconds: convert to seconds
            y[dimIdx] = Float64(time)/1.0e9
        end
        if arch == archs[1]
            refTime = copy(y)
        end
        p = plot!(x,y, marker=mark[archIdx], size=(800,800), xlabel="# of cells", ylabel="Time (s)", label=arch, legend=(0.2, 0.9), grid=true)
    end

    @show SpeedUp = refTime./y
    p = plot!(x,SpeedUp, marker=:diamond, size=(800,800), xlabel="# of cells", ylabel="Time (s)", label="Architectural speed up", legend=(0.2, 0.9), grid=true)

    #display(p)
    figureTitle = "performance_plot_" * modelName * ".png"
    #savefig(p, figureTitle)

end


text:
# Autogenerated wrapper script for HGSCVRP_jll for i686-w64-mingw32-cxx03
export hgs, libhgscvrp

JLLWrappers.@generate_wrapper_header("HGSCVRP")
JLLWrappers.@declare_library_product(libhgscvrp, "libhgscvrp.dll")
JLLWrappers.@declare_executable_product(hgs)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libhgscvrp,
        "bin\\libhgscvrp.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        hgs,
        "bin\\hgs.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 7.
#

fl = Float64

function f(x::fl)::fl
    return fl(sin(x) + cos(3*x))
end

function f_derivative(x::fl)::fl
    return fl(cos(x) - 3 * sin(3*x))
end

# break out this part of the approximate function for later use
function _approximate_derivative_internal(func, x::fl, h::fl)::fl
    return fl(func(x + h) - func(x))
end
function approximate_derivative_at_given_point(func, x::fl, h::fl)::fl
    return fl( _approximate_derivative_internal(func, x, h) / h )
end

# calculate the results and print them out
for n in -25:-1:-32
    h = fl(2.0^n)
    x = fl(1.0)
    appr = approximate_derivative_at_given_point(f, x, h)
    exact = f_derivative(x)

    println(
        "n = ", n, "\n",
        "h                         = ", h, "\n",
        "f(x + h) - f(x)           = ", _approximate_derivative_internal(f, x, h), "\n",
        "\\tilde{f’}(x)             = ", appr, "\n",
        "f’(x)                     = ", exact, "\n",
        "|f’(x) - \\tilde{f’}(x)|   = ", abs( fl(exact - appr) ),
        "\n"
    )
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# It implements the Actor-model
#

# --------------------------------
# registry API
# --------------------------------

"""
    register(name::Symbol, lk::Link)

Register the actor `lk` with `name`. Returns `true` if the
registration succeeds, `false` if `name` is already in use.
"""
function register(name::Symbol, lk::Link)
    res = myid() == 1 ?
        call!(_REG, register, name, lk) :
        call!(_REG, register, name, Link(RemoteChannel(()->lk.chn),myid(),:remote))
    update!(lk, name, s=:name)
    return res
end

"""
    unregister(name::Symbol)

Remove any registrations associated with `name`.
"""
function unregister(name::Symbol)
    lk = whereis(name)
    if !ismissing(lk)
        call!(_REG, unregister, name)
        update!(lk, nothing, s=:name)
    end
end

"""
    whereis(name::Symbol)

Find out whether name is registered. Return the actor link
`lk` or `missing` if not found.
"""
whereis(name::Symbol) = call!(_REG, whereis, name)

"""
    registered()

Return an Array of all registered actors in the system.
"""
registered() = call!(_REG, registered, myid())

# behavior functions
function _reg(d::Dict{Symbol,Link}, ::typeof(register), name::Symbol, lk::L) where L<:Link
    if !haskey(d, name)
        d[name] = lk
        return true
    else
        return false
    end
end
_reg(d::Dict{Symbol,Link}, ::typeof(unregister), name::Symbol) = delete!(d, name)
_reg(d::Dict{Symbol,Link}, ::typeof(whereis), name::Symbol) = get(d, name, missing)
_reg(d::Dict{Symbol,Link}, ::typeof(registered), id::Int) =
    id == 1 ?   collect(pairs(d)) :
                [Pair(i[1], _rlink(i[2])) for i in pairs(d)]


text:
"""
    getTRMMlsm(
        geo :: GeoRegion = GeoRegion("GLB");
	    sroot :: AbstractString = homedir()
    ) -> nothing

Downloads the TRMM-TMPA Land-Sea Mask for a GeoRegion specified by `geo` and saves it into `sroot`

Arguments
=========
- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat

Keyword Arguments
=================
- `sroot` : the directory to save the TRMM-TMPA Land-Sea Mask, default is the home directoy called by `homedir()`
"""
function getTRMMlsm(
	geo   :: GeoRegion = GeoRegion("GLB");
    sroot :: AbstractString = homedir(),
)

    @info "$(modulelog()) - Retrieving Land-Sea Mask for the TRMM-TMPA dataset"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the TRMM LandSea Mask GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("TRMMLSM")
	else
		isinGeoRegion(geo,GeoRegion("TRMMLSM"))
	end

    lon,lat = trmmlonlat(full=true); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary array for extraction of TRMM-TMPA Land-Sea Mask data for the $(geo.name) GeoRegion from the original global gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp  = zeros(Float32,nlat,nlon)
	var  = zeros(Float32,nglon,nglat)

	@info "$(modulelog()) - Retrieving the original TRMM-TMPA Land-Sea Mask data from NASA's EOSDIS OPeNDAP servers"
	hroot = "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/AUXILIARY"
    npdnc = "TRMM_TMPA_LandSeaMask.2/TRMM_TMPA_LandSeaMask.2.nc4"
    npdds = NCDataset(joinpath(hroot,npdnc))
	NCDatasets.load!(npdds["landseamask"].var,tmp,:,:)
	close(npdds)

	@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"

	extractregionlsm!(var,tmp,ginfo)

	saveTRMMlsm(var,geo,ginfo,sroot)

end

function saveTRMMlsm(
	lsm   :: AbstractArray{<:Real,2},
	geo   :: GeoRegion,
	ginfo :: RegionGrid,
	sroot :: AbstractString
)

	@info "$(modulelog()) - Saving TRMM-TMPA Land-Sea Mask data in the $(geo.name) GeoRegion"

	if !isdir(sroot); mkpath(sroot) end
	fnc = joinpath(sroot,"trmmlsm-$(geo.regID).nc")
	if isfile(fnc)
		@info "$(modulelog()) - Overwrite stale NetCDF file $(fnc) ..."
		rm(fnc);
	end

	@info "$(modulelog()) - Creating NetCDF file $(fnc) ..."
	ds = NCDataset(fnc,"c",attrib = Dict(
		"Conventions" => "CF-1.4",
		"NCO"         => "4.5.1",
		"Authors"     => "Bill Olson, Dave Bolvin, George Huffman",
		"Description" => "This TMPA land sea mask contains values ranging from 0% to 100% with 0% representing all land and 100% representing all ocean",
		"Title"       => "Land/Sea static mask relevant to TMPA precipitation 0.25x0.25 degree",
		"Version"     => "2",
		"History"     => "This land sea mask originated from the NOAA group at SSEC in the 1980s. It was originally produced at 1/6 deg resolution, and then regridded for the purposes of GPCP, TMPA, and IMERG precipitation products. NASA code 610.2, Global Change Data Center, restructured this TMPA land sea mask to match the TMPA grid, and converted the file to CF-compliant netCDF4. Version 2 was created in May, 2019 to resolve detected inaccuracies in coastal regions.",
		"DOI"         => "10.5067/K1VR7BN7YDWN",
	))

	ds.dim["longitude"] = length(ginfo.glon)
	ds.dim["latitude"]  = length(ginfo.glat)

	nclon = defVar(ds,"longitude",Float32,("longitude",),attrib = Dict(
	    "units"     => "degrees_east",
	    "long_name" => "longitude",
	))

	nclat = defVar(ds,"latitude",Float32,("latitude",),attrib = Dict(
	    "units"     => "degrees_north",
	    "long_name" => "latitude",
	))

	nclsm = defVar(ds,"lsm",Float32,("longitude","latitude",),attrib = Dict(
	    "units"         => "0-1",
	    "long_name"     => "land_sea_mask",
		"full_name"     => "Land-Sea Mask",
	))

	nclon[:] = ginfo.glon
	nclat[:] = ginfo.glat
	nclsm[:] = 1 .- lsm/100

	close(ds)

	@info "$(modulelog()) - TRMM-TMPA Land-Sea Mask data in the $(geo.name) GeoRegion has been saved into $(fnc)"

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using Test, CUDA, Random, Knet, AutoGrad
using CUDA.CUDNN:
    cudnnConvolutionForward,
    cudnnConvolutionForward!,
    cudnnConvolutionBackwardFilter,
    cudnnConvolutionBackwardData,
    cudnnGetConvolutionNdForwardOutputDim,
    cudnnSetConvolutionMathType,
    cudnnSetConvolutionReorderType,
    cudnnSetConvolutionGroupCount,
    cudnnFindConvolutionForwardAlgorithmEx,
        cudnnConvolutionFwdAlgoPerf_t,
    cudnnFindConvolutionBackwardFilterAlgorithmEx,
        cudnnConvolutionBwdFilterAlgoPerf_t,
    cudnnFindConvolutionBackwardDataAlgorithmEx,
        cudnnConvolutionBwdDataAlgoPerf_t,
    cudnnConvolutionDescriptor,
        cudnnConvolutionDescriptor_t,
        cudnnCreateConvolutionDescriptor,
        cudnnSetConvolutionNdDescriptor,
        cudnnDestroyConvolutionDescriptor,
    cudnnConvolutionMode_t,
        CUDNN_CONVOLUTION,       # 0
        CUDNN_CROSS_CORRELATION, # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnMathType_t,
        CUDNN_DEFAULT_MATH,                    # 0
        CUDNN_TENSOR_OP_MATH,                  # 1
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION, # 2
        CUDNN_FMA_MATH,                        # 3
    cudnnReorderType_t,
        CUDNN_DEFAULT_REORDER, # 0
        CUDNN_NO_REORDER,      # 1
    cudnnConvolutionFwdAlgo_t,
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,         # 0
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM, # 1
        CUDNN_CONVOLUTION_FWD_ALGO_GEMM,                  # 2
        CUDNN_CONVOLUTION_FWD_ALGO_DIRECT,                # 3
        CUDNN_CONVOLUTION_FWD_ALGO_FFT,                   # 4
        CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING,            # 5
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD,              # 6
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED,     # 7
        CUDNN_CONVOLUTION_FWD_ALGO_COUNT,                 # 8
    cudnnConvolutionBwdFilterAlgo_t,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3,                 # 3, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD,          # 4, /* not implemented */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING,        # 6,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT,             # 7
    cudnnConvolutionBwdDataAlgo_t,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING,        # 3,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD,          # 4,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT,             # 6
    cudnnDataType,
    convdims,
    math_mode,
    handle

@testset "cudnn/convolution" begin
    T = Float64
    ax,aw,ag,ab = randn(T,8,8,4,4),randn(T,3,3,4,4),randn(T,3,3,2,4),randn(T,1,1,4,1)
    cx,cw,cg,cb = Param.(CuArray.((ax,aw,ag,ab)))
    cz = Param(cudnnConvolutionForward(cw,cx))

    # These call cudnnConvolutionForward
    @test @gcheck cudnnConvolutionForward(cw,cx)
    @test @gcheck cudnnConvolutionForward(cw,cx;padding=1)
    @test @gcheck cudnnConvolutionForward(cw,cx;stride=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;dilation=2)
    @test @gcheck cudnnConvolutionForward(cg,cx;group=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;mathType=CUDNN_DEFAULT_MATH)
    @test @gcheck cudnnConvolutionForward(cw,cx;mathType=CUDNN_TENSOR_OP_MATH)
    @test @gcheck cudnnConvolutionForward(cw,cx;mathType=CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION)
    @test @gcheck cudnnConvolutionForward(cw,cx;reorderType=CUDNN_NO_REORDER)
    @test @gcheck cudnnConvolutionForward(cw,cx;alpha=2)

    # These call cudnnConvolutionBiasActivationForward
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb)
    @test @gcheck cudnnConvolutionForward(cw,cx;z=cz)
    @test @gcheck cudnnConvolutionForward(cw,cx;z=cz,beta=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;activation=CUDNN_ACTIVATION_RELU)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,z=cz)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,activation=CUDNN_ACTIVATION_RELU)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,padding=1)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,stride=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,dilation=2)
    @test @gcheck cudnnConvolutionForward(cg,cx;bias=cb,group=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,mathType=CUDNN_DEFAULT_MATH)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,mathType=CUDNN_TENSOR_OP_MATH)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,mathType=CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,reorderType=CUDNN_NO_REORDER)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,alpha=2)
    @test @gcheck cudnnConvolutionForward(cw,cx;bias=cb,beta=2,z=cz)
end


text:
x = range(0, stop=1, length=20)
plot(x, x.^2, x.^3, lw=0.05, ls=["--", "-"], marker=["wo", "none"],
     keys=["quadratic", "cubic"],)
scatter!(x, sin.(x), key=["sine"])
legend(; position="top-left")


text:
using JpegTurbo
using JpegTurbo.LibJpeg
using Test

@testset "JpegTurbo.jl" begin
    @testset "config" begin
        @test_nowarn JpegTurbo.versioninfo()

        # ensure we're using SIMD-built version for performance
        @test LibJpeg.WITH_SIMD == 1

        # use the 8-bit version
        @test LibJpeg.JSAMPLE == UInt8
        @test LibJpeg.BITS_IN_JSAMPLE == 8
        @test LibJpeg.MAXJSAMPLE == 255
        @test LibJpeg.CENTERJSAMPLE == 128
    end
end


text:
using Documenter, ImageBase, ImageShow

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true")

makedocs(modules  = [ImageBase, ],
         format   = format,
         sitename = "ImageBase",
         pages    = ["index.md", "reference.md"])

deploydocs(repo   = "github.com/JuliaImages/ImageBase.jl.git")


text:
"""
    Recombinations

Basis recombination module.

Defines [`RecombinedBSplineBasis`](@ref) and [`RecombineMatrix`](@ref) types.
"""
module Recombinations

export RecombinedBSplineBasis, RecombineMatrix
export recombination_matrix, nzrows, constraints

import Base: @propagate_inbounds
import LinearAlgebra

using LinearAlgebra: UniformScaling, ldiv!, dot
using StaticArrays

using ..BSplines
using ..DifferentialOps
using ..Splines

# These are redefined for RecombinedBSplineBasis.
import ..BSplines:
    boundaries, knots, order, evaluate, nonzero_in_segment, support, summary_basis

include("matrices.jl")
include("bases.jl")
include("splines.jl")

end


text:
module SimpleADMM

import Printf: @sprintf
import SimpleNL: Expression, Model, variable, parameter, constraint, objective, num_variables, num_constraints, optimize!, instantiate!, index,  non_caching_eval
import SimpleNLUtils: get_terms, get_entries_expr, sparsity, KKTErrorEvaluator
import Requires: @require

default_subproblem_optimizer() = @isdefined(DEFAULT_SUBPROBLEM_OPTIMIZER) ? DEFAULT_SUBPROBLEM_OPTIMIZER : error("DEFAULT_SUBPROBLEM_OPTIMIZER is not defined. To use Ipopt as a default subproblem optimizer, do: using Ipopt")
default_option() = Dict(
    :rho=>1.,
    :maxiter=>400,
    :tol=>1e-6,
    :subproblem_optimizer=>default_subproblem_optimizer(),
    :subproblem_option=>Dict(:print_level=>0),
    :save_output=>false,
    :optional=>admm->nothing
)

mutable struct SubModel
    model::Model

    x_V_orig
    x_V_sub
    l_V_orig
    l_V_sub

    z_orig
    z_sub
    x_sub
    l_sub

    function SubModel(m::Model,optimizer,V,rho,
                      objs,cons,obj_sparsity,con_sparsity;opt...)
        msub = Model(optimizer;opt...)

        V_bdry = Int[]

        V_bool = falses(num_variables(m))
        V_bool[V] .= true

        for vs in con_sparsity
            if examine(V_bool,vs) == 1
                union!(V_bdry,vs)
            end
        end

        sort!(V_bdry)

        V_inner = setdiff(V,V_bdry)
        V_con = [i for i in eachindex(con_sparsity) if examine(V_bool,con_sparsity[i]) >= 1]

        x = Dict{Int,Expression}()
        z = Dict{Int,Expression}()
        l = Dict{Int,Expression}()

        for i in V_inner
            x[i]= variable(msub;lb=m.xl[i],ub=m.xu[i],start=m.x[i])
        end
        for i in V_bdry
            x[i]= variable(msub;lb=m.xl[i],ub=m.xu[i],start=m.x[i])
            z[i]= parameter(msub)
        end
        for i in V_bdry
            l[i]= parameter(msub)
        end
        for i in V_con
             constraint(msub,non_caching_eval(cons[i],x,m.p);lb=m.gl[i],ub=m.gu[i])
        end

        objective(msub,sum(non_caching_eval(objs[i],x,m.p) for i in eachindex(objs) if examine(V_bool,obj_sparsity[i]) >= 2) + sum((x[i]-z[i]) * l[i] + 0.5 * rho * (x[i]-z[i])^2 for i in V_bdry))

        instantiate!(msub)

        x_V_orig = view(m.x,V)
        x_V_sub = view(msub.x,([x[i].index for i in V]))
        l_V_orig = view(m.l,V_con)
        l_V_sub = msub.l

        z_orig = view(m[:z],V_bdry)
        z_sub = view(msub.p,1:length(V_bdry))
        x_sub = view(msub.x,length(V_inner)+1:length(V_inner)+length(V_bdry))
        l_sub = view(msub.p,1+length(V_bdry):2*length(V_bdry))

        return new(msub, x_V_orig, x_V_sub, l_V_orig, l_V_sub, z_orig, z_sub, x_sub, l_sub)
    end
end

mutable struct Optimizer
    model::Model
    submodels::Vector{SubModel}
    kkt_error_evaluator
    opt::Dict{Symbol,Any}
end

function primal_update!(z,z_counter)
    @inbounds @simd for i in eachindex(z)
        z_counter[i] != 0 && (z[i] /= z_counter[i])
    end
end

function optimize!(admm::Optimizer)
    save_output = admm.opt[:save_output]
    if save_output
        output = Tuple{Float64,Float64}[]
        start = time()
    end

    iter = 0
    while (err=admm.kkt_error_evaluator(admm.model.x,admm.model.l,admm.model.gl)) > admm.opt[:tol] && iter < admm.opt[:maxiter]
        save_output && push!(output,(err,time()-start))
        println(@sprintf "%4i %4.2e" iter+=1 err)

        Threads.@threads for sm in admm.submodels
            sm.z_sub .= sm.z_orig
            sm.l_sub .+= admm.opt[:rho] .* (sm.x_sub .- sm.z_sub)
            optimize!(sm)
            sm.x_V_orig .= sm.x_V_sub
        end
        admm.model.l.=0
        for sm in admm.submodels
            sm.l_V_orig .+= sm.l_V_sub
        end
        admm.model[:z] .= 0
        for sm in admm.submodels
            sm.z_orig .+= sm.x_sub + sm.l_sub ./ 2 ./ admm.opt[:rho]
        end
        primal_update!(admm.model[:z],admm.model[:z_counter])
    end
    save_output && push!(output,(err,time()-start))
    println(@sprintf "%4i %4.2e" iter+=1 err)
    save_output && (admm.model.ext[:output]=output)
    admm.opt[:optional](admm)
end

function difference(a,b)
    diff = .0
    @inbounds @simd for i in eachindex(a)
        diff = max(diff,abs(a[i] - b[i]))
    end
    return diff
end

function examine(W_bool,keys)
    numtrue = 0
    for k in keys
        W_bool[k] && (numtrue += 1)
    end
    return numtrue == 0 ? 0 : numtrue == length(keys) ? 2 : 1
end

optimize!(sm::SubModel) = optimize!(sm.model)
instantiate!(sm::SubModel) = instantiate!(sm.model)

function Optimizer(m::Model)

    opt = default_option()
    for (sym,val) in m.opt
        opt[sym] = val
    end

    objs = get_terms(m.obj)
    cons = get_entries_expr(m.con)
    obj_sparsity = [sparsity(e) for e in objs]
    con_sparsity = [sparsity(e) for e in cons]

    m[:z] = zeros(num_variables(m))
    m[:z_counter] = zeros(Int,num_variables(m))

    sms = Vector{SubModel}(undef,length(m[:Vs]))

    K = [k for k in keys(m[:Vs])]
    Threads.@threads for k in K
        sms[k] = SubModel(m,opt[:subproblem_optimizer],m[:Vs][k],opt[:rho],
                          objs,cons,obj_sparsity,con_sparsity;
                          opt[:subproblem_option]...)
    end

    for k in K
        m[:z_counter][sms[k].z_orig.indices[1]] .+= 1
    end

    Optimizer(m,sms,KKTErrorEvaluator(m),opt)
end

function __init__()
    @require Ipopt="b6b21f68-93f8-5de0-b562-5493be1d77c9" @eval begin
        import ..Ipopt
        DEFAULT_SUBPROBLEM_OPTIMIZER = Ipopt.Optimizer
    end
end

end # module


text:
function getFavorites(config::Config, key::String)
    favorites = get(config, key, Array{String, 1}())
    config[key] = favorites

    return favorites
end

function addFavorite(config::Config, favoritesKey::String, favorite::String, favorites::Array{T, 1}=getFavorites(config, favoritesKey)) where T
    if !(favorite in favorites)
        push!(favorites, favorite)
    end
end

function removeFavorite(config::Config, favoritesKey::String, favorite::String, favorites::Array{T, 1}=getFavorites(config, favoritesKey)) where T
    filter!(f -> f != favorite, favorites)
end

function toggleFavorite(config::Config, favoritesKey::String, favorite::String)
    favorites = getFavorites(config, favoritesKey)

    if favorite in favorites
        removeFavorite(config, favoritesKey, favorite, favorites)

    else
        addFavorite(config, favoritesKey, favorite, favorites)
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
using DelimitedFiles

const data = readdlm("day11.txt", String)

# floor will be 0, empty will be 1 and occupied will be 2

function create_matrix(data::Array{<:AbstractString,2})
    # We will pad the data with 0s around
    height = length(data) + 2
    width = length(data[1]) + 2
    mat = zeros(Int, (height, width))
    for (i, line) in enumerate(data)
        for (j, c) in enumerate(data[i])
            if c == 'L'
                # +1 because padding
                mat[i + 1 , j + 1] = 1
            end
        end
    end
    return mat
end

function process(mat::Array{Int,2}, i::Int, j::Int)
    state = mat[i, j]
    # empty
    if state == 1
        # if any is occupied, stay empty, else occupied
        new_state = any(mat[xi, xj] == 2 for xi ∈ i - 1:i + 1, xj ∈ j - 1:j + 1) ? 1 : 2
    # occupied
    elseif state == 2
        # If 4 adjacent (i.e. 5 including itself) are occupied, becomes empty
        new_state = sum(mat[xi, xj] == 2 for xi ∈ i - 1:i + 1, xj ∈ j - 1:j + 1) >= 5 ? 1 : 2
    else
        new_state = state
    end
    return new_state
end

function one_pass(mat::Array{Int,2})
    h, w = size(mat)
    new_mat = zeros(Int, (h, w))
    for i ∈ 2:h - 1, j ∈ 2:w - 1
        new_mat[i, j] = process(mat, i, j)
    end
    return new_mat
end

function wait_for_stable(mat::Array{Int,2})
    new_mat = one_pass(mat)
    while new_mat != mat
        mat = new_mat
        new_mat = one_pass(mat)
    end
    return count(==(2), mat)
end

const mat = create_matrix(data)

println("Part 1: $(wait_for_stable(mat))")

# Now we want the first non empty seat in each of the 8 directions

function look_in_direction(mat::Array{Int,2}, i::Int, j::Int, xi::Int, xj::Int)
    maxi, maxj = size(mat)
    while true
        new_i, new_j = i + xi, j + xj
        if new_i ∈ [1, maxi] || new_j ∈ [1, maxj]
            # we reached the end
            return 0
        end
        if mat[new_i, new_j] != 0
            return mat[new_i, new_j]
        end
        i, j = new_i, new_j
    end
end

function process_p2(mat::Array{Int,2}, i::Int, j::Int)
    state = mat[i, j]
    # empty
    if state == 1
        # if any is occupied, stay empty, else occupied
        new_state = any(look_in_direction(mat, i, j, xi, xj) == 2 for xi ∈ - 1:1, xj ∈ - 1:1) ? 1 : 2
    # occupied
    elseif state == 2
        # If 5 adjacent (i.e. 6 including itself) are occupied, becomes empty
        new_state = sum(look_in_direction(mat, i, j, xi, xj) == 2 for xi ∈ - 1:1, xj ∈ - 1:1) >= 6 ? 1 : 2
    else
        new_state = state
    end
    return new_state
end

function one_pass_p2(mat::Array{Int,2})
    h, w = size(mat)
    new_mat = zeros(Int, (h, w))
    for i ∈ 2:h - 1, j ∈ 2:w - 1
        new_mat[i, j] = process_p2(mat, i, j)
    end
    return new_mat
end

function wait_for_stable_p2(mat::Array{Int,2})
    new_mat = one_pass(mat)
    while new_mat != mat
        mat = new_mat
        new_mat = one_pass_p2(mat)
    end
    return count(==(2), mat)
end

println("Part 2: $(wait_for_stable_p2(mat))")


text:
immutable Spec2D{T<:Real,D<:PlaneDomains,
        G<:PlaneGroups,W<:RedundantWaveletClass} <: AbstractSpec{T,D,G,W}
    ɛ::Float64
    ϕmeta::ΦMeta
    ψmetas::Array{ΨMeta2D,3}
    class::W
    domain::D
    log2_size::Int
    max_aspectratio::Float64
    max_qualityfactor::Float64
    max_scale::Float64
    motherfrequency::Float64
    n_filters_per_octave::Int
    n_octaves::Int
    pointgroup::G
    signaltype::Type{T}

    function (::Type{Spec2D}){T,D,W}( ;
            class::W = Morlet(),
            signaltype::Type{T} = Float32,
            domain::D = FourierDomain(2),
            ɛ = default_ɛ(signaltype),
            log2_size = 5,
            max_aspectratio = nothing,
            max_qualityfactor = nothing,
            max_scale = Inf,
            n_filters_per_octave = nothing,
            n_octaves = nothing,
            n_orientations = nothing)
        n_orientations = default_n_orientations(class, n_orientations)
        """Infer point group from wavelet."""
        pointgroup = issteerable(class) ?
            RotationGroup(n_orientations) : TrivialGroup()
        """Maximum aspect ratio (length-to-width) of the wavelets."""
        max_aspectratio = default_max_aspectratio(class, max_aspectratio)
        """Maximum quality factor and number of filters per octave are
        set simultaneously, so that they default to each other if either
        is present."""
        max_qualityfactor, n_filters_per_octave =
            default_max_qualityfactor(max_qualityfactor, n_filters_per_octave),
            default_n_filters_per_octave(n_filters_per_octave, max_qualityfactor)
        """Center frequency of the mother wavelet (`motherfrequency`)."""
        motherfrequency = default_motherfrequency(class, n_filters_per_octave)
        """By default, the number of octaves in the filter bank is such that
        all wavelet scales do not exceed the signal length, nor the
        user-defined maximum scale `max_scale` if present."""
        n_octaves = default_n_octaves(n_octaves, class, log2_size,
            max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave)
        """The number of orientations of a one-dimensional filter bank is
        equal to 1 for non-steerable real wavelets (`pointgroup::TrivialGroup`),
        equal to n_orientations otherwise (`pointgroup::RotationGroup`)."""
        nΘs = get_n_orientations(pointgroup)
        ħ = uncertainty(class)
        """The meta-information of the wavelets is computed as a 3-d array,
        whose dimensions respectively correspond to
        * the orientation variable `θ`,
        * the chroma variable `χ`, and
        * the octave variable `j`."""
        ψmetas = Array{ΨMeta2D}(nΘs, n_filters_per_octave, n_octaves)
        for j in 0:(n_octaves-1), χ in 0:(n_filters_per_octave-1)
            γ = j * n_filters_per_octave + χ
            resolution = exp2(-γ / n_filters_per_octave)
            centerfrequency = motherfrequency * resolution
            unbounded_scale =
                ħ * max_aspectratio * max_qualityfactor / centerfrequency
            scale = min(unbounded_scale, max_scale)
            unbounded_q = scale * centerfrequency / (ħ * max_aspectratio)
            qualityfactor = clamp(unbounded_q, 1.0, max_qualityfactor)
            unbounded_scale =
                ħ * max_aspectratio * qualityfactor / centerfrequency
            scale = min(unbounded_scale, max_scale)
            aspectratio = scale * centerfrequency / (ħ * scale)
            bandwidth = centerfrequency / qualityfactor
            for θ in 0:(nΘs-1)
                ψmetas[1+θ, 1+χ, 1+j] = ΨMeta2D(γ, θ, χ,
                    aspectratio, bandwidth, centerfrequency, j,
                    qualityfactor, scale)
            end
        end
        """The bandwidth of the lowpass filter `ϕ` is such that the Fourier
        transform of `ϕ` will cover the low-frequency part of the spectrum
        that is not covered by the wavelets `ψs`."""
        ϕbandwidth = motherfrequency * exp2(-n_octaves)
        ϕscale = ħ / ϕbandwidth
        ϕmeta = ΦMeta(ϕbandwidth, ϕscale)
        G = typeof(pointgroup)
        spec = new{T,D,G,W}(ɛ, ϕmeta, ψmetas, class, domain, log2_size,
            max_aspectratio, max_qualityfactor, max_scale, motherfrequency,
            n_filters_per_octave, n_octaves, pointgroup, signaltype)
        """Before returning the `spec`, we call the function `checkspec`
        which enforces properties of the wavelet filter bank to satisfy null
        mean, limited spatial support, and Littlewood-Paley inequality."""
        checkspec(spec) && return spec
    end
end

"""Enforces properties of two-dimensional wavelets.

* For steerable wavelets, such as `Morlet`, `max_aspectratio` must be `≧1.0` and
`n_orientations` must be `≧2`.

* For non-steerable wavelets, such as `MexicanHat`, `max_aspectratio` must be
equal to `1.0` and `n_orientations` must be equal to `1`."""
function checkspec(spec::Spec2D)
    if issteerable(spec.class)
        if spec.max_aspectratio < 1.0
            error("`max_aspectratio` must be `≧1.0`",
                "for steerable wavelets.")
        end
        if get_n_orientations(spec.pointgroup) < 2
            error("`n_orientations` must be `≧2` for steerable wavelets")
        end
    else
        if spec.max_aspectratio != 1.0
            error("`max_aspectratio` must be equal to `1.0`",
                "for non-steerable wavelets.")
        end
        if get_n_orientations(spec.pointgroup) != 1
            error("`n_orientations` must be equal to `1`",
                "for non-steerable wavelets")
        end
    end
    checkspec_super(spec) && return true
end

"""By default, the maximum aspect ratio is set to `2.0` for steerable
wavelets, such as `Morlet`. For non-steerable wavelets, such as `MexicanHat`,
the maximum aspect ratio must be equal to `1.0`."""
default_max_aspectratio(class::RedundantWaveletClass, max_aspectratio::Void) =
    issteerable(class) ? 2.0 : 1.0
default_max_aspectratio(class::RedundantWaveletClass, max_aspectratio::Any) =
    Float64(max_aspectratio)

"""By default, the number of orientations is set to `4` for steerable
wavelets, such as `Morlet`. For non-steerable wavelets, such as `MexicanHat`,
the number of orientations must be equal to `1`."""
default_n_orientations(class::RedundantWaveletClass, n_orientations::Void) =
    issteerable(class) ? 4 : 1

function default_n_orientations(class::RedundantWaveletClass, n_orientations::Any)
    if issteerable(class)
        return Int(n_orientations)
    else
        error("`n_orientations` must be equal to `1` for steerable wavelets")
    end
end


text:
"""
`DEFAULT_PARAMS` - Create dictionary of default parameters for blockmatching
"""
function default_params()
  return Dict(  "mesh_dist" => 750,
                "block_size" => 200,
                "search_r" => 150,
                "min_r" => 0.25)
end

"""
`BLOCKMATCH` - Build mesh to produce nodes, and blockmatch between two images.
"""
function blockmatch(src_img, dst_img; src_offset=[0,0], dst_offset=[0,0],
                                                      params=default_params())
  mesh = create_mesh(src_img, src_offset, params["mesh_dist"])
  matches = blockmatch(mesh.src_nodes, src_img, dst_img,
                                        src_offset, dst_offset, params)
  return mesh, matches
end

"""
`MATCHES_TO_MESH` - Using matches and old mesh, create new mesh with just nodes
  related to the matches.

```
new_mesh = matches_to_mesh(matches::Matches, mesh::Mesh)
```

* matches: set of Matches
* mesh: Mesh from which the Matches were generated
* new_mesh: adjusted Mesh using matches.mesh_indices to filter only the nodes
    where matches were found. The edge-node incidence matrix will also be
    adjusted.
"""
function matches_to_mesh(matches::Matches, mesh::Mesh)
  assert(length(matches.mesh_indices) == size(matches.dst_points, 1))
  new_mesh = Mesh()
  mask = matches.mesh_indices
  new_mesh.src_nodes = mesh.src_nodes[mask, :]
  new_mesh.dst_nodes = matches.dst_points
  edges = mesh.edges[mask, :]
  non_orphans = find_zero_indices(sum(edges, 1))
  new_mesh.edges = edges[:, non_orphans]
  return new_mesh
end

function meshwarp(img, mesh::Mesh, offset=[0,0])
  triangles = incidence_to_triangles(mesh.edges)
  return @time meshwarp(img, mesh.src_nodes, mesh.dst_nodes, triangles, offset)
end

"""
Transform Nx2 pts by 3x3 tform matrix
"""
function warp_pts(tform, pts)
    pts = hcat(pts, ones(size(pts,1)))
    tpts = pts * tform
    return tpts[:,1:2]
end

"""
Create array of UInt8 type from from image path
"""
function load_uint8_image(path::AbstractString)
  img = Images.load(path)
  return reinterpret(UInt8, Images.data(img)[:,:,1])'
end

"""
Load test images
"""
function load_test_images()
  dir = Pkg.dir("ImageRegistration")
  pathA = joinpath(dir, "test", "test_images", "imgA.tif")
  pathB = joinpath(dir, "test", "test_images", "imgB.tif")
  return load_uint8_image(pathA), load_uint8_image(pathB)
end


text:
using Documenter, OMEinsum

makedocs(;
    modules=[OMEinsum],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Parsing" => "parsing.md",
        "Implementations" => "implementation.md",
        "Extending OMEinsum" => "extending.md",
        "DocStrings" => "docstrings.md"
    ],
    repo="https://github.com/under-Peter/OMEinsum.jl/blob/{commit}{path}#L{line}",
    sitename="OMEinsum.jl",
    authors="Andreas Peter",
)

deploydocs(;
    repo="github.com/under-Peter/OMEinsum.jl",
)


text:
push!(LOAD_PATH, Sys.STDLIB)
using SlurmAddAllocatedProcs
using Test
addprocs_slurm_allocated()
using Distributed

SLURM_TASKS_PER_NODE = ENV["SLURM_TASKS_PER_NODE"]
@test nworkers() == parse_slurm_tasks_per_node(SLURM_TASKS_PER_NODE)
hosts = [@fetchfrom w Libc.gethostname() for w in workers()]
@test length(unique(hosts)) == parse(Int, ENV["SLURM_JOB_NUM_NODES"])


text:
# Plotting Options
using Printf
using Plots
using LaTeXStrings
using BSON
pgfplotsx()

includet("./models/dense_rren.jl")

theme(:default)

ms = 10
lw = 2.0
line_alpha = 0.8
c = palette(:default);

plot_args = (xaxis = :log10, seriestype = :scatter)
xaxis = (L"\gamma", (0, 1000), :log, font(20, "Courier"))

function plot_gamma_vs_perf!(data; kwargs...)
    plot!([data.gamma_lower], [mean.([data.val[1]["nrmse"]])]; kwargs...)
end

function plot_gamma_vs_perf!(data, gamma_upper; kwargs...)
    plot!([data.gamma_lower], [mean.([data.val[1]["nrmse"]])]; kwargs..., plot_args...)
    plot!([gamma_upper, gamma_upper], [0, 1.0];  color=kwargs[:color], label=nothing, ls=:dashdot)
end

function format!()
    tfont = font(22, "Courier")
    lfont = font(26, "Courier")

    plot!(;xlabel=L"\underline{\gamma}", xtickfont=tfont, labelfontsize=24, legend=(0.61, 0.99))
    plot!(;ylabel="NRMSE", ytickfont=tfont, legendfontsize=14, framestyle=:box)
end

# ## Plot data F16 Vibration Test
# mf16_ren_4p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_4.0")[2]
# mf16_ren_6p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_6.0")[2]
# mf16_ren_10p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_10.0")[2]
# mf16_ren_20p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_20.0")[2]
# mf16_ren_40p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_40.0")[2]
# mf16_ren_60p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_60.0")[2]
# mf16_ren_100p0 = load_sys_id_test("./results/sys_id/f16/bounded_ren_75_150_100.0")[2]
# mf16_ren_stable = load_sys_id_test("./results/sys_id/f16/stable_ren_75_150")[2]

# # ## Plot data F16 Vibration Test
# mf16_bounded_4p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_4.0")[2]
# mf16_bounded_6p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_6.0")[2]
# mf16_bounded_10p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_10.0")[2]
# mf16_bounded_20p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_20.0")[2]
# mf16_bounded_40p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_40.0")[2]
# mf16_bounded_60p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_60.0")[2]
# mf16_bounded_100p0 = load_sys_id_test("./results/sys_id_cdc2021/f16/bounded_rren_75_150_100.0")[2]
# mf16_stable = load_sys_id_test("./results/sys_id_cdc2021/f16/stable_rren_75_150")[2]

# mf16_ff_4p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_4.0")[2]
# mf16_ff_6p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_6.0")[2]
mf16_ff_10p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_10.0")[2]
mf16_ff_20p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_20.0")[2]
mf16_ff_40p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_40.0")[2]
# mf16_ff_60p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_60.0")[2]
# mf16_ff_100p0 = load_sys_id_test("./results/sys_id/f16/bounded_ffren_75_150_100.0")[2]
mf16_ff_stable = load_sys_id_test("./results/sys_id/f16/stable_ffren_75_150")[2]

# mf16_rrnn_4p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_4.0")[2]
# mf16_rrnn_6p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_6.0")[2]
# mf16_rrnn_10p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_10.0")[2]
# mf16_rrnn_20p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_20.0")[2]
# mf16_rrnn_40p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_40.0")[2]
# mf16_rrnn_60p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_60.0")[2]
# mf16_rrnn_100p0 = load_sys_id_test("./results/sys_id/f16/robust_rnn_lipschitz75_150_100.0")[2]
# mf16_rrnn_stable = load_sys_id_test("./results/sys_id/f16/robust_rnn_stable75_150")[2]

mf16_lstm = load_sys_id_test("./results/sys_id/f16/lstm_170")[2]
mf16_rnn = load_sys_id_test("./results/sys_id/f16/rnn_340")[2]

# gr()
pgfplotsx()
c = palette(:default);
lw = 1.5
plot_font = "Computer Modern"

default(fontfamily=plot_font, labelfontsize=22, legendfontsize=14, tickfontsize=16,
        linewidth=lw, framestyle=:box, label=nothing, grid=true)

p = plot()

plot_gamma_vs_perf!(mf16_ff_10p0, 10.0; label=L"\textrm{R-aREN } \gamma < 10", color=c[2], m=:diamond, msize=ms)
plot_gamma_vs_perf!(mf16_ff_20p0, 20.0; label=L"\textrm{R-aREN }{\gamma <20}", color=c[3], m=:diamond, msize=ms)
plot_gamma_vs_perf!(mf16_ff_40p0, 40.0; label=L"\textrm{R-aREN }{\gamma <40}", color=c[4], m=:diamond, msize=ms)
plot_gamma_vs_perf!(mf16_ff_stable; label=L"\textrm{C-aREN }", color=c[7], m=:diamond, msize=ms)

# Data comes from python impletmentation of Robust RNN
# Data can be found at https://github.com/imanchester/RobustRNN
# branch f16_exp - results/f16
plot!([9.045], [0.5610889]; label=L"\textrm{Robust RNN }{\gamma<10}", color=c[2], m=:dtriangle, msize=ms)
plot!([18.546], [0.475547]; label=L"\textrm{Robust RNN }{\gamma<20}", color=c[3], m=:dtriangle, msize=ms)
plot!([38.846], [0.325848]; label=L"\textrm{Robust RNN }{\gamma<40}", color=c[4], m=:dtriangle, msize=ms)
plot!([148.340], [0.243539]; label=L"\textrm{Robust RNN }{\gamma<\infty}", color=c[7], m=:dtriangle, msize=ms)
# plot!([49.896], [0.335553]; label=L"\textrm{Robust RNN }{\gamma<\infty}", color=c[7], m=:dtriangle, msize=ms)

plot_gamma_vs_perf!(mf16_lstm; label="LSTM", color=c[1], m=:rect, msize=ms)
plot_gamma_vs_perf!(mf16_rnn; label="RNN", color=c[1], m=:circle, msize=ms)

plot!(;xlabel=L"\underline{\gamma}", ylabel="NRMSE", legend=(0.59, 0.98))
ylims!(0.1, 0.60)
minor = [j*10.0^i for i=(0-1):(3+1) for j=2:9]
plot!(; xticks=((10, 100, 1000), ("10", "100", "1000")), minorticks=minor)

savefig(p, "./results/sys_id/f16.pdf")

# plot_gamma_vs_perf!(mf16_lstm; label="LSTM", color=c[1], m=:rect, msize=ms)
# plot_gamma_vs_perf!(mf16_rnn; label="RNN", color=c[1], m=:circle, msize=ms)

# ylims!(0.18, 0.75)
# savefig(p, "./results/sys_id/f16.png")

mean(mf16_ff_stable.val[1]["nrmse"][:])
mean(mf16_bounded_100p0.val[1]["nrmse"][:])

# mf16_bounded_10p0.gamma_lower
# mf16_bounded_20p0.gamma_lower
# mf16_bounded_40p0.gamma_lower
# mf16_bounded_60p0.gamma_lower
# mf16_bounded_100p0.gamma_lower

# mf16_ff_10p0.gamma_lower
# mf16_ff_20p0.gamma_lower
# mf16_ff_40p0.gamma_lower
# mf16_ff_60p0.gamma_lower
# mf16_ff_100p0.gamma_lower
# mf16_ff_stable.gamma_lower

# # Example outputs for f16 dataset & 91.0
# function make_plots()
#     plot_array = Any[]
#     range = 3500:3750
#     for output = 1:3
#         for batch = 1:3
#             y = extract(mf16_bounded_20p0.val[1]["y"], output, batch)[range]
#             yest1 = extract(mf16_bounded_4p0.val[1]["yest"], output, batch)[range]
#             # yest2 = extract(mf16_bounded_10p0.val[1]["yest"], output, batch)[range]
#             yest3 = extract(mf16_bounded_20p0.val[1]["yest"], output, batch)[range]
#             yest6 = extract(mf16_bounded_100p0.val[1]["yest"], output, batch)[range]
#             yest7 = extract(mf16_stable.val[1]["yest"], output, batch)[range]
#             yest8 = extract(mf16_lstm.val[1]["yest"], output, batch)[range]

#             push!(plot_array, plot())
#             plot!(y; c=:black, lw=3)
#             plot!(yest1; c=c[1],lw=1.5)
#             # plot!(yest2; c=c[2],lw=1.5)
#             plot!(yest3; c=c[2], lw=1.5)
#             plot!(yest7; c=c[3], lw=1.5)
#             plot!(yest8; c=c[4], lw=1.5,  legend=nothing)
#         end
#     end
#     return plot(plot_array..., layouts=(3, 3))
# end
# p = make_plots()
# ylims!(-4, 4)
# savefig(p, "./results/sys_id/f16_outputs.png")

# Plot training loss
# c1 = cgrad([:orange, :blue], LinRange(0, 1, 6))
c1 = palette(:thermal, 5)

p = plot(mf16_ff_4p0.tloss, label=L"\textrm{R-aREN }{\gamma < 4}", c=c1[4])
plot!(mf16_ff_10p0.tloss, label=L"\textrm{R-aREN }{\gamma < 10}", c=c1[3])
plot!(mf16_ff_40p0.tloss, label=L"\textrm{R-aREN }{\gamma < 40}", c=c1[2])
plot!(mf16_ff_stable.tloss, label=L"\textrm{C-aREN }", c=c1[1])

plot!(mf16_lstm.tloss, label="LSTM", lw=1.5, c=c[1])
plot!(mf16_rnn.tloss, label="RNN", lw=1.5, c=c[3])
plot!(;xlabel="Epochs", ylabel="Loss" ,yaxis=:log10, legend = (0.66, 0.98))
ylims!(0.5, 11.1)

plot!(;yticks=((1, 10), ("1", "10")))
savefig(p, "./results/sys_id/f16_training.pdf")

# Plot Wiener Hammerstein Results
wh_ff_1p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_1.0")[2]
wh_ff_1p5 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_1.5")[2]
wh_ff_2p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_2.0")[2]
wh_ff_2p5 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_2.5")[2]
wh_ff_3p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_3.0")[2]
wh_ff_3p5 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_3.5")[2]
wh_ff_5p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_5.0")[2]
wh_ff_8p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_8.0")[2]
wh_ff_15p0 = load_sys_id_test("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_15.0")[2]
wh_ff_stable = load_sys_id_test("./results/sys_id/wiener_hammerstein/stable_ffren_40_100")[2]
wh_lstm = load_sys_id_test("./results/sys_id/wiener_hammerstein/lstm_100")[2]
wh_rnn = load_sys_id_test("./results/sys_id/wiener_hammerstein/rnn_200")[2]

# gr()
pgfplotsx()

p = plot()
plot_gamma_vs_perf!(wh_ff_1p5, 1.5; label=L"\textrm{R-aREN } {\gamma < 1.5}", color=c[2], m=:diamond, msize=ms)
plot_gamma_vs_perf!(wh_ff_2p5, 2.5; label=L"\textrm{R-aREN } {\gamma < 2.5}", color=c[3], m=:diamond, msize=ms)
plot_gamma_vs_perf!(wh_ff_3p5, 3.5; label=L"\textrm{R-aREN } {\gamma < 3.5}", color=c[4], m=:diamond, msize=ms)
plot_gamma_vs_perf!(wh_ff_stable; label=L"\textrm{C-aREN }", color=c[7], m=:diamond, msize=ms)
plot_gamma_vs_perf!(wh_lstm; label="LSTM", color=c[1], m=:utriangle, msize=ms)
plot_gamma_vs_perf!(wh_rnn; label="RNN", color=c[1], m=:circle, msize=ms)

ylims!(0.25, 0.5)
xlims!(1.3, 1000)
plot!(;xlabel=L"\underline{\gamma}", ylabel="NRMSE", legend=(0.64, 0.43))
plot!(; xticks=((10, 100, 1000), ("10", "100", "1000")))
savefig(p, "./results/sys_id/wiener_hammerstein.pdf")

# Plot adversarial perterbations for f16 example
function calc_pert(data; pert_size=0.05)
    model = data.model
    u1 = data.u1
    u2 = data.u2

    du = u2 - u1
    Δu = pert_size * du / norm(du)

    x0 = init_state(model, 1)
    y1 = model(x0, u1)[1]
    y2 = model(x0, u1 + Δu)[1]

    Δy = y2 - y1
    println("Sensitivity: ", norm(Δy) / norm(Δu))

    return Δy, Δu
end

pgfplotsx()
c = palette(:default);
lw = 1.5
plot_font = "Computer Modern"

default(fontfamily=plot_font, labelfontsize=22, legendfontsize=14, tickfontsize=16,
        linewidth=lw, framestyle=:box, label=nothing, grid=true)

ffren_10_dat = BSON.load("./results/sys_id/f16_adversarial_perturbations/ffren_10.bson")[:data]
ffren_40_dat = BSON.load("./results/sys_id/f16_adversarial_perturbations/ffren_40.bson")[:data]
rnn_dat = BSON.load("./results/sys_id/f16_adversarial_perturbations/rnn.bson")[:data]
lstm_dat = BSON.load("./results/sys_id/f16_adversarial_perturbations/lstm.bson")[:data]

rnn_dy, rnn_du = calc_pert(rnn_dat)
lstm_dy, lstm_du = calc_pert(lstm_dat)
ffren_40_dy, ffren_40_du = calc_pert(ffren_40_dat)
ffren_10_dy, ffren_10_du = calc_pert(ffren_10_dat)

k = 1
p=plot(;ylabel=L"\Delta y_1")
plot!(extract(rnn_dy, k, 1); label=L"\textrm{RNN}", c=c[4])
plot!(extract(lstm_dy, k, 1); label=L"\textrm{LSTM}", c=c[3])
plot!(extract(ffren_40_dy, k, 1); label=L"\textrm{R-aREN} \gamma<40", c=c[2])
plot!(extract(ffren_10_dy, k, 1); label=L"\textrm{R-aREN} \gamma<10", c=c[1])
plot!(;xlabel="Time Steps", legend=(0.02, 0.98))
savefig(p, "./results/sys_id/adv_pert.pdf")

# Zoomed in version of the above plot
plot_range = 1000:1500
p=plot(;ylabel=L"\Delta y_1")
plot!(plot_range, extract(rnn_dy[plot_range], k, 1); label=L"\textrm{RNN}", c=c[4])
plot!(plot_range, extract(lstm_dy[plot_range], k, 1); label=L"\textrm{LSTM}", c=c[3])
plot!(plot_range, extract(ffren_40_dy[plot_range], k, 1); label=L"\textrm{R-aREN} \gamma<40", c=c[2])
plot!(plot_range, extract(ffren_10_dy[plot_range], k, 1); label=L"\textrm{R-aREN} \gamma<10", c=c[1])
plot!(;xlabel="Time Steps", legend=(0.01, 0.98))
savefig(p, "./results/sys_id/adv_pert_zoomed.pdf")

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
"""
# HDFSException
Thrown by HDFS APIs.
"""
struct HDFSException <: Exception
    message::AbstractString
end

"""
HDFSProtocol: Hadoop RPC client for namenode protocol.
"""
const HDFSProtocol = HadoopRpcProtocol{ClientNamenodeProtocolBlockingStub}

for fn in (:getListing, :getFileInfo, :getBlockLocations, :getServerDefaults, :getFsStats, :setReplication, :getContentSummary, :mkdirs, :rename2, :rename, :delete, :create, :complete, :addBlock, :renewLease, :setTimes)
    @eval begin
        (hadoop.hdfs.$fn)(p::HDFSProtocol, inp) = (hadoop.hdfs.$fn)(p.stub, p.controller, inp)
    end
end

"""
# HDFSClient
A client to the namenode in a HDFS cluster and holds a connection to it.
It also stores the folder context for using relative paths in APIs that use the client.
"""
mutable struct HDFSClient
    nn_conn::HDFSProtocol
    wd::AbstractString
    server_defaults::Union{Nothing,FsServerDefaultsProto}

    function HDFSClient(host::AbstractString, port::Integer, ugi::UserGroupInformation=UserGroupInformation())
        nn_conn = HDFSProtocol(host, port, ugi)
        new(nn_conn, "/", nothing)
    end
end

function show(io::IO, client::HDFSClient)
    show(io, client.nn_conn)
    println(io, "    pwd: $(client.wd)")
    nothing
end

set_debug(client::HDFSClient, debug::Bool) = set_debug(client.nn_conn, debug)

abspath(client::HDFSClient, path::AbstractString) = abspath(client.wd, path)
function cd(client::HDFSClient, path::AbstractString)
    wd = abspath(client, path)
    !isdir(client, wd) && throw(HDFSException("Directory not found: $wd"))
    client.wd = wd
end
pwd(client::HDFSClient) = client.wd

"""
# HDFSFile
Holds a HDFSClient and path pointing to a file on HDFS
"""
struct HDFSFile
    client::HDFSClient
    path::AbstractString
end

function HDFSFile(uristr::AbstractString; ugi::Union{UserGroupInformation,Nothing}=nothing, proxy::Bool=false)
    uri = URI(uristr)
    (uri.scheme == "hdfs") || throw(HDFSException("Not a HDFS URI: $uristr"))
    if ugi === nothing
        ugi = isempty(uri.userinfo) ? UserGroupInformation() : UserGroupInformation(uri.userinfo)
    end
    client = HDFSClient(uri.host, uri.port, ugi)
    HDFSFile(client, uri.path)
end

function show(io::IO, file::HDFSFile)
    println(io, "HDFSFile: ", string(convert(URI, file)))
    nothing
end

function convert(::Type{URI}, file::HDFSFile)
    client = file.client
    ch = client.nn_conn.channel
    user = username(ch.ugi)
    user_spec = isempty(user) ? user : "$(user)@"
    URI("hdfs://$(user_spec)$(ch.host):$(ch.port)$(abspath(client, file.path))")
end

"""
# HDFSFileInfo
Stat structure for file/foledr in HDFS
"""
mutable struct HDFSFileInfo
    kind::Int32
    name::AbstractString
    size::UInt64
    replications::UInt32
    block_sz::UInt64
    owner::AbstractString
    grp::AbstractString
    permissions::UInt32
    last_mod::UInt64
    last_access::UInt64

    HDFSFileInfo(fs::HdfsFileStatusProto) = new(fs.fileType, String(copy(fs.path)),
                    fs.length, fs.block_replication, fs.blocksize,
                    fs.owner, fs.group, fs.permission.perm,
                    fs.modification_time, fs.access_time)
end

function show(io::IO, st::HDFSFileInfo)
    println(io, "HDFSFileInfo: $(st.name)")
    if isdir(st)
        println(io, "    type: folder")
    elseif islink(st)
        println(io, "    type: link")
    else
        println(io, "    type: file")
        println(io, "    size: $(st.size)")
        println(io, "    block_sz: $(st.block_sz)")
    end
    println(io, "    owner: $(st.owner)")
    println(io, "    group: $(st.grp)")
    nothing
end

function _as_dict(obj, d=Dict{Symbol,Any}())
    for name in propertynames(obj)
        if hasproperty(obj, name)
            d[name] = getproperty(obj, name)
        end
    end
    d
end

function _walkdir(client::HDFSClient, path::AbstractString, process_entry::Function)
    path = abspath(client, path)
    cont = true
    start_after = UInt8[]
    while cont
        inp = GetListingRequestProto(src=path, startAfter=start_after, needLocation=false)
        resp = getListing(client.nn_conn, inp)

        if hasproperty(resp, :dirList)
            dir_list = resp.dirList
            if hasproperty(dir_list, :partialListing)
                partial_listing = dir_list.partialListing
                for filestatus in partial_listing
                    cont = process_entry(filestatus)
                    (cont == false) && break
                end
                (cont == false) && break
            end
            if hasproperty(dir_list, :remainingEntries) && (dir_list.remainingEntries > 0)
                start_after = dir_list.partialListing[end].path
            else
                cont = false
            end
        else
            cont = false
        end
    end
    nothing
end

function _get_file_info(client::HDFSClient, path::AbstractString)
    path = abspath(client, path)
    inp = GetFileInfoRequestProto(src=path)
    resp = getFileInfo(client.nn_conn, inp)
    hasproperty(resp, :fs) ? resp.fs : nothing
end

function _get_block_locations(client::HDFSClient, path::AbstractString, offset::UInt64=zero(UInt64), length::UInt64=zero(UInt64))
    path = abspath(client, path)
    if length == 0
        length = UInt64(1024)
    end
    inp = GetBlockLocationsRequestProto(src=path, offset=offset, length=length)
    resp = getBlockLocations(client.nn_conn, inp)
    hasproperty(resp, :locations) ? resp.locations : nothing
end

#
# Server defaults
function _get_server_defaults(client::HDFSClient)
    if client.server_defaults === nothing
        inp = GetServerDefaultsRequestProto()
        resp = getServerDefaults(client.nn_conn, inp)
        client.server_defaults = resp.serverDefaults
    end
    client.server_defaults
end

hdfs_server_defaults(client::HDFSClient) = _as_dict(_get_server_defaults(client))
hdfs_default_block_size(client::HDFSClient) = _get_server_defaults(client).blockSize
hdfs_default_replication(client::HDFSClient) = _get_server_defaults(client).replication

#
# File System Status
_get_fs_status(client::HDFSClient) = getFsStats(client.nn_conn, GetFsStatusRequestProto())

hdfs_status(client::HDFSClient) = _as_dict(_get_fs_status(client))

hdfs_capacity(client::HDFSClient) = _get_fs_status(client).capacity
hdfs_capacity_used(client::HDFSClient) = _get_fs_status(client).used
hdfs_capacity_remaining(client::HDFSClient) = _get_fs_status(client).remaining

#
# File Status
stat(file::HDFSFile) = stat(file.client, file.path)
function stat(client::HDFSClient, path::AbstractString)
    fileinfo = _get_file_info(client, path)
    (fileinfo === nothing) && throw(HDFSException("Path not found $path"))
    hdfs_file_info = HDFSFileInfo(fileinfo)
    hdfs_file_info.name = path
    hdfs_file_info
end

exists(file::HDFSFile) = exists(file.client, file.path)
exists(client::HDFSClient, path::AbstractString) = (_get_file_info(client, path) !== nothing)

isdir(file::HDFSFile) = isdir(file.client, file.path)
isdir(client::HDFSClient, path::AbstractString) = isdir(stat(client, path))
isdir(fileinfo::HDFSFileInfo) = (fileinfo.kind == HdfsFileStatusProto_FileType.IS_DIR)

isfile(file::HDFSFile) = isfile(file.client, file.path)
isfile(client::HDFSClient, path::AbstractString) = isfile(stat(client, path))
isfile(fileinfo::HDFSFileInfo) = (fileinfo.kind == HdfsFileStatusProto_FileType.IS_FILE)

islink(file::HDFSFile) = islink(file.client, file.path)
islink(client::HDFSClient, path::AbstractString) = islink(stat(client, path))
islink(fileinfo::HDFSFileInfo) = (fileinfo.kind == HdfsFileStatusProto_FileType.IS_SYMLINK)

filesize(file::HDFSFile) = filesize(file.client, file.path)
filesize(client::HDFSClient, path::AbstractString) = filesize(stat(client, path))
filesize(fileinfo::HDFSFileInfo) = fileinfo.size

filemode(file::HDFSFile) = filemode(file.client, file.path)
filemode(client::HDFSClient, path::AbstractString) = filemode(stat(client, path))
filemode(fileinfo::HDFSFileInfo) = fileinfo.permissions

mtime(file::HDFSFile) = mtime(file.client, file.path)
mtime(client::HDFSClient, path::AbstractString) = mtime(stat(client, path))
mtime(fileinfo::HDFSFileInfo) = fileinfo.last_mod

atime(file::HDFSFile) = atime(file.client, file.path)
atime(client::HDFSClient, path::AbstractString) = atime(stat(client, path))
atime(fileinfo::HDFSFileInfo) = fileinfo.last_access

hdfs_blocks(file::HDFSFile, offset::Integer=0, length::Integer=typemax(Int)) = hdfs_blocks(file.client, file.path, offset, length)
function hdfs_blocks(client::HDFSClient, path::AbstractString, offset::Integer=0, length::Integer=typemax(Int))
    blocks = Tuple{UInt64,Array}[]
    locations = _get_block_locations(client, path, UInt64(offset), UInt64(length))
    (locations === nothing) && (return blocks)
    for block in locations.blocks
        block.corrupt && throw(HDFSException("Corrupt block found at offset $(block.offset)"))
        node_ips = AbstractString[]
        for loc in block.locs
            (loc.adminState == DatanodeInfoProto_AdminState.NORMAL) || continue
            node_id = loc.id
            push!(node_ips, node_id.ipAddr)
        end
        push!(blocks, (block.offset, node_ips))
    end
    blocks
end

hdfs_set_replication(file::HDFSFile, replication::Integer) = hdfs_set_replication(file.client, file.path, replication)
function hdfs_set_replication(client::HDFSClient, path::AbstractString, replication::Integer)
    path = abspath(client, path)
    inp = SetReplicationRequestProto(src=path, replication=replication)

    resp = setReplication(client.nn_conn, inp)
    resp.result
end

#
# Disk Usage
function _get_content_summary(client::HDFSClient, path::AbstractString)
    path = abspath(client, path)
    inp = GetContentSummaryRequestProto(path=path)

    resp = getContentSummary(client.nn_conn, inp)
    resp.summary
end

du(file::HDFSFile, details::Union{Nothing,Dict{Symbol,Any}}=nothing) = du(file.client, file.path, details)
function du(client::HDFSClient, path::AbstractString=".", details::Union{Nothing,Dict{Symbol,Any}}=nothing)
    summary = _get_content_summary(client, path)
    (details === nothing) || _as_dict(summary, details)
    summary.length
end

#
# File create, delete, list
readdir(file::HDFSFile, limit::Int=typemax(Int)) = readdir(file.client, file.path, limit)
function readdir(client::HDFSClient, path::AbstractString=".", limit::Int=typemax(Int))
    result = String[]
    _walkdir(client, path, (filestatus)->begin
                push!(result, String(copy(filestatus.path)))
                (length(result) < limit)
            end)
    result
end

const DEFAULT_FOLDER_MODE = UInt32(0o755)
const DEFAULT_FILE_MODE = UInt32(0o644)
mkdir(file::HDFSFile, mode::UInt32) = mkdir(file, false, mode)
mkdir(file::HDFSFile, create_parents::Bool=false, mode::UInt32=DEFAULT_FOLDER_MODE) = mkdir(file.client, file.path, create_parents, mode)
function mkdir(client::HDFSClient, path::AbstractString, create_parents::Bool=false, mode::UInt32=DEFAULT_FOLDER_MODE)
    path = abspath(client, path)
    perms = FsPermissionProto(perm=mode)
    inp = MkdirsRequestProto(src=path, createParent=create_parents, masked=perms)

    resp = mkdirs(client.nn_conn, inp)
    resp.result
end

mkpath(file::HDFSFile, mode::UInt32=DEFAULT_FOLDER_MODE) = mkdir(file, true, mode)
dirname(file::HDFSFile) = HDFSFile(file.client, dirname(file.path))

mv(src::HDFSFile, dst::AbstractString, overwrite::Bool) = mv(src.client, src.path, dst, overwrite)
function mv(client::HDFSClient, src::AbstractString, dst::AbstractString, overwrite::Bool)
    src = abspath(client, src)
    dst = abspath(client, dst)
    inp = Rename2RequestProto(src=src, dst=dst, overwriteDest=overwrite)

    rename2(client.nn_conn, inp)
    true
end

mv(src::HDFSFile, dst::AbstractString) = mv(src.client, src.path, dst)
function mv(client::HDFSClient, src::AbstractString, dst::AbstractString)
    src = abspath(client, src)
    dst = abspath(client, dst)
    inp = RenameRequestProto(src=src, dst=dst)

    resp = rename(client.nn_conn, inp)
    resp.result
end

rm(file::HDFSFile, recursive::Bool=false) = rm(file.client, file.path, recursive)
function rm(client::HDFSClient, path::AbstractString, recursive::Bool=false)
    path = abspath(client, path)
    inp = DeleteRequestProto(src=path, recursive=recursive)

    resp = delete(client.nn_conn, inp)
    resp.result
end

joinpath(hdfsfile::HDFSFile, parts...) = HDFSFile(hdfsfile.client, joinpath(hdfsfile.path, parts...))

function _create_file(client::HDFSClient, path::AbstractString, overwrite::Bool=false, replication::UInt32=zero(UInt32), blocksz::UInt64=zero(UInt64), mode::UInt32=DEFAULT_FILE_MODE, docomplete::Bool=true)
    path = abspath(client, path)

    if (blocksz == 0) || (replication == 0)
        defaults = _get_server_defaults(client)
        (blocksz == 0) && (blocksz = defaults.blockSize)
        (replication == 0) && (replication = defaults.replication)
    end

    perms = FsPermissionProto(perm=mode)
    createFlag = UInt32(overwrite ? CreateFlagProto.OVERWRITE : CreateFlagProto.CREATE)

    inp = CreateRequestProto(src=path, masked=perms, clientName=ELLY_CLIENTNAME, createFlag=createFlag, createParent=false, replication=replication, blockSize=blocksz)

    resp = create(client.nn_conn, inp)
    hasproperty(resp, :fs) || (return nothing)

    if docomplete
        _complete_file(client, path) || (return nothing)
    end

    return resp.fs
end

function _complete_file(client::HDFSClient, path::AbstractString, last::Union{Nothing,ExtendedBlockProto}=nothing)
    path = abspath(client, path)

    endinp = CompleteRequestProto(src=path, clientName=ELLY_CLIENTNAME)
    if last !== nothing
        setproperty!(endinp, :last, last)
        @debug("setting last block as", last)
    end

    endresp = complete(client.nn_conn, endinp)
    endresp.result
end

function _add_block(::Type{T}, client::HDFSClient, path::AbstractString, previous::Union{Nothing,T}=nothing) where T<:LocatedBlockProto
    (previous === nothing) && (return _add_block(ExtendedBlockProto, client, path))
    @debug("adding next block to", previous.b)
    _add_block(ExtendedBlockProto, client, path, previous.b)
end
function _add_block(::Type{T}, client::HDFSClient, path::AbstractString, previous::Union{Nothing,T}=nothing) where T<:ExtendedBlockProto
    path = abspath(client, path)

    inp = AddBlockRequestProto(src=path, clientName=ELLY_CLIENTNAME)
    (previous === nothing) || setproperty!(inp, :previous, previous)

    resp = addBlock(client.nn_conn, inp)
    return resp.block
end

"""
Applications that write infrequently/slowly must call renewlease periodically to prevent
the namenode from assuming the client from having abandoned the file or some other client
from recovering the lease.
"""
function hdfs_renewlease(client::HDFSClient)
    inp = RenewLeaseRequestProto(clientName=ELLY_CLIENTNAME)
    renewLease(client.nn_conn, inp)
    nothing
end

touch(file::HDFSFile, replication::UInt32=zero(UInt32), blocksz::UInt64=zero(UInt64), mode::UInt32=DEFAULT_FILE_MODE) = touch(file.client, file.path, replication, blocksz, mode)
function touch(client::HDFSClient, path::AbstractString, replication::UInt32=zero(UInt32), blocksz::UInt64=zero(UInt64), mode::UInt32=DEFAULT_FILE_MODE)
    if exists(client, path)
        path = abspath(client, path)
        t = UInt64(datetime2unix(now(UTC))*1000)
        inp = SetTimesRequestProto(src=path, mtime=t, atime=t)

        setTimes(client.nn_conn, inp)
    else
        fs = _create_file(client, path, false, replication, blocksz, mode)
        (fs === nothing) && throw(HDFSException("Could not create file $path"))
    end
    nothing
end


text:
import Base.Iterators: flatten
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph, Edge

graph_attrs = Attributes(:rankdir=>"LR")
node_attrs  = Attributes(:shape=>"plain", :style=>"filled", :color=>"white")
edge_attrs  = Attributes(:splines=>"splines")

function edgify(δ, transition, reverse::Bool)
    return [Edge(reverse ? ["T_$transition", "S_$k"] :
                           ["S_$k", "T_$transition"],
                 Attributes(:label=>"$(δ[k])", :labelfontsize=>"6"))
             for k in collect(keys(δ)) if δ[k] != 0]
end

"""
    Graph(model::Model)

convert a Model into a GraphViz Graph. Transition are green boxes and states are blue circles. Arrows go from the input states to the output states for each transition.
"""
function Graph(model::Model)
    ks = collect(keys(model.Δ))
    statenodes = [Node(string("S_$s"), Attributes(:shape=>"circle", :color=>"#6C9AC3")) for s in model.S]
    transnodes = [Node(string("T_$k"), Attributes(:shape=>"square", :color=>"#E28F41")) for k in ks]

    stmts = vcat(statenodes, transnodes)
    edges = map(ks) do k
      vcat(edgify(first(model.Δ[k]), k, false),
           edgify(last(model.Δ[k]), k, true))
    end |> flatten |> collect
    stmts = vcat(stmts, edges)
    g = Graphviz.Digraph("G", stmts; graph_attrs=graph_attrs, node_attrs=node_attrs, edge_attrs=edge_attrs)
    return g
end


text:
using Test
import DPFEHM
import Zygote

doplot = false
if doplot == true
	import PyPlot
end

function thiem(R, r, T, Q)
	return Q * log(R / r) / (2 * pi * T)
end
#utitlity for checking gradients
function checkgradientquickly(f, x0, gradf, n; delta::Float64=1e-8, kwargs...)
	indicestocheck = sort(collect(1:length(x0)), by=i->abs(gradf[i]), rev=true)[1:n]
	#indicestocheck = [indicestocheck; rand(1:length(x0), n)]
	f0 = f(x0)
	for i in indicestocheck
		x = copy(x0)
		x[i] += delta
		fval = f(x)
		grad_f_i = (fval - f0) / delta
		@test isapprox(gradf[i], grad_f_i; kwargs...)
	end
end

#test Thiem solution against groundwater model
steadyhead = 1e3
sidelength = 50.0
thickness = 10.0
R = sidelength
n = 101
ns = [n, n]
coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid2d([-sidelength, -sidelength], [sidelength, sidelength], ns, thickness)
k = 1e-5
Ks = fill(k, length(neighbors))
Q = 1e-3
Qs = zeros(size(coords, 2))
Qs[ns[2] * (div(ns[1] + 1, 2) - 1) + div(ns[2] + 1, 2)] = -Q#put a fluid source in the middle
dirichletnodes = Int[]
dirichleths = zeros(size(coords, 2))
for i = 1:size(coords, 2)
	if sqrt(coords[1, i] ^ 2 + coords[2, i] ^ 2) >= R
		push!(dirichletnodes, i)
		dirichleths[i] = steadyhead
	end
end
print("steady state groundwater forward")
@time h_gw = DPFEHM.groundwater_steadystate(Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs; reltol=1e-12)
r0 = 0.1
goodnodes = collect(filter(i->coords[2, i] == 0 && coords[1, i] > r0, 1:size(coords, 2)))
rs = coords[1, goodnodes]
T = thickness * k
thiem_drawdowns = thiem.(R, rs, T, Q)
gw_drawdowns = -h_gw[goodnodes] .+ steadyhead
if doplot
	fig, ax = PyPlot.subplots()
	ax.plot(rs, thiem_drawdowns, "r.", ms=20, label="Thiem")
	ax.plot(rs, gw_drawdowns, "k", linewidth=3, label="DPFEHM groundwater")
	ax.set_xlabel("x [m]")
	ax.set_ylabel("drawdown [m]")
	ax.legend()
	display(fig)
	println()
	PyPlot.close(fig)
end
@test isapprox(thiem_drawdowns, gw_drawdowns; rtol=1e-1)
g_gw(Ks) = DPFEHM.groundwater_steadystate(Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs; reltol=1e-12)[goodnodes[round(Int, 0.25 * end)], end]
print("groundwater gradient")
@time grad_gw_zygote = Zygote.gradient(g_gw, Ks)[1]
checkgradientquickly(g_gw, Ks, grad_gw_zygote, 3; delta=1e-8, rtol=1e-1)


text:
# Implementation of
#  "Table-driven Implementation of the Logarithm Function in IEEE Floating-point Arithmetic"
#  Tang, Ping-Tak Peter
#  ACM Trans. Math. Softw. (1990), 16(4):378--400
#  http://dx.doi.org/10.1145/98267.98294

# Other implementations:
#  Apple OS X libm (APSL v2)
#   http://www.opensource.apple.com/source/Libm/Libm-2026/Source/Intel/xmm_log.c
#  libacml_mv (GPL v2)
#   http://svn.open64.net/svnroot/open64/trunk/osprey/libacml_mv/src/log1p.c
#  libclc (BSD): modified to use different table size
#   https://github.com/llvm-mirror/libclc/blob/master/generic/lib/math/log1p.cl

# At the moment we do not correctly handle floating point flags (inexact, div-by-zero, etc).

import Base.unsafe_trunc

const t_log_64 = Array((Float64,Float64),129)
N=39 # (can be up to N=42, which appears to be what Apple's libm uses).
sN = 2.0^N
isN = 1.0/sN
s7 = 2.0^7
is7 = 1.0/s7

for j=0:128
    l_big = Base.log(big(1.0+j*is7))
    l_hi = isN*float64(round(sN*l_big))
    l_lo = float64(l_big-l_hi)
    t_log_64[j+1] = (l_hi,l_lo)
    # c_invF[j+1] = 1.0+is7*j
end

const t_log_32 = Array((Float32,Float32),129)
N=16
sN = 2f0^N
isN = 1f0/sN
for j=0:128
    l_big = Base.log(big(1.0+j*is7))
    l_hi = isN*float32(round(sN*l_big))
    l_lo = float32(l_big-l_hi)
    t_log_32[j+1] = (l_hi,l_lo)
end

const tt_log_32 = Array(Float64,129)
N=16
sN = 2f0^N
isN = 1f0/sN
for j=0:128
    tt_log_32[j+1] = float64(Base.log(big(1.0+j*is7)))
end

# determine if hardware FMA is available
# should probably check with LLVM, see #9855.
const FMA_NATIVE = muladd(nextfloat(1.0),nextfloat(1.0),-nextfloat(1.0,2)) == -4.930380657631324e-32

# truncate lower order bits (up to 26)
# ideally, this should be able to use ANDPD instructions, see #9868.
@inline function truncbits(x::Float64)
    reinterpret(Float64, reinterpret(UInt64,x) & 0xffff_ffff_f800_0000)
end

# Procedure 1
@inline function log_proc1(y::Float64,mf::Float64,F::Float64,f::Float64,jp::Int)
    ## Steps 1 and 2
    @inbounds hi,lo = libm.t_log_64[jp]
    l_hi = mf* 0x1.62e42fefa4p-1 + hi
    l_lo = mf*-0x1.8432a1b0e2634p-43 + lo

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # u = f/F
    # q = u*u*@horner(u,
    #                 -0x1.0_0000_0000_0001p-1,
    #                 +0x1.5_5555_5550_9ba5p-2,
    #                 -0x1.f_ffff_ffeb_6526p-3,
    #                 +0x1.9_99b4_dfed_6fe4p-3,
    #                 -0x1.5_5576_6647_2e04p-3)

    ## Step 3' (alternative)
    u = (2.0f)/(y+F)
    v = u*u
    q = u*v*@horner(v,
                    0x1.5_5555_5555_0286p-4,
                    0x1.9_99a0_bc71_2416p-7)

    ## Step 4
    l_hi + (u + (q + l_lo))
end

# Procedure 2
@inline function log_proc2(f::Float64)
    ## Step 1
    g = 1.0/(2.0+f)
    u = 2.0*f*g
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    0x1.5_5555_5555_54e6p-4,
                    0x1.9_9999_99ba_c6d4p-7,
                    0x1.2_4923_07f1_519fp-9,
                    0x1.c_8034_c85d_fff0p-12)

    ## Step 3
    # based on:
    #   2(f-u) = 2(f(2+f)-2f)/(2+f) = 2f^2/(2+f) = fu
    #   2(f-u1-u2) - f*(u1+u2) = 0
    #   2(f-u1) - f*u1 = (2+f)u2
    #   u2 = (2(f-u1) - f*u1)/(2+f)
    if FMA_NATIVE
        return u + fma(fma(-u,f,2(f-u)), g, q)
    else
        u1 = truncbits(u) # round to 24 bits
        f1 = truncbits(f)
        f2 = f-f1
        u2 = ((2.0*(f-u1)-u1*f1)-u1*f2)*g
        ## Step 4
        return u1 + (u2 + q)
    end
end

@inline function log_proc1(y::Float32,mf::Float32,F::Float32,f::Float32,jp::Int)
    ## Steps 1 and 2
    @inbounds hi = tt_log_32[jp]
    l = mf*0.6931471805599453 + hi

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # q = u*u*@horner(u,
    #                 Float32(-0x1.00006p-1),
    #                 Float32(0x1.55546cp-2))

    ## Step 3' (alternative)
    u = (2f0f)/(y+F)
    v = u*u
    q = u*v*Float32(0x1.555584p-4)

    ## Step 4
    float32(l + (u + q))
end

@inline function log_proc2(f::Float32)
    ## Step 1
    # compute in higher precision
    u64 = Float64(2f0*f)/(2.0+f)
    u = Float32(u64)
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    Float32(0x1.555552p-4),
                    Float32(0x1.9a012ap-7))

    ## Step 3: not required

    ## Step 4
    Float32(u64 + q)
end

function log_tang(x::Float64)
    if x > 0.0
        x == Inf && return x

        # Step 2
        if 0x1.e_0fab_fbc7_02a3p-1 < x < 0x1.1_082b_577d_34eep0
            f = x-1.0
            return log_proc2(f)
        end

        # Step 3
        xu = reinterpret(UInt64,x)
        m = int(xu >> 52) & 0x07ff
        if m == 0 # x is subnormal
            x *= 0x1p54 # normalise significand
            xu = reinterpret(UInt64,x)
            m = int(xu >> 52) & 0x07ff - 54
        end
        m -= 1023
        y = reinterpret(Float64,(xu & 0x000f_ffff_ffff_ffff) | 0x3ff0_0000_0000_0000)

        mf = Float64(m)
        F = (y + 0x1p45) - 0x1p45 # 0x1p-7*round(0x1p7*y)
        f = y-F
        jp = unsafe_trunc(Int,0x1p7*F)-127

        return log_proc1(y,mf,F,f,jp)
    elseif x == 0.0
        -Inf
    elseif isnan(x)
        NaN
    else
        throw(DomainError())
    end
end

function log_tang(x::Float32)
    if x > 0f0
        x == Inf32 && return x

        # Step 2
        if Float32(0x1.e0fabep-1) < x < Float32(0x1.1082b6p+0)
            f = x-1f0
            return log_proc2(f)
        end

        # Step 3
        xu = reinterpret(UInt32,x)
        m = int(xu >> 23) & 0x00ff
        if m == 0 # x is subnormal
            x *= Float32(0x1p25) # normalise significand
            xu = reinterpret(UInt32,x)
            m = int(xu >> 23) & 0x00ff - 25
        end
        m -= 127
        y = reinterpret(Float32,(xu & 0x007f_ffff) | 0x3f80_0000)

        mf = Float32(m)
        F = (y + Float32(0x1p16)) - Float32(0x1p16) # 0x1p-7*round(0x1p7*y)
        f = y-F
        jp = unsafe_trunc(Int,Float32(0x1p7)*F)-127

        log_proc1(y,mf,F,f,jp)
    elseif x == 0f0
        -Inf32
    elseif isnan(x)
        NaN32
    else
        throw(DomainError())
    end
end

function log1p_tang(x::Float64)
    if x > -1.0
        x == Inf && return x
        if -0x1p-53 < x < 0x1p-53
            return x

        # Step 2
        elseif -0x1.f_0540_438f_d5c4p-5 < x < 0x1.0_82b5_77d3_4ed8p-4
            return log_proc2(x)
        end

        # Step 3
        z = 1.0 + x
        zu = reinterpret(UInt64,z)
        s = reinterpret(Float64,0x7fe0_0000_0000_0000 - (zu & 0xfff0_0000_0000_0000)) # 2^-m
        m = int(zu >> 52) & 0x07ff - 1023 # z cannot be subnormal
        c = m > 0 ? 1.0-(z-x) : x-(z-1.0) # 1+x = z+c exactly
        y = reinterpret(Float64,(zu & 0x000f_ffff_ffff_ffff) | 0x3ff0_0000_0000_0000)

        mf = Float64(m)
        F = (y + 0x1p45) - 0x1p45 # 0x1p-7*round(0x1p7*y)
        f = (y - F) + c*s #2^m(F+f) = 1+x = z+c
        jp = unsafe_trunc(Int,0x1p7*F)-127

        log_proc1(y,mf,F,f,jp)
    elseif x == -1.0
        -Inf
    elseif isnan(x)
        NaN
    else
        throw(DomainError())
    end
end

function log1p_tang(x::Float32)
    if x > -1f0
        x == Inf32 && return x
        if -0x1p-53 < x < 0x1p-53
            return x # Inexact

        # Step 2
        elseif Float32(-0x1.f05406p-5) < x < Float32(0x1.082b58p-4)
            return log_proc2(x)
        end

        # Step 3
        z = 1f0 + x
        zu = reinterpret(UInt32,z)
        s = reinterpret(Float32,0x7f000000 - (zu & 0xff80_0000)) # 2^-m
        m = int(zu >> 23) & 0x00ff - 127 # z cannot be subnormal
        c = m > 0 ? 1f0-(z-x) : x-(z-1f0) # 1+x = z+c
        y = reinterpret(Float32,(zu & 0x007f_ffff) | 0x3f80_0000)

        mf = Float32(m)
        F = (y + Float32(0x1p16)) - Float32(0x1p16) # 0x1p-7*round(0x1p7*y)
        f = (y - F) + s*c #2^m(F+f) = 1+x = z+c
        jp = unsafe_trunc(Int,Float32(0x1p7)*F)-127

        log_proc1(y,mf,F,f,jp)
    elseif x == -1f0
        -Inf32
    elseif isnan(x)
        NaN32
    else
        throw(DomainError())
    end
end

@vectorize_1arg Real log_tang
@vectorize_1arg Real log1p_tang


text:
@inline function _mutable_struct_equality(a::T, b::T)::Bool where T
    if isstructtype(T) && T.mutable
        T_fieldnames = fieldnames(T)
        T_num_fieldnames = length(T_fieldnames)
        temp_vector = Vector{Bool}(undef, T_num_fieldnames)
        for i in 1:T_num_fieldnames
            name_i = T_fieldnames[i]
            temp_vector[i] = _mutable_struct_equality(getfield(a, name_i), getfield(b, name_i))
        end
        return all(temp_vector)
    else
        return a == b
    end
end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ e7a4bf7b-f5f5-4564-9bca-7c69b794f8ce
md"# Generate a Mandelbrot set with Julia"

# ╔═╡ eb23f2c6-61f4-4a82-b795-033f6f2a0674
md"Notebook from [HelloJulia.jl](https://github.com/ablaom/HelloJulia.jl)"

# ╔═╡ ec6ad8e5-c854-41db-935f-ddf6a6bfbbdd
md"Instantiate package environment:"

# ╔═╡ 754d60ac-4cb0-443a-af2a-25235e544a31
begin
  using Pkg
  Pkg.activate(@__DIR__)
  Pkg.instantiate()
end

# ╔═╡ 045ab7e2-eb64-4301-8af5-148d95ea2900
begin
  using Plots
  plotly() # choose plotting backend
end

# ╔═╡ e6ca23dc-5cc4-4ec0-a6bf-47194d7e8e12
function mandel(z)
    c = z     # starting value and constant shift
    max_iterations = 20
    for n = 1:max_iterations
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return max_iterations
end

# ╔═╡ 9c765694-8292-4f26-828a-4cc485149963
begin
  xs = -2.5:0.005:0.75
  ys = -1.5:0.005:1.5
  z = [mandel(x + im*y) for y in ys, x in xs]
  heatmap(xs, ys, z)
end

# ╔═╡ da7f8fbf-e71a-48d6-9e54-3f7b3ca87ecb
savefig("mandelbrot.html")

# ╔═╡ 135dac9b-0bd9-4e1d-ba1f-6363f43ec697
md"""
---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
"""

# ╔═╡ Cell order:
# ╟─e7a4bf7b-f5f5-4564-9bca-7c69b794f8ce
# ╟─eb23f2c6-61f4-4a82-b795-033f6f2a0674
# ╟─ec6ad8e5-c854-41db-935f-ddf6a6bfbbdd
# ╠═754d60ac-4cb0-443a-af2a-25235e544a31
# ╠═045ab7e2-eb64-4301-8af5-148d95ea2900
# ╠═e6ca23dc-5cc4-4ec0-a6bf-47194d7e8e12
# ╠═9c765694-8292-4f26-828a-4cc485149963
# ╠═da7f8fbf-e71a-48d6-9e54-3f7b3ca87ecb
# ╟─135dac9b-0bd9-4e1d-ba1f-6363f43ec697


text:
module TestAccessors
import PerformanceTestTools
import Accessors

using Documenter: doctest
if VERSION == v"1.6"
    # ⨟ needs to be defined
    doctest(Accessors)
else
    @info "Skipping doctests, on VERSION = $VERSION"
end

include("test_optics.jl")
include("test_examples.jl")
include("test_staticarrays.jl")
include("test_quicktypes.jl")
include("test_setmacro.jl")
include("test_setindex.jl")
include("test_core.jl")
include("test_functionlenses.jl")
PerformanceTestTools.@include("perf.jl")

end  # module


text:
ChainRulesCore.@scalar_rule(
    betalogpdf(α::Real, β::Real, x::Number),
    @setup(z = digamma(α + β)),
    (
        log(x) + z - digamma(α),
        log1p(-x) + z - digamma(β),
        (α - 1) / x + (1 - β) / (1 - x),
    ),
)

ChainRulesCore.@scalar_rule(
    binomlogpdf(n::Real, p::Real, k::Real),
    @setup(z = digamma(n - k + 1)),
    (
        ChainRulesCore.NoTangent(),
        (k / p - n) / (1 - p),
        ChainRulesCore.NoTangent(),
    ),
)

ChainRulesCore.@scalar_rule(
    chisqlogpdf(k::Real, x::Number),
    @setup(hk = k / 2),
    (
        (log(x) - logtwo - digamma(hk)) / 2,
        (hk - 1) / x - one(hk) / 2,
    ),
)

ChainRulesCore.@scalar_rule(
    fdistlogpdf(ν1::Real, ν2::Real, x::Number),
    @setup(
        xν1 = x * ν1,
        temp1 = xν1 + ν2,
        a = (x - 1) / temp1,
        νsum = ν1 + ν2,
        di = digamma(νsum / 2),
    ),
    (
        (-log1p(ν2 / xν1) - ν2 * a + di - digamma(ν1 / 2)) / 2,
        (-log1p(xν1 / ν2) + ν1 * a + di - digamma(ν2 / 2)) / 2,
        ((ν1 - 2) / x - ν1 * νsum / temp1) / 2,
    ),
)

ChainRulesCore.@scalar_rule(
    gammalogpdf(k::Real, θ::Real, x::Number),
    @setup(
        invθ = inv(θ),
        xoθ = invθ * x,
        z = xoθ - k,
    ),
    (
        log(xoθ) - digamma(k),
        invθ * z,
        - (1 + z) / x,
    ),
)

ChainRulesCore.@scalar_rule(
    poislogpdf(λ::Number, x::Number),
    ((iszero(x) && iszero(λ) ? zero(x / λ) : x / λ) - 1, ChainRulesCore.NoTangent()),
)

ChainRulesCore.@scalar_rule(
    tdistlogpdf(ν::Real, x::Number),
    @setup(
        νp1 = ν + 1,
        xsq = x^2,
        invν = inv(ν),
        a = xsq * invν,
        b = νp1 / (ν + xsq),
    ),
    (
        (digamma(νp1 / 2) - digamma(ν / 2) + a * b - log1p(a) - invν) / 2,
        - x * b,
    ),
)


text:
using StaticArrays
using MieSeries
using Plots

n_terms = 80   # number of terms in Mie series
a = 2          # sphere radius

# normalized monostatic RCS (σ/(πa^2))
obs_point = SVector(1e-7, 1e-7, -1.)          # observation point
radius_in_λ = range(0.01, 1.6, length=100)
λ_list = a ./ radius_in_λ                     # wavelengths
σ_mono_list = Float64[]
for λ in λ_list
    k = 2π/λ       # wavenumber
    E_far = MieSeries.mieseries(obs_point; k, a, n_terms, farfield=true)
    E_far_norm2 = sum(abs2.(E_far))
    σ_mono = 4*π*E_far_norm2
    push!(σ_mono_list, σ_mono)
end
σ_mono_list /= π*a^2;  # normalize

## plot
plot(radius_in_λ, σ_mono_list, yscale=:log10, legend=false)
hline!([1], linestyle=:dash)
xlabel!("Sphere radius in wavelengths")
ylabel!("Normalized monostatic RCS")

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
@testset "LeNet5" begin
    # Test that networks can be constructed, have the right length and that
    # outputs have the right shape
    batchsize = 5

    # Forward net:
    lenet = lenet5()
    input = rand(28, 28, 1, batchsize)
    @test length(lenet) == 9
    @test size(lenet(input)) == (10, batchsize)
    lenet = lenet5(pad=0)
    input = rand(32, 32, 1, batchsize)
    @test length(lenet) == 9
    @test size(lenet(input)) == (10, batchsize)

    # Feedback net:
    lenet = lenet5_fb()
    @test length(lenet) == 13
    input = rand(28, 28, 1, batchsize)
    net = wrapfb_lenet5(lenet, batchsize)
    @test net.init isa Dict
    @test net.state isa Dict
    @test net.init["conv2"] == zeros(10, 10, 16, batchsize)
    @test net.init["fc1"] == zeros(1, 1, 120, batchsize)
    @test size(net(input)) == (10, batchsize)
end # @testset "LeNet5"


text:
Base.@kwdef struct ParallelQuickSortAlg{Alg,SmallSize,BaseSize} <: ParallelSortAlgorithm
    smallsort::Alg = Base.Sort.DEFAULT_UNSTABLE
    smallsize::SmallSize = nothing  # lazily determined
    basesize::BaseSize = 10_000
end
# `basesize` is tuned using `Float64`.  Make it `eltype`-aware?

function Base.sort!(
    v::AbstractVector,
    lo::Integer,
    hi::Integer,
    a::ParallelQuickSortAlg,
    o::Ordering,
)
    if a.basesize === nothing
        a = @set a.basesize = default_basesize(hi - lo + 1)
    end
    if a.smallsize === nothing
        a = @set a.smallsize = a.basesize
    end
    ys = view(v, lo:hi)
    _quicksort!(
        similar(ys),
        ys,
        a,
        o,
        Vector{Int8}(undef, length(ys)),
        false,  # ys_is_result
        true,   # mutable_xs
    )
    return v
end

function _quicksort!(
    ys,
    xs,
    alg,
    order,
    cs = Vector{Int8}(undef, length(ys)),
    ys_is_result = true,
    mutable_xs = false,
)
    @check length(ys) == length(xs)
    if length(ys) <= max(8, alg.smallsize)
        if ys_is_result
            zs = copyto!(ys, xs)
        else
            zs = xs
        end
        return sort!(zs, alg.smallsort, order)
    end
    pivot = _median(
        order,
        (
            xs[1],
            xs[end÷8],
            xs[end÷4],
            xs[3*(end÷8)],
            xs[end÷2],
            xs[5*(end÷8)],
            xs[3*(end÷4)],
            xs[7*(end÷8)],
            xs[end],
        ),
    )

    # TODO: Calculate extrema during the first pass if it's possible
    # to use counting sort.
    # TODO: When recursing, fuse copying _from_ `ys` to `xs` with the
    # first pass.

    # Compute sizes of each partition for each chunks.
    chunks = zip(_partition(xs, alg.basesize), _partition(cs, alg.basesize))
    results = maptasks(partition_sizes!(pivot, order), chunks)
    nbelows::Vector{Int} = map(first, results)
    nequals::Vector{Int} = map(last, results)
    naboves::Vector{Int} =
        [length(c) - (b + e) for (b, e, (c, _)) in zip(nbelows, nequals, chunks)]
    @check length(chunks) == length(nbelows) == length(nequals) == length(naboves)
    @check all(>=(0), naboves)
    singleton_chunkid = map(nbelows, nequals, naboves) do nb, ne, na
        if (nb > 0) + (ne > 0) + (na > 0) == 1
            return 1 * (nb > 0) + 2 * (ne > 0) + 3 * (na > 0)
        else
            return 0
        end
    end

    below_offsets = copy(nbelows)
    equal_offsets = copy(nequals)
    above_offsets = copy(naboves)
    acc = exclusive_cumsum!(below_offsets)
    acc = exclusive_cumsum!(equal_offsets, acc)
    acc = exclusive_cumsum!(above_offsets, acc)
    @check acc == length(xs)

    @sync begin
        for (i, (xs_chunk, cs_chunk)) in enumerate(chunks)
            singleton_chunkid[i] > 0 && continue
            @spawn unsafe_quicksort_scatter!(
                ys,
                xs_chunk,
                cs_chunk,
                below_offsets[i],
                equal_offsets[i],
                above_offsets[i],
            )
        end
        for (i, (xs_chunk, _)) in enumerate(chunks)
            singleton_chunkid[i] > 0 || continue
            idx = (
                below_offsets[i]+1:get(below_offsets, i + 1, equal_offsets[1]),
                equal_offsets[i]+1:get(equal_offsets, i + 1, above_offsets[1]),
                above_offsets[i]+1:get(above_offsets, i + 1, length(ys)),
            )[singleton_chunkid[i]]
            # There is only one partition. Short-circuit scattering.
            ys_chunk = view(ys, idx)
            copyto!(ys_chunk, xs_chunk)
            # Is it better to multi-thread this?
        end
    end

    partitions = (1:equal_offsets[1], above_offsets[1]+1:length(xs))
    @sync begin
        for idx in partitions
            length(idx) <= alg.smallsize && continue
            ys_new = view(ys, idx)
            xs_new = view(xs, idx)
            cs_new = view(cs, idx)
            @spawn let zs
                if mutable_xs
                    zs = xs_new
                else
                    zs = similar(ys_new)
                end
                _quicksort!(zs, ys_new, alg, order, cs_new, !ys_is_result, true)
            end
        end
        for idx in partitions
            length(idx) <= alg.smallsize || continue
            if ys_is_result
                ys_new = view(ys, idx)
            else
                ys_new = copyto!(view(xs, idx), view(ys, idx))
            end
            sort!(ys_new, alg.smallsort, order)
        end
        if !ys_is_result
            let idx = equal_offsets[1]+1:above_offsets[1]
                copyto!(view(xs, idx), view(ys, idx))
            end
        end
    end

    return ys_is_result ? ys : xs
end

partition_sizes!(pivot, order) = ((xs, cs),) -> partition_sizes!(xs, cs, pivot, order)

function partition_sizes!(xs, cs, pivot, order)
    nbelows = 0
    nequals = 0
    @inbounds for i in eachindex(xs, cs)
        x = xs[i]
        b = Base.lt(order, x, pivot)
        a = Base.lt(order, pivot, x)
        cs[i] = ifelse(b, -Int8(1), ifelse(a, Int8(1), Int8(0)))
        nbelows += Int(b)
        nequals += Int(!(a | b))
    end
    return (nbelows, nequals)
end

function unsafe_quicksort_scatter!(
    ys,
    xs_chunk,
    cs_chunk,
    below_offset,
    equal_offset,
    above_offset,
)
    b = below_offset
    e = equal_offset
    a = above_offset
    _foldl((b, a, e), Unroll{4}(eachindex(xs_chunk, cs_chunk))) do (b, a, e), i
        @inbounds x = xs_chunk[i]
        @inbounds c = cs_chunk[i]
        is_equal = c == 0
        is_above = c > 0
        is_below = c < 0
        e += Int(is_equal)
        a += Int(is_above)
        b += Int(is_below)
        @inbounds ys[ifelse(is_equal, e, ifelse(is_above, a, b))] = x
        (b, a, e)
    end
    return
end


text:
# Alignment Operations
# ====================
#
# Alignment operation type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# An alignment operation type
bitstype 8 Operation

# Conversion to and from integers
# -------------------------------

Base.convert(::Type{Operation}, op::UInt8) = reinterpret(Operation, op)
Base.convert(::Type{UInt8}, op::Operation) = reinterpret(UInt8, op)

Base.convert{T<:Unsigned}(::Type{Operation}, op::T) = Operation(UInt8(op))
Base.convert{T<:Unsigned}(::Type{T}, op::Operation) = T(UInt8(op))

# Operation encoding definitions
# ------------------------------

const OP_MATCH        = convert(Operation, 0x00) # M
const OP_INSERT       = convert(Operation, 0x01) # I
const OP_DELETE       = convert(Operation, 0x02) # D
const OP_SKIP         = convert(Operation, 0x03) # N
const OP_SOFT_CLIP    = convert(Operation, 0x04) # S
const OP_HARD_CLIP    = convert(Operation, 0x05) # H
const OP_PAD          = convert(Operation, 0x06) # P
const OP_SEQ_MATCH    = convert(Operation, 0x07) # =
const OP_SEQ_MISMATCH = convert(Operation, 0x08) # X
const OP_BACK         = convert(Operation, 0x09) # B
const OP_START        = convert(Operation, 0x0a) # 0 (non-standard)
const OP_INVALID      = convert(Operation, 0xff)

const OP_MAX_VALID = OP_START

# classify operations
function ismatchop(op::Operation)
    return op == OP_MATCH || op == OP_SEQ_MATCH || op == OP_SEQ_MISMATCH
end

function isinsertop(op::Operation)
    return op == OP_INSERT || op == OP_SOFT_CLIP || op == OP_HARD_CLIP
end

function isdeleteop(op::Operation)
    return op == OP_DELETE || op == OP_SKIP
end

# Conversion from characters to operations
# ----------------------------------------

# Lookup table for conversion from Char to Operation
const char_to_op = [
    OP_SEQ_MATCH, OP_INVALID, OP_INVALID,   OP_INVALID,
    OP_INVALID,   OP_BACK,    OP_INVALID,   OP_DELETE,
    OP_INVALID,   OP_INVALID, OP_INVALID,   OP_HARD_CLIP,
    OP_INSERT,    OP_INVALID, OP_INVALID,   OP_INVALID,
    OP_MATCH,     OP_SKIP,    OP_INVALID,   OP_PAD,
    OP_INVALID,   OP_INVALID, OP_SOFT_CLIP, OP_INVALID,
    OP_INVALID,   OP_INVALID, OP_INVALID,   OP_SEQ_MISMATCH ]

function Base.convert(::Type{Operation}, c::Char)
    @inbounds op = '=' <= c <= 'X' ? char_to_op[c - '=' + 1] :
                          c == '0' ? OP_START                : OP_INVALID
    if op == OP_INVALID
        error("Invalid alignment operation '$c'")
    end
    return op
end

# Conversion from characters to operations
# ----------------------------------------

const op_to_char = ['M', 'I', 'D', 'N', 'S', 'H', 'P', '=', 'X', 'B', '0']

function Base.convert(::Type{Char}, op::Operation)
    @assert op != OP_INVALID error("Alignment operation is not valid.")
    return op_to_char[convert(UInt8, op) + 1]
end

function Base.show(io::IO, op::Operation)
    write(io, convert(Char, op))
    return
end