text:
export iauApio13
"""
For a terrestrial observer, prepare star-independent astrometry
parameters for transformations between CIRS and observed
coordinates.  The caller supplies UTC, site coordinates, ambient air
conditions and observing wavelength.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
utc1   double      UTC as a 2-part...
utc2   double      ...quasi Julian Date (Notes 1,2)
dut1   double      UT1-UTC (seconds)
elong  double      longitude (radians, east +ve, Note 3)
phi    double      geodetic latitude (radians, Note 3)
hm     double      height above ellipsoid (m, geodetic Notes 4,6)
xp,yp  double      polar motion coordinates (radians, Note 5)
phpa   double      pressure at the observer (hPa = mB, Note 6)
tc     double      ambient temperature at the observer (deg C)
rh     double      relative humidity at the observer (range 0-1)
wl     double      wavelength (micrometers, Note 7)

Returned:
astrom iauASTROM*  star-independent astrometry parameters:
pmt    double       unchanged
eb     double[3]    unchanged
eh     double[3]    unchanged
em     double       unchanged
v      double[3]    unchanged
bm1    double       unchanged
bpn    double[3][3] unchanged
along  double       longitude + s' (radians)
xpl    double       polar motion xp wrt local meridian (radians)
ypl    double       polar motion yp wrt local meridian (radians)
sphi   double       sine of geodetic latitude
cphi   double       cosine of geodetic latitude
diurab double       magnitude of diurnal aberration vector
eral   double       "local" Earth rotation angle (radians)
refa   double       refraction constant A (radians)
refb   double       refraction constant B (radians)

Returned (function value):
        int         status: +1 = dubious year (Note 2)
                                0 = OK
                        -1 = unacceptable date

Notes:

1.  utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
convenient way between the two arguments, for example where utc1
is the Julian Day Number and utc2 is the fraction of a day.

However, JD cannot unambiguously represent UTC during a leap
second unless special measures are taken.  The convention in the
present function is that the JD day represents UTC days whether
the length is 86399, 86400 or 86401 SI seconds.

Applications should use the function iauDtf2d to convert from
calendar date and time of day into 2-part quasi Julian Date, as
it implements the leap-second-ambiguity convention just
described.

2.  The warning status "dubious year" flags UTCs that predate the
introduction of the time scale or that are too far in the future
to be trusted.  See iauDat for further details.

3)  UT1-UTC is tabulated in IERS bulletins.  It increases by exactly
one second at the end of each positive UTC leap second,
introduced in order to keep UT1-UTC within +/- 0.9s.  n.b. This
practice is under review, and in the future UT1-UTC may grow
essentially without limit.

4.  The geographical coordinates are with respect to the WGS84
reference ellipsoid.  TAKE CARE WITH THE LONGITUDE SIGN:  the
longitude required by the present function is east-positive
(i.e. right-handed), in accordance with geographical convention.

5.  The polar motion xp,yp can be obtained from IERS bulletins.  The
values are the coordinates (in radians) of the Celestial
Intermediate Pole with respect to the International Terrestrial
Reference System (see IERS Conventions 2003), measured along the
meridians 0 and 90 deg west respectively.  For many applications,
xp and yp can be set to zero.

Internally, the polar motion is stored in a form rotated onto
the local meridian.

6.  If hm, the height above the ellipsoid of the observing station
in meters, is not known but phpa, the pressure in hPa (=mB), is
available, an adequate estimate of hm can be obtained from the
expression

        hm = -29.3 * tsl * log ( phpa / 1013.25 );

where tsl is the approximate sea-level air temperature in K
(See Astrophysical Quantities, C.W.Allen, 3rd edition, section
52).  Similarly, if the pressure phpa is not known, it can be
estimated from the height of the observing station, hm, as
follows:

        phpa = 1013.25 * exp ( -hm / ( 29.3 * tsl ) );

Note, however, that the refraction is nearly proportional to the
pressure and that an accurate phpa value is important for
precise work.

7.  The argument wl specifies the observing wavelength in
micrometers.  The transition from optical to radio is assumed to
occur at 100 micrometers (about 3000 GHz).

8.  It is advisable to take great care with units, as even unlikely
values of the input parameters are accepted and processed in
accordance with the models used.

9.  In cases where the caller wishes to supply his own Earth
rotation information and refraction constants, the function
iauApc can be used instead of the present function.

10. This is one of several functions that inserts into the astrom
structure star-independent parameters needed for the chain of
astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

The various functions support different classes of observer and
portions of the transformation chain:

        functions         observer        transformation

iauApcg iauApcg13    geocentric      ICRS <-> GCRS
iauApci iauApci13    terrestrial     ICRS <-> CIRS
iauApco iauApco13    terrestrial     ICRS <-> observed
iauApcs iauApcs13    space           ICRS <-> GCRS
iauAper iauAper13    terrestrial     update Earth rotation
iauApio iauApio13    terrestrial     CIRS <-> observed

Those with names ending in "13" use contemporary SOFA models to
compute the various ephemerides.  The others accept ephemerides
supplied by the caller.

The transformation from ICRS to GCRS covers space motion,
parallax, light deflection, and aberration.  From GCRS to CIRS
comprises frame bias and precession-nutation.  From CIRS to
observed takes account of Earth rotation, polar motion, diurnal
aberration and parallax (unless subsumed into the ICRS <-> GCRS
transformation), and atmospheric refraction.

11. The context structure astrom produced by this function is used
by iauAtioq and iauAtoiq.

Called:
iauUtctai    UTC to TAI
iauTaitt     TAI to TT
iauUtcut1    UTC to UT1
iauSp00      the TIO locator s', IERS 2000
iauEra00     Earth rotation angle, IAU 2000
iauRefco     refraction constants for given ambient conditions
iauApio      astrometry parameters, CIRS-observed

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauApio13(utc1::Real, utc2::Real, dut1::Real,
               elong::Real, phi::Real, hm::Real,
               xp::Real, yp::Real,
               phpa::Real, tc::Real,
               rh::Real, wl::Real)
    # Allocate return value
    ref_astrom = Ref{iauASTROM}(iauASTROM())

    status = ccall((:iauApio13, libsofa_c), Cint,
            (Cdouble, Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
            Ref{iauASTROM}),
            convert(Float64, utc1),
            convert(Float64, utc2),
            convert(Float64, dut1),
            convert(Float64, elong),
            convert(Float64, phi),
            convert(Float64, hm),
            convert(Float64, xp),
            convert(Float64, yp),
            convert(Float64, phpa),
            convert(Float64, tc),
            convert(Float64, rh),
            convert(Float64, wl),
            ref_astrom)

    return status, ref_astrom[]
end

ext: jl


text:
## random number generation

using Random
import RandomNumbers

# global state

# shared memory with the actual seed, per warp, loaded lazily or overridden by calling `seed!`
@eval @inline function global_random_keys()
    ptr = Base.llvmcall(
        $("""@global_random_keys = weak addrspace($(AS.Shared)) global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Shared))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Shared))* @global_random_keys, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Shared))* %ptr to i8 addrspace($(AS.Shared))*
                 ret i8 addrspace($(AS.Shared))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Shared}, Tuple{})
    CuDeviceArray((32,), ptr)
end

# shared memory with per-warp counters, incremented when generating numbers
@eval @inline function global_random_counters()
    ptr = Base.llvmcall(
        $("""@global_random_counters = weak addrspace($(AS.Shared)) global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Shared))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Shared))* @global_random_counters, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Shared))* %ptr to i8 addrspace($(AS.Shared))*
                 ret i8 addrspace($(AS.Shared))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Shared}, Tuple{})
    CuDeviceArray((32,), ptr)
end

@device_override Random.make_seed() = clock(UInt32)

# generators

using Random123: philox2x_round, philox2x_bumpkey

# GPU-compatible/optimized version of the generator from Random123.jl
struct Philox2x32{R} <: RandomNumbers.AbstractRNG{UInt64}
    @inline function Philox2x32{R}() where R
        rng = new{R}()
        if rng.key == 0
            # initialize the key. this happens when first accessing the (0-initialized)
            # shared memory key from each block. if we ever want to make the device seed
            # controlable from the host, this would be the place to read a global seed.
            #
            # note however that it is undefined how shared memory persists across e.g.
            # launches, so we may not be able to rely on the zero initalization then.
            rng.key = Random.make_seed()
        end
        return rng
    end
end

# default to 7 rounds; enough to pass SmallCrush
@inline Philox2x32() = Philox2x32{7}()

@inline function Base.getproperty(rng::Philox2x32, field::Symbol)
    threadId = threadIdx().x + (threadIdx().y - 0x1) * blockDim().x +
                               (threadIdx().z - 0x1) * blockDim().x * blockDim().y
    warpId = (threadId - 0x1) >> 0x5 + 0x1  # fld1

    if field === :seed
        @inbounds global_random_seed()[1]
    elseif field === :key
        @inbounds global_random_keys()[warpId]
    elseif field === :ctr1
        @inbounds global_random_counters()[warpId]
    elseif field === :ctr2
        blockId = blockIdx().x + (blockIdx().y - 0x1) * gridDim().x +
                                 (blockIdx().z - 0x1) * gridDim().x * gridDim().y
        globalId = threadId + (blockId - 0x1) * (blockDim().x * blockDim().y * blockDim().z)
        globalId%UInt32
    end::UInt32
end

@inline function Base.setproperty!(rng::Philox2x32, field::Symbol, x)
    threadId = threadIdx().x + (threadIdx().y - 0x1) * blockDim().x +
                               (threadIdx().z - 0x1) * blockDim().x * blockDim().y
    warpId = (threadId - 0x1) >> 0x5 + 0x1  # fld1

    if field === :key
        @inbounds global_random_keys()[warpId] = x
    elseif field === :ctr1
        @inbounds global_random_counters()[warpId] = x
    end
end

@device_override @inline Random.default_rng() = Philox2x32()

"""
    Random.seed!(rng::Philox2x32, seed::Integer, [counter::Integer=0])

Seed the on-device Philox2x32 generator with an UInt32 number.
Should be called by at least one thread per warp.
"""
function Random.seed!(rng::Philox2x32, seed::Integer, counter::Integer=0)
    rng.key = seed % UInt32
    rng.ctr1 = counter
    return
end

if VERSION >= v"1.7-"
@device_override Random.seed!(::Random._GLOBAL_RNG, seed) =
    Random.seed!(Random.default_rng(), seed)
end

"""
    Random.rand(rng::Philox2x32, UInt32)

Generate a byte of random data using the on-device Tausworthe generator.
"""
function Random.rand(rng::Philox2x32{R},::Type{UInt64}) where {R}
    ctr1, ctr2, key = rng.ctr1, rng.ctr2, rng.key

    if R > 0                               ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 1  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 2  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 3  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 4  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 5  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 6  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 7  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 8  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 9  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 10 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 11 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 12 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 13 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 14 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 15 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end

    # update the warp counter
    # NOTE: this performs the same update on every thread in the warp, but each warp writes
    #       to a unique location so the duplicate writes are innocuous
    # XXX: what if this overflows? we can't increment ctr2. bump the key?
    rng.ctr1 += 0x1

    # NOTE: it's too expensive to keep both numbers around in case the user only wanted one,
    #       so just make our 2x32 generator return 64-bit numbers by default.
    return (ctr1 % UInt64) << 32 | (ctr2 % UInt64)
end


text:
alphas = [4.0] .* pi/180
ploti = 1:length(alphas)
wing_b = 640e-3 * 2
wing_TR = 1.0 # hershey bar
wing_c = 240e-3
wing_AR = wing_b / wing_c
wing_θroot = 0.0
wing_θtip = 0.0
wing_ϕ = 0.0
simulationdata = AS.lift_distribution_template(ploti, alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_ϕ; surfacenames = ["PROWIM, w/o props"])
objective = AS.runsimulation!(simulationdata...)

y2b_cn_0deg_noprops = [
    0.09694654312188172 0.00459016393442635;
    0.17520171062009982 0.00721311475409836;
    0.2525958660014256 0.005901639344262133;
    0.31520456165359956 0.005901639344262355;
    0.37782466143977206 0.0006557377049181135;
    0.4091290092658591 0.0006557377049181135;
    0.5308624376336422 0.003278688524590123;
    0.5621667854597292 0.003278688524590123;
    0.6247811831789025 0.0006557377049181135;
    0.6865174625801855 0.0019672131147541183;
    0.8108652886671419 0.0019672131147541183;
    0.8717405559515325 -0.0006557377049178914;
    0.9186970776906629 -0.0006557377049181135
]

y2b_cn_4deg_noprops = [
    0.09710049893086248 0.3337704918032788;
    0.17448610121168928 0.3363934426229509;
    0.25274982181040634 0.3350819672131151;
    0.3153670705630795 0.33114754098360666;
    0.37714041339985754 0.31540983606557405;
    0.40933143264433364 0.30754098360655724;
    0.5311190306486102 0.28524590163934427;
    0.5615424091233072 0.2904918032786885;
    0.6250292230933714 0.2865573770491804;
    0.6867883107626516 0.277377049180328;
    0.7485759087669281 0.25508196721311494;
    0.8112330719885961 0.23278688524590174;
    0.8730434782608698 0.20000000000000018;
    0.9192045616535995 0.16590163934426239;
]

y2b_cn_10deg_noprops = [
    0.09708624376336422 0.740327868852459;
    0.1744803991446902 0.739016393442623;
    0.2518859586600142 0.7324590163934428;
    0.314520313613685 0.7206557377049181;
    0.3771746258018532 0.6996721311475411;
    0.4085217391304348 0.6800000000000002;
    0.5303292943692088 0.6485245901639345;
    0.5616136849607984 0.657704918032787;
    0.62424233784747 0.6485245901639345;
    0.686041339985745 0.6209836065573772;
    0.8114240912330721 0.5449180327868852;
    0.8724533143264435 0.47147540983606573;
    0.919552387740556 0.4059016393442625;
]

prowim_liftdistribution = [
    y2b_cn_0deg_noprops,
    y2b_cn_4deg_noprops,
    y2b_cn_10deg_noprops
]

fig = plt.figure("liftdistribution")
axs = fig.get_axes()
alpha_labels = [0,4,10]
for (idata, data) in enumerate(prowim_liftdistribution)
    cratio = idata / length(alpha_labels)
    axs[3].scatter(data[:,1], data[:,2], marker = "x", color = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label=LS.L"PROWIM, \alpha = " * "$(alpha_labels[idata])" * LS.L"^\circ")
end
axs[3].legend(loc="upper left", bbox_to_anchor=(1.01,1))
# @test objective == 0

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
nums = [parse(Int64, x) for x in readlines("inputs/01.txt")]
increases = 0
for i in 2:length(nums)
  if nums[i] > nums[i-1]
    increases += 1
  end
end
print(increases)


text:
using Revise
using ZigZagBoomerang
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using Makie, AbstractPlotting
cd(@__DIR__)
function gridlaplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1.
                    S[linear[i2, j2], linear[i, j]] -= 1.

                    S[linear[i, j], linear[i, j]] += 1.
                    S[linear[i2, j2], linear[i2, j2]] += 1.
                end
            end
        end
    end
    S
end

# Define precision operator of a Gaussian random field (sparse matrix operating on `vec`s of `n*n` matrices)
#n = 100
n = 60
const σ2 = 0.5
Γ0 = 2gridlaplacian(Float64, n, n)
Γ = 0.1I + Γ0
mat(x) = reshape(x, (n, n)) # vector to matrix
function mat0(y)
 #   mat(y + 0.1*sign.(y) .- 0.0)
    mat(y  .- 0.1)
end
# Γ is very sparse
@show nnz(Γ)/length(Γ) # 0.000496

# Corresponding Gaussian potential
# ϕ(x', Γ) = 0.5*x'*Γ*x  # not needed

# Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i]
#∇ϕ(x, i, Γ, y) = ZigZagBoomerang.idot(Γ, i, x) - ZigZagBoomerang.idot(Γ, i, y)    # more efficient that dot(Γ[:, i], x)
∇ϕ(x, i, Γ, y) = ZigZagBoomerang.idot(Γ, i, x)  + (x[i]-y[i])/σ2 # more efficient that dot(Γ[:, i], x)

# Random initial values
t0 = 0.0
h(x, y) = x^2+(5y/4-sqrt(abs(x)))^2
heart = [ max.(1 - h(x, y), 0).^0.3 for x in range(-1.5,1.5, length=n), y   in range(-1.5,1.5, length=n)]

μ0 = 7.0*vec(heart)
y = μ = μ0 + randn(n*n)
μpost = yhat = (I + Γ)\y
Γpost = (Γ + I)/σ2

#x0 = 0.2μ + 0*randn(n*n)
θ0 = rand([-1.0,1.0], n*n)

# Rejection bounds
c = [norm(Γpost[:, i], 2) for i in 1:n*n]

# Define ZigZag
Z = ZigZag(Γpost, μpost)
# or try the FactBoomerang
#Z = FactBoomerang(Γ, x0*0, 0.1)

κ = 0.01
# Run sparse ZigZag for T time units and collect trajectory
T = 200.0
@time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ, μ; adapt = false)

@time traj0 = collect(discretize(trace, 0.2))

su = false
adapt = false
trace, (t, x, θ), (acc, num), c = @time sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ, Γ, μ;
                                                strong_upperbounds = su ,
                                                adapt = adapt)

@time traj = collect(discretize(trace, 0.2))

# Prepare surface plot
M = Node(mat0(traj[end].second))
scene = Scene()
surface!(scene, M; shading=false, show_axis=false, colormap = :oleron, colorrange = (-3,3))
scale!(scene, 1.0, 1.0, 2.0)
zlims!(scene, -5.0, 5.0)
#display(scene)

scene1, _ = heatmap(abs.([mat(mean(last.(traj0[end÷2:end]))) mat(mean(last.(traj[end÷2:end]))) ]))
scene2, _ = heatmap(abs.([mat(μ0 - mean(last.(traj0[end÷2:end]))) mat(μ0 - mean(last.(traj[end÷2:end]))) ]))

error("stop")
# Movie frames
mkpath(joinpath(@__DIR__, "output"))
for i in eachindex(traj)
    M[] = mat0(traj[i].second)
    FileIO.save(joinpath(@__DIR__, "output", "surfs$i.png"), scene)
end

# Make video
dir = joinpath(@__DIR__, "output")
run(`ffmpeg -y -r 40 -f image2 -i $dir/surfs%d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p $(typeof(Z).name)field.mp4`)


text:
using Logging
logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

using DataStructures

#>> Sort by char <<#

d = SortedDict{Char,Int}(Base.Order.Forward) # Base.Order.Reverse

d['c'] = 3
d['d'] = 4
d['a'] = 1

col = collect(d)
@info typeof(col)
@info col

for p in col
    println("$(p[1]) $(p[2])")
end

#>> Sort by DateTime <<#

using Dates
d2 = SortedDict{DateTime, Int}(Base.Order.Forward)
d2[now()] = 13
sleep(0.1)
d2[now()] = 17
sleep(0.2)
d2[now()] = 19

col2 = collect(d2)
@info typeof(col2)
@info col2

#>> Sort by Date <<#

# SortedDict of Dict
d3 = SortedDict{Date, Dict}(Base.Order.Forward)

d3[Date(2020, 12, 24)] = Dict("A"=>1)
d3[Date(2020, 12, 23)] = Dict("B"=>1)
d3[Date(2020, 12, 22)] = Dict("C"=>3)
@info d3

# Add item existing Dict
# Indeed we are modifying a ref not a copy
d3[Date(2020, 12, 23)]["D"] = 4
@info d3

@info "Values are in order: $(values(d3))"

@info in(Date(2020, 12, 23), keys(d3))
@info in(Date(2019, 1, 1), keys(d3))
@info haskey(d3, Date(2020, 12, 23))
@info haskey(d3, Date(2019, 1, 1))

nothing # suppress printing the last line


text:
# example taken from Aragón, J. L., R. A. Barrio, T. E. Woolley, R. E. Baker, and P. K. Maini. “Nonlinear Effects on Turing Patterns: Time Oscillations and Chaos.” Physical Review E 86, no. 2 (August 8, 2012): 026201. https://doi.org/10.1103/PhysRevE.86.026201.
using Revise
	using DiffEqOperators, ForwardDiff, DifferentialEquations, Sundials
	using PseudoArcLengthContinuation, LinearAlgebra, Plots, SparseArrays, Parameters, Setfield
	const PALC = PseudoArcLengthContinuation

norminf = x -> norm(x, Inf)
f(u, v, p) = p.η * (      u + p.a * v - p.C * u * v - u * v^2)
g(u, v, p) = p.η * (p.H * u + p.b * v + p.C * u * v + u * v^2)

function Laplacian(N, lx, bc = :Dirichlet)
	hx = 2lx/N
	D2x = CenteredDifference(2, 2, hx, N)
	if bc == :Neumann
		Qx = Neumann0BC(hx)
	elseif bc == :Dirichlet
		Qx = Dirichlet0BC(typeof(hx))
	end
	D2xsp = sparse(D2x * Qx)[1]
end

function NL!(dest, u, p, t = 0.)
	N = div(length(u), 2)
	u1 =  @view (u[1:N])
	u2 =  @view (u[N+1:end])

	dest[1:N]     .= f.(u1, u2, Ref(p))
	dest[N+1:end] .= g.(u1, u2, Ref(p))

	return dest
end

function Fbr!(f, u, p, t = 0.)
	mul!(f, p.Δ, u)
	f .= f .+ NL(u, p)
end

function NL(u, p)
	out = similar(u)
	NL!(out, u, p)
	out
end

function Fbr(x, p, t = 0.)
	f = similar(x)
	Fbr!(f, x, p)
end

dNL(x, p, dx) = ForwardDiff.derivative(t -> NL(x .+ t .* dx, p), 0.)

function dFbr(x, p, dx)
	f = similar(dx)
	mul!(f, p.Δ, dx)

	nl = dNL(x, p, dx)
	f .= f .+ nl
end

Jbr(x, p) = sparse(ForwardDiff.jacobian(x -> Fbr(x, p), x))
# Jbr(sol0, par_br)

####################################################################################################
N = 100
	n = 2N
	lx = 3pi /2
	X = LinRange(-lx,lx, N)
	Δ = Laplacian(N, lx, :Neumann)
	D = 0.08
	par_br = (η = 1.0, a = -1., b = -3/2., H = 3.0, D = D, C = -0.6, Δ = blockdiag(D*Δ, Δ))

	u0 = cos.(2X)
	solc0 = vcat(u0, u0)
####################################################################################################
# eigls = DefaultEig()
eigls = EigArpack(0.5, :LM)
	optnewton = NewtonPar(eigsolver = eigls, verbose=true, maxIter = 3200, tol=1e-9)
	out, _, _ = @time newton(
		x -> Fbr(x, par_br),
		x -> Jbr(x, par_br),
		solc0, optnewton,
		normN = norminf)

		plot();plot!(X,out[1:N]);plot!(X,solc0[1:N], label = "sol0",line=:dash)

optcont = ContinuationPar(dsmax = 0.01, ds = -0.001, pMin = -1.8, detectBifurcation = 1, nev = 21, plotEveryNsteps = 50, newtonOptions = optnewton, maxSteps = 400)

	br, _ = @time continuation(
		(x, p) -> Fbr(x, @set par_br.C = p),
		(x, p) -> Jbr(x, @set par_br.C = p),
		solc0, -0.2,
		optcont;
		plot = true, verbosity = 2,
		printSolution = (x, p) -> norm(x, Inf),
		plotSolution = (x; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...))
####################################################################################################
ind_hopf = 1
	hopfpt = HopfPoint(br, ind_hopf)

# we get the parameters from the Hopf point
hopfpoint = HopfPoint(br, ind_hopf)
C_hopf = hopfpoint.p[1]
ωH = hopfpoint.p[end] |> abs
# number of time slices
M = 61

# we get the parameters from the Hopf point
orbitguess = zeros(2N, M)
vec_hopf = geteigenvector(optnewton.eigsolver,
			br.eig[br.bifpoint[ind_hopf].idx][2],
			br.bifpoint[ind_hopf].ind_bif-1)
vec_hopf ./=  norm(vec_hopf)
# vec_hopf is the eigenvector for the eigenvalues iω
phase = []
	for ii=1:M
		t = (ii-1)/(M-1)
		orbitguess[:, ii] .= real.(hopfpoint.u .+
				2.1 * vec_hopf * exp(-2pi * complex(0, 1) * (t - .2750)))
		push!(phase, dot(orbitguess[:, ii] - hopfpoint.u,real.(vec_hopf)))
	end
plot(phase)
orbitguess_f = vcat(vec(orbitguess), 2pi/ωH) |> vec

poTrap = p -> PeriodicOrbitTrapProblem(
	x ->  Fbr(x, @set par_br.C = p),
	x ->  Jbr(x, @set par_br.C = p),
	real.(vec_hopf),
	hopfpoint.u,
	M)

poTrap(-0.9)(orbitguess_f) |> plot
PALC.plotPeriodicPOTrap(orbitguess_f, N, M)

deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zero(orbitguess_f)])
# deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo_f])

opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 120)

outpo_f, _, flag = @time PALC.newton(
		poTrap(-0.9),
		orbitguess_f, opt_po, :FullLU;
		normN = norminf)
	flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", PALC.amplitude(outpo_f,2*N, M),"\n")
	PALC.plotPeriodicPOTrap(outpo_f, N, M)

eig = EigKrylovKit(tol= 1e-10, x₀ = rand(2N), verbose = 2, dim = 40)
# eig = EigArpack()
eig = DefaultEig()
optcontpo = ContinuationPar(dsmin = 0.001, dsmax = 0.015, ds= 0.01, pMin = -1.8, maxSteps = 140, newtonOptions = (@set opt_po.eigsolver = eig), nev = 25, precisionStability = 1e-7, detectBifurcation = 0)
	br_po, _ , _ = @time continuationPOTrap(poTrap,
		outpo_f, -0.87,
		optcontpo; verbosity = 2,
		plot = true,
		# callbackN = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", amplitude(x, n, M));true),
		plotSolution = (x;kwargs...) ->  heatmap!(reshape(x[1:end-1], 2*N, M)'; ylabel="time", color=:viridis, kwargs...),
		printSolution = (u, p) -> PALC.maximumPOTrap(u, N, M; ratio = 2),
		normC = norminf)

# branches = [br_po]
push!(branches, br_po)
plotBranch(vcat(branches, br), label = "");title!("")
####################################################################################################
# Period doubling
ind_pd = 1
vec_pd = geteigenvector(eig,
		br_po.eig[br_po.bifpoint[ind_pd].idx][2],
		br_po.bifpoint[ind_pd].ind_bif)

# orbitguess_f = br_po.bifpoint[1].u .+ 0.1 * real.(vec_pd)
# orbitguess_f[end] *= 2
#
# outpo_pd, _, _ = @time PALC.newton(
# 		poTrap(br_po.bifpoint[1].param),
# 		orbitguess_f, opt_po, :FullLU;
# 		normN = norminf)
# 	printstyled(color=:red, "--> T = ", outpo_pd[end], ", amplitude = ", PALC.amplitude(outpo_pd, 2N, M),"\n")
# 	PALC.plotPeriodicPOTrap(outpo_pd, N, M)
####################################################################################################
# shooting
orbitsection = Array(sol[:,[end]])
# orbitsection = orbitguess[:, 1]

initpo = vcat(vec(orbitsection), 3.)

PALC.plotPeriodicShooting(initpo[1:end-1], 1);title!("")

probSh = p -> ShootingProblem(u -> Fbr(u, p), p, prob_sp, ETDRK2(krylov=true),
		1, x -> PALC.sectionShooting(x, Array(sol[:,[end]]), p, Fbr); atol = 1e-14, rtol = 1e-14, dt = 0.1)

probSh = p -> ShootingProblem(u -> Fbr(u, p), p, prob, Rodas4P(),
		1, x -> sectionArray(x, Array(sol[:,[end]]), p))#; atol = 1e-12, rtol = 1e-10)

ls = GMRESIterativeSolvers(tol = 1e-7, N = length(initpo), maxiter = 50, verbose = false)
	# ls = GMRESKrylovKit{Float64}(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 120, linsolver = ls)
	# deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outposh, _, flag = @time newton(probSh(@set par_br.C = -0.86),
		initpo, optn;
		callbackN = (x, f, J, res, iteration; kwargs...) -> (@show x[end];true),
		normN = norminf)
	flag && printstyled(color=:red, "--> T = ", outposh[end], ", amplitude = ", PALC.getAmplitude(probSh(@set par_br.C = -0.86), outposh; ratio = 2),"\n")

	plot(initpo[1:end-1], label = "Init guess")
	plot!(outpo[1:end-1], label = "sol")

eig = EigKrylovKit(tol= 1e-12, x₀ = rand(2N), verbose = 2, dim = 40)
eig = DefaultEig()
optcontpo = ContinuationPar(dsmin = 0.0001, dsmax = 0.02, ds= 0.005, pMin = -1.8, maxSteps = 70, newtonOptions = (@set optn.eigsolver = eig), nev = 5, precisionStability = 1e-3, detectBifurcation = 1)
	br_po_sh, _ , _ = @time continuationPOShooting(
		p -> probSh(@set par_br.C = p),
		outposh, -0.86,
		optcontpo; verbosity = 2,
		plot = true,
		plotSolution = (x; kwargs...) -> PALC.plotPeriodicShooting!(x[1:end-1], 1; kwargs...),
		printSolution = (u, p) -> PALC.getMaximum(probSh(@set par_br.C = p), u; ratio = 2), normC = norminf)

# branches = [br_po_sh]
# push!(branches, br_po_sh)
# plotBranch(branches);title!("")

plotBranch(vcat(br_po_sh, br_po, br), label = "");title!("")

####################################################################################################
# shooting PD
prob_sp = SplitODEProblem(f1, f2, solc0, (0.0, 300.0), @set par_br.C = -1.32)

solpd = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1)
	# heatmap(sol.t, X, sol[N,:], color=:viridis, xlim=(20,280.0))

plot(solpd.t, solpd[N÷2, :], xlim=(290,296.2))

solpd.t[end-100]

orbitsectionpd = Array(solpd[:,end-100])
initpo_pd = vcat(vec(orbitsectionpd), 6.2)
PALC.plotPeriodicShooting(initpo_pd[1:end-1], 1);title!("")

ls = GMRESIterativeSolvers(tol = 1e-7, N = length(initpo), maxiter = 50, verbose = false)
	# ls = GMRESKrylovKit{Float64}(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 120, linsolver = ls)
	# deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outposh_pd, _, flag = @time newton(x -> probSh(@set par_br.C = -1.32)(x),
		x -> (dx -> probSh(@set par_br.C = -1.32)(x, dx)),
		initpo_pd, optn;
		callbackN = (x, f, J, res, iteration; kwargs...) -> (@show x[end];true),
		normN = norminf)
	flag && printstyled(color=:red, "--> T = ", outposh_pd[end], ", amplitude = ", PALC.getAmplitude(probSh(@set par_br.C = -0.86), outposh_pd; ratio = 2),"\n")

	plot(initpo[1:end-1], label = "Init guess")
	plot!(outposh_pd[1:end-1], label = "sol")

optcontpo = ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= 0.001, pMin = -1.8, maxSteps = 100, newtonOptions = (@set optn.eigsolver = eig), nev = 5, precisionStability = 1e-3, detectBifurcation = 1)
	br_po_sh_pd, _ , _ = @time continuationPOShooting(
		p -> probSh(@set par_br.C = p),
		outposh_pd, -1.32,
		optcontpo; verbosity = 2,
		plot = true,
		# callbackN = cb_ss,
		plotSolution = (x; kwargs...) -> PALC.plotPeriodicShooting!(x[1:end-1], 1; kwargs...),
		printSolution = (u, p) -> PALC.getMaximum(probSh(@set par_br.C = p), u; ratio = 2), normC = norminf)

plotBranch(vcat(br_po_sh, br, br_po_sh_pd), label = "");title!("")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, SparseArrays, Test

#=
In [36]: C = Matrix([[0], [30], [15]]) # node coordinates
In [37]: A = Matrix([P.subs({x: C[i,0]}).T for i in range(len(P))])
In [38]: N = P.T*A.inv()
In [39]: Me = integrate(N.T*N, (x, 0, 30))
In [40]: De = diag(*integrate(N, (x, 0, 30)))
In [41]: Me
Out[41]:
Matrix([
[ 4, -1,  2],
[-1,  4,  2],
[ 2,  2, 16]])
In [42]: De
Out[42]:
Matrix([
[5, 0,  0],
[0, 5,  0],
[0, 0, 20]])
=#

@testset "dirichlet problem in 1 dimension" begin

    time = 0.0

    element = Element(Seg2, (1, 2))
    X = Dict(1 => [0.0, 0.0], 2 => [6.0, 0.0])
    update!(element, "geometry", X)
    update!(element, "temperature 1", 0.0)

    problem1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
    problem1.properties.variational = true
    problem1.properties.dual_basis = false
    add_element!(problem1, element)
    assemble!(problem1, time)
    C1 = problem1.assembly.C1
    C2 = problem1.assembly.C2
    @test isapprox(C1, C2)
    @test isapprox(C1, [2.0 1.0; 1.0 2.0])

    problem2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
    problem2.properties.variational = true
    problem2.properties.dual_basis = true
    add_element!(problem2, element)
    assemble!(problem2, time)
    C1 = problem2.assembly.C1
    C2 = problem2.assembly.C2
    @test isapprox(C1, C2)
    @test isapprox(C1, [3.0 0.0; 0.0 3.0])

    element = Element(Seg3, (1, 2, 3))
    X = Dict(1 => [0.0, 0.0], 2 => [30.0, 0.0], 3 => [15.0, 0.0])
    update!(element, "geometry", X)
    update!(element, "temperature 1", 0.0)
    problem3 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
    problem3.properties.variational = true
    problem3.properties.dual_basis = false
    add_element!(problem3, element)
    assemble!(problem3, time)
    C1 = problem3.assembly.C1
    C2 = problem3.assembly.C2
    @test isapprox(C1, C2)
    @test isapprox(C1, [4.0 -1.0 2.0; -1.0 4.0 2.0; 2.0 2.0 16.0])

    problem4 = Problem(Dirichlet, "quadratic 2", 1, "temperature")
    problem4.properties.variational = true
    problem4.properties.dual_basis = true
    add_element!(problem4, element)
    assemble!(problem4, time)
    C1 = problem4.assembly.C1
    C2 = problem4.assembly.C2
    @test isapprox(C1, C2)
    @test isapprox(C1, [5.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 20.0])
end

#=
@testset "dirichlet problem using tri3 surface element" begin
    element = Tri3([1, 2, 3])
    element["geometry"] = Node[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
    element["temperature"] = 0.0
    problem = Problem(Dirichlet, "test problem", 1, "temperature")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, 1/24*[2 1 1; 1 2 1; 1 1 2])
end

@testset "dirichlet problem in 2 dimensions" begin
    element = Seg2([1, 2])
    element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
    element["displacement 1"] = 0.0
    element["displacement 2"] = 0.0
    problem = Problem(Dirichlet, "test problem", 2, "displacement")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    g = full(problem.assembly.g)
    @test isapprox(C1, C2)
    C1_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
    @test isapprox(C1, C1_expected)
    @test isapprox(g, [0.0, 0.0, 0.0, 0.0])
end

@testset "dirichlet problem in 2 dimensions, with 1 dof fixed" begin
    element = Seg2([1, 2])
    element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
    element["displacement 2"] = 0.0
    problem = Problem(Dirichlet, "test problem", 2, "displacement")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    g = full(problem.assembly.g)
    @test isapprox(C1, C2)
    C1_expected = 1/6*[
        0 0 0 0
        0 2 0 1
        0 0 0 0
        0 1 0 2]
    @test isapprox(C1, C1_expected)
    @test isapprox(g, [0.0, 0.0, 0.0, 0.0])
end
=#

@testset "test analytical boundary condition" begin
    X = Dict(1 => [0.0, 0.0],
             2 => [1.0, 0.0])
    element = Element(Seg2, (1, 2))
    update!(element, "geometry", X)
    function f(element, ip, time)
        x, y = element("geometry", ip, time)
        val = x*time
        @debug("analytical function called", ip, time, x, y, val)
        return val
    end
    update!(element, "displacement 1", 0.0)
    update!(element, "displacement 2", f)
    problem = Problem(Dirichlet, "test boundary", 2, "displacement")
    add_element!(problem, element)
    time = 0.0
    assemble!(problem, time)
    @test isapprox(problem.assembly.g, [0.0, 0.0, 0.0, 0.0])
    empty!(problem.assembly)
    time = 1.0
    assemble!(problem, time)
    g2 = Vector(problem.assembly.g, 4)
    C2 = Matrix(problem.assembly.C2, 4, 4)
    u = C2 \ g2
    @debug("displacement vector", u)
    @test isapprox(u, [0.0, 0.0, 0.0, 1.0])
end


text:
using GLM, DataFrames, StatsFuns
glm(@formula(y ~ 1), DataFrame(y = float(bitrand(10))), Bernoulli())

n = 2_500_000
srand(1234321)
const df2 = DataFrame(x1 = rand(Normal(), n),
                x2 = rand(Exponential(), n),
                ss = pool(rand(DiscreteUniform(50), n)));
const mf = ModelFrame(@formula(x1 ~ x1 + x2 + ss), df2)
const mm = ModelMatrix(mf)
const β = unshift!(rand(Normal(),52), 0.5); # "true" parameter values

## Create linear predictor and mean response

const y = [float(rand() < logistic(η)) for η in mm.m * β];        # simulate observed responses

gm6 = glm(mm.m, y, Bernoulli())
@time glm(mm.m, y, Bernoulli());
#@profile glm(mm.m, y, Binomial());
#using ProfileView
#ProfileView.view()


text:
# Automatically generated using Clang.jl wrap_c

# Skipping MacroDefinition: const RMW_SAFE_FWRITE_TO_STDERR = RCUTILS_SAFE_FWRITE_TO_STDERR
const rmw_initialize_error_handling_thread_local_storage = rcutils_initialize_error_handling_thread_local_storage
const rmw_set_error_state = rcutils_set_error_state

# Skipping MacroDefinition: RMW_CHECK_ARGUMENT_FOR_NULL ( argument , error_return_type ) RCUTILS_CHECK_ARGUMENT_FOR_NULL ( argument , error_return_type )
# Skipping MacroDefinition: RMW_CHECK_FOR_NULL_WITH_MSG ( value , msg , error_statement ) RCUTILS_CHECK_FOR_NULL_WITH_MSG ( value , msg , error_statement )
# Skipping MacroDefinition: RMW_SET_ERROR_MSG ( msg ) RCUTILS_SET_ERROR_MSG ( msg )
# Skipping MacroDefinition: RMW_SET_ERROR_MSG_WITH_FORMAT_STRING ( format_string , ... ) RCUTILS_SET_ERROR_MSG_WITH_FORMAT_STRING ( format_string , __VA_ARGS__ )

const rmw_error_is_set = rcutils_error_is_set
const rmw_get_error_state = rcutils_get_error_state
const rmw_get_error_string = rcutils_get_error_string
const rmw_reset_error = rcutils_reset_error
const rmw_error_string_t = rcutils_error_string_t
const rmw_error_state_t = rcutils_error_state_t
const rmw_context_impl_t = Cvoid

struct rmw_context_t
    instance_id::UInt64
    implementation_identifier::Cstring
    impl::Ptr{rmw_context_impl_t}
end

const rmw_init_options_impl_t = Cvoid

struct rmw_init_options_t
    instance_id::UInt64
    implementation_identifier::Cstring
    allocator::rcutils_allocator_t
    impl::Ptr{rmw_init_options_impl_t}
end

# Skipping MacroDefinition: RMW_STRINGIFY ( x ) RCUTILS_STRINGIFY ( x )

# Skipping MacroDefinition: const RMW_WARN_UNUSED = RCUTILS_WARN_UNUSED

struct rmw_names_and_types_t
    names::rcutils_string_array_t
    types::Ptr{rcutils_string_array_t}
end

const RMW_RET_OK = 0
const RMW_RET_ERROR = 1
const RMW_RET_TIMEOUT = 2
const RMW_RET_BAD_ALLOC = 10
const RMW_RET_INVALID_ARGUMENT = 11
const RMW_RET_INCORRECT_RMW_IMPLEMENTATION = 12
const rmw_ret_t = Int32
const rmw_get_zero_initialized_serialized_message = rcutils_get_zero_initialized_uint8_array
const rmw_serialized_message_init = rcutils_uint8_array_init
const rmw_serialized_message_fini = rcutils_uint8_array_fini
const rmw_serialized_message_resize = rcutils_uint8_array_resize
const rmw_serialized_message_t = rcutils_uint8_array_t
const RMW_GID_STORAGE_SIZE = 24

struct rmw_node_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    name::Cstring
    namespace_::Cstring
    context::Ptr{rmw_context_t}
end

struct rmw_publisher_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    topic_name::Cstring
end

struct rmw_subscription_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    topic_name::Cstring
end

struct rmw_service_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    service_name::Cstring
end

struct rmw_client_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    service_name::Cstring
end

struct rmw_guard_condition_t
    implementation_identifier::Cstring
    data::Ptr{Cvoid}
    context::Ptr{rmw_context_t}
end

struct rmw_subscriptions_t
    subscriber_count::Csize_t
    subscribers::Ptr{Ptr{Cvoid}}
end

struct rmw_services_t
    service_count::Csize_t
    services::Ptr{Ptr{Cvoid}}
end

struct rmw_clients_t
    client_count::Csize_t
    clients::Ptr{Ptr{Cvoid}}
end

struct rmw_guard_conditions_t
    guard_condition_count::Csize_t
    guard_conditions::Ptr{Ptr{Cvoid}}
end

struct rmw_wait_set_t
    implementation_identifier::Cstring
    guard_conditions::Ptr{rmw_guard_conditions_t}
    data::Ptr{Cvoid}
end

struct rmw_request_id_t
    writer_guid::NTuple{16, Int8}
    sequence_number::Int64
end

struct rmw_time_t
    sec::UInt64
    nsec::UInt64
end

@cenum(rmw_security_enforcement_policy_t,
    RMW_SECURITY_ENFORCEMENT_PERMISSIVE = 0,
    RMW_SECURITY_ENFORCEMENT_ENFORCE = 1,
)

struct rmw_node_security_options_t
    enforce_security::rmw_security_enforcement_policy_t
    security_root_path::Cstring
end

@cenum(rmw_qos_reliability_policy_t,
    RMW_QOS_POLICY_RELIABILITY_SYSTEM_DEFAULT = 0,
    RMW_QOS_POLICY_RELIABILITY_RELIABLE = 1,
    RMW_QOS_POLICY_RELIABILITY_BEST_EFFORT = 2,
)
@cenum(rmw_qos_history_policy_t,
    RMW_QOS_POLICY_HISTORY_SYSTEM_DEFAULT = 0,
    RMW_QOS_POLICY_HISTORY_KEEP_LAST = 1,
    RMW_QOS_POLICY_HISTORY_KEEP_ALL = 2,
)
@cenum(rmw_qos_durability_policy_t,
    RMW_QOS_POLICY_DURABILITY_SYSTEM_DEFAULT = 0,
    RMW_QOS_POLICY_DURABILITY_TRANSIENT_LOCAL = 1,
    RMW_QOS_POLICY_DURABILITY_VOLATILE = 2,
)

struct rmw_qos_profile_t
    history::rmw_qos_history_policy_t
    depth::Csize_t
    reliability::rmw_qos_reliability_policy_t
    durability::rmw_qos_durability_policy_t
    avoid_ros_namespace_conventions::Bool
end

struct rmw_gid_t
    implementation_identifier::Cstring
    data::NTuple{24, UInt8}
end

struct rmw_message_info_t
    publisher_gid::rmw_gid_t
    from_intra_process::Bool
end

@cenum(RWM_PUBLIC_TYPE,
    RMW_LOG_SEVERITY_DEBUG = 10,
    RMW_LOG_SEVERITY_INFO = 20,
    RMW_LOG_SEVERITY_WARN = 30,
    RMW_LOG_SEVERITY_ERROR = 40,
    RMW_LOG_SEVERITY_FATAL = 50,
)

const rmw_log_severity_t = RWM_PUBLIC_TYPE
const RMW_TOPIC_VALID = 0
const RMW_TOPIC_INVALID_IS_EMPTY_STRING = 1
const RMW_TOPIC_INVALID_NOT_ABSOLUTE = 2
const RMW_TOPIC_INVALID_ENDS_WITH_FORWARD_SLASH = 3
const RMW_TOPIC_INVALID_CONTAINS_UNALLOWED_CHARACTERS = 4
const RMW_TOPIC_INVALID_CONTAINS_REPEATED_FORWARD_SLASH = 5
const RMW_TOPIC_INVALID_NAME_TOKEN_STARTS_WITH_NUMBER = 6
const RMW_TOPIC_INVALID_TOO_LONG = 7
const RMW_TOPIC_MAX_NAME_LENGTH = (Cuint(255)    # impl constraint
                                 - Cuint(8))     # reserved for prefixes
const RMW_NAMESPACE_VALID = 0
const RMW_NAMESPACE_INVALID_IS_EMPTY_STRING = 1
const RMW_NAMESPACE_INVALID_NOT_ABSOLUTE = 2
const RMW_NAMESPACE_INVALID_ENDS_WITH_FORWARD_SLASH = 3
const RMW_NAMESPACE_INVALID_CONTAINS_UNALLOWED_CHARACTERS = 4
const RMW_NAMESPACE_INVALID_CONTAINS_REPEATED_FORWARD_SLASH = 5
const RMW_NAMESPACE_INVALID_NAME_TOKEN_STARTS_WITH_NUMBER = 6
const RMW_NAMESPACE_INVALID_TOO_LONG = 7
const RMW_NAMESPACE_MAX_LENGTH = RMW_TOPIC_MAX_NAME_LENGTH - UInt32(2)
const RMW_NODE_NAME_VALID = 0
const RMW_NODE_NAME_INVALID_IS_EMPTY_STRING = 1
const RMW_NODE_NAME_INVALID_CONTAINS_UNALLOWED_CHARACTERS = 2
const RMW_NODE_NAME_INVALID_STARTS_WITH_NUMBER = 3
const RMW_NODE_NAME_INVALID_TOO_LONG = 4
const RMW_NODE_NAME_MAX_NAME_LENGTH = 255

# Skipping MacroDefinition: RMW_EXPORT __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: RMW_PUBLIC __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: RMW_LOCAL __attribute__ ( ( visibility ( "hidden" ) ) )

const RMW_AVOID_MEMORY_ALLOCATION = 0

# Skipping MacroDefinition: RMW_TRY_PLACEMENT_NEW ( Destination , BufferForNew , FailureAction , Type , ... ) try { Destination = new ( BufferForNew ) Type ( __VA_ARGS__ ) ; } catch ( const std :: exception & exception ) { RMW_SET_ERROR_MSG ( ( std :: string ( "caught C++ exception " ) + rmw :: impl :: cpp :: demangle ( exception ) + " constructing " # Type ": " + exception . what ( ) ) . c_str ( ) ) ; FailureAction ; } catch ( ... ) { RMW_SET_ERROR_MSG ( "caught unknown C++ exception constructing " # Type ) ; FailureAction ; }
# Skipping MacroDefinition: RMW_TRY_DESTRUCTOR ( Statement , Type , FailureAction ) try { Statement ; } catch ( const std :: exception & exception ) { RMW_SET_ERROR_MSG ( ( std :: string ( "caught C++ exception in destructor of " # Type ": " ) + rmw :: impl :: cpp :: demangle ( exception ) + ": " + exception . what ( ) ) . c_str ( ) ) ; FailureAction ; } catch ( ... ) { RMW_SET_ERROR_MSG ( "caught unknown C++ exception in destructor of " # Type ) ; FailureAction ; }
# Skipping MacroDefinition: RMW_TRY_DESTRUCTOR_FROM_WITHIN_FAILURE ( Statement , Type ) try { Statement ; } catch ( const std :: exception & exception ) { std :: stringstream ss ; ss << "caught C++ exception in destructor of " # Type " while handling a failure: " << rmw :: impl :: cpp :: demangle ( exception ) << ": " << exception . what ( ) << ", at: " << __FILE__ << ":" << __LINE__ << '\n' ; ( std :: cerr << ss . str ( ) ) . flush ( ) ; } catch ( ... ) { std :: stringstream ss ; ss << "caught unknown C++ exception in destructor of " # Type << " while handling a failure at: " << __FILE__ << ":" << __LINE__ << '\n' ; ( std :: cerr << ss . str ( ) ) . flush ( ) ; }
# Skipping MacroDefinition: RMW_CHECK_TYPE_IDENTIFIERS_MATCH ( ElementName , ElementTypeID , ExpectedTypeID , OnFailure ) { if ( ElementTypeID != ExpectedTypeID ) { int __bytes_that_would_have_been_written = rcutils_snprintf ( NULL , 0 , # ElementName " implementation '%s'(%p) does not match rmw implementation '%s'(%p)" , ElementTypeID , reinterpret_cast < const void * > ( ElementTypeID ) , ExpectedTypeID , reinterpret_cast < const void * > ( ExpectedTypeID ) ) ; if ( __bytes_that_would_have_been_written < 0 ) { RMW_SET_ERROR_MSG ( "RMW_CHECK_TYPE_IDENTIFIERS_MATCH(): rcutils_snprintf(NULL, 0, ...) failed" ) ; OnFailure ; } else { char * __msg = reinterpret_cast < char * > ( rmw_allocate ( __bytes_that_would_have_been_written + 1 ) ) ; if ( NULL == __msg ) { RMW_SET_ERROR_MSG ( "RMW_CHECK_TYPE_IDENTIFIERS_MATCH(): rmw_allocate() failed" ) ; } else { int ret = rcutils_snprintf ( __msg , __bytes_that_would_have_been_written + 1 , # ElementName " implementation '%s'(%p) does not match rmw implementation '%s'(%p)" , ElementTypeID , reinterpret_cast < const void * > ( ElementTypeID ) , ExpectedTypeID , reinterpret_cast < const void * > ( ExpectedTypeID ) ) ; if ( ret < 0 ) { RMW_SET_ERROR_MSG ( "RMW_CHECK_TYPE_IDENTIFIERS_MATCH(): rcutils_snprintf() failed" ) ; } else { RMW_SET_ERROR_MSG ( __msg ) ; } } rmw_free ( __msg ) ; OnFailure ; } } }


text:
# TODO: DEBUG!
# - [ ] Is it perhaps needed to sample the other parameters multiple times
#       after updating beta? Especially when the value of beta has changed?
# - [ ] Is it necessary to fit the model with fixed beta first?
# - [ ] Perhaps fit two models; one with beta=0, one with beta=1,
#       then sample beta after somehow?
# - [ ] What happens if I don't marginalize over lambda when updating beta?
# - NOTE: nu is a trouble-maker. Fixing it solves most of the problems!!!

ENV["GKSwstype"] = "nul"  # For StatsPlots

import Pkg; Pkg.activate("../../../")

using CytofDensityEstimation
using LaTeXStrings
using Distributions
using StatsPlots
using HypothesisTests
import Random
using BSON

const CDE = CytofDensityEstimation

resultsdir = "results"
imgdir = "$(resultsdir)/img"
mkpath(imgdir)

Random.seed!(2);
simdata = CDE.Model.generate_samples(NC=1000, NT=1000, gammaC=0.3, gammaT=0.2,
                                     etaC=[.5, .5, 0], etaT=[.5, .2, .3],
                                     # etaC=[.5, .5, 0], etaT=[.5, .5, .0],
                                     loc=[-1, 1, 3.], scale=[1, 1, 1]/10,
                                     df=[15, 30, 10], skew=[-20, -5, 0.])
yC, yT = simdata[:yC], simdata[:yT]

# Plot data.
density(yC[isfinite.(yC)],  lw=3, label=L"y_C", color=:blue)
density!(yT[isfinite.(yT)], lw=3, label=L"y_T", color=:red)
savefig(joinpath(imgdir, "data.pdf"))
closeall()

Random.seed!(1)
K = 6
data = CDE.Model.Data(yC, yT)
prior_mu = let
  yfinite = [data.yC_finite; data.yT_finite]
  Normal(mean(yfinite), std(yfinite)*3)
end
prior = CDE.Model.Prior(K, mu=prior_mu, nu=LogNormal(3, .01), p=Beta(.1, .9))
state = CDE.Model.State(data, prior)
tuners = CDE.Model.Tuners(K)

println("Priors:\n", prior)

# state.beta=0
# chain0, laststate, summarystats = CDE.Model.fit(state, data, prior, tuners,
#                                                 nsamps=[1000], nburn=2000,
#                                                 fix=[:beta])
#
# state.beta=1
# chain1, laststate, summarystats = CDE.Model.fit(state, data, prior, tuners,
#                                                 nsamps=[1000], nburn=2000,
#                                                 fix=[:beta])

# state.beta = true
# state.gammaC = simdata[:gammaC]
# state.gammaT = simdata[:gammaT]
# state.etaC .= [simdata[:etaC]; zeros(2)]
# state.etaT .= [simdata[:etaT]; zeros(2)]
# state.mu .= [simdata[:loc]; rand(prior.mu, 2)]
# state.nu .= [simdata[:df]; rand(prior.nu, 2)]
# state.vC .= 1
# state.vT .= 1
# state.psi .= let
#   _psi = CDE.Util.toaltskew.(simdata[:scale], simdata[:skew])
#   [_psi; rand(prior.psi, 2)]
# end
# state.omega .= let
#   _omega = CDE.Util.toaltscale.(simdata[:scale], simdata[:skew]) .^ 2
#   [_omega; rand(prior.omega, 2)]
# end

# fix=[:gammaC, :gammaT, :mu, :nu, :omega, :psi, :etaC, :etaT]
# fix=[:gammaC, :gammaT, :mu, :nu, :omega, :psi, :v]
# fix=Symbol[]  # TODO: Test this case.
# fix=[:gammaC, :gammaT, :nu, :omega, :psi, :etaC, :etaT]
# fix=[:gammaC, :gammaT, :etaC, :etaT, :nu]
# fix=[:gammaC, :gammaT, :etaC, :etaT, :v]
# fix=[:v]
# fix=Symbol[]
fix=Symbol[:nu]

init = deepcopy(state)
chain, laststate, summarystats = CDE.Model.fit(init, data, prior, tuners,
                                               nsamps=[1000], nburn=1000,
                                               fix=fix)

# Save results
BSON.bson("$(resultsdir)/results.bson",
          Dict(:chain=>chain, :laststate=>laststate, :summarystats=>summarystats))

# Load via:
# using BSON, CytofDensityEstimation
out = BSON.load("$(resultsdir)/results.bson")
chain = out[:chain]
laststate = out[:laststate]
summarystats = out[:summarystats]

ks_fit = ApproximateTwoSampleKSTest(yC, yT)
println(ks_fit)

CDE.Model.printsummary(chain, summarystats)

@time CDE.Model.plotpostsummary(chain, summarystats, yC, yT, imgdir,
                                simdata=simdata, bw_postpred=.1)

# Hard to recover (for K≥3):
# simdata = CDE.Model.generate_samples(NC=1000, NT=1000, gammaC=0.3, gammaT=0.2,
#                                      etaC=[.5, .5, 0], etaT=[.5, .2, .3],
#                                      loc=[-1, 1, 3.], scale=[1, 1, 1]/10,
#                                      df=[15, 30, 10], skew=[-20, -5, 0.])
# - nu, phi (psi). All other parameters are easily recovered.

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# yearmonthday is the opposite of totaldays
# taking Rata Die Day # and returning proleptic Gregorian date
@test Dates.yearmonthday(-306) == (0,2,29)
@test Dates.yearmonth(-306) == (0,2)
@test Dates.monthday(-306) == (2,29)
@test Dates.yearmonthday(-305) == (0,3,1)
@test Dates.yearmonth(-305) == (0,3)
@test Dates.monthday(-305) == (3,1)
@test Dates.yearmonthday(-2) == (0,12,29)
@test Dates.yearmonth(-2) == (0,12)
@test Dates.monthday(-2) == (12,29)
@test Dates.yearmonthday(-1) == (0,12,30)
@test Dates.yearmonth(-1) == (0,12)
@test Dates.monthday(-1) == (12,30)
@test Dates.yearmonthday(0) == (0,12,31)
@test Dates.yearmonth(-0) == (0,12)
@test Dates.monthday(-0) == (12,31)
@test Dates.yearmonthday(1) == (1,1,1)
@test Dates.yearmonth(1) == (1,1)
@test Dates.monthday(1) == (1,1)
# year, month, and day return the indivial components
# of yearmonthday, avoiding additional calculations when possible
@test Dates.year(-1) == 0
@test Dates.month(-1) == 12
@test Dates.day(-1) == 30
@test Dates.year(0) == 0
@test Dates.month(0) == 12
@test Dates.day(0) == 31
@test Dates.year(1) == 1
@test Dates.month(1) == 1
@test Dates.day(1) == 1
@test Dates.yearmonthday(730120) == (2000,1,1)
@test Dates.year(730120) == 2000
@test Dates.month(730120) == 1
@test Dates.day(730120) == 1

# Test totaldays and yearmonthday from January 1st of "from" to December 31st of "to"
# test_dates(-10000,10000) takes about 15 seconds
# test_dates(year(typemin(Date)),year(typemax(Date))) is full range
# and would take.......a really long time
function test_dates1(from,to)
    y = m = d = 0
    test_day = Dates.totaldays(from,1,1)
    for y in from:to
        for m = 1:12
            for d = 1:Dates.daysinmonth(y,m)
                days = Dates.totaldays(y,m,d)
                @test days == test_day
                @test (y,m,d) == Dates.yearmonthday(days)
                test_day += 1
            end
        end
    end
end
test_dates1(0,2100)

# Test year, month, day, hour, minute
function test_dates1()
    y = m = d = h = mi = 0
    for m = 1:12
        for d = 1:Dates.daysinmonth(y,m)
            for h = 0:23
                for mi = 0:59
                    dt = Dates.DateTime(y,m,d,h,mi)
                    @test y == Dates.year(dt)
                    @test m == Dates.month(dt)
                    @test d == Dates.day(dt)
                    @test d == Dates.dayofmonth(dt)
                    @test h == Dates.hour(dt)
                    @test mi == Dates.minute(dt)
                    @test (m,d) == Dates.monthday(dt)
                    #@test s == Dates.second(dt)
                    #@test ms == Dates.millisecond(dt)
                end
            end
        end
    end
end
test_dates1()

# Test second, millisecond
function test_dates2()
    y = m = d = h = mi = s = ms = 0
    for y in [-2013,-1,0,1,2013]
        for m in [1,6,12]
            for d in [1,15,Dates.daysinmonth(y,m)]
                for h in [0,12,23]
                    for s = 0:59
                        for ms in [0,1,500,999]
                            dt = Dates.DateTime(y,m,d,h,mi,s,ms)
                            @test y == Dates.year(dt)
                            @test m == Dates.month(dt)
                            @test d == Dates.day(dt)
                            @test h == Dates.hour(dt)
                            @test s == Dates.second(dt)
                            @test ms == Dates.millisecond(dt)
                        end
                    end
                end
            end
        end
    end
end
test_dates2()

function test_dates2(from,to)
    y = m = d = 0
    for y in from:to
        for m = 1:12
            for d = 1:Dates.daysinmonth(y,m)
                dt = Dates.Date(y,m,d)
                @test y == Dates.year(dt)
                @test m == Dates.month(dt)
                @test d == Dates.day(dt)
            end
        end
    end
end
test_dates2(0,2100)

# week function
# Tests from https://en.wikipedia.org/wiki/ISO_week_date
@test Dates.week(Dates.Date(2005,1,1)) == 53
@test Dates.week(Dates.Date(2005,1,2)) == 53
@test Dates.week(Dates.Date(2005,12,31)) == 52
@test Dates.week(Dates.Date(2007,1,1)) == 1
@test Dates.week(Dates.Date(2007,12,30)) == 52
@test Dates.week(Dates.Date(2007,12,31)) == 1
@test Dates.week(Dates.Date(2008,1,1)) == 1
@test Dates.week(Dates.Date(2008,12,28)) == 52
@test Dates.week(Dates.Date(2008,12,29)) == 1
@test Dates.week(Dates.Date(2008,12,30)) == 1
@test Dates.week(Dates.Date(2008,12,31)) == 1
@test Dates.week(Dates.Date(2009,1,1)) == 1
@test Dates.week(Dates.Date(2009,12,31)) == 53
@test Dates.week(Dates.Date(2010,1,1)) == 53
@test Dates.week(Dates.Date(2010,1,2)) == 53
@test Dates.week(Dates.Date(2010,1,2)) == 53
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=1999
dt = Dates.DateTime(1999,12,27)
dt1 = Dates.Date(1999,12,27)
check = (52,52,52,52,52,52,52,1,1,1,1,1,1,1,2,2,2,2,2,2,2)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2000
dt = Dates.DateTime(2000,12,25)
dt1 = Dates.Date(2000,12,25)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Test from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2030
dt = Dates.DateTime(2030,12,23)
dt1 = Dates.Date(2030,12,23)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2004
dt = Dates.DateTime(2004,12,20)
dt1 = Dates.Date(2004,12,20)
check = (52,52,52,52,52,52,52,53,53,53,53,53,53,53,1,1,1,1,1,1,1)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end

# Vectorized accessors
a = Dates.Date(2014,1,1)
dr = [a,a,a,a,a,a,a,a,a,a]
@test Dates.year.(dr) == repmat([2014],10)
@test Dates.month.(dr) == repmat([1],10)
@test Dates.day.(dr) == repmat([1],10)

a = Dates.DateTime(2014,1,1)
dr = [a,a,a,a,a,a,a,a,a,a]
@test Dates.year.(dr) == repmat([2014],10)
@test Dates.month.(dr) == repmat([1],10)
@test Dates.day.(dr) == repmat([1],10)
@test Dates.hour.(dr) == repmat([0],10)
@test Dates.minute.(dr) == repmat([0],10)
@test Dates.second.(dr) == repmat([0],10)
@test Dates.millisecond.(dr) == repmat([0],10)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module TestAqua

using Aqua
using TestFunctionRunner

test() = Aqua.test_all(TestFunctionRunner)

end  # module


text:
@inline function _getZetas(smcio::SMCIO{Particle}, p::Int64) where Particle
  @assert 1 <= p <= smcio.n
  p == smcio.n && return smcio.zetas
  return smcio.allZetas[p]
end

@inline function _getWs(smcio::SMCIO, hat::Bool, p::Int64)
  @assert 1 <= p <= smcio.n
  if hat
    p == smcio.n && return smcio.ws
    return smcio.allWs[p]
  else
    return smcio.internal.vecOnes
  end
end

"""
    eta(smcio::SMCIO{Particle}, f::F, hat::Bool, p::Int64) where {Particle, F<:Function}

Compute:
  - ```!hat```: ``\\eta^N_p(f)``
  - ```hat```:  ``\\hat{\\eta}_p^N(f)``
"""
function eta(smcio::SMCIO{Particle}, f::F, hat::Bool, p::Int64) where
  {Particle, F<:Function}
  zetas::Vector{Particle} = _getZetas(smcio, p)
  ws::Vector{Float64} = _getWs(smcio, hat, p)
  v = f(zetas[1]) * ws[1]
  for i = 2:smcio.N
    @inbounds v += f(zetas[i]) * ws[i]
  end
  return v / smcio.N
end

"""
    allEtas(smcio::SMCIO, f::F, hat::Bool) where F<:Function

Compute ```eta(smcio::SMCIO, f::F, hat::Bool, p)``` for p in {1, …, smcio.n}
"""
function allEtas(smcio::SMCIO, f::F, hat::Bool) where F<:Function
  T = typeof(f(smcio.zetas[1]) / smcio.N)
  result::Vector{T} = Vector{T}(undef, smcio.n)
  for p = 1:smcio.n
    @inbounds result[p] = eta(smcio, f, hat, p)
  end
  return result
end

"""
    slgamma(smcio::SMCIO, f::F, hat::Bool, p::Int64) where {Particle, F<:Function}

Compute:
- ```!hat```: ``(\\eta^N_p(f) \\geq 0, \\log |\\gamma^N_p(f)|)``
- ```hat```:  ``(\\hat{\\eta}^N_p(f) \\geq 0, \\log |\\hat{\\gamma}_p^N(f)|)``
The result is returned as a ```Tuple{Bool, Float64}```: the first component
represents whether the returned value is non-negative, the second is the
logarithm of the absolute value of the approximation.
"""
function slgamma(smcio::SMCIO, f::F, hat::Bool, p::Int64) where
  F<:Function
  @assert 1 <= p <= smcio.n "p must be between 1 and smcio.n"
  idx::Int64 = p - 1 + hat
  logval = idx == 0 ? 0.0 : smcio.logZhats[idx]
  v::Float64 = eta(smcio, f, hat, p)
  logval += log(abs(v))
  return (v >= 0, logval)
end

"""
    allGammas(smcio::SMCIO, f::F, hat::Bool) where F<:Function

Compute ```slgamma(smcio::SMCIO, f::F, hat::Bool, p)``` for p in {1, …, smcio.n}
"""
function allGammas(smcio::SMCIO, f::F, hat::Bool) where F<:Function
  result::Vector{Tuple{Bool, Float64}} =
    Vector{Tuple{Bool, Float64}}(undef, smcio.n)
  for p = 1:smcio.n
    @inbounds result[p] = slgamma(smcio, f, hat, p)
  end
  return result
end


text:
function solve()
    a = parse(Int, readline())
    a^3 + a^2 + a
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
"""
`Prog = run_speedy(NF,kwargs...)`
`Prog = run_speedy(kwargs...)`

Runs SpeedyWeather.jl with number format `NF` and any additional parameters in the keyword arguments
`kwargs...`. Any unspeficied parameters will use the default values as defined in `src/parameters.jl`."""
function run_speedy(::Type{NF}=Float64;         # number format, use Float64 as default
                    kwargs...                   # all additional non-default parameters
                    ) where {NF<:AbstractFloat}

    P = Params(NF=NF,kwargs...)
    C = Constants{NF}(P)
    G = GeoSpectral{NF}(P)
    B = Boundaries{NF}(P,G)

    Prog = initial_conditions(P,B,G)

    # TODO
    # Prog = PrognosticVars{T}()
    # Diag = DiagnosticVars{T}()

    M = ModelSetup(P,C,G,B)

    #time_stepping!(Prog,Diag,M)

    return Prog
end


text:
# ---
# title: 962. Maximum Width Ramp
# id: problem962
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/maximum-width-ramp/description/>
# hidden: true
# ---
#
# Given an array `A` of integers, a _ramp_  is a tuple `(i, j)` for which `i <
# j` and `A[i] <= A[j]`.  The width of such a ramp is `j - i`.
#
# Find the maximum width of a ramp in `A`.  If one doesn't exist, return 0.
#
#
#
# **Example 1:**
#
#
#
#     Input: [6,0,8,2,1,5]
#     Output: 4
#     Explanation:
#     The maximum width ramp is achieved at (i, j) = (1, 5): A[1] = 0 and A[5] = 5.
#
#
# **Example 2:**
#
#
#
#     Input: [9,8,1,0,1,9,4,0,4,1]
#     Output: 7
#     Explanation:
#     The maximum width ramp is achieved at (i, j) = (2, 9): A[2] = 1 and A[9] = 1.
#
#
#
#
# **Note:**
#
#   1. `2 <= A.length <= 50000`
#   2. `0 <= A[i] <= 50000`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
module OCIM1
#=
This module serves to load the OCIM1 matrix and grid.
These are loaded from the public, persistant URL in FigShare.
The Julia BSON format version of the OCIM1 was created with the code
in the GitHub repository https://github.com/briochemc/OceanCirculations.
=#

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using BSON                  # For saving circulation as BSON format
using Unitful               # for units
using Reexport
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

# OCIM1 URLs
function url(; version="CTL") # TODO add other OCIM1 circulations
    url_start = "https://files.figshare.com/"
    version == "CTL" ? "$(url_start)23420027/OCIM1_CTL.bson" :
    OCIM1versionerror(version)
end

OCIM1versionerror(version) = error("""`$version` is not a valid OCIM1 version name (from within AIBECS).

                                   Only valid version is `CTL`.

                                   Start an issue on GitHub if you want me to add other OCIM1 circulations.""")

# OCIM1 Hashes
function sha(; version="CTL")
    version == "CTL" ? "28d4290c04844180191e71d5e1bc832ff43c1127dd621b719c6622a2288d5e24" :
    OCIM1versionerror(version)
end

# Create registry entry for OCIM1 in BSON format
function register_OCIM1(; version="CTL")
    register(
        DataDep(
            "AIBECS-OCIM1_$version",
            """
            References:
            $(citation())
            """,
            url(version=version),
            sha(version=version),
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load

Returns the grid and the transport matrix (in that order).

### Usage

```
julia> grd, T = OCIM1.load()
```

See *DeVries and Primeau* (2011) and *DeVries* (2014) for more details.
"""
function load(; version="CTL")
    register_OCIM1(version=version)
    bson_file = @datadep_str string("AIBECS-OCIM1_$version/", "OCIM1_$version.bson")
    BSON.@load bson_file grid T
    @info """You are about to use the OCIM1_$version model.
          If you use it for research, please cite:

          $(citation())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "DeVries_Primeau_2011" and "DeVries_2014" keys.)
          """
    return grid, ustrip.(T)
end

citation() = """
             - DeVries, T., 2014: The oceanic anthropogenic CO2 sink: Storage, air‐sea fluxes, and transports over the industrial era, Global Biogeochem. Cycles, 28, 631–647, doi:10.1002/2013GB004739.
             - DeVries, T. and F. Primeau, 2011: Dynamically and Observationally Constrained Estimates of Water-Mass Distributions and Ages in the Global Ocean. J. Phys. Oceanogr., 41, 2381–2401, doi:10.1175/JPO-D-10-05011.1
             """

versions() = ["CTL"]

end # end module

export OCIM1


text:
module FunctorTest

using Funs, Functor
using Base.Test

function main()
    # Array
    let
        xs = [1,2,4,8,16]
        ys = [1,2,3,4,5]
        @test fmap(x->2*x, xs) == [2,4,8,16,32]
        @test fmap(+, xs, ys) == [2,4,7,12,21]
        x2s = [1 2 4; 8 16 32]
        y2s = [1 2 3; 4 5 6]
        @test fmap(x->2*x, x2s) == [2 4 8; 16 32 64]
        @test fmap(+, x2s, y2s) == [2 4 7; 12 21 38]
    end

    # Fun (Callable)
    let
        f = Fun{Float64}(x::Int -> float(x)+0.5)
        g = Fun{Complex{Float64}}(x::Float64 -> x+im)
        @test fmap(g, f)(12) == 12.5+im
        f2 = Fun{Int}(x::Int -> x+1)
        g2 = Fun{Complex{Float64}}((x::Float64, y::Int) -> x+float(y)*im)
        @test fmap(g2, f, f2)(12,13) == 12.5+14.0im
    end

    # Set
    let
        xs = Set([1,2,4,8,16])
        @test fmap(x->2*x, xs) == Set([2,4,8,16,32])
    end

    # Tuple
    let
        xs = (1,2,4,8,16)
        ys = (1,2,3,4,5)
        @test fmap(x->2*x, xs) == (2,4,8,16,32)
        @test fmap(+, xs, ys) == (2,4,7,12,21)
    end
end

end


text:
using SDELab2 # SDELab2.jl should be findable via LOAD_PATH
println("run.jl")

function drift_f(u,r)
  return (r*u)
end
#
function diff_g_mat(u,dw,alpha)
  return alpha*(u.*dw)
end
#
function diff_g_vecs(alpha)
  m=2
  A=Array(Function,m)
  g1(u)=[u[1]*alpha,0]
  g2(u)=[0,u[2]*alpha]
  A[1]=g1
  A[2]=g2
  return A
end
#
function gbm_exact(y0,r,sig,t,w,ito_flag)
  if ito_flag==true
    out=exp((r-0.5*sig^2)*t+sig*w).*y0
  else
    out=exp(r*t+sig*w).*y0
  end
end
#
function get_error(dt,solver,test_case)
  r=test_case["r"]
  sig=test_case["sig"]
  noise_type=test_case["noise_type"]
  d=2
  m=2
  tf=1
  u=test_case["initial_data"]
  fcn=set_fcn(d,m,x->drift_f(x,r),(x,dw)->diff_g_mat(x,dw,sig),diff_g_vecs(sig), noise_type)
  opt=set_opt(dt,solver)
  t,y,W=sde_strong_solution(fcn, [0,tf], u, opt)
  exact=gbm_exact(u,r,sig,t[end],W,opt["is_ito"])
  yf=reshape(y[end,:],fcn.d)
  return norm(exact-yf)/norm(exact)# relative error
end
test_case1=Dict("r"=>-0.1, "sig"=>0.5, "noise_type"=>2, "initial_data"=>[1.,1.])
#

#fcn=set_fcn(2,2,x->drift_f(x,r),(x,dw)->diff_g_mat(x,dw,sig),diff_g_vecs(1), 2)
#opt=set_opt(0.01,"SIBDF")
#opt["nlsolve_param"].xtol=1e-5

#t,y,W=sde_strong_solution(fcn, linspace(0.0,1.0,10000), u, opt)
#exact=gbm_exact(u,r,sig,t[end],W,opt["is_ito"])
#
#plot(t,y)
#print(W,"u=",u)
dt_vals=[0.01,0.005,0.002,0.001,0.0005,1e-4]
error=similar(dt_vals)
all_error=Dict()
M=10
for key in keys(SDELab2.method_code)
  print("Integrator:", method_code[key])
  for i in 1:length(dt_vals)
    sum=0
    for m in 1:M
      err=get_error(dt_vals[i],key,test_case1)
      sum+=err^2
   end
  error[i]=sqrt(sum/M)
  end
  all_error[key]=deepcopy(error)
end
using PyPlot
for key in keys(method_code)
  loglog(dt_vals,all_error[key],label=method_code[key])
  println(method_code[key],all_error[key])
end
legend()


text:
export FinVolMesh, add_point!, add_face!, add_volume!, process_geometry!, mesh_smoothing!

using LinearAlgebra

"""
Mutable struct to define a mesh

* `ndim_space`, `ndim_mesh`: numbers of dimensions for the space in which the mesh is inserted and for the mesh
    itself, respectively. A 2D, surface mesh defined as the surface of a sphere, for example, would yield `ndim_mesh=2, ndim_space=3`
* `points`: array of arrays of floats for corner points
* `point_connectivity`: array of arrays of ints, indicating points with which each point shares
    a face. Used for mesh smoothing
* `point_fixed`: array of bools. Indicates whether or not each point in the mesh should be fixed when
    performing mesh smoothing
* `point_weight`: weights for each point when performing elliptic mesh smoothing. Array of floats
* `faces`: array of arrays of ints for faces (point indices)
* `face_centers`: array of arrays of floats for face centers
* `face_normals`: array of arrays of floats indicating face normals (with norm equal to the face's area)
* `face_normal_unit`: the same as `face_normals`, but with unit norm
* `face_connectivity`: array of arrays of ints indicating, respectively, the volume behind the face and the volume ahead of it
    (or more specifically their indices)
* `face_dists`: distances between the face and its adjacent volumes, in the same format as `face_connectivity`
* `volumes`: array of arrays of ints for faces (face indices)
* `volume_centers`: array of arrays of floats with positions for volume centers
* `volume_volumes`: volume of each finite volume (array of floats)
"""
mutable struct FinVolMesh{N}
    ndim_space::Int64
    ndim_mesh::Int64
    points::Vector{Vector{Float64}}
    point_connectivity::Vector{Vector{Int64}}
    point_fixed::Vector{Bool}
    point_weight::Vector{Float64}
    faces::Vector{Vector{Int64}}
    face_centers::Vector{Vector{Float64}}
    face_normals::Vector{Vector{Float64}}
    face_normal_unit::Vector{Vector{Float64}}
    face_connectivity::Vector{Vector{Int64}}
    face_dists::Vector{Vector{Float64}}
    volumes::Vector{Vector{Int64}}
    volume_centers::Vector{Vector{Float64}}
    volume_volumes::Vector{Float64} # duh
end

"""
Jacobian space volume
"""
function Jac_volume(Jac::Matrix{Float64})
    v=1.0
    a=1.0
    d=1.0

    for i=1:size(Jac, 2)
        for j=1:(i-1)
            Jac[:, i].-=dot(Jac[:, i], Jac[:, j]).*Jac[:, j]
        end

        locn=norm(Jac[:, i])
        Jac[:, i]./=(locn+1e-18)

        if i==size(Jac, 2)
            a=v
            d=locn
        end

        v*=locn/i
    end

    Jac[:, end]*a, v, d
end

"""
Function to process the geometry of a finite volume mesh and obtain pertinent geometrical properties

* `msh`: the mesh at hand
"""
function process_geometry!(msh::FinVolMesh{N}) where N
    for i=1:length(msh.points)
        msh.point_connectivity[i]=[]
    end

    for i=1:length(msh.faces)
        msh.face_centers[i].=0.0
        msh.face_normals[i].=0.0
        msh.face_normal_unit[i].=0.0
        msh.face_dists[i].=0.0

        msh.face_connectivity[i].=0
    end

    msh.volume_volumes.=0.0

    for i=1:length(msh.faces)
        for ptind in msh.faces[i]
            msh.face_centers[i].+=msh.points[ptind]
        end

        msh.face_centers[i]./=length(msh.faces[i])

        for j=1:length(msh.faces[i])
            p1=msh.faces[i][j]
            p2=msh.faces[i][j%length(msh.faces[i])+1]

            if !(p1 in msh.point_connectivity[p2])
                push!(msh.point_connectivity[p2], p1)
            end

            if !(p2 in msh.point_connectivity[p1])
                push!(msh.point_connectivity[p1], p2)
            end
        end
    end

    for i=1:length(msh.volumes)
        for fcind in msh.volumes[i]
            msh.volume_centers[i].+=msh.face_centers[fcind]
        end

        msh.volume_centers[i]./=length(msh.volumes[i])
    end

    Jac=Matrix{Float64}(undef, N, msh.ndim_mesh)

    for i=1:length(msh.volumes)
        for fc in msh.volumes[i]
            fcpts=msh.faces[fc]
            nfc=length(fcpts)

            freeside=1
            if msh.face_connectivity[fc][freeside]!=0
                freeside=2
            end

            if msh.face_connectivity[fc][freeside]!=0
                throw(
                    error(
                        "process_geometry!:ERROR:face $(fc) connected to more than two volumes"
                    )
                )
            end

            for j=1:(nfc-msh.ndim_mesh+1)
                for k=1:(msh.ndim_mesh-1)
                    kind=k+j

                    Jac[:, k].=msh.points[fcpts[kind]].-msh.points[fcpts[1]]
                end

                if freeside==1
                    Jac[:, msh.ndim_mesh].=msh.face_centers[fc].-msh.volume_centers[i]
                else
                    Jac[:, msh.ndim_mesh].=msh.volume_centers[i].-msh.face_centers[fc]
                end

                shat, v, d=Jac_volume(Jac)

                msh.face_normals[fc].+=shat
                msh.face_dists[fc][freeside]=d
                msh.face_connectivity[fc][freeside]=i

                msh.volume_volumes[i]+=v
            end

            if freeside==2
                msh.face_normals[fc]./=2
            end
        end
    end

    for fc=1:length(msh.faces)
        msh.face_normal_unit[fc].=msh.face_normals[fc]/norm(msh.face_normals[fc])
    end
end

"""
Constructor for an empty mesh

* `ndim_mesh`: number of dimensions for the mesh, if differnet from space dimensionality.
   A 2D, surface mesh defined as the surface of a sphere, for example, would yield `ndim_mesh=2, ndim_space=3`
"""
function FinVolMesh{N}(; ndim_mesh::Union{Int64, Nothing}=nothing) where N
    FinVolMesh{N}(
        N,
        something(ndim_mesh, N),
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        [],
        []
    )
end

"""
Add a point to a mesh

* `msh`: the mesh
* `pt`: a vector with the point at hand
* `fixed`: a boolean identifying whether the point should be fixed during mesh smoothing
* `weight`: a weight for mesh smoothing
* `geps`: if set to a non-zero value, is set as a minimum threshold for merging with other pre-existing points

* return: the new point's index
"""
function add_point!(msh::FinVolMesh{N}, pt::Vector{Float64}; fixed::Bool=false, weight::Float64=1.0, geps::Float64=0.0) where N
    @assert N==length(pt)

    if geps>0.0
        for (i, ept) in enumerate(msh.points)
            if norm(ept.-pt)<geps
                return i
            end
        end
    end

    push!(msh.points, pt)

    push!(msh.point_fixed, fixed)
    push!(msh.point_weight, weight)
    push!(msh.point_connectivity, [])

    length(msh.points)
end

function _is_face_eq(f1, f2)
    if length(f1)!=length(f2)
        return false
    end

    for p1 in f1
        if !(p1 in f2)
            return false
        end
    end

    true
end

"""
Add a face to a mesh

* `msh`: the mesh at hand
* `pts`: a vector of point indices
* `univocal`: if true (default), allows for merging with equal, previously set faces

* return: the new face's index
"""
function add_face!(msh::FinVolMesh{N}, pts::Vector{Int64}; univocal::Bool=true) where N
    for (i, fc) in enumerate(msh.faces)
        if _is_face_eq(fc, pts)
            return i
        end
    end

    push!(msh.faces, pts)

    push!(msh.face_centers, zeros(Float64, N))
    push!(msh.face_normals, zeros(Float64, N))
    push!(msh.face_normal_unit, zeros(Float64, N))
    push!(msh.face_connectivity, zeros(Int64, 2))
    push!(msh.face_dists, zeros(Float64, 2))

    length(msh.faces)
end

"""
Add volume to a mesh

* `msh`: the mesh
* `fcs`: vector of face indices

* return: an index for the volume
"""
function add_volume!(msh::FinVolMesh{N}, fcs::Vector{Int64}) where N
    push!(msh.volumes, fcs)

    push!(msh.volume_centers, zeros(Float64, N))
    push!(msh.volume_volumes, 0.0)

    length(msh.volumes)
end

"""
Function to apply elliptic equation mesh smoothing. Must be used after `process_geometry!`, which is
also automatically re-ran after changes are made

* `msh`: the mesh at hand
* `niter`: number of iterations (defaults to 100)
* `relaxation_parameter`: relaxation parameter for iteration
"""
function mesh_smoothing!(msh::FinVolMesh{N}; niter::Int64=100, relaxation_parameter::Float64=1.0) where N
    pos=zeros(Float64, N)
    w=0.0

    for nit=1:niter
        for (i, (pt, cons)) in enumerate(zip(msh.points, msh.point_connectivity))
            if !msh.point_fixed[i]
                w=msh.point_weight[i]
                pos.=pt.*w

                for c in cons
                    pos.+=msh.points[c].*msh.point_weight[c]
                    w+=msh.point_weight[c]
                end

                pos./=w

                msh.points[i].=relaxation_parameter.*pos.+(1.0-relaxation_parameter).*msh.points[i]

                w/=(length(cons)+1)
                msh.point_weight[i]=(1.0-relaxation_parameter)*msh.point_weight[i]+relaxation_parameter*w
            end
        end
    end

    process_geometry!(msh)
end


text:
function set_∂ψ!(∂ψ, ϕ, d::Domain, p::Params, ψ0::Function)
    Δx² = d.Δx^2
    Δy² = d.Δy^2
    c = CartesianIndices(d)
    l = LinearIndices(d)
    @. ∂ψ = 0
    # ∂x²
    for I = c[1,:,:]
        x, y, z = d[I-dx]
        ∂ψ[l[I]] += ψ0(z)/Δx²
    end
    for I = c[end,:,:]
        x, y, z = d[I+dx]
        ∂ψ[l[I]] += ψ0(z)/Δx²
    end
    # + ∂y²
    for I = c[:,1,:]
        x, y, z = d[I-dy]
        ∂ψ[l[I]] += ψ0(z)/Δy²
    end
    for I = c[:,end,:]
        x, y, z = d[I+dy]
        ∂ψ[l[I]] += ψ0(z)/Δy²
    end
    # * (1 + ∂z²ϕ)
    for I = c
        ∂ψ[l[I]] *= (1 + ∂z²(I,ϕ,d))
    end
    return ∂ψ
end

function set_linear_∂ψ!(∂ψ′, ϕ0, d::Domain, p::Params)
    return set_∂ψ!(∂ψ′, ϕ0, d, p, (z) -> 0.0)
end

function set_∂ψ!(∂ψ, ϕ, d::Domain, p::Params)
    return set_∂ψ!(∂ψ, ϕ, d, p, (z) -> background_ψ(z,p))
end

function set_∂ϕ!(∂ϕ, ψ, d::Domain, p::Params, ϕ0::Function)
    Δx² = d.Δx^2
    Δy² = d.Δy^2
    Δz² = d.Δz^2
    c = CartesianIndices(d)
    l = LinearIndices(d)
    @. ∂ϕ = 0
    # ∂z²
    for I = c[:,:,1]
        x′, y′, z′ = d[I]
        x, y, z = d[I-dz]
        ∂ϕ[l[I]] += (ϕ0(z) - ϕ0(z′))/Δz²
    end
    for I = c[:,:,end]
        x′, y′, z′ = d[I]
        x, y, z = d[I+dz]
        ∂ϕ[l[I]] += (ϕ0(z) - ϕ0(z′))/Δz²
    end
    # * (1 + ∇²ψ)
    for I = c
        ∂ϕ[l[I]] *= (1 + ∂x²(I,ψ,d) + ∂y²(I,ψ,d))
    end
    # + ∂x²
    for I = c[1,:,:]
        x, y, z = d[I-dx]
        ∂ϕ[l[I]] += ϕ0(z)/Δx²
    end
    for I = c[end,:,:]
        x, y, z = d[I+dx]
        ∂ϕ[l[I]] += ϕ0(z)/Δx²
    end
    # + ∂y²
    for I = c[:,1,:]
        x, y, z = d[I-dy]
        ∂ϕ[l[I]] += ϕ0(z)/Δy²
    end
    for I = c[:,end,:]
        x, y, z = d[I+dy]
        ∂ϕ[l[I]] += ϕ0(z)/Δy²
    end
    return ∂ϕ
end

function set_linear_∂ϕ!(∂ϕ′, ψ0, d::Domain, p::Params)
    return set_∂ϕ!(∂ϕ′, ψ0, d, p, (z) -> 0.0)
end

function set_∂ϕ!(∂ϕ, ψ, d::Domain, p::Params)
    return set_∂ϕ!(∂ϕ, ψ, d, p, (z) -> background_ϕ(z,p))
end

function set_∇²_∂!(∂, d::Domain, ϕ0::Function)
    Δx² = d.Δx^2
    Δy² = d.Δy^2
    c = CartesianIndices(d)
    l = LinearIndices(d)
    @. ∂ = 0
    # ∂x²
    for I = c[1,:,:]
        x, y, z = d[I-dx]
        ∂[l[I]] += ϕ0(z)/Δx²
    end
    for I = c[end,:,:]
        x, y, z = d[I+dx]
        ∂[l[I]] += ϕ0(z)/Δx²
    end
    # + ∂y²
    for I = c[:,1,:]
        x, y, z = d[I-dy]
        ∂[l[I]] += ϕ0(z)/Δy²
    end
    for I = c[:,end,:]
        x, y, z = d[I+dy]
        ∂[l[I]] += ϕ0(z)/Δy²
    end
end

function set_cv_∂ψ!(∂ψ, d::Domain, p::Params)
    return set_∇²_∂!(∂ψ, d, (z) -> background_ψ(z,p))
end

function set_cv_∂ϕ!(∂ϕ, d::Domain, p::Params)
    return set_∇²_∂!(∂ϕ, d, (z) -> background_ϕ(z,p))
end


text:
module Variational

using ..Core, ..Utilities
using DocStringExtensions: TYPEDEF, TYPEDFIELDS
using Distributions, Bijectors, DynamicPPL
using LinearAlgebra
using ..Turing: PROGRESS, Turing
using DynamicPPL: Model, SampleFromPrior, SampleFromUniform
using Random: AbstractRNG

using ForwardDiff
using Tracker

import ..Core: getchunksize, getADbackend

import AbstractMCMC
import ProgressLogging

using Requires
function __init__()
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        apply!(o, x, Δ) = Flux.Optimise.apply!(o, x, Δ)
        Flux.Optimise.apply!(o::TruncatedADAGrad, x, Δ) = apply!(o, x, Δ)
        Flux.Optimise.apply!(o::DecayedADAGrad, x, Δ) = apply!(o, x, Δ)
    end
    @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
        function Variational.grad!(
            vo,
            alg::VariationalInference{<:Turing.ZygoteAD},
            q,
            model,
            θ::AbstractVector{<:Real},
            out::DiffResults.MutableDiffResult,
            args...
        )
            f(θ) = if (q isa VariationalPosterior)
                - vo(alg, update(q, θ), model, args...)
            else
                - vo(alg, q(θ), model, args...)
            end
            y, back = Tracker.pullback(f, θ)
            dy = back(1.0)
            DiffResults.value!(out, y)
            DiffResults.gradient!(out, dy)
            return out
        end
    end
    @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        function Variational.grad!(
            vo,
            alg::VariationalInference{<:Turing.ReverseDiffAD{false}},
            q,
            model,
            θ::AbstractVector{<:Real},
            out::DiffResults.MutableDiffResult,
            args...
        )
            f(θ) = if (q isa VariationalPosterior)
                - vo(alg, update(q, θ), model, args...)
            else
                - vo(alg, q(θ), model, args...)
            end
            tp = Turing.Core.tape(f, θ)
            ReverseDiff.gradient!(out, tp, θ)
            return out
        end
        @require Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73" begin
            function Variational.grad!(
                vo,
                alg::VariationalInference{<:Turing.ReverseDiffAD{true}},
                q,
                model,
                θ::AbstractVector{<:Real},
                out::DiffResults.MutableDiffResult,
                args...
            )
                f(θ) = if (q isa VariationalPosterior)
                    - vo(alg, update(q, θ), model, args...)
                else
                    - vo(alg, q(θ), model, args...)
                end
                ctp = Turing.Core.memoized_tape(f, θ)
                ReverseDiff.gradient!(out, ctp, θ)
                return out
            end
        end
    end
end

export
    vi,
    ADVI,
    ELBO,
    elbo,
    TruncatedADAGrad,
    DecayedADAGrad

abstract type VariationalInference{AD} end

getchunksize(::Type{<:VariationalInference{AD}}) where AD = getchunksize(AD)
getADbackend(::VariationalInference{AD}) where AD = AD()

abstract type VariationalObjective end

const VariationalPosterior = Distribution{Multivariate, Continuous}

"""
    grad!(vo, alg::VariationalInference, q, model::Model, θ, out, args...)

Computes the gradients used in `optimize!`. Default implementation is provided for
`VariationalInference{AD}` where `AD` is either `ForwardDiffAD` or `TrackerAD`.
This implicitly also gives a default implementation of `optimize!`.

Variance reduction techniques, e.g. control variates, should be implemented in this function.
"""
function grad! end

"""
    vi(model, alg::VariationalInference)
    vi(model, alg::VariationalInference, q::VariationalPosterior)
    vi(model, alg::VariationalInference, getq::Function, θ::AbstractArray)

Constructs the variational posterior from the `model` and performs the optimization
following the configuration of the given `VariationalInference` instance.

# Arguments
- `model`: `Turing.Model` or `Function` z ↦ log p(x, z) where `x` denotes the observations
- `alg`: the VI algorithm used
- `q`: a `VariationalPosterior` for which it is assumed a specialized implementation of the variational objective used exists.
- `getq`: function taking parameters `θ` as input and returns a `VariationalPosterior`
- `θ`: only required if `getq` is used, in which case it is the initial parameters for the variational posterior
"""
function vi end

# default implementations
function grad!(
    vo,
    alg::VariationalInference{<:ForwardDiffAD},
    q,
    model,
    θ::AbstractVector{<:Real},
    out::DiffResults.MutableDiffResult,
    args...
)
    f(θ_) = if (q isa VariationalPosterior)
        - vo(alg, update(q, θ_), model, args...)
    else
        - vo(alg, q(θ_), model, args...)
    end

    chunk_size = getchunksize(typeof(alg))
    # Set chunk size and do ForwardMode.
    chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
    config = ForwardDiff.GradientConfig(f, θ, chunk)
    ForwardDiff.gradient!(out, f, θ, config)
end

function grad!(
    vo,
    alg::VariationalInference{<:TrackerAD},
    q,
    model,
    θ::AbstractVector{<:Real},
    out::DiffResults.MutableDiffResult,
    args...
)
    θ_tracked = Tracker.param(θ)
    y = if (q isa VariationalPosterior)
        - vo(alg, update(q, θ_tracked), model, args...)
    else
        - vo(alg, q(θ_tracked), model, args...)
    end
    Tracker.back!(y, 1.0)

    DiffResults.value!(out, Tracker.data(y))
    DiffResults.gradient!(out, Tracker.grad(θ_tracked))
end

"""
    optimize!(vo, alg::VariationalInference{AD}, q::VariationalPosterior, model, θ; optimizer = TruncatedADAGrad())

Iteratively updates parameters by calling `grad!` and using the given `optimizer` to compute
the steps.
"""
function optimize!(
    vo,
    alg::VariationalInference,
    q,
    model,
    θ::AbstractVector{<:Real};
    optimizer = TruncatedADAGrad(),
    progress = Turing.PROGRESS[],
    progressname = "[$(alg_str(alg))] Optimizing..."
)
    # TODO: should we always assume `samples_per_step` and `max_iters` for all algos?
    samples_per_step = alg.samples_per_step
    max_iters = alg.max_iters

    # TODO: really need a better way to warn the user about potentially
    # not using the correct accumulator
    if optimizer isa TruncatedADAGrad && θ ∉ keys(optimizer.acc)
        # this message should only occurr once in the optimization process
        @info "[$(alg_str(alg))] Should only be seen once: optimizer created for θ" objectid(θ)
    end

    diff_result = DiffResults.GradientResult(θ)

    # Create the progress bar.
    AbstractMCMC.@ifwithprogresslogger progress name=progressname begin
        # add criterion? A running mean maybe?
        for i in 1:max_iters
            grad!(vo, alg, q, model, θ, diff_result, samples_per_step)

            # apply update rule
            Δ = DiffResults.gradient(diff_result)
            Δ = apply!(optimizer, θ, Δ)
            @. θ = θ - Δ

            Turing.DEBUG && @debug "Step $i" Δ DiffResults.value(diff_result)

            # Update the progress bar.
            progress && ProgressLogging.@logprogress i/max_iters
        end
    end

    return θ
end

"""
    make_logjoint(model::Model; weight = 1.0)

Constructs the logjoint as a function of latent variables, i.e. the map z → p(x ∣ z) p(z).

The weight used to scale the likelihood, e.g. when doing stochastic gradient descent one needs to
use `DynamicPPL.MiniBatch` context to run the `Model` with a weight `num_total_obs / batch_size`.

## Notes
- For sake of efficiency, the returned function is closes over an instance of `VarInfo`. This means that you *might* run into some weird behaviour if you call this method sequentially using different types; if that's the case, just generate a new one for each type using `make_logjoint`.
"""
function make_logjoint(model::Model; weight = 1.0)
    # setup
    ctx = DynamicPPL.MiniBatchContext(
        DynamicPPL.DefaultContext(),
        weight
    )
    varinfo_init = Turing.VarInfo(model, ctx)

    function logπ(z)
        varinfo = VarInfo(varinfo_init, SampleFromUniform(), z)
        model(varinfo)

        return getlogp(varinfo)
    end

    return logπ
end

function logjoint(model::Model, varinfo, z)
    varinfo = VarInfo(varinfo, SampleFromUniform(), z)
    model(varinfo)

    return getlogp(varinfo)
end

# objectives
include("objectives.jl")

# optimisers
include("optimisers.jl")

# VI algorithms
include("advi.jl")

end


text:
using Knet
@show hidden = 100
@show winit = Gaussian(0,.01)
@show lr = 1.0
@show gclip = 5.0
@show batchsize = 5
maxnorm = zeros(2)
losscnt = zeros(2)

sgen1=SketchEngine(`xzcat sdfkl32KCsd_enTenTen12.vert.xz`; dict="tenten10k.dict", column=1)
sgen2=SketchEngine(`xzcat sdfkl32KCsd_enTenTen12.vert.xz`; dict="tentenpos58.dict", column=2)
data=TagData(sgen1,sgen2; batchsize=batchsize)
@show maxtoken(data,1)
@show vocab = maxtoken(data,2)

fnet = Net(lstm; out=hidden)
bnet = Net(lstm; out=hidden)
pnet = Net(add2; out=vocab, f=soft, ninputs=2)

model = Tagger(fnet, bnet, pnet)
setopt!(model; lr=lr)
train(model, data, softloss; gclip=gclip, maxnorm=maxnorm, losscnt=losscnt, lossreport=1000)

:ok


text:
using FlightSimulator
using Base.Test

path = abspath(joinpath(dirname(@__FILE__),"..","resources"))

#Create download path if required.
if !isdir(path)
    mkdir(path)
end

#download Kernels if they do not exist
if !isfile("$path/de430.bsp")
    download("https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp", "$path/de430.bsp")
end

if !isfile("$path/pck00010.tpc")
    download("https://naif.jpl.nasa.gov/pub/naif/pds/data/lro-l-spice-6-v1.0/lrosp_1000/data/pck/pck00010.tpc", "$path/pck00010.tpc")
end

#Load Kernels
LoadKernels("$path/de430.bsp")
LoadKernels("$path/pck00010.tpc")

#Set time after J2000
ep = 0.0

#Define a standard system
ListOfBodies = [Earth, Moon]
StandardSystem = PlanetarySystem(ListOfBodies)
StandardPos = [6378.0, 0.0, 0.0]
scpos = GetInertialFramePos(StandardSystem, Earth, StandardPos, ep)

#Another StandardSystem
ListOfBodies2 = [Sun,Earth,Moon]
StandardSystem2 = PlanetarySystem(ListOfBodies2)
scpos2 = GetInertialFramePos(StandardSystem2, Earth, StandardPos, ep)

#Perform tests
@testset "Environment Tests" begin
    @test FlightSimulator.parent(Earth) == FlightSimulator.parent(Moon)
    @test round.(Earth.Environment.GravityModel(Earth, [6378.0, 0.0, 0.0])*100000000)/100000 == [0.0, 0.0, 9.79871]
    @test PlanetarySystem(ListOfBodies) == StandardSystem
    @test PlanetarySystem(ListOfBodies2) == StandardSystem2
    @test StandardSystem.InertialFrame_Body == EarthMoonBarycenter
    @test StandardSystem2.InertialFrame_Body == SSB
    @test round.(GetGravAccel(StandardSystem, scpos, ep)*100000000)/100000 == [-1.72622, -9.64544, -1.0e-5]
    @test round.(GetGravAccel(StandardSystem2, scpos2, ep)*100000000)/100000 == [-1.76066, -9.47293, 0.06882]
end


text:
include("../src/LarSurf.jl")
using Plasm
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
# using LarSurf

block_size = [5, 5, 5]
margin_size = 0

## Artifical data
data3d = ones(Int16, 3, 3, 3)
data3d[2, 2, 2] = 10
# data3d

voxelsize_mm = [0.5, 1.0, 2.]
threshold=0

verts, trifaces = LarSurf.import_data3d(data3d, voxelsize_mm, 1)

println(verts)

println(trifaces)

println("======= To LAR ========")

V, EV, FE = LarSurf.to_lar(verts, trifaces)

println(V)
println(EV)
println(FE)

println("V: ", typeof(V), size(V))
println("EV: ", typeof(EV))
print("FV: ", typeof(FE))

@show EV;

# convert to list of lists
# trifaces_list = [trifaces[k, :] for k=1:size(trifaces,1)]
# this is not general
trifaces_list = [trifaces[k, 3:-1:1] for k=1:size(trifaces,1)]
Vt = permutedims(V, [2,1])

# Plasm.view(Lar.cuboid([1,1,1]))

Plasm.view(Vt, trifaces_list)
Plasm.viewexploded(Vt, trifaces_list)(2,2,2)


text:
# Tokenizer
# =========

struct Tokenizer
    machine::Machine
    actions_code::Vector{Tuple{Symbol,Expr}}
end

function compile(tokens::Pair{RegExp.RE,Expr}...; optimize::Bool=true)
    start = NFANode()
    final = NFANode()
    actions = Dict{Symbol,Action}()
    for i in 1:lastindex(tokens)
        # HACK: place token exit actions after any other actions
        action = Action(Symbol(:__token, i), 10000 - i)
        actions[action.name] = action
    end
    actions_code = Tuple{Symbol,Expr}[]
    for (i, (re, code)) in enumerate(tokens)
        re′ = RegExp.shallow_desugar(re)
        push!(re′.actions[:enter], :__token_start)
        name = Symbol(:__token, i)
        push!(re′.actions[:final], name)
        nfa = re2nfa(re′, actions)
        push!(start.edges, (Edge(eps), nfa.start))
        push!(nfa.final.edges, (Edge(eps), final))
        push!(actions_code, (name, code))
    end
    nfa = NFA(start, final)
    dfa = nfa2dfa(remove_dead_nodes(nfa))
    if optimize
        dfa = remove_dead_nodes(reduce_nodes(dfa))
    end
    return Tokenizer(dfa2machine(dfa), actions_code)
end

function generate_init_code(tokenizer::Tokenizer)
    # TODO: deprecate this?
    return generate_init_code(CodeGenContext(), tokenizer)
end

function generate_init_code(ctx::CodeGenContext, tokenizer::Tokenizer)
    quote
        $(ctx.vars.p)::Int = 1
        $(ctx.vars.p_end)::Int = 0
        $(ctx.vars.p_eof)::Int = -1
        $(ctx.vars.ts)::Int = 0
        $(ctx.vars.te)::Int = 0
        $(ctx.vars.cs)::Int = $(tokenizer.machine.start_state)
    end
end

function generate_exec_code(ctx::CodeGenContext, tokenizer::Tokenizer, actions=nothing)
    if actions == nothing
        actions = Dict{Symbol,Expr}()
    elseif actions == :debug
        actions = debug_actions(tokenizer.machine)
    elseif isa(actions, AbstractDict{Symbol,Expr})
        actions = copy(actions)
    else
        throw(ArgumentError("invalid actions argument"))
    end
    actions[:__token_start] = :($(ctx.vars.ts) = $(ctx.vars.p))
    for (i, (name, _)) in enumerate(tokenizer.actions_code)
        actions[name] = :(t = $(i); $(ctx.vars.te) = $(ctx.vars.p))
    end
    return generate_table_code(ctx, tokenizer, actions, true)
end

function generate_table_code(ctx::CodeGenContext, tokenizer::Tokenizer, actions::AbstractDict{Symbol,Expr}, check::Bool)
    action_dispatch_code, set_act_code = generate_action_dispatch_code(ctx, tokenizer.machine, actions)
    trans_table = generate_transition_table(tokenizer.machine)
    getbyte_code = generate_geybyte_code(ctx)
    cs_code = :($(ctx.vars.cs) = $(trans_table)[($(ctx.vars.cs) - 1) << 8 + $(ctx.vars.byte) + 1])
    eof_action_code = generate_eof_action_code(ctx, tokenizer.machine, actions)
    token_exit_code = generate_token_exit_code(tokenizer)
    return quote
        $(ctx.vars.mem) = $(SizedMemory)($(ctx.vars.data))
        t = 0
        ts = 0
        while p ≤ p_end && cs > 0
            $(getbyte_code)
            $(set_act_code)
            $(cs_code)
            $(action_dispatch_code)
            p += 1
        end
        if p > p_eof ≥ 0 && cs ∈ $(tokenizer.machine.final_states)
            $(eof_action_code)
            cs = 0
        elseif cs < 0
            p -= 1
        end
        if t > 0 && (cs ≤ 0 || p > p_end ≥ 0)
            $(token_exit_code)
            p = te + 1
            if cs != 0
                cs = $(tokenizer.machine.start_state)
            end
        end
    end
end

function generate_token_exit_code(tokenizer::Tokenizer)
    i = 0
    default = :()
    return foldr(default, reverse(tokenizer.actions_code)) do name_code, els
        _, code = name_code
        i += 1
        Expr(:if, :(t == $(i)), code, els)
    end
end