text:
# ---
# title: 1450. Number of Students Doing Homework at a Given Time
# id: problem1450
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/number-of-students-doing-homework-at-a-given-time/description/>
# hidden: true
# ---
#
# Given two integer arrays `startTime` and `endTime` and given an integer
# `queryTime`.
#
# The `ith` student started doing their homework at the time `startTime[i]` and
# finished it at time `endTime[i]`.
#
# Return _the number of students_ doing their homework at time `queryTime`. More
# formally, return the number of students where `queryTime` lays in the interval
# `[startTime[i], endTime[i]]` inclusive.
#
#
#
# **Example 1:**
#
#
#
#     Input: startTime = [1,2,3], endTime = [3,2,7], queryTime = 4
#     Output: 1
#     Explanation: We have 3 students where:
#     The first student started doing homework at time 1 and finished at time 3 and wasn't doing anything at time 4.
#     The second student started doing homework at time 2 and finished at time 2 and also wasn't doing anything at time 4.
#     The third student started doing homework at time 3 and finished at time 7 and was the only student doing homework at time 4.
#
#
# **Example 2:**
#
#
#
#     Input: startTime = [4], endTime = [4], queryTime = 4
#     Output: 1
#     Explanation: The only student was doing their homework at the queryTime.
#
#
# **Example 3:**
#
#
#
#     Input: startTime = [4], endTime = [4], queryTime = 5
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: startTime = [1,1,1,1], endTime = [1,3,2,4], queryTime = 7
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: startTime = [9,8,7,6,5,4,3,2,1], endTime = [10,10,10,10,10,10,10,10,10], queryTime = 5
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `startTime.length == endTime.length`
#   * `1 <= startTime.length <= 100`
#   * `1 <= startTime[i] <= endTime[i] <= 1000`
#   * `1 <= queryTime <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Estuaries
using DataTables
using DataUtils
using Feather

# const filename = "littletest.feather"
const filename = "bigtest.feather"

function maketestfeather(filename::String; nrows::Integer=10^8)
    data = randomData(Int64, Float64, String, DateTime, nrows=nrows)
    featherWrite(filename, data, overwrite=true)
end

# maketestfeather(filename, nrows=100)

s = Feather.Source(filename)
E = Estuaries.Source(s)

ext: jl


text:
# One iteration of Value Iteration
function valueiterationsolver(mdp::MDP, stage::Int64, util)
    next_stage_value = util               # maximum value in each row
    stage_q = fill(0., (no_states(mdp), length(actions(mdp)))) # q_matrix preinitialization # XXX: use length(stage_states(mdp)) instead of mdp.no_states?

    for s in stage_states(mdp, stage)
        isterminal(mdp, s) && continue

        for a in actions(mdp, s)
            si = stage_stateindex(mdp, s)
            ai = actionindex(mdp, a)
            for (sp, p) in weighted_iterator(transition(mdp, s, a))
                spi = stage_stateindex(mdp, sp)
                stage_q[si, ai] += p * (reward(mdp, s, a, sp) + discount(mdp) * next_stage_value[spi])
            end
        end
    end

    util = maximum(stage_q; dims=2)
    pol = [i[2] for i in findmax(stage_q; dims=2)[2]][:]

    return stage_q, util, pol
end


text:
using Test
using LinearAlgebra
using StaticArrays
using IFGF

K(x,y)   = 1/norm(x-y)
IFGF.wavenumber(::typeof(K)) = 0

@testset "Near field" begin
    @testset "Single leaf" begin
        nx,ny = 100, 200
        nz = 5
        Xpts  = rand(SVector{2,Float64},nx)
        Ypts  = rand(SVector{2,Float64},ny)
        p   = (4,4)
        ds_func  = (x) -> (1/4,2π/4)
        A_mat = [K(x,y) for x in Xpts, y in Ypts]
        B     = rand(ny)
        B_mat = rand(ny,nz)
        C     = zeros(nx)
        A     = assemble_ifgf(K,Xpts,Ypts;order=p,nmax=250)
        mul!(C,A,B,1,0,threads=false)
        @test size(A) == (nx,ny)
        @test C ≈ A_mat*B
        # @test A*B_mat ≈ A_mat*B_mat
    end
    @testset "Tree" begin
        nx,ny = 100, 200
        p   = (4,4)
        nz    = 5
        Xpts  = rand(SVector{2,Float64},nx)
        Ypts  = rand(SVector{2,Float64},ny)
        A_mat = [K(x,y) for x in Xpts, y in Ypts]
        B     = rand(ny)
        B_mat = rand(ny,nz)
        C     = zeros(nx)
        A     = assemble_ifgf(K,Xpts,Ypts;order=p,nmax=100)
        mul!(C,A,B,1,0,threads=false)
        @test size(A) == (nx,ny)
        @test C ≈ A_mat*B
        # @test A*B_mat ≈ A_mat*B_mat
    end
end

@testset "Far field leaf" begin
    @testset "Single leaf" begin
        nx,ny = 100, 200
        nz = 5
        p   = (4,4)
        Xpts  = rand(SVector{2,Float64},nx)
        Ypts  = [SVector(10,10)+rand(SVector{2,Float64}) for _ in 1:ny]
        A_mat = [K(x,y) for x in Xpts, y in Ypts]
        B     = rand(ny)
        B_mat = rand(ny,nz)
        C     = zeros(nx)
        A     = assemble_ifgf(K,Xpts,Ypts;order=p,nmax=250)
        mul!(C,A,B,1,0,threads=false)
        @test size(A) == (nx,ny)
        @test C ≈ A_mat*B
        # @test A*B_mat ≈ A_mat*B_mat
    end
    @testset "Tree" begin
        nx,ny = 100, 200
        nz    = 5
        p   = (4,4)
        Xpts  = rand(SVector{2,Float64},nx)
        Ypts  = [SVector(10,10)+rand(SVector{2,Float64}) for _ in 1:ny]
        A_mat = [K(x,y) for x in Xpts, y in Ypts]
        B     = rand(ny)
        B_mat = rand(ny,nz)
        C     = zeros(nx)
        A     = assemble_ifgf(K,Xpts,Ypts;order=p,nmax=100)
        mul!(C,A,B,1,0;threads=false)
        @test size(A) == (nx,ny)
        @test C ≈ A_mat*B
        # @test A*B_mat ≈ A_mat*B_mat
    end
end

@testset "Near and far field" begin
    nx,ny  = 1000, 1000
    nz = 3
    Xpts   = rand(SVector{3,Float64},nx)
    Ypts   = [SVector(1,1,1)+rand(SVector{3,Float64}) for _ in 1:ny]
    p     = (7,7,7)
    A_mat = [K(x,y) for x in Xpts, y in Ypts]
    B     = rand(ny)
    B_mat = rand(ny,nz)
    C     = zeros(nx)
    A     = assemble_ifgf(K,Xpts,Ypts;order=p,nmax=50)
    bytes_before = Base.summarysize(A)
    mul!(C,A,B,1,0;threads=false)
    @test size(A) == (nx,ny)
    @test C ≈ A_mat*B
    @test C == mul!(C,A,B,1,0;threads=false) # recompute to verify that the result does not change
    @test bytes_before == Base.summarysize(A)
    #@test A*B_mat ≈ A_mat*B_mat
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
export Standardize
export standardize

struct Standardize
    scale::Var
    bias::Var
    runmean
    runvar
end

function Standardize(::Type{T}, insize::Tuple) where T
    dim = 2
    dims = ntuple(i -> i == dim ? 1 : insize[i], length(insize))
    scale = zerograd(ones(T,dims))
    bias = zerograd(zeros(T,dims))
    runmean = zeros(T, dims)
    runvar = ones(T, dims)
    Standardize(scale, bias, runmean, runvar)
end
(f::Standardize)(x, train) = standardize(x, train, f.scale, f.bias, f.runmean, f.runvar)

function standardize(x::Var, train::Bool, scale::Var, bias::Var, runmean, runvar; eps=1e-4, decay=0.9)
    T = eltype(x.data)
    if train
        xmean = mean(x.data, 2)
        xvar = varm(x.data, xmean, 2, corrected = size(x.data,2) > 1)
        n = length(xmean)
        @. runmean = T(decay) * runmean + T(1-decay) * xmean
        @. runvar = T(decay) * runvar + T(1-decay) * xvar
        invstd = T(1) ./ sqrt.(xvar + T(eps))
        xhat = (x.data .- xmean) .* invstd
        data = xhat .* scale.data .+ bias.data
        Var(data, (standardize,x,scale,bias,invstd,xhat))
    else
        data = (x.data .- runmean) ./ sqrt.(runvar + T(eps)) .* scale.data .+ bias.data
        Var(data, (standardize,x,scale,bias))
    end
end

function addgrad!(y::Var, ::typeof(standardize), x::Var, scale::Var, bias::Var, invstd, xhat)
    T = eltype(y.data)
    gscale = sum(y.grad .* xhat, 2)
    gbias = sum(y.grad, 2)

    if !isvoid(x.grad)
        n = size(x.data, 2)
        g = scale.data .* (y.grad .- (xhat .* gscale .+ gbias) / n) .* invstd
        BLAS.axpy!(T(1), g, x.grad)
    end
    isvoid(scale.grad) || BLAS.axpy!(T(1), gscale, scale.grad)
    isvoid(bias.grad) || BLAS.axpy!(T(1), gbias, bias.grad)
end


text:
function _curr_remotes(repodir::String)
    !_check_gitdir(repodir) && return String[]
    out = _run("git -C $(repodir) remote 2>&1"; verbose = false, ignorestatus = true)
    return split(out, "\n"; keepempty = false)
end


text:
using ApproxFun, Base.Test
import Compat: view

## Try constructor variants

ff=(x,y)->exp(-10(x+.2)^2-20(y-.1)^2)*cos(x*y)
gg=x->exp(-10(x[1]+.2)^2-20(x[1]-.1)^2)*cos(x[1]*x[2])
f=Fun(ff,Interval()^2,10000)
@test_approx_eq f(0.,0.) ff(0.,0.)

f=Fun(gg,Interval()^2,10000)
@test_approx_eq f(0.,0.) ff(0.,0.)

f=Fun(ff,Interval()^2)
@test_approx_eq f(0.,0.) ff(0.,0.)
f=Fun(gg,Interval()^2)
@test_approx_eq f(0.,0.) ff(0.,0.)

f=Fun(ff)
@test_approx_eq f(0.,0.) ff(0.,0.)
f=Fun(gg)
@test_approx_eq f(0.,0.) ff(0.,0.)

# Fun +-* constant
f=Fun((x,y)->exp(x)*cos(y))

@test_approx_eq f(0.1,0.2)+2 (f+2)(0.1,0.2)
@test_approx_eq f(0.1,0.2)-2 (f-2)(0.1,0.2)
@test_approx_eq f(0.1,0.2)*2 (f*2)(0.1,0.2)

## ProductFun
u0   = ProductFun((x,y)->cos(x)+sin(y) +exp(-50x.^2-40(y-.1).^2)+.5exp(-30(x+.5).^2-40(y+.2).^2))

@test values(u0)-values(u0|>LowRankFun)|>norm < 1000eps()
@test chebyshevtransform(values(u0))-coefficients(u0)|>norm < 100eps()

##TODO: need to do adaptive to get better accuracy
@test sin(u0)(.1,.2)-sin(u0(.1,.2))|>abs < 10e-4

## LowRankFun

F = LowRankFun((x,y)->besselj0(10(y-x)),Chebyshev(),Chebyshev())

@test_approx_eq F(.123,.456) besselj0(10(.456-.123))

G = LowRankFun((x,y)->besselj0(10(y-x));method=:Cholesky)

@test_approx_eq G(.357,.246) besselj0(10(.246-.357))

F = LowRankFun((x,y)->hankelh1(0,10abs(y-x)),Chebyshev([1.0,2.0]),Chebyshev([1.0im,2.0im]))

@test_approx_eq F(1.5,1.5im) hankelh1(0,10abs(1.5im-1.5))

## Periodic
f=LowRankFun((x,y)->cos(x)*sin(y),PeriodicInterval(),PeriodicInterval())
@test_approx_eq f(.1,.2) cos(.1)*sin(.2)

f=LowRankFun((x,y)->cos(cos(x)+sin(y)),PeriodicInterval(),PeriodicInterval())
@test_approx_eq f(.1,.2) cos(cos(.1)+sin(.2))
@test norm(Float64[cos(cos(x)+sin(y)) for x=ApproxFun.vecpoints(f,1),y=ApproxFun.vecpoints(f,2)]-values(f))<10000eps()

f=ProductFun((x,y)->cos(cos(x)+sin(y)),PeriodicInterval()^2)
@test_approx_eq f(.1,.2) cos(cos(.1)+sin(.2))
x,y=points(f)
@test norm(Float64[cos(cos(x[k,j])+sin(y[k,j])) for k=1:size(f,1),j=1:size(f,2)]-values(f))<10000eps()

d=PeriodicInterval()^2
f=ProductFun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
@test (f.'-f|>coefficients|>norm)< 1000eps()

## Functional*Fun

d=Interval()
B=ldirichlet(d)
f=ProductFun((x,y)->cos(cos(x)*sin(y)),d^2)
@test norm(B*f-Fun(y->cos(cos(-1)*sin(y)),d))<20000eps()
@test norm(f*B-Fun(x->cos(cos(x)*sin(-1)),d))<20000eps()

## 1D in 2D

d=Interval((0.,0.),(1.,1.))
f=Fun(xy->exp(-xy[1]-2cos(xy[2])),d)
@test_approx_eq f(0.5,0.5) exp(-0.5-2cos(0.5))
@test_approx_eq f(FixedSizeArrays.Vec(0.5,0.5)) exp(-0.5-2cos(0.5))

f=Fun(xy->exp(-xy[1]-2cos(xy[2])),d,20)
@test_approx_eq f(0.5,0.5) exp(-0.5-2cos(0.5))

f=Fun((x,y)->exp(-x-2cos(y)),d)
@test_approx_eq f(0.5,0.5) exp(-0.5-2cos(0.5))

f=Fun((x,y)->exp(-x-2cos(y)),d,20)
@test_approx_eq f(0.5,0.5) exp(-0.5-2cos(0.5))

d=Circle((0.,0.),1.)
f=Fun(xy->exp(-xy[1]-2cos(xy[2])),Fourier(d),40)
@test_approx_eq f(cos(0.1),sin(0.1)) exp(-cos(0.1)-2cos(sin(0.1)))
@test_approx_eq f(FixedSizeArrays.Vec(cos(0.1),sin(0.1))) exp(-cos(0.1)-2cos(sin(0.1)))

f=Fun((x,y)->exp(-x-2cos(y)),Fourier(d),40)
@test_approx_eq f(cos(0.1),sin(0.1)) exp(-cos(0.1)-2cos(sin(0.1)))

f=Fun((x,y)->exp(-x-2cos(y)),Fourier(d))
@test_approx_eq f(cos(0.1),sin(0.1)) exp(-cos(0.1)-2cos(sin(0.1)))

println("    Calculus tests")

## Sum

ff=(x,y)->(x-y)^2*exp(-x^2/2.-y^2/2)
f=Fun(ff,[-4.,4.],[-4.,4.])

@test_approx_eq sum(f,1)(0.1) 2.5162377980828357
f=LowRankFun(f)
@test_approx_eq evaluate(f.A,0.1) map(f->f(0.1),f.A)

## Kron operator
Mx=Multiplication(Fun(cos),Chebyshev())
My=Multiplication(Fun(sin),Chebyshev())
K=Mx⊗My
@test_approx_eq ApproxFun.BandedBlockBandedMatrix(view(K,1:10,1:10)) K[1:10,1:10]
C=Conversion(Chebyshev()⊗Chebyshev(),Ultraspherical(1)⊗Ultraspherical(1))
@test_approx_eq C[1:100,1:100] Float64[C[k,j] for k=1:100,j=1:100]

# 2d derivative (issue #346)
let d = Chebyshev()^2
    f = Fun((x,y) -> sin(x) * cos(y), d)
    C=Conversion(Chebyshev()⊗Chebyshev(),Ultraspherical(1)⊗Ultraspherical(1))
    @test_approx_eq (C*f)(0.1,0.2) f(0.1,0.2)
    Dx = Derivative(d, [1,0])
    f = Fun((x,y) -> sin(x) * cos(y), d)
    fx = Fun((x,y) -> cos(x) * cos(y), d)
    @test (Dx*f)(0.2,0.3) ≈ fx(0.2,0.3)
    Dy = Derivative(d, [0,1])
    fy = Fun((x,y) -> -sin(x) * sin(y), d)
    @test (Dy*f)(0.2,0.3) ≈ fy(0.2,0.3)
    L=Dx+Dy
    @test_approx_eq (L*f)(0.2,0.3) (fx(0.2,0.3)+fy(0.2,0.3))

    B=ldirichlet(d[1])⊗ldirichlet(d[2])
    @test_approx_eq B*f f(-1.,-1.)

    B=Evaluation(d[1],0.1)⊗ldirichlet(d[2])
    @test_approx_eq B*f f(0.1,-1.)

    B=Evaluation(d[1],0.1)⊗Evaluation(d[2],0.3)
    @test_approx_eq B*f f(0.1,0.3)

    B=Evaluation(d,(0.1,0.3))
    @test_approx_eq B*f f(0.1,0.3)
end

let d = Space([0,1]) * Space([0,2])
    Dx = Derivative(d, [1,0])
    f = Fun((x,y) -> sin(x) * cos(y), d)
    fx = Fun((x,y) -> cos(x) * cos(y), d)
    @test (Dx*f)(0.2,0.3) ≈ fx(0.2,0.3)
    Dy = Derivative(d, [0,1])
    fy = Fun((x,y) -> -sin(x) * sin(y), d)
    @test (Dy*f)(0.2,0.3) ≈ fy(0.2,0.3)
    L=Dx+Dy
    @test_approx_eq (L*f)(0.2,0.3) (fx(0.2,0.3)+fy(0.2,0.3))

    B=ldirichlet(d[1])⊗ldirichlet(d[2])
    @test abs(B*f-f(0.,0.)) ≤ 10eps()

    B=Evaluation(d[1],0.1)⊗ldirichlet(d[2])
    @test_approx_eq B*f f(0.1,0.)

    B=Evaluation(d[1],0.1)⊗Evaluation(d[2],0.3)
    @test_approx_eq B*f f(0.1,0.3)

    B=Evaluation(d,(0.1,0.3))
    @test_approx_eq B*f f(0.1,0.3)
end


text:
# push!(LOAD_PATH,"../src")
using MPI
using Test
using StructJuMP

case = "./test/PowerGrid/data/case9"
prof = false
contingencies = ["2", "8"]

# Compute reference solution using Ipopt and JuMP

include("./PowerGrid/scopf.jl")

@testset "Testing PIPS Solver Interface" begin

    opfdata = opf_loaddata(case)

    lines_off=Array{Line}(undef, length(contingencies))
    for l in 1:length(lines_off)
        lines_off[l] = opfdata.lines[parse(Int,contingencies[l])]
    end
    @show lines_off
    scopfdata = SCOPFData(opfdata,lines_off)
    scopfmodel = scopf_model(scopfdata,case)
    opfmodel,status = scopf_solve(scopfmodel,scopfdata,case)
    @test status==:Optimal
    if status==:Optimal
        scopf_outputAll(opfmodel,scopfdata)
    end
    refobjective = getobjectivevalue(opfmodel)
    refVm=getvalue(getindex(opfmodel,:Vm))
    refVa=getvalue(getindex(opfmodel,:Va))
    refPg=getvalue(getindex(opfmodel,:Pg))
    refQg=getvalue(getindex(opfmodel,:Qg))

    include("./PowerGrid/scopf_structjump.jl")

    solver = "Ipopt"

    raw = loadrawdata(case)

    scopfmodel, scopfdata = scopf_model(raw, contingencies)

    scopf_init_x(scopfmodel,scopfdata)

    model,status = scopf_solve(scopfmodel,scopfdata, solver, prof)
    @test status == :Solve_Succeeded

    @test refobjective ≈ getObjectiveVal(model)
    @test refVa[:,0] ≈ getvalue(getindex(model,:Va)) atol=1e-5
    @test refVm[:,0] ≈ getvalue(getindex(model,:Vm)) atol=1e-5
    @show refPg
    @show getvalue(getindex(model,:Pg))
    @test refPg ≈ getvalue(getindex(model,:Pg)) atol=1e-5
    @test refQg ≈ getvalue(getindex(model,:Qg)) atol=1e-5
    if status==:Optimal || status == :Solve_Succeeded
        scopf_structjump_outputAll(scopfmodel,scopfdata)
    end

    solver = "PipsNlp"

    raw = loadrawdata(case)

    scopfmodel, scopfdata = scopf_model(raw, contingencies)

    scopf_init_x(scopfmodel,scopfdata)

    model,status = scopf_solve(scopfmodel,scopfdata, solver, prof)
    @test status == :Solve_Succeeded

    @test refobjective ≈ getObjectiveVal(model)
    @test refVa[:,0] ≈ getvalue(getindex(model,:Va)) atol=1e-5
    @test refVm[:,0] ≈ getvalue(getindex(model,:Vm)) atol=1e-5
    @test refPg ≈ getvalue(getindex(model,:Pg)) atol=1e-5
    @test refQg ≈ getvalue(getindex(model,:Qg)) atol=1e-5
    if status==:Optimal || status == :Solve_Succeeded
        scopf_structjump_outputAll(scopfmodel,scopfdata)
    end
end


text:
immutable RLSModel
    R::Float64
    β::Float64
    Δt::Float64
    T::Int64
    nESS::Int64
end

immutable CostParameters
    k1::Float64
    k2::Float64
    C::LinSpace # Capital C refers to the potential C, not a realization
    Cmin::Float64
    nC::Int64
    p::Vector{Float64}
    π::Vector{Float64}
end

type EQ
    P1::Float64
    vN1::Float64
    vI1::Float64
    P2::Float64
    vN2::Float64
    vI2::Float64
end

type Stage
    EQs::Array{EQ,3}
end

type ESS
    number::Array{Int64,1}
    bases::Array{Int64,1}
    nC::Int64
    nESS::Int64
    active::Int64
    ic_to_ess::Array{Int64,3}
    ess_to_ic::Dict{Int64, Array{Int64,1}}
end

function ESS(nC::Int64)
    a,b = index_dict(nC)
    ESS(zeros(Int64,number_of_strings(nC)), ones(Int64,number_of_strings(nC)),nC,number_of_strings(nC),1,a,b)
end

type Statistics
    P1::Float64
    P2::Float64
end

function getindex(x::ESS, ic1::Int64, ic2::Int64, ic::Int64)
    return x.number[x.ic_to_ess[ic1,ic2,ic]]+1
end

function index(nC::Int64, ic1::Int64, ic2::Int64, ic::Int64)
    if (ic1,ic2) == (ic,ic)
        return 1 + div(nC*(nC+1)*(2*nC+1),6) - div(ic*(ic+1)*(2*ic+1),6)
    elseif ic2 == ic
        return 1 + div(nC*(nC+1)*(2*nC+1),6) - div(ic*(ic+1)*(2*ic+1),6) + ic1
    elseif ic1 == ic
        return 1 + div(nC*(nC+1)*(2*nC+1),6) - div(ic*(ic+1)*(2*ic+1),6) + ic - 1 + ic2
    else
        return 1 + div(nC*(nC+1)*(2*nC+1),6) - div(ic*(ic+1)*(2*ic+1),6) + 2*(ic - 1) + sub2ind((ic-1,ic-1),ic1,ic2)
    end
end

function index_dict(nC::Int64)
    count = 1
    nESS = div(nC*(nC+1)*(2*nC+1),6)
    A = Dict{Int64, Array{Int64,1}}()
    B = Array{Int64}(nC, nC, nC)
    for ic = 1:nC
        for ic1 = 1:ic
            for ic2 = 1:ic
                A[index(nC, ic1, ic2, ic)] = [ic1,ic2,ic]
                B[ic1,ic2,ic] = index(nC, ic1, ic2, ic)
                count += 1
            end
        end
    end
    return B,A
end

length(x::ESS) = length(x.number)


text:
import StaticArrays: FieldVector, MVector

"""
    resample!(v::AbstractArray{T, 1}, x::AbstractUncertainValue)
    resample!(v::MVector{N, T}, x::AbstractUncertainValue) where {N, T}
    resample!(v::FieldVector{N, T}, x::AbstractUncertainValue) where {N, T}

    resample!(v::MVector{N, T}, x::Vararg{AbstractUncertainValue, N}) where {N, T}
    resample!(v::FieldVector{N, T}, x::Vararg{AbstractUncertainValue, N}) where {N, T}

    resample!(v::AbstractArray{T, 1}, x::UVAL_COLLECTION_TYPES) where T
    resample!(v::AbstractArray{T, 2}, x::UVAL_COLLECTION_TYPES) where T

    resample!(idxs::AbstractArray{T, 1}, vals::AbstractArray{T, 1},
        x::AbstractUncertainIndexValueDataset) where T
    resample!(idxs::AbstractArray{T, 2}, vals::AbstractArray{T, 2},
        x::AbstractUncertainIndexValueDataset) where T

Resample a uncertain value `x`, or a collection of uncertain values `x`, into a
pre-allocated container `v`.

## Uncertain values

- If `x` is a single uncertain value, and `v` is vector-like, then fill
    `v` with `N` draws of `x`. Works with vectors of length `N`,
    `MVector{N, T}`s or `FieldVector{N, T}`s.

## Uncertain collections

Uncertain collections may be a `Vector{AbstractUncertainValue}` of length `N`,
an  `AbstractUncertainValueDataset`, or an `NTuple{N, AbstractUncertainValue}`.
See also [`UVAL_COLLECTION_TYPES`](@ref).

- If `x` is a collection of uncertain values and `v` is vector-like, then
    fill `v[i]` with a draw of `x[i]` for `i = 1:N`.
- If `x` is a collection of uncertain values and `v` is a 2D-array, then
    fill the `i`-th column of `v` with `length(x)` draws of the `i`-th
    uncertain value in `x`.

## Uncertain index-value collections

- If two mutable vector-like containers, `idxs` and `vals`, are provided along
    with an uncertain index-value dataset `x`, then fill `idxs[i]` with a
    random draw from `x.indices[i]` and fill `vals[i]` with a random draw
    from `x.values[i]`.
- If two mutable matrix-like containers, `idxs` and `vals` are provided along
    with an uncertain index-value dataset `x` (where the number of
    columns in both `idxs` and `vals` matches `length(x)`), then fill the
    `i`-th column of `idxs` with `size(idxs, 1)` draws from `x.indices[i]`,
    and fill the `i`-th column of `vals` with `size(idxs, 1)` draws
    from `x.values[i]`.
"""
function resample! end

function resample!(v, x::AbstractUncertainValue)
    v[:] = rand(x, length(v))

    return v
end

function resample!(v, x::UVAL_COLLECTION_TYPES) where T
    for i in eachindex(v)
        @inbounds v[i] = rand(x[i])
    end

    return v
end

function resample!(v::AbstractArray{T, 2}, x::UVAL_COLLECTION_TYPES) where T
    # The i-th column is filled with random values from the ith uncertain
    # value in the collection.
    n_draws = size(v, 1)
    n_vals = length(x)
    for i in 1:n_vals
        v[:, i] = rand(x[i], n_draws)
    end

    return v
end

function resample!(idxs::Vector{T}, vals::Vector{T},
        x::AbstractUncertainIndexValueDataset) where T

    if !(length(idxs) == length(vals) == length(x))
        error("`length(idxs) == length(vals) == length(x)` evaluated to false")
    end

    for i in eachindex(idxs)
        @inbounds idxs[i] = rand(x.indices[i])
        @inbounds vals[i] = rand(x.values[i])
    end

    return idxs, vals
end

function resample!(idxs::AbstractArray{T, 2}, vals::AbstractArray{T, 2},
        x::AbstractUncertainIndexValueDataset) where T

    if !(size(idxs, 2) == size(vals, 2) == length(x))
        error("`length(idxs) == length(vals) == length(x)` evaluated to false")
    end
    n_draws = size(idxs, 1)
    n_uvals = length(x)

    # The i-th column in `idxs` is filled with random values from the ith uncertain
    # index in the collection, and vice versa for `vals`.
    for i in 1:n_uvals
        @inbounds idxs[:, i] = rand(x.indices[i], n_draws)
        @inbounds vals[:, i] = rand(x.values[i], n_draws)
    end

    return idxs, vals
end

#################################################################################
# Multiple draws of a single uncertain value into mutable vector-like containers.
#################################################################################
function resample!(v::MVector{N, T}, x::AbstractUncertainValue) where {N, T}
    v[:] = resample(x, N)

    return v
end

function resample!(v::FieldVector{N, T}, x::AbstractUncertainValue) where {N, T}
    v[:] = resample(x, N)

    return v
end

##########################################################################################
# A single draw of `N` uncertain values into a length-`N` a mutable vector-like container.
##########################################################################################
function resample!(v::Vector{T}, x::NTuple{N, AbstractUncertainValue}) where {N, T}
    if length(v) != N
        throw(ArgumentError("length(v) == $(length(v)), has to match the number of uncertain values ($N)"))
    end

    for i = 1:N
        @inbounds v[i] = resample(x[i])
    end

    return v
end

function resample!(v::MVector{N, T}, x::Vararg{AbstractUncertainValue, N}) where {N, T}
    @inbounds for i = 1:N
        v[i] = resample(x[i])
    end

    return v
end
function resample!(v::MVector{N, T}, x::NTuple{N, AbstractUncertainValue}) where {N, T}
    for i = 1:N
        @inbounds v[i] = resample(x[i])
    end

    return v
end

function resample!(v::FieldVector{N, T}, x::Vararg{AbstractUncertainValue, N}) where {N, T}
    @inbounds for i = 1:N
        v[i] = resample(x[i])
    end

    return v
end

function resample!(v::FieldVector{N, T}, x::NTuple{N, AbstractUncertainValue}) where {N, T}
    for i = 1:N
        @inbounds v[i] = resample(x[i])
    end

    return v
end


text:
using PyPlot
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../test/useful_routines.jl")
include("../elasticity-solver-uniform-pressure.jl")

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

K1, K2 = 247.0, 247.0    # Pa
mu1, mu2 = 126.0, 126.0   # Pa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
penaltyfactor = 1e2
eta = 1

polyorder = 2
numqp = required_quadrature_order(polyorder) + 2
nelmts = 17
basis = TensorProductBasis(2, polyorder)
interface_center = [0.5, 0.5]
interface_radius = 1.0/3.0

distancefunc(x) = -circle_distance_function(x,interface_center,interface_radius)

nodaldisplacement, mergedmesh, cellquads, facequads, interfacequads =
    nodal_displacement(
        distancefunc,
        stiffness,
        nelmts,
        basis,
        numqp,
        penaltyfactor,
        eta,
    )

refseedpoints, spatialseedpoints, seedcellids =
    CutCellDG.seed_zero_levelset_with_interfacequads(interfacequads, mergedmesh)
interfacenormals =
    CutCellDG.collect_interface_normals(interfacequads, mergedmesh)

spatialpoints = spatialseedpoints[1, :, :]
referencepoints = refseedpoints
referencecellids = seedcellids

relspatialpoints = spatialpoints .- interface_center
angularposition = angular_position(relspatialpoints)
sortidx = sortperm(angularposition)

angularposition = angularposition[sortidx]
spatialpoints = spatialpoints[:, sortidx]
referencepoints = referencepoints[:, :, sortidx]
referencecellids = referencecellids[:, sortidx]
interfacenormals = interfacenormals[:, sortidx]

productdisplacement = CutCellDG.displacement_at_reference_points(
    nodaldisplacement,
    basis,
    referencepoints,
    referencecellids,
    +1,
    mergedmesh,
)
parentdisplacement = CutCellDG.displacement_at_reference_points(
    nodaldisplacement,
    basis,
    referencepoints,
    referencecellids,
    -1,
    mergedmesh,
)
# exactdisplacement =
#     mapslices(x -> displacement(0.1, x), spatialpoints, dims = 1)

fig,ax = PyPlot.subplots(2,1)
ax[1].plot(angularposition,productdisplacement[1,:],label="product u1")
ax[1].plot(angularposition,parentdisplacement[1,:],label="parent u1")
# ax[1].plot(angularposition,exactdisplacement[1,:],"--",label="exact")
ax[2].plot(angularposition,productdisplacement[2,:],label="product u1")
ax[2].plot(angularposition,parentdisplacement[2,:],label="parent u1")
# ax[2].plot(angularposition,exactdisplacement[2,:],"--",label="exact")
ax[1].grid()
ax[2].grid()
ax[1].legend()
ax[2].legend()
fig

productdisplacementmagnitude = vec(mapslices(norm,productdisplacement,dims=1))
parentdisplacementmagnitude = vec(mapslices(norm,parentdisplacement,dims=1))

diffproduct = maximum(productdisplacementmagnitude) - minimum(productdisplacementmagnitude)
diffparent = maximum(parentdisplacementmagnitude) - minimum(parentdisplacementmagnitude)

fig,ax = PyPlot.subplots()
ax.plot(angularposition,productdisplacementmagnitude,label="product")
ax.plot(angularposition,parentdisplacementmagnitude,label="parent")
ax.legend()
ax.grid()
fig

productstress = CutCellDG.stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    +1,
    mergedmesh,
)
producttraction = CutCellDG.traction_force_at_points(productstress,interfacenormals)

parentstress = CutCellDG.stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    -1,
    mergedmesh,
)
parenttraction = CutCellDG.traction_force_at_points(parentstress,interfacenormals)
#
# exactstress = mapslices(x->stress_field(lambda1,mu1,0.1,x),spatialpoints,dims=1)
# exacttraction = CutCellDG.traction_force_at_points(exactstress,interfacenormals)
#
#
# fig,ax = PyPlot.subplots(2,1)
#
# ax[1].plot(angularposition,producttraction[1,:],label="product t1")
# ax[1].plot(angularposition,parenttraction[1,:],label="parent t1")
# ax[1].plot(angularposition,exacttraction[1,:],"--",label="exact t1")
# ax[2].plot(angularposition,producttraction[2,:],label="product t2")
# ax[2].plot(angularposition,parenttraction[2,:],label="parent t2")
# ax[2].plot(angularposition,exacttraction[2,:],"--",label="exact t2")
#
# ax[1].grid()
# ax[2].grid()
# ax[1].legend()
# ax[2].legend()
#
# fig

producttractionmagnitude = vec(mapslices(norm,producttraction,dims=1))
parenttractionmagnitude = vec(mapslices(norm,parenttraction,dims=1))

fig,ax = PyPlot.subplots()
ax.plot(angularposition,producttractionmagnitude)
ax.plot(angularposition,parenttractionmagnitude)
ax.grid()
fig


text:
using Documenter
using Images, TestImages
using Convolution2D

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

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

deploydocs(repo   = "github.com/johnnychen94/Convolution2d.jl.git")


text:
"""
    augment([input], pipeline) -> out

Apply the operations of the given `pipeline` sequentially to the
given image `input` and return the resulting image `out`.

```julia-repl
julia> input = testpattern();

julia> out = augment(input, FlipX() |> FlipY())
3×2 Array{Gray{N0f8},2}:
[...]
```

The parameter `input` can either be a single image, or a tuple of
multiple images. In case `input` is a tuple of images, its elements
will be assumed to be conceptually connected. Consequently, all
images in the tuple will take the exact same path through the
pipeline; even when randomness is involved. This is useful for the
purpose of image segmentation, for which the input and output are
both images that need to be transformed exactly the same way.

```julia
input1 = testpattern()
input2 = Gray.(testpattern())
out1, out2 = augment((input1, input2), FlipX() |> FlipY())
```

The parameter `pipeline` can be a `AugmentorBase.Pipeline`, a tuple
of `AugmentorBase.Operation`, or a single `AugmentorBase.Operation`.

```julia
input = testpattern()
augment(input, FlipX() |> FlipY())
augment(input, (FlipX(), FlipY()))
augment(input, FlipX())
```

If `input` is omitted, Augmentor will use the augmentation test
image provided by the function [`testpattern`](@ref) as the input
image.

```julia
augment(FlipX())
```
"""
function augment(input, pipeline::AbstractPipeline)
    plain_array(_augment(input, pipeline))
end

function augment(input, pipeline::Union{ImmutablePipeline{1},NTuple{1,Operation}})
    augment(input, first(operations(pipeline)))
end

function augment(input, op::Operation)
    plain_array(applyeager(op, input))
end

function augment(op::Union{AbstractPipeline,Operation})
    augment(use_testpattern(), op)
end

@inline function _augment(input, pipeline::AbstractPipeline)
    _augment(input, operations(pipeline)...)
end

@generated function _augment(input, pipeline::Vararg{Operation})
    Expr(:block, Expr(:meta, :inline), augment_impl(:input, pipeline, false))
end

# --------------------------------------------------------------------

"""
    augment!(out, input, pipeline) -> out

Apply the operations of the given `pipeline` sequentially to the
image `input` and write the resulting image into the preallocated
parameter `out`. For convenience `out` is also the function's
return-value.

```julia
input = testpattern()
out = similar(input)
augment!(out, input, FlipX() |> FlipY())
```

The parameter `input` can either be a single image, or a tuple of
multiple images. In case `input` is a tuple of images, the
parameter `out` has to be a tuple of the same length and
ordering. See [`augment`](@ref) for more information.

```julia
inputs = (testpattern(), Gray.(testpattern()))
outs = (similar(inputs[1]), similar(inputs[2]))
augment!(outs, inputs, FlipX() |> FlipY())
```

The parameter `pipeline` can be a `AugmentorBase.Pipeline`, a tuple
of `AugmentorBase.Operation`, or a single `AugmentorBase.Operation`.

```julia
input = testpattern()
out = similar(input)
augment!(out, input, FlipX() |> FlipY())
augment!(out, input, (FlipX(), FlipY()))
augment!(out, input, FlipX())
```
"""
augment!(out, input, op::Operation) = augment!(out, input, (op,))

function augment!(out::AbstractArray, input::AbstractArray, pipeline::AbstractPipeline)
    out_lazy = _augment_avoid_eager(input, pipeline)
    copy!(match_idx(out, axes(out_lazy)), out_lazy)
    out
end

function augment!(outs::NTuple{N,AbstractArray}, inputs::NTuple{N,AbstractArray}, pipeline::AbstractPipeline) where N
    outs_lazy = _augment_avoid_eager(inputs, pipeline)
    map(outs, outs_lazy) do out, out_lazy
        copy!(match_idx(out, axes(out_lazy)), out_lazy)
    end
    outs
end

@inline function _augment_avoid_eager(input, pipeline::AbstractPipeline)
    _augment_avoid_eager(input, operations(pipeline)...)
end

@generated function _augment_avoid_eager(input, pipeline::Vararg{Operation})
    Expr(:block, Expr(:meta, :inline), augment_impl(:input, pipeline, true))
end


text:
function get_removal_parameters(tree)
    components_to_remove_params = []
    parallel_detection = false
    components = listcomponents(tree)
    if length(tree) == 3
        return components_to_remove_params
    end
    for component in components
        component_type = component.Type
        Parent = tree[component.ParentIndex]
        parentoperation = Parent.Type
        if !(component.Parameter in components_to_remove_params)
            sibling = Parent.Child1Index == component.Index ? tree[Parent.Child2Index] : tree[Parent.Child1Index]
            sibling_decendents = descendents(tree,sibling)
            descendant_counter = 1
            while !parallel_detection && descendant_counter <= length(sibling_decendents)
                if sibling_decendents[descendant_counter].Type == '-'
                    parallel_detection = true
                elseif sibling_decendents[descendant_counter].Type == component_type
                    push!(components_to_remove_params,sibling_decendents[descendant_counter].Parameter)
                end
                descendant_counter += 1
            end
        end
        parallel_detection = false
    end
    if length(components_to_remove_params) + 1 == length(components)
        return []
    else
        return components_to_remove_params
    end
end

function simplifycircuit(tree)
    if length(tree) == 3
        return tree
    end
    simplified = false
    while !simplified
        removal_parameters = get_removal_parameters(tree)
        if isempty(removal_parameters)
            simplified = true
        else
            for P in removal_parameters
                if length(tree) > 3
                component = tree[findfirst(x->x.Parameter == P,tree)]
                tree = removecomponent(tree,component)
                end
            end
        end
    end
    return tree
end

function simplifycircuit(circuit::Circuit,terminals="RCPL")
    if count(isoperation,circuit.karva[1:3]) == 1
        return circuit
    end
    circuit_coding_length = length(circuit.karva)
    initialtree = karva_to_tree(circuit.karva,circuit.parameters)
    tree = simplifycircuit(initialtree)
    if tree == initialtree
        return circuit
    end
    rest_karva = circuit.karva[length(initialtree):end]
    rest_parameters = circuit.parameters[length(initialtree):end]
    extra_karva = join(rand(terminals,circuit_coding_length-length(tree)-length(rest_karva)))
    extra_parameters = karva_parameters(extra_karva)
    return Circuit(join([node.Type for node in tree])*rest_karva*extra_karva,vcat(get_tree_parameters(tree),rest_parameters,extra_parameters),nothing)
end

function simplifycircuit!(circuit::Circuit,terminals="RCPL")
    if count(isoperation,circuit.karva[1:3]) != 1
        circuit_coding_length = length(circuit.karva)
        initialtree = karva_to_tree(circuit.karva,circuit.parameters)
        tree = simplifycircuit(initialtree)
        if tree != initialtree
            rest_karva = circuit.karva[length(initialtree):end]
            rest_parameters = circuit.parameters[length(initialtree):end]
            extra_karva = join(rand(terminals,circuit_coding_length-length(tree)-length(rest_karva)))
            extra_parameters = karva_parameters(extra_karva)
            circuit.karva, circuit.parameters =  join([node.Type for node in tree])*rest_karva*extra_karva, vcat(get_tree_parameters(tree),rest_parameters,extra_parameters)
        end
    end
end

function replace_redundant_cpes!(circuit::Circuit)
    cpe_positions = findall(x->x=='P',circuit.karva)
    mutable_karva = split(circuit.karva,"")
    for position in cpe_positions
        if circuit.parameters[position][2] > 0.99 && (1/(circuit.parameters[position][1]) < 10)
            mutable_karva[position] = "C"
            circuit.parameters[position] = 1/(circuit.parameters[position][1])
        elseif circuit.parameters[position][2] < 0.01
            mutable_karva[position] = "R"
            circuit.parameters[position] = circuit.parameters[position][1]
        end
    end
    circuit.karva = join(mutable_karva);
end

#Extra simplification to be included: if the second CPE parameter is sufficiently close to 0.5: it can be replaced by a Warburg element.

function replace_redundant_cpes!(population::Array{Circuit,1})
    for circuit in population
        replace_redundant_cpes!(circuit)
    end
end

function removeduplicates(population::Array{Circuit,1})
    fitnesses = [p.fitness for p in population]
    unique_inds = findfirst.(isequal.(unique(fitnesses)), [fitnesses])
    unique_population = population[unique_inds]
    return unique_population
end

function Terminalnode_replacement(oldtree, oldnode::TreeNode, new_type::Char)
    @assert occursin(new_type,"RCLP") "invalid terminal provided."
    tree = copy(oldtree)
    new_parameter = karva_parameters(new_type)[1]
    tree[oldnode.Index] = TreeNode(oldnode.Index, oldnode.ParentIndex,
     oldnode.Child1Index, oldnode.Child2Index, new_type,new_parameter)
    return tree
end

function CPE_replacements(circuit::Circuit)
    tree = get_circuit_tree(circuit)
    CPE_nodes = tree[findall(x->x.Type=='P',tree)]
    if isempty(CPE_nodes)
        return circuit
    end
    replacement_tree = ""
    replacement_circuits = []
    for new_type in ['R','C']
        for CPE_node in CPE_nodes
        replacement_tree = Terminalnode_replacement(tree, CPE_node, new_type)
        push!(replacement_circuits,Circuit(tree_to_karva(replacement_tree),
        get_tree_parameters(replacement_tree),nothing))
        end
    end

    return replacement_circuits
end

function CPE_testing(circuit::Circuit,measurements,frequencies,cutoff = 0.8)
    replacement_circuits = CPE_replacements(circuit)
    evaluate_fitness!(replacement_circuits,measurements,frequencies)
    if circuit.fitness/minimum(replacement_circuits).fitness >= cutoff
        return minimum(replacement_circuits)
    else
        return circuit
    end
end


text:
#!/usr/bin/env julia
using Test

struct Ellipsoid
    a::Float64  # semi-major axis [m]
    f::Float64  # flattening
    b::Float64  # semi-minor axis [m]
end

function ecef2geodetic(x::Real, y::Real, z::Real, ell::Ellipsoid)
    #=
    convert ECEF (meters) to geodetic coordinates

    input
    -----
    x,y,z  [meters] target ECEF location                             [0,Infinity)
    ell    reference ellipsoid

    output
    ------
    lat,lon   (degrees/radians)
    alt  (meters)

    based on:
    You, Rey-Jer. (2000). Transformation of Cartesian to Geodetic Coordinates without Iterations.
    Journal of Surveying Engineering. doi: 10.1061/(ASCE)0733-9453
    =#

    r = sqrt(x^2 + y^2 + z^2)

    E = sqrt(ell.a^2 - ell.b^2)

    # eqn. 4a
    u = sqrt(0.5 * (r^2 - E^2) + 0.5 * sqrt((r^2 - E^2)^2 + 4 * E^2 * z^2))

    Q = hypot(x, y)

    huE = hypot(u, E)

    # eqn. 4b
    Beta = atan(huE / u * z / hypot(x, y))

    # eqn. 13
    eps = ((ell.b * u - ell.a * huE + E^2) * sin(Beta)) / (ell.a * huE * 1 / cos(Beta) - E^2 * cos(Beta))

    Beta += eps
# %% final output
    lat = atan(ell.a / ell.b * tan(Beta))

    lon = atan(y, x)

    # eqn. 7
    alt = sqrt((z - ell.b * sin(Beta))^2 + (Q - ell.a * cos(Beta))^2)

    lat > -pi / 2 && lat < pi / 2 || throw(DomainError("-90 <= lat <= 90"))
    lon > -pi && lon < 2 * pi || throw(DomainError("-180 <= lat <= 360"))

    return (rad2deg(lat), rad2deg(lon), alt)

end

function geodetic2ecef(lat::Real, lon::Real, alt::Real, ell::Ellipsoid)
    #=
    Point
    input:
    -----
    lat, lon (degrees)
    alt (altitude, meters)    [0, Infinity)
    ell    reference ellipsoid

    output: ECEF x,y,z (meters)
    =#

    lat = deg2rad(lat)
    lon = deg2rad(lon)

    lat > -pi / 2 && lat < pi / 2 || throw(DomainError("-90 <= lat <= 90"))

    lon > -pi && lon < 2 * pi || throw(DomainError("-180 <= lat <= 360"))

    alt > 0 || throw(DomainError("altitude in  [0, Infinity)"))
    # radius of curvature of the prime vertical section
    N = get_radius_normal(lat, ell)
    # Compute cartesian (geocentric) coordinates given  (curvilinear) geodetic
    # coordinates.
    x = (N + alt) * cos(lat) * cos(lon)
    y = (N + alt) * cos(lat) * sin(lon)
    z = (N * (ell.b / ell.a)^2 + alt) * sin(lat)

    return (x, y, z)
end

function enu2uvw(east::Real, north::Real, up::Real,
                 lat0::Real, lon0::Real)

    lat0 = deg2rad(lat0)
    lon0 = deg2rad(lon0)

    t = cos(lat0) * up - sin(lat0) * north
    w = sin(lat0) * up + cos(lat0) * north

    u = cos(lon0) * t - sin(lon0) * east
    v = sin(lon0) * t + cos(lon0) * east

    return (u, v, w)
end

function get_radius_normal(lat_radians, ell)
    # Compute normal radius of planetary body

    a = ell.a
    b = ell.b

    return a^2 / sqrt(a^2 * cos(lat_radians)^2 + b^2 * sin(lat_radians)^2)
end

if basename(PROGRAM_FILE) == basename(@__FILE__)

  if length(ARGS) == 3

  else
    lla0 = (42., -82., 200.)
    ell = Ellipsoid(6378137., 1 / 298.2572235630, 6378137. * (1 - 1 / 298.2572235630))

    x, y, z = geodetic2ecef(lla0..., ell)

    @test isapprox(x, 660675.2518247)
    @test isapprox(y, -4700948.68316)
    @test isapprox(z, 4245737.66222)

    lat1, lon1, alt1 = ecef2geodetic(x,y,z, ell)

    @test isapprox(lat1, lla0[1])
    @test isapprox(lon1, lla0[2])
    @test isapprox(alt1, lla0[3])

  end
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot,Random
Random.seed!(124)
n = 30
x, y, z = randn(n), randn(n), randn(n)
@gp x y 5z z "with points notit pt 4 ps var lc palette"
@gp :- xrange = (-2.2,2.5) yrange = (-2.2,2.2)
save(term="pngcairo size 600,400", output="./code/burbuja2.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
function drop!(db::DB, table::AbstractString; ifexists::Bool = false)
    exists = ifexists ? "IF EXISTS" : ""
    return execute(db, "DROP TABLE $exists $(esc_id(table))")
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
module SummitCloud

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Summit Cloud" => Ahorn.EntityPlacement(
        Maple.SummitCloud
    )
)

function Ahorn.selection(entity::Maple.SummitCloud)
    x, y = Ahorn.position(entity)

    rng = Ahorn.getSimpleEntityRng(entity)

    sprite = rand(rng, textures)
    scaleX = rand(rng, Int[-1, 1])

    return Ahorn.getSpriteRectangle(sprite, x, y, sx=scaleX, sy=1)
end

const textures = String[
    "scenery/summitclouds/cloud00.png",
    "scenery/summitclouds/cloud01.png",
    "scenery/summitclouds/cloud03.png"
]

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SummitCloud, room::Maple.Room)
    Ahorn.Cairo.save(ctx)

    rng = Ahorn.getSimpleEntityRng(entity)

    sprite = rand(rng, textures)
    scaleX = rand(rng, Int[-1, 1])

    Ahorn.Cairo.scale(ctx, scaleX, 1)
    Ahorn.drawSprite(ctx, sprite, 0, 0)

    Ahorn.Cairo.restore(ctx)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
@discord_object struct AlowedMentions
    parse::Vector{String}
    roles::Vector{Snowflake}
    users::Vector{Snowflake}
end


text:
addprocs(3)
using Primes.factor
@assert nprocs() == 4
future = remotecall(factor, 4, 21883298135690819)
@assert isa(future, Future)
factors = fetch(future)
@assert factors == Dict(234711901=>1,93234719=>1)


text:
function _create_rngs(num_rngs::Integer; seed::Integer = 0)
    # Create the RNG to create seeds
    if seed == 0
        # Without a given seed
        rng_master = PCG.PCGStateOneseq()
    else
        # With a given seed
        rng_master = PCG.PCGStateOneseq(seed)
    end
    # From this RNG, create `num_rngs` seeds
    seed_list = rand(rng_master, UInt64, num_rngs)
    # With these seeds, seed the new RNG's
    rng_list = map(Xorshifts.Xorshift1024Star, seed_list)

    # Convert to static vector for faster performance
    seed_list = SVector{num_rngs}(seed_list)
    rng_list = SVector{num_rngs}(rng_list)

    return (seed_list, rng_list)
end

function rng_matrix!(rnd_matrix::AbstractArray, rng_list::AbstractArray)
    rnd_type = eltype(rnd_matrix)
    for i in axes(rnd_matrix, 2)
        rnd_matrix[:, i] = randn(rng_list[i], rnd_type, size(rnd_matrix, 1))
    end
end

function todataframe(positions::AbstractArray, forces::AbstractArray)
    # Save positions
    pos_df = DataFrame()
    pos_df.x = positions[:, 1]
    pos_df.y = positions[:, 2]
    pos_df.z = positions[:, 3]

    # Save forces
    forces_df = DataFrame()
    forces_df.x = forces[:, 1]
    forces_df.y = forces[:, 2]
    forces_df.z = forces[:, 3]

    return pos_df, forces_df
end


text:
# using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood: DrillingCost_TimeFE
using JLD2

JLDFILES = [
    # "E:/projects/ShaleDrillingResults/2019-12/16481962/estimation-results-16481962-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16482113/estimation-results-16482113-no-rigs-data_first_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16482515/estimation-results-16482515-no-rigs-data_first_lease_only_all_leased.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16482550/estimation-results-16482550-no-rigs-data_last_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16486124/estimation-results-16486124-WITH-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16486177/estimation-results-16486177-WITH-rigs-data_first_lease_only_all_leased.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16486227/estimation-results-16486227-WITH-rigs-data_first_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16486272/estimation-results-16486272-WITH-rigs-data_last_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16544476/estimation-results-16544476-WITH-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16627219/estimation-results-16627219-no-rigs-data_first_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16627762/estimation-results-16627762-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16653809/estimation-results-16653809-WITH-rigs-data_first_lease_only.jld2",
    # "E:/projects/ShaleDrillingResults/2019-12/16677572/estimation-results-16677572-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17477499/estimation-results-17477499-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2020-01/17489762/estimation-results-17489762-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2020-01/17491570/estimation-results-17491570-no-rigs-data_all_leases.jld2",

    # "E:/projects/ShaleDrillingResults/2020-01/17491625/estimation-results-17491625-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2020-01/17491863/estimation-results-17491863-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2020-01/17491864/estimation-results-17491864-no-rigs-data_all_leases.jld2",
    # "E:/projects/ShaleDrillingResults/2020-01/17493674/estimation-results-17493674-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17493679/estimation-results-17493679-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17525858/estimation-results-17525858-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17529657/estimation-results-17529657-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17608936/estimation-results-17608936-no-rigs-data_all_leases.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17608940/estimation-results-17608940-no-rigs-data_first_lease_only_all_leased.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17608948/estimation-results-17608948-no-rigs-data_first_lease_only.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17608951/estimation-results-17608951-no-rigs-data_last_lease_only.jld2",
    "E:/projects/ShaleDrillingResults/2020-01/17608965/estimation-results-17608965-no-rigs-data_all_leases.jld2",
]

for JLD2FILE in JLDFILES

    println_time_flush("Loading results from $JLD2FILE")
    file = jldopen(JLD2FILE, "r")
        DATAPATH = file["DATAPATH"]
        M        = file["M"]
        ddmnovf  = file["ddm_novf"]
        theta    = file["theta1"]
    close(file)
end


text:
const MAXUISECT  = 10      # maximum number of sections
const MAXUIITEM  = 80      # maximum number of items per section
const MAXUITEXT  = 500     # maximum number of chars in edittext and other
const MAXUINAME  = 40      # maximum number of chars in name
const MAXUIMULTI = 20      # maximum number of radio/select items in group
const MAXUIEDIT  = 5       # maximum number of elements in edit list
const MAXUIRECT  = 15      # maximum number of rectangles

# key codes matching GLFW (user must remap for other frameworks)
const KEY_ESCAPE    = 256
const KEY_ENTER     = 257
const KEY_TAB       = 258
const KEY_BACKSPACE = 259
const KEY_INSERT    = 260
const KEY_DELETE    = 261
const KEY_RIGHT     = 262
const KEY_LEFT      = 263
const KEY_DOWN      = 264
const KEY_UP        = 265
const KEY_PAGE_UP   = 266
const KEY_PAGE_DOWN = 267
const KEY_HOME      = 268
const KEY_END       = 269
const KEY_F1        = 290
const KEY_F2        = 291
const KEY_F3        = 292
const KEY_F4        = 293
const KEY_F5        = 294
const KEY_F6        = 295
const KEY_F7        = 296
const KEY_F8        = 297
const KEY_F9        = 298
const KEY_F10       = 299
const KEY_F11       = 300
const KEY_F12       = 301

@enum mjtButton::Cint begin        # mouse button
   BUTTON_NONE = 0          # no button
   BUTTON_LEFT              # left button
   BUTTON_RIGHT             # right button
   BUTTON_MIDDLE             # middle button
end

@enum mjtEvent::Cint begin         # mouse and keyboard event type
   EVENT_NONE = 0           # no event
   EVENT_MOVE               # mouse move
   EVENT_PRESS              # mouse button press
   EVENT_RELEASE            # mouse button release
   EVENT_SCROLL             # scroll
   EVENT_KEY                # key press
   EVENT_RESIZE              # resize
end

@enum mjtItem::Cint begin          # UI item type
   ITEM_END = -2            # end of definition list (not an item)
   ITEM_SECTION = -1        # section (not an item)
   ITEM_SEPARATOR = 0       # separator
   ITEM_STATIC              # static text
   ITEM_BUTTON              # button

   # the rest have data pointer
   ITEM_CHECKINT            # check box,value::Cint
   ITEM_CHECKBYTE           # check box, XXXByte value
   ITEM_RADIO               # radio group
   ITEM_SELECT              # selection box
   ITEM_SLIDERINT           # slider,value::Cint
   ITEM_SLIDERNUM           # slider, XXXNum value
   ITEM_EDITINT             # editable array,values::Cint
   ITEM_EDITNUM             # editable array, XXXNum values
   ITEM_EDITTXT             # editable text

   NITEM                     # number of item types
end

# predicate function: set enable/disable based on item category
#typedef int (*mjfItemEnable)(int category, void* data)

struct mjuiState               # mouse and keyboard state
   # constants set by user
   nrect::Cint                  # number of rectangles used
   rect::SVector{MAXUIRECT, mjrRect}  # rectangles (index 0: entire window)
   userdata::Ptr{Cvoid}             # pointer to user data (for callbacks)

   # event type
   _type::Cint                   # (type XXXEvent)

   # mouse buttons
   left::Cint                   # is left button down
   right::Cint                  # is right button down
   middle::Cint                 # is middle button down
   doubleclick::Cint            # is last press a double click
   button::Cint                 # which button was pressed (XXXButton)
   buttontime::Cdouble          # time of last button press

   # mouse position
   x::Cdouble                   # x position
   y::Cdouble                   # y position
   dx::Cdouble                  # x displacement
   dy::Cdouble                  # y displacement
   sx::Cdouble                  # x scroll
   sy::Cdouble                  # y scroll

   # keyboard
   control::Cint                # is control down
   shift::Cint                  # is shift down
   alt::Cint                    # is alt down
   key::Cint                    # which key was pressed
   keytime::Cdouble             # time of last key press

   # rectangle ownership and dragging
   mouserect::Cint              # which rectangle contains mouse
   dragrect::Cint               # which rectangle is dragged with mouse
   dragbutton::Cint             # which button started drag (XXXButton)
end

struct mjuiThemeSpacing        # UI visualization theme spacing
   total::Cint              # total width
   scroll::Cint             # scrollbar width
   label::Cint              # label width
   section::Cint            # section gap
   itemside::Cint           # item side gap
   itemmid::Cint            # item middle gap
   itemver::Cint            # item vertical gap
   texthor::Cint            # text horizontal gap
   textver::Cint            # text vertical gap
   linescroll::Cint         # number of pixels to scroll
   samples::Cint            # number of multisamples
end

struct mjuiThemeColor          # UI visualization theme color
   master::SVector{3, Cfloat}            # master background
   thumb::SVector{3, Cfloat}             # scrollbar thumb
   secttitle::SVector{3, Cfloat}         # section title
   sectfont::SVector{3, Cfloat}          # section font
   sectsymbol::SVector{3, Cfloat}        # section symbol
   sectpane::SVector{3, Cfloat}          # section pane
   shortcut::SVector{3, Cfloat}          # shortcut background
   fontactive::SVector{3, Cfloat}        # font active
   fontinactive::SVector{3, Cfloat}      # font inactive
   decorinactive::SVector{3, Cfloat}     # decor inactive
   decorinactive2::SVector{3, Cfloat}    # inactive slider color 2
   button::SVector{3, Cfloat}            # button
   check::SVector{3, Cfloat}             # check
   radio::SVector{3, Cfloat}             # radio
   select::SVector{3, Cfloat}            # select
   select2::SVector{3, Cfloat}           # select pane
   slider::SVector{3, Cfloat}            # slider
   slider2::SVector{3, Cfloat}           # slider color 2
   edit::SVector{3, Cfloat}              # edit
   edit2::SVector{3, Cfloat}             # edit invalid
   cursor::SVector{3, Cfloat}            # edit cursor
end

struct mjuiItem                # UI item
   # common properties
   _type::Cint                   # type (XXXItem)
   name::SVector{MAXUINAME, Cchar}     # name
   state::Cint                  # 0: disable, 1: enable, 2+: use predicate
   pdata::Ptr{Cvoid}                # data pointer (type-specific)
   sectionid::Cint              # id of section containing item
   itemid::Cint                 # id of item within section

   # JULIA cannot do union declaration as follows, so using largest size
   nelem::Cint
   _name::SVector{MAXUIMULTI, SVector{MAXUINAME, Cchar}}
   # type-specific properties
   #union
   #begin
   #    # check and button-related
   #    struct
   #    begin
   #       modifier::Cint       # 0: none, 1: control, 2: shift; 4: alt
   #       shortcut::Cint       # shortcut key; 0: undefined
   #    end single

   #    # static, radio and select-related
   #    struct
   #    begin
   #       nelem::Cint          # number of elements in group
   #        char name[mjMAXUIMULTI][mjMAXUINAME] # element names
   #    end multi

   #    # slider-related
   #    struct
   #    begin
   #        double range[2]    # slider range
   #        double divisions   # number of range divisions
   #    end slider

   #    # edit-related
   #    struct
   #    begin
   #       nelem::Cint          # number of elements in list
   #        double range[mjMAXUIEDIT][2] # element range (min>=max: ignore)
   #    end edit
   #end

   # internal
   rect::mjrRect                # rectangle occupied by item
end

struct mjuiSection             # UI section
   # properties
   name::SVector{MAXUINAME, Cchar}     # name
   state::Cint                  # 0: closed, 1: open
   modifier::Cint               # 0: none, 1: control, 2: shift; 4: alt
   shortcut::Cint               # shortcut key; 0: undefined
   nitem::Cint                  # number of items in use
   item::SVector{MAXUIITEM, mjuiItem} # preallocated array of items

   # internal
   rtitle::mjrRect              # rectangle occupied by title
   rcontent::mjrRect            # rectangle occupied by content
end

struct mjUI                    # entire UI
   # constants set by user
   spacing::mjuiThemeSpacing    # UI theme spacing
   color::mjuiThemeColor        # UI theme color
   predicate::Function #mjfItemEnable     # callback to set item state programmatically
   userdata::Ptr{Cvoid}             # pointer to user data (passed to predicate)
   rectid::Cint                 # index of this ui rectangle in mjuiState
   auxid::Cint                  # aux buffer index of this ui
   radiocol::Cint               # number of radio columns (0 defaults to 2)

   # UI sizes (framebuffer units)
   width::Cint                  # width
   height::Cint                 # current heigth
   maxheight::Cint              # height when all sections open
   scroll::Cint                 # scroll from top of UI

   # mouse focus
   mousesect::Cint              # 0: none, -1: scroll, otherwise 1+section
   mouseitem::Cint              # item within section
   mousehelp::Cint              # help button down: print shortcuts

   # keyboard focus and edit
   editsect::Cint               # 0: none, otherwise 1+section
   edititem::Cint               # item within section
   editcursor::Cint             # cursor position
   editscroll::Cint             # horizontal scroll
   edittext::SVector{MAXUITEXT, Cchar} # current text
   editchanged::Ptr{mjuiItem}      # pointer to changed edit in last mjui_event

   # sections
   nsect::Cint                  # number of sections in use
   sect::SVector{MAXUISECT, mjuiSection}  # preallocated array of sections
end

struct mjuiDef
   _type::Cint                   # type (XXXItem); -1: section
   name::SVector{MAXUINAME, Cchar}     # name
   state::Cint                  # state
   pdata::Ptr{Cvoid}                # pointer to data
   other::SVector{MAXUITEXT, Cchar}     # string with type-specific properties
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # `MeshArrays.Interpolate` To Arbitrary Coordinates
#
# A `MeshArray` is a collection of arrays that collectively form a global domain grid. Here we interpolate from the global grid to a set of arbitary locations as is commonly done e.g. to compare climate models to sparse field observations.
#
# ```
# (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
# Depth_interpolated=Interpolate(Γ.Depth,f,i,j,w)
# ```
#
# This notebook breaks down the `MeshArray.Interpolate` function, normally called as shown here. In brief, the program finds a grid point quadrilateral (4 grid points) that encloses each chosen location. Computation is chuncked in subdomains (tiles) to allow for parallelism. `MeshArrays.InterpolationFactors` outputs interpolation coefficients -- reusing those is easy and fast, as shown here.

# + [markdown] {"slideshow": {"slide_type": "skip"}}
# ## Initialize Framework
#
# 1. import `MeshArrays` and plotting tools
# 2. choose e.g. a standard `MITgcm` grid
# 3. download the grid if needed
#

# + {"cell_style": "center", "slideshow": {"slide_type": "skip"}}
using MeshArrays, MITgcmTools, Plots

pth=MeshArrays.GRID_LLC90
γ=GridSpec("LatLonCap",pth)
Γ=GridLoad(γ;option="full");
# -

# ## Choose Target Coordinates And Interpolate

# +
lon=collect(0.1:0.5:2.1); lat=collect(0.1:0.5:2.1);
#lon=[66.75]; lat=[-69.75];

(f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
Depth_int=Interpolate(Γ.Depth,f,i,j,w)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## The Interpolation Code? Let's Break It Down:
#
# - find nearest neighbor (`MeshArray` & `set`)
# - define subdomain tiles (`MeshArray` -> `tiles`)
# - exchange and start loop (`tile` & `subset`)
#     - local stereographic projection
#     - define array of quadrilaterals
#     - find enclosing quadrilaterals
#     - compute interpolation coefficients

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# ### Find Nearest Neighbor

# + {"slideshow": {"slide_type": "-"}}
(f,i,j,c)=knn(Γ.XC,Γ.YC,lon,lat)
[write(Γ.XC)[c] write(Γ.YC)[c]]

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# ### Define Subdomain Tiles
#
# _(as done in `InterpolationFactors`)_

# + {"slideshow": {"slide_type": "-"}}
ni=30; nj=30;
τ=Tiles(γ,ni,nj)

tiles=MeshArray(γ,Int);
[tiles[τ[i].face][τ[i].i,τ[i].j].=i for i in 1:length(τ)];

import Plots: heatmap
function heatmap(x::MeshArray; args...)
    n=x.grid.nFaces
    p=()
    for i=1:n; p=(p...,heatmap(x[i]; args...)); end
    plot(p...)
end

heatmap(tiles,title="Tile ID",clims=(0,length(τ)))

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# ### Example Setup
#
# Hereafter we depict the interpolation problem setup:
#
# 1. in geographic coordinates.
# 2. in the local stereographic coordinates used for interpolation.

# + {"slideshow": {"slide_type": "subslide"}}
XCtiles=Tiles(τ,exchange(Γ.XC))
YCtiles=Tiles(τ,exchange(Γ.YC))

iiTile=tiles[f[1]][i[1],j[1]]; iiFace=τ[iiTile].face

ii0=minimum(τ[iiTile].i)+Int(ni/2); jj0=minimum(τ[iiTile].j)+Int(nj/2)
XC0=Γ.XG.f[iiFace][ii0,jj0]; YC0=Γ.YG.f[iiFace][ii0,jj0]
#XC0=66.5000; YC0=-64.4201

scatter(XCtiles[iiTile],YCtiles[iiTile],marker=:+,c=:blue,leg=false,xlabel="longitude",ylabel="latitude")
scatter!([XC0],[YC0],c=:red); scatter!([lon],[lat],c=:green)

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# ### Projection, Quadrilaterals, & Coefficients
#
# _(as done in `InterpolationFactors`)_

# + {"slideshow": {"slide_type": "-"}}
#Local Stereographic Projection
(x_grid,y_grid)=StereographicProjection(XC0,YC0,XCtiles[iiTile],YCtiles[iiTile])
(x_trgt,y_trgt)=StereographicProjection(XC0,YC0,lon,lat)
~isa(x_trgt,Array) ? x_trgt=[x_trgt] : nothing
~isa(y_trgt,Array) ? y_trgt=[y_trgt] : nothing

#Identify Enclosing Quadrilaterals
(x_quad,y_quad,i_quad,j_quad)=MeshArrays.QuadArrays(x_grid,y_grid)
angsum=zeros(size(x_quad,1),size(x_trgt,1))
MeshArrays.PolygonAngle(x_quad,y_quad,x_trgt,y_trgt,angsum)
ii=findall(angsum.>180.)
ii=[ii[j].I[1] for j in 1:length(ii)]

#Interpolation Coefficients
px=permutedims(x_quad[ii[1],:]); py=permutedims(y_quad[ii[1],:])
ox=x_trgt[1]; oy=y_trgt[1]; ow=zeros(4)
MeshArrays.QuadCoeffs(px,py,ox,oy,ow);

# + {"slideshow": {"slide_type": "subslide"}}
scatter(x_grid,y_grid,marker=:+,c=:blue,leg=false,xlabel="x",ylabel="y")
scatter!([0.],[0.],c=:red); scatter!(x_quad[ii,:],y_quad[ii,:],c=:orange)
scatter!(x_trgt,y_trgt,c=:green)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Putting It All Together
#
# Here we just call `InterpolationFactors`. We then interpolate longitude and latitude using the pre-computed coefficients and compare with initial `XC,YC` as a verification.

# +
(f,i,j,w)=InterpolationFactors(Γ,lon,lat)

lon_a=NaN*similar(lon)
lat_a=NaN*similar(lat)
for jj=1:length(lon)
    if !isnan(sum(w[jj,:]))
        x=[Γ.XC[f[jj,ii]][i[jj,ii],j[jj,ii]] for ii=1:4]
        y=[Γ.YC[f[jj,ii]][i[jj,ii],j[jj,ii]] for ii=1:4]
        lon_a[jj]=sum(w[jj,:].*x)
        lat_a[jj]=sum(w[jj,:].*y)
    end
end

#or equivalently:
lon_b=Interpolate(Γ.XC,f,i,j,w)
lat_b=Interpolate(Γ.YC,f,i,j,w)

# + {"slideshow": {"slide_type": "subslide"}}
scatter(XCtiles[iiTile],YCtiles[iiTile],marker=:+,c=:blue,leg=false,xlabel="longitude",ylabel="latitude")
scatter!([XC0],[YC0],c=:red); scatter!(lon_a,lat_a,c=:yellow,marker=:square)
scatter!(lon,lat,c=:green,marker=:star4)

# + {"slideshow": {"slide_type": "skip"}}
#or equivalently:

lon_c=similar(lon)
lat_c=similar(lat)
for jj=1:length(lon)
    x=[XCtiles[iiTile][i_quad[ii[jj],i]+1,j_quad[ii[jj],i]+1] for i=1:4]
    y=[YCtiles[iiTile][i_quad[ii[jj],i]+1,j_quad[ii[jj],i]+1] for i=1:4]
    lon_c[jj]=sum(x.*w[jj,:])
    lat_c[jj]=sum(y.*w[jj,:])
end

scatter(XCtiles[iiTile],YCtiles[iiTile],marker=:+,c=:blue,leg=false,xlabel="longitude",ylabel="latitude")
scatter!([XC0],[YC0],c=:red); scatter!(lon_c,lat_c,c=:yellow,marker=:square)
scatter!(lon,lat,c=:green,marker=:star4)
# -

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
module Trade

import ..Clarus

function cashflows(;params...)
  return Clarus.Services.api_request("Trade","Cashflows";params...)
end

function convert(;params...)
  return Clarus.Services.api_request("Trade","Convert";params...)
end

function price(;params...)
  return Clarus.Services.api_request("Trade","Price";params...)
end

function validate(;params...)
  return Clarus.Services.api_request("Trade","Validate";params...)
end

end #MODULE-END


text:
include("EmbeddingLags.jl")
include("EmbeddingPositions.jl")

"""
    EmbeddingData{N, T}

The data and parameters that goes into a embedding.

## Fields
If `rd` is a `EmbeddingData` instance, then

- **`rd.dataseries`**: gives a vector of dataseries (each itself being
a vector) that was used for the embedding.
- **`rd.in_which_pos`**: indicates which of the dataseries appear in which
columns of the embedding. For example, `rd.in_which_pos[2] = 3` means
that the 3rd timeseries (which can accessed by `rd.dataseries[3]`)
appears in the second column of the embedding.
- **`rd.at_what_lags[i]`**: gives the lag in the i-th column of the embedding.
- **`rd.labels[i]` returns the label of the i-th dataseries
- **`rd.descriptions[i]`**: returns the description of the i-th dataseries.
"""
struct EmbeddingData{D, T}
    dataseries::Vector{Union{AbstractVector{T}, SVector{D, T}}}
    labels::AbstractVector{String}
    descriptions::AbstractVector{String}
    in_which_pos::EmbeddingPositions
    at_what_lags::EmbeddingLags
end

####################################
# Indexing.
####################################
t = Union{Int, Colon, UnitRange, Vector{Int}}
Base.getindex(rl::EmbeddingLags, i::t) = Base.getindex(rl.lags, i)
Base.getindex(rp::EmbeddingPositions, i::t) = Base.getindex(rp.positions, i)

####################################
# Pretty printing.
####################################
function summarise(r::EmbeddingData)
    n  = length(r.dataseries)
    _type = typeof(r)
    """$_type consisting of $n data series."""
end

function description(r::EmbeddingData)
    S = String[]
    lags_s = string.(r.at_what_lags)
    positions_s = join(string.(r.in_which_pos))
    indexed_s = join(string.(1:length(r.dataseries)), " ")
    if all(r.labels .== "")

        push!(S, " The time series are both unnamed")
    else
        lbls = join([" The time series have labels", r.labels], ", ")
        push!(S, lbls)
    end
    push!(S, ", are indexed $indexed_s, and the data series appear")
    push!(S, " as columns of the embedding in $positions_s at $lags_s.")
    S = join([summarise(r); S], "")
    return S
end

Base.show(io::IO, r::EmbeddingData) = println(io, summarise(r))

export EmbeddingData


text:
# TODO: remove
struct PlanarContactMode{T}
    frame::CartesianFrame3D
    projected_points::Vector{SVector{2, T}}
    hull::DConvexHull{T}
end

function PlanarContactMode{T}(frame::CartesianFrame3D) where T
    projected_points = SVector{2, T}[]
    hull = DConvexHull{T}()
    PlanarContactMode{T}(frame, projected_points, hull)
end

_transform_to_root(state::MechanismState, bodyid::BodyID) = transform_to_root(state, bodyid)
_transform_to_root(transforms_to_root::AbstractDict{BodyID, <:Transform3D}, bodyid::BodyID) = transforms_to_root[bodyid]

function update!(mode::PlanarContactMode{T}, active_points::Dict{BodyID, Vector{SPoint3D{T}}},
        transformsource::Union{MechanismState, AbstractDict{BodyID, <:Transform3D}};
        root_to_frame::Transform3D=one(Transform3D{T}, mode.frame)) where T
    empty!(mode.projected_points)
    for (bodyid, points_body) in active_points
        to_root = _transform_to_root(transformsource, bodyid)
        to_frame = root_to_frame * to_root
        for point_body in points_body
            point_frame = to_frame * point_body
            push!(mode.projected_points, horizontal_projection(point_frame.v))
        end
    end
    update!(mode)
    mode
end

function update!(mode::PlanarContactMode)
    jarvis_march!(mode.hull, mode.projected_points)
end


text:
using ZeroBasedIndexing, Base.Test
zerobased()

x = [1, 2, 3]
@test x[0] == 1
x[1] = 5
@test x == [1, 5, 3]


text:
# Solve -u'' + V u + α u^3 = f

using BenchmarkTools
using ChainRulesCore
using FFTW
using FiniteDiff
using IterativeSolvers: gmres # gmres
using KrylovKit: linsolve # linsolve
using LinearAlgebra
using LinearMaps
using NLsolve
using Zygote

N = 100
V(x) = 2 + sin(x)
f(x) = cos(x)
const α = 1.0

# frequencies, avoiding aliasing: eg if N=5, this is 0, 1, 2, -2, 1
ω(i, N) = (i <= div(N,2)+1) ? (i-1) : -(N-i+1)

const x = range(0, 2π, length=N+1)[1:end-1]
const ωs = ω.(1:N, N)
const Vx = V.(x)
const fx = f.(x)

# Zygote specific: need to wrap input u in real(.), otherwise one gets
# implicit promotion to complex numbers in the backwards pass arriving at u.
# real(.) has a custom rrule that ensures real-valued cotangents
# https://github.com/JuliaDiff/ChainRules.jl/blob/65833a19629ee890d30edc96b61bbdbc4e0da72f/src/rulesets/Base/base.jl#L23
lap(u) = real(ifft(ωs.^2 .* fft(real(u)))) # computes u''
precond(u) = real(ifft(fft(real(u)) ./ (1 .+ ωs .^ 2)))

# Goal: find e.g. d(sum(u*)) / dVx via AD

function blackbox_nlsolve_linsolve(Vx)
    residual(u) = -lap(u) .+ @. Vx * u + α * u^3 - fx
    lin_op(u, u_frozen) = -lap(u) .+ @. Vx * u + α * u_frozen^2 * u
    # we solve the fixed-point u_frozen = g(u_frozen) where g(u_frozen) is the solution to the linear equation -u'' + Vu + α ufrozen^2 u = f
    function fp_frozen(u_frozen)
        f_inner(u) = precond(lin_op(u, u_frozen))
        linsolve(f_inner, precond(fx), tol=1e-14)[1]
    end
    f_outer(u) = fp_frozen(u) - u
    u_sol2 = nlsolve(f_outer, zeros(N), method=:anderson).zero
    sum(u_sol2)
end

@time blackbox_nlsolve_linsolve(Vx)

# Approach 1: Direct black-box AD with Zygote (unrolling solvers in time)

# Zygote.gradient(blackbox_nlsolve_linsolve, Vx) # ERROR
### the following error was extracted from the vs code debugger (otherwise caught by Zygote in interface2.jl)
# MethodError: no method matching (::var"#3#4")(::Vector{Float64}, ::Vector{Float64})
# Closest candidates are:
#   (::var"#3#4")(::Any) at ~/Documents/cloud/sandbox/julia/dftk/DifferentiableDFTK/sandbox/1D/1D.jl:42 (method too new to be called from this world context.)
# Stacktrace:
#   [1] value!!(obj::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, F::Vector{Float64}, x::Vector{Float64})
#     @ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:166
#   [2] value!!(obj::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
#     @ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:163
#   [3] anderson_(df::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, beta::Int64, aa_start::Int64, droptol::Float64, cache::NLsolve.AndersonCache{Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}})
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/anderson.jl:73
#   [4] anderson(df::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, beta::Int64, aa_start::Int64, droptol::Float64, cache::NLsolve.AndersonCache{Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}})
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/anderson.jl:203
#   [5] anderson(df::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/anderson.jl:188
#   [6] nlsolve(df::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#29#31", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:30
#   [7] (::NLsolve.var"#nlsolve##kw")(::NamedTuple{(:method, :show_trace), Tuple{Symbol, Bool}}, ::typeof(nlsolve), df::NLSolversBase.NonDifferentiable{Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64})
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:18
#   [8] nlsolve(f::var"#3#4", initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:show_trace,), Tuple{Bool}}})
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
#   [9] (::NLsolve.var"#nlsolve##kw")(::NamedTuple{(:method, :show_trace), Tuple{Symbol, Bool}}, ::typeof(nlsolve), f::var"#3#4", initial_x::Vector{Float64})
#     @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:46
#  [10] top-level scope

#===#
# Approach 2: Add the general rrule of nlsolve

function ChainRulesCore.rrule(config::RuleConfig{>:HasReverseMode}, ::typeof(nlsolve), f, x0; kwargs...)
    result = nlsolve(f, x0; kwargs...)
    function nlsolve_pullback(Δresult)
        Δx = Δresult[].zero
        x = result.zero
        _, f_pullback = rrule_via_ad(config, f, x)
        JT(v) = f_pullback(v)[2] # w.r.t. x
        # solve JT*Δfx = -Δx
        L = LinearMap(JT, length(x0))
        Δfx = gmres(L, -Δx)
        ∂f = f_pullback(Δfx)[1] # w.r.t. f itself (implicitly closed-over variables)
        return (NoTangent(), ∂f, ZeroTangent())
    end
    return result, nlsolve_pullback
end

# 2.1 Try with nesting NLsolve only (since we already defined its rrule above)

function blackbox_nlsolve_nlsolve(Vx)
    residual(u) = -lap(u) .+ @. Vx * u + α * u^3 - fx
    lin_op(u, u_frozen) = -lap(u) .+ @. Vx * u + α * u_frozen^2 * u
    # we solve the fixed-point u_frozen = g(u_frozen) where g(u_frozen) is the solution to the linear equation -u'' + Vu + α ufrozen^2 u = f
    function fp_frozen(u_frozen)
        nlsolve(u -> precond(lin_op(u, u_frozen)) - precond(fx), zeros(N), method=:anderson).zero
    end
    f_outer(u) = fp_frozen(u) - u
    u_sol2 = nlsolve(f_outer, zeros(N), method=:anderson).zero
    sum(u_sol2)
end

let
    @time blackbox_nlsolve_nlsolve(Vx) # 0.605172 seconds (1.47 M allocations: 105.368 MiB, 10.58% gc time, 98.51% compilation time)
    @time g1 = Zygote.gradient(blackbox_nlsolve_nlsolve, Vx)[1] # 10.130100 seconds (16.95 M allocations: 982.691 MiB, 5.67% gc time, 99.35% compilation time)
    @time g2 = FiniteDiff.finite_difference_gradient(blackbox_nlsolve_nlsolve, Vx) # 3.132538 seconds (4.94 M allocations: 846.075 MiB, 10.13% gc time, 14.61% compilation time)
    @show g1[1:5]
    @show g2[1:5]
    @show sum(abs, g1 - g2) / length(g1) # 0.014476332143572825
    @btime blackbox_nlsolve_nlsolve($x) # 21.492 ms (41011 allocations: 8.33 MiB)
    @btime Zygote.gradient(blackbox_nlsolve_nlsolve, $Vx)[1] # 28.739 ms (116532 allocations: 15.70 MiB)
    nothing
end # looks good

# next step: implement an rrule for linsolve

function ChainRulesCore.rrule(config::RuleConfig{>:HasReverseMode}, ::typeof(linsolve), A, b; kwargs...)
    result = linsolve(A, b; kwargs...)
    x = result[1]
    function linsolve_pullback(Δresult)
        Δx = Δresult[1]
        _, A_pullback = rrule_via_ad(config, A, x)
        Aᵀ(v) = A_pullback(v)[2]
        b̄, _ = linsolve(Aᵀ, -Δx; kwargs...) # discard convergence info
        Ā = A_pullback(b̄)[1]  # w.r.t. A itself (implicitly closed-over variables)
        return (NoTangent(), Ā, b̄)
    end
    return result, linsolve_pullback
end

let
    @time blackbox_nlsolve_linsolve(Vx) # 0.636815 seconds (1.73 M allocations: 108.934 MiB, 7.03% gc time, 98.51% compilation time)
    @time g1 = Zygote.gradient(blackbox_nlsolve_linsolve, Vx)[1] # 4.534159 seconds (12.70 M allocations: 809.877 MiB, 6.34% gc time, 92.63% compilation time)
    @time g2 = FiniteDiff.finite_difference_gradient(blackbox_nlsolve_linsolve, Vx) # 2.504600 seconds (5.37 M allocations: 791.759 MiB, 7.66% gc time, 6.94% compilation time)
    @show g1[1:5]
    @show g2[1:5]
    @show sum(abs, g1 - g2) / length(g1) # 1.8373434323529963e-6
    @btime blackbox_nlsolve_linsolve($x) # 12.168 ms (37263 allocations: 5.48 MiB)
    @btime Zygote.gradient(blackbox_nlsolve_linsolve, $Vx)[1] # 293.429 ms (1954880 allocations: 132.05 MiB)
    nothing
end
# works! TODO: figure out speed improvement

# TODO alternatively without KrylovKit: implement rrule for IterativeSolvers.gmres

function ChainRulesCore.rrule(config::RuleConfig{>:HasReverseMode}, ::typeof(gmres), A, b; kwargs...)
    x = gmres(A, b; kwargs...)
    function linsolve_pullback(Δx)
        _, A_pullback = rrule_via_ad(config, A, x)
        Aᵀ(v) = A_pullback(v)[2]
        ∂b = gmres(Aᵀ, -Δx; kwargs...)
        ∂A = A_pullback(∂b)[1]  # w.r.t. A itself (implicitly closed-over variables)
        return (NoTangent(), ∂A, ∂b)
    end
    return x, linsolve_pullback
end

# function blackbox_nlsolve_gmres(Vx)
#     residual(u) = -lap(u) .+ @. Vx * u + α * u^3 - fx
#     lin_op(u, u_frozen) = -lap(u) .+ @. Vx * u + α * u_frozen^2 * u
#     # we solve the fixed-point u_frozen = g(u_frozen) where g(u_frozen) is the solution to the linear equation -u'' + Vu + α ufrozen^2 u = f
#     function fp_frozen(u_frozen)
#         L = LinearMap(u -> precond(lin_op(u, u_frozen)), length(precond(fx)))
#         gmres(L, precond(fx))
#     end
#     u_sol = nlsolve(u -> fp_frozen(u) - u, zeros(N), method = :anderson).zero
#     sum(u_sol)
# end

# let
#     @time blackbox_nlsolve_gmres(Vx) # 0.636815 seconds (1.73 M allocations: 108.934 MiB, 7.03% gc time, 98.51% compilation time)
#     @time g1 = Zygote.gradient(blackbox_nlsolve_gmres, Vx)[1]
# end

max_stars_repo_name: niklasschmitz/DifferentiableDFTK

max_issues_repo_name: niklasschmitz/DifferentiableDFTK

max_forks_repo_name: niklasschmitz/DifferentiableDFTK


text:
using SpecialFunctions: logabsgamma

@doc raw"""
Represents the base-class of all continuous wavelets.
"""
abstract type GenericContinuousWavelet; end

@doc raw"""
    eval_analysis(wav::ContinuousWavelet.GenericContinuousWavelet, t::Float64)

Evaluates the *mother wavelet* at timepoint `t`.
"""
function eval_analysis end

@doc raw"""
    eval_synthesis(wav::ContinuousWavelet.GenericContinuousWavelet, t::Float64)

Evaluates the *mother synthesis wavelet* at timepoint `t`.
"""
function eval_synthesis end

@doc raw"""
    eval_repkern(wav::ContinuousWavelet.GenericContinuousWavelet, a::Float64, b::Float64)

Evaluates the reproducing kernel at scale `a` and time-point `b`.
"""
function eval_repkern end

@doc raw"""
    cutoff_time(wav::ContinuousWavelet.GenericContinuousWavelet)

Returns the cut-off time-point for the given wavelet at scale 1.
"""
function cutoff_time end

@doc raw"""
    cutoff_freq(wav::ContinuousWavelet.GenericContinuousWavelet)

Returns the cut-off frequency for the given wavelet at scale 1.
"""
function cutoff_freq end

include("morlet.jl")
include("cauchy.jl")


text:
# TODO: Check. Document.
"""
x: the running means
y: the new values
"""
function update_running_mean!(x::AbstractVector{<:Real}, y::AbstractVector{<:Real},
                              iter::Integer, nburn::Real)
  iter == nburn && (x .= 0)
  if iter <= nburn
    x .= x * (iter - 1) / iter + y / iter
  else
    update_running_mean!(x, y, iter - nburn, Inf)
  end
end

"""
- `init::State`: Initial state.
- `data::Data`: Data.
- `prior::Prior`: Priors.
- `tuners::Tuners`: Tuners.
- `nsamps::Vector{Int}`: Number of posterior samples for each monitor. Needs to
   be in descending order. Defaults to `[1000]`.
- `nburn::Int`: Number of burn-in iterations. Defaults to `1000`.
- `thin::Int`: Thinning factor for first monitor. Defaults to `1`.
- `monitors::Vector{Vector{Symbol}}`: Parameters to keep samples of. Defaults
   to `default_monitors()`.
- `fix::Vector{Symbol}`: Parameters to hold constant. Defaults to `Symbol[]`.
- `flags::Vector{Flag}`: Available options include `UpdateBetaWithSkewT()`
   and `UpdateLambdaWithSkewT()`.
   Defaults to `default_flags()`.
- `verbose::Int`: How much info to print. Defaults to 1.
- `rep_aux::Int`: How many times to update `zeta, v, lambda` each iteration.
   Defaul=1.
- `callback_fn`: A function with signature `(state::State, data::Data,
   prior::Prior, iter::Int, pbar::ProgressBar)` that returns a dicionary of
   summary statistics. Enables one to add metrics to progress bar.
"""
function fit_via_pseudo_prior(init::State, data::Data, prior::Prior,
                              tuners::Tuners;
                              nsamps::Vector{Int}=[1000], nburn::Int=1000,
                              thin::Int=1, warmup::Integer=10,
                              monitors::Vector{Vector{Symbol}}=default_monitors(),
                              callback_fn::Function=default_callback_fn,
                              fix::Vector{Symbol}=Symbol[],
                              rep_aux::Integer=1, seed=nothing,
                              flags::Vector{Flag}=default_flags(),
                              verbose::Int=1)
  seed == nothing || Random.seed!(seed)

  fix = unique(fix)
  flags = unique(flags)

  println("seed: ", seed)
  println("flags: ", flags)
  println("fix: ", fix)
  println("monitors: ", monitors)

  isfixed(sym::Symbol) = sym in fix

  state0 = deepcopy(init); state0.beta = false
  state1 = deepcopy(init); state1.beta = true
  tuners0 = deepcopy(tuners)
  tuners1 = deepcopy(tuners)

  # Warmup phase for state0, state1.
  for _ in MCMC.ProgressBars.ProgressBar(1:warmup)
    update_theta!(state1, data, prior, tuners1, rep_aux=rep_aux, fix=fix,
                  flags=flags)
    update_theta!(state0, data, prior, tuners0, rep_aux=rep_aux, fix=fix,
                  flags=flags)
    flush(stdout)
  end

  # Define update function.
  function update!(state)
    isfixed(:p) || update_p!(state, data, prior)
    curr_p = state.p
    update_state_via_pseudo_prior!(state, state0, state1, data, prior, tuners,
                                   tuners0, tuners1, rep_aux=rep_aux, fix=fix,
                                   flags=flags)
    state.p = curr_p
  end

  # Initialize beta_hat tracker.
  beta_hat = [0.0]

  # Define callback.
  function _callback_fn(state::State, iter::Int,
                        pbar::MCMC.ProgressBars.ProgressBar)
    update_running_mean!(beta_hat, [state.beta], iter, nburn)
    cb_out = callback_fn(state, data, prior, iter, pbar)
    pbar.postfix = (beta_hat=round(beta_hat[1], digits=3), pbar.postfix...)
    return cb_out
  end

  # Run Gibbs sampler.
  chain, laststate, summarystats = MCMC.gibbs(
      init, update!, nsamps=nsamps, nburn=nburn, thin=thin, monitors=monitors,
      callback_fn=_callback_fn)

  return chain, laststate, summarystats
end

ppfit = fit_via_pseudo_prior

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
import Base.Math.@horner

# BinomialRmathSampler is non-exported, which is mainly used
# for benchmarking or testing purpose
#
immutable BinomialRmathSampler <: Sampleable{Univariate,Discrete}
    n::Int
    prob::Float64
end

rand(s::BinomialRmathSampler) =
    int(ccall((:rbinom, "libRmath-julia"), Float64, (Float64, Float64), s.n, s.prob))

# compute probability vector of a Binomial distribution
function binompvec(n::Int, p::Float64)
    pv = Array(Float64, n+1)
    if p == 0.0
        fill!(pv, 0.0)
        pv[1] = 1.0
    elseif p == 1.0
        fill!(pv, 0.0)
        pv[n+1] = 1.0
    else
        q = 1.0 - p
        a = p / q
        @inbounds pv[1] = pk = q ^ n
        for k = 1:n
            @inbounds pv[k+1] = (pk *= ((n - k + 1) / k) * a)
        end
    end
    return pv
end

# Remainder term after Stirling's approximation to the log-gamma function
# lstirling(x) = lgamma(x) + x - (x-0.5)*log(x) - 0.5*log2π
#              = 1/(12x) - 1/(360x^3) + 1/(1260x^5) + ...
# Asymptotic expansion from:
#   Temme, N. (1996) Special functions: An introduction to the classical
#   functions of mathematical physics, Wiley, New York, ISBN: 0-471-11313-1,
#   Chapter 3.6, pp 61-65.
# Relative error of approximation is bounded by
#   (174611/125400 x^-19) / (1/12 x^-1 - 1/360 x^-3)
# which is < 1/2 ulp for x >= 10.0
# total numeric error appears to be < 2 ulps
#
lstirling_asym(x::Integer) = lstirling_asym(float(x))

function lstirling_asym(x::Float64)
    t = 1.0/(x*x)
    @horner(t,
             8.33333333333333333e-2, #  1/12 x^-1
            -2.77777777777777778e-3, # -1/360 x^-3
             7.93650793650793651e-4, #  1/1260 x^-5
            -5.95238095238095238e-4, # -1/1680 x^-7
             8.41750841750841751e-4, #  1/1188 x^-9
            -1.91752691752691753e-3, # -691/360360 x^-11
             6.41025641025641026e-3, #  1/156 x^-13
            -2.95506535947712418e-2, # -3617/122400 x^-15
             1.79644372368830573e-1)/x #  43867/244188 x^-17
end

function lstirling_asym(x::Float32)
    t = 1f0/(x*x)
    @horner(t,
             8.333333333333f-2, #  1/12 x^-1
            -2.777777777777f-3, # -1/360 x^-3
             7.936507936508f-4, #  1/1260 x^-5
            -5.952380952381f-4, # -1/1680 x^-7
             8.417508417508f-4)/x #  1/1188 x^-9
end

# Geometric method:
#
#   Devroye. L.
#   "Generating the maximum of independent identically  distributed random variables"
#   Computers and Marhemafics with Applicalions 6, 1960, 305-315.
#
immutable BinomialGeomSampler <: Sampleable{Univariate,Discrete}
    comp::Bool
    n::Int
    scale::Float64
end

BinomialGeomSampler() = BinomialGeomSampler(false, 0, 0.0)

function BinomialGeomSampler(n::Int, prob::Float64)
    if prob <= 0.5
        comp = false
        scale = -1.0/log1p(-prob)
    else
        comp = true
        scale = prob < 1.0 ? -1.0/log(prob) : Inf
    end
    BinomialGeomSampler(comp, n, scale)
end

function rand(s::BinomialGeomSampler)
    y = 0
    x = 0
    n = s.n
    while true
        er = Base.Random.randmtzig_exprnd()
        v = er * s.scale
        if v > n  # in case when v is very large or infinity
            break
        end
        y += iceil(v)
        if y > n
            break
        end
        x += 1
    end
    (s.comp ? s.n - x : x)::Int
end

# BTPE algorithm from:
#
#   Kachitvichyanukul, V.; Schmeiser, B. W.
#   "Binomial random variate generation."
#   Comm. ACM 31 (1988), no. 2, 216–222.
#
# Note: only use this sampler when n * min(p, 1-p) is large enough
#       e.g., it is greater than 20.
#
immutable BinomialTPESampler <: Sampleable{Univariate,Discrete}
    comp::Bool
    n::Int
    r::Float64
    q::Float64
    nrq::Float64
    M::Float64
    Mi::Int
    p1::Float64
    p2::Float64
    p3::Float64
    p4::Float64
    xM::Float64
    xL::Float64
    xR::Float64
    c::Float64
    λL::Float64
    λR::Float64
end

BinomialTPESampler() =
    BinomialTPESampler(false, 0, 0., 0., 0., 0., 0,
                       0., 0., 0., 0., 0., 0., 0., 0., 0., 0.)

function BinomialTPESampler(n::Int, prob::Float64)
    if prob <= 0.5
        comp = false
        r = prob
        q = 1.0 - prob
    else
        comp = true
        r = 1.0 - prob
        q = prob
    end

    nrq = n*r*q
    fM = (n+1)*r #
    M = floor(fM)
    Mi = integer(M)
    p1 = floor(2.195*sqrt(nrq)-4.6*q) + 0.5
    xM = M+0.5
    xL = xM-p1
    xR = xM+p1
    c = 0.134 + 20.5/(15.3+M)
    a = (fM-xL)/(fM-xL*r) #
    λL = a*(1.0 + 0.5*a)
    a = (xR-fM)/(xR*q) #
    λR = a*(1.0 + 0.5*a)
    p2 = p1*(1.0 + 2.0*c)
    p3 = p2 + c/λL
    p4 = p3 + c/λR

    BinomialTPESampler(comp,n,r,q,nrq,M,Mi,p1,p2,p3,p4,
                        xM,xL,xR,c,λL,λR)
end

function rand(s::BinomialTPESampler)
    y = 0
    while true
        # Step 1
        u = s.p4*rand()
        v = rand()
        if u <= s.p1
            y = ifloor(s.xM-s.p1*v+u)
            # Goto 6
            break
        elseif u <= s.p2 # Step 2
            x = s.xL + (u-s.p1)/s.c
            v = v*s.c+1.0-abs(s.M-x+0.5)/s.p1
            if v > 1
                # Goto 1
                continue
            end
            y = ifloor(x)
            # Goto 5
        elseif u <= s.p3 # Step 3
            y = ifloor(s.xL + log(v)/s.λL)
            if y < 0
                # Goto 1
                continue
            end
            v *= (u-s.p2)*s.λL
            # Goto 5
        else # Step 4
            y = ifloor(s.xR-log(v)/s.λR)
            if y > s.n
                # Goto 1
                continue
            end
            v *= (u-s.p3)*s.λR
            # Goto 5
        end

        # Step 5
        # 5.0
        k = abs(y-s.Mi)
        if (k <= 20) || (k >= 0.5*s.nrq-1)
            # 5.1
            S = s.r/s.q
            a = S*(s.n+1)
            F = 1.0
            if s.Mi < y
                for i = (s.Mi+1):y
                    F *= a/i-S
                end
            elseif s.Mi > y
                for i = (y+1):s.Mi
                    F /= a/i-S
                end
            end
            if v > F
                # Goto 1
                continue
            end
            # Goto 6
            break
        else
            # 5.2
            ρ = (k/s.nrq)*((k*(k/3.0+0.625)+1.0/6.0)/s.nrq+0.5)
            t = -k^2/(2.0*s.nrq)
            A = log(v)
            if A < t - ρ
                # Goto 6
                break
            elseif A > t + ρ
                # Goto 1
                continue
            end

            # 5.3
            x1 = float64(y+1)
            f1 = float64(s.Mi+1)
            z = float64(s.n+1-s.Mi)
            w = float64(s.n-y+1)

            if A > (s.xM*log(f1/x1) + ((s.n-s.Mi)+0.5)*log(z/w) + (y-s.Mi)*log(w*s.r/(x1*s.q)) +
                    lstirling_asym(f1) + lstirling_asym(z) + lstirling_asym(x1) + lstirling_asym(w))
                # Goto 1
                continue
            end

            # Goto 6
            break
        end
    end
    # 6
    (s.comp ? s.n - y : y)::Int
end

# Constructing an alias table by directly computing the probability vector
#
immutable BinomialAliasSampler <: Sampleable{Univariate,Discrete}
    table::AliasTable
end

BinomialAliasSampler(n::Int, p::Float64) =
    BinomialAliasSampler(make_alias_table!(binompvec(n, p)))

rand(s::BinomialAliasSampler) = rand(s.table) - 1

# Integrated Polyalgorithm sampler that automatically chooses the proper one
#
# It is important for type-stability
#
type BinomialPolySampler <: Sampleable{Univariate,Discrete}
    use_btpe::Bool
    geom_sampler::BinomialGeomSampler
    btpe_sampler::BinomialTPESampler
end

function BinomialPolySampler(n::Int, p::Float64)
    q = 1.0 - p
    if n * min(p, q) > 20
        use_btpe = true
        geom_sampler = BinomialGeomSampler()
        btpe_sampler = BinomialTPESampler(n, p)
    else
        use_btpe = false
        geom_sampler = BinomialGeomSampler(n, p)
        btpe_sampler = BinomialTPESampler()
    end
    BinomialPolySampler(use_btpe, geom_sampler, btpe_sampler)
end

BinomialPolySampler(n::Real, p::Real) = BinomialPolySampler(int(n), float64(p))

rand(s::BinomialPolySampler) = s.use_btpe ? rand(s.btpe_sampler) : rand(s.geom_sampler)


text:
function get_executable_path(package::AbstractString, exec::AbstractString)
    base_path = Base.find_package(package)
    sep = joinpath("src", "$package.jl")
    exec_path = joinpath(split(base_path, sep)[1], exec)

    return exec_path
end

function symlink_user_bin(path_::AbstractString)
    exec = splitdir(path_)[end]
    if Sys.iswindows()
        bin_path = split(path_, "\\$exec")[1]
        @warn bin_path
        run(`setx path "%path%;$bin_path"`)
    else
        ln = "/usr/local/bin/$exec"
        rm(ln, force=true)
        symlink(path_, ln)
        @warn "Created symlink: $ln -> $path_"
    end
end

path_ = get_executable_path("VariantVisualization", "viva")
symlink_user_bin(path_)