text:
"""
Selection of random objects on the screen is realized by rendering an
object id + plus an arbitrary index into the framebuffer.
The index can be used for e.g. instanced geometries.
"""
struct SelectionID{T <: Integer} <: FieldVector{2, T}
    id::T
    index::T
end

mutable struct GLFramebuffer
    resolution::Observable{NTuple{2, Int}}
    id::NTuple{2, GLuint}

    buffers::Dict{Symbol, Texture}
    render_buffer_ids::Vector{GLuint}
end

# it's guaranteed, that they all have the same size
Base.size(fb::GLFramebuffer) = size(fb.buffers[:color])

function attach_framebuffer(t::Texture{T, 2}, attachment) where T
    glFramebufferTexture2D(GL_FRAMEBUFFER, attachment, GL_TEXTURE_2D, t.id, 0)
end

function GLFramebuffer(fb_size::NTuple{2, Int})
    # First Framebuffer
    render_framebuffer = glGenFramebuffers()
    glBindFramebuffer(GL_FRAMEBUFFER, render_framebuffer)

    color_buffer = Texture(RGBA{N0f8}, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)
    objectid_buffer = Texture(Vec{2, GLuint}, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)

    depth_buffer = Texture(
        Ptr{GLAbstraction.DepthStencil_24_8}(C_NULL), fb_size,
        minfilter = :nearest, x_repeat = :clamp_to_edge,
        internalformat = GL_DEPTH24_STENCIL8,
        format = GL_DEPTH_STENCIL
    )

    attach_framebuffer(color_buffer, GL_COLOR_ATTACHMENT0)
    attach_framebuffer(objectid_buffer, GL_COLOR_ATTACHMENT1)
    attach_framebuffer(depth_buffer, GL_DEPTH_ATTACHMENT)
    attach_framebuffer(depth_buffer, GL_STENCIL_ATTACHMENT)

    status = glCheckFramebufferStatus(GL_FRAMEBUFFER)
    @assert status == GL_FRAMEBUFFER_COMPLETE

    # Second Framebuffer
    # postprocessor adds buffers here
    color_luma_framebuffer = glGenFramebuffers()
    glBindFramebuffer(GL_FRAMEBUFFER, color_luma_framebuffer)

    @assert status == GL_FRAMEBUFFER_COMPLETE

    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    fb_size_node = Observable(fb_size)

    buffers = Dict(
        :color => color_buffer,
        :objectid => objectid_buffer,
        :depth => depth_buffer
    )

    return GLFramebuffer(
        fb_size_node,
        (render_framebuffer, color_luma_framebuffer),
        buffers,
        [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1]
    )
end

function Base.resize!(fb::GLFramebuffer, window_size)
    ws = Int.((window_size[1], window_size[2]))
    if ws != size(fb) && all(x-> x > 0, window_size)
        for (name, buffer) in fb.buffers
            resize_nocopy!(buffer, ws)
        end
        fb.resolution[] = ws
    end
    nothing
end

struct MonitorProperties
    name::String
    isprimary::Bool
    position::Vec{2, Int}
    physicalsize::Vec{2, Int}
    videomode::GLFW.VidMode
    videomode_supported::Vector{GLFW.VidMode}
    dpi::Vec{2, Float64}
    monitor::GLFW.Monitor
end

function MonitorProperties(monitor::GLFW.Monitor)
    name = GLFW.GetMonitorName(monitor)
    isprimary = GLFW.GetPrimaryMonitor() == monitor
    position = Vec{2, Int}(GLFW.GetMonitorPos(monitor)...)
    physicalsize = Vec{2, Int}(GLFW.GetMonitorPhysicalSize(monitor)...)
    videomode = GLFW.GetVideoMode(monitor)
    sfactor = Sys.isapple() ? 2.0 : 1.0
    dpi = Vec(videomode.width * 25.4, videomode.height * 25.4) * sfactor ./ Vec{2, Float64}(physicalsize)
    videomode_supported = GLFW.GetVideoModes(monitor)

    MonitorProperties(name, isprimary, position, physicalsize, videomode, videomode_supported, dpi, monitor)
end

was_destroyed(nw::GLFW.Window) = nw.handle == C_NULL

function GLContext()
    context = GLFW.GetCurrentContext()
    version = opengl_version_number()
    glsl_version = glsl_version_number()
    return GLContext(context, version, glsl_version, unique_context_counter())
end

function ShaderAbstractions.native_switch_context!(x::GLFW.Window)
    GLFW.MakeContextCurrent(x)
end

function ShaderAbstractions.native_context_alive(x::GLFW.Window)
    GLFW.is_initialized() && !was_destroyed(x)
end

function destroy!(nw::GLFW.Window)
    was_current = ShaderAbstractions.is_current_context(nw)
    if !was_destroyed(nw)
        GLFW.SetWindowShouldClose(nw, true)
        GLFW.PollEvents()
        GLFW.DestroyWindow(nw)
        nw.handle = C_NULL
    end
    was_current && ShaderAbstractions.switch_context!()
end

function windowsize(nw::GLFW.Window)
    was_destroyed(nw) && return (0, 0)
    size = GLFW.GetFramebufferSize(nw)
    return (size.width, size.height)
end

function Base.isopen(window::GLFW.Window)
    was_destroyed(window) && return false
    try
        return !GLFW.WindowShouldClose(window)
    catch e
        # can't be open if GLFW is already terminated
        e.code == GLFW.NOT_INITIALIZED && return false
        rethrow(e)
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings <: SwaggerModel
    domainNameLabel::Any # spec type: Union{ Nothing, String } # spec name: domainNameLabel

    function VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings(;domainNameLabel=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings, Symbol("domainNameLabel"), domainNameLabel)
        setfield!(o, Symbol("domainNameLabel"), domainNameLabel)
        o
    end
end # type VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings

const _property_map_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings = Dict{Symbol,Symbol}(Symbol("domainNameLabel")=>Symbol("domainNameLabel"))
const _property_types_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings = Dict{Symbol,String}(Symbol("domainNameLabel")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings }) = collect(keys(_property_map_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings))
Swagger.property_type(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings[property_name]

function check_required(o::VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings)
    (getproperty(o, Symbol("domainNameLabel")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings }, name::Symbol, val)
end

ext: jl


text:
using jInv.Mesh;
using Multigrid

# function getNodalLaplacianMatrixDirichlet(Minv::RegularMesh)
# newDomain = copy(Minv.domain);
# newDomain[1] -= Minv.h[1];
# newDomain[2] += Minv.h[1];
# newDomain[3] -= Minv.h[2];
# newDomain[4] += Minv.h[2];
# MinvPadded = getRegularMesh(newDomain,Minv.n+2);
# mask = zeros(Bool,tuple(MinvPadded.n+1...));
# mask[2:end-1,2:end-1] = true;
# A = getNodalLaplacianMatrix(MinvPadded);
# A = A[mask[:],mask[:]];
# return A;
# end

# A = getNodalLaplacianMatrixDirichlet(getRegularMesh([0.0,1.0,0.0,1.0],[4,4]))

m = ones(129,129);

Minv = getRegularMesh([0.0,1.0,0.0,1.0],collect(size(m))-1);

levels      = 4;
numCores 	= 2;
maxIter     = 5;
relativeTol = 1e-2;
relaxType   = "Jac";
relaxParam  = 0.8;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "NoMUMPS";

MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,
				relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0)

## TEST DIRICHLET
# L = getNodalLaplacianMatrixDirichlet(Minv);
# s = zeros(tuple(Minv.n+1...));
# s[2:end-1,2:end-1] = rand(tuple(Minv.n-1...));
# s = s[:];
# MGsetup(getNodalLaplacianMatrixDirichlet,Minv,MG,Float64,1,true)
# MGsetup(L,Minv,MG,Float64,1,true)

## TEST NEUMANN
L = getNodalLaplacianMatrix(Minv);
s = rand(prod(Minv.n+1));
b = L*s;
x0 = 0.0*s;

MGsetup(getMultilevelOperatorConstructor([],getNodalLaplacianMatrix,[]),Minv,MG,Float64,1,true);
solveMG(MG,b,x0,true);

xn = getCellCenteredGrid(Minv);
sig = 3*xn[:,1].*(1.0-xn[:,1]) + 2*xn[:,2].*(1.0-xn[:,2]);
sig = sig[:];
# sig[:]=1.0;

## The following lines also appear in a different version in jInv.Mesh.
## The difference is the averaging: in jInv.Mesh sigma is averaged with 0.0 outside the domain (there's a difference in how av is defined).
## geometric MG does not work with the jInv version. only with this one.

function av(n)
# A = av(n), 1D average operator
	av = spdiagm((fill(.5,n),fill(.5,n)),(-1,0),n+1,n);
	av[1,1] = 1.0;
	av[end,end] = 1.0;
	return av;
end
function getEdgeAverageMatrix(Mesh::AbstractTensorMesh)
	if Mesh.dim==3
		A1 = kron(av(Mesh.n[3]),kron(av(Mesh.n[2]),speye(Mesh.n[1])))
		A2 = kron(av(Mesh.n[3]),kron(speye(Mesh.n[2]),av(Mesh.n[1])))
		A3 = kron(speye(Mesh.n[3]),kron(av(Mesh.n[2]),av(Mesh.n[1])))
		Mesh.Ae = [A1; A2; A3]
	elseif Mesh.dim==2
		A1 = kron(av(Mesh.n[2]),speye(Mesh.n[1]))
		A2 = kron(speye(Mesh.n[2]),av(Mesh.n[1]))
		Ae = [A1; A2]
	else
		error("getEdgeAverageMatrix not implemented fot $(Mesh.dim)D Meshes")
	end
	return Ae
end
function getDivSigGrad(M::AbstractTensorMesh,sig::Vector{Float64})
 	G       = getNodalGradientMatrix(M)
	Ae      = getEdgeAverageMatrix(M);
    A       = (G'*spdiagm(Ae*(vec(sig))))*G
	return A
end

Ar = getDivSigGrad(Minv,sig);
b = Ar*s;
MGsetup(Ar,Minv,MG,Float64,1,true);
x0 = 0.0*s;
solveMG(MG,b,x0,true);

# getDivSigGrad(mesh::RegularMesh,param::Array{Float64,1}) = getDivSigGradMatrix(param,mesh);
restrictSigma(mesh_fine,mesh_coarse,param_fine,level) = restrictCellCenteredVariables(param_fine,mesh_fine.n);
MGsetup(getMultilevelOperatorConstructor(sig,getDivSigGrad,restrictSigma),Minv,MG,Float64,1,true);
x0 = 0.0*s;
solveMG(MG,b,x0,true);
println("DONE!")


text:
module TestSamples

using DataFrames
using Orihime
using Test

function plot_lines(n = 100)
    df = DataFrame(
        a = repeat([1, 2, 3, 4], outer = [n]),
        b = repeat(1:n, outer = [4]),
        c = cumsum(randn(n * 4)),
    )
    plt = Orihime.@plot(
        :line,
        x = {:b, type = :quantitative},
        y = {:c, type = :quantitative},
        extra_facets = [:a],
        data = df
    )
    return plt
end

function plot_hists(n = 100)
    df = DataFrame(
        a = repeat([1, 2, 3, 4], outer = [n]),
        b = repeat(1:2, outer = [2n]),
        c = randn(n * 4),
    )
    plt = Orihime.@plot(
        :bar,
        x = {:c, bin = true},
        y = {aggregate = :count},
        extra_facets = [:a, :b],
        data = df
    )
    return plt
end

function check_smoke(f)
    plt = f()
    io = IOBuffer()
    show(io, "image/png", plt)
    @test position(io) > 0
end

function test_smoke()
    @testset for f in [plot_lines, plot_hists]
        check_smoke(f)
    end
end

end  # module


text:
import FractionalCalculus.FracDiffAlg

"""
Riesz sense fractional derivative
"""
abstract type Riesz <: FracDiffAlg end

"""
# Riesz sense symmetric fractional derivative algorithm.

    fracdiff(f, α, end_point, h, Riesz_Symmetric())

Compute fractional derivative of Riesz sense using Triangular Strip Matrix algorithm.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.01, Riesz_Symmetric())
```
"""
struct Riesz_Symmetric <: Riesz end

################################################################
###                    Type defination done                  ###
################################################################

function fracdiff(f, α, end_point, h, ::Riesz_Symmetric)
    N=Int(floor(end_point/h))

    mat = RieszMatrix(α, N+1, h)
    return mat*f.(collect(0:h:end_point))
end
function RieszMatrix(α, N, h)
    caputo = B(N+1, α)
    caputo = caputo[2:(N+1), 1:N]
    result = 1/2*(caputo+caputo')
    result = h^(-α)*result

    return result
end

function B(N, p)
    result=zeros(N, N)
    temp=omega(N, p)

    @inbounds @simd for i ∈ 1:N
        @views result[i, 1:i]=reverse(temp[1:i])
    end

    return result
end

function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    @fastmath @inbounds @simd for i ∈ 1:n
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega
end


text:
using PyCall, Pkg.TOML

export era5download, monthlyera5download, saveconfig, download_and_convert_era5

getconfig(key) = getconfig()[key]

function getconfig()
    configfile = joinpath(homedir(), ".GlobalEnergyGIS_config")
    if !isfile(configfile)
        error("Configuration file missing, please run saveconfig(datafolder, uid, api_key) first. See GlobalEnergyGIS README.")
    end
    return TOML.parsefile(configfile)
end

function saveconfig(datafolder::AbstractString, uid::Int, api_key::AbstractString; agree_terms=false)
    !agree_terms && error("You must agree to the terms of use of all datasets to proceed. See GlobalEnergyGIS README.")
    downloadsfolder = joinpath(datafolder, "downloads")
    mkpath(downloadsfolder)
    configfile = joinpath(homedir(), ".GlobalEnergyGIS_config")
    open(configfile, "w") do io
        d = Dict("datafolder"=>datafolder, "agree_terms"=>agree_terms)
        TOML.print(io, d)
        println("Configuration file written to $configfile.")
    end
    cds_id(uid, api_key)
end

function cds_id(uid::Int, api_key::AbstractString)
    filename = joinpath(homedir(), ".cdsapirc")
    isfile(filename) && error("$filename already exists, no changes made. Please check its contents manually.")
    open(filename, "w") do file
        write(file, "url: https://cds.climate.copernicus.eu/api/v2\n")
        write(file, "key: $uid:$api_key\n")
        println("Copernicus credentials written to $filename.")
    end
end

function download_and_convert_era5(year=2018; datasets=["wind", "solar", "temp"])
    for dataset in datasets
        println("\nDownloading ERA5 $dataset data from Copernicus...")
        era5download(year; datasets=[dataset])
        println("\nConverting downloaded $dataset data to HDF5 and recompressing...")
        if dataset == "solar"
            makesolarera5(; year)
        elseif dataset == "wind"
            makewindera5(; year)
        elseif dataset == "temp"
            maketempera5(; year)
        end
        println("\nCleanup: deleting downloaded $dataset data...")
        clearvars_era5(; year, datasets=[dataset])
    end
    println("\nERA5 datasets $datasets downloaded and converted. Temporary files cleaned up.")
end

function era5download(year=2018; datasets=["wind", "solar", "temp"])
    mkpath(in_datafolder("downloads"))
    count = 0
    for dataset in datasets, month = 1:12, monthhalf = 1:2
        if dataset == "wind"
            vars = ["100m_u_component_of_wind", "100m_v_component_of_wind"]
        elseif dataset == "solar"
            vars = ["surface_solar_radiation_downwards", "total_sky_direct_solar_radiation_at_surface"]
        else
            vars = ["2m_temperature"]
        end
        if monthhalf == 1
            firstday, lastday = "01", "15"
        else
            firstday = "16"
            lastday = Dates.daysinmonth(Date("$year-$month"))
        end
        monthstr = lpad(month,2,'0')
        date1, date2 = "$year-$monthstr-$firstday", "$year-$monthstr-$lastday"
        outfile = in_datafolder("downloads", "$dataset$year-$monthstr$firstday-$monthstr$lastday.nc")
        count += 1
        println("\nFile $count of $(24*length(datasets)):")
        request_era5_vars(outfile, vars, date1, date2)
    end
end

function monthlyera5download(; datasets=["wind", "solar", "temp"])
    mkpath(in_datafolder("downloads"))
    for dataset in datasets
        if dataset == "wind"
            vars = ["100m_u_component_of_wind", "100m_v_component_of_wind"]
        elseif dataset == "solar"
            vars = ["surface_solar_radiation_downwards", "total_sky_direct_solar_radiation_at_surface"]
        else
            vars = ["2m_temperature"]
        end
        # avoid 2020 to avoid partial-year problems and near real-time ERA5T data
        # (which includes an additional dimension to indicate ERA5/ERA5T)
        years = 1979:2019
        outfile = in_datafolder("downloads", "monthly$(dataset)_$(years[1])-$(years[end]).nc")
        request_monthly_era5_vars(outfile, vars, collect(years))
    end
end

function request_era5_vars(outfile::String, vars::Vector{String}, firstdate::String, lastdate::String)
    datestring = "$firstdate/$lastdate"
    py"""
    import cdsapi

    c = cdsapi.Client()
    c.retrieve(
        'reanalysis-era5-single-levels',
        {
            'product_type': 'reanalysis',
            'format': 'netcdf',
            'variable': $vars,
            'grid': '0.28125/0.28125',
            'area': '89.859375/-179.859375/-89.859375/179.859375',
            'date': $datestring,
            'time': '00/to/23/by/1'
        },
        $outfile)
    """
end

function request_monthly_era5_vars(outfile::String, vars::Vector{String}, years::Vector{Int})
    stryears = string.(years)
    months = string.(1:12)
    py"""
    import cdsapi

    c = cdsapi.Client()
    c.retrieve(
        'reanalysis-era5-single-levels-monthly-means',
        {
            'product_type': 'monthly_averaged_reanalysis',
            'format': 'netcdf',
            'variable': $vars,
            'grid': '0.28125/0.28125',
            'area': '89.859375/-179.859375/-89.859375/179.859375',
            'year': $stryears,
            'month': $months,
            'time': '00:00'
        },
        $outfile)
    """
end

# For some reason the delivered NetCDF files are unreadable unless they are limited to 15-16 days each.

# default resolution:  'grid': '0.25/0.25'  (but the max ERA5 resolution is 0.28125 degrees, so use that instead)
# default area:        'area': '90/0/-90/179.75'  (note max latitude is 90)

# The default grid size using a resolution of 0.25 is 721x1440 pixels.
# The default grid size using a resolution of 0.28125 is 641x1280 pixels.

# 'area': '89.859375/-179.859375/-89.859375/179.859375': get correct center points for a 640x1280 grid (north/west/south/east)

# The CDS server doesn't use high enough precision to evaluate the area coordinates exactly. Each lon & lat is off
# by about 0.000625 degrees, but this error is very small compared to the resolution (0.28125 degrees). Otherwise
# using the default grid we would have to interpolate to get correct center points, and this would also lose some
# precision. This way we avoid doing all that computational work for a 641x1280x8760 grid.

# Note that using the default grid with a resolution of 0.28125 also makes longitude errors of up to 0.000244 degrees
# (the numbers come from inspecting latitude & longitude coordinates of ERA5 test data downloads).

# https://confluence.ecmwf.int/display/CKB/How+to+install+and+use+CDS+API+on+Windows
# https://cds.climate.copernicus.eu/api-how-to#install-the-cds-api-key
# https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5
# https://confluence.ecmwf.int/display/CKB/ERA5+data+documentation
# https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=form
# https://cds.climate.copernicus.eu/cdsapp#!/yourrequests

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(GraphIdx.Grid.compute_dirs), Int64})
    precompile(Tuple{Type{GraphIdx.Grid.GridGraph}, Int64, Int64, Int64})
    precompile(Tuple{typeof(GraphIdx.Utils.setlocale), String})
    precompile(Tuple{typeof(GraphIdx.Utils.fmt_thousands), Int64})
    precompile(Tuple{Type{GraphIdx.PrimMstMem}, GraphIdx.IncidenceIndex})
    precompile(Tuple{typeof(GraphIdx.prim_mst_edges), Array{Float64, 1}, Int64, GraphIdx.PrimMstMem})
    precompile(Tuple{typeof(GraphIdx.Tree.reset!), GraphIdx.Tree.ChildrenIndex, Array{Int64, 1}, Int64})
    precompile(Tuple{typeof(Base.sizehint!), GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, Int64})
    precompile(Tuple{typeof(GraphIdx.Utils.Heap.percolate_down!), GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, Int64})
    precompile(Tuple{typeof(GraphIdx.Utils.Heap.dequeue!), GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}})
    precompile(Tuple{typeof(GraphIdx.Utils.Heap.percolate_up!), GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, Int64})
    precompile(Tuple{typeof(GraphIdx.Utils.Heap.enqueue!), GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, Base.Pair{Int64, Float64}})
    precompile(Tuple{typeof(GraphIdx.prim_mst_edges), Array{Float64, 1}, Array{Bool, 1}, Array{Float64, 1}, Array{Int64, 1}, GraphIdx.IncidenceIndex, Array{Int64, 1}, GraphIdx.Utils.Heap.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, Int64})
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# This program creates a simple 3D vel model for FD seismic modeling
#
# by GRAM | 8 Apr 2016

using Seismic

# def size
nx = 50
ny = 50
nz = 200
dz = 10 # depth increment
dx = 20
dy = 20

vel = zeros(Float32,nz,nx,ny)

# build model
for i = 1:nz, j = 1:nx, k = 1:ny
	if i < nz / 3
		vel[i,j,k] = 2000
	elseif ((i < 2 * nz / 3) & (i > nz / 3))
		vel[i,j,k] = 3000
	elseif ((i > (2 * nz / 3) - 25) & (i < (2 * nz / 3) + 25))
		vel[i,j,k] = 3500
	else
		vel[i,j,k] = 4000
	end
end

#= display
print(size(vel), "\n")
print(vel,"\n")
=#

#= write out text file for QC
out = open("../dat/vel_1.txt","w")
for k = 1:ny
	for i = 1:nz
		for j = 1:nx
			showcompact(out, vel[i,j,k])
			write(out," ")
		end
		write(out, "\r\n")
	end
	write(out, "\r\n")
end
close(out)
=#

# write out vel model for Seismic.jl
h = Array(Header,nx*ny)
ex = Seismic.Extent(convert(Int32,nz), convert(Int32,nx), convert(Int32,ny),
	1, 1, 0, 0, 0, 0, 0, convert(Float32,dz), convert(Float32,dx),
	convert(Float32,dy), 1, 1, "Depth", "mx", "my", "", "", "", "", "",
	"", "", "")

for i = 1 : ny
	for j = 1 : nx
		h[j + nx * (i - 1)] = Seismic.InitSeisHeader()
		h[j + nx * (i - 1)].tracenum = convert(Int32, j + nx * (i - 1))
		h[j + nx * (i - 1)].o1 = convert(Float32, 0)
		h[j + nx * (i - 1)].n1 = convert(Int32, nz)
		h[j + nx * (i - 1)].d1 = convert(Float32, dz)
		h[j + nx * (i - 1)].mx = convert(Float32, j)
		h[j + nx * (i - 1)].my = convert(Float32, i)
	end
end

SeisWrite("../dat/vel_incorrect",vel,h,ex)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url)
    return GitRemote(rmt_ptr_ptr[])
end

function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString, fetch_spec::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_with_fetchspec, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url, fetch_spec)
    return GitRemote(rmt_ptr_ptr[])
end

function GitRemoteAnon(repo::GitRepo, url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_anonymous, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, url)
    return GitRemote(rmt_ptr_ptr[])
end

function get(::Type{GitRemote}, repo::GitRepo, rmt_name::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_lookup, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name)
    return GitRemote(rmt_ptr_ptr[])
end

function url(rmt::GitRemote)
    url_ptr = ccall((:git_remote_url, :libgit2), Cstring, (Ptr{Void}, ), rmt.ptr)
    url_ptr == Cstring_NULL && return ""
    return  unsafe_string(url_ptr)
end

function fetch_refspecs(rmt::GitRemote)
    sa_ref = Ref{StrArrayStruct}()
    try
        @check ccall((:git_remote_get_fetch_refspecs, :libgit2), Cint,
                      (Ptr{LibGit2.StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
        convert(Vector{AbstractString}, sa_ref[])
    finally
        finalize(sa_ref[])
    end
end

function push_refspecs(rmt::GitRemote)
    sa_ref = Ref{StrArrayStruct}()
    try
        @check ccall((:git_remote_get_push_refspecs, :libgit2), Cint,
                      (Ptr{LibGit2.StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
        convert(Vector{AbstractString}, sa_ref[])
    finally
        finalize(sa_ref[])
    end
end

function fetch{T<:AbstractString}(rmt::GitRemote, refspecs::Vector{T};
               options::FetchOptions = FetchOptions(),
               msg::AbstractString="")
    msg = "libgit2.fetch: $msg"
    no_refs = isempty(refspecs)
    !no_refs && (sa = StrArrayStruct(refspecs...))
    try
        @check ccall((:git_remote_fetch, :libgit2), Cint,
            (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{FetchOptions}, Cstring),
            rmt.ptr, no_refs ? C_NULL : Ref(sa), Ref(options), msg)
    finally
        !no_refs && finalize(sa)
    end
end

function push{T<:AbstractString}(rmt::GitRemote, refspecs::Vector{T};
                                 force::Bool = false,
                                 options::PushOptions = PushOptions())
    no_refs = isempty(refspecs)
    !no_refs && (sa = StrArrayStruct(refspecs...))
    try
        @check ccall((:git_remote_push, :libgit2), Cint,
                      (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{PushOptions}),
                       rmt.ptr, no_refs ? C_NULL : Ref(sa), Ref(options))
    finally
        !no_refs && finalize(sa)
    end
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Query

import OldPkg
import ..PkgError
using ..Types

function init_resolve_backtrace(reqs::Requires, fix::Dict{String,Fixed} = Dict{String,Fixed}())
    bktrc = ResolveBacktrace()
    for (p,f) in fix
        bktrc[p] = ResolveBacktraceItem(:fixed, f.version)
    end
    for (p,vs) in reqs
        bktrcp = get!(bktrc, p) do; ResolveBacktraceItem() end
        push!(bktrcp, :required, vs)
    end
    return bktrc
end

function check_fixed(reqs::Requires, fix::Dict{String,Fixed}, avail::Dict)
    for (p1,f1) in fix
        for p2 in keys(f1.requires)
            haskey(avail, p2) || haskey(fix, p2) || throw(PkgError("unknown package $p2 required by $p1"))
        end
        satisfies(p1, f1.version, reqs) ||
            @warn "$p1 is fixed at $(f1.version) conflicting with top-level requirement: $(reqs[p1])"
        for (p2,f2) in fix
            satisfies(p1, f1.version, f2.requires) ||
                @warn "$p1 is fixed at $(f1.version) conflicting with requirement for $p2: $(f2.requires[p1])"
        end
    end
end

function propagate_fixed!(reqs::Requires, bktrc::ResolveBacktrace, fix::Dict{String,Fixed})
    for (p,f) in fix
        merge_requires!(reqs, f.requires)
        for (rp,rvs) in f.requires
            bktrcp = get!(bktrc, rp) do; ResolveBacktraceItem() end
            push!(bktrcp, p=>bktrc[p], rvs)
        end
    end
    for (p,f) in fix
        delete!(reqs, p)
    end
    reqs
end

# Specialized copy for the avail argument below because the deepcopy is slow
function availcopy(avail)
    new_avail = empty(avail)
    for (pkg, vers_avail) in avail
        new_vers_avail = empty(vers_avail)
        for (version, pkg_avail) in vers_avail
            new_vers_avail[version] = copy(pkg_avail)
        end
        new_avail[pkg] = new_vers_avail
    end
    return new_avail
end

# Generate a reverse dependency graph (package names only)
function gen_backdeps(avail::Dict)
    backdeps = Dict{String,Set{String}}()
    for (ap,av) in avail, (v,a) in av, rp in keys(a.requires)
        s = get!(backdeps, rp) do; Set{String}() end
        push!(s, ap)
    end
    return backdeps
end

function dependencies(avail::Dict, fix::Dict = Dict{String,Fixed}("julia"=>Fixed(VERSION)))
    avail = availcopy(avail)
    conflicts = Dict{String,Set{String}}()
    to_expunge = VersionNumber[]
    emptied = String[]
    backdeps = gen_backdeps(avail)

    for (fp,fx) in fix
        delete!(avail, fp)
        haskey(backdeps, fp) || continue
        # for (ap,av) in avail
        for ap in backdeps[fp]
            haskey(avail, ap) || continue
            av = avail[ap]
            empty!(to_expunge)
            for (v,a) in av
                if satisfies(fp, fx.version, a.requires)
                    delete!(a.requires, fp)
                else
                    conflicts_ap = get!(conflicts, ap) do; Set{String}() end
                    push!(conflicts_ap, fp)
                    # don't delete v from av right away so as not to screw up iteration
                    push!(to_expunge, v)
                end
            end
            for v in to_expunge
                delete!(av, v)
            end
            isempty(av) && push!(emptied, ap)
        end
    end
    while !isempty(emptied)
        deleted_pkgs = String[]
        for ap in emptied
            delete!(avail, ap)
            push!(deleted_pkgs, ap)
        end
        empty!(emptied)

        for dp in deleted_pkgs
            haskey(backdeps, dp) || continue
            for ap in backdeps[dp]
                haskey(avail, ap) || continue
                av = avail[ap]
                empty!(to_expunge)
                for (v,a) in av
                    haskey(a.requires, dp) || continue
                    conflicts_ap = get!(conflicts, ap) do; Set{String}() end
                    union!(conflicts_ap, conflicts[dp])
                    push!(to_expunge, v)
                end
                for v in to_expunge
                    delete!(av, v)
                end
                isempty(av) && push!(emptied, ap)
            end
        end
    end
    avail, conflicts
end

function partial_update_mask(instd::Dict{String,Tuple{VersionNumber,Bool}},
        avail::Dict{String,Dict{VersionNumber,Available}}, upkgs::Set{String})
    dont_update = Set{String}()
    isempty(upkgs) && return dont_update
    avail_new = deepcopy(avail)
    for p in upkgs
        haskey(instd, p) || throw(PkgError("Package $p is not installed"))
        v = instd[p][1]
        if haskey(avail, p)
            for vn in keys(avail[p])
                vn < v && delete!(avail_new[p], vn)
            end
        end
    end
    avail_new = dependencies_subset(avail_new, upkgs)

    for p in keys(avail)
        !haskey(avail_new, p) && push!(dont_update, p)
    end
    for p in keys(instd)
        !haskey(avail_new, p) && p ∉ upkgs && push!(dont_update, p)
    end
    return dont_update
end

# Try to produce some helpful message in case of a partial update which does not go all the way
# (Does not do a full analysis, it only checks requirements and direct dependents.)
function check_partial_updates(reqs::Requires,
                               deps::Dict{String,Dict{VersionNumber,Available}},
                               want::Dict{String,VersionNumber},
                               fixed::Dict{String,Fixed},
                               upkgs::Set{String})
    for p in upkgs
        if !haskey(want, p)
            if !haskey(fixed, p)
                @warn "Something went wrong with the update of package $p, please submit a bug report"
                continue
            end
            v = fixed[p].version
        else
            v = want[p]
            if haskey(fixed, p) && v != fixed[p].version
                @warn "Something went wrong with the update of package $p, please submit a bug report"
                continue
            end
        end
        haskey(deps, p) || continue
        vers = sort!(collect(keys(deps[p])))
        higher_vers = vers[vers .> v]
        isempty(higher_vers) && continue # package p has been set to the highest available version

        # Determine if there are packages which depend on `p` and somehow prevent its update to
        # the latest version
        blocking_parents = Set{String}()
        for (p1,d1) in deps
            p1 in upkgs && continue # package `p1` is among the ones to be updated, skip the check
            haskey(fixed, p1) || continue # if package `p1` is not fixed, it can't be blocking
            r1 = fixed[p1].requires # get `p1` requirements
            haskey(r1, p) || continue # check if package `p1` requires `p`
            vs1 = r1[p] # get the versions of `p` allowed by `p1` requirements
            any(hv in vs1 for hv in higher_vers) && continue # package `p1` would allow some of the higher versions,
                                                             # therefore it's not responsible for blocking `p`
            push!(blocking_parents, p1) # package `p1` is blocking the update of `p`
        end

        # Determine if the update of `p` is prevented by explicit user-provided requirements
        blocking_reqs = (haskey(reqs, p) && all(hv ∉ reqs[p] for hv in higher_vers))

        # Determine if the update of `p` is prevented by it being fixed (e.g. it's dirty, or pinned...)
        isfixed = haskey(fixed, p)

        @info begin
            "Package $p was set to version $v, but a higher version $(vers[end]) exists.\n" *
            if isfixed
                "The package is fixed. You can try using `OldPkg.free(\"$p\")` to update it."
            elseif blocking_reqs
                "The update is prevented by explicit requirements constraints. Edit your REQUIRE file to change this."
            elseif !isempty(blocking_parents)
                string("To install the latest version, you could try updating these packages as well: ", join(blocking_parents, ", ", " and "), ".")
            else
                "To install the latest version, you could try doing a full update with `OldPkg.update()`."
            end
        end
    end
end

const PackageState = Union{Nothing,VersionNumber}

function diff(have::Dict, want::Dict, avail::Dict, fixed::Dict)
    change = Vector{Tuple{String,Tuple{PackageState,PackageState}}}()
    remove = Vector{Tuple{String,Tuple{PackageState,PackageState}}}()

    for pkg in collect(union(keys(have),keys(want)))
        h, w = haskey(have,pkg), haskey(want,pkg)
        if h && w
            if have[pkg] != want[pkg]
                push!(change, (pkg,(have[pkg], want[pkg])))
            end
        elseif h
            push!(remove, (pkg,(have[pkg],nothing)))
        elseif w
            push!(change, (pkg,(nothing,want[pkg])))
        end
    end
    append!(sort!(change), sort!(remove))
end

function check_requirements(reqs::Requires, deps::Dict{String,Dict{VersionNumber,Available}}, fix::Dict)
    for (p,vs) in reqs
        if !any(vn->(vn in vs), keys(deps[p]))
            remaining_vs = VersionSet()
            err_msg = "fixed packages introduce conflicting requirements for $p: \n"
            available_list = sort!(collect(keys(deps[p])))
            for (p1,f1) in fix
                f1r = f1.requires
                haskey(f1r, p) || continue
                err_msg *= "         $p1 requires versions $(f1r[p])"
                if !any([vn in f1r[p] for vn in available_list])
                    err_msg *= " [none of the available versions can satisfy this requirement]"
                end
                err_msg *= "\n"
                remaining_vs = intersect(remaining_vs, f1r[p])
            end
            if isempty(remaining_vs)
                err_msg *= "       the requirements are unsatisfiable because their intersection is empty"
            else
                err_msg *= "       available versions are $(join(available_list, ", ", " and "))"
            end
            throw(PkgError(err_msg))
        end
    end
end

# If there are explicitly required packages, dicards all versions outside
# the allowed range.
# It also propagates requirements: when all allowed versions of a required package
# require some other package, this creates a new implicit requirement.
# The propagation is tracked so that in case a contradiction is detected the error
# message allows to determine the cause.
# This is a pre-pruning step, so it also creates some structures which are later used by pruning
function filter_versions(reqs::Requires, deps::Dict{String,Dict{VersionNumber,Available}}, bktrc::ResolveBacktrace)
    allowed = Dict{String,Dict{VersionNumber,Bool}}()
    staged = copy(reqs)
    while !isempty(staged)
        staged_next = Requires()
        for (p,vs) in staged
            # Parse requirements and store allowed versions.
            depsp = deps[p]
            if !haskey(allowed, p)
                allowedp = Dict{VersionNumber,Bool}(vn=>true for vn in keys(depsp))
                allowed[p] = allowedp
                seen = false
            else
                allowedp = allowed[p]
                oldallowedp = copy(allowedp)
                seen = true
            end
            for vn in keys(depsp)
                allowedp[vn] &= vn ∈ vs
            end
            @assert !isempty(allowedp)
            if !any(values(allowedp))
                err_msg = "Unsatisfiable requirements detected for package $p:\n"
                err_msg *= string(bktrc[p])
                err_msg *= """The intersection of the requirements is $(bktrc[p].versionreq).
                              None of the available versions can satisfy this requirement."""
                throw(PkgError(err_msg))
            end

            # If we've seen this package already and nothing has changed since
            # the last time, we stop here.
            seen && allowedp == oldallowedp && continue

            # Propagate requirements:
            # if all allowed versions of a required package require some other package,
            # then compute the union of the allowed versions for that other package, and
            # treat that as a new requirement.
            # Start by filtering out the non-allowed versions
            fdepsp = Dict{VersionNumber,Available}(vn=>depsp[vn] for vn in keys(depsp) if allowedp[vn])
            # Collect all required packages
            isreq = Dict{String,Bool}(rp=>true for a in values(fdepsp) for rp in keys(a.requires))
            # Compute whether a required package appears in all requirements
            for rp in keys(isreq)
                isreq[rp] = all(haskey(a.requires, rp) for a in values(fdepsp))
            end

            # Create a list of candidates for new implicit requirements
            staged_new = Set{String}()
            for a in values(fdepsp), (rp,rvs) in a.requires
                # Skip packages that may not be required
                isreq[rp] || continue
                # Compute the union of the version sets
                if haskey(staged_next, rp)
                    snvs = staged_next[rp]
                    union!(snvs, rvs)
                else
                    snvs = copy(rvs)
                    staged_next[rp] = snvs
                end
                push!(staged_new, rp)
            end
            for rp in staged_new
                @assert isreq[rp]
                srvs = staged_next[rp]
                bktrcp = get!(bktrc, rp) do; ResolveBacktraceItem(); end
                push!(bktrcp, p=>bktrc[p], srvs)
                if isa(bktrcp.versionreq, VersionSet) && isempty(bktrcp.versionreq)
                    err_msg = "Unsatisfiable requirements detected for package $rp:\n"
                    err_msg *= string(bktrcp)
                    err_msg *= "The intersection of the requirements is empty."
                    throw(PkgError(err_msg))
                end
            end
        end
        staged = staged_next
    end

    filtered_deps = Dict{String,Dict{VersionNumber,Available}}()
    for (p,depsp) in deps
        filtered_deps[p] = Dict{VersionNumber,Available}()
        allowedp = get(allowed, p) do; Dict{VersionNumber,Bool}() end
        fdepsp = filtered_deps[p]
        for (vn,a) in depsp
            get(allowedp, vn, true) || continue
            fdepsp[vn] = a
        end
    end

    return filtered_deps, allowed
end

# Reduce the number of versions by creating equivalence classes, and retaining
# only the highest version for each equivalence class.
# Two versions are equivalent if:
#   1) They appear together as dependecies of another package (i.e. for each
#      dependency relation, they are both required or both not required)
#   2) They have the same dependencies
# Preliminarily calls filter_versions.
function prune_versions(reqs::Requires, deps::Dict{String,Dict{VersionNumber,Available}}, bktrc::ResolveBacktrace)
    filtered_deps, allowed = filter_versions(reqs, deps, bktrc)
    if !isempty(reqs)
        filtered_deps = dependencies_subset(filtered_deps, Set{String}(keys(reqs)))
    end

    # To each version in each package, we associate a BitVector.
    # It is going to hold a pattern such that all versions with
    # the same pattern are equivalent.
    vmask = Dict{String,Dict{VersionNumber, BitVector}}()

    # For each package, we examine the dependencies of its versions
    # and put together those which are equal.
    # While we're at it, we also collect all dependencies into alldeps
    alldeps = Dict{String,Set{VersionSet}}()
    for (p,fdepsp) in filtered_deps
        # Extract unique dependencies lists (aka classes), thereby
        # assigning an index to each class.
        uniqdepssets = unique(a.requires for a in values(fdepsp))

        # Store all dependencies seen so far for later use
        for r in uniqdepssets, (rp,rvs) in r
            get!(alldeps, rp) do; Set{VersionSet}() end
            push!(alldeps[rp], rvs)
        end

        # If the package has just one version, it's uninteresting
        length(deps[p]) == 1 && continue

        # Grow the pattern by the number of classes
        luds = length(uniqdepssets)
        @assert !haskey(vmask, p)
        vmask[p] = Dict{VersionNumber,BitVector}()
        vmaskp = vmask[p]
        for vn in keys(fdepsp)
            vmaskp[vn] = falses(luds)
        end
        for (vn,a) in fdepsp
            vmind = findfirst(isequal(a.requires), uniqdepssets)
            @assert vmind !== nothing
            vm = vmaskp[vn]
            vm[vmind] = true
        end
    end

    # Produce dependency patterns.
    for (p,vss) in alldeps, vs in vss
        # packages with just one version, or dependencies
        # which do not distiguish between versions, are not
        # interesting
        (length(deps[p]) == 1 || vs == VersionSet()) && continue

        # Store the dependency info in the patterns
        @assert haskey(vmask, p)
        for (vn,vm) in vmask[p]
            push!(vm, vn in vs)
        end
    end

    # At this point, the vmask patterns are computed. We divide them into
    # classes so that we can keep just one version for each class.
    pruned_vers = Dict{String,Vector{VersionNumber}}()
    eq_classes = Dict{String,Dict{VersionNumber,Vector{VersionNumber}}}()
    for (p, vmaskp) in vmask
        vmask0_uniq = unique(values(vmaskp))
        nc = length(vmask0_uniq)
        classes = [VersionNumber[] for c0 = 1:nc]
        for (vn,vm) in vmaskp
            c0 = findfirst(isequal(vm), vmask0_uniq)::Int
            push!(classes[c0], vn)
        end
        map(sort!, classes)

        # For each nonempty class, we store only the highest version)
        pruned_vers[p] = VersionNumber[]
        prunedp = pruned_vers[p]
        eq_classes[p] = Dict{VersionNumber,Vector{VersionNumber}}()
        eqclassp = eq_classes[p]
        for cl in classes
            if !isempty(cl)
                vtop = maximum(cl)
                push!(prunedp, vtop)
                @assert !haskey(eqclassp, vtop)
                eqclassp[vtop] = cl
            end
        end
        sort!(prunedp)
    end
    # Put non-allowed versions into eq_classes
    for (p, allowedp) in allowed
        haskey(eq_classes, p) || continue
        eqclassp = eq_classes[p]
        for (vn, a) in allowedp
            a && continue
            eqclassp[vn] = [vn]
        end
    end
    # Put all remaining packages into eq_classes
    for (p, depsp) in deps
        haskey(eq_classes, p) && continue
        eq_classes[p] = Dict{VersionNumber,Vector{VersionNumber}}()
        eqclassp = eq_classes[p]
        for vn in keys(depsp)
            eqclassp[vn] = [vn]
        end
    end

    # Recompute deps. We could simplify them, but it's not worth it
    new_deps = Dict{String,Dict{VersionNumber,Available}}()

    for (p,depsp) in filtered_deps
        @assert !haskey(new_deps, p)
        if !haskey(pruned_vers, p)
            new_deps[p] = depsp
            continue
        end
        new_deps[p] = Dict{VersionNumber,Available}()
        pruned_versp = pruned_vers[p]
        for (vn,a) in depsp
            vn ∈ pruned_versp || continue
            new_deps[p][vn] = a
        end
    end

    #println("pruning stats:")
    #numvers = 0
    #numdeps = 0
    #for (p,d) in deps, (vn,a) in d
    #    numvers += 1
    #    for r in a.requires
    #        numdeps += 1
    #    end
    #end
    #numnewvers = 0
    #numnewdeps = 0
    #for (p,d) in new_deps, (vn,a) in d
    #    numnewvers += 1
    #    for r in a.requires
    #        numnewdeps += 1
    #    end
    #end
    #println("  before: vers=$numvers deps=$numdeps")
    #println("  after: vers=$numnewvers deps=$numnewdeps")
    #println()

    return new_deps, eq_classes
end
prune_versions(deps::Dict{String,Dict{VersionNumber,Available}}) =
    prune_versions(Dict{String,VersionSet}(), deps, ResolveBacktrace())
prune_versions(deps::Dict{String,Dict{VersionNumber,Available}}, bktrc::ResolveBacktrace) =
    prune_versions(Dict{String,VersionSet}(), deps, bktrc)

# Build a graph restricted to a subset of the packages
function subdeps(deps::Dict{String,Dict{VersionNumber,Available}}, pkgs::Set{String})
    sub_deps = Dict{String,Dict{VersionNumber,Available}}()
    for p in pkgs
        haskey(sub_deps, p) || (sub_deps[p] = Dict{VersionNumber,Available}())
        sub_depsp = sub_deps[p]
        for (vn,a) in deps[p]
            sub_depsp[vn] = a
        end
    end

    return sub_deps
end

# Build a subgraph incuding only the (direct and indirect) dependencies
# of a given package set
function dependencies_subset(deps::Dict{String,Dict{VersionNumber,Available}}, pkgs::Set{String})
    staged::Set{String} = filter(p->p in keys(deps), pkgs)
    allpkgs = copy(staged)
    while !isempty(staged)
        staged_next = Set{String}()
        for p in staged, a in values(get(deps, p, Dict{VersionNumber,Available}())), rp in keys(a.requires)
            rp ∉ allpkgs && rp ≠ "julia" && push!(staged_next, rp)
        end
        union!(allpkgs, staged_next)
        staged = staged_next
    end

    return subdeps(deps, allpkgs)
end

# Build a subgraph incuding only the (direct and indirect) dependencies and dependants
# of a given package set
function undirected_dependencies_subset(deps::Dict{String,Dict{VersionNumber,Available}}, pkgs::Set{String})
    graph = Dict{String, Set{String}}()

    for (p,d) in deps
        haskey(graph, p) || (graph[p] = Set{String}())
        for a in values(d), rp in keys(a.requires)
            push!(graph[p], rp)
            haskey(graph, rp) || (graph[rp] = Set{String}())
            push!(graph[rp], p)
        end
    end

    staged = pkgs
    allpkgs = copy(pkgs)
    while !isempty(staged)
        staged_next = Set{String}()
        for p in staged, rp in graph[p]
            rp ∉ allpkgs && push!(staged_next, rp)
        end
        union!(allpkgs, staged_next)
        staged = staged_next
    end

    return subdeps(deps, allpkgs)
end

function prune_dependencies(reqs::Requires,
                            deps::Dict{String,Dict{VersionNumber,Available}},
                            bktrc::ResolveBacktrace = init_resolve_backtrace(reqs))
    deps, _ = prune_versions(reqs, deps, bktrc)
    return deps
end

end # module

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
# Usage:
#   julia run_benchmarks.jl [<#runs> <#threads>]

using Statistics
using Serialization
using Printf

const RUNS = isempty(ARGS) ? 10 : parse(Int, ARGS[1])
const THREADS = isempty(ARGS) ? 1 : parse(Int, ARGS[2])
const JULIAVER = Base.julia_cmd()[1]
dir = joinpath(@__DIR__, "benches")
cd(dir)

gctime(stat) = stat.total_time

for category in readdir()
    @show category
    cd(category)
    for test in readdir()
        endswith(test, ".jl") || continue
        @show test
        value = []
        times = []
        stats = []
        for _ in 1:RUNS
            r = open(deserialize, `$JULIAVER --project=. --threads=$THREADS $test SERIALIZE`)
            push!(value, r.value)
            push!(times, r.times)
            push!(stats, r.stats)
        end
        @printf("run time: %0.0fms min, %0.0fms max %0.0fms median\n",
           minimum(times) / 1_000_000,
           maximum(times) / 1_000_000,
           median(times)  / 1_000_000)
        time = map(gctime, stats)
        @printf("gc time: %0.0fms min, %0.0fms max, %0.0fms median\n",
           minimum(time) / 1_000_000,
           maximum(time) / 1_000_000,
           median(time)  / 1_000_000)
    end
    cd("..")
end

max_stars_repo_name: oscardssmith/GCBenchmarks

max_issues_repo_name: oscardssmith/GCBenchmarks

max_forks_repo_name: oscardssmith/GCBenchmarks


text:
using Xpress, JuMP, MathProgBase

include(Pkg.dir("JuMP","test","solvers.jl"))
lp_solvers = [XpressSolver()]
ip_solvers = [XpressSolver()]

include(Pkg.dir("JuMP","test","model.jl"))
include(Pkg.dir("JuMP","test","probmod.jl"))

FactCheck.exitstatus()


text:
export BotClient, BearerClient

struct Token
    token::String
end

Base.show(io::IO, ::Token) = print(io, "<token>")

struct BotClient
    token::Token
    rate_limiter::RateLimiter

    BotClient(token; rate_limiter=RateLimiter()) = new(Token(token), rate_limiter)
end

struct BearerClient
    token::Token
    rate_limiter::RateLimiter

    BearerClient(token; rate_limiter=RateLimiter()) = new(Token(token), rate_limiter)
end

auth_header(c::BotClient) = "Bot $(c.token.token)"
auth_header(c::BearerClient) = "Bearer $(c.token.token)"

rate_limiter(c::BotClient) = c.rate_limiter
rate_limiter(c::BearerClient) = c.rate_limiter


text:
using TreeSitterHighlight
using Markdown

const treesitter_javascript =
    "~/Projects/tree-sitter-javascript/build_so/libtreesitter_javascript.so" |> expanduser

const treesitter_julia =
    "~/Projects/tree-sitter-julia/build_so/libtreesitter_julia.so" |> expanduser

captures = [
    "symbol",
    "include",
    "variable",
    "comment",
    "tag",
    "function",
    "string",
    "keyword",
    "punctuation",
]

highlighter = Highlighter(Dict(capture => "class=$capture" for capture in captures))

js_scope = "source.js"
js_injection_regex = "^javascript"

js = Language(
    :javascript,
    @ccall treesitter_javascript.tree_sitter_javascript()::Ptr{Nothing}
)

js_highlights_query = read("./queries/javascript/highlights.scm", String)
js_injections_query = read("./queries/javascript/injections.scm", String)
js_locals_query = read("./queries/javascript/locals.scm", String)

add_language!(
    highlighter,
    js;
    scope = js_scope,
    injection_regex = js_injection_regex,
    highlights_query = js_highlights_query,
    injections_query = js_injections_query,
    locals_query = js_locals_query,
)

jl_scope = "source.jl"
jl_injection_regex = "^julia"

jl = Language(:julia, @ccall treesitter_julia.tree_sitter_julia()::Ptr{Nothing})

jl_highlights_query = read("./queries/julia/highlights.scm", String)
jl_injection_query = read("./queries/julia/injections.scm", String)
jl_locals_query = read("./queries/julia/locals.scm", String)

add_language!(
    highlighter,
    jl;
    scope = jl_scope,
    injection_regex = jl_injection_regex,
    highlights_query = jl_highlights_query,
    injections_query = jl_injection_query,
    locals_query = jl_locals_query,
)

scopes = Dict{String,String}("javascript" => js_scope, "julia" => jl_scope)

function Markdown.html(io::IO, code::Markdown.Code)
    if code.language ∈ ("julia", "javascript")
        write(io, "<pre>")
        write(io, highlight(highlighter, code.code; scope = scopes[code.language]))
        write(io, "</pre>")
    else
        write(io, code.code)
    end
end

readme = read("./README.md", String) |> Markdown.parse

generate_index() =
    open("index.html", "w") do f
        write(
            f,
            """
   <html>
   <head>
   <style>
   body {
       margin: 0;
       padding: 0;
       font-family: -apple-system,BlinkMacSystemFont,"Segoe UI",Helvetica,Arial,sans-serif;
   }

   h1, h2, h3, p {
       margin: 0;
       padding-top: 5px;
       padding-left: 20px;
       padding-right: 20px;
   }

   pre {
       padding: 20px;
       background: hsla(46, 90%, 98%, 1);
       color: #41323f;
   }

   pre p {
       white-space: pre-wrap;
       margin: 0;
   }

   .symbol, .symbol * {
       color: #815ba4 !important;
   }

   .comment {
       color: #e96ba8;
   }

   .keyword, .include {
       color: #ef6155;
   }

   .string {
       color: #da5616;
   }

   .function, .function * {
       color: #cc80ac !important;
   }

   .number {
       color: #815ba4;
   }

   .punctuation {
       color: #41323f;
   }

   .variable {
       color: #5668a4;
   }
   </style>
   </head>
   <body>
   """,
        )

        show(f, MIME"text/html"(), readme)

        write(
            f,
            """
   </body>
   </html>
   """,
        )
    end

generate_index()


text:
using NiLang, NiLang.AD

@i function iexp(out!, x::T; atol::Float64=1e-14) where T
    @anc anc1 = zero(T)
    @anc anc2 = zero(T)
    @anc anc3 = zero(T)
    @anc iplus = 0
    @anc expout = zero(T)

    out! += identity(1.0)
    @routine r1 begin
        anc1 += identity(1.0)
        while (value(anc1) > atol, iplus != 0)
            iplus += identity(1)
            anc2 += anc1 * x
            anc3 += anc2 / iplus
            expout += identity(anc3)
            # speudo inverse
            anc1 -= anc2 / x
            anc2 -= anc3 * iplus
            SWAP(anc1, anc3)
        end
    end

    out! += identity(expout)

    ~@routine r1
end

using Test
# NOTE: to allow high performance use of f += T(x).
# Now, this kind of use is a performance killer.
@testset "iexp" begin
    out = 0.0
    x = 1.3
    @instr iexp(out, x)
    res = exp(x)
    @test check_inv(iexp,(out, x))
    @test out ≈ res

    out = 0.0
    x = 1e-9
    @instr iexp(out, x)
    res = exp(x)
    @test check_inv(iexp,(out, x))
    @test out ≈ res

    out = 0.0
    x = 1.0

    @instr iexp(out, x)
    res = exp(x)
    @test check_inv(iexp,(out, x))
    @test out ≈ res
end

@testset "iexp grad" begin
    NiLangCore.GLOBAL_ATOL[] = 1e-2
    out = 0.0
    x = 1.6
    gres = exp(x)
    @test check_inv(iexp, (out, x); verbose=true)
    @test check_grad(iexp, (Loss(out), x); verbose=true)

    out = 0.0
    x = 1.6
    @instr iexp'(Loss(out), x)
    @test grad(x) ≈ gres

    h1 = (iexp''(Loss(0.0), 1.6); collect_hessian())
    h2 = simple_hessian(iexp, (Loss(0.0), 1.6))
    nh = nhessian(iexp, (Loss(0.0), 1.6))
    @show h1, h2
    @test isapprox(h1, nh, atol=1e-3)
    @test isapprox(h2, nh, atol=1e-3)
end


text:
# the structure is mutable for the sake of optimization and internal needs. don't mutate your phantoms
mutable struct AetherlangPhantomCollection{T}
    """A phantom structure for basic collections"""
    ref::T
    addindx::Dict{Integer, Any}
    length::Int64
end

function Base.getindex(p::AetherlangPhantomCollection, idx::Integer)
    try
        return p.addindx[idx]
    catch
        return p.ref[idx]
    end
end

Base.length(p::AetherlangPhantomCollection) = p.length

Base.copy(p::AetherlangPhantomCollection) = AetherlangPhantomCollection(p.ref, copy(p.addindx), p.length)

function aetherlang_show(io::IO, v::AetherlangPhantomCollection)
    print(io, "(phantomarray "*join([aetherlang_repr(v[i]) for i in 1:length(v)], ' ')*")")
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
module Dircrypt

function rand_modulus(ix, modulus)
    ix = 16807*(ix % 127773) - 2836*trunc(Int64, ix / 127773) & 0xFFFFFFFF
    return ix, ix % modulus
end

function gen_one(seed="1DBA8930", tld="")
    ix = trunc(Int64, rand()*parse(Int64, seed, base=16))
    domain = ""
    ix, len_domain = rand_modulus(ix, 13)
    len_domain += 8
    for _ in 1:len_domain
        ix, tmp = rand_modulus(ix, 26)
        domain *= 'a'+tmp
    end
    domain *= tld

    return domain
end

function gen_many(cnt=100)
    res = []
    for _ in 1:cnt
        domain = gen_one()
        res = [res..., domain]
    end
    return res
end

end  # module Dircrypt


text:
export EuclideanVectorState, DiscreteState, ACEConfig

@doc raw"""
`struct EuclideanVectorState` : a $\mathbb{R}^3$-vector, which transforms under
the rotation group as
```math
g_Q \cdot {\bm r} = Q {\bm r}.
```
It typically defines a position or a force.
"""
struct EuclideanVectorState{T} <: AbstractContinuousState
   rr::SVector{3, T}
   label::String
end

EuclideanVectorState(rr::SVector{3}) = EuclideanVectorState(rr, "𝒓")
EuclideanVectorState{T}(label::String = "𝒓") where {T} = EuclideanVectorState(zero(SVector{3, T}), label)
EuclideanVectorState(label::String = "𝒓") = EuclideanVectorState(zero(SVector{3, Float64}), label)

Base.length(X::EuclideanVectorState) = 3

Base.show(io::IO, s::EuclideanVectorState) = print(io, "$(s.label)$(s.rr)")

import Base: *, -, +
*(A::Union{Number, AbstractMatrix}, X::EuclideanVectorState) =
      EuclideanVectorState(A * X.rr, X.label)
-(X::T) where {T <: EuclideanVectorState} = T( -X.rr, X.label)
+(X1::T, X2::T) where {T <: EuclideanVectorState} =
      T(X1.rr + X2.rr, X1.label)
+(X1::T, u::SVector{3}) where {T <: EuclideanVectorState} =
      T(X1.rr + u, X1.label)

@doc raw"""
`struct DiscreteState` : a state ``\mu`` specified by a discrete number of
possible values, e.g. ranging through ``\mathbb{Z}`` or ``\mathbb{Z}_p``.
Discrete states cannot possibly have a non-trivial transformation under
a continuoue group action, hence it is assumed to be invariant, i.e.,
```math
g_Q \cdot \mu = \mu
```
"""
struct DiscreteState{T, SYM} <: AbstractDiscreteState
   val::T
   _valsym::Val{SYM}
end

DiscreteState(sym::Symbol) = DiscreteState(0, Val(sym))
DiscreteState{T, SYM}(val::T) where {T, SYM} = DiscreteState(val, Val{SYM}())

Base.show(io::IO, s::DiscreteState{T, SYM}) where {T, SYM} =
         print(io, "$(SYM)[$(s.val)]")

Base.getproperty(s::DiscreteState{T, SYM}, sym) where {T, SYM} =
      sym == SYM ? getfield(s, :val) : getfield(s, sym)

struct ACEConfig{STT} <: AbstractConfiguration
   Xs::Vector{STT}   # list of states
end

# --- iterator to go through all states in an abstract configuration assuming
#     that the states are stored in cfg.Xs

Base.iterate(cfg::AbstractConfiguration) =
   length(cfg.Xs) == 0 ? nothing : (cfg.Xs[1], 1)

Base.iterate(cfg::AbstractConfiguration, i::Integer) =
   length(cfg.Xs) == i ? nothing : (cfg.Xs[i+1], i+1)

Base.length(cfg::AbstractConfiguration) = length(cfg.Xs)


text:
module SnakePlot

using DataFrames
using Statistics
using GLMakie

include("readFasta.jl")
include("fasta2df.jl")
include("assignStructure.jl")
include("positionTMs!.jl")
include("loopLayout!.jl")
include("positionLoops!.jl")
include("bundleTMs!.jl")

export readFasta, fasta2df, assignStructure, positionTMs!, positionLoops!, bundleTMs!
end


text:
using RecursiveLeastSquares
using Test

@testset "RecursiveLeastSquares.jl" begin
    # Write your tests here.
end


text:
using ComputationalResources
using Test

@testset "ComputationalResources" begin
    @testset "Settings" begin
        for r in (CPU1(5), CPUThreads(5), CPUProcesses(5), ArrayFireLibs(5), CUDALibs(5), OpenCLLibs(5))
            for T in (CPU1, CPUThreads, CPUProcesses, ArrayFireLibs, CUDALibs, OpenCLLibs)
                r1 = T(r)
                @test r1.settings === 5
            end
        end

        @test isa(CPUThreads(), AbstractCPU{Nothing})
        @test isa(CPUProcesses(), AbstractCPU{Nothing})
        @test isa(CUDALibs(), AbstractResource{Nothing})
        @test isa(OpenCLLibs(), AbstractResource{Nothing})
    end

    @testset "Dummy1" begin
        push!(LOAD_PATH, joinpath(dirname(@__FILE__), "packages"))

        @eval import Dummy1

        @test Dummy1.foo(CPU1())
        @test_throws MethodError Dummy1.foo(ArrayFireLibs())

        addresource(ArrayFireLibs)
        @test haveresource(ArrayFireLibs)

        Main.Dummy1.__init__()

        @test Main.Dummy1.foo(CPU1())
        @test Main.Dummy1.foo(ArrayFireLibs())

        rmresource(ArrayFireLibs)

        include("packages/Dummy1/src/Dummy1.jl")

        @test Main.Dummy1.foo(CPU1())
        @test_throws MethodError Main.Dummy1.foo(ArrayFireLibs())

        pop!(LOAD_PATH)
    end
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

path = "examples/examples_files/doublependulum.urdf"
mech, shapes = Mechanism(path)

storage = simulate!(mech, 10., record = true)
visualize(mech, storage, shapes)


text:
# Kalman Filter tests

# Test construction with absolute state
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = DiscreteState([1.0, 2.0])
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
@test typeof(kf.estimate) <: UncertainDiscreteState
@test kf.estimate.P == zeros(2,2)

# Test consider constructors
kf = KalmanFilter(linear_sys, linear_obs, initial_est, [2])
kf = KalmanFilter(linear_sys, linear_obs, make_uncertain(initial_est), [2])
@test_throws DimensionMismatch KalmanFilter(linear_sys, linear_obs, initial_est, [3])
@test_throws ArgumentError KalmanFilter(linear_sys, linear_obs, initial_est, [3,3])

# Test discrete kalman filter
srand(2);
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0], 0.1*eye(2))
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
for i = 1:10
    true_state, measurement = simulate(simulator, i)
    process!(kf, measurement)
    @test mahalanobis(true_state, kf.estimate) < 3
end

# Test discrete consider kalman filter
srand(2);
A = [1.0 0.5 0.01; -0.5 1.0 0.01; 0.0 0.0 1.0]
linear_sys = LinearSystem(A, 0.001*eye(3))
linear_obs = LinearObserver(eye(2,3), 0.1*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0, 3.0], 0.1*eye(3))
consider_states = [3]
kf = KalmanFilter(linear_sys, linear_obs, initial_est, consider_states)
simulator = make_simulator(kf)
for i = 1:10
    true_state, measurement = simulate(simulator, i)
    process!(kf, measurement)
    @test mahalanobis(true_state, kf.estimate) < 3
end
@test kf.estimate.P[3,3] == initial_est.P[3,3]

# Test discrete kalman filter with archiving
srand(2);
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.01*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0], 0.1*eye(2))
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
archive = EstimatorHistory(kf)
for i = 1:10
    true_state, measurement = simulate(simulator, i)
    process!(kf, measurement, archive)
    @test mahalanobis(true_state, kf.estimate) < 3
end

# Test continuous kalman filter
srand(2);
linear_sys = LinearSystem(-eye(3), 0.001*eye(3))
linear_obs = LinearObserver(eye(3), 0.01*eye(3))
initial_est = UncertainContinuousState([1.0, 2.0, 3.0], 0.1*eye(3))
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
for i = 1:10
    true_state, measurement = simulate(simulator, kf.estimate.t+0.1)
    process!(kf, measurement)
    @test mahalanobis(true_state, kf.estimate) < 3
end

# Test continuous kalman filter with archiving
srand(2);
linear_sys = LinearSystem(-eye(3), 0.001*eye(3))
linear_obs = LinearObserver(eye(3), 0.01*eye(3))
initial_est = UncertainContinuousState([1.0, 2.0, 3.0], 0.1*eye(3))
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
archive = EstimatorHistory(kf)
for i = 1:10
    true_state, measurement = simulate(simulator, kf.estimate.t+0.1)
    process!(kf, measurement, archive)
    @test mahalanobis(true_state, kf.estimate) < 3
end

# Test continuous kalman filter convenience function
srand(2);
linear_sys = LinearSystem(-eye(3), 0.001*eye(3))
linear_obs = LinearObserver(eye(3), 0.001*eye(3))
estimate = UncertainContinuousState([1.0, 2.0, 3.0], 0.1*eye(3))
simulator = make_simulator(KalmanFilter(linear_sys, linear_obs, estimate))
for i = 1:5
    true_state, measurement = simulate(simulator, kf.estimate.t+0.1)
    estimate = kalman_filter(estimate, linear_sys, linear_obs, measurement)
end

# Test continuous kalman filter with archiving
srand(2);
linear_sys = LinearSystem(-eye(3), 0.001*eye(3))
linear_obs = LinearObserver(eye(3), 0.001*eye(3))
estimate = UncertainContinuousState([1.0, 2.0, 3.0], 0.1*eye(3))
kf = KalmanFilter(linear_sys, linear_obs, estimate)
simulator = make_simulator(kf)
archive = EstimatorHistory(kf)
for i = 1:5
    true_state, measurement = simulate(simulator, kf.estimate.t+0.1)
    estimate =
        kalman_filter(estimate, linear_sys, linear_obs, measurement, archive)
end


text:
# susing Revise
# using Profile
# using JLD2
using CSV
# using DataFrames
using DecisionTrees

# include(joinpath("..","test\\runtests.jl"))

datafile = string("data\\freMTPL2\\freMTPL2.csv")
@assert isfile(datafile);
@time fullData = CSV.read(datafile, rows_for_type_detect=100000, allowmissing=:none, categorical=false);

areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]

dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, trnvalcol="trnTest", numcol="ClaimNb", denomcol="Exposure", weightcol="Exposure", independent_vars=selected_explanatory_vars);

updateSettingsMod!(sett,minWeight=-0.03,model_type="boosted_tree",iterations=40,learningRate=0.025,subsampling_features_prop=.7,calculatePoissonError=true)

resultingFiles, resM = dtm(dtmtable, sett) # 31s on notebook (x260 with profiling disabled but Revise on!)

##
# consider a boosting model
##
import DecisionTrees: calculateSplitValue,increasing_subsets,build_listOfMeanResponse,_split_feature,some_tree_settings,build_tree!,global_statsperiter_header,createScorebandsUTF8List,initBootstrapSample,_split,sample_data_and_build_tree!,find_max_type,emptyNode,ExcelData,Tree,Leaf,get_feature_pools,build_tree_iteration!,ExcelSheet,Chart

trnidx = dtmtable.trnidx
validx = dtmtable.validx
features = dtmtable.features
weight = dtmtable.weight
actualNumerator = dtmtable.numerator
denominatorDTM = dtmtable.denominator
mappings = dtmtable.mappings
candMatWOMaxValues = dtmtable.candMatWOMaxValues
obstrn, obsval, trn_meanobservedvalue, val_meanobservedvalue, trn_numtot, val_numtot, trn_denomtot, val_denomtot, empty_rows_after_iteration_stats, showTimeUsedByEachIteration, chosen_apply_tree_fn, startAtMean, adaptiveLearningRate, moderationvector, iterations, nameOflistofpredictorsSheet, nameOfpredictorsSheet, nameOfModelStatisticsSheet, nameOfScoresSheet, nameOfOverviewSheet, nameOfSettingsSheet, nameOfValidationSheet, xlData, showProgressBar_time, prroduceEstAndLeafMatrices = DecisionTrees.initSettingsWhichAreTheSameForBoostingAndBagging(trnidx, validx, actualNumerator, denominatorDTM, sett)
    # trn_meanobservedvalue is the mean observed RATIO, i.e. sum(numerator)/sum(denominatorDTM) for the training data
obs = obstrn + obsval
current_error = 0.0
moderationfactor = moderationvector[1]
if (size(moderationvector, 1) != iterations && size(moderationvector, 1) != 1)
    error("DTM: This should not have happend: size(mdf_vector) is neither 1 nor #iterations")
end
if size(moderationvector, 1) == 1
    actual_moderationvector = Array{Float64}(undef, iterations)
    fill!(actual_moderationvector, deepcopy(moderationvector[1]))
else
    actual_moderationvector = deepcopy(moderationvector)
end

if startAtMean
    estimatedRatio = ones(length(weight)) .* trn_meanobservedvalue
else
    error("DTM: This is currently not working: please set startAtMean=true)")
end

estimatedNumerator = estimatedRatio .* denominatorDTM # these are for instance the estimated losses for a LR model

indicatedRelativityForApplyTree_reused = zeros(obs)
reused_fitted_leafnr_vector = zeros(Int, obs)
sortvec_reused_trn_only = zeros(Int, length(trnidx))
intVarsUsed = Array{Array{Array{Int,1},1}}(undef, 0);sizehint!(intVarsUsed, iterations)
inds_considered = Array{Array{Int,1}}(undef, 0);sizehint!(inds_considered, iterations)

res = Array{Union{Leaf,Node{UInt8},Node{UInt16}}}(undef, iterations)
if prroduceEstAndLeafMatrices
    est_matrix = Array{Float64}(undef, obs, iterations + 1)
    est_matrixFromScores = deepcopy(est_matrix)
    MatrixOfLeafNumbers = Array{Int}(undef, obs, iterations + 1)
    MatrixOfLeafNumbers[:,1] .= 0
    est_matrix[:,1] = copy(transpose(estimatedRatio))
    est_matrixFromScores[:,1] = copy(transpose(estimatedRatio))
else
    est_matrix = Array{Float64}(undef, 0, 0)
    est_matrixFromScores = deepcopy(est_matrix)
    MatrixOfLeafNumbers = Array{Int}(undef, 0, 0)
end

scores = zeros(Int, obs)
estFromScores = zeros(obs)
    # vectorOfRulePathsToLeavesArrays=Array{Array{Array{Rulepath,1},1}}(iterations+1)
    # vectorOfRulePathsToLeavesArrays[1]=Array{Array{Rulepath,1}}(undef,0) #the first entry is not defined
vectorOfLeafArrays = Array{Array{Leaf,1}}(undef, iterations + 1)
vectorOfLeafArrays[1] = Array{Leaf}(undef, 0) # the first entry is not defined

p = ProgressMeter.Progress(iterations, 2, "Progress of Boosting Model:") # minimum update interval: x seconds (2)
((adaptiveLearningRate >= 1.0) || (adaptiveLearningRate <= 0.0)) ? (BoolAdaptiveLearningRate = false) : (BoolAdaptiveLearningRate = true)
    # currently not supported
BoolAdaptiveLearningRate && error("DTM: this is not yet updated to reflect numerator/denominatorDTM objectives!!")
    # this header needs to be in line with the output of the function createTrnValStatsForThisIteration
statsPerIteration = global_statsperiter_header
currentRelativity = Array{Float64}(undef, obs)
scoreBandLabels = createScorebandsUTF8List(sett.scorebandsstartingpoints, sett.nScores, addtotal=true)
    # this row is only here such that the variables are initialized outside the for loop
local fristRowOfThisTable, obsPerScore, vectorWeightPerScore, numPerScore, denomPerScore, scores, scoresVAL, cumulativeStatsPerScoreBand, maxRawRelativityPerScoreSorted, rawObservedRatioPerScore, MAPPINGSmoothedEstimatePerScore

    # the following two lines are for the initialization in case any subsampling (boostrapping) of data is done
    # the idea is to preallocate many of the arrays which are used for sampling (and to save allocation/memory space)
sampleSizeCanBeNEGATIVE = convert(Int, round(sett.subsampling_prop * obstrn))
if abs(sampleSizeCanBeNEGATIVE) < 1
    sampleSizeCanBeNEGATIVE = convert(Int, sign(sett.subsampling_prop))
end
abssampleSize, sampleVector = initBootstrapSample(sampleSizeCanBeNEGATIVE) # ,trnidx,validx)

    # for iter=1:iterations
        # @timeConditional(showTimeUsedByEachIteration, begin
        # Build Iteration iter
            # showTimeUsedByEachIteration&&println("Calculating iteration $(iter). Time: $(Dates.now())")
iter = 1

current_mdf = moderationvector[min(iter, size(moderationvector, 1))] # mdf remains constant if a vector of size 1 was provided (instead of a vector of size iterations)

T_Uint8_or_UInt16 = find_max_type(features)

tmpTree = sample_data_and_build_tree!(trnidx, validx, indicatedRelativityForApplyTree_reused, candMatWOMaxValues, mappings, deepcopy(sett), actualNumerator, estimatedNumerator, weight, features, sampleSizeCanBeNEGATIVE, abssampleSize, sampleVector, T_Uint8_or_UInt16)

using BenchmarkTools

@code_warntype sample_data_and_build_tree!(trnidx, validx, indicatedRelativityForApplyTree_reused, candMatWOMaxValues, mappings, deepcopy(sett), actualNumerator, estimatedNumerator, weight, features, sampleSizeCanBeNEGATIVE, abssampleSize, sampleVector, T_Uint8_or_UInt16)

# if no subsampling is done:
nnnode = build_tree!(trnidx, validx, candMatWOMaxValues, mappings, sett, actualNumerator, estimatedNumerator, weight, features, indicatedRelativityForApplyTree_reused, T_Uint8_or_UInt16)

# if no subsampling is done:
@code_warntype build_tree!(trnidx, validx, candMatWOMaxValues, mappings, sett, actualNumerator, estimatedNumerator, weight, features, indicatedRelativityForApplyTree_reused, T_Uint8_or_UInt16)

settings = sett
fitted_values = indicatedRelativityForApplyTree_reused

intVarsUsed, inds, minweightcalculated = some_tree_settings(trnidx, validx, settings.fixedinds, candMatWOMaxValues, mappings, settings.minWeight, weight, settings.subsampling_features_prop, size(features, 2))

	settings.minWeight = minweightcalculated # update minWeight
	if !(length(inds) > 0)
    throw(ErrorException(string("Error: no features were selected length(inds)=", length(inds))))
end
	empty_xl_data = ExcelData(Array{ExcelSheet}(undef, 0), Array{Chart}(undef, 0))
	fp = get_feature_pools(features)
	resultingTree = Tree(emptyNode, intVarsUsed, candMatWOMaxValues, mappings, inds, settings, empty_xl_data, fp)
	resultingTree.rootnode = build_tree_iteration!(trnidx, validx, settings, resultingTree, actualNumerator, estimatedNumerator, weight, features, 0, settings.randomw, Array{Rulepath}(undef, 0), Distributed.myid(), fitted_values, T_Uint8_or_UInt16)

@time resultingTree.rootnode = build_tree_iteration!(trnidx, validx, settings, resultingTree, actualNumerator, estimatedNumerator, weight, features, 0, settings.randomw, Array{Rulepath}(undef, 0), Distributed.myid(), fitted_values, T_Uint8_or_UInt16)

@code_warntype build_tree_iteration!(trnidx, validx, settings, resultingTree, actualNumerator, estimatedNumerator, weight, features, 0, settings.randomw, Array{Rulepath}(undef, 0), Distributed.myid(), fitted_values, T_Uint8_or_UInt16)

    ##

thisTree = resultingTree

    	# !!!! the current concept forsees that features is always the FULL DataFrame
	boolRandomizeOnlySplitAtTopNode = settings.boolRandomizeOnlySplitAtTopNode
	local inds
	# if fixedinds are provided we never change the selection of features for the whole tree! (bagged boosting model)
	if boolRandomizeOnlySplitAtTopNode || (length(settings.fixedinds) > 0)
		# these settings are randomly chosen by some_tree_settings in this case
		inds = thisTree.inds_considered
	else
		# inds need to be randomly picked for each iteration
		inds = randomFeatureSelection(size(features, 2), settings.subsampling_features_prop)
	end
	intVarsUsed = thisTree.intVarsUsed

	minweight = settings.minWeight
	crit = settings.crit
	spawnsmaller = settings.spawnsmaller
	catSortByThreshold = settings.catSortByThreshold
	catSortBy = settings.catSortBy

	nobs = size(actualNumerator, 1)
	fnames = names(features)
	# @code_warntype _split(settings.number_of_num_features,trnidx,validx,numerator,denominator,weight,fnames,features, minweight, depth,randomweight,crit,parallel_level_threshold,parallel_weight_threshold,inds,catSortByThreshold,catSortBy)
	# T<:Union{UInt8,UInt16}=find_max_type(features)<:Union{UInt8,UInt16}
    # T=find_max_type(features)::DataType #Union{UInt8,UInt16}

depth = 0
randomweight = sett.randomw

best_split = _split(one(T_Uint8_or_UInt16), settings.number_of_num_features, trnidx, validx, actualNumerator, estimatedNumerator, weight, fnames, features, minweight, depth, randomweight, crit, inds, catSortByThreshold, catSortBy)

@code_warntype _split(one(T_Uint8_or_UInt16), settings.number_of_num_features, trnidx, validx, actualNumerator, estimatedNumerator, weight, fnames, features, minweight, depth, randomweight, crit, inds, catSortByThreshold, catSortBy)

T = T_Uint8_or_UInt16
val_of_some_UInt_type = one(T)
number_of_num_features = sett.number_of_num_features

    # This function selects the maximal possible split defined by crit (thus depending on the impurity function, we need to put a minus sign in front of it)
		tmpsz = 0
		if sum(view(weight, trnidx)) < 2 * minweight;
		#	return const_default_splitdef;
		end
		tmp_splitlist = Vector{Splitdef{T}}(undef, 0)
		# for i in inds
i = inds[1]
			# ATTENTION: for char variables we pass the variable i with a negative sing!!
			# this allows us to distinguish whether we are working on a char or num variable later on
			modified_i = eltype(features[i]) <: AbstractString ? number_of_num_features - i  :  i
			tmplist = _split_feature(val_of_some_UInt_type, number_of_num_features, trnidx, validx, actualNumerator, estimatedNumerator, weight, fnames[i], features[i], minweight, crit, modified_i, randomweight, catSortByThreshold, catSortBy)

@warn("be careful here, features is overwritten")
features = features[i]

crit_type = typeof(crit)
# This function is now for numeric and character variables!
# feature_column_id is negative in case of character variables
best_value = -Inf
best_thresh = best_wl = best_wr = NaN
best_subset = Array{UInt8}(undef, 0)
trnfeatures = view(features, trnidx)
elt = T# elt=eltype(trnfeatures.parent.refs)
labellist_sorted = collect(one(elt):convert(elt, length(trnfeatures.parent.pool))) # this used to be levels(features) #this also contains the val feature levels here! It is considerably faster than levels(view) \factor 100 or so
	# THIS IS CRITICAL
	# THE WAY A PooledArray IS CONSTRUCTED, WE WILL ALWAYS HAVE
	# pda.pool[2] == "some string" -> any pda.refs[x].==0x02 will be "some string"
	# todo/tbd countsort here might be obsolete: we should check if levels is always sorted by construction
  # also we will later sort the labels in a different order anyway (then again the list probably needs to be sorted in the natural manner such that build_listOfMeanResponse is working properly)
if size(labellist_sorted, 1) <= 1
    if T == UInt8
        return UInt8VECTORemptySplitDef
    else
        return UInt16VECTORemptySplitDef # collect(Vector{Splitdef{T}}(undef,0))::Vector{Splitdef{T}}
    end
else
	# countsort!(labellist_sorted)
    # here I am (somewhat) misusing multiple dispatch since I was too lazy to parametrize this at the moment (todo/tbd in the future)
	if (crit_type == DifferenceSplit || crit_type == PoissonDevianceSplit || crit_type == MaxValueSplit || crit_type == MaxMinusValueSplit)
		labellist, sumnumerator, sumdenominator, sumweight, countlistfloat = build_listOfMeanResponse(crit, trnidx, validx, actualNumerator, estimatedNumerator, weight, trnfeatures, labellist_sorted, minweight)

        # @code_warntype build_listOfMeanResponse(crit,trnidx,validx,actualNumerator,estimatedNumerator,weight,trnfeatures,labellist_sorted,minweight)

	elseif (crit_type == msePointwiseSplit)
		labellist, sumnumerator, sumdenominator, sumweight, countlistfloat, moments_per_pdaclass = build_listOfMeanResponse(crit, numerator, denominator, weight, trnfeatures, labellist_sorted, minweight)
	# else #we catch this possibility earlier when checking the settings
	#	throw(ErrorException(string("Invalid Splitting criterion $(crit)")))
	end
  # todo/tbd
  # here we can introduce the possiblity to sort the labellist (e.g. by meanobserved or median (to be calculated)).
  # then we could only loop through the "increasing list" of sorted labels (instead of doing the 2^ncategories exhaustive search (bitflip_graycode_subsets))
    if feature_column_id > 0 # id>0 -> we are working on a numeric column
		# only consider to split at the candidate split points
		subs = increasing_subsets(labellist)
    else # id<0 -> we are working on a character column
		# distinguish between exhaustive and "increasing" search for split point
	if size(labellist_sorted, 1) > catSortByThreshold
		if (crit_type == DifferenceSplit || crit_type == MaxValueSplit || crit_type == MaxMinusValueSplit)
			sortlists!(catSortBy, labellist, sumnumerator, sumdenominator, sumweight, countlistfloat) # catSortBy::SortBy=SORTBYMEAN
		elseif (crit_type == msePointwiseSplit)
			sortlists!(catSortBy, labellist, sumnumerator, sumdenominator, sumweight, countlistfloat, moments_per_pdaclass)
		end
		subs = increasing_subsets(labellist)
	else
	# perform exhaustive search
		subs = bitflip_graycode_subsetsHALF(labellist)
	end
	end

    fname = :hopefully_not_important

	if randomweight == 0.0
	if (crit_type == DifferenceSplit || crit_type == MaxValueSplit || crit_type == MaxMinusValueSplit)
		tmp_result = calculateSplitValue(crit, fname, number_of_num_features, labellist, sumnumerator, sumdenominator, sumweight, countlistfloat, minweight, subs)


text:
@testset "reported issues" begin
    @testset "barplot nl in text (#102)" begin
        p = barplot(["I need a\nbreak", "I don't but could", :Hello, :Again], [30, 40, 20, 10])
        test_ref("references/issues/barplot_nl.txt", @show_col(p))
    end

    @testset "plot hang (tan → ∞) (#129)" begin
        p = lineplot([cos, sin, tan], -π/2, 2π)
        test_ref("references/issues/cos_sin_tan.txt", @show_col(p))
    end

    @testset "integer edges (#139)" begin
        p_int = histogram(fit(Histogram, rand(10), 0:3))
        p_float = histogram(fit(Histogram, rand(10), 0.0:1.0:3.0))
        @test @print_col(p_int) == @print_col(p_float)
    end

    @testset "isolated colorbar (#169)" begin
        p = @inferred Plot([0], [0], colorbar=true, colormap=:cividis, width=0, min_width=0)
        test_ref("references/issues/isolated_colorbar.txt", @show_col(p))
    end
end


text:
module StanModels

using Reexport

@reexport using StanSample

end # module


text:
module ARFF

# Simply uses the CSV.jl parser for fast ARFF loading with the header parsing logic from ARFFFiles.jl
# Once there is a sufficiently fast ARFF library available in Julia, this module will become obsolete

# Doesn't support:
# - multi-instance attributes
# - weighted data
# - sparse data

using CSV
using Dates:DateTime

export ARFFType,
       ARFFNumericType,
       ARFFStringType,
       ARFFDateType,
       ARFFNominalType,
       ARFFRelation,
       ARFFAttribute,
       ARFFData,
       ARFFHeader,
       header,
       read

include("parsing.jl")

# The unique name of the relation in the corresponding file.
struct ARFFRelation
    name::String
end

# Abstract type of ARFF types
abstract type ARFFType end

struct ARFFNumericType <: ARFFType
    native::Type
end
ARFFNumericType() = ARFFNumericType(Float64)

struct ARFFStringType <: ARFFType
        native::Type
end
ARFFStringType() = ARFFStringType(String)

struct ARFFDateType <: ARFFType
    native::Type
    format::String
end
ARFFDateType() = ARFFDateType(DateTime, "yyyy-MM-dd'T'HH:mm:ss")

struct ARFFNominalType <: ARFFType
    native::Type
    classes::Vector{String}
end
ARFFNominalType(classes::Vector{String}) = ARFFNominalType(String, classes)

# Represents a single ARFF @attribute.
struct ARFFAttribute
    name::String
    type::ARFFType
end

# Indicates the start of the data rows (to be parsed as CSV).
struct ARFFData end

# Represents the header information in an ARFF file.
struct ARFFHeader
    relation::String
    attributes::Vector{ARFFAttribute}
end

function header(path)
    @debug "Reading header metadata"
    open(path, "r") do io
        return parse_header(io)
    end
end

function read(path; kwargs...)
    h, lineno = header(path)

    @debug "Reading data"
    colnames = map(a -> a.name, h.attributes)
    coltypes = map(a -> a.type.native, h.attributes)
    formats = [a.format for a in h.attributes if isa(a, ARFFDateType)]
    @assert all(y -> y==formats[1], formats) "Expected all date formats to be equal (see constraints in the README)"

    f = CSV.File(path;
        skipto=lineno+1,
        header=colnames,
        types=coltypes,
        missingstring="?",
        comment="%",
        dateformat=length(formats) > 0 ? formats[1] : nothing,
        kwargs...)
    return f
end

end


text:
using Distributions, LinearAlgebra, Random, Measures, Plots; pyplot()
Random.seed!(1)

a, varXi, varZeta, = 0.8, 0.36, 1.0
X0, spikeTime, spikeSize = 10.0, 50, -20.0
Tsmall, Tlarge, warmTime = 100, 10^6, 10^4
kKalman = 0.3

function luenbergerTrack(k, T, spikeTime = Inf)
    X, Xhat = X0, 0.0
    xTraj, xHatTraj, yTraj = [X], [Xhat,Xhat], [X0]

    for t in 1:T-1
        X = a*X + rand(Normal(0,sqrt(varXi)))
        Y = X + rand(Normal(0,sqrt(varZeta)))
        Xhat =a*Xhat - k*(Xhat - Y)

        push!(xTraj,X)
        push!(xHatTraj,Xhat)
        push!(yTraj,Y)

        if t == spikeTime
            X += spikeSize
        end
    end
    deleteat!(xHatTraj,length(xHatTraj))
    xTraj, xHatTraj, yTraj
end

smallTraj, smallHat, smallY = luenbergerTrack(kKalman, Tsmall, spikeTime)

p1 = scatter(smallTraj, c=:blue,
        ms=3, msw=0, label="System trajectory")
p1 = scatter!(smallY, c = :black,
        ms=3, msw=0, label="Measurments")
p1 = scatter!(smallHat,c=:red,
        ms=3, msw=0, label="Kalman filter tracking",
        xlabel = "Time", ylabel = "Temperature",
        xlims=(0, Tsmall))

kRange = 0.2:0.005:0.4
errs = []
for k in kRange
    xTraj, xHatTraj, _ = luenbergerTrack(k, Tlarge)
    mse = norm(xTraj[warmTime:end] - xHatTraj[warmTime:end])^2/(Tlarge-warmTime)
    push!(errs, mse)
end

analyticErr(k) =(varXi + k^2*varZeta) / (1-(a-k)^2)

p2 = scatter(kRange,errs, c=:black, ms=3, msw=0,
    xlabel="k", ylabel="MSE", label = "Monte Carlo")
p2 = plot!(kRange,analyticErr.(kRange), c = :red,
    xlabel="k", ylabel="MSE", label = "Analytic", ylim =(0.58,0.64))

plot(p1, p2, size=(1000,400), margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
struct InputField
    label::String
    μ::Widget
    units::Union{Unitful.FreeUnits, Nothing}
    field::Node
    dumby::Bool
end

"""
    InputField( label::String, μ::Float64, units::Union{ Unitful.FreeUnits,Nothing })

Creates a labelled textbox input field in Blink.

"""
function InputField( label::String, μ::Float64, units::Union{ Unitful.FreeUnits,Nothing })
    mu = Widgets.textbox(; value = string( μ ) )
    if !isa(units, Nothing)
        units_str = string(units)
        field = Interact.hbox( label * " ($(units_str)): ", mu )
    else
        field = Interact.hbox( label * ":", mu )
    end
    return InputField( label, mu, units, field, true )
end

"""
    value(uif::InputField)

Retrieves a numeric value from a labelled textbox InputField.

"""
function value(uif::InputField)
    if isa(uif.units, Nothing)
        return parse(Float64, uif.μ[] )
    else
        type = string( typeof( uif.units ).parameters[1][1] )
        return uparse(uif.μ[] * type )
    end
end

struct UncertainInputField
    label::String
    μ::Widget
    σ::Widget
    units::Unitful.FreeUnits
    field::Node
end

"""
    UncertainInputField( label::String, μ::Float64, σ::Float64, units::Unitful.FreeUnits)

Creates a labelled textbox input field with an uncertainty textbox in Blink.

"""
function UncertainInputField( label::String, μ::Float64, σ::Float64, units::Unitful.FreeUnits)
    ( σ > μ ) && @warn("Uncertainty greater then value. Are you sure you understand the parameters?")
    mu = Widgets.textbox(; value = string( μ ) )
    sigma = Widgets.textbox(; value = string( σ ) )
    units_str = string(units)
    field = Interact.hbox( label * " ($(units_str)): ", mu, " ± ", sigma  )
    return UncertainInputField( label, mu, sigma, units, field )
end

(uif::UncertainInputField)() = uif.field

"""
    value(uif::UncertainInputField)

Retrieves a numeric value with uncertainty from a labelled textbox UncertainInputField.

"""
function value(uif::UncertainInputField)
    type = string( typeof( uif.units ).parameters[1][1] )
    return uparse(uif.μ[] * type ) ± uparse(uif.σ[] * type )
end

mutable struct CalibrationData
    steps::Vector{Float64}
    displacement::Vector{Float64}
    disp_uncertainty::Vector{Float64}
    beta::Float64
    beta_uncertainty::Float64
end

function (cd::CalibrationData)(step, disp, disp_u)
    push!(cd.steps, step)
    push!(cd.displacement, disp)
    push!(cd.disp_uncertainty, disp_u)
    if length(cd.steps) > 1
        disp        = cd.displacement .± cd.disp_uncertainty
        steps       = cumsum( cd.steps )
        b           = inv(disp'*disp)*disp'*steps
        cd.beta = b.val
        cd.beta_uncertainty = b.err
    end
    return cd
end

struct Calibration
    steps_and_disps::Observable
    widget::Node
end

function Calibration( units::String )
    units_str = units

    steps = Widgets.textbox()
    addbtn = Widgets.button("Add")
    stepsentry = hbox( "Steps(#): ", steps )

    disp = Widgets.textbox()
    unc_disp = Widgets.textbox()
    dispentry = hbox( "Displacement($units_str): ", disp, " ± ", unc_disp )

    calib_data = Observable( CalibrationData(Vector{Float64}(undef,0), Vector{Float64}(undef,0),
                                             Vector{Float64}(undef,0), 0.0, 0.0) )

    map!( x -> calib_data[]( parse(Float64, steps[]), parse(Float64, disp[]),
                             parse(Float64, unc_disp[]) ),
                             calib_data, addbtn )

    return Calibration( calib_data, Interact.vbox( stepsentry, dispentry, addbtn ) )
end

cumdiff( a ) = [ 0; -cumsum( reverse( diff( reverse( a ) ) ) ) ]

function makeplot(calib_data)
    disp        = calib_data.displacement .± calib_data.disp_uncertainty
    steps       = cumsum( calib_data.steps )
    retplot     = scatter(  steps, map(x -> x.val, disp), legend = :topleft,
                            yerr = map(x -> x.err, disp), label  = "Calibration",
                            xlabel = "Steps (#)", ylabel = "Displacement (mm)")
    retplot     = Plots.abline!(retplot, 1.0/calib_data.beta, 0.0,
                                line=:dash, color = :red, label = "Trendline")
    return retplot
end

"""
    launch_pump_interface(com_port::String; baud_rate::Int64 = 9600)

Launches a Blink/Electron WebUI for calibrating and driving the pump.

"""
function launch_pump_interface(com_port::String; baud_rate::Int64 = 9600)
    mcu = nothing
    try
        mcu = open( "/dev/ttyACM0", baud_rate )
    catch
        @error("Failed to establish COM connection at: $com_port.")
    end

    driver          = InputField( "Move N Steps", 100.0,  nothing);
    driver_btn      = Widgets.button("Drive Pump");
    c               = Calibration("mm");
    calibration_plt = Observable{Any}(nothing);
    drive_status    = Observable{Any}(nothing);
    map!( x ->  drive_pump( mcu, convert( Int64, value(driver) ),
                convert( Int64, 10000) ), drive_status, driver_btn)

    map!( x -> makeplot( x ), calibration_plt, c.steps_and_disps)

    calibrate_txt = "Drive the stepper motor N steps, and measure the displacement of the plunger with a caliper." *
                    "Add the experimental values to the calibration model." *
                    "After numerous runs, you will be given a calibration + uncertainty for dispensing fluids based on a linear model."

    calibrationbox = vbox( "Calibrate:", calibrate_txt, hbox(driver.field, driver_btn),
                            calibration_plt, c.widget
                    );

    units       = UncertainInputField( "Marking vs Volume", 50.0, 0.02, u"mm / mL");
    volume      = InputField( "Volume", 5.0, u"mL");
    duration    = InputField( "Duration", 3.0, u"s");
    driver_btn2 = Widgets.button("Drive Pump");
    stop_btn    = Widgets.button("Stop Pump");

    map!( x -> drive_physical( mcu, c.steps_and_disps[], units, volume, duration ), drive_status, driver_btn2)

    driverbox   = vbox( units.field, volume.field, duration.field,
                        hbox(driver_btn2, stop_btn) );
    tabbed      = Widgets.tabulator(Observable( Dict( "Calibrate" => calibrationbox, "Driver" => driverbox ) ));

    w = Window(); body!(w, tabbed)
    return w
end

max_stars_repo_name: JuliaMakers/MakerProjects

max_issues_repo_name: JuliaMakers/MakerProjects

max_forks_repo_name: JuliaMakers/JuliaMakerProjects


text:
using SMC
import JSON

function load_states(fname::String)
    states = JSON.parsefile(fname)
    return convert(Vector{Int}, states)
end

function load_observations(fname::String)
    observations = JSON.parsefile(fname)
    return convert(Vector{Int}, observations)
end

"""
Load matrix from column-major format
"""
function to_float_matrix(vec::Vector, num_rows::Int, num_cols::Int)
    if length(vec) != num_cols
        error("Dimension mismatch: vector has length $(length(vec)) but expected $num_cols")
    end
    mat = fill(NaN, (num_rows, num_cols))
    for (i, column) in enumerate(vec)
        if length(column) != num_rows
            error("Dimension mismatch")
        end
        mat[:,i] = column
    end
    if any(isnan.(mat)[:])
        error("Dimension mismatch")
    end
    return mat
end

function load_hmm(fname::String)
    dict = JSON.parsefile(fname)
    num_states = dict["num_states"]
    num_obs = dict["num_obs"]
    initial_state_prior = dict["initial_state_prior"]
    initial_state_prior = convert(Vector{Float64}, initial_state_prior)
    if length(initial_state_prior) != num_states
        error("Dimension mismatch")
    end
    transition_model = dict["transition_model"]
    transition_model = to_float_matrix(transition_model, num_states, num_states)
    observation_model = dict["observation_model"]
    observation_model = to_float_matrix(observation_model, num_states, num_obs)
    log_initial_state_prior = dict["log_initial_state_prior"]
    log_transition_model = dict["log_transition_model"]
    log_observation_model = dict["log_observation_model"]
    return HiddenMarkovModel(initial_state_prior, transition_model, observation_model)
end


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Base.Test, DiffEqBase
prob = prob_ode_linear

sol = solve(prob,BS3();dt=1//2^(4),tstops=[0.5],saveat=0:0.01:1)
sol(0.9)

integrator = init(prob,BS3();dt=1//2^(4),tstops=[0.5],saveat=0:0.01:1)
step!(integrator)
@test integrator.iter == 1
solve!(integrator)
@test integrator.t == 1.0
integrator(0.95)
integrator.tprev
integrator.t

push!(integrator.opts.tstops,5.0)
integrator.opts.advance_to_tstop=true
step!(integrator)
@test integrator.t == 5.0
integrator.opts.advance_to_tstop=false
step!(integrator)
@test integrator.t > 5

@test integrator.sol(0.9) == sol(0.9)

integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5],advance_to_tstop=true)
tupint = tuples(integrator)
for (t,u) in tuples(integrator)
  @test t ∈ [0.5,1.0]
end

integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5],advance_to_tstop=true,stop_at_next_tstop=true)
for (t,u) in tuples(integrator)
  @test t == 0.5
end
integrator([1.0;2.0])

integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5])
for (tprev,uprev,t,u) in intervals(integrator)
  @show tprev,t
end
integrator([1.0;2.0])

integrator = init(prob,RK4();dt=1//2^(9))
for i in take(integrator,12) end
@test integrator.iter == 12
for i in take(integrator,12) end
@test integrator.iter == 24

integrator = init(prob_ode_2Dlinear,Tsit5();dt=1//2^(4),tstops=[0.5],advance_to_tstop=true,stop_at_next_tstop=true)
for (t,u) in tuples(integrator)
  @test t == 0.5
end
A = integrator([1.0;2.0])
B = integrator([1.0;2.0],idxs=1:2:5)

@test minimum([A[i][1:2:5] == B[i] for i in 1:length(A)])

integrator(A[1],0.5)
@test A[1] == integrator(0.5)

A = zeros(3)
integrator(A,0.6,idxs=1:2:5)
@test A == integrator(0.6;idxs=1:2:5)