text:
# remove deps.jl if it exists, in case build.jl fails
isfile("deps.jl") && rm("deps.jl")

using Libdl

println("Pardiso library")
println("===============")

const LIBPARDISONAMES =
if Sys.iswindows()
[
    "libpardiso600-WIN-X86-64.dll",
]
elseif Sys.isapple()
[
    "libpardiso600-MACOS-X86-64.dylib",
]
elseif Sys.islinux()
[
    "libpardiso600-GNU800-X86-64.so",
    "libpardiso600-GNU720-X86-64.so",
]
else
    error("unhandled OS")
end

println("Looking for libraries with name: ", join(LIBPARDISONAMES, ", "), ".")

PATH_PREFIXES = [@__DIR__; get(ENV, "JULIA_PARDISO", [])]

if !haskey(ENV, "JULIA_PARDISO")
    println("INFO: use the `JULIA_PARDISO` environment variable to set a path to " *
            "the folder where the Pardiso library is located")
end

function find_paradisolib()
    found_lib = false
    for prefix in PATH_PREFIXES
        println("Looking in \"$(abspath(prefix))\" for libraries")
        for libname in LIBPARDISONAMES
            local path
            try
                path = joinpath(prefix, libname)
                if isfile(path)
                    println("    found \"$(abspath(path))\", attempting to load it...")
                    Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
                    println("    loaded successfully!")
                    global PARDISO_LIB_FOUND = true
                    return path, true
                end
            catch e
                println("    failed to load due to:")
                Base.showerror(stderr, e)
            end
        end
    end
    println("did not find libpardiso, assuming PARDISO 5/6 is not installed")
    return "", false
end

pardisopath, found_pardisolib = find_paradisolib()

#################################################

println("\nMKL Pardiso")
println("=============")
function find_mklparadiso()
    if haskey(ENV, "MKLROOT")
        println("found MKLROOT environment varaible, enabling local MKL")
        return true
    end
    println("did not find MKLROOT environment variable, using provided MKL")
    return false
end

found_mklpardiso = find_mklparadiso()

open("deps.jl", "w") do f
    print(f,
"""
const LOCAL_MKL_FOUND = $found_mklpardiso
const PARDISO_LIB_FOUND = $found_pardisolib
const PARDISO_PATH = raw"$pardisopath"
"""
)

end


text:
using GRIB
using Test
using Statistics

@testset "GRIB" begin

    GribFile(joinpath(dirname(@__FILE__), "samples", "regular_latlon_surface.grib2")) do f

        # Test that it starts at 0
        @test position(f) == 0

        # Get the message
        msg = Message(f)

        # Test that position has changed
        @test position(f) == 1

        # Test haskey
        @test haskey(msg, "shortName")
        @test !haskey(msg, "gibberish")

        # Test getindex
        @test msg["shortName"] == "2t"
        @test msg["level"] == 2
        @test missingvalue(msg) == 9999
        @test typeof(maskedvalues(msg)) == Array{Union{Missing, Float64}, 2}

        # Test data
        lons, lats, vals = data(msg)
        validlats = collect(0.:2:60.)
        validlons = collect(0.:2:30.)
        @test all([l in validlats for l in lats])
        @test all([l in validlons for l in lons])
        @test size(vals) == (16, 31)
        @test mean(vals) == msg["average"]
        @test vals == msg["values"]
        nx, ny = size(lons)
        @test all([lons[1, y] - lons[1, 1] == 0 for y in 1:ny]) || all([lons[x, 1] - lons[1, 1] == 0 for x in 1:nx])
        @test all([lats[1, y] - lats[1, 1] == 0 for y in 1:ny]) || all([lats[x, 1] - lats[1, 1] == 0 for x in 1:nx])

        # Test clone, setindex!, and write
        msg2 = clone(msg)
        @test data(msg2) == data(msg)
        msg2["shortName"] = "t"
        msg2["level"] = 10
        msg2["latitudeOfFirstGridPointInDegrees"] = msg["latitudeOfFirstGridPointInDegrees"] + 10
        msg2["values"] = msg["values"] .+ 20
        writemessage(msg2, joinpath(dirname(@__FILE__), "samples", "modified_latlon_surface.grib2"))
        GribFile(joinpath(dirname(@__FILE__), "samples", "modified_latlon_surface.grib2")) do f
            msg3 = Message(f)
            @test msg3["shortName"] == "t"
            @test msg3["level"] == 10
        end

        # Just testing to make sure these don't throw
        bytes = getbytes(msg)
        msg3 = Message(bytes)
        @test msg3["shortName"] == msg["shortName"]
        bytesarr = [bytes, getbytes(msg2)]
        msg4 = Message(bytesarr)
        msg5 = Message(msg3)
        msg6 = deepcopy(msg5)

        # Test GribIterate
        lons1d = Vector{Float64}(undef, length(lons))
        lats1d = similar(lons1d)
        vals1d = similar(lons1d)
        @test length(eachpoint(msg)) == length(vals)
        for (i, pt) in enumerate(eachpoint(msg))
            lons1d[i] = pt[1]
            lats1d[i] = pt[2]
            vals1d[i] = pt[3]
        end
        @test vec(lons) == lons1d
        @test vec(lats) == lats1d
        @test vec(vals) == vals1d

        # Test keys and values
        ks = Set{String}()
        for k in keys(msg)
            push!(ks, k)
        end
        @test eltype(keys(msg)) == String
        @test "skewness" in ks && "latitudes" in ks && "centre" in ks && "maximum" in ks
        vs = Set()
        for v in values(msg)
            push!(vs, v)
        end
        @test "2t" in vs && 2 in vs && 9999 in vs
        ks2 = Set{String}()
        vs2 = Set()
        for (key, value) in msg
            push!(ks2, key)
            push!(vs2, value)
        end
        @test issetequal(ks, ks2)
        @test issetequal(vs, vs2)

        # Test Nearest
        lons, lats, vals, dists = Nearest(msg) do n
            find(n, msg, 25.75, 25.25)
        end
        @test lons == [26, 26, 24, 24]
        @test lats == [26, 24, 26, 24]
        @test all(isapprox.(dists, [87.08, 141.28, 194.27, 224.97], atol=0.01))

        inlons = collect(2.:2:10.)
        inlats = collect(2.:2:10.)
        lons, lats, vals, dists = findmultiple(msg, inlons, inlats)
        @test all([d == 0 for d in dists])
        @test_throws DomainError findmultiple(msg, inlons[1:2], inlats)
    end

    # Test Index
    Index(joinpath(dirname(@__FILE__), "samples", "flux.grb"), "edition") do index
        select!(index, "edition", 2)
        count = 0
        for msg in index
            count += 1
        end
        @test count == 4
    end

    Index(joinpath(dirname(@__FILE__), "samples", "flux.grb"), "typeOfLevel") do index
        select!(index, "typeOfLevel", "surface")
        @test keycount(index, "typeOfLevel") == 2
        names = Set{String}()
        for msg in index
            push!(names, msg["shortName"])
        end
        @test issetequal(names, ["prate", "sp"])
    end

    Index(joinpath(dirname(@__FILE__), "samples", "flux.grb"), "latitudeOfFirstGridPointInDegrees") do index
        addfile!(index, joinpath(dirname(@__FILE__), "samples", "regular_latlon_surface.grib2"))
        select!(index, "latitudeOfFirstGridPointInDegrees", 88.542)
    end

    GribFile(joinpath(dirname(@__FILE__), "samples", "flux.grb")) do f
        count = 0
        for msg in f
            count += 1
        end
        @test count == 4
    end

    # Test the seeking methods
    GribFile(joinpath(dirname(@__FILE__), "samples", "flux.grb")) do f
        @test eltype(f) == Message
        msgs = read(f, 2)
        @test length(msgs) == 2
        @test position(f) == 2
        seekstart(f)
        @test position(f) == 0
        seek(f, 3)
        @test position(f) == 3
        skip(f, -2)
        @test position(f) == 1
        skip(f, 1)
        @test position(f) == 2
        seek(f, 1)
        @test position(f) == 1
        @test_throws DomainError seek(f, -1)
        @test_throws DomainError skip(f, -15)
        @test length(read(f, 8)) == 3
        @test position(f) == 9
    end

    # Test that it doesn't try to open a non-existent file
    if isfile(joinpath(dirname(@__FILE__), "samples", "nofilehere.grb"))
        rm(joinpath(dirname(@__FILE__), "samples", "nofilehere.grb"))
    end
    @test_throws SystemError GribFile(joinpath(dirname(@__FILE__), "samples", "nofilehere.grb"))

    # Test multi-field
    GribFile(joinpath(dirname(@__FILE__), "samples", "isobaricsmaller.grb2")) do f
        shortnames = AbstractString[]
        for msg in f
            push!(shortnames, msg["shortName"])
        end
        @test "v" in shortnames
        @test length(shortnames) == 18
    end

    GribFile(joinpath(dirname(@__FILE__), "samples", "isobaricsmaller.grb2")) do f
        @test eltype(f) == Message
        msgs = read(f, 2)
        @test length(msgs) == 2
        @test position(f) == 2
        seekstart(f)
        @test position(f) == 0
        seek(f, 3)
        @test position(f) == 3
        skip(f, -2)
        @test position(f) == 1
        skip(f, 1)
        @test position(f) == 2
        seek(f, 1)
        @test position(f) == 1
    end

    nomultisupport()
    GribFile(joinpath(dirname(@__FILE__), "samples", "isobaricsmaller.grb2")) do f
        shortnames = AbstractString[]
        for msg in f
            push!(shortnames, msg["shortName"])
        end
        @test !("v" in shortnames)
        @test length(shortnames) == 12
    end

    GribFile(joinpath(dirname(@__FILE__), "samples", "isobaricsmaller.grb2")) do f
        @test eltype(f) == Message
        msgs = read(f, 2)
        @test length(msgs) == 2
        @test position(f) == 2
        seekstart(f)
        @test position(f) == 0
        seek(f, 3)
        @test position(f) == 3
        skip(f, -2)
        @test position(f) == 1
        skip(f, 1)
        @test position(f) == 2
        seek(f, 1)
        @test position(f) == 1
    end

    # Test that iterator methods are defined because Julia assumes weird defaults
    GribFile(joinpath(dirname(@__FILE__), "samples", "isobaricsmaller.grb2")) do f
       @test Base.IteratorSize(f) == Base.SizeUnknown()
       msg = read(f)
       @test Base.IteratorSize(keys(msg)) == Base.SizeUnknown()
       @test Base.IteratorSize(values(msg)) == Base.SizeUnknown()
    end

    # Test that the local definitions are working
    GribFile(joinpath(dirname(@__FILE__), "samples", "hrrr.echotop.grib2")) do f
        msg = read(f)
        @test msg["name"] == "Echo Top"
    end
end

ext: jl


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

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

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

const _property_map_VirtualHardDisk = Dict{Symbol,Symbol}(Symbol("uri")=>Symbol("uri"))
const _property_types_VirtualHardDisk = Dict{Symbol,String}(Symbol("uri")=>"String")
Base.propertynames(::Type{ VirtualHardDisk }) = collect(keys(_property_map_VirtualHardDisk))
Swagger.property_type(::Type{ VirtualHardDisk }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualHardDisk[name]))}
Swagger.field_name(::Type{ VirtualHardDisk }, property_name::Symbol) =  _property_map_VirtualHardDisk[property_name]

function check_required(o::VirtualHardDisk)
    true
end

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


text:
using GeometryBasics;
using FileIO;
using LinearAlgebra;
using GaussianProcesses;
using MovingFrame;
import Makie;
using OwnTime;
using Profile;
const mv = MovingFrame
##11
meshfile = Base.download("https://github.com/libigl/libigl-tutorial-data/raw/master/bumpy-cube.obj")
shape = FileIO.load(FileIO.File{FileIO.format"OBJ"}(meshfile))
verts = decompose(Point3f0, shape)
faces = decompose(TriangleFace{Int}, shape)
vmat = MovingFrame.array_of_array_to_matrix(verts)
fmat = MovingFrame.array_of_array_to_matrix(faces)
##
mv.surface_area2(vmat,fmat)
mv.surface_area(vmat,fmat)
##
mv.mean_curvature_flow(verts,faces)
mv.mean_curvature_flow(vmat,fmat)


text:
module Turf
    import GeoJSON

    export
        inside,
        within,
        tag,
        remove,
        filter,
        distance,
        area,
        nearest,
        bboxpolygon,
        envelope,
        extent,
        square,
        # size,
        center,
        centroid,
        midpoint,
        bearing,
        destination,
        linedistance,
        along,
        quantile,
        reclass,
        convex,
        flip,
        explode,
        combine,
        isclockwise,
        kinks

    include("joins.jl")
    include("data.jl")
    include("measurement.jl")
    # include("interpolation.jl")
    include("classification.jl")
    # include("aggregation.jl")
    # include("transformation.jl")
    include("misc.jl")
    # include("wgs84constants.jl")
end


text:
# a bunch of overlays for debugging
mutable struct BeliefOverlay <: SceneOverlay
    verbosity::Int
    color::Colorant
    font_size::Int
    belief::Any
    policy::DecomposedPolicy

    function BeliefOverlay(belief, policy::DecomposedPolicy,
        verbosity::Int=1;
        color::Colorant=colorant"white",
        font_size::Int=20
        )

        new(verbosity, color,font_size, belief, policy)
    end
end

function AutoViz.render!(rendermodel::RenderModel, overlay::BeliefOverlay, scene::Scene, env::OccludedEnv)
    AutoViz.render!(rendermodel, overlay.belief, overlay, scene, env)
    return rendermodel
end

function AutoViz.render!(rendermodel::RenderModel, belief::PreviousObsDecBelief, overlay::BeliefOverlay, Scene::Scene, env::OccludedEnv)
    for (prob_key, b) in belief
        n_obstacles = 0
        if prob_key == :obsped || prob_key == :obscar
            n_obstacles = 3
        end
        problem = overlay.policy.problem_map[prob_key]
        AutoViz.render!(rendermodel, overlay,n_obstacles, problem, b, env)
    end
end

function AutoViz.render!(rendermodel::RenderModel, belief::KMarkovDecBelief, overlay::BeliefOverlay, Scene::Scene, env::OccludedEnv)
    for (prob_key, b) in belief
        n_obstacles = 0
        if prob_key == :obsped || prob_key == :obscar
            n_obstacles = 3
        end
        problem = overlay.policy.problem_map[prob_key]
        for i=1:size(b, 2)
            AutoViz.render!(rendermodel, overlay,n_obstacles, problem, b[:,i][:], env)
        end
    end
end

function AutoViz.render!(rendermodel::RenderModel, overlay::BeliefOverlay, n_obstacles::Int64, problem::UrbanPOMDP, b::Vector{Float64}, env::OccludedEnv, transparency = 0.1)
    o = deepcopy(b)
    o = unrescale!(o, problem)
    ego, car_map, ped_map, obs_map = AutomotivePOMDPs.split_o(o, problem, 4, n_obstacles)
    for (key, state) in car_map
        x, y, θ, v = state
        color = color=RGBA(75.0/255, 66.0/255, 244.0/255, transparency)
        add_instruction!(rendermodel, render_vehicle, (x, y, θ, problem.car_type.length, problem.car_type.width, color, color, RGBA(1.,1.,1.,transparency)))
    end
    for (key, state) in ped_map
        x, y, θ, v = state
        color = RGBA(75.0/255, 66.0/255, 244.0/255, transparency)
        add_instruction!(rendermodel, render_vehicle, (x, y, θ, problem.ped_type.length, problem.ped_type.width, color, color, RGBA(1.,1.,1.,transparency)))
    end
end

function AutomotivePOMDPs.animate_scenes(scenes::Vector{Scene},
                        actions::Vector{Float64},
                        beliefs,
                        policy::DecomposedPolicy,
                        env;
                        overlays::Vector{SceneOverlay} = SceneOverlay[],
                        cam=FitToContentCamera(0.),
                        sim_dt=0.1)
    duration =length(scenes)*sim_dt
    fps = Int(1/sim_dt)
    function render_rec(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        return AutoViz.render(scenes[frame_index],
                              env,
                              cat(1, overlays,
                                     TextOverlay(text = ["Acc: $(actions[frame_index])"]),
                                     BeliefOverlay(beliefs[frame_index], policy)),
                              cam=cam,
                              car_colors=get_colors(scenes[frame_index]))
    end
    return duration, fps, render_rec
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
"""
	mmconv(sig,imp)
Convolve two signals and cut the result to valid time interval only (cut out edges)
Cutting of edges is the only difference between Julias `conv` and this function

**Input:**
* sig: signal to be filtered (vector)
* imp: filter = impulse response (vector)
> length of `imp` must be an odd number

**Output**
* outsig: filtered 'valid' input signal, i.e., NaNs where filter/edge effect
  (same length as `sig`)

**Example**
```
t = collect(1:1:400.);
s = sin.(2*pi*1/50.*t) .+ randn(length(t))./6;
sf = mmconv(s,ones(7)./7);
```
"""
function mmconv(sig::Vector{Float64},imp::Vector{Float64})
	outsig = DSP.conv(sig,imp);
	return correctconv(outsig,length(imp));
end

"""
	correctconv(sig,impl)
Correct phase shift and edge effect (replace by NaNs)
**Input**
* sig: signal to be filtered (vector)
* impl: length of impulse response (=filter length)
> `impl` must be odd number

**Output**
* valid output signal

**Example**
```
sig = [1.,2.,3.,4.,3.,2.,1.,0.,-1];
imp = [1/3,1/3,1/3];
outsig = correctconv(sig,length(imp));
```
"""
function correctconv(sig::Vector{Float64},impl::Int)
	iseven(impl) ? error("Filter must have odd number of coefficients") : nothing
	temp = div(impl,2);
	# Remove phase shift
	out = sig[1+temp:end-temp];
	# Remove filter/edge effect
	out[1:temp] .= NaN; # left
	out[end-temp+1:end] .= NaN; # right
	return out;
end

"""
	findblocks(invec)
Find blocks without NaNs

**Input**
* invec: input vector to be examined. Regular sampling is assumed! (use time2regular)

**Output**
* start,stop indices of blocks without NaNs

**Example**
```
invec = collect(1.:1:17.);
invec[[6,10,11,12,14,17]] .= NaN;
istart,istop = findblocks(invec);
```
"""
function findblocks(invec::Vector{Float64})
	idall = findall(isnan.(invec));
	findblocks_decide(idall,length(invec));
end

function findblocks_decide(idall::Vector{Int},length_invec::Int)
	if !isempty(idall)
		return findblocks_main(idall,length_invec);
	else
		return [1],[length_invec]
	end
end
function findblocks_main(idall::Vector{Int},length_invec::Int)
	idstart = Vector{Int}();
	idstop = Vector{Int}();
	for i in 1:length(idall)
		if i == 1 && i != idall[1]
			push!(idstart,1);
			push!(idstop,idall[1]-1);
		elseif i == 1
			continue;
		elseif idall[i]-idall[i-1] != 1
			push!(idstart,idall[i-1]+1)
			push!(idstop,idall[i]-1)
		end
	end
	if idall[end] != length_invec
		push!(idstart,idall[end]+1);
		push!(idstop,length_invec);
	end
	return idstart,idstop
end

"""
	findnanblocks(invec)
Find blocks of NaNs

**Input**
* invec: input vector to be examined. Regular sampling is assumed! (use time2regular)

**Output**
* start,stop indices of blocks of NaNs

**Example**
```
invec = collect(1.:1:17.);
invec[[6,10,11,12,14,17]] .= NaN;
nstart,nstop = findnanblocks(invec);
```
"""
function findnanblocks(invec)
	inans = map(isnan,invec);
	ostart,ostop = Vector{Int64}(), Vector{Int64}();
	if any(inans)
		istart,istop = findblocks(invec);
		ostart = istop[1:end-1] .+ 1;
		ostop = istart[2:end] .- 1;
		if istart[1] != 1
			ostart = vcat(1,ostart);
			ostop = vcat(istart[1]-1,ostop);
		end
		if istop[end] != length(invec)
			ostart = vcat(ostart,istop[end]+1);
			ostop = vcat(ostop,length(invec));
		end
	end
	return ostart,ostop
end

"""
	filtblocks(sig,imp)
Filter signal assuming input time series contains NaNs. Thus, piecewise filtering
of the input singal will be applied (`mmconv` will be utilized).

See `mmconv` help for **Input** and **Output**

**Example**
```
t = collect(1:1:100.);
s = sin.(2*pi*1/20.*t) .+ randn(length(t))./6;
s[40:49] .= NaN;
sf = filtblocks(s,ones(5)./5);
```
"""
function filtblocks(sig::Vector{Float64},imp::Vector{Float64})
	id1,id2 = findblocks(sig);
	out = Vector{Float64}(undef,length(sig)) .+ NaN;
	for i in 1:1:length(id1)
		if id2[i]-id1[i] > length(imp)*2.
			out[id1[i]:id2[i]] = mmconv(sig[id1[i]:id2[i]],imp);
		end
	end
	return out
end

"""
	demean(sig)
subtract mean value from input Vector

**Input:**
* sig: data to be reduced (can contain NaNs)

**Output**
* outsig: data after subtraction of mean value

**Example**
```
s = [-1.,NaN,0.,1.] .+ 1.234
sf = demean(s);
```
"""
function demean(sig::Vector{Float64})::Vector{Float64}
	c = 0;s = 0.;
	for i in sig
		if !isnan(i)
			c += 1;
			s += i;
		end
	end
	return sig .- s/c
end

"""
	detrend(x,y;deg)
subtract polynomial from input Vector

**Input:**
* x: x coordinates
* y: y coordinates
* deg: polynomial degree to be subtracted (0=mean,1=linear,etc)

**Output**
* outsig: data after subtraction of given polynomial

**Example**
```
x = collect(1.:1:10);
y = ones(length(x));
out = detrend(x,y,deg=1)
```
"""
function detrend(x,y;deg::Int=1)
	if deg==0
		return demean(y)
	else
		xc,yc = copy(x),copy(y);
		# remove NaNs
		xc = xc[.!isnan.(y)];
		yc = yc[.!isnan.(y)];
		# fit
		fit,err = ResampleAndFit.fitpoly(xc,yc,deg=deg);
		return y - ResampleAndFit.evalpoly(x,fit);
	end
end

"""
	defirst(sig)
subtract first not NaN value from input Vector

**Input:**
* sig: data to be reduced (can contain NaNs)

**Output**
* outsig: data after subtraction of fist not NaN value

**Example**
```
s = [-1.,NaN,0.,1.];
sf = defirst(s);
```
"""
function defirst(sig::Vector{Float64})::Vector{Float64}
	s = 0.;
	for i in sig
		if !isnan(i)
			s = i;
			break;
		end
	end
	return sig .- s
end


text:
# Graph analysis tools

function rangeErrMaxPoint2(fgl1::FactorGraph, id1, fgl2::FactorGraph ,id2)
  mv1 = getKDEMax(getVertKDE(fgl1,id1))
  mv2 = getKDEMax(getVertKDE(fgl2,id2))
  return norm(mv1[1:2]-mv2[1:2])
end

function rangeCompAllPoses(fgl1::FactorGraph, fgl2::FactorGraph)
  ranges = Float64[]
  xx,ll = ls(fgl1)
  for x in xx
    push!(ranges, rangeErrMaxPoint2(fgl1,x,fgl2,x))
  end
  return ranges
end

function rangeCompAllPoses(
    valsbaseline::Dict{Int,Array{Float64,1}},
    fglbaseline::FactorGraph,
    fgltest::FactorGraph)

  ranges = Float64[]
  xx,ll = ls(fgltest)
  for x in xx
    mv1 = valsbaseline[fglbaseline.IDs[x]]
    mv2 = getKDEMax(getVertKDE(fgltest,x))
    push!(ranges, norm(mv1[1:2]-mv2[1:2]))
  end
  return ranges
end


text:
function count_nucleotides(strand::AbstractString)

end


text:
using Gadfly

function main()
	csv = readdlm("ml-latest/fixed/ratings.csv", ',', header=true)
	users = csv[1][:,1]
	movies = csv[1][:,2]
	colours = csv[1][:,3]
	csv = 0

	width = length(union(users))
	height = length(union(movies))

	graph = plot(x=users, y=movies, color=colours, Geom.point,
		Theme(default_point_size=0.1mm, highlight_width=0mm, grid_line_width=0mm),
		Guide.colorkey("Rating"), Guide.xlabel("User"), Guide.ylabel("Movie"),
		Guide.xticks(ticks=collect(1:width), label=false), Guide.yticks(label=false))

	draw(PNG("movie_user_graph.png", width * 0.01mm, height * 0.01mm), graph)
end

main()


text:
function LinearAlgebra.qr!(a::IndicesArray, args...; kwargs...)
    return _qr(qr!(parent(a), args...; kwargs...), axes(a))
end

function _qr(inner::QR, axs::Tuple)
    return QR(
        IndicesArray(inner.factors, axs, AllUnique, LengthChecked),
        inner.τ
       )
end

function Base.parent(fact::QR{<:Any, <:IndicesArray})
    return QR(parent(getfield(fact, :factors)), getfield(fact, :τ))
end

function _qr(inner::LinearAlgebra.QRCompactWY, inds::Tuple)
    return LinearAlgebra.QRCompactWY(
        IndicesArray(inner.factors, inds, AllUnique, LengthChecked),
        inner.T
       )
end
function Base.parent(fact::LinearAlgebra.QRCompactWY{<:Any, <:IndicesArray})
    return LinearAlgebra.QRCompactWY(
        parent(getfield(fact, :factors)),
        getfield(fact, :T)
       )
end

function _qr(inner::QRPivoted, axs::Tuple)
    return QRPivoted(
        IndicesArray(inner.factors, axs, AllUnique, LengthChecked),
        inner.τ,
        inner.jpvt
       )
end
function Base.parent(fact::QRPivoted{<:Any, <:IndicesArray})
    return QRPivoted(
        parent(getfield(fact, :factors)),
        getfield(fact, :τ),
        getfield(fact, :jpvt)
       )
end

function Base.getproperty(fact::IndicesQR, d::Symbol)
    inner = getproperty(parent(fact), d)
    axs = axes(fact)
    if d === :Q
        return IndicesArray(
            inner,
            (unsafe_axes(axs, 1), nothing),
            AllUnique,
            LengthChecked
           )
    elseif d === :R
        return IndicesArray(
            inner,
            (nothing, unsafe_axes(axs, 2)),
            AllUnique,
            LengthChecked
           )
    elseif fact isa QRPivoted && d === :P
        return IndicesArray(inner,
            (unsafe_axes(axs, 1), unsafe_axes(axs, 1)),
            AllUnique,
            LengthChecked
           )
    elseif fact isa QRPivoted && d === :p
        return IndicesArray(
            inner,
            (unsafe_axes(axs, 1),),
            AllUnique,
            LengthChecked
           )
    else
        return inner
    end
end


text:
using LiterateDollop, Test

ϵ = 1e-5

@testset "optfun" begin
    f1(x) = x^2
    f2(x) = (x-1)^2*10 + 10
    f3(x) = (x+10)^2 - 10

    @test isapprox(optfun(f1, -3, 10), 0, atol=ϵ)
    @test isapprox(optfun(f2, -3, 10), 1, atol=ϵ)
    @test isapprox(optfun(f3, -30, 30), -10, atol=ϵ)
end


text:
include("utils.jl")
include("rbm.jl")
include("conditional.jl")
include("nets.jl")
include("io.jl")
include("reporter.jl")
include("test.jl")


text:
module CrossValidators

using Statistics: mean, std

# standard included modules
using DataFrames
using Random
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils

export crossvalidate

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
		       pfunc::Function,nfolds=10)
  ## flatten arrays
  @assert size(X)[1] == length(Y)
  ppl = deepcopy(pl)
  input = X
  output = Y |> Vector{String}
  rowsize = size(input)[1]
  folds = kfold(rowsize,nfolds)
  pacc = Float64[]
  for (fold,trainndx) in enumerate(folds)
    testndx = setdiff(1:rowsize,trainndx)
    trX = input[trainndx,:]
    trY = output[trainndx] |> collect
    tstX = input[testndx,:]
    tstY = output[testndx] |> collect
    res = pipe_accuracy(ppl,pfunc,trX,trY,tstX,tstY)
    push!(pacc,res)
    println("fold: ",fold,", ",res)
  end
  (mean=mean(pacc),std=std(pacc),folds=nfolds)
end

function pipe_accuracy(plearner,perf,trainX,trainY,testX,testY)
  learner = deepcopy(plearner)
  fit!(learner,trainX,trainY)
  pred = transform!(learner,testX)
  perf(pred,testY)
end

end


text:
module CellValuesOperations

using Gridap
using Gridap.Helpers

using Gridap.Kernels: CellSumKernel
using Gridap.Kernels: CellNewAxisKernel

export cellsum
export cellmean
export cellnewaxis
export compress
import Base: +,-,*,/,\, ==, ≈, transpose
import LinearAlgebra: inv, det
import TensorValues: inner, outer, meas

for op in (:+,:-,:inv,:det,:meas,:transpose)
  @eval begin

    function ($op)(m::CellNumber)
      apply($op,m)
    end

    function ($op)(m::CellArray)
      apply($op,m,broadcast=true)
    end

    function ($op)(m::CellMap)
      apply($op,m,broadcast=true)
    end

    function ($op)(m::Map)
      apply($op,m,broadcast=true)
    end

  end
end

for op in (:+,:-,:*,:/,:\,:inner,:outer)
  @eval begin

    function ($op)(a::CellNumber,b::CellNumber)
      apply($op,a,b)
    end

    function ($op)(a::CellArray,b::CellArray)
      apply($op,a,b,broadcast=true)
    end

    function ($op)(a::CellArray,b::CellNumber)
      apply($op,a,b,broadcast=true)
    end

    function ($op)(a::CellNumber,b::CellArray)
      apply($op,a,b,broadcast=true)
    end

    function ($op)(a::CellMap,b::CellMap)
      apply($op,a,b,broadcast=true)
    end

    function ($op)(a::CellMap,b::Number)
      l = length(a)
      apply($op,a,ConstantCellValue(b,l),broadcast=true)
    end

    function ($op)(a::Number,b::CellMap)
      l = length(b)
      apply($op,ConstantCellValue(a,l),b,broadcast=true)
    end

    function ($op)(a::Map,b::Map)
      apply($op,a,b,broadcast=true)
    end

  end
end

for op in (:(==),:≈)
  @eval begin

    function ($op)(a::CellNumber,b::CellNumber)
      _eq_kernel($op,a,b)
    end

    function ($op)(a::CellArray,b::CellArray)
      _eq_kernel($op,a,b)
    end

  end
end

@inline function _eq_kernel(op,a,b)
  length(a) != length(b) && return false
  for (ai,bi) in zip(a,b)
    !(op(ai,bi)) && return false
  end
  return true
end

function cellsum(self::CellArray{T,N};dim::Int) where {T,N}
  k = CellSumKernel{dim}()
  apply(k,self)
end

function cellsum(self::CellArray{T,1};dim::Int) where T
  apply(sum,self)
end

_mean(a) = sum(a)/length(a)

function cellmean(self::CellArray)
  apply(_mean,self)
end

function cellnewaxis(ca::CellArray;dim::Int)
  k = CellNewAxisKernel{dim}()
  apply(k,ca)
end

function cellnewaxis(ca::CellMap;dim::Int)
  k = CellNewAxisKernel{dim}()
  apply(k,ca)
end

function compress(cell_to_vector::CellVector{T}) where T
  n = length(cell_to_vector)
  cell_to_vector_ptrs = zeros(Int,n+1)
  cell_to_vector_data = T[]
  _compress_kernel!(
    cell_to_vector_data,
    cell_to_vector_ptrs,
    cell_to_vector)
  (cell_to_vector_data, cell_to_vector_ptrs)
end

function  _compress_kernel!(
  cell_to_vector_data,
  cell_to_vector_ptrs,
  cell_to_vector)
  for (cell,vector) in enumerate(cell_to_vector)
    cell_to_vector_ptrs[cell+1] = length(vector)
    for v in vector
      push!(cell_to_vector_data,v)
    end
  end
  length_to_ptrs!(cell_to_vector_ptrs)
end

end # module


text:
# If you aren't trying to add to transD_GP but want to use what's
# already available.
module BarPhysics
# completely independent of transD_GP
# you can use whatever packages you want here

# this is probably where you'll store all variables needed for your physics
# e.g., wavenumbers and so on
mutable struct Bar
    x         :: Int
    somefield :: Vector{Float64}
end

# this is what will return the computation from our physics engine
function returnphysics!(F::Bar, m::AbstractArray)
    n = length(F.somefield)
    F.somefield .= F.x*rand(n) .+ sum(m[:]) # pre-allocated assignment
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
"""
    bubble_sorting(array::Array{Int64,1})

The Bubble sorting algorithm (BSA) is a primitive sorting algorithm that repeatedly steps
through the array by using a double for-loop with n and n-1 size. During the walkthrough,
the BSA compares adjacent elements and swaps wrong ordered elements until the array is
sorted. For more information see: [https://en.wikipedia.org/wiki/Bubble_sort](https://en.wikipedia.org/wiki/Bubble_sort)

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [64, 34, 25, 12, 22, 11, 90]
julia> ClassicAlgorithmsCollections.bubble_sorting(arr)
[11, 12, 22, 25, 34, 64, 90]
```
"""
function bubble_sorting(array::Array{Int64,1})
    n = length(array)

    # Traverse through all array elements
    for i in 1:n
        # Last i elements are already in place
        for j in 1:n-i-1
            # traverse the array from 1 to n-i-1
            if array[j] > array[j+1]
                # Swap if the element found is greater than the next element
                array[j], array[j+1] = array[j+1], array[j]
            end
        end
    end
    return array
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
module GaussianVCBenchmark

using InteractiveUtils, LinearAlgebra, Profile, Random
using BenchmarkTools, GLMCopula

Random.seed!(123)

n = 100  # number of observations
ns = rand(100:300, n) # ni in each observation
p = 3   # number of mean parameters
m = 2   # number of variance components
gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
# true parameter values
βtruth = ones(p)
σ2truth = collect(1.:m)
σ02truth = 1.0
for i in 1:n
    ni = ns[i]
    # set up covariance matrix
    V1 = convert(Matrix, Symmetric([Float64(i * (ni - j + 1)) for i in 1:ni, j in 1:ni])) # a pd matrix
    V1 ./= norm(V1) / sqrt(ni) # scale to have Frobenius norm sqrt(n)
    prob = fill(1/ni, ni)
    V2 = ni .* (Diagonal(prob) - prob * transpose(prob))
    V2 ./= norm(V2) / sqrt(ni) # scale to have Frobenious norm sqrt(n)
    Ω = σ2truth[1] * V1 + σ2truth[2] * V2 + σ02truth * I
    Ωchol = cholesky(Symmetric(Ω))
    # simulate design matrix
    X = [ones(ni) randn(ni, p-1)]
    # generate responses
    y = X * βtruth + Ωchol.L * randn(ni)
    # add to data
    gcs[i] = GaussianCopulaVCObs(y, X, [V1, V2])
end

gcm = GaussianCopulaVCModel(gcs)

@info "Initial point:"
init_β!(gcm)
@show gcm.β
fill!(gcm.Σ, 1)
update_Σ!(gcm)
@show gcm.τ
@show gcm.Σ
# @btime update_Σ!(gcm) setup=(fill!(gcm.Σ, 1))

@show loglikelihood!(gcm.data[1], gcm.β, gcm.τ[1], gcm.Σ, true, false)
# @code_warntype loglikelihood!(gcm.data[1], gcm.β, gcm.τ[1], gcm.Σ, true, false)
@btime loglikelihood!(gcm.data[1], gcm.β, gcm.τ[1], gcm.Σ, true, false)

@show loglikelihood!(gcm, true, false)
@show [gcm.∇β; gcm.∇τ; gcm.∇Σ]
# @code_warntype loglikelihood!(gcm, false, false)
# @code_llvm loglikelihood!(gcm, false, false)
@btime loglikelihood!(gcm, true, false)

# Solvers:
# Ipopt.IpoptSolver(print_level=0)
# NLopt.NLoptSolver(algorithm=:LD_MMA, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_CCSAQ, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_SLSQP, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_LBFGS, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_TNEWTON_PRECOND_RESTART, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_TNEWTON_PRECOND, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_TNEWTON_RESTART, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_TNEWTON, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_VAR1, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LD_VAR2, maxeval=4000)
# NLopt.NLoptSolver(algorithm=:LN_COBYLA, maxeval=10000)
# NLopt.NLoptSolver(algorithm=:LN_BOBYQA, maxeval=10000)

@info "MLE:"
solver = Ipopt.IpoptSolver(print_level=0)
@time fit!(gcm, solver)
@show [gcm.β; gcm.τ; gcm.Σ]
@show [gcm.∇β; gcm.∇τ; gcm.∇Σ]
@show loglikelihood!(gcm)
# @btime fit!(gcm, solver) setup=(init_β!(gcm);
#     standardize_res!(gcm); update_quadform!(gcm, true); fill!(gcm.Σ, 1);
#     update_Σ!(gcm))

# Profile.clear()
# @profile begin
for solver in [
    NLopt.NLoptSolver(algorithm=:LN_BOBYQA, maxeval=10000),
    NLopt.NLoptSolver(algorithm=:LN_COBYLA, maxeval=10000),
    NLopt.NLoptSolver(algorithm=:LD_MMA, maxeval=4000),
    NLopt.NLoptSolver(algorithm=:LD_CCSAQ, maxeval=4000),
    NLopt.NLoptSolver(algorithm=:LD_SLSQP, maxeval=4000),
    NLopt.NLoptSolver(algorithm=:LD_LBFGS, maxeval=4000),
    # NLopt.NLoptSolver(algorithm=:LD_TNEWTON_PRECOND_RESTART, maxeval=4000),
    #NLopt.NLoptSolver(algorithm=:LD_TNEWTON_PRECOND, maxeval=4000),
    #NLopt.NLoptSolver(algorithm=:LD_TNEWTON_RESTART, maxeval=4000),
    #NLopt.NLoptSolver(algorithm=:LD_TNEWTON, maxeval=4000),
    # NLopt.NLoptSolver(algorithm=:LD_VAR1, maxeval=4000),
    # NLopt.NLoptSolver(algorithm=:LD_VAR2, maxeval=4000),
    # Ipopt.IpoptSolver(print_level=0)
    ]
    println()
    @show solver
    # re-set starting point
    init_β!(gcm)
    fill!(gcm.Σ, 1)
    update_Σ!(gcm)
    # fit
    fit!(gcm, solver)
    @time fit!(gcm, solver)
    @show loglikelihood!(gcm)
    @show [gcm.β; gcm.τ; gcm.Σ]
    @show [gcm.∇β; gcm.∇τ; gcm.∇Σ]
    println()
end
# end
# Profile.print(format=:flat)

end


text:
#=
Each complex should have collection of cells per dimension:
- cells::Dict{Int,Vector{C}} or Vector{Vector{C}}
=#
abstract type AbstractComplex end

#
# AbstractComplex Public Interface
#
"""Return a complex boundary given element and dimension"""
boundary(cplx::AbstractComplex, i::Integer, d::Int, ::Type{PID}) where {PID} =
    throw(MethodError(boundary, (typeof(cplx),Int,Int,PID)))
boundary(cplx::AbstractComplex, i::Integer, d::Int) = boundary(cplx, i, d, Int)
boundary(cplx::AbstractComplex, cell::AbstractCell, ::Type{PID}) where {PID} =
    boundary(cplx, hash(cell), dim(cell), PID)
boundary(cplx::AbstractComplex, cell::AbstractCell) = boundary(cplx, hash(cell), dim(cell), Int)

"""Return a complex coboundary given element and dimension"""
coboundary(cplx::AbstractComplex, i::Integer, d::Int, ::Type{PID}) where {PID} =
    throw(MethodError(coboundary, (typeof(cplx),Int,Int,PID)))
coboundary(cplx::AbstractComplex, i::Integer, d::Int) = coboundary(cplx, i, d, Int)
coboundary(cplx::AbstractComplex, cell::AbstractCell, ::Type{PID}) where {PID} =
    coboundary(cplx, hash(cell), dim(cell), PID)
coboundary(cplx::AbstractComplex, cell::AbstractCell) = coboundary(cplx, hash(cell), dim(cell), Int)

"""
    faces(cplx::AbstractComplex, cellidx::Integer) -> Vector{Integer}

Return an index collection of faces of a cell element with an indentifier `cellidx` in the complex `cplx`.
"""
faces(cplx::AbstractComplex, cellidx::Integer) = throw(MethodError(faces, (typeof(cplx),Integer)))

"""
    cofaces(cplx::AbstractComplex, cellidx::Integer) -> Vector{Integer}

Return an index collection of cofaces of a cell element with an indentifier `cellidx` in the complex `cplx`.
"""
cofaces(cplx::AbstractComplex, cellidx::Integer) = throw(MethodError(cofaces, (typeof(cplx),Integer)))

"""
Return a complex cell type
"""
eltype(cplx::AbstractComplex) = throw(MethodError(eltype, (typeof(cplx),)))

"""
    cells(cplx::AbstractComplex) -> Dict{Int, AbstractVector{AbstractCell}}

Return a cell collection per dimension (increasing)
"""
cells(cplx::AbstractComplex) = throw(MethodError(cells, (typeof(cplx),)))

"""
    cells(cplx::AbstractComplex, d::Int) -> AbstractVector{AbstractCell}

Return a cell collection for the dimenion `d`.
"""
cells(cplx::AbstractComplex, d::Int) = throw(MethodError(cells, (typeof(cplx),Int)))

"""
    push!(cplx::AbstractComplex, cell::AbstractCell; recursive=false) -> Vector{AbstractCell}

Insert a `cell` to a complex `cplx`, and returns an array of inserted cell(s). If `recursive=true` is passed then all faces of the `cell` are also added to the complex `cplx`.
"""
push!(cplx::AbstractComplex, c::AbstractCell; recursive=false) =
    throw(MethodError(push!, (typeof(cplx),typeof(c))))

#
# Public Methods
#
"""Return a number of cells per dimension"""
size(cplx::AbstractComplex) = (map(length, cells(cplx))...,)

"""Return a dimension of the complex"""
dim(cplx::AbstractComplex) = length(size(cplx))-1

"""Return a total number of cells in the complex"""
length(cplx::AbstractComplex) = sum(size(cplx))

"""Return a number of the cell in the complex of a dimension `d` (0-based)"""
function size(cplx::AbstractComplex, d::Int)
    sz = size(cplx)
    szlen = length(sz)
    (d < 0 || d >= szlen) && return 0
    return sz[d+1]
end

"""
    position(complex, index, dimenion)

Return a position of the cell in an order of cells of the same dimenion of the `complex` given its `index` and `dimenion`.
"""
function position(cplx::AbstractComplex, idx::Integer, d::Int)
    dcells = cells(cplx, d)
    length(dcells) == 0 && return 0
    cidx = findfirst(c->hash(c) == idx, dcells)
    return cidx === nothing ? 0 : cidx
end

"""
    position(complex, cell)

Return a position of the `cell` in an order of cells of the same dimenion of the `complex`.
"""
position(cplx::AbstractComplex, c::AbstractCell) = position(cplx, hash(c), dim(c))

"""
    cplx[idx, d]

Return a `d`-dimensional cell given its index `idx` and dimenion `d`.
"""
function getindex(cplx::AbstractComplex, idx::Integer, d::Int)
    cidx = position(cplx, idx, d)
    cidx == 0 && return nothing
    return cells(cplx, d)[cidx]
end

"""
    boundary(cplx, ch)

Return the chain `ch` boundary in the complex `cplx`.
"""
function boundary(cplx::AbstractComplex, ch::Chain{IX,R}) where {R, IX<:Integer}
    d = dim(ch)
    cc = Chain(d-1, IX, R)
    for (elem, coef) in ch
        append!(cc, coef * boundary(R, cplx[elem, d]))
    end
    return simplify(cc)
end

"""
    coboundary(cplx, ch)

Return the chain `ch` coboundary in the complex `cplx`.
"""
function coboundary(cplx::AbstractComplex, ch::Chain{IX,R}) where {R, IX<:Integer}
    d = dim(ch)
    cc = Chain(d+1, IX, R)
    # δₙ₋₁(c)(σ) = c(∂ₙ(σ))
    for (elem, coef) in ch
        append!(cc, coef * coboundary(cplx, elem, d, R))
    end
    return simplify(cc)
end

"""
    boundary(complex, d, PID)

Generate a boundary matrix from the cell `complex` of the dimension `d` in using coefficients of `PID`.
"""
function boundary(cplx::AbstractComplex, d::Int, ::Type{PID}) where {PID}
    csize = size(cplx)
    rows = d > 0 ? csize[d] : 0
    cols = d <= dim(cplx) ? csize[d+1] : 0
    bm = spzeros(PID, rows, cols)
    if d>=0 && d <= dim(cplx)
        for c in cells(cplx, d)
            i = position(cplx, hash(c), d)
            for (elem, coef) in boundary(PID, c)
                j = position(cplx, elem, d-1)
                bm[j, i] = coef
            end
        end
    end
    return bm
end

"""
    cell in cplx

Checks if the `cell` is in the complex `cplx`
"""
in(c::AbstractCell, cplx::AbstractComplex) = position(cplx, c) > 0

"""
    cochain(complex, d, coefficients)

Return a cochain of the dimension `d` for the `complex` with PID `coefficients`.
"""
function cochain(cplx::AbstractComplex, d::Int, coefs::Vector{PID}) where {PID}
    cs = cells(cplx, d)
    @assert length(cs) == length(coefs) "Number of coefficients must match number of cells of dimension $d"
    return Chain(d, coefs, map(c->hash(c), cs))
end
cochain(cplx::AbstractComplex, d::Int, coef::PID) where {PID} =
    cochain(cplx, d, fill(coef, size(cplx, d)))

"""
    adjacency_matrix(cplx, [T=Int])

Construct an adjacency matrix of type `T` from a 1-skeleton (1D subcomplex) of the complex `cplx`.
"""
function adjacency_matrix(cplx::AbstractComplex, ::Type{T}) where {T<:Real}
    C0 = map(hash, cells(cplx, 0))
    N = length(C0)
    adj = spzeros(T,N,N)
    for c in cells(cplx, 1)
        i, j = map(h->findfirst(isequal(h), C0), vertices(c))
        adj[i, j] = adj[j, i] = one(T)
    end
    return adj
end
adjacency_matrix(cplx::AbstractComplex) = adjacency_matrix(cplx, Int)

"""
    showchain(cplx::AbstractComplex, ch::AbstractChain)

Display the chain `ch` with elements from the complex `cplx`.
"""
function showchain(cplx::AbstractComplex, ch::AbstractChain)
    d = dim(ch)
    R = valtype(ch)
    pos = one(R)
    print("[$d]: ")
    if iszero(ch)
        print("0")
    else
        for (i,(id,v)) in enumerate(ch.cells)
            val = abs(v)
            if sign(v) == pos
                i != 1 && print(" + ")
            else
                print(" - ")
            end
            splx = strip(repr("text/plain", cplx[id, d]), ['\"'])
            print("$val⋅$splx")
        end
    end
end


text:
mutable struct Problem{T<:Real}
    head::Symbol
    objective::AbstractExpr
    constraints::Array{Constraint}
    status::MOI.TerminationStatusCode
    model::Union{MOI.ModelLike,Nothing}

    function Problem{T}(
        head::Symbol,
        objective::AbstractExpr,
        constraints::Array = Constraint[],
    ) where {T<:Real}
        if sign(objective) == Convex.ComplexSign()
            error("Objective cannot be a complex expression")
        else
            return new(
                head,
                objective,
                constraints,
                MOI.OPTIMIZE_NOT_CALLED,
                nothing,
            )
        end
    end
end

function Base.getproperty(p::Problem, s::Symbol)
    if s === :optval
        if getfield(p, :status) == MOI.OPTIMIZE_NOT_CALLED
            return nothing
        else
            return objective_value(p)
        end
    else
        return getfield(p, s)
    end
end

dual_status(p::Problem) = MOI.get(p.model, MOI.DualStatus())
primal_status(p::Problem) = MOI.get(p.model, MOI.PrimalStatus())
termination_status(p::Problem) = MOI.get(p.model, MOI.TerminationStatus())
objective_value(p::Problem) = MOI.get(p.model, MOI.ObjectiveValue())

Problem(args...) = Problem{Float64}(args...)

function vexity(p::Problem)
    bad_vex = [ConcaveVexity, NotDcp]

    obj_vex = vexity(p.objective)
    if p.head == :maximize
        obj_vex = -obj_vex
    end
    typeof(obj_vex) in bad_vex &&
        @warn "Problem not DCP compliant: objective is not DCP"

    constr_vex = ConstVexity()
    for i in 1:length(p.constraints)
        vex = vexity(p.constraints[i])
        typeof(vex) in bad_vex &&
            @warn "Problem not DCP compliant: constraint $i is not DCP"
        constr_vex += vex
    end
    problem_vex = obj_vex + constr_vex
    # this check is redundant
    # typeof(problem_vex) in bad_vex && warn("Problem not DCP compliant")
    return problem_vex
end

function conic_form!(p::Problem, unique_conic_forms::UniqueConicForms)
    objective_var = Variable()
    objective = conic_form!(objective_var, unique_conic_forms)
    conic_form!(p.objective - objective_var == 0, unique_conic_forms)
    for constraint in p.constraints
        conic_form!(constraint, unique_conic_forms)
    end
    return objective, objective_var.id_hash
end

function Problem{T}(
    head::Symbol,
    objective::AbstractExpr,
    constraints::Constraint...,
) where {T<:Real}
    return Problem{T}(head, objective, [constraints...])
end

# Allow users to simply type minimize
function minimize(
    objective::AbstractExpr,
    constraints::Constraint...;
    numeric_type = Float64,
)
    return Problem{numeric_type}(:minimize, objective, collect(constraints))
end
function minimize(
    objective::AbstractExpr,
    constraints::Array{<:Constraint} = Constraint[];
    numeric_type = Float64,
)
    return Problem{numeric_type}(:minimize, objective, constraints)
end
function minimize(
    objective::Value,
    constraints::Constraint...;
    numeric_type = Float64,
)
    return minimize(
        convert(AbstractExpr, objective),
        collect(constraints);
        numeric_type = numeric_type,
    )
end
function minimize(
    objective::Value,
    constraints::Array{<:Constraint} = Constraint[];
    numeric_type = Float64,
)
    return minimize(
        convert(AbstractExpr, objective),
        constraints;
        numeric_type = numeric_type,
    )
end

# Allow users to simply type maximize
function maximize(
    objective::AbstractExpr,
    constraints::Constraint...;
    numeric_type = Float64,
)
    return Problem{numeric_type}(:maximize, objective, collect(constraints))
end
function maximize(
    objective::AbstractExpr,
    constraints::Array{<:Constraint} = Constraint[];
    numeric_type = Float64,
)
    return Problem{numeric_type}(:maximize, objective, constraints)
end
function maximize(
    objective::Value,
    constraints::Constraint...;
    numeric_type = Float64,
)
    return maximize(
        convert(AbstractExpr, objective),
        collect(constraints);
        numeric_type = numeric_type,
    )
end
function maximize(
    objective::Value,
    constraints::Array{<:Constraint} = Constraint[];
    numeric_type = Float64,
)
    return maximize(
        convert(AbstractExpr, objective),
        constraints;
        numeric_type = numeric_type,
    )
end

# Allow users to simply type satisfy (if there is no objective)
function satisfy(constraints::Constraint...; numeric_type = Float64)
    return Problem{numeric_type}(:minimize, Constant(0), [constraints...])
end
function satisfy(
    constraints::Array{<:Constraint} = Constraint[];
    numeric_type = Float64,
)
    return Problem{numeric_type}(:minimize, Constant(0), constraints)
end
function satisfy(constraint::Constraint; numeric_type = Float64)
    return satisfy([constraint]; numeric_type = numeric_type)
end

# +(constraints, constraints) is defined in constraints.jl
function add_constraints!(p::Problem, constraints::Array{<:Constraint})
    return append!(p.constraints, constraints)
end
function add_constraints!(p::Problem, constraint::Constraint)
    return add_constraints!(p, [constraint])
end

function add_constraint!(p::Problem, constraints::Array{<:Constraint})
    return add_constraints!(p, constraints)
end
function add_constraint!(p::Problem, constraint::Constraint)
    return add_constraints!(p, constraint)
end

# caches conic form of x when x is the solution to the optimization problem p
function cache_conic_form!(
    conic_forms::UniqueConicForms,
    x::AbstractExpr,
    p::Problem,
)
    objective = conic_form!(p.objective, conic_forms)
    for c in p.constraints
        conic_form!(c, conic_forms)
    end
    return cache_conic_form!(conic_forms, x, objective)
end


text:
function classify(cur_string::AbstractString)
  cur_string = capitalize(cur_string)

  cur_string = join(map(x -> ucfirst(x), split(cur_string, r"[- _]")))

  cur_string
end


text:
function input()
    ccall(:);
end


text:
"""
    truncated(d0::UnivariateDistribution; [lower::Real], [upper::Real])
    truncated(d0::UnivariateDistribution, lower::Real, upper::Real)

A _truncated distribution_ `d` of a distribution `d0` to the interval
``[l, u]=```[lower, upper]` has the probability density (mass) function:

```math
f(x; d_0, l, u) = \\frac{f_{d_0}(x)}{P_{Z \\sim d_0}(l \\le Z \\le u)}, \\quad x \\in [l, u],
```
where ``f_{d_0}(x)`` is the probability density (mass) function of ``d_0``.

The function throws an error if ``l > u``.

```julia
truncated(d0; lower=l)           # d0 left-truncated to the interval [l, Inf)
truncated(d0; upper=u)           # d0 right-truncated to the interval (-Inf, u]
truncated(d0; lower=l, upper=u)  # d0 truncated to the interval [l, u]
truncated(d0, l, u)              # d0 truncated to the interval [l, u]
```

The function falls back to constructing a [`Truncated`](@ref) wrapper.

# Implementation

To implement a specialized truncated form for distributions of type `D`, one or more of the
following methods should be implemented:
- `truncated(d0::D, l::T, u::T) where {T <: Real}`: interval-truncated
- `truncated(d0::D, ::Nothing, u::Real)`: right-truncated
- `truncated(d0::D, l::Real, u::Nothing)`: left-truncated
"""
function truncated end
function truncated(d::UnivariateDistribution, l::Real, u::Real)
    return truncated(d, promote(l, u)...)
end
function truncated(
    d::UnivariateDistribution;
    lower::Union{Real,Nothing}=nothing,
    upper::Union{Real,Nothing}=nothing,
)
    return truncated(d, lower, upper)
end
function truncated(d::UnivariateDistribution, ::Nothing, u::Real)
    # (log)ucdf = (log)tp = (log) P(X ≤ u) where X ~ d
    logucdf = logtp = logcdf(d, u)
    ucdf = tp = exp(logucdf)

    Truncated(d, promote(oftype(float(u), -Inf), u, zero(ucdf), ucdf, tp, logtp)...)
end
function truncated(d::UnivariateDistribution, l::Real, ::Nothing)
    # (log)lcdf = (log) P(X < l) where X ~ d
    loglcdf = if value_support(typeof(d)) === Discrete
        logsubexp(logcdf(d, l), logpdf(d, l))
    else
        logcdf(d, l)
    end
    lcdf = exp(loglcdf)

    # (log)tp = (log) P(l ≤ X) where X ∼ d
    logtp = log1mexp(loglcdf)
    tp = exp(logtp)

    Truncated(d, promote(l, oftype(float(l), Inf), lcdf, one(lcdf), tp, logtp)...)
end
truncated(d::UnivariateDistribution, ::Nothing, ::Nothing) = d
function truncated(d::UnivariateDistribution, l::T, u::T) where {T <: Real}
    l <= u || error("the lower bound must be less or equal than the upper bound")

    # (log)lcdf = (log) P(X < l) where X ~ d
    loglcdf = if value_support(typeof(d)) === Discrete
        logsubexp(logcdf(d, l), logpdf(d, l))
    else
        logcdf(d, l)
    end
    lcdf = exp(loglcdf)

    # (log)ucdf = (log) P(X ≤ u) where X ~ d
    logucdf = logcdf(d, u)
    ucdf = exp(logucdf)

    # (log)tp = (log) P(l ≤ X ≤ u) where X ∼ d
    logtp = logsubexp(loglcdf, logucdf)
    tp = exp(logtp)

    Truncated(d, promote(l, u, lcdf, ucdf, tp, logtp)...)
end

"""
    Truncated

Generic wrapper for a truncated distribution.
"""
struct Truncated{D<:UnivariateDistribution, S<:ValueSupport, T <: Real} <: UnivariateDistribution{S}
    untruncated::D      # the original distribution (untruncated)
    lower::T      # lower bound
    upper::T      # upper bound
    lcdf::T       # cdf of lower bound (exclusive): P(X < lower)
    ucdf::T       # cdf of upper bound (inclusive): P(X ≤ upper)

    tp::T         # the probability of the truncated part, i.e. ucdf - lcdf
    logtp::T      # log(tp), i.e. log(ucdf - lcdf)
    function Truncated(d::UnivariateDistribution, l::T, u::T, lcdf::T, ucdf::T, tp::T, logtp::T) where {T <: Real}
        new{typeof(d), value_support(typeof(d)), T}(d, l, u, lcdf, ucdf, tp, logtp)
    end
end

### Constructors of `Truncated` are deprecated - users should call `truncated`
@deprecate Truncated(d::UnivariateDistribution, l::Real, u::Real) truncated(d, l, u)

params(d::Truncated) = tuple(params(d.untruncated)..., d.lower, d.upper)
partype(d::Truncated) = partype(d.untruncated)
Base.eltype(::Type{Truncated{D, S, T} } ) where {D, S, T} = T

### range and support

islowerbounded(d::Truncated) = islowerbounded(d.untruncated) || isfinite(d.lower)
isupperbounded(d::Truncated) = isupperbounded(d.untruncated) || isfinite(d.upper)

minimum(d::Truncated) = max(minimum(d.untruncated), d.lower)
maximum(d::Truncated) = min(maximum(d.untruncated), d.upper)

function insupport(d::Truncated{D,<:Union{Discrete,Continuous}}, x::Real) where {D<:UnivariateDistribution}
    return d.lower <= x <= d.upper && insupport(d.untruncated, x)
end

### evaluation

quantile(d::Truncated, p::Real) = quantile(d.untruncated, d.lcdf + p * d.tp)

function pdf(d::Truncated, x::Real)
    result = pdf(d.untruncated, x) / d.tp
    return d.lower <= x <= d.upper ? result : zero(result)
end

function logpdf(d::Truncated, x::Real)
    result = logpdf(d.untruncated, x) - d.logtp
    return d.lower <= x <= d.upper ? result : oftype(result, -Inf)
end

function cdf(d::Truncated, x::Real)
    result = (cdf(d.untruncated, x) - d.lcdf) / d.tp
    return if x < d.lower
        zero(result)
    elseif x >= d.upper
        one(result)
    else
        result
    end
end

function logcdf(d::Truncated, x::Real)
    result = logsubexp(logcdf(d.untruncated, x), log(d.lcdf)) - d.logtp
    return if x < d.lower
        oftype(result, -Inf)
    elseif x >= d.upper
        zero(result)
    else
        result
    end
end

function ccdf(d::Truncated, x::Real)
    result = (d.ucdf - cdf(d.untruncated, x)) / d.tp
    return if x <= d.lower
        one(result)
    elseif x > d.upper
        zero(result)
    else
        result
    end
end

function logccdf(d::Truncated, x::Real)
    result = logsubexp(logccdf(d.untruncated, x), log1p(-d.ucdf)) - d.logtp
    return if x <= d.lower
        zero(result)
    elseif x > d.upper
        oftype(result, -Inf)
    else
        result
    end
end

## random number generation

function _rand!(rng::AbstractRNG, d::Truncated)
    d0 = d.untruncated
    if d.tp > 0.25
        while true
            r = _rand!(rng, d0)
            if d.lower <= r <= d.upper
                return r
            end
        end
    else
        return quantile(d0, d.lcdf + rand(rng) * d.tp)
    end
end

## show

function show(io::IO, d::Truncated)
    print(io, "Truncated(")
    d0 = d.untruncated
    uml, namevals = _use_multline_show(d0)
    uml ? show_multline(io, d0, namevals) :
          show_oneline(io, d0, namevals)
    if d.lower > -Inf
        if d.upper < Inf
            print(io, "; lower=$(d.lower), upper=$(d.upper))")
        else
            print(io, "; lower=$(d.lower))")
        end
    elseif d.upper < Inf
        print(io, "; upper=$(d.upper))")
    else
        print(io, ")")
    end
    uml && println(io)
end

_use_multline_show(d::Truncated) = _use_multline_show(d.untruncated)

### specialized truncated distributions

include(joinpath("truncated", "normal.jl"))
include(joinpath("truncated", "exponential.jl"))
include(joinpath("truncated", "uniform.jl"))
include(joinpath("truncated", "loguniform.jl"))


text:
module ParallelRun

import FITSIO
import JLD

import ..Log
using ..Model
import ..SDSSIO
import ..Infer
import ..SDSSIO: RunCamcolField

#set this to false to use source-division parallelism
const SKY_DIVISION_PARALLELISM=true

const TILE_WIDTH = 20
const MIN_FLUX = 2.0

# Use distributed parallelism (with Dtree)
if haskey(ENV, "USE_DTREE") && ENV["USE_DTREE"] != ""
    const Distributed = true
    using Dtree
else
    const Distributed = false
    const dt_nodeid = 1
    const dt_nnodes = 1
    DtreeScheduler(n, f) = ()
    initwork(dt) = 0, (1, 0)
    getwork(dt) = 0, (1, 0)
    runtree(dt) = 0
    cpu_pause() = ()
end

# Use threads (on the loop over sources)
const Threaded = true
if Threaded && VERSION > v"0.5.0-dev"
    using Base.Threads
else
    # Pre-Julia 0.5 there are no threads
    nthreads() = 1
    threadid() = 1
    macro threads(x)
        x
    end
    SpinLock() = 1
    lock(l) = ()
    unlock(l) = ()
end

# A workitem is of this ra / dec size
const wira = 0.025
const widec = 0.025

"""
Timing information.
"""
type InferTiming
    query_fids::Float64
    num_infers::Int64
    read_photoobj::Float64
    read_img::Float64
    fit_psf::Float64
    opt_srcs::Float64
    num_srcs::Int64
    write_results::Float64
    wait_done::Float64

    InferTiming() = new(0.0, 0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0)
end

function add_timing!(i::InferTiming, j::InferTiming)
    i.query_fids = i.query_fids + j.query_fids
    i.num_infers = i.num_infers + j.num_infers
    i.read_photoobj = i.read_photoobj + j.read_photoobj
    i.read_img = i.read_img + j.read_img
    i.fit_psf = i.fit_psf + j.fit_psf
    i.opt_srcs = i.opt_srcs + j.opt_srcs
    i.num_srcs = i.num_srcs + j.num_srcs
    i.write_results = i.write_results + j.write_results
    i.wait_done = i.wait_done + j.wait_done
end

immutable BoundingBox
    ramin::Float64
    ramax::Float64
    decmin::Float64
    decmax::Float64
end

@inline nputs(nid, s) = ccall(:puts, Cint, (Ptr{Int8},), string("[$nid] ", s))
@inline phalse(b) = b[] = false

include("source_division_inference.jl")

"""
Given a BoundingBox that is to be divided into `nra` x `ndec` subareas,
return the `i`th subarea. `i` is a linear index between 1 and
`nra * ndec`.

This function assumes a cartesian (rather than spherical) coordinate system!
"""
function divide_skyarea(box, nra, ndec, i)
    global wira, widec
    ix, iy = ind2sub((nra, ndec), i)

    return (box.ramin + (ix - 1) * wira,
            min(box.ramin + ix * wira, box.ramax),
            box.decmin + (iy - 1) * widec,
            min(box.decmin + iy * widec, box.decmax))
end

function time_puts(elapsedtime, bytes, gctime, allocs)
    s = @sprintf("%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units),
                            Int64(1024))
        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units),
                            Int64(1000))
        if ma == 1
            s = string(s, @sprintf(" (%d%s allocation%s: ", allocs,
                                   Base._cnt_units[ma], allocs==1 ? "" : "s"))
        else
            s = string(s, @sprintf(" (%.2f%s allocations: ", allocs,
                                   Base._cnt_units[ma]))
        end
        if mb == 1
            s = string(s, @sprintf("%d %s%s", bytes,
                                   Base._mem_units[mb], bytes==1 ? "" : "s"))
        else
            s = string(s, @sprintf("%.3f %s", bytes, Base._mem_units[mb]))
        end
        if gctime > 0
            s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
        end
        s = string(s, ")")
    elseif gctime > 0
        s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
    end
    nputs(dt_nodeid, s)
end

"""
Divide `N` into `np` parts as evenly as possible, returning `my` part as
a (first, last) tuple.
"""
function divparts(N, np, my)
    len, rem = divrem(N, np)
    if len == 0
        if my > rem
            return 1, 0
        end
        len, rem = 1, 0
    end
    # compute my part
    f = 1 + ((my-1) * len)
    l = f + len - 1
    # distribute remaining evenly
    if rem > 0
        if my <= rem
            f = f + (my-1)
            l = l + my
        else
            f = f + rem
            l = l + rem
        end
    end
    return f, l
end

"""
Divide the given ra, dec range into sky areas of `wira`x`widec` and
use Dtree to distribute these sky areas to nodes. Within each node
use `one_node_infer()` to fit the Celeste model to sources in each sky area.
"""
function divide_sky_and_infer(
                box::BoundingBox,
                stagedir::String;
                timing=InferTiming(),
                outdir=".")
    if dt_nodeid == 1
        nputs(dt_nodeid, "running on $dt_nnodes nodes")
    end

    # how many `wira` X `widec` sky areas (work items)?
    global wira, widec
    nra = ceil(Int64, (box.ramax - box.ramin) / wira)
    ndec = ceil(Int64, (box.decmax - box.decmin) / widec)

    num_work_items = nra * ndec
    each = ceil(Int64, num_work_items / dt_nnodes)

    if dt_nodeid == 1
        nputs(dt_nodeid, "work item dimensions: $wira X $widec")
        nputs(dt_nodeid, "$num_work_items work items, ~$each per node")
    end

    # create Dtree and get the initial allocation
    dt, is_parent = DtreeScheduler(num_work_items, 0.4)
    ni, (ci, li) = initwork(dt)
    rundt = Ref(runtree(dt))
    @inline function rundtree(again)
        if again[]
            again[] = runtree(dt)
            cpu_pause()
        end
        again[]
    end

    # work item processing loop
    nputs(dt_nodeid, "initially $ni work items ($ci to $li)")
    itimes = InferTiming()
    while ni > 0
        li == 0 && break
        if ci > li
            nputs(dt_nodeid, "consumed allocation (last was $li)")
            ni, (ci, li) = getwork(dt)
            nputs(dt_nodeid, "got $ni work items ($ci to $li)")
            continue
        end
        item = ci
        ci = ci + 1

        # map item to subarea
        iramin, iramax, idecmin, idecmax = divide_skyarea(box, nra, ndec, item)

        # Get vector of (run, camcol, field) triplets overlapping this patch
        tic()
        box = BoundingBox(iramin, iramax, idecmin, idecmax)
        rcfs = get_overlapping_fields(box, stagedir)
        itimes.query_fids = toq()

        # run inference for this subarea
        results = one_node_infer(rcfs, stagedir;
                        box=BoundingBox(iramin, iramax, idecmin, idecmax),
                        reserve_thread=rundt,
                        thread_fun=rundtree,
                        timing=itimes)
        tic()
        save_results(outdir, iramin, iramax, idecmin, idecmax, results)
        itimes.write_results = toq()

        timing.num_infers = timing.num_infers+1
        add_timing!(timing, itimes)
        rundtree(rundt)
    end
    nputs(dt_nodeid, "out of work")
    tic()
    while rundt[]
        rundtree(rundt)
    end
    finalize(dt)
    timing.wait_done = toq()
end

function load_images(rcfs, stagedir)
    images = TiledImage[]
    image_names = String[]
    image_count = 0

    for i in 1:length(rcfs)
        Log.info("reading field $(rcfs[i])")
        rcf = rcfs[i]
        field_images = SDSSIO.load_field_images(rcf, stagedir)
        for b=1:length(field_images)
            image_count += 1
            push!(image_names,
                "$image_count run=$(rcf.run) camcol=$(rcf.camcol) field=$(rcf.field) b=$b")
            tiled_image = TiledImage(field_images[b])
            push!(images, tiled_image)
        end
    end
    gc()

    Log.debug("Image names:")
    Log.debug(string(image_names))

    images
end

"""
Use mulitple threads on one node to
fit the Celeste model to sources in a given bounding box.

- rcfs: Array of run, camcol, field triplets that the source occurs in.
- box: a bounding box specifying a region of sky

Returns:

- Dictionary of results, keyed by SDSS thing_id.
"""
function one_node_infer(
               rcfs::Vector{RunCamcolField},
               stagedir::String;
               objid="",
               box=BoundingBox(-1000., 1000., -1000., 1000.),
               primary_initialization=true,
               reserve_thread=Ref(false),
               thread_fun=phalse,
               timing=InferTiming())
    nprocthreads = nthreads()
    if reserve_thread[]
        nprocthreads = nprocthreads-1
    end
    Log.info("Running with $(nprocthreads) threads")

    # Read all primary objects in these fields.
    tic()
    duplicate_policy = primary_initialization ? :primary : :first
    catalog = SDSSIO.read_photoobj_files(rcfs, stagedir,
                        duplicate_policy=duplicate_policy)
    timing.read_photoobj = toq()
    Log.info("$(length(catalog)) primary sources")

    reserve_thread[] && thread_fun(reserve_thread)

    # Filter out low-flux objects in the catalog.
    catalog = filter(entry->(maximum(entry.star_fluxes) >= MIN_FLUX), catalog)
    Log.info("$(length(catalog)) primary sources after MIN_FLUX cut")

    # Filter any object not specified, if an objid is specified
    if objid != ""
        Log.info(catalog[1].objid)
        catalog = filter(entry->(entry.objid == objid), catalog)
    end

    # Get indicies of entries in the  RA/Dec range of interest.
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_sources = find(entry_in_range, catalog)

    nputs(dt_nodeid, string("processing $(length(target_sources)) sources in ",
          "$(box.ramin), $(box.ramax), $(box.decmin), $(box.decmax)"))

    # If there are no objects of interest, return early.
    if length(target_sources) == 0
        return Dict{Int, Dict}()
    end

    reserve_thread[] && thread_fun(reserve_thread)

    # Read in images for all (run, camcol, field).
    tic()

    images = load_images(rcfs, stagedir)
    timing.read_img = toq()

    reserve_thread[] && thread_fun(reserve_thread)

    Log.info("finding neighbors")
    tic()
    neighbor_map = Infer.find_neighbors(target_sources, catalog, images)
    Log.info("neighbors found in $(toq()) seconds")

    reserve_thread[] && thread_fun(reserve_thread)

    # iterate over sources
    curr_source = 1
    last_source = length(target_sources)
    sources_lock = SpinLock()
    results = Dict[]
    results_lock = SpinLock()
    function process_sources()
        tid = threadid()

        if reserve_thread[] && tid == 1
            while reserve_thread[]
                thread_fun(reserve_thread)
                cpu_pause()
            end
        else
            while true
                lock(sources_lock)
                ts = curr_source
                curr_source += 1
                unlock(sources_lock)

                if ts > last_source
                    break
                end
#                try
                    s = target_sources[ts]
                    entry = catalog[s]
                    nputs(dt_nodeid, "processing source $s: objid = $(entry.objid)")

                    t0 = time()
                    # TODO: subset images to images_local too.
                    vs_opt = Infer.infer_source(images,
                                                catalog[neighbor_map[ts]],
                                                entry)
                    runtime = time() - t0
#                catch ex
#                    Log.error(ex)
#                end

                lock(results_lock)
                push!(results, Dict(
                    "thing_id"=>entry.thing_id,
                    "objid"=>entry.objid,
                    "ra"=>entry.pos[1],
                    "dec"=>entry.pos[2],
                    "vs"=>vs_opt,
                    "runtime"=>runtime))
                unlock(results_lock)
            end
        end
    end

    tic()
    ccall(:jl_threading_run, Void, (Any,), Core.svec(process_sources))
    ccall(:jl_threading_profile, Void, ())
    timing.opt_srcs = toq()
    timing.num_srcs = length(target_sources)

    results
end

"""
Query the SDSS database for all fields that overlap the given RA, Dec range.
"""
function get_overlapping_field_extents(query::BoundingBox, stagedir::String)
    f = FITSIO.FITS("$stagedir/field_extents.fits")

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}
    all_ramin = read(hdu, "ramin")::Vector{Float64}
    all_ramax = read(hdu, "ramax")::Vector{Float64}
    all_decmin = read(hdu, "decmin")::Vector{Float64}
    all_decmax = read(hdu, "decmax")::Vector{Float64}

    close(f)

    ret = Tuple{RunCamcolField, BoundingBox}[]

    # The ramin, ramax, etc is a bit unintuitive because we're looking
    # for any overlap.
    for i in eachindex(all_ramin)
        if (all_ramax[i] > query.ramin && all_ramin[i] < query.ramax &&
                all_decmax[i] > query.decmin && all_decmin[i] < query.decmax)
            cur_box = BoundingBox(all_ramin[i], all_ramax[i],
                                  all_decmin[i], all_decmax[i])
            cur_fe = RunCamcolField(all_run[i], all_camcol[i], all_field[i])
            push!(ret, (cur_fe, cur_box))
        end
    end

    return ret
end

"""
Like `get_overlapping_fields`, but return a Vector of
(run, camcol, field) triplets.
"""
function get_overlapping_fields(query::BoundingBox, stagedir::String)
    fes = get_overlapping_field_extents(query, stagedir)
    [fe[1] for fe in fes]
end

"""
called from main entry point for inference for one field
(used for accuracy assessment, infer-box is the primary inference
entry point)
"""
function infer_field(rcf::RunCamcolField,
                     stagedir::String,
                     outdir::String;
                     objid="")
    results = one_node_infer([rcf,], stagedir; objid=objid, primary_initialization=false)
    fname = if objid == ""
        @sprintf "%s/celeste-%06d-%d-%04d.jld" outdir rcf.run rcf.camcol rcf.field
    else
        @sprintf "%s/celeste-objid-%s.jld" outdir objid
    end
    JLD.save(fname, "results", results)
    Log.debug("infer_field finished successfully")
end

"""
Save provided results to a JLD file.
"""
function save_results(outdir, ramin, ramax, decmin, decmax, results)
    fname = @sprintf("%s/celeste-%.4f-%.4f-%.4f-%.4f.jld",
                     outdir, ramin, ramax, decmin, decmax)
    JLD.save(fname, "results", results)
end

save_results(outdir, box, results) =
    save_results(outdir, box.ramin, box.ramax, box.decmin, box.decmax, results)

"""
called from main entry point.
"""
function infer_box(box::BoundingBox, stagedir::String, outdir::String)
    # Base.@time hack for distributed environment
    gc_stats = ()
    gc_diff_stats = ()
    elapsed_time = 0.0
    gc_stats = Base.gc_num()
    elapsed_time = time_ns()

    times = InferTiming()
    if !SKY_DIVISION_PARALLELISM
        Log.debug("source division parallelism")
        divide_sources_and_infer(box, stagedir; timing=times, outdir=outdir)
    elseif dt_nnodes > 1
        Log.debug("sky division parallelism")
        divide_sky_and_infer(box, stagedir; timing=times, outdir=outdir)
    else
        Log.debug("multithreaded parallelism only")
        tic()
        # Get vector of (run, camcol, field) triplets overlapping this patch
        rcfs = get_overlapping_fields(box, stagedir)
        times.query_fids = toq()

        results = one_node_infer(rcfs, stagedir; box=box, timing=times)

        tic()
        save_results(outdir, box, results)
        times.write_results = toq()
    end

    # Base.@time hack for distributed environment
    elapsed_time = time_ns() - elapsed_time
    gc_diff_stats = Base.GC_Diff(Base.gc_num(), gc_stats)
    time_puts(elapsed_time, gc_diff_stats.allocd, gc_diff_stats.total_time,
              Base.gc_alloc_count(gc_diff_stats))

    times.num_srcs = max(1, times.num_srcs)
    nputs(dt_nodeid, "timing: query_fids=$(times.query_fids)")
    nputs(dt_nodeid, "timing: num_infers=$(times.num_infers)")
    nputs(dt_nodeid, "timing: read_photoobj=$(times.read_photoobj)")
    nputs(dt_nodeid, "timing: read_img=$(times.read_img)")
    nputs(dt_nodeid, "timing: fit_psf=$(times.fit_psf)")
    nputs(dt_nodeid, "timing: opt_srcs=$(times.opt_srcs)")
    nputs(dt_nodeid, "timing: num_srcs=$(times.num_srcs)")
    nputs(dt_nodeid, "timing: average opt_srcs=$(times.opt_srcs/times.num_srcs)")
    nputs(dt_nodeid, "timing: write_results=$(times.write_results)")
    nputs(dt_nodeid, "timing: wait_done=$(times.wait_done)")
end

end


text:
__precompile__()

module FastStyleTransfer

using Flux, Metalhead, Images, CUDA, NNlib, BSON

CUDA.allowscalar(false)

export train, stylize

include("utils.jl")
include("layers.jl")
include("vgg.jl")
include("transformer_net.jl")
include("neural_style.jl")

end # module


text:
#----------- Basic Interaction --------------#

"Returns the flag with given symbol."
function Base.getindex(input::DFInput, n::Symbol)
    if haskey(input.flags, n)
        return input.flags[n]
    else
        error("Flag:$n Not found in flags of Input $(name(input))")
    end
end
"Sets the flag"
Base.setindex!(input::DFInput, dat, key) = setflags!(input, key => dat)

"""
    setflags!(input::DFInput, flags...; print=true)

Sets the specified flags in the input.
"""
function setflags!(input::DFInput{T}, flags...; print=true) where T
    found_keys = Symbol[]
    for (flag, value) in flags
        flag_type = flagtype(input, flag)
        if flag_type != Nothing
            !(flag in found_keys) && push!(found_keys, flag)
            try
                if isa(value, AbstractVector{<:AbstractVector}) && flag_type <: AbstractVector
                    value = [convert.(eltype(flag_type), v) for v in value]
                else
                    value = convert(flag_type, value)
                end
            catch
                print && ( @warn "Filename '$(name(input))':\n  Could not convert '$value' into '$flag_type'.\n    Flag '$flag' not set.\n")
                continue
            end
            old_data = haskey(input.flags, flag) ? input.flags[flag] : ""
            input.flags[flag] = value
            print && (@info "$(name(input)):\n  -> $flag:\n      $old_data set to: $value\n")
        end
    end
    return found_keys, input
end

"""
    rmflags!(input::DFInput, flags...)

Remove the specified flags.
"""
function rmflags!(input::DFInput, flags...; print=true)
    for flag in flags
        if haskey(input.flags, flag)
            pop!(input.flags, flag, false)
            print && (@info "Removed flag '$flag' from input '$(name(input))'")
        end
    end
    return input
end

data(input::DFInput, n) = getfirst(x-> name(x) == n, data(input))

"""
    setdata!(input::DFInput, block_name::Symbol, new_block_data; option=nothing, print=true)

sets the data of the specified 'InputData' to the new data. Optionally also sets the 'InputData' option.
"""
function setdata!(input::DFInput, block_name::Symbol, new_block_data; option=nothing, print=true)
    setd = false
    for data_block in input.data
        if data_block.name == block_name
            if typeof(data_block.data) != typeof(new_block_data)
                if print @warn "Overwritten data of type '$(typeof(data_block.data))' with type '$(typeof(new_block_data))'." end
            end
            old_data        = data_block.data
            data_block.data = new_block_data
            data_block.option = option == nothing ? data_block.option : option
            if print
                @info "Block data '$(data_block.name)' in input  '$(name(input))' is now:\n\t$(string(data_block.data)) \n\toption: $(data_block.option)\n"
            end
            setd = true
        end
    end
    if !setd
        setdata!(input, InputData(block_name, option, new_block_data))
        setd = true
    end
    return setd, input
end

"""
    setdataoption!(input::DFInput, name::Symbol, option::Symbol; print=true)

Sets the option of specified data.
"""
function setdataoption!(input::DFInput, name::Symbol, option::Symbol; print=true)
    for data in input.data
        if data.name == name
            old_option  = data.option
            data.option = option
            print && (@info "Option of InputData '$(data.name)' in input '$(name(input))' set from '$old_option' to '$option'")
        end
    end
    return input
end

execs(input::DFInput, exec::String) =
	filter(x -> occursin(exec, x.exec), input.execs)

exec(input::DFInput, exec::String) =
	getfirst(x -> occursin(exec, x.exec), input.execs)

execflags(input::DFInput, exec::String) =
	[x.exec => x.flags for x in execs(input, exec)]

function setexecflags!(input::DFInput, exec::String, flags...)
	for e in execs(input, exec)
		setflags!(e, flags...)
	end
end

setexecdir!(input::DFInput, exec, dir) =
	setexecdir!.(execs(input, exec), dir)

rmexecflags!(input::DFInput, exec::String, flags...) =
	rmflags!.(execs(input, exec), flags...)

runcommand(input::DFInput) = input.execs[1]

"Returns the outputdata for the input."
function outputdata(input::DFInput; print=true, overwrite=true)
    if hasoutput(input) && !overwrite
        return outdata(input)
    end
    if hasoutfile(input)
        input.outdata = readoutput(input)
        return input.outdata
    end
    print && (@warn "No output data or output file found for input: $(name(input)).")
    return SymAnyDict()
end
#---------- Extended Interaction ------------#

"""
    setkpoints!(input::DFInput, k_grid)

Sets the kpoints of the input. Will automatically generate the kgrid values if necessary.
"""
function setkpoints!(input::DFInput{Wannier90}, k_grid::NTuple{3, Int}; print=true)
    setflags!(input, :mp_grid => [k_grid...], print=print)
    setdata!(input, :kpoints, kgrid(k_grid..., :wan), print=print)
    return input
end

function setkpoints!(input::DFInput{QE}, k_grid::NTuple{3, Int}; print=true) #nscf

    calc = flag(input, :calculation)
    print && calc != "nscf" && (@warn "Expected calculation to be 'nscf'.\nGot $calc.")
    setdata!(input, :k_points, kgrid(k_grid..., :nscf), option = :crystal, print=print)
    prod(k_grid) > 100 && setflags!(input, :verbosity => "high", print=print)
    return input
end

function setkpoints!(input::DFInput{QE}, k_grid::NTuple{6, Int}; print=true) #scf
    calc = flag(input, :calculation)
    print && calc != "scf" && !occursin("relax", calc) && (@warn "Expected calculation to be scf, vc-relax, relax.\nGot $calc.")
    setdata!(input, :k_points, [k_grid...], option = :automatic, print=print)
    prod(k_grid[1:3]) > 100 && setflags!(input, :verbosity => "high", print=print)
    return input
end

function setkpoints!(input::DFInput{QE}, k_grid::Vector{NTuple{4, T}}; print=true, k_option=:crystal_b) where T<:AbstractFloat
    calc = flag(input, :calculation)
    print && calc != "bands" && (@warn "Expected calculation to be bands, got $calc.")
    @assert in(k_option, [:tpiba_b, :crystal_b, :tpiba_c, :crystal_c]) error("Only $([:tpiba_b, :crystal_b, :tpiba_c, :crystal_c]...) are allowed as a k_option, got $k_option.")
    if k_option in [:tpiba_c, :crystal_c]
        @assert length(k_grid) == 3 error("If $([:tpiba_c, :crystal_c]...) is selected the length of the k_points needs to be 3, got length: $(length(k_grid)).")
    end
    num_k = 0.0
    for k in k_grid
        num_k += k[4]
    end
    if num_k > 100.
        setflags!(input, :verbosity => "high", print=print)
        if print
            @info "Verbosity is set to high because num_kpoints > 100,\n
                       otherwise bands won't get printed."
        end
    end
    setdata!(input, :k_points, k_grid, option=k_option, print=print)
    return input
end

"""
    readbands(input::DFInput)

Tries to read bands from the outputfile of the input.
"""
function readbands(input::DFInput)
    hasoutput_assert(input)
    to = readoutput(input)
    if haskey(to, :bands)
        return to[:bands]
    else
        error("No bands found in $(name(input)).")
    end
end

"""
    readfermi(input::DFInput)

Tries to read the fermi level from the outputfile of the input.
"""
function readfermi(input::DFInput)
    hasoutput_assert(input)
    to = readoutput(input)
    if haskey(to, :fermi)
        return to[:fermi]
    else
        error("No fermi found in $(name(input)).")
    end
end

"""
    setwanenergies!(waninput::DFInput{Wannier90}, structure::AbstractStructure, nscf::DFInput , Emin::Real; Epad=5.0)

Automatically calculates and sets the wannier energies. This uses the projections,
`Emin` and the output of the nscf calculation to infer the other limits.
`Epad` allows one to specify the padding around the inner and outer energy windows
"""
function setwanenergies!(input::DFInput{Wannier90}, structure::AbstractStructure, nscf::DFInput, Emin::Real; Epad=5.0)
    hasoutput_assert(nscf)
    iscalc_assert(nscf, "nscf")
    nbnd = isnoncolincalc(nscf) ? 2 * nprojections(structure) : nprojections(structure)
    print && (@info "num_bands=$nbnd (inferred from provided projections).")
    bands = readbands(nscf)
    winmin, frozmin, frozmax, winmax = wanenergyranges(Emin, nbnd, bands, Epad)
    setflags!(input, :dis_win_min => winmin, :dis_froz_min => frozmin, :dis_froz_max => frozmax, :dis_win_max => winmax, :num_wann => nbnd, :num_bands=>length(bands))
    return input
end

#-------- Generating new DFInputs ---------- #

function input_from_kpoints(template::DFInput, newname, kpoints, newflags...)
    newcalc = DFInput(template, newname, newflags...)
    setkpoints!(newcalc, kpoints, print=false)
    return newcalc
end

"""
    gencalc_scf(template::DFInput, kpoints::NTuple{6, Int}, newflags...; name="scf")

Uses the information from the template and supplied kpoints to generate an scf input.
Extra flags can be supplied which will be set for the generated input.
"""
gencalc_scf(template::DFInput, kpoints::NTuple{6, Int}, newflags...; name="scf") =
    input_from_kpoints(template, name, kpoints, :calculation => "scf", newflags...)

"""
    gencalc_bands(template::DFInput, kpoints::Vector{NTuple{4}}, newflags...; name="bands")

Uses the information from the template and supplied kpoints to generate a bands input.
Extra flags can be supplied which will be set for the generated input.
"""
gencalc_bands(template::DFInput, kpoints::Vector{<:NTuple{4}}, newflags...; name="bands") =
    input_from_kpoints(template, name, kpoints, :calculation => "bands", newflags...)

"""
    gencalc_nscf(template::DFInput, kpoints::NTuple{3, Int}, newflags...; name="nscf")

Uses the information from the template and supplied kpoints to generate an nscf input.
Extra flags can be supplied which will be set for the generated input.
"""
gencalc_nscf(template::DFInput, kpoints::NTuple{3, Int}, newflags...; name="nscf") =
    input_from_kpoints(template, name, kpoints, :calculation => "nscf", newflags...)

"""
    gencalc_projwfc(template::DFInput, Emin, Emax, DeltaE, newflags...; name="projwfc")

Uses the information from the template and supplied kpoints to generate a projwfc.x input.
Extra flags can be supplied which will be set for the generated input.
"""
function gencalc_projwfc(template::DFInput, Emin, Emax, DeltaE, extraflags...; name="projwfc")
    occflag = flag(template, :occupations)
    ngauss  = 0
    if occflag == "smearing"
        smearingflag = flag(template, :smearing)
        if smearingflag ∈ ("methfessel-paxton", "m-p", "mp")
            ngauss = 1
        elseif smearingflag ∈ ("marzari-vanderbilt", "cold", "m-v", "mv")
            ngauss = -1
        elseif smearingflag ∈ ("fermi-dirac", "f-d", "fd")
            ngauss = -99
        end
    end
    tdegaussflag = flag(template, :degauss)
    degauss = tdegaussflag != nothing ? tdegaussflag : 0.0
    if length(execs(template)) == 2
        excs = [execs(template)[1], Exec("projwfc.x", execs(template)[end].dir)]
    else
        excs = [Exec("projwfc.x", execs(template)[end].dir)]
    end

    out = DFInput(template, name, excs=excs)
    empty!(out.flags)
    setflags!(out, :Emin => Emin, :Emax => Emax, :DeltaE => DeltaE,
              :ngauss => ngauss, :degauss => degauss, print=false)
    setflags!(out, extraflags...)
    return out
end

"""
    gencalc_wan(nscf::DFInput{QE}, structure::AbstractStructure, Emin;
                Epad     = 5.0,
                wanflags = nothing,
                wanexec  = Exec("wannier90.x", ""))

Generates a Wannier90 input to follow on the supplied `nscf` calculation. It uses the projections defined in the `structure`, and starts counting the required amount of bands from `Emin`.
The `nscf` needs to have a valid output since it will be used in conjunction with `Emin` to find the required amount of bands and energy window for the Wannier90 calculation.
"""
function gencalc_wan(nscf::DFInput{QE}, structure::AbstractStructure, Emin;
                     Epad     = 5.0,
                     wanflags = nothing,
                     wanexec  = Exec("wannier90.x", ""))

    hasoutput_assert(nscf)
    iscalc_assert(nscf, "nscf")
    hasprojections_assert(structure)
    if iscolincalc(nscf)
        wannames = ["wanup", "wandn"]
        @info "Spin polarized calculation found (inferred from nscf input)."
    else
        wannames = ["wan"]
    end

    if flag(nscf, :nosym) != true
        @info "'nosym' flag was not set in the nscf calculation.
                If this was not intended please set it and rerun the nscf calculation.
                This generally gives errors because of omitted kpoints, needed for pw2wannier90.x"
    end

    nbnd = isnoncolincalc(nscf) ? 2 * nprojections(structure) : nprojections(structure)
    @info "num_bands=$nbnd (inferred from provided projections)."

    bands = readbands(nscf)
    wanflags = wanflags != nothing ? SymAnyDict(wanflags) : SymAnyDict()
    wanflags[:dis_win_min], wanflags[:dis_froz_min], wanflags[:dis_froz_max], wanflags[:dis_win_max] = wanenergyranges(Emin, nbnd, bands, Epad)

    wanflags[:num_bands] = length(bands)
    wanflags[:num_wann]  = nbnd
    kpoints = data(nscf, :k_points).data
    wanflags[:mp_grid] = kakbkc(kpoints)
    wanflags[:preprocess] = true
    @info "mp_grid=$(join(wanflags[:mp_grid]," ")) (inferred from nscf input)."

	isnoncolincalc(nscf) && (wanflags[:spinors] = true)

    kdata = InputData(:kpoints, :none, [k[1:3] for k in kpoints])

    waninputs = DFInput{Wannier90}[]
    for  wanfil in wannames
        push!(waninputs, DFInput{Wannier90}(wanfil, dir(nscf), copy(wanflags), [kdata], [Exec(), wanexec], true))
    end

    if length(waninputs) > 1
        setflags!(waninputs[1], :spin => "up")
        setflags!(waninputs[2], :spin => "down")
    end
    return waninputs
end

"""
    gencalc_wan(nscf::DFInput{QE}, structure::AbstractStructure, projwfc::DFInput{QE}, threshold::Real; kwargs...)

Generates a wannier calculation, that follows on the `nscf` calculation. Instead of passing Emin manually, the output of a projwfc.x run
can be used together with a `threshold` to determine the minimum energy such that the contribution of the
projections to the DOS is above the `threshold`.
"""
function gencalc_wan(nscf::DFInput{QE}, structure::AbstractStructure, projwfc::DFInput{QE}, threshold::Real; kwargs...)
    hasexec_assert(projwfc, "projwfc.x")
    hasoutput_assert(projwfc)
    Emin = Emin_from_projwfc(structure, projwfc, threshold)
    gencalc_wan(nscf, structure, Emin; kwargs...)
end

"""
    isconverged(input::DFInput{QE})

Returns whether an `scf` calculation was converged.
"""
function isconverged(input::DFInput{QE})
    hasoutput_assert(input)
    iscalc_assert(input, "scf")
    return outputdata(input)[:converged]
end


text:
# Multivariate Lazy Kernel Density Estimator ====================================================================

module MultiKDE

export KDEMulti, KDEUniv, pdf, gpke
export DimensionType, ContinuousDim, CategoricalDim, UnorderedCategoricalDim

using Distributions

abstract type  DimensionType end
struct ContinuousDim <: DimensionType end
struct CategoricalDim <: DimensionType
    level::Int
end
struct UnorderedCategoricalDim <: DimensionType
    level::Int
end

include("kernel.jl")
include("kde.jl")

end


text:
using GRIB

keylist = Vector{String}()
grib_file = "/home/neo/Projekte/Predly/weather/gfs_workflow/download/data/gfs.20210109/06/gfs.t06z.pgrb2.0p25.f002"

parameters = [7,11,586,422,423,424,426,428,430,433,434,435,437,438,439,440,441,442,443,451,453,455,456,457,458,472,473,474,476,479,480,481,482,483,495,496,498,588]

unneeded_identifiers = ["values", "codedValues", "latLonValues", "latitudes", "longitudes", "distinctLatitudes", "distinctLongitudes"]

f = GribFile(grib_file)

function getValue(f, band::Int64, lat::Float64, lon::Float64)

    start_time = time()

    for (index, message) in enumerate(f)

        band_time = time()
        data_array = []
        parameter_name = message["parameterName"]
        parameter_unit = message["parameterUnits"]
        parameter_ecmf = message["nameECMF"]
        units_ecmf     = message["unitsECMF"]
        name           = message["name"]
        units          = message["units"]
        level          = message["level"]
        type_of_level  = message["typeOfLevel"]
        short_name     = message["shortName"]

        if index == band

            Nearest(message) do near
                lons, lats, values, distances = find(near, message, 8, 50)
                seek(f, 0)
                println("Took ", time()-start_time)
                println(lons,lats,values,distances)
                return lons, lats, values, distances
            end
        end
    end

end

getValue(f, 450, 50.0, 8.25)

destroy(f)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

_hmax(γ::Variogram) = 3range(γ)
_hmax(γ::PowerVariogram) = 3.0
_hmax(γ::NuggetEffect) = 3.0

@recipe function f(γ::Variogram, minlag=0, maxlag=_hmax(γ); nlags=100)
  # start at 1e-6 instead of 0 to avoid nugget artifact in plot
  h = range(minlag+1e-6, stop=maxlag, length=nlags)

  seriestype --> :path
  label --> string(nameof(typeof(γ)))
  xlims --> (0, maxlag)
  ylims --> (0, Inf)
  xguide --> "h"
  yguide --> "γ(h)"

  h, γ.(h)
end