text:
# Counts of discrete values

#################################################
#
#  counts on given levels
#
#################################################

typealias IntRange1{T<:Integer} Range1{T}

#### functions for counting a single list of integers (1D)

function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntRange1)
	# add counts of integers from x to r

	k = length(levels)
	length(r) == k || raise_dimerror()

	m0 = levels[1]
	m1 = levels[end]
	b = m0 - 1

	@inbounds for i in 1 : length(x)
		xi = x[i]
		if m0 <= xi <= m1
			r[xi - b] += 1
		end
	end
	return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntRange1, wv::WeightVec)
	k = length(levels)
	length(r) == k || raise_dimerror()

	m0 = levels[1]
	m1 = levels[end]
	b = m0 - 1
	w = values(wv)

	@inbounds for i in 1 : length(x)
		xi = x[i]
		if m0 <= xi <= m1
			r[xi - b] += w[i]
		end
	end
	return r
end

counts(x::IntegerArray, levels::IntRange1) = addcounts!(zeros(Int, length(levels)), x, levels)
counts(x::IntegerArray, levels::IntRange1, wv::WeightVec) = addcounts!(zeros(eltype(wv), length(levels)), x, levels, wv)
counts(x::IntegerArray, k::Integer) = counts(x, 1:k)
counts(x::IntegerArray, k::Integer, wv::WeightVec) = counts(x, 1:k, wv)

proportions(x::IntegerArray, levels::IntRange1) = counts(x, levels) .* inv(length(x))
proportions(x::IntegerArray, levels::IntRange1, wv::WeightVec) = counts(x, levels, wv) .* inv(sum(wv))
proportions(x::IntegerArray, k::Integer) = proportions(x, 1:k)
proportions(x::IntegerArray, k::Integer, wv::WeightVec) = proportions(x, 1:k, wv)

#### functions for counting a single list of integers (2D)

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1))
	# add counts of integers from x to r

	n = length(x)
	length(y) == n || raise_dimerror()

	xlevels, ylevels = levels

	kx = length(xlevels)
	ky = length(ylevels)
	size(r) == (kx, ky) || raise_dimerror()

	mx0 = xlevels[1]
	mx1 = xlevels[end]
	my0 = ylevels[1]
	my1 = ylevels[end]

	bx = mx0 - 1
	by = my0 - 1

	for i = 1:n
		xi = x[i]
		yi = y[i]
		if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
			r[xi - bx, yi - by] += 1
		end
	end
	return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1), wv::WeightVec)
	# add counts of integers from x to r

	n = length(x)
	length(y) == length(wv) == n || raise_dimerror()

	xlevels, ylevels = levels

	kx = length(xlevels)
	ky = length(ylevels)
	size(r) == (kx, ky) || raise_dimerror()

	mx0 = xlevels[1]
	mx1 = xlevels[end]
	my0 = ylevels[1]
	my1 = ylevels[end]

	bx = mx0 - 1
	by = my0 - 1
	w = values(wv)

	for i = 1:n
		xi = x[i]
		yi = y[i]
		if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
			r[xi - bx, yi - by] += w[i]
		end
	end
	return r
end

# facet functions

function counts(x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1))
	addcounts!(zeros(Int, length(levels[1]), length(levels[2])), x, y, levels)
end

function counts(x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1), wv::WeightVec)
	addcounts!(zeros(eltype(wv), length(levels[1]), length(levels[2])), x, y, levels, wv)
end

counts(x::IntegerArray, y::IntegerArray, levels::IntRange1) = counts(x, y, (levels, levels))
counts(x::IntegerArray, y::IntegerArray, levels::IntRange1, wv::WeightVec) = counts(x, y, (levels, levels), wv)

counts(x::IntegerArray, y::IntegerArray, ks::(Integer, Integer)) = counts(x, y, (1:ks[1], 1:ks[2]))
counts(x::IntegerArray, y::IntegerArray, ks::(Integer, Integer), wv::WeightVec) = counts(x, y, (1:ks[1], 1:ks[2]), wv)
counts(x::IntegerArray, y::IntegerArray, k::Integer) = counts(x, y, (1:k, 1:k))
counts(x::IntegerArray, y::IntegerArray, k::Integer, wv::WeightVec) = counts(x, y, (1:k, 1:k), wv)

proportions(x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1)) = counts(x, y, levels) .* inv(length(x))
proportions(x::IntegerArray, y::IntegerArray, levels::(IntRange1, IntRange1), wv::WeightVec) = counts(x, y, levels, wv) .* inv(sum(wv))

proportions(x::IntegerArray, y::IntegerArray, ks::(Integer, Integer)) = proportions(x, y, (1:ks[1], 1:ks[2]))
proportions(x::IntegerArray, y::IntegerArray, ks::(Integer, Integer), wv::WeightVec) = proportions(x, y, (1:ks[1], 1:ks[2]), wv)
proportions(x::IntegerArray, y::IntegerArray, k::Integer) = proportions(x, y, (1:k, 1:k))
proportions(x::IntegerArray, y::IntegerArray, k::Integer, wv::WeightVec) = proportions(x, y, (1:k, 1:k), wv)

#################################################
#
#  countmap on unknown levels
#
#  These methods are based on dictionaries, and
#  can be used on any kind of hashable values.
#
#################################################

## auxiliary functions

function _normalize_countmap{T}(cm::Dict{T}, s::Real)
	r = Dict{T,Float64}()
	for (k, c) in cm
		r[k] = c / s
	end
	return r
end

## 1D

function addcounts!{T}(cm::Dict{T}, x::AbstractArray{T})
	for v in x
		cm[v] = get(cm, v, 0) + 1
	end
	return cm
end

function addcounts!{T,W}(cm::Dict{T}, x::AbstractArray{T}, wv::WeightVec{W})
	n = length(x)
	length(wv) == n || raise_dimerror()
	w = values(wv)
	z = zero(W)

	for i = 1 : n
		@inbounds xi = x[i]
		@inbounds wi = w[i]
		cm[xi] = get(cm, xi, z) + wi
	end
	return cm
end

countmap{T}(x::AbstractArray{T}) = addcounts!(Dict{T,Int}(), x)
countmap{T,W}(x::AbstractArray{T}, wv::WeightVec{W}) = addcounts!(Dict{T,W}(), x, wv)

proportionmap(x::AbstractArray) = _normalize_countmap(countmap(x), length(x))
proportionmap(x::AbstractArray, wv::WeightVec) = _normalize_countmap(countmap(x, wv), sum(wv))


text:
struct GPModel{TGP,Tx,Ty} <: AbstractMCMC.AbstractModel
    gp::TGP
    x::Tx
    y::Ty
end

struct GPSampler{TS} <: AbstractMCMC.AbstractSampler
    sampler::TS
end

function StatsBase.sample(model::MCGP, N::Int; kwargs...)
    return sample(Random.GLOBAL_RNG, model, N; kwargs...)
end
function StatsBase.sample(rng::Random.AbstractRNG, model::MCGP, N::Int; kwargs...)
    sampler = inference(model)
    thinning = get(kwargs, :thinning, sampler.thinning)
    discard_initial = get(kwargs, :discard_initial, sampler.nBurnin)
    progressname = get(kwargs, :progressname, "Sampling with $sampler")
    return sample(
        rng,
        GPModel(model, input(model.data), output(model.data)),
        GPSampler(sampler),
        N;
        thinning,
        discard_initial,
        progressname,
        kwargs...,
    )
end
function AbstractMCMC.step(
    rng::AbstractRNG, model::GPModel, sampler::GPSampler{<:GibbsSampling}; kwargs...
)
    state = compute_kernel_matrices(model.gp, (;), model.x, true)
    local_vars = init_local_vars(likelihood(model.gp), length(model.x))
    local_vars = sample_local!(local_vars, likelihood(model.gp), model.y, means(model.gp))
    state = merge(state, (; local_vars))
    f =
        sample_global!.(
            ∇E_μ(model.gp, model.y, state.local_vars),
            ∇E_Σ(model.gp, model.y, state.local_vars),
            Zviews(model.gp),
            getf(model.gp),
            state.kernel_matrices,
        )
    sampler.sampler.n_iter += 1
    return f, state
end

function AbstractMCMC.step(
    rng::AbstractRNG, model::GPModel, sampler::GPSampler{<:GibbsSampling}, state; kwargs...
)
    local_vars = sample_local!(
        state.local_vars, likelihood(model.gp), model.y, means(model.gp)
    )
    state = merge(state, (; local_vars))
    f =
        sample_global!.(
            ∇E_μ(model.gp, model.y, state.local_vars),
            ∇E_Σ(model.gp, model.y, state.local_vars),
            Zviews(model.gp),
            getf(model.gp),
            state.kernel_matrices,
        )
    sampler.sampler.n_iter += 1
    return f, state
end

function AbstractMCMC.bundle_samples(
    samples, model::GPModel, sampler::GPSampler, state, chain_type::Type; kwargs...
)
    resume = get(kwargs, :cat, true)
    if !resume || isempty(sampler.sampler.sample_store) # check if either one should restart sampling or if no samples was ever taken
        sampler.sampler.sample_store = samples
    else
        sampler.sampler.sample_store = vcat(sampler.sampler.sample_store, samples)
    end
    return samples
end

ext: jl


text:
using JSON

#####################
# Work in progress!
#####################

# set of tools for listing packages based on their popularity (by github stars)
# and exporting all info we get hold of about the items in the package to JSON
# before generating a big "glossary" html file

type PkgInfo
	name::String
	url::String
	stars::Int
end

# using Requests
function findstars()

	function tourl(ssh::String)
		url = replace(ssh, r"^git", "https")
		replace(url, r".git$", "")
	end

	function geturls()
		pkgs = PkgInfo[]
		metadir = Pkg.dir("METADATA")
		plist = readdir(metadir)
		for p in plist
			urlfname = joinpath(metadir, p, "url")
			if isfile(urlfname)
				ssh = strip(open(readall, urlfname, "r"), ['\n', ' ', '\t'])
				url = tourl(ssh)
				push!(pkgs, PkgInfo(p, url, 0))
			end
	    end
	    return pkgs
	end

	function getstars(pkg::PkgInfo, attempt=0)
		if attempt > 10
			println("max attempts reached")
			return
		end
		m = nothing
		response = nothing
		try
			response = get(pkg.url)
			m = match(r"stargazers\">.*\n *(\d+)", response.data)
		catch
			println(pkg.url, ", error on get request!, trying again")
			return getstars(pkg, attempt + 1)
		end
		if response.headers["status_code"] == "301" && haskey(response.headers, "Location")
			newurl = response.headers["Location"]
			println("redirection to $newurl")
			pkg.url = newurl
			return getstars(pkg, attempt + 1)
		end
		if m == nothing
			println("\n url: $(pkg.url)")
			@show response
			println("no Stars found!")
			return
		end
		pkg.stars = int(m.captures[1])
		println(pkg.stars, " stars")
	end

	#> old and very naive way of getting list of packages:
	# url = "http://julia.readthedocs.org/en/latest/packages/packagelist/"
	# println("Getting package list from $url...")
	# html = get(url).data
	# pkgs = PkgInfo[]
	# offset = 1
	# i = 0
	# while true
	# 	m = match(r"
[\n ]*(.*?)", html, offset)
	# 	m == nothing ? break : ()
	# 	offset = m.offset + 1
	# 	push!(pkgs, PkgInfo(m.captures[2], m.captures[1], 0))
	# end
	pkgs = geturls()
	println("looking up $(length(pkgs)) packages...")
	for (i, pkg) in enumerate(pkgs)
		print("$i looking up $(pkg.name)... ")
		getstars(pkg)
		# i > 20 ? break : ()
	end
	sort!(pkgs, by=i -> i.stars, rev=true)
	println("\nPackage Stars:")
	for pkg in pkgs
		println(rpad("$(pkg.name):", 30), pkg.stars)
	end
	pkgs
end

function listpopular()
	jsonfile = "pkgs.json"
	if isfile(jsonfile)
		pkgs = JSON.parsefile(jsonfile)
	else
		println("finding popular packages")
		pkgs = findstars()
		f=open(jsonfile, "w")
		JSON.print(f, pkgs, 2)
		close(f)
	end
end
# listpopular()

type ItemInfo
	name::String
	summary::String
	help::String
	methods::Array{String, 1}
end

function moduleinfo(m::Module, existingmethods::Set)
	# we want to ignore methods already defined in other modules
	# but still show them if they've only be been previously seen
	# in this module
	function checkexists(s::String)
		h = hash(s)
		in(h, existingmethods) ? (return false) : ()
		push!(existingmethods_new, h)
		return true
	end
	existingmethods_new = copy(existingmethods)
	info = ItemInfo[]
	for v in sort(names(m))
		s = string(v)
		if isdefined(m,v)
			i = eval(m,v)
			sumstr = summary(i)
			helpstr = ""
			sio = IOBuffer()
			help(sio, s)
			helpstr = takebuf_string(sio)
			close(sio)
			mlist = String[]
			if isa(i, Function)
				try
					mlist = map(string, methods(eval(m,v)))
					mlist = filter(checkexists, mlist)
				catch
				end
			end
			push!(info, ItemInfo(s, sumstr, helpstr, mlist))
		end
	end
	existingmethods = existingmethods_new
	return info
end

function getinfo(infofile)
	pkgs = [
		"Base",
		"Gadfly",
		"DataFrames",
		"PyCall",
		"GeneticAlgorithms",
		"Optim",
		"JuMP",
		"Morsel",
		"GLM",
		"Winston"
	]
	for pkg in pkgs
		pkg == "Base" ? continue: ()
		if Pkg.installed(pkg) == nothing
			println("Adding packages $pkg ...")
			Pkg.add(pkg)
		end
	end
	allinfo = Any[]
	methods = Set()
	for mod in pkgs
		println("processing $mod ...")
		eval(parse("using $mod"))
		modinfo = ["name"=>string(mod), "items"=> moduleinfo(eval(parse(mod)), methods)]
		push!(allinfo, modinfo)
	end
	f=open(infofile, "w")
	JSON.print(f, allinfo, 2)
	close(f)
end
# getinfo("allinfo.json")

using Mustache
using HttpCommon
function create_page()
	infofile = "allinfo.json"
	if !isfile(infofile)
		getinfo(infofile)
	end
	info = JSON.parsefile(infofile)
	template = open(readall, "glossary.template.html", "r")
	i = 1
	for mod in info
		mod["modname"] = pop!(mod, "name")
		for item in mod["items"]
			if length(item["methods"]) == 0
				delete!(item, "methods")
			else
				mesc = map(escapeHTML, item["methods"])
				mstr = join(mesc, "
\n
")
				item["methods"] = "
$mstr
"
			end
			item["id"] = i
			i += 1
		end
	end
	context = ["modules"=>info]
	page = render(template, context)
	open("glossary.html","w") do f
		print(f, page)
	end
	open("gmodules.js", "w") do f
		print(f, "modules = ")
		JSON.print(f, info)
	end
end
create_page()

max_stars_repo_name: lkuper/JuliaByExample

max_issues_repo_name: svaksha/JuliaByExample

max_forks_repo_name: svaksha/JuliaByExample


text:
InputFolder = './Images/EyeFundus/';
OutputFolder = './Results/Images/EyeFundusEnhance/';
Lbl = 1;
Fill = -1;

@iA = '*.tif';

@fxg_gEnhanceVessels2D [iA] > [F];
params.Scales = [1 1.5 2];
params.Spacings = [1 1];
params.Tau = 1;
params.Polarity = 0;
params.Gamma = 1;
/endf

/show iA>;
/show F >;
/%keep F > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
import LightXML
export peaks, shift, PeakStruct, ToXML

mutable struct PeakStruct
    positions::Array{Float64,1}
    heights::Array{Float64,1}
    widths::Array{Float64,1}
    intensities::Array{Float64,1}
    deconvolution::Data1D
end

"""
	shift(p::PeakStruct,delta::Real)

Shift the positions of all peaks in `p` by delta
"""
function shift(p::PeakStruct, delta::Real)
	np=deepcopy(p)
	np.positions.+=delta
	return np
end

"""
`smooth(d::Data1D,n=10)`
returns a smoothed version of the input data set `d`, using
convolution with a Gaussian peak.
"""
function smooth(d::Data1D,n=10)
      l=length(d.dat);
      kern=[exp(-5*j^2/n^2) for j=-n:n]
      kern /= sum(kern)
      return(Data1D( [sum([d.dat[k+j]*kern[j+n+1] for j=-n:n]) for k=(n+1):(l-n)],pos2ind(d,n+1),pos2ind(d,l-n-1)));
end

lorentzian(x0::Float64,σ::Float64,x::Float64) = 1.0/π*sqrt(σ)/(1.0+σ*(x-x0)^2)

@doc raw"""
    clorentzian(x0::Float64,σ::Float64,x::Float64)

returns the complex Lorentzian

``L(x) = \sqrt{\sigma}\dfrac{1+i\sqrt{\sigma}(x-x_0)}{1+\sigma(x-x_0)^2}``
"""
clorentzian(x0::Float64,σ::Float64,x::Float64) =  1.0/π*sqrt(σ)*(1.0+sqrt(σ)*im*(x-x0))/(1.0+σ*(x-x0)^2)

"""
`peaks(dinput::Data1D;threshold=1,athresh=1,regions=128)`
identifies the peaks in `dinput`. Peak positions are located roughly
by searching for the downward zero crossings of the first derivative.
The data points surrounding each peak are then used for a quadratic regression,
from which the position, curvature, and height of the maximum are interpolated.
The signal `dinput` is then decomposed into a linear superposition of Lorentzian
peaks at the determined positions and with the determined curvatures.

`peaks()` returns a data structure
```
struct PeakStruct
    positions::Array{Float64,1}
    heights::Array{Float64,1}
    widths::Array{Float64,1}
    intensities::Array{Float64,1}
    deconvolution::Data1D
  end
```

where `positions` and `heights` are the interpolated peak positions and heights,
respectively; `widths` are the peak curvatures, which correspond
to 1/FWHM line widths of the Lorentzian peaks. `intensities` are the
peak intensities resulting from deconvolution. These values are
proportional to the underlying concentrations in an NMR spectrum.
Finally, `deconvolution` contains the deconvoluted approximation to the
original input signal.
"""
function peaks(dinput::Data1D;threshold=1,athresh=1,regions=128)
    d=smooth(dinput)
    dr=derivative(d);
    n=length(d.dat)
    δ=(d.istop-d.istart)/n
    lreg=floor(Int64,n/regions)

    # find smallest standard deviation in regions
    mindev=minimum(abs.([std(dinput.dat[k:(k+lreg)]) for k=1:lreg:(n-lreg)]))
    mindevd=minimum(abs.([std(dr.dat[k:(k+lreg)]) for k=1:lreg:(n-lreg)]))

    # locate indices of peak positions roughly on the basis of the derivative
    # and original data set
    ind=10:1:(length(dr.dat)-10);
    peaks=filter(k->dr.dat[k-1]>0 && dr.dat[k]<0 && dr.dat[k-1]-dr.dat[k] > athresh*mindevd && dinput.dat[k]>mindev*threshold, ind);

    # For each peak, perform a second order polynomial fit to the
    # peak tip. More accurate locations and equivalent peak width
    # information is extracted from the resulting coefficients

    dpeaks=[ind2pos(dinput,pos2ind(dr,k)) for k in peaks] # a hack to eliminate aligment differences
    coeffs = [NMR.polyfit(δ*collect(-9:9),dinput.dat[(p-9):(p+9)],2)
                    for p in dpeaks]

    # compute updated peak positions, intensities, and lorentzian equivalent
    # line widths from the expansion coefficients

    σ=[abs(C[3]) for C in coeffs] ;    # inverse line widths
    cpks=[pos2ind(dr,k) for k in peaks]+[-0.5*C[2]/C[3] for C in coeffs]; # corrected peak positions
    hpks=[C[1]-0.25*C[2]*C[2]/C[3] for C in coeffs]; # corrected peak heights

    # deconvolute the spectrum with Lorentzian peaks of the determined
    # positions and widths

    B=collect(lorentzian(cpks[k],0.66*σ[k],x) for x in NMR.ind(dinput), k in 1:length(cpks))  ;
    c=pinv(B)* NMR.val(dinput);

    return(PeakStruct(cpks,hpks,σ,c,Data1D(B*c,dinput.istart,dinput.istop)))
end
"""

	ToXML(p)

generates an XML structure containing the peak intensities, positions, and widths
contained in the `PeakStruct` `p`.
"""
function ToXML(p::PeakStruct)
	root=LightXML.new_element("PeakStruct")
	for k=1:length(p.positions)
		e=LightXML.new_child(root, "Peak")
		LightXML.add_text( LightXML.new_child(e,"Intensity"), "$(p.intensities[k])" )
		LightXML.add_text( LightXML.new_child(e,"Position"), "$(p.positions[k])" )
		LightXML.add_text( LightXML.new_child(e,"Width"), "$(p.widths[k])" )
	end
	return root
end

"""

	 FromXML(e::XMLElement)

converts XML representation to `PeakStruct`
"""
function FromXML(e::LightXML.XMLElement)
	int=Array(Float64,1)([])
	pos=similar(int)
	wid=similar(int)

# we assume that the name of e is "PeakStruct", we don't check it here anymore
	for k in e["Peak"]
		push!(pos,content(k["Position"]));
		push!(int,content(k["Intensity"]));
		push!(wid,content(k["Width"]));
	end

	return PeakStruct(pos,[],wid,int,Data1D([],[]))
end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
# Test Trace

using Turing

import Turing: Trace, current_trace, fork, VarName, Sampler

if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :n)
  n[] = 0
else
  const n = Ref(0)
end

alg = PG(5, 1)
spl = Turing.Sampler(alg)
dist = Normal(0, 1)

function f2()
  t = TArray(Int, 1);
  t[1] = 0;
  while true
    ct = current_trace()
    vn = VarName(gensym(), :x, "[$n]", 1)
    Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
    produce(t[1]);
    vn = VarName(gensym(), :x, "[$n]", 1)
    Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
    t[1] = 1 + t[1]
  end
end

# Test task copy version of trace
t = Trace(f2)

consume(t); consume(t)
a = fork(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4


text:
###############################################################################
#
# Optimize leaf invstment Vcmax+Jmax and Leaf area
#
###############################################################################
"""
    optimize_hs!(
                node::SPACSimple{FT},
                weather::Array{FT,2},
                Δt::FT = FT(1)
    ) where {FT<:AbstractFloat}

Optimize hydraulic conductance and leaf investment, given
- `node` [`SPACSimple`] type struct
- `weather` Weather profile in a growing season
- `Δt` Time period in `[h]`
"""
function optimize_hs!(
            node::SPACSimple{FT},
            weather::Array{FT,2},
            Δt::FT = FT(1)
) where {FT<:AbstractFloat}
    # 1. use the opt_laba and opt_vmax to initialize
    @inline f(x) = (
        tmp_node = deepcopy(node);
        (tmp_node).hs.root.k_max = (tmp_node).containerKS[1] * x[1];
        (tmp_node).hs.stem.k_max = (tmp_node).containerKS[2] * x[1];
        (tmp_node).hs.leaf.k_sla = (tmp_node).containerKS[3] * x[1];
        (tmp_node).hs.root.k_element .= (tmp_node).hs.root.k_max * (tmp_node).hs.root.N;
        (tmp_node).hs.stem.k_element .= (tmp_node).hs.stem.k_max * (tmp_node).hs.stem.N;
        (tmp_node).hs.leaf.k_element .= (tmp_node).hs.leaf.k_sla * (tmp_node).hs.leaf.N;
        leaf_allocation!(tmp_node, x[2], x[3]);
        tmp_prof = annual_profit(tmp_node, weather, Δt);
        return tmp_prof
    );

    ms = ReduceStepMethodND{FT}(
                x_mins=FT[0.001,1,5],
                x_maxs=FT[1e10, node.gaba*10, 200],
                x_inis=FT[node.hs.root.k_max/node.containerKS[1],
                          node.opt_laba,
                          node.opt_vmax],
                Δ_inis=FT[0.1,100,10]);
    st = SolutionToleranceND{FT}(FT[0.0009, 0.9, 0.09], 50);
    klv = find_peak(f, ms, st);

    # 2. update the opt_laba and opt_vmax to node
    node.hs.root.k_max = node.containerKS[1] * klv[1];
    node.hs.stem.k_max = node.containerKS[2] * klv[1];
    node.hs.leaf.k_sla = node.containerKS[3] * klv[1];
    node.hs.root.k_element .= node.hs.root.k_max * node.hs.root.N;
    node.hs.stem.k_element .= node.hs.stem.k_max * node.hs.stem.N;
    node.hs.leaf.k_element .= node.hs.leaf.k_sla * node.hs.leaf.N;
    leaf_allocation!(node, klv[2], klv[3]);
    node.opt_laba = klv[2];
    node.opt_vmax = klv[3];

    return nothing
end


text:
function waveAnim(path::ASCIIString, pathout="NULL"; inter=5, attribute="p", wbox=6, hbox=6, ox=0.0, dx=5.0, oz=0.0, dz=5.0, clip=1.0, cmap="seismic", aspect="auto", interval=40)
    fid = open(path, "r")
    isz = read(fid, Int64); isx   = read(fid, Int64);
    nz  = read(fid, Int64); nx    = read(fid, Int64);
    ext = read(fid, Int64); iflag = read(fid, Int64);
    ot  = read(fid, Float64); dt  = read(fid, Float64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx  = nx + 2*ext
    head    = sizeof(Int64) * 8
    cptSize = sizeof(Float64) * Nz * Nx
    sptSize = cptSize * 4
    nt =  round(Int64, (filesize(fid)-head)/sptSize)
    fig = plt.figure(figsize=(wbox, hbox))
    ims = PyCall.PyObject[]
    for it = 1 : inter : nt
        if attribute == "vz"
           position = head + (it-1)*sptSize; seek(fid, position);
           d = reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        elseif attribute == "vx"
           position = head + (it-1)*sptSize + cptSize; seek(fid, position);
           d = reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        elseif attribute == "p"
           position = head + (it-1)*sptSize + cptSize*2; seek(fid, position);
           d = reshape(read(fid, Float64, Nz*Nx), Nz, Nx) + reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        else
           error("wrong attribute")
        end
        vmax = quantile(vec(abs(d)), clip)
        vmin = -vmax
        im = plt.imshow(d, cmap=cmap, vmin=vmin,vmax=vmax,extent=[ox, ox+(size(d,2)-1)*dx, oz+(size(d,1)-1)*dz,oz], aspect=aspect)
        push!(ims, PyCall.PyObject[im])
    end
    close(fid)
    ani = anim.ArtistAnimation(fig, ims, interval=interval, blit=true)
    pathout = join([pathout ".mp4"])
    ani[:save](pathout)
    return nothing
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using Impero, Printf, MPI, LinearAlgebra, Statistics, GaussQuadrature
using ImperoHooks
using Makie
using Observables

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
const FT = Float64

N = 4
Ne_vert  = 1
Ne_horz  = 3
Rrange = range(0.5; length = Ne_vert + 1, stop = FT(1))

topl = StackedCubedSphereTopology(mpicomm, Ne_horz, Rrange)
grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
)

scene = visualize(grid)
## save interaction
fps = 10
record(scene, pwd() * "/sphere.mp4"; framerate = fps) do io
    for i = 1:200
        sleep(1/fps)
        recordframe!(io)
    end
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using NeidSolarScripts    # For continuum, solar rotaiton, differential rotation,...
 using RvSpectMLBase      # For spectra data structures
 using EchelleInstruments # For reading spectra
 using JLD2, FileIO       # For outputs
 using ArgParse           # For command line arguments

 # using CSV, DataFrames, Query, Dates, SHA
 #=
 println("# Parsing arguments...")
    function parse_commandline()
      s = ArgParseSettings( description = "Calculate continuum normalization from NEID L1 FITS files.")
      add_arg_group!(s, "Files", :argg_files)
      @add_arg_table! s begin
          "manifest"
              help = "Manifest file (CVS) containing FITS files to analyze.\nExpects columns named Filename, bjd, and target."
              arg_type = String
              default = "manifest.csv"
              #required = true
          "output"
              help = "Filename for output CCFs (jld2)"
              arg_type = String
              default = "daily_ccfs.jld2"
          "--overwrite"
             help = "Specify it's ok to overwrite the output file."
             action = :store_true
       end
       add_arg_group!(s, "Normalization parameters", :argg_norm_param)
       @add_arg_table! s begin
         "--orders_to_use"
            help = "First and last order _index_ to compute CCFs for"
            nargs = 2
            arg_type = Int64
            #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
            #default = [ first(orders_to_use_default(NEID2D())), last(orders_to_use_default(NEID2D())) ]
       end
       add_arg_group!(s, "Filter manifest for ", :argg_filter_param)
       @add_arg_table! s begin
          "--target"
             help = "Target field"
             arg_type = String
             default = "Sun"
          "--datestr"
             help = "Filenames that contain date string."
             arg_type = String
         "--max_solar_hour_angle"
            help = "Maximum absolute value of solar hour angle."
            arg_type = Float64
            default = 6.0
          "--max_airmass"
             help = "Maximum airmass."
             arg_type = Float64
             #default = 10.0
          "--min_snr_factor"
             help = "Minimum SNR relative to max_snr."
             arg_type = Float64
             #default = 0.0
          "--max_snr"
             help = "Specify max_snr manually."
             arg_type = Float64
          "--start_time"
             help = "Specify daily start time for CCFs to be processed (HH MM)"
             nargs = 2
             arg_type = Int64
             #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
             default = [0, 0]
         "--stop_time"
            help = "Specify daily stop time for CCFs to be processed (HH MM)"
            nargs = 2
            arg_type = Int64
            #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
            default = [23, 59]
          "--max_num_spectra"
             help = "Maximum number of spectra to process."
             arg_type = Int
             default = 250
       end

      return parse_args(s)
  end
  args = parse_commandline()
  =#

sed_path = "/home/eford/Code/RvSpectMLEcoSystem/NeidSolarScripts/data"
   sed = Continuum.read_master_sed_neid(path=sed_path)

fits_path = "/mnt/data_simons/NEID/DRPv0.7-fixedflatfielding2/"
   fits_fn = "neidL1_20210104T182500.fits"
   continuum_afs_path = joinpath(fits_path,"output","continuum")

   continuum_afs_fn = "neidL1_20210104T182500_continuum=afs.jld2"
   continuum_ras_path = joinpath(fits_path,"output","continuum")
   continuum_ras_fn = "neidL1_20210104T182500_continuum=ras1.jld2"
   write_output = false
   spec = NEID.read_data(joinpath(fits_path,fits_fn))

min_order_for_continuum = min_order(NEID2D())
  max_order_for_continuum = max_order(NEID2D())
  orders_to_use = min_order_for_continuum:max_order_for_continuum
  (f_norm, var_norm) = Continuum.normalize_by_sed(spec.λ,spec.flux,spec.var, sed; orders_to_use=orders_to_use)

using Plots
order_idx = 110
 plt = scatter(spec.λ[:,order_idx],f_norm[:,order_idx],ms=1.5,legend=:none)
 #plot!(plt, lambda,poly, color=:blue, lw=2)
 #plt = scatter(lambda,f_obs./(sed_view.*coeff.(lambda)),ms=1.5,legend=:none)
 #plt
#savefig("fig.png")
 ylims!(0.,2.0)
true

 #f_norm = spec.flux ./sed
 #xlabel!("λ")
 #ylabel!("Flux/SED")
 #savefig("fig.png")

get_pixel_range_for_continuum = GetPixelRangeFunctor(NEID2D())

#f_norm = spec.flux
@time (anchors, continuum_output, f_filtered) = Continuum.calc_continuum(spec.λ,f_norm,var_norm,stretch_factor=5.0, ν=1.0,
               merging_threshold=0.25, orders_to_use=orders_to_use, get_pixel_range=get_pixel_range_for_continuum,
               verbose=false)

order_idx = 90
  lambda_obs = view(spec.λ,:,order_idx)
  f_obs = view(f_norm,:,order_idx)
  var_obs = view(var_norm,:,order_idx)
  (anchors, continuum, f_filtered_internal) = Continuum.calc_continuum(lambda_obs,f_obs,var_obs,
       stretch_factor=5.0, ν=1.0, merging_threshold=0.25,  min_R_factor = 100.0, smoothing_half_width = 40, verbose=false)
  f_filtered = f_filtered_internal
  scatter(lambda_obs,f_obs, color=:black, ms=2, legend=:none)
  plot!(lambda_obs,f_filtered,color=:black)
  plot!(lambda_obs,continuum,color=:cyan)
  scatter!(lambda_obs[anchors], f_filtered[anchors], color=:blue, ms=4)

if write_output
  jldopen(continuum_ras_fn, "w") do file
    file["order_idx_with_continuum"] = orders_to_use
    file["continuum_anchors_idx"] = anchors
    file["continuum"] = continuum_output
    file["flux_filtered"] = f_filtered
  end
end

using Plots
order_idx = 50
 println("# Order Index: ", order_idx)
 lambda_obs = view(spec.λ,:,order_idx)
 f_obs = view(f_norm,:,order_idx)
 var_obs = view(var_norm,:,order_idx)
 (anchors, continuum, f_filtered_internal) = Continuum.calc_continuum(lambda_obs,f_obs,var_obs,
      stretch_factor=10.0, ν=1.0, merging_threshold=0.25, smoothing_half_width = 40, min_R_factor = 200.0, verbose=false)
  f_filtered = f_filtered_internal
  println(" size(anchors) = ", size(anchors), " size(continuum) = ", size(continuum), " size(f_filtered) = ", size(f_filtered) )

  #plt = scatter(lambda_obs,f_obs, color=:black, ms=2, legend=:none)
  plot!(plt,lambda_obs,f_filtered,color=:black)
  plot!(plt,lambda_obs,continuum,color=:cyan)
  #=
  if length(anchors) >=1
    scatter!(plt,lambda_obs[anchors], f_filtered[anchors], color=:blue, ms=4)
  end
  plt
  =#
  #=
  (anchors2, continuum2, f_filtered_internal2) = Continuum.calc_continuum(lambda,f_norm[:,order_idx],
       stretch_factor=5.0, ν=1.0, merging_threshold=0.25, smoothing_half_width = 40, verbose=false)
  println("# num_anchors1 = ", length(anchors), "  num_anchors2 = ", length(anchors2))
  plot!(lambda,continuum2,color=:magenta)
  scatter!(lambda[anchors2], f_filtered_internal2[anchors2], color=:red, ms=4)
  =#

xmin = first(lambda)-0.05; xmax = xmin+20;
 xlims!(xmin,xmax)
 idx_pix = searchsortedfirst(lambda,xmin):searchsortedlast(lambda,xmax)
 ylims!(NaNMath.extrema(f_obs[idx_pix]))

xmax = last(lambda)+0.05; xmin = xmax-20;
  xlims!(xmin,xmax)
  idx_pix = searchsortedfirst(lambda,xmin):searchsortedlast(lambda,xmax)
  ylims!(NaNMath.extrema(f_obs[idx_pix]))

  #xmin = 4570; xmax = 4580;
  xmin = 4945; xmax = 4965;
    xlims!(xmin,xmax)
    idx_pix = searchsortedfirst(lambda,xmin):searchsortedlast(lambda,xmax)
    ylims!(NaNMath.extrema(f_obs[idx_pix]))

xlims!(10840,10860)
xlims!(5880,5915)
  #xlims!(xmin,xmax)
xlims!(5350,5385)
xlims!(5385,5400)

xlims!(4950,4960)
true

  xlims!(5875,5905)
  #ylims!(0.95,1.15)
  #xmin = minimum(lambda_plt)
  #xmax = minimum(lambda_plt)
  #xlims!(xmin,xmax)
  #ylims!(extrema(f_filtered[xmin.<=lambda.<=xmax]))

  lambda = spec.λ[:,order_idx]
  f_obs = spec.flux[:,order_idx]
  findall(isnan.(f_obs))

 (anchors, continuum, f_filtered_internal) = Continuum.calc_continuum(lambda,f_obs,
       stretch_factor=5.0, ν=1.0, merging_threshold=0.25, verbose=true)

@time continuum = Continuum.calc_continuum(lambda,f_obs,stretch_factor=5.0, ν=1.0, merging_threshold=0.25, verbose=false)

#using Plots
continuum = zeros(size(spec.flux))
  for order_idx in 4:118
    #order_idx =30
    #possible_pix = get_pixel_range(get_inst(spec),order_idx)
    possible_pix = min_pixel_in_order(get_inst(spec)):max_pixel_in_order(get_inst(spec))
    bad_pix = bad_col_ranges(get_inst(spec),order_idx)
    pix_rng = EchelleInstruments.calc_complement_index_ranges(possible_pix,bad_pix)
    pix = mapreduce(p->collect(p),vcat,pix_rng)
    lambda = spec.λ[pix,order_idx]
    f_obs = spec.flux[pix,order_idx]
    (anchors, continuum_order, f_filtered) = Continuum.calc_continuum(lambda,f_obs,λout=spec.λ[possible_pix,order_idx],
         stretch_factor=10.0, merging_threshold=1.0, verbose=true)
    continuum[possible_pix,order_idx] .= continuum_order
    lambda_plt = spec.λ[possible_pix,order_idx]
    plot(lambda_plt,continuum_order)
    scatter!(lambda[anchors], f_filtered[anchors], color=:blue, ms=4)
  end

plot(spec.λ[:,60],continuum[:,60])
true

plot(lambda, spec.flux[:,order_idx]./continuum, legend=:none)
 xlabel!("λ")
 ylabel!("Flux (continuum-normalized)")
 ylims!(0.9,1.05)
 #ylims!(extrema( f_obs./continuum))
 #xlims!(minimum(lambda),5340)
true
=#

#=
  #plot!(lambda,f_clean,ms=1.5, color=:green)
  #plot!(lambda,f_clip_threshold,color=:blue)
  #ylims!(0.9,1.2)
  #plot!(lambda,f_median_filtered.-f_clean)

 #plot!(lambda,rollingpin_radius./rollingpin_R_max)
 #plot(lambda,f_obs, legend=:none)
 #plot(lambda,f_obs, legend=:none)
 #plot!(lambda,f_interp)
 anchor_vals = f_median_filtered[anch]
 #scatter!(lambda[anch],anchor_vals, color=:blue)
 scatter!(lambda[anch],anchor_vals, color=:red)
 scatter!(lambda[anch[anch_mask]],anchor_vals[anch_mask], color=:blue)
 #interp_cubic = CubicSpline(f_median_filtered[anch],anchor_vals )
 #plot!(lambda,interp_cubic.(lambda), color=:green)
 ylims!(0.9,1.2)
 #xlims!(5006.5,5008)

  scatter!(lambda[anchors_merged],f_median_filtered[anchors_merged], color=:magenta)
 plot!(lambda,calc_continuum_from_anchors(lambda,f_median_filtered,anchors_merged),color=:red)

 #=
 plot!(λ_med_filtered,f_median_filtered)
 plot!(λ_med_filtered,f_median_filtered.+iqr)
 plot!(λ_med_filtered,f_median_filtered.-iqr)
 plot!(λ_med_filtered,f_interp, ms=2.0)
 scatter!(λ_med_filtered[anch],f_median_filtered[anch], ms=4.0)
 =#
 #ylims!(0.9,1.03)

#interp_bsp = BSplineApprox(lambda, f_obs,3,4,:ArcLen,:Average)
=#


text:
module search_utils

function save_settings(save_loc, settings_vec)
    if split(basename(save_loc), ".")[end] == "txt"
        open(save_loc, "w") do f
            for v in settings_vec
                write(f, string(v)*"\n")
            end
        end
    else
        @save save_loc Dict("settings"=>settings_vec)
    end
end

"""
    best_settings

This function takes an experiment directory and finds the best setting for the product of arguments
with keys specified by product_args. To see a list of viable arguments use
    `ic = ItemCollection(exp_loc); diff(ic.items)`

If a save_loc is provided, this will save to the file specified. The fmt must be supported by FileIO and be able to take dicts.

Additional kwargs are passed to order_settings.

"""

function best_settings(exp_loc, product_args::Vector{String};
                       save_loc="", kwargs...)

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)

    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    settings_dict = Dict()
    for (arg_idx, arg) in enumerate(args)
        search_dict = Dict([product_args[idx]=>key for (idx, key) in enumerate(arg)]...)
        ret = order_settings(exp_loc; set_args=search_dict, ic=ic, kwargs...)
        settings_dict[search_dict] = ret[1]
    end

    if save_loc != ""
        save(save_loc, Dict("best_settings"=>settings_dict))
    else
        return settings_dict
    end

end

"""
    order_settings

    This provides a mechanism to order the settings of an experiment.

    kwargs:
        `set_args(=Dict{String, Any}())`: narrowing the search parameters. See best_settings for an example of use.

        `clean_func(=identity)`: The function used to clean the loaded data
        `runs_func(=mean)`: The function which takes a vector of floats and produces statistics. Must return either a Float64 or Dict{String, Float64}. (WIP, any container/primitive which implements get_index).

        `lt(=<)`: The less than comparator.
        `sort_idx(=1)`: The idx of the returned `runs_func` structure used for sorting.
        `run_key(=run)`: The key used to specify an ind run for an experiment.

        `results_file(=\"results.jld2\")`: The string of the file containing experimental results.
        `save_loc(=\"\")`: The save location (returns settings_vec if not provided).
        `ic(=ItemCollection([])`: Optional item_collection, not needed in normal use.
"""

function order_settings(exp_loc;
                        results_file="results.jld2",
                        clean_func=identity, runs_func=mean,
                        lt=<, sort_idx=1, run_key="run",
                        set_args=Dict{String, Any}(),
                        ic=ItemCollection([]), save_loc="")

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end

    exp_path = dirname(exp_loc)
    if length(ic.items) == 0
        ic = ItemCollection(exp_loc)
    end
    diff_dict = diff(ic.items)
    product_args = collect(filter((k)->(k!=run_key && k∉keys(set_args)), keys(diff_dict)))

    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    settings_vec =
        Vector{Tuple{Union{Float64, Vector{Float64}, Dict{String, Float64}}, Dict{String, Any}}}(undef, length(args))

    #####
    # Populate settings Vector
    #####
    @showprogress 0.1 "Setting: " for (arg_idx, arg) in enumerate(args)

        search_dict = merge(
            Dict([product_args[idx]=>key for (idx, key) in enumerate(arg)]...),
            set_args)
        _, hashes, _ = search(ic, search_dict)
        μ_runs = zeros(length(hashes))
        for (idx_d, d) in enumerate(hashes)
            if isfile(joinpath(exp_path, d, results_file))
                results = load(joinpath(exp_path, d, results_file))
                μ_runs[idx_d] = clean_func(results)
            else
                μ_runs[idx_d] = Inf
            end
        end
        settings_vec[arg_idx] = (runs_func(μ_runs), search_dict)
    end

    #####
    # Sort settings vector
    #####
    sort!(settings_vec; lt=lt, by=(tup)->tup[1][sort_idx])

    #####
    # Save
    #####
    if save_loc != ""
        save(save_loc, Dict("settings"=>settings_vec))
    else
        return settings_vec
    end
end

function collect_data(exp_loc;
                      run_arg="run",
                      results_file="results.jld2", settings_file="settings.jld2", clean_func=identity,
                      save_dir="collected")

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end
    head_dir = dirname(exp_loc)

    if !isdir(joinpath(exp_loc, save_dir))
        mkdir(joinpath(exp_loc, save_dir))
    end

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)
    # args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    search_dict = Dict(run_arg=>diff_dict[run_arg][1])

    _, hashes, _ = search(ic, search_dict)

    settings_vec = Vector{Dict}(undef, length(hashes))

    # collect the parameter settings run
    for (idx, h) in enumerate(hashes)
        sett = load(joinpath(head_dir, h, settings_file))["parsed_args"]
        settings_vec[idx] = Dict(k=>sett[k] for k in filter(v -> v != run_arg, keys(diff_dict)))
    end

    @showprogress for (idx, stngs) in enumerate(settings_vec)
        # println(length(search(ic, stngs)[2]))
        hashes = search(ic, stngs)[2]

        v = Vector{Any}(undef, length(hashes))
        for (idx, h) in enumerate(hashes)
            v[idx] = clean_func(load(joinpath(head_dir, h, results_file)))
        end
        save(joinpath(exp_loc, save_dir, join(["$(k)_$(stngs[k])" for k in keys(stngs)], '_')*".jld2"), Dict("results"=>v, "settings"=>stngs))
    end

end

function collect_sens_data(exp_loc, sens_param, product_args;
                           run_arg="run",
                           results_file="results.jld2", settings_file="settings.jld2", clean_func=identity,
                           save_dir="collected_sens", ignore_nans=false, ignore_sens=nothing)

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end
    head_dir = dirname(exp_loc)

    if !isdir(joinpath(exp_loc, save_dir))
        mkdir(joinpath(exp_loc, save_dir))
    end

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)
    if ignore_sens != nothing
        diff_dict[sens_param] = filter(x->x!=ignore_sens, diff_dict[sens_param])
        println(diff_dict[sens_param])
    end
    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    println(collect(args))

    for arg in collect(args)

        println([k=>arg[k_idx] for (k_idx, k) in enumerate(product_args)])

        search_dict = Dict(run_arg=>diff_dict[run_arg][1], [k=>arg[k_idx] for (k_idx, k) in enumerate(product_args)]...)

        _, hashes, _ = search(ic, search_dict)

        settings_vec = Vector{Dict}(undef, length(hashes))

        # collect the parameter settings run
        for (idx, h) in enumerate(hashes)
            sett = load(joinpath(head_dir, h, settings_file))["parsed_args"]
            settings_vec[idx] = Dict([k=>sett[k] for k in filter(v -> v ∉ keys(search_dict), keys(diff_dict))]..., [k=>arg[k_idx] for (k_idx, k) in enumerate(product_args)]...)
        end

        avg_res = zeros(length(diff_dict[sens_param]))
        std_err = zeros(length(diff_dict[sens_param]))

        for (idx, stngs) in enumerate(settings_vec)
            # println(length(search(ic, stngs)[2]))
            hashes = search(ic, stngs)[2]

            v = zeros(length(hashes))
            for (idx, h) in enumerate(hashes)
                v[idx] = clean_func(load(joinpath(head_dir, h, results_file)))
            end

            sens_idx = findfirst(x->x==stngs[sens_param], diff_dict[sens_param])
            if sens_idx != nothing
                filtered = filter(x->!isnan(x), v)
                println(stngs, ": ", filtered)
                avg_res[sens_idx] = mean(filtered)
                std_err[sens_idx] = std(filter(x->!isnan(x), v))/sqrt(length(filter(x->!isnan(x), v)))
            end

        end

        save(joinpath(exp_loc, save_dir, "collect_"*join(["$(k)_$(arg[k_idx])" for (k_idx, k) in enumerate(product_args)], '_')*".jld2"), Dict("avg"=>avg_res, "std_err"=>std_err, "sens"=>diff_dict[sens_param], "settings"=>Dict([k=>arg[k_idx] for (k_idx, k) in enumerate(product_args)])))
    end
end

end # module search_utils


text:
using JET, Lux, Random, Test

include("../utils.jl")

rng = Random.default_rng()
Random.seed!(rng, 0)

@testset "RNNCell" begin for rnncell in (RNNCell(3 => 5, identity),
                                         RNNCell(3 => 5, tanh),
                                         RNNCell(3 => 5, tanh; bias=false),
                                         RNNCell(3 => 5, identity; bias=false))
    println(rnncell)
    ps, st = Lux.setup(rng, rnncell)
    x = randn(rng, Float32, 3, 2)
    h, st_ = Lux.apply(rnncell, x, ps, st)

    run_JET_tests(rnncell, x, ps, st)
    run_JET_tests(rnncell, (x, h), ps, st_)

    function loss_loop_rnncell(p)
        h, st_ = rnncell(x, p, st)
        for i in 1:10
            h, st_ = rnncell((x, h), p, st_)
        end
        return sum(abs2, h)
    end

    test_gradient_correctness_fdm(loss_loop_rnncell, ps, atol=1e-3, rtol=1e-3)
end end

@testset "LSTMCell" begin
    lstmcell = LSTMCell(3 => 5)
    println(lstmcell)
    ps, st = Lux.setup(rng, lstmcell)
    x = randn(rng, Float32, 3, 2)
    (h, c), st_ = Lux.apply(lstmcell, x, ps, st)

    run_JET_tests(lstmcell, x, ps, st)
    run_JET_tests(lstmcell, (x, h, c), ps, st_)

    function loss_loop_lstmcell(p)
        (h, c), st_ = lstmcell(x, p, st)
        for i in 1:10
            (h, c), st_ = lstmcell((x, h, c), p, st_)
        end
        return sum(abs2, h)
    end

    test_gradient_correctness_fdm(loss_loop_lstmcell, ps, atol=1e-3, rtol=1e-3)
end

@testset "GRUCell" begin
    grucell = GRUCell(3 => 5)
    println(grucell)
    ps, st = Lux.setup(rng, grucell)
    x = randn(rng, Float32, 3, 2)
    h, st_ = Lux.apply(grucell, x, ps, st)

    run_JET_tests(grucell, x, ps, st)
    run_JET_tests(grucell, (x, h), ps, st_)

    function loss_loop_grucell(p)
        h, st_ = grucell(x, p, st)
        for i in 1:10
            h, st_ = grucell((x, h), p, st_)
        end
        return sum(abs2, h)
    end

    test_gradient_correctness_fdm(loss_loop_grucell, ps, atol=1e-3, rtol=1e-3)
end

@testset "multigate" begin
    x = rand(6, 5)
    res, (dx,) = Zygote.withgradient(x) do x
        x1, _, x3 = Lux.multigate(x, Val(3))
        sum(x1) + sum(x3 .* 2)
    end
    @test res == sum(x[1:2, :]) + 2sum(x[5:6, :])
    @test dx == [ones(2, 5); zeros(2, 5); fill(2, 2, 5)]
end


text:
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables
using CategoricalArrays
using Tables
using StatsModels
cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS2-optimization-intro")

function fullquestions()
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# question 1
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
	minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
	startval = rand(1)   # random starting value
	result = optimize(minusf, startval, BFGS())

	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# question 2
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
	df = CSV.read(HTTP.get(url).body)
	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = df.married.==1

	function ols(beta, X, y)
		ssr = (y.-X*beta)'*(y.-X*beta)
		return ssr
	end

	beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
	println(vec(beta_hat_ols.minimizer))

	bols = inv(X'*X)*X'*y
	df.white = df.race.==1
	bols_lm = lm(@formula(married ~ age + white + collgrad), df)

	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# question 3
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	function logit(alpha, X, d)
		# your turn
		loglike = sum(X*alpha.*d .- log.(1 .+exp.(X*alpha)))
		return -loglike
	end

	beta_hat_logit = optimize(b -> logit(b, X, y),beta_hat_ols.minimizer , LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
	println(beta_hat_logit.minimizer)

	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# question 4
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# check it
	df.raceind = df.race .==1
	beta_glm = glm(@formula(married ~ age + raceind+collgrad), df, Binomial(), LogitLink())
	if !all(round.(coef(beta_glm) .- beta_hat_logit.minimizer,digits=5).==0)
		error()
	end

	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	# question 5
	#:::::::::::::::::::::::::::::::::::::::::::::::::::
	freqtable(df, :occupation) # note small number of obs in some occupations
	df = dropmissing(df, :occupation)
	df[df.occupation.==8 ,:occupation] .= 7
	df[df.occupation.==9 ,:occupation] .= 7
	df[df.occupation.==10,:occupation] .= 7
	df[df.occupation.==11,:occupation] .= 7
	df[df.occupation.==12,:occupation] .= 7
	df[df.occupation.==13,:occupation] .= 7
	freqtable(df, :occupation) # problem solved

	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = df.occupation

	# set base occ to be first one.
	function mlogit(alpha, X, d)
		u = reduce(hcat,[X*alpha[(4*i-3):(4*i)] for i in 1:6])
		help = reduce(hcat,[d.==i for i in 2:7])
		loglike = sum(help .*u) - sum(log.(1 .+sum(exp.(u), dims=2)))
		return -loglike
	end

	# starting values - set intercepts to fraction present in data.
	starting = Vector(prop(freqtable(df, :occupation)))
	starting = repeat(starting[Not(1)], inner = 4)
	starting[Not((1:6) .*4 .-3)] .= 0

	beta_hat_mlogit = optimize(b -> mlogit(b, X, y),starting, BFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
	println(vec(beta_hat_mlogit.minimizer))
end

fullquestions()


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
omib_file_dir = joinpath(dirname(@__FILE__), "OMIB.raw")
omib_sys =
    System(PowerModelsData(omib_file_dir), runchecks = false, unit_system = "device_base")
add_source_to_ref(omib_sys)
res = solve_powerflow!(omib_sys)
############### Data Dynamic devices ########################
function dyn_gen_classic(generator)
    return PSY.DynamicGenerator(
        generator,
        1.0, #ω_ref
        machine_classic(), #machine
        shaft_damping(), #shaft
        avr_none(), #avr
        tg_none(), #tg
        pss_none(),
    ) #pss
end

#Attach dynamic generator. Currently use PSS/e format based on bus #.
gen = [g for g in get_components(Generator, omib_sys)][1]
case_gen = dyn_gen_classic(gen)
add_component!(omib_sys, case_gen)

#Compute Y_bus after fault
fault_branch = deepcopy(collect(get_components(Branch, omib_sys))[1])
fault_branch.r = 0.02;
fault_branch.x = 0.1;
Ybus_fault = PSY.Ybus([fault_branch], get_components(Bus, omib_sys))[:, :]


text:
using Revise

using JDF

using CSV, DataFrames, Blosc, JLSO, Base.GC

# use 12 threads
Blosc.set_num_threads(6)
# GC.gc()
@time a = CSV.read(
    "C:/data/Performance_All/Performance_2016Q1.txt",
    delim = '|',
    header = false
);

GC.gc()

t = time()
@time metadatas = psavejdf(a, "c:/data/large.jdf");
time() - t
GC.gc()

GC.gc()
@time metadatas = savejdf(a, "c:/data/large.jdf");
GC.gc()

@time JLSO.save("c:/data/large.meta.jlso", metadatas)
GC.gc()

using Revise, JLSO, JDF
@time metadatas = JLSO.load("c:/data/large.meta.jlso")["data"];
GC.gc()

@time a2 = loadjdf("c:/data/large.jdf", metadatas);
GC.gc()

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))

# CSVFiles vs CSV
using Revise
using CSV, DataFrames, JDF, Base.GC
import Base.Threads: @spawn

# a = gf()
fp = [f for f in readdir("c:/data/AirOnTimeCSV/") if endswith(f, ".csv")];

using CSVFiles

@time a= load(joinpath("c:/data/AirOnTimeCSV", fp[1]), type_detect_rows = 2000)

#psavejdf(a, joinpath("c:/data/hehe/", fp[1]))
sort!(fp, by = x->filesize(joinpath("c:/data/AirOnTimeCSV/", x)), rev=true)

fn4(fp) = begin
    #@time a = load(joinpath("C:/data/AirOnTimeCSV", fp), type_detect_rows = 2000) |> DataFrame
    @time a = CSV.read(joinpath("C:/data/AirOnTimeCSV", fp))
    #mkdir(joinpath("c:/data/hehe/", fp))
    res = psavejdf(a, joinpath("c:/data/hehe/", fp))
    a = nothing
    res
end

res = Vector{Any}(undef, length(fp))

@time for (i, f) in enumerate(fp[1:6])
    print(i)
    res[i] = fn4(f)
end

#@time fn4(fp[1])
#@time fn4.(fp)
2+2


text:
#=
JoiceEstimationFunctions:
- Julia version: 1.6
- Author: connorforsythe
- Date: 2021-11-07
=#

using JuMP, NLopt, LinearAlgebra, Statistics, PyCall, ProgressMeter, Distributions, FiniteDiff

export EstimateModel, testFun

function estimateHessianUsingObjective(modelInfo::ModelInfo, theta::T...) where {T<:Real}

    logitVals::Union{Vector{T}, Matrix{T}} = logit(modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta...)

    if typeof(logitVals)==Matrix{T}

        parameterDraws::Union{Matrix{T}, Nothing} = getSimulatedParameters(modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, theta...)
        normalDraws::Union{Matrix{T}, Nothing} = getStandardNormalDraws(modelInfo.RandomDraws)
        df::Union{DataFrame, Nothing} = DataFrame(logitVals, :auto)

    else
        parameterDraws = nothing
        normalDraws = nothing
        df = nothing
    end

    obj(x::Vector{T}) = getObjLM(logitVals, df, parameterDraws, normalDraws, modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ChoiceColumn, modelInfo.ModelSpace, modelInfo.DataStructure, x...)

    hess = FiniteDiff.finite_difference_hessian(obj, [theta...])

    return hess
end

function getObjLM(logitVals::AbstractMatrix{T}, df::DataFrame, parameterDraws::AbstractMatrix{T}, normalDraws::AbstractMatrix{T}, data::DataFrame, parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, randomDraws::Matrix{T}, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, choiceColumn::Symbol, modelSpace::Preference, dataStructure::Panel, theta::T...) where {T<:Real}
    logit!(logitVals, parameterDraws, normalDraws, data, parameters, randomParameters, randomDraws, choiceSetIDColumn, observationIDColumn, modelSpace, dataStructure, theta...)
    return logLikelihoodLM(logitVals, data[:, choiceColumn], data[:, observationIDColumn], df)
end

function getObjLM(logitVals::AbstractVector{T}, df::Nothing, parameterDraws::Nothing, normalDraws::Nothing, data::DataFrame, parameters::Vector{Symbol}, randomParameters::Nothing, randomDraws::Nothing, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, choiceColumn::Symbol, modelSpace::Preference, dataStructure::Panel, theta::T...) where {T<:Real}
    logit!(logitVals, parameterDraws, normalDraws, data, parameters, randomParameters, randomDraws, choiceSetIDColumn, observationIDColumn, modelSpace, dataStructure, theta...)
    return logLikelihoodLM(logitVals, data[:, choiceColumn], data[:, observationIDColumn], df)
end
function EstimateModel(modelInfo::ModelInfo; numMultistarts::Int = 1, post::Bool = true)

    startingVals::Vector{T} where {T<:Real} = getStartingValues(modelInfo)
    res = EstimateModel(modelInfo, startingVals...; numMultistarts=numMultistarts, post=post)

    return res

end

function EstimateModel(modelInfo::ModelInfo, theta0::T...; numMultistarts::Int = 1, post::Bool = true) where {T<:Real}
    println("Beginning estimation procedure...")
    res = buildAndEstimateJuMPModel(modelInfo, theta0...; numMultistarts=numMultistarts, post=post)
    println("Estimation procedure completed!")
    return res

end

function EstimateModelLM(modelInfo::ModelInfo; numMultistarts::Int = 1, post::Bool = true)

    startingVals::Vector{T} where {T<:Real} = getStartingValues(modelInfo)
    println("Starting values: $startingVals")
    res = EstimateModelLM(modelInfo, startingVals...; numMultistarts=numMultistarts, post=post)

    return res

end

function EstimateModelLM(modelInfo::ModelInfo, theta0::T...; numMultistarts::Int = 1, post::Bool = true) where {T<:Real}
    @time logitVals::Union{Vector{T}, Matrix{T}} = logit(modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta0...)

    logitVals = logitVals.*0

    if typeof(logitVals)==Matrix{T}
        println("This is true")
        logitVals = Matrix{T}(undef, ncol(modelInfo.Data), modelInfo.InterDrawCount)
        parameterDraws::Matrix{T} = getSimulatedParameters(modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, theta0...)
        normalDraws::Matrix{T} = getStandardNormalDraws(modelInfo.RandomDraws)

        @time logit!(logitVals, parameterDraws, normalDraws, modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta0...)
    else
        @time logit!(logitVals, modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta0...)
    end

    @time logit!(logitVals, modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta0...)
    println("-------")

end

function createMemoObjGrad(objGradFun::Function)
    prevX, prevObjGrad = nothing, nothing
    function memoGetObjGrad(i, x::T...) where{T<:Real}
        if(x!=prevX)
            prevX, prevObjGrad = x, objGradFun(x...)
        end
        return prevObjGrad[i]
    end

    return [(x...)-> memoGetObjGrad(i, x...) for i in 1:2]

end

function buildAndEstimateJuMPModel(modelInfo::ModelInfo, theta0::T...; numMultistarts::Int = 1, post::Bool=true) where{T<:Real}
    println("Constructing JuMP Model...")
    objGradToMemo(x...) = getObjGrad(modelInfo, x...)

    memoObjGrad = createMemoObjGrad(objGradToMemo)

    numMultistarts = 1

    numParams::Int = length(theta0)

    obj(x...) = memoObjGrad[1](x...)
    grad(x...) = memoObjGrad[2](x...)

    scaledGrad(x...) = memoObjGrad[2](x...)./10
    jumpGrad(g::AbstractVector{T}, x::T...) where{T<:Real} = placeGradInVector(g, scaledGrad, x...)

    rModel = Model(NLopt.Optimizer)
    set_optimizer_attribute(rModel, "algorithm", :LD_LBFGS)
    set_optimizer_attribute(rModel, "xtol_rel", 1e-4)
    set_optimizer_attribute(rModel, "ftol_rel", 1e-4)
    set_optimizer_attribute(rModel, "xtol_abs", 1e-4)
    set_optimizer_attribute(rModel, "ftol_abs", 1e-4)
    set_optimizer_attribute(rModel, "maxeval", 1000)
    register(rModel, :LL, numParams, obj, jumpGrad; autodiff=false)
    @variable(rModel, theta[1:numParams])

    for i in 1:length(theta0)

        JuMP.set_start_value(theta[i], theta0[i])

    end

    @NLobjective(rModel, Max, LL(theta...))
    println("JuMP Model Constructed...")
    println("JuMP Model Initial Estimation...")
    maxLL = -Inf
    bestParams = [theta0...]
    try
        JuMP.optimize!(rModel)

        maxLL = objective_value(rModel)
        bestParams = value.(theta)
    catch e
    end
    println("JuMP Model Initial Estimation Complete...")

    if(numMultistarts>1)
        randomDraws = getDraws(numParams, numMultistarts-1)
        parameterDraws = min(modelInfo.PreferenceSearchBounds...).+randomDraws.*(max(modelInfo.PreferenceSearchBounds...)-min(modelInfo.PreferenceSearchBounds...))

        pMeter = Progress(numMultistarts-1, desc="Running Multistarts:")
        for j in 1:(numMultistarts-1), i in 1:numParams
            JuMP.set_start_value(theta[i], parameterDraws[i, j])
            if(i==numParams)
                try
                    JuMP.optimize!(rModel)

                    maxLL = objective_value(rModel)
                    bestParams = value.(theta)

                    if(objective_value(rModel)>maxLL)
                        maxLL = objective_value(rModel)
                        bestParams = value.(theta)
                    end
                catch e
                end
                next!(pMeter)
            end
        end

    end

    println("Refining Model Estimates...")

    set_optimizer_attribute(rModel, "xtol_rel", 1e-15)
    set_optimizer_attribute(rModel, "ftol_rel", 1e-15)
    set_optimizer_attribute(rModel, "xtol_abs", 1e-15)
    set_optimizer_attribute(rModel, "ftol_abs", 1e-15)
    set_optimizer_attribute(rModel, "maxeval", 10000)

    for i in 1:length(bestParams)
        JuMP.set_start_value(theta[i], bestParams[i])
    end

    try
        JuMP.optimize!(rModel)

        maxLL = objective_value(rModel)
        bestParams = value.(theta)
    catch e
        println("Issue with refining estimates - try a different starting point.")
        println("Returning nothing.")
        return nothing
    end
    println("Model Estimates Refined...")
    println("Post-Estimation Calculations Beginning...")
    r = Dict{Any, Any}()
    coefNames::Vector{Any} = getParameterNames(modelInfo.Parameters, modelInfo.RandomParameters)
    r["LL"] = maxLL
    if post
        r["coef"] = NamedTuple{Tuple(Symbol.(coefNames))}(Tuple(bestParams))
    else
        r["coef"] = [bestParams...]
    end

    r["grad"] = grad(bestParams...)

    if post
        obsCount::Int = length(unique(modelInfo.Data[:, modelInfo.ObservationIDColumn]))

        r["hessian"] = estimateHessianUsingObjective(modelInfo, bestParams...)

        r["cov"] = -inv(r["hessian"])
        r["se"] = NamedTuple{Tuple(Symbol.(coefNames))}(Tuple(sqrt.(diag(r["cov"]))))

        r["obsGrad"] = getObservationGradients(modelInfo, value.(theta)...)

        meanGradMat = repeat(mean(r["obsGrad"], dims=1), size(r["obsGrad"])[1])

        r["robustCov"] = r["cov"]*(transpose(r["obsGrad"].-meanGradMat)*(r["obsGrad"].-meanGradMat))*r["cov"]
        r["robustCov"] = (obsCount/(obsCount-1)).*r["robustCov"]
        r["robustSE"] = NamedTuple{Tuple(Symbol.(coefNames))}(Tuple(sqrt.(diag(r["robustCov"]))))
        println("Post-Estimation Calculations Completed...")
    end
    println("Model Estimated!")
    return r

end

function getParameterNames(parameters::Vector{Symbol}, randomParameters::Nothing) return string.(parameters) end
function getParameterNames(parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String})
    r = string.(deepcopy(parameters))
    for coef in keys(randomParameters)
        push!(r, string(coef, "-μ"))
    end

    for coef in keys(randomParameters)
        push!(r, string(coef, "-σ"))
    end

    return r
end

function getHessianNumericallyFromObjective(obj::Function, x::T...) where {T<:Real}
    passFun(x::Vector{T}) = obj(x...)
    hess = FiniteDiff.finite_difference_hessian(passFun, [x...])

    return hess

end

function placeGradInVector(g::AbstractVector{T}, fun::Function, x::T...) where{T}
    r = fun(x...)
    for i in 1:length(x)
        g[i] = r[i]
    end
end

function getObjGrad(modelInfo::ModelInfo, theta::T...) where {T<:Real}

    logitVals::Union{Vector{T}, Matrix{T}} = logit(modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta...)
    parameterCols = getParameterColumns(modelInfo.Parameters, modelInfo.RandomParameters)
    logLikelihoodVal::T = logLikelihood(logitVals, modelInfo.Data[:, modelInfo.ChoiceColumn], modelInfo.Data[:, modelInfo.ObservationIDColumn])
    logLikelihoodGradVals::Vector{T} = logLikelihoodGrad(logitVals, modelInfo.Data[:, modelInfo.ChoiceColumn], Matrix(modelInfo.Data[:, parameterCols]), modelInfo.Data[:, modelInfo.ObservationIDColumn], modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, theta...)

    return logLikelihoodVal, logLikelihoodGradVals
end

function getObjGradLM(logitVals::AbstractVecOrMat{T}, modelInfo::ModelInfo, theta::T...) where {T<:Real}

    logit!(logitVals, modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta...)
    parameterCols = getParameterColumns(modelInfo.Parameters, modelInfo.RandomParameters)
    logLikelihoodVal::T = logLikelihood(logitVals, modelInfo.Data[:, modelInfo.ChoiceColumn], modelInfo.Data[:, modelInfo.ObservationIDColumn])
    logLikelihoodGradVals::Vector{T} = logLikelihoodGrad(logitVals, modelInfo.Data[:, modelInfo.ChoiceColumn], Matrix(modelInfo.Data[:, parameterCols]), modelInfo.Data[:, modelInfo.ObservationIDColumn], modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, theta...)

    return logLikelihoodVal, logLikelihoodGradVals
end

function getParameterColumns(parameters::Vector{Symbol}, randomParameters::Nothing)

    return parameters

end

function getParameterColumns(parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String})
    r = deepcopy(parameters)

    push!(r, keys(randomParameters)...)

    return r

end

function getObservationGradients(modelInfo::ModelInfo, theta::T...) where {T<:Real}

    obsCol = modelInfo.ObservationIDColumn

    logitVals::Union{Vector{T}, Matrix{T}} = logit(modelInfo.Data, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, modelInfo.ChoiceSetIDColumn, modelInfo.ObservationIDColumn, modelInfo.ModelSpace, modelInfo.DataStructure, theta...)
    uniqueObsIDs = unique(modelInfo.Data[:, obsCol])

    obsGrads = Matrix{T}(undef, length(uniqueObsIDs), length(theta))

    i::Int = 1

    pMeter = Progress(length(uniqueObsIDs); desc="Computing observation gradients for covariance estimation: " )

    for tempID in uniqueObsIDs

        tempObsIndices = modelInfo.Data[:, obsCol].==repeat([tempID], nrow(modelInfo.Data))
        tempObs = modelInfo.Data[tempObsIndices, obsCol]
        parameterCols = getParameterColumns(modelInfo.Parameters, modelInfo.RandomParameters)
        tempXMat = Matrix(modelInfo.Data[tempObsIndices, parameterCols])
        tempChoices = modelInfo.Data[tempObsIndices, modelInfo.ChoiceColumn]
        if(modelInfo.RandomParameters!=nothing)
            tempLogit = logitVals[tempObsIndices, :]
        else
            tempLogit = logitVals[tempObsIndices]
        end

        tempGrad::Vector{T} = logLikelihoodGrad(tempLogit, tempChoices, tempXMat, tempObs, modelInfo.Parameters, modelInfo.RandomParameters, modelInfo.RandomDraws, theta...)

        obsGrads[i, 1:length(theta)] = tempGrad

        i+=1
        next!(pMeter)

    end

    return obsGrads

end

function logit(data::DataFrame, parameters::Vector{Symbol}, randomParameters::Nothing, randomDraws::Nothing, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, modelSpace::Preference, dataStructure::Union{Cross, Panel}, theta::T...) where {T<:Real}

    tempDF::DataFrame = DataFrame(expU = exp.(Matrix(data[:, parameters])*[theta...]))
    tempDF[:, choiceSetIDColumn]=data[:, choiceSetIDColumn]

    gData::GroupedDataFrame = groupby(tempDF, choiceSetIDColumn)
    gDataRes::DataFrame = combine(gData, :expU => sum)

    tempDF = innerjoin(tempDF, gDataRes, on=choiceSetIDColumn) #TODO fix so that obsID is forced to be unique

    tempDF[:, :p] .= tempDF[:, :expU]./tempDF[:, :expU_sum]
    r::Vector{T} = tempDF[:, :p]

    return r
end

function logit!(logitVals::AbstractMatrix{T}, index::Int, data::DataFrame, parameters::Vector{Symbol}, randomParameters::Nothing, randomDraws::Nothing, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, modelSpace::Preference, dataStructure::Union{Cross, Panel}, theta::T...) where {T<:Real}

    tempDF::DataFrame = DataFrame(expU = exp.(Matrix(data[:, parameters])*[theta...]))
    tempDF[:, choiceSetIDColumn]=data[:, choiceSetIDColumn]

    gData::GroupedDataFrame = groupby(tempDF, choiceSetIDColumn)
    gDataRes::DataFrame = combine(gData, :expU => sum)

    tempDF = innerjoin(tempDF, gDataRes, on=choiceSetIDColumn) #TODO fix so that obsID is forced to be unique

    tempDF[:, :p] .= tempDF[:, :expU]./tempDF[:, :expU_sum]
    logitVals[:, index] .= tempDF[:, :p]

end

function logit!(logitVals::AbstractVector{T}, parameterDraws::Nothing, normalDraws::Nothing, data::DataFrame, parameters::Vector{Symbol}, randomParameters::Nothing, randomDraws::Nothing, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, modelSpace::Preference, dataStructure::Union{Cross, Panel}, theta::T...) where {T<:Real}

    tempDF::DataFrame = DataFrame(expU = exp.(Matrix(data[:, parameters])*[theta...]))
    tempDF[:, choiceSetIDColumn]=data[:, choiceSetIDColumn]

    gData::GroupedDataFrame = groupby(tempDF, choiceSetIDColumn)
    gDataRes::DataFrame = combine(gData, :expU => sum)

    tempDF = innerjoin(tempDF, gDataRes, on=choiceSetIDColumn) #TODO fix so that obsID is forced to be unique

    tempDF[:, :p] .= tempDF[:, :expU]./tempDF[:, :expU_sum]
    logitVals .= tempDF[:, :p]

end

function logit(data::DataFrame, parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, randomDraws::Matrix{T}, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, modelSpace::Preference, dataStructure::Panel, theta::T...) where {T<:Real}

    parameterDraws::Matrix{T} = getSimulatedParameters(parameters, randomParameters, randomDraws, theta...)
    normalDraws::Matrix{T} = getStandardNormalDraws(randomDraws)

    paramsInOrder = [parameters..., keys(randomParameters)...]

    numDraws = size(normalDraws)[2]

    logitVals = Matrix{T}(undef, size(data)[1], numDraws)

    for i in 1:numDraws
        logitVals[:, i] = logit(data, paramsInOrder, nothing, nothing, choiceSetIDColumn, observationIDColumn, modelSpace, dataStructure, parameterDraws[:, i]...)
    end

    return logitVals

end

function logit!(logitVals::AbstractMatrix{T}, parameterDraws::AbstractMatrix{T}, normalDraws::AbstractMatrix{T}, data::DataFrame, parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, randomDraws::Matrix{T}, choiceSetIDColumn::Symbol, observationIDColumn::Symbol, modelSpace::Preference, dataStructure::Panel, theta::T...) where {T<:Real}

    parameterDraws = getSimulatedParameters(parameters, randomParameters, randomDraws, theta...)
    normalDraws = getStandardNormalDraws(randomDraws)

    paramsInOrder = [parameters..., keys(randomParameters)...]

    numDraws = size(normalDraws)[2]

    for i in 1:numDraws
        logit!(logitVals, i, data, paramsInOrder, nothing, nothing, choiceSetIDColumn, observationIDColumn, modelSpace, dataStructure, parameterDraws[:, i]...)
    end

end

function logLikelihood(logit::Vector{L}, choices::Vector{C}, observationIDs::Vector{O}) where {L<:Real, C<:Real, O<:Any}

    return sum(log.(logit).*choices)

end

function logLikelihood(logit::Matrix{L}, choices::Vector{C}, observationIDs::Vector{O}) where {L<:Real, C<:Real, O<:Any}
    logitDF::DataFrame = getIndividualLogits(logit, choices, observationIDs)
    meanLogit::Matrix{L} = mean(Matrix(logitDF[:, 2:size(logitDF)[2]]), dims=2)

    return sum(log.(meanLogit))

end

function logLikelihoodLM(logit::Vector{L}, choices::Vector{C}, observationIDs::Vector{O}, df::Nothing) where {L<:Real, C<:Real, O<:Any}

    return logLikelihood(logit, choices, observationIDs)

end

function logLikelihoodLM(logit::Matrix{L}, choices::Vector{C}, observationIDs::Vector{O}, df::DataFrame) where {L<:Real, C<:Real, O<:Any}

    meanLogit::Matrix{L} = mean(Matrix(getIndividualLogitsLM(logit, choices, observationIDs, df)[:, 2:size(df)[2]-1]), dims=2)

    return sum(log.(meanLogit))

end

function getIndividualLogits(logit::Matrix{L}, choices::Vector{C}, observationIDs::Vector{O}) where {L<:Real, C<:Real, O<:Any}

    logitCopy = deepcopy(logit)

    for i in 1:size(logitCopy)[2]
        logitCopy[:, i] = logitCopy[:, i].^choices
    end

    df::DataFrame = DataFrame(logitCopy, :auto)
    obsNames = Symbol.([names(df)...])
    df[:, :obs] .= observationIDs

    gData::GroupedDataFrame = groupby(df, :obs)
    df = combine(gData, obsNames .=> prod)

    return df

end

function getIndividualLogitsLM(logit::Matrix{L}, choices::Vector{C}, observationIDs::Vector{O}, df::DataFrame) where {L<:Real, C<:Real, O<:Any}

    for i in 1:size(logit)[2]
        logit[:, i] = logit[:, i].^choices
        df[:, Symbol(string("x", i))] .= logit[:, i]
    end

    obsNames = Symbol.([names(df)...])

    df[:, :obs] .= observationIDs

    gData::GroupedDataFrame = groupby(df, :obs)

    return combine(gData, obsNames[2:end] .=> prod)

end

function logLikelihoodGrad(logit::Vector{L}, choices::Vector{C}, xMat::Matrix{M}, observationIDs::Vector{O}, parameters::Vector{Symbol}, randomParameters::Nothing, randomDraws::Nothing, theta::L...) where{L<:Real, C<:Real, M<:Real, O<:Any}

    res::Vector{L} = choices.-logit

    gradient::Vector{L} = transpose(xMat)*res

    return gradient

end

function logLikelihoodGrad(logit::Matrix{L}, choices::Vector{C}, xMat::Matrix{M}, observationIDs::Vector{O}, parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, randomDraws::Matrix{L}, theta::L...) where{L<:Real, C<:Real, M<:Real, O<:Any}

    res::Matrix{L} = Matrix{L}(undef, size(logit)...)

    for i in 1:size(logit)[2]
        res[:, i] = choices.-logit[:, i]
    end

    fixedInds = 1:length(parameters)
    randomInds = (length(parameters)+1):size(xMat)[2]

    individualLogits::DataFrame = getIndividualLogits(logit, choices, observationIDs)
    meanLogit::Matrix{L} = mean(Matrix(individualLogits[:, 2:size(individualLogits)[2]]), dims=2)
    ll::L = sum(log.(meanLogit))

    normalDraws::Matrix{L} = getStandardNormalDraws(randomDraws)
    grad = repeat([0.], length(theta))

    #Perform calculations once in order to pre-alloc memory for temporary variables in loops below
    #outer loop
    tempObs::O = individualLogits[1,:obs]
    tempInds::BitVector = observationIDs.==tempObs
    tempFixedXMat::Matrix{L} = xMat[tempInds, :]
    indGrad = repeat([0.], length(theta))
    meanIndLogit::L = 0.

    #inner loop
    i::Int = 1
    tempRes::Vector{L} = res[tempInds, i]
    tempLogit::L = individualLogits[1, Symbol(string("x", i, "_prod"))]
    tempNormalDrawMat::Matrix{L} = transpose(reshape(repeat(normalDraws[:, i], size(tempFixedXMat)[1]), length(randomParameters), size(tempFixedXMat)[1]))
    partialMat::Matrix{L} = hcat(tempFixedXMat, tempFixedXMat[:, randomInds].*tempNormalDrawMat)
    for tempRow in eachrow(individualLogits)
        tempObs = tempRow[:obs]
        tempInds = observationIDs.==tempObs
        tempFixedXMat = xMat[tempInds, :]
        indGrad = repeat([0.], length(theta))
        meanIndLogit = 0.
        for i in 1:size(logit)[2]
            tempRes = res[tempInds, i]
            tempLogit = tempRow[Symbol(string("x", i, "_prod"))]
            meanIndLogit+=tempLogit
            tempNormalDrawMat= transpose(reshape(repeat(normalDraws[:, i], size(tempFixedXMat)[1]), length(randomParameters), size(tempFixedXMat)[1]))

            partialMat = hcat(tempFixedXMat, tempFixedXMat[:, randomInds].*tempNormalDrawMat)

            indGrad .= indGrad.+(tempLogit.*(transpose(partialMat)*tempRes))

        end
        grad = grad.+(indGrad./(meanIndLogit/size(logit)[2]))
    end

    grad = grad./(size(logit)[2])

    return grad

end

function constructSubGradients(res::Vector{L}, xMat::Matrix{M}, observationIDs::Vector{O}, parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, normalDraws::Vector{L}) where {L<:Real, M<:Real, O<:Any}

    for i in unique(observationIDs)
        tempInds = observationIDs.==i
        tempRes = res[tempInds, :]
        tempXMat = xMat[tempInds, :]
        tempNormalDrawMat = reshape(repeat(normalDraws, size(tempXMat)[1]), :, 2)
    end

end
function getDraws(k::Int, n::Int; RandomDraws::Bool=false)

    qmc = pyimport("scipy.stats.qmc")

    sobol = qmc.Sobol(k, seed=1995)
    draws = sobol.random(n)

    return transpose(draws)

end

function getDraws(randomParameters::Dict{Symbol, String}, nInterDraws::Int; RandomDraws::Bool=false)

    return Matrix(getDraws(length(randomParameters), nInterDraws; RandomDraws=RandomDraws))

end

function getDraws(randomParameters::Nothing, nInterDraws::Int; RandomDraws::Bool=false)

    return nothing

end

function getSimulatedParameters(parameters::Vector{Symbol}, randomParameters::Dict{Symbol, String}, randomDraws::Matrix{T}, theta::T...) where{T<:Real}

    numDraws = size(randomDraws)[2]
    numFixed = length(parameters)
    numRandom = length(randomParameters)*2

    fixedParams = theta[1:numFixed]

    params = reshape(repeat([fixedParams...], numDraws), numFixed, :)

    i = 1
    tempParams::Vector{T} = quantile.(Normal(), randomDraws[i, :])
    tempDist::String = "N"
    tempMean::T = theta[numFixed+i]
    tempSD::T = theta[numFixed+i+length(randomParameters)]
    for param in keys(randomParameters)
        tempDist = randomParameters[param]
        tempMean = theta[numFixed+i]
        tempSD = theta[numFixed+i+length(randomParameters)]

        tempParams = quantile.(Normal(tempMean, abs(tempSD)), randomDraws[i, :])
        params = vcat(params, transpose(tempParams))

        i+=1
    end

    return params

end

function getStandardNormalDraws(randomDraws::Matrix{T}) where{T<:Real}

    return quantile.(Normal(), randomDraws)

end

function getStandardNormalDraws(randomDraws::Nothing)

    return nothing

end

function getStartingValues(modelInfo::ModelInfo)

    if(modelInfo.RandomParameters==nothing)

        return repeat([0.], length(modelInfo.Parameters))
    else
        tempFixedParameters = deepcopy(modelInfo.Parameters)

        for tempVar in keys(modelInfo.RandomParameters)
            push!(tempFixedParameters, tempVar)
        end

        simpleModelInfo = ModelInfo(modelInfo.Data, modelInfo.ChoiceColumn, tempFixedParameters, modelInfo.ChoiceSetIDColumn, nothing, nothing, modelInfo.PanelIDColumn, modelInfo.ObservationIDColumn, modelInfo.DataStructure, modelInfo.ModelSpace, modelInfo.Name, modelInfo.InterDrawCount, modelInfo.WeightColumn, modelInfo.TimeCreated, modelInfo.PreferenceSearchBounds)
        tempRes = EstimateModel(simpleModelInfo; post=false)
        numRandParams::Int = length(modelInfo.RandomParameters)

        startSigmas = repeat([0.3], numRandParams)

        startValues = tempRes["coef"]

        push!(startValues, startSigmas...)

        return startValues

    end
end

max_stars_repo_name: crforsythe/Joice-Releases

max_issues_repo_name: crforsythe/Joice-Releases

max_forks_repo_name: crforsythe/Joice-Releases


text:
"""
    Torsion(i, j, k, l, f1, f2, f3, f4)

A dihedral torsion angle between four atoms.
"""
struct Torsion{T} <: SpecificInteraction
    i::Int
    j::Int
    k::Int
    l::Int
    f1::T
    f2::T
    f3::T
    f4::T
end

@inline @inbounds function force(d::Torsion,
                                  coords,
                                  s::Simulation)
    ba = vector(coords[d.j], coords[d.i], s.box_size)
    bc = vector(coords[d.j], coords[d.k], s.box_size)
    dc = vector(coords[d.l], coords[d.k], s.box_size)
    p1 = normalize(ba × bc)
    p2 = normalize(-dc × -bc)
    θ = atan(dot((-ba × bc) × (bc × -dc), normalize(bc)), dot(-ba × bc, bc × -dc))
    angle_term = (d.f1*sin(θ) - 2*d.f2*sin(2*θ) + 3*d.f3*sin(3*θ)) / 2
    fa = (angle_term / (norm(ba) * sin(acosbound(dot(ba, bc) / (norm(ba) * norm(bc)))))) * p1
    # fd clashes with a function name
    f_d = (angle_term / (norm(dc) * sin(acosbound(dot(bc, dc) / (norm(bc) * norm(dc)))))) * p2
    oc = bc / 2
    tc = -(oc × f_d + (-dc × f_d) / 2 + (ba × fa) / 2)
    fc = (1 / dot(oc, oc)) * (tc × oc)
    fb = -fa - fc - f_d
    return [d.i, d.j, d.k, d.l], [fa, fb, fc, f_d]
end

@inline @inbounds function potential_energy(d::Torsion,
                                            s::Simulation)
    ba = vector(s.coords[d.j], s.coords[d.i], s.box_size)
    bc = vector(s.coords[d.j], s.coords[d.k], s.box_size)
    dc = vector(s.coords[d.l], s.coords[d.k], s.box_size)

    θ = atan(dot((-ba × bc) × (bc × -dc), normalize(bc)), dot(-ba × bc, bc × -dc))

    (d.f1 * (1 + cos(θ)) + d.f2 * (1 - cos(2θ)) + d.f3 * (1 + cos(3θ)) + d.f4) / 2
end


text:
export distance

"""
    distance(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle;
             [p]::Real=2)

Compute the standard distance between two hyperrectangular sets, defined as

```math
    \\inf_{x \\in H_1, y \\in H_2} \\{ d(x, y) \\}.
```

### Input

- `H1` -- hyperrectangular set
- `H2` -- hyperrectangular set
- `p`  -- (optional; default: `2`) value of the ``p``-norm

### Output

The distance, which is zero if the sets intersect and otherwise the ``p``-norm
of the shortest line segment between any pair of points.

### Notes

See also [`hausdorff_distance`](@ref) for an alternative distance notion.
"""
function distance(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle;
                  p::Real=2)
    n = dim(H1)
    @assert n == dim(H2) "incompatible set dimensions $n and $(dim(H2))"

    N = promote_type(eltype(H1), eltype(H2))
    d = Vector{N}(undef, n)
    @inbounds for i in 1:n
        # find relative position in dimension i
        # (if c1 == c2, the results are equivalent independent of the branch)
        if center(H1, i) >= center(H2, i)
            lhs = low(H1, i)
            rhs = high(H2, i)
        else
            lhs = low(H2, i)
            rhs = high(H1, i)
        end
        if _leq(lhs, rhs)
            d[i] = zero(N)
        else
            d[i] = rhs - lhs
        end
    end
    return norm(d, p)
end

"""
    distance(S::AbstractSingleton, X::LazySet; [p]::Real=2.0)

Compute the distance between the singleton `S` and the set `X` with respect to
the given `p`-norm.

### Input

- `S` -- singleton, i.e., a set with one element
- `X` -- set
- `p` -- (optional, default: `2.0`) the `p`-norm used; `p = 2.0` corresponds to
         the usual Euclidean norm

### Output

A scalar representing the distance between the element wrapped by `S` and the
set `X`.
"""
function distance(S::AbstractSingleton, X::LazySet; p::Real=2.0)
    return distance(element(S), X; p=p)
end

distance(X::LazySet, S::AbstractSingleton; p::Real=2.0) = distance(S, X; p=p)

distance(S1::AbstractSingleton, S2::AbstractSingleton; p::Real=2.0) =
    distance(element(S1), element(S2); p=p)

# disambiguation
distance(S::AbstractSingleton, H::AbstractHyperrectangle; p::Real=2.0) =
    distance(element(S), H; p=p)
distance(H::AbstractHyperrectangle, S::AbstractSingleton; p::Real=2.0) =
    distance(element(S), H; p=p)


text:
function reflect(@nospecialize(func), @nospecialize(types);
                 optimize::Bool=true, run_enzyme::Bool=true, second_stage::Bool=true, split::Bool=false)
    primal, adjoint, rt = fspec(func, types)

    target = Compiler.EnzymeTarget()
    params = Compiler.EnzymeCompilerParams()
    job    = Compiler.CompilerJob(target, primal, params)

    # Codegen the primal function and all its dependency in one module
    mod, primalf = Compiler.codegen(:llvm, job, optimize=false, #= validate=false =#)

    # Run pipeline and Enzyme pass
    if optimize
        optimize!(mod)
    end

    if run_enzyme
        annotate!(mod)
        adjointf, augmented_primalf = enzyme!(mod, primalf, adjoint, rt, split)

        if second_stage
            post_optimze!(mod)
        end
        llvmf = adjointf
    else
        llvmf = primalf
    end

    return llvmf, mod
end

function enzyme_code_llvm(io::IO, @nospecialize(func), @nospecialize(types);
                          optimize::Bool=true, run_enzyme::Bool=true, second_stage::Bool=true,
                          raw::Bool=false, debuginfo::Symbol=:default, dump_module::Bool=false)
    llvmf, mod = reflect(func, types, optimize=optimize, run_enzyme=run_enzyme, second_stage=second_stage)

    str = ccall(:jl_dump_function_ir, Ref{String},
                (LLVM.API.LLVMValueRef, Bool, Bool, Ptr{UInt8}),
                llvmf, !raw, dump_module, debuginfo)
    print(io, str)
end
enzyme_code_llvm(@nospecialize(func), @nospecialize(types); kwargs...) = enzyme_code_llvm(stdout, func, types; kwargs...)

function enzyme_code_native(io::IO, @nospecialize(func), @nospecialize(types))
    llvmf, mod = reflect(func, types)
    str = String(LLVM.emit(tm[], mod, LLVM.API.LLVMAssemblyFile))
    print(io, str)
end
enzyme_code_native(@nospecialize(func), @nospecialize(types); kwargs...) = enzyme_code_native(stdout, func, types; kwargs...)


text:
export QuadDirect

struct QuadDirect end

function __map_optimizer_args(prob::OptimizationProblem, opt::QuadDirect;
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)
    if !isnothing(maxtime)
        @warn "common maxtime is currently not used by $(opt)"
    end

    if !isnothing(cb)
        @warn "callbacks are currently not used by $(opt)"
    end

    mapped_args = (; kwargs...)

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., maxevals=maxiters)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., atol=abstol)
    end

    if !isnothing(reltol)
        mapped_args = (; mapped_args..., rtol=reltol)
    end

    return mapped_args
end

function __solve(prob::OptimizationProblem, opt::QuadDirect;
    splits = nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    local x, _loss

    maxiters = _check_and_convert_maxiters(maxiters)

    if splits === nothing
        error("You must provide the initial locations at which to evaluate the function in `splits` (a list of 3-vectors with values in strictly increasing order and within the specified bounds).")
    end

    _loss = function(θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    opt_arg = __map_optimizer_args(prob, opt; maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol, kwargs...)
    t0 = time()
    # root, x0 = !(isnothing(maxiters)) ? QuadDIRECT.analyze(_loss, splits, prob.lb, prob.ub; maxevals = maxiters, kwargs...) : QuadDIRECT.analyze(_loss, splits, prob.lb, prob.ub; kwargs...)
    root, x0 = QuadDIRECT.analyze(_loss, splits, prob.lb, prob.ub; opt_arg...)
    box = minimum(root)
    t1 = time()

    SciMLBase.build_solution(prob, opt, QuadDIRECT.position(box, x0), QuadDIRECT.value(box); original=root)
end


text:
function autorange(img, tform)
    R = CartesianIndices(axes(img))
    autorange(R, tform)
end

function autorange(R::CartesianIndices, tform)
    tform = _round(tform)
    mn = mx = tform(SVector(first(R).I))
    for I in CornerIterator(R)
        x = tform(SVector(I.I))
        # we map min and max to prevent type-inference issues
        # (because min(::SVector,::SVector) -> Vector)
        mn = map(min, x, mn)
        mx = map(max, x, mx)
    end
    _autorange(Tuple(mn), Tuple(mx))
end

@noinline _autorange(mn,mx) = map((a,b)->floor(Int,a):ceil(Int,b), mn, mx)

## Iterate over the corner-indices of a rectangular region
struct CornerIterator{I<:CartesianIndex}
    start::I
    stop::I
end
CornerIterator(R::CartesianIndices) = CornerIterator(first(R), last(R))

eltype(::Type{CornerIterator{I}}) where {I} = I
Base.Iterators.IteratorSize(::Type{CornerIterator{I}}) where {I<:CartesianIndex{N}} where {N} = Base.Iterators.HasShape{N}()

# in 0.6 we could write: 1 .+ (iter.stop.I .- iter.start.I .!= 0)
size(iter::CornerIterator{CartesianIndex{N}}) where {N} = ntuple(d->iter.stop.I[d]-iter.start.I[d]==0 ? 1 : 2, Val(N))::NTuple{N,Int}
length(iter::CornerIterator) = prod(size(iter))

@inline function Base.iterate(iter::CornerIterator{<:CartesianIndex})
    if any(map(>, iter.start.I, iter.stop.I))
        return nothing
    end
    iterate(iter, iter.start)
end

@inline function Base.iterate(iter::CornerIterator{I}, state) where I<:CartesianIndex
    if state.I[end] > iter.stop.I[end]
        return nothing
    end
    state, I(inc(state.I, iter.start.I, iter.stop.I))
end

# increment & carry
@inline inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
# the max(1, ...) makes sure that the code doesn't break
# for corner cases where if stop == start
# (e.g. CornerIterator(CartesianRange((1,1))))
@inline inc(state::Tuple{Int}, start::Tuple{Int}, stop::Tuple{Int}) = (state[1]+(max(1,stop[1]-start[1])),)
@inline function inc(state, start, stop)
    if state[1] < stop[1]
        return (stop[1],tail(state)...)
    end
    newtail = inc(tail(state), tail(start), tail(stop))
    (start[1], newtail...)
end

# 0-d is special-cased to iterate once and only once
Base.iterate(iter::CornerIterator{<:CartesianIndex{0}}) = (iter.start, nothing)
Base.iterate(iter::CornerIterator{<:CartesianIndex{0}}, state) = nothing

# Ensure that we use StaticArrays for linear transformations
try_static(tfm, img) = tfm   # fallback
try_static(tfm::AffineMap{<:SMatrix, <:SVector}, img::AbstractArray{T,N}) where {T,N} = tfm
try_static(tfm::AffineMap{<:AbstractMatrix, <:AbstractVector}, img::AbstractArray{T,N}) where {T,N} =
    AffineMap(SMatrix{N,N}(tfm.linear), SVector{N}(tfm.translation))
try_static(tfm::LinearMap{<:SMatrix}, img::AbstractArray{T,N}) where {T,N} = tfm
try_static(tfm::LinearMap{<:AbstractMatrix}, img::AbstractArray{T,N}) where {T,N} =
    LinearMap(SMatrix{N,N}(tfm.linear))
try_static(tfm::Translation{<:SVector}, img::AbstractArray{T,N}) where {T,N} = tfm
try_static(tfm::Translation{<:AbstractVector}, img::AbstractArray{T,N}) where {T,N} =
    Translation(SVector{N}(tfm.translation))

# Slightly round/discretize the transformation so that the warpped image size isn't affected by
# numerical stability
# https://github.com/JuliaImages/ImageTransformations.jl/issues/104
_default_digits(::Type{T}) where T<:Number = _default_digits(floattype(T))
# these constants come from eps() digits
_default_digits(::Type{<:AbstractFloat}) = 15
_default_digits(::Type{Float64}) = 15
_default_digits(::Type{Float32}) = 7

function _round(tform::T; kwargs...) where T<:CoordinateTransformations.Transformation
    rounded_fields = map(Base.OneTo(fieldcount(T))) do i
        __round(getfield(tform, i); kwargs...)
    end
    T(rounded_fields...)
end
if isdefined(Base, :ComposedFunction)
    _round(tform::ComposedFunction; kwargs...) = _round(tform.outer; kwargs...) ∘ _round(tform.inner; kwargs...)
end
_round(tform; kwargs...) = tform

__round(x; kwargs...) = x
__round(x::AbstractArray; digits=_default_digits(eltype(x)), kwargs...) = round.(x; digits=digits, kwargs...)
__round(x::T; digits=_default_digits(T), kwargs...) where T<:Number = round(x; digits=digits, kwargs...)


text:
@testset "Print" begin
    buf = IOBuffer()
    print(buf, ReferenceSequence(dna""))
    @test String(take!(buf)) == ""

    print(buf, ReferenceSequence(dna"ACGTN"))
    @test String(take!(buf)) == "ACGTN"

    print(buf, ReferenceSequence(dna"A"^100))
    @test String(take!(buf)) == "A"^100

    print(buf, ReferenceSequence(dna"A"^100), width=80)
    @test String(take!(buf)) == "A"^80 * '\n' * "A"^20

    print(buf, ReferenceSequence(dna"A"^20), width=20)
    @test String(take!(buf)) == "A"^20

    # Check more than 4096
    print(buf, ReferenceSequence(dna"A"^4101), width=100)
    @test String(take!(buf)) == repeat("A"^100 * '\n', 41) * "A"
end

@testset "Show" begin
    buf = IOBuffer()
    show(buf, ReferenceSequence(dna""))
    @test String(take!(buf)) == "< EMPTY SEQUENCE >"

    show(buf, ReferenceSequence(dna"ACGTN"))
    @test String(take!(buf)) == "ACGTN"
end


text:
typealias BigFloats Union{BigFloat,Complex{BigFloat}}

real(x...)=Base.real(x...)
real{T<:Real}(::Type{T})=T
real{T<:Real}(::Type{Complex{T}})=T

# The following implements Bluestein's algorithm, following http://www.dsprelated.com/dspbooks/mdft/Bluestein_s_FFT_Algorithm.html
# To add more types, add them in the union of the function's signature.
function Base.fft{T<:BigFloats}(x::Vector{T})
    n = length(x)
    ispow2(n) && return fft_pow2(x)
    ks = linspace(zero(real(T)),n-one(real(T)),n)
    Wks = exp(-im*convert(T,π)*ks.^2/n)
    xq,wq = x.*Wks,conj([exp(-im*convert(T,π)*n);reverse(Wks);Wks[2:end]])
    return Wks.*conv(xq,wq)[n+1:2n]
end

function Base.fft!{T<:BigFloats}(x::Vector{T})
    x[:] = fft(x)
    return x
end

# add rfft for BigFloat, by calling fft
#  this creates ToeplitzMatrices.rfft, so avoids changing Base.rfft
Base.rfft{T<:BigFloats}(v::Vector{T})=fft(v)[1:div(length(v),2)+1]
function Base.irfft{T<:BigFloats}(v::Vector{T},n::Integer)
    @assert n==2length(v)-1
    r=Array(Complex{BigFloat},n)
    r[1:length(v)]=v
    r[length(v)+1:end]=reverse(conj(v[2:end]))
    real(ifft(r))
end

Base.ifft{T<:BigFloats}(x::Vector{T}) = conj!(fft(conj(x)))/length(x)
function Base.ifft!{T<:BigFloats}(x::Vector{T})
    x[:] = ifft(x)
    return x
end

function Base.conv{T<:BigFloats}(u::StridedVector{T}, v::StridedVector{T})
    nu,nv = length(u),length(v)
    n = nu + nv - 1
    np2 = nextpow2(n)
    append!(u,zeros(T,np2-nu)),append!(v,zeros(T,np2-nv))
    y = ifft_pow2(fft_pow2(u).*fft_pow2(v))
    #TODO This would not handle Dual/ComplexDual numbers correctly
    y = T<:Real ? real(y[1:n]) : y[1:n]
end

# This is a Cooley-Tukey FFT algorithm inspired by many widely available algorithms including:
# c_radix2.c in the GNU Scientific Library and four1 in the Numerical Recipes in C.
# However, the trigonometric recurrence is improved for greater efficiency.
# The algorithm starts with bit-reversal, then divides and conquers in-place.
function fft_pow2!{T<:BigFloat}(x::Vector{T})
    n,big2=length(x),2one(T)
    nn,j=n÷2,1
    for i=1:2:n-1
        if j>i
            x[j], x[i] = x[i], x[j]
            x[j+1], x[i+1] = x[i+1], x[j+1]
        end
        m = nn
        while m ≥ 2 && j > m
            j -= m
            m = m÷2
        end
        j += m
    end
    logn = 2
    while logn < n
        θ=-big2/logn
        wtemp = sinpi(θ/2)
        wpr, wpi = -2wtemp^2, sinpi(θ)
        wr, wi = one(T), zero(T)
        for m=1:2:logn-1
            for i=m:2logn:n
                j=i+logn
                mixr, mixi = wr*x[j]-wi*x[j+1], wr*x[j+1]+wi*x[j]
                x[j], x[j+1] = x[i]-mixr, x[i+1]-mixi
                x[i], x[i+1] = x[i]+mixr, x[i+1]+mixi
            end
            wr = (wtemp=wr)*wpr-wi*wpi+wr
            wi = wi*wpr+wtemp*wpi+wi
        end
        logn = logn << 1
    end
    return x
end

function fft_pow2{T<:BigFloat}(x::Vector{Complex{T}})
    y = interlace(real(x),imag(x))
    fft_pow2!(y)
    return complex(y[1:2:end],y[2:2:end])
end
fft_pow2{T<:BigFloat}(x::Vector{T}) = fft_pow2(complex(x))

function ifft_pow2{T<:BigFloat}(x::Vector{Complex{T}})
    y = interlace(real(x),-imag(x))
    fft_pow2!(y)
    return complex(y[1:2:end],-y[2:2:end])/length(x)
end

function Base.dct(a::AbstractArray{Complex{BigFloat}})
	N = big(length(a))
    c = fft([a; flipdim(a,1)])
    d = c[1:N] .* exp(-im*big(pi)*(0:N-1)/(2*N))
    d[1] = d[1] / sqrt(big(2))
    scale!(inv(sqrt(2N)), d)
end

Base.dct(a::AbstractArray{BigFloat}) = real(dct(complex(a)))

function Base.idct(a::AbstractArray{Complex{BigFloat}})
    N = big(length(a))
    b = a * sqrt(2*N)
    b[1] = b[1] * sqrt(big(2))
    b = b .* exp(im*big(pi)*(0:N-1)/(2*N))
    b = [b; 0; conj(flipdim(b[2:end],1))]
    c = ifft(b)
    c[1:N]
end

Base.idct(a::AbstractArray{BigFloat}) = real(idct(complex(a)))

Base.dct!{T<:BigFloats}(a::AbstractArray{T}) = (b = dct(a); a[:] = b)
Base.idct!{T<:BigFloats}(a::AbstractArray{T}) = (b = idct(a); a[:] = b)

# dummy plans
type DummyFFTPlan{T,inplace} <: Base.DFT.Plan{T} end
type DummyiFFTPlan{T,inplace} <: Base.DFT.Plan{T} end
type DummyrFFTPlan{T,inplace} <: Base.DFT.Plan{T} end
type DummyirFFTPlan{T,inplace} <: Base.DFT.Plan{T} end
type DummyDCTPlan{T,inplace} <: Base.DFT.Plan{T} end
type DummyiDCTPlan{T,inplace} <: Base.DFT.Plan{T} end

for (Plan,iPlan) in ((:DummyFFTPlan,:DummyiFFTPlan),
                     (:DummyrFFTPlan,:DummyirFFTPlan),
                     (:DummyDCTPlan,:DummyiDCTPlan))
   @eval begin
       Base.inv{T,inplace}(::$Plan{T,inplace})=$iPlan{T,inplace}()
       Base.inv{T,inplace}(::$iPlan{T,inplace})=$Plan{T,inplace}()
    end
end

for (Plan,ff,ff!) in ((:DummyFFTPlan,:fft,:fft!),
                      (:DummyiFFTPlan,:ifft,:ifft!),
                      (:DummyrFFTPlan,:rfft,:rfft!),
                      (:DummyirFFTPlan,:irfft,:irfft!),
                      (:DummyDCTPlan,:dct,:dct!),
                      (:DummyiDCTPlan,:idct,:idct!))
    @eval begin
        *{T,N}(p::$Plan{T,true}, x::StridedArray{T,N})=$ff!(x)
        *{T,N}(p::$Plan{T,false}, x::StridedArray{T,N})=$ff(x)
        function Base.A_mul_B!(C::StridedVector,p::$Plan,x::StridedVector)
            C[:]=$ff(x)
            C
        end
    end
end

Base.plan_fft!{T<:BigFloats}(x::Vector{T}) = DummyFFTPlan{Complex{BigFloat},true}()
Base.plan_ifft!{T<:BigFloats}(x::Vector{T}) = DummyiFFTPlan{Complex{BigFloat},true}()
#Base.plan_rfft!{T<:BigFloats}(x::Vector{T}) = DummyrFFTPlan{Complex{BigFloat},true}()
#Base.plan_irfft!{T<:BigFloats}(x::Vector{T},n::Integer) = DummyirFFTPlan{Complex{BigFloat},true}()
Base.plan_dct!{T<:BigFloats}(x::Vector{T}) = DummyDCTPlan{T,true}()
Base.plan_idct!{T<:BigFloats}(x::Vector{T}) = DummyiDCTPlan{T,true}()

Base.plan_fft{T<:BigFloats}(x::Vector{T}) = DummyFFTPlan{Complex{BigFloat},false}()
Base.plan_ifft{T<:BigFloats}(x::Vector{T}) = DummyiFFTPlan{Complex{BigFloat},false}()
Base.plan_rfft{T<:BigFloats}(x::Vector{T}) = DummyrFFTPlan{Complex{BigFloat},false}()
Base.plan_irfft{T<:BigFloats}(x::Vector{T},n::Integer) = DummyirFFTPlan{Complex{BigFloat},false}()
Base.plan_dct{T<:BigFloats}(x::Vector{T}) = DummyDCTPlan{T,false}()
Base.plan_idct{T<:BigFloats}(x::Vector{T}) = DummyiDCTPlan{T,false}()

function interlace{S<:Number,V<:Number}(a::Vector{S},b::Vector{V})
    na=length(a);nb=length(b)
    T=promote_type(S,V)
    if nb≥na
        ret=zeros(T,2nb)
        ret[1:2:1+2*(na-1)]=a
        ret[2:2:end]=b
        ret
    else
        ret=zeros(T,2na-1)
        ret[1:2:end]=a
        if !isempty(b)
            ret[2:2:2+2*(nb-1)]=b
        end
        ret
    end
end