text:
module LaleOps

using PyCall

# standard included modules
using DataFrames
using Random
using ..AbsTypes
using ..LaleAbsTypes
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!, lalepropertynames
import ..LaleAbsTypes: predict, pretty_print
import ..LaleLibOps: LalePipe

export fit!, transform!, fit, transform, predict
export LaleOp, skops, autogenops, lalelibops, pretty_print, lalesummary

const sk_dict = Dict{String,PyObject}()
const ag_dict = Dict{String,PyObject}()
const ll_dict = Dict{String,PyObject}()
const mb_dict = Dict{String,PyObject}()
const SK   = PyNULL()
const AG   = PyNULL()
const LL   = PyNULL()
const MB   = PyNULL()

include("sklearn-const.jl") # SKVec
include("autogen-const.jl") # AGVec
include("liblale-const.jl") # LLVec
include("imblearn-const.jl") # MBVec

function __init__()
   copy!(SK, pyimport("lale.lib.sklearn"))
   copy!(AG, pyimport("lale.lib.autogen"))
   copy!(LL, pyimport("lale.lib.lale"))
   copy!(MB, pyimport("lale.lib.imblearn"))

   for lr in SKVec
      sk_dict[lr] = SK
   end
   for lr in AGVec
      ag_dict[lr] = AG
   end
   for lr in LLVec
      ll_dict[lr] = LL
   end
   for lr in MBVec
      mb_dict[lr] = MB
   end
end

"""
    LaleOp(learner::String, args::Dict=Dict())

A Scikitlearn wrapper to load the different machine learning models.
Invoking `skops()` or `autogenops()` will list the available learners. Please
consult Scikitlearn documentation for arguments to pass.

Implements `fit!` and `transform!`.
"""
mutable struct LaleOp <: LaleOperator
   name::String
   model::Dict{Symbol,Any}

   function LaleOp(args=Dict{Symbol,Any}())
      default_args=Dict{Symbol,Any}(
         :name => "laleop",
         :type => "sklearn", # or "autogen" or "lale"
         :output => :class,
         :learner => "LinearSVC",
         :impl_args => Dict{Symbol,Any}()
      )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      lr = cargs[:learner]
      type = cargs[:type]
      if type == "sklearn"
         if !(lr in keys(sk_dict))
            println("$lr is not supported.")
            println()
            skops()
            throw(ArgumentError("Argument keyword error"))
         end
      elseif type == "autogen"
         if !(lr in keys(ag_dict))
            println("$lr is not supported.")
            println()
            autogenops()
            throw(ArgumentError("Argument keyword error"))
         end
      elseif type == "lale"
         if !(lr in keys(ll_dict))
            println("$lr is not supported.")
            println()
            lalelibops()
            throw(ArgumentError("Argument keyword error"))
         end
      elseif type == "imblearn"
         if !(lr in keys(mb_dict))
            println("$lr is not supported.")
            println()
            lalelibops()
            throw(ArgumentError("Argument keyword error"))
         end
      else
         throw(ArgumentError("Argument keyword error"))
      end

      impl_args = cargs[:impl_args]
      learner = cargs[:learner]
      if type == "sklearn"
         py_learner = getproperty(sk_dict[learner],learner)
      elseif type == "autogen"
         py_learner = getproperty(ag_dict[learner],learner)
      elseif type == "lale"
         py_learner = getproperty(ll_dict[learner],learner)
      elseif type == "imblearn"
         py_learner = getproperty(mb_dict[learner],learner)
      else
         throw(ArgumentError("$learner not found"))
      end
      cargs[:laleobj]   = py_learner
      new(cargs[:name],cargs)
   end
end

# initialized laleobject and pass the args for fit
function (x::LaleOp)(;args...)
   pyobj = x.model[:laleobj]
   pr = pyobj(;args...)
   x.model[:laleobj] = pr
   x.model[:impl_args]=Dict(pairs(args))
   return x
end

# check if operator is Hyperopt or SMOTE,etc
function (x::LaleOp)(pipe::LalePipe; args...)
   pyobj = x.model[:laleobj]
   lpipe = pipe.model[:laleobj]
   if x.model[:learner] == "Hyperopt" # use estimator arg
      pr = pyobj(estimator=lpipe;args...)
      x.model[:laleobj] = pr
      x.model[:impl_args]=Dict(pairs(args))
      return x
   else #imbalance operator, use operator arg
      pr = pyobj(operator=lpipe;args...) |> LalePipe
      return pr
   end
end

function LaleOp(learner::String, args::Dict)
   LaleOp(Dict(:learner => learner,:name=>learner, args...))
end

function LaleOp(learner::String,type::String="sklearn"; args...)
   LaleOp(Dict(:learner => learner,:name=>learner,:type=>type,:impl_args=>Dict(pairs(args))))
end

"""
    function skops()

List the available scikitlearn machine learners.
"""
function skops()
  learners = keys(sk_dict) |> collect |> x-> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: LaleOp(name::String, type=\"sklearn\"; args...)")
  println("where 'name can be one of:")
  println()
  [print(learner," ") for learner in learners]
  println()
  println()
  println("and 'args' are the corresponding learner's initial parameters.")
  println("Note: Consult Scikitlearn's online help for more details about the learner's arguments.")
end

"""
    function autogenops()

List the available scikitlearn machine learners.
"""
function autogenops()
  learners = keys(ag_dict) |> collect |> x-> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: LaleOp(name::String, type=\"autogen\"; args...)")
  println("where 'name can be one of:")
  println()
  [print(learner," ") for learner in learners]
  println()
  println()
  println("and 'args' are the corresponding learner's initial parameters.")
  println("Note: Consult Scikitlearn's online help for more details about the learner's arguments.")
end

"""
    function lalelibops()

List the available scikitlearn machine learners.
"""
function lalelibops()
  learners = keys(ll_dict) |> collect |> x-> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: LaleOp(name::String, type=\"lale\"; args...)")
  println("where 'name can be one of:")
  println()
  [print(learner," ") for learner in learners]
  println()
  println()
  println("and 'args' are the corresponding learner's initial parameters.")
  println("Note: Consult Scikitlearn's online help for more details about the learner's arguments.")
end

function fit!(lale::LaleOp, xx::DataFrame, y::Vector=Vector())::Nothing
  x = xx |> Array
  impl_args = copy(lale.model[:impl_args])
  learner = lale.model[:learner]
  type = lale.model[:type]
  py_learner=PyObject
  if type == "sklearn"
     py_learner = getproperty(sk_dict[learner],learner)
  elseif type == "autogen"
     py_learner = getproperty(ag_dict[learner],learner)
  elseif type == "lale"
     py_learner = getproperty(ll_dict[learner],learner)
  else
     throw(ArgumentError("$learner not found"))
  end

  # Train
  modelobj = py_learner(;impl_args...)
  trained = PyNULL
  if :predict ∈ lalepropertynames(modelobj)
     # learner
     trained=modelobj.fit(x,y)
  else # transformer
     trained=modelobj.fit(x)
  end
  lale.model[:laleobj]   = trained
  lale.model[:impl_args] = impl_args
  return nothing
end

function fit(lale::LaleOp, xx::DataFrame, y::Vector=Vector())::LaleOp
   fit!(lale,xx,y)
   lcopy =deepcopy(lale)
   return lcopy
end

function transform!(lale::LaleOp, xx::DataFrame)::Union{DataFrame, Vector}
   x = deepcopy(xx)|> Array
   laleobj = lale.model[:laleobj]
   # transform is predict for learners
   if :predict ∈ lalepropertynames(laleobj)
      return collect(laleobj.predict(x))
   elseif :transform ∈ lalepropertynames(laleobj)
      return collect(laleobj.transform(x)) |> x -> DataFrame(x,:auto)
   else
      throw(KeyError("predict/transform function not available"))
   end
end

function transform(lale::LaleOp, xx::DataFrame)::Union{DataFrame, Vector}
   return transform!(lale,xx)
end

predict(lale::LaleOp, xx::DataFrame)  = transform!(lale,xx)

function pretty_print(lop::LaleOp;args...)
~    lop.model[:laleobj].pretty_print(;args...)
end

function lalesummary(op::LaleOp)
   pyobject = op.model[:laleobj]
   if :summary ∈ lalepropertynames(pyobject)
      pyobject.summary()
   else
      nothing
   end
end

end


text:
#####
##### Blasius Smagorinsky: the version of the Smagorinsky turbulence closure
##### used in the UK Met Office code "Blasius" (see Polton and Belcher 2007).
#####

"""
    BlasiusSmagorinsky{ML, FT}

Parameters for the version of the Smagorinsky closure used in the UK Met Office code
Blasius, according to Polton and Belcher (2007).
"""
struct BlasiusSmagorinsky{ML, FT, P, K} <: AbstractSmagorinsky{FT}
               Pr :: P
                ν :: FT
                κ :: K
    mixing_length :: ML

    function BlasiusSmagorinsky{FT}(Pr, ν, κ, mixing_length) where FT
        Pr = convert_diffusivity(FT, Pr)
         κ = convert_diffusivity(FT, κ)
        return new{typeof(mixing_length), FT, typeof(Pr), typeof(κ)}(Pr, ν, κ, mixing_length)
    end
end

"""
    BlasiusSmagorinsky(FT=Float64; Pr=1.0, ν=1.05e-6, κ=1.46e-7)

Returns a `BlasiusSmagorinsky` closure object of type `FT`.

Keyword arguments
=================
    - `Pr` : Turbulent Prandtl numbers for each tracer. Either a constant applied to every
             tracer, or a `NamedTuple` with fields for each tracer individually.
    - `ν`  : Constant background viscosity for momentum
    - `κ`  : Constant background diffusivity for tracer. Can either be a single number
             applied to all tracers, or `NamedTuple` of diffusivities corresponding to each
             tracer.

References
==========
Polton, J. A., and Belcher, S. E. (2007), "Langmuir turbulence and deeply penetrating jets
    in an unstratified mixed layer." Journal of Geophysical Research: Oceans.
"""
BlasiusSmagorinsky(FT=Float64; Pr=1.0, ν=ν₀, κ=κ₀, mixing_length=nothing) =
    BlasiusSmagorinsky{FT}(Pr, ν, κ, mixing_length)

function with_tracers(tracers, closure::BlasiusSmagorinsky{ML, FT}) where {ML, FT}
    Pr = tracer_diffusivities(tracers, closure.Pr)
     κ = tracer_diffusivities(tracers, closure.κ)
    return BlasiusSmagorinsky{FT}(Pr, closure.ν, κ, closure.mixing_length)
end

const BS = BlasiusSmagorinsky

@inline Lm²(i, j, k, grid, closure::BS{<:AbstractFloat}, args...) =
    closure.mixing_length^2

@inline Lm²(i, j, k, grid, closure::BS{<:Nothing}, args...) =
    geo_mean_Δᶠ(i, j, k, grid)^2

#####
##### Mixing length via Monin-Obukhov theory
#####

"""
    MoninObukhovMixingLength(FT=Float64; kwargs...)

Returns a `MoninObukhovMixingLength closure object of type `FT` with

    * Cκ : Von Karman constant
    * z₀ : roughness length
    * L₀ : 'base' mixing length
    * Qu : surface velocity flux (Qu = τ/ρ₀)
    * Qb : surface buoyancy flux

The surface velocity flux and buoyancy flux are restricted to constants for now.
"""
struct MoninObukhovMixingLength{L, FT, U, B}
    Cκ :: FT
    z₀ :: FT
    L₀ :: L
    Qu :: U
    Qb :: U
end

MoninObukhovMixingLength(FT=Float64; Qu, Qb, Cκ=0.4, z₀=0.0, L₀=nothing) =
    MoninObukhovMixingLength(FT(Cκ), FT(z₀), L₀, Qu, Qb)

const MO = MoninObukhovMixingLength

@inline ϕm(z, Lm) = ϕm(z, Lm.Qu, Lm.Qb, Lm.Cκ)

"""
    ϕm(z, Qu, Qb, Cκ)

Return the stability function for monentum at height
`z` in terms of the velocity flux `Qu`, buoyancy flux
`Qb`, and von Karman constant `Cκ`.
"""
@inline function ϕm(z::FT, Qu, Qb, Cκ) where FT
    if Qu == 0
        return zero(FT)
    elseif Qb > 0 # unstable
        return ( 1 - 15Cκ * Qb * z / Qu^FT(1.5) )^FT(-0.25)
    else # neutral or stable
        return 1 + 5Cκ * Qb * z / Qu^FT(1.5)
    end
end

@inline function Lm²(i, j, k, grid, clo::BlasiusSmagorinsky{<:MO, FT}, Σ², N²) where FT
    Lm = clo.mixing_length
    z = @inbounds grid.zC[i, j, k]
    Lm⁻² = (
             ( ϕm(z + Lm.z₀, Lm) * buoyancy_factor(Σ², N²)^FT(0.25) / (Lm.Cκ * (z + Lm.z₀)) )^2
            + 1 / L₀(i, j, k, grid, Lm)^2
           )
    return 1 / Lm⁻²
end

@inline L₀(i, j, k, grid, mixing_length::MO{<:Number})  = mixing_length.L₀
@inline L₀(i, j, k, grid, mixing_length::MO{<:Nothing}) = geo_mean_Δᶠ(i, j, k, grid)
@inline buoyancy_factor(Σ², N²::FT) where FT = ifelse(Σ²==0, zero(FT), max(zero(FT), one(FT) - N² / Σ²))

"""
    νₑ_blasius(Lm², Σ², N²)

Return the eddy viscosity for the BLASIUS version of constant Smagorinsky
(Polton and Belcher, 2007) given the squared mixing length scale `Lm²`,
strain tensor dot product `Σ²`, and buoyancy gradient / squared buoyancy
frequency `N²`.
"""
@inline νₑ_blasius(Lm², Σ², N²::FT) where FT = Lm² * sqrt(Σ²/2 * buoyancy_factor(Σ², N²))

@inline function νᶜᶜᶜ(i, j, k, grid, clo::BlasiusSmagorinsky{ML, FT}, buoyancy, U, C) where {ML, FT}
    Σ² = ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w)
    N² = max(zero(FT), ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, C))
    Lm_sq = Lm²(i, j, k, grid, clo, Σ², N²)

    return νₑ_blasius(Lm_sq, Σ², N²) + clo.ν
end

ext: jl


text:
#=##############################################################################
# DESCRIPTION
    Post-processing acoustic methods.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez and Tyler Critchfield
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Nov 2019
  * License   : MIT
=###############################################################################

"Adds two sound pressure levels together."
function addSPL(x, y)

    #* For better understanding of what's going on, the following commented
    #* steps are equivalent to the returned line
    # pref = 2e-5

    # rms1 = pref^2 * 10 ^ (oaspl1 / 10) # convert to root mean square pressure
    # rms2 = pref^2 * 10 ^ (oaspl2 / 10)
    # rms_combined = rms1 + rms2 # add together

    # oaspl_combined = 10 * log10(rms_combined / (pref^2)) # convert back to decibels

    # return oaspl_combined

    return 10 * log10.(10 .^ (x/10) .+ 10 .^ (y/10))
end

"Adds a vector of sound pressure levels."
function addSPL(x::AbstractArray)

    return 10 * log10.(sum(10 .^ (x/10)))
end

"SPL spectrum to OASPL - works for A-weighted and non-A-weighted
    spectrums"
function SPL2OASPL(spls)                # vector of sound pressure levels over the spectrum

    # don't need a reference pressure when converting back and forth
    ps = 10 .^ (spls / 20)
    ps_sum = sum(ps.^2)

    oaspl = 10 * log10(ps_sum)

    return oaspl
end

"Takes a sound pressure level spectrum in the frequency domain and
    converts it to pressure in the time domain through an inverse
    fast fourier transform."
function spl2pressure(spls)             # vector of sound pressure levels over the spectrum

    pref = 20e-6 # reference pressure
    p2s = pref^2 * 10 .^ (spls / 10.0)

    # inverse fourier transform
    ps = FFTW.ifft(p2s)

    return ps
end

"Takes in an unweighted SPL in decibels and frequency and returns
    the A-weighted SPL in decibels"
function aWeight(freq,                    # frequency
                 spl)                     # unweighted sound pressure level at each frequency

    if size(freq, 1) != size(spl, 1) # ensure vectors have same dimensions
        spl = spl'
    end

    wa = aWeight(freq)

    pref = 2e-5 #? may not need?

    p = pref^2 * 10.0 .^ (spl / 10) # convert to Pa
    pa = wa .* p # a-weight

    spla = 10 * log10.(pa / pref^2) # back to dB

    return spla
end

function aWeight(freq)

    # unchanging constants defined in wopwop user manual
    K1 = 2.243e16 #s^-4
    K3 = 1.562
    f1 = 20.599 # Hz
    f2 = 107.653 # Hz
    f3 = 737.862 # Hz
    f4 = 12194.22 # Hz

    wc = K1 * freq .^ 4 ./ ((freq .^ 2 .+ f1^2) .^ 2 .* (freq .^ 2 .+ f4^2) .^ 2) # C-weighting function
    wa = wc .* K3 .* freq .^ 4 ./ ((freq .^ 2 .+ f2^2) .* (freq .^ 2 .+ f3^2)) # A-weighting function
end

"Takes pressure output from wopwop and converts it to sound pressure
    level spectrum or a scalar overall sound pressure level if
    oasplFlag == true (for both cases, it is A-weighted if
    AweightingFlag == true).

Note 2: Currently does not include any windowing or filtering methods.
"
function pressure2SPL(time,                     # time history
                        p;                      # pressure (Pa) at each time step
                        AweightingFlag=false,   # SPL values are A-weighted
                        oasplFlag=false)        # returns scalar OASPL en lieu of SPL spectrum

    pref = 2e-5

    deltat = time[2] - time[1]
    T = time[end]
    deltaf = 1/T
    N = length(time)
    M = convert(Int64, floor(N/2 + 1))

    x = 0:M-1
    f = x * deltaf # frequencies

    Y = FFTW.fft(p)
    Y = Y[1:M] # cutoff frequencies higher than Nyquist frequency

    Pc = Y * 2 * deltat / T
    Pc[1] *= 0.5 # complex pressure

    p2local = 0.5 * abs2.(Pc) #mean square pressure
    spls = 10 * log10.(p2local / pref^2) #sound pressure level

    if AweightingFlag == true
        spls = aWeight(f, spls)
        if oasplFlag
            p2local = pref^2 * 10 .^ (spls / 10)
        end
    end

    if oasplFlag
        p2 = 1.0 * sum(p2local)
        oaspl = 10 * log10(p2 / pref^2)

        return oaspl
    else
        return f, spls
    end
end

"Same as function above, but weights the SPLs based on the frequency bin width.

Note: Intuitively this seems like it should work better. It will not match OASPL output from PSU-WOPWOP."
function pressure2SPL_weighted(time,                     # time history
                                p;                      # pressure (Pa) at each time step
                                AweightingFlag=false,   # SPL values are A-weighted
                                oasplFlag=false)        # returns scalar OASPL en lieu of SPL spectrum

    pref = 2e-5

    deltat = time[2] - time[1]
    T = time[end]
    deltaf = 1/T
    N = length(time)
    M = convert(Int64, floor(N/2 + 1))

    x = 0:M-1
    f = x * deltaf # frequencies

    Y = FFTW.fft(p)
    Y = Y[1:M] # cutoff frequencies higher than Nyquist frequency

    Pc = Y * 2 * deltat / T
    Pc[1] *= 0.5 # complex pressure

    p2local = 0.5 * abs2.(Pc) * deltaf #mean square pressure
    spls = 10 * log10.(p2local / pref^2) #sound pressure level

    if AweightingFlag == true
        spls = aWeight(f, spls)
        if oasplFlag
            p2local = pref^2 * 10 .^ (spls / 10)
        end
    end

    if oasplFlag
        p2 = 1.0 * sum(p2local)
        oaspl = 10 * log10(p2 / pref^2)

        return oaspl
    else
        return f, spls
    end
end

# """
#     `pressure_time2frequency(time, pressure, windowtype, SPLtype)`

# Converts pressure from time domain to frequency domain, calculates overall sound pressure level (OASPL)

# # ARGUMENTS
# * `time::Array{Real, 1}`        : Vector of time segments associated with each pressure.
# * `pressure::Array{Real, 1}`    : Pressure as a function of time.

# # OPTIONAL ARGUMENTS
# * `windowtype::String`          : Data windowing options:
#                                     "None", "Hann"
# * `windowtype::String`          : Specifies wether non-weighted or A-weighted SPL should be calculated. Options:
#                                     "Normal", "A-weighted"

# """
# function pressure_time2frequency(time::Array{Float64, 1}, pressure::Array{Float64, 1},
#                                     windowtype::String, SPLtype::String)

#     #Define constants
#     K1 = 2.243*10^16
#     K3 = 1.562
#     f1 = 20.599
#     f2 = 107.653
#     f3 = 737.862
#     f4 = 12194.22

#     #Time & frequency parameters
#     deltat = time[2] - time[1] #time spacing
#     T = time[end] #total time
#     deltaf = 1/T #frequency spacing
#     N = length(time) # number of time points
#     M = convert(Int64,floor(N/2 + 1)) # number of frequency bins of type Int64 for indexing later

#     # println("Pressure before windowing:")
#     # println(pressure)
#     # Window data if specified - treats data as periodic for fourier transform
#     if windowtype != "None"
#         pressure = FLOWNoise.windowdata!(time, pressure, windowtype)
#     end

#     # println("Pressure after windowing:")
#     # println(pressure)

#     # Use FFTW package to perform fast fourier transform on pressure data
#     Y = FFTW.fft(pressure)

#     # Initialize variables
#     Pc = zeros(Complex{Float64},M,1)
#     p2local = zeros(Float64,M,1)
#     SPLlocal_dB = zeros(Float64,M,1)
#     p2Alocal = zeros(Float64,M,1)
#     SPLlocal_dBA = zeros(Float64,M,1)
#     wc = ones(Float64,M,1)
#     wa = ones(Float64,M,1)
#     p2 = 0.0

#     x = 1:M
#     f = x*deltaf

#     for m in 1:M

#         if m > 1
#             Pc[m] = Y[m] * 2 * deltat / T #Complex pressure (magnitude and phase) at each frequency bin
#         else
#             Pc[m] = Y[m] * deltat / T
#         end

#         if SPLtype == "A-weighted" #if "Normal" then wa[m] = 1.0 and does not affect output

#             #A-weighted prep, from PSU-WOPWOP manual
#             wc[m] = K1 * f[m]^4 / ((f[m]^2 + f1^2)^2*(f[m]^2 + f4^2)^2)
#             wa[m] = wc[m] * K3 * f[m]^4 / ((f[m]^2 + f2^2)*(f[m]^2 + f3^2))

#         elseif SPLtype != "Normal"
#             error("Invalid SPLtype specified.")
#         end

#         # global p2 += wa[m] * abs2(Pc[m]) #overall mean square pressure
#         p2 += wa[m] * abs2(Pc[m]) #overall mean square pressure
#         p2local[m] = 0.5 * wa[m] * abs2(Pc[m]) #local mean square pressure for specific frequency
#         SPLlocal_dB[m] = 10.0 * log(10,p2local[m]/((20*10.0^-6)^2)) #local SPL

#     end

#     p2 *= 0.5

#     OASPL_dB = 10.0 * log(10,p2/(20*10.0^-6)^2) # overall sound pressure level in dB

#     return SPLlocal_dB, OASPL_dB
# end

"""
    `windowdata!(time, data, windowtype)`

Windows data, or "tapers" the ends to zero. This effectively makes the data periodic to avoid discontinuities in a fourier transform.

# ARGUMENTS
* `time::Array{Real, 1}`        : Vector of time segments associated with each data point in data.
* `data::Array{Real, 1}`        : Data (later to be transformed) as a function of time.

# OPTIONAL ARGUMENTS
* `windowtype::String`          : Data windowing options:
                                    "Hann", "Flat Top", "Blackman"

NOTES:
* Windowing changes the amplitude of the data in the time history, so a scaling factor must also be used, depending on which type of windowing is being performed.
* Windowing adds certain frequencies to the spectrum. Do not trust windoewed time signals which have frequencies closer than five bins away from both ends of the spectrum's frequency limits.
"""
function windowdata!(time, data, windowtype::String="Hann")
    T = time[end]
    tindex = 1
    # Apply windowing
    for t in time
        if windowtype == "Hann"
            wh = 1/2 - 1/2 * cos(2*pi*t/T)
            scalingfactor = 1.633
        elseif windowtype == "Flat Top"
            wh = 1 - 1.93 * cos(2*pi*t/T) + 1.29 * cos(4*pi*t/T) - 0.388 * cos(6*pi*t/T) + 0.0322 * cos(8*pi*t/T)
            scalingfactor = 0.515
        elseif windowtype == "Blackman"
            wh = 0.42 - 0.50 * cos(2*pi*t/T) + 0.08 * cos(4*pi*t/T)
            scalingfactor = 1.8119
        else
            error("Invalid windowtype defined.")
        end

        data[tindex] *= wh * scalingfactor
        # global tindex += 1
        tindex += 1
    end # TODO: broadcast this in the future

    return data

#TODO: currently not solving problem in case 2 not matching up
end

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
using Distributions: MvNormal
using LinearAlgebra: norm, I
using Statistics: mean
using Random

"""
Implementation of hard version of dpmm. DPMeans for small variance asymptotics
approach taken in this paper: 'Revisiting k-means: New Algorithms via Bayesian
Nonparametrics'.
"""

mutable struct Mixture
    components::Array{Distribution}
    weights::Array{Float64}
end

function likelihood(mixture::Mixture, point)
    sum([pdf(mixture.components[i], point) * mixture.weights[i] for i in 1:length(mixture.components)])
end

function sample(mixture::Mixture, n_samples=1)
    # Sample how many times each component is randomly selected.
    sampled_components = rand(Multinomial(n_samples, mixture.weights))
    samples = []
    for (ci, count) in enumerate(sampled_components)
        # sample from the component the number of times it was selected.
        for col in eachcol(rand(mixture.components[ci], count))
            push!(samples, col)
        end
    end
    return samples
end

mutable struct DPMeans
    lambda::Float64
    kernels::Array{MvNormal, 1}
    weights::Array{Float64, 1}
end

function DPMeans(lambda=0.5)
    lambda = lambda
    kernels = Array{MvNormal, 1}()
    DPMeans(lambda, kernels, Array{Float64, 1}())
end

"""
Train a DPMeans on the provided data. TODO: make this efficient (sub O(n)).
"""
function fit!(model::DPMeans, data, n_iterations=5, lambda_c=20, min_cluster_size=1)
    # If this is the first time fit was called, set the kernel to a random
    update_lambda=false   # point in the data.
    if size(model.kernels, 1) == 0
        push!(model.kernels, MvNormal(data[rand(1:size(data, 1)), :], model.lambda))
        push!(model.weights, 1.0)
    end

    labels = [0 for i in 1:size(data, 1)]
    for iteration in 1:n_iterations
        # For all observations.
        for i in 1:size(data, 1)
            # Calculate distances to all cluster centers.
            #   There has got to be a nicer way to do this.
            row_distances = [data[i, :] .- model.kernels[ki].μ for ki in 1:size(model.kernels, 1)]
            row_distances = [norm(row_distances[i, :]) for i in 1:size(row_distances, 1)]

            # Find the min and check if it is > lambda.
            min_index = argmin(row_distances)
            min_value = row_distances[min_index]
            if min_value > model.lambda
                # Add a new cluster to model.kernels.
                push!(model.kernels, MvNormal(vec(data[i, :]), model.lambda))
                push!(model.weights, 0.0)
            end
            labels[i] = min_index
        end

        # Move the cluster means to reflect the labels.
        # Remove any clusters that don't have points assigned to them.
        keep_inds = []
        for label in 1:maximum(labels)
            cluster_data = data[labels .== label, :]
            covariance = model.lambda
            if size(cluster_data, 1) == 1 || sum(covariance) == 0.0
                covariance = I * model.lambda
            end
            model.kernels[label] = MvNormal(vec(mean(cluster_data, dims=1)), covariance)
            model.weights[label] = size(cluster_data, 1)

            # We can't have clusters with one thing in them.
            if length(cluster_data) > min_cluster_size
                push!(keep_inds, label)
            end

        end
        model.kernels = model.kernels[keep_inds]
        model.weights = model.weights[keep_inds] / sum(model.weights[keep_inds])

        # Now make a mixture
        mixture = Mixture(model.kernels, model.weights)
    end
end

function get_labels(model::DPMeans, data)
    labels = [0 for i in 1:size(data, 1)]
    for i in 1:size(data, 1)
       # Calculate distances to all cluster centers.
       #   There has got to be a nicer way to do this.
        row_distances = [data[i, :] .- model.kernels[ki].μ for ki in 1:size(model.kernels, 1)]
        row_distances = [norm(row_distances[i, :]) for i in 1:size(row_distances, 1)]

        # Find the min and check if it is > lambda.
        min_index = argmin(row_distances)
        min_value = row_distances[min_index]
        labels[i] = min_index
    end
    labels
end

function get_mixture(model::DPMeans)
    return Mixture(model.kernels, model.weights)
end

max_stars_repo_name: arider/DPMeans

max_issues_repo_name: arider/DPMeans

max_forks_repo_name: arider/DPMeans


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -16.424759279035918 <= q <= 20.23452089871623)

                     add_NL_constraint(m, :(log(1 + exp(0.8420310557117703 + 0.34198141401583015*log(1 + exp(-0.5389217244927917 + 0.5798025520034438*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + 0.4428632765655749*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + 0.9033224385695577*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + -0.1887519896617813*log(1 + exp(-0.23921485861490321 + 0.5030319708107833*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.2081481501163327*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.2939216839670169*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + -0.2882490310085508*log(1 + exp(-0.6545847926105979 + 0.9774628615021363*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.21180887804098658*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.7013494059895882*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + log(1 + exp(-0.35188345522271236 + -0.11415509244934263*log(1 + exp(-0.5389217244927917 + 0.5798025520034438*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + 0.4428632765655749*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + 0.9033224385695577*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + -0.8088728193024375*log(1 + exp(-0.23921485861490321 + 0.5030319708107833*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.2081481501163327*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.2939216839670169*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + 0.8271426964720727*log(1 + exp(-0.6545847926105979 + 0.9774628615021363*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.21180887804098658*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.7013494059895882*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + log(1 + exp(0.4798469380896484 + 0.8134725364327324*log(1 + exp(-0.5389217244927917 + 0.5798025520034438*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + 0.4428632765655749*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + 0.9033224385695577*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + -0.12585836968848696*log(1 + exp(-0.23921485861490321 + 0.5030319708107833*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.2081481501163327*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.2939216839670169*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))) + 0.613011911779838*log(1 + exp(-0.6545847926105979 + 0.9774628615021363*log(1 + exp(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))) + -0.21180887804098658*log(1 + exp(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))) + -0.7013494059895882*log(1 + exp(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#!/home/kimura/Julia1027/julia/julia

const MECHAJYO_VERSION = "0.7"

function setupIO()
    stdinx::Ptr{Uint8} = 0
    stdoutx::Ptr{Uint8} = 0
    ret = ccall((:setunbuffering, "../lib/libMJ.so.1"), Int32, ())
end

setupIO()
s = readline(STDIN)
s = chomp(s)
bookfile = "./Joseki.db"
usebook  = true

if s == "usi"
    println("id name Mecha Jyoshi Shogi ",MECHAJYO_VERSION)
    println("id author Sayuri TAKEBE, Mio WATANABE, Rieko TSUJI and Takeshi KIMURA")
    println("option name BookFile type string default $(bookfile)")
    println("option name UseBook type check default $(usebook)")
    println("usiok");
end

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
using ForwardDiff: derivative

"""
    DerivativeMean{T∇<:MeanFunction} <: MeanFunction

"""
struct DerivativeMean{T∇<:MeanFunction} <: MeanFunction
    m::T∇
end
(μ::DerivativeMean)(x::Real) = derivative(x->μ.m(x), x)

"""
    DerivativeKernel{Tk<:Kernel} <: Kernel

"""
struct DerivativeKernel{Tk<:Kernel} <: Kernel
    k::Tk
end
function (k::DerivativeKernel)(x::Real, x′::Real)
    return derivative(x′->derivative(x->k.k(x, x′), x), x′)
end
(k::DerivativeKernel)(x::Real) = k(x, x)

"""
    DerivativeLhsCross{Tk<:CrossKernel} <: CrossKernel

"""
struct DerivativeLhsCross{Tk<:CrossKernel} <: CrossKernel
    k::Tk
end
function (k::DerivativeLhsCross)(x::Real, x′::Real)
    return derivative(x->k.k(x, x′), x)
end

"""
    DerivativeRhsCross{Tk<:CrossKernel} <: CrossKernel

"""
struct DerivativeRhsCross{Tk<:CrossKernel} <: CrossKernel
    k::Tk
end
function (k::DerivativeRhsCross)(x::Real, x′::Real)
    return derivative(x′->k.k(x, x′), x′)
end

"""
    DerivativeCross{Tk<:CrossKernel} <: CrossKernel

"""
struct DerivativeCross{Tk<:CrossKernel} <: CrossKernel
    k::Tk
end
function (k::DerivativeCross)(x::Real, x′::Real)
    return derivative(x′->derivative(x->k.k(x, x′), x), x′)
end
(k::DerivativeCross)(x::Real) = k(x, x)


text:
module TabularTMazeEvalSet

using Curiosity
using MinimalRLCore
const TTMU = Curiosity.TabularTMazeUtils
include("../experiments/tabular_tmaze.jl")
using Random
using StatsBase
using GVFHordes
using Statistics
using JLD2

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(s, a) for a in actions]))

function gen_dataset()
    parsed = TabularTMazeExperiment.default_args()
    parsed["cumulant_schedule"] = "Constant"
    parsed["cumulant"] = 1.0

    horde = TabularTMazeExperiment.get_horde(parsed, 21, 4, nothing)

    cumulant_schedule = TTMU.get_cumulant_schedule(parsed)

    exploring_starts = parsed["exploring_starts"]
    env = TabularTMaze(exploring_starts, cumulant_schedule)
    start_states = Curiosity.valid_state_mask()

    function del!(arr, ind_to_del)
        arr = arr[1:end .!= ind_to_del]
        return arr
    end
    #NOTE: This lets the agent start in terminal states. This is not normally possible.
    # Work backwards with index so you dont have to reindex due to shifting
    ind_to_delete = [21, 17, 5, 1]
    for i in ind_to_delete
        start_states = del!(start_states, i)
    end
    observations = []
    start_states_all = []
    for s in start_states
        MinimalRLCore.reset!(env,s)
        push!(observations, MinimalRLCore.get_state(env))
        push!(start_states_all, s)
    end
    total_observations = []
    total_start_states = []
    total_actions = []
    for (obs_ind,obs) in enumerate(observations)
            unique_actions = unique([StatsBase.sample(GVFHordes.policy(gvf), obs, get_actions(env)) for gvf in horde.gvfs])
            for ua in unique_actions
                push!(total_actions, ua)
                push!(total_observations, obs)
                push!(total_start_states, start_states[obs_ind])
            end
    end

    num_returns = 2
    γ_thresh=1e-6
    actions = zeros(Int,length(horde.gvfs), length(total_observations))
    horde_rets = zeros(length(horde.gvfs), length(total_observations))
    for (gvf_i,gvf) in enumerate(horde.gvfs)
        rets = monte_carlo_returns(env, gvf, total_start_states, total_actions, num_returns, γ_thresh)
        rets_avg = mean(hcat(rets...),dims=1) # stack horizontally and then average over runs
        horde_rets[gvf_i,:] = rets_avg
    end

    return horde_rets, total_observations, total_actions
end

function save_data(rets,obs,actions)
    TTMazeEvalSet = Dict()
    TTMazeEvalSet["ests"] = rets
    TTMazeEvalSet["actions"] = actions
    TTMazeEvalSet["states"] = obs
    @save "./src/data/TTMazeEvalSet.jld2" TTMazeEvalSet
end

function gen_roundrobin_dataset(num_start_states = 100)
    parsed = TabularTMazeExperiment.default_args()
    parsed["cumulant_schedule"] = "Constant"
    parsed["cumulant"] = 1.0

    horde = TabularTMazeExperiment.get_horde(parsed, 21, 4, nothing)

    cumulant_schedule = TTMU.get_cumulant_schedule(parsed)

    exploring_starts = true
    env = TabularTMaze(exploring_starts, cumulant_schedule)

    all_states = []
    all_actions = []
    for gvf_i in 1:4
        # policy = gvf.policy
        policy = TTMU.GoalPolicy(gvf_i)
        for i in 1:num_start_states
            s = MinimalRLCore.start!(env)
            a = policy(s)
            push!(all_states, s)
            push!(all_actions, a)
            term = false
            while term == false
                s, r, term = MinimalRLCore.step!(env, a)
                if term
                    break
                end
                a = policy(s)
                push!(all_states, s)
                push!(all_actions, a)
            end
        end
    end
    eval_set_size = 500
    indices = sample(1:size(all_states)[1], eval_set_size, replace=false)
    ss_eval = all_states[indices]
    as_eval = all_actions[indices]
    num_returns = 2
    γ_thresh=1e-6

    horde_rets = zeros(length(horde.gvfs), length(ss_eval))

    for (gvf_i,gvf) in enumerate(horde.gvfs)
        ss_modified = [Int(s[1]) for s in ss_eval]
        rets = monte_carlo_returns(env, gvf, ss_modified, as_eval, num_returns, γ_thresh)
        rets_avg = mean(hcat(rets...),dims=1) # stack horizontally and then average over runs
        horde_rets[gvf_i,:] = rets_avg
    end
    return horde_rets, ss_eval, as_eval

end

function save_data_roundrobin(rets,obs,actions)
    TTMazeRoundRobinEvalSet = Dict()
    TTMazeRoundRobinEvalSet["ests"] = rets
    TTMazeRoundRobinEvalSet["actions"] = actions
    TTMazeRoundRobinEvalSet["states"] = obs
    @save "./src/data/TTMazeRoundRobinEvalSet.jld2" TTMazeRoundRobinEvalSet
end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function lint_helper(ex::Expr, ctx)
    if ex.head == :macrocall
        if ex.args[1] == Symbol("@defcomp")
            push!(ctx.callstack[end].types, ex.args[2])
            return true
        end
    end
    return false
end


text:
@testset "VTKOutput" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Float64}(mesh)

    u = addVariable!(diffusion_system, "u")

    diffusion_kernel = Diffusion(u)
    addKernel!(diffusion_system, diffusion_kernel)

    left_boundary = DirichletBC(u, [4], 0.0)
    addBC!(diffusion_system, left_boundary)

    right_boundary = DirichletBC(u, [2], 1.0)
    addBC!(diffusion_system, right_boundary)

    initialize!(diffusion_system)

    solver = JuliaDenseImplicitSolver(diffusion_system)

    solve!(solver)

    out = VTKOutput()

    rm("test_out.vtu", force=true)

    output(out, solver, "test_out")

    @test isfile("test_out.vtu")
end

@testset "VTKOutputFD" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Dual{4,Float64}}(mesh)

    u = addVariable!(diffusion_system, "u")

    diffusion_kernel = Diffusion(u)
    addKernel!(diffusion_system, diffusion_kernel)

    left_boundary = DirichletBC(u, [4], 0.0)
    addBC!(diffusion_system, left_boundary)

    right_boundary = DirichletBC(u, [2], 1.0)
    addBC!(diffusion_system, right_boundary)

    initialize!(diffusion_system)

    solver = JuliaDenseImplicitSolver(diffusion_system)

    solve!(solver)

    out = VTKOutput()

    rm("test_out.vtu", force=true)

    output(out, solver, "test_out")

    @test isfile("test_out.vtu")
end


text:
channel = RemoteChannel()

function asynccall(f)
	@async begin # if we use multiple processes
	put!(channel, true)
	f()
	take!(channel)
	end
end

println("START CHANNEL")

asynccall(()->println("1"))
asynccall(()->println("2"))
asynccall(()->println("3"))
asynccall(()->println("4"))

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
"""
	import_project_information(config::REDCap.Config, data; format::String="json")

Update basic attributes of given REDCap project.
NOTE: Only for projects in development

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `format` - "json", "xml", "csv", or "odm". declares format of imported data

#### Returns:
Number of successfully imported values
"""
function import_project_information(config::REDCap.Config, data; format::String="json", returnFormat::String="json")
	return api_pusher("import", "project_settings", config, data = import_file_checker(data, format), format=format, returnFormat=returnFormat)
end

"""
	import_metadata(config::REDCap.Config, data; format::String="json", returnFormat::String="json")

Import metadata (i.e., Data Dictionary) into a project.
NOTE: Only for projects in development

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `returnFormat` - Error message format

#### Returns:
Number of successfully imported fields
"""
###BROKEN(?)###
function import_metadata(config::REDCap.Config, data; format::String="json", returnFormat::String="json")
	return api_pusher("import", "metadata", config, data = import_file_checker(data, format), format=format, returnFormat=returnFormat)
end
#=
Breaks on JSON, XML,
Works for CSV

newmeta = Dict("required_field"=>"",
  "section_header"=>"",
  "matrix_ranking"=>"",
  "select_choices_or_calculations"=>"",
  "field_type"=>"file",
  "field_note"=>"",
  "form_name"=>"demographics",
  "matrix_group_name"=>"",
  "field_label"=>"File Upload",
  "custom_alignment"=>"",
  "question_number"=>"",
  "text_validation_max"=>"",
  "text_validation_type_or_show_slider_number"=>"",
  "branching_logic"=>"",
  "field_annotation"=>"",
  "identifier"=>"",
  "text_validation_min"=>"",
  "field_name"=>"file_upload")

push!(meta, newmeta)

ERROR: HTTP.ExceptionRequest.StatusError(400, HTTP.Messages.Response:
"""
HTTP/1.1 400 Bad Request
Date: Wed, 25 Jul 2018 18:12:00 GMT
Server: Apache/2.4.6 (Red Hat Enterprise Linux)
X-Powered-By: PHP/5.4.16
Expires: 0
cache-control: no-store, no-cache, must-revalidate
Pragma: no-cache
Access-Control-Allow-Origin: *
REDCap-Random-Text: VxGV9VFPeaU9yzkMsFXFKUQuhMYCygnBXEgejwoi2wXnX
Content-Length: 343
Connection: close
Content-Type: application/json; charset=utf-8

{"error":"Each row MUST have a variable\/field name, but the following cells have a field name missing: A2.\nEach row MUST have a form name, but the following cells have a form name missing: B2.\nA field validation type is required in order to have a minimum or maximum validation value. The following cells are missing a validation type: I2"}""")
=#

"""
	import_users(config::REDCap.Config, data; format::String="json", returnFormat::String="json")

Update/import new users into a project.

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `returnFormat` - Error message format

#### Returns:
Number of succesfully added/modified users.
"""
###BROKEN###
function import_users(config::REDCap.Config, data; format::String="json", returnFormat::String="json")
	return api_pusher("import", "user", config, data = import_file_checker(data, format), format=format, returnFormat=returnFormat)
end

#=
No longer modifies users! Can add them wholesale, not modify them after???

=#

"""
	import_arms(config::REDCap.Config, data; override::Int=0, format::String="json", returnFormat::String="json")

Update/import Arms into a project.

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `override` - 0 (false) 1 (true) - overwrites existing arms
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `returnFormat` - Error message format

#### Returns:
Number of successfully imported arms
"""
###BROKEN###
function import_arms(config::REDCap.Config, data; override::Int=0, format::String="json", returnFormat::String="json")
	return api_pusher("import", "arm", config, data = import_file_checker(data, format), override=override, format=format, returnFormat=returnFormat)
end

#=
julia> arms = export_arms(config)
POSTing
POSTd
1-element Array{Any,1}:
 Dict{String,Any}(Pair{String,Any}("name", "Arm 1"),Pair{String,Any}("arm_num", 1))

julia> newarm
Dict{String,Any} with 2 entries:
  "name"    => "Arm 1"
  "arm_num" => 1

julia> push!(arms, newarm)
2-element Array{Any,1}:
 Dict{String,Any}(Pair{String,Any}("name", "Arm 1"),Pair{String,Any}("arm_num", "1"))
 Dict{String,Any}(Pair{String,Any}("name", "Arm 2"),Pair{String,Any}("arm_num", "2"))

julia> import_arms(config, arms)
POSTing
POSTd
2	<- this should indicate 2 arms added.

julia> export_arms(config)
POSTing
POSTd
1-element Array{Any,1}:
 Dict{String,Any}(Pair{String,Any}("name", "Arm 1"),Pair{String,Any}("arm_num", 1))

=#

"""
	import_events(config::REDCap.Config, data; override::Int=0, format::String="json", returnFormat::String="json")

Update/import Events into a project.

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `override` - 0 (false) 1 (true) - overwrites existing events
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `returnFormat` - Error message format

#### Returns:
Number of successfully imported events
"""
function import_events(config::REDCap.Config, data; override::Int=0, format::String="json", returnFormat::String="json")
	return api_pusher("import", "event", config, data = import_file_checker(data, format), override=override, format=format, returnFormat=returnFormat)
end

"""
	import_records(config::REDCap.Config, data::Any; format::String="json", dtype::String="flat", overwriteBehavior::String="normal", forceAutoNumber::Bool=false, dateFormat::String="YMD", returnContent::String="count", returnFormat::String="json")

Import a set of records for a project.

#### Parameters:
* `config` - Struct containing url and api-key
* `recordData` - Array of record data to be imported - pass as a file location to import from disk
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `dtype` - "flat" (one record per row) or "eav" (one data point per row)
* `overwriteBehavior` - "normal" - will not overwrite, "overwrite" - will
* `forceAutoNumber` - Force auto-numbering and overwrite given id number
* `dateFormat` - "YMD", "MDY", or "DMY"
* `returnContent` - "count" (number of successfully uploaded records),
						"ids" (list of record numbers imported),
						"auto-ids" (pair of assigned id and given id)
* `returnFormat` - Error message format

#### Returns:
Specified by returnContent
"""
function import_records(config::REDCap.Config, data; format::String="json", dtype::String="flat", overwriteBehavior::String="normal", forceAutoNumber::Bool=false, dateFormat::String="YMD", returnContent::String="count", returnFormat::String="json")
	return api_pusher("import", "record", config, data = import_file_checker(data, format), format=format, dtype=dtype, overwriteBehavior=overwriteBehavior, forceAutoNumber=forceAutoNumber, dateFormat=dateFormat, returnContent=returnContent, returnFormat=returnFormat)
end

"""
	import_instrument_event_mappings(config::REDCap.Config, data; format::String="json", returnFormat::String="json")

Import Instrument-Event Mappings into a project

#### NOTE: This only works for longitudinal projects.

#### Parameters:
* `config` - Struct containing url and api-key
* `data` - Data to be imported - pass as a file location to import from disk
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `returnFormat` - Error message format

#### Returns:
Number of successfully imported inst-event mappings
"""
###BROKEN(?)###
function import_instrument_event_mappings(config::REDCap.Config, data; format::String="json", returnFormat::String="json")
	return api_pusher("import", "formEventMapping", config, data = import_file_checker(data, format), format=format, returnFormat=returnFormat)
end

"""
	import_file(config::REDCap.Config, record::String, field::String, event::String, file::String; repeat_instance::Int=1, returnFormat::String="json")

Upload a document to specific record to the designated uploading field.

#### Parameters:
* `config` - Struct containing url and api-key
* `record` - Destination record id
* `field` - Destination file upload field
* `event` - Destination event
* `repeat_instance` - Number of repeated instances (long project)
* `file` - File to be imported
* `returnFormat` - Error message format

#### Returns:
Nothing/errors
"""
###BROKEN###
function import_file(config::REDCap.Config, record::String, field::String, event::String, file::String; repeat_instance::Int=1, returnFormat::String="json")
	return api_pusher("import", "file", config, record=record, field=field, event=event, file=open(file), repeat_instance=repeat_instance, returnFormat=returnFormat)
end
#=
No longer properly passes files- REDCap returns an error of invalid file.
=#

"""
	create_project(config::REDCap.Config, project_title::String, purpose::Integer; format::String="json", returnFormat::String="json", odm="", purpose_other::String="", project_notes::String="", is_longitudinal::Integer=0, surveys_enabled::Integer=0, record_autonumbering_enabled::Integer=1)

Creates a project with the given parameters

#### Parameters:
* `config` - Struct containing url and super-api-key
* `format` - "json", "xml", "csv", or "odm". declares format of imported data
* `data` - Attributes of project to create- only project_title and purpose are required (* for default)
	* `project_title`: Title
	* `purpose`: Must be numerical (0 - test, 1 - other, 2 - research, 3 - Qual+, 4 - OpSupport)
	* `purpose_other`: If purpose == 1, string of purpose
	* `project_notes`: Notes
	* `is_longitudinal`: 0 - false*, 1 - true
	* `surveys_enabled`: 0 - false*, 1 - true
	* `record_autonumbering_enabled`: 0 - false, 1 - true*
* `returnFormat` - Error message format
* `odm` - XML string containing metadata

#### Returns:
The standard config for that project.
"""
function create_project(config::REDCap.Config, project_title::String, purpose::Integer; format::String="json", returnFormat::String="json", odm="", purpose_other::String="", project_notes::String="", is_longitudinal::Integer=0, surveys_enabled::Integer=0, record_autonumbering_enabled::Integer=1)
	if length(config.key)==64
		data = json_formatter([Dict{String, Any}("project_title" => project_title,
													"purpose" => purpose,
													"purpose_other" => purpose_other,
													"project_notes" => project_notes,
													"is_longitudinal" => is_longitudinal,
													"surveys_enabled" => surveys_enabled,
													"record_autonumbering_enabled" => record_autonumbering_enabled)], "import")
		#Send through api_pusher NOT poster
		response = api_pusher("import", "project", config, format=format, data=data, returnFormat=returnFormat, odm=odm)
		return Config(config.url, response; ssl=config.ssl) #inherit all settings except the newly generated key
	else
		@error("Please use a config object that contains a properly entered Super API key.\n$(config.key) is an invalid Super-API key.")
	end
end


text:
include("faber.jl")
include("types.jl")

"""
    acc_rej(n::Int64, S::System, X::AbstractModel, f1::Regular, u, v)

For Regular sampling_scheme: always accept
"""
function acc_rej(n::Int64, S::System, X::AbstractModel, f1::Regular, u::Float64, v::Float64, t::Float64)
    return true
end
"""
    acc_rej(n::Int64, S::System, X::AbstractModel, f1::SubSampling, u, v)

Subsampling scheme: accept the event time as coming from the real Poisson rate.
"""
function acc_rej(n::Int64, S::System, X::AbstractModel, ::SubSampling, u::Float64, v::Float64, t::Float64)
    if λratio(n, S, X, u, v, t)> rand()
        return true
    else
        return false
    end
end

"""
    update_events!(n, S::System, X::AbstractModel, ::SubSampling, ::FullIndependence, acc)

renovate time when the nth coefficient rings (i.e. τ_n = 0). In the full independece case
we do not need to simulate any other time, but just rescale them.
"""
function update_events!(n::Int64, τ0::Float64, S::System, X::AbstractModel, ::SubSampling, ::FullIndependence, acc::Bool, u::Float64, v::Float64, t::Float64)
    S.τ[n] = λbar(n, S, X , u, v, t)
    for i in 1:(n-1)
        S.τ[i] -= τ0
    end
    for i in  (n + 1): length(S.ϕ)
        S.τ[i] -=  τ0
    end
end

function update_events!(n::Int64, τ0::Float64, S::System, X::AbstractModel, f1::Regular, f2::PartialIndependence, acc::Bool, u::Float64, v::Float64, t::Float64)
    for i in S.ϕ[n].nhb
        S.τ[i] = λbar(i, S, X , u, v, t)
    end
    for i in S.ϕ[n].notnhb
        S.τ[i] -= τ0
    end
end

"""
    update_events!(n, S::System, X::AbstractModel, f1::Regular, f2::PartialIndependence, acc)

renovate time when the `n`_th coefficient rings (i.e. τ_n = 0). In the partial independece case
we need to simulate all the time relative to the functions sharing the same support and rescale the others.
"""
function update_events!(n::Int64, τ0::Float64, S::System, X::AbstractModel, f1::Regular, f2::PartialIndependence, acc::Bool, u::Float64, v::Float64, t::Float64)
    for i in S.ϕ[n].nhb
        S.τ[i] = λbar(i, S, X , u, v, t)
    end
    for i in S.ϕ[n].notnhb
        S.τ[i] -= τ0
    end
end

"""
    update_events!(n, S::System, X::AbstractModel, f1::SubSampling, f2::PartialIndependence, acc)

renovate time when the `n`_th coefficient rings (i.e. τ_n = 0). In the partial independece case
we need to simulate all the time relative to the functions sharing the same support and rescale the others.
if `acc` is false, means that in the previous step we rejected the event and we did not change velocity.
This implies that we need to renovate only the `n`th time like in the Full independece case.
"""
function update_events!(n::Int64, τ0::Float64, S::System, X::AbstractModel, ::SubSampling, ::PartialIndependence, acc::Bool, u::Float64, v::Float64, t::Float64)
    if acc == true
        update_events!(n, τ0, S, X, Regular(), PartialIndependence(), acc, u, v, t)
    else
        update_events!(n, τ0, S, X, SubSampling(), FullIndependence(), acc, u, v, t)
    end
end

"""
    zigzagsampler(X::AbstractModel, T, L, u, v, TT)

run the ZigZag sampler for diffusion `X` conditioned to start at `u` and
end at `v` at time `T`. The infinite summation is truncated at level `L` and
the ZigZag process will run up to time `TT`. Optionally set velocities θ
"""
function zz_sampler(X::AbstractModel, T::Float64, L::Int64, u::Float64, v::Float64, clock::Float64, ξ =fill(0.0, 2<<L - 1) , θ = fill(1.0, 2<<L - 1))
    t = 0.0
    S = System(L, T, ξ, θ)
    τ0 = 0.0
    n0 = 0
    Ξ = Skeleton[]
    while t < clock
        τ0 , n0 = findmin(S.τ)
        S.ξ .+=  S.θ*τ0
        t += τ0
        if acc_rej(n0, S, X, sampling_scheme(X), u, v, t)
            S.θ[n0] = -S.θ[n0]
            push!(Ξ, (Skeleton(copy(S.ξ), t)))
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), true, u, v, t)
        else
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), false, u, v, t)
        end
    end
    return Ξ
end

function zz_sampler_count(X::AbstractModel, T::Float64, L::Int64, u::Float64, v::Float64, clock::Float64, ξ =fill(0.0, 2<<L - 1), θ = fill(1.0, 2<<L - 1))
    t = 0.0
    S = System(L, T, ξ, θ)
    τ0 = 0.0
    n0 = 0
    Ξ = Skeleton[]
    num_events = fill(0, 2^(L+1) -1)
    while t < clock
        τ0 , n0 = findmin(S.τ)
        S.ξ .+=  S.θ*τ0
        t += τ0
        if acc_rej(n0, S, X, sampling_scheme(X), u, v, t)
            S.θ[n0] = -S.θ[n0]
            num_events[n0] += 1
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), true, u, v, t)
        else
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), false, u, v, t)
        end
    end
    return num_events
end

function zz_sampler_count_ub(X::AbstractModel, T::Float64, L::Int64, u::Float64, v::Float64, clock::Float64, ξ =fill(0.0, 2<<L - 1), θ = fill(1.0, 2<<L - 1))
    t = 0.0
    S = System(L, T, ξ, θ)
    τ0 = 0.0
    n0 = 0
    Ξ = Skeleton[]
    num_events = fill(0, 2^(L+1) -1)
    while t < clock
        τ0 , n0 = findmin(S.τ)
        S.ξ .+=  S.θ*τ0
        t += τ0
        num_events[n0] += 1
        if acc_rej(n0, S, X, sampling_scheme(X), u, v, t)
            S.θ[n0] = -S.θ[n0]
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), true, u, v, t)
        else
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), false, u, v, t)

        end
    end
    return num_events
end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
module LegibleLambdas

using MacroTools: postwalk
export @λ, @lambda, LegibleLambda

"""
    @λ <lambda definition>

Create legible lambdas.

# Example

```julia
julia> @λ(x -> g(x)/3)
(x -> g(x)/3)
```
"""
macro λ(ex)
    return parse_lambda(ex)
end

eval(:(const $(Symbol("@lambda")) = $(Symbol("@λ"))))

"""
    parse_lambda(ex)

Parse a lambda expression to `LegibleLambda` constructor expression.
"""
function parse_lambda(ex)
    if ex.head == :(->)
        :(LegibleLambda($(QuoteNode(ex)), $(esc(ex))))
    else
        throw("Must be called on a Lambda expression")
    end
end

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

struct LegibleLambda{F <: Function} <: Function
    ex::Expr
    λ::F
end

(f::LegibleLambda)(args...) = (f.λ)(args...)

function Base.show(io::IO, f::LegibleLambda)
    show(io, MIME"text/plain"(), f)
end

function Base.show(io::IO, ::MIME"text/plain", f::LegibleLambda)
    l_names = propertynames(f.λ)
    ex = f.ex
    foreach(l_names) do name
        ex = postwalk(ex) do s
            s === name ? getfield(f.λ, name) : s
        end
    end
    ex_cut = ex |> (ex -> postwalk(cutlnn, ex)) |> (ex -> postwalk(cutblock, ex))
    name = replace((repr(ex_cut |> pretty_kwargs)[2:end]), "->" => " -> ")
    print(io, name)
end

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

cutlnn(x) = x
function cutlnn(ex::Expr)
    Expr(ex.head, ex.args[[!(arg isa LineNumberNode) for arg in ex.args]]...)
end

cutblock(x) = x
function cutblock(ex::Expr)
    if (ex.head == :block) && (length(ex.args)==1)
        ex.args[1]
    else
        ex
    end
end

tupleargs(x::Symbol) = (x,)
function tupleargs(ex::Expr)
    if ex.head == :tuple
        Tuple(ex.args)
    else
        throw("expression must have head `:tuple`")
    end
end

pretty_kwargs(x::Symbol) = x
function pretty_kwargs(ex)
    if ex.args[1] isa Expr
        if ex.args[1].head == :block
            exout = copy(ex)
            exout.args[1].head = :tuple
            exout
        else
            ex
        end
    else
        ex
    end
end

end


text:
# Test classified collections

using DataStructures
using Base.Test

# classified lists

c = classified_lists(ASCIIString, Int)

push!(c, "low", 1)
push!(c, "low", 2)
push!(c, "low", 3)
push!(c, "high", 4)
push!(c, "high", 5)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

@test c["low"] == [1, 2, 3]
@test c["high"] == [4, 5]

# classified sets

c = classified_sets(ASCIIString, Int)

push!(c, "low", 1)
push!(c, "low", 2)
push!(c, "low", 3)
push!(c, "low", 1)
push!(c, "low", 2)

push!(c, "high", 4)
push!(c, "high", 5)
push!(c, "high", 5)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

@test isa(c["low"], Set{Int})
@test isa(c["high"], Set{Int})

@test sort(collect(c["low"])) == [1, 2, 3]
@test sort(collect(c["high"])) == [4, 5]

# classified counters

c = classified_counters(ASCIIString, Float64)

push!(c, "low", 1.)
push!(c, "low", 2.)
push!(c, "low", 3.)
push!(c, "low", 1.)
push!(c, "low", 2.)
push!(c, "low", 2.)

push!(c, "high", 4.)
push!(c, "high", 5.)
push!(c, "high", 5.)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

cl = c["low"]
ch = c["high"]

@test isa(cl, Accumulator{Float64, Int})
@test isa(ch, Accumulator{Float64, Int})

@test cl[1.] == 2
@test cl[2.] == 3
@test cl[3.] == 1
@test ch[4.] == 1
@test ch[5.] == 2


text:
push!(LOAD_PATH,"../src/")

using Documenter, JournaledJets

makedocs(
    sitename="JournaledJets.jl",
    modules=[JournaledJets],
    pages = [
        "Home" => "index.md",
        "reference.md",
        ]
)

deploydocs(
    repo = "github.com/ChevronETC/JournaledJets.jl"
)


text:
# Autogenerated wrapper script for liblinear_jll for x86_64-apple-darwin14
export predict, liblinear, train

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `liblinear`
const liblinear_splitpath = ["lib", "liblinear.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
liblinear_path = ""

# liblinear-specific global declaration
# This will be filled out by __init__()
liblinear_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const liblinear = "@rpath/liblinear.dylib"

# Relative path to `predict`
const predict_splitpath = ["bin", "predict"]

# This will be filled out by __init__() for all products, as it must be done at runtime
predict_path = ""

# predict-specific global declaration
function predict(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(predict_path)
    end
end

# Relative path to `train`
const train_splitpath = ["bin", "train"]

# This will be filled out by __init__() for all products, as it must be done at runtime
train_path = ""

# train-specific global declaration
function train(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(train_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"liblinear")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global liblinear_path = normpath(joinpath(artifact_dir, liblinear_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global liblinear_handle = dlopen(liblinear_path)
    push!(LIBPATH_list, dirname(liblinear_path))

    global predict_path = normpath(joinpath(artifact_dir, predict_splitpath...))

    push!(PATH_list, dirname(predict_path))
    global train_path = normpath(joinpath(artifact_dir, train_splitpath...))

    push!(PATH_list, dirname(train_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# Copyright (c) 2013: Cbc.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Cbc

import Cbc_jll
import MathOptInterface
import SparseArrays

function __init__()
    global libcbcsolver = Cbc_jll.libcbcsolver
    version_str = unsafe_string(Cbc_getVersion())
    version = if version_str == "devel"
        # Support un-released versions of Cbc. These may differ in C API
        # compatibility! Use at your own peril.
        v"2.10.5"
    else
        VersionNumber(version_str)
    end
    if !(v"2.10.0" <= version <= v"2.10.5")
        error("""
        You have installed version $version of Cbc, which is not supported by
        Cbc.jl We require Cbc version 2.10. After installing Cbc 2.10, run:

            import Pkg
            Pkg.rm("Cbc")
            Pkg.add("Cbc")

        If you have a newer version of Cbc installed, changes may need to be made
        to the Julia code. Please open an issue at
        https://github.com/jump-dev/Cbc.jl.
        """)
    end
    return
end

include("gen/libcbc_common.jl")
include("gen/libcbc_api.jl")
include("MOI_wrapper/MOI_wrapper.jl")

# Cbc exports all `Cbc_xxx` symbols. If you don't want all of these symbols in
# your environment, then use `import Cbc` instead of `using Cbc`.

for sym in filter(s -> startswith("$s", "Cbc_"), names(@__MODULE__, all = true))
    @eval export $sym
end

end


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using ParticleFilters

frames = Frames(MIME("image/png"), fps=2)

pomdp = VDPTagPOMDP()

filter = SIRParticleFilter(pomdp, 10000, rng=MersenneTwister(100))

hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(1), show_progress=true)
hist = simulate(hr, pomdp, RandomPolicy(pomdp, rng=MersenneTwister(2)), filter)

gr()
@showprogress "Creating gif..." for i in 1:length(hist)
    push!(frames, plot(pomdp, view(hist, 1:i)))
end

filename = string(tempname(), "_vdprun.gif")
write(filename, frames)
println(filename)
run(`setsid gifview $filename`)


text:
#####
##### First order fluxes
#####

# Mass
eq_tends(::Mass, ::AtmosLinearModel, ::Flux{FirstOrder}) = (Advect(),)

# Momentum
eq_tends(::Momentum, ::AtmosLinearModel, ::Flux{FirstOrder}) =
    (LinearPressureGradient(),)

# Energy
eq_tends(::Energy, m::AtmosLinearModel, tt::Flux{FirstOrder}) =
    (LinearEnergyFlux(),)

# AbstractMoistureVariable
# TODO: Is this right?
eq_tends(::AbstractMoistureVariable, ::AtmosLinearModel, ::Flux{FirstOrder}) =
    ()

# Tracers
eq_tends(::Tracers{N}, ::AtmosLinearModel, ::Flux{FirstOrder}) where {N} = ()

#####
##### Second order fluxes
#####

eq_tends(pv::PV, ::AtmosLinearModel, ::Flux{SecondOrder}) where {PV} = ()

#####
##### Sources
#####
eq_tends(pv::PV, ::AtmosLinearModel, ::Source) where {PV} = ()

eq_tends(pv::Momentum, ::AtmosAcousticGravityLinearModel, ::Source) =
    (Gravity(),)


text:
using Documenter,Analytical

makedocs(
	modules   = [Analytical],
	format    = Documenter.HTML(),
	sitename  = "ABC-MK",
	authors   = "Jesús Murga-Moreno, Lawrence Uricchio, David Enard",
	pages     = [
		"Home" => "index.md",
		"Package Overview" => [
			"Analytical estimations" => "analytical.md",
			"Processing data" => "data.md",
			"MK approaches" => "mk.md"
		],
		"Infering the rate and strength of adaptation" =>[
			"Empirical estimation" => "empirical.md",
			"Rates" => "rates.md",
			"Input data" => "input.md",
			"Summary statistics" => "summstat.md",
			"ABC inference" => "abc.md",
			"Multiple dataset" => "multiple.md",
		],
		"Command-Line interface" => "cli.md",
		"Reference" => "reference.md"
	]
)

deploydocs(
	repo      = "github.com/jmurga/Analytical.jl.git"
)


text:
# Note: we define type aliases outside of the @muladd block to avoid Revise breaking when code
# inside the @muladd block is edited. See https://github.com/trixi-framework/Trixi.jl/issues/801
# for more details.

# `DGMulti` refers to both multiple DG types (polynomial/SBP, simplices/quads/hexes) as well as
# the use of multi-dimensional operators in the solver.
const DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} =
  DG{<:RefElemData{NDIMS, ElemType, ApproxType}, Mortar, SurfaceIntegral, VolumeIntegral} where {Mortar}

# Type aliases. The first parameter is `ApproxType` since it is more commonly used for dispatch.
const DGMultiWeakForm{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralWeakForm} where {NDIMS}

const DGMultiFluxDiff{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralFluxDifferencing} where {NDIMS}

const DGMultiFluxDiffSBP{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralFluxDifferencing} where {NDIMS, ApproxType<:Union{SBP, AbstractDerivativeOperator}}

const DGMultiSBP{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} where {NDIMS, ElemType, ApproxType<:Union{SBP, AbstractDerivativeOperator}, SurfaceIntegral, VolumeIntegral}

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# these are necessary for pretty printing
polydeg(dg::DGMulti) = dg.basis.N
Base.summary(io::IO, dg::DG) where {DG <: DGMulti} = print(io, "DGMulti(polydeg=$(polydeg(dg)))")
Base.real(rd::RefElemData{NDIMS, Elem, ApproxType, Nfaces, RealT}) where {NDIMS, Elem, ApproxType, Nfaces, RealT} = RealT

"""
    DGMulti(; polydeg::Integer,
              element_type::AbstractElemShape,
              approximation_type=Polynomial(),
              surface_flux=flux_central,
              surface_integral=SurfaceIntegralWeakForm(surface_flux),
              volume_integral=VolumeIntegralWeakForm(),
              RefElemData_kwargs...)

Create a discontinuous Galerkin method which uses
- approximations of polynomial degree `polydeg`
- element type `element_type` (`Tri()`, `Quad()`, `Tet()`, and `Hex()` currently supported)

Optional:
- `approximation_type` (default is `Polynomial()`; `SBP()` also supported for `Tri()`, `Quad()`,
  and `Hex()` element types).
- `RefElemData_kwargs` are additional keyword arguments for `RefElemData`, such as `quad_rule_vol`.
  For more info, see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/).
"""
function DGMulti(; polydeg=nothing,
                   element_type::AbstractElemShape,
                   approximation_type=Polynomial(),
                   surface_flux=flux_central,
                   surface_integral=SurfaceIntegralWeakForm(surface_flux),
                   volume_integral=VolumeIntegralWeakForm(),
                   kwargs...)

  # call dispatchable constructor
  DGMulti(element_type, approximation_type, volume_integral, surface_integral;
          polydeg=polydeg, kwargs...)
end

# dispatchable constructor for DGMulti to allow for specialization
function DGMulti(element_type::AbstractElemShape,
                 approximation_type,
                 volume_integral,
                 surface_integral;
                 polydeg::Integer,
                 kwargs...)

  rd = RefElemData(element_type, approximation_type, polydeg; kwargs...)
  return DG(rd, nothing #= mortar =#, surface_integral, volume_integral)
end

# now that DGMulti is defined, we can define constructors for VertexMappedMesh which use dg::DGMulti
"""
    VertexMappedMesh(vertex_coordinates, EToV, dg::DGMulti;
                     is_on_boundary = nothing,
                     is_periodic::NTuple{NDIMS, Bool} = ntuple(_->false, NDIMS)) where {NDIMS, Tv}

Constructor which uses `dg::DGMulti` instead of `rd::RefElemData`.
"""
VertexMappedMesh(vertex_coordinates, EToV, dg::DGMulti; kwargs...) =
  VertexMappedMesh(vertex_coordinates, EToV, dg.basis; kwargs...)

"""
    VertexMappedMesh(triangulateIO, dg::DGMulti, boundary_dict::Dict{Symbol, Int})

Constructor which uses `dg::DGMulti` instead of `rd::RefElemData`.
"""
VertexMappedMesh(triangulateIO, dg::DGMulti, boundary_dict::Dict{Symbol, Int}) =
  VertexMappedMesh(triangulateIO, dg.basis, boundary_dict)

# Todo: DGMulti. Add traits for dispatch on affine/curved meshes here.

# Matrix type for lazy construction of physical differentiation matrices
# Constructs a lazy linear combination of B = ∑_i coeffs[i] * A[i]
struct LazyMatrixLinearCombo{Tcoeffs, N, Tv, TA <: AbstractMatrix{Tv}} <: AbstractMatrix{Tv}
  matrices::NTuple{N, TA}
  coeffs::NTuple{N, Tcoeffs}
  function LazyMatrixLinearCombo(matrices, coeffs)
    @assert all(matrix -> size(matrix) == size(first(matrices)), matrices)
    new{typeof(first(coeffs)), length(matrices), eltype(first(matrices)), typeof(first(matrices))}(matrices, coeffs)
  end
end
Base.eltype(A::LazyMatrixLinearCombo) = eltype(first(A.matrices))
Base.IndexStyle(A::LazyMatrixLinearCombo) = IndexCartesian()
Base.size(A::LazyMatrixLinearCombo) = size(first(A.matrices))

@inline function Base.getindex(A::LazyMatrixLinearCombo{<:Real, N}, i, j) where {N}
  val = zero(eltype(A))
  for k in Base.OneTo(N)
    val = val + A.coeffs[k] * getindex(A.matrices[k], i, j)
  end
  return val
end

# `SimpleKronecker` lazily stores a Kronecker product `kron(ntuple(A, NDIMS)...)`.
# This object also allocates some temporary storage to enable the fast computation
# of matrix-vector products.
struct SimpleKronecker{NDIMS, TA, Ttmp}
  A::TA
  tmp_storage::Ttmp # temporary array used for Kronecker multiplication
end

# constructor for SimpleKronecker which requires specifying only `NDIMS` and
# the 1D matrix `A`.
function SimpleKronecker(NDIMS, A, eltype_A=eltype(A))
  @assert size(A, 1) == size(A, 2) # check if square
  tmp_storage=[zeros(eltype_A, ntuple(_ -> size(A, 2), NDIMS)...) for _ in 1:Threads.nthreads()]
  return SimpleKronecker{NDIMS, typeof(A), typeof(tmp_storage)}(A, tmp_storage)
end

# Computes `b = kron(A, A) * x` in an optimized fashion
function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{2}, x_in)

  @unpack A = A_kronecker
  tmp_storage = A_kronecker.tmp_storage[Threads.threadid()]
  n = size(A, 2)

  # copy `x_in` to `tmp_storage` to avoid mutating the input
  @assert length(tmp_storage) == length(x_in)
  for i in eachindex(tmp_storage)
    tmp_storage[i] = x_in[i]
  end
  x = reshape(tmp_storage, n, n)
  b = reshape(b_in, n, n)

  @turbo for j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for ii in 1:n
      tmp = tmp + A[i, ii] * x[ii, j]
    end
    b[i, j] = tmp
  end

  @turbo for j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for jj in 1:n
      tmp = tmp + A[j, jj] * b[i, jj]
    end
    x[i, j] = tmp
  end

  @turbo for i in eachindex(b_in)
    b_in[i] = x[i]
  end

  return nothing
end

# Computes `b = kron(A, A, A) * x` in an optimized fashion
function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{3}, x_in)

  @unpack A = A_kronecker
  tmp_storage = A_kronecker.tmp_storage[Threads.threadid()]
  n = size(A, 2)

  # copy `x_in` to `tmp_storage` to avoid mutating the input
  for i in eachindex(tmp_storage)
    tmp_storage[i] = x_in[i]
  end
  x = reshape(tmp_storage, n, n, n)
  b = reshape(b_in, n, n, n)

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for ii in 1:n
      tmp = tmp + A[i, ii] * x[ii, j, k]
    end
    b[i, j, k] = tmp
  end

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for jj in 1:n
      tmp = tmp + A[j, jj] * b[i, jj, k]
    end
    x[i, j, k] = tmp
  end

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for kk in 1:n
      tmp = tmp + A[k, kk] * x[i, j, kk]
    end
    b[i, j, k] = tmp
  end

  return nothing
end

end # @muladd


text:
using Flux
using Flux: onehot, chunk, batchseq, throttle, crossentropy
using StatsBase: wsample
using Base.Iterators: partition

cd(@__DIR__)

isfile("input.txt") ||
  download("https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt",
           "input.txt")

text = collect(String(read("input.txt")))
alphabet = [unique(text)..., '_']
text = map(ch -> onehot(ch, alphabet), text)
stop = onehot('_', alphabet)

N = length(alphabet)
seqlen = 50
nbatch = 50

Xs = collect(partition(batchseq(chunk(text, nbatch), stop), seqlen))
Ys = collect(partition(batchseq(chunk(text[2:end], nbatch), stop), seqlen))

m = Chain(
  LSTM(N, 128),
  LSTM(128, 128),
  Dense(128, N),
  softmax)

m = gpu(m)

function loss(xs, ys)
  l = sum(crossentropy.(m.(gpu.(xs)), gpu.(ys)))
  Flux.truncate!(m)
  return l
end

opt = ADAM(0.01)
tx, ty = (Xs[5], Ys[5])
evalcb = () -> @show loss(tx, ty)

Flux.train!(loss, params(m), zip(Xs, Ys), opt,
            cb = throttle(evalcb, 30))

# Sampling

function sample(m, alphabet, len; temp = 1)
  m = cpu(m)
  Flux.reset!(m)
  buf = IOBuffer()
  c = rand(alphabet)
  for i = 1:len
    write(buf, c)
    c = wsample(alphabet, m(onehot(c, alphabet)).data)
  end
  return String(take!(buf))
end

sample(m, alphabet, 1000) |> println

# evalcb = function ()
#   @show loss(Xs[5], Ys[5])
#   println(sample(deepcopy(m), alphabet, 100))
# end

max_stars_repo_name: shashikdm/model-zoo

max_issues_repo_name: shashikdm/model-zoo

max_forks_repo_name: shashikdm/model-zoo


text:
function plotFunc_timecourse(
        Sim::Module,
        n_file::Vector{Int},
        viz_type::String,
        show_all::Bool,
        stdev::Bool,
        simulations_all::Array{Float64,4})
    if !isdir(strip(MODEL_PATH, '/') * "/figure/simulation/$viz_type")
        mkpath(strip(MODEL_PATH, '/') * "/figure/simulation/$viz_type")
    end

    Viz.set_rcParams()

    for (i,name) in enumerate(observables)
        gca().spines["right"].set_visible(false)
        gca().spines["top"].set_visible(false)
        gca().yaxis.set_ticks_position("left")
        gca().xaxis.set_ticks_position("bottom")
        if viz_type != "experiment"
            if show_all
                for j in eachindex(n_file)
                    for (l,condition) in enumerate(Sim.conditions)
                        if !(condition in Viz.options[i]["dont_show"])
                            plot(
                                Sim.t ./ Viz.options[i]["divided_by"],
                                simulations_all[i,j,:,l] ./ ifelse(
                                    Sim.normalization,
                                    maximum(simulations_all[i,j,:,:]),1.0),
                                color=Viz.options[i]["cmap"][l],
                                lw=0.5,alpha=0.35
                            )
                        end
                    end
                end
            end
            if viz_type == "average"
                normalized = Array{Float64,4}(
                    undef,length(observables),length(n_file),length(Sim.t),length(Sim.conditions)
                )
                @inbounds for j in eachindex(n_file)
                    @simd for l in eachindex(Sim.conditions)
                        normalized[i,j,:,l] = (
                            simulations_all[i,j,:,l] ./ ifelse(
                                Sim.normalization,maximum(simulations_all[i,j,:,:]),1.0
                            )
                        )
                    end
                end
                norm_max::Float64 = maximum(
                    vcat(
                        [
                            [
                                mean(
                                    filter(
                                        !isnan,normalized[i,:,k,l]
                                    )
                                ) for k in eachindex(Sim.t)
                            ] for l in eachindex(Sim.conditions)
                        ]...
                    )
                )
                for j in eachindex(n_file)
                    for k in eachindex(Sim.t)
                        for l in eachindex(Sim.conditions)
                            @inbounds normalized[i,j,k,l] /= ifelse(Sim.normalization,norm_max,1.0)
                        end
                    end
                end
                for (l,condition) in enumerate(Sim.conditions)
                    if !(condition in Viz.options[i]["dont_show"])
                        plot(
                            Sim.t ./ Viz.options[i]["divided_by"],[
                                mean(
                                    filter(
                                        !isnan,normalized[i,:,k,l]
                                    )
                                ) for k in eachindex(Sim.t)
                            ],color=Viz.options[i]["cmap"][l],
                            label=condition
                        )
                    end
                end
                if stdev
                    for (l,condition) in enumerate(Sim.conditions)
                        if !(condition in Viz.options[i]["dont_show"])
                            y_mean = [
                                mean(
                                    filter(
                                        !isnan,normalized[i,:,k,l]
                                    )
                                ) for k in eachindex(Sim.t)
                            ]
                            y_std = [
                                std(
                                    filter(
                                        !isnan,normalized[i,:,k,l]
                                    )
                                ) for k in eachindex(Sim.t)
                            ]
                            fill_between(
                                Sim.t ./ Viz.options[i]["divided_by"],
                                y_mean-y_std,y_mean+y_std,
                                lw=0,color=Viz.options[i]["cmap"][l],alpha=0.1
                            )
                        end
                    end
                end
            else
                for (l,condition) in enumerate(Sim.conditions)
                    if !(condition in Viz.options[i]["dont_show"])
                        plot(
                            Sim.t ./ Viz.options[i]["divided_by"],
                            Sim.simulations[i,:,l] / ifelse(
                                Sim.normalization,
                                maximum(Sim.simulations[i,:,:]),1.0),
                            color=Viz.options[i]["cmap"][l],
                            label=condition
                        )
                    end
                end
            end
        end

        if Viz.options[i]["exp_data"] && isassigned(Exp.experiments,i)
            exp_t = Exp.get_timepoint(name)
            if isassigned(Exp.error_bars,i)
                for (l,condition) in enumerate(Sim.conditions)
                    if condition in keys(Exp.experiments[i]) && !(condition in Viz.options[i]["dont_show"])
                        if viz_type == "experiment"
                            Viz.options[i]["shape"][l] *= "-"
                        end
                        exp_data = errorbar(
                            exp_t ./ Viz.options[i]["divided_by"],
                            Exp.experiments[i][condition],
                            yerr=Exp.error_bars[i][condition],
                            lw=1,markerfacecolor="None",
                            color=Viz.options[i]["cmap"][l],
                            markeredgecolor=Viz.options[i]["cmap"][l],
                            ecolor=Viz.options[i]["cmap"][l],
                            fmt=Viz.options[i]["shape"][l],capsize=8,
                            clip_on=false
                        )
                        for capline in exp_data[2]
                            capline.set_clip_on(false)
                        end
                        for barlinecol in exp_data[3]
                            barlinecol.set_clip_on(false)
                        end
                    end
                end
            else
                for (l,condition) in enumerate(Sim.conditions)
                    if condition in keys(Exp.experiments[i]) && !(condition in Viz.options[i]["dont_show"])
                        if viz_type == "experiment"
                            Viz.options[i]["shape"][l] *= "-"
                        end
                        plot(
                            exp_t ./ Viz.options[i]["divided_by"],
                            Exp.experiments[i][condition],
                            Viz.options[i]["shape"][l],
                            color=Viz.options[i]["cmap"][l],
                            markerfacecolor="None",
                            markeredgecolor=Viz.options[i]["cmap"][l],
                            clip_on=false
                        )
                    end
                end
            end
        end

        if length(Viz.options[i]["xlim"]) > 0
            xlim(Viz.options[i]["xlim"]...)
        end
        if Viz.options[i]["xticks"] !== nothing
            xticks(Viz.options[i]["xticks"])
        end
        xlabel(Viz.options[i]["xlabel"])
        if length(Viz.options[i]["ylim"]) > 0
            ylim(Viz.options[i]["ylim"]...)
        end
        if Viz.options[i]["yticks"] !== nothing
            yticks(Viz.options[i]["yticks"])
        end
        ylabel(Viz.options[i]["ylabel"])
        if Viz.options[i]["legend_loc"] !== nothing
            legend(
                loc=Viz.options[i]["legend_loc"], frameon=false, fontsize=15
            )
        end
        savefig(
            strip(MODEL_PATH, '/') * "/figure/simulation/$viz_type/$name.pdf", bbox_inches="tight"
        )
        close()
    end
end

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim