text:
using ToyFHE
using ToyFHE: degree
using ToyFHE.NTT
using GaloisFields
using OffsetArrays
using Test
using Primes
using BitIntegers

params = BFVParams(
    65537, # plaintext modulus
    ; eval_mult_count = 1
)
kp = keygen(params)

plain = SlotEncoding(zero(plaintext_space(params)))
plain[0] = 1
plain[1] = 1

plain2 = SlotEncoding(zero(plaintext_space(params)))
plain2[:] .= 10
plain2[0] = 5

c1 = encrypt(kp, plain)
c2 = encrypt(kp, plain2)
y = c1*c2
let data = SlotEncoding(decrypt(kp, y))
    @test data[0] == 5
    @test data[1] == 10
    @test all(iszero, data[2:end])
end


text:
"""
	las2pointcloud(fnames::String...) -> PointCloud

Read more than one file `.las` and extrapolate the LAR model and the color of each point.
"""
function las2pointcloud(fnames::String...)::PointCloud
	Vtot = Array{Float64,2}(undef, 3, 0)
	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
	l = Threads.ReentrantLock()

	Threads.@threads for fname in fnames
		partialV = las2larpoints(fname)
		partialRGB = las2color(fname)
		Threads.lock(l)
		Vtot = hcat(Vtot, partialV)
		rgbtot = hcat(rgbtot, partialRGB)
		Threads.unlock(l)
	end
	return PointCloud(Vtot,rgbtot)
end

"""
	las2larpoints(file::String) -> Lar.Points

Return coordinates of points in LAS file.
"""
function las2larpoints(file::String)::Lar.Points
	header, laspoints = read_LAS_LAZ(file)
	npoints = length(laspoints)
	x = [LasIO.xcoord(laspoints[k], header) for k in 1:npoints]
	y = [LasIO.ycoord(laspoints[k], header) for k in 1:npoints]
	z = [LasIO.zcoord(laspoints[k], header) for k in 1:npoints]
	return vcat(x',y',z')
end

"""
	las2aabb(file::String) -> AABB

Return LAS file's bounding box.
"""
function las2aabb(file::String)::AABB
	header, p = read_LAS_LAZ(file)
	aabb = LasIO.boundingbox(header)
	return AABB(aabb.xmax, aabb.xmin, aabb.ymax, aabb.ymin, aabb.zmax, aabb.zmin)
end

"""
	las2color(file::String)::Lar.Points

Return color, rgb, associated to each point in LAS file.
"""
function las2color(file::String)::Lar.Points
	header, laspoints =  read_LAS_LAZ(file)
	npoints = length(laspoints)
	type = LasIO.pointformat(header)
	if type != LasPoint0 && type != LasPoint1
		r = LasIO.ColorTypes.red.(laspoints)
		g = LasIO.ColorTypes.green.(laspoints)
		b = LasIO.ColorTypes.blue.(laspoints)
		return vcat(r',g',b')
	end
	return rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
end

"""
	color(p::LasPoint, header::LasHeader)

Return color of one point in LAS file.
"""
function color(p::LasPoint, header::LasHeader)
	type = LasIO.pointformat(header)
	if type != LasPoint0 && type != LasPoint1
		r = LasIO.ColorTypes.red(p)
		g = LasIO.ColorTypes.green(p)
		b = LasIO.ColorTypes.blue(p)
		return vcat(r',g',b')
	end
	return rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
end

"""
	 xyz(p::LasPoint, h::LasHeader)

Return coords of this point `p`.
"""
function xyz(p::LasPoint, h::LasHeader)
	return [LasIO.xcoord(p, h); LasIO.ycoord(p, h); LasIO.zcoord(p, h)]
end

"""
	read_LAS_LAZ(fname::String)

Read point cloud files: LAS or LAZ.
"""
function read_LAS_LAZ(fname::String)
	if endswith(fname,".las")
		header, laspoints = LasIO.FileIO.load(fname)
	elseif endswith(fname,".laz")
		header, laspoints = LazIO.load(fname)
	else endswith(fname,".bin")
		header, laspoints = load_bin(fname)
	end
	return header,laspoints
end

ext: jl


text:
# src/runtests.jl
# Test Suite for AbstractRationals

using AbstractRationals


text:
using Flux, DiffEqFlux, StochasticDiffEq

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lotka_volterra_noise(du,u,p,t)
  du[1] = 0.01u[1]
  du[2] = 0.01u[2]
end
prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,10.0))
p = param([2.2, 1.0, 2.0, 0.4])
function predict_fd_sde()
  diffeq_fd(p,sol->sol[1,:],101,prob,SOSRI(),saveat=0.1)
end
loss_fd_sde() = sum(abs2,x-1 for x in predict_fd_sde())
loss_fd_sde()
Flux.back!(loss_fd_sde())

prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,0.5))
function predict_rd_sde()
  Tracker.collect(diffeq_rd(p,prob,SOSRI(),saveat=0.0:0.1:0.5))
end
loss_rd_sde() = sum(abs2,x-1 for x in predict_rd_sde())
loss_rd_sde()
Flux.back!(loss_rd_sde())


text:
function division(a, b)
	k=0
    while a*k<b
        k = k+1
    end
    return k
end

function syra(n::Int64)
	if n==1
        return 0
    elseif n%2 == 0
        return 1 + syra(division(2,n))
    else
        return 1 + syra(3*n+1)
    end
end

println(syra(42))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# ---
# title: 901. Online Stock Span
# id: problem901
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/online-stock-span/description/>
# hidden: true
# ---
#
# Write a class `StockSpanner` which collects daily price quotes for some stock,
# and returns the _span_  of that stock's price for the current day.
#
# The span of the stock's price today is defined as the maximum number of
# consecutive days (starting from today and going backwards) for which the price
# of the stock was less than or equal to today's price.
#
# For example, if the price of a stock over the next 7 days were `[100, 80, 60,
# 70, 60, 75, 85]`, then the stock spans would be `[1, 1, 1, 2, 1, 4, 6]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["StockSpanner","next","next","next","next","next","next","next"], [[],[100],[80],[60],[70],[60],[75],[85]]
#     Output: [null,1,1,1,2,1,4,6]
#     Explanation:
#     First, S = StockSpanner() is initialized.  Then:
#     S.next(100) is called and returns 1,
#     S.next(80) is called and returns 1,
#     S.next(60) is called and returns 1,
#     S.next(70) is called and returns 2,
#     S.next(60) is called and returns 1,
#     S.next(75) is called and returns 4,
#     S.next(85) is called and returns 6.
#
#     Note that (for example) S.next(75) returned 4, because the last 4 prices
#     (including today's price of 75) were less than or equal to today's price.
#
#
#
#
# **Note:**
#
#   1. Calls to `StockSpanner.next(int price)` will have `1 <= price <= 10^5`.
#   2. There will be at most `10000` calls to `StockSpanner.next` per test case.
#   3. There will be at most `150000` calls to `StockSpanner.next` across all test cases.
#   4. The total time limit for this problem has been reduced by 75% for C++, and 50% for all other languages.
#
#
## @lc code=start
using LeetCode

Base.@kwdef struct StockSpanner
    history::Vector{Pair{Int,Int}} = []
end

function next!(s::StockSpanner, price::Int)
    span = 1
    while !isempty(s.history) && s.history[end][1] <= price
        span += s.history[end][2]
        pop!(s.history)
    end
    push!(s.history, price => span)
    return span
end
## @lc code=end


text:
export StateCachedEnv

"""
Cache the state so that `state(env)` will always return the same result before
the next interaction with `env`. This function is useful because some
environments are stateful during each `state(env)`. For example:
`StateOverriddenEnv(StackFrames(...))`.
"""
mutable struct StateCachedEnv{S,E<:AbstractEnv} <: AbstractEnvWrapper
    s::S
    env::E
    is_state_cached::Bool
end

StateCachedEnv(env) = StateCachedEnv(state(env), env, true)

function (env::StateCachedEnv)(args...; kwargs...)
    env.env(args...; kwargs...)
    env.is_state_cached = false
end

function RLBase.state(env::StateCachedEnv, args...; kwargs...)
    if env.is_state_cached
        env.s
    else
        env.s = state(env.env, args...; kwargs...)
        env.is_state_cached = true
        env.s
    end
end


text:
function generate_symmetric_distance_matrix(xs::Vector{T},
                                            ys::Vector{T}) where T <: Real
    @assert length(xs) == length(ys)
    n = length(xs)
    matrix = Matrix{Float64}(undef, n, n)
    for j in 1:n
        matrix[j,j] = 0.0
        for i in j+1:n
            matrix[i,j] = matrix[j,i] = sqrt((xs[j]-xs[i])^2 + (ys[j]-ys[i])^2)
        end
    end
    return matrix
end

function intersect(tw1::Range, tw2::Range)
    if max(tw1.lb,tw2.lb) > min(tw1.ub,tw2.ub)
        return nothing
    else
        return Range(max(tw1.lb,tw2.lb), min(tw1.ub,tw2.ub))
    end
end

function ranges_intersection(tws_1::Vector{Range}, tws_2::Vector{Range})
    intersection_vec = Range[]
    for tw1 in tws_1
        for tw2 in tws_2
            intersection = intersect(tw1, tw2)
            if intersection != nothing
                push!(intersection_vec, intersection)
            end
        end
    end
    return intersection_vec
end

function set_indices(data::RvrpInstance)
    for loc_idx in 1:length(data.locations)
        data.locations[loc_idx].index = loc_idx
    end
end

function build_computed_data(data::RvrpInstance)
    location_id_2_index = Dict{String,Int}(
        data.locations[i].id => i for i in 1:length(data.locations)
    )
    location_group_id_2_index = Dict{String,Int}(
        data.location_groups[i].id => i for i in 1:length(data.location_groups)
    )
    product_specification_class_id_2_index = Dict{String,Int}(
        data.product_specification_classes[i].id =>
        i for i in 1:length(data.product_specification_classes)
    )
    request_id_2_index = Dict{String,Int}(
        data.requests[i].id =>
        i for i in 1:length(data.requests)
    )
    vehicle_category_id_2_index = Dict{String,Int}(
        data.vehicle_categories[i].id =>
        i for i in 1:length(data.vehicle_categories)
    )
    travel_specification_id_2_index = Dict{String,Int}(
        data.travel_specifications[i].id =>
        i for i in 1:length(data.travel_specifications)
    )
    vehicle_set_id_2_index = Dict{String,Int}(
        data.vehicle_sets[i].id =>
        i for i in 1:length(data.vehicle_sets)
    )
    capacity_id_2_index = Dict{String, Int}()
    property_id_2_index = Dict{String, Int}()
    cap_idx = 1
    prop_idx = 1
    for v_cat in data.vehicle_categories
        for (cap_id,v) in v_cat.capacities.of_vehicle
            if !haskey(capacity_id_2_index, cap_id)
                capacity_id_2_index[cap_id] = cap_idx
                cap_idx += 1
            end
        end
        for (prop_id,v) in v_cat.properties.of_vehicle
            if !haskey(property_id_2_index, prop_id)
                property_id_2_index[prop_id] = prop_idx
                prop_idx += 1
            end
        end
    end
    for p_spec_class in data.product_specification_classes
        for (cap_id,v) in p_spec_class.capacity_consumptions
            if !haskey(capacity_id_2_index, cap_id)
                capacity_id_2_index[cap_id] = cap_idx
                cap_idx += 1
            end
        end
        for (prop_id,v) in p_spec_class.property_requirements
            if !haskey(property_id_2_index, prop_id)
                property_id_2_index[prop_id] = prop_idx
                prop_idx += 1
            end
        end
    end
    return RvrpComputedData(
        location_id_2_index, location_group_id_2_index,
        product_specification_class_id_2_index, request_id_2_index,
        vehicle_category_id_2_index, vehicle_set_id_2_index,
        capacity_id_2_index, property_id_2_index,
        travel_specification_id_2_index, BitSet(), false
    )
end

function create_singleton_location_groups(locations::Vector{Location})
    return [LocationGroup(
        id = string(l.id, "_loc_group"),
        location_ids = [l.id]
    ) for l in locations]
end

function get_capacity_consumptions(req::Request, product_specification_classes::Vector{ProductSpecificationClass}, computed_data::RvrpComputedData)
    quantity = req.product_quantity_range.ub
    product_specification_class = product_specification_classes[computed_data.product_specification_class_id_2_index[req.product_specification_class_id]]
    c = Dict{String,Float64}(
        k => ceil(quantity/v[2]) * v[1]
        for (k,v) in product_specification_class.capacity_consumptions
    )
    return c
end

function preprocess_instance(data::RvrpInstance)
    # push!(data.locations, Location(id = "default_id"))
    # push!(data.location_groups, LocationGroup(id = "default_id"))
    push!(data.product_compatibility_classes,
          ProductCompatibilityClass(id = "default_id"))
    push!(data.product_sharing_classes,
          ProductSharingClass(id = "default_id"))
    push!(data.product_specification_classes,
          ProductSpecificationClass(id = "default_id"))
    # push!(data.requests, Request(id = "default_id"))
    push!(data.vehicle_categories, VehicleCategory(id = "default_id"))
    # push!(data.vehicle_sets, HomogeneousVehicleSet(id = "default_id"))
end

function check_initial_feasibility(data::RvrpInstance)
    v_sets = data.vehicle_sets
    nb_vehicles = 0
    for v_set in v_sets
        nb_vehicles += v_set.nb_of_vehicles_range.soft_range.ub
    end
    if (nb_vehicles == 0
        || length(data.requests) == 0
        || length(data.travel_specifications) == 0)
        return false
    else
        return true
    end
end

function build_empty_sol(data::RvrpInstance, id::String)
    unassigned = [req.id for req in data.requests]
    problem_id = data.id
    sol = RvrpSolution(id, problem_id, 0.0, Route[], unassigned)
    return sol
end


text:
"""
    gammaHSIC(X::Vector{T}, Y::Vector{T})

"""

function gammaHSIC(X::Vector{T}, Y::Vector{T}; randomSubSet = 100, kernelSize = -1) where {T<:Real}

	M = length(X)

	# get kernel sizes
	if kernelSize == -1
		sigx = estimateKernelSize(X, sampleSize = randomSubSet)
		sigy = estimateKernelSize(Y, sampleSize = randomSubSet)
	else
		sigx = kernelSize
		sigy = kernelSize
	end

	bone = ones(T, M, 1)
	H = eye(M) - 1/M * ones(T, M,M)

	K = rbfDotProduct(X, X, sigx)
	L = rbfDotProduct(Y, Y, sigy)

	# NOTE: these are slightly biased estimates of centred Gram matrices
	Kc = H*K*H
	Lc = H*L*H

	# NOTE: we fit Gamma to testStat*M
	testStat = 1/M * sum(sum(Kc'.*Lc))

	varHSIC = (1/6 * Kc.*Lc).^2
	varHSIC = 1/M/(M-1)* (  sum(sum(varHSIC)) - sum(diag(varHSIC))  ) # second subtracted term is bias correction
	varHSIC = 72*(M-4)*(M-5)/M/(M-1)/(M-2)/(M-3) * varHSIC

	K = K-diagm(0 => diag(K))
	L = L-diagm(0 => diag(L))

	muX = 1/M/(M-1)*bone'*(K*bone)
	muY = 1/M/(M-1)*bone'*(L*bone)

	mHSIC  = 1/M * ( 1 .+muX*muY - muX - muY )

	al = mHSIC^2 / varHSIC
	bet = varHSIC*M ./ mHSIC

	return 1 - gammacdf(al[1], bet[1], testStat)

end


text:
type DummyKernel <: Kernel
end

@testset "Kernel" begin
    residual = Vector{Float64}(0)

    dummy_kernel = DummyKernel()

    # Test the Kernel Interface
    @test_throws MethodError MOOSE.computeQpResidual(dummy_kernel, 1, 1)
end


text:
@inline function density_integral(bound_low::Real, bound_up::Real,
                                  norm::Real, slope::Real, ρ::Real)

    # density integral (eq. 9 M01)
    nb = 4π * norm * bound_low^3 / ρ
    density = nb * ( (bound_up/bound_low)^(3.0 - slope ) - 1.0 ) / ( 3.0 - slope )

    if ( 3.0 - slope_soft ) < slope < ( 3.0 + slope_soft )
        slope_var = (slope - 3.0)/slope_soft
        density2 = nb * log10(bound_up/bound_low)
        if slope_var != 0.0
            density = density * slope_var + density2 * ( 1.0 - slope_var )
        else
            density = density2
        end
    end

    return density
end

# function calculateCRNumber(CR_N, CR_S, CR_Cut, ρ;pmin::Real=10.0, pmax::Real=1.e7,
#                            mc::Real=0.0, SelectBin::Int=-1,
#                            units::GadgetPhysical=GadgetPhysical(), verbose=true)
#     # calculates engergy per bin in cgs units.

#     CR_Num = 0.0
#     cnst_c = 2.9979e10
#     ρ *= units.m_unit / units.l_unit^3  # convert rho in cgs units

#     CR_N .*= 1.e20  # compensate for io

#     Nbins = length(CR_N)
#     bin_width = log10(pmax/pmin)/Nbins
#     bins = 10.0.^collect(log10(pmin):bin_width:log10(pmax)) .* mc
#     above_cut = findall( bins .> CR_Cut/mc )

#     if length(above_cut) > 0
#         bins[above_cut] .= CR_Cut/mc
#     end

#     if SelectBin == -1
#         CR_Num = zeros(length(CR_N))
#         for i = 1:length(CR_N)
#             CR_Num[i] = density_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
#         end
#     elseif length(SelectBin) == 1
#         CR_Num = density_integral(bin[SelectBin], bin[SelectBin+1], CR_N[SelectBin], CR_S[SelectBin], ρ)
#     elseif length(SelectBin) > 1
#         CR_Num = zeros(length(SelectBin))
#         if i ∈ SelectBin
#             CR_Num[i-SelectBin[1]+1] = density_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
#         end
#     end
#     return CR_Num
# end


text:
#=
tester:
- Julia version: 1.3.1
- Author: erathorn
- Date: 2020-10-07
OPENBLAS_NUM_THREADS=5 JULIA_NUM_THREADS=5 /home/jo/Julia14/julia-1.4.2/bin/julia -O3
=#
using Distributed
addprocs(2)
@everywhere include("./src/MCPhylo.jl")
@everywhere using .MCPhylo
@everywhere using Random

#@everywhere Random.seed!(42)
#
#@everywhere using LinearAlgebra
#@everywhere BLAS.set_num_threads(10)
#include("./src/MCPhylo.jl")
#using .MCPhylo
#using Random
#Random.seed!(42)

#mt, df = make_tree_with_data("Example.nex", binary=true); # load your own nexus file
mt, df = make_tree_with_data("notebook/Dravidian.cc.phy.nex"); # load your own nexus file

mt2 = deepcopy(mt)
randomize!(mt2)
my_data = Dict{Symbol, Any}(
  :mtree => mt,
  :df => df,
  :nnodes => size(df)[3],
  :nbase => size(df)[1],
  :nsites => size(df)[2],
);

# model setup
model =  Model(
    df = Stochastic(3, (mtree, mypi) ->  PhyloDist(mtree, mypi, [1.0], [1.0], Restriction), false, false),
    mypi = Stochastic(1, () -> Dirichlet(2,1)),
    mtree = Stochastic(Node(), () -> CompoundDirichlet(1.0,1.0,0.100,1.0), true)
     )
# intial model values
inits = [ Dict{Symbol, Union{Any, Real}}(
    :mtree => mt,
    :mypi=> rand(Dirichlet(2,1)),
    :df => my_data[:df],
    :nnodes => my_data[:nnodes],
    :nbase => my_data[:nbase],
    :nsites => my_data[:nsites],
    :a => rand(),
    ),
    Dict{Symbol, Union{Any, Real}}(
        :mtree => mt2,
        :mypi=> rand(Dirichlet(2,1)),
        :df => my_data[:df],
        :nnodes => my_data[:nnodes],
        :nbase => my_data[:nbase],
        :nsites => my_data[:nsites],
        :a => rand()
        )
    ]

scheme = [PNUTS(:mtree, target=0.7, targetNNI=4),
           SliceSimplex(:mypi),
          ]

setsamplers!(model, scheme);

# do the mcmc simmulation. if trees=true the trees are stored and can later be
# flushed ot a file output.
sim = mcmc(model, my_data, inits, 500000, burnin=250000,thin=50, chains=2, trees=true)

# request more runs
sim = mcmc(sim, 2500, trees=true)

# write the output to a path specified as the second argument
to_file(sim, "example_run")


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

#=
Jiménez, J., Segura, R. and Feito, F. 2009.
(https://www.sciencedirect.com/science/article/pii/S0925772109001448?via%3Dihub)
=#
function intersection(f, s::Segment{3,T}, t::Triangle{3,T}) where {T}
  vₛ = vertices(s)
  vₜ = vertices(t)

  A = vₛ[1] - vₜ[3]
  B = vₜ[1] - vₜ[3]
  C = vₜ[2] - vₜ[3]

  W₁ = B × C
  wᵥ = A ⋅ W₁

  D = vₛ[2] - vₜ[3]
  sᵥ = D ⋅ W₁

  if wᵥ > atol(T)
    # rejection 2
    if sᵥ > atol(T)
      return @IT NoIntersection nothing f
    end

    W₂ = A × D
    tᵥ = W₂ ⋅ C

    # rejection 3
    if tᵥ < -atol(T)
      return @IT NoIntersection nothing f
    end

    uᵥ = -(W₂ ⋅ B)

    # rejection 4
    if uᵥ < -atol(T)
      return @IT NoIntersection nothing f
    end

    # rejection 5
    if wᵥ < (sᵥ + tᵥ + uᵥ)
      return @IT NoIntersection nothing f
    end
  elseif wᵥ < -atol(T)
    # rejection 2
    if sᵥ < -atol(T)
      return @IT NoIntersection nothing f
    end

    W₂ = A × D
    tᵥ = W₂ ⋅ C

    # rejection 3
    if tᵥ > atol(T)
      return @IT NoIntersection nothing f
    end

    uᵥ = -(W₂ ⋅ B)

    # rejection 4
    if uᵥ > atol(T)
      return @IT NoIntersection nothing f
    end

    # rejection 5
    if wᵥ > (sᵥ + tᵥ + uᵥ)
      return @IT NoIntersection nothing f
    end
  else
    if sᵥ > atol(T)
      W₂ = D × A
      tᵥ = W₂ ⋅ C

      # rejection 3
      if tᵥ < -atol(T)
        return @IT NoIntersection nothing f
      end

      uᵥ = -(W₂ ⋅ B)

      # rejection 4
      if uᵥ < -atol(T)
        return @IT NoIntersection nothing f
      end
      # rejection 5
      if -sᵥ < (tᵥ + uᵥ)
        return @IT NoIntersection nothing f
      end
    elseif sᵥ < -atol(T)
      W₂ = D × A
      tᵥ = W₂ ⋅ C

      # rejection 3
      if tᵥ > atol(T)
        return @IT NoIntersection nothing f
      end

      uᵥ = -(W₂ ⋅ B)

      # rejection 4
      if uᵥ > atol(T)
        return @IT NoIntersection nothing f
      end

      # rejection 5
      if -sᵥ > (tᵥ + uᵥ)
        return @IT NoIntersection nothing f
      end
    else
      # rejection 1, coplanar segment
      return @IT NoIntersection nothing f
    end
  end

  λ = clamp(wᵥ / (wᵥ - sᵥ), zero(T), one(T))

  return @IT IntersectingSegmentTriangle s(λ) f
end


text:
export j1900, j1950, j2000, j2100, jyear

"""
    j1900()

Returns the Julian Date of 1899 DEC 31 12:00:00 (1900 JAN 0.5).

[https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j1900_c.html]
(https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j1900_c.html)
"""
function j1900()
    ccall((:j1900_c, libcspice), Cdouble, ())
end

"""
    j1950()

Returns the Julian Date of 1950 JAN 01 00:00:00 (1950 JAN 1.0).

[https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j1950_c.html]
(https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j1950_c.html)
"""
function j1950()
    ccall((:j1950_c, libcspice), Cdouble, ())
end

"""
    j2000()

Returns the Julian Date of 2000 JAN 01 12:00:00 (2000 JAN 1.5).

[https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j2000_c.html]
(https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j2000_c.html)
"""
function j2000()
    ccall((:j2000_c, libcspice), Cdouble, ())
end

"""
    j2100()

Returns the Julian Date of 2100 JAN 01 12:00:00 (2100 JAN 1.5).

[https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j2100_c.html]
(https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/j2100_c.html)
"""
function j2100()
    ccall((:j2100_c, libcspice), Cdouble, ())
end

"""
    jyear()

Returns the number of seconds per Julian year.

[https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/jyear_c.html]
(https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/jyear_c.html)
"""
function jyear()
    ccall((:jyear_c, libcspice), Cdouble, ())
end


text:
using Match

include("debug_println.jl")

# An AstWalker walks over and visits all elements of the AST.
abstract AstWalker

visit(w::AstWalker, s::Symbol) = visit_symbol(w, s)
visit(w::AstWalker, i::Int) = visit_int(w, i)
visit(w::AstWalker, l::LineNumberNode) = visit_line_node(w, l)
visit(w::AstWalker, q::QuoteNode) = visit_quote_node(w, q)
visit(w::AstWalker, e::Expr) = visit_expr(w, e)
visit(w::AstWalker, n) = error("Cannot visit AST node $n of type $(typeof(n))")

visit_expr(w::AstWalker, e::Expr) = dispatch_on_expr_head(w, e)

function dispatch_on_expr_head(w::AstWalker, e::Expr)
	dprintln("Dispatching on $e")
	@match e begin
		Expr(:function,  [signature, body], _)  => visit_expr_function(w, signature, body, e)
		Expr(:block,  stmts, _)      			=> visit_expr_block(w, stmts, e)
		Expr(:call,  [name, args...], _)      	=> visit_expr_call(w, name, args, e)
      	_   									=> error("Cannot visit Expr $e")
  	end
end

# Default behavior is to just visit nodes and return the node itself.
# Override for more specific behavior.
visit_int(w::AstWalker, i::Int) = i
visit_symbol(w::AstWalker, s::Symbol) = s
visit_quote_node(w::AstWalker, q::QuoteNode) = q
visit_line_node(w::AstWalker, l::LineNumberNode) = l
function visit_expr_function(w::AstWalker, sig, body, funcex)
	visit_expr_function_sig(w, sig, funcex)
	visit_expr_function_body(w, body, funcex)
	funcex
end
visit_expr_function_sig(w::AstWalker, sig, funcex) = sig
function visit_expr_function_body(w::AstWalker, body, funcex)
	visit(w, body)
	body
end
function visit_expr_block(w::AstWalker, stmts, blockex)
	dprintln("Visiting block:\n$e")
	map(s -> visit(w, s), stmts)
	blockex
end
function visit_expr_call(w::AstWalker, name, args, callex)
	map(s -> visit(w, s), args)
	callex
end


text:
@doc raw"""
    PSFModels.AiryDisk([T=Float64]; fwhm, x=0, y=0, amp=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.AiryDisk([T=Float64]; fwhm, pos, amp=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.AiryDisk([T=Float64]; fwhm, r, theta, amp=1, maxsize=3, extent=maxsize .* fwhm)

An unnormalized Airy disk. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. By default the model is placed at the origin. The position can also be given as a polar coordinate using `r`/`ρ` and `theta`/`θ`, optionally centered around `origin`.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). For efficient calculations, we recommend using [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl). Here, `maxsize` is a multiple of the fwhm, and can be given as a scalar or as a tuple for each axis. The `extent` defines the bounding box for the model and is used for the default rendering size.

# Functional form

The Airy disk is a distribution over the radius `r` (the square-Euclidean distance)

```
f(x | x̂, FWHM) = [ 2J₁(q) / q ]^2
```
where `J₁` is the first-order Bessel function of the first kind and
```
q ≈ π * r / (0.973 * FWHM)
```
"""
struct AiryDisk{T,FT,VT<:AbstractVector,IT<:Tuple} <: PSFModel{T}
    fwhm::FT
    pos::VT
    amp::T
    indices::IT

    function AiryDisk(::Type{T}, fwhm::FT, pos::VT, amp::T, indices::IT) where {T,VT<:AbstractVector,FT,IT<:Tuple}
        new{T,FT,VT,IT}(fwhm, pos, amp, indices)
    end
end

## constructors
# default type is Float64
AiryDisk(; kwargs...) = AiryDisk(Float64; kwargs...)
function AiryDisk(T; fwhm, amp=one(T), maxsize=3, extent = maxsize .* fwhm, position...)
    # get the position from keyword distpatch
    pos = _position(; position...)
    return AiryDisk(T, fwhm, pos, convert(T, amp), indices_from_extent(pos, extent))
end

Base.size(a::AiryDisk) = map(length, a.indices)
Base.axes(a::AiryDisk) = a.indices

const rz = 3.8317059702075125 / π

function (a::AiryDisk{T})(point::AbstractVector) where T
    radius = a.fwhm * 1.18677
    Δ = euclidean(point, a.pos)
    r = Δ / (radius / rz)
    val = ifelse(iszero(r), a.amp, a.amp * (2 * besselj1(π * r) / (π * r))^2)
    return convert(T, val)
end

function (a::AiryDisk{T,<:Union{AbstractVector,Tuple}})(point::AbstractVector) where T
    weights = SA[(rz / (first(a.fwhm) * 1.18677))^2, (rz / (last(a.fwhm) * 1.18677))^2]
    r = weuclidean(point, a.pos, weights)
    val = ifelse(iszero(r), a.amp, a.amp * (2 * besselj1(π * r) / (π * r))^2)
    return convert(T, val)
end


text:
using Statistics
using Random
using Distributions
using JSON
using POMDPs
using MCTS
using NearestNeighbors
using GeometryTypes
using JLD2

using SCoBA
using SCoBA.SCoBADomains
using SCoBA.SCoBADomains: LatLonCoords
using SCoBA.SCoBASolver

# Hardcode the depot locations to get good coverage
const DEPOT1 = LatLonCoords((lat = 37.774524, lon = -122.473656))
const DEPOT2 = LatLonCoords((lat = 37.751751, lon = -122.410654))
const DEPOT3 = LatLonCoords((lat = 37.718779, lon = -122.462401))
const DEPOT4 = LatLonCoords((lat = 37.789290, lon = -122.426797))
const DEPOT5 = LatLonCoords((lat = 37.739611, lon = -122.492203))
const DEPOTLOCS = [DEPOT1, DEPOT2, DEPOT3, DEPOT4, DEPOT5]
const TRAVELTIME_EST = "./param_files/sf_halton_tt_estimates_scoba.jld2"
const PARAMS_FN = "./param_files/sf_bb_params.toml"
const rng = MersenneTwister(1345)
const mcts_params = (trials=100, explore=0.1, depth=20)

function main(n_depots::Int64, n_drones::Int64, n_requests::Int64, outfile::String,
              n_trials::Int64=10, method::String="scoba")

    # To store the computation times
    comp_times = Float64[]

    # Load up travel time estimates
    @load TRAVELTIME_EST  scoba_halton_tree travel_time_estimates

    city_params = SCoBADomains.parse_city_params(PARAMS_FN)
    lat_dist = Uniform(city_params.lat_start, city_params.lat_end)
    lon_dist = Uniform(city_params.lon_start, city_params.lon_end)

    # Iterate over trials
    for t = 1:n_trials+1

        # Define set of drones
        drones_per_depot = div(n_drones, n_depots)
        drone_ordering = String[]
        drone_set = Dict{String,Drone}()

        # Distribute drones equally across depots
        n = 1
        for nd = 1:n_depots
            for dd = 1:drones_per_depot

                drone_nm = string("dn", n)
                n += 1
                new_drone = Drone(depot_number=nd, depot_loc=DEPOTLOCS[nd])

                push!(drone_ordering, drone_nm)
                drone_set[drone_nm] = new_drone
            end
        end # for in n_depots

        # Boilerplate code since drones don't have any intrinsic properties here
        drone_prop_set = Dict{String,DroneProperties}()
        for drone_nm in drone_ordering
            drone_prop_set[drone_nm] = DroneProperties()
        end

        drone_server = RoutingAllocation(agent_set=drone_set, agent_prop_set=drone_prop_set,
                                         agent_ordering=drone_ordering,
                                         max_tasks_to_consider=20,
                                         conflict_threshold=20)

        sim = setup_routing_sim(PARAMS_FN, scoba_halton_tree, travel_time_estimates,
                                num_init_requests=n_requests,
                                new_request_prob=1.0, rng=rng)

        if method == "scoba"

            update_time_windows!(sim, drone_server)
            t = @elapsed scoba_routing!(drone_server, sim, rng)

            push!(comp_times, t)

        elseif method == "hungarian"

            update_time_windows!(sim, drone_server)
            t = @elapsed expected_hungarian!(drone_server, sim, rng)

            push!(comp_times, t)

        elseif method == "mcts"

            mdp = RoutingMCTSMDP(sim, drone_server, 100, Int64[], Tuple{Float64,Float64}[])
            nearestpol = NearestPkgPolicy(mdp)

            mcts = MCTSSolver(n_iterations=mcts_params.trials, depth=mcts_params.depth, rng=rng,
                              exploration_constant=mcts_params.explore,
                              reuse_tree=false,
                              estimate_value = RolloutEstimator(nearestpol))
            mcts_policy = POMDPs.solve(mcts, mdp)

            update_time_windows!(sim, drone_server)
            update_routing_mcts_fullstate!(mdp, rng)

            curr_action = Int64[]

            sum_t = 0.0

            for idx = 1:n_drones

                drone_state = get_current_routing_mcts_state(mdp, idx)
                t = @elapsed drone_action = POMDPs.action(mcts_policy, drone_state)
                mdp.drone_pkg_assignment[idx] = drone_action
                push!(curr_action, drone_action)
                sum_t += t
            end

            push!(comp_times, sum_t)
        end
    end # t = 1:n_Trials + 1

    meanval =  mean(comp_times[2:end]) # Ignore first time
    @show meanval

    timing_results = Dict("mean"=>meanval,
                          "sem"=>std(comp_times[2:end])/sqrt(n_trials),
                          "median"=>median(comp_times[2:end]))

    open(outfile, "w") do f
        JSON.print(f, timing_results, 2)
    end
end # main

# Example call to main - for 5 depots, 15 drones, 50 requests, and 50 trials
# main(5, 15, 50, "routing_time_trial50_dr15dep5req50.json", 50)


text:
using Test
using TestPackage

@test plusTwo(3) == 5


text:
mutable struct Filler
    x::Integer
    y::Integer
    w::Integer
    h::Integer
end

function Base.Dict(filler::Filler)
    return Dict{String, Any}(
        "__name" => "rect",

        "x" => filler.x,
        "y" => filler.y,

        "w" => filler.w,
        "h" => filler.h
    )
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
export rl_unplugged_bsuite_dataset

using TFRecord

struct BSuiteRLTransition <: RLTransition
    state
    action
    reward
    terminal
    next_state
end

function BSuiteRLTransition(example::TFRecord.Example, game::String)
    f = example.features.feature
    if game == "catch"
        s = reshape(f["observation"].float_list.value[1:end÷2], 10, 5)
        s′ = reshape(f["observation"].float_list.value[end÷2+1:end], 10, 5)
    else
        s = f["observation"].float_list.value[1:end÷2]
        s′ = f["observation"].float_list.value[end÷2+1:end]
    end
    a = f["action"].int64_list.value[1]
    r = f["reward"].float_list.value[]
    t = f["step_type"].int64_list.value[1] == 2
    BSuiteRLTransition(s, a, r, t, s′)
end
"""
    rl_unplugged_bsuite_dataset(game, shards, type; <keyword arguments>)

Return a `RingBuffer`(@ref) of [`BSuiteRLTransition`](@ref) batches which supports
multi threaded loading. Check out `bsuite_params()` for more info on arguments.

# Arguments

- `game::String`: name of the dataset. available datasets: `cartpole`, `mountain_car` and  `catch`.
- `shards::Vector{Int}`: the shards that are to be loaded.
- `type::String`: can be `full`, `full_train` and `full_valid`.
- `is_shuffle::Bool`
- `stochasticity::Float32`: represents the stochasticity of the dataset. can be
in the range: `0.0:0.1:0.5`.
- `shuffle_buffer_size::Int=10_000`: size of the shuffle_buffer used in loading AtariRLTransitions.
- `tf_reader_bufsize::Int=10_000`: the size of the buffer `bufsize` that is used internally
in `TFRecord.read`.
- `tf_reader_sz::Int=10_000`: the size of the `Channel`, `channel_size` that is returned by
`TFRecord.read`.
- `batch_size::Int=256`: The number of samples within the batches that are returned by the `Channel`.
- `n_preallocations::Int=nthreads()*12`: the size of the buffer in the `Channel` that is returned.

!!! note

    To enable reading records from multiple files concurrently, remember to set the number of
    threads correctly (See [JULIA_NUM_THREADS](https://docs.julialang.org/en/v1/manual/environment-variables/#JULIA_NUM_THREADS)).
"""
function rl_unplugged_bsuite_dataset(
    game::String,
    shards::Vector{Int},
    type::String;
    is_shuffle::Bool=true,
    stochasticity::Float64=0.0,
    shuffle_buffer_size::Int=10_000,
    tf_reader_bufsize::Int=10_000,
    tf_reader_sz::Int=10_000,
    batch_size::Int=256,
    n_preallocations::Int=nthreads()*12
)
    n = nthreads()

    repo = "rl-unplugged-bsuite"

    folders= [
        @datadep_str "$repo-$game-$stochasticity-$shard-$type"
        for shard in shards
    ]

    ch_files = Channel{String}(length(folders)) do ch
        for folder in cycle(folders)
            file = folder * "/$(readdir(folder)[1])"
            put!(ch, file)
        end
    end

    if is_shuffle
        files = buffered_shuffle(ch_files, length(folders))
    else
        files = ch_files
    end

    ch_src = Channel{BSuiteRLTransition}(n * tf_reader_sz) do ch
        for fs in partition(files, n)
            Threads.foreach(
                TFRecord.read(
                    fs;
                    compression=:gzip,
                    bufsize=tf_reader_bufsize,
                    channel_size=tf_reader_sz,
                );
                schedule=Threads.StaticSchedule()
            ) do x
                put!(ch, BSuiteRLTransition(x, game))
            end
        end
    end

    if is_shuffle
        transitions = buffered_shuffle(
        ch_src,
        shuffle_buffer_size
        )
    else
        transitions = ch_src
    end

    taskref = Ref{Task}()

    ob_size = game=="mountain_car" ? 3 : 6

    if game == "catch"
        obs_template = Array{Float32, 3}(undef, 10, 5, batch_size)
    else
        obs_template = Array{Float32, 2}(undef, ob_size, batch_size)
    end

    buffer = BSuiteRLTransition(
        copy(obs_template),
        Array{Int, 1}(undef, batch_size),
        Array{Float32, 1}(undef, batch_size),
        Array{Bool, 1}(undef, batch_size),
        copy(obs_template),
    )

    res = RingBuffer(buffer;taskref=taskref, sz=n_preallocations) do buff
        Threads.@threads for i in 1:batch_size
            batch!(buff, take!(transitions), i)
        end
    end

    bind(ch_src, taskref[])
    bind(ch_files, taskref[])
    res
end


text:
export process_raw, process_raw!, process_fft, compute_fft
export onebit!, onebit, remove_response!, remove_response, amplitude!, amplitude
export clip, clip!, clamp, clamp!, mute!, mute
import Base: clamp, clamp!
"""
    process_raw!(S,fs)

Pre-process raw seismic data.

- Removes mean from each channel in `S`.
- Detrends each channel in `S`.
- Downsamples data to sampling rate `fs`
- Phase-shifts data to begin at 00:00:00.0

# Arguments
- `S::SeisData`: SeisData structure.
- `fs::Float64`: Sampling rate to downsample `S`.
"""
function process_raw!(S::SeisData, fs::Float64; ϕshift::Bool=true)
    merge!(S)
    ungap!(S)

    for ii = 1:S.n
        SeisIO.detrend!(S[ii])         # remove mean & trend from channel
        SeisNoise.taper!(S[ii].x,S[ii].fs)         # taper channel ends
        if fs ∉ S.fs
            lowpass!(S[ii].x,fs/2,S[ii].fs,zerophase=true)    # lowpass filter before downsampling
        end
        resample!(S,chans=ii,fs=fs) # downsample to lower fs
        SeisNoise.taper!(S[ii].x,S[ii].fs)
        phase_shift!(S[ii], ϕshift=ϕshift) # timing offset from sampling period
    end
    return nothing
end
process_raw(S::SeisData, fs::Float64;
           ϕshift::Bool=true) = (U = deepcopy(S);
           process_raw!(U,fs, ϕshift=ϕshift); return U)

"""

  compute_fft(R)

Computes windowed rfft of ambient noise data. Returns FFTData structure.

# Arguments
- `R::RawData`: RawData structure
"""
function compute_fft(R::RawData)
    FFT = rfft(R.x,1)
    return FFTData(R.name, R.id,R.loc, R.fs, R.gain, R.freqmin, R.freqmax,
                   R.cc_len, R.cc_step, R.whitened, R.time_norm, R.resp,
                   R.misc, R.notes, R.t, FFT)
end

"""

  mute(A,factor=factor)

Set high amplitudes in array `A` to zero.
Uses median of envelope of `A` to find outliers.
"""
function mute!(A::AbstractArray;factor::Int=3)
    T = eltype(A)
    envelope = abs.(hilbert(A))
    levels = mean(envelope,dims=1)
    level = factor * median(levels)
    A[envelope .> level] .= T(0)
    return nothing
end
mute(A::AbstractArray;factor::Int=3) = (U = deepcopy(A); mute!(U,factor=factor);
     return U)
mute!(R::RawData;factor::Int=3) = mute!(R.x,factor=factor)
mute(R::RawData;factor::Int=3) = (U = deepcopy(R); mute!(U,factor=factor);
     return U)

"""
  clip(A,factor)

Truncate array A at `factor` times the root mean square of each column.

#Arguments
- `A::AbstractArray`: N-d time series array
- `factor::Real`:
- `f::Function`: Input statistical function (e.g. rms, var, std, mad)
- `dims`: Dimension of `A` to apply clipping (defaults to 1)
"""
function clip!(A::AbstractArray{T,N}, factor::Real; f::Function=rms,dims=1) where {T,N}
    if N == 1
        high = f(A) .* factor
        clamp!(@view(A[:]),-high,high)
    else
        for ii = 1:size(A,2)
            high = f(@view(A[:,ii])) * factor
            clamp!(@view(A[:,ii]),-high,high)
        end
    end
    return nothing
end
clip(A::AbstractArray, factor::Real; f::Function=rms, dims=1) = (U = deepcopy(A);
     clip!(U,factor,f=f,dims=dims);return U)
clip!(R::RawData,factor::Real;f::Function=rms,dims=1) = clip!(R.x,factor,f=f,dims=dims)
clip(R::RawData,factor::Real;f::Function=rms,dims=1) = (U = deepcopy(R);
     clip!(U.x,factor,f=f,dims=dims); return U)
clamp!(R::RawData,val::Real) = clamp!(R.x,-abs(val),abs(val))
clamp!(R::RawData,lo::Real,hi::Real) = clamp!(R.x,lo,hi)
clamp(R::RawData,val::Real) = (U = deepcopy(R); clamp!(U,val);return U)
clamp(R::RawData,lo::Real,hi::Real) = (U = deepcopy(R); clamp!(U,lo,hi);return U)

"""
  amplitude!(R)

Filter raw data based on amplitude.
"""
function amplitude!(R::RawData; max_std::Float64=10.)
    # remove nonzero columns
    zeroind = nonzero(R.x)
    if length(zeroind) == 0
        R = nothing
    elseif size(R.x,2) != length(zeroind)
        R.x = R.x[:,zeroind]
        starts = starts[zeroind]
    end

    # amplitude threshold indices
    stdind = std_threshold(R.x,max_std)
    if length(stdind) == 0
        R = nothing
    elseif size(R.x,2) != length(stdind)
        R.x = R.x[:,stdind]
        starts = starts[stdind]
        ends = ends[stdind]
    end
    return nothing
end
amplitude(R::RawData; max_std::Float64=10.) = (U = deepcopy(R);amplitude!(U,
          max_std=max_std); return U)

"""

  onebit!(R)

One-bit amplitude modification of RawData `R`.
"""
function onebit!(R::RawData)
    R.x .= sign.(R.x)
    return nothing
end
onebit(R::RawData) = (U = deepcopy(R); sign!(U);return U)

"""
  remove_response!(S, stationXML, freqmin, freqmax)

Loads instrument response from stationXML and removes response from `S`.

# Arguments
- `S::SeisData`: SeisData structure.
- `stationXML::String`: Path to stationXML file, e.g. "/path/to/file.xml"
- `freqmin::Float64`: minimum frequency for pre-filtering.
- `freqmax::Float64`: maximum frequency for pre-filtering.
- `np::Int`: number of poles in pre-filter.
- `t_max::Float64`: Length of taper in seconds.
- `wl::Float32`: Water level for instrument response.
"""
function remove_response!(S::SeisData, stationXML::String, freqmin::Float64,
                          freqmax::Float64;np::Int=2, t_max::Float64=20.,
                          wl::Float32=eps(Float32))

    # read response file
    if !isfile(stationXML)
        error("$stationXML does not exist. Instrument response not removed.")
    end

    s,e = string.(start_end(S[1]))
    R = read_sxml(stationXML,s=s,t=e)
    # loop through responses
    Rid = R.id
    Sid = S.id

    # remove trend, taper and prefilter
    SeisIO.detrend!(S)
    SeisIO.taper!(S,t_max=t_max)
    filtfilt!(S,fl=freqmin,fh=freqmax,np=np,rt="Bandpass")

    # remove instrument response for each channel in S
    @inbounds for ii = 1:S.n
        id = S[ii].id
        ind = findfirst(x -> x == id,Rid)
        LOC = R[ind].loc
        GAIN = R[ind].gain
        RESP = R[ind].resp
        UNITS = R[ind].units
        setindex!(S.loc,LOC,ii)
        setindex!(S.gain,GAIN,ii)
        setindex!(S.units,UNITS,ii)
        translate_resp!(S,RESP,chans=ii,wl=wl)
    end
    return nothing
end
remove_response(S::SeisData, stationXML::String, freqmin::Float64,
                          freqmax::Float64;np::Int=2, t_max::Float64=20.,
                          wl::Float32=eps(Float32)) = (U = deepcopy(S);
            remove_response!(U,stationXML,freqmin,freqmax,np=np,t_max=t_max,
                             wl=wl); return U)

"""
nonzero(A)

Find indices of all nonzero columns in array `A`.
"""
function nonzero(A::AbstractArray)
    Nrows, Ncols = size(A)
    ind = Int64[]
    sizehint!(ind,Ncols)
    for ii = 1:Ncols
        for jj = 1:Nrows
            if !iszero(A[jj,ii])
                append!(ind,ii)
                break
            end
        end
    end
    return ind
end

"""
    compute_fft(S, freqmin, freqmax, fs, cc_step, cc_len;
                time_norm=false, to_whiten=false, max_std=5.)
Computes windowed fft of ambient noise data.
# Arguments
- `S::SeisData`: SeisData structure.
- `freqmin::Float64`: minimum frequency for filtering/whitening.
- `freqmax::Float64`: maximum frequency for filtering/whitening.
- `fs::Float64`: Sampling rate to downsample `S`.
- `cc_step::Int`: time, in seconds, between successive cross-correlation windows.
- `cc_len::Int`: length of noise data window, in seconds, to cross-correlate.
- `time_norm::Union{Bool,String}`: time domain normalization to perform.
- `to_whiten::Bool`: Apply whitening in frequency domain.
- `max_std::Float64=5.`: Number of standard deviations above mean to reject windowed data.
"""
function compute_fft(S::SeisData,freqmin::Float64,freqmax::Float64,fs::Float64,
                     cc_step::Int, cc_len::Int;
                     time_norm::Union{Bool,String}=false,
                     to_whiten::Bool=false,
                     max_std::Float64=5.,
                     ϕshift::Bool=true)

    # sync!(S,s=starttime,t=endtime)
    merge!(S)
    ungap!(S)
    process_raw!(S,fs,ϕshift=ϕshift)  # demean, detrend, taper, lowpass, downsample

    # subset by time
    starttime, endtime = u2d.(nearest_start_end(S[1],cc_len, cc_step))
    # check if waveform length is < cc_len
    if Int(floor((endtime - starttime).value / 1000)) < cc_len
        return nothing
    end

    sync!(S,s=starttime,t=endtime)
    A, starts, ends = slide(S[1], cc_len, cc_step)

    # remove nonzero columns
    zeroind = nonzero(A)
    if length(zeroind) == 0
        return nothing
    elseif size(A,2) != length(zeroind)
        A = A[:,zeroind]
        starts = starts[zeroind]
        ends = ends[zeroind]
    end

    # amplitude threshold indices
    stdind = std_threshold(A,max_std)
    if length(stdind) == 0
        return nothing
    elseif size(A,2) != length(stdind)
        A = A[:,stdind]
        starts = starts[stdind]
        ends = ends[stdind]
    end

    FFT = process_fft(A, freqmin, freqmax, fs, time_norm=time_norm,
                      to_whiten=to_whiten)
    return F = FFTData(S[1].id, Dates.format(u2d(starts[1]),"Y-mm-dd"),
                       S[1].loc, S[1].fs, S[1].gain, freqmin, freqmax,
                       cc_len, cc_step, to_whiten, time_norm, S[1].resp,
                       S[1].misc, S[1].notes, starts, FFT)
end


text:
#######
# det #
#######

function LinearAlgebra.det(x::TrackedArray{V,D}) where {V,D}
    tp = tape(x)
    x_value = value(x)
    det_x_value = det(x_value)
    out = track(det_x_value, D, tp)
    record!(tp, SpecialInstruction, det, x, out)
    return out
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(det)})
    output = instruction.output
    if output != 0
        input = instruction.input
        input_deriv = deriv(input)
        inv_input_value = inv(value(input))
        k = deriv(output) * value(output)
        for i in 1:size(input_deriv, 1), j in 1:size(input_deriv, 2)
            input_deriv[i, j] += k * inv_input_value[j, i]
        end
    end
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(det)})
    output, input = instruction.output, instruction.input
    value!(output, det(value(input)))
    return nothing
end

#######
# inv #
#######

function LinearAlgebra.inv(x::TrackedArray{V,D}) where {V,D}
    tp = tape(x)
    out_value = inv(value(x))
    out = track(out_value, D, tp)
    cache = (similar(out_value, D), similar(out_value, D))
    record!(tp, SpecialInstruction, inv, x, out, cache)
    return out
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(inv)})
    output = instruction.output
    output_value, output_deriv = value(output), deriv(output)
    output_tmp1, output_tmp2 = instruction.cache
    mul!(output_tmp1, output_deriv, adjoint(output_value))
    mul!(output_tmp2, adjoint(output_value), output_tmp1)
    decrement_deriv!(instruction.input, output_tmp2)
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(inv)})
    input = instruction.input
    value!(instruction.output, inv(value(input)))
    return nothing
end


text:
using Gadfly

set_default_plot_size(6inch, 8inch)

plot(x=rand(40), y=rand(40),
     xgroup=vcat(fill("A", 20), fill("B", 20)),
     ygroup=vcat(fill("X", 10), fill("Y", 10), fill("X", 10), fill("Y", 10)),
     Geom.subplot_grid(Geom.smooth))


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("mcmc_weighting.jl")
include("proposaldist.jl")
include("mcmc_sampleid.jl")
include("mcmc_algorithm.jl")
include("mcmc_noop_tuner.jl")
include("mcmc_stats.jl")
include("mcmc_convergence.jl")
include("chain_pool_init.jl")
include("multi_cycle_burnin.jl")
include("mcmc_sample.jl")
include("mh/mh.jl")
include("ahmc/ahmc.jl")


text:
"""
    setdiff(d1,d2)

Find set difference between two MuttsDict{K,V}.  Returns an unordered Vector{K,V}
enumerating the (key,value) pairs that are in d1 but not d2.

Internal nodes that are common to both d1 and d2 are skipped.  If d2 is obtained
by branching d1 and performing δ insert/delete operations, the cost of setdiff
is ϴ(δ n^(1/7)).  (This also applies if d1 derives from a branch of d2.)
"""
function Base.setdiff(d1::MuttsDict{K,V}, d2::MuttsDict{K,V}) where {K,V}
    # Recursively traverse d2. At each node of d2, check whether the current hash path
    # reaches the same node in d1; if so, skip the subtree. When get to a leaf that is
    # in d1 but not d2, iterate the entries and compare to d2.
    diffs = Vector{Tuple{K,V}}()
    height = length(d1.config.inode_capacity)
    pathmask = UInt64(0)
    pathhash = UInt64(0)
    _collect_diffs(diffs, d1.root, d2, height, pathmask, pathhash)
    return diffs
end

# Collect diffs for an inode subtree in d1
function _collect_diffs(diffs::Vector{Tuple{K,V}},
    inode::INode{K,V},
    d2::MuttsDict{K,V},
    height::Int,
    pathmask::UInt64,
    pathhash::UInt64
) where {K,V}
    if _has_node(d2, objectid(inode), pathmask, pathhash)
        # This path reaches the same node object in both dictionaries, so there
        # will be no diffs in this subtree.
        return
    end

    if isa(inode.entries, Vector{INode{K,V}})
        capacity = length(inode.entries)
        mask = reinterpret(UInt64, capacity-1)
        shift = reinterpret(UInt64, height+1)*UInt64(8)
        subtree_pathmask = pathmask | (mask << shift)
        for i=1:capacity
            subtree_pathhash = pathhash | ((i-1) << shift)
            child_inode = @inbounds inode.entries[i]
            _collect_diffs(diffs, child_inode, d2, height-1, subtree_pathmask,
                subtree_pathhash)
        end
    else
        # Branch specialized for isa(inode.entries,Vector{Leaf{K,V}})
        capacity = length(inode.entries)
        mask = reinterpret(UInt64, capacity-1)
        shift = reinterpret(UInt64,height+1)*UInt64(8)
        subtree_pathmask = pathmask | (mask << shift)
        for i=1:capacity
            subtree_pathhash = pathhash | ((i-1) << shift)
            leaf = @inbounds inode.entries[i]
            _collect_diffs(diffs, leaf, d2, height-1, subtree_pathmask, subtree_pathhash)
        end
    end
end

# Collect diffs for a leaf in d1
function _collect_diffs(diffs::Vector{Tuple{K,V}},
    leaf::Leaf{K,V},
    d2::MuttsDict{K,V},
    height::Int,
    pathmask::UInt64,
    pathhash::UInt64
) where {K,V}
    if _has_node(d2, objectid(leaf), pathmask, pathhash)
        # This path reaches the same node object in both dictionaries, so there
        # will be no diffs in this subtree.
        return
    end

    # Collect diffs from individual entries in this leaf
    num_slots = length(leaf.entries)
    for i=1:num_slots
        slot = reinterpret(UInt64,i)
        if _slot_occupied(leaf, slot, num_slots)
            (ki,vi) = @inbounds leaf.entries[i]
            h = hash(ki)
            if (h & pathmask) == pathhash
                if !in((ki,vi),d2)
                    push!(diffs, (ki,vi))
                end
            end
        end
    end
end

# Returns true if dict contains a node where objectid(node)==id, reachable along
# the specified path.
function _has_node(dict::MuttsDict{K,V}, id::UInt64, pathmask, pathhash) where {K,V}
    node = dict.root
    (objectid(node) == id) && return true
    height = length(dict.config.inode_capacity)
    while isa(node, INode{K,V})
        entries = node.entries
        mask = reinterpret(UInt64, length(entries)-1)
        shift = reinterpret(UInt64, height+1)*UInt64(8)
        slot = 1 + reinterpret(Int64, ((pathhash >> shift) & mask))
        node = @inbounds entries[slot]
        if objectid(node) == id
            return true
        end
        height -= 1
    end

    return false
end

function Base.:(==)(dict1::MuttsDict{K,V}, dict2::MuttsDict{K,V}) where {K,V}
    (dict1 === dict2) && return true
    (length(dict1) != length(dict2)) && return false
    if length(dict1) < 20
        # Small dictionary - iterate & lookup, avoid overhead of allocation from
        # setdiff()
        for (k,v) in dict1
            !((k,v) in dict2) && return false
        end
        return true
    else
        return isempty(setdiff(dict1,dict2))
    end
end


text:
mutable struct PRW
    w # weight
end

function PRW(weight::Int, scale=0.01)
    w = weight
    return PRW(w)
end

function (dataset::PRW)(x)
    return l.w[:, x]
end

function load_queries(l::PRW)
    queries = []
    # TODO replicate from seqnet code
    return queries
end

function load_split_img_names(l::PRW)
    images = []
    # TODO replicate from seqnet code
    return images
end

function load_annotations(l::PRW)
    annotations = []
    # TODO replicate from seqnet code
    return annotations
end

max_stars_repo_name: ucatbas/seqnet-julia

max_issues_repo_name: ucatbas/seqnet-julia

max_forks_repo_name: ucatbas/seqnet-julia


text:
export filterFrequencies

"""
  filterFrequencies(f; kargs...) => Vector{Int64}

Supported keyword arguments:
* SNRThresh
* minFreq
* maxFreq
* recChannels
* sortBySNR
* numUsedFreqs
* stepsize
* maxMixingOrder
* sortByMixFactors
"""
function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
               maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f),
               sortBySNR=false, numUsedFreqs=-1, stepsize=1, maxMixingOrder=-1,
               sortByMixFactors=false, numPeriodAverages=1, numPeriodGrouping=1)

  nFreq = rxNumFrequencies(f, numPeriodGrouping)
  nReceivers = rxNumChannels(f)
  nPeriods = 1 #acqNumPeriodsPerFrame(f)

  minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
  maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1

  freqMask = zeros(Bool,nFreq,nReceivers,nPeriods)

  freqMask[:,recChannels,:] .= true

  if measIsFrequencySelection(f)
    freqMask[:,recChannels,:] .= false
    idx = measFrequencySelection(f)
    freqMask[idx,recChannels,:] .= true
  else
    freqMask[:,recChannels,:] .= true
  end

  if minFreq > 0
    freqMask[1:(minIdx),:,:] .= false
  end

  if maxIdx < nFreq
    freqMask[(maxIdx):end,:,:] .= false
  end

  if maxMixingOrder > 0
    if numPeriodGrouping == 1
      mf = mixingFactors(f)
      for l=1:size(mf,1)
        if mf[l,4] > maxMixingOrder || mf[l,4] > maxMixingOrder
          freqMask[(l-1)+1,recChannels] .= false
        end
      end
    else # This is a hack until we have a general solution

      numPatches = div(acqNumPeriodsPerFrame(f), numPeriodAverages)
      fBands = (collect(0:(rxNumFrequencies(f)-1)).-1) .* numPatches
      freqMaskMO = zeros(Bool,nFreq,nReceivers,nPeriods)

      for i=1:length(fBands)
        for y=-maxMixingOrder:maxMixingOrder
          index = fBands[i]+y
          if 1 <= index <= size(freqMask,1)
            freqMaskMO[index,recChannels] .= true
          end
        end
      end
      freqMask .&= freqMaskMO
    end
  end

  SNRAll = nothing

  if SNRThresh > 0 || numUsedFreqs > 0 || sortBySNR
    SNR = zeros(nFreq, nReceivers)
    idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq

    SNRAll = calibSNR(f)
    if SNRAll != nothing
      SNR[idx,:] = SNRAll[:,:,1]
    end
  end

  if SNRThresh > 0 && numUsedFreqs > 0
    error("It is not possible to use SNRThresh and SNRFactorUsedFreq similtaneously")
  end

  if SNRThresh > 0 && SNRAll != nothing
    freqMask[ findall(vec(SNR) .< SNRThresh) ] .= false
  end

  if numUsedFreqs > 0 && SNRAll != nothing
    numFreqsAlreadyFalse = sum(!freqMask)
    numFreqsFalse = round(Int,length(freqMask)* (1-numUsedFreqs))
    S = sortperm(vec(SNR))

    l = 1
    j = 1
    while j<  (numFreqsFalse-numFreqsAlreadyFalse)
      if freqMask[S[l]] == true
        freqMask[S[l]] = false
        j += 1
      end
      l += 1
    end

  end

  if stepsize > 1
    freqStepsizeMask = zeros(Bool,nFreq, nReceivers, nPatches)
    freqStepsizeMask[1:stepsize:nFreq,:,:] = freqMask[1:stepsize:nFreq,:,:]
    freqMask = freqStepsizeMask
  end

  freq = findall( vec(freqMask) )

  if sortBySNR && !sortByMixFactors && SNRAll != nothing
    SNR = vec(SNR[1:stepsize:nFreq,:,:])

    freq = freq[reverse(sortperm(SNR[freq]),dims=1)]
  end

  if !sortBySNR && sortByMixFactors
    mfNorm = zeros(nFreq,nReceivers,nPeriods)
    mf = mixingFactors(f)
    for k=1:nFreq
      mfNorm[k,:,:] = norm(mf[k,1:3])
    end

    freq = freq[sortperm(mfNorm[freq])]
  end

  freq
end

function rowsToSubsampledRows(f::MPIFile, rows)
  if measIsFrequencySelection(f)
    # In this case we need to convert indices
    tmp = zeros(Int64, rxNumFrequencies(f), rxNumChannels(f) )
    idxAvailable = measFrequencySelection(f)
    for d=1:rxNumChannels(f)
      tmp[idxAvailable, d] = (1:length(idxAvailable)) .+ (d-1)*length(idxAvailable)
    end
    rows_ = vec(tmp)[rows]
    if findfirst(x -> x == 0, rows_) != nothing
      @error "Indices applied to systemMatrix are not available in the file"
    end
  else
    rows_ = rows
  end
  return rows_
end


text:
"""
Symbolic simplifications of +, -, *, /, ^

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

If possible, the operation is performed. Special care about adding zero and multiplying with one.
"""
module Simplify

export add, sub, mult, divide, power, simplify

using Base.Meta: isexpr
using Unitful

isUnit(x) = typeof(x) <: Unitful.FreeUnits || isexpr(x, :macrocall)

function add(x, y)
    if typeof(x) in [Float64, Int64] && typeof(y) in [Float64, Int64]
        x + y
    elseif x == 0
        y
    elseif y == 0
        x
    else
        :($x + $y)
    end
end

function sub(x, y)
    if typeof(x) in [Float64, Int64] && typeof(y) in [Float64, Int64]
        x - y
    elseif x == 0 && isexpr(y, :call) && length(y.args) == 2 && y.args[1] == :-
        y.args[2]
    elseif x == 0
        :(- $y)
    elseif y == 0
        x
    else
        :($x - $y)
    end
end

function mult(x, y)
    if typeof(x) in [Float64, Int64] && typeof(y) in [Float64, Int64]
        x * y
    elseif (x == 0 && ! isUnit(y)) || (y == 0 && ! isUnit(x))  # Does not remove unit even if 0
        0
    elseif x == 1 && ! isUnit(y)
        y
    elseif y == 1 && ! isUnit(x)
        x
    else
        :($x * $y)
    end
end

function divide(x, y)
    if typeof(x) in [Float64, Int64] && typeof(y) in [Float64, Int64]
        x / y
    elseif x === 0
        0
    elseif y === 1
        x
    elseif y === -1
        sub(0, x)
    else
        :($x / $y)
    end
end

function power(x, y)
    if typeof(x) in [Float64, Int64] && typeof(y) in [Float64, Int64]
        x ^ y
    elseif y == 1
        x
    else
        :($x ^ $y)
    end
end

simplify(ex) = ex

function simplify(ex::Expr)
    args = [simplify(arg) for arg in ex.args]
    if ex.head == :call && ex.args[1] in [:+, :-, :*] && all([typeof(a) in [Float64, Int64] for a in args[2:end]])
        if ex.args[1] in [:+]
            return sum(args[2:end])
        elseif ex.args[1] in [:-]
            return args[2] - sum(args[3:end])
        else
            return prod(args[2:end])
        end
    elseif ex.head == :call && ex.args[1] == :* && any([a == 0 for a in args[2:end]])
        # x*0*y = 0
        return 0
    elseif ex.head == :call && ex.args[1] == :^ && length(args) == 3 && args[3] == 1
        # x^1 = x
        return args[2]
    elseif ex.head == :call && ex.args[1] == :^ && length(args) == 3 && args[3] == 0
        # x^0 = 1
        return 1
    else
        return Expr(ex.head, args...)
    end
end

function testSimplify()
    @show simplify(:(1+3))
    @show simplify(:(3-1))
    @show simplify(:(4*3))
    @show simplify(:(4*3+x))
    @show simplify(:(1+3+x))
    @show simplify(:(1+x+3))
    @show simplify(:(x^0))
    @show simplify(:(x^1))
    @show simplify(:(x^(3-2)))
end

end

max_stars_repo_name: ModiaSim/ModiaBase

max_issues_repo_name: ModiaSim/ModiaBase

max_forks_repo_name: ModiaSim/ModiaBase


text:
abstract type AbstractPlanarSurfacePrimitive{T} <: AbstractSurfacePrimitive{T} end
abstract type AbstractCurvedSurfacePrimitive{T} <: AbstractSurfacePrimitive{T} end

include("Plane.jl")
include("ConeMantle.jl")
include("EllipticalSurface.jl")
include("EllipsoidMantle.jl")
include("Polygon.jl")
include("TorusMantle.jl")

# Normal vectors of planar surface primitives do not depend on the point,
# but the normal method is generally called with a point.
normal(p::AbstractPlanarSurfacePrimitive, ::AbstractCoordinatePoint) = normal(p)

function _get_n_points_in_arc_φ(p::AbstractSurfacePrimitive{T}, n_arc::Int64)::Int64 where {T}
    φMin, φMax = get_φ_limits(p)
    f = (φMax - φMin)/T(2π)
    Int(ceil(n_arc*f))
end

function _get_n_points_in_arc_θ(p::AbstractSurfacePrimitive{T}, n_arc::Int64)::Int64 where {T}
    θMin, θMax = get_θ_limits(p)
    f = (θMax - θMin)/T(2π)
    Int(ceil(n_arc*f))
end

function _get_vert_lines_range(p::AbstractSurfacePrimitive{T}, n_arc::Int64, n_vert_lines::Int64)::Vector{Int64} where {T}
    φMin, φMax = get_φ_limits(p)
    if mod(φMax-φMin, T(2π)) == 0
        if n_vert_lines == 0
            []
        elseif n_vert_lines == 1
            [1]
        else
            Int.(ceil.(range(1, (n_arc+1)*(n_vert_lines-1)/n_vert_lines, length = n_vert_lines)))
        end
    else
        Int.(ceil.(range(1, n_arc+1, length = max(n_vert_lines, 2))))
    end
end


text:
using PlanktonIndividuals

grid = RegularRectilinearGrid(size = (16, 16, 16), spacing = (2, 2, 2), halo = (2, 2, 2))

model = PlanktonModel(CPU(), grid; mode = CarbonMode())

TC = sum((interior(model.nutrients.DIC.data, grid) .+
          interior(model.nutrients.DOC.data, grid) .+
          interior(model.nutrients.POC.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TC = TC + sum(model.individuals.phytos.sp1.data.Bm)

uvel = zeros(16,16,16,11)
vvel = zeros(16,16,16,11)
wvel = zeros(16,16,17,11)

for i in 1:11
    uvel[:,:,:,i] .= randn(16,16,16) .* 1e-4
    vvel[:,:,:,i] .= randn(16,16,16) .* 1e-4
    wvel[:,:,:,i] .= randn(16,16,17) .* 1e-4
end

# add boundary conditions for DOC
model.nutrients.DON.bc.west  = 1.0e-3 # west boundary condition of 0.1 mmol/m^2/second
model.nutrients.DON.bc.east  = randn(20,20) .* 1e-3 # east boundary condition of -0.1 mmol/m^2/second
model.nutrients.DON.bc.south = randn(20,20,10) .* 1e-3 # south boundary condition of 0.1 mmol/m^2/second
model.nutrients.DON.bc.north = randn(20,20,10) .* 1e-3 # north boundary condition of -0.1 mmol/m^2/second

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10, vels=(u=uvel, v=vvel, w=wvel))

update!(sim)

TCt = sum((interior(model.nutrients.DIC.data, grid) .+
           interior(model.nutrients.DOC.data, grid) .+
           interior(model.nutrients.POC.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TCt = TCt + sum(model.individuals.phytos.sp1.data.Bm)


text:
#=
Find toplevel items (bind / call)

- downstreams: modules.jl, outline.jl, goto.jl
=#

abstract type ToplevelItem end

struct ToplevelBinding <: ToplevelItem
    expr::CSTParser.EXPR
    bind::CSTParser.Binding
    lines::UnitRange{Int}
end

struct ToplevelCall <: ToplevelItem
    expr::CSTParser.EXPR
    lines::UnitRange{Int}
    callstr::String
end

struct ToplevelTupleH <: ToplevelItem
    expr::CSTParser.EXPR
    lines::UnitRange{Int}
end

"""
    toplevelitems(text; kwargs...)::Vector{ToplevelItem}

Finds and returns toplevel "item"s (call and binding) in `text`.

keyword arguments:
- `mod::Union{Nothing, String}`: if not `nothing` don't return items within modules
    other than `mod`, otherwise enter into every module.
- `inmod::Bool`: if `true`, don't include toplevel items until it enters into `mod`.
"""
function toplevelitems(text; kwargs...)
    parsed = CSTParser.parse(text, true)
    _toplevelitems(text, parsed; kwargs...)
end

function _toplevelitems(
    text, expr,
    items::Vector{ToplevelItem} = Vector{ToplevelItem}(), line = 1, pos = 1;
    mod::Union{Nothing, String} = nothing,
    inmod::Bool = false,
)
    # add items if `mod` isn't specified or in a target modle
    if mod === nothing || inmod
        # binding
        bind = CSTParser.bindingof(expr)
        if bind !== nothing
            lines = line:line+countlines(expr, text, pos, false)
            push!(items, ToplevelBinding(expr, bind, lines))
        end

        lines = line:line+countlines(expr, text, pos, false)

        # toplevel call
        if iscallexpr(expr)
            push!(items, ToplevelCall(expr, lines, str_value_as_is(expr, text, pos)))
        end

        # destructure multiple returns
        if ismultiplereturn(expr)
            push!(items, ToplevelTupleH(expr, lines))
        end
    end

    # look for more toplevel items in expr:
    if shouldenter(expr, mod)
        if expr.args !== nothing
            if ismodule(expr) && shouldentermodule(expr, mod)
                inmod = true
            end
            for arg in expr.args
                _toplevelitems(text, arg, items, line, pos; mod = mod, inmod = inmod)
                line += countlines(arg, text, pos)
                pos += arg.fullspan
            end
        end
    end
    return items
end

function shouldenter(expr::CSTParser.EXPR, mod::Union{Nothing, String})
    !(scopeof(expr) !== nothing && !(
        expr.typ === CSTParser.FileH ||
        (ismodule(expr) && shouldentermodule(expr, mod)) ||
        isdoc(expr)
    ))
end

shouldentermodule(expr::CSTParser.EXPR, mod::Nothing) = true
shouldentermodule(expr::CSTParser.EXPR, mod::String) = expr.binding.name == mod


text:
export random_init!, layered_init!
# Initializations for models

function random_init!(m::Model, nin::Int64, nout::Int64;
                      nhidden::Int64=nin, nreward::Int64=0)
    nc = m.cfg["n_cell_params"]
    for i in 1:nin
        cell = Cell(m.cfg, rand(m.cfg["n_cell_params"]); interface=true)
        cell.params[nc] = 0.25
        push!(m.cells, cell)
        push!(m.inputs, Ref(cell))
    end
    for i in 1:nout
        cell = Cell(m.cfg, rand(m.cfg["n_cell_params"]); interface=true)
        cell.params[nc] = 0.5
        push!(m.cells, cell)
        push!(m.outputs, Ref(cell))
    end
    for i in 1:nreward
        cell = Cell(m.cfg, rand(m.cfg["n_cell_params"]); interface=true)
        cell.params[nc] = 0.75
        push!(m.cells, cell)
        push!(m.rewards, Ref(cell))
    end
    for i in 1:nhidden
        cell = Cell(m.cfg, rand(m.cfg["n_cell_params"]))
        cell.params[nc] = 1.0
        push!(m.cells, cell)
    end
end

function layer_xy(i::Int64, N::Int64)
    N = ceil(Int64, sqrt(N))
    x = fld(i - 1, N) + 1
    y = i - ((x - 1) * N)
    x / (N + 1), y / (N + 1)
end

function layered_init!(m::Model, nin::Int64, nout::Int64;
                       nhidden::Int64=nin, nreward::Int64=1)
    for i in 1:nin
        params = zeros(m.cfg["n_cell_params"])
        params[1], params[2] = layer_xy(i, nin)
        params[3] = 0.0
        params[end] = 0.25
        cell = Cell(m.cfg, params; interface=true)
        push!(m.cells, cell)
        push!(m.inputs, Ref(cell))
    end
    for i in 1:nhidden
        params = zeros(m.cfg["n_cell_params"])
        params[1], params[2] = layer_xy(i, nhidden)
        params[3] = 0.5
        params[end] = 0.5
        cell = Cell(m.cfg, params)
        push!(m.cells, cell)
    end
    for i in 1:nout
        params = zeros(m.cfg["n_cell_params"])
        params[1], params[2] = layer_xy(i, nout)
        params[3] = 1.0
        params[end] = 0.75
        cell = Cell(m.cfg, params; interface=true)
        push!(m.cells, cell)
        push!(m.outputs, Ref(cell))
    end
    for i in 1:nreward
        params = zeros(m.cfg["n_cell_params"])
        params[1], params[2] = layer_xy(i, nreward)
        params[3] = 0.0
        params[end] = 1.0
        cell = Cell(m.cfg, params; interface=true)
        push!(m.cells, cell)
        push!(m.rewards, Ref(cell))
    end
    all_new_conns = Array{Array{Conn, 1}}(undef, 0)
    for source in m.cells
        new_conns = Array{Conn}(undef, 0)
        for dest in m.cells
            if source.params[3] < dest.params[3] && dest.params[end] != 1.0
                inputs = vcat(source.params, get_params(dest, m.cfg["n_cell_params"]))
                params = m.cont.new_conn_params(inputs)
                conn = Conn(m.cfg, source, dest, params)
                push!(new_conns, conn)
            end
        end
        push!(all_new_conns, new_conns)
    end
    for i in eachindex(m.cells)
        append!(m.cells[i].conns, all_new_conns[i])
    end
    reward = m.rewards[1][]
    for cell in m.cells
        if cell != reward
            for dest in cell.conns
                inputs = vcat(reward.params, get_params(dest, m.cfg["n_cell_params"]))
                params = m.cont.new_conn_params(inputs)
                conn = Conn(m.cfg, reward, dest, params)
                push!(reward.conns, conn)
            end
        end
    end
end


text:
# OUTDATED...

using MathProgBase
using GLPKMathProgInterface
using Presolve

function make_lp(m::Int, n::Int, s::Float64)
    #m = rand((1:100))
    #n = rand((1:100))
    #m = 1000
    #n = 1000
    #enron data set contained 36,692 nodes and 183,831 nodes at sparsity of 0.00013
    #m = 36692
    #n = 36692
    c = ones(n)
    x = rand((1:1000),n)
    A = sprand(m,n,s)
    b = A*x
    rowlb = b
    rowub = b
    collb = float(copy(x))
    #ub = 1000000*ones(n)
    #lb = fill(-Inf,n)
    colub = fill(Inf,n)
    #return m,n,c,A,b,lb,ub,x
    return m,n,A,collb,colub,c,rowlb,rowub,x
end

function make_lp(m::Int, n::Int, s::Float64, nice::Bool)
    if(!nice)
        return make_lp(m,n,s)
    end
    c = ones(n)
    x = rand((1:10),n)
    a = randn(m,n)
    for i in 1:m
       for j in 1:n
           if(a[i,j] > 0.5)
               a[i,j] = float(rand((1:10)))
           else
               a[i,j] = 0.
           end
       end
     end
    A = sparse(a)
    b = A*x
    #lb = fill(-Inf,n)
    colub = fill(Inf,n)

    collb = float(copy(x))
    rowlb = b
    rowub = b
    #    ub = 10.0*ones(n)
    #return m,n,c,A,b,lb,ub,x
    return m,n,A,collb,colub,c,rowlb,rowub,x
end

# this is only needed because of the way I have generated sample LP instances. Not significant time cost anyway
function trim(x::Array{Float64,1}, y::Array{Float64,1})
    for i in 1:length(x)
        if(x[i] < y[i])
            x[i] = y[i]
        end
    end
    return x
end

#simplest test, fetch one nice instance and compare.
function correctness_test(in1::Int, in2::Int, in3::Float64, in4::Bool)
    i=1
    j=1
    tol = 1e-3

    while(i<= 5)
        @show i,j
        println("---------STARTING ITERATION $i---------")
        m,n,A,collb,colub,c,rowlb,rowub,x = make_lp(in1,in2,in3,true)
        tolerance = tol * ones(n)
        in4 && @show x
        in4 && @show full(A)
        in4 && @show collb
        in4 && @show colub
        in4 && @show c
        in4 && @show rowlb
        in4 && @show rowub

        # p = Presolve_Problem(true,m,n)

        newA,newcollb,newcolub,newc,newrowlb,newrowub,independentvar,active_constr,pstack = presolver!(false,A,collb,colub,c,rowlb,rowub)

        println("AFTER PRESOLVER CALL -------------")
        in4 && @show size(newA), typeof(newA)
        in4 && @show full(newA)
        in4 && @show size(newcollb),size(newcolub)
        in4 && @show newcollb
        in4 && @show newcolub
        in4 && @show size(newc)
        in4 && @show newc
        in4 && @show size(newrowlb),size(newrowub)
        in4 && @show newrowlb
        in4 && @show newrowub
        in4 && @show size(x), typeof(x)
        in4 && @show independentvar
        ans = Array{Float64,1}()
        #ans = fill(0.0,length(x))
        if(length(find(independentvar))!=0)
            #presol = linprog(newc, newA, '=', newb,newlb,newub, GLPKSolverLP(presolve=false))
            presol = linprog(newc, newA, newrowlb, newrowub, newcollb,newcolub)
            #presol = linprog(newc, newA, '=', newb, -Inf, Inf)
            in4 && println(newc)
            in4 && println(newA)
            in4 && println(newrowlb)
            in4 && println(newrowub)
            in4 && println(newcollb)
            in4 && println(newcolub)

            in4 && @show ans = presol.sol
            in4 && @show independentvar
            if(presol.status!= :Optimal)
                println("$(presol.status)")
                i += 1
                break
            end
            presol.status != :Optimal && error("ERROR - Input feasible problem with an optimal solution but after presolving solver status is not optimal")
            #ans = presol.sol
        end
        #@show pstack
        finalsol = return_postsolved(ans,independentvar,active_constr,pstack)
        #@show finalsol
        finalsol = trim(finalsol,collb)

        in4 && @show x
        in4 && @show finalsol

        if( ((x - finalsol) .< tolerance) == trues(n) )
            println("PASS!")
            j+=1
        else
            @show x - finalsol
            @show A*x
            @show A*finalsol
            if(is_equal(A*x,A*finalsol))
                println("PASS!")
                j+=1
            else
                error("DIDNT PASS!!!")
            end
        end
        i+=1
    end
    if(i == j)
        @show i,j
        println("------Presolve works subject to randomized testing-----")
    end
end

function time_test(in1::Int, in2::Int, in3::Float64, in4::Bool)
    tol = 1e-3
    println("---------STARTING TIME PROFILE TEST--------")
    m,n,c,A,b,lb,ub,x = make_lp(in1,in2,in3,in4)
    @show m,n

    tolerance = tol * ones(n)
    println("GLPK Presolve")
    @time begin
        answer = linprog(c,A,'=',b,lb,ub,GLPKSolverLP(presolve=true))
        answer.status != :Optimal && error("Input feasible problem with an optimal solution but solver status is not optimal")
    end
println("My Presolve")

@time begin
        newc,newA,newb,newlb,newub,independentvar,pstack = presolver!(in4,c,A,b,lb,ub)
        ans = Array{Float64,1}()
    #    @show size(newA)
        #(newA == A) && error("No Presolving Done")
            if(length(find(independentvar))!=0)
            presol = linprog(newc, newA, '=', newb,newlb,newub, GLPKSolverLP(presolve=false))
            #presol = linprog(newc, newA, '=', newb,newlb,newub)
            presol.status != :Optimal && error("Input feasible problem with an optimal solution but after presolving solver status is not optimal")
            ans = presol.sol
            end
        finalsol = return_postsolved(ans,independentvar,pstack)
        finalsol = trim(finalsol,lb)
    end
    @show size(newA)

    #@show answer.sol
    #@show finalsol
end

function noob_test()
    m = 4
    n = 4
    A = sprand(m,n,0.5)
    for i in 1:4
       for j in 1:4
           if(A[i,j] > 0.5)
               A[i,j] = float(rand((1:10)))
           else
               A[i,j] = 0.
           end
       end
   end
   @show full(A)
   lb = fill(-Inf,4)
   ub = fill(Inf,4)
   x = [1,2,3,4]
   b = A*x
   @show b
   c = [1.0,1,1,1]
   newc,newA,newb,newlb,newub,independentvar,pstack = presolver!(true,c,A,b,lb,ub)
   @show size(newA)
   #(newA == A) && error("No Presolving Done")
   ans = Array{Float64,1}()
       if(length(find(independentvar))!=0)
       presol = linprog(newc, newA, '=', newb,newlb,newub, GLPKSolverLP(presolve=false))
       #presol = linprog(newc, newA, '=', newb,newlb,newub)
       presol.status != :Optimal && error("Input feasible problem with an optimal solution but after presolving solver status is not optimal")
       ans = presol.sol
       end
   finalsol = return_postsolved(ans,independentvar,pstack)
   finalsol = trim(finalsol,float(x))
#end
#@show answer.sol
@show finalsol
@show x
end

function do_tests(correctness::Bool, time::Bool)
    if(correctness)
        #correctness_tests
        correctness_test(5,5,0.3,false)
    end

    if(time)
        # Time-Profile tests
        #time_test(10,10,0.5,false)
        #println("AGAIN")
        #time_test(100,100,0.01,false)
        #println("AGAIN")
        time_test(1000,1000,0.001,false)
        #println("AGAIN")
        #time_test(10000,10000,0.0001,false)
        #println("AGAIN")
        #time_test(100000,100000,0.00001,false)
    end
end

`
function checkapi(in1::Int, in2::Int, in3::Float64, in4::Bool)
    tol = 1e-3
    println("-------Testing the API--------")
    m,n,c,A,b,lb,ub,x = make_lp(in1,in2,in3,in4)
    @show m,n

    tolerance = tol * ones(n)
    println("GLPK Presolve")
    @time begin
        answer = linprog(c,A,'=',b,lb,ub,GLPKSolverLP(presolve=true))
        answer.status != :Optimal && error("Input feasible problem with an optimal solution but solver status is not optimal")
    end
println("My Presolve")

@time begin
        answer = linprog(c,A,'=',b,lb,ub,GLPKSolverLP(presolve=false),presolve=true)
        answer.status != :Optimal && error("Input feasible problem with an optimal solution but solver status is not optimal")
    end
    @show size(newA)

end`

println("-------------------RANDOMIZED CORRECTNESS TESTS-----------------")
do_tests(true,false)
#noob_test()

#time_test(1,1,0.3,false)
#Profile.clear()

#println("-------------------RANDOMIZED TIME TESTS---------------------")
#do_tests(false,true)
#@profile do_tests(false,true)

#Profile.print(format=:flat)
`using ProfileView
ProfileView.view()
ProfileView.svgwrite("profile_results2.svg")`

#checkapi(10,10,0.5,false)


text:
fixture_num_states() = 5

fixture_num_actions() = 3

fixture_transition() = cat(3,
    [0.4  0.4  0.1  0.05 0.05;
     0.1  0.35 0.1  0.1  0.35;
     0.25 0.2  0.25 0.05 0.25;
     0.05 0.25 0.3  0.25 0.15;
     0.2  0.2  0.15 0.05 0.40]',
    [0.3  0.15 0.2  0.05 0.3 ;
     0.25 0.1  0.25 0.3  0.1 ;
     0.25 0.2  0.25 0.1  0.2 ;
     0.2  0.15 0    0.4  0.25;
     0.4  0    0.45 0.15 0   ]',
    [0.1  0.2  0.3  0.4  0   ;
     0.15 0.3  0.4  0.1  0.05;
     0.35 0.15 0.3  0.1  0.1 ;
     0.3  0.2  0.3  0.1  0.1 ;
     0.3  0.3  0.2  0.05 0.15]'
)

fixture_reward(shape::Symbol) = shape == :vector ?
    [0.82, -0.45, -0.42, -0.68,  0.37] :
    [0.82 -0.45 -0.42 -0.68  0.37;
     0.49 -0.95  0.46  0.66 -0.57;
    -0.73 -0.31 -0.07 -0.54 -0.29]'

fixture_reward() = fixture_reward(:array)

fixture_initial_value() = [0.04, 0.49, 0.24, 0.45, 0.82]

fixture_discount() = 0.9

fixture_expected_value(reward_shape::Symbol) = reward_shape == :vector ?
    [1.1818, 0.0283499999999999, -0.06405, -0.26015, 0.81325] :
    [1.08955, 0.02835, 0.7993, 1.07985, 0.81325]

fixture_expected_value() = fixture_expected_value(:array)

fixture_expected_policy(reward_shape::Symbol) = reward_shape == :vector ?
    [2, 1, 1, 2, 1] :
    [1, 1, 2, 2, 1]

fixture_expected_policy() = fixture_expected_policy(:array)


text:
using LSH
using Distances
using Base.Test

using MNIST

data = Array(Vector{Float64},60000)
for i = 1:60000
    data[i] = trainfeatures(i)
    data[i] = data[i]/norm(data[i])
end

#println("Building data structure p = $p, N=$(length(D))")
T = @time LSHTable(0.6,LSH.createUHashesAM04(784,4.0,26,2346,0.6,69),data; progress=true)


text:
function hex2binary(x::String)
    # reduce(*, bin.(parse.(Int, collect(x)), 4))
    reduce(*, string.(hex2bytes(x), base=2, pad=8))
end

function puzzle14(input::String="jzgqcdpd")
    input = "jzgqcdpd"
    m = Array{Any,2}(nothing, 128, 128);
    for p in 1:128
        bytes = Int.(collect("$input-$(p-1)"))
        append!(bytes, [17, 31, 73, 47, 23])
        list = knothash(bytes, iterations=64)
        sparsehash = reshape(list, 16, 16)
        densehash = [reduce(xor, sparsehash[:, i]) for i in 1:size(sparsehash, 2)]
        hashstring = reduce(*, lpad.(string.(densehash, base=16), 2, "0"))
        binarystring = hex2binary(hashstring)
        m[p, :] = parse.(Int, collect(binarystring))
    end
    partone = sum(m)

    ref = reshape(1:16384, 128, 128)'
    graph = Dict()
    for i in 1:16384
        if m[i] == 1
            # find coordinates
            coord = findall(ref .== i)
            for c in coord
                x, y = c[1], c[2]
                nb_coords = [[(_x+x, y) for _x=-1:2:1]; [(x, _x+y) for _x=-1:2:1]]
                nb_coords = filter(x -> (0 < x[1][1] < 129) & (0 < x[2][1] < 129), nb_coords)
                neighbors = [ref[x[1], x[2]] for x in collect.(flatten.(nb_coords))]
                # if neighbor is also 1, add to edge list
                graph[i] = [nb for nb in neighbors if m[nb]==1]
            end
        end
    end

    nodes = [first(keys(graph))]
    connected = nodes
    groups = Dict()

    while length(graph) > 0
        while true
            g = []
            for node in nodes
                append!(g, [get(graph, node, nothing)])
            end
            g = reduce(vcat, g)
            nodes = setdiff(g, connected)
            if length(nodes) == 0
                break
            else
                for n in nodes
                    append!(connected, n)
                end
            end
        end
        groups[first(connected)] = connected
        for connection in connected
            delete!(graph, connection)
        end
        length(graph) == 0 ? break : nothing
        nodes = [first(keys(graph))]
        connected = nodes
    end
    parttwo = length(groups)

    return [partone, parttwo]
end


text:
module TestModelSearch

# using Revise
using Test
using MLJModels
import MLJBase.info

pca = info("PCA", pkg="MultivariateStats")
cnst = info("ConstantRegressor", pkg="MLJModels")

@test_throws ArgumentError info("Julia")

@info(ConstantRegressor) == cnst
@test info(Standardizer()) == info("Standardizer", pkg="MLJModels")

@testset "localmodels" begin
    tree = info("DecisionTreeRegressor")
    @test cnst in localmodels(modl=TestModelSearch)
    @test !(tree in localmodels(modl=TestModelSearch))
    import MLJModels
    import DecisionTree
    import MLJModels.DecisionTree_.DecisionTreeRegressor
    @test tree in localmodels(modl=TestModelSearch)
end

@testset "models() and localmodels()" begin
    t(model) = model.is_pure_julia
    mods = models(t)
    @test pca in mods
    @test cnst in mods
    @test !(info("SVC") in mods)
    mods = localmodels(t, modl=TestModelSearch)
    @test cnst in mods
    @test !(pca in mods)
    u(model) = !(model.is_supervised)
    @test pca in models(u, t)
    @test !(cnst in models(u, t))
end

end
true


text:
"""
REWRITE!!
forw(f::Net, input...; kwargs...) applies the compiled knet function f
to the given input.

If f has a return statement, its result is returned, otherwise nothing
is returned.  Following forw(), the program registers can be queried
using get(f,regname).  The special register get(f,:return) contains
the return value, if any.  Note that arrays returned by get(f,regname)
are internal to the Net and will be overwritten in the next forw call.

The input consists of zero or more arrays representing a single
minibatch.  For sequence models, the input typically corresponds to a
minibatch for a single time step.  The input arrays can be any Julia
array type (gpu/cpu, dense/sparse, Float64/32/16).  They get copied
before processing, they are never overwritten.

The keyword arguments are used as boolean condition variables during
function execution.

The internal registers are not cleared between calls.  This should not
be a problem for feed forward nets and is necessary for RNNs which
have read-before-write registers.  An explicit call to reset!(f)
clears all internal registers, which may be necessary at the beginning
of a sequence for RNNs.

forw() saves the state useful for the gradient calculation on f.stack.
This is only necessary during training where the stack will be popped
by back().  During testing apply() should be used instead, which
performs the same computation but does not touch the stack.

"""
forw(f::Net, input...; kwargs...) = _forw(f,false,input...; kwargs...)
sforw(f::Net, input...; kwargs...) = _forw(f,true,input...; kwargs...)

function _forw(f::Net, seq::Bool, input...; kwargs...)
    initforw(f, input...; kwargs...)
    lastinput = 0
    for y in regs(f)
        if getp(y,:forw)
            if isa(y.op, Input)
                y.out = copysync!(y.out0, input[lastinput += 1])
            else
                y.out = forw(y.op, inputs(f,y)..., y.out0; kwargs...)
            end
            # println(:forw, (findfirst(regs(f), y), typeof(y.op), y.argv, size(y.out0), vecnorm0(y.out)))
        end
        seq && push!(f, y)
    end
    # error(:ok)
    r = reg(f,:return)
    return (r==nothing ? r : r.out)
end

### DEAD CODE:
    # yout != nothing && copysync!(yout, r.out[N])
    # loss1 = 0.0
    # if ygold != nothing
    #     r.dif[N] = copysync!(r.dif0[N], ygold)                                      # t:6/628
    #     loss1 = loss(r.op[N], r.dif[N], r.out[N])                               # t:442/628
    # end
    # return loss1

# """
# forw(r::Net,x::Vector) for sequences.
# x can be a Vector of Arrays representing items.
# x can be a Vector of Tuples representing multiple inputs.
# x cannot be a Vector of scalars (DONE:think this over)
# """
# function forw(r::Net, x::Vector, yout=nothing; ygold=nothing, a...)
#     # display((:forwseq0,length(x),vecnorm0(r.out),vecnorm0(r.stack[1:r.sp])))
#     isbits(eltype(x)) && error("forw expects a minibatch")
#     x1 = (isa(x[1], Tuple) ? x[1] : (x[1],))
#     y1 = (ygold == nothing ? nothing : ygold[1])
#     initforw(r, x1...; seq=true, ygold=y1, a...)
#     loss = 0.0
#     for i=1:length(x)
#         xi = (isa(x[i], Tuple) ? x[i] : (x[i],))
#         yi = (yout == nothing ? nothing : yout[i])
#         yg = (ygold == nothing ? nothing : ygold[i])
#         loss += forw(r, xi..., yi; seq=true, ygold=yg, a...)
#     end
#     # display((:forwseq1,length(x),vecnorm0(r.out),vecnorm0(r.stack[1:r.sp])))
#     return loss
# end

# The
# output is always a dense array (on gpu if available) that is internal
# to the net, and will be overwritten in the next call, so it should be
# copied by the caller if needed long term.  initforw allocates out0 if
# necessary but does not initialize out or out0, so for rnn's that have
# forward references a call to reset at the beginning of a sequence is
# necessary.

# # If y.out has been saved on stack, find new storage so we do not overwrite
# function copy_on_write(f::Net,y::Reg)
#     if !ispersistent(y) && haskey(f.sdict,y.out0)
#         y.out0 = nothing
#         initout0(f,y)
#     end
# end


text:
mutable struct sNES_state{T}
    ημ::T
    ησ::T
    samples::Int
    termination::Termination{T}
    function sNES_state{T}(d::Integer;P...) where T
        ημ=T(1)
        ησ=T( (3+log(d))/(5*√d) )
        samples=4 + ceil(Int, log(3*d))
        Term=Termination{T}()
        S=new{T}(ημ,ησ,samples,Term)
        for k in keys(P)
            hasfield(sNES_state{T},k) && setfield!(S,k,P[k])
            hasfield(Termination{T},k) && setfield!(Term,k,P[k])
        end
        S
    end
end

struct sNES <: OptMethod end

function separable_nes(f,x0::AbstractVector{T},σ::AbstractVector{T},params::sNES_state{T}) where T
    samples=params.samples
    samples>3 || error("atleast 3 samples are required for sNES to work.")
    ημ = params.ημ
    ησ = params.ησ
    N=length(x0)
    length(σ) == N || error("the length of 'σ' must be equal to the length of 'x'.")
    x = copy(x0)
    Z = zero(T)
    O = one(T)
    F = fill(f(x0),samples)
    ∇f = fill(Z,N)
    ∇fσ = fill(Z,N)
    ϵ = map(i->zeros(T,N), 1:samples)
    u = utility_function{T}(samples)
    idx = collect(1:samples)
    tmp_x = copy(x0)
    for i in 1:samples
        randn!(ϵ[i])
        tmp_x .= x .+ σ .* ϵ[i]
        F[i] = f(tmp_x)
    end
    sortperm!(idx,F)
    iterations=0
    while !termination_criteria(params.termination,geo_mean(σ),iterations)
        for i in 1:samples
            j = idx[i]
            if i >= (samples ÷ 2)
                randn!(ϵ[j])
            else
                rescaling = sqrt(randchisq(T, N)) / norm(ϵ[j])
                ϵ[j] .*= rescaling
            end
            tmp_x .= x .+ σ .* ϵ[j]
            F[j] = f(tmp_x)
        end
        sortperm!(idx,F)
        ∇f .= Z
        ∇fσ .= Z
        for i in 1:samples
            j = idx[i]
            ∇f .+= u[i] .* ϵ[j]
            ∇fσ .+= u[i] .* (ϵ[j] .* ϵ[j] .- O)
        end
        x .+= ημ .* σ .* ∇f
        σ .*= exp.(ησ ./ 2 .* ∇fσ)
        iterations+=1
    end
    (sol = x, cost = f(x))
end

function separable_nes(f,μ::AbstractVector{T},σ::T,params::sNES_state{T}) where T
    separable_nes(f,μ,fill(σ,length(μ)),params)
end

function separable_nes(f,μ::AbstractVector{T},σ; P...) where T
    separable_nes(f,μ,σ,sNES_state{T}(length(μ);P...))
end

pickmethod(::Type{sNES}) = separable_nes


text:
"""
"""
macro law(fundef)
  s = "The @law macro is only allowed in function definitions"
  @assert isa(fundef,Expr) s
  @assert fundef.head in (:(=), :function) s
  funname = fundef.args[1].args[1]
  nargs = length(fundef.args[1].args)-1
  @assert nargs >= 1 "The @law macro is only allowed in functions with at leats one argument"

  y = fundef.args[1].args[2]
  x =  fundef.args[1].args[3:end]
  a = fundef.args[1].args[2:end]
  q = quote
    function $(funname)($(y)::GridapType,$(x...))
      operate($(funname),$(a...))
    end
    $(fundef)
  end
  :($(esc(q)))
end


text:
import Optim
export joint_svd, joint_svd_gradient_descent

"""
Calculate the joint SVD of two real matrices using gradient descent on the
manifold St(n,n)⊗St(n,n).  This is pretty slow... and not particularly reliable

Reference:

G. Hori. "A general framework for SVD flows and joint SVD flows"
IEEE ICASSP'03 (2003)

"""
function joint_svd_gradient_descent(A::Matrix{T}, B::Matrix{T},
                U0::Union{Matrix{T}, Nothing}, V0::Union{Matrix{T}, Nothing};
                tol=sqrt(eps(T))) where T<:Real

    (N, M) = size(A)

    D = N*M

    to_vector(U, V) = cat(reshape(U, :, 1), reshape(U, :, 1), dims=1)
    ϕ(X) = sum(abs.(diag(X)).^2)
    function to_matrices(X)
        U = reshape(X[1:D], N, M)
        V = reshape(X[D+1:end], N, M)
        return U, V
    end
    function f(X)
        U, V = to_matrices(X)
        return -(ϕ(U'*A*V) + ϕ(U'*B*V))
    end
    U0 = nothing
    V0 = nothing
    if U0 == nothing
        U0 = Float64.(Matrix(I, (n,nn)))
    end
    if V0 == nothing
        V0 = Float64.(Matrix(I, (n,nn)))
    end
    X0 = to_vector(U0, V0)
    #manifold over which to optimize: Steifel manifold of NxN orthogonal matrices
    #Julia is amazing
    manif = Optim.ProductManifold(Optim.Stiefel(), Optim.Stiefel(), (N,M), (N,M))
    result = Optim.optimize(f, X0, Optim.GradientDescent(manifold=manif))

    U, V = to_matrices(Optim.minimizer(result))

    Λa = U'*A*V
    Λb = U'*B*V

    cost(X,Y) = offdiag(X)+offdiag(Y)
    start_cost = cost(A, B)
    final_cost = cost(Λa, Λb)

    success = true
    if (!Optim.converged(result) || (final_cost > tol))
         @warn "Was unable to find a simultaneous SVD using gradient descent:
                    start cost function: $(start_cost), final: $(final_cost)."
         success = false
    end

    U = Matrix(U') #convention shuffling
    V = Matrix(V')

     return (U, V, Diagonal(Λa), Diagonal(Λb), success)
end

"""
Calculatesd the approximate joint singular value decomposition of two real
matrices A, B such that:

A = U' Λa V and B = U' Λb V

Algorithm is from:

J. Miao, G. Cheng, Y. Cai and J. Xia. "Approximate Joint Singular Value
Decomposition Based on Givens-Like Rotation"
IEEE Signal Processing Letter, 25(5) 620-624 (2018).

"""
function joint_svd(A::Matrix{T}, B::Matrix{T}; tol=sqrt(eps(T)),
            maxiter=1000) where T<:Real

    (n,nn) = size(A)

    U = Float64.(Matrix(I, (n,nn)))
    V = Float64.(Matrix(I, (n,nn)))

    Λa = A
    Λb = B

    function MM(X, Y, i, j)
        α1 = X[i,i]*X[j,j] + Y[i,i]*Y[j,j]
        α2 = X[i,i]*X[j,i] + Y[i,i]*Y[j,i]
        α3 = X[j,j]*X[i,j] + Y[j,j]*Y[i,j]
        α4 = X[i,i]*X[i,j] + Y[i,i]*Y[i,j]
        α5 = X[j,j]*X[j,i] + Y[j,j]*Y[j,i]
        β = X[i,i]^2 + X[j,j]^2 + Y[i,i]^2 + Y[j,j]^2
        m = [[β -2*α1]; [-2*α1 β]]
        f = [α2-α3; α4-α5]
        (x,y) = qr(m, Val(true))\f
        return (x,y)
    end

    function R(x, i, j)
        λ = 1/sqrt(1 + x^2)
        r = Float64.(Matrix(I, (n,nn)))
        r[i,i] = λ
        r[j,j] = λ
        r[i,j] = λ*x
        r[j,i] = -λ*x
        return r
    end

    cost(X,Y) = offdiag(X)+offdiag(Y)

    start_cost = cost(A, B)

    target = tol*(norm(A,2)+norm(B,2))
    k = 1

    while ((cost(Λa, Λb) > tol) && (k < maxiter))
        for i=1:n-1
            for j=i+1:n
                (x,y) = MM(Λa, Λb, i, j)
                R1 = R(x, i, j)
                R2 = R(y, i, j)
                U = R1*U
                V = R2*V
                Λa = R1*Λa*R2'
                Λb = R1*Λb*R2'
                k+=1
            end
        end
    end

    final_cost = cost(Λa,Λb)

    success = true
    if k >= maxiter
        @warn "Was unable to find a simultaneous SVD: start cost function:
                    $(start_cost), final: $(final_cost)."
        success = false
    end

    return (U, V, Diagonal(Λa), Diagonal(Λb), success)
end


text:
function log_keyval(lg::TBLogger, key, val, step)
    # Preprocess the as if it was a key=>value pair
    data = Vector{Pair{String,Any}}()
    preprocess("$key", val, data)

    summ    = SummaryCollection()
    for (name,val) in data
        push!(summ.value, summary_impl(name, val))
    end

    # Serialize to file
    write_event(lg.file, make_event(lg, summ, step=step))
end


text:
import Distributions: Exponential, product_distribution

"""
$(TYPEDEF)

A symetric random telegraph noise with switch rate `γ/2` magnitude `b`

$(FIELDS)
"""
struct SymetricRTN
    "Magnitude"
    b::Any
    "Two times the switching probability"
    γ::Any
end

correlation(τ, R::SymetricRTN) = R.b^2 * exp(-R.γ * τ)
spectrum(ω, R::SymetricRTN) = 2 * R.b^2 * R.γ / (ω^2 + R.γ^2)
construct_distribution(R::SymetricRTN) = Exponential(1 / R.γ)

"""
$(TYPEDEF)

An ensemble of random telegraph noise.

$(FIELDS)
"""
struct EnsembleFluctuator{T} <: StochasticBath
    "A list of RTNs"
    f::Vector{T}
end

"""
$(SIGNATURES)

Build the `EnsembleFluctuator` object from a list of amplitudes `b` and a list of switch rates `ω`.

# Examples
```julia-repl
julia> EnsembleFluctuator([1, 1], [1, 2])
Fluctuator ensemble with 2 fluctuators
```
"""
EnsembleFluctuator(b::AbstractArray{T}, ω::AbstractArray{T}) where {T <: Number} =
    EnsembleFluctuator([SymetricRTN(x, y) for (x, y) in zip(b, ω)])
correlation(τ, E::EnsembleFluctuator) = sum((x) -> correlation(τ, x), E.f)
spectrum(ω, E::EnsembleFluctuator) = sum((x) -> spectrum(ω, x), E.f)
construct_distribution(E::EnsembleFluctuator) =
    product_distribution([construct_distribution(x) for x in E.f])

Base.length(E::EnsembleFluctuator) = Base.length(E.f)
Base.show(io::IO, ::MIME"text/plain", E::EnsembleFluctuator) =
    print(io, "Fluctuator ensemble with ", length(E), " fluctuators")
Base.show(io::IO, E::EnsembleFluctuator) =
    print(io, "Fluctuator ensemble with ", length(E), " fluctuators")


text:
# Day 11: Dumbo Octopus
load(::Day{11}, path) = parse.(Int, hcat(split.(readlines(path), "")...)) # transposed!

# Part 1 – count flashes in cellular automaton
function moore_neighborhood(I::CartesianIndex; xmax=10, ymax=10)
    r, c = Tuple(I)
    r₋, r₊ = clamp.([r - 1, r + 1], 1, xmax)
    c₋, c₊ = clamp.([c - 1, c + 1], 1, ymax)
    return CartesianIndices((r₋:r₊, c₋:c₊))
end

function solve1(::Day{11}, data; steps=100)
    data = deepcopy(data) # in-place `solve1!` could comment out this line for performance
    count = 0
    for i in 1:steps
        data .+= 1
        to_spread = Set(findall((==)(10), data))

        while !isempty(to_spread)
            neighbors = moore_neighborhood(pop!(to_spread))
            data[neighbors] .+= 1
            new_flashes = filter(n -> data[n] == 10, neighbors)
            union!(to_spread, new_flashes)
        end

        flashed = data .> 9
        count += sum(flashed)
        data[flashed] .= 0
    end

    return count
end

# Part 2 - return step in which octopuses synchronize
function solve2(::Day{11}, data)
    data = deepcopy(data) # in-place `solve1!` could comment out this line for performance
    step = 1
    while true # iterate until value is returned on synchronization
        data .+= 1
        to_spread = Set(findall((==)(10), data))

        while !isempty(to_spread)
            neighbors = moore_neighborhood(pop!(to_spread))
            data[neighbors] .+= 1
            new_flashes = filter(n -> data[n] == 10, neighbors)
            union!(to_spread, new_flashes)
        end

        flashed = data .> 9
        data[flashed] .= 0

        if sum(flashed) == 100
            return step
        end
        step += 1
    end
end

testresult1(::Day{11}) = 1656
testresult2(::Day{11}) = 195


text:
import TrajectoryOptimization
urdf_folder = joinpath(@__DIR__, "urdfs")
urdf_furuta = joinpath(urdf_folder, "FurutaPendulum.urdf")

furuta_model = TrajectoryOptimization.Model(urdf_furuta)