text:
#=
    MCMC smapler is defined by a sequence of block updates and a master plan
    with a schedule on how to choose and modify those blocks
=#

#===============================================================================
                        Trackers of acceptance probability
===============================================================================#
"""
    mutable struct AccptTracker{S}

container to track the acceptance probability of the proposal
"""
mutable struct AccptTracker{S}
    accpt::S
    prop::S
    AccptTracker(::S) where {S <: Integer} = new{S}(zero(S), zero(S))
end

"""
    reset!(at::AccptTracker{S}) where S

Reset the acceptance probabliity of the proposal
"""
function reset!(at::AccptTracker{S}) where S
    at.accpt = zero(S)
    at.prop = zero(S)
end

"""
    reset!(at::Vector{<:AccptTracker})

Reset the acceptance probabliities
"""
function reset!(at::Vector{<:AccptTracker})
    for a in at
        reset!(a)
    end
end

"""
    acceptance_rate(at::AccptTracker{S}) where S

Compute the acceptance rate
"""
function acceptance_rate(at::AccptTracker{S}) where S
    at.prop == zero(S) && return 0.0
    at.accpt/at.prop
end
"""
    register_accpt!(at::AccptTracker, accepted)

Update acceptance probability tracker. `accepted` can be a vector
"""
function register_accpt!(at::AccptTracker{S}, accepted) where S
    at.prop += one(S)
    at.accpt += one(S)*accepted
end

function register_accpt!(at::Vector{AccptTracker{S}}, accepted::Vector) where S
    @assert length(at) == length(accepted)
    for i in 1:length(at)
        register_accpt!(at[i], accepted[i])
    end
end

function register_accpt!(at::Vector{AccptTracker{S}}, accepted::Bool) where S
    @assert length(at) == 1
    register_accpt!(at[1], accepted)
end

"""
    display(at::AccptTracker)

Show the acceptance rates
"""
function display(at::AccptTracker)
    print("Acceptance rate: ", round(accpetance_rate(at), digits=3), ".\n")
end

#===============================================================================
                        Definition of a single Gibbs step
===============================================================================#

abstract type AuxiliaryInfo end
"""
    struct UpdtAuxiliary{ST} <: AuxiliaryInfo

Info for auxiliary process. Tells which solver you use and whether to recompute the odes
"""
struct UpdtAuxiliary{ST} <: AuxiliaryInfo
    solver::ST
    recompute_ODEs::Bool
    function UpdtAuxiliary(solver::ST, recompute_ODEs::Bool) where ST<:ODESolverType
        new{ST}(solver, recompute_ODEs)
    end
end

"""
    struct ParamUpdate{UpdtType,...} <: MCMCUpdate

Needs to have a prior, transition kernel, acceptance rate tracker, information
about the which coordinates are being updated and an object with auxiliary,
use-specific information say whether to recompute_ODEs
"""
struct ParamUpdate{UpdtType,S,T,U,V} <: MCMCUpdate
    updt_coord::S
    t_kernel::T
    priors::U
    accpt_tracker::AccptTracker{Int64}
    accpt_history::Vector{Bool}
    aux::V
    readjust_param::ReadjustT
    readjust_history::Vector{Float64}

    function ParamUpdate(::UpdtType, updt_coord, θ, t_kernel::T, priors,
                         aux::V, readjust_param=(100, 0.1, 0.0, 999.9, 0.234, 50)
                         ) where {UpdtType<:ParamUpdateType,T,V<:AuxiliaryInfo}
        updt_coord = reformat_updt_coord(updt_coord, θ)
        priors = Priors(priors, updt_coord)
        S, U = typeof(updt_coord), typeof(priors)
        new{UpdtType,S,T,U,V}(updt_coord, t_kernel, priors, AccptTracker(0),
                              Bool[], aux, named_readjust(readjust_param),
                              [get_dispersion(t_kernel)])
    end

    function ParamUpdate(::UpdtType, updt_coord, θ, t_kernel::T, priors::U,
                         aux::V, readjust_param=(100, 0.1, 0.0, 999.9, 0.234, 50)
                         ) where {UpdtType<:ParamUpdateType,T,U<:Priors,V<:AuxiliaryInfo}
        updt_coord = reformat_updt_coord(updt_coord, θ)
        S = typeof(updt_coord)
        new{UpdtType,S,T,U,V}(updt_coord, t_kernel, priors, AccptTracker(0),
                              Bool[], aux, named_readjust(readjust_param),
                              [get_dispersion(t_kernel)])
    end
end

get_dispersion(::Any) = 0.0

IntContainer = Union{Number,NTuple{N,<:Integer},Vector{<:Integer}} where N

"""
    reformat_updt_coord(updt_coord::IntContainer, θ)

Single joint update of multiple parameters at once
"""
function reformat_updt_coord(updt_coord::S, θ) where S<:IntContainer
    @assert all([1 <= uc <= length(θ) for uc in updt_coord])
    Val{Tuple([i in updt_coord for i in 1:length(θ)])}()
end
#TODO
#WARNING the above might have been commented out whereas it is lines 27--29 in
#the file coordinate_access that should have been commented out...

"""
    readjust!(pu::ParamUpdate, cov_mat, mcmc_iter)

Readjust the proposal dispersion
"""
function readjust!(pu::ParamUpdate, cov_mat, mcmc_iter)
    at = pu.accpt_tracker
    p = pu.readjust_param
    disp = readjust!(pu.t_kernel, at, p, cov_mat, mcmc_iter, pu.updt_coord)
    push!(pu.readjust_history, disp)
    reset!(at)
end

readjust!(::Nothing, params...) = 0.0

aux_params(::Any, aux) = aux

function register_accpt!(pu::ParamUpdate, acc)
    register_accpt!(pu.accpt_tracker, acc)
    push!(pu.accpt_history, acc)
end
#===============================================================================
                        Definition of an imputation step
===============================================================================#

"""
    struct Imputation{B,ST,T} <: MCMCUpdate

Cointainer for various aspect of the imputation step
"""
struct Imputation{B,ST,T} <: MCMCUpdate
    accpt_tracker::Vector{AccptTracker{Int64}}
    accpt_history::Vector{T}
    ρs::Vector{Float64}
    blocking::B
    solver::ST
    readjust_param::ReadjustT
    ρs_history::Vector{Vector{Float64}}

    function Imputation(blocking::B, ρ::S, solver::ST=Ralston3(),
                            readjust_param=(100, 0.1, 0.00001, 0.99999, 0.234, 50)
                            ) where {B,S<:Number,ST<:ODESolverType}
        accpt_tracker = [AccptTracker(0) for _ in 1:length(blocking)]
        T = NTuple{length(blocking),Bool}
        accpt_history = T[]
        readjust_param = named_readjust(readjust_param)
        ρs = fill(ρ, length(blocking))
        new{B,ST,T}(accpt_tracker, accpt_history, ρs, blocking, solver,
                    readjust_param, [copy(ρs)])
    end
end
"""
    readjust!(pu::Imputation, ::Any, mcmc_iter)

Adaptively changes the  Crank-Nickolson parameter  `ρ`
"""
function readjust!(pu::Imputation, ::Any, mcmc_iter)
    at = pu.accpt_tracker
    p = pu.readjust_param
    δ = compute_δ(p, mcmc_iter)
    ar = [acceptance_rate(tracker) for tracker in at]
    for i in 1:length(pu.ρs)
        pu.ρs[i] = compute_ϵ(pu.ρs[i], p, ar[i], δ, -1.0, logit, sigmoid)
    end
    push!(pu.ρs_history, copy(pu.ρs))
    display_new_ρ(ar, pu.ρs)
    reset!(at)
end

function display_new_ρ(accpt_rates, ρs)
    print("imputation acceptance rate: ")
    for i in 1:length(accpt_rates) print(round(accpt_rates[i], digits=2), "  | ") end
    print("\nnew ρs: ")
    for i in 1:length(ρs) print(round(ρs[i], digits=3), "  | ") end
    print("\n")
end

aux_params(pu::Imputation, aux) = pu.blocking

function register_accpt!(pu::Imputation, acc)
    register_accpt!(pu.accpt_tracker, acc)
    push!(pu.accpt_history, tuple(acc...))
end

#===============================================================================
                        Overall schedule defining MCMC
===============================================================================#

"""
    MCMCSchedule{T}

Special iterator for the MCMC schedule.
Defines an iterator going step by step throught the update steps of the mcmc scheme
defined as numbered actions. Each group of sequential actions
is counted as single step for until the total `num_mcmc_steps` is reached.

Provides the functionality to change the update action after a certain amount of steps
through the `fuse` entry of the named tuple argument.

    MCMCSchedule(100, [[1], [2,3]], (save=10^3, verbose=10^4, warm_up=100,
                                readjust=(x->x%100==0), fuse=(x->false)))

In this example, first step 1 is performed, then steps 2, 3 after each other.
Will run 50 times step 1, 50 times step 2 followed by 3, saving the path every 10^3 times,
printing summaries every 10^4, burning the first 100 iteration
readjusting with the adaptive method every 100 iteration.
"""
struct MCMCSchedule{T}
    num_mcmc_steps::Int64
    updt_idx::Vector{Vector{Int64}}
    start::Int64
    actions::T

    function MCMCSchedule(num_mcmc_steps, updt_idx, actions)
        try
            names = (:save, :verbose, :warm_up, :readjust, :fuse)
            for name in names
                @assert name in keys(actions)
            end
        catch e
            if isa(e, AssertionError)
                @assert length(actions) == 5
                actions = (save=actions[1], verbose=actions[2], warm_up=actions[3],
                           readjust=actions[4], fuse=actions[5])

                print("`actions` has been passed without names. The names are ",
                      " set automatically:\n",
                      "save imputation every ", actions.save, " iterations,\n",
                      "print to console every ", actions.verbose, " iterations,\n",
                      "warm up the sampler for ", actions.warm_up, " iterations,\n",
                      "readjust proposals using function ", actions.readjust, "\n",
                      ",\nfuse the most correlated param updates at iterations ",
                      "determined by the function ", actions.fuse, ".\n")
            else
                rethrow(e)
            end
        end
        T = typeof(actions)
        new{T}(num_mcmc_steps, updt_idx, 1, actions)
    end
end

function Base.iterate(iter::MCMCSchedule, state=(iter.start, 1))
    element, i = state

    i > iter.num_mcmc_steps && return nothing

    at = iter.actions
    actions = (idx=iter.updt_idx[element],
               save=(i > at.warm_up && i % at.save == 0),
               verbose=(i % at.verbose == 0),
               param_updt=(i > at.warm_up),
               readjust=at.readjust(i),
               fuse=at.fuse(i),
               iter=i)

    return (actions, (transition(iter, element), i + 1))
end

Base.length(iter::MCMCSchedule) = iter.length

function Base.eltype(iter::MCMCSchedule)
    NamedTuple{(:idx, :save, :verbose, :param_updt, :readjust, :fuse, :iter),
               Tuple{Array{Int64,1},Bool,Bool,Bool,Bool,Bool,Int64}}
end
transition(schedule::MCMCSchedule, elem) = mod1(elem + 1, length(schedule.updt_idx))
"""
    reschedule!(schedule::MCMCSchedule, idx_to_remove, new_idx)
changes the schedule of the MCMC
"""
function reschedule!(schedule::MCMCSchedule, idx_to_remove, new_idx)
    for updt_indices in schedule.updt_idx
        updated = false
        N = length(updt_indices)
        for i in N:-1:1
            if updt_indices[i] in idx_to_remove
                if !updated
                    updt_indices[i] = new_idx
                    updated = true
                else
                    splice!(updt_indices, i)
                end
            end
        end
    end
end


text:
"""
    k,wk,Tkx = makescatteringtrans(M;n=4,ω=1.0,basis="Hermite",Nk=M,Na=M)

Construct transforms and associated arrays for precise numerical quadrature
evaluation of the scattering potential ``V_\\epsilon(x,t)`` of the energy-damped SPGPE.

 `M`: number of modes in the c-field.

 `n`: order of field product that will be integrated to machine precision.

 `ω`: trap frequency relative to the chosen reference frequency.

 `basis`: name of orthonormal eigenbasis representing c-field state. Currently `"Hermite"` is implemented.

 `Na`: number of auxiliary oscillator states.

 External links

 [Numerical method for the stochastic projected Gross-Pitaevskii equation, Physical Review E 89, 013302 (2014)](https://arxiv.org/abs/1310.0161)

"""

function makescatteringtrans(M;n=3,ω=1.0,basis="Hermite",Na=M)

    #n-field to x-space
    x,wx,Tx = nfieldtrans(M,n,ω=ω)

    #n-field to k-space. To include Nk need low level call to eigmat
    k,wk,Tk = nfieldtrans(M,n,ω=1/ω)

    #n-field aux transforms to x
    xa,wax,Tax = nfieldtrans(Na,n,ω=ω)

    #n-field aux transform to k
    ka,wak,Tak = nfieldtrans(Na,n,ω = 1/ω)
    Tak = Tak*diagm((-im).^(0:Na-1))

    #n-field transform from x --> k via auxiliary states.
    Txk = Tak*Tax'*diagm(wax)

return x,wx,Tx,k,wk,Txk
end

ext: jl


text:
module RatingCurves

using CSV, DataFrames, Distributions, Optim, Random, Statistics, StatsBase

import Base.sort, Base.minimum, Base.maximum, Statistics.var, StatsBase.bic

include("data.jl")
include("structures.jl")
include("parameterestimation.jl")

export

# Datatype
Gauging,
RatingCurve,
CompoundRatingCurve,

bic,
cint,
crcfit,
level,
discharge,
logdischarge,
minimum,
maximum,
pint,
pintlog,
rcfit,
sort

end


text:
using Flux
using Flux: Chain, Dense, params, train!, glorot_uniform
using Flux.Tracker: Params, gradient, update!

f = x -> (1/(0.82+x))^12 -2 * (1/(.82+x))^6
#f = x -> sin.(x.*7.0)

nsamples = 200
XS = rand(nsamples)
YS = f.(XS)
samples = zip(XS,YS)
#samples = [(x, f(x)) for x in XS]
#samples= [(XS, YS)]

xerr = range(0, 1, length=1000)
eval_nn(model, x) = model(x).data
maxfit_nn(model) = maximum(abs, f.(XS) - model(XS)')

struct Linear{F,S,T}
  W::S
  b::T
  σ::F
end

Linear(W, b) = Linear(W, b, identity)

function Linear(in::Integer, out::Integer, σ = identity)
  return Linear(param(randn(out, in)), param(randn(out)), σ)
end

@Flux.treelike Linear

function (a::Linear)(x::AbstractArray)
    W, b, σ = a.W, a.b, a.σ
    σ.(W * x' .+ b)
end

# Try to avoid hitting generic matmul in some simple cases
# Base's matmul is so slow that it's worth the extra conversion to hit BLAS
(a::Linear{<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  invoke(a, Tuple{AbstractArray}, x)

(a::Linear{<:Any,W})(x::AbstractArray{<:Real}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  a(T.(x))

function train_nn(nnodes, η)
    model = Chain(Linear(1, nnodes, Flux.tanh),   # first hidden layer
                  #Dense(nnodes, nnodes), # second hidden layer
                  Dense(nnodes, 1))      # output layer

    loss(x, y) = Flux.mse(model(x), y)

    θ = Flux.params(model)
    #println(θ)

    for i=1:5000
        #Flux.train!(loss, model_weights, samples, Descent(0.1))
        for j=1:length(XS)
            g = gradient(() -> loss([XS[j]], YS[j]), θ)
        #g = gradient(() -> loss(XS, YS), θ)
            for x in θ
                update!(x, -g[x]*η)
            end
        end
        if i%100==1
            @show loss(XS, YS)
            #println(size(model(XS)),size(YS))
            @show maxfit_nn(model)
        end
        if i%200==1
            η = η * 1.1
        end
    end
    return model
end

η = 0.000005 # Learning rate (should be smaller corresponding to the complex regression times the # of nodes)
NN_nodes = 5:5:15
errs_nn = Float64[]

for N in NN_nodes
    @show N
    nn = train_nn(N,η)
    println(Flux.params(nn))
    #@show e = err_nn(nn)
    #push!(errs_nn, e)
end
@show errs_nn

max_stars_repo_name: berkonat/ann-examples

max_issues_repo_name: berkonat/ann-examples

max_forks_repo_name: berkonat/ann-examples


text:
export BrickletIndustrialDual020mACurrentCallbackThreshold
struct BrickletIndustrialDual020mACurrentCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletIndustrialDual020mAIdentity
struct BrickletIndustrialDual020mAIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDual020mA
"""
Measures two DC currents between 0mA and 20mA (IEC 60381-1)
"""
mutable struct BrickletIndustrialDual020mA <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDual020mA(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_dual_0_20ma")
        deviceInternal = package.BrickletIndustrialDual020mA(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_current
"""
    $(SIGNATURES)

Returns the current of the specified sensor.

It is possible to detect if an IEC 60381-1 compatible sensor is connected
and if it works properly.

If the returned current is below 4mA, there is likely no sensor connected
or the sensor may be defect. If the returned current is over 20mA, there might
be a short circuit or the sensor may be defect.

If you want to get the current periodically, it is recommended to use the
:cb:`Current` callback and set the period with
:func:`Set Current Callback Period`.
"""
function get_current(device::BrickletIndustrialDual020mA, sensor)
    return device.deviceInternal.get_current(sensor)
end

export set_current_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Current` callback is triggered
periodically for the given sensor. A value of 0 turns the callback off.

The :cb:`Current` callback is only triggered if the current has changed since the
last triggering.
"""
function set_current_callback_period(device::BrickletIndustrialDual020mA, sensor, period)
    device.deviceInternal.set_current_callback_period(sensor, period)
end

export get_current_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Current Callback Period`.
"""
function get_current_callback_period(device::BrickletIndustrialDual020mA, sensor)
    return device.deviceInternal.get_current_callback_period(sensor)
end

export set_current_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Current Reached` callback for the given
sensor.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the current is *outside* the min and max values"
 "'i'",    "Callback is triggered when the current is *inside* the min and max values"
 "'<'",    "Callback is triggered when the current is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the current is greater than the min value (max is ignored)"
"""
function set_current_callback_threshold(device::BrickletIndustrialDual020mA, sensor, option, min, max)
    device.deviceInternal.set_current_callback_threshold(sensor, option, min, max)
end

export get_current_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Current Callback Threshold`.
"""
function get_current_callback_threshold(device::BrickletIndustrialDual020mA, sensor)
    return device.deviceInternal.get_current_callback_threshold(sensor)
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Current Reached`

is triggered, if the threshold

* :func:`Set Current Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletIndustrialDual020mA, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletIndustrialDual020mA)
    return device.deviceInternal.get_debounce_period()
end

export set_sample_rate
"""
    $(SIGNATURES)

Sets the sample rate to either 240, 60, 15 or 4 samples per second.
The resolution for the rates is 12, 14, 16 and 18 bit respectively.

.. csv-table::
 :header: "Value", "Description"
 :widths: 10, 100

 "0",    "240 samples per second, 12 bit resolution"
 "1",    "60 samples per second, 14 bit resolution"
 "2",    "15 samples per second, 16 bit resolution"
 "3",    "4 samples per second, 18 bit resolution"
"""
function set_sample_rate(device::BrickletIndustrialDual020mA, rate)
    device.deviceInternal.set_sample_rate(rate)
end

export get_sample_rate
"""
    $(SIGNATURES)

Returns the sample rate as set by :func:`Set Sample Rate`.
"""
function get_sample_rate(device::BrickletIndustrialDual020mA)
    return device.deviceInternal.get_sample_rate()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDual020mA)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDual020mA, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using OrthogonalPolynomialsQuasi, QuasiArrays, ContinuumArrays, BandedMatrices, LazyArrays, FastTransforms, ArrayLayouts, Test, FillArrays, Base64
import OrthogonalPolynomialsQuasi: Clenshaw, recurrencecoefficients, clenshaw, paddeddata, jacobimatrix
import LazyArrays: ApplyStyle
import QuasiArrays: MulQuasiMatrix
import Base: OneTo
import ContinuumArrays: MappedWeightedBasisLayout

@testset "Chebyshev" begin
    @testset "ChebyshevGrid" begin
        for kind in (1,2)
            @test all(ChebyshevGrid{kind}(10) .=== chebyshevpoints(Float64,10, Val(kind)))
            for T in (Float16, Float32, Float64)
                @test all(ChebyshevGrid{kind,T}(10) .=== chebyshevpoints(T,10, Val(kind)))
            end
            @test ChebyshevGrid{kind,BigFloat}(10) == chebyshevpoints(BigFloat,10, Val(kind))
        end
    end

    @testset "Evaluation" begin
        @testset "T" begin
            T = Chebyshev()
            @test @inferred(T[0.1,Base.OneTo(0)]) == Float64[]
            @test @inferred(T[0.1,Base.OneTo(1)]) == [1.0]
            @test @inferred(T[0.1,Base.OneTo(2)]) == [1.0,0.1]
            for N = 1:10
                @test @inferred(T[0.1,Base.OneTo(N)]) ≈ @inferred(T[0.1,1:N]) ≈ cos.((0:N-1)*acos(0.1))
                @test @inferred(T[0.1,N]) ≈ cos((N-1)*acos(0.1))
            end
            @test T[0.1,[2,5,10]] ≈ [0.1,cos(4acos(0.1)),cos(9acos(0.1))]

            @test axes(T[1:1,:]) === (Base.OneTo(1), Base.OneTo(∞))
            @test T[1:1,:][:,1:5] == ones(1,5)
            @test T[0.1,:][1:10] ≈ T[0.1,1:10] ≈ (T')[1:10,0.1]
        end

        @testset "U" begin
            U = ChebyshevU()
            @test @inferred(U[0.1,Base.OneTo(0)]) == Float64[]
            @test @inferred(U[0.1,Base.OneTo(1)]) == [1.0]
            @test @inferred(U[0.1,Base.OneTo(2)]) == [1.0,0.2]
            for N = 1:10
                @test @inferred(U[0.1,Base.OneTo(N)]) ≈ @inferred(U[0.1,1:N]) ≈ [sin((n+1)*acos(0.1))/sin(acos(0.1)) for n = 0:N-1]
                @test @inferred(U[0.1,N]) ≈ sin(N*acos(0.1))/sin(acos(0.1))
            end
            @test U[0.1,[2,5,10]] ≈ [0.2,sin(5acos(0.1))/sin(acos(0.1)),sin(10acos(0.1))/sin(acos(0.1))]
        end

        @testset "Function eval" begin
            T = Chebyshev()
            f = T * [1; 2; 3; zeros(∞)]
            @test f[0.1] ≈ 1 + 2*0.1 + 3*cos(2acos(0.1))
            x = [0.1,0.2]
            @test f[x] ≈ @.(1 + 2*x + 3*cos(2acos(x)))
        end
    end

    @testset "Transform" begin
        @testset "ChebyshevT" begin
            T = Chebyshev()
            @test T == ChebyshevT() == chebyshevt()
            Tn = @inferred(T[:,OneTo(100)])
            @test grid(Tn) == chebyshevpoints(100, Val(1))
            P = factorize(Tn)
            u = T*[P.plan * exp.(grid(Tn)); zeros(∞)]
            @test u[0.1] ≈ exp(0.1)

            # auto-transform
            x = axes(T,1)
            u = Tn * (P \ exp.(x))
            @test u[0.1] ≈ exp(0.1)

            u = Tn * (Tn \ exp.(x))
            @test u[0.1] ≈ exp(0.1)

            Tn = T[:,2:100]
            @test factorize(Tn) isa ContinuumArrays.ProjectionFactorization
            @test grid(Tn) == chebyshevpoints(100, Val(1))
            @test (Tn \ (exp.(x) .- 1.26606587775201)) ≈ (Tn \ u) ≈ (T\u)[2:100]

            u = T * (T \ exp.(x))
            @test u[0.1] ≈ exp(0.1)

            v = T[:,2:end] \ (exp.(x) .- 1.26606587775201)
            @test v[1:10] ≈ (T\u)[2:11]
        end
        @testset "ChebyshevU" begin
            U = ChebyshevU()
            @test U == chebyshevu()
            @test U ≠ ChebyshevT()
            x = axes(U,1)
            F = factorize(U[:,Base.OneTo(5)])
            @test @inferred(F \ x) ≈ [0,0.5,0,0,0]
            v = (x -> (3/20 + x + (2/5) * x^2)*exp(x)).(x)
            @inferred(U[:,Base.OneTo(5)]\v)
            w = @inferred(U\v)
            @test U[0.1,:]'w ≈ v[0.1]
        end
    end

    @testset "Mapped Chebyshev" begin
        x = Inclusion(0..1)
        T = Chebyshev()[2x .- 1,:]
        @test (T*(T\x))[0.1] ≈ 0.1
        @test (T* (T \ exp.(x)))[0.1] ≈ exp(0.1)
        @test chebyshevt(0..1) == T

        U = chebyshevu(0..1)
        @test (U*(U\x))[0.1] ≈ 0.1
        @test (U* (U \ exp.(x)))[0.1] ≈ exp(0.1)
    end

    @testset "weighted" begin
        T = ChebyshevT()
        w = ChebyshevTWeight()
        wT = WeightedChebyshevT()
        x = axes(wT,1)
        @test (x .* wT).args[2] isa BandedMatrix

        c = [1; 2; 3; zeros(∞)]
        a = wT * c;
        @test a[0.1] ≈ (T * c)[0.1]/sqrt(1-0.1^2)
        u = wT * (wT \ @.(exp(x)/sqrt(1-x^2)))
        @test u[0.1] ≈ exp(0.1)/sqrt(1-0.1^2)

        @testset "mapped" begin
            x = Inclusion(0..1)
            wT̃ = wT[2x .- 1, :]
            @test MemoryLayout(wT̃) isa MappedWeightedBasisLayout
            v = wT̃ * (wT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
            @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end

            WT̃ = w[2x .- 1] .* T[2x .- 1, :]
            @test MemoryLayout(wT̃) isa MappedWeightedBasisLayout
            v = wT̃ * (wT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
            @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end
        end
    end

    @testset "operators" begin
        T = ChebyshevT()
        U = ChebyshevU()
        @test axes(T) == axes(U) == (Inclusion(ChebyshevInterval()),Base.OneTo(∞))
        D = Derivative(axes(T,1))

        @test T\T === pinv(T)*T === Eye(∞)
        @test U\U === pinv(U)*U === Eye(∞)

        @test D*T isa MulQuasiMatrix
        D₀ = U\(D*T)
        @test D₀ isa BandedMatrix
        @test D₀[1:10,1:10] isa BandedMatrix{Float64}
        @test D₀[1:10,1:10] == diagm(1 => 1:9)
        @test colsupport(D₀,1) == 1:0

        S₀ = (U\T)[1:10,1:10]
        @test S₀ isa BandedMatrix{Float64}
        @test S₀ == diagm(0 => [1.0; fill(0.5,9)], 2=> fill(-0.5,8))

        x = axes(T,1)
        J = T\(x.*T)
        @test J isa BandedMatrix
        @test J[1:10,1:10] == jacobimatrix(T)[1:10,1:10]
    end

    @testset "test on functions" begin
        T = ChebyshevT()
        U = ChebyshevU()
        D = Derivative(axes(T,1))
        f = T*Vcat(randn(10), Zeros(∞))
        @test (U*(U\f)).args[1] isa Chebyshev{2}
        @test (U*(U\f))[0.1] ≈ f[0.1]
        @test (D*f)[0.1] ≈ ForwardDiff.derivative(x -> (ChebyshevT{eltype(x)}()*f.args[2])[x],0.1)
    end

    @testset "U->T lowering"  begin
        wT = ChebyshevWeight() .* Chebyshev()
        wU = ChebyshevUWeight() .*  ChebyshevU()
        @test (wT \ wU)[1:10,1:10] == diagm(0 => fill(0.5,10), -2 => fill(-0.5,8))
    end

    @testset "sub-of-sub" begin
        T = Chebyshev()
        V = T[:,2:end]
        @test copy(V) == V
        @test view(V,0.1:0.1:1,:) isa SubArray
        @test V[0.1:0.1:1,:] isa SubArray
        @test V[0.1:0.1:1,:][:,1:5] == T[0.1:0.1:1,2:6]
        @test parentindices(V[:,OneTo(5)])[1] isa Inclusion
    end

    @testset "==" begin
        @test Chebyshev() == ChebyshevT() == ChebyshevT{Float32}()
        @test ChebyshevU() == ChebyshevU{Float32}()
        @test Chebyshev{3}() == Chebyshev{3,Float32}()
        @test Chebyshev() ≠ ChebyshevU()
    end

    @testset "sum" begin
        wT = ChebyshevWeight() .* Chebyshev()
        @test sum(wT; dims=1)[1,1:10] == [π; zeros(9)]
        @test sum(wT * [[1,2,3]; zeros(∞)]) == 1.0π
        wU = ChebyshevUWeight() .* ChebyshevU()
        @test sum(wU; dims=1)[1,1:10] == [π/2; zeros(9)]
        @test sum(wU * [[1,2,3]; zeros(∞)]) == π/2

        x = Inclusion(0..1)
        @test sum(wT[2x .- 1, :]; dims=1)[1,1:10] == [π/2; zeros(9)]
        @test sum(wT[2x .- 1, :] * [[1,2,3]; zeros(∞)]) == π/2
    end

    @testset "algebra" begin
        T = ChebyshevT()
        U = ChebyshevU()

        a = T * [1; zeros(∞)]
        b = T * [1; 2; 3; zeros(∞)]
        c = U * [4; 5; zeros(∞)]

        @test a + b == b + a == T * [2;  2;  3; zeros(∞)]
        @test a - b == T * [0; -2; -3; zeros(∞)]
        @test a + c == c + a == U * [5; 5; zeros(∞)]
    end

    @testset "multiplication" begin
        U = ChebyshevU()
        x = axes(U,1)
        a = U * [[0.25; 0.5; 0.1]; Zeros(∞)];
        M = Clenshaw(a, U)

        @test colsupport(M,4:10) == rowsupport(M,4:10) == 2:12
        @test colsupport(M',4:10) == rowsupport(M',4:10) == 2:12

        A,B,C = recurrencecoefficients(U)
        c = paddeddata(a.args[2])
        m = 10
        X = (U \ (x .* U))[1:m,1:m]
        @test clenshaw(c, A, B, C, X) ≈ 3/20*I + X + (2/5) * X^2

        @test bandwidths(M) == (2,2)
        @test M[2:10,2:10] == M[1:10,1:10][2:10,2:10]
        @test M[3:10,5:20] == M[1:10,1:20][3:10,5:20]
        @test M[100_000:101_000,100_000:101_000] == M[2:1002,2:1002]

        @test Eye(∞) * M isa Clenshaw

        M = U \ (a .* U)
        @test M isa Clenshaw

        f = U \ exp.(x)
        @test M*f ≈ U \ (x -> (3/20 + x + (2/5) * x^2)*exp(x)).(x)

        @testset "weighted" begin
            T = ChebyshevT()
            wT = WeightedChebyshevT()
            a = wT * [1; 2; 3; zeros(∞)];
            @test (a .* T)[0.1,1:10] ≈ a[0.1] * T[0.1,1:10]
        end
    end

    @testset "show" begin
        T = Chebyshev()
        @test stringmime("text/plain", T * [1; 2; Zeros(∞)]) == "Chebyshev{1,Float64} * [1.0, 2.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, …]"
    end

    @testset "Complex eltype" begin
        @test axes(ChebyshevT{ComplexF64}(),1) ≡ Inclusion{ComplexF64}(ChebyshevInterval())
        @test ChebyshevT{ComplexF64}()[0.1+0im,1:5] isa Vector{ComplexF64}
        @test ChebyshevT{ComplexF64}()[0.1+0im,1:5] ≈ ChebyshevT()[0.1,1:5]
    end

    @testset "extrapolation" begin
        @test Base.unsafe_getindex(Chebyshev(), 5.0, Base.OneTo(5)) ≈ [cos(k*acos(5+0im)) for k=0:4]
        @test Base.unsafe_getindex(Chebyshev(), 5.0, 5)  ≈ cos(4*acos(5+0im))
        @test Base.unsafe_getindex(ChebyshevT{ComplexF64}(), 5.0+im, 5)  ≈ cos(4*acos(5+im))
        @test Base.unsafe_getindex(Chebyshev(), 5.0, 1:5) ≈ [cos(n*acos(5+0im)) for n=0:4]
        @test Base.unsafe_getindex(Chebyshev(), 5.0, Base.OneTo(5)) ≈ [cos(n*acos(5+0im)) for n=0:4]
        @test Base.unsafe_getindex(Chebyshev(), [1.0,5.0], 5)  ≈ [cos(4*acos(x+0im)) for x in [1,5]]
        @test Base.unsafe_getindex(Chebyshev(), [1.0,5.0], 1:5)  ≈ [cos(n*acos(x+0im)) for x in [1,5], n in 0:4]
        @test Base.unsafe_getindex(Chebyshev(), 5.0, [1,5])  ≈ [cos(n*acos(5.0+0im)) for n in [0,4]]
        @test Base.unsafe_getindex(Chebyshev(), [1.0,5.0], [1,5])  ≈ [cos(n*acos(x+0im)) for x in [1,5], n in [0,4]]
        @test Base.unsafe_getindex(Chebyshev(), [1.0,5.0], Base.OneTo(5))  ≈ [cos(n*acos(x+0im)) for x in [1,5], n in 0:4]

        v = Base.unsafe_getindex(Chebyshev(), 2, 2:∞)
        @test eltype(v) == Float64
        @test v[5] == Base.unsafe_getindex(Chebyshev(), 2, 6)
    end

    @testset "Christoffel–Darboux" begin
        T = Chebyshev()
        X = T\ (axes(T,1) .* T)
        Mi = inv(T'*(ChebyshevWeight() .* T))
        x,y = 0.1,0.2
        n = 10
        Pn = Diagonal([Ones(n); Zeros(∞)])
        Min = Pn * Mi
        @test (X*Min - Min*X')[1:n,1:n] ≈ zeros(n,n)
        @test (x-y) * T[x,1:n]'Mi[1:n,1:n]*T[y,1:n] ≈ T[x,n:n+1]' * (X*Min - Min*X')[n:n+1,n:n+1] * T[y,n:n+1]

        @testset "extrapolation" begin
            x,y = 0.1,3.4
            @test (x-y) * T[x,1:n]'Mi[1:n,1:n]*Base.unsafe_getindex(T,y,1:n) ≈ T[x,n:n+1]' * (X*Min - Min*X')[n:n+1,n:n+1] * Base.unsafe_getindex(T,y,n:n+1)
        end
    end

    @testset "pointwise evaluation" begin
        T = Chebyshev()
        v = T[0.1,:]
        @test v[1:100_000] == T[0.1,1:100_000]
        @test Base.BroadcastStyle(typeof(v)) isa LazyArrays.LazyArrayStyle{1}
        @test (v./v)[1:10] == ones(10)
        v = T[0.1,2:∞]
        @test v[1:100_000] == T[0.1,2:100_001]
        @test Base.BroadcastStyle(typeof(v)) isa LazyArrays.LazyArrayStyle{1}
        @test (v./v)[1:10] == ones(10)
        v = Base.unsafe_getindex(T, 2.0, :)
        @test isequal(v[1:10_000], Base.unsafe_getindex(T, 2.0, 1:10_000))
        @test Base.BroadcastStyle(typeof(v)) isa LazyArrays.LazyArrayStyle{1}
        @test (v./v)[1:10] == ones(10)
    end
end


text:
"""
    space(::SiteType"S=1";
          conserve_qns = false,
          conserve_sz = conserve_qns,
          qnname_sz = "Sz")

Create the Hilbert space for a site of type "S=1".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function space(
  ::SiteType"S=1"; conserve_qns=false, conserve_sz=conserve_qns, qnname_sz="Sz"
)
  if conserve_sz
    return [QN(qnname_sz, +2) => 1, QN(qnname_sz, 0) => 1, QN(qnname_sz, -2) => 1]
  end
  return 3
end

state(::SiteType"S=1", ::StateName"Up") = 1
state(::SiteType"S=1", ::StateName"Z0") = 2
state(::SiteType"S=1", ::StateName"Dn") = 3

state(st::SiteType"S=1", ::StateName"↑") = state(st, StateName("Up"))
state(st::SiteType"S=1", ::StateName"0") = state(st, StateName("Z0"))
state(st::SiteType"S=1", ::StateName"↓") = state(st, StateName("Dn"))

function op!(Op::ITensor, ::OpName"Sz", ::SiteType"S=1", s::Index)
  Op[s' => 1, s => 1] = +1.0
  return Op[s' => 3, s => 3] = -1.0
end

op!(Op::ITensor, ::OpName"Sᶻ", t::SiteType"S=1", s::Index) = op!(Op, OpName("Sz"), t, s)

function op!(Op::ITensor, ::OpName"S+", ::SiteType"S=1", s::Index)
  Op[s' => 2, s => 3] = sqrt(2)
  return Op[s' => 1, s => 2] = sqrt(2)
end

op!(Op::ITensor, ::OpName"S⁺", t::SiteType"S=1", s::Index) = op!(Op, OpName("S+"), t, s)

op!(Op::ITensor, ::OpName"Splus", t::SiteType"S=1", s::Index) = op!(Op, OpName("S+"), t, s)

function op!(Op::ITensor, ::OpName"S-", ::SiteType"S=1", s::Index)
  Op[s' => 3, s => 2] = sqrt(2)
  return Op[s' => 2, s => 1] = sqrt(2)
end

op!(Op::ITensor, ::OpName"S⁻", t::SiteType"S=1", s::Index) = op!(Op, OpName("S-"), t, s)

op!(Op::ITensor, ::OpName"Sminus", t::SiteType"S=1", s::Index) = op!(Op, OpName("S-"), t, s)

function op!(Op::ITensor, ::OpName"Sx", ::SiteType"S=1", s::Index)
  Op[s' => 2, s => 1] = 1 / sqrt(2)
  Op[s' => 1, s => 2] = 1 / sqrt(2)
  Op[s' => 3, s => 2] = 1 / sqrt(2)
  return Op[s' => 2, s => 3] = 1 / sqrt(2)
end

op!(Op::ITensor, ::OpName"Sˣ", t::SiteType"S=1", s::Index) = op!(Op, OpName("Sx"), t, s)

function op!(Op::ITensor, ::OpName"iSy", ::SiteType"S=1", s::Index)
  Op[s' => 2, s => 1] = -1 / sqrt(2)
  Op[s' => 1, s => 2] = +1 / sqrt(2)
  Op[s' => 3, s => 2] = -1 / sqrt(2)
  return Op[s' => 2, s => 3] = +1 / sqrt(2)
end

op!(Op::ITensor, ::OpName"iSʸ", t::SiteType"S=1", s::Index) = op!(Op, OpName("iSy"), t, s)

function op!(Op::ITensor, ::OpName"Sy", ::SiteType"S=1", s::Index)
  Op[s' => 2, s => 1] = -1im / sqrt(2)
  Op[s' => 1, s => 2] = +1im / sqrt(2)
  Op[s' => 3, s => 2] = -1im / sqrt(2)
  return Op[s' => 2, s => 3] = +1im / sqrt(2)
end

op!(Op::ITensor, ::OpName"Sʸ", t::SiteType"S=1", s::Index) = op!(Op, OpName("Sy"), t, s)

function op!(Op::ITensor, ::OpName"Sz2", ::SiteType"S=1", s::Index)
  Op[s' => 1, s => 1] = +1.0
  return Op[s' => 3, s => 3] = +1.0
end

function op!(Op::ITensor, ::OpName"Sx2", ::SiteType"S=1", s::Index)
  Op[s' => 1, s => 1] = 0.5
  Op[s' => 3, s => 1] = 0.5
  Op[s' => 2, s => 2] = 1.0
  Op[s' => 1, s => 3] = 0.5
  return Op[s' => 3, s => 3] = 0.5
end

function op!(Op::ITensor, ::OpName"Sy2", ::SiteType"S=1", s::Index)
  Op[s' => 1, s => 1] = +0.5
  Op[s' => 3, s => 1] = -0.5
  Op[s' => 2, s => 2] = +1.0
  Op[s' => 1, s => 3] = -0.5
  return Op[s' => 3, s => 3] = +0.5
end

space(::SiteType"SpinOne"; kwargs...) = space(SiteType("S=1"); kwargs...)

state(::SiteType"SpinOne", st::AbstractString) = state(SiteType("S=1"), st)

op!(Op::ITensor, o::OpName, ::SiteType"SpinOne", s::Index) = op!(Op, o, SiteType("S=1"), s)


text:
"""
Star fileformat

Output dictionary with for each Airport a list of STAR names.

See EUROCONTROL NEST Manual for Star fileformat description

Star fileformat example
BIKF THORI METIL ASTAM KEF CELLO
BIRK THORI METIL ASTAM KEF CELLO
DTMB MON.A
DTNH NBA.A
"""
module Star

export read

function read(file)
    airportstars = Dict{String, Vector{AbstractString}}()
    for line in eachline(file)
        line_elements = split(line, ' ')
        airport_name = line_elements[1]
        airportstars[airport_name] = [star for star in line_elements[2:end]]
    end
    return airportstars
end

end # module


text:
# This script is a minimal example of how to use a `FeedbackChain` to build a
# convolutional network with feedback operations and to apply it to data.
using FeedbackNets
using Flux

# some parameters
batchsize = 10
usegpu = true

if usegpu
    using CuArrays
end

# We generate the FeedbackChain much like a normal Flux.Chain. The difference is
# that we set Splitters where feedback branches off and Mergers where the feedback
# is folded back in.
# The basis for this architecture is LeNet5, but with ReLUs instead of sigmoids.
chain = FeedbackChain(
    Conv((5,5), 1=>6, relu),
    Merger("conv2", ConvTranspose((10,10), 16=>6, relu, stride=2), +),
    MaxPool((2,2), stride=(2,2)),
    Conv((5,5), 6=>16, relu),
    Splitter("conv2"),
    Merger("fc1", ConvTranspose((10,10), 120=>16, relu), +),
    MaxPool((2,2), stride=(2,2)),
    Conv((5,5), 16=>120, relu),
    Splitter("fc1"),
    x -> reshape(x, 120, size(x, 4)),
    Dense(120, 84, relu),
    Dense(84, 10)
)

# In order to apply this chain, we need to pass it a dictionary with the initial
# states of all Splitters. Here, we simply set all initial states to zero.
h = Dict(
    "conv2" => zeros(10, 10, 16, batchsize),
    "fc1" => zeros(1, 1, 120, batchsize)
)

# If we want to calculate the model on a GPU, we also need to shift the initial
# states.
if usegpu
    chain = chain |> gpu
    h = Dict(key => gpu(val) for (key, val) in pairs(h))
end

# Now we can turn out model into a Flux.Recur, which handles its state
# internally.
model = Flux.Recur(chain, h)

# At this point, we can treat out model like any other recurrent network in Flux.
# We can apply it to a single input batch ...
batch = randn(32, 32, 1, batchsize)
if usegpu
    batch = gpu(batch)
end
out1 = model(batch)
# ... or broadcast it over a (time)sequence of batches.
seq = [batch, batch, batch]
out2 = model.(seq)


text:
using Test
using HMatrices
using LinearAlgebra
using Random
using StaticArrays
using HMatrices: RkMatrix

include(joinpath(HMatrices.PROJECT_ROOT,"test","testutils.jl"))

N,r     = 500,3
atol    = 1e-6
X  = Y     = rand(SVector{3,Float64},N)
splitter = CardinalitySplitter(nmax=50)
Xclt = Yclt  = ClusterTree(X,splitter)
adm       = StrongAdmissibilityStd(eta=3)
rtol      = 1e-5
comp      = PartialACA(rtol=rtol)
K         = helmholtz_matrix(X,Y,1.3)
H         = assemble_hmat(K,Xclt,Yclt;adm,comp,threads=false,distributed=false)
T         = eltype(H)
_H   = Matrix(H;global_index=false)
R    = RkMatrix(rand(T,N,r),rand(T,N,r))
_R   = Matrix(R)
M    = rand(T,N,N)
_M   = Matrix(M)
a    = rand()

@testset "axpy!" begin
    @test axpy!(a,M,deepcopy(R)) ≈ a*_M + _R
    tmp = axpy!(a,M,deepcopy(H))
    @test Matrix(tmp;global_index=false) ≈ a*_M + _H
    @test axpy!(a,R,deepcopy(M)) ≈ a*_R + _M
    @test axpy!(a,R,deepcopy(R)) ≈ a*_R + _R
    tmp = axpy!(a,R,deepcopy(H))
    Matrix(tmp;global_index=false) ≈ a*_R + _H
    @test axpy!(a,H,deepcopy(M)) ≈ a*_H + _M
    @test axpy!(a,H,deepcopy(R)) ≈ a*_H + _R
    tmp = axpy!(a,H,deepcopy(H))
    @test Matrix(tmp;global_index=false) ≈ a*_H + _H
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
# # Granger method vs. Phase-Slope-Index
#
# Here we attempt at reproducing the *Figure 1* of [Nolte et al. 2008](http://link.aps.org/abstract/PRL/v100/e234101)
#
#md # The notebook can be viewed here:
#md # * [![binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/benchmark.ipynb)
#md # * [![nbviewer](https://img.shields.io/badge/render-nbviewer-orange.svg)](@__NBVIEWER_ROOT_URL__/generated/benchmark.ipynb)
#
#nb # ### Acknowledgement
#nb # This work was funded by the German Federal Ministry of Education and Research ([BMBF](https://www.bmbf.de/)) in the project ALICE III under grant ref. 01IS18049B.
#
#nb # ### Reference
# *Robustly Estimating the Flow Direction of Information in Complex Physical Systems* paper, by *Guido Nolte, Andreas Ziehe, Vadim V. Nikulin, Alois Schlögl, Nicole Krämer, Tom Brismar, and Klaus-Robert Müller*, ([Nolte et al. 2008](http://link.aps.org/abstract/PRL/v100/e234101)).
#
# ### Load packages

using TimeSeriesCausality
using Distributions: MvNormal
using Plots: plot, bar

# ### Unidirectional flux with white noise
# - channel 1 is the causal (independent) time series
# - channel 2 is the effect (dependent) time series

# parameters
order = 1
n_samples = 1024  # number of data points per channel
segment_length = 128
epoch_length = 128
noise_scale = 0.25
plot_slice=128

## data generation
rand_data = randn(Float64, (n_samples + order))
signal = zeros(Float64, (n_samples, 2))
white_noise = noise_scale * randn(n_samples, 2)
signal[:, 1] = rand_data[order+1:end]
signal[:, 2] = rand_data[1:n_samples]
signal += white_noise

## Causal estimation
psi, psi_std = psi_est(signal, segment_length; method="jackknife", eplen=epoch_length)
granger_idx, err_std = granger_est(signal, segment_length; method="jackknife")

p1 = plot(
    1:plot_slice,
    signal[1:plot_slice, :];
    title="Unidirectional Flux with white noise",
    label=["Channel 1" "Channel 2"],
    linecolor=["red" "blue"],
    xlabel="Time [bins]",
    lw = 2,
);
p2 = bar(
    ["Granger" "PSI"],
    [psi[1, 2] granger_idx];
    title="Granger vs PSI",
    yerror=[psi_std[1, 2] err_std],
    legend=false
);
plot(p1, p2; layout=(1, 2), size=(800, 300))

# ### Pink (correlated) noise

## data generation
signal = zeros(Float64, (n_samples, 2))
signal[:, 1] = TimeSeriesCausality.pink_noise(n_samples, 1.0)
signal[:, 2] = TimeSeriesCausality.pink_noise(n_samples, 1.0)

## Causal estimation
psi, psi_std = psi_est(signal, segment_length; method="jackknife", eplen=epoch_length)
granger_idx, err_std = granger_est(signal, segment_length; method="jackknife")

p1 = plot(
    1:plot_slice,
    signal[1:plot_slice, :];
    title="Correlated noise",
    label=["Channel 1" "Channel 2"],
    linecolor=["red" "blue"],
    xlabel="Time [bins]",
    lw = 2,
);
p2 = bar(
    ["Granger" "PSI"],
    [psi[1, 2] granger_idx];
    title="Granger vs PSI",
    yerror=[psi_std[1, 2] err_std],
    legend=false
);
plot(p1, p2; layout=(1, 2), size=(800, 300))


text:
using DataFrames
using JFSM2
using PyPlot
using IcanProj
using CSV
using Statistics

function load_results()

    res_all = Dict()

    for spaceres in ["5km", "10km", "25km", "50km"]

        filename = joinpath(dirname(pathof(IcanProj)), "..", "data", "table_errors_$(spaceres).txt")

        res_all[spaceres] = CSV.File(filename, delim = ",") |> DataFrame

    end

    return res_all

end

function error_matrix(res_all, variable, measure)

    spaceres = ["5km", "10km", "25km", "50km"]

    cfgs = 1:32

    data = fill(0.0, length(cfgs), length(spaceres))

    for j in eachindex(spaceres)

        df_res = res_all[spaceres[j]]

        for i in cfgs

            tmp = df_res[Symbol("$(measure)_$(variable)_cfg$(cfgs[i])")]

            tmp = tmp[map(x -> !isnan(x), tmp)]

            data[i,j] = median(tmp) #df_res[Symbol("$(measure)_$(variable)_cfg$(cfgs[i])")])

        end

    end

    return data, spaceres

end

function plot_error_scales(data, df_cfg, plottitle, ylimits)

    data = permutedims(data)

    exchng_off = convert(Array{Bool}, df_cfg[:exchng] .== 0)
    exchng_on  = convert(Array{Bool}, df_cfg[:exchng] .== 1)

    figure(figsize = (5, 4))

    plot(collect(1:size(data,1)), data, color = "gray")

    fill_between(collect(1:size(data,1)),
                 maximum(data[:, exchng_off], dims = 2)[:],
                 minimum(data[:, exchng_off], dims = 2)[:],
                 facecolor = "red", edgecolor = "red", alpha = 0.5,
                 label = "Exchng = 0")

    fill_between(collect(1:size(data,1)),
                 maximum(data[:, exchng_on], dims = 2)[:],
                 minimum(data[:, exchng_on], dims = 2)[:],
                 facecolor = "blue", edgecolor = "blue", alpha = 0.5,
                 label = "Exchng = 1")

    xticks(collect(1:size(data,1)), ["5km", "10km", "25km", "50km"])

    title(plottitle)

    # ylim(ylimits)

    # ylabel("NSE (-)")

    legend()

end

# Plot results

res_all = load_results()

df_cfg = cfg_table()

figpath = joinpath(dirname(pathof(IcanProj)), "..", "plots", "error_scales")

tmp = [("swe", "SWE"),
       ("snowdepth", "Snowdepth"),
       ("rnet", "Net radiation")]

for (v, t) in tmp

    global nseres, spaceres

    nseres, spaceres = error_matrix(res_all, v, "bias")

    plot_error_scales(nseres, df_cfg, t, (0.95, 1.0))

    savefig(joinpath(figpath, "$(v)_error.png"), dpi = 200)

end

tmp = [("hatmo", "Sensible heat flux"),
       ("latmo", "Latent heat flux"),
       ("melt", "Surface melt")]

for (v, t) in tmp

    global nseres, spaceres

    nseres, spaceres = error_matrix(res_all, v, "bias")

    plot_error_scales(nseres, df_cfg, t, (0.3, 1.0))

    savefig(joinpath(figpath, "$(v)_error.png"), dpi = 200)

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMIFlux

@debug "Debugging messages enabled for FMIFlux ..."

using Requires

using FMIImport
using FMIImport: fmi2ValueReference, FMU, FMU2, FMU2Component

# ToDo: this can not be imported correctly from FMIImport.jl
fmi2Struct = Union{FMU2, FMU2Component}

using FMIImport: fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Reset, fmi2Terminate
using FMIImport: fmi2NewDiscreteStates, fmi2SetContinuousStates, fmi2GetContinuousStates, fmi2GetNominalsOfContinuousStates
using FMIImport: fmi2SetTime, fmi2CompletedIntegratorStep, fmi2GetEventIndicators, fmi2GetDerivatives, fmi2GetReal
using FMIImport: fmi2SampleDirectionalDerivative, fmi2GetDirectionalDerivative, fmi2GetJacobian, fmi2GetJacobian!
using FMIImport: fmi2True, fmi2False

include("component_shadow.jl")
include("FMI2_neural.jl")
include("FMI_neural.jl")
include("misc.jl")

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        import .Plots

        include("FMI_plot.jl")

        Plots.plot(nfmu::NeuralFMU) = fmiPlot(nfmu)
    end

    @require FMI="14a09403-18e3-468f-ad8a-74f8dda2d9ac" begin
        @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            import .FMI

            FMI.fmiPlot(nfmu::NeuralFMU) = fmiPlot(nfmu)
        end
    end
end

# FMI2_neural.jl
export fmi2EvaluateME, fmi2DoStepCS
export fmi2InputDoStepCSOutput
export ME_NeuralFMU, CS_NeuralFMU, NeuralFMU, NeuralFMU_TrainingModeConfig

# FMI_neural.jl
export fmiEvaluateME, fmiDoStepCS
export fmiInputDoStepCSOutput

# misc.jl
export mse_interpolate, transferParams!, transferFlatParams!, lin_interp

# debugging only
# export _build_jac_dx_x_slow

end # module


text:
using LabJack

lsdev()

dev = open(dtT7, ctANY)

write_digital("EIO0", 0, dev)

write_analog("TDAC1", 5.0, dev)

read_digital("EIO0")


text:
println("################################################################################")
# https://docs.julialang.org/en/v1/stdlib/Test/#Basic-Unit-Tests
using Test, Collatz

# Test function collatz_function(n::Integer; P::Integer=2, a::Integer=3, b::Integer=1)
@testset verbose = true "collatz_function" begin
    # Default/Any (P,a,b); 0 trap
    @test Collatz.collatz_function(0) == 0
    # Default/Any (P,a,b); 1 cycle; positives
    @test Collatz.collatz_function(1) == 4
    @test Collatz.collatz_function(4) == 2
    @test Collatz.collatz_function(2) == 1
    # Default/Any (P,a,b); -1 cycle; negatives
    @test Collatz.collatz_function(-1) == -2
    @test Collatz.collatz_function(-2) == -1
    # Test a wider modulo sweep by upping P to 5, a to 2, and b to 3.
    @test Collatz.collatz_function(1, P=5, a=2, b=3) == 5
    @test Collatz.collatz_function(2, P=5, a=2, b=3) == 7
    @test Collatz.collatz_function(3, P=5, a=2, b=3) == 9
    @test Collatz.collatz_function(4, P=5, a=2, b=3) == 11
    @test Collatz.collatz_function(5, P=5, a=2, b=3) == 1
    # Test negative P, a and b. %, used in the function, is "floor" in python
    # rather than the more reasonable euclidean, but we only use it's (0 mod P)
    # conjugacy class to determine functionality, so the flooring for negative P
    # doesn't cause any issue.
    @test Collatz.collatz_function(1, P=-3, a=-2, b=-5) == -7
    @test Collatz.collatz_function(2, P=-3, a=-2, b=-5) == -9
    @test Collatz.collatz_function(3, P=-3, a=-2, b=-5) == -1
    # Set P and a to 0 to @test on __assert_sane_parameterisation
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.collatz_function(1, P=0, a=2, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.collatz_function(1, P=0, a=0, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_A) Collatz.collatz_function(1, P=1, a=0, b=3)
end

# Test function reverse_collatz_function(n::Integer; P::Integer=2, a::Integer=3, b::Integer=1)
@testset verbose = true "reverse_collatz_function" begin
    # Default (P,a,b); 0 trap [as b is not a multiple of a]
    @test Collatz.reverse_collatz_function(0) == [0]
    # Default (P,a,b); 1 cycle; positives
    @test Collatz.reverse_collatz_function(1) == [2]
    @test Collatz.reverse_collatz_function(4) == [1, 8]
    @test Collatz.reverse_collatz_function(2) == [4]
    # Default (P,a,b); -1 cycle; negatives
    @test Collatz.reverse_collatz_function(-1) == [-2]
    @test Collatz.reverse_collatz_function(-2) == [-4, -1]
    # Test a wider modulo sweep by upping P to 5, a to 2, and b to 3.
    @test Collatz.reverse_collatz_function(1, P=5, a=2, b=3) == [-1, 5]
    @test Collatz.reverse_collatz_function(2, P=5, a=2, b=3) == [10]
    @test Collatz.reverse_collatz_function(3, P=5, a=2, b=3) == [15]  # also tests !0
    @test Collatz.reverse_collatz_function(4, P=5, a=2, b=3) == [20]
    @test Collatz.reverse_collatz_function(5, P=5, a=2, b=3) == [1, 25]
    # Test negative P, a and b. %, used in the function, is "floor" in python
    # rather than the more reasonable euclidean, but we only use it's (0 mod P)
    # conjugacy class to determine functionality, so the flooring for negative P
    # doesn't cause any issue.
    @test Collatz.reverse_collatz_function(1, P=-3, a=-2, b=-5) == [-3]  # != [-3, -3]
    @test Collatz.reverse_collatz_function(2, P=-3, a=-2, b=-5) == [-6]
    @test Collatz.reverse_collatz_function(3, P=-3, a=-2, b=-5) == [-9, -4]
    # Set P and a to 0 to @test on __assert_sane_parameterisation
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.reverse_collatz_function(1, P=0, a=2, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.reverse_collatz_function(1, P=0, a=0, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_A) Collatz.reverse_collatz_function(1, P=1, a=0, b=3)
    # If b is a multiple of a, but not of Pa, then 0 can have a reverse.
    @test Collatz.reverse_collatz_function(0, P=17, a=2, b=-6) == [0, 3]
    @test Collatz.reverse_collatz_function(0, P=17, a=2, b=102) == [0]
end

# Test function hailstone_sequence(initial_value::Integer; P::Integer=2, a::Integer=3, b::Integer=1, max_total_stopping_time::Integer=1000, total_stopping_time::Bool=true, verbose::Bool=true)
@testset verbose = true "hailstone_sequence" begin
    # Test 0's immediated termination.
    @test Collatz.hailstone_sequence(0) == [[_CC.ZERO_STOP, 0]]
    @test Collatz.hailstone_sequence(0, verbose=false) == [0]
    # The cycle containing 1 wont yield a cycle termination, as 1 is considered
    # the "total stop" that is the special case termination.
    @test Collatz.hailstone_sequence(1) == [[_CC.TOTAL_STOPPING_TIME, 0]]
    @test Collatz.hailstone_sequence(1, verbose=false) == [1]
    # Test the 3 known default parameter's cycles (ignoring [1,4,2])
    # If not verbose, then the result will be the cycle plus the final value
    # being the start of the cycle. If verbose, then the output will be the
    # values not in the cycle, a CC flag, then the cycle and another CC flag.
    for kc in [kc for kc in _KNOWN_CYCLES if !(1 in kc)]
        _temp_kc = deepcopy(kc); push!(_temp_kc, _temp_kc[1])
        @test Collatz.hailstone_sequence(kc[1], verbose=false) == _temp_kc
        @test Collatz.hailstone_sequence(kc[1]) == [_CC.CYCLE_INIT, kc, [_CC.CYCLE_LENGTH, length(kc)]]
    end
    # Test the lead into a cycle by entering two of the cycles.
    seq = [kc for kc in _KNOWN_CYCLES if -5 in kc][1]
    seq = append!(append!([seq[2]*4, seq[2]*2], seq[2:length(seq)]), [seq[1]])
    @test Collatz.hailstone_sequence(-56, verbose=false) == push!(deepcopy(seq), seq[3])
    @test Collatz.hailstone_sequence(-56) == [seq[1], seq[2], _CC.CYCLE_INIT, seq[3:length(seq)], [_CC.CYCLE_LENGTH, length(seq)-2]]
    seq = [kc for kc in _KNOWN_CYCLES if -17 in kc][1]
    seq = append!(append!([seq[2]*4, seq[2]*2], seq[2:length(seq)]), [seq[1]])
    @test Collatz.hailstone_sequence(-200, verbose=false) == push!(deepcopy(seq), seq[3])
    @test Collatz.hailstone_sequence(-200) == [seq[1], seq[2], _CC.CYCLE_INIT, seq[3:length(seq)], [_CC.CYCLE_LENGTH, length(seq)-2]]
    # 1's cycle wont yield a description of it being a "cycle" as far as the
    # hailstones are concerned, which is to be expected, so..
    @test Collatz.hailstone_sequence(4, verbose=false) == [4, 2, 1]
    @test Collatz.hailstone_sequence(4) == [4, 2, 1, [_CC.TOTAL_STOPPING_TIME, 2]]
    @test Collatz.hailstone_sequence(16, verbose=false) == [16, 8, 4, 2, 1]
    @test Collatz.hailstone_sequence(16) == [16, 8, 4, 2, 1, [_CC.TOTAL_STOPPING_TIME, 4]]
    # Test the regular stopping time check.
    @test Collatz.hailstone_sequence(4, total_stopping_time=false) == [4, 2, [_CC.STOPPING_TIME, 1]]
    @test Collatz.hailstone_sequence(5, total_stopping_time=false) == [5, 16, 8, 4, [_CC.STOPPING_TIME, 3]]
    # Test small max_total_stopping_time: (minimum internal value is one)
    @test Collatz.hailstone_sequence(4, max_total_stopping_time=-100) == [4, 2, [_CC.MAX_STOP_OOB, 1]]
    # Test the zero stop mid hailing. This wont happen with default params tho.
    @test Collatz.hailstone_sequence(3, P=2, a=3, b=-9) == [3, 0, [_CC.ZERO_STOP, -1]]
    # Lastly, while the function wont let you use a P value of 0, 1 and -1 are
    # still allowed, although they will generate immediate 1 or 2 length cycles
    # respectively, so confirm the behaviour of each of these hailstones.
    @test Collatz.hailstone_sequence(3, P=1, verbose=false) == [3, 3]
    @test Collatz.hailstone_sequence(3, P=-1, verbose=false) == [3, -3, 3]
    # Set P and a to 0 to @test on __assert_sane_parameterisation
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.hailstone_sequence(1, P=0, a=2, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.hailstone_sequence(1, P=0, a=0, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_A) Collatz.hailstone_sequence(1, P=1, a=0, b=3)
end

# Test function stopping_time(initial_value::Integer; P::Integer=2, a::Integer=3, b::Integer=1, max_stopping_time::Integer=1000, total_stopping_time::Bool=false)
@testset verbose = true "stopping_time" begin
    # Test 0's immediated termination.
    @test Collatz.stopping_time(0) == 0
    # The cycle containing 1 wont yield a cycle termination, as 1 is considered
    # the "total stop" that is the special case termination.
    @test Collatz.stopping_time(1) == 0
    # Test the 3 known default parameter's cycles (ignoring [1,4,2])
    # If not verbose, then the result will be the cycle plus the final value
    # being the start of the cycle. If verbose, then the output will be the
    # values not in the cycle, a CC flag, then the cycle and another CC flag.
    for c in [kc for kc in _KNOWN_CYCLES if !(1 in kc)]
        @test Collatz.stopping_time(c[1], total_stopping_time=true) == Base.Inf
    end
    # Test the lead into a cycle by entering two of the cycles. -56;-5, -200;-17
    @test Collatz.stopping_time(-56, total_stopping_time=true) == Base.Inf
    @test Collatz.stopping_time(-200, total_stopping_time=true) == Base.Inf
    # 1's cycle wont yield a description of it being a "cycle" as far as the
    # hailstones are concerned, which is to be expected, so..
    @test Collatz.stopping_time(4, total_stopping_time=true) == 2
    @test Collatz.stopping_time(16, total_stopping_time=true) == 4
    # Test the regular stopping time check.
    @test Collatz.stopping_time(4) == 1
    @test Collatz.stopping_time(5) == 3
    # Test small max_total_stopping_time: (minimum internal value is one)
    @test Collatz.stopping_time(5, max_stopping_time=-100) == nothing
    # Test the zero stop mid hailing. This wont happen with default params tho.
    @test Collatz.stopping_time(3, P=2, a=3, b=-9) == -1
    # Lastly, while the function wont let you use a P value of 0, 1 and -1 are
    # still allowed, although they will generate immediate 1 or 2 length cycles
    # respectively, so confirm the behaviour of each of these stopping times.
    @test Collatz.stopping_time(3, P=1) == Base.Inf
    @test Collatz.stopping_time(3, P=-1) == Base.Inf
    # One last one for the fun of it..
    @test Collatz.stopping_time(27, total_stopping_time=true) == 111
    # And for a bit more fun, common trajectories on
    for x in 0:0  #TODO: Extend this up to 0:4 when hailstone is made arbitrary integer safe.
        @test Collatz.stopping_time(27+576460752303423488*x) == 96
    end
    # Set P and a to 0 to @test on __assert_sane_parameterisation
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.stopping_time(1, P=0, a=2, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.stopping_time(1, P=0, a=0, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_A) Collatz.stopping_time(1, P=1, a=0, b=3)
end

# Test function tree_graph(initial_value::Integer, max_orbit_distance::Integer; P::Integer=2, a::Integer=3, b::Integer=1, __cycle_prevention::Union{Set{Integer},Nothing}=nothing)
@testset verbose = true "tree_graph" begin
    C = _CC.CYCLE_INIT  # Shorthand the cycle terminus
    D = Dict()  # Just to colourise the below in the editor..
    # The default zero trap
    @test Collatz.tree_graph(0, 0) == Dict(0=>D)
    @test Collatz.tree_graph(0, 1) == Dict(0=>Dict(C=>0))
    @test Collatz.tree_graph(0, 2) == Dict(0=>Dict(C=>0))
    # The roundings of the 1 cycle.
    @test Collatz.tree_graph(1, 1) == Dict(1=>Dict(2=>D))
    @test Collatz.tree_graph(1, 0) == Dict(1=>D)
    @test Collatz.tree_graph(1, 1) == Dict(1=>Dict(2=>D))
    @test Collatz.tree_graph(1, 2) == Dict(1=>Dict(2=>Dict(4=>D)))
    @test Collatz.tree_graph(1, 3) == Dict(1=>Dict(2=>Dict(4=>Dict(C=>1,8=>D))))
    @test Collatz.tree_graph(2, 3) == Dict(2=>Dict(4=>Dict(1=>Dict(C=>2),8=>Dict(16=>D))))
    @test Collatz.tree_graph(4, 3) == Dict(4=>Dict(1=>Dict(2=>Dict(C=>4)),8=>Dict(16=>Dict(5=>D,32=>D))))
    # The roundings of the -1 cycle
    @test Collatz.tree_graph(-1, 1) == Dict(-1=>Dict(-2=>D))
    @test Collatz.tree_graph(-1, 2) == Dict(-1=>Dict(-2=>Dict(-4=>D,C=>-1)))
    # Test a wider modulo sweep by upping P to 5, a to 2, and b to 3.
    T = function (x,y) Dict(1=>Dict(-1=>x,5=>y)) end
    orb_1 = T(D,D)
    orb_2 = T(Dict(-5=>D,-2=>D),Dict(C=>1,25=>D))
    T = function (x,y,z) Dict(1=>Dict(-1=>Dict(-5=>x,-2=>y),5=>Dict(C=>1,25=>z))) end
    orb_3 = T(Dict(-25=>D,-4=>D),Dict(-10=>D),Dict(11=>D,125=>D))
    @test Collatz.tree_graph(1, 1, P=5, a=2, b=3) == orb_1
    @test Collatz.tree_graph(1, 2, P=5, a=2, b=3) == orb_2
    @test Collatz.tree_graph(1, 3, P=5, a=2, b=3) == orb_3
    # Test negative P, a and b.
    orb_1 = Dict(1=>Dict(-3=>D))
    T = function (x,y) Dict(1=>Dict(-3=>Dict(-1=>x,9=>y))) end
    orb_2 = T(D,D)
    orb_3 = T(Dict(-2=>D,3=>D),Dict(-27=>D,-7=>D))
    @test Collatz.tree_graph(1, 1, P=-3, a=-2, b=-5) == orb_1
    @test Collatz.tree_graph(1, 2, P=-3, a=-2, b=-5) == orb_2
    @test Collatz.tree_graph(1, 3, P=-3, a=-2, b=-5) == orb_3
    # Set P and a to 0 to @test on __assert_sane_parameterisation
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.tree_graph(1, 1, P=0, a=2, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_P) Collatz.tree_graph(1, 1, P=0, a=0, b=3)
    @test_throws AssertionError(_ErrMsg.SANE_PARAMS_A) Collatz.tree_graph(1, 1, P=1, a=0, b=3)
    # If b is a multiple of a, but not of Pa, then 0 can have a reverse.
    @test Collatz.tree_graph(0, 1, P=17, a=2, b=-6) == Dict(0=>Dict(C=>0,3=>D))
    @test Collatz.tree_graph(0, 1, P=17, a=2, b=102) == Dict(0=>Dict(C=>0))
end

println("################################################################################")

max_stars_repo_name: Skenvy/Collatz

max_issues_repo_name: Skenvy/Collatz

max_forks_repo_name: Skenvy/Collatz


text:
#! format: off

@testset "filter" begin
    n = filter(<(3), n1)
    @test _eval(n) == Block([
        DateTime(2000, 1, 1) => 1,
        DateTime(2000, 1, 2) => 2,
    ])

    n = filter(>=(3), n1)
    @test _eval(n) == Block([
        DateTime(2000, 1, 3) => 3,
        DateTime(2000, 1, 4) => 4,
    ])
end

@testset "skipmissing" begin
    # If no missing inputs, should be a no-op.
    n = empty_node(Int64)
    @test value_type(n) == Int64
    @test n === skipmissing(n)

    # Constant semantics.
    @test constant(1) === skipmissing(constant(1))
    @test skipmissing(constant(missing)).op isa TimeDag.SkipMissing

    # Standard case.
    n = block_node(
        Block([
            DateTime(2000, 1, 1) => missing,
            DateTime(2000, 1, 2) => 2,
            DateTime(2000, 1, 3) => 3,
            DateTime(2000, 1, 4) => missing,
        ]),
    )
    @test value_type(n) == Union{Missing,Int64}
    n2 = skipmissing(n)
    @test value_type(n2) == Int64

    result = _evaluate(n2, DateTime(2000, 1, 1), DateTime(2000, 1, 5))
    @test value_type(result) == Int64
    @test result == Block([
        DateTime(2000, 1, 2) => 2,
        DateTime(2000, 1, 3) => 3,
    ])
end

#! format: on


text:
module SCRIPT

using GLFW, ModernGL, SharedArrays

using ..GraphicsManager
using ..DefaultModelData
using ..CameraManager
using ..FrustumManager
using ..ChunkManager
using ..MeshManager
using ..TextureManager
using ..ShaderManager
using ..ScriptManager

using ..TimeManager
using ..LogManager
using ..Math

const GPU = GraphicsManager

WINDOW = nothing
prevTime = Ref(0.0)
FRAMES = 0
MAX_FRAMES = 0
ITERATION = 0

function showFrames()
  currenttime = time() #GetTimer("FRAME_TIMER")

  global FRAMES, MAX_FRAMES, ITERATION

  ITERATION +=1
  if !OnTime(1.0, prevTime; time=currenttime) FRAMES += 1; return end

  if MAX_FRAMES < FRAMES MAX_FRAMES = FRAMES end
  fps = FRAMES
  max_fps = MAX_FRAMES
  fpms = FRAMES > 0 ? (1000.0 / FRAMES) : 0
  max_fmps = MAX_FRAMES > 0 ? (1000.0 / MAX_FRAMES) : 0
  norm_fps = FRAMES / MAX_FRAMES

  GLFW.SetWindowTitle(WINDOW, "Test - FPS $(round(fps; digits=2))[$(round(max_fps; digits=2))] | FMPS $(round(fpms; digits=2))[$(round(max_fmps; digits=2))] - IT $ITERATION")
  FRAMES = 0
end

function rezizeWindow(window, width, height)
  GLFW.SetWindowSize(window, width, height)
  glViewport(0, 0, width, height)
end

WINDOW_FOCUS = true
function OnFocus(window, focus)
  global WINDOW_FOCUS = focus > 0 ? true : false
end

function OnUpdatedKeys()
  keyValue, keyPressed = getKey()
  if keyPressed
    resetKeys()

    if keyValue == 82 #r
      ScriptManager.reload()
      return false
    end
  end
  true
end

BUFFER = nothing
function init(inputs::Dict{Symbol,Any})
  global WINDOW = inputs[:WINDOW]
  rezizeWindow(WINDOW, 800, 600)
  global BUFFER = zeros(Float32,128^3*6*40) #createBuffer(zeros(Float32,6*128^3),40)
end

function clean()
  println("I got cleaned")
  global BUFFER = nothing
end

CLOSE = true
function update()
  UpdateTimers()
  if CLOSE return end
end

function render()
  showFrames()
  if !WINDOW_FOCUS
    sleep(0.1)
  else
    if !OnUpdatedKeys() return end
  end
end

end

max_stars_repo_name: Gilga/JLGLAssets

max_issues_repo_name: Gilga/JLGLAssets

max_forks_repo_name: Gilga/JLGLAssets


text:
using Dates

lines = readlines("04.input")

function parseevent(event::AbstractString)
    num = match(r"(\d+)", event)
    if num isa Nothing
        return event == "falls asleep" ? :sleep : :wake
    else
        return parse(Int, num.captures[1])
    end
end

function parseline(line)
    date_string, event = match(r"\[(.*)\] (.*)", line).captures
    datetime = Dates.DateTime(date_string, dateformat"y-m-d H:M")
    event = parseevent(event)
    datetime, event
end

input = sort!(parseline.(lines), by=first)

function collect_sleeps(input)
    id = -1
    last_sleep = DateTime(0)
    sleeps = Dict{Int,Vector{Tuple{DateTime,DateTime}}}()
    for (time, event) in input
        if event isa Int
            id = event
        elseif event == :sleep
            last_sleep = time
        else # event == :wake
            push!(get!(sleeps, id, Vector{Tuple{DateTime,DateTime}}()),
                  (last_sleep, time))
        end
    end
    sleeps
end

sleeps = collect_sleeps(input)

function find_sleepiest_min(sleeps::Vector{Tuple{DateTime,DateTime}})
    mins_asleep = zeros(Int, 60)
    for (sleep, wake) in sleeps
        sleep_min = minute(sleep)+1 # one-based indexing ;)
        wake_min = minute(wake)
        mins_asleep[sleep_min:wake_min] .+= 1
    end
    sleep_freq, sleepiest_min = findmax(mins_asleep) .- (0,1)
end

total_sleeps = Dict(id => mapreduce(((x,y),) -> Minute(y-x), +, sleeps)
                    for (id, sleeps) in sleeps)
mins, id = findmax(total_sleeps)
_, sleepiest_min = find_sleepiest_min(sleeps[id])

@show id * sleepiest_min

(freq, mins2), id2 = findmax(Dict(id => find_sleepiest_min(sleeps)
                                  for (id, sleeps) in sleeps))
@show id2 * mins2


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ImageOSDisk <: SwaggerModel
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    osState::Any # spec type: Union{ Nothing, String } # spec name: osState
    snapshot::Any # spec type: Union{ Nothing, SubResource } # spec name: snapshot
    managedDisk::Any # spec type: Union{ Nothing, SubResource } # spec name: managedDisk
    blobUri::Any # spec type: Union{ Nothing, String } # spec name: blobUri
    caching::Any # spec type: Union{ Nothing, String } # spec name: caching
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    storageAccountType::Any # spec type: Union{ Nothing, StorageAccountType } # spec name: storageAccountType

    function ImageOSDisk(;osType=nothing, osState=nothing, snapshot=nothing, managedDisk=nothing, blobUri=nothing, caching=nothing, diskSizeGB=nothing, storageAccountType=nothing)
        o = new()
        validate_property(ImageOSDisk, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(ImageOSDisk, Symbol("osState"), osState)
        setfield!(o, Symbol("osState"), osState)
        validate_property(ImageOSDisk, Symbol("snapshot"), snapshot)
        setfield!(o, Symbol("snapshot"), snapshot)
        validate_property(ImageOSDisk, Symbol("managedDisk"), managedDisk)
        setfield!(o, Symbol("managedDisk"), managedDisk)
        validate_property(ImageOSDisk, Symbol("blobUri"), blobUri)
        setfield!(o, Symbol("blobUri"), blobUri)
        validate_property(ImageOSDisk, Symbol("caching"), caching)
        setfield!(o, Symbol("caching"), caching)
        validate_property(ImageOSDisk, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(ImageOSDisk, Symbol("storageAccountType"), storageAccountType)
        setfield!(o, Symbol("storageAccountType"), storageAccountType)
        o
    end
end # type ImageOSDisk

const _property_map_ImageOSDisk = Dict{Symbol,Symbol}(Symbol("osType")=>Symbol("osType"), Symbol("osState")=>Symbol("osState"), Symbol("snapshot")=>Symbol("snapshot"), Symbol("managedDisk")=>Symbol("managedDisk"), Symbol("blobUri")=>Symbol("blobUri"), Symbol("caching")=>Symbol("caching"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("storageAccountType")=>Symbol("storageAccountType"))
const _property_types_ImageOSDisk = Dict{Symbol,String}(Symbol("osType")=>"String", Symbol("osState")=>"String", Symbol("snapshot")=>"SubResource", Symbol("managedDisk")=>"SubResource", Symbol("blobUri")=>"String", Symbol("caching")=>"String", Symbol("diskSizeGB")=>"Int32", Symbol("storageAccountType")=>"StorageAccountType")
Base.propertynames(::Type{ ImageOSDisk }) = collect(keys(_property_map_ImageOSDisk))
Swagger.property_type(::Type{ ImageOSDisk }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ImageOSDisk[name]))}
Swagger.field_name(::Type{ ImageOSDisk }, property_name::Symbol) =  _property_map_ImageOSDisk[property_name]

const _allowed_ImageOSDisk_osType = ["Windows", "Linux"]

const _allowed_ImageOSDisk_osState = ["Generalized", "Specialized"]

const _allowed_ImageOSDisk_caching = ["None", "ReadOnly", "ReadWrite"]

function check_required(o::ImageOSDisk)
    true
end

function validate_property(::Type{ ImageOSDisk }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "ImageOSDisk", :enum, val, _allowed_ImageOSDisk_osType)
    end
    if name === Symbol("osState")
        Swagger.validate_param(name, "ImageOSDisk", :enum, val, _allowed_ImageOSDisk_osState)
    end
    if name === Symbol("caching")
        Swagger.validate_param(name, "ImageOSDisk", :enum, val, _allowed_ImageOSDisk_caching)
    end
end


text:
################################################################################
#
#  Change of ring
#
################################################################################

@doc Markdown.doc"""
    restrict_scalars(A::AlgAss{nf_elem}, Q::FlintRationalField)
    restrict_scalars(A::AlgAss{fq_nmod}, Fp::GaloisField)
    restrict_scalars(A::AlgAss{fq}, Fp::Generic.ResField{fmpz})
      -> AlgAss, Map

Given an algebra $A$ over a field $L$ and the prime field $K$ of $L$, this
function returns the restriction $B$ of $A$ to $K$ and the morphism $B \to A$.
"""
# Top level functions to avoid "type mix-ups" (like AlgAss{fq_nmod} with FlintQQ)
function restrict_scalars(A::AbsAlgAss{nf_elem}, Q::FlintRationalField)
  return _restrict_scalars(A, Q)
end

function restrict_scalars(A::AbsAlgAss{fq_nmod}, Fp::GaloisField)
  return _restrict_scalars(A, Fp)
end

function restrict_scalars(A::AbsAlgAss{fq}, Fp::GaloisFmpzField)
  return _restrict_scalars(A, Fp)
end

#function restrict_scalars(A::AbsAlgAss{gfp_elem}, Fp::GaloisField)
#  function AtoA(x::AlgAssElem)
#    return x
#  end
#  return A, AtoA, AtoA
#end

@doc Markdown.doc"""
    restrict_scalars(A::AlgAss{nf_elem}, KtoL::NfToNfMor)
      -> AlgAss, Map

Given an algebra $A$ over a number field $L$ and an inclusion map `KtoL` from
a number field $K$ to $L$, this function returns the restriction $B$ of $A$
to $K$ and the morphism $B \to A$.
"""
function restrict_scalars(A::AbsAlgAss{nf_elem}, KtoL::NfToNfMor)
  return _restrict_scalars(A, KtoL)
end

#function restrict_scalars(A::AbsAlgAss{Nemo.gfp_fmpz_elem}, Fp::Nemo.GaloisFmpzField)
#  function AtoA(x::AlgAssElem)
#    return x
#  end
#  return A, AtoA, AtoA
#end

mutable struct AlgAssResMor{S, T, U, V} <: Map{S, T, HeckeMap, AlgAssResMor}
  header::MapHeader{S, T}
  domain::S
  codomain::T
  f::U
  B::V

  function AlgAssResMor(B::S, A::T, f::U, Bas::V) where {S, T, U, V}
    z = new{S, T, U, V}()
    z.domain = B
    z.codomain = A
    z.B = Bas
    z.f = f
    z.header = MapHeader(B, A)
    return z
  end
end

degree(::FlintRationalField) = 1

domain(f::AlgAssResMor) = f.domain

codomain(f::AlgAssResMor) = f.codomain

function preimage(f::AlgAssResMor, a)
  @assert parent(a) == codomain(f)
  B = domain(f)
  A = codomain(f)
  m = div(dim(B), dim(A))
  y = coefficients(a, copy = false)
  yy = Vector{elem_type(base_ring(B))}(undef, dim(B))
  for i in 1:length(y)
    ee = image(f.f, y[i])
    for j in 1:length(ee)
      yy[(i - 1)*m + j] = ee[j]
    end
  end
  return domain(f)(yy)
end

function image(f::AlgAssResMor, a)
  @assert parent(a) == domain(f)
  A = codomain(f)
  B = domain(f)
  d = div(dim(B), dim(A))
  y = coefficients(a, copy = false)
  yy = Vector{elem_type(base_ring(A))}(undef, dim(A))
  for i in 1:dim(A)
    yy[i] = preimage(f.f, y[(i - 1)*d + 1:i*d])
  end
  return A(yy)
end

#function _restrict_scalars_to_prime_field(A::AlgAss{T}, prime_field::Union{FlintRationalField, GaloisField, Generic.ResField{fmpz}}) where { T <: Union{nf_elem, fq_nmod, fq} }
# TODO: fix the type
function _restrict_scalars(A::AbsAlgAss{T}, prime_field) where { T }
  K = base_ring(A)
  n = dim(A)
  # We use b * A[i], b running through a basis of A over prime_field
  # the basis for A over the prime field.
  # Precompute the powers a^k:

  f = FldToVecMor(K, prime_field)
  Bas = f.B
  F = f.K

  m = div(degree(K), degree(F))
  nm = n*m

  absbasis = elem_type(A)[]
  BA = basis(A)
  for i in 1:n
    for j in 1:length(Bas)
      push!(absbasis, Bas[j] * BA[i])
    end
  end

  m1 = m - 1
  mult_table = zeros(F, nm, nm, nm)

  Aij = A()
  t = A()
  tt = K()
  firstindex = 1
  for i in 1:n
    for i2 in 1:m
      secondindex = 1
      for j in 1:n
        Aij = mul!(Aij, A[i], A[j])
        if iszero(Aij)
          secondindex += m
          continue
        end

        for j2 in 1:m
          e = Bas[i2] * Bas[j2] * Aij
          y = Vector{elem_type(F)}(undef, nm)
          yy = coefficients(e, copy = false)
          for i in 1:n
            ee = image(f, yy[i])
            @assert length(ee) == m
            for j in 1:m
              y[(i - 1)*m + j] = ee[j]
            end
          end
          mult_table[firstindex, secondindex, :] = y
          secondindex += 1
        end
      end
      firstindex += 1
    end
  end

  e = one(A)
  y = Vector{elem_type(F)}(undef, nm)
  yy = coefficients(e, copy = false)
  for i in 1:n
    ee = image(f, yy[i])
    for j in 1:m
      y[(i - 1)*m + j] = ee[j]
    end
  end

  B = AlgAss(F, mult_table, y)
  B.iscommutative = iscommutative(A)

  return B, AlgAssResMor(B, A, f, absbasis)
end

mutable struct AlgAssExtMor{S, T, U, V, W, X, Y} <: Map{S, T, HeckeMap, AlgAssExtMor}
  header::MapHeader{S, T}
  domain::S
  codomain::T
  f::U
  B::V
  BB::W
  MMinv::X
  BAoverC::Y

  function AlgAssExtMor(B::S, A::T, f::U, BC::V, BL::W, MMinv::X, BAoverC::Y) where {S, T, U, V, W, X, Y}
    z = new{S, T, U, V, W, X, Y}()
    z.domain = B
    z.codomain = A
    z.f = f
    z.B = BC
    z.BB = BL
    z.BAoverC = BAoverC
    z.MMinv = MMinv
    z.header = MapHeader(B, A)
    return z
  end
end

domain(f::AlgAssExtMor) = f.domain

codomain(f::AlgAssExtMor) = f.codomain

function image(f::AlgAssExtMor, a)
  @assert parent(a) == domain(f)
  B = domain(f)
  A = codomain(f)
  m = div(dim(A), dim(B))
  y = coefficients(a, copy = false)
  z = zero(A)
  for i in 1:length(y)
    z = z + dot(f.B, (f.f\y[i]).coeffs) * f.BAoverC[i]
  end
  return z
end

function preimage(f::AlgAssExtMor, a)
  @assert parent(a) == codomain(f)
  A = codomain(f)
  B = domain(f)
  d = div(dim(A), dim(B))
  _y = matrix(base_ring(A), 1, dim(A), coefficients(a, copy = false)) * f.MMinv
  y = elem_type(base_ring(A))[_y[1, i] for i in 1:dim(A)]
  yy = Vector{elem_type(base_ring(B))}(undef, dim(B))
  for i in 1:dim(B)
    yy[i] = f.f(dot(basis(domain(f.f)), y[(i - 1)*d + 1:i*d]))
  end
  return B(yy)
end

function _as_algebra_over_center(A::AlgAss{T}) where { T <: Union{nf_elem, fmpq}}
  @assert !iszero(A)

  K = base_ring(A)
  C, CtoA = center(A)
  fields = as_number_fields(C)
  @assert length(fields) == 1
  L, CtoL = fields[1]
  # Maybe do something more efficient A iscentral
  return __as_algebra_over_center(A, K, L, CtoA, CtoL)
end

function _as_algebra_over_center(A::AlgAss{T}) where { T } #<: Union{fmpq, gfp_elem, Generic.ResF{fmpz}, fq, fq_nmod} }
  @assert !iszero(A)

  K = base_ring(A)
  C, CtoA = center(A)

  # Maybe do something more efficient A iscentral

  L, CtoL = _as_field_with_isomorphism(C)
  return __as_algebra_over_center(A, K, L, CtoA, CtoL)
end

function __as_algebra_over_center(A, K, L, CtoA, CtoL)
  C = domain(CtoA)

  iscentral = ( dim(C) == 1 )

  basisC = basis(C)
  basisCinA = Vector{elem_type(A)}(undef, dim(C))
  basisCinL = Vector{elem_type(L)}(undef, dim(C))
  for i = 1:dim(C)
    basisCinA[i] = CtoA(basisC[i])
    basisCinL[i] = CtoL(basisC[i])
  end

  # We construct a basis of A over C (respectively L) by using the following fact:
  # A subset M of basis(A) is a C-basis of A if and only if |M| = dim(A)/dim(C)
  # and all possible products of elements of M and basisCinA form a K-basis of A,
  # with K := base_ring(A).
  AoverK = basis(A)
  AoverC = Vector{Int}()
  M = zero_matrix(K, dim(C), dim(A))
  MM = zero_matrix(K, 0, dim(A))
  r = 0
  for i = 1:dim(A)
    b = AoverK[i]

    for j = 1:dim(C)
      elem_to_mat_row!(M, j, b*basisCinA[j])
    end

    N = vcat(MM, M)
    s = rank(N)
    if s > r
      push!(AoverC, i)
      MM = N
      r = s
    end
    if r == dim(A)
      break
    end
  end

  m = div(dim(A), dim(C))

  @assert length(AoverC) == m
  @assert nrows(MM) == dim(A)

  iMM = inv(MM)

  dC = dim(C)

  @assert dC * m == dim(A)

  mult_table = zeros(L, m, m, m)
  Aij = A()
  for i = 1:m
    for j = 1:m
      Aij = mul!(Aij, A[AoverC[i]], A[AoverC[j]])
      if iszero(Aij)
        continue
      end

      xx = matrix(K, 1, dim(A), coefficients(Aij, copy = false))
      xxM = xx * iMM

      y = Vector{elem_type(L)}(undef, m)

      for k in 1:m
        # This linear indexing is fine, since it is a row/vector
        y[k] = CtoL(dot(basis(domain(CtoL)), elem_type(K)[xxM[(k - 1) * dC + l] for l in 1:dC]))
      end

      mult_table[i, j, :] = y
    end
  end

  xx = matrix(K, 1, dim(A), coefficients(one(A), copy = false))
  xxM = xx * iMM

  y = Vector{elem_type(L)}(undef, m)

  for k in 1:m
    # This linear indexing is fine, since it is a row/vector
    y[k] = CtoL(dot(basis(domain(CtoL)), elem_type(K)[xxM[(k - 1) * dC + l] for l in 1:dC]))
  end

  B = AlgAss(L, mult_table, y)
  B.iscommutative = A.iscommutative

  BtoA = AlgAssExtMor(B, A, CtoL, basisCinA, basisCinL, iMM, elem_type(A)[A[i] for i in AoverC])

  return B, BtoA
end

mutable struct PrimeFieldEmbedStub{S, T}
  domain::S
  codomain::T
end

function FldToVecMor(f::PrimeFieldEmbedStub)
  L = codomain(f)
  K = domain(f)
  return FldToVecMor(L, K)
end

function _embed_center_into_field(m::AbsAlgAssToNfAbsMor{AlgAss{fmpq}})
  return PrimeFieldEmbedStub(base_ring(domain(m)), codomain(m))
end

function _embed_center_into_field(m::AbsAlgAssToNfAbsMor{AlgAss{nf_elem}})
  return PrimeFieldEmbedStub(base_ring(domain(m)), codomain(m))
end

function _embed_center_into_field(m::AbsAlgAssToFqMor{AlgAss{gfp_elem}})
  return PrimeFieldEmbedStub(base_ring(domain(m)), codomain(m))
end

function _embed_center_into_field(m::AbsAlgAssToFqMor{AlgAss{gfp_fmpz_elem}})
  return PrimeFieldEmbedStub(base_ring(domain(m)), codomain(m))
end

domain(f::PrimeFieldEmbedStub) = f.domain

codomain(f::PrimeFieldEmbedStub) = f.codomain


text:
using Images, MXNet

### LOADING THE MODEL
const MODEL_NAME = "weights/mobilenet-v2/mobilenet_v2"
const MODEL_CLASS_NAMES = "weights/mobilenet-v2/synset.txt"

nnet = mx.load_checkpoint(MODEL_NAME, 0, mx.FeedForward; context = mx.cpu());
synset = readlines(MODEL_CLASS_NAMES);

### SEARCH FOR A LAYER OF INTERESTS
layers = mx.get_internals(nnet.arch);
layers_flatten = nothing
layers_to_remove = Symbol[]

# We iterate over all layers until we find the one matching our requirements
# and remove the ones to follow after
for i = 1:2000

    layer = layers[i];
    layer_name = mx.get_name(layer)

    if layers_flatten == nothing && layer_name == :pool6
        layers_flatten = layer
    elseif layers_flatten != nothing
        push!(layers_to_remove, layer_name)
        if layer_name in [:softmax, :label, :prob] break end
    end
end

# UPDATE NETWORK STRUCTURE WITH NEW
nnet.arch = @mx.chain layers_flatten => Flatten()

### REMOVE SOFTMAX OUTPUT AND FULLY CONNECTED LAYER
map(x -> delete!(nnet.arg_params, x), layers_to_remove);
map(x -> delete!(nnet.aux_params, x), layers_to_remove);

mx.save_checkpoint(nnet, "weights/mobilenet-v2/MobiletNet-FE", mx.OptimizationState(1, 0, 0, 0))

### PREPARING THE INPUT

img = imresize(load("sample-images/bird-3183441_640.jpg"), (224, 224));
img = permutedims(Float16.(channelview(img)), (3, 2, 1));

img[:, :, :] *= 256.0;
img[:, :, :] -= 128.;
img[:, :, :] /= 128.;

img = reshape(img, 224, 224, 3, 1);

mx_data = mx.zeros((224, 224, 3, 1));
mx_data[1:1] = img;
data_provider = mx.ArrayDataProvider(:data => mx_data);

### PREDICTING

@time pred = mx.predict(nnet, data_provider)


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

using LaTeXStrings
using Plots
pyplot()

# pick a pixel (xprime, yprime), and trace a full ray through the disk

function plot_ray(xprime, yprime)

    vlos = 1.0 # km/s

    # Get a decent estimate of DeltaVmax
    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # DeltaVmax = 0.5 # km/s
    rmax = 600 * AU
    nzp = 512
    zps = linspace(-600 * AU, 600 * AU, nzp)

    rcyls = Array{Float64}(nzp)
    zs = Array{Float64}(nzp)
    Ss = Array{Float64}(nzp)
    Us = Array{Float64}(nzp)
    Deltavs = Array{Float64}(nzp)
    DeltaVs = Array{Float64}(nzp)
    Deltas = Array{Float64}(nzp)
    alphas = Array{Float64}(nzp)

    for i=1:nzp
        zprime = zps[i]
        rcyl = DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)
        z = DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

        rcyls[i] = rcyl
        zs[i] = z

        Ss[i] = DiskTracer.model.S_nu(rcyl, z, mol.nu_0, pars)
        Us[i] = DiskTracer.model.Upsilon_nu(rcyl, z, mol.nu_0, pars, mol)

        # Calculate Delta v at this position
        Deltav = vlos * 1e5 - DiskTracer.geometry.get_vlos(xprime, yprime, zprime, pars)
        Deltavs[i] = Deltav

        DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)
        DeltaVs[i] = DeltaV * 1e-5 # km/s

        Deltas[i] = exp(-Deltav^2/DeltaV^2)

        alphas[i] = Us[i] * exp(-Deltav^2/DeltaV^2)

    end

    # get bounding zps
    zbounds = DiskTracer.geometry.get_bounding_zps(xprime, yprime, pars, vlos, DeltaVmax, rmax)
    if length(zbounds) == 2
        zstart, zend = zbounds
        println("H bound verified merged: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, zstart, zend, pars, 0.5))
    else
        z1start, z1end, z2start, z2end = zbounds
        println("H bound verified front-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z1start, z1end, pars, 0.5))
        println("H bound verified rear-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z2start, z2end, pars, 0.5))
    end

    # hbounds, derivs = DiskTracer.geometry.get_bounding_scale_heights(xprime, yprime, pars, 1.0)

    p1 = plot(zps./AU, rcyls./AU, ylabel=L"$r_\mathrm{cyl}$ [AU]")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p2 = plot(zps./AU, zs./AU, ylabel=L"$z$ [AU]")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p25 = plot(zps./AU, abs.(zs./rcyls), ylabel=L"$|z|/r_\mathrm{cyl}$", ylims=(0, 3))
    # vline!((hbounds./AU)', line=([:black :black]))

    p26 = plot(zps./AU, Deltavs .* 1e-5, ylabel=L"$\Delta v$")

    p3 = plot(zps./AU, Ss, ylabel=L"$S_\nu$", yaxis=(scale=:log10))
    # p3 = plot(zps./AU, Ss, ylabel=L"$S_\nu$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p4 = plot(zps./AU, Us, xlabel=L"$z^\prime$", ylabel=L"$\Upsilon$")# , yaxis=(scale=:log10))
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p45 = plot(zps./AU, DeltaVs, xlabel=L"$z^\prime$", ylabel=L"$\Delta V$")

    p5 = plot(zps./AU, Deltas, xlabel=L"$z^\prime$", ylabel=L"$\exp(-\Delta v^2/\Delta V^2)$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p6 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p = plot(p1,p2,p25,p26,p3,p4,p45,p5,p6,layout=(9,1), legend=false, size=(800,1200))

    savefig(p, "ray_plot.png")

end

plot_ray(100.*AU, 200.0*AU)


text:
#####
##### Time-stepping functionality that is independent of the TimeStepper
#####

"""
    calculate_explicit_substep!(tendencies, velocities, tracers, pressures, diffusivities, model)

Calculate the initial and explicit substep of the two-step fractional step method with pressure correction.
"""
function calculate_explicit_substep!(tendencies, velocities, tracers, pressures, diffusivities, model)
    time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)
    calculate_tendencies!(tendencies, velocities, tracers, pressures, diffusivities, model)
    return nothing
end

"""
    time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)

Perform precomputations necessary for an explicit timestep or substep.
"""
function time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)

    fill_halo_regions!(merge(model.velocities, model.tracers), model.architecture,
                       model.clock, state(model))

    calculate_diffusivities!(diffusivities, model.architecture, model.grid, model.closure,
                             model.buoyancy, velocities, tracers)

    fill_halo_regions!(model.diffusivities, model.architecture, model.clock, state(model))

    @launch(device(model.architecture), config=launch_config(model.grid, :xy),
            update_hydrostatic_pressure!(pressures.pHY′, model.grid, model.buoyancy, tracers))

    fill_halo_regions!(model.pressures.pHY′, model.architecture)

    return nothing
end

"""
    calculate_tendencies!(diffusivities, pressures, velocities, tracers, model)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(tendencies, velocities, tracers, pressures, diffusivities, model)

    # Note:
    #
    # "tendencies" is a NamedTuple of OffsetArrays corresponding to the tendency data for use
    # in GPU computations.
    #
    # "model.timestepper.Gⁿ" is a NamedTuple of Fields, whose data also corresponds to
    # tendency data.

    # Arguments needed to calculate tendencies for momentum and tracers

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_interior_tendency_contributions!(tendencies, model.architecture, model.grid,
                                               model.coriolis, model.buoyancy, model.surface_waves, model.closure,
                                               velocities, tracers, pressures.pHY′, diffusivities, model.forcing,
                                               model.clock)

    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_boundary_tendency_contributions!(model.timestepper.Gⁿ, model.architecture, model.velocities,
                                               model.tracers, model.clock, state(model))

    # Calculate momentum tendencies on boundaries in `Bounded` directions.
    #calculate_velocity_tendencies_on_boundaries!(tendencies, model.architecture, model.grid,
    #                                             model.coriolis, model.buoyancy, model.surface_waves, model.closure,
    #                                             velocities, tracers, pressures.pHY′, diffusivities, model.forcing,
    #                                             model.clock)

    return nothing
end

"""
    calculate_pressure_correction!(nonhydrostatic_pressure, Δt, tendencies, velocities, model)

Calculate the (nonhydrostatic) pressure correction associated `tendencies`, `velocities`, and step size `Δt`.
"""
function calculate_pressure_correction!(nonhydrostatic_pressure, Δt, predictor_velocities, model)
    fill_halo_regions!(model.timestepper.predictor_velocities, model.architecture,
                       model.clock, state(model))

    solve_for_pressure!(nonhydrostatic_pressure, model.pressure_solver,
                        model.architecture, model.grid, Δt, predictor_velocities)

    fill_halo_regions!(model.pressures.pNHS, model.architecture)

    return nothing
end

#####
##### Fractional and time stepping
#####

"""
Update the horizontal velocities u and v via

    `u^{n+1} = u^n + (Gu^{n+½} - δₓp_{NH} / Δx) Δt`

Note that the vertical velocity is not explicitly time stepped.
"""
function _fractional_step_velocities!(U, grid, Δt, pNHS)
    @loop_xyz i j k grid begin
        @inbounds U.u[i, j, k] -= ∂xᶠᵃᵃ(i, j, k, grid, pNHS) * Δt
        @inbounds U.v[i, j, k] -= ∂yᵃᶠᵃ(i, j, k, grid, pNHS) * Δt
    end
    return nothing
end

"Update the solution variables (velocities and tracers)."
function fractional_step_velocities!(U, C, arch, grid, Δt, pNHS)
    @launch device(arch) config=launch_config(grid, :xyz) _fractional_step_velocities!(U, grid, Δt, pNHS)
    return nothing
end


text:
# Linking of different PTX modules

import Base: unsafe_convert, show, cconvert

export
    CuLink, complete, destroy,
    addData, addFile

typealias CuLinkState_t Ptr{Void}

immutable CuLink
    handle::CuLinkState_t

    options::Dict{CUjit_option,Any}
    optionKeys::Vector{CUjit_option}
    optionVals::Vector{Ptr{Void}}

    function CuLink()
        handle_ref = Ref{CuLinkState_t}()

        options = Dict{CUjit_option,Any}()
        options[ERROR_LOG_BUFFER] = Array(UInt8, 1024*1024)
        @static if DEBUG
            options[GENERATE_LINE_INFO] = true
            options[GENERATE_DEBUG_INFO] = true

            options[INFO_LOG_BUFFER] = Array(UInt8, 1024*1024)
            options[LOG_VERBOSE] = true
        end
        optionKeys, optionVals = encode(options)

        @apicall(:cuLinkCreate, (Cuint, Ptr{CUjit_option}, Ptr{Ptr{Void}}, Ptr{CuLinkState_t}),
                                length(optionKeys), optionKeys, optionVals, handle_ref)

        new(handle_ref[], options, optionKeys, optionVals)
    end
end

unsafe_convert(::Type{CuLinkState_t}, link::CuLink) = link.handle
show(io::IO,link::CuLink) = print(io, typeof(link), "(", link.handle, ")")

"""
Complete a pending linker invocation, returning an output image.

Note that this image is still owned by the JIT, meaning it will be destroyed upon calling `destroy` on the linker.
"""
function complete(link::CuLink)
    cubin_ref = Ref{Ptr{Void}}()
    size_ref = Ref{Csize_t}()

    try
        @apicall(:cuLinkComplete, (CuLinkState_t, Ptr{Ptr{Void}}, Ptr{Csize_t}),
                                  link, cubin_ref, size_ref)
    catch err
        (err == ERROR_NO_BINARY_FOR_GPU || err == ERROR_INVALID_IMAGE) || rethrow(err)
        options = decode(link.optionKeys, link.optionVals)
        rethrow(CuError(err.code, options[ERROR_LOG_BUFFER]))
    end

    @static if DEBUG
        options = decode(link.optionKeys, link.optionVals)
        if isempty(options[INFO_LOG_BUFFER])
            debug("JIT info log is empty")
        else
            debug("JIT info log: ", repr_indented(options[INFO_LOG_BUFFER]))
        end
    end

    return unsafe_wrap(Array, convert(Ptr{UInt8}, cubin_ref[]), size_ref[])
end

"""
Destroys state for a JIT linker invocation.

Note that this destroys any existing references to images generated by this JIT, despite having live Julia objects pointing to them.
"""
function destroy(link::CuLink)
    @apicall(:cuLinkDestroy, (CuLinkState_t,), link)
end

function addData(link::CuLink, name::String, data::Union{Vector{UInt8},String}, typ::CUjit_input)
    # NOTE: ccall can't directly convert String to Ptr{Void}, so step through a typed Ptr
    if typ == PTX
        # additionally, in the case of PTX there shouldn't be any embedded NULLs
        typed_ptr = pointer(unsafe_convert(Cstring, cconvert(Cstring, String(data))))
    else
        typed_ptr = pointer(unsafe_convert(Vector{UInt8}, cconvert(Vector{UInt8}, data)))
    end
    untyped_ptr = convert(Ptr{Void}, typed_ptr)

    @apicall(:cuLinkAddData,
             (CuLinkState_t, CUjit_input, Ptr{Void}, Csize_t, Cstring, Cuint, Ptr{CUjit_option}, Ptr{Ptr{Void}}),
             link, typ, untyped_ptr, length(data), name, 0, C_NULL, C_NULL)

    return nothing
end

function addFile(link::CuLink, path::String, typ::CUjit_input)
    @apicall(:cuLinkAddFile,
             (CuLinkState_t, CUjit_input, Cstring, Cuint, Ptr{CUjit_option}, Ptr{Ptr{Void}}),
             link, typ, path, 0, C_NULL, C_NULL)

    return nothing
end


text:
function compute_estimated_iif!(ata_model)
    T = ata_model.settings.T
    n_items = ata_model.settings.n_items
    irt_parameters = ata_model.settings.irt.parameters
    irt_model = ata_model.settings.irt.model
    irt_D = ata_model.settings.irt.D
    irt_parametrization = ata_model.settings.irt.parametrization
    IIF = Vector{Array{Float64,2}}(undef, T)
    K = zeros(Int, T)
    for t = 1:T
        K[t] = size(ata_model.obj.cores[t].points, 1)
        IIF[t] = zeros(K[t], n_items)
        for k = 1:K[t]
            IIF[t][k, :] = item_info(
                irt_parameters,
                ata_model.obj.cores[t].points[k],
                model = irt_model,
                parametrization = irt_parametrization,
                D = irt_D,
            )# K[t] x I
        end
        ata_model.obj.cores[t].IIF = IIF[t]
    end
    JLD2.@save "opt/IIF.jld2" IIF
    return nothing
end


text:
d = rand(2:4)
A, B = rand(d, d), rand(d, d)
α = rand()

@test isapprox(det(A * B), det(A) * det(B))
@test isapprox(det(A * B), det(B) * det(A))
@test isapprox(det(A * B), det(B * A))
@test isapprox(det(convert(Array{Float64}, transpose(A))), det(A))
@test isapprox(det(α .* A), α^d * det(A))

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta