text:
mutable struct Validation <: IFeature
    mVkInstance::vk.VkInstance

    mPaused::Bool
    mMessageCount::UInt32
    mCallback::vk.VkDebugReportCallbackEXT

    function Validation()
        this = new()
        this.mPaused = false
        this.mMessageCount = 0
        return this
    end
end

@class Validation

function LavaCore.:destroy(this::Validation)
    if this.mCallback != 0
        VkExt.destroyDebugReportCallbackEXT(this.mVkInstance, this.mCallback)
        this.mCallback = vk.VkDebugReportCallbackEXT(0)
    end
end

function create(::Type{Validation})
    return Validation()
end

function paused(this::Validation)
    return this.mPaused
end

function bumpMessageCount(this::Validation)
    this.mMessageCount += 1
end

function debugCallback(flags::vk.VkDebugReportFlagsEXT,
                     objType::vk.VkDebugReportObjectTypeEXT,
                         obj,
                    location::UInt64,
                        code::Int32,
                 layerPrefix,
                         msg,
                    userData)

    validationPtr = Base.unsafe_convert(Ptr{Validation}, userData)
    validation = Base.unsafe_load(validationPtr)
    if (paused(validation))
        return VkExt.VK_FALSE
    end

    bumpMessageCount(validation)
    println("--------", validation.mMessageCount , ".Debug Report--------")
    println("flags:", flags)
    println("objType:", objType)
    println("obj:", obj)
    println("location:", location)
    println("code:", code)
    println("layerPrefix:", unsafe_string(layerPrefix))
    println("msg:", unsafe_string(msg))

    return VkExt.VK_FALSE
end

function beforeInstanceDestruction(this::Validation)
    println("Before instance destruct in Validation: ", this.mCallback)
    if this.mCallback != 0
        VkExt.destroyDebugReportCallbackEXT(this.mVkInstance, this.mCallback)
        this.mCallback = vk.VkDebugReportCallbackEXT(0)
    end
end

########## override begin ##########
function LavaCore.:layers(this::Validation, available::Vector{String})::Vector{String}
    return ["VK_LAYER_LUNARG_standard_validation"]
end

function LavaCore.:instanceExtensions(this::Validation, available::Vector{String})::Vector{String}
    return [vk.VK_EXT_DEBUG_REPORT_EXTENSION_NAME]
end

function LavaCore.:onInstanceCreated(this::Validation, vkInstance::vk.VkInstance)
    this.mVkInstance = vkInstance
    pfnCallback = @cfunction(debugCallback, vk.VkBool32, (vk.VkDebugReportFlagsEXT, vk.VkDebugReportObjectTypeEXT, Culonglong, Csize_t, Cint, Cstring, Cstring, Ptr{Cvoid}))
    validationRef = Ref(this)
    debug = vk.VkDebugReportCallbackCreateInfoEXT(
        vk.VK_STRUCTURE_TYPE_DEBUG_REPORT_CALLBACK_CREATE_INFO_EXT, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        vk.VK_DEBUG_REPORT_ERROR_BIT_EXT | vk.VK_DEBUG_REPORT_WARNING_BIT_EXT, #flags::VkDebugReportFlagsEXT
        pfnCallback, #pfnCallback::PFN_vkDebugReportCallbackEXT
        Base.unsafe_convert(Ptr{Cvoid}, validationRef) #pUserData::Ptr{Cvoid}
    )
    this.mCallback = VkExt.createDebugReportCallbackEXT(vkInstance, debug);
end
########## override end ##########


text:
using Distributions, PDMats
using Test, LinearAlgebra
import Distributions: ispossemdef

# RealInterval
r = RealInterval(1.5, 4.0)
@test minimum(r) == 1.5
@test maximum(r) == 4.0
@test extrema(r) == (1.5, 4.0)

@test partype(Gamma(1, 2)) == Float64
@test partype(Gamma(1.1, 2)) == Float64
@test partype(Normal(1//1, 2//1)) == Rational{Int}
@test partype(MvNormal(rand(Float32, 5), Matrix{Float32}(I, 5, 5))) == Float32

# special cases
@test partype(Kolmogorov()) == Float64
@test partype(Hypergeometric(2, 2, 2)) == Float64
@test partype(DiscreteUniform(0, 4)) == Float64

A = rand(1:10, 5, 5)
B = rand(Float32, 4)
C = 1//2
L = rand(Float32, 4, 4)
D = PDMats.PDMat(L * L')

# Ensure that utilities functions works with abstract arrays

@test isprobvec(GenericArray([1, 1, 1])) == false
@test isprobvec(GenericArray([1/3, 1/3, 1/3]))

# Positive definite matrix
M = GenericArray([1.0 0.0; 0.0 1.0])
# Non-invertible matrix
N = GenericArray([1.0 0.0; 1.0 0.0])

@test Distributions.isApproxSymmmetric(N) == false
@test Distributions.isApproxSymmmetric(M)

n = 10
areal = randn(n,n)/2
aimg  = randn(n,n)/2
@testset "For A containing $eltya" for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int)
    ainit = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    @testset "Positive semi-definiteness" begin
        notsymmetric = ainit
        notsquare    = [ainit ainit]
        @test !ispossemdef(notsymmetric)
        @test !ispossemdef(notsquare)
        for truerank in 0:n
            X = ainit[:, 1:truerank]
            A = truerank == 0 ? zeros(eltya, n, n) : X * X'
            @test ispossemdef(A)
            for testrank in 0:n
                if testrank == truerank
                    @test ispossemdef(A, testrank)
                else
                    @test !ispossemdef(A, testrank)
                end
            end
            @test !ispossemdef(notsymmetric, truerank)
            @test !ispossemdef(notsquare, truerank)
            @test_throws ArgumentError ispossemdef(A, -1)
            @test_throws ArgumentError ispossemdef(A, n + 1)
        end
    end
end

ext: jl


text:
#using TensorToolbox, Test, LinearAlgebra

println("\n\n**** Testing tensor.jl")

X=rand(20,10,50,5)
N=ndims(X)
println("\n**Test tensor X of size: ", size(X))

println("\n...Testing functions matten and tenmat (by mode).")
for n=1:N
  Xn=tenmat(X,n);
  M=setdiff(1:N,n);
  println("Size of $n-mode matricization: ", [size(Xn)...])
  @test size(Xn) == (size(X,n),prod([size(X,M[k]) for k=1:length(M)]))
  println("Check if it folds back correctly: ",matten(Xn,n,[size(X)...]) == X)
  @test matten(Xn,n,[size(X)...]) == X
end

println("\n...Testing functions matten and tenmat (by rows and columns).")
R=[2,1];C=[4,3];
Xmat=tenmat(X,row=R,col=C);
println("Size of R=$R and C=$C matricization: ", [size(Xmat)...])
println("Check if it folds back correctly: ",matten(Xmat,R,C,[size(X)...]) == X)
@test matten(Xmat,R,C,[size(X)...]) == X

println("\n...Testing function ttm.")
M=MatrixCell(undef,N)
for n=1:N
  M[n]=rand(5,size(X,n))
end
println("Created $N matrices with 5 rows and appropriate number of columns.")
Xprod=ttm(X,M)
println("Size of tensor Y=ttm(X,M): ",size(Xprod))
err=norm(tenmat(Xprod,1) - M[1]*tenmat(X,1)*kron(M[end:-1:2])')
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing function ttv.")
Xk=reshape(collect(1:24),(3,4,2))
mode=2
v=collect(1:4)
println("Multiplying a tensor X by a vector v in mode $mode.")
Xprod=ttv(Xk,v,mode)
println("Size of tensor Y=ttv(X,v): ",size(Xprod))
res=[70 190;80 200;90 210]
@test Xprod==res

println("\n...Testing function krontm.")
X=rand(5,4,3)
Y=rand(2,5,4)
println("Created two tensors X and Y of order ",ndims(X)," and sizes ",size(X)," and ",size(Y),".")
mode=3
M1=rand(20,10)
M2=rand(20,20)
M3=rand(20,12)
println("Multiplying tkron(X,Y) by random matrix in mode $mode.")
Z=krontm(X,Y,M3,mode)
err= norm(Z-ttm(tkron(X,Y),M3,mode))
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10
mode=[3,2]
M=[M3,M2]
println("Multiplying tkron(X,Y) by random matrices in modes $mode.")
Z=krontm(X,Y,M,mode)
err = norm(Z-ttm(tkron(X,Y),M,mode))
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10
M=[M1,M2,M3]
println("Multiplying tkron(X,Y) by random matrices in all modes.")
Z=krontm(X,Y,M)
err = norm(Z-ttm(tkron(X,Y),M))
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing function mkrontv.")
v=rand(240)
mode=1
println("Multiplying mode-$mode matricized tkron(X,Y) by a random vector.")
Z=mkrontv(X,Y,v,mode)
err = norm(Z-tenmat(tkron(X,Y),mode)*v)
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10
v=rand(10)
Z=mkrontv(X,Y,v,mode,'t')
err = norm(Z-tenmat(tkron(X,Y),mode)'*v)
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing function mttkrp.")
X=rand(5,4,3)
mode=1
A1=rand(2,5);
A2=rand(4,5);
A3=rand(3,5);
A=[A1,A2,A3]
println("Multiplying mode-$mode matricized tensor X by Khatri-Rao product of matrices.")
Z=mttkrp(X,A,mode)
err = norm(Z-tenmat(X,mode)*khatrirao(A3,A2))
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing function dropdims.")
X=rand(5,4,1,3,6,1)
Xsq=dropdims(X)
println("Tensor X of size :",size(X)," squeezed to size :",size(Xsq),".")

println("\n...Testing contraction (contract).")
X=rand(5,4,3);
Y=rand(3,4,2);
Z=contract(X,Y)
for i1=1:5, i2=1:4, j2=1:4,j3=1:2
    s=0;
    for i=1:3
        s+=X[i1,i2,i]*Y[i,j2,j3];
    end
    @test Z[i1,i2,j2,j3] ≈ s atol=1e-10
end
X=rand(5,4,3,2);
Y=rand(5,7,3,6);
Z=contract(X,[1,3],Y,[1,3])
err=norm(Z-reshape(tenmat(X,row=[1,3])'*tenmat(Y,row=[1,3]),(4,2,7,6)))
println("Error: ",err)
@test err ≈ 0 atol=1e-10


text:
## ExoplanetsSysSim/src/mission_constants.jl
## (c) 2015 Eric B. Ford

#module KeplerConstants
 const num_channels = 84
 const num_modules = 42
 const num_quarters = 17              # QUERY:  I'm favoring leaving out quarter 0, since that was engineering data.  Agree?
 const cdpp_durations = [1.5,2.,2.5,3.,3.5,4.5,5.,6.,7.5,9.,10.5,12.,12.5,15.]
 const duration_symbols = [:rrmscdpp01p5, :rrmscdpp02p0,:rrmscdpp02p5,:rrmscdpp03p0,:rrmscdpp03p5,:rrmscdpp04p5,:rrmscdpp05p0,:rrmscdpp06p0,:rrmscdpp07p5,:rrmscdpp09p0,:rrmscdpp10p5,:rrmscdpp12p0,:rrmscdpp12p5,:rrmscdpp15p0 ]
 const num_cdpp_timescales = 14
 @assert num_cdpp_timescales == length(cdpp_durations) == length(duration_symbols)
 const mission_data_span = 1459.789   # maximum(ExoplanetsSysSim.StellarTable.df[:dataspan])
 const mission_duty_cycle = 0.8751    # median(ExoplanetsSysSim.StellarTable.df[:dutycycle])

 const kepler_exp_time_internal  =  6.019802903/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const kepler_read_time_internal = 0.5189485261/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const num_exposures_per_kepler_LC = 270
 const num_exposures_per_kepler_SC = 9
 const kepler_LC_integration_time = kepler_exp_time_internal*num_exposures_per_kepler_LC
 const kepler_SC_integration_time = kepler_exp_time_internal*num_exposures_per_kepler_SC
 const kepler_LC_read_time = kepler_read_time_internal*num_exposures_per_kepler_LC
 const kepler_SC_read_time = kepler_read_time_internal*num_exposures_per_kepler_SC
 const kepler_LC_duration = kepler_LC_integration_time +  kepler_LC_read_time
 const kepler_SC_duration = kepler_SC_integration_time +  kepler_SC_read_time
 const kepler_LC_rate = 1.0/kepler_LC_duration
 const kepler_SC_rate = 1.0/kepler_SC_duration

 const tess_sector_duration = 27.4 # days
 const tess_exp_time_internal = 96.0 / (24*60*60) # https://archive.stsci.edu/files/live/sites/mast/files/home/missions-and-data/active-missions/tess/_documents/TESS_Instrument_Handbook_v0.1.pdf
 const tess_read_time_internal = 24.0 / (24*60*60)
 const num_exposures_per_tess_LC = 15
 const num_exposures_per_tess_SC = 1
 const tess_LC_integration_time = tess_exp_time_internal * num_exposures_per_tess_LC
 const tess_SC_integration_time = tess_exp_time_internal * num_exposures_per_tess_SC
 const tess_LC_read_time = tess_read_time_internal*num_exposures_per_tess_LC
 const tess_SC_read_time = tess_read_time_internal*num_exposures_per_tess_SC
 const tess_LC_duration = tess_LC_integration_time +  tess_LC_read_time
 const tess_SC_duration = tess_SC_integration_time +  tess_SC_read_time
 const tess_LC_rate = 1.0/tess_LC_duration
 const tess_SC_rate = 1.0/tess_SC_duration

# export num_channels,num_modules,num_quarters,num_cdpp_timescales,mission_data_span,mission_duty_cycle
# export kepler_exp_time_internal,kepler_read_time_internal,num_exposures_per_kepler_LC,num_exposures_per_kepler_SC
# export kepler_LC_integration_time,kepler_SC_integration_time,kepler_LC_read_time,kepler_SC_read_time,kepler_LC_duration,kepler_SC_duration,kepler_LC_rate,kepler_SC_rate
#end

#using KeplerConstants

# Standard conversion factors on which unit system is based
const global AU_in_m_IAU2012 = 149597870700.0
 const global G_in_mks_IAU2015 = 6.67384e-11
 const global G_mass_sun_in_mks = 1.3271244e20
 const global G_mass_earth_in_mks = 3.986004e14
 const global sun_radius_in_m_IAU2015 = 6.9566e8
 const global earth_radius_eq_in_m_IAU2015 = 6.3781e6
 const global sun_mass_in_kg_IAU2010 = 1.988547e30

# Constants used by this code
const global sun_mass = 1.0
 const global earth_mass = G_mass_earth_in_mks/G_mass_sun_in_mks  # about 3.0024584e-6
 const global earth_radius = earth_radius_eq_in_m_IAU2015 / sun_radius_in_m_IAU2015 # about 0.0091705248
 const global rsol_in_au = sun_radius_in_m_IAU2015 / AU_in_m_IAU2012  # about 0.00464913034
 const global sec_in_day = 24*60*60
 const global grav_const = G_in_mks_IAU2015 * sec_in_day^2 * sun_mass_in_kg_IAU2010 / AU_in_m_IAU2012^3 # about 2.9591220363e-4 in AU^3/(day^2 Msol)
 const global day_in_year = 365.25


text:
# helper scripts for event handlung
"""
copy field-info from source to closest target

`match_fun` can be "closest", "s before t" (eq. "t after s") or "t before s" (eq. "s after t")

## Example
julia> # copy reaction time values from button press to closest stimulus immediately before button press
julia> copy_eventinfo!(evts,"button"=>"stimulus",:reaction_time;match_fun="s after t")
"""
function copy_eventinfo!(evts,fromTo,field;match_fun="closest")

    source = fromTo.first
    target = fromTo.second

    if isa(field,Pair)
        source_field = field.first
        target_field = field.second
    else
        source_field = field
        target_field = field
    end
      source_ix = findall(isequal(source),evts.trial_type)
      target_ix = findall(isequal(target),evts.trial_type)

    isempty(source_ix) && error("couldnt find source entries ($source) in evts.trial_type" )
    isempty(target_ix) && error("couldnt find target entries ($target) in evts.trial_type" )

    # for some matching functions we want to find the minimum, but no negative number
    filter_greaterzero = x-> x >= 0 ? x : Inf

    if match_fun == "closest"
      match_fun = (a,b)-> argmin(abs.(a .- b)) # find closest before or after
    elseif match_fun == "s after t" || match_fun == "t before s"
      match_fun = (s,t)-> argmin(filter_greaterzero.(s .- t)) # source after target
    elseif match_fun == "s before t" || match_fun == "t after s"
      match_fun = (s,t)-> argmin(filter_greaterzero.(t .- s)) # source before target
    elseif !callable(match_fun)
        error("bad match_fun input")
    end
      ix = [match_fun(evts[s,:latency], evts[target_ix,:latency]) for s in source_ix]
    ix = target_ix[ix]
    payload = evts[source_ix,source_field]

    # if targetfield does not exist, create one with Union(missing,target_field_datatype)
    if !any(names(evts) .== target_field)
        evts[!,target_field] = Array{Union{Missing,typeof(payload[1])}}(missing,nrow(evts))
    end

    # fill it in
    evts[ix,target_field] .= payload
    return evts
end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
# Timing tests of transit_poly_g!:

include("transit_poly_struct.jl")
#trans = transit_init(0.1,0.5,[0.5,0.5],true)
#u_n = [0.5,0.25,0.125,0.0625,0.03125]
u_n = [0.5,0.25]
#trans = transit_init(0.1,0.5,u_n,true)
#n = 6
#u_n = 2.^(-linspace(1,n,n))
#trans = transit_init(0.1,0.5,[0.5,0.25,0.125,0.0625,0.03125],false)
trans = transit_init(0.1,0.5,u_n,true)
#Transit_Struct{Float64}(0.1, 0.5, [0.5, 0.5], 2, 3, [-0.25, 1.5, -0.125], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], false, [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0 0.0; 0.0 0.0; 0.0 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0])

function transit_b!(trans::Transit_Struct{T},b::Array{T,1},flux::Array{T,1}) where {T <: Real}
#  flux = zeros(length(b))
#  gradient = zeros(length(b),2+trans.n)
  for i=1:length(b)
    trans.b = b[i]
#    flux[i]=transit_poly_g!(trans)
#    gradient[i,1:2]=trans.dfdrb
#    gradient[i,3:2+trans.n]=trans.dfdu
    if trans.grad
      flux[i]=transit_poly_g!(trans)
    else
      flux[i]=transit_poly_g(trans)
    end
  end
#  return flux,gradient
  return
end

nb = 10000000
b = abs.(linspace(-1.2,1.2,nb))
flux = zeros(nb)
@time transit_b!(trans,b,flux);
@time transit_b!(trans,b,flux);

#=
Profile.clear()
Profile.init(n=10^7,delay=0.01)
@profile transit_b(trans,b)
Profile.print()
=#


text:
"""
Abstract type for [derivative contract](https://en.wikipedia.org/wiki/Derivative_(finance)).
"""
abstract type Derivative end

"""
Abstract type for [forward contract](https://en.wikipedia.org/wiki/Forward_contract).
"""
abstract type Forward <: Derivative end

"""
Abstract type for [futures contract](https://en.wikipedia.org/wiki/Futures_contract).
"""
abstract type Future <: Derivative end

"""
Abstract type for [option](https://en.wikipedia.org/wiki/Option_(finance)).
"""
abstract type Option <: Derivative end

"""
Abstract type for [swap](https://en.wikipedia.org/wiki/Swap_(finance)).
"""
abstract type Swap <: Derivative end

"""
Abstract type for [swaption](https://en.wikipedia.org/wiki/Swaption).
"""
abstract type Swaption <: Derivative end

"""
    InterestRateDerivative(k, r, σ, θ, t)

[Interest rate derivative](https://en.wikipedia.org/wiki/Interest_rate_derivative).

# Arguments
- `k`: speed of reversion
- `r`: initial interest rate
- `σ`: instantaneous volatility
- `θ`: long term mean level
- `t`: time interval
"""
struct InterestRateDerivative{T<:Number}
    k::T
    r::T
    σ::T
    θ::T
    t::T
end

"""
    AmericanOption(s, k, r, σ, t, call)

[American option](https://en.wikipedia.org/wiki/Option_style#American_and_European_options).

# Arguments
- `s`: underlying price
- `k`: strike price
- `r`: risk-free interest rate
- `σ`: volatility
- `t`: time to expiration
- `call`: 1 if call, -1 if put
"""
struct AmericanOption{T<:Number} <: Option
    s::T
    k::T
    r::T
    σ::T
    t::T
    call::Int64
end

"""
    AsianOption(s, k, r, σ, t, call)

[Asian option](https://en.wikipedia.org/wiki/Option_style#Asian_option).

# Arguments
- `s`: underlying price
- `k`: strike price
- `r`: risk-free interest rate
- `σ`: volatility
- `t`: time to expiration
- `call`: 1 if call, -1 if put
"""
struct AsianOption{T<:Number} <: Option
    s::T
    k::T
    r::T
    σ::T
    t::T
    call::Int64
end

"""
    EuropeanOption(s, k, r, σ, t, call)

[European option](https://en.wikipedia.org/wiki/Option_style#American_and_European_options).

# Arguments
- `s`: underlying price
- `k`: strike price
- `r`: risk-free interest rate
- `σ`: volatility
- `t`: time to expiration
- `call`: 1 if call, -1 if put
"""
struct EuropeanOption{T<:Number} <: Option
    s::T
    k::T
    r::T
    σ::T
    t::T
    call::Int64
end

"""
    FXOption(s, k, r_d, r_f, σ, t, call)

[Foreign exchange option](https://en.wikipedia.org/wiki/Foreign_exchange_option).

# Arguments
- `s`: underlying price
- `k`: strike price
- `r_d`: domestic risk-free interest rate
- `r_f`: foreign risk-free interest rate
- `σ`: volatility
- `t`: time to expiration
- `call`: 1 if call, -1 if put
"""
struct FXOption{T<:Number} <: Option
    s::T
    k::T
    r_d::T
    r_f::T
    σ::T
    t::T
    call::Int64
end


text:
module StenoGraphs

    export AbstractNode
    export Node
    export MetaNode
    export SimpleNode
    export AbstractEdge
    export Edge
    export MetaEdge
    export DirectedEdge
    export UndirectedEdge
    export Modifier
    export EdgeModifier
    export NodeModifier
    export NodeOrEdgeModifier
    export ModifiedEdge
    export ModifyingNode
    export ModifiedNode
    export @StenoGraph
    export *
    export convert
    export promote_rule
    export →, ←, ↔, ⇒, ⇐, ⇔
    export ==
    export unmeta
    export unarrow
    export show

    include("graph_types.jl")
    include("syntax/symbol.jl")
    include("syntax/arrow.jl")
    include("syntax/addition.jl")
    include("broadcast.jl")
    include("macro_communative.jl")
    include("syntax/macro.jl")
    include("modification.jl")
    include("symbol_as_node.jl")
    include("compare.jl")
    include("keep.jl")
    include("show.jl")

end


text:
##################################################
# Search index structures and associated methods #
##################################################
"""
IVFADC index type for storing vectors. It is a wrapper around a
`IVFADCIndex` (inverted file system with asymmetric distance computation)
structure and performs a billion-scale search using a distance-based
similarity between vectors.

# References
 * [Jègou et al. "Product quantization for nearest neighbor search"](https://hal.inria.fr/file/index/docid/514462/filename/paper_hal.pdf)
 * [Baranchuk et al. "Revisiting the inverted indices for billion-scale approximate nearest neighbors"](http://openaccess.thecvf.com/content_ECCV_2018/papers/Dmitry_Baranchuk_Revisiting_the_Inverted_ECCV_2018_paper.pdf)
"""
struct IVFIndex{U,I,Dc,Dr,T,Q} <: AbstractIndex
    index::IVFADCIndex{U,I,Dc,Dr,T,Q}
end

IVFIndex(data, args...; kwargs...) = IVFIndex(IVFADCIndex(densify(data); kwargs...))  # args are ignored

# Nearest neighbor search method
function knn_search(index::IVFIndex{U,I,Dc,Dr,T,Q},
                    point::AbstractVector{T},
                    k::Int,
                    keep::Vector{Int}=collect(1:length(index));
                    w::Int=1
                   ) where {U,I,Dc,Dr,T,Q}
    # Uses Euclidean distance by default
    _idxs, scores = knn_search(index.index, Vector(point), k; w=w)
    idxs = Int.(_idxs) .+ 1
    if length(keep) == length(index)
        # all data points are valid
        return idxs, scores
    else
        # this bit is slow if 'keep' is large
        mask = map(idx->in(idx, keep), idxs)
        return idxs[mask], scores[mask]
    end
end

# Length method
length(index::IVFIndex) = length(index.index)

# push!, pushfirst!, pop!, popfirst!, delete_from_index!
Base.push!(index::IVFIndex, point) = push!(index.index, densify(point))

Base.pushfirst!(index::IVFIndex, point) = pushfirst!(index.index, densify(point))

Base.pop!(index::IVFIndex) = pop!(index.index)

Base.popfirst!(index::IVFIndex) = popfirst!(index.index)

delete_from_index!(index::IVFIndex, pos) = begin
    IVFADC.delete_from_index!(index.index, pos)
    nothing
end


text:
function mergesort(arr::Vector)
    if length(arr) ≤ 1 return arr end
    mid = length(arr) ÷ 2
    lpart = mergesort(arr[1:mid])
    rpart = mergesort(arr[mid+1:end])
    rst = similar(arr)
    i = ri = li = 1
    @inbounds while li ≤ length(lpart) && ri ≤ length(rpart)
        if lpart[li] ≤ rpart[ri]
            rst[i] = lpart[li]
            li += 1
        else
            rst[i] = rpart[ri]
            ri += 1
        end
        i += 1
    end
    if li ≤ length(lpart)
        copy!(rst, i, lpart, li)
    else
        copy!(rst, i, rpart, ri)
    end
    return rst
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", mergesort(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
exports fmt and @fmt, functions to build
format specifiers in an alternate
way to @fmt_str
"""
module Cons
export fmt, @fmt
import Base: @pure
import ..GenPrintf: PrintArgument, PrintValue,
                    PrintBoth, FormatSpecifier

"""
builds format specifiers.
examples:
  @fmt("abc")
  @fmt("hello, ", {})
"""
macro fmt(args...)
  fmt(args...)
end

"""
builds format specifiers.
examples:
  fmt("abc")

  a = "hello, "
  fmt(a, :{})
"""
@pure function fmt(literal::Union{Number, String})
  PrintValue{Symbol(String(literal))}()
end

@pure function fmt(ex::Expr)
  if ex.head == :tuple
    fmt(ex.args...)
  else
    @assert ex == :{} "invalid formatting expression"
    PrintArgument()
  end
end

@pure function fmt(ex::Tuple)
  fmt(args...)
end

@pure function fmt(ex::FormatSpecifier)
  ex
end

@pure function fmt(args...)
  PrintBoth(fmt(args[1]), fmt(args[2:end]...))
end

end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  G'(t) = -(p1 + X(t)) * G(t) + p1 * Gb(t) + v * R(t),
  X'(t) = -p2 * X(t) + p3 * (u(t) - Ib(t)),
  R'(t) = k,
  Ib'(t) = 0,
  Gb'(t) = 0,
  y1(t) = G(t),
  y2(t) = Ib(t),
  y3(t) = Gb(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using Documenter, SymbolicRegression
using SymbolicRegression: Node, PopMember, Population, evalTreeArray, Dataset, HallOfFame, CONST_TYPE, stringTree

makedocs(sitename="SymbolicRegression.jl",
         authors="Miles Cranmer",
         doctest=false, clean=true,
         format=Documenter.HTML(canonical="https://astroautomata.com/SymbolicRegression.jl/stable")
         )

deploydocs(repo="github.com/MilesCranmer/SymbolicRegression.jl.git")


text:
function check_sol(T0, P, phi, p; i_exp=0)
    l_plt = []
    ts, pred = get_idt(T0, P, phi, p; dT=dT)
    ts .+= 1.e-6
    plt = plot(ts, pred[end, :], lw=2, label="Train")
    ylabel!(plt, "Temperature [K]")
    xlabel!(plt, "Time [s]")
    title!(plt, @sprintf("%s, IDT=%.2e [s] \n @%.1f K, %.1f atm, phi=%.1f",
                fuel, ts[end], T0, P / one_atm, phi))
    push!(l_plt, plt)
    for s in [fuel, oxygen]
        plt =
            plot(ts,
            pred[species_index(gas, "$s"), :],
            lw=2, label="Train")
        ylabel!(plt, "Y $s")
        xlabel!(plt, "Time [s]")
        push!(l_plt, plt)
    end
    pltsum = plot(l_plt..., legend=false, framestyle=:box, xscale=:log10)

    png(pltsum, string(fig_path, "/conditions/sol_$i_exp"))
end

function plot_sol(sol)
    plt = Plots.scatter(sol.t, sol[end, :]);
    xlabel!(plt, "Time [s]");
    ylabel!(plt, "T [K]");
    ind_ign = get_ind_ign(sol; dT=dT)
    idt = sol.t[ind_ign]
    xlims!(plt, (0.0, idt * 2.0))
    png(plt, string(fig_path, "/conditions/sol_"));
end

max_stars_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_issues_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_forks_repo_name: DENG-MIT/ArrheniusActiveSubspace


text:
module UnitQuaternions

import Base: +, -, ^, ==, ×, angle, inv, log, mean, show, getindex

export UnitQuaternion, ⊕, ⊞, ⊟, axis, covariance, rotatevector, rotateframe, rotationmatrix, slerp, vector

const EPS = 1e-9

####################################################################################################
# Type Definition
####################################################################################################

immutable UnitQuaternion
    ϵ::Vector{Float64}
    η::Float64

    function UnitQuaternion(ϵ, η)
        length(ϵ) == 3 || error("Must be a 3-vector.")
        η = abs(η) < EPS ? zero(η) : η # small scalars are exactly zero
        # normalization algorithm taken from http://stackoverflow.com/a/12934750/1053656
        squared_mag = η * η + sum(abs2(ϵ))
        if abs(one(squared_mag) - squared_mag) < 2.107342e-08
            scale = 2.0 / (1.0 + squared_mag)
        else
            scale = 1.0 / sqrt(squared_mag)
        end
        η >= zero(η) ? new(ϵ * scale, η * scale) : new(-ϵ * scale, -η * scale)
    end
end

####################################################################################################
# Constructors
####################################################################################################

"""
`UnitQuaternion()`

Constructs the identity unit quaternion.
"""
UnitQuaternion() = UnitQuaternion([0, 0, 0], 1)

"""
`UnitQuaternion(ϵ1, ϵ2, ϵ3, η)`

Constructs a unit quaternion where the vector part is `[ϵ1, ϵ2, ϵ3]` and the scalar part `η`.
"""
UnitQuaternion(ϵ1::Real, ϵ2::Real, ϵ3::Real, η::Real) = UnitQuaternion([ϵ1, ϵ2, ϵ3], η)

"""
`UnitQuaternion(v)`

Constructs a unit quaternion from the vector `v`. If `v` is a 4-vector, its first three elements are
taken as the vector part and its fourth element is taken as the scalar part. If `v` is a 3-vector, it
is assumed to be a rotation vector and the rotation is converted to a unit quaternion.
"""
function UnitQuaternion(v::AbstractVector)
    if length(v) == 3
        θ = norm(v)
        return θ > EPS ? UnitQuaternion(sin(θ/2.0) * v / θ, cos(θ/2.0)) : UnitQuaternion()
    elseif length(v) == 4
        return UnitQuaternion([v[1], v[2], v[3]], v[4])
    else
        error("Must be a 3-vector (rotation vector) or a 4-vector (quaternion).")
    end
end

####################################################################################################
# Operators
####################################################################################################

"""
`+(q)`

Applies the left-hand compound operator to the unit quaternion `q`, which returns a 4x4 matrix.
"""
+(q::UnitQuaternion) = [q.η * eye(3) - ×(q.ϵ) q.ϵ; -q.ϵ' q.η]

"""
`+(p, q)`

Calculates the quaternion product `p + q` using the left-hand compound operator.
"""
+(p::UnitQuaternion, q::UnitQuaternion) = UnitQuaternion(+(p) * vector(q))

"""
`^(q, n)`

Raises the unit quaternion `q` to the power of `n`.
"""
^(q::UnitQuaternion, n::Integer) = UnitQuaternion(2n * log(q))
^(q::UnitQuaternion, n::Real) = UnitQuaternion(2n * log(q))

"""
`⊕(q)`

Applies the right-hand compound operator to the unit quaternion `q`, which returns a 4x4 matrix.
"""
⊕(q::UnitQuaternion) = [q.η * eye(3) + ×(q.ϵ) q.ϵ; -q.ϵ' q.η]

"""
`⊕(p, q)`

Calculates the quaternion product `q + p` using the right-hand compound operator (i.e., `p ⊕ q`).
"""
⊕(p::UnitQuaternion, q::UnitQuaternion) = UnitQuaternion(⊕(p) * vector(q))

"""
`⊟(p, q)`

Calculates the rotation from unit quaternion `p` to unit quaternion `q` (i.e., their difference)
and converts the result to a rotation vector.
"""
⊟(p::UnitQuaternion, q::UnitQuaternion) = 2.0 * log(inv(q) + p)

"""
`⊞(q, rotationvector)`

Perturbs the unit quaternion `q` by the rotation vector `rotationvector`. Equivalent to
`q + UnitQuaternion(rotationvector)`.
"""
function ⊞(q::UnitQuaternion, rotationvector::AbstractVector)
    length(rotationvector) == 3 || error("Must be 3-vector.")
    q + UnitQuaternion(rotationvector)
end

"""
`==(p, q)`

Checks if unit quaterions `p` and `q` are equal.
"""
==(p::UnitQuaternion, q::UnitQuaternion) = abs(p.η - q.η) < EPS && all(i->(abs(i) < EPS), p.ϵ - q.ϵ)

"""
`getindex(q, i)`

Returns the `i`-th element of the unit quaternion `q`, where `i` = 1:3 form the vector part of `q`
and `i` = 4 is the scalar part of `q`.
"""
getindex(q::UnitQuaternion, i::Integer) = i == 4 ? q.η : q.ϵ[i]

####################################################################################################
# Methods
####################################################################################################

"""
`angle(q)`

Returns the angle of rotation of the unit quaternion `q`.
"""
angle(q::UnitQuaternion) = 2.0 * atan2(norm(q.ϵ), q.η)

"""
`axis(q)`

Returns the axis of rotation of the unit quaternion `q`.
"""
axis(q::UnitQuaternion) = norm(q.ϵ) > EPS ? q.ϵ / norm(q.ϵ) : [0.0, 0.0, 1.0]

"""
`covariance(qs, qmean, weights)`

Returns the 3x3 covariance matrix of a vector of unit quaternions `qs` given the unit quaternion
mean `qmean`. The optional argument `weights` is a vector of weights used to calculate a weighted
covariance of the `qs`. If no weights are given, all unit quaternions are weighted equally.
"""
function covariance{T<:Real}(v::AbstractVector{UnitQuaternion}, qmean::UnitQuaternion,
                             weights::Vector{T} = ones(length(v)))
    normweights = map(w -> w/sum(weights), weights)
    cov = zeros(3, 3)
    for i = 1:length(v)
        diff = v[i] ⊟ qmean
        cov += normweights[i] * diff * diff'
    end
    return cov
end

"""
`covariance(qs, weights)`

Returns the 3x3 covariance matrix of a vector of unit quaternions `qs`. The optional argument `weights` is
a vector of weights used to calculate a weighted covariance of the `qs`. If no weights are given, all unit
quaternions are weighted equally.
"""
function covariance{T<:Real}(qs::AbstractVector{UnitQuaternion}, weights::Vector{T} = ones(length(qs)))
    qmean = mean(qs, weights)
    return covariance(qs, qmean, weights)
end

"""
`inv(q)`

Calculates the inverse of the unit quaternion `q`.
"""
inv(q::UnitQuaternion) = UnitQuaternion(-q.ϵ, q.η)

"""
`mean(qs, weights)`

Returns the quaternion mean of a vector of unit quaternions `qs`. The optional argument `weights` is
a vector of weights used to calculate a weighted mean of the `qs`. If no weights are given, all unit
quaternions are weighted equally.
"""
function mean{T<:Real}(qs::AbstractVector{UnitQuaternion}, weights::Vector{T} = ones(length(qs)))
    normweights = map(w -> w/sum(weights), weights)
    M = sum([w * vector(q) * vector(q)' for (w, q) in zip(normweights, qs)])
    evals, evecs = eig(M)
    largestevalindex = sortperm(evals)[end]
    return UnitQuaternion(evecs[:,largestevalindex])
end

"""
`rotatevector(q, v)`

Given a unit quaternion `q` and a 3-vector `v`, performs a vector rotation and
returns the rotated vector `v`.
"""
function rotatevector(q::UnitQuaternion, v::AbstractVector)
    length(v) == 3 || error("Must be a 3-vector.")
    (2.0 * q.η^2 - 1.0) * v + (2.0 * q.η * ×(q.ϵ)) * v + (2.0 * q.ϵ * q.ϵ') * v
end

"""
`rotateframe(q, v)`

Given a unit quaternion `q` and a 3-vector `v`, performs a coordinate frame rotation and
returns `v` expressed in the rotated coordinate frame.
"""
function rotateframe(q::UnitQuaternion, v::AbstractVector)
    length(v) == 3 || error("Must be a 3-vector.")
    (2.0 * q.η^2 - 1.0) * v - (2.0 * q.η * ×(q.ϵ)) * v + (2.0 * q.ϵ * q.ϵ') * v
end

"""
`rotationmatrix(q)`

Creates a rotation matrix (direction cosine matrix) of the rotation parameterized by the unit quaternion `q`.
"""
function rotationmatrix(q::UnitQuaternion)
    # Note, the "natural order" is used here (see Trawny or Shuster), which is different
    # from the textbook by Kuipers, who uses the "historical" (Hamiltonian) order.
    # Returned rotation matrix implements rotateframe (i.e., rotateframe(q, r) = rotationmatrix(q) * r)
    return (+(q) * ⊕(inv(q)))[1:3, 1:3]
end

"""
`log(q)`

Takes the log of the unit quaternion `q`.
"""
function log(q::UnitQuaternion)
    if q == UnitQuaternion()
        return [0.0, 0.0, 0.0]
    elseif abs(q.η) < EPS
        return π/2 * axis(q)
    else
        return angle(q)/2.0 * axis(q)
    end
end

"""
`show(io, q)`

Prints the unit quaternion `q` to the stream `io`.
"""
function show(io::IO, q::UnitQuaternion)
    out = @sprintf("ϵ = [%.3f, %.3f, %.3f], η = %.3f", q.ϵ[1], q.ϵ[2], q.ϵ[3], q.η)
    print(io, out)
end

"""
`slerp(p, q, t)`

Performs spherical linear interpolation from unit quaternion `p` to unit quaternion `q`. The
interpolation parameter `t` specifies the fraction of arc to traverse. The returned unit
quaternion parameterizes a rotation from `p` to the end point specified by `t`.
"""
slerp(p::UnitQuaternion, q::UnitQuaternion, t::Real) = (q + inv(p))^t + p

"""
`vector(q)`

Returns the elements of the unit quaternion `q` as a 4-vector.
"""
vector(q::UnitQuaternion) = [q.ϵ[1], q.ϵ[2], q.ϵ[3], q.η]

####################################################################################################
# Utility Functions
####################################################################################################

"""
`×(a)`

Returns the skew-symmetric cross product matrix of the 3-vector `a`.
"""
function ×(a::AbstractVector)
    length(a) == 3 || error("Must be a 3-vector.")
    return [0.0 -a[3] a[2]; a[3] 0.0 -a[1]; -a[2] a[1] 0.0]
end

end # module


text:
# Struct to test AbstractMvNormal methods
struct CholeskyMvNormal{M,T} <: Distributions.AbstractMvNormal
    m::M
    L::T
end

# Constructor for diagonal covariance matrices used in the tests belows
function CholeskyMvNormal(m::Vector, Σ::Diagonal)
    L = Diagonal(map(sqrt, Σ.diag))
    return CholeskyMvNormal{typeof(m),typeof(L)}(m, L)
end

Distributions.length(p::CholeskyMvNormal) = length(p.m)
Distributions.mean(p::CholeskyMvNormal) = p.m
Distributions.cov(p::CholeskyMvNormal) = p.L * p.L'
Distributions.logdetcov(p::CholeskyMvNormal) = 2 * logdet(p.L)
function Distributions.sqmahal(p::CholeskyMvNormal, x::AbstractVector)
    return sum(abs2, p.L \ (mean(p) - x))
end
function Distributions._rand!(rng::AbstractRNG, p::CholeskyMvNormal, x::Vector)
    return x .= p.m .+ p.L * randn!(rng, x)
end

@testset "Expectations" begin
    # univariate distributions
    for d in (Normal(), Poisson(2.0), Binomial(10, 0.4))
        @test Distributions.expectation(d, identity) ≈ mean(d) atol=1e-3
        @test @test_deprecated(Distributions.expectation(d, identity, 1e-10)) ≈ mean(d) atol=1e-3
    end

    # multivariate distribution
    d = MvNormal([1.5, -0.5], I)
    @test Distributions.expectation(d, identity; nsamples=10_000) ≈ mean(d) atol=1e-2
end

@testset "KL divergences" begin
    function test_kl(p, q)
        @test kldivergence(p, q) >= 0
        @test kldivergence(p, p) ≈ 0 atol=1e-1
        @test kldivergence(q, q) ≈ 0 atol=1e-1
        if p isa UnivariateDistribution
            @test kldivergence(p, q) ≈ invoke(kldivergence, Tuple{UnivariateDistribution,UnivariateDistribution}, p, q) atol=1e-1
        elseif p isa MultivariateDistribution
            @test kldivergence(p, q) ≈ invoke(kldivergence, Tuple{MultivariateDistribution,MultivariateDistribution}, p, q; nsamples=10000) atol=1e-1
        end
    end

    @testset "univariate" begin
        @testset "Beta" begin
            p = Beta(2, 10)
            q = Beta(3, 5)
            test_kl(p, q)
        end
        @testset "Categorical" begin
            @test kldivergence(Categorical([0.0, 0.1, 0.9]), Categorical([0.1, 0.1, 0.8])) ≥ 0
            @test kldivergence(Categorical([0.0, 0.1, 0.9]), Categorical([0.1, 0.1, 0.8])) ≈
                kldivergence([0.0, 0.1, 0.9], [0.1, 0.1, 0.8])
        end
        @testset "Exponential" begin
            p = Exponential(2.0)
            q = Exponential(3.0)
            test_kl(p, q)
        end
        @testset "Gamma" begin
            p = Gamma(2.0, 1.0)
            q = Gamma(3.0, 2.0)
            test_kl(p, q)
        end
        @testset "InverseGamma" begin
            p = InverseGamma(2.0, 1.0)
            q = InverseGamma(3.0, 2.0)
            test_kl(p, q)
        end
        @testset "Normal" begin
            p = Normal(0, 1)
            q = Normal(0.5, 0.5)
            test_kl(p, q)
        end
        @testset "Poisson" begin
            p = Poisson(4.0)
            q = Poisson(3.0)
            test_kl(p, q)

            # special case (test function also checks `kldivergence(p0, p0)`)
            p0 = Poisson(0.0)
            test_kl(p0, p)
        end
    end

    @testset "multivariate" begin
        @testset "AbstractMvNormal" begin
            p_mvnormal = MvNormal([0.2, -0.8], Diagonal([0.5, 0.75]))
            q_mvnormal = MvNormal([1.5, 0.5], Diagonal([1.0, 0.2]))
            test_kl(p_mvnormal, q_mvnormal)

            p_cholesky = CholeskyMvNormal([0.2, -0.8], Diagonal([0.5, 0.75]))
            q_cholesky = CholeskyMvNormal([1.5, 0.5], Diagonal([1.0, 0.2]))
            test_kl(p_cholesky, q_cholesky)

            # check consistency and mixed computations
            v = kldivergence(p_mvnormal, q_mvnormal)
            @test kldivergence(p_mvnormal, q_cholesky) ≈ v
            @test kldivergence(p_cholesky, q_mvnormal) ≈ v
            @test kldivergence(p_cholesky, q_cholesky) ≈ v
        end
    end
end


text:
# helper functions
value(z::DiffCheck) = z.value
value(x::Number) = x

derivatives(z::DiffCheck) = z.der
derivatives(x::Number) = Set{Symbol}()


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to ECI to ECEF transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  :IAU1980)
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", :IAU2000A)

# File: ./src/transformations/eci_to_ecef.jl
# ==========================================

# Functions: rECItoECEF
# ---------------------

################################################################################
#                                  IAU-76/FK5
################################################################################

## GCRF to ITRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF GCRF => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    D_ITRF_GCRF = rECItoECEF(GCRF(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = D_ITRF_GCRF*r_gcrf

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_GCRF = rECItoECEF(Quaternion, GCRF(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = vect(conj(q_ITRF_GCRF)*r_gcrf*q_ITRF_GCRF)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## J2000 to ITRF
## =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF J2000 => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    D_ITRF_J2000 = rECItoECEF(J2000(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = D_ITRF_J2000*r_j2000

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_J2000 = rECItoECEF(Quaternion, J2000(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = vect(conj(q_ITRF_J2000)*r_j2000*q_ITRF_J2000)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## TOD to ITRF
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF TOD => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]

    D_ITRF_TOD = rECItoECEF(TOD(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = D_ITRF_TOD*r_tod

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_TOD = rECItoECEF(Quaternion, TOD(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = vect(conj(q_ITRF_TOD)*r_tod*q_ITRF_TOD)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## MOD to ITRF
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF MOD => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_mod  = [5094.02837450; 6127.87081640; 6380.24851640]

    D_ITRF_MOD = rECItoECEF(MOD(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = D_ITRF_MOD*r_mod

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_MOD = rECItoECEF(Quaternion, MOD(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = vect(conj(q_ITRF_MOD)*r_mod*q_ITRF_MOD)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## TEME to ITRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF TEME => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    D_ITRF_TEME = rECItoECEF(TEME(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = D_ITRF_TEME*r_teme

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_TEME = rECItoECEF(Quaternion, TEME(), ITRF(), JD_UTC, eop_iau1980)
    r_itrf = vect(conj(q_ITRF_TEME)*r_teme*q_ITRF_TEME)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## GCRF to PEF
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECItoECEF GCRF => PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    D_PEF_GCRF = rECItoECEF(GCRF(), PEF(), JD_UTC, eop_iau1980)
    r_pef = D_PEF_GCRF*r_gcrf

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    q_PEF_GCRF = rECItoECEF(Quaternion, GCRF(), PEF(), JD_UTC, eop_iau1980)
    r_pef = vect(conj(q_PEF_GCRF)*r_gcrf*q_PEF_GCRF)

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4
end

## J2000 to PEF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECItoECEF J2000 => PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    D_PEF_J2000 = rECItoECEF(J2000(), PEF(), JD_UTC, eop_iau1980)
    r_pef = D_PEF_J2000*r_j2000

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    q_PEF_J2000 = rECItoECEF(Quaternion, J2000(), PEF(), JD_UTC, eop_iau1980)
    r_pef = vect(conj(q_PEF_J2000)*r_j2000*q_PEF_J2000)

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4
end

## TOD to PEF
## ==========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get the same result by using:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#
# Notice, however, that this result is obtained by using the correct UT1.
#
################################################################################

@testset "Function rECItoECEF TOD => PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_tod  = [5094.51620300; 6127.36527840; 6380.34453270]

    D_PEF_TOD = rECItoECEF(TOD(), PEF(), JD_UTC, eop_iau1980)
    r_pef = D_PEF_TOD*r_tod

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    q_PEF_TOD = rECItoECEF(Quaternion, TOD(), PEF(), JD_UTC, eop_iau1980)
    r_pef = vect(conj(q_PEF_TOD)*r_tod*q_PEF_TOD)

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    # No EOP corrections
    # ==================

    # Notice that if we use JD_UT1, then the correction to TT will be wrong.
    # However, this will lead to a much smaller error than assuming that UTC =
    # UT1.

    JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400
    r_tod  = [5094.51478040; 6127.36646120; 6380.34453270]

    D_PEF_TOD = rECItoECEF(TOD(), PEF(), JD_UT1)
    r_pef = D_PEF_TOD*r_tod

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7

    q_PEF_TOD = rECItoECEF(Quaternion, TOD(), PEF(), JD_UT1)
    r_pef = vect(conj(q_PEF_TOD)*r_tod*q_PEF_TOD)

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7
end

## PEF to MOD
## ==========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get the same result by using:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#
# Notice, however, that this result is obtained by using the correct UT1.
#
################################################################################

@testset "Function rECItoECEF MOD => PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_mod  = [5094.02837450; 6127.87081640; 6380.24851640]

    D_PEF_MOD = rECItoECEF(MOD(), PEF(), JD_UTC, eop_iau1980)
    r_pef = D_PEF_MOD*r_mod

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    q_PEF_MOD = rECItoECEF(Quaternion, MOD(), PEF(), JD_UTC, eop_iau1980)
    r_pef = vect(conj(q_PEF_MOD)*r_mod*q_PEF_MOD)

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    # No EOP corrections
    # ==================

    # Notice that if we use JD_UT1, then the correction to TT will be wrong.
    # However, this will lead to a much smaller error than assuming that UTC =
    # UT1.

    JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400
    r_mod  = [5094.02901670; 6127.87093630; 6380.24788850]

    D_PEF_MOD = rECItoECEF(MOD(), PEF(), JD_UT1)
    r_pef = D_PEF_MOD*r_mod

    @test r_pef[1] ≈ -1033.47503130 atol=8e-6
    @test r_pef[2] ≈ +7901.30558560 atol=8e-6
    @test r_pef[3] ≈ +6380.34453270 atol=8e-6

    q_PEF_MOD = rECItoECEF(Quaternion, MOD(), PEF(), JD_UT1)
    r_pef = vect(conj(q_PEF_MOD)*r_mod*q_PEF_MOD)

    @test r_pef[1] ≈ -1033.47503130 atol=8e-6
    @test r_pef[2] ≈ +7901.30558560 atol=8e-6
    @test r_pef[3] ≈ +6380.34453270 atol=8e-6
end

## TEME to PEF
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECItoECEF TEME => PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    D_PEF_TEME = rECItoECEF(TEME(), PEF(), JD_UTC, eop_iau1980)
    r_pef = D_PEF_TEME*r_teme

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4

    q_PEF_TEME = rECItoECEF(Quaternion, TEME(), PEF(), JD_UTC, eop_iau1980)
    r_pef = vect(conj(q_PEF_TEME)*r_teme*q_PEF_TEME)

    @test r_pef[1] ≈ -1033.47503130 atol=3e-4
    @test r_pef[2] ≈ +7901.30558560 atol=3e-4
    @test r_pef[3] ≈ +6380.34453270 atol=3e-4
end

################################################################################
#                                IAU-2006/2010
################################################################################

## GCRF => TIRS
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#
# one gets:
#
#   r_tirs = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECItoECEF GCRF => TIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]

    D_TIRS_GCRF = rECItoECEF(GCRF(), TIRS(), JD_UTC, eop_iau2000a)
    r_tirs = D_TIRS_GCRF*r_gcrf

    @test r_tirs[1] ≈ -1033.47503120 atol=3e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=3e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=3e-4

    q_TIRS_GCRF = rECItoECEF(Quaternion, GCRF(), TIRS(), JD_UTC, eop_iau2000a)
    r_tirs = vect(q_TIRS_GCRF\r_gcrf*q_TIRS_GCRF)

    @test r_tirs[1] ≈ -1033.47503120 atol=3e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=3e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=3e-4
end

## GCRF => ITRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF GCRF => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]

    D_ITRF_GCRF = rECItoECEF(GCRF(), ITRF(), JD_UTC, eop_iau2000a)
    r_itrf = D_ITRF_GCRF*r_gcrf

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_GCRF = rECItoECEF(Quaternion, GCRF(), ITRF(), JD_UTC, eop_iau2000a)
    r_itrf = vect(q_ITRF_GCRF\r_gcrf*q_ITRF_GCRF)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end

## CIRS => TIRS
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_tirs = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECItoECEF CIRS => TIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]

    D_TIRS_CIRS = rECItoECEF(CIRS(), TIRS(), JD_UTC, eop_iau2000a)
    r_tirs = D_TIRS_CIRS*r_cirs

    @test r_tirs[1] ≈ -1033.47503120 atol=3e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=3e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=3e-4

    q_TIRS_CIRS = rECItoECEF(Quaternion, CIRS(), TIRS(), JD_UTC, eop_iau2000a)
    r_tirs = vect(q_TIRS_CIRS\r_cirs*q_TIRS_CIRS)

    @test r_tirs[1] ≈ -1033.47503120 atol=3e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=3e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=3e-4
end

## CIRS => ITRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
################################################################################

@testset "Function rECItoECEF CIRS => ITRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]

    D_ITRF_CIRS = rECItoECEF(CIRS(), ITRF(), JD_UTC, eop_iau2000a)
    r_itrf = D_ITRF_CIRS*r_cirs

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4

    q_ITRF_CIRS = rECItoECEF(Quaternion, CIRS(), ITRF(), JD_UTC, eop_iau2000a)
    r_itrf = vect(q_ITRF_CIRS\r_cirs*q_ITRF_CIRS)

    @test r_itrf[1] ≈ -1033.4793830 atol=3e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=3e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=3e-4
end


text:
# ./src/lif_current.jl

############
# Current-based (classic) LIF

# threshold-linear input-output function
struct NTLIF <: NeuronType
  τ::Float64 # time constant
	v_threshold::Float64 # fixed spiking threshold
	v_reset::Float64 # reset after spike
	τ_refractory::Float64 # refractory time
	τ_output_decay::Float64 # decay of postsynaptic currents
end

struct PSLIF{NT} <: PSGeneralCurrentIFType{NT}
  neurontype::NT
  n::Int64 # pop size
  state_now::Vector{Float64}
  input::Vector{Float64}
	alloc_dv::Vector{Float64}
	last_fired::Vector{Float64}
	isfiring::BitArray{1}
	isrefractory::BitArray{1}
end
function PSLIF(p::NTLIF,n::Int64)
  PSLIF(p,n,ntuple(_-> zeros(Float64,n),4)...,ntuple(_-> falses(n),2)...)
end
function PSLIF(τ::Float64,v_threshold::Float64,
	  v_reset::Float64 ,
	  τ_refractory::Float64 ,
	  τ_output_decay::Float64, n::Int64)
  nt=NTLIF(τ,v_threshold,v_reset,τ_refractory,τ_output_decay)
  return PSLIF(nt,n)
end
function reset!(ps::PSLIF)
  fill!(ps.state_now,0.0)
  fill!(ps.input,0.0)
  fill!(ps.last_fired,-Inf)
  fill!(ps.isfiring,false)
  fill!(ps.isrefractory,false)
  return nothing
end

# Connection : ConnGeneralIF, defined in firingneruons_shared.jl

function local_update!(t_now::Float64,dt::Float64,ps::PSLIF)
	# computes the update to internal voltage, given the total input
  # dv =  (v_leak - v ) dt / τ + I dt / Cap
  copy!(ps.alloc_dv,ps.state_now)  # v
  ps.alloc_dv .-= ps.input  # (v-I)
  lmul!(-dt/ps.neurontype.τ,ps.alloc_dv) # du =  dt/τ (-v+I)
  ps.state_now .+= ps.alloc_dv # v_{t+1} = v_t + dv
	# update spikes and refractoriness, and end
  # function defined in firingneruons_shared.jl
  return _spiking_state_update!(ps.state_now,ps.isfiring,ps.isrefractory,ps.last_fired,
    t_now,ps.neurontype.τ_refractory,ps.neurontype.v_threshold,ps.neurontype.v_reset)
end

# analytic stuff

function expected_period_norefr(neu::NTLIF,input)
	return expected_period_norefr(neu.τ,neu.v_reset,neu.v_threshold,input)
end

function expected_period_norefr(τ::Real,vr::Real,vth::Real,input::Real)
	return τ*log( (input-vr)/(input-vth) )
end

# forward_signal(...) defined in firingneruons_shared.jl

#=
# threshold-linear input-output function
struct PopLIF <: Population
  τ::Float64 # time constant
	v_threshold::Float64 # fixed spiking threshold
	v_reset::Float64 # reset after spike
	τ_refractory::Float64 # refractory time
	τ_output_decay::Float64 # decay of postsynaptic currents
end

struct PSLIF{P} <: PopulationState
  population::P
  state_now::Vector{Float64}
  input::Vector{Float64}
	alloc_dv::Vector{Float64}
	last_fired::Vector{Float64}
	isfiring::BitArray{1}
	isrefractory::BitArray{1}
end
function PSLIF(p::Population)
  PSLIF(p, ntuple(_-> zeros(Float64,p.n),4)...,ntuple(_-> falses(p.n),2)...)
end

struct ConnectionLIF <: Connection
  postps::PSLIF # postsynaptic population state
  preps::PSLIF # presynaptic population state
  adjR::Vector{Int64} # adjacency matrix, rows
  adjC::Vector{Int64} # adjacency matrix, columns
  weights::SparseMatrixCSC{Float64,Int64}
  post_current::Vector{Float64}
end

function ConnectionLIF(post::PSLIF,weights::SparseMatrixCSC,pre::PSLIF)
  aR,aC,_ = findnz(weights)
	npost = post.population.n
  ConnectionLIF(post,pre,aR,aC,weights,zeros(Float64,npost))
end

function dynamics_step!(t_now::Float64,dt::Float64,ps::PSLIF)
	# computes the update to internal voltage, given the total input
  copy!(ps.alloc_dv,ps.state_now)  # v
  ps.alloc_dv .-= ps.input  # (v-I)
  lmul!(-dt/ps.population.τ,ps.alloc_dv) # du =  dt/τ (-v+I)
  ps.state_now .+= ps.alloc_dv # v_{t+1} = v_t + dv
	# now spikes and refractoriness
	reset_spikes!(ps)
	@inbounds @simd for i in eachindex(ps.state_now)
		if ps.state_now[i] > ps.population.v_threshold
			ps.isfiring[i] = true
			ps.state_now[i] =  ps.population.v_reset
			ps.last_fired[i] = t_now
			ps.isrefractory[i] = true
		# check only when refractory
		elseif ps.isrefractory[i] &&
				( (t_now-ps.last_fired[i]) >= ps.population.τ_refractory)
			ps.isrefractory[i] = false
		end
	end
	#
  return nothing
end

function send_signal!(t_now::Float64,conn::ConnectionLIF)
	post_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
	τ_decay = conn.preps.population.τ_output_decay
	for _pre in findall(conn.preps.isfiring)
		_posts_nz = nzrange(conn.weights,_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
			post_idx = post_idxs[_pnz]
			# update the post currents even when refractory
				conn.post_current[post_idx] +=
						weightsnz[_pnz] / τ_decay
		end
	end
	# finally , add the currents to postsynaptic input
	# ONLY non-refractory neurons
	post_refr = findall(conn.postps.isrefractory) # refractory ones
	@inbounds @simd for i in eachindex(conn.postps.input)
		if !(i in post_refr)
			conn.postps.input[i] += conn.post_current[i]
		end
	end
  return nothing
end

# postsynaptic currents decay in time
function dynamics_step!(t_now::Real,dt::Real,conn::ConnectionLIF)
	τ_decay = conn.preps.population.τ_output_decay
	@inbounds @simd for i in eachindex(conn.post_current)
		conn.post_current[i] -= dt*conn.post_current[i] / τ_decay
	end
	return nothing
end

@inline function send_signal!(t::Real,input::PopInputStatic{P}) where P<:PSLIF
	refr = findall(input.population_state.isrefractory) # refractory ones
  @inbounds @simd for i in eachindex(input.population_state.input)
		if !(i in refr)
  		input.population_state.input[i] += input.h[i]
		end
	end
  return nothing
end

		for cc = 1:Ncells
    ###plasticity detectors, forward euler###
		  x[cc]  += -dt*x[cc] /p.tau_i #inhibitory traces
			o1[cc] += -dt*o1[cc]/p.tau_minus #triplet terms
			o2[cc] += -dt*o2[cc]/p.tau_y
			r1[cc] += -dt*r1[cc]/p.tau_plus
			r2[cc] += -dt*r2[cc]/p.tau_x

      ## NEURON DYNAMICS
			if t > (lastSpike[cc] + p.taurefrac) #only after refractory period
        ###CONDUCTANCE DYNAMICS####
			  #forward euler
				Ie[cc] += -dt*Ie[cc]/p.tauedecay
				Ip[cc] += -dt*Ip[cc]/p.taupdecay
				Is[cc] += -dt*Is[cc]/p.tausdecay
				Iv[cc] += -dt*Iv[cc]/p.tauvdecay
        ###CONDUCTANCE DYNAMICS END###

        ###LIF NEURONS###
				#forward euler
				if cc <= p.Ne #excitatory neuron
					dv =  -v[cc]/p.taue + Ie[cc] - Ip[cc] - Is[cc] +  I_stim[cc]+ I_E0[cc]
					v[cc] += dt*dv
					if  v[cc] > p.vth_e
						spiked[cc] = true
					end
				elseif p.Ne < cc <= (p.Ne+p.Np) #pv neurons
					dv = -v[cc]/p.taup + Ie[cc] - Ip[cc] - Is[cc] + I_stim[cc] + I_P0[cc]
					v[cc] += dt*dv
					if v[cc] > p.vth_i
						spiked[cc] = true
					end
				elseif (p.Ne+p.Np) < cc <= (p.Ne+p.Np+p.Ns)#sst neurons
					dv = -v[cc]/p.taus + Ie[cc] -Iv[cc]+ I_S0[cc] + I_stim[cc]
					v[cc] += p.dt*dv
					if v[cc] > p.vth_i
						spiked[cc] = true
					end
				else #vip neurons
					dv = -v[cc]/p.tauv + Ie[cc] -Is[cc] -Ip[cc] + I_V0[cc] + I_stim[cc]
					v[cc] += dt*dv
					if v[cc] > p.vth_i
						spiked[cc] = true
					end
				end
        ###LIF NEURONS END###

        ###UPDATE WHEN SPIKE OCCURS
				if spiked[cc] #spike occurred
					v[cc] = p.vre #voltage back to reset potential
					lastSpike[cc] = t #record last spike time
					ns[cc] += 1 #lists number of spikes per neuron
					if ns[cc] <= p.Nspikes #spike time are only record for ns < Nspikes
						times[cc,ns[cc]] = t #recording spiking times
					end
					# for synaptic updates other than input upon every presynaptic spike
					for dd = 1:Ncells
						if cc <= p.Ne
							Ie[dd]+= weights[cc,dd]/p.tauedecay
						elseif p.Ne < cc <= (p.Np+p.Ne)
							Ip[dd]+= weights[cc,dd]/p.taupdecay
						elseif (p.Ne+p.Np) < cc <=(p.Np+p.Ne+p.Ns)
							Is[dd]+= weights[cc,dd]/p.tausdecay
						else
							Iv[dd]+= weights[cc,dd]/p.tauvdecay
						end
					end
				end #end if(spiked)
			end #end if(not refractory)

=#


text:
@symbol_func function G_CR(cur_reactor::AbstractReactor)
  cur_epsilon = epsilon_b(cur_reactor)

  ( isa(cur_epsilon, AbstractFloat) && cur_epsilon >= 1 ) && return NaN

  cur_left_part = safe_symbol(cur_reactor, :B_0) ^ 2

  cur_left_part /= (8/5) * cur_reactor.pi

  cur_left_part /= cur_reactor.sigma_TF

  cur_left_part /= ( 1 - cur_epsilon )

  cur_right_part = 4 * cur_epsilon

  cur_right_part /= 1 + cur_epsilon

  cur_right_part += log(
    ( 1 + cur_epsilon ) / ( 1 - cur_epsilon )
  )

  cur_G = cur_left_part

  cur_G *= cur_right_part

  cur_G
end


text:
import Base: start, next, done

"""
Data generators should implement the Iterator interface to produce minibatches:

* `start(iter) => state`
* `next(iter,state) => (item,state)`
* `done(iter,state) => Bool`

ItemTensor is a data generator that is constructed from a tuple of
arrays (x[d...,i], y[d...,i],...) where the last dimension is
interpreted as the item index.  Produces tuples of minibatches. For
non-sequential data.

"""
type ItemTensor; data; rng; datasize; epochsize; batchsize; bootstrap; shuffle; batch;
    function ItemTensor(x...; rng=MersenneTwister(), epoch=ccount(x[1]), batch=128, bootstrap=false, shuffle=false)
        nx = ccount(x[1])
        all(xi->ccount(xi)==nx, x) || error("Item count mismatch")
        idx = (shuffle ? shuffle!(rng,[1:nx;]) : nothing)
        buf = map(xi->itembatch(xi,batch), x)
        new(x, rng, nx, epoch, batch, bootstrap, idx, buf)
    end
end

start(d::ItemTensor)=(d.shuffle != nothing && shuffle!(d.rng, d.shuffle); 0)

done(d::ItemTensor, n)=(n + d.batchsize > d.epochsize)

function next(d::ItemTensor, n)
    idx = nextidx(d,n)
    for i=1:length(d.data)
        cslice!(d.batch[i], d.data[i], idx) # 549/556
    end
    (d.batch, n+length(idx))
end

function nextidx(d, n)
    nx = d.datasize
    nb = min(d.batchsize, d.epochsize-n)
    if d.bootstrap
        ix = rand(d.rng, 1:nx, nb)
    elseif d.shuffle != nothing
        i1 = mod1(n+1, nx)
        i2 = min(i1+nb-1, nx)
        ix = d.shuffle[i1:i2]
        while length(ix) < nb
            shuffle!(d.rng, d.shuffle)
            i2 = min(nb - length(ix), nx)
            ix = [ix; d.shuffle[1:i2]]
        end
    else
        i1 = mod1(n+1, nx)
        i2 = min(i1+nb-1, nx)
        ix = (i1:i2)
        while length(ix) < nb
            i2 = min(nb - length(ix), nx)
            ix = [ix; 1:i2]
        end
    end
    length(ix) == nb || error()
    return ix
end

# DynamicArrayCPU{T}(::Type{T}, d::Dims)=KUdense(Array,T,d)
# DynamicArrayGPU{T}(::Type{T}, d::Dims)=KUdense(CudaArray,T,d)
# SparseArrayCPU{T}(::Type{T}, d::Dims)=spzeros(T,d...)
# SparseArrayGPU{T}(::Type{T}, d::Dims)=CudaSparseMatrixCSR(spzeros(T,d...))

# SparseArrayCPU{T}(::Type{T}, d::Dims)=spzeros(T,d...)
# DynamicArrayCPU{T}(::Type{T}, d::Dims)=KUdense(Array,T,d)
# itembatch(x,n)=(issparse(x)?SparseArrayCPU:DynamicArrayCPU)(eltype(x),csize(x,n))

itembatch(x,n)=(issparse(x) ? spzeros(eltype(x),csize(x,n)...) : Array(eltype(x),csize(x,n)))


text:
load("src/neural.jl")

# Force JIT'ing
net = init_network([2, 3, 2])
train(net, [0.15, 0.7],[0.1, 0.9])
net_eval(net, [0.15, 0.7])

@elapsed for i = 1:10_000
  net = init_network([2, 3, 2])
  train(net, [0.15, 0.7],[0.1, 0.9])
  net_eval(net, [0.15, 0.7])
end


text:
using Lale
import Pandas
using DataFrames
using PyCall

fetch = pyimport("lale.datasets.openml").fetch
(train_X, train_y), (test_X, test_y) = fetch("breast-cancer", "classification", preprocess=true)
train_X = DataFrame(train_X,:auto)
test_X = DataFrame(test_X,:auto)

MinMaxScaler = laleoperator("MinMaxScaler")
PCA = laleoperator("PCA")
NoOp = laleoperator("NoOp","lale")
LR= laleoperator("LogisticRegression")
RF= laleoperator("RandomForestClassifier")
SMOTE = laleoperator("SMOTE","imblearn")
CondensedNearestNeighbour = laleoperator("CondensedNearestNeighbour","imblearn")
SMOTEENN = laleoperator("SMOTEENN","imblearn")
accuracy_score=pyimport("sklearn.metrics").accuracy_score
Hyperopt = laleoperator("Hyperopt","lale")

pipeline_without_correction =  MinMaxScaler() >> PCA() >> RF()
optimizer = Hyperopt(pipeline_without_correction, max_evals = 10, scoring="roc_auc")
trained_optimizer = fit(optimizer,train_X, train_y)
predictions = predict(trained_optimizer,test_X)
lalesummary(trained_optimizer)
score(:accuracy,test_y,predictions)

# SMOTE
pipeline_with_correction =  SMOTE(MinMaxScaler() >> PCA() >> RF())
optimizer = Hyperopt(pipeline_with_correction, max_evals = 10, scoring="roc_auc")
trained_optimizer = fit(optimizer,train_X, train_y)
predictions = predict(trained_optimizer,test_X)
lalesummary(trained_optimizer)
score(:accuracy,test_y,predictions)

# CondensedNearestNeighbour
pipeline_with_correction =  CondensedNearestNeighbour(MinMaxScaler() >> PCA() >> RF())
optimizer = Hyperopt(pipeline_with_correction, max_evals = 10, scoring="roc_auc")
trained_optimizer = fit(optimizer,train_X, train_y)
predictions = predict(trained_optimizer,test_X)
lalesummary(trained_optimizer)
score(:accuracy,test_y,predictions)

# SMOTEENN
pipeline_with_correction =  SMOTEENN(MinMaxScaler() >> PCA() >> RF())
optimizer = Hyperopt(pipeline_with_correction, max_evals = 10, scoring="roc_auc")
trained_optimizer = fit(optimizer,train_X, train_y)
predictions = predict(trained_optimizer,test_X)
lalesummary(trained_optimizer)
score(:accuracy,test_y,predictions)


text:
# Remove `*` as infix operator between tensors
function Base.:*(S1::AbstractTensor, S2::AbstractTensor)
    error("use `⋅` (`\\cdot`) for single contraction and `⊡` (`\\boxdot`) for double contraction instead of `*`")
end

for f in (:transpose, :adjoint)
    @eval function LinearAlgebra.$f(::Vec)
        throw(ArgumentError("the (no-op) $($f) is discontinued for `Tensors.Vec`"))
    end
end


text:
using Dierckx
using DelimitedFiles

function opal_test_spl()
    # function opal_test_spl()
    #     opal_all = myarray=(open(readdlm,"data/opal_all.txt"))
    #     i = 78
    #     table = opal_all[182+ (i-1)*69 + (i-1)*4:182+i*(69)+(i-1)*4,:]
    #     replace!(table, 9.999 => 0);
    #     # replace!(table, "" => 0);
    #     logT = table[:,1]
    table = myarray=(open(readdlm,"data/solar.txt"))
    replace!(table, "" => 0);
    logT = table[:,1]
    logR = Array(range(-8, 1., step=0.5))
    logkappa = table[:,2:end]

    spline = Spline2D(logT, logR, logkappa)
    return spline
end

function opacity(spl, rho, T)
    R = rho ./ (T / 1e6) .^ 3
    # @info(log10(T), log10(R))
    return 10 .^ spl(log10.(T), log10.(R))
end


text:
function pinvertible(A::Symmetric, tol = eps(real(float(one(eltype(A))))))
    eigval, eigvect = eigen(A)
    small = eigval .<= tol
    if any(small)
        @warn "estimated covariance matrix of moment conditions not of full rank.
                 model tests should be interpreted with caution."
        eigval[small] .= 0
        return Symmetric(eigvect' * Diagonal(eigval) * eigvect)
    else
        return A
    end
end


text:
using KitBase

D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d2f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "maxwell"
    D[:cfl] = 0.8
    D[:maxTime] = 0.01

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 45
    D[:y0] = 0.0
    D[:y1] = 1.0
    D[:ny] = 45
    D[:pMeshType] = "uniform"
    D[:nxg] = 0
    D[:nyg] = 0

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 28
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 28
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.075
    D[:mach] = 2.0
    D[:prandtl] = 1.0
    D[:inK] = 1.0
    D[:omega] = 0.72
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5

    D[:uLid] = 0.15
    D[:vLid] = 0.0
    D[:tLid] = 1.0
end

ks = SolverSet(D)

ctr, a1face, a2face = init_fvm(ks, ks.pSpace; structarray = true)
ctr0, a1face0, a2face0 = init_fvm(ks, ks.pSpace; structarray = false)

using BenchmarkTools

@btime KitBase.solve!(ks, ctr, a1face, a2face, 0.0)
```~0.010509362```

@btime KitBase.solve!(ks, ctr0, a1face0, a2face0, 0.0)
```~0.010509326```


text: using Pomodorousing Base.Test@testset "All Tests" begin  if ( endswith(pwd(), "/test") ) ; cd("..") ; end  Julz.include_all_files("$(pwd())/test")endreturn

max_stars_repo_name: marcelovmaciel/Pomodoro

max_issues_repo_name: marcelovmaciel/Pomodoro

max_forks_repo_name: marcelovmaciel/Pomodoro


text:
module CommunalHelperShieldedRefill

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/ShieldedRefill" ShieldedRefill(
    x::Integer,
    y::Integer,
    twoDashes::Bool=false,
    oneUse::Bool=false,
    bubbleRepel::Bool=true,
)

const placements = Ahorn.PlacementDict(
    "Shielded Refill (Communal Helper)" => Ahorn.EntityPlacement(
        ShieldedRefill,
    ),
    "Shielded Refill (Two Dashes) (Communal Helper)" => Ahorn.EntityPlacement(
        ShieldedRefill,
        "point",
        Dict{String,Any}(
            "twoDashes" => true
        ),
    ),
)

const spriteOneDash = "objects/refill/idle00"
const spriteTwoDash = "objects/refillTwo/idle00"

getSprite(entity::ShieldedRefill) = get(entity.data, "twoDashes", false) ? spriteTwoDash : spriteOneDash

function Ahorn.selection(entity::ShieldedRefill)
    x, y = Ahorn.position(entity)
    sprite = getSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ShieldedRefill)
    Ahorn.Cairo.save(ctx)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    Ahorn.drawCircle(ctx, 0, 0, 8, (1.0, 1.0, 1.0, 1.0))
    Ahorn.Cairo.restore(ctx)

    sprite = getSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
import SimpleMolecule.AbstractMolecularContainer
import SimpleMolecule.obtain
import SimpleMolecule

"""
A test for obtain the value of
a field in an object
of type Atom.

Arguments
----------
obj:Atom
    object of type Atom

field:Symbol
    the requested field name

solution:Any
    answer key
"""
function test_obtain{T<:AbstractMolecularContainer}(
    obj::T,
    field::Symbol,
    solution::Any)

    println("----------------------------------------")
    print_with_color(:blue, "Test obtain(obj,:$(field))\n\n")

    answer = obtain(obj, field)

    if issubtype(typeof(answer), AbstractFloat)
        tol = 1e-3
        @test abs(answer - solution) < tol
    else
        @test answer == solution
    end

    print_with_color(:green, "\nVERIFIED! obtain(obj,:$(field)) for type Atom\n")
    println("----------------------------------------\n\n")

    return true
end


text:
module DPC

include("model.jl")
include("simulate.jl")
include("utils.jl")
include("plotting.jl")

end


text:
module LocalPolynomialFitting

export lpparams, lpfit,
       optbandn,
       epanwin

"""Epanechnikov's window"""
epanwin(n::Integer) = 3/4 * (1 - linspace(-1,1,n+2).^2)[2:end-1]

include("lanshammar.jl")

end


text:
# initalize an array
function Init!(X::Array{Float64,2},L)
    siz=size(X)[1]
    h=L/siz
    for j=0:siz-1
        for i=0:siz-1
            if i>floor(siz//8) && i<floor(siz//2)+floor(siz//8)
                X[i+1,j+1]=1.-2*(i-floor(siz/8))*h/L;
            else
                X[i+1,j+1]=0.0
            end
        end
    end
end

# proc[1-5]! functions make the same computation, with different programing
# styles.

function proc1!(In::Array{Float64,2},Out::Array{Float64,2},niter::Int64)
    const siz=size(In)[1]
    h2::Float64= (1./siz)^2
    for it=1:niter
        Out[2:siz-1,2:siz-1]=
        h2*(In[1:siz-2,2:siz-1 ] + In[2:siz-1,1:siz-2]-
            4.0*In[2:siz-1,2:siz-1]+
            In[3:siz,2:siz-1]+In[2:siz-1,3:siz])
        In,Out=Out,In
    end
end
function proc2!(In::Array{Float64,2},Out::Array{Float64,2},niter::Int64)
    const siz=size(In)[1]
    h2::Float64= (1./siz)^2
    for it=1:niter
       @. Out[2:siz-1,2:siz-1]=
        h2*(In[1:siz-2,2:siz-1 ] + In[2:siz-1,1:siz-2]-
            4.0*In[2:siz-1,2:siz-1]+
            In[3:siz,2:siz-1]+In[2:siz-1,3:siz])
        In,Out=Out,In
    end
end
function proc3!(In::Array{Float64,2},Out::Array{Float64,2},niter::Int64)
    const siz=size(In)[1]
    h2::Float64= (1./siz)^2
    for it=1:niter
        for j=2:siz-1
            for i=2:siz-1
                Out[i,j]=
                h2*(In[i-1,j] + In[i,j-1]-
                    4.0*In[i,j]+
                    In[i+1,j]+In[i,j+1])
            end
        end
        In,Out=Out,In
    end
end
function proc4!(In::Array{Float64,2},Out::Array{Float64,2},niter::Int64)
    const siz=size(In)[1]
    h2::Float64= (1./siz)^2
    for it=1:niter
        for j=2:siz-1
            @simd for i=2:siz-1
                Out[i,j]=
                h2*(In[i-1,j] + In[i,j-1]-
                    4.0*In[i,j]+
                    In[i+1,j]+In[i,j+1])
            end
        end
        In,Out=Out,In
    end
end
function proc5!(In::Array{Float64,2},Out::Array{Float64,2},niter::Int64)
    const siz=size(In)[1]
    h2::Float64= (1./siz)^2
    for it=1:niter
        for i=2:siz-1
            @simd for j=2:siz-1
                Out[i,j]=
                h2*(In[i-1,j] + In[i,j-1]-
                    4.0*In[i,j]+
                    In[i+1,j]+In[i,j+1])
            end
        end
        In,Out=Out,In
    end
end
# measure as precisely as possible the execution time for a given
# p= proc[1-5]!
function test(p,In::Array{Float64,2},Out::Array{Float64,2},nit::Int)

    niter=nit
    # be sure to run once before actually running the benchmark!
    Init!(In,1.)
    Init!(Out,1.)
    p(In,Out,niter)
    #
    T=0.
    while true
        Init!(In,1.)
        Init!(Out,1.)
        t1 = time_ns()
        p(In,Out,niter)
        t = (time_ns() -t1)/niter
        if abs(t-T)/t<0.025
            break
    else
            T=t
            if niter> 150000
                error("too many iterations",niter)
            else
                niter*=2
            end
    end

end
    T,niter
end
DD=Dict("proc1!"=>"Vectorisation with braces             ",
        "proc2!"=>"Vectorisation with @. and braces      ",
        "proc3!"=>"Unrolled loop (good order)            ",
        "proc4!"=>"Unrolled loop with @simd (good order) ",
        "proc5!"=>"Unrolled loop with @simd (bad order)  "
        )

fw=open("RunningOn"*gethostname()*"_lapl_2","w")

# computation starts here:
siz=32
sizemax=2049
const niter=2
while siz<sizemax
    println("size: ",siz)
    In=Array{Float64,2}(siz,siz)
    Out=Array{Float64,2}(siz,siz)
    tbest=10.^29
    best=0
    t=0.0
    for p= [proc1!,proc2!,proc3!,proc4!,proc5!]
        t,it=test(p,In,Out,niter)
        if t<tbest
            tbest=t
            best=p
        end
        t*=10.0^(-9)
        println(DD[string(p)]," : t= ",t," seconds ")
    end

    nflops=6*(siz-2)^2
    flops=nflops/tbest
    println("\nbest: ",DD[string(best)])

    write(fw,string(siz)," ",string(tbest*10.0^(-9)),"\n")

    println("nb. flops: ",nflops, ", Gflops/s (best): ",flops)
    println("-------")

    siz*=2
    println()
end
close(fw)

max_stars_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_issues_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_forks_repo_name: wolfv/BenchmarksPythonJuliaAndCo


text:
using DataAssimilation
using Test

@testset "DataAssimilation.jl" begin
    # Write your tests here.
end


text:
using Glob
using FileIO
using CSV
using DataFrames
using Statistics
using Colors
using Plots
using Plots.PlotMeasures
using LaTeXStrings

#using PackageCompiler
#create_sysimage(:Plots, sysimage_path="sys_plots.so")

pyplot()

include("../config.jl")
include("../utils/utils.jl")

function get_data(fix_size ; fix)
    if fix == "size"
        size_str =  get_tensorinfo(fix_size)
        result_path = "../../result/synthetic/size$size_str/"
    elseif fix == "rank"
        rank_str =  get_tensorinfo(fix_size)
        result_path = "../../result/synthetic/reqrank$rank_str/"
    else
        error("invalid fix $fix")
    end
    horizon_list = glob("$result_path*")

    horizontal_ax = []
    for horizon in horizon_list
        tmp = split(horizon, "/")[end]
        rank_str = split(tmp,"_")[2:end]
        rank = parse.(Int, rank_str)
        append!(horizontal_ax, rank[1])
        #append!(horizontal_ax, prod(rank))
    end

    KLerror = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )
    LSerror = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )
    runtime = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )

    KLerror_std = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )
    LSerror_std = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )
    runtime_std = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )

    for horizon in horizon_list
        for method in methods
            target_path = horizon*"/"*method
            datalists = glob("$target_path*")
            tmp_for_runtime = []
            tmp_for_LSerror = []
            tmp_for_KLerror = []
            for datalist in datalists
                data = load(datalist)
                append!(tmp_for_runtime, data["runtime"])
                append!(tmp_for_LSerror, data["LSerror"])
                append!(tmp_for_KLerror, data["KLerror"])
            end
            append!(LSerror[method], mean( tmp_for_LSerror ))
            append!(KLerror[method], mean( tmp_for_KLerror ))
            append!(runtime[method], mean( tmp_for_runtime ))
            append!(LSerror_std[method], std( tmp_for_LSerror ))
            append!(KLerror_std[method], std( tmp_for_KLerror ))
            append!(runtime_std[method], std( tmp_for_runtime ))
        end
    end
    return horizontal_ax, KLerror, LSerror, runtime, KLerror_std, LSerror_std, runtime_std
end

function plot_error(horizontal_ax, error, error_std, fixed_target ;fix, cost="LS", rel=false, xlog=false)
    ylabel = "$cost error"
    if rel == true
        rel_error = Dict( method for method = zip(methods, [ [] for i=1:length(methods)] ) )
        for method in methods
            rel_error[method] = ( error[method] - error["TLrRR"] ) ./ error[method]
            error_std[method] .= 0.0
        end
        ylabel = "Relative $cost error"
        error = rel_error
    end

    plt_time = plot()
    if xlog
        plt_time = plot!(plt_time, xaxis=:log)
    end

    if fix == "size"
        tensor_size = fixed_target
        title     = "Random $tensor_size tensor"
        xlabel    = "Tensor rank"
    elseif fix == "rank"
        tensor_rank = fixed_target
        title     = "Fixed $tensor_rank rank tensor"
        xlabel    = "Tensor size"
    else
        error("invalid fix $fix")
    end

    #ylim_min = minimum(minimum.(values(error)))
    #ylim_max = maximum(maximum.(values(error)))
    ylim_min = -0.05
    ylim_max = +0.05
    for method in methods
        #if method == "TLrRR" && rel == true
        #    continue
        #end
        plot!(plt_time, horizontal_ax[2:end], error[method][2:end],
            title     = title,
            xlabel    = xlabel,
            ylabel    = ylabel,
            grid      = "on",
            legend    = legend,
            legendfont = legendfont,
            size       = img_size,

            #ylim = (ylim_min, ylim_max),

            label       = label_dict[method],
            linestyle   = linetype_dict[method],
            markershapes= markershapes_dict[method],
            linewidth   = linewidth,

            yguidefont  = fnt3,
            xguidefont  = fnt3,
            xtickfont   = fnt1,
            ytickfont   = fnt1,
            #yerror = error_std[method],
            #xticks = [2,4,6,8,10,12],
            xticks = ([20,100,300,1000], [20, 100, 300, 1000]),

            markersize  = markersize,
            markercolor = :transparent,
            markerstrokewidth = markerstrokewidth,
            markerstrokecolor = colors_dict[method],
            linecolor = colors_dict[method],
            )
    end
    fixval = get_tensorinfo(fixed_target)
    save_path = "../../png/synthetic/$cost"*"error_$fix$fixval.pdf"
    savefig(plt_time, save_path)
    println("$save_path has been saved")
end

function plot_time(horizontal_ax, runtime, runtime_std, fixed_target; fix, xlog)
    if fix == "size"
        tensor_size = fixed_target
        title     = "Random $tensor_size tensor"
        xlabel    = "Tensor rank"
    elseif fix == "rank"
        tensor_rank = fixed_target
        title     = "Fixed $tensor_rank rank tensor"
        xlabel    = "Tensor size"
    else
        error("invalid fix $fix")
    end

    plt_time = plot()
    if xlog
        plt_time = plot!(plt_time, xaxis=:log)
    end

    ylim_min = minimum(minimum.(values(runtime)))
    ylim_max = maximum(maximum.(values(runtime)))
    for method in methods
        plot!(plt_time, horizontal_ax[2:end], runtime[method][2:end],
            yaxis     = :log,
            title     = title,
            xlabel    = xlabel,
            ylabel    = "Running time(sec.)",
            grid      = "on",
            legend    = legend,
            legendfont = legendfont,
            size       = img_size,
            linewidth  = linewidth,

            ylim = (ylim_min, ylim_max),
            label       = label_dict[method],
            linestyle   = linetype_dict[method],
            markershapes= markershapes_dict[method],

            yguidefont  = fnt3,
            xguidefont  = fnt3,
            xtickfont   = fnt1,
            ytickfont   = fnt1,
            yerror = runtime_std[method],
            #xticks = ([20,100,300,1000], [20, 100, 300, 1000]),
            xticks = [2,4,6,8,10,12],

            markersize  = markersize,
            markercolor = :transparent,
            markerstrokewidth = markerstrokewidth,
            markerstrokecolor = colors_dict[method],
            linecolor = colors_dict[method],
             )
    end
    fixval = get_tensorinfo(fixed_target)
    save_path = "../../png/synthetic/time_$fix$fixval.pdf"
    savefig(plt_time, save_path)
    println("$save_path has been saved")
end

function plots_fix_size(tensor_size)
    horizontal_ax, KLerror, LSerror, runtime, KLerror_std, LSerror_std, runtime_std = get_data(tensor_size, fix="size")
    println("Data has been collected")
    plot_time(horizontal_ax, runtime, runtime_std, tensor_size,  fix="size", xlog=false)
    plot_error(horizontal_ax, LSerror, LSerror_std, tensor_size, fix="size", cost="LS", rel=false, xlog=false)
    plot_error(horizontal_ax, KLerror, KLerror_std, tensor_size, fix="size", cost="KL", rel=false, xlog=false)
end

function plots_fix_rank(reqrank)
    horizontal_ax, KLerror, LSerror, runtime, KLerror_std, LSerror_std, runtime_std = get_data(reqrank, fix="rank")
    println("Data has been collected")
    plot_time(horizontal_ax,  runtime, runtime_std, reqrank, fix="rank", xlog=true)
    plot_error(horizontal_ax, LSerror, LSerror_std, reqrank, fix="rank", cost="LS", rel=false, xlog=true)
    plot_error(horizontal_ax, KLerror, KLerror_std, reqrank, fix="rank", cost="KL", rel=false, xlog=true)
end

#plots_fix_size([25,25,25,25,25])
#plots_fix_size([30,30,30,30,30])
#plots_fix_size([200,200,200])
#plots_fix_rank([3,3,3,3,3])
plots_fix_rank([10,10,10])

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
# See examples/Finiteness_conjecture_counterexample.ipynb

using HybridSystems
using StaticArrays
using SwitchOnSafety

@testset "Finiteness conjecture counterexample with $lp_optimizer_constructor" for lp_optimizer_constructor in lp_optimizer_constructors
    A0 = @SMatrix [1 1; 0 1]
    A1 = @SMatrix [1 0; 1 1]
    A(α) = discreteswitchedsystem([A0, α * A1])
    s1 = A(1.0)

    s = A(0.8)
    s10 = SwitchOnSafety.periodicswitching(s, [2, 1])
    psw, done, polys = invariant_polytopes(s, lp_optimizer_constructor,
        s10, tol=0.0, verbose=0)
    @test psw == s10
    @test done
    psw, done, polys = invariant_polytopes(s, lp_optimizer_constructor,
        s10, tol=1e-15, verbose=0)
    @test psw == s10
    @test !done
    psw, done, polys = invariant_polytopes(s, lp_optimizer_constructor,
        s10, tol=1e-16, verbose=0)
    @test psw == s10
    @test done
    psw, done, polys = invariant_polytopes(s, lp_optimizer_constructor,
        s10, verbose=0)
    @test psw == s10
    @test done
end