text:
using DataFrames, FreqTables, StatsBase, BenchmarkTools
srand(1); x = rand(1:100, 10^6); y = categorical(x); z = string.(x);

using FastGroupBy

@benchmark freqtable($x)
@benchmark fastby(sum, $x, $x |> length |> fcollect)
@benchmark sumby(x, x |> length |> fcollect)
@benchmark countmap($x)

@benchmark freqtable($y)
# @benchmark fastby(sum, $y, $y |> length |> fcollect)
@benchmark sumby(y, y |> length |> fcollect)
@benchmark countmap($y)

@benchmark freqtable($z)
@benchmark fastby(sum, $z, $z |> length |> fcollect)
@benchmark sumby(z, z |> length |> fcollect)
@benchmark countmap($z)

@benchmark by(DataFrame(x = $x), :x, nrow)

@benchmark fastby(sum, DataFrame(x = $x, weight = $x |> length |> fcollect), :x, :weight)

df1 = DataFrame(grps = rand(1:100, 1_000_000), val = rand(1_000_000))
# compute the difference between the number rows in that group and the mean of `val` in that group
res = fastby(val_grouped -> length(val_grouped) - mean(val_grouped), df1, :grps, :val)
# convert to dataframe
resdf = DataFrame(grps = keys(res) |> collect, len_minus_mean_val = values(res) |> collect)

# or use fastby


text:
struct DiscreteMeasure{X<:AbstractVector,P<:AbstractVector}
    xs::X
    ps::P

    function DiscreteMeasure{X,P}(xs::X, ps::P) where {X,P}
        length(xs) == length(ps) ||
            error("length of support `xs` and probabilities `ps` must be equal")
        return new{X,P}(xs, ps)
    end
end

"""
    DiscreteMeasure(xs::AbstractVector, ps::AbstractVector)

Construct a discrete measure with support `xs` and corresponding weights `ps`.
"""
function DiscreteMeasure(xs::AbstractVector, ps::AbstractVector)
    return DiscreteMeasure{typeof(xs),typeof(ps)}(xs, ps)
end

Base.length(μ::DiscreteMeasure) = length(μ.xs)

# initial gradient step (unregularized subgradient)
function gradient_step(c, τ, ν::DiscreteMeasure, x, ::Nothing)
    tmp = @. c((x,), ν.xs)
    i = argmin(tmp)
    z = τ .* ν.ps
    z[i] -= τ
    return z
end

# initial gradient step (regularized gradient)
function gradient_step(c, τ, ν::DiscreteMeasure, x, ε::Real)
    tmp = @. -c((x,), ν.xs) / ε
    LogExpFunctions.softmax!(tmp)
    z = @. τ * (ν.ps - tmp)
    return z
end

# gradient step (unregularized subgradient)
function gradient_step!(
    c,
    z::AbstractVector,
    τ,
    v::AbstractVector,
    ν::DiscreteMeasure,
    x,
    tmp::AbstractVector,
    ::Nothing;
    reset::Bool=false,
)
    @. tmp = c((x,), ν.xs) - v
    if reset
        @. z = τ * ν.ps
    else
        @. z += τ * ν.ps
    end
    z[argmin(tmp)] -= τ
    return z
end

# gradient step (regularized gradient)
function gradient_step!(
    c,
    z::AbstractVector,
    τ,
    v::AbstractVector,
    ν::DiscreteMeasure,
    x,
    tmp::AbstractVector,
    ε::Real;
    reset::Bool=false,
)
    @. tmp = (v - c((x,), ν.xs)) / ε
    LogExpFunctions.softmax!(tmp)
    if reset
        @. z = τ * (ν.ps - tmp)
    else
        @. z += τ * (ν.ps - tmp)
    end
    return z
end

@doc raw"""
    ctransform(c, v, x, ν, ε)

Compute the c-transform
```math
v^{c,ε}(x) = \begin{cases}
- ε \log\bigg(\int \exp{\Big(\frac{v_y - c(x, y)}{ε}\Big)} \, ν(\mathrm{d}y)\bigg) & \text{if } ε > 0,\\
\min_{y} c(x, y) - v_y & \text{otherwise}.
\end{cases}
```
"""
function ctransform(c, v::AbstractVector, x, ν::DiscreteMeasure, ::Nothing)
    return minimum(c(x, yᵢ) - vᵢ for (vᵢ, yᵢ) in zip(v, ν.xs))
end
function ctransform(c, v::AbstractVector, x, ν::DiscreteMeasure, ε::Real)
    t = LogExpFunctions.logsumexp(
        (vᵢ - c(x, yᵢ)) / ε + log(νᵢ) for (vᵢ, yᵢ, νᵢ) in zip(v, ν.xs, ν.ps)
    )
    return -ε * (t + 1)
end

ext: jl


text:
#hiJulia config file

include(Pkg.dir("JulieTest/src/watch.jl"))

# Configure the watcher

watch = ["src", "test"]
interval = 0.25 # in seconds

# Configure the test reporter
reporter = "spec"

# start hiJulia
julieTest()


text:
using RandomMatrices
using Base.Test

srand(1)
include("GaussianEnsembles.jl")
include("FormalPowerSeries.jl")
include("Haar.jl")
include("StochasticProcess.jl")

include("densities/Semicircle.jl")
include("densities/TracyWidom.jl")


text:
using NFFT

### performance test ###
function nfft_performance()
  let m = 3, sigma = 2.0
    @info "NFFT Performance Test 1D"
    let N = 2^19, M = N, x = rand(M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LUT, NFFT.FULL]
        @info "* precomputation = " pre
        @info "* initialization"
        @time p = plan_nfft(x, N, m, sigma, precompute=pre)

        @info "* adjoint"
        @time fApprox = nfft_adjoint(p, fHat, true)

        @info "* trafo"
        @time nfft(p, fApprox, true)
      end
    end

    @info "NFFT Performance Test 2D"
    let N = 1024, M = N*N, x2 = rand(2,M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LUT, NFFT.FULL]
        @info "* precomputation = " pre
        @info "* initialization"
        @time p = plan_nfft(x2, (N,N), m, sigma, precompute=pre)

        @info "* adjoint"
        @time fApprox = nfft_adjoint(p, fHat, true)

        @info "* trafo"
        @time nfft(p, fApprox, true)
      end
    end

    @info "NFFT Performance Test 3D"
    let N = 32, M = N*N*N, x3 = rand(3,M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LUT, NFFT.FULL]
        @info "* precomputation = " pre
        @info "* initialization"
        @time p = plan_nfft(x3, (N,N,N), m, sigma, precompute=pre)

        @info "* adjoint"
        @time fApprox = nfft_adjoint(p, fHat, true)

        @info "* trafo"
        @time nfft(p, fApprox, true)
      end
    end
  end
  return nothing
end

nfft_performance()


text:
# Unfortunately, while the standard ImageMagick TIFF reader can read the image data from ASPEX TIFF files,
# it can't read the spectrum data.  So I've written a lean TIFF tag reader to access the tags and
# to read the spectral data into a Spectrum object.

struct _ATField
    tagId::UInt16
    tagType::UInt16
    tagCount::UInt32
    tagData::Any

    function _ATField(ios, order)
        _TIFF_TYPES = (
            UInt8, # byte
            UInt8, # ASCII
            UInt16, # short
            UInt32, # long
            Rational{UInt32}, # rational
            Int8, # signed byte
            Int8, # undefined
            Int16, # signed short
            Int32, # signed long
            Rational{Int32}, # signed rational
            Float32, # float
            Float64, # double
        )
        tagId = order(read(ios, UInt16))
        tagType = order(read(ios, UInt16))
        tagCount = order(read(ios, UInt32))
        if tagCount * sizeof(_TIFF_TYPES[tagType]) <= 4
            tagData = [ order(read(ios, _TIFF_TYPES[tagType])) for _ in Base.OneTo(tagCount)]
            extra = 4 - sizeof(_TIFF_TYPES[tagType]) * tagCount
            seek(ios, position(ios) + extra) # make sure that we have read 4 bytes
        else
            tagOffset = order(read(ios, UInt32))
            ret = position(ios)
            seek(ios, tagOffset)
            tagData = [order(read(ios, _TIFF_TYPES[tagType])) for _ in Base.OneTo(tagCount)]
            seek(ios, ret)
        end
        return new(tagId, tagType, tagCount, tagData)
    end
end

struct _TIFFIFD
    ifdTags::Vector{_ATField}
    ifdNext::UInt
    function _TIFFIFD(ios, order::Function)
        sz = order(read(ios, UInt16))
        ifdTags = [_ATField(ios, order) for _ in Base.OneTo(sz)]
        ifdNext = order(read(ios, UInt32))
        return new(ifdTags, ifdNext)
    end
end

function Base.get(ifd::_TIFFIFD, id::Integer, def)
    idx = findfirst(atf -> atf.tagId == id, ifd.ifdTags)
    return isnothing(idx) ? def : ifd.ifdTags[idx]
end

struct _TIFFInternals
    ifds::Vector{_TIFFIFD}
    function _TIFFInternals(ios)
        magic = read(ios, UInt16)
        order = magic == 0x4D4D ? ntoh : (magic == 0x4949 ? ltoh : nothing)
        fortytwo = order(read(ios, UInt16))
        if isnothing(order) || (fortytwo ≠ 42)
            @error "This file is not a valid TIFF file."
        end
        if order ≠ ltoh
            @error "ASPEX Spectrum TIFF files are always little-endian ordered."
        end
        offset = order(read(ios, UInt32))
        ifds = _TIFFIFD[]
        while offset ≠ 0
            seek(ios, offset)
            ifd = _TIFFIFD(ios, order)
            push!(ifds, ifd)
            offset = ifd.ifdNext
        end
        return new(ifds)
    end
end

function _parseDesc(desc::AbstractString)
    number(v) = parse(Float64, match(r"([+-]?[0-9]+[.]?[0-9]*)", v)[1])
    date, time = missing, missing
    res, stagePos = Dict{Symbol,Any}(), Dict{Symbol,Float64}()
    for item in split(desc, "\n")
        (k, v) = split(item, "=")
        k = uppercase(k)
        try
            if k == "MAG"
                res[:ImageMag] = number(v)
                res[:FieldOfView] = (3.5 * 25.4 * 1.0e3) / number(v) # Assumed 3.5" image width in microns
            elseif k == "ZOOM"
                res[:ImageZoom] = number(v)
            elseif k == "ANALYSIS_DATE"
                date = Date(DateTime(v, "m/d/Y"))
            elseif k == "ANALYSIS_TIME"
                try
                    time = Time(DateTime(v, "H:M:S \\P\\M")) + Hour(12)
                catch
                    time = Time(DateTime(v, "H:M:S \\A\\M"))
                end # try
            elseif k == "OPERATOR"
                res[:Operator] = string(v)
            elseif k == "ACCELERATING_VOLTAGE"
                res[:BeamEnergy] = 1000.0 * number(v)
            elseif k == "WORKING_DISTANCE"
                res[:WorkingDistance] = 0.1 * number(v)
            elseif k == "COMMENT"
                res[:Name] = string(v)
            elseif k == "SAMPLE_DESCRIPTION"
                res[:Comment] = string(v)
            elseif k == "LIVE_TIME"
                res[:LiveTime] = number(v)
            elseif k == "ACQUISITION_TIME"
                res[:RealTime] = number(v)
            elseif k == "DETECTOR_TILT"
                res[:TakeOffAngle] = deg2rad(90.0 - number(v))
            elseif k == "STAGE_X"
                stagePos[:X] = 0.1 * number(v)
            elseif k == "STAGE_Y"
                stagePos[:Y] = 0.1 * number(v)
            elseif k == "STAGE_Z"
                stagePos[:Z] = 0.1 * number(v)
            elseif k == "STAGE_R"
                stagePos[:R] = deg2rad(number(v))
            elseif k == "STAGE_T"
                stagePos[:T] = deg2rad(number(v))
            elseif k == "STAGE_B"
                stagePos[:T] = deg2rad(number(v))
            end
        catch
            @info "Failed to parse $(k) => $(v)"
        end
    end
    if !isempty(stagePos)
        res[:StagePosition] = stagePos
    end
    if !ismissing(date)
        res[:AcquisitionTime] = ismissing(time) ? DateTime(date) : DateTime(date, time)
    end
    return res
end

function isAspexTIFF(filename)
    try
        open(filename, read = true) do ios
            return detectAspexTIFF(ios)
        end
    catch
        return false
    end
end

function detectAspexTIFF(ios)
    TIFF_SPECTRAL_DATA = 0x8352
    try
        seekstart(ios)
        ti = _TIFFInternals(ios)
        return any(ifd -> any(fld -> fld.tagId == TIFF_SPECTRAL_DATA, ifd.ifdTags), ti.ifds)
    catch
        return false
    end
end

function readAspexTIFF(
    file::AbstractString;
    withImgs = false,
    astype::Type{<:Real} = Float64,
)
    open(file, read = true) do ios
        return readAspexTIFF(ios, withImgs = withImgs, astype = astype)
    end
end

function readAspexTIFF(ios::IO; withImgs = false, astype::Type{<:Real} = Float64)
    # Special TIFF tags
    TIFF_SPECTRAL_DATA = 0x8352
    TIFF_SPECTRAL_XRES = 0x8353
    TIFF_SPECTRAL_XOFF = 0x8354
    TIFF_SPECTRAL_YRES = 0x8355
    TIFF_SPECTRAL_YOFF = 0x8356
    TIFF_IMAGE_DESCRIPTION = 270 # ASCII
    TIFF_SOFTWARE = 305 # ASCII
    floatonly(v) = parse(Float64, match(r"([+-]?[0-9]+[.]?[0-9]*)", v)[1])
    number(v) = parse(
        astype,
        match(
            astype isa Type{<:Integer} ? r"([+-]?[0-9]+)" : r"([+-]?[0-9]+[.]?[0-9]*)",
            v,
        )[1],
    )
    res = missing
    ti = _TIFFInternals(ios)
    for ifd in ti.ifds
        id = get(ifd, TIFF_IMAGE_DESCRIPTION, missing)
        sw = get(ifd, TIFF_SOFTWARE, missing)
        sp = get(ifd, TIFF_SPECTRAL_DATA, missing)
        sxr = get(ifd, TIFF_SPECTRAL_XRES, missing)
        sxo = get(ifd, TIFF_SPECTRAL_XOFF, missing)
        syr = get(ifd, TIFF_SPECTRAL_YRES, missing)
        syo = get(ifd, TIFF_SPECTRAL_YOFF, missing)
        if !ismissing(sp)
            @assert !(ismissing(sxr) || ismissing(sxo)) "X gain and offset data is missing from ASPEX TIFF spectrum"
            evperch = floatonly(ismissing(sxr) ? "10 eV/ch" : String(sxr.tagData))
            offset = floatonly(ismissing(sxo) ? "0 eV" : String(sxo.tagData))
            yoff = number(ismissing(syo) ? "0 counts" : String(syo.tagData))
            yres = number(ismissing(syr) ? "1 counts" : String(syr.tagData))
            energy = LinearEnergyScale(offset, evperch)
            data = map(i -> yoff + yres * convert(astype, i), sp.tagData)
            props = _parseDesc(String(id.tagData))
            if !ismissing(sw)
                props[:Instrument] = String(sw.tagData)
            end
            res = Spectrum(energy, data, props)
            break
        end
    end
    if withImgs && (!ismissing(res))
        try
            seekstart(ios)
            imgs = FileIO.load(Stream{format"TIFF"}(ios))
            nimgs = ndims(imgs)>2 ? size(imgs,3) : 1
            if haskey(res, :ImageMag) || haskey(res, :FieldOfView)
                fov = haskey(res, :ImageMag) ? (3.5*25.4)/res[:ImageMag] : 1.0e3*res[:FieldOfView] # X field-of-view in mm
                off = haskey(res, :StagePosition) ? (res[:StagePosition][:Y], res[:StagePosition][:X]) : (0.0, 0.0)
                ratio, pix = size(imgs,1) / size(imgs,2), fov / (size(imgs, 2)-1)
                ax = Axis{:x}((off[2]-0.5*fov)*mm:pix*mm:(off[2]-0.5*fov+pix*(size(imgs,2)-1))*mm)
                ay = Axis{:y}((off[1]+0.5*fov)*ratio*mm:-pix*mm:(off[1]+0.5*fov-pix*(size(imgs,1)-1))*ratio*mm)
                @assert length(ay)==size(imgs,1) && length(ax) == size(imgs,2)
                if nimgs == 2
                    # Macro image
                    res[Symbol("Image2")] = AxisArray(imgs[:, :, 2], ay, ax)
                    # Micro image
                    imgZoom = get(res, :ImageZoom, 1.0)  # >= 1.0
                    rfov, rpix = fov / imgZoom, fov / (imgZoom * (size(imgs, 2)-1))
                    ay = Axis{:y}(0.5*rfov*ratio*mm:-rpix*mm:(0.5*rfov-rpix*(size(imgs,1)-1))*ratio*mm)
                    ax = Axis{:x}(-0.5*rfov*mm:rpix*mm:(-0.5*rfov+rpix*(size(imgs,2)-1))*mm)
                    @assert length(ay)==size(imgs,1) && length(ax) == size(imgs,2)
                    res[Symbol("Image1")] = AxisArray(imgs[:, :, 1], ay, ax)
                else
                    # All images same FOV
                    if nimgs == 1
                        res[Symbol("Image1")] = AxisArray(imgs, ay, ax)
                    else
                        foreach(i -> res[Symbol("Image$i")] = AxisArray(imgs[:, :, i], ay, ax), 1:nimgs)
                    end
                end
            else
                # No scale data
                if nimgs == 1
                    res[Symbol("Image1")] = imgs
                else
                    foreach(i -> res[Symbol("Image$i")] = imgs[:, :, i], 1:nimgs)
                end
            end
        catch err
            @info err
            @info "Unable to read images from $(ios)."
        end
    end
    return res
end

"""
    LinearAlgebra.norm(aa::AxisArray, pt1, pt2, p::Real=2)

Returns the distance between `pt1` and `pt2` using the scaled axes in an AxisArray.
"""
LinearAlgebra.norm(aa::AxisArray, pt1, pt2, p::Real=2) =
    norm( ( AxisArrays.axes(aa,i)[pt1[i]] for i in eachindex(pt1)) .- (AxisArrays.axes(aa,i)[pt2[i]] for i in eachindex(pt2)), p)
LinearAlgebra.norm(aa::AxisArray, pt1::CartesianIndex, pt2::CartesianIndex, p::Real=2) =
    norm( ( AxisArrays.axes(aa,i)[pt1[i]] for i in Base.OneTo(length(pt1))) .- (AxisArrays.axes(aa,i)[pt2[i]] for i in Base.OneTo(length(pt1))), p)

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using Nettle
using Formatting
using Iterators

_regex = r"\w+"
# default in py : r"[\w\u4e00-\u9fcc]+"
# going with the more straightforward get the words

type Simhash
    # this, unfortunately, does not hash on init
    value
    f::Int

    function Simhash(value::BigInt, f::Int)
        new(value, f)
    end
    function Simhash(value::Simhash, f::Int)
        new(value.value, f)
    end
    function Simhash(value::AbstractString, f::Int)
        v = build_by_text(value, f)
        new(v, f)
    end
    function Simhash(value::Dict, f::Int)
        v = build_by_features(value, f)
        new(v, f)
    end
    function Simhash(value::AbstractArray, f::Int)
        v = build_by_features(value, f)
        new(v, f)
    end
    function Simhash(s::Simhash)
        new(s.value, s.f)
    end
end

function _hashfunc(x)
    h = Hasher("md5")
    update!(h, x)
    return parse(BigInt, hexdigest!(h), 16)
end

function distance(one::Simhash, another::Simhash)
    # this is ridiculous, julia
    x =  BigInt((one.value $ another.value) & ((BigInt(1) << BigInt(one.f)) - BigInt(1)))
    ans = 0
    while x > 0
        ans += 1
        x &= x - 1
    end
    return ans
end

function _tokenize(content::String)
    content = lowercase(content)
    content = join(matchall(_regex, content), "")
    return _slide(content)
end

function _slide(content::String, width=4)
    return [content[i:min(i + width-1, length(content))] for i in 1:max(length(content) - width, 1) + 1]
end

function build_by_text(content::String, f)
    features = _tokenize(content)
    fdict = Dict{String, Integer}(g[1] => length(g) for g = groupby(x -> x, sort(features)))
    return build_by_features(fdict, f)
end

function build_by_features(features::Any, f)
    v = zeros(f)
    masks = [BigInt(1) << BigInt(i-1) for i in 1:f]

    if isa(features, Dict)
        # to list of tuples
        features = collect(features)
    end

    for feat in features
        # strings are by default unicode
        if isa(feat, String)
            h = _hashfunc(feat)
            w = 1
        else
            # without assertion for the iterables
            h = _hashfunc(feat[1])
            w = feat[2]
        end

        for i in 1:f
            v[i] += h & masks[i] != 0 ? w : -w
        end
    end

    ans = 0
    for i in 1:f
        if v[i] >= 0
            ans |= masks[i]
        end
    end

    # numeric rep of the string/iterable
    return ans
end


text:
using NumCME
using Test
using Sundials: CVODE_BDF

𝕊 = [[-1, 1, 0] [1, -1, 0] [0, 0, 1] [0, 0, -1]]
x₀ = [1, 0, 0]
k₀₁ = 0.05
k₁₀ = 0.1
λ = 5.0
γ = 0.5
L = 20.0
θ = [k₀₁, k₁₀, λ, γ, L]
a1 = propensity() do x, p
    p[1] * x[1]
end
a2 = propensity() do x, p
    p[2] * x[2]
end
a3 = propensity() do x, p
    p[3] * x[2]
end
a4 = propensity((t, p) -> max(0.0, 1.0 - sin(π * t / p[5]))) do x, p
    p[4] * x[3]
end
propensities = [a1, a2, a3, a4]
sensmodel = CmeModelWithSensitivity(CmeModel(𝕊, propensities, θ))
init_cond = forwardsens_initial_condition([x₀], [1.0], [[0.0] for i in 1:5])
sensfspalgorithm = AdaptiveForwardSensFspSparse(
    ode_method = CVODE_BDF(linear_solver = :GMRES),
    space_adapter = ForwardSensRStepAdapter(10, 10, true)
)
senssol = solve(sensmodel,
    init_cond,
    (0.0, 400.0),
    sensfspalgorithm;
    saveat = [],
    fsptol = 1.0E-6,
    odeatol = 1.0E-10,
    odertol = 1.0E-4)


text:
function gensys_uncertain_altpol(m::AbstractDSGEModel, prob_vec::AbstractVector{S},
                                 altpolicies::Vector{AltPolicy} = [get_setting(m, :alternative_policy)];
                                 apply_altpolicy::Bool = false, regime_switching::Bool = false,
                                 regimes::Union{Int, Vector{Int}, UnitRange{Int}} = 1,
                                 TTT::Matrix{S} = Matrix{S}(undef, 0, 0),
                                 Γ0s::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Γ1s::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Cs::Vector{Vector{S}} = Vector{Vector{S}}(undef, 0),
                                 Ψs::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Πs::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0)) where {S <: Real}

    @assert sum(prob_vec) == 1. "The vector of probabilities must sum to 1"

    @assert length(altpolicies) == length(prob_vec) - 1

    if regime_switching
        error("Regime switching has not been completed for gensys_uncertain_altpol")
        if isempty(Γ0s) || isempty(Γ1s) || isempty(Cs) || isempty(Ψs) || isempty(Πs)
            # If any of these are empty, we recompute the relevant matrices
            for fcast_reg in regimes
                Γ0s[fcast_reg], Γ1s[fcast_reg], Cs[fcast_reg], Ψs[fcast_reg], Πs[fcast_reg] = eqcond(m, fcast_reg)
            end

            # ALSO NEED TO ADD THE APPLY ALTPOLICY OPTION
        end

        # Solve model for creating the "default" matrices in each period
        TTTs_gensys, RRRs_gensys, CCCs_gensys, eus =
            map((Γ0, Γ1, C, Ψ, Π) -> gensys(Γ0, Γ1, C, Ψ, Π, 1 + 1e-6, verbose = verbose), Γ0s, Γ1s, Cs, Ψs, Πs)

        # Check for LAPACK exception, existence and uniqueness
        for eu in eus
            if eu[1] != 1 || eu[2] != 1
                throw(GensysError())
            end
        end

        TTTs_gensys = map(TTT_gensys -> real(TTT_gensys), TTTs_gensys)
        RRRs_gensys = map(RRR_gensys -> real(RRR_gensys), RRRs_gensys)
        CCCs_gensys = map(CCC_gensys -> real(CCC_gensys), CCCs_gensys)

        Γ0_tils, Γ1_tils, Γ2_tils, C_tils, Ψ_tils =
            map((Γ0, Γ1, C, Ψ, Π) -> gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π), Γ0s, Γ1s, Cs, Ψs, Πs)

        # Rows are different altpolicies in each regime, columns are different regimes
        Td = Matrix{Matrix{S}}(undef, length(prob_vec))
        Cd = Matrix{Vector{S}}(undef, length(prob_vec))
        Td[1, :] = TTTs_gensys
        Cd[1, :] = CCCs_gensys
        for fcast_reg in regimes
            for (i, altpolicy) in enumerate(altpolicies)
                j = i + 1 # may refactor to directly use solve(m; apply_altpolicy = true) by updating the alternative policy
                tmpT, _, tmpC = altpolicy.solve(m) # no regime switching for now ; regime_switching = regime_switching, regimes = regimes)
                Td[j, fcast_reg] = tmpT
                Cd[j, fcast_reg] = tmpC
            end
        end

        T̄s = map(reg -> sum([p .* Td_reg[j, reg] for (j, p) in enumerate(prob_vec)]), regimes)
        C̄s = map(reg -> sum([p .* Cd_reg[j, reg] for (j, p) in enumerate(prob_vec)]), regimes)

        Tcal = Vector{Matrix{S}}(undef, length(regimes))
        Rcal = Vector{Matrix{S}}(undef, length(regimes))
        Ccal = Vector{Vector{S}}(undef, length(regimes))

        for i in length(regimes)
            Lmat = (Γ2_tils[i] * T̄s[i] + Γ0_tils[i])
            Tcals[i] = Lmat \ Γ1_tils[i]
            Rcals[i] = Lmat \ Ψ_tils[i]
            Ccals[i] = Lmat \ (C_tils[i] - Γ2_tils[i] * C̄s[i])
        end

        return Tcals, Rcals, Ccals
    else
        Γ0, Γ1, C, Ψ, Π = if apply_altpolicy
            get_setting(m, :alternative_policy).eqcond(m)
        else
            eqcond(m)
        end

        if isempty(TTT)
            assert_cond = haskey(get_settings(m), :uncertain_altpolicy) ? !get_setting(m, :uncertain_altpolicy) : true
            @assert assert_cond "The Setting :uncertain_altpolicy must be false if TTT is empty."
            TTT, _, _ = solve(m; apply_altpolicy = apply_altpolicy)
        end

        Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π)

        inds = 1:n_states(m) # Don't want augmented states
        Td = Vector{Matrix{S}}(undef, length(prob_vec))
        Cd = Vector{Vector{S}}(undef, length(prob_vec))
        Td[1] = TTT[inds, inds]
        Cd[1] = C[inds]
        for (i, altpolicy) in enumerate(altpolicies)
            j = i + 1 # may refactor to directly use solve(m; apply_altpolicy = true) by updating the alternative policy
            tmpT, _, tmpC = altpolicy.solve(m) # no regime switching for now ; regime_switching = regime_switching, regimes = regimes)
            Td[j] = tmpT[inds, inds]
            Cd[j] = tmpC[inds]
        end

        T̄ = sum([p .* Td[i] for (i, p) in enumerate(prob_vec)])
        C̄ = sum([p .* Cd[i] for (i, p) in enumerate(prob_vec)])

        Lmat = (Γ2_til * T̄ + Γ0_til)
        Tcal = Lmat \ Γ1_til
        Rcal = Lmat \ Ψ_til
        Ccal = Lmat \ (C_til - Γ2_til * C̄)

        return Tcal, Rcal, Ccal
    end
end


text:
module RegularizedRegression

include(joinpath("utils", "load_R_library.jl"))
include(joinpath("utils", "install_SLOPE_R_code.jl"))

include("loss_functions.jl")

abstract DataSet
num_indepvars(ds::DataSet) = length(indepvars(ds))
include("matrix_dataset.jl")

# Different Regression algorithms all implement the abstract RegressionLearner interface.
abstract RegressionLearner
abstract RegressionResult
learn(l::RegressionLearner, ds::DataSet) = error("Not implemented!")

include(joinpath("learners", "ncvreg_scad_and_mcp_regression.jl"))
include(joinpath("learners", "slope_regression.jl"))
include(joinpath("learners", "glmnet_lasso_and_ridge_regression.jl"))

export mape, mase, r2, r2adjusted
#export ScadRegression, McpRegression, learn, coefficients

end # RegularizedRegression


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

using Test, Tensors

# mat = Material(Chaboche, tuple())
# props = mat.properties
# props.youngs_modulus = 200.0e3
# props.poissons_ratio = 0.3
# props.yield_stress = 100.0
# props.K_n = 100.0
# props.n_n = 3.0
# props.C_1 = 0.0
# props.D_1 = 100.0
# props.C_2 = 0.0
# props.D_2 = 1000.0
# props.Q = 0.0
# props.b = 0.1
E = 200.0e3
nu = 0.3
parameters = ChabocheParameterState(E = E,
                                    nu = nu,
                                    R0 = 100.0,
                                    Kn = 100.0,
                                    nn = 3.0,
                                    C1 = 0.0,
                                    D1 = 100.0,
                                    C2 = 0.0,
                                    D2 = 1000.0,
                                    Q = 0.0,
                                    b = 0.1)
mat = Chaboche(parameters = parameters)
times = [0.0]
loads = [0.0]
dt = 0.5
G = 0.5*E/(1+nu)
syield = 100.0
# vm = sqrt(3)*G*ga | ea = ga
ea = 2*syield/(sqrt(3)*G)
# Go to elastic border
push!(times, times[end]+dt)
push!(loads, loads[end] + ea*dt)
 # Proceed to plastic flow
push!(times, times[end]+dt)
push!(loads, loads[end] + ea*dt)
 # Reverse direction
push!(times, times[end]+dt)
push!(loads, loads[end] - ea*dt)
 # Continue and pass yield criterion
push!(times, times[end]+dt)
push!(loads, loads[end] - ea*dt)
push!(times, times[end]+dt)
push!(loads, loads[end] - ea*dt)

eeqs = [mat.variables.cumeq]
stresses = [copy(tovoigt(mat.variables.stress))]
for i=2:length(times)
    dtime = times[i]-times[i-1]
    dstrain31 = loads[i]-loads[i-1]
    dstrain = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain31]
    # mat.dtime = dtime
    # mat.dstrain = dstrain
    dstrain_ = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0)
    mat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain_)
    integrate_material!(mat)
    # mat.time += mat.dtime
    # mat.strain .+= mat.dstrain
    # mat.stress .+= mat.dstress
    # mat.properties.plastic_strain .+= mat.properties.dplastic_strain
    # mat.properties.backstress1 .+= mat.properties.dbackstress1
    # mat.properties.backstress2 .+= mat.properties.dbackstress2
    # mat.properties.R += mat.properties.dR
    update_material!(mat)
    push!(stresses, copy(tovoigt(mat.variables.stress)))
    push!(eeqs, mat.variables.cumeq)
    # @info "time = $(mat.time), stress = $(mat.stress), cumeq = $(mat.properties.cumulative_equivalent_plastic_strain))"
end

for i in 1:length(times)
    @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6)
end
s31 = [s[6] for s in stresses]

@test isapprox(s31[2], syield/sqrt(3.0))
@test isapprox(s31[3]*sqrt(3.0), syield + 100.0*((eeqs[3]-eeqs[2])/dt)^(1.0/3.0); rtol=1e-2)
@test isapprox(s31[4], s31[3]-G*ea*dt)
@test isapprox(s31[6]*sqrt(3.0), -(syield + 100.0*((eeqs[6]-eeqs[5])/dt)^(1.0/3.0)); rtol=1e-2)


text:
# # A functional tendency specification layer

# In the balance law (mutating) functions, where we specify fluxes and sources,

# - [`flux_first_order!`](@ref ClimateMachine.BalanceLaws.flux_first_order!)
# - [`flux_second_order!`](@ref ClimateMachine.BalanceLaws.flux_second_order!)
# - and [`source!`](@ref ClimateMachine.BalanceLaws.source!),

# an additional (functional) tendency specification
# layer can be placed on-top that has several nice
# properties. The functional layer:

# - Separates tendency definitions from which tendencies are included in a particular model.
# - Reduces duplicate implementations of tendency definitions (e.g., in optional submodel variants)
# - Allows a more flexible combination of tendencies
# - Allows a simple way to loop over all tendencies for all prognostic variables and recover
#   _each_ flux / source term. This will allow us a simple way to evaluate, for example, the energy budget.

# ## Used modules / imports

# Make running locally easier from ClimateMachine.jl/:
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
    nothing
end

# First, using necessary modules:
using ClimateMachine.BalanceLaws
using ClimateMachine.VariableTemplates
using StaticArrays, Test

# Import methods to overload
import ClimateMachine.BalanceLaws: prognostic_vars, eq_tends, flux

# ## Define a balance law

# Here, we define a simple balance law:

struct MyBalanceLaw <: BalanceLaw end

# ## Define prognostic variable types

# Here, we'll define some prognostic variable types,
# by sub-typing [`AbstractPrognosticVariable`](@ref ClimateMachine.BalanceLaws.AbstractPrognosticVariable),
# for mass and energy:
struct Mass <: AbstractPrognosticVariable end
struct Energy <: AbstractPrognosticVariable end

# Define [`prognostic_vars`](@ref ClimateMachine.BalanceLaws.prognostic_vars),
# which returns _all_ prognostic variables
prognostic_vars(::MyBalanceLaw) = (Mass(), Energy());

# ## Define tendency definition types

# Tendency definitions types are made by subtyping
# [`TendencyDef`](@ref ClimateMachine.BalanceLaws.TendencyDef).
# `TendencyDef` has one type parameters: the
# `AbstractTendencyType`, which can be either
# `Flux{FirstOrder}`, `Flux{SecondOrder}`, or `Source`.
struct Advection <: TendencyDef{Flux{FirstOrder}} end
struct Source1 <: TendencyDef{Source} end
struct Source2 <: TendencyDef{Source} end
struct Diffusion <: TendencyDef{Flux{SecondOrder}} end

# Define [`eq_tends`](@ref ClimateMachine.BalanceLaws.eq_tends),
# which returns a tuple of tendency definitions (those sub-typed
# by [`TendencyDef`](@ref ClimateMachine.BalanceLaws.TendencyDef)),
# given
#  - the prognostic variable
#  - the model (balance law)
#  - the tendency type ([`Flux`](@ref ClimateMachine.BalanceLaws.Flux) or
#    [`Source`](@ref ClimateMachine.BalanceLaws.Source))
eq_tends(::Mass, ::MyBalanceLaw, ::Flux{FirstOrder}) = (Advection(),);
eq_tends(::Energy, ::MyBalanceLaw, ::Flux{FirstOrder}) = (Advection(),);
eq_tends(::Mass, ::MyBalanceLaw, ::Flux{SecondOrder}) = ();
eq_tends(::Energy, ::MyBalanceLaw, ::Flux{SecondOrder}) = (Diffusion(),);
eq_tends(::Mass, ::MyBalanceLaw, ::Source) = (Source1(), Source2());
eq_tends(::Energy, ::MyBalanceLaw, ::Source) = (Source1(), Source2());

# ## Testing `prognostic_vars` `eq_tends`

# To test that `prognostic_vars` and `eq_tends` were
# implemented correctly, we'll create a balance law
# instance and call [`show_tendencies`](@ref ClimateMachine.BalanceLaws.show_tendencies),
# to make sure that the tendency table is accurate.

bl = MyBalanceLaw()
show_tendencies(bl; table_complete = true)

# The table looks correct. Now we're ready to
# add the specification layer.

# ## Adding the tendency specification layer

# For the purpose of this tutorial, we'll only focus
# on adding the layer to the first order flux, since
# doing so for the second order flux and source
# functions follow the same exact pattern. In other words,
# we'll add a layer that tests the `Flux{FirstOrder}` column
# in the table above. First, we'll define individual
# [`flux`](@ref ClimateMachine.BalanceLaws.flux) kernels:
flux(::Mass, ::Advection, bl::MyBalanceLaw, args) =
    args.state.ρ * SVector(1, 1, 1);
flux(::Energy, ::Advection, bl::MyBalanceLaw, args) =
    args.state.ρe * SVector(1, 1, 1);

# !!! note
#     - `flux` should return a 3-componet vector for scalar equations
#     - `flux` should return a 3xN-componet tensor for N-component vector equations
#     - `source` should return a scalar for scalar equations
#     - `source` should return a N-componet vector for N-component vector equations

# Define `flux_first_order!` and utilize `eq_tends`
function flux_first_order!(
    bl::MyBalanceLaw,
    flx::Grad,
    state::Vars,
    aux,
    t,
    direction,
)

    tend_type = Flux{FirstOrder}()
    args = (; state, aux, t, direction)

    ## `Σfluxes(Mass(), eq_tends(Mass(), bl, tend_type), bl, args)` calls
    ## `flux(::Mass, ::Advection, ...)` defined above:
    eqt_ρ = eq_tends(Mass(), bl, tend_type)
    flx.ρ = Σfluxes(Mass(), eqt_ρ, bl, args)

    ## `Σfluxes(Energy(), eq_tends(Energy(), bl, tend_type), bl, args)` calls
    ## `flux(::Energy, ::Advection, ...)` defined above:
    eqt_ρe = eq_tends(Energy(), bl, tend_type)
    flx.ρe = Σfluxes(Energy(), eqt_ρe, bl, args)
    return nothing
end;

# ## Testing the tendency specification layer

# Now, let's test `flux_first_order!` we need to initialize
# some dummy data to call it first:

FT = Float64; # float type
aux = (); # auxiliary fields
t = 0.0; # time
direction = nothing; # Direction

state = Vars{@vars(ρ::FT, ρe::FT)}([1, 2]);
flx = Grad{@vars(ρ::FT, ρe::FT)}(zeros(MArray{Tuple{3, 2}, FT}));

# call `flux_first_order!`
flux_first_order!(bl, flx, state, aux, t, direction);

# Test that `flx` has been properly mutated:
@testset "Test results" begin
    @test flx.ρ == [1, 1, 1]
    @test flx.ρe == [2, 2, 2]
end

nothing


text:
#!/usr/bin/env julia --project

include("../shared_source/boilerplate.jl")
include("TwoDimensionalCompressibleNavierStokesEquations.jl")

ClimateMachine.init()

setups = [
    (;
        name = "rusanov_periodic",
        flux = RusanovNumericalFlux(),
        periodicity = true,
        Nover = 0,
    ),
    (;
        name = "roeflux_periodic",
        flux = RoeNumericalFlux(),
        periodicity = true,
        Nover = 0,
    ),
    (;
        name = "rusanov",
        flux = RusanovNumericalFlux(),
        periodicity = false,
        Nover = 0,
    ),
    (;
        name = "roeflux",
        flux = RoeNumericalFlux(),
        periodicity = false,
        Nover = 0,
    ),
    (;
        name = "rusanov_overintegration",
        flux = RusanovNumericalFlux(),
        periodicity = false,
        Nover = 1,
    ),
    (;
        name = "roeflux_overintegration",
        flux = RoeNumericalFlux(),
        periodicity = false,
        Nover = 1,
    ),
]

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("refvals_bickley_jet.jl")

    ########
    # Define timestepping parameters
    ########
    start_time = 0
    end_time = 200.0
    Δt = 0.02
    method = LSRK54CarpenterKennedy

    timestepper = TimeStepper(method = method, timestep = Δt)

    callbacks = (Info(), StateCheck(10))

    ########
    # Define physical parameters and parameterizations
    ########
    parameters = (
        ϵ = 0.1,  # perturbation size for initial condition
        l = 0.5, # Gaussian width
        k = 0.5, # Sinusoidal wavenumber
        ρₒ = 1.0, # reference density
        c = 2,
        g = 10,
    )

    physics = FluidPhysics(;
        advection = NonLinearAdvectionTerm(),
        dissipation = ConstantViscosity{Float64}(μ = 0, ν = 0, κ = 0),
        coriolis = nothing,
        buoyancy = nothing,
    )

    ########
    # Define boundary conditions
    ########
    ρu_bc = Impenetrable(FreeSlip())
    ρθ_bc = Insulating()
    ρu_bcs = (south = ρu_bc, north = ρu_bc)
    ρθ_bcs = (south = ρθ_bc, north = ρθ_bc)
    BC = (ρθ = ρθ_bcs, ρu = ρu_bcs)

    ########
    # Define initial conditions
    ########

    # The Bickley jet
    U₀(p, x, y, z) = cosh(y)^(-2)

    # Slightly off-center vortical perturbations
    Ψ₀(p, x, y, z) =
        exp(-(y + p.l / 10)^2 / (2 * (p.l^2))) * cos(p.k * x) * cos(p.k * y)

    # Vortical velocity fields (ũ, ṽ) = (-∂ʸ, +∂ˣ) ψ̃
    u₀(p, x, y, z) = Ψ₀(p, x, y, z) * (p.k * tan(p.k * y) + y / (p.l^2))
    v₀(p, x, y, z) = -Ψ₀(p, x, y, z) * p.k * tan(p.k * x)
    θ₀(p, x, y, z) = sin(p.k * y)

    ρ₀(p, x, y, z) = p.ρₒ
    ρu₀(p, x...) = ρ₀(p, x...) * (p.ϵ * u₀(p, x...) + U₀(p, x...))
    ρv₀(p, x...) = ρ₀(p, x...) * p.ϵ * v₀(p, x...)
    ρθ₀(p, x...) = ρ₀(p, x...) * θ₀(p, x...)

    ρu⃗₀(p, x...) = @SVector [ρu₀(p, x...), ρv₀(p, x...)]
    initial_conditions = (ρ = ρ₀, ρu = ρu⃗₀, ρθ = ρθ₀)

    for setup in setups
        @testset "$(setup.name)" begin
            Ωˣ = Periodic(-2π, 2π)
            Ωʸ = IntervalDomain(-2π, 2π, periodic = setup.periodicity)

            grid = DiscretizedDomain(
                Ωˣ × Ωʸ,
                elements = 16,
                polynomial_order = 3,
                overintegration_order = setup.Nover,
            )

            model = SpatialModel(
                balance_law = Fluid2D(),
                physics = physics,
                numerics = (flux = setup.flux,),
                grid = grid,
                boundary_conditions = BC,
                parameters = parameters,
            )

            simulation = Simulation(
                model = model,
                initial_conditions = initial_conditions,
                timestepper = timestepper,
                callbacks = callbacks,
                time = (; start = start_time, finish = end_time),
            )

            ########
            # Run the model
            ########
            evolve!(
                simulation,
                model;
                refDat = getproperty(refVals, Symbol(setup.name)),
            )
        end
    end
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, NetCDF
using YAXArrays.Datasets: open_mfdataset
GLMakie.activate!() # HIDE
let
    dataset = open_mfdataset("./data/ETOPO1_halfdegree.nc")
    lon = dataset.lon.values
    lat = dataset.lat.values
    data = dataset.ETOPO1avg.data[:,:]

    fig = Figure(resolution = (1250,700),fontsize = 25, backgroundcolor = :black)
    ax = Axis3(fig, aspect = (1,0.6,0.1),azimuth = -0.65π, elevation = 0.225π)
    pltobj=surface!(ax,lon,lat, data, colormap = :starrynight, backlight = 1.0f0,
        colorrange =(-6500,5500))
    cbar  =Colorbar(fig,pltobj,label = "ETOPO1 [m]",ticklabelsize = 20,
        tickalign =1, vertical = true, ticksize=15, labelcolor = :white,
        ticklabelcolor = :white, tickcolor = :grey90, height = Relative(0.5))
    hidespines!(ax)
    hidedecorations!(ax)
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "topography.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["GLMakie","NetCDF", "YAXArrays"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# Network Learning
__precompile__(true)

"""
NetworkLearning implements a generic framework for network classification. It could in theory be used for other functionality such as regression and density estimation,
provided that appropriate methods for relational learning (i.e. relational variable generation) and collective inference are added. The framework is designed to make as little assumptions as possible on the elements involved in the process.

# References
[1] S.A. Macskassy, F. Provost "Classification in networked data: A toolkit and a univariate case study", Journal of Machine learning Research 8, 2007, 935-983

[2] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Gallagher, T. Eliassi-Rad "Collective classification in network data", AI Magazine 29(3), 2008
"""
module NetworkLearning

	using LinearAlgebra, SparseArrays
	using Distances
	using LearnBase, MLDataPattern, MLLabelUtils
	using LightGraphs, SimpleWeightedGraphs
	using Random: shuffle!
	using Statistics: mean
	using DelimitedFiles: readdlm
	using Future: copy!

	# Verbosity level, used for debugging:
	# 0 - off, 1 - minimal (i.e. warnings, convergence information), 2 - maximum (i.e. iteration details)
	global const VERBOSE = 0

	# Exports
	export  # Adjacencies
		AbstractAdjacency,
		MatrixAdjacency,
		GraphAdjacency,
		ComputableAdjacency,
		PartialAdjacency,
		EmptyAdjacency,

		# Relational learners
		AbstractRelationalLearner,
		SimpleRN,
		WeightedRN,
		BayesRN,
		ClassDistributionRN,

		# Collective inference
		AbstractCollectiveInferer,
		RelaxationLabelingInferer,
		IterativeClassificationInferer,
		GibbsSamplingInferer,

		# Network learners
	 	AbstractNetworkLearner,
		NetworkLearnerObs,
		NetworkLearnerEnt,

		# Functionality
		@print_verbose,
		fit,
		predict,
		predict!,
		infer!,
		transform,
		transform!,
		adjacency,
		add_adjacency!,
		update_adjacency!,
		strip_adjacency,
		adjacency_matrix,
		adjacency_graph,
		adjacency_obs,
		intdim,
		oppdim,
		matrix_prealloc

	abstract type AbstractNetworkLearner end

	include("utils.jl")									# Small utility functions
	include("adjacency.jl") 								# Adjacency-related structures
	include("rlearners.jl")									# Relational learners
	include("cinference.jl")								# Collective inference algorithms
	include("obslearning.jl")								# Observation-based learning
	include("entlearning.jl")								# Entity-based learning

end


text:
function get_real_uav()
    opts = SimulationOptions()
    uavs = build_uavs(opts)
    return uavs[1]
end

# some functions and constants common to multiple test files
const SHOW_PLOTS = false


text:
abstract type CycloneDetectionAlgorithm end

struct PressureMinima <: CycloneDetectionAlgorithm end

struct VorticityMaxima <: CycloneDetectionAlgorithm end

"""
smooth_and_mask(surf_pres_anomaly, threshold = -9)

Takes a 2d array and smooths it using a rolling median filter.
It then returns the elements of the filtered array whose values are less than the threshold.
"""
function smooth_and_mask(surf_pres_anomaly,threshold=-9, resolution = 2000; target_lengthscale = 30000)
buf1 = similar(surf_pres_anomaly)
buf2 = similar(buf1)
smooth_and_mask!(buf1,buf2,surf_pres_anomaly,threshold,resolution; target_lengthscale)
return buf2
end

"""
smooth_and_mask!(buf,surf_pres_anomaly, threshold = -9)

Takes a 2d array and smooths it using a rolling median filter.
It then returns the elements of the filtered array whose values are less than the threshold.
"""
function smooth_and_mask!(buf,buf2,target_variable,threshold=-9, resolution = 2000; target_lengthscale = 30000)
windowsize = target_lengthscale ÷ (2*resolution) # Length order of magnitude of the eye
windowsize % 2 == 0 ? windowsize = windowsize + 1 : windowsize
mapwindow!(median!,buf,target_variable,[windowsize,windowsize],border="circular");
imfilter!(buf2,buf,Kernel.gaussian(3),"circular");
@inbounds @simd for ind in eachindex(buf2)
    if buf2[ind] >= threshold
         buf2[ind] = -9999
    end
end
return buf2

end
"""
findcyclonecenters_aspressureminima(surf_pres_anomaly,detection_threshold)

Takes a surface pressure anomaly array surf_pres_anomaly[x,y] = surf_pres[x,y] .- mean(surf_pres,dims(1,2))
a detection threshold for the anomaly, and return an array of tuples (x,y)
where each tuple represents the centers of cyclones identified as the minima of the anomaly.
"""
function findcyclonecenters(cyclone_detection_algorithm,target_variable,detection_threshold;grid_spacing=2000, target_lengthscale = 30000)
buf1 = similar(target_variable)
buf2 = similar(buf1)
peaks = findcyclonecenters!(cyclone_detection_algorithm,buf1,buf2,target_variable,detection_threshold;grid_spacing, target_lengthscale)
return peaks
end

"""
findcyclonecenters_aspressureminima!(buf1,buf2,surf_pres_anomaly,detection_threshold)

Takes a surface pressure anomaly array surf_pres_anomaly[x,y] = surf_pres[x,y] .- mean(surf_pres,dims(1,2))
a detection threshold for the anomaly, and return an array of tuples (x,y)
where each tuple represents the centers of cyclones identified as the minima of the anomaly.
"""
function findcyclonecenters!(::VorticityMaxima,buf1,buf2,target_variable,detection_threshold;grid_spacing=2000, target_lengthscale = 30000)
filtered_variable = smooth_and_mask!(buf1,buf2,target_variable,detection_threshold,grid_spacing; target_lengthscale);
peaks              = findlocalmaxima(filtered_variable,[1,2],false);
return peaks
end

"""
findcyclonecenters_aspressureminima!(buf1,buf2,surf_pres_anomaly,detection_threshold)

Takes a surface pressure anomaly array surf_pres_anomaly[x,y] = surf_pres[x,y] .- mean(surf_pres,dims(1,2))
a detection threshold for the anomaly, and return an array of tuples (x,y)
where each tuple represents the centers of cyclones identified as the minima of the anomaly.
"""
function findcyclonecenters!(::PressureMinima,buf1,buf2,target_variable,detection_threshold;grid_spacing=2000, target_lengthscale = 30000)
filtered_variable = smooth_and_mask!(buf1,buf2,target_variable,detection_threshold,grid_spacing; target_lengthscale);
peaks              = findlocalminima(filtered_variable,[1,2],false);
return peaks
end

"""
neighbours_2d(arraysize,indicesCartesian)
Compute de neighboring points of an index in a 2-d array considering periodic boundaries
"""
function neighbours_2d(arraysize,indicesCartesian)
decrement_mod1(num, n) = num == 1 ? n : num - 1
increment_mod1(num, n) = num == n ? one(num) : num + 1
x_ind = indicesCartesian[1]
y_ind = indicesCartesian[2]
size_x = arraysize[1]
size_y = arraysize[2]
@inbounds neighbours = Iterators.product(
    (decrement_mod1(x_ind,size_x), x_ind, increment_mod1(x_ind,size_x)),
    (decrement_mod1(y_ind,size_y), y_ind, increment_mod1(y_ind,size_y)) )
return CartesianIndex.(neighbours)
end

"""
function detect_cyclones!(buf1, buf2,pressure_anomaly,pressure_threshold,resolution)
receive 2 buffers and a pressure anomaly and returns a segmented image with the cyclones

"""
function detect_cyclones(cyclone_detection_algorithm,target_variable,detection_threshold,resolution; target_lengthscale = 30000)
buf1 = similar(target_variable)
buf2 = similar(buf1)
frame_with_detected_cyclones = detect_cyclones!(cyclone_detection_algorithm,buf1,buf2,target_variable,detection_threshold,resolution; target_lengthscale)
return frame_with_detected_cyclones
end

"""
function detect_cyclones!(buf1, buf2,pressure_anomaly,pressure_threshold,resolution)
receive 2 buffers and a pressure anomaly and returns a segmented image with the cyclones

"""
function detect_cyclones!(cyclone_detection_algorithm,buf1,buf2,target_variable,detection_threshold,resolution; target_lengthscale = 30000)
neighbourhood_gen(arraysize) = point -> AvailablePotentialEnergyFramework.neighbours_2d(arraysize,point)
centers = findcyclonecenters!(cyclone_detection_algorithm,buf1,buf2,target_variable,detection_threshold; grid_spacing = resolution, target_lengthscale)#buf2 is the masked array
#@info centers
if length(centers) == 0
    return (nothing,nothing)
end
centers_and_labels = [(centers[i],i) for i in 1:length(centers)]
#mask = pres_anomaly .<= pressure_threshold
push!(centers_and_labels,(findfirst(==(-9999),buf2),1000))
cyclones = seeded_region_growing(buf2,centers_and_labels,neighbourhood_gen(size(target_variable))) #Many allocations? this may be the culprit
return (centers_and_labels,cyclones)
end

"""
    isinteracting(adjacencyMatrix :: SparseMatrixCSC,regionnumber)
Computes the adjacency of identified cyclones and returns true if the cylone-number is adjacent to a region other than the background.
"""
function isinteracting(adjacencyMatrix :: SparseMatrixCSC{T,N},regionnumber) where {T,N}
    @inbounds for i in 1:(adjacencyMatrix.m - 1)
        @views adjacencyMatrix[regionnumber,i] == 1  && (return true)
    end
    return false
end

"""
    isinteracting(cyclones :: SegmentedImage,regionnumber)
Computes the adjacency of identified cyclones and returns true if the cylone-number is adjacent to a region other than the background.
"""
function isinteracting(cyclones :: SegmentedImage,regionnumber)
    adjacency, vert_map = region_adjacency_graph(cyclones,(i,j)->1)
    @inbounds for i in 1:(adjacency.weights.m - 1)
        adjacency.weights[regionnumber,i] == 1 ? (return true) : nothing
    end
    return false
end

"""
add_allcyclones(array :: Array{T,2},radius_bins,array :: Array{T,3},segmentedcyclones,cyclonescenters,gridspacing)

Compute the azimuthal average of some quantity around a center. Repeats the process and averages about all the tropical cyclones detected on the array.
It receives an array with the radius bins to use,the field to average, called `array`, each cyclone as a SegmentedImage,the centers of the cyclones and the gridspacing.
"""
function add_allcyclones(array,segmentedcyclones,cyclonescenters;maskcyclones = true)
addition = zeros(eltype(array),size(array))
buf1 = similar(addition)
buf2 = similar(addition)
cyclonecount = add_allcyclones!(addition,buf1,buf2,array,segmentedcyclones,cyclonescenters;maskcyclones)
return (cyclonecount,addition)
end

"""
add_allcyclones!(addition :: Array{T,2}, buf :: Array{T,2},array :: Array{T,2},radius_bins,array :: Array{T,3},segmentedcyclones,cyclonescenters,gridspacing)

Compute the azimuthal average of some quantity around a center. Repeats the process and averages about all the tropical cyclones detected on the array.
It receives an array with the radius bins to use,the field to average, called `array`, each cyclone as a SegmentedImage,the centers of the cyclones and the gridspacing.
"""
function add_allcyclones!(addition,buf1,buf2,array,segmentedcyclones,cyclonescenters;maskcyclones = true)
if !(size(addition) == size(buf1) == size(array))
    DimensionMismatch("Addition, buffer and array must all have the same dimensions")
end
    center = size(array)[1:2] .÷ 2
    #@info center
adjacency, vert_map = region_adjacency_graph(segmentedcyclones, (i,j)->1)
cyclonecount = 0
if maskcyclones
    @inbounds for cyclone in 1:(length(segmentedcyclones.segment_labels)-1)
        labelsmap = labels_map(segmentedcyclones)
        @inbounds for ind in CartesianIndices(array)
            if labelsmap[ind[1],ind[2]] == cyclone
                buf2[ind] = array[ind]
            else
                buf2[ind] = 0.0
            end
        end
        if !isinteracting(adjacency.weights,cyclone)
            cyclonecount += 1
            addition .+=  shifter!(buf1,buf2,center,cyclonescenters[cyclone][1])
        end
    end
else #no cyclone masking
    @inbounds for cyclone in 1:(length(segmentedcyclones.segment_labels)-1)
        if !isinteracting(adjacency.weights,cyclone)
            cyclonecount += 1
            addition .+=  shifter!(buf1,array,center,cyclonescenters[cyclone][1])
        else
            @debug "Cyclone is too close to another one"
        end
    end
end
return cyclonecount

end

"""
    isindexindistancebin(binlimits,index,center = (0,0),gridspacing=1) = (binlimits[1] < distance(index,center,gridspacing) <= binlimits[2]) ? true : false
computes the distance of one index to the origin and returns true if that distance is inside a bin
"""
isindexindistancebin(binlimits,index,center,gridspacing=1) = (binlimits[1] < distance(index,center,gridspacing) <= binlimits[2]) ? true : false

"""
    averageallindistance(radiusbin,array :: Array{T,2},mask,center,gridspacing = 1)
Create an average of the quantity in array at all the points located between radiusbin[1] and radiusbin[2] from a center.
The points should be masked by a boolean array. It assumes a uniform gridspacing.
"""
function averageallindistance(radiusbin,array :: Array{T,2},center,gridspacing = one(T)) where T
    average = zero(T)
    count = 0
    @inbounds for index in CartesianIndices(array)
        if isindexindistancebin(radiusbin,index,center,gridspacing)
                count += 1
                average += array[index]
        end
    end
    iszero(count) ? average :  average/count
end

"""
    averageallindistance(radiusbin,array :: Array{T,3},mask,center,gridspacing = 1)
Create an average of the quantity in array at all the points located between radiusbin[1] and radiusbin[2] from a center.
The points should be masked by a boolean array. It assumes a uniform gridspacing.
"""
function averageallindistance(radiusbin,array :: Array{T,3},center,gridspacing = one{T} ) where T
    sx,sy,sz = size(array)
    average = zeros(T,sz)
    averageallindistance!(average,radiusbin,array,center,gridspacing)
    return average
end

"""
    averageallindistance!(radiusbin,array :: Array{T,3},center,gridspacing = 1)
Create an average of the quantity in array at all the points located between radiusbin[1] and radiusbin[2] from a center.
The points should be masked by a boolean array. It assumes a uniform gridspacing.
"""
function averageallindistance!(average,radiusbin,array :: Array{T,3},center,gridspacing = 1) where T
    sx,sy,sz = size(array)
    count = 0
    @inbounds for index in CartesianIndices((1:sx,1:sy))
        if isindexindistancebin(radiusbin,index,center,gridspacing)
            #@info "In distance bin"
            count += 1
            for idz in 1:sz
                average[idz] += array[index[1],index[2],idz]
            end
        end
    end
    if !iszero(count)
        average .= average./count
    end
    return average
end

"""
    azimuthalaverage_allcyclones(radius_bins,array :: Array{T,3},segmentedcyclones,cyclonescenters,gridspacing)

Compute the azimuthal average of some quantity around a center. Repeats the process and averages about all the tropical cyclones detected on the array.
It receives an array with the radius bins to use,the field to average, called `array`, each cyclone as a SegmentedImage,the centers of the cyclones and the gridspacing.
"""
function  azimuthalaverage_allcyclones(radius_bins,array :: Array{T,3},segmentedcyclones,cyclonescenters,gridspacing) where T
    G, vert_map = region_adjacency_graph(segmentedcyclones, (i,j)->1)
    labelsmap = labels_map(segmentedcyclones)
    adjacencymatrix = G.weights
    average = zeros(size(array,3),length(radius_bins)-1)
    cyclonecount = 0
    for cyclone in 1:(length(segmentedcyclones.segment_labels)-1)
        if !isinteracting(segmentedcyclones,cyclone)
            cyclonecount += 1
            for bin in 1:(length(radius_bins) - 1)
                average[:,bin] .+= averageallindistance((radius_bins[bin],radius_bins[bin+1]),array,cyclonescenters[cyclone][1],gridspacing)
            end
        end
    end
    if !iszero(cyclonecount)
        return    average./cyclonecount
    else
        return average
    end
end

"""
    azimuthalaverage_allcyclones(radius_bins,array :: Array{T,3},segmentedcyclones,cyclonescenters,gridspacing)

Compute the azimuthal average of some quantity around a center. Repeats the process and averages about all the tropical cyclones detected on the array.
It receives an array with the radius bins to use,the field to average, called `array`, each cyclone as a SegmentedImage,the centers of the cyclones and the gridspacing.
"""
function  azimuthalaverage_allcyclones(radius_bins,array :: Array{T,2},segmentedcyclones,cyclonescenters,gridspacing) where T
    average = zeros(length(radius_bins)-1)
    G, vert_map = region_adjacency_graph(segmentedcyclones, (i,j)->1)
    labelsmap = labels_map(segmentedcyclones)
    adjacencymatrix = G.weights
    cyclonecount = 0
    for cyclone in 1:(length(segmentedcyclones.segment_labels)-1)
        if !isinteracting(segmentedcyclones,cyclone)
            cyclonecount += 1
            for bin in 1:(length(radius_bins) - 1)
                average[bin] += averageallindistance((radius_bins[bin],radius_bins[bin+1]),array,cyclonescenters[cyclone][1],gridspacing)
            end
        end
    end
    if !iszero(cyclonecount)
        return    average./cyclonecount
    else
        return average
    end
end

"""
    shifter(array,domain_center,peak)

Returns an array in which a pressure perturbation center is displaced to the center of the domain using circshift.
Using it may assume periodic domain.
Receives and SAM 3D+time or 2D+time array and two tuples, the (x,y) indices of the domain center,
and the (x,y) indices of the location of the pressure perturbation peak.
"""
function shifter(array,domain_center,peak)
    if ndims(array)==2
      return  circshift(array[:,:],[domain_center[1]-peak[1],domain_center[2]-peak[2]]);
    elseif ndims(array)==3
      return  circshift(array[:,:,:],[domain_center[1]-peak[1],domain_center[2]-peak[2],0]);
    end
end
"""
    shifter!(dest,array,time,domain_center,peak)

Stores in `dest` an array in which a pressure perturbation center is displaced to the center of the domain using circshift.
Using it may assume periodic domain.
Receives and SAM 3D+time or 2D+time array and two tuples, the (x,y) indices of the domain center,
and the (x,y) indices of the location of the pressure perturbation peak.
"""
function shifter!(dest,array,domain_center,peak)
        return  circshift!(dest,array,[domain_center[1]-peak[1],domain_center[2]-peak[2]]);
end

"""
    timemean_nofalseframe(input)

Takes a 3-d or 3-d array and computes the average along the third or fourth dimension, skipping the slices input[:,:,:,i] for which the maxiumum value is 0.0

"""
function timemean_nofalseframe(input::Array{T,4}) where {T<:Real}
        sx,sy,sz,st = size(input)
        copy = zeros(T,sx,sy,sz)
        trueframes = 1
        for t in 1:size(input,4)
            if maximum(view(input,:,:,:,t)) != 0.0
               copy[:,:,:] .+= @views  input[:,:,:,t]
               trueframes = trueframes + 1
            end
        end
        return copy .= copy ./ trueframes
end

function timemean_nofalseframe(input::Array{T,3}) where {T<:Real}
        sx,sy,st = size(input)
        copy = zeros(T,sx,sy)
        trueframes = 1
        for t in 1:size(input,4)
            if maximum(view(input,:,:,t)) != 0.0
               copy[:,:] .+= @views  input[:,:,t]
               trueframes = trueframes + 1
            end
        end
        return copy .= copy ./ trueframes
end

function removefalseframes(input::Array{T,4},peaktimes) where {T<:Real}
        sizet = length(unique(peaktimes))
        copy = zeros(T,size(input,1),size(input,2),size(input,3),sizet)
        trueframes = 1
        Threads.@threads for t in 1:size(input,4)
            if maximum(view(input,:,:,:,t)) != 0
               @views copy[:,:,:,trueframes] .= input[:,:,:,t]
               trueframes = trueframes + 1
            end
        end
        return copy
end

function removefalseframes(input::Array{T,3},peaktimes) where {T<:Real}
        sizet = length(unique(peaktimes))
        copy = zeros(T,size(input,1),size(input,2),sizet)
        trueframes = 1
        Threads.@threads for t in 1:size(input,3)
            if maximum(view(input,:,:,t)) != 0
               @views copy[:,:,trueframes] .= input[:,:,t]
               trueframes = trueframes + 1
            end
        end
        return copy
end


text:
mutable struct NodeType2
    Point :: Array{Float64,1}
    LB :: Float64
    UB :: Float64
end
mutable struct Node
    X :: Array{Float64,1}
    XLB :: Array{Float64}
    XUB :: Array{Float64}
    LB :: Float64
    UB :: Float64
end
mutable struct CutType
    Coeff :: Array{Float64,1}
    RHS :: Float64
end
mutable struct VerticeType
    Point :: Array{Float64,1}
    adj ::  Array{Int64,1}
    J :: Array{Int64,1}
    T :: Float64
end


text:
module CHIP8

export Chip,
       emulate

include("chip.jl")
include("main.jl")

end # module


text:
### Responders are functions which take a req::HTTP.Request and some sort of
### object, and use that object to modify req.response and return it

export json_responder
"""
    json_responder(req::HTTP.Request, o)

Forms a response in a standard format.  `o` can be a `Number`, `String`,
`Dict`, or an `Array` of one of the other types.  The response will look like:

```
{
    "error": false,
    "result": <JSON serialization of o>
}
```
"""
function json_responder(req::HTTP.Request, o::T) where {T<:Union{Number, AbstractString, Array, Dict}}
    b = """{"error": false, "result": $(JSON.json(o))}"""
    req.response.body = bytes(b)
    return req.response
end

export simple_json_responder
"""
    simple_json_responder(req::HTTP.Request, o)

Respond simply with the JSON serialization of an object.  `o` can be a `Number`, `String`,
`Dict`, or an `Array` of one of the other types.
"""
function simple_json_responder(req::HTTP.Request, o::T) where {T<:Union{Number, AbstractString, Array, Dict}}
    b = JSON.json(o)
    req.response.body = bytes(b)
    return req.response
end

export error_responder
"""
    error_responder(req::HTTP.Request, e::String)
    error_responder(req::HTTP.Request, e::Exception)
    error_responder(e::Exception)

Returns a response representing an error in a standard format.  Errors captured
within endpoints are considered to be user errors and return 200 responses.
Errors not explicitly captured by a `try ... catch` block return 500 responses.
The response will look like:

```
{
    "error": true,
    "result": "You need to pass variables named x and y in the query string"
}
```
"""
function error_responder(req::HTTP.Request, e::String)

    # Endpoints can do their own error checking.  Errors caught by them
    # return a 200-code JSON response that has a key 'error' set to true

    # Make sure we don't have unescaped quotes in the error string since
    # we are going to return it as a json
    ee = replace(e, r"(?<!\\)(?:\\{2})*\K\"" => "\\\"")
    b = """{"error": true, "message": "$(ee)"}"""
    req.response.body = bytes(b)
    return req.response
end

# If we think the julia error will be clear to the user then we can just use
# that.
error_responder(req::HTTP.Request, e::Exception) =
    error_responder(req, string(e))

"""
    unhandled_error_responder(req::HTTP.Request,e::Exception)

The default error handler, not intended for use by users.  It includes a gentle reminder that errors should typically be caught elsewhere.
"""
function unhandled_error_responder(req::HTTP.Request, e::Exception)
    @info string("Unhandled error!  You probably want to catch these: ", string(e))
    error_responder(req, string(e))
end

export critical_error_responder
"""
    critical_error_responder(req::HTTP.Request,e::Exception)

Respond with a 500 error response.
"""
function critical_error_responder(req::HTTP.Request, e::Exception)
    @warn string("Critical error!  Returning a 500 response to the client: ", string(e))
    setstatus(req.response, 500)
    return req.response
end


text:
using XLSXReader, RDatasets
using Test, Dates

#Test reference datasets and options
@test readxlsx("datasets/mtcars.xlsx") == dataset("datasets", "mtcars")
@test readxlsx("datasets/mtcars.xlsx") == readxlsx("datasets/mtcars.xlsx", 1)
@test readxlsx("datasets/mtcars.xlsx") == readxlsx("datasets/mtcars.xlsx", "Sheet1")

tdf = DataFrame(x1 = Array(Any[2, missing ,3.2, missing]),
                   x2 = Array(Any[missing, 3.2 ,missing, missing]),
                   x3 = Array(Any[missing, 7.5 ,missing, 4]))

#Test missing values
rdf = readxlsx("datasets/testsets.xlsx", header = false)
@test all(skipmissing(rdf[:x1]) .== skipmissing(tdf[:x1]))
@test all(skipmissing(rdf[:x2]) .== skipmissing(tdf[:x2]))
@test all(skipmissing(rdf[:x3]) .== skipmissing(tdf[:x3]))

#Test dates
ddf = readxlsx("datasets/testsets.xlsx", "dates")
@test names(ddf) == Symbol[:Dates, :Datetimes]
@test ddf[:Dates] == DateTime("2017-01-01T00:00:00") + Dates.Day.(0:8)
@test ddf[:Datetimes] == DateTime("2017-03-01T10:20:15") + Dates.Day.(0:8)

#Test shared string with format
sdf = readxlsx("datasets/testsets.xlsx", "sharedstrings")
@test sdf == readxlsx("datasets/testsets.xlsx", 4)
@test names(sdf) == [:x1col, :x2col, :x3col]
@test Vector(sdf[1,:]) == ["string", "style", "text"]

#Test booleans, empty cols and NA header
bdf = readxlsx("datasets/testsets.xlsx", "booleans")
@test names(bdf) == Symbol[:boolcol, :x1, :moredata]
@test bdf[:boolcol] == [true, false, true, false, false]


text:
try
    cd("C:/Users/francis.smart.ctr/GitDir/SoftDates.jl")
catch
   cd("C:/Users/Francis Smart/Documents/GitHub/SoftDates.jl")
end
using Pkg
Pkg.activate(".")
using Pkg, Dates, DataFrames, Distributions, CSV, StatsBase, IterTools

##############################################################################
# Read in Dates

# Load the dates without errors first
dtrange = CSV.read("sampledata/SampleDateRanges.csv")

DTStart = dtrange[:Start]
DTEnd   = dtrange[:End]
DTtext  = dtrange[:Text]

x = softdate(DTtext[1], Date(DTStart[1]), Date(DTEnd[1]))
[println(y) for y in x[:txt]]

function trunk(x,n=60)
 x = replace(x, "\n\r"=>"]")
 join([x[1:min(length(x), n)], repeat(" ", n-min(length(x), n))])
end
trunk(x[:txt][1])

#presenting problems
# 459
# 463
# 465 # range not Interprettable! What happenned to "17 12 16 Maecenas tellus metus, venenatis ut faucibus in, eleifend at ipsum Nullam enim nibh, rhoncus eu odio nec, condimentum ele"

i = 465
while i <= length(DTtext)
  println("\n------Input $i -- Start Date $(Date(DTStart[i])) till $(Date(DTEnd[i]))------")
  [println(trunk(tx,130)) for tx in split(DTtext[i], r"\n\r|\r\n|\r") if strip(tx) != ""]
  println("------SoftDates-------")
   x = softdate(DTtext[i], dtstart = Date(DTStart[i]), dtend = Date(DTEnd[i]))
  sort!(x, :date)
  [println("$(x[:date][j]) || $(trunk(x[:txt][j], 130)) || $(x[:indt][j]) || $(x[:score][j])") for j in 1:size(x)[1]]
  println()
  print("Type \"break\" or enter to continue:")
  userinput = readline(stdin)
  (length(userinput)>0) && (userinput[1:1] ∈ ["b"]) && break
  global i
  i += 1
end

x = "10 14 16 through 16 12 16"
((x,I...)->replace(x, r"([0-9]{2}).*?([0-9]{2}).*?(1[0-9])(.*)([0-9]{2}).*?([0-9]{2}).*?\3"=>s"\1 \2 20\3 \4 \5 \6 20\3"))(x)


text:
#Simulate data series for inflation
using DSGE, HDF5, StateSpaceRoutines

m = KrusellSmithCT()

#Are these variables defined differently?
n_vars = get_setting(m, :n_vars)
n_expectation_errors = get_setting(m, :n_expectation_errors)
n_shocks = get_setting(m, :n_shocks)
x = zeros(Float64, n_vars)

I = get_setting(m, :I)
n_v = get_setting(m, :n_jump_vars)
n_g = get_setting(m, :n_state_vars)

#How is x set up?
x[1:n_v-1] = reshape(m[:V_ss].value, 2*I, 1)
x[n_v] = m[:inflation_ss].value
x[n_v + 1 : n_v + n_g - 1] = m[:g_ss].value[1:end-1]
x[n_v + n_g + 1] = m[:w_ss].value
x[n_v + n_g + 2] = m[:N_ss].value
x[n_v + n_g + 3] = m[:C_ss].value
x[n_v + n_g + 4] = m[:Y_ss].value
x[n_v + n_g + 5] = m[:B_ss].value

T, R, ~, inverse_basis = solvect(m)

nst = size(T,1) ###
C = zeros(nst)
TTT, RRR, C = transform_transition_matrices(m, T, R, C, true)

n = 40#100
s = zeros(nst,n) #s = zeros(n,406)
s[:,1] = inverse_basis'*x
e = randn(n)

#NOTE: this is dependent on the number of states between you want
dt = 1/2

for i = 2:n
  s[:,i] = T*s[:,i-1]*dt + R*e[i]*(sqrt(dt))
end

data = zeros(1,Int(n/2))
data[1,:] = [s[64,2*i] for i = 1:Int(n/2)]

#Z is 1x130

#write("output.h5", "data", data)

#Create error covariance matrix, matrices in measurement equation

Q = eye(2,2)#
E = zeros(1,1)#
   #steady state output
D = [mean(data)]
Z = zeros(1, nst*2)#(1, nst) # to allow for reduction * inverse_basis #
   #where inflation series shows up in matrix of states
Z[nst*2-1] = 1

s_0 = zeros(nst*2)
   #Exceeds max iterations if we use discrete lyapunov solver
P_0 = 1e6*eye(nst*2) #solve_discrete_lyapunov(T_tilde, R*Q*R')

#@assert false

#Call standard kalman filter
output = kalman_filter(data, TTT, RRR, C, Q, Z, D, E, s_0, P_0)

###
#=

=#


text:
module Weilianqi
export newgame, save, loadgame, expandboard!, checkharvest, center, setcolorset, getgroup, newunit, units, sync!, string2game, loadgame!, placeunit!, Unit, Game, Group
using Gtk, Graphics

dir=joinpath(homedir(),"weilianqi","saves")
if !ispath(dir)
	mkpath(dir)
end

include("qi.jl")

end # module


text:
export runTTC, ttc, solve_ttc, readImage, uniformBlur, rgb2gray

##############################################
# Generic convn function using direct method for computing convolutions:
# Accelerated Convolutions for Efficient Multi-Scale Time to Contact Computation in Julia
# Alexander Amini, Alan Edelman, Berthold Horn
##############################################

function readImage(f)
    img1 = restrict(read(f,Image))
end
function rgb2gray(img)
    img = convert(Array{Float64,3}, separate(img))
    (mean(img,3)[:,:,1])
end
function uniformBlur(img, n)
    I = convn(img, ones(n,n)./(n^2))
    I[n:size(I,1)-n+1,n:size(I,2)-n+1]
end

function krestrict(k,img)
    for i=k:-1:1
        img=restrict(img);
    end
    img
end

function solve_ttc( img1, img2 )
    #SOLVE_TTC Computes the time to contact and focus of expansion
    # using two images as input
    ## Inputs:
    #   img1 = mxn matrix double image
    #   img2 = mxn matrix doulbe image (after img1)
    ## Outputs:
    #   out  = 1x3 vector with the following components:
    #          - T  = time to contact
    #          - x0 = x coordinate of FOE
    #          - y0 = y coordinate of FOE

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

    stack = cat(3,img1,img2);
    (m,n) = size(img1);

    ########################################
    ####  DEFINE Kernel Scaling Factors ####
    DIVIDE_EX = 16/(n); ####################
    DIVIDE_EY = 16/(m); ####################
    DIVIDE_ET = 4/(min(n,m)); ##############
    DIVIDE_GG = 1; #########################
    ########################################

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

    ## setup

    xv = linspace(-n/2,n/2,n-4); yv=linspace(-m/2,m/2,m-4);
    x = [j for i in yv, j in xv];
    y = [i for i in yv, j in xv];

    ## create kernels that average both images
    sobel = [1;2;1]*[-1 0 1]; #sobel operator is 3x3
    sobel3_x = (1/DIVIDE_EX) * repeat(sobel,outer=[1,1,2]);
    sobel3_y = (1/DIVIDE_EY) * repeat(sobel',outer=[1,1,2]);
    prewitt2_t = (1/DIVIDE_ET) * reshape([1;1]*[-1 -1 1 1],2,2,2); #prewitt operator is 2x2

    EX = convn(stack,sobel3_x); #x dir grad
    EY = convn(stack,sobel3_y); #y dir grad
    ET = convn(stack,prewitt2_t); #t dir grad

    EX = EX[3:m-2,3:n-2,2];
    EY = EY[3:m-2,3:n-2,2];
    ET = ET[3:m-2,3:n-2,2];
    GG = EX.*x./(DIVIDE_GG) + EY.*y./(DIVIDE_GG); #radial grad

    ## solve
    E = [EX[:] EY[:] GG[:]];
    v = E\(-ET[:]); #v=3x1 vector

    ## compile results
    x0 = -v[1]/v[3]; #foe x
    y0 = -v[2]/v[3]; #foe y
    T = 1/v[3];      #ttc
    return T,x0,y0,EX,EY,ET,v; #all results

end

function ttc(img1,img2)

    (T_n,x0_n,y0_n,EX_n,EY_n,ET_n,v_n)=solve_ttc(img1,img2)
    (T_p,x0_p,y0_p,EX_p,EY_p,ET_p,v_p) = (T_n,x0_n,y0_n,EX_n,EY_n,ET_n,v_n)

    T_p = T_n+1

    k=0
    while T_p > T_n && T_n>0
        (T_p,x0_p,y0_p,EX_p,EY_p,ET_p,v_p) = (T_n,x0_n,y0_n,EX_n,EY_n,ET_n,v_n)

        img1_ = restrict(uniformBlur(img1,2));
        img2_ = restrict(uniformBlur(img2,2));

        if minimum(size(img1_))<=6
            break
        end

        (T_n,x0_n,y0_n,EX_n,EY_n,ET_n,v_n)=solve_ttc(img1_,img2_)

        img1=img1_; img2=img2_
        k+=1
    end
    #println(k)
    return T_p, x0_p*k, y0_p*k, EX_p, EY_p, ET_p, v_p*k, k

end

runTTC(f) = runTTC(f, "search");
function runTTC(f, method)
    #f = VideoIO.openvideo("./test_vids/Newman4.avi")
    #f = VideoIO.openvideo("./test_vids/Tipper3.avi")
    #f = VideoIO.openvideo("./test_vids/ttc_car_clip.mp4")
    #f = VideoIO.openvideo("./test_vids/toyota.mp4")

    img1_rgb = readImage(f)
    img1 = rgb2gray(img1_rgb)

    if method=="search"
        img1 = uniformBlur(img1,5);
    elseif isa(method,Int)
        img1 = uniformBlur(img1,method)
    end

    (T,x0,y0,EX,EY,ET,v)=solve_ttc(img1,img1)

    (m,n) = size(img1)

    c = canvasgrid(2,2)
    ops = [:pixelspacing => [1,1]]

    (canvas,cEX,cEY,cET) = (c[1,1],c[1,2],c[2,2],c[2,1])

    canvas, ptnr = ImageView.view(canvas, img1_rgb; ops...)
    cEX, ptnrX = ImageView.view(cEX, EX;ops...)
    cEY, ptnrY = ImageView.view(cEY, EY;ops...)
    cET, ptnrT = ImageView.view(cET, ET;ops...)

    FOE = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(0,0, "X"))
    TTC = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(0,0, "X"))

    A = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(0,0, "A"))
    B = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(0,0, "B"))
    C = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(0,0, "C"))

    ImageView.annotate!(cEX, ptnrX, ImageView.AnnotationText(n/2,30, "EX", color=RGB(1,0,0), fontsize=35))
    ImageView.annotate!(cEY, ptnrY, ImageView.AnnotationText(n/2,30, "EY", color=RGB(1,0,0), fontsize=35))
    ImageView.annotate!(cET, ptnrT, ImageView.AnnotationText(n/2,30, "ET", color=RGB(1,0,0), fontsize=35))

    records = Matrix(0,4);
    iter = 0;
    while !eof(f)
        try
            #read and view new image
            img2_rgb = readImage(f)
            img2 = rgb2gray(img2_rgb)

            if method=="search"
                img2 = uniformBlur(img2,5);
                #solve the least squares time to contact
                (T,x0,y0,EX,EY,ET,v,k)=ttc(img1,img2)
            elseif method=="none"
                (T,x0,y0,EX,EY,ET,v)=solve_ttc(img1,img2)
            elseif isa(method,Int)
                img2 = uniformBlur(img2,method)
                (T,x0,y0,EX,EY,ET,v)=solve_ttc(img1,img2)
            end

            iter += 1

            if T>10^5 || T<0
                continue
            end

            records = [records; [iter T x0 y0]];

            _, ptnr = ImageView.view(canvas, img2_rgb;ops...)
            _, ptnrX = ImageView.view(cEX, EX;ops...)
            _, ptnrY = ImageView.view(cEY, EY;ops...)
            _, ptnrT = ImageView.view(cET, ET;ops...)

            #remove old annotations
            ImageView.delete!(canvas, FOE)
            ImageView.delete!(canvas, TTC)
            ImageView.delete!(canvas, A)
            ImageView.delete!(canvas, B)
            ImageView.delete!(canvas, C)

            #add new ones
            FOE = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(x0+(n)/2, y0+(m)/2, "X", color=RGB(0,1,0), fontsize=20))
            TTC = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(n/2,30, string(round(T*10)/10), color=RGB(0,0.5,0.5), fontsize=35))
            A = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(n-60,m-60, string("A:    ",round(v[1]*100)/100), color=RGB(0,0.5,0.5), fontsize=15))
            B = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(n-60,m-40, string("B:    ",round(v[2]*100)/100), color=RGB(0,0.5,0.5), fontsize=15))
            C = ImageView.annotate!(canvas, ptnr, ImageView.AnnotationText(n-60,m-20, string("C: ",round(v[3]*10000)/10000), color=RGB(0,0.5,0.5), fontsize=15))

            #increment moving image window
            img1_rgb = img2_rgb;
            img1 = img2;

            #sleep to slow down
            #sleep(0)

        catch e
            if isa(e, InterruptException) || (isa(e,ErrorException) && e.msg=="invalid window")
                return  records
            else
                throw(e)
            end
        end
    end
    return records
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

const CUDNN_MAJOR = 6
const CUDNN_MINOR = 0
const CUDNN_PATCHLEVEL = 20
const CUDNN_VERSION = CUDNN_MAJOR * 1000 + CUDNN_MINOR * 100 + CUDNN_PATCHLEVEL
const CUDNN_DIM_MAX = 8
const CUDNN_LRN_MIN_N = 1
const CUDNN_LRN_MAX_N = 16
const CUDNN_LRN_MIN_K = 1.0e-5
const CUDNN_LRN_MIN_BETA = 0.01
const CUDNN_BN_MIN_EPSILON = 1.0e-5

typealias cudnnContext Void
typealias cudnnHandle_t Ptr{cudnnContext}

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const CUDNN_STATUS_SUCCESS = (UInt32)(0)
const CUDNN_STATUS_NOT_INITIALIZED = (UInt32)(1)
const CUDNN_STATUS_ALLOC_FAILED = (UInt32)(2)
const CUDNN_STATUS_BAD_PARAM = (UInt32)(3)
const CUDNN_STATUS_INTERNAL_ERROR = (UInt32)(4)
const CUDNN_STATUS_INVALID_VALUE = (UInt32)(5)
const CUDNN_STATUS_ARCH_MISMATCH = (UInt32)(6)
const CUDNN_STATUS_MAPPING_ERROR = (UInt32)(7)
const CUDNN_STATUS_EXECUTION_FAILED = (UInt32)(8)
const CUDNN_STATUS_NOT_SUPPORTED = (UInt32)(9)
const CUDNN_STATUS_LICENSE_ERROR = (UInt32)(10)
const CUDNN_STATUS_RUNTIME_PREREQUISITE_MISSING = (UInt32)(11)
# end enum ANONYMOUS_1

# begin enum cudnnStatus_t
typealias cudnnStatus_t UInt32
const CUDNN_STATUS_SUCCESS = (UInt32)(0)
const CUDNN_STATUS_NOT_INITIALIZED = (UInt32)(1)
const CUDNN_STATUS_ALLOC_FAILED = (UInt32)(2)
const CUDNN_STATUS_BAD_PARAM = (UInt32)(3)
const CUDNN_STATUS_INTERNAL_ERROR = (UInt32)(4)
const CUDNN_STATUS_INVALID_VALUE = (UInt32)(5)
const CUDNN_STATUS_ARCH_MISMATCH = (UInt32)(6)
const CUDNN_STATUS_MAPPING_ERROR = (UInt32)(7)
const CUDNN_STATUS_EXECUTION_FAILED = (UInt32)(8)
const CUDNN_STATUS_NOT_SUPPORTED = (UInt32)(9)
const CUDNN_STATUS_LICENSE_ERROR = (UInt32)(10)
const CUDNN_STATUS_RUNTIME_PREREQUISITE_MISSING = (UInt32)(11)
# end enum cudnnStatus_t

typealias cudnnTensorStruct Void
typealias cudnnTensorDescriptor_t Ptr{cudnnTensorStruct}
typealias cudnnConvolutionStruct Void
typealias cudnnConvolutionDescriptor_t Ptr{cudnnConvolutionStruct}
typealias cudnnPoolingStruct Void
typealias cudnnPoolingDescriptor_t Ptr{cudnnPoolingStruct}
typealias cudnnFilterStruct Void
typealias cudnnFilterDescriptor_t Ptr{cudnnFilterStruct}
typealias cudnnLRNStruct Void
typealias cudnnLRNDescriptor_t Ptr{cudnnLRNStruct}
typealias cudnnActivationStruct Void
typealias cudnnActivationDescriptor_t Ptr{cudnnActivationStruct}
typealias cudnnSpatialTransformerStruct Void
typealias cudnnSpatialTransformerDescriptor_t Ptr{cudnnSpatialTransformerStruct}
typealias cudnnOpTensorStruct Void
typealias cudnnOpTensorDescriptor_t Ptr{cudnnOpTensorStruct}
typealias cudnnReduceTensorStruct Void
typealias cudnnReduceTensorDescriptor_t Ptr{cudnnReduceTensorStruct}

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const CUDNN_DATA_FLOAT = (UInt32)(0)
const CUDNN_DATA_DOUBLE = (UInt32)(1)
const CUDNN_DATA_HALF = (UInt32)(2)
const CUDNN_DATA_INT8 = (UInt32)(3)
const CUDNN_DATA_INT32 = (UInt32)(4)
const CUDNN_DATA_INT8x4 = (UInt32)(5)
# end enum ANONYMOUS_2

# begin enum cudnnDataType_t
typealias cudnnDataType_t UInt32
const CUDNN_DATA_FLOAT = (UInt32)(0)
const CUDNN_DATA_DOUBLE = (UInt32)(1)
const CUDNN_DATA_HALF = (UInt32)(2)
const CUDNN_DATA_INT8 = (UInt32)(3)
const CUDNN_DATA_INT32 = (UInt32)(4)
const CUDNN_DATA_INT8x4 = (UInt32)(5)
# end enum cudnnDataType_t

# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const CUDNN_NOT_PROPAGATE_NAN = (UInt32)(0)
const CUDNN_PROPAGATE_NAN = (UInt32)(1)
# end enum ANONYMOUS_3

# begin enum cudnnNanPropagation_t
typealias cudnnNanPropagation_t UInt32
const CUDNN_NOT_PROPAGATE_NAN = (UInt32)(0)
const CUDNN_PROPAGATE_NAN = (UInt32)(1)
# end enum cudnnNanPropagation_t

# begin enum ANONYMOUS_4
typealias ANONYMOUS_4 UInt32
const CUDNN_NON_DETERMINISTIC = (UInt32)(0)
const CUDNN_DETERMINISTIC = (UInt32)(1)
# end enum ANONYMOUS_4

# begin enum cudnnDeterminism_t
typealias cudnnDeterminism_t UInt32
const CUDNN_NON_DETERMINISTIC = (UInt32)(0)
const CUDNN_DETERMINISTIC = (UInt32)(1)
# end enum cudnnDeterminism_t

# begin enum ANONYMOUS_5
typealias ANONYMOUS_5 UInt32
const CUDNN_TENSOR_NCHW = (UInt32)(0)
const CUDNN_TENSOR_NHWC = (UInt32)(1)
const CUDNN_TENSOR_NCHW_VECT_C = (UInt32)(2)
# end enum ANONYMOUS_5

# begin enum cudnnTensorFormat_t
typealias cudnnTensorFormat_t UInt32
const CUDNN_TENSOR_NCHW = (UInt32)(0)
const CUDNN_TENSOR_NHWC = (UInt32)(1)
const CUDNN_TENSOR_NCHW_VECT_C = (UInt32)(2)
# end enum cudnnTensorFormat_t

# begin enum ANONYMOUS_6
typealias ANONYMOUS_6 UInt32
const CUDNN_OP_TENSOR_ADD = (UInt32)(0)
const CUDNN_OP_TENSOR_MUL = (UInt32)(1)
const CUDNN_OP_TENSOR_MIN = (UInt32)(2)
const CUDNN_OP_TENSOR_MAX = (UInt32)(3)
const CUDNN_OP_TENSOR_SQRT = (UInt32)(4)
# end enum ANONYMOUS_6

# begin enum cudnnOpTensorOp_t
typealias cudnnOpTensorOp_t UInt32
const CUDNN_OP_TENSOR_ADD = (UInt32)(0)
const CUDNN_OP_TENSOR_MUL = (UInt32)(1)
const CUDNN_OP_TENSOR_MIN = (UInt32)(2)
const CUDNN_OP_TENSOR_MAX = (UInt32)(3)
const CUDNN_OP_TENSOR_SQRT = (UInt32)(4)
# end enum cudnnOpTensorOp_t

# begin enum ANONYMOUS_7
typealias ANONYMOUS_7 UInt32
const CUDNN_REDUCE_TENSOR_ADD = (UInt32)(0)
const CUDNN_REDUCE_TENSOR_MUL = (UInt32)(1)
const CUDNN_REDUCE_TENSOR_MIN = (UInt32)(2)
const CUDNN_REDUCE_TENSOR_MAX = (UInt32)(3)
const CUDNN_REDUCE_TENSOR_AMAX = (UInt32)(4)
const CUDNN_REDUCE_TENSOR_AVG = (UInt32)(5)
const CUDNN_REDUCE_TENSOR_NORM1 = (UInt32)(6)
const CUDNN_REDUCE_TENSOR_NORM2 = (UInt32)(7)
# end enum ANONYMOUS_7

# begin enum cudnnReduceTensorOp_t
typealias cudnnReduceTensorOp_t UInt32
const CUDNN_REDUCE_TENSOR_ADD = (UInt32)(0)
const CUDNN_REDUCE_TENSOR_MUL = (UInt32)(1)
const CUDNN_REDUCE_TENSOR_MIN = (UInt32)(2)
const CUDNN_REDUCE_TENSOR_MAX = (UInt32)(3)
const CUDNN_REDUCE_TENSOR_AMAX = (UInt32)(4)
const CUDNN_REDUCE_TENSOR_AVG = (UInt32)(5)
const CUDNN_REDUCE_TENSOR_NORM1 = (UInt32)(6)
const CUDNN_REDUCE_TENSOR_NORM2 = (UInt32)(7)
# end enum cudnnReduceTensorOp_t

# begin enum ANONYMOUS_8
typealias ANONYMOUS_8 UInt32
const CUDNN_REDUCE_TENSOR_NO_INDICES = (UInt32)(0)
const CUDNN_REDUCE_TENSOR_FLATTENED_INDICES = (UInt32)(1)
# end enum ANONYMOUS_8

# begin enum cudnnReduceTensorIndices_t
typealias cudnnReduceTensorIndices_t UInt32
const CUDNN_REDUCE_TENSOR_NO_INDICES = (UInt32)(0)
const CUDNN_REDUCE_TENSOR_FLATTENED_INDICES = (UInt32)(1)
# end enum cudnnReduceTensorIndices_t

# begin enum ANONYMOUS_9
typealias ANONYMOUS_9 UInt32
const CUDNN_32BIT_INDICES = (UInt32)(0)
const CUDNN_64BIT_INDICES = (UInt32)(1)
const CUDNN_16BIT_INDICES = (UInt32)(2)
const CUDNN_8BIT_INDICES = (UInt32)(3)
# end enum ANONYMOUS_9

# begin enum cudnnIndicesType_t
typealias cudnnIndicesType_t UInt32
const CUDNN_32BIT_INDICES = (UInt32)(0)
const CUDNN_64BIT_INDICES = (UInt32)(1)
const CUDNN_16BIT_INDICES = (UInt32)(2)
const CUDNN_8BIT_INDICES = (UInt32)(3)
# end enum cudnnIndicesType_t

# begin enum ANONYMOUS_10
typealias ANONYMOUS_10 UInt32
const CUDNN_CONVOLUTION = (UInt32)(0)
const CUDNN_CROSS_CORRELATION = (UInt32)(1)
# end enum ANONYMOUS_10

# begin enum cudnnConvolutionMode_t
typealias cudnnConvolutionMode_t UInt32
const CUDNN_CONVOLUTION = (UInt32)(0)
const CUDNN_CROSS_CORRELATION = (UInt32)(1)
# end enum cudnnConvolutionMode_t

# begin enum ANONYMOUS_11
typealias ANONYMOUS_11 UInt32
const CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum ANONYMOUS_11

# begin enum cudnnConvolutionFwdPreference_t
typealias cudnnConvolutionFwdPreference_t UInt32
const CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionFwdPreference_t

# begin enum ANONYMOUS_12
typealias ANONYMOUS_12 UInt32
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_ALGO_GEMM = (UInt32)(2)
const CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = (UInt32)(3)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT = (UInt32)(4)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING = (UInt32)(5)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD = (UInt32)(6)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED = (UInt32)(7)
const CUDNN_CONVOLUTION_FWD_ALGO_COUNT = (UInt32)(8)
# end enum ANONYMOUS_12

# begin enum cudnnConvolutionFwdAlgo_t
typealias cudnnConvolutionFwdAlgo_t UInt32
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_ALGO_GEMM = (UInt32)(2)
const CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = (UInt32)(3)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT = (UInt32)(4)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING = (UInt32)(5)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD = (UInt32)(6)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED = (UInt32)(7)
const CUDNN_CONVOLUTION_FWD_ALGO_COUNT = (UInt32)(8)
# end enum cudnnConvolutionFwdAlgo_t

immutable Array_4_Cint
    d1::Cint
    d2::Cint
    d3::Cint
    d4::Cint
end

zero(::Type{Array_4_Cint}) = begin  # /home/shindo/local-lemon/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
        Array_4_Cint(fill(zero(Cint),4)...)
    end

immutable cudnnConvolutionFwdAlgoPerf_t
    algo::cudnnConvolutionFwdAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    reserved::Array_4_Cint
end

# begin enum ANONYMOUS_13
typealias ANONYMOUS_13 UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum ANONYMOUS_13

# begin enum cudnnConvolutionBwdFilterPreference_t
typealias cudnnConvolutionBwdFilterPreference_t UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionBwdFilterPreference_t

# begin enum ANONYMOUS_14
typealias ANONYMOUS_14 UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3 = (UInt32)(3)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD = (UInt32)(4)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED = (UInt32)(5)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING = (UInt32)(6)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT = (UInt32)(7)
# end enum ANONYMOUS_14

# begin enum cudnnConvolutionBwdFilterAlgo_t
typealias cudnnConvolutionBwdFilterAlgo_t UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3 = (UInt32)(3)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD = (UInt32)(4)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED = (UInt32)(5)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING = (UInt32)(6)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT = (UInt32)(7)
# end enum cudnnConvolutionBwdFilterAlgo_t

immutable cudnnConvolutionBwdFilterAlgoPerf_t
    algo::cudnnConvolutionBwdFilterAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    reserved::Array_4_Cint
end

# begin enum ANONYMOUS_15
typealias ANONYMOUS_15 UInt32
const CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum ANONYMOUS_15

# begin enum cudnnConvolutionBwdDataPreference_t
typealias cudnnConvolutionBwdDataPreference_t UInt32
const CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionBwdDataPreference_t

# begin enum ANONYMOUS_16
typealias ANONYMOUS_16 UInt32
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING = (UInt32)(3)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD = (UInt32)(4)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED = (UInt32)(5)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT = (UInt32)(6)
# end enum ANONYMOUS_16

# begin enum cudnnConvolutionBwdDataAlgo_t
typealias cudnnConvolutionBwdDataAlgo_t UInt32
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING = (UInt32)(3)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD = (UInt32)(4)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED = (UInt32)(5)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT = (UInt32)(6)
# end enum cudnnConvolutionBwdDataAlgo_t

immutable cudnnConvolutionBwdDataAlgoPerf_t
    algo::cudnnConvolutionBwdDataAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    reserved::Array_4_Cint
end

# begin enum ANONYMOUS_17
typealias ANONYMOUS_17 UInt32
const CUDNN_SOFTMAX_FAST = (UInt32)(0)
const CUDNN_SOFTMAX_ACCURATE = (UInt32)(1)
const CUDNN_SOFTMAX_LOG = (UInt32)(2)
# end enum ANONYMOUS_17

# begin enum cudnnSoftmaxAlgorithm_t
typealias cudnnSoftmaxAlgorithm_t UInt32
const CUDNN_SOFTMAX_FAST = (UInt32)(0)
const CUDNN_SOFTMAX_ACCURATE = (UInt32)(1)
const CUDNN_SOFTMAX_LOG = (UInt32)(2)
# end enum cudnnSoftmaxAlgorithm_t

# begin enum ANONYMOUS_18
typealias ANONYMOUS_18 UInt32
const CUDNN_SOFTMAX_MODE_INSTANCE = (UInt32)(0)
const CUDNN_SOFTMAX_MODE_CHANNEL = (UInt32)(1)
# end enum ANONYMOUS_18

# begin enum cudnnSoftmaxMode_t
typealias cudnnSoftmaxMode_t UInt32
const CUDNN_SOFTMAX_MODE_INSTANCE = (UInt32)(0)
const CUDNN_SOFTMAX_MODE_CHANNEL = (UInt32)(1)
# end enum cudnnSoftmaxMode_t

# begin enum ANONYMOUS_19
typealias ANONYMOUS_19 UInt32
const CUDNN_POOLING_MAX = (UInt32)(0)
const CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = (UInt32)(1)
const CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = (UInt32)(2)
const CUDNN_POOLING_MAX_DETERMINISTIC = (UInt32)(3)
# end enum ANONYMOUS_19

# begin enum cudnnPoolingMode_t
typealias cudnnPoolingMode_t UInt32
const CUDNN_POOLING_MAX = (UInt32)(0)
const CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = (UInt32)(1)
const CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = (UInt32)(2)
const CUDNN_POOLING_MAX_DETERMINISTIC = (UInt32)(3)
# end enum cudnnPoolingMode_t

# begin enum ANONYMOUS_20
typealias ANONYMOUS_20 UInt32
const CUDNN_ACTIVATION_SIGMOID = (UInt32)(0)
const CUDNN_ACTIVATION_RELU = (UInt32)(1)
const CUDNN_ACTIVATION_TANH = (UInt32)(2)
const CUDNN_ACTIVATION_CLIPPED_RELU = (UInt32)(3)
const CUDNN_ACTIVATION_ELU = (UInt32)(4)
# end enum ANONYMOUS_20

# begin enum cudnnActivationMode_t
typealias cudnnActivationMode_t UInt32
const CUDNN_ACTIVATION_SIGMOID = (UInt32)(0)
const CUDNN_ACTIVATION_RELU = (UInt32)(1)
const CUDNN_ACTIVATION_TANH = (UInt32)(2)
const CUDNN_ACTIVATION_CLIPPED_RELU = (UInt32)(3)
const CUDNN_ACTIVATION_ELU = (UInt32)(4)
# end enum cudnnActivationMode_t

# begin enum ANONYMOUS_21
typealias ANONYMOUS_21 UInt32
const CUDNN_LRN_CROSS_CHANNEL_DIM1 = (UInt32)(0)
# end enum ANONYMOUS_21

# begin enum cudnnLRNMode_t
typealias cudnnLRNMode_t UInt32
const CUDNN_LRN_CROSS_CHANNEL_DIM1 = (UInt32)(0)
# end enum cudnnLRNMode_t

# begin enum ANONYMOUS_22
typealias ANONYMOUS_22 UInt32
const CUDNN_DIVNORM_PRECOMPUTED_MEANS = (UInt32)(0)
# end enum ANONYMOUS_22

# begin enum cudnnDivNormMode_t
typealias cudnnDivNormMode_t UInt32
const CUDNN_DIVNORM_PRECOMPUTED_MEANS = (UInt32)(0)
# end enum cudnnDivNormMode_t

# begin enum ANONYMOUS_23
typealias ANONYMOUS_23 UInt32
const CUDNN_BATCHNORM_PER_ACTIVATION = (UInt32)(0)
const CUDNN_BATCHNORM_SPATIAL = (UInt32)(1)
# end enum ANONYMOUS_23

# begin enum cudnnBatchNormMode_t
typealias cudnnBatchNormMode_t UInt32
const CUDNN_BATCHNORM_PER_ACTIVATION = (UInt32)(0)
const CUDNN_BATCHNORM_SPATIAL = (UInt32)(1)
# end enum cudnnBatchNormMode_t

# begin enum ANONYMOUS_24
typealias ANONYMOUS_24 UInt32
const CUDNN_SAMPLER_BILINEAR = (UInt32)(0)
# end enum ANONYMOUS_24

# begin enum cudnnSamplerType_t
typealias cudnnSamplerType_t UInt32
const CUDNN_SAMPLER_BILINEAR = (UInt32)(0)
# end enum cudnnSamplerType_t

typealias cudnnDropoutStruct Void
typealias cudnnDropoutDescriptor_t Ptr{cudnnDropoutStruct}

# begin enum ANONYMOUS_25
typealias ANONYMOUS_25 UInt32
const CUDNN_RNN_RELU = (UInt32)(0)
const CUDNN_RNN_TANH = (UInt32)(1)
const CUDNN_LSTM = (UInt32)(2)
const CUDNN_GRU = (UInt32)(3)
# end enum ANONYMOUS_25

# begin enum cudnnRNNMode_t
typealias cudnnRNNMode_t UInt32
const CUDNN_RNN_RELU = (UInt32)(0)
const CUDNN_RNN_TANH = (UInt32)(1)
const CUDNN_LSTM = (UInt32)(2)
const CUDNN_GRU = (UInt32)(3)
# end enum cudnnRNNMode_t

# begin enum ANONYMOUS_26
typealias ANONYMOUS_26 UInt32
const CUDNN_UNIDIRECTIONAL = (UInt32)(0)
const CUDNN_BIDIRECTIONAL = (UInt32)(1)
# end enum ANONYMOUS_26

# begin enum cudnnDirectionMode_t
typealias cudnnDirectionMode_t UInt32
const CUDNN_UNIDIRECTIONAL = (UInt32)(0)
const CUDNN_BIDIRECTIONAL = (UInt32)(1)
# end enum cudnnDirectionMode_t

# begin enum ANONYMOUS_27
typealias ANONYMOUS_27 UInt32
const CUDNN_LINEAR_INPUT = (UInt32)(0)
const CUDNN_SKIP_INPUT = (UInt32)(1)
# end enum ANONYMOUS_27

# begin enum cudnnRNNInputMode_t
typealias cudnnRNNInputMode_t UInt32
const CUDNN_LINEAR_INPUT = (UInt32)(0)
const CUDNN_SKIP_INPUT = (UInt32)(1)
# end enum cudnnRNNInputMode_t

# begin enum ANONYMOUS_28
typealias ANONYMOUS_28 UInt32
const CUDNN_RNN_ALGO_STANDARD = (UInt32)(0)
const CUDNN_RNN_ALGO_PERSIST_STATIC = (UInt32)(1)
const CUDNN_RNN_ALGO_PERSIST_DYNAMIC = (UInt32)(2)
# end enum ANONYMOUS_28

# begin enum cudnnRNNAlgo_t
typealias cudnnRNNAlgo_t UInt32
const CUDNN_RNN_ALGO_STANDARD = (UInt32)(0)
const CUDNN_RNN_ALGO_PERSIST_STATIC = (UInt32)(1)
const CUDNN_RNN_ALGO_PERSIST_DYNAMIC = (UInt32)(2)
# end enum cudnnRNNAlgo_t

typealias cudnnRNNStruct Void
typealias cudnnRNNDescriptor_t Ptr{cudnnRNNStruct}
typealias cudnnPersistentRNNPlan Void
typealias cudnnPersistentRNNPlan_t Ptr{cudnnPersistentRNNPlan}