text:
# Generate observed linearized data for the Sigsbee2A velocity model
# Author: Philipp Witte, pwitte.slim@gmail.com
# Date: May 2018
#

# TO DO:
# Replace w/ full path the directory where the observed data will be stored
file_path="/path/to/data/"
file_name="sigsbee2A_marine"

using JUDI, PyPlot, JLD, SegyIO

# Load Sigsbee model
if !isfile("sigsbee2A_model.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/sigsbee2A_model.jld`)
end
M = load("sigsbee2A_model.jld")

# Setup info and model structure
model0 = Model(M["n"], M["d"], M["o"], M["m0"])
dm = vec(M["dm"])

# Source interval [m]
source_spacing = 25

## Set up receiver geometry
domain_x = (model0.n[1] - 1)*model0.d[1]
nrec = 1200     # no. of receivers
xmin = 500f0
xmax = domain_x - 500f0
min_offset = 100f0  # [m]
max_offset = 12000f0
xmid = domain_x / 2

# Source/receivers
nsrc = Int(length(xmin:source_spacing:xmax))
xrec = Array{Any}(undef, nsrc)
yrec = Array{Any}(undef, nsrc)
zrec = Array{Any}(undef, nsrc)
xsrc = Array{Any}(undef, nsrc)
ysrc = Array{Any}(undef, nsrc)
zsrc = Array{Any}(undef, nsrc)

# Vessel goes from left to right in right-hand side of model
nsrc_half = Int(nsrc/2)
for j=1:nsrc_half
    xloc = xmid + (j-1)*source_spacing
    xrec[j] = range(xloc - max_offset, xloc - min_offset, length=nrec)
    yrec[j] = 0.
    zrec[j] = range(50f0, 50f0, length=nrec)
    xsrc[j] = xloc
    ysrc[j] = 0f0
    zsrc[j] = 20f0
end

# Vessel goes from right to left in left-hand side of model
for j=1:nsrc_half
    xloc = xmid - (j-1)*source_spacing
    xrec[nsrc_half + j] = range(xloc + min_offset, xloc + max_offset, length=nrec)
    yrec[nsrc_half + j] = 0f0
    zrec[nsrc_half + j] = range(50f0, 50f0, length=nrec)
    xsrc[nsrc_half + j] = xloc
    ysrc[nsrc_half + j] = 0f0
    zsrc[nsrc_half + j] = 20f0
end

# receiver sampling and recording time
timeR = 10000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR)

# Source sampling and number of time steps
timeS = 10000f0
dtS = 2f0

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# Setup wavelet
f0 = 0.015  # dominant frequency in [kHz]
wavelet = ricker_wavelet(timeS,dtS,f0)
q = judiVector(srcGeometry,wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model0)
info = Info(prod(model0.n),nsrc,ntComp)

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

opt = Options(isic=true,    # impedance modeling
              save_data_to_disk=true,
              file_path = file_path,  # directory for saving generated shots
              file_name = file_name
              )

# Setup operators
Pr = judiProjection(info, recGeometry)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*Ps', q)

# Linearized modeling (shots written to disk as SEG-Y files automatically)
J*dm

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
#=
Implements the Kalman filter for a linear Gaussian state space model.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-07-29

References
----------

https://julia.quantecon.org/tools_and_techniques/kalman.html

TODO: Do docstrings here after implementing LinerStateSpace
=#

mutable struct Kalman
    A
    G
    Q
    R
    k
    n
    cur_x_hat
    cur_sigma
end

# Initializes current mean and cov to zeros
function Kalman(A, G, Q, R)
    k = size(G, 1)
    n = size(G, 2)
    xhat = n == 1 ? zero(eltype(A)) : zeros(n)
    Sigma = n == 1 ? zero(eltype(A)) : zeros(n, n)
    return Kalman(A, G, Q, R, k, n, xhat, Sigma)
end

function set_state!(k::Kalman, x_hat, Sigma)
    k.cur_x_hat = x_hat
    k.cur_sigma = Sigma
    Nothing
end

@doc doc"""
Updates the moments (`cur_x_hat`, `cur_sigma`) of the time ``t`` prior to the
time ``t`` filtering distribution, using current measurement ``y_t``.
The updates are according to

```math
    \hat{x}^F = \hat{x} + \Sigma G' (G \Sigma G' + R)^{-1}
                    (y - G \hat{x}) \\

    \Sigma^F = \Sigma - \Sigma G' (G \Sigma G' + R)^{-1} G
               \Sigma
```

#### Arguments

- `k::Kalman` An instance of the Kalman filter
- `y` The current measurement

"""
function prior_to_filtered!(k::Kalman, y)
    # simplify notation
    G, R = k.G, k.R
    x_hat, Sigma = k.cur_x_hat, k.cur_sigma

    # and then update
    if k.k > 1
        reshape(y, k.k, 1)
    end
    A = Sigma * G'
    B = G * Sigma * G' + R
    M = A / B
    k.cur_x_hat = x_hat + M * (y .- G * x_hat)
    k.cur_sigma = Sigma - M * G * Sigma
    Nothing
end

"""
Updates the moments of the time ``t`` filtering distribution to the
moments of the predictive distribution, which becomes the time
``t+1`` prior

#### Arguments
- `k::Kalman` An instance of the Kalman filter

"""
function filtered_to_forecast!(k::Kalman)
    # simplify notation
    A, Q = k.A, k.Q
    x_hat, Sigma = k.cur_x_hat, k.cur_sigma

    # and then update
    k.cur_x_hat = A * x_hat
    k.cur_sigma = A * Sigma * A' + Q
    Nothing
end

"""
Updates `cur_x_hat` and `cur_sigma` given array `y` of length `k`.  The full
update, from one period to the next

#### Arguments

- `k::Kalman` An instance of the Kalman filter
- `y` An array representing the current measurement

"""
function update!(k::Kalman, y)
    prior_to_filtered!(k, y)
    filtered_to_forecast!(k)
    Nothing
end

function stationary_values(k::Kalman)
    # simplify notation

    A, Q, G, R = k.A, k.Q, k.G, k.R

    # solve Riccati equation, obtain Kalman gain
    Sigma_inf = solve_discrete_riccati(A', G', Q, R)
    K_inf = A * Sigma_inf * G' * inv(G * Sigma_inf * G' .+ R)
    return Sigma_inf, K_inf
end

"""
computes log-likelihood of period ``t``

##### Arguments
- `kn::Kalman`: `Kalman` specifying the model. Current values must be the
                forecast for period ``t`` observation conditional on ``t-1``
                observation.
- `y::AbstractVector`: Respondentbservations at period ``t``

##### Returns
- `logL::Real`: log-likelihood of observations at period ``t``
"""
function log_likelihood(k::Kalman, y::AbstractVector)
    eta = y - k.G*k.cur_x_hat # forecast error
    P = k.G*k.cur_sigma*k.G' + k.R # covariance matrix of forecast error
    logL = - (length(y)*log(2pi) + logdet(P) .+ eta'/P*eta)[1]/2
    return logL
end

"""
computes log-likelihood of entire observations

##### Arguments
- `kn::Kalman`: `Kalman` specifying the model. Initial value must be the prior
                for t=1 period observation, i.e. ``x_{1|0}``.
- `y::AbstractMatrix`: `n x T` matrix of observed data.
                       `n` is the number of observed variables in one period.
                       Each column is a vector of observations at each period.

##### Returns
- `logL::Real`: log-likelihood of all observations
"""
function compute_loglikelihood(kn::Kalman, y::AbstractMatrix)
    T = size(y, 2)
    logL = 0
    # forecast and update
    for t in 1:T
        logL = logL + log_likelihood(kn, y[:, t])
        update!(kn, y[:, t])
    end
    return logL
end

"""
##### Arguments
- `kn::Kalman`: `Kalman` specifying the model. Initial value must be the prior
                for t=1 period observation, i.e. ``x_{1|0}``.
- `y::AbstractMatrix`: `n x T` matrix of observed data.
                       `n` is the number of observed variables in one period.
                       Each column is a vector of observations at each period.

##### Returns
- `x_smoothed::AbstractMatrix`: `k x T` matrix of smoothed mean of states.
                                `k` is the number of states.
- `logL::Real`: log-likelihood of all observations
- `sigma_smoothed::AbstractArray` `k x k x T` array of smoothed covariance matrix of states.
"""
function smooth(kn::Kalman, y::AbstractMatrix)
    G, R = kn.G, kn.R

    T = size(y, 2)
    n = kn.n
    x_filtered = Matrix{Float64}(undef, n, T)
    sigma_filtered = Array{Float64}(undef, n, n, T)
    sigma_forecast = Array{Float64}(undef, n, n, T)
    logL = 0
    # forecast and update
    for t in 1:T
        logL = logL + log_likelihood(kn, y[:, t])
        prior_to_filtered!(kn, y[:, t])
        x_filtered[:, t], sigma_filtered[:, :, t] = kn.cur_x_hat, kn.cur_sigma
        filtered_to_forecast!(kn)
        sigma_forecast[:, :, t] = kn.cur_sigma
    end
    # smoothing
    x_smoothed = copy(x_filtered)
    sigma_smoothed = copy(sigma_filtered)
    for t in (T-1):-1:1
        x_smoothed[:, t], sigma_smoothed[:, :, t] =
            go_backward(kn, x_filtered[:, t], sigma_filtered[:, :, t],
                        sigma_forecast[:, :, t], x_smoothed[:, t+1],
                        sigma_smoothed[:, :, t+1])
    end

    return x_smoothed, logL, sigma_smoothed
end

"""
##### Arguments
- `kn::Kalman`: `Kalman` specifying the model.
- `x_fi::Vector`: filtered mean of state for period ``t``
- `sigma_fi::Matrix`: filtered covariance matrix of state for period ``t``
- `sigma_fo::Matrix`: forecast of covariance matrix of state for period ``t+1``
                      conditional on period ``t`` observations
- `x_s1::Vector`: smoothed mean of state for period ``t+1``
- `sigma_s1::Matrix`: smoothed covariance of state for period ``t+1``

##### Returns
- `x_s1::Vector`: smoothed mean of state for period ``t``
- `sigma_s1::Matrix`: smoothed covariance of state for period ``t``
"""
function go_backward(k::Kalman, x_fi::Vector,
                     sigma_fi::Matrix, sigma_fo::Matrix,
                     x_s1::Vector, sigma_s1::Matrix)
    A = k.A
    temp = sigma_fi*A'/sigma_fo
    x_s = x_fi + temp*(x_s1-A*x_fi)
    sigma_s = sigma_fi + temp*(sigma_s1-sigma_fo)*temp'
    return x_s, sigma_s
end

ext: jl


text:
const N_CLEAR_FILES = (0x3f3f3f3f3f3f3f3f, 0x7f7f7f7f7f7f7f7f,
                       0xfefefefefefefefe, 0xfcfcfcfcfcfcfcfc,
                       0xfcfcfcfcfcfcfcfc, 0xfefefefefefefefe,
                       0x7f7f7f7f7f7f7f7f, 0x3f3f3f3f3f3f3f3f)
const N_SHIFTS = (-10, -17, -15, -6, 10, 17, 15, 6)

function knightMovesGen()
    knight_moves = Dict{UInt64,UInt64}()

    for square in values(PGN2UINT)
        targets = EMPTY
        for (clear, shift) in zip(N_CLEAR_FILES, N_SHIFTS)
            if (square & clear) >> shift != EMPTY
                targets |= (square & clear) >> shift
            end
        end
        push!(knight_moves, square=>targets)
    end
    return knight_moves
end
const KNIGHT = knightMovesGen()


text:
@with_kw mutable struct L1normParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    λ::T = 1
end

struct L1norm <: AbstractProx
    params::L1normParameters{Float64}

    function (::Type{L1norm})(; kwargs...)
        return new(L1normParameters(; kwargs...))
    end
end

function prox!(l1norm::L1norm,step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    @unpack λ = l1norm.params
    xval = xprev - step*gcurr
    temp = max.(abs.(xval) .- λ*step, zeros(length(xprev)))
    xcurr .= temp
    xcurr[xval .< 0] *= -1
end


text:
"""
Pinched cylinder with free ends and a concentrated force.
Mostly inextensional deformation.

Asymptotic deformation modes of benchmark
problems suitable for evaluating shell elements
Demetres Briassoulis
Comput. Methods Appl. Mech. Engrg. 194 (2005) 2385–2405
"""
module pinched_cylinder_free_ends_examples

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute_dsg_model(formul, n = 2, visualize = true)
    E = 10.5e6;
    nu = 0.3125
    R = 4.953
    L = 10.35
    P = 100.0

    # [8] D. Briassoulis, Testing the asymptotic behaviour of shell elements.
    # Part I. The classical benchmark tests, Int. J. Numer. Meth. Engrg. 54(3)
    # (2002) 421–452.

    # Analytical: series solution; inextensible [8,17]
    thickness = R/1.0e6
    analyt_sol = -1.907*P*R/E/thickness^3;
    # thickness = R/1.0e4
    # analyt_sol = -1.912*P*R/E/thickness^3;
    # Finite element solution (25 · 35 mesh) [8]
    # thickness = R/1.0e2
    # analyt_sol = -1.997*P*R/E/thickness^3;

    @info "Mesh: $n elements per side"

    # Mesh

    tolerance = R/n/1000
    fens, fes = T3block(90/360*2*pi,L/2,n,n,:b);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a), y, R*cos(a))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    # femm.drilling_stiffness_scale = 1.0e0
    associategeometry! = formul.associategeometry!
    stiffness = formul.stiffness

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # free end
    # l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    # for i in [1,3]
    #     setebc!(dchi, l1, true, i)
    # end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Z
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0.0 0.0], inflate = tolerance)
    for i in [3,4,5]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[0 0 L/2 L/2 -Inf Inf], inflate = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[0, 0, -P/4, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    @show n,  dchi.values[nl, 3]/analyt_sol*100

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (L/8)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return true
end

function test_convergence(formul = FEMMShellT3FFModule)
    @info "Pinched cylinder, formulation=$(formul)"
    for n in [2, 4, 8, 16, 32] # 3:64 #
        _execute_dsg_model(formul, n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
@testset "$(rpad("Lagrange Polynomials Tests",80))" begin

    b = LagrangeBasis([0.0, 1.0])
    l = LagrangePolynomial([0.0, 1.0], ones(2))

    @test l.b == b
    @test l.c == ones(2)

    xp = [0.0, 0.5, 1.0]
    yp = zeros(3)
    yr = [1.0, 0.5, 0.0]

    evaluate!(l, xp, yp)
    @test yp == ones(3)

    y = zeros(2)
    y[1] = 1

    evaluate!(similar(l, y), xp, yp)
    @test yp == yr

    evaluate!(b, y, xp, yp)
    @test yp == yr

    y = zeros(2)
    y[2] = 1

    evaluate!(similar(l, y), xp, yp)
    @test yp == xp

    evaluate!(b, y, xp, yp)
    @test yp == xp

end


text:
function ENDF_interp(nbt, int, x, y, q)
    # https://t2.lanl.gov/nis/endf/intro09.html
    if q >= x[end]
        return y[end]
    elseif q < x[1]
        return y[1]
    end
    idx_1 = findfirst(X -> X > q, x)
    idx_0 = idx_1 - 1
    x_1 = x[idx_1]
    x_0 = x[idx_0]
    y_1 = y[idx_1]
    y_0 = y[idx_0]
    interp_idx = findfirst(X -> X > idx_0, nbt)
    interp = int[interp_idx]
    #println(x_1, y_1, x_0, y_0, interp)
    if interp == 1
        return y_0
    elseif interp == 2
        return y_0 + (y_1 - y_0) * (q - x_0) / (x_1 - x_0)
    elseif interp == 3
        return y_0 + (y_1 - y_0) * log(q/x_0) / log(x_1/x_0)
    elseif interp == 4
        return exp(log(y_0) + log(y_1/y_0) * (q - x_0) / (x_1 - x_0))
    elseif interp == 5
        return exp(log(y_0) + log(y_1/y_0) * log(q/x_0) / log(x_1/x_0))
    elseif interp == 6
        B = log((y_1*x_1)/(y_0*x_0)) / (1/sqrt(x_0-0) - 1/sqrt(x_1-0))
        A = exp(B/sqrt(x_0-0))*y_0*x_0
        return (A/B) * exp(-B/sqrt(q-0))
    end
end

function ENDF_interp_func(NBT, INT, ES, P)
    return E -> ENDF_interp(NBT, INT, ES, P, E)
end

#filename ="C:/Users/warner/Downloads/Lib80x/Lib80x/H/1001.800nc"
filename = "C:/Users/warner/Desktop/My Stuff/Personal Learning/Monte Carlo/Lib80x/Lib80x/C/6012.805nc"
file_lines = readlines(filename)
pos = 1

if length(split(file_lines[1])) == 3
    pos = pos + 2
end
HZ, AW, TZ, HD = split(file_lines[pos], " ", keepempty = false)
HK = file_lines[pos+1][1:70]
HM = file_lines[pos+1][71:80]
pos = pos + 2

IZAW = split(join(file_lines[pos:pos+3]))
IZ = parse.(Int, IZAW[1:2:end])
AW = parse.(Float64, IZAW[2:2:end])
pos = pos + 4

NXS = parse.(Int, split(join(file_lines[pos:pos+1])))
LEN_XSS = NXS[ 1] # Length of second block of data (XSS array)
ZA      = NXS[ 2] # 1000 * Z + A
NES     = NXS[ 3] # Number of energies
NTR     = NXS[ 4] # Number of reactions excluding elastic scattering
NR      = NXS[ 5] # Number of reactions having secondary neutrons excluding elastic scattering
NTRP    = NXS[ 6] # Number of photon production reactions
NTYPE   = NXS[ 7] # Number of particle types for which production data is given
NPCR    = NXS[ 8] # Number of delayed neutron precurser families
S       = NXS[ 9] # Excited state
Z       = NXS[10] # Atomic number
A       = NXS[11] # Atomic mass number
pos = pos + 2

JXS = parse.(Int, split(join(file_lines[pos:pos+3])))
pos = pos + 4
XSS = parse.(Float64, split(join(file_lines[pos:end])))

ESZ     = JXS[ 1] # Energy table
NU      = JXS[ 2] # Fission yield data
MTR     = JXS[ 3] # MT array
LQR     = JXS[ 4] # Q-value array
TYR     = JXS[ 5] # Reaction type array
LSIG    = JXS[ 6] # Table of cross section locators
SIG     = JXS[ 7] # Cross sections
LAND    = JXS[ 8] # Table of angular distribution locators
AND     = JXS[ 9] # Angular distributions
LDLW    = JXS[10] # Table of energy distribution locators
DLW     = JXS[11] # Energy distributions
GPD     = JXS[12] # Photon production data
MTRP    = JXS[13] # Photon production MT array
LSIGP   = JXS[14] # Table of photon production cross section locators
SIGP    = JXS[15] # Photon production cross sections
LANDP   = JXS[16] # Table of photon production angular distribution locators
ANDP    = JXS[17] # Photon production angular distributions
LDLWP   = JXS[18] # Table of photon production energy distribution locators
DLWP    = JXS[19] # Photon production energy distributions
YP      = JXS[20] # Table of yield multipliers
FIS     = JXS[21] # Total fission cross section
END     = JXS[22] # Last word of the conventional table (last word of photon production data)
LUNR    = JXS[23] # Probability tables
DNU     = JXS[24] # Delayed neutron yield data
BDD     = JXS[25] # Basic delayed neutron precursor data (decay constant’s, probabilities)
DNEDL   = JXS[26] # Table of delayed neutron energy distribution locators
DNED    = JXS[27] # Delayed neutron energy distributions
PTYPE   = JXS[30] # Particle type array
NTRO    = JXS[31] # Array containing the number of particle production reactions
NEXT    = JXS[32] # Table of particle production locators (IXS array)

# Read ESZ section
Energy      = XSS[ESZ : ESZ + NES - 1]
sigma_t     = XSS[ESZ + NES : ESZ + 2*NES - 1]
sigma_a     = XSS[ESZ + 2*NES : ESZ + 3*NES - 1]
sigma_el    = XSS[ESZ + 3*NES : ESZ + 4*NES - 1]
H_ave       = XSS[ESZ + 4*NES : ESZ + 5*NES - 1]

# Read NU section
function nu_function(KNU)
    if XSS[KNU] == 1
        NC = convert(Int, XSS[KNU + 1])
        C = XSS[KNU + 2 : KNU + 1 + NC]
        l = collect(1:1:NC)
        return E -> sum(C.*E^(l-1))
    elseif XSS[KNU] == 2
        NR = convert(Int, XSS[KNU + 1])
        NBT = convert.(Int, XSS[KNU + 2 : KNU + 1 + NR])
        INT = convert.(Int, XSS[KNU + 2 + NR : KNU + 1 + 2*NR])
        NE = convert(Int, XSS[KNU + 2 + 2*NR])
        ES = XSS[KNU + 3 + 2*NR : KNU + 2 + 2*NR + NE]
        nubar = XSS[KNU + 3 + 2*NR + NE : KNU + 2 + 2*NR + 2*NE]
        if NR == 0
            NBT = [NE + 1]
            INT = [2]
        end
        return E -> ENDF_interp(NBT, INT, ES, nubar, E)
    end
end

# If XSS[NU] > 1, then there is only one array given that starts at NU[2]
if NU > 0
    if XSS[NU] > 0
        if DNU == 0
            nu_total = nu_function(NU)
        else
            nu_prompt = nu_function(NU)
        end
    # Otherwise, prompt and total arrays are given
    else
        nu_prompt = nu_function(NU + 1)
        nu_total = nu_function(NU + 1 + convert(Int, abs(XSS[NU])))
    end
end
if DNU > 0
    nu_delayed = nu_function(DNU)
end

#=
function dnu_function(start)

end

if BDD > 0
    DNP_decay_constants = []
    DNP_probabilities = Dict()
    idx = BDD
    for group = 1:NPCR
        DEC = XSS[idx]
        NR = convert(Int, XSS[idx + 1])
        NBT = convert.(Int, XSS[idx + 2 : idx + 1 + NR])
        INT = convert.(Int, XSS[idx + 2 + NR : idx + 1 + 2*NR])
        NE = convert(Int, XSS[idx + 2 + 2*NR])
        ES = XSS[idx + 3 + 2*NR : idx + 2 + 2*NR + NE]
        P = XSS[idx + 3 + 2*NR + NE : idx + 2 + 2*NR + 2*NE]
        if NR == 0
            NBT = [NE + 1]
            INT = [2]
        end
        append!(DNP_decay_constants, DEC)
        DNP_probabilities[group] = ENDF_interp_func(NBT, INT, ES, P)
        idx = idx + 3 + 2*NR + 2*NE
    end
end
=#

#=
using PyPlot
loglog(Energy, sigma_el)
=#

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct AllocationRecord <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AllocationRecord(; kwargs...)
        obj = new(meta(AllocationRecord), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AllocationRecord
const __meta_AllocationRecord = Ref{ProtoMeta}()
function meta(::Type{AllocationRecord})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AllocationRecord)
            __meta_AllocationRecord[] = target = ProtoMeta(AllocationRecord)
            allflds = Pair{Symbol,Union{Type,String}}[:alloc_micros => Int64, :alloc_bytes => Int64]
            meta(target, AllocationRecord, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AllocationRecord[]
    end
end
function Base.getproperty(obj::AllocationRecord, name::Symbol)
    if name === :alloc_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :alloc_bytes
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct AllocatorMemoryUsed <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AllocatorMemoryUsed(; kwargs...)
        obj = new(meta(AllocatorMemoryUsed), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AllocatorMemoryUsed
const __meta_AllocatorMemoryUsed = Ref{ProtoMeta}()
function meta(::Type{AllocatorMemoryUsed})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AllocatorMemoryUsed)
            __meta_AllocatorMemoryUsed[] = target = ProtoMeta(AllocatorMemoryUsed)
            fnum = Int[1,2,3,4,6,5]
            allflds = Pair{Symbol,Union{Type,String}}[:allocator_name => AbstractString, :total_bytes => Int64, :peak_bytes => Int64, :live_bytes => Int64, :allocation_records => Base.Vector{AllocationRecord}, :allocator_bytes_in_use => Int64]
            meta(target, AllocatorMemoryUsed, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AllocatorMemoryUsed[]
    end
end
function Base.getproperty(obj::AllocatorMemoryUsed, name::Symbol)
    if name === :allocator_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :total_bytes
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :peak_bytes
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :live_bytes
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :allocation_records
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AllocationRecord}
    elseif name === :allocator_bytes_in_use
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct NodeOutput <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NodeOutput(; kwargs...)
        obj = new(meta(NodeOutput), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NodeOutput
const __meta_NodeOutput = Ref{ProtoMeta}()
function meta(::Type{NodeOutput})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NodeOutput)
            __meta_NodeOutput[] = target = ProtoMeta(NodeOutput)
            fnum = Int[1,3]
            allflds = Pair{Symbol,Union{Type,String}}[:slot => Int32, :tensor_description => TensorDescription]
            meta(target, NodeOutput, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NodeOutput[]
    end
end
function Base.getproperty(obj::NodeOutput, name::Symbol)
    if name === :slot
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :tensor_description
        return (obj.__protobuf_jl_internal_values[name])::TensorDescription
    else
        getfield(obj, name)
    end
end

mutable struct MemoryStats <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function MemoryStats(; kwargs...)
        obj = new(meta(MemoryStats), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct MemoryStats
const __meta_MemoryStats = Ref{ProtoMeta}()
function meta(::Type{MemoryStats})
    ProtoBuf.metalock() do
        if !isassigned(__meta_MemoryStats)
            __meta_MemoryStats[] = target = ProtoMeta(MemoryStats)
            fnum = Int[1,3,5,2,4,6]
            pack = Symbol[:persistent_tensor_alloc_ids,:device_persistent_tensor_alloc_ids]
            allflds = Pair{Symbol,Union{Type,String}}[:temp_memory_size => Int64, :persistent_memory_size => Int64, :persistent_tensor_alloc_ids => Base.Vector{Int64}, :device_temp_memory_size => Int64, :device_persistent_memory_size => Int64, :device_persistent_tensor_alloc_ids => Base.Vector{Int64}]
            meta(target, MemoryStats, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_MemoryStats[]
    end
end
function Base.getproperty(obj::MemoryStats, name::Symbol)
    if name === :temp_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :persistent_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :persistent_tensor_alloc_ids
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :device_temp_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :device_persistent_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :device_persistent_tensor_alloc_ids
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct NodeExecStats <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NodeExecStats(; kwargs...)
        obj = new(meta(NodeExecStats), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NodeExecStats
const __meta_NodeExecStats = Ref{ProtoMeta}()
function meta(::Type{NodeExecStats})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NodeExecStats)
            __meta_NodeExecStats[] = target = ProtoMeta(NodeExecStats)
            allflds = Pair{Symbol,Union{Type,String}}[:node_name => AbstractString, :all_start_micros => Int64, :op_start_rel_micros => Int64, :op_end_rel_micros => Int64, :all_end_rel_micros => Int64, :memory => Base.Vector{AllocatorMemoryUsed}, :output => Base.Vector{NodeOutput}, :timeline_label => AbstractString, :scheduled_micros => Int64, :thread_id => UInt32, :referenced_tensor => Base.Vector{AllocationDescription}, :memory_stats => MemoryStats, :all_start_nanos => Int64, :op_start_rel_nanos => Int64, :op_end_rel_nanos => Int64, :all_end_rel_nanos => Int64, :scheduled_nanos => Int64]
            meta(target, NodeExecStats, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NodeExecStats[]
    end
end
function Base.getproperty(obj::NodeExecStats, name::Symbol)
    if name === :node_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :all_start_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :op_start_rel_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :op_end_rel_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :all_end_rel_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :memory
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AllocatorMemoryUsed}
    elseif name === :output
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{NodeOutput}
    elseif name === :timeline_label
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :scheduled_micros
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :thread_id
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :referenced_tensor
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AllocationDescription}
    elseif name === :memory_stats
        return (obj.__protobuf_jl_internal_values[name])::MemoryStats
    elseif name === :all_start_nanos
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :op_start_rel_nanos
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :op_end_rel_nanos
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :all_end_rel_nanos
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :scheduled_nanos
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct DeviceStepStats_ThreadNamesEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function DeviceStepStats_ThreadNamesEntry(; kwargs...)
        obj = new(meta(DeviceStepStats_ThreadNamesEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct DeviceStepStats_ThreadNamesEntry (mapentry)
const __meta_DeviceStepStats_ThreadNamesEntry = Ref{ProtoMeta}()
function meta(::Type{DeviceStepStats_ThreadNamesEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DeviceStepStats_ThreadNamesEntry)
            __meta_DeviceStepStats_ThreadNamesEntry[] = target = ProtoMeta(DeviceStepStats_ThreadNamesEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => UInt32, :value => AbstractString]
            meta(target, DeviceStepStats_ThreadNamesEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_DeviceStepStats_ThreadNamesEntry[]
    end
end
function Base.getproperty(obj::DeviceStepStats_ThreadNamesEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct DeviceStepStats <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function DeviceStepStats(; kwargs...)
        obj = new(meta(DeviceStepStats), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct DeviceStepStats
const __meta_DeviceStepStats = Ref{ProtoMeta}()
function meta(::Type{DeviceStepStats})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DeviceStepStats)
            __meta_DeviceStepStats[] = target = ProtoMeta(DeviceStepStats)
            allflds = Pair{Symbol,Union{Type,String}}[:device => AbstractString, :node_stats => Base.Vector{NodeExecStats}, :thread_names => Base.Dict{UInt32,AbstractString}]
            meta(target, DeviceStepStats, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_DeviceStepStats[]
    end
end
function Base.getproperty(obj::DeviceStepStats, name::Symbol)
    if name === :device
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :node_stats
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{NodeExecStats}
    elseif name === :thread_names
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{UInt32,AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct StepStats <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StepStats(; kwargs...)
        obj = new(meta(StepStats), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct StepStats
const __meta_StepStats = Ref{ProtoMeta}()
function meta(::Type{StepStats})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StepStats)
            __meta_StepStats[] = target = ProtoMeta(StepStats)
            allflds = Pair{Symbol,Union{Type,String}}[:dev_stats => Base.Vector{DeviceStepStats}]
            meta(target, StepStats, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StepStats[]
    end
end
function Base.getproperty(obj::StepStats, name::Symbol)
    if name === :dev_stats
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{DeviceStepStats}
    else
        getfield(obj, name)
    end
end

export AllocationRecord, AllocatorMemoryUsed, NodeOutput, MemoryStats, NodeExecStats, DeviceStepStats_ThreadNamesEntry, DeviceStepStats, StepStats
# mapentries: "DeviceStepStats_ThreadNamesEntry" => ("UInt32", "AbstractString")


text:
# convert Char {N,V} to cusolverEigMode_t
function cusolverjob(jobz::Char)
    if jobz == 'N'
        return CUSOLVER_EIG_MODE_NOVECTOR
    end
    if jobz == 'V'
        return CUSOLVER_EIG_MODE_VECTOR
    end
    throw(ArgumentError("unknown cusolver eigmode $jobz."))
end


text:
using GLPK
const GLPKOPT = with_optimizer(GLPK.Optimizer)

@testset "GLPK Linear Problems" begin
    list_of_linear_problems = [
        lp1_test,
        lp2_test,
        lp3_test,
        lp4_test,
        lp5_test,
        lp6_test,
        lp7_test,
        # lp8_test, Int64 problem, does not work
        # lp9_test, Interval is not implemented
        lp10_test,
        # lp11_test, Feasibility not supported
        lp12_test,
        lp13_test
    ]
    test_strong_duality(list_of_linear_problems, GLPKOPT)
end


text:
function POMDPLinter.requirements_info(solver::AbstractPOMCPSolver, problem::POMDP)
    println("""
    Since POMCP is an online solver, most of the computation occurs in `action(planner, state)`. In order to view the requirements for this function, please, supply an initial beleif to `requirements_info`, e.g.

            @requirements_info $(typeof(solver))() $(typeof(problem))() initialstate(pomdp)

        """)
end

function POMDPLinter.requirements_info(solver::AbstractPOMCPSolver, problem::POMDP, b)
    policy = solve(solver, problem)
    requirements_info(policy, b)
end

POMDPs.requirements_info(policy::POMCPPlanner, b) = @show_requirements action(policy, b)

@POMDP_require action(p::POMCPPlanner, b) begin
    tree = POMCPTree(p.problem, b, p.solver.tree_queries)
    @subreq search(p, b, tree)
end

@POMDP_require search(p::POMCPPlanner, b, t::POMCPTree) begin
    P = typeof(p.problem)
    @req rand(::typeof(p.rng), ::typeof(b))
    s = rand(p.rng, b)
    @req isterminal(::P, ::statetype(P))
    @subreq simulate(p, s, POMCPObsNode(t, 1), p.solver.max_depth)
end

@POMDP_require simulate(p::POMCPPlanner, s, hnode::POMCPObsNode, steps::Int) begin
    P = typeof(p.problem)
    S = statetype(P)
    A = actiontype(P)
    O = obstype(P)
    @req gen(::P, ::S, ::A, ::typeof(p.rng))
    @req isequal(::O, ::O)
    @req hash(::O)
    # from insert_obs_node!
    @req actions(::P)
    AS = typeof(actions(p.problem))
    @req length(::AS)
    @subreq estimate_value(p.solved_estimator, p.problem, s, hnode, steps)
    @req discount(::P)
end

@POMDP_require estimate_value(f::Function, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) begin
    @req f(::typeof(pomdp), ::typeof(start_state), ::typeof(h), ::typeof(steps))
end


text:
include("kitchen.jl")
function leafset(t)
    if isa(t, Int)
        return Multiset([t])
    else
        return reduce(+, map(leafset, t))
    end
end

using DocOpt

doc = """validate.jl

Usage:
    validate.jl <trees>

Options:
    -h --help   Show this screen.
"""
args = docopt(doc)
inpath = args["<trees>"]

newicks = eval.(Meta.parse.(readlines(inpath)))
@show issingly.(leafset.(newicks))
function derr(l, r)
    "$l has a clade that is smaller than $(r)!" |> println
    "  by $(leafset(l)) being smaller than $(leafset(r))" |> println
end

for (l, r) in combinations(newicks, 2)
    if leafset(l) < leafset(r)
        derr(l, r)
    end

    if leafset(l) > leafset(r)
        derr(r, l)
    end
end

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
using Base.Test
using Stella
@test 1 == 1


text:
using Test
using AlgamesDriving
using Algames
using BenchmarkTools
using Colors
using CoordinateTransformations
using FileIO
using GeometryBasics
using LinearAlgebra
using MeshCat
using MeshIO
using Parameters
using REPL
using Random
using Rotations
using StaticArrays

#Struct
include("struct/vehicle_state.jl")
include("struct/starting_area.jl")
include("struct/lane.jl")
include("struct/roadway.jl")
include("struct/player.jl")

# Scenario
include("scenario/scenario.jl")

# Visualization
include("visualization/camera.jl")

# MPC
include("mpc/mpc_struct.jl")
include("mpc/mpc_methods.jl")


text:
export CuBLAS

module CuBLAS
using ..CUDA

# cublasStatus_t
const CUBLAS_STATUS_SUCCESS         = 0
const CUBLAS_STATUS_NOT_INITIALIZED = 1
const CUBLAS_STATUS_ALLOC_FAILED    = 3
const CUBLAS_STATUS_INVALID_VALUE   = 7
const CUBLAS_STATUS_ARCH_MISMATCH   = 8
const CUBLAS_STATUS_MAPPING_ERROR   = 11
const CUBLAS_STATUS_EXECUTION_FAILED= 13
const CUBLAS_STATUS_INTERNAL_ERROR  = 14
const CUBLAS_STATUS_NOT_SUPPORTED   = 15
const CUBLAS_STATUS_LICENSE_ERROR   = 16

immutable CuBLASError <: Exception
  code :: Int
end
using Compat
const cublas_error_description = @compat(Dict(
  CUBLAS_STATUS_SUCCESS         => "Success",
  CUBLAS_STATUS_NOT_INITIALIZED => "Not initialized",
  CUBLAS_STATUS_ALLOC_FAILED    => "Alloc failed",
  CUBLAS_STATUS_INVALID_VALUE   => "Invalid value",
  CUBLAS_STATUS_ARCH_MISMATCH   => "Arch mismatch",
  CUBLAS_STATUS_MAPPING_ERROR   => "Mapping error",
  CUBLAS_STATUS_EXECUTION_FAILED=> "Execution failed",
  CUBLAS_STATUS_INTERNAL_ERROR  => "Internal error",
  CUBLAS_STATUS_NOT_SUPPORTED   => "Not supported",
  CUBLAS_STATUS_LICENSE_ERROR   => "License error"
))
import Base.show
show(io::IO, error::CuBLASError) = print(io, cublas_error_description[error.code])

if is_windows()
  const libcublas = Libdl.find_library(["cublas64_80.dll", "cublas64_70.dll", "cublas64_65.dll",
      "cublas32_70.dll", "cublas32_65.dll", "cublas64_75.dll"], [""])
  @assert (libcublas != "") "Could not find cuBLAS DLL [cublas64_70.dll, cublas64_65.dll, cublas32_70.dll, cublas32_65.dll, cublas64_75.dll]. See: http://mochajl.readthedocs.io/en/latest/user-guide/backend.html#cuda-backend"
else
  const libcublas = Libdl.find_library(["libcublas"], [""])
  @assert (libcublas != "") "Could not find cuBLAS DLL [libcublas]. See http://mochajl.readthedocs.io/en/latest/user-guide/backend.html#cuda-backend"
end

macro cublascall(fv, argtypes, args...)
  args = map(esc, args)
  f = eval(fv)
  quote
    _curet = ccall( ($(Meta.quot(f)), $libcublas), Cint, $argtypes, $(args...)  )
    if round(Int, _curet) != CUBLAS_STATUS_SUCCESS
      throw(CuBLASError(round(Int, _curet)))
    end
  end
end

const Handle = Ptr{Void}
const StreamHandle = Ptr{Void}

function create()
  handle = Handle[0]
  @cublascall(:cublasCreate_v2, (Ptr{Handle},), handle)
  return handle[1]
end
function destroy(handle :: Handle)
  @cublascall(:cublasDestroy_v2, (Handle,), handle)
end
function set_stream(handle::Handle, stream::StreamHandle)
  @cublascall(:cublasSetStream_v2, (Handle, StreamHandle), handle, stream)
end
function get_stream(handle::Handle)
  s_handle = StreamHandle[0]
  @cublascall(:cublasGetStream_v2, (Handle, Ptr{StreamHandle}), handle, s_handle)
  return s_handle[1]
end

############################################################
# Copy a vector from host to device
############################################################
function set_vector(n::Int, elem_size::Int, src::Ptr{Void}, incx::Int, dest::Ptr{Void}, incy::Int)
  @cublascall(:cublasSetVector, (Cint, Cint, Ptr{Void}, Cint, Ptr{Void}, Cint),
      n, elem_size, src, incx, dest, incy)
end
function set_vector(n::Int, elem_size::Int, src::Ptr{Void}, incx::Int, dest::CuPtr, incy::Int)
  set_vector(n, elem_size, src, incx, Base.unsafe_convert(Ptr{Void}, dest.p), incy)
end
function set_vector{T}(src::Array{T}, incx::Int, dest::CuPtr, incy::Int)
  elem_size = sizeof(T)
  n = length(src)
  src_buf = convert(Ptr{Void}, pointer(src))
  set_vector(n, elem_size, src_buf, incx, dest, incy)
end
set_vector{T}(src::Array{T}, dest::CuPtr) = set_vector(src, 1, dest, 1)

############################################################
# Copy a vector from device to host
############################################################
function get_vector(n::Int, elem_size::Int, src::CuPtr, incx::Int, dest::Ptr{Void}, incy::Int)
  @cublascall(:cublasGetVector, (Cint, Cint, Ptr{Void}, Cint, Ptr{Void}, Cint),
      n, elem_size, Base.unsafe_convert(Ptr{Void}, src.p), incx, dest, incy)
end
function get_vector{T}(src::CuPtr, incx::Int, dest::Array{T}, incy::Int)
  elem_size = sizeof(T)
  n = length(dest)
  dest_buf = convert(Ptr{Void}, pointer(dest))
  get_vector(n, elem_size, src, incx, dest_buf, incy)
end
get_vector{T}(src::CuPtr, dest::Array{T}) = get_vector(src, 1, dest, 1)

############################################################
# y = α y
############################################################
for (fname, elty) in ((:cublasSscal_v2, :Float32),
                      (:cublasDscal_v2, :Float64))
  @eval begin
    function scal(handle::Handle, n::Int, alpha::$elty, x, incx::Int)
      x = Base.unsafe_convert(Ptr{Void}, x)
      alpha_box = $elty[alpha]
      @cublascall($(string(fname)), (Handle, Cint, Ptr{Void}, Ptr{Void}, Cint),
                  handle, n, alpha_box, x, incx)
    end
    function scal(handle::Handle, n::Int, alpha::$elty, x::CuPtr, incx::Int)
      scal(handle, n, alpha, x.p, incx)
    end
  end
end

############################################################
# y = α x + y
############################################################
for (fname, elty) in ((:cublasSaxpy_v2, :Float32),
                      (:cublasDaxpy_v2, :Float64))
  @eval begin
    function axpy(handle::Handle, n::Int, alpha::$elty, x, incx::Int, y, incy::Int)
      x = Base.unsafe_convert(Ptr{Void}, x)
      y = Base.unsafe_convert(Ptr{Void}, y)
      alpha_box = $elty[alpha]
      @cublascall($(string(fname)), (Handle, Cint, Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}, Cint),
                  handle, n, alpha_box, x, incx, y, incy)
    end
    function axpy(handle::Handle, n::Int, alpha::$elty, x::CuPtr, incx::Int, y::CuPtr, incy::Int)
      axpy(handle, n, alpha, x.p, incx, y.p, incy)
    end
  end
end

############################################################
# vector dot product
############################################################
for (fname, elty) in ((:cublasSdot_v2, :Float32),
                      (:cublasDdot_v2, :Float64))
  @eval begin
    function dot(handle::Handle, ::Type{$elty}, n::Int, x::CuPtr, incx::Int, y::CuPtr, incy::Int)
      result = $elty[0]
      @cublascall($(string(fname)), (Handle, Cint, Ptr{Void}, Cint, Ptr{Void}, Cint, Ptr{Void}),
                  handle, n, x.p, incx, y.p, incy, result)
      return result[1]
    end
  end
end

############################################################
# blas copy
# Note blascopy is copying from x to y, while most of the
# copy! functions in julia and also for blobs are copying
# from y to x.
############################################################
for (fname, elty) in ((:cublasScopy_v2, :Float32),
                      (:cublasDcopy_v2, :Float64))
  @eval begin
    function copy(handle::Handle, ::Type{$elty}, n::Int, x, incx::Int, y, incy::Int)
      x = Base.unsafe_convert(Ptr{Void}, (x))
      y = Base.unsafe_convert(Ptr{Void}, (y))
      @cublascall($(string(fname)), (Handle, Cint, Ptr{Void}, Cint, Ptr{Void}, Cint),
                  handle, n, x, incx, y, incy)
    end
  end
end

############################################################
# cublasOperation_t
############################################################
const OP_N=0
const OP_T=1
const OP_C=2

############################################################
# C = α A * B + β C
############################################################
function gemm{T}(handle::Handle, trans_a::Int, trans_b::Int, m::Int, n::Int, k::Int,
    alpha::T, A::CuPtr, lda::Int, B::CuPtr, ldb::Int, beta::T, C::CuPtr, ldc::Int)
  @assert OP_N <= trans_a <= OP_C
  @assert OP_N <= trans_b <+ OP_C
  alpha_box = T[alpha]
  beta_box = T[beta]
  gemm_impl(handle, trans_a, trans_b, m, n, k, alpha_box, A, lda, B, ldb, beta_box, C, ldc)
end

for (fname, elty) in ((:cublasSgemm_v2, :Float32),
                      (:cublasDgemm_v2, :Float64))
  @eval begin
    function gemm_impl(handle::Handle, trans_a::Int, trans_b::Int, m::Int, n::Int, k::Int,
        alpha_box::Array{$elty}, A::CuPtr, lda::Int, B::CuPtr, ldb::Int, beta_box::Array{$elty}, C::CuPtr, ldc::Int)
      @cublascall($(string(fname)), (Handle, Cint,Cint, Cint,Cint,Cint, Ptr{Void},
                  Ptr{Void},Cint, Ptr{Void},Cint, Ptr{Void}, Ptr{Void},Cint),
                  handle, trans_a, trans_b, m, n, k, alpha_box, A.p, lda, B.p, ldb, beta_box, C.p, ldc)
    end
  end
end

end # module


text:
# test boundary forces
@doc """
Advection Equation -- test_adjoint

The function tests for the correctness functional
computation and then tests if the adjoint vector is being computed correctly.
This is a serial test and uses the input file called

`input_vals_functional_DG.jl`

"""->
function test_adjoint()
  #= TODO: Uncomment when CG meshes support mesh.bndry_geo_nums

  @testset "--- Testing Boundary Functional Computation on CG Mesh ---" begin
    clean_dict(arg_dict)
    fname = "input_vals_functional_CG.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    println("use_DG = ", arg_dict["use_DG"])

    @test ( mesh.isDG )== false
    @test ( opts["functional_name1"] )== "qflux"
    @test isapprox( opts["analytical_functional_val"], 2*(exp(1) - 1)) atol=1e-12
    @test ( opts["geom_edges_functional1"] )== [2,3]

    fname = "./functional_error1.dat"
    error = readdlm(fname)

    @test isapprox( error[1], 0.0060826244541961885) atol=1e-6

    rm("./functional_error1.dat") # Delete the file
  end
  =#

  @testset "--- Testing Boundary Functional & Adjoint Computation On DG Mesh ---" begin

    fname = "input_vals_functional_DG.jl"
    mesh, sbp, eqn, opts, pmesh = createObjects(fname)

    @assert mesh.isDG == true
    @assert opts["jac_method"] == 2
    @assert opts["run_type"] == 5

    functional = createFunctional(mesh, sbp, eqn,
                                  opts, opts["num_functionals"])

    @testset "Checking Functional Object Creation" begin

      @test ( functional.bcnums )== [2,3]
      @test ( functional.val )== zero(Complex{Float64})
      @test ( functional.target_qflux )== zero(Complex{Float64})

    end # End context("Checking Functional Object Creation")

    solvePDE(mesh, sbp, eqn, opts, pmesh)
    evalFunctional(mesh, sbp, eqn, opts, functional)

    @testset "Checking Functional Computation" begin

      analytical_val = 3.0
      functional_error = norm(real(functional.val) - analytical_val,2)
      @test isapprox( functional_error, 0.0) atol=1e-12

      # test another functional
      func = AdvectionEquationMod.IntegralQDataConstructor(Complex128, mesh, sbp, eqn, opts, opts["functional_bcs1"])
      AdvectionEquationMod.calcBndryFunctional(mesh, sbp, eqn, opts, func)
      @test isapprox( func.val, analytical_val) atol=1e-13

    end # End context("Checking Functional Computation")

    @testset "Checking Adjoint Computation on DG mesh" begin

      adjoint_vec = zeros(Complex{Float64}, mesh.numDof)
      pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
      ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
      calcAdjoint(mesh, sbp, eqn, opts, ls, functional, adjoint_vec, recalc_jac=true, recalc_pc=true)

      for i = 1:length(adjoint_vec)
        @test isapprox( real(adjoint_vec[i]), 1.0) atol=1e-10
      end

    end # End context("Checking Adjoint Computation on DG mesh")

  end # End testset("--- Testing Boundary Functional Computation On DG Mesh ---")

end # End function test_adjoint

add_func1!(AdvectionTests, test_adjoint, [TAG_ADJOINT, TAG_SHORTTEST])


text:
using DisplayStructure
const DS = DisplayStructure
using Terming
const T = Terming
using Test

@testset "DisplayStructure.jl" begin

    T.set_term!(T.PseudoTerminal(
        Base.BufferStream(), Base.BufferStream(), Base.BufferStream()
    ))

    TERM_SIZE = T.displaysize()

    include("util.jl")
    include("display_row.jl")
    include("display_array.jl")
    include("feature/feature.jl")

end


text:
module AndTraits

include("andtrait.jl")
include("tag.jl")

end # module


text:
###################################
## Demo script for Bayesian Dimensionality Reduction

using PyPlot, PyCall
@pyimport sklearn.datasets as datasets

push!(LOAD_PATH,".")
import DimensionalityReduction

function load_facedata(skip::Int)
    face = datasets.fetch_olivetti_faces()
    Y_raw = face["images"]
    N, S_raw, _ = size(Y_raw)

    L = round(Int, S_raw / skip)
    Y_tmp = Y_raw[:,1:skip:end, 1:skip:end]
    Y = convert(Array{Float64, 2}, reshape(Y_tmp, N, size(Y_tmp,2)*size(Y_tmp,3))')
    D = size(Y, 1)

    return Y, D, L
end

function visualize(Y::Array{Float64,2}, L::Int)
    D, N = size(Y)
    base = round(Int, sqrt(N))
    v = round(Int, (L*ceil(N / base)))
    h = L * base
    pic = zeros(v, h)

    for n in 1 : N
        i = round(Int, (L*ceil(n / base)))
        idx1 = i - L + 1 : i
        idx2 = L*mod(n-1, base)+1 : L*(mod(n-1, base) + 1)
        pic[idx1,idx2] = reshape(Y[:,n], L, L)
    end
    imshow(pic, cmap=ColorMap("gray"))
end

function visualize(Y::Array{Float64,2}, L::Int, mask::BitArray{2})
    # for missing
    D, N = size(Y)
    base = round(Int, sqrt(N))
    v = round(Int, (L*ceil(N / base)))
    h = L * base
    pic = zeros(v, h, 3)

    Y_3dim = zeros(D, N, 3)
    for i in 1 : 3
        if i == 2
            Y_tmp = deepcopy(Y)
            Y_tmp[mask] = 1
            Y_3dim[:,:,i] = Y_tmp
        else
            Y_tmp = deepcopy(Y)
            Y_tmp[mask] = 0
            Y_3dim[:,:,i] = Y_tmp
        end
    end

    for n in 1 : N
        i = round(Int, (L*ceil(n / base)))
        idx1 = i - L + 1 : i
        idx2 = L*mod(n-1, base)+1 : L*(mod(n-1, base) + 1)
        for i in 1 : 3
            pic[idx1,idx2,i] = reshape(Y_3dim[:,n,i], L, L)
        end
    end
    imshow(pic, cmap=ColorMap("gray"))
end

"""
Run a demo script of missing data interpolation for face dataset.
"""
function test_face_missing()
    # load data
    skip = 2
    Y, D, L = load_facedata(skip)

    # mask
    missing_rate = 0.50
    mask = rand(size(Y)) .< missing_rate
    Y_obs = deepcopy(Y)
    Y_obs[mask] = NaN

    # known parames
    M = 16
    sigma2_y = 0.001
    Sigma_W = zeros(M,M,D)
    Sigma_mu = 1.0 * eye(D)
    for d in 1 : D
        Sigma_W[:,:,d] = 0.1 * eye(M)
    end
    prior = DimensionalityReduction.DRModel(D, M, sigma2_y, zeros(M, D), Sigma_W, zeros(D), Sigma_mu)

    # learn & generate
    max_iter = 100
    posterior, X_est = DimensionalityReduction.VI(deepcopy(Y_obs), prior, max_iter)
    Y_est = posterior.m_W'*X_est + repmat(posterior.m_mu, 1, size(X_est, 2))
    Y_itp = deepcopy(Y_obs)
    Y_itp[mask] = Y_est[mask]

    #visualize
    N_show = 4^2

    figure("Observation")
    clf()
    visualize(Y_obs[:,1:N_show], L, mask[:,1:N_show])
    title("Observation")

    #figure("Estimation")
    #clf()
    #visualize(Y_est[:,1:N_show], L)
    #title("Estimation")

    figure("Interpolation")
    clf()
    visualize(Y_itp[:,1:N_show], L)
    title("Interpolation")

    figure("Truth")
    clf()
    visualize(Y[:,1:N_show], L)
    title("Truth")
    show()
end

"""
Run a dimensionality reduction demo using Iris dataset.
"""
function test_iris()
    ##################
    # load data
    iris = datasets.load_iris()
    Y_obs = iris["data"]'
    label_list = [iris["target_names"][elem+1] for elem in iris["target"]]
    D, N = size(Y_obs)

    ##################
    # 2D compression

    # model
    M = 2
    sigma2_y = 0.001
    Sigma_W = zeros(M,M,D)
    Sigma_mu = 1.0 * eye(D)
    for d in 1 : D
        Sigma_W[:,:,d] = 0.1 * eye(M)
    end
    prior = DimensionalityReduction.DRModel(D, M, sigma2_y, zeros(M, D), Sigma_W, zeros(D), Sigma_mu)

    # learn & generate
    max_iter = 100
    posterior, X_est = DimensionalityReduction.VI(deepcopy(Y_obs), prior, max_iter)

    # visualize
    figure("2D plot")
    clf()
    scatter(X_est[1,1:50], X_est[2,1:50], color="r")
    scatter(X_est[1,51:100], X_est[2,51:100], color="g")
    scatter(X_est[1,101:end], X_est[2,101:end], color="b")
    xlabel("\$x_1\$", fontsize=20)
    ylabel("\$x_2\$", fontsize=20)
    legend([label_list[1], label_list[51], label_list[101]], fontsize=16)

    ##################
    # 3D compression

    # model
    M = 3
    sigma2_y = 0.001
    Sigma_W = zeros(M,M,D)
    Sigma_mu = 1.0 * eye(D)
    for d in 1 : D
        Sigma_W[:,:,d] = 0.1 * eye(M)
    end
    prior = DimensionalityReduction.DRModel(D, M, sigma2_y, zeros(M, D), Sigma_W, zeros(D), Sigma_mu)

    # learn & generate
    max_iter = 100
    posterior, X_est = DimensionalityReduction.VI(deepcopy(Y_obs), prior, max_iter)

    # visualize
    figure("3D plot")
    clf()
    scatter3D(X_est[1,1:50], X_est[2,1:50], X_est[3,1:50], c="r")
    scatter3D(X_est[1,51:100], X_est[2,51:100], X_est[3,51:100], c="g")
    scatter3D(X_est[1,101:end], X_est[2,101:end], X_est[3,101:end], c="b")
    legend([label_list[1], label_list[51], label_list[101]], fontsize=16)
    xlabel("\$x_1\$", fontsize=20)
    ylabel("\$x_2\$", fontsize=20)
    zlabel("\$x_3\$", fontsize=20)
    show()
end

#test_face_missing()
test_iris()


text:
module GPBounding

using GaussianProcesses
using LinearAlgebra
using Random
using Distributions

include("squared_exponential.jl")

export compute_μ_bounds_bnb, compute_σ_ub_bounds, compute_σ_ub_bounds_auto, compute_σ_ub_bounds_approx, compute_σ_ub_bounds_from_gp

function compute_μ_bounds_bnb(gp, x_L, x_U; max_iterations=100, bound_epsilon=1e-2, max_flag=false)
    # By default, it calculates bounds on the minimum.
    # Calculate this vector outside of the loop to save computation.
    theta_vec_train_squared = zeros(gp.nobs);
    theta_vec = ones(gp.dim) * 1 ./ (2*gp.kernel.ℓ2)
    for i = 1:gp.nobs
        theta_vec_train_squared[i] = transpose(theta_vec) * (gp.x[:, i].^2)
    end

    x_best, lbest, ubest = compute_μ_lower_bound(gp, x_L, x_U, theta_vec_train_squared, upper_flag=max_flag)
    if max_flag
        temp = lbest
        lbest = -ubest
        ubest = -temp
    end

    candidates = [[(x_L, x_U), lbest, ubest]]
    iterations = 0
    lb_list = []

    while !isempty(candidates) && iterations < max_iterations
        new_candidates = []
        for candidate in candidates

            extent = candidate[1]
            lb_can = candidate[2]
            ub_can = candidate[3]
            bound_pairs = split_region(extent[1], extent[2])

            for pair in bound_pairs
                x_lb1, lb1, ub1 = compute_μ_lower_bound(gp, pair[1], pair[2], theta_vec_train_squared, upper_flag=max_flag)

                if max_flag
                    temp = lb1
                    lb1 = -ub1
                    ub1 = -temp
                end

                if ub1 <= ubest
                    ubest = ub1
                    lbest = lb1
                    x_best = x_lb1
                    push!(new_candidates, hcat([pair, lb1, ub1]))
                elseif lb1 < ubest
                    push!(new_candidates, hcat([pair, lb1, ub1]))
                end

            end

        end

        if norm(ubest - lbest) < bound_epsilon
            if max_flag
                temp = lbest
                lbest = -ubest
                ubest = -temp
            end
            return x_best, lbest, ubest
        end
        candidates = new_candidates
        iterations += 1
    end
    if max_flag
        temp = lbest
        lbest = -ubest
        ubest = -temp
    end
    return x_best, lbest, ubest
end

function compute_σ_ub_bounds(gp, K_inv, x_L, x_U; max_iterations=10, bound_epsilon=1e-4)

    lbest = -Inf
    ubest = 1.
    x_best = nothing

    candidates = [[(x_L, x_U), lbest, ubest]]
    iterations = 0

    while !isempty(candidates) && iterations < max_iterations
        new_candidates = []
        for candidate in candidates

            extent = candidate[1]
            lb_can = candidate[2]
            ub_can = candidate[3]

            if ub_can < lbest
                continue
            end

            round_lbest = lbest
            round_ubest = ubest

            bound_pairs = split_region(extent[1], extent[2])

            for pair in bound_pairs
                x_ub1, lb1, ub1 = compute_σ_upper_bound(gp, pair[1], pair[2], K_inv)

                if lb1 > lbest
                    lbest = lb1
                    ubest = ub1
                    x_best = x_ub1
                    push!(new_candidates, hcat([pair, lb1, ub1]))
                elseif ub1 < ubest && ub1 > lbest
                    ubest = ub1
                    push!(new_candidates, hcat([pair, lb1, ub1]))
                elseif lb1 > lb_can && ub1 < ub_can && ub1 > lbest
                    push!(new_candidates, hcat([pair, lb1, ub1]))
                end

                if norm(ubest - lbest) < bound_epsilon
                    @debug ubest, lbest
                    return x_best, lbest, ubest
                end
            end

        end

        if (length(new_candidates) == 0 && norm(ubest - lbest) > bound_epsilon)
            delta_x = 10. ^(-iterations)
            if delta_x < 1e-6
                break
            end
            x_L_n = [x_best[i] - delta_x for i=1:length(x_best)]
            x_U_n = [x_best[i] + delta_x for i=1:length(x_best)]
            new_candidates = [[(x_L_n, x_U_n), lbest, ubest]]
        end
        candidates = new_candidates
        iterations += 1
    end

    while norm(ubest - lbest) > bound_epsilon
        delta_x = 10. ^(-iterations)
        if delta_x < 1e-6
            break
        end
        x_L_n = [x_best[i] - delta_x for i=1:length(x_best)]
        x_U_n = [x_best[i] + delta_x for i=1:length(x_best)]
        x_best, lbest, ubest = compute_σ_upper_bound(gp, x_L_n, x_U_n, K_inv)
        iterations += 1
    end

    return x_best, lbest, ubest

end

function compute_σ_ub_bounds_approx(gp, x_L, x_U)
    N = 100
    mt = MersenneTwister(11)
    σ2_best = -Inf
    # Get N samples uniformly dist.
    x_samp = zeros(2,1)
    # @info x_samp
    for i=1:N
        # x_samp = vcat([rand(mt, Uniform(x_L[1], x_U[1]), 1, 1)[1], rand(mt, Uniform(x_L[2], x_U[1]), 1, 1)[1]])
        x_samp[1] = rand(mt, Uniform(x_L[1], x_U[1]), 1, 1)[1]
        x_samp[2] = rand(mt, Uniform(x_L[2], x_U[2]), 1, 1)[1]
        # @info x_samp
        _, σ2 = predict_f(gp, x_samp)
        if σ2[1] > σ2_best
            σ2_best = σ2[1]
        end
    end
    return 0., 0., sqrt(σ2_best[1])
end

function create_x_matrix(xL, xU, N)
    dim = length(xL)
    xds = [(xU[i]-xL[i])/N for i=1:dim]
    # Assume 2D for now
    xd_prod = collect(Iterators.product(xL[1]:xds[1]:xU[1], xL[2]:xds[2]:xU[2]))
    return reshape(collect(Iterators.flatten(xd_prod)), dim, length(xd_prod))
end

function prepare_σ_gp(gp, x_L, x_U, ub; N=10, ll=-2.0)
    xn = create_x_matrix(x_L, x_U, N)
    yn = sqrt.(predict_f(gp, xn)[2])

    meanfcn = MeanConst(ub)
    kernel = SE(ll, 0.0)

    gp_σ = GP(xn, yn, meanfcn, kernel, 0.0)
    return gp_σ
end

function compute_σ_ub_bounds_from_gp(gp, x_L, x_U; ub=1.0)
    σ_gp = prepare_σ_gp(gp, x_L, x_U, ub)
    res_test = GPBounding.compute_μ_bounds_bnb(σ_gp, x_L, x_U, max_flag=true, max_iterations=4)
    return res_test[1], 0., res_test[3][1]+ub
end

end # module


text:
### A Pluto.jl notebook ###
# v0.12.11

using Markdown
using InteractiveUtils

# ╔═╡ 8dbcf370-1dec-11eb-0a2d-e967ec1f2fde
md" ### Nomenclatura"

# ╔═╡ c4a8db60-1dec-11eb-0c39-73a8aac79092
md"Seja $X$ uma matriz com dimensões $n x m$ com $n$ linhas e $m$ colunas que contém nossos exemplos. $x_i^{(j)}$ corresponde ao elemento de $X$ em que $i$ é uma linha da matriz (exemplo) e $j$ uma coluna da matriz (atributo), cujo acesso na linguagem representamos por x[i,j]. Para pegarmos um vetor linha desta matriz (exemplo) representamos formalmente por $X_i$ e na linguagem por x[i,:]."

# ╔═╡ cd532e00-1dec-11eb-06d9-9b7c76e9c1ba
md"Seja $\theta$ um vetor linha com $m$ elementos tal que $m$ corresponde ao número de atributos mais um. $\theta = [\theta^{(0)},\theta^{(1)}, ...,\theta^{(m)}]$ em que cada elemento corresponde ao j-ésimo parâmetro associado ao j-ésimo atributo. Na nossa linguagem acessamos um elemento por $\theta[i]$ e o vetor por $\theta$"

# ╔═╡ dc365e10-1dec-11eb-0e4f-95624ee468ea
md"Seja $y$ nosso vetor alvo em que representamos o i-ésimo elemento por $y_i$. Na nossa linguagem representamos por $y[i]$."

# ╔═╡ 571a5ea0-1dee-11eb-0eaf-bf746a6faa86
md"### Algumas libs importantes usadas neste experimento"

# ╔═╡ 317bbe60-1ded-11eb-2faa-c3e47c0df176
begin
   import Statistics: mean,std
   import BenchmarkTools: @benchmark, @btime
   import StaticArrays: @SMatrix
   import Plots: plot,plot!,scatter,scatter!
   import StatsBase: zscore
   μ = mean
end

# ╔═╡ f9af3a70-1dec-11eb-34f0-91b7fc591b87
md" ### Hipótese e modelo linear"

# ╔═╡ 78d2c0e0-1f98-11eb-06bc-b51c3b4b9f35
md"A seguir descrevemos a hipótese que consiste da aplicação de um modelo (conjunto de parâmetros) em um exemplo. Note que a hipótese sempre se refere à aplicação de um único exemplo e não ao conjunto de dados todo. "

# ╔═╡ 2e7709e0-1ded-11eb-2abd-d770bf1afe8a
md"$H_{\theta}(X_i) = \sum_{j}^{m}{x_{i}^{(j)} \theta^{(j)}} = X_i\theta^{T}$"

# ╔═╡ b05f9930-1f97-11eb-0083-07ebd21313e0
md"Exemplificando o produto de vetores, nossa hipótese faria o seguinte cálculo para um exemplo $X_i = \begin{bmatrix} 1 & 500 \end{bmatrix}$ e um vetor de parâmetros $\theta^T = \begin{bmatrix}
    0.5 \\
    1.0
\end{bmatrix}$:"

# ╔═╡ cbf998e0-1f96-11eb-0a8b-a3bf631bfe02
md"e.g. $H_{\theta}(X_i) = \begin{bmatrix}
    1 & 500
  \end{bmatrix}\begin{bmatrix}
    0.5 \\
    1.0
\end{bmatrix} = 500.5$"

# ╔═╡ 4cc76c80-1f98-11eb-0c70-29aec0b92ac4
md"Note que o resultado é um escalar, pois estamos estimando um valor contínuo (problema de regressão)."

# ╔═╡ c0b094a0-1f98-11eb-2efb-07f680447b12
md"### Implementação a hipótese"

# ╔═╡ be0deb50-1ded-11eb-06bb-33e94217e815
H(χ,θ) = χ*Θ'

# ╔═╡ 2f5a6500-1ded-11eb-11dd-bd88bf102b98
md" ### Função de custo"

# ╔═╡ 30972ac0-1ded-11eb-1007-1d3833dfc264

md"$J(\theta)= \frac{1}{2n}\sum_{i=1}^{n}{(H_{\theta}(X_i) - y_i)²}$"

# ╔═╡ bba6d590-1e1d-11eb-2df6-67bd6c5d4777
md" A seguir estou calculando a função de custo sem loop, apenas por matrizes.Desta maneira, geramos em uma multiplicação de matrizes, um vetor de predições com a mesma dimensão da saída. Após isto, subtraímos elemento a elemento. Após isto, elevamos cada elemento ao quadrado e depois tiramos uma média. Explicado na aula anterior."

# ╔═╡ 9837a3d0-1ded-11eb-3c49-17a45387e46e
J(Θ,χ,γ) = .5μ((χ*Θ' .- γ).^2)

# ╔═╡ 2ebaad22-1dee-11eb-2fb1-f51c49ec13fb
md"### Gradiente descendente"

# ╔═╡ d1cf31b0-1e17-11eb-22d1-916c9a331700
md"A seguir, mostramos o cálculo do gradiente. Como podemos notar, o expoente 2 desce e é eliminado pela constante 2 do denominador. Adicionalemente, vemos que é uma subtração do gradiente, pois estamos minimizando a função de custo. Por isto, estamos caminhando no sentido contrário ao do gradiente. É importante também ressaltar que o viés já está embutido na matriz $X$ na forma de uma coluna com valores 1. Desta forma, o cálculo do gradiente é o mesmo para o parâmetro do viés ($θ^{(0)}$) e os demais."

# ╔═╡ cbf654c0-1dfa-11eb-0994-85a67ded79b0
md"$θ^{(j)} = θ^{(j)} - α \frac{1}{n} \sum_{i=1}^{n}{(H_θ(X_i) - y_i)x_{i}^{(j)}}$"

# ╔═╡ 46484b80-1e1d-11eb-2bdc-e50b24f6b530
md"Entretanto, podemos condensar esta fórmula usando notação de vetores e matrizes. Lembrando que $\theta$ é um vetor agora e os cálculos são feitos simultaneamente para cada eleemnto deste. Podemos observar que o índice da coluna some, pois são processadas automaticamente. A seguir, veja a fórmula:"

# ╔═╡ f1089b10-1e1d-11eb-120e-cdbd993d72dc
md"$θ = θ - α \frac{1}{n} \sum_{i=1}^{n}{(H_θ(X_i) - y_i) \circ X_{i}}$"

# ╔═╡ 79bc38d0-1e1f-11eb-066c-2929340c479e
md"Cabe observar que $\circ$ é o produto de hadamard."

# ╔═╡ f5d33090-1e1f-11eb-3f25-d10f2616803d
md"### Algoritmo "

# ╔═╡ 7a22eab0-1def-11eb-2672-495fb9cf19f8
function train(θ,χ,γ;verbose=false,α = 1e-3,ϵ=1e-5,τ=1e3)

	e₂  = J(θ,χ,γ)
	e₁  = e₂ + .1
	#Δ₁,Δ₂ = e₁ - e₂ + 2ϵ,e₁ - e₂
	💻 = [[e₂ θ]]
	i   = 1
	while e₂ < e₁ && e₂ > ϵ && i < τ
	#while Δ₁ > Δ₂  && e₂ > ϵ && i < τ

		∇     = μ((χ*θ' .- γ).* χ,dims=1)
		θ    -= α .* ∇

		e₁,e₂ = e₂,J(θ,χ,γ)
		#Δ₁,Δ₂ = Δ₂,e₁ - e₂
		verbose && if (i%10==0) append!(💻,[[e₂ θ]]) end

		i += 1
	end
	θ,vcat(💻...)
end

# ╔═╡ 1f0e51f0-1df4-11eb-0feb-b53552fa63c9
md"### Testes de unitários"

# ╔═╡ 58dd7720-1e09-11eb-1904-af5f9991ece5
md" Sempre que for desenvolver suas funções, realize uma bateria de testes para cada uma antes mesmo de executar todo o seu programa. Isso vai lhe poupar muito tempo para uma futura depuração de código. Assim, para cada função, faça pelo menos um teste. Elabore casos bons de entrada para suas funções de forma que cubra bem o espaço de possibilidades. Seus casos devem ser bem simples para que num futuro, se der erro, você consiga rapidamente resolver."

# ╔═╡ 259ac040-1e39-11eb-3670-f5c5519a1978
md"### Padronizando os Dados "

# ╔═╡ 246cc792-1e39-11eb-3c20-f1623d8ea737
md"O objetivo de padronizar os dados é acelerar a convergência do método. Aqui utilizaremos a padronização zscore. Entretanto, para amostras pequenas, a padronização não funciona muito bem."

# ╔═╡ c467ca00-1e08-11eb-3d41-75a38532652b
md"### Verificando a curva de erro"

# ╔═╡ d1a88a00-1f8f-11eb-36b7-f39dabee0657
md"### Veificando o comportamento dos parâmetros"

# ╔═╡ b405694e-1f8f-11eb-0fad-5d0c3527d206
md"### Verificando o modelo"

# ╔═╡ d497051e-1f94-11eb-2ccb-a7d3e52570d1
md"Como podemos observar, o gradiente jogou o θ₁ para zero (viés, responsável pela translação) e ajustou o θ₂ (responsável pela rotação). Cabe ressaltar que estamos trabalhando com os dados padronizados. É possível, no entanto, despadronizar e obter predições no campo dos números originais."

# ╔═╡ 37c52e40-1e48-11eb-1e13-7ba1bebfabc4
md"### Depuração "

# ╔═╡ 44997b82-1e48-11eb-1821-f7f3aecd8093
md" Costumo colocar alguma variáveis soltas para que eu veja se está tudo correto. "

# ╔═╡ a419b220-1e4b-11eb-1299-295bf39df3f4
md"### Solução direta por álgebra linear"

# ╔═╡ a39cb9f0-1e4b-11eb-296e-bf19e5438416
train(X,Y) = inv(X'X)X'Y

# ╔═╡ 628624c0-1dfa-11eb-2afb-137bf77669e1
begin
	Θ1 = [0. 4.]
	χ =  [1 50.;1 60.;1 100.; 1 200.]
	γ =  [200  ;   240  ;   400; 800]
	@assert J(Θ1,χ,γ) == 0.0
	@assert train(Θ1,χ,γ)[1] == [0. 4.] # iniciando com a solução, o θ deve continuar igual
end

# ╔═╡ 4c632fee-1e39-11eb-00c6-c3f2a0774005
begin
	χₚ = hcat(ones((4,1)), zscore(χ[:,2:end])) # padronizamos só o que não é viés
	γₚ = zscore(γ)
end

# ╔═╡ d8ab5330-1e47-11eb-06b3-af3d2ee97c7b
χₚ

# ╔═╡ 86be2920-1e48-11eb-10f1-6740a2459951
γₚ

# ╔═╡ 74604cb2-1e05-11eb-19ef-09fd04dc52e3
(θᵏ,r) = train([4 5],χₚ,γₚ,verbose=true,α = 1e-3,ϵ=1e-5,τ=1e5);

# ╔═╡ b7913600-1e0c-11eb-362e-0b4b9ef9d1aa
plot(1:length(r[:,1]),r[:,1],xlabel="iterações",ylabel="J(θ)",title="Análise de erros",label="Erro")

# ╔═╡ fa0bd460-1e23-11eb-131f-f31bec373c07
begin
	p = plot(1:length(r[:,2]),r[:,2],xlabel="iterações",ylabel="θ₁,θ₂",title="Análise de parâmetros",label="θ₁")
	plot!(p,1:length(r[:,3]),r[:,3],xlabel="iterações",label="θ₂")
end

# ╔═╡ cc73d620-1f8f-11eb-357d-b19b50c19280
begin
	f(x) = θᵏ[1] + θᵏ[2]*x # o mesmo que H(X) lá do início
	plot(χₚ[:,2],f,xlabel="m² (padronizado)",ylabel="Preço (padronizado)",title="Análise do da predição do modelo",label="θ₁ + θ₂x")
	scatter!(χₚ[:,2],γₚ,label="Dados Coletados")
end

# ╔═╡ 41505b10-1e34-11eb-06db-750eb3ab11cf
r[end,:]

# ╔═╡ 29b04780-1e44-11eb-3bcf-8b5f121d157f
J(θᵏ,χₚ,γₚ)

# ╔═╡ d0de29e0-1e36-11eb-398c-5f0b55ea281e
θᵏ

# ╔═╡ 70775a70-1e4c-11eb-3cec-11fd54e9b637
train(χₚ,γₚ)

# ╔═╡ Cell order:
# ╟─8dbcf370-1dec-11eb-0a2d-e967ec1f2fde
# ╟─c4a8db60-1dec-11eb-0c39-73a8aac79092
# ╟─cd532e00-1dec-11eb-06d9-9b7c76e9c1ba
# ╟─dc365e10-1dec-11eb-0e4f-95624ee468ea
# ╟─571a5ea0-1dee-11eb-0eaf-bf746a6faa86
# ╠═317bbe60-1ded-11eb-2faa-c3e47c0df176
# ╟─f9af3a70-1dec-11eb-34f0-91b7fc591b87
# ╟─78d2c0e0-1f98-11eb-06bc-b51c3b4b9f35
# ╟─2e7709e0-1ded-11eb-2abd-d770bf1afe8a
# ╟─b05f9930-1f97-11eb-0083-07ebd21313e0
# ╟─cbf998e0-1f96-11eb-0a8b-a3bf631bfe02
# ╟─4cc76c80-1f98-11eb-0c70-29aec0b92ac4
# ╟─c0b094a0-1f98-11eb-2efb-07f680447b12
# ╠═be0deb50-1ded-11eb-06bb-33e94217e815
# ╟─2f5a6500-1ded-11eb-11dd-bd88bf102b98
# ╟─30972ac0-1ded-11eb-1007-1d3833dfc264
# ╟─bba6d590-1e1d-11eb-2df6-67bd6c5d4777
# ╠═9837a3d0-1ded-11eb-3c49-17a45387e46e
# ╟─2ebaad22-1dee-11eb-2fb1-f51c49ec13fb
# ╟─d1cf31b0-1e17-11eb-22d1-916c9a331700
# ╟─cbf654c0-1dfa-11eb-0994-85a67ded79b0
# ╟─46484b80-1e1d-11eb-2bdc-e50b24f6b530
# ╟─f1089b10-1e1d-11eb-120e-cdbd993d72dc
# ╟─79bc38d0-1e1f-11eb-066c-2929340c479e
# ╟─f5d33090-1e1f-11eb-3f25-d10f2616803d
# ╠═7a22eab0-1def-11eb-2672-495fb9cf19f8
# ╟─1f0e51f0-1df4-11eb-0feb-b53552fa63c9
# ╟─58dd7720-1e09-11eb-1904-af5f9991ece5
# ╠═628624c0-1dfa-11eb-2afb-137bf77669e1
# ╟─259ac040-1e39-11eb-3670-f5c5519a1978
# ╟─246cc792-1e39-11eb-3c20-f1623d8ea737
# ╠═4c632fee-1e39-11eb-00c6-c3f2a0774005
# ╟─c467ca00-1e08-11eb-3d41-75a38532652b
# ╠═74604cb2-1e05-11eb-19ef-09fd04dc52e3
# ╟─b7913600-1e0c-11eb-362e-0b4b9ef9d1aa
# ╟─d1a88a00-1f8f-11eb-36b7-f39dabee0657
# ╠═fa0bd460-1e23-11eb-131f-f31bec373c07
# ╟─b405694e-1f8f-11eb-0fad-5d0c3527d206
# ╟─d497051e-1f94-11eb-2ccb-a7d3e52570d1
# ╠═cc73d620-1f8f-11eb-357d-b19b50c19280
# ╟─37c52e40-1e48-11eb-1e13-7ba1bebfabc4
# ╟─44997b82-1e48-11eb-1821-f7f3aecd8093
# ╠═41505b10-1e34-11eb-06db-750eb3ab11cf
# ╠═29b04780-1e44-11eb-3bcf-8b5f121d157f
# ╠═d0de29e0-1e36-11eb-398c-5f0b55ea281e
# ╠═d8ab5330-1e47-11eb-06b3-af3d2ee97c7b
# ╠═86be2920-1e48-11eb-10f1-6740a2459951
# ╟─a419b220-1e4b-11eb-1299-295bf39df3f4
# ╠═a39cb9f0-1e4b-11eb-296e-bf19e5438416
# ╠═70775a70-1e4c-11eb-3cec-11fd54e9b637

max_stars_repo_name: lalvim/machinelearning

max_issues_repo_name: lalvim/machinelearning

max_forks_repo_name: lalvim/machinelearning


text:
module Yelp

using Nettle, Requests

type OAuth
    consumer_secret::String
    oauth_secret::String
    consumer_key::String
    oauth_token::String
end

type Query
    cred::OAuth
    baseurl::String
    params::Dict
end

encode(s::Dict)   = tostring(s,x->Requests.encodeURI(string(x)))
encode(s::String) = Requests.encodeURI(s)

function tostring(params::Dict, func=identity)
    pkeys = collect(keys(params))
    sort!(pkeys)
    q = ""
    for k in pkeys
        k == "oauth_signature" && continue
        v = func(params[k])
        q *= "$k=$v&"
    end
    return chop(q)
end

function oauth_sign!(q::Query)
    signing_key = "$(encode(q.cred.consumer_secret))&$(encode(q.cred.oauth_secret))"
    h = Nettle.HMACState("SHA1", signing_key)
    sig = "GET&$(encode(q.baseurl))&$(encode(tostring(q.params)))"
    # for whatever reason, Yelp doesn't like a few values
    # that Requests.encodeURI puts in
    sig = replace(sig,"%2C","%252C")
    sig = replace(sig,"%20","%2520")
    sig = replace(sig,"%7","%257")
    println(sig)
    Nettle.update!(h, sig)
    q.params["oauth_signature"] = encode(base64(Nettle.digest!(h)))
    return "OAuth realm= oauth_consumer_key=\"$(q.params["oauth_consumer_key"])\", oauth_nonce=\"$(q.params["oauth_nonce"])\", oauth_signature=\"$(q.params["oauth_signature"])\", oauth_signature_method=\"$(q.params["oauth_signature_method"])\", oauth_timestamp=\"$(q.params["oauth_timestamp"])\", oauth_token=\"$(q.params["oauth_token"])\""
end

SEARCH_API = "http://api.yelp.com/v2/search"
BUSINESS_API = "http://api.yelp.com/v2/business"

BEST_MATCHED = 0
DISTANCE = 1
HIGHEST_RATED = 2
SORT_OPTIONS = [BEST_MATCHED,DISTANCE,HIGHEST_RATED]

MAX_RADIUS_FILTER = 40000

typealias LatLon (Float64,Float64)

function Base.search(cred::OAuth;
                     term::String = "",
                     location::String = "San Francisco, CA",
                     cll::LatLon = (0.0,0.0),
                     ll::LatLon = (0.0,0.0),
                     sw_bounds::LatLon = (0.0,0.0),
                     ne_bounds::LatLon = (0.0,0.0),
                     limit::Int = 5,
                     offset::Int = 0,
                     sort::Int = BEST_MATCHED,
                     category_filter::String = "",
                     radius_filter::Int = 40000,
                     deals_filter::Bool = false,
                     cc::String = "US",
                     lang::String = "en")
    # input validation
    radius_filter <= MAX_RADIUS_FILTER || throw(ArgumentError("radius_filter = $radius_filter; must be <= $MAX_RADIUS_FILTER"))
    sort in SORT_OPTIONS || throw(ArgumentError("sort = $sort; must be in $SORT_OPTIONS"))

    # build options dict
    params = ["limit"         => limit,
              "offset"        => offset,
              "sort"          => sort,
              "radius_filter" => radius_filter,
              "deals_filter"  => deals_filter,
              "cc"            => cc,
              "lang"          => lang]

    # only include these parameters if explicitly set
    term            == "" || (params["term"] = term)
    category_filter == "" || (params["category_filter"] = category_filter)

    # add a few OAuth parameters
    params["oauth_consumer_key"] = cred.consumer_key
    params["oauth_token"] = cred.oauth_token
    params["oauth_signature_method"] = "HMAC-SHA1"

    # figure out location parameter
    if ll != (0.0,0.0)
        params["ll"] = string(ll)[2:end-1]
    elseif sw_bounds != (0.0,0.0) && ne_bounds != (0.0,0.0)
        s = string(string(sw_bounds)[2:end-1],"|",string(ne_bounds)[2:end-1])
        params["bounds"] = s
    else
        params["location"] = location
        if cll != (0.0,0.0)
            params["cll"] = string(cll)[2:end-1]
        end
    end

    # construct and return a Query instance
    return Query(cred,SEARCH_API,params)
end

function business(cred, business_id; cc="US",lang="en",lang_filter::Bool=false)
    params = ["cc" => cc, "lang" => lang, "lang_filter" => lang_filter]
    return Query(cred,BUSINESS_API*business_id,params)
end

function Base.get(q::Query)
    q.params["oauth_nonce"] = randstring(32)
    q.params["oauth_timestamp"] = int(time())
    oauth_header = oauth_sign!(q)
    url = URI("$(q.baseurl)?$(encode(q.params))")
    return Requests.get(url;
                headers = {"Content-Type" => "application/x-www-form-urlencoded",
                "Authorization" => oauth_header,
                "Accept" => "*/*"})
end

end # module


text:
@userplot Plotoptimmodel

@recipe function f(p::Plotoptimmodel; rvec=nothing, δt=1.0, attrib=:obs)
	pa=p.args[1]
	(rvec===nothing) && (rvec=1:pa.nr)

	#=

	gxobs=zeros(2*size(pa.obs.g,1)-1, size(pa.obs.g,2))
	gxcal=similar(gxobs)
	scobs=vecnorm(pa.obs.g[:,1])^2
	sccal=vecnorm(pa.cal.g[:,1])^2
	for ir in 1:pa.nr
		gxobs[:,ir] = xcorr(pa.obs.g[:,1], pa.obs.g[:, ir])/scobs
		gxcal[:,ir] = xcorr(pa.cal.g[:,1], pa.cal.g[:, ir])/sccal
	end
	#
	g=collect(1:pa.ntg)*δt
	gx=collect(-pa.ntg+1:1:pa.ntg-1)*δt

	# time vectors
	# autocorr s
	asobs=autocor(pa.obs.s[:,1], 1:pa.nt-1, demean=true)
	as=autocor(pa.cal.s[:,1], 1:pa.nt-1, demean=true)
	sli=max(maximum(abs,asobs), maximum(abs,as))
	# autocorr g
	agobs=autocor(pa.obs.g,1:pa.ntg-1, demean=true)
	ag=autocor(pa.cal.g,1:pa.ntg-1, demean=true)
	gli=max(maximum(abs,agobs), maximum(abs,ag))

#	asobs=asobs[1:fact:ns] # resample
#	as=as[1:fact:ns] # resample

	#fact=(pa.nt*pa.nr>1000) ? round(Int,pa.nt*pa.nr/1000) : 1
	# cut receivers
#	dcal=pa.cal.d[1:fact:pa.nt*pa.nr]
#	dobs=pa.obs.d[1:fact:pa.nt*pa.nr]

=#

	if((attrib==:obs) || (attrib==:cal))
		layout := (3,1)
		g=getfield(pa,attrib).g
		s=getfield(pa,attrib).s
		temp=getfield(pa,attrib)
		dd=getfield(temp, fieldnames(temp)[5])

		ns=length(s)
		fact=(ns>1000) ? round(Int,ns/1000) : 1
		fnrp=(pa.nr>10) ? round(Int,pa.nr/10) : 1

		@series begin
			subplot := 1
	#		aspect_ratio := :auto
			legend := false
	#		l := :plot
			title := string("\$g_i\$ ", attrib)
			w := 1
			g[:,1:fnrp:end]
		end
		@series begin
			subplot := 2
	#		aspect_ratio := :auto
			legend := false
			title := string("\$s\$ ", attrib)
			w := 1
			s[1:fact:end,1]
		end
		@series begin
			subplot := 3
	#		aspect_ratio := :auto
			legend := false
			title := string("\$d_i\$ ", attrib)
			w := 1
			dd[1:fact:end,1:fnrp:end]
		end
	end

	if(attrib==:x)
		ns=length(pa.nt)
		fact=(ns>1000) ? round(Int,ns/1000) : 1

		fnrp=(pa.nr>10) ? round(Int,pa.nr/10) : 1

		xsobs=pa.obs.s[1:fact:end,1]
		xscal=pa.cal.s[1:fact:end,1]
		xdcal=getfield(pa.cal,fieldnames(pa.cal)[5])[1:fact*pa.nr:end]
	        xdobs=getfield(pa.obs,fieldnames(pa.obs)[5])[1:fact*pa.nr:end]
		xgcal=pa.cal.g[:,1:fnrp:end][:]
		xgobs=pa.obs.g[:,1:fnrp:end][:]
		layout := (1,3)
		@series begin
			subplot := 1
			aspect_ratio := :equal
			seriestype := :scatter
			title := "scatter s"
			legend := false
			normalize(xsobs), normalize(xscal)
		end
		@series begin
			subplot := 2
			aspect_ratio := :equal
			seriestype := :scatter
			title := "scatter g"
			legend := false
			normalize(xgobs), normalize(xgcal)
		end

		@series begin
			subplot := 3
			aspect_ratio := :equal
			seriestype := :scatter
			title := "scatter d"
			legend := false
			normalize(xdobs), normalize(xdcal)
		end
	end

	if(attrib==:precon)
		layout := (3,1)
		@series begin
			subplot := 1
	#		aspect_ratio := :auto
			legend := false
			title := "gprecon"
			w := 1
			pa.gprecon
		end
		@series begin
			subplot := 2
	#		aspect_ratio := :auto
			legend := false
			title := "sprecon"
			w := 1
			pa.sprecon
		end
		@series begin
			subplot := 3
	#		aspect_ratio := :auto
			legend := false
			title := "gweights"
			w := 1
			pa.gweights
		end

	end
end

@userplot Plotobsmodel

@recipe function f(p::Plotobsmodel; rvec=nothing, δt=1.0)
	pa=p.args[1]
	layout := (3,1)

	tg=collect(1:pa.ntg)*δt
	td=collect(1:pa.nt)*δt

	@series begin
		subplot := 1
#		aspect_ratio := :auto
		legend := false
#		l := :plot
		title := "\$g_i\$"
		w := 1
		tg, pa.g
	end
	@series begin
		subplot := 2
#		aspect_ratio := :auto
		legend := false
		title := "\$s\$"
		w := 1
		td, pa.s
	end
	@series begin
		subplot := 3
#		aspect_ratio := :auto
		legend := false
		title := "\$d_i\$"
		w := 1
		td, getfield(pa,fieldnames(pa)[6])
	end
end

@userplot Plotprmodel

@recipe function f(p::Plotprmodel; gobs=nothing, rvec=nothing, δt=1.0, attrib=:obs)
	pa=p.args[1]
	nr=size(pa.g,2)
	ntg=size(pa.g,1)
	(rvec===nothing) && (rvec=1:nr)

	layout := (3,2)
	g1=pa.g
	if(gobs===nothing)
		gobs=zeros(g1)
	end
	dd0=hcat(pa.p_misfit_xcorr.cy...)
	dd=hcat(pa.p_misfit_xcorr.pxcorr.cg...)

	nr1=size(dd,2)
	ns=size(dd,1)
	fact=(ns>1000) ? round(Int,ns/1000) : 1
	fnrp=(nr1>10) ? round(Int,nr1/10) : 1

	@series begin
		subplot := 1
#		aspect_ratio := :auto
		legend := false
#		l := :plot
		title := string("\$g_i\$ ", "obs")
		w := 1
		gobs
	end
	@series begin
		subplot := 2
#		aspect_ratio := :auto
		legend := false
#		l := :plot
		title := string("\$g_i\$ ", "cal")
		w := 1
		g1
	end
	@series begin
		subplot := 3
#		aspect_ratio := :auto
		legend := false
		title := string("\$d_i\$ ", "obs")
		w := 1
		dd0[1:fact:end,1:fnrp:end]
	end
	@series begin
		subplot := 4
#		aspect_ratio := :auto
		legend := false
		title := string("\$d_i\$ ", "cal")
		w := 1
		dd[1:fact:end,1:fnrp:end]
	end

	xdcal=dd[1:fact*nr1:end]
	xdobs=dd0[1:fact*nr1:end]
	@series begin
		subplot := 5
		aspect_ratio := :equal
		seriestype := :scatter
		title := "scatter d"
		legend := false
		normalize(xdobs), normalize(xdcal)
	end
end


text:
using Gnuplot
x = LinRange(-10,10,200)
@gp x sin.(x) sin.(x) .+ 1 "with filledcu lc '#56B4E9' fs transparent solid 0.3"
@gp :- x cos.(x) 1 .+ cos.(x) "with filledcu lc 'red' fs transparent solid 0.5"
save(term="pngcairo size 600,400", output="./code/plt_ex15.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples