text:
using Profile
using LinearMaps
using BenchmarkTools
#using KrylovKit
using IterativeSolvers

"""
Represents a state in an set of abitrary (yet dense) subspaces of a set of FockConfigs.

v[FockConfig][TuckerConfig] => Dense Tensor
E.g., used in n-body Tucker

    clusters::Vector{Cluster}
    data::OrderedDict{FockConfig,OrderedDict{TuckerConfig,Array}}
    p_spaces::Vector{ClusterSubspace}
    q_spaces::Vector{ClusterSubspace}
"""
struct TuckerState{T,N} <: AbstractState
    clusters::Vector{Cluster}
    data::OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Array{T}}}
    p_spaces::Vector{ClusterSubspace}
    q_spaces::Vector{ClusterSubspace}
end
Base.haskey(ts::TuckerState, i) = return haskey(ts.data,i)
Base.getindex(ts::TuckerState, i) = return ts.data[i]
Base.setindex!(ts::TuckerState, i, j) = return ts.data[j] = i
Base.iterate(ts::TuckerState, state=1) = iterate(ts.data, state)

# Conversions
#Base.convert(::Type{TuckerConfig}, input::Vector{UnitRange{T}}) where T<:Integer = TuckerConfig(input)

"""
    TuckerState(clusters, p_spaces, q_spaces, foi; nroots=1)

Constructor
- `clusters::Vector{Cluster}`
- `p_spaces::Vector{ClusterSubspace}`
- `q_spaces::Vector{ClusterSubspace}`
- `na::Int` Number of alpha
- `nb::Int` Number of beta
"""
function TuckerState(clusters::Vector{Cluster}, p_spaces::Vector{ClusterSubspace}, q_spaces::Vector{ClusterSubspace}, na, nb; nroots=1)
   #={{{=#
    length(p_spaces) == length(clusters) || error("# of clusters don't match # of subspaces")
    length(p_spaces) == length(q_spaces) || error(" p_spaces/q_spaces don't have same length", length(p_space), length(q_space))
    for ci in clusters
        for fock in p_spaces[ci.idx].data
            if haskey(q_spaces[ci.idx].data, fock)
                all(collect(p_spaces[ci.idx][fock]) .!= collect(q_spaces[ci.idx][fock])) || error(" Not orthogonal")
            end
        end
    end

    N = length(clusters)
    #s = TuckerState(clusters)
    data = OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Array{Float64}} }()
    ns = []
    for cssi in p_spaces
        nsi = []
        for (fock,range) in cssi.data
            push!(nsi,fock)
        end
        push!(ns,nsi)
    end

    for newfock in product(ns...)
        nacurr = 0
        nbcurr = 0
        for c in newfock
            nacurr += c[1]
            nbcurr += c[2]
        end
        if (nacurr == na) && (nbcurr == nb)

            fockconfig = FockConfig(collect(newfock))

            tuckconfig = TuckerConfig([p_spaces[ci.idx].data[newfock[ci.idx]] for ci in clusters])
#            tuckconfig = TuckerConfig()
#            for ci in clusters
#                cssi = p_spaces[ci.idx]
#                push!(tuckconfig, cssi.data[newfock[ci.idx]])
#            end

            haskey(data, fockconfig) == false || error(" here:", fockconfig)
            data[fockconfig] = OrderedDict(tuckconfig => zeros((size(tuckconfig)...,nroots)) )
            #add_fockconfig!(data,fockconfig)
            #data[fockconfig][tuckconfig] = zeros((size(tuckconfig)...,nroots))  # todo - finish this
        end
    end
    return TuckerState(clusters, data, p_spaces, q_spaces)
#=}}}=#
end

"""
    TuckerState(clusters, p_spaces, q_spaces, foi; nroots=1)

Constructor to build state directly from a definition of a first order interacting space (or more generic even i suppose)
- `clusters::Vector{Cluster}`
- `p_spaces::Vector{ClusterSubspace}`
- `q_spaces::Vector{ClusterSubspace}`
- `foi::OrderedDict{FockConfig,Vector{TuckerConfig}}`
- `nroots`
"""
function TuckerState(clusters::Vector{Cluster}, p_spaces::Vector{ClusterSubspace}, q_spaces::Vector{ClusterSubspace},
        foi::OrderedDict{FockConfig,Vector{TuckerConfig}}; nroots=1)
   #={{{=#
    length(p_spaces) == length(clusters) || error("# of clusters don't match # of subspaces")
    length(p_spaces) == length(q_spaces) || error(" p_spaces/q_spaces don't have same length", length(p_space), length(q_space))
    for ci in clusters
        for fock in p_spaces[ci.idx].data
            if haskey(q_spaces[ci.idx].data, fock)
                all(collect(p_spaces[ci.idx][fock]) .!= collect(q_spaces[ci.idx][fock])) || error(" Not orthogonal")
            end
        end
    end

    #s = TuckerState(clusters)
    data = OrderedDict{FockConfig,OrderedDict{TuckerConfig,Array}}()
    for (fock,tconfig_list) in foi
        data2 = OrderedDict{TuckerConfig,Array}()
        for tconfig in tconfig_list
            data2[tconfig] = zeros((size(tconfig)...,nroots))
        end

        data[fock] = data2
    end
    return TuckerState(clusters, data, p_spaces, q_spaces)
#=}}}=#
end
"""
    +(ts1::FermiCG.TuckerState, ts2::FermiCG.TuckerState)
"""
function Base.:+(ts0::TuckerState, ts2::TuckerState)
#={{{=#
    ts1 = deepcopy(ts0)
    for (fock,configs) in ts2
        if haskey(ts1, fock)
            for (config,coeffs) in configs
                if haskey(ts1[fock], config)
                    ts1[fock][config] .+= ts2[fock][config]
                else
                    ts1[fock][config] = ts2[fock][config]
                end
            end
        else
            ts1[fock] = ts2[fock]
        end
    end
    return ts1
#=}}}=#
end

"""
    add!(ts1::FermiCG.TuckerState, ts2::FermiCG.TuckerState)

Add coeffs in `ts2` to `ts1`
"""
function add!(ts1::TuckerState, ts2::TuckerState)
#={{{=#
    for (fock,configs) in ts2
        if haskey(ts1, fock)
            for (config,coeffs) in configs
                if haskey(ts1[fock], config)
                    ts1[fock][config] .+= ts2[fock][config]
                else
                    ts1[fock][config] = ts2[fock][config]
                end
            end
        else
            ts1[fock] = ts2[fock]
        end
    end
#=}}}=#
end
"""
    dot(ts1::FermiCG.TuckerState, ts2::FermiCG.TuckerState)

Dot product between in `ts2` to `ts1`
"""
function dot(ts1::TuckerState, ts2::TuckerState)
#={{{=#

    unfold!(ts1)
    unfold!(ts2)
    nroots1 = nothing
    for (fock,configs) in ts1
        for (config,coeffs) in configs
            if nroots1 == nothing
                nroots1 = last(size(coeffs))
            else
                nroots1 == last(size(coeffs)) || error(" mismatch in number of roots")
            end
        end
    end
    nroots2 = nothing
    for (fock,configs) in ts2
        for (config,coeffs) in configs
            if nroots2 == nothing
                nroots2 = last(size(coeffs))
            else
                nroots2 == last(size(coeffs)) || error(" mismatch in number of roots")
            end
        end
    end

    overlap = zeros(nroots1,nroots2)
    for (fock,configs) in ts2
        haskey(ts1, fock) || continue
        for (config,coeffs) in configs
            haskey(ts1[fock], config) || continue
            overlap .+= ts1[fock][config]' * coeffs
        end
    end
    fold!(ts1)
    fold!(ts2)
    return overlap
#=}}}=#
end

"""
    scale!(ts::FermiCG.TuckerState, a)

Scale `ts` by a constant `a`
"""
function scale!(ts::TuckerState, a)
    #={{{=#
    for (fock,configs) in ts
        for (config,tuck) in configs
            ts[fock][config] .*= a
        end
    end
    #=}}}=#
end

"""
    unfold!(ts::TuckerState)
"""
function unfold!(ts::TuckerState)
#={{{=#
    for (fock,configs) in ts.data
        #display(fock)
        for (config,coeffs) in configs
            #length(size(coeffs)) == length(ts.clusters)+1 || error(" Can only unfold a folded vector")
            nroots = last(size(coeffs))

            #display(config)
            #display(size(coeffs))
            #display(size(config))
            #display(nroots)

            #display(fock)
            #display(config)
            ts[fock][config] = reshape(coeffs, (prod(size(config)),nroots))
        end
    end
#=}}}=#
end
"""
    fold!(ts::TuckerState)
"""
function fold!(ts::TuckerState)
#={{{=#
    for (fock,configs) in ts.data
        for (config,coeffs) in configs
            #length(size(coeffs)) == 2 || error(" Can only fold an unfolded vector")
            nroots = last(size(coeffs))

            ts[fock][config] = reshape(coeffs, (size(config)...,nroots))
        end
    end
#=}}}=#
end

"""
    randomize!(ts::TuckerState; scale=1)

Add some random noise to the vector
"""
function randomize!(ts::TuckerState; scale=1)
    #={{{=#
    for (fock,configs) in ts
        for (config,coeffs) in configs
            ts[fock][config] += scale .* (rand(size(coeffs)...) .- .5)
        end
    end
    #=}}}=#
end

"""
    mult!(ts::TuckerState, A)

Multiple `ts` by a matrix A. This is a multiplication over global state index
"""
function mult!(ts::TuckerState, A)
    #={{{=#
    unfold!(ts)
    for (fock,configs) in ts
        for (config,coeffs) in configs
            ts[fock][config] = coeffs * A
        end
    end
    fold!(ts)
    #=}}}=#
end

"""
    orthogonalize!(ts::TuckerState)

Symmetric Orthogonalization of vectors in ts
"""
function orthogonalize!(ts::TuckerState)
    #={{{=#
    S = dot(ts,ts)
    F = eigen(S)
    l = sqrt.(inv.(F.values))
    U = F.vectors
    for li in 1:length(l)
        U[:,li] .*= l[li]
    end
    mult!(ts,U)
    #=}}}=#
end

"""
    add_fockconfig!(s::ClusteredState, fock::FockConfig)
"""
function add_fockconfig!(s::TuckerState, fock::FockConfig)
    s.data[fock] = OrderedDict{TuckerConfig, Array}()
    #s.data[fock] = OrderedDict{TuckerConfig, Array}(TuckerConfig([1:1 for i in 1:length(s.clusters)])=>[0.0])
end

"""
    Base.length(s::TuckerState)
"""
function Base.length(s::TuckerState)
    l = 0
    for (fock,configs) in s.data
        for (config,vector) in configs
            l += dim(config)
        end
    end
    return l
end
"""
    prune_empty_fock_spaces!(s::TuckerState)

remove fock_spaces that don't have any configurations
"""
function prune_empty_fock_spaces!(s::TuckerState)
    focklist = keys(s.data)
    for fock in focklist
        if length(s.data[fock]) == 0
            delete!(s.data, fock)
        end
    end
    focklist = keys(s.data)
    #println(" Here's what's left")
    #for fock in keys(s.data)
    #end
    #    display(fock)
    #    display(s[fock])
    #end
end
"""
    get_vector(s::TuckerState)
"""
function get_vector(ts::TuckerState)
    nroots = nothing
    for (fock,configs) in ts
        for (config,coeffs) in configs
            if nroots == nothing
                nroots = last(size(coeffs))
            else
                nroots == last(size(coeffs)) || error(" mismatch in number of roots")
            end
        end
    end

    nroots != nothing || error(" Couldn't find nroots? length:", length(ts))

    v = zeros(length(ts), nroots)
    #println(length(ts), nroots)
    idx = 1
    for (fock, configs) in ts
        for (config, coeffs) in configs
            dims = size(coeffs)

            dim1 = prod(dims[1:end-1])
            v[idx:idx+dim1-1,:] = copy(reshape(coeffs,dim1,nroots))
            idx += dim1
        end
    end
    return v
end
"""
    set_vector!(s::TuckerState)
"""
function set_vector!(ts::TuckerState, v)

    length(size(v)) <= 2 || error(" Only takes matrices", size(v))
    nbasis = size(v)[1]
    nroots = 0
    if length(size(v)) == 1
        nroots = 1
    elseif length(size(v)) == 2
        nroots = size(v)[2]
    end

    #println(length(ts), nroots)
    idx = 1
    for (fock, tconfigs) in ts
        for (tconfig, coeffs) in tconfigs
            dims = size(tconfig)

            dim1 = prod(dims)
            #display(size(v[idx:idx+dim1-1,:]))
            ts[fock][tconfig] = copy(v[idx:idx+dim1-1,:])
            idx += dim1
        end
    end
    nbasis == idx-1 || error("huh?", nbasis, " ", idx)
    fold!(ts)
    return
end
"""
    zero!(s::TuckerState)
"""
function zero!(s::TuckerState)
    for (fock, configs) in s
        for (config, coeffs) in configs
            fill!(s[fock][config], 0.0)
        end
    end
end
"""
    eye!(s::TuckerState)
"""
function eye!(s::TuckerState)
    idx1 = 1
    idx2 = 1
    for (fock, configs) in s
        for (config, coeffs) in configs
            for config in product(config.config)
                coeffs[config,idx] = 1
                idx += 1
            end
        end
    end
end

"""
    Base.display(s::TuckerState; thresh=1e-3)

Pretty print
"""
function Base.display(s::TuckerState; root=1, thresh=1e-3)
#={{{=#
    println()
    @printf(" --------------------------------------------------\n")
    @printf(" ---------- State        -------------------: %5i  \n",root)
    @printf(" ---------- # Fockspaces -------------------: %5i  \n",length(keys(s.data)))
    @printf(" ---------- # Configs    -------------------: %5i  \n",length(s))
    @printf(" --------------------------------------------------\n")
    @printf(" Printing contributions greater than: %f", thresh)
    @printf("\n")
    @printf(" %-20s%-20s%-20s\n", "Weight", "# Configs", "(α,β)...")
    @printf(" %-20s%-20s%-20s\n", "-------", "---------", "----------")
    unfold!(s)
    for (fock,configs) in s.data
        prob = 0
        len = 0
        for (config, coeffs) in configs
            prob += coeffs[:,root]' * coeffs[:,root]
            len += length(coeffs[:,root])
        end
        if prob > thresh
            @printf(" %-20.3f%-20i", prob,len)
            for sector in fock
                @printf("(%2i,%-2i)", sector[1],sector[2])
            end
            println()

            @printf("     %-16s%-20s%-20s\n", "Weight", "", "Subspaces")
            @printf("     %-16s%-20s%-20s\n", "-------", "", "----------")
            for (config, coeffs) in configs
                probi = coeffs[:,root]' * coeffs[:,root]
                @printf("     %-16.3f%-20i", probi,length(coeffs))
                for range in config
                    @printf("%4s", range)
                end
                println()
            end
            println()
        end
    end
    fold!(s)
    print(" --------------------------------------------------\n")
    println()
#=}}}=#
end
"""
    print_fock_occupations(s::TuckerState; root=1, thresh=1e-3)

Pretty print
"""
function print_fock_occupations(s::TuckerState; root=1, thresh=1e-3)
#={{{=#

    println()
    unfold!(s)
    @printf(" --------------------------------------------------\n")
    @printf(" ---------- State        -------------------: %5i  \n",root)
    @printf(" ---------- # Fockspaces -------------------: %5i  \n",length(keys(s.data)))
    @printf(" ---------- # Configs    -------------------: %5i  \n",length(s))
    @printf(" --------------------------------------------------\n")
    @printf(" Printing contributions greater than: %f", thresh)
    @printf("\n")
    @printf(" %-20s%-20s%-20s\n", "Weight", "# Configs", "(α,β)...")
    @printf(" %-20s%-20s%-20s\n", "-------", "---------", "----------")
    sum = 0
    for (fock,configs) in s.data
        prob = 0
        len = 0
        for (config, coeffs) in configs
            prob += coeffs[:,root]' * coeffs[:,root]
            len += length(coeffs[:,root])
        end
        if prob > thresh
            @printf(" %-20.3f%-20i", prob,len)
            for sector in fock
                @printf("(%2i,%-2i)", sector[1],sector[2])
            end
            println()
        end
        sum += prob
    end
    print(" --------------------------------------------------\n")
    fold!(s)
    println()
#=}}}=#
end

"""
    expand_each_fock_space!(s::TuckerState, bases)

For each fock space sector defined, add all possible basis states
- `bases::Vector{ClusterBasis}`
"""
function expand_each_fock_space!(s::TuckerState, bases::Vector{ClusterBasis}; nroots=1)
    # {{{
    println("\n Make each Fock-Block the full space")
    # create full space for each fock block defined
    for (fblock,configs) in s.data
        #println(fblock)
        dims = TuckerConfig()
        #display(fblock)
        for c in s.clusters
            # get number of vectors for current fock space
            dim = size(bases[c.idx][fblock[c.idx]], 2)
            push!(dims, 1:dim)
        end

        s.data[fblock][dims] = zeros((length.(dims)...,nroots))

    end
end
# }}}

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script creating monthly data files
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, Dates

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

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

# initia DataFrame
df_VLAIO = DataFrame(PERIOD=Int64[], PRODUCT_NC=String[], PRODUCT_NC_LAB1=String[], PRODUCT_NC_LAB2=String[], DECLARANT_ISO=String[],
                TRADE_TYPE=String[], PARTNER_ISO=String[], FLOW=String[], VALUE_IN_EUROS=Union{Missing, Float64}[], QUANTITY_IN_KG=Union{Missing, Float64}[])

# EU27 in Dutch
EU27 = ["Roemenië", "Griekenland", "Oostenrijk", "Polen", "Duitsland", "Spanje", "Hongarije", "Slovakije", "Italië", "Nederland",
       "Frankrijk", "Letland", "Kroatië", "Cyprus", "Malta", "Litouwen", "Slovenië", "Estland", "Portugal", "Finland", "Tsjechië",
       "Luxemburg", "Zweden", "Denemarken", "Bulgarije", "Ierland", "België"]

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

for i in string.(2014:2021)

    path = dir_dropbox * "rawdata/VLAIO/" * "dbo.FactProductFlows_PD8_" * i * ".csv"
    df = CSV.read(path, DataFrame, header=false, delim=",", quoted=true)

    # same formatting as comext
    cols_names = ["PERIOD", "PRODUCT_NC", "PRODUCT_NC_LAB1","PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE",
                    "PARTNER_ISO", "FLOW", "VALUE_IN_EUROS", "QUANTITY_IN_KG", "SUP_QUANTITY", "SUP_UNIT"]
    rename!(df, cols_names)
    df = df[:, Not([:SUP_QUANTITY, :SUP_UNIT])]
    transform!(df, [:PARTNER_ISO, :DECLARANT_ISO, :TRADE_TYPE] .=> ByRow(string), renamecols=false)

    # missing values for zeros
    # Notes:
    #   - a lot of observations have entries for VALUES but only zeros for QUANTITY
    #   - for example in 2014: 25_304 of 97_708 (~25%) observations lost
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(x == 0.0, missing, x)), renamecols=false)

    g(x) = x[1:6]
    df.PERIOD = g.(string.(df.PERIOD))
    df.PERIOD = parse.(Int64, df.PERIOD)

    df.PRODUCT_NC = lpad.(string.(df.PRODUCT_NC), 8, '0')

    df.FLOW = ifelse.(df.FLOW .== "X", "exports", "imports")

    # append
    append!(df_VLAIO, df) # missing values are kept in output, essentially just formatted data

end

# export
CSV.write(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", df_VLAIO)

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

# # english names

# df_VLAIO = CSV.read(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", DataFrame)
# transform!(df_VLAIO, ["PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE", "PARTNER_ISO", "FLOW"] .=> ByRow(string), renamecols=false)
# transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> convert(Union{Missing, Float64}, x)), renamecols=false)
# df_VLAIO.PRODUCT_NC = lpad.(string.(df_VLAIO.PRODUCT_NC), 8, '0') # needs to be done again
# df_VLAIO.PARTNER_ISO .= "VL"

# df_names = DataFrame(XLSX.readtable(dir_home * "ctry_english_dutch" * ".xlsx", "Sheet1")...)
# transform!(df_ctry, names(df_ctry) .=> ByRow(string), renamecols=false)
# rename!(df_ctry, :lab_dutch => :PARTNER_ISO)

# df_VLAIO = leftjoin(df_VLAIO, df_names, on=:PARTNER_ISO)

# CSV.write(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", df_VLAIO)

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

ext: jl


text:
@testset "parse single subject" begin
    sub_path = joinpath(
        @__DIR__,
        "data",
        "bids_root",
        "sub-subtest"
    )

    my_sub = Subject(sub_path)

    @test my_sub.path == sub_path
    @test my_sub.identifier == "subtest"
    @test total_sessions(my_sub) == 1
    @test total_files(my_sub) == 1
    @test sprint(show, my_sub) == """
        Subject:
            identifier = subtest
            total session = 1
            total files = 1"""
end


text:
using OrdinaryDiffEq, DiffEqSensitivity, Zygote, Test

function loss1(p;sensealg=nothing)
  f(x,p,t) = [p[1]]
  prob = DiscreteProblem(f, [0.0], (1,10), p)
  sol = solve(prob, FunctionMap(scale_by_time = true), saveat=[1,2,3])
  return sum(sol)
end
dp1 = Zygote.gradient(loss1,[1.0])[1]
dp2 = Zygote.gradient(x->loss1(x,sensealg=TrackerAdjoint()),[1.0])[1]
dp3 = Zygote.gradient(x->loss1(x,sensealg=ReverseDiffAdjoint()),[1.0])[1]
dp4 = Zygote.gradient(x->loss1(x,sensealg=ForwardDiffSensitivity()),[1.0])[1]
@test dp1 == dp2
@test dp1 == dp3
@test dp1 == dp4


text:
# TODO Check faster implementations in computational methods for local regression
# William S. Cleveland & E. Grosse https://link.springer.com/article/10.1007/BF01890836
function ghat(x::T;
           A::AbstractMatrix{T},
           b::AbstractVector{T},
           d::Integer=2,
           q::Integer,
           rho::AbstractVector{T}) where T<:Real

    xv = A[:,d]
    yv = b

    ## λ_q
    n = length(xv)
    q = min(q,n)
    xvx = @. abs(xv-x)
    qidx = sortperm(xvx)[1:q]
    qdist = abs(xv[last(qidx)]-x)*max(1.0,q/n)

    ## upsilon
    w = zeros(n)
    for wi in qidx
        aq = abs(xv[wi]-x)/qdist
        w[wi] = max((1-aq^3)^3,0.0)
    end

    A = @. A*(w*rho)
    b = @. b*(w*rho)

    lsq_x = A\b

    d == 1 ? [x,1.0]'*lsq_x : [x^2.0,x,1.0]'*lsq_x

end

"""
Package: Forecast

    loess(xv, yv;
          d = 2,
          q = Int64(round(3/4*length(xv))),
          rho = repeat([1.0],inner=length(xv)),
          predict = xv)

Smooth a vector of observations using locally weighted regressions.

Although loess can be used to smooth observations for any given number of independent variables, this implementation is univariate. The speed of loess can be greatly increased by using fast aproximations for the linear fitting calculations, however this implementation calculates only exact results.

The loess functionality and nomenclature follows the descriptions in:

"STL: A Seasonal, Trend Decomposition Procedure Based on Loess"
Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning.
Journal of Official Statistics Vol. 6. No. 1, 1990, pp. 3-73 (c) Statistics Sweden.

# Arguments
- `xv`: Observations' support.
- `yv`: Observation values.
- `d`: Degree of the linear fit, it accepts values 1 or 2.
- `q`: As q increases loess becomes smoother, when q tends to infinity loess tends to an ordinary least square poynomial fit of degree `d`. It defaults to the rounding of 3/4 of xv's length.
    - `rho`: Weights expressing the reliability of the observations (e.g. if yi had variances σ^2*ki where ki where known, the rhoi could be 1/ki). It defaults to 1.0.
- `predict`: Vector containing the real values to be predicted, by default predicts xv.

# Returns
The loess values for the values contained in `predict`.

# Examples
```julia-repl
julia> loess(rand(5), rand(5); predict=rand(10))
10-element Array{Float64,1}:
[...]
```
"""
function loess(xv::AbstractVector{R},
               yv::AbstractVector{<:Union{Missing,S}};
               d::Integer=2,
               q::Integer=Int64(round(3/4*length(xv))),
               rho::AbstractVector{<:Union{Missing,T}}=fill(1.0,length(xv)),
               predict::AbstractVector{V} = xv) where {R<:Real, S<:Real, T<:Real, V<:Real}

    @assert (d==1) | (d==2) "Linear Regression must be of degree 1 or 2"
    @assert length(xv) == length(yv)

    myi = findall(x -> !ismissing(x),yv)
    xv = xv[myi]
    yv = yv[myi]
    rho = rho[myi]

    # Promote to same type
    P = promote_type(R,S,T,V)
    xv = R == P ? xv :  P.(xv)
    yv =  P.(yv)
    rho = P.(rho)
    predict = V == P ? predict : P.(predict)

    res = Vector{P}(undef, length(predict))

    ## Ax = b
    A = hcat(xv,fill(P(1),length(xv)))
    b = yv
    d == 2 && (A = hcat(xv .^ 2, A))

    for (i,xi) in enumerate(predict)
        res[i] = ghat(xi;A,b,d,q,rho)
    end

    res
end


text:
using TopOpt

nels = (20, 10, 10)
sizes = (1.0,1.0,1.0)
E = 1.0
ν = 0.3
force = -1.0
problem = PointLoadCantilever(Val{:Linear}, nels, sizes, E, ν, force)
solver1 = FEASolver(Displacement, Direct, problem)
# Takes forever to compile
# solver2 = FEASolver(Displacement, CG, MatrixFree, problem)
solver3 = FEASolver(Displacement, CG, Assembly, problem)

x0 = rand(length(solver1.vars))
solver1.vars .= x0
# solver2.vars .= x0
solver3.vars .= x0

solver1()
# solver2()
solver3()

#@test solver1.u ≈ solver2.u
@test solver1.u ≈ solver3.u


text:
import ImageView.display

type RasterImage{T,N,A<:AbstractArray} <: AbstractImageDirect{T,N}
    data::A
    properties::Dict
end

RasterImage(data::AbstractArray,props::Dict) = RasterImage{eltype(data),ndims(data),typeof(data)}(data,props)
RasterImage(raster::Raster) = RasterImage{eltype(raster.data),ndims(raster.data),typeof(raster.data)}(raster.data,["limits"=>(minimum(raster.data),maximum(raster.data)),"colorspace"=>"Gray","spatialorder"=>["y","x"],"transform"=>raster.transform,"projection"=>raster.projection])

display(raster::Raster) = display(RasterImage(raster),pixelspacing=[1,1])


text: S=SharedArray( Float64, (100, 100, 5), pids=[2,3,4,5]);

max_stars_repo_path: Module 2/Chapter07/shared_array.jl

max_issues_repo_path: Module 2/Chapter07/shared_array.jl

max_forks_repo_path: Module 2/Chapter07/shared_array.jl


text:
# Hold a topology in a format that should be close to what is available from
# network operators.
struct Topology
    name::AbstractString # Name (optional).
    nodes::AbstractMatrix{Any} # n×2 matrix: first column is a numerical index
    # (Int), second column is a node name.
    edges::AbstractMatrix{Float64} # m×3 matrix: first column is the index of
    # the edge source, second column is the index of the edge destination,
    # third column is the capacity of this link (no unit mandated; the user is
    # responsible for consistency).
    traffic::AbstractMatrix{Float64} # m×3 matrix: first column is the index
    # of the source node for this demand, second column is the index of the
    # destination node for this demand, third column is the actual demand
    # (no unit mandated; the user is responsible for consistency).
end

Topology(nodes::AbstractMatrix{Any}, edges::AbstractMatrix{Float64},
         traffic::AbstractMatrix{Float64}) = Topology("", nodes, edges, traffic)

function topology_to_graphs(t::Topology; make_network_undirected::Bool=true)
    # make_network_undirected: whether edges should be duplicated to make the
    # network effectively undirected.

    # Topology.
    g = MetaDiGraph()
    for i in 1:size(t.nodes, 1)
        add_vertex!(g, :name, t.nodes[i, 2])
    end
    for i in 1:size(t.edges, 1)
        add_edge!(g, Int(t.edges[i, 1]), Int(t.edges[i, 2]), :capacity, t.edges[i, 3])
        if make_network_undirected
            add_edge!(g, Int(t.edges[i, 2]), Int(t.edges[i, 1]), :capacity, t.edges[i, 3])
        end
    end

    # Demands.
    k = MetaDiGraph()
    dm = Dict{Edge{Int}, Float64}()
    sizehint!(dm, size(t.nodes, 1))

    for i in 1:size(t.nodes, 1)
        add_vertex!(k, :name, t.nodes[i, 2])
    end
    for i in 1:size(t.traffic, 1)
        ds = Int(t.traffic[i, 1])
        dt = Int(t.traffic[i, 2])
        add_edge!(k, ds, dt, :demand, t.traffic[i, 3]) # TODO: deprecate :demand here? So that all metadata from the graph is ignored. This should ease implementation of stochastic uncertainty.
        dm[Edge(ds, dt)] = t.traffic[i, 3]
    end

    return g, k, dm
end


text:
module IndexedDims

using AcceleratedArrays

export IndexedDimsArray, bypass

include("array.jl")  # IndexedDimsArray type
include("base.jl")
include("autoaccelerate.jl")
include("stdindexing.jl")  # a wrapper for non-value (regular) indexing (hypothetical)
include("dimindexing.jl")  # value indexing

end # module


text:
using Test
using MPI

MPI.Init()

function gather(A, root)
    comm = MPI.COMM_WORLD
    C = MPI.Gather(A, root, comm)
end

comm = MPI.COMM_WORLD
root = 0

for typ in Base.uniontypes(MPI.MPIDatatype)
    A = typ[MPI.Comm_rank(comm) + 1]
    C = gather(A, root)
    if MPI.Comm_rank(comm) == root
        @test typeof(C) === Vector{typ}
        @test C == map(typ, collect(1:MPI.Comm_size(comm)))
    end
    A = convert(typ,MPI.Comm_rank(comm) + 1)
    C = gather(A, root)
    if MPI.Comm_rank(comm) == root
        @test typeof(C) === Vector{typ}
        @test C == map(typ, collect(1:MPI.Comm_size(comm)))
    end
end

MPI.Finalize()


text:
"""
This is the main class for block ciphers.
"""

"""
    Cipher

A cipher is an algorithm used to translate plaintext into ciphertext.

# Attack models
Attack model - assumptions about what the attacker can and cannot do.

1. Ciphertext-only attackers (COA) observe ciphertexts but don't know
the associated plaintexts, and don't know how the plaintexts were
selected. Attackers in the COA model are passive and can't perform
encryption or decryption queries.
2. Known-plaintext attackers (KPA) observe ciphertexts and do know the
associated plaintexts. Attackers in the KPA model thus get a list of
plaintext–ciphertext pairs, where plaintexts are assumed to be
randomly selected. Again, KPA is a passive attacker model.
3. Chosen-plaintext attackers (CPA) can perform encryption queries for
plaintexts of their choice and observe the resulting ciphertexts. This
model captures situations where attackers can choose all or part of
the plaintexts that are encrypted and then get to see the ciphertexts.
Unlike COA or KPA, which are passive models, CPA are active
attackers, because they influence the encryption processes rather
than passively eavesdropping.
4. Chosen-ciphertext attackers (CCA) can both encrypt and decrypt; that
is, they get to perform encryption queries and decryption queries. CDs
could use this method since attackers can encrypt and decrypt but they want
the key so that they can distribute CDs.

# Security models
Security model - states that we consider a successful attack.

Intuitively, a cipher is secure if, even given a large number of plaintext–
ciphertext pairs, nothing can be learned about the cipher’s behavior when
applied to other plaintexts or ciphertexts.

In mathematical terms this means
1. Indistinguishability (IND) - Ciphertexts should be indistinguishable
from random strings. This is usually illustrated with this hypothetical
game: if an attacker picks two plaintexts and then receives a ciphertext
of one of the two (chosen at random), they shouldn’t be able to tell
which plaintext was encrypted, even by performing encryption queries
with the two plaintexts
2. Non-malleability (NM) - Given a ciphertext C1 = E(K, P1), it should be
impossible to create another ciphertext, C2, whose corresponding
plaintext, P2, is related to P1 in a meaningful way.

# Combining the two
Security notion - the combination of some security goal with some
attack model. We'll say that a cipher achieves a certain security notion if
any attacker working in a given model can't achieve the security goal.
The notion is described as GOAL-MODEL.

Semantic security - IND-CPA. That is, ciphertexts shouldn't leak any
information about plaintexts as long as the key is secret. To achieve IND-CPA
security, encryption must return different ciphertexts if called twice
on the same plaintext; otherwise, an attacker could identify duplicate
plaintexts from their ciphertexts, contradicting the definition that
ciphertexts shouldn't reveal any information.

Some relations
- IND-CCA implies IND-CPA
- NM-CCA implies NM-CPA
- NM-CPA implies IND-CPA

So encryption often is randomised with R, some fresh random bits
each time the system is called.

Kerckhoffs's principle - the security of a cipher should rely only on the secrecy
of the key and not on the secrecy of the cipher.

Entropy = - sum over x of P(X=x)*lb(P(X=x))
The binary logarithm expresses the information in bits and yields integer
values when probabilities are powers of two. Entropy is the level of surprise
of the outcome.

# Mersenne Twister
Why Mersenne Twister (MT) algorithm is insecure.
Its internal state is an array, S, consisting of 624 32-bit words.
This array is initially set to S1, S2, . . . , S624 and evolves to S2, . . . , S625,
then S3, . . . , S626, and so on, according to this equation:
S(k + 624) := S(k + 397) ⊕ A((Sk ∧ 0x80000000) ∨ (S(k + 1) ∧ 0xfffffff))
all instructions are bitwise and A is a function that transforms some
32-bit word, x, to (x >> 1), if x's most significant bit is 0,
or to (x >> 1) ⊕ 0x9908b0df otherwise.

Notice in this equation that bits of S interact with each other only
through XORs. The operators ∧ and ∨ never combine two bits of S
together, but just bits of S with bits from the constants 0x80000000 and
0x7fffffff. This way, any bit from S625 can be expressed as an XOR of bits
from S398, S1, and S2, and any bit from any future state can be expressed
as an XOR combination of bits from the initial state S1, . . . , S624.

Because there are exactly 624 × 32 = 19,968 bits in the initial state (or
624 32-bit words), any output bit can be expressed as an equation with at
most 19,969 terms (19,968 bits plus one constant bit). That’s just about
2.5 kilobytes of data. The converse is also true: bits from the initial state
can be expressed as an XOR of output bits.

A cipher is informationally
secure only if, even given unlimited computation time and memory, it
cannot be broken. Computational security views a cipher as
secure if it cannot be broken within a reasonable amount of time, and with
reasonable resources such as memory, hardware, budget, energy, and so
on.

Computational security is sometimes expressed in terms of two values
t and epsilon. We then say that a cryptographic scheme is (t, ε)-secure if an attacker
performing at most t operations - whatever those operations are - has a
probability of success that is no higher than ε. We can conclude that
a cipher with a key of n bits is at best (t, t/2^n)-secure, for any t between 1 and 2^n.

t-secure = (t, 1)-secure
n-bit security = 2^n-secure
Does not provide much information on the cost of the attack
because the attacker can break the cipher with less than expected
number of operations.

Confusion means that the input
(plaintext and encryption key) undergoes complex transformations, and
diffusion means that these transformations depend equally on all bits of
the input.
"""
abstract type Cipher end
encrypt(cipher, plaintext::String) = String(encrypt(cipher, Vector{UInt8}(plaintext)))
decrypt(cipher, ciphertext::String) = String(encrypt(cipher, Vector{UInt8}(ciphertext)))

"""
    OneTimePad

The one time pad is a symmetric block cipher and it is the most secure cipher
but the most pointless. It requires a key the same size as the plaintext.
It then produces a ciphertext the same length as the plaintext as well by
XORing the key and plaintext. Decryption works exactly the same.

Keys should only be used one time otherwise an attacker will be able to
determine `XOR(plaintext1, plaintext2)` by calculating `XOR(ciphertext1, ciphertext2)`.
Further, if an attacker can identify `plaintext1` then they can calculate `plaintext2`.

As long as the key is random, the ciphertext will appear completely random
because the XOR of a random key and a nonrandom plaintext appears random.
The attacker can only identify the plaintext's length.
"""
struct OneTimePad <: Cipher
    key::Vector{UInt8}
end

#using Base.xor
#xor(a::Array{Uint8, 1}, b::Array{Uint8, 1}) = xor()
encrypt(cipher::OneTimePad, plaintext::Vector{UInt8}) = xor.(cipher.key, plaintext)
decrypt(cipher::OneTimePad, ciphertext::Vector{UInt8}) = xor.(cipher.key, ciphertext)

struct FeiselScheme <: Cipher
    key::Vector{UInt8}
end

function encrypt(cipher::FeiselScheme, plaintext::Vector{UInt8}, block_size=64, rounds=15)
    left = plaintext[1:(block_size ÷ 2)]
    right = plaintext[(block_size ÷ 2 + 1):end]
    subkeys = rijndael_key_schedule(cipher.key, rounds)

    for i in 1:rounds
        temp = xor.(left, round_function(right, subkeys[i]))
        left = right
        right = temp
    end
end

"""
The cipher key used for encryption is 128, 192 or 256 bits long.

If you want to encrypt only a single block of exactly 16 bytes, use padding=false.

Without KeyExpansion, all rounds would use the same key, K, and
AES would be vulnerable to slide attacks.
Without AddRoundKey, encryption wouldn’t depend on the key;
hence, anyone could decrypt any ciphertext without the key.
SubBytes brings nonlinear operations, which add cryptographic
strength. Without it, AES would just be a large system of linear
equations that is solvable using high-school algebra.
Without ShiftRows, changes in a given column would never affect
the other columns, meaning you could break AES by building four
232-element codebooks for each column. (Remember that in a secure
block cipher, flipping a bit in the input should affect all the output
bits.)
Without MixColumns, changes in a byte would not affect any other
bytes of the state. A chosen-plaintext attacker could then decrypt any
ciphertext after storing 16 lookup tables of 256 bytes each that hold
the encrypted values of each possible value of a byte.

You should never encrypt blocks independently. This reveals patterns and
identical plaintext values because the ciphertext will be the same. This is
called electronic codebook and should be avoided. One should rather use
cipher block chaining.
"""
struct AES <: Cipher
    key::Vector{UInt8}
    rounds::Int64
end

function AES(key)
    n = length(key) * 8

    # Determine the number of rounds to run
    n == 128 ? rounds = 10 : nothing
    n == 192 ? rounds = 12 : nothing
    n == 256 ? rounds = 14 : nothing
    return AES(key, rounds)
end

function encrypt(cipher::AES, plaintext::Vector{UInt8}, padding=true)
    padding ? plaintext = pad(plaintext, 16) : nothing

    if length(plaintext) > 16
        result = encrypt(cipher, plaintext[1:16], false)  # should be xored with nonce
        for i in 17:16:16fld(length(plaintext), 16)
            input = xor.(plaintext[i:(i+15)], result[(i-16):(i-1)])
            append!(result, encrypt(cipher, input, false))
        end
        return result
    end

    rounds = cipher.rounds

    # Generate the roundkeys from the key
    subkeys = rijndael_key_schedule(cipher.key, rounds)
    roundkeys = [hcat(subkeys[i:(i+3)]...) for i in 1:4:length(subkeys)]  # transposable

    # Define the 4 central functions
    xorkeys(s, i) = xor.(s, roundkeys[i])
    subbytes(s) = sbox(s)
    function shiftrows(s)
        result = []
        n = size(s, 1)
        for (i, row) in enumerate(eachrow(s))
            # push!(result, vcat(row[(n-i+2):end], row[1:(n-i+1)]))
            push!(result, vcat(row[i:end], row[1:(i-1)]))
        end
        return hcat(result...)'
    end
    function mixcolumns(s)  # https://en.wikipedia.org/wiki/Rijndael_MixColumns
        ss = copy(s)
        for j in 1:size(s, 2)
            ss[1, j] = xor(gmul(0x02, s[1, j]), gmul(0x03, s[2, j]), s[3, j], s[4, j])
            ss[2, j] = xor(s[1, j], gmul(0x02, s[2, j]), gmul(0x03, s[3, j]), s[4, j])
            ss[3, j] = xor(s[1, j], s[2, j], gmul(0x02, s[3, j]), gmul(0x03, s[4, j]))
            ss[4, j] = xor(gmul(0x03, s[1, j]), s[2, j], s[3, j], gmul(0x02, s[4, j]))
        end
        return ss
    end

    # main algo
    s = reshape(plaintext, (4, 4))  # transposable
    s = xorkeys(s, 1)
    for i in 2:rounds
        s = subbytes(s)
        s = shiftrows(s)
        s = mixcolumns(s)
        s = xorkeys(s, i)
    end
    s = subbytes(s)
    s = shiftrows(s)
    s = xorkeys(s, rounds+1)

    return reshape(s, 16)  # transposable
end

"""
If you encrypted only a single block of exactly 16 bytes, use padding=false.
"""
function decrypt(cipher::AES, ciphertext::Vector{UInt8}, padding=true)
    if length(ciphertext) > 16
        result = decrypt(cipher, ciphertext[1:16], false)
        for i in 17:16:16fld(length(ciphertext), 16)
            input = decrypt(cipher, ciphertext[i:(i+15)], false)
            input = xor.(input, ciphertext[(i-16):(i-1)])
            append!(result, input)
        end

        padding ? result = unpad(result) : nothing
        return result
    end

    rounds = cipher.rounds

    # Generate the roundkeys from the key
    subkeys = rijndael_key_schedule(cipher.key, rounds)
    roundkeys = [hcat(subkeys[i:(i+3)]...) for i in 1:4:length(subkeys)]  # transposable

    # Define the 4 central functions
    xorkeys(s, i) = xor.(s, roundkeys[i])
    invsubbytes(s) = invsbox(s)
    function invshiftrows(s)
        result = []
        n = size(s, 1)
        for (i, row) in enumerate(eachrow(s))
            push!(result, vcat(row[(n-i+2):end], row[1:(n-i+1)]))
        end
        return hcat(result...)'
    end
    function invmixcolumns(s)  # https://en.wikipedia.org/wiki/Rijndael_MixColumns
        ss = copy(s)
        for j in 1:size(s, 2)
            ss[1, j] = xor(gmul(0x0e, s[1, j]), gmul(0x0b, s[2, j]), gmul(0x0d, s[3, j]), gmul(0x09, s[4, j]))
            ss[2, j] = xor(gmul(0x09, s[1, j]), gmul(0x0e, s[2, j]), gmul(0x0b, s[3, j]), gmul(0x0d, s[4, j]))
            ss[3, j] = xor(gmul(0x0d, s[1, j]), gmul(0x09, s[2, j]), gmul(0x0e, s[3, j]), gmul(0x0b, s[4, j]))
            ss[4, j] = xor(gmul(0x0b, s[1, j]), gmul(0x0d, s[2, j]), gmul(0x09, s[3, j]), gmul(0x0e, s[4, j]))
        end
        return ss
    end

    # main algo
    s = reshape(ciphertext, (4, 4))  # transposable
    s = xorkeys(s, rounds+1)
    s = invshiftrows(s)
    s = invsubbytes(s)
    for i in rounds:-1:2
        s = xorkeys(s, i)
        s = invmixcolumns(s)
        s = invshiftrows(s)
        s = invsubbytes(s)
    end
    s = xorkeys(s, 1)
    s = reshape(s, 16) # transposable

    padding ? s = unpad(s) : nothing
    return s
end


text:
function abc_West!(r)
    w=0
    while (!isborder(r,Nord) && isborder(r, West))
        move!(r, Nord)
        w = w + 1
    end
    if (!isborder(r, West))
        move!(r, West)
        while isborder(r, Sud)
            move!(r, West)
        end
    end
    while w > 0
        move!(r, Sud)
        w = w - 1
    end
end
function abc_Nord!(r)
    n=0
    while (!isborder(r,Ost) && isborder(r, Nord))
        move!(r, Ost)
        n = n + 1
    end
    if (!isborder(r, Nord))
        move!(r, Nord)
        while isborder(r, West)
            move!(r, Nord)
        end
    end
    while n > 0
        move!(r, West)
        n = n - 1
    end
end
function abc_Ost!(r)
    o=0
    while (!isborder(r,Sud) && isborder(r, Ost))
        move!(r, Sud)
        o = o + 1
    end
    if (!isborder(r, Ost))
        move!(r, Ost)
        while isborder(r, Nord)
            move!(r, Ost)
        end
    end
    while o > 0
        move!(r, Nord)
        o = o - 1
    end
end
function abc_Sud!(r)
    s=0
    while (!isborder(r,West) && isborder(r, Sud))
        move!(r, West)
        s = s + 1
    end
    if (!isborder(r, Sud))
        move!(r, Sud)
        while isborder(r, Ost)
            move!(r, Sud)
        end
    end
    while s > 0
        move!(r, Ost)
        s = s - 1
    end
end
function mark_cross!(r::Robot)
    a=0
    b=0
    while !isborder(r, Nord)
        move!(r, Nord)
        a = a + 1
    end
    abc_Nord!(r)
    if (isborder(r, Sud))
        putmarker!(r)
        while !isborder(r, Nord)
            move!(r, Nord)
            putmarker!(r)
        end
        while !isborder(r,Sud)
            move!(r,Sud)
        end
        abc_Sud!(r)
    end
    while !isborder(r, Sud)
        putmarker!(r)
        move!(r, Sud)
        b = b + 1
    end
    putmarker!(r)
    abc_Sud!(r)
    if isborder(r, Nord)
        putmarker!(r)
        while !isborder(r, Sud)
            move!(r, Sud)
            putmarker!(r)
        end
        while !isborder(r, Nord)
            move!(r, Nord)
        end
        abc_Nord!(r)
    end
    b = b - a
    while b > 0
        move!(r, Nord)
        b = b - 1
    end
    a = 0
    b = 0
    while !isborder(r, West)
        move!(r, West)
        a = a + 1
    end
    abc_West!
    if isborder(r,Ost)
        putmarker!(r)
        while !isborder(r,West)
            move!(r,West)
            putmarker!(r)
        end
        while !isborder(r, Ost)
            move!(r, Ost)
        end
        abc_Ost!(r)
    end
    while !isborder(r, Ost)
        putmarker!(r)
        move!(r, Ost)
        b = b + 1
    end
    putmarker!(r)
    abc_Ost!(r)
    if isborder(r,West)
        putmarker!(r)
        while !isborder(r, Ost)
            move!(r, Ost)
            putmarker!(r)
        end
        while !isborder(r, West)
            move!(r, West)
        end
        abc_West!(r)
    end
    b = b - a
    while b > 0
        move!(r, West)
        b = b - 1
    end
end


text:
module Algorithm101

    filename = [
        "merge_sort",
        "quick_sort",
        "binary_search",
        "tree_traverse",
        "binary_search_tree",
        "bfs",
        "dfs",
        "topological_sort",
        "monte_carlo_integral",
        "mcm",
        "mst",
    ]

    for f in filename
        include("$(f).jl")
    end

end # module


text:
using CSV, DataFrames, StatsPlots, StanSample, MCMCChains
using Distributions, Statistics, Random

ProjDir = @__DIR__

# Simulate the data

#Random.seed!(123)

N = 15
obspairs = ((30.0, 53), (35.0, 45), (40.0, 28), (45.0, 26), (50.0, 25))

ppu = []; quantity = []
for i in 1:N in
	obs = rand(obspairs, 1)
	append!(ppu, [rand(Normal(obs[1][1], 1.0))])
	append!(quantity, [rand(Normal(obs[1][2], 1.0))])
end

df = DataFrame()
df[!, :ppu] = ppu
df[!, :quantity] = quantity
df[!, :ppu_l] = log.(ppu)
df[!, :quantity_l] = log.(quantity)

CSV.write(joinpath(ProjDir, "data", "data.csv"), df)

opt_01 = "
data {
	int<lower=0> N;
	vector[N] q; // Outcome
	vector[N] p; // Predictor
}
parameters {
	real<lower=0> a;
	real b;
	real<lower=0> sigma;
}
model {
	vector[N] mu;            		// mu is a vector
	a ~ normal(2.0, 10.0);			// intercept positive
	b ~ normal(-2.0, 10.0);			// demand drops as price increases
	sigma ~ uniform(0, 50);
	mu = a + b * p;
	q ~ normal(mu, sigma);
}
generated quantities {
	vector[N] mu_pred;
	mu_pred = exp(a + b * p);
}
";

data = Dict(
  :N => size(df, 1),
  :q => df[:, :quantity_l],
  :p => df[:, :ppu_l]
)

p1 = scatter(df[:, :ppu], df[:, :quantity],
	xlabel="ppu", ylabel="quantity sold", leg=false)
p2 = scatter(data[:p], data[:q],
	xlabel="log(ppu)", ylabel="log(quantity sold)", leg=false)
plot(p1, p2, layout=(2, 1))
savefig(joinpath(ProjDir, "plots", "data.png"))

sm = SampleModel("opt_01", opt_01)

rc = stan_sample(sm, data=data)

if success(rc)
	chn = read_samples(sm; output_format=:mcmcchains)

	chns = set_section(chn,
	    Dict(
	      :parameters => ["a", "b", "sigma"],
	      :mu => ["mu_pred.$i" for i in 1:length(ppu)],
	      :internals => names(chn, [:internals])
	    )
	 )

	plot(chns)
	savefig(joinpath(ProjDir, "plots", "chains.png"))
	dfsa = DataFrame(chn)
	density(dfsa[:, Symbol("mu_pred.1")], label="mu_pred.1")
	for i in 2:3:N
		density!(dfsa[:, Symbol("mu_pred.$i")], label="mu_pred.$i")
	end
	savefig(joinpath(ProjDir, "plots", "density.png"))

	dfs = DataFrame(chns)
	describe(chns, sections=[:parameters, :mu]) |> display

	p = 3.4:0.01:4.0
	ps = collect(p)
	p1 = scatter(data[:p], data[:q],
		xlims=(3.3, 3.95), ylims=(3.1, 4.05),
		xlabel="ppu_l", ylabel="quantity_l sold", leg=false)
	a_mean = mean(dfs[:, :a])
	b_mean = mean(dfs[:, :b])
	mu_q = a_mean .+ b_mean .* p
	plot!(p1, p, mu_q)
	savefig(joinpath(ProjDir, "plots", "regression.png"))

	price_points_l = data[:p]
	plot(xlabel="Price", ylabel="Quantity", leg=false,
		xlims=(0, 100), ylims=(0, 100)
	)
	selected_rows = rand(1:size(dfs, 1), 75)
	for row in eachrow(dfs[selected_rows,:])
		qs = exp.(row[:a] .+ row[:b] .* (ps .- mean(ps)))
		plot!(exp.(p), qs, color=:gray)
	end
	qqs = exp.(a_mean .+ b_mean .* (price_points_l .- mean(ps)))
	scatter!(exp.(price_points_l), qqs)
	savefig(joinpath(ProjDir, "plots", "sampling.png"))

	scatter(dfs[:, :a], dfs[:, :b],
		xlabel="a", ylabel="b", leg=false)
	savefig(joinpath(ProjDir, "plots", "correlations.png"))

end


text:
# Author: Ashwani Rathee
#  Program to Check Whether Number is Even or Odd
println("Enter number to check if number is even or odd  : ")
number=readline()
number=parse(Int64,number)
check=number%2
if check == 1
    println("The Number is odd")
elseif check == 0
    println("The Number is even")
else
    println("Don't Know what this stuff is")
end

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
module TestExamples3DAdvection

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_3d_dgsem")

@testset "Linear scalar advection" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [0.00016263963870641478],
      linf = [0.0014537194925779984])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.00016017848135651983],
      linf = [0.0014175368788298393])
  end

  @trixi_testset "elixir_advection_extended.jl with initial_condition_sin" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [0.002647730309275237],
      linf = [0.02114324070353557],
      initial_condition=Trixi.initial_condition_sin)
  end

  @trixi_testset "elixir_advection_extended.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [7.728011630010656e-16],
      linf = [3.9968028886505635e-15],
      initial_condition=initial_condition_constant)
  end

  @trixi_testset "elixir_advection_extended.jl with initial_condition_linear_z and periodicity=false" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [3.007995700405795e-16],
      linf = [2.886579864025407e-15],
      initial_condition=Trixi.initial_condition_linear_z,
      boundary_conditions=Trixi.boundary_condition_linear_z, periodicity=false)
  end

  @trixi_testset "elixir_advection_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_mortar.jl"),
      l2   = [0.001810141301577316],
      linf = [0.017848192256602058])
  end

  @trixi_testset "elixir_advection_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"),
      l2   = [9.773852895157622e-6],
      linf = [0.0005853874124926162])
  end
end

end # module


text:
export BestResponsePolicy

struct BestResponsePolicy{E,S,A,X,P<:AbstractPolicy} <: AbstractCFRPolicy
    cfr_reach_prob::Dict{S,Vector{Pair{E,Float64}}}
    best_response_action_cache::Dict{S,A}
    best_response_value_cache::Dict{E,Float64}
    best_responder::X
    policy::P
end

"""
    BestResponsePolicy(policy, env, best_responder)

- `policy`, the original policy to be wrapped in the best response policy.
- `env`, the environment to handle.
- `best_responder`, the player to choose best response action.
"""
function BestResponsePolicy(policy, env, best_responder;)
    S = eltype(state_space(env))
    A = eltype(action_space(env))
    E = typeof(env)

    p = BestResponsePolicy(
        Dict{S,Vector{Pair{E,Float64}}}(),
        Dict{S,A}(),
        Dict{E,Float64}(),
        best_responder,
        policy,
    )

    e = copy(env)
    @assert e == env "The copy method doesn't seem to be implemented for environment: $env"
    @assert hash(e) == hash(env) "The hash method doesn't seem to be implemented for environment: $env"
    RLBase.reset!(e)  # start from the root!
    init_cfr_reach_prob!(p, e)
    p
end

function (p::BestResponsePolicy)(env::AbstractEnv)
    if current_player(env) == p.best_responder
        best_response_action(p, env)
    else
        p.policy(env)
    end
end

function init_cfr_reach_prob!(p, env, reach_prob = 1.0)
    if !is_terminated(env)
        if current_player(env) == p.best_responder
            push!(get!(p.cfr_reach_prob, state(env), []), env => reach_prob)

            for a in legal_action_space(env)
                init_cfr_reach_prob!(p, child(env, a), reach_prob)
            end
        elseif current_player(env) == chance_player(env)
            for (a, pₐ) in zip(action_space(env), prob(env))
                if pₐ > 0
                    init_cfr_reach_prob!(p, child(env, a), reach_prob * pₐ)
                end
            end
        else  # opponents
            for a in legal_action_space(env)
                init_cfr_reach_prob!(p, child(env, a), reach_prob * prob(p.policy, env, a))
            end
        end
    end
end

function best_response_value(p, env)
    get!(p.best_response_value_cache, env) do
        if is_terminated(env)
            reward(env, p.best_responder)
        elseif current_player(env) == p.best_responder
            a = best_response_action(p, env)
            best_response_value(p, child(env, a))
        elseif current_player(env) == chance_player(env)
            v = 0.0
            A, P = action_space(env), prob(env)
            @assert length(A) == length(P)
            for (a, pₐ) in zip(A, P)
                if pₐ > 0
                    v += pₐ * best_response_value(p, child(env, a))
                end
            end
            v
        else
            v = 0.0
            for a in legal_action_space(env)
                v += prob(p.policy, env, a) * best_response_value(p, child(env, a))
            end
            v
        end
    end
end

function best_response_action(p, env)
    get!(p.best_response_action_cache, state(env)) do
        best_action, best_action_value = nothing, typemin(Float64)
        for a in legal_action_space(env)
            # for each information set (`state(env)` here), we may have several paths to reach it
            # here we sum the cfr reach prob weighted value to find out the best action
            v = sum(p.cfr_reach_prob[state(env)]) do (e, reach_prob)
                reach_prob * best_response_value(p, child(e, a))
            end
            if v > best_action_value
                best_action, best_action_value = a, v
            end
        end
        best_action
    end
end

RLBase.update!(p::BestResponsePolicy, args...) = nothing

function RLBase.prob(p::BestResponsePolicy, env::AbstractEnv)
    if current_player(env) == p.best_responder
        onehot(p(env), action_space(env))
    else
        prob(p.policy, env)
    end
end

function RLBase.prob(p::BestResponsePolicy, env::AbstractEnv, action)
    if current_player(env) == p.best_responder
        action == best_response_action(p, env)
    else
        prob(p.policy, env, action)
    end
end


text:
# Main file for solving the kNN-based reweighted SAA problem of Bertsimas and Kallus

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi

const gurobi_env = Gurobi.Env()

# read these as inputs from the corresponding "instance file"
modRepNum = parse(ARGS[1])
degNum = parse(ARGS[2])
dataRepNum = parse(ARGS[3])
sampSizeNum = parse(ARGS[4])

#*========= MODEL INFORMATION ==========
const caseNum = 1
const paramsFile = "params_case" * string(caseNum) * "_knn.jl"
#*======================================

#*========= INCLUDE FILES ====================
include("maxParameters.jl")
include("randomSeeds.jl")

include(paramsFile)
# the next two data files are defined in paramsFile
include(modelFile)
include(dataFile)

include("genkNNScenarios.jl")
include("solveModel.jl")
include("estimateSolnCost.jl")
include("evaluateSoln.jl")
#*======================================

# set seed for reproducibility
srand(startingSeed)

# directory name for storing results
const baseDirName = "case" * string(caseNum) * "_knn/" * "mod_" * string(modRepNum) * "/" * "deg_" * string(degNum) * "/"
const subDirName = baseDirName * "rep_" * string(dataRepNum) * "/"
const subDirName2 = subDirName * "samp_" * string(sampSizeNum) * "/"
mkpath(subDirName2)

#*========= PRINT OPTIONS ====================
const storeResults = true

const infoFile = "ddsp.txt"
const modelDataFile = "model_data.txt"
const numSampFile = "num_samples.txt"

const knnObjFile = "knn_obj.txt"
const knnDDObjFile = "knn_ddobj.txt"
const knnTimeFile = "knn_time.txt"
const koptFile = "k_opt.txt"

const covRealFile = "covariate_obs.txt"
#*======================================

#*========= GENERATE MODEL PARAMETERS ====================

srand(randomSeeds_models[modRepNum])

covariate_mean, covariate_covMat, coeff_true = generateBaseDemandModel(numCovariates)

#*========= STORE RESULTS ====================
if(storeResults)

	# write details to text file, including some key details about the test instance
	details_file = baseDirName * infoFile
	open(details_file, "w") do f
		write(f,"case number: $caseNum \n")
		write(f,"numResources: $numResources \n")
		write(f,"numCustomers: $numCustomers \n")
		write(f,"model replicate number: $modRepNum \n")
		write(f,"degree: $(degree[degNum]) \n")
		write(f,"sample sizes: $numDataSamples \n")
		write(f,"demand_errors_scaling: $demand_errors_scaling \n")
		write(f,"numFolds: $numFolds \n")
		write(f,"minExponent: $minExponent \n")
		write(f,"maxExponent: $maxExponent \n\n")

		write(f,"randomSeeds_MC: $randomSeeds_MC \n")
		write(f,"randomSeeds_models: $randomSeeds_models \n")
		write(f,"randomSeeds_data: $randomSeeds_data \n")
	end

	mod_data_file = baseDirName * modelDataFile
	open(mod_data_file, "w") do f
		write(f,"covariate_mean = $covariate_mean \n")
		write(f,"covariate_covMat = $covariate_covMat \n")
		write(f,"coeff_true = $coeff_true \n")
	end

	samp_size_file = baseDirName * numSampFile
	open(samp_size_file, "w") do f
		for i = 1:length(numDataSamples)
			write(f,"$(numDataSamples[i]) \n")
		end
	end

end
#*============================================

#*========= GENERATE DATA REPLICATE ====================

srand(randomSeeds_data[dataRepNum])

covariate_obs = generateCovariateReal(numCovariates,covariate_mean,covariate_covMat)

demand_data, covariate_data = generateDemandData(numCovariates,numDataSamples[sampSizeNum],degree[degNum],covariate_mean,covariate_covMat,coeff_true)

#*========= STORE RESULTS ====================
if(storeResults)

	cov_real_file = subDirName * covRealFile
	open(cov_real_file, "w") do f
		write(f,"covariate_obs = $covariate_obs")
	end

end
#*============================================

#*========= SOLVE kNN-based SAA MODEL ====================

tic()

demand_scen_knn = generatekNNScenarios(demand_data,covariate_data,numFolds,minExponent,maxExponent,covariate_obs)

z_soln_knn, objDDknn = solveSAAModel(demand_scen_knn)

knnObjEstimates = estimateSolnQuality(z_soln_knn,covariate_obs,degree[degNum],numMCScenarios,numMCReplicates,coeff_true)

knnTime = toq()

#*========= STORE RESULTS ====================
if(storeResults)

	# write data to text file
	obj_est_file = subDirName2 * knnObjFile
	open(obj_est_file, "w") do f
		for i = 1:length(knnObjEstimates)
			write(f,"$(knnObjEstimates[i]) \n")
		end
	end

	ddobj_est_file = subDirName2 * knnDDObjFile
	open(ddobj_est_file, "w") do f
		write(f,"$objDDknn")
	end

	time_obj_file = subDirName2 * knnTimeFile
	open(time_obj_file, "w") do f
		write(f,"$knnTime")
	end

	kopt_file = subDirName2 * koptFile
	open(kopt_file, "w") do f
		write(f,"$(size(demand_scen_knn,1)) \n")
	end

end
#*============================================

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
function reshapesingleton(x::AbstractArray, dim::Integer)
	@assert maximum(size(x))==length(x) "Only one nonsingleton dimension allowed"
	sz = ones(Int,dim)
	sz[dim] = length(x)
	reshape(x,sz...)
end

#function loadswarm(id::AbstractString,freqFolder,consensusFolder="")
function loadswarm(id::AbstractString, freqPath::AbstractString, consensusPath::AbstractString="")
	d = load(freqPath)
	@assert ["codonFreqs", "positions", "coverage"] ⊆ keys(d) "$id.jld does not have the required fields."
	N = length(d["codonFreqs"])

	pos = pop!(d,"positions")
	cf = pop!(d,"codonFreqs")
	cov = pop!(d,"coverage")
	segmentInfo = pop!(d,"segmentInfo")

	# create annotated array - putting segments after each other
	# 4x4x4xP - codons x total number of positions in all segments
	# annotations
		# nucleotide in pos 3, 2, 1
		# codon
		# amino acid
		# position in segment
		# segment
		# coverage
		# consensus
		# annotations from codon frequency computations

	data = cat(cf...; dims=4)

	A = AnnotatedArray(data)

	# id
	# annotate!(A,:id,[id])
	annotate!(A,:id,[id])

	# nucleotides
	nucs = ['A','C','G','T']
	annotate!(A,:n3,copy(nucs))
	annotate!(A,:n2,reshapesingleton(copy(nucs),2))
	annotate!(A,:n1,reshapesingleton(copy(nucs),3))

	# codons
	codons = [ Symbol(String([n1,n2,n3])) for n3 in nucs, n2 in nucs, n1 in nucs ]
	annotate!(A,:codon,codons)

	# amino acids
	aa = codon2aa(codons)
	annotate!(A,:aminoacid,aa)

	# nucleotide position
	pos = map( x->reshapesingleton(x,4), pos )
	pos = cat(pos...; dims=4)
	annotate!(A,:position,pos)

	if length(segmentInfo)==1
		annotate!(A,:segment,[Symbol(segmentInfo[1][1])])
	else
		# if multiple segments, we need to annotate each position
		segments = map( x->repeat([Symbol(x[1])], x[2]), segmentInfo )
		segments = vcat(segments...)
		annotate!(A,:segment,reshapesingleton(segments,4))
	end

	# coverage
	cov = map( x->reshapesingleton(x,4), cov )
	cov = cat(cov...; dims=4)
	annotate!(A,:coverage,cov)

	# load consensus if supplied
	if !isempty(consensusPath)
		# consensus = loadfasta(joinpath(consensusFolder,"$(id)_consensus.fasta"))
		consensus = loadfasta(consensusPath)

		# check that the consensus sequence names and length agrees with the .jld info
		consensusInfo = [(x[1],length(x[2])) for x in consensus]
		@assert segmentInfo==consensusInfo "Mismatch between consensus and codon frequency files"

		# codon consensus
		seqs = map(x->x[2], consensus)
		codonConsensus = map( x->x*"NN", seqs ) # add two NN to the end (we are beyond the reference genome)
		codonConsensus = map(seq2codons, codonConsensus)
		codonConsensus = vcat(codonConsensus...)
		annotate!(A,:consensus,reshapesingleton(codonConsensus,4))
	end

	# other annoations from dict
	for (k,v) in d
		annotate!(A, Symbol(k), [v])
	end

	A
end

function loadswarm(ids::AbstractArray, freqPaths::AbstractArray, consensusPaths::AbstractArray=String[])
	@assert length(ids)==length(freqPaths)
	@assert isempty(consensusPaths) || length(freqPaths)==length(consensusPaths)
	isempty(consensusPaths) && (consensusPaths = repeat([""],length(freqPaths)))
	cat(map((id,f,c)->loadswarm(id,f,c),ids,freqPaths,consensusPaths)...; dims=5)
end

function loadswarm(ids::AbstractArray,freqFolder::AbstractString,consensusFolder::AbstractString="")
	freqPaths = map(id->joinpath(freqFolder,"$id.jld"), ids)
	if isempty(consensusFolder)
		loadswarm(ids, freqPaths)
	else
		loadswarm(ids, freqPaths, map(id->joinpath(consensusFolder,"$(id)_consensus.fasta"), ids))
	end
end


text:
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.BalanceLaws: vars_state, Prognostic, Auxiliary
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates: flattenednames
using ClimateMachine.Ocean.SplitExplicit01
using ClimateMachine.GenericCallbacks
using ClimateMachine.VTK
using ClimateMachine.Checkpoint

using Test
using MPI
using Random
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

## using FileIO
using JLD2
using QuickVizExample
using GLMakie

function spex01_get_checkpoints(idir="jld-inputs")
        jp(x)=joinpath(idir,x)
        bcfiles=map( jp, filter(x -> occursin("baroclinic_checkpoint"  ,x), readdir( idir ) )
                   )
        btfiles=map( jp, filter(x -> occursin("barotropic_checkpoint"  ,x), readdir( idir ) )
                   )
        gfile=map(   jp, filter(  x -> occursin("spex_01_grid",x), readdir( idir ) )
                 )
        if length(gfile) == 1
         gfile=gfile[1]
        end
 return gfile,bcfiles,btfiles
end

function spex01_lookat(
                       file_list=nothing,
                       grid_file=nothing;
                       is_baroclinic=false,
                       is_barotropic=false,
                       plot_points=nothing,
                       np_def=10
                       )

         if !is_baroclinic && !is_barotropic
           # Figure out from file names
           if any(map( x->occursin("baroclinic_checkpoint",x),file_list ))
             is_baroclinic=true
           end
           if any(map( x->occursin("barotropic_checkpoint",x),file_list ))
             is_barotropic=true
           end
         end
         # Check we are only doing 1 of barotropic and baorclinic
         if is_baroclinic && is_barotropic
           # Can't do both barotropic and baroclinic rules are same time
           println("Can't plot both barotropic and baroclinic at once")
           return
         end

         # Read and process grid file - create grids for 3d xyz (baroclinic) and 2d xy
         # (barotropic).
         jld_grid=jldopen(grid_file)
         xr,yr,zr=jld_grid["xrange"],jld_grid["yrange"],jld_grid["zrange"]
         Nel=jld_grid["N"]
         close(jld_grid)
         t3d=StackedBrickTopology( MPI.COMM_WORLD,(xr,yr,zr) )
         g3d=DiscontinuousSpectralElementGrid(t3d,
                                              FloatType=Float64,
                                              DeviceArray=ClimateMachine.array_type(),
                                              polynomialorder=Nel )
         t2d=BrickTopology( MPI.COMM_WORLD,(xr,yr) )
         g2d=DiscontinuousSpectralElementGrid(t2d,
                                              FloatType=Float64,
                                              DeviceArray=ClimateMachine.array_type(),
                                             polynomialorder=Nel )
         gh3=GridHelper(g3d)
         gh2=GridHelper(g2d)
         x3,y3,z3=coordinates(g3d)
         x2,y2,z2=coordinates(g2d)

         # Create interpolation grid
         if plot_points == nothing
          xnew=range(xr[1],xr[end],length=np_def)
          ynew=range(yr[1],yr[end],length=np_def)
          znew=range(zr[1],zr[end],length=np_def)
         else
          xnew=plot_points[1]
          ynew=plot_points[2]
          znew=plot_points[3]
         end
         if is_barotropic
          znew=range(0,0,length=1)
         end

         # Load pickup data, interpolate and put into states for plotting
         npf=length(file_list)
         if is_barotropic
          fld1 = zeros(length(xnew), length(ynew), npf)
          fld2 = zeros(length(xnew), length(ynew), npf)
          fld3 = zeros(length(xnew), length(ynew), npf)
          ϕ = ScalarField(copy(x2), gh2)
          nf=1
          for fn in file_list
           jlf=jldopen(fn)
           # U[1],U[2],η
           ϕ .= view(jlf["h_Q"],:,1,:)
           fld1[:,:,nf] .= view(ϕ(xnew,ynew,znew),:,:,1)
           ϕ .= view(jlf["h_Q"],:,2,:)
           fld2[:,:,nf] .= view(ϕ(xnew,ynew,znew),:,:,1)
           ϕ .= view(jlf["h_Q"],:,3,:)
           fld3[:,:,nf] .= view(ϕ(xnew,ynew,znew),:,:,1)
           close(jlf)
           nf=nf+1
          end
          states=[fld1, fld2, fld3]
          statenames=["U", "V", "η"]
          plot_head="Barotropic "
         end

         if is_baroclinic
          fld1 = zeros(length(xnew), length(ynew), length(znew))
          fld2 = zeros(length(xnew), length(ynew), length(znew))
          fld3 = zeros(length(xnew), length(ynew), length(znew))
          fld4 = zeros(length(xnew), length(ynew), length(znew))
          ϕ = ScalarField(copy(x3), gh3)
          fn=file_list[1]
          jlf=jldopen(fn)
          # u[1],u[2],η,θ
          ϕ .= view(jlf["h_Q"],:,1,:)
          fld1[:,:,:] .= view(ϕ(xnew,ynew,znew),:,:,:)
          ϕ .= view(jlf["h_Q"],:,2,:)
          fld2[:,:,:] .= view(ϕ(xnew,ynew,znew),:,:,:)
          ϕ .= view(jlf["h_Q"],:,3,:)
          fld3[:,:,:] .= view(ϕ(xnew,ynew,znew),:,:,:)
          ϕ .= view(jlf["h_Q"],:,4,:)
          fld4[:,:,:] .= view(ϕ(xnew,ynew,znew),:,:,:)
          plot_tim=@sprintf("%g",jlf["t"])
          close(jlf)
          states=[fld1, fld2, fld3, fld4 ]
          statenames=["u, t="*plot_tim, "v, t="*plot_tim, "η, t="*plot_tim, "θ, t="*plot_tim ]
          plot_head="Baroclinic "
         end

         volumeslice(states, statenames = statenames, title=plot_head)

         return
end

gfile,bcfiles,btfiles=spex01_get_checkpoints("/Users/chrishill/projects/QuickVizExample.jl/jld-inputs")
spex01_lookat(btfiles,gfile;np_def=100)
# spex01_lookat([bcfiles[30]],gfile;np_def=100)

max_stars_repo_name: christophernhill/jmc-channel-plot

max_issues_repo_name: christophernhill/jmc-channel-plot

max_forks_repo_name: christophernhill/jmc-channel-plot


text:
import Base

abstract type Layer end

function (layer::Layer)(x...)
    outputs = forward(layer, x...)
    return outputs
end

function cleargrads!(layer::Layer)
    for param in parameters(layer)
        cleargrad!(param)
    end
end


text:
import Pkg;

offline_registry_name = Pkg.TOML.parsefile(
    joinpath(splitpath(@__DIR__)..., "Registry.toml")
    )["name"]
offline_registry_uuid = Pkg.TOML.parsefile(
    joinpath(splitpath(@__DIR__)..., "Registry.toml")
    )["uuid"]
pushfirst!(Base.DEPOT_PATH,joinpath(splitpath(@__DIR__)..., "depot",),);
pushfirst!(
    Base.DEPOT_PATH,
    joinpath(
        homedir(),
        ".julia.isolated",
        offline_registry_name,
        offline_registry_uuid,
        ),
    );
filter!((x) -> !(Base.Filesystem.samefile(expanduser(x),expanduser("~/.julia")) || lowercase(strip(abspath(expanduser(x)))) == lowercase(strip(abspath(expanduser("~/.julia")))) || Base.Filesystem.samefile(x,"~/.julia") || lowercase(strip(abspath(x))) == lowercase(strip(abspath("~/.julia"))) || Base.Filesystem.samefile(expanduser(x),expanduser(joinpath(homedir(), ".julia"))) || lowercase(strip(abspath(expanduser(x)))) == lowercase(strip(abspath(expanduser(joinpath(homedir(), ".julia"))))) || Base.Filesystem.samefile(x,joinpath(homedir(), ".julia")) || lowercase(strip(abspath(x))) == lowercase(strip(abspath(joinpath(homedir(), ".julia"))))), Base.DEPOT_PATH);
unique!(Base.DEPOT_PATH);

if !isfile(
        joinpath(
            homedir(),
            ".julia.isolated",
            offline_registry_name,
            offline_registry_uuid,
            "environments",
            string("v", VERSION.major, ".", VERSION.minor,),
            "Project.toml",
            )
        )
    mkpath(
        joinpath(
            homedir(),
            ".julia.isolated",
            offline_registry_name,
            offline_registry_uuid,
            "environments",
            string("v", VERSION.major, ".", VERSION.minor,),
            )
        );
    touch(
        joinpath(
            homedir(),
            ".julia.isolated",
            offline_registry_name,
            offline_registry_uuid,
            "environments",
            string("v", VERSION.major, ".", VERSION.minor,),
            "Project.toml",
            )
        );
end

Pkg.activate(
    joinpath(
        homedir(),
        ".julia.isolated",
        offline_registry_name,
        offline_registry_uuid,
        "environments",
        string("v", VERSION.major, ".", VERSION.minor,),
        )
    );

max_forks_repo_name: aluthge/OfflineRegistry


text:
using ImageTransformations: center
using Photometry
using Distributions
using Statistics
using HCIToolbox: get_annulus_segments
using StaticArrays
using Rotations

"""
    detectionmap([method=snr], data, fwhm; fill=0)

Parallel implementation of arbitrary detection mapping applied to each pixel in the input image. Any invalid values will be set to `fill`.

The following methods are provided in the [`Metrics`](@ref) module:
* [`snr`](@ref) - signal-to-noise ratio (S/N) using student-t statistics to account for small sample penalty.
* [`significance`](@ref) - Gaussian signifance using student-t statistics to account for samll sample penalty.

!!! tip
    This code is automatically multi-threaded, so be sure to set `JULIA_NUM_THREADS` before loading your runtime to take advantage of it!
"""
function detectionmap(method, data::AbstractMatrix{T}, fwhm; fill=zero(T)) where T
    out = fill!(similar(data), fill)
    width = minimum(size(data)) / 2 - 1.5 * fwhm

    masked = get_annulus_segments(data, fwhm/2 + 2, width, mode=:apply)
    coords = findall(!iszero, masked)

    Threads.@threads for coord in coords
        val = method(data, coord, fwhm)
        @inbounds out[coord] = isfinite(val) ? val : fill
    end

    return out
end

detectionmap(data, fwhm) = detectionmap(snr, data, fwhm)

"""
    snr(data, position, fwhm)

Calculate the signal to noise ratio (SNR, S/N) for a test point at `position` using apertures of diameter `fwhm` in a residual frame.

Uses the method of Mawet et al. 2014 which includes penalties for small sample statistics. These are encoded by using a student's t-test for calculating the SNR.

!!! note
    SNR is not equivalent to significance, use [`significance`](@ref) instead
"""
function snr(data::AbstractMatrix, position, fwhm)
    x, y = position
    cy, cx = center(data)
    separation = sqrt((x - cx)^2 + (y - cy)^2)
    separation > fwhm / 2 + 1 || return NaN
    r = fwhm / 2

    N = floor(Int, 2 * π  * separation / fwhm)
    dθ = 2π / N
    R = RotMatrix{2}(dθ)

    # initial points
    rx = x - cx
    ry = y - cy

    fluxes = similar(data, N)

    for idx in eachindex(fluxes)
        ap_x = rx + cx
        ap_y = ry + cy
        rx, ry = R * SA[rx, ry]
        ap = CircularAperture(ap_x, ap_y, r)
        fluxes[idx] = photometry(ap, data).aperture_sum
    end

    other_elements = @view fluxes[2:end]
    bkg_σ = std(other_elements) # ddof = 1 by default
    return (fluxes[1] - mean(other_elements)) / (bkg_σ * sqrt(1 + 1/(N - 1)))
end

snr(data::AbstractMatrix, idx::CartesianIndex, fwhm) = snr(data, reverse(idx.I), fwhm)

"""
    significance(data, position, fwhm)

Calculates the Gaussian significance from the signal-to-noise ratio (SNR, S/N) for a test point at `position` using apertures of diameter `fwhm` in a residual frame.

The Gaussian signifiance is calculated from converting the SNR confidence limit from a student t distribution to a Gaussian via

``\\text{sig}(\\text{SNR}) = \\Phi^{-1}\\left[\\int_0^\\text{SNR}{t_\\nu(x)dx}\\right]``

where the degrees of freedom ``\\nu`` is given as ``2\\pi r / \\Gamma - 2`` where r is the radial distance of each pixel from the center of the frame.

# See Also
[`snr`](@ref)
"""
function significance(data::AbstractMatrix, position, fwhm)
    x, y = position
    cy, cx = center(data)
    separation = sqrt((x - cx)^2 + (y - cy)^2)
    _snr = snr(data, position, fwhm)
    # put in one line to allow broadcast fusion
    return snr_to_sig(_snr, separation, fwhm)
end
significance(data::AbstractMatrix, idx::CartesianIndex, fwhm) = significance(data, reverse(idx.I), fwhm)

function snr_to_sig(snr, separation, fwhm)
    dof = 2 * π * separation / fwhm - 2
    dof > 0 || return NaN
    return quantile(Normal(), cdf(TDist(dof), float(snr)))
end
function sig_to_snr(sig, separation, fwhm)
    dof = 2 * π * separation / fwhm - 2
    dof > 0 || return NaN
    return quantile(TDist(dof), cdf(Normal(), float(sig)))
end

"""
    noise(data, position, fwhm)

Calculate the statistical noise for a test point at `position` using apertures of diameter `fwhm` in a residual frame.

Uses the standard deviation of the apertures in the entire annulus. This is distinct from the [`snr`](@ref) noise calculation, which defines a confidence interval using student-t statistics. This means you cannot simply create a noise map and divide it from the signal to create an equivalent S/N map.
"""
function noise(data::AbstractMatrix, position, fwhm)
    x, y = position
    cy, cx = center(data)
    separation = sqrt((x - cx)^2 + (y - cy)^2)
    separation > fwhm / 2 + 1 || return NaN
    r = fwhm / 2

    N = floor(Int, 2 * π  * separation / fwhm)
    dθ = 2π / N
    R = RotMatrix{2}(dθ)

    # initial points
    rx = x - cx
    ry = y - cy

    fluxes = similar(data, N)

    for idx in eachindex(fluxes)
        ap_x = rx + cx
        ap_y = ry + cy
        rx, ry = R * SA[rx, ry]
        ap = CircularAperture(ap_x, ap_y, r)
        fluxes[idx] = photometry(ap, data).aperture_sum
    end

    return std(fluxes)
end

noise(data::AbstractMatrix, idx::CartesianIndex, fwhm) = noise(data, reverse(idx.I), fwhm)


text:
using Codify, Test, FieldDefaults, FieldMetadata, Unitful

import FieldMetadata: flattenable, @flattenable

@flattenable @default_kw struct SubTest
    x::Float64 | 3.0 | false
    y::Float64 | 7.0 | true
end

@default_kw struct TestStruct{T}
    a::Float64 | 2.0
    b::Any     | 5.0
    c::T       | nothing
end

t = TestStruct(1.0, 1.0u"m*s^-1",SubTest(8.0,9.0))
s = Codify.codify_inner(typeof(t))
s = codify(t)
@test s == "TestStruct(\n    a = 1.0,\n    b = 1.0u\"m*s^-1\",\n    c = SubTest(\n        y = 9.0,\n    ),\n)"

t2 = eval(Meta.parse(s))

@test typeof(t) == typeof(t2)
@test t.a == t2.a
@test t.b == t2.b
@test t.c.y == 9.0
@test t2.c.x == 3.0


text:
"""
This is a sample input file for the CloudScat code.  To use it run julia
and type

julia> include("cloud_geometry.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia cloud_geometry.jl

- - -
2019 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat
using StaticArrays
using LinearAlgebra
using CoordinateTransformations

const co = CloudScat.constants

function run()
    params = Params(
        # NUmber of simulated photons
        N = 600000,

        # The wavelength of the scattered light
        λ = 337 * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 8 * co.kilo],
        source_b = [0, 0, 13 * co.kilo])

    # World geometry
    # Describe the cloud geometry here.
    # The geometry is specified using several elementary shapes
    # that can be deformed and combined.  The geometry primitives are defined
    # in src/geometry.jl
    #
    # ELEMENTARY SHAPES
    # =================
    #
    # Currently the supported elementary shapes are:
    # - A sphere:
    #   Sphere(xcenter, ycenter, zcenter, radius)
    # - A vertical, finite cylinder:
    #   Cylinder(zbottom, ztop, xcenter, ycenter, radius)
    # - A vertical, finite cone (possibly a double-cone):
    #   Cone(zbottom, ztop, xvertex, yvertex, zvertex, m).
    #   m is the slope of the cone, i.e. the cone surface satisfies
    #   sqrt((x - xv)^2 + (y - yv)^2) = m (z - zv).
    # - A half-space delimited by a horizontal plane:
    #   Plane(z, up)
    #   If up is true, the half-space is above the plane, else it is below the
    #   plane
    #
    # COMBINATIONS
    # ============
    #
    # - Union of several shapes:
    #   union(shape1, shape2, ...)
    #   Alternatively, (shape1 + shape2 + ...) would also work.
    # - Intersection:
    #   intersect(shape1, shape2, ...)
    # - Difference
    #   shapediff(shape1, shape2)
    #   shape1 - shape2 would also work.
    #
    # TRANSFORMATIONS
    # ===============
    #
    # Any linear transformation in the package CoordinateTransformations is
    # supported.  See https://github.com/FugroRoames/CoordinateTransformations.jl
    # To use these transformations use ]add CoordinateTransformations in the julia
    # prompt.
    # Rotations from the Rotations package (]add CoordinateTransformations)
    # are also supported when embedded in a LinearMap.
    # For performance reasons, use StaticArrays whenever possible.

    # As an example we will build a geometry composed by a wide cylinder topped
    # by a conical "turret" which is in turn topped by an ellipsoid and has
    # a small semi-spherical "hole" on top.

    # First a 20-km wide cylinder:
    cyl = Cylinder(7 * co.kilo,  10 * co.kilo, 0, 0, 20 * co.kilo)

    # Now a cone on the axis
    cone = Cone(9 * co.kilo, 15 * co.kilo, 0, 0, 25 * co.kilo,
                0.75)

    # Finally to define the ellipsiod we start with a sphere centered on the
    # origin with radius 1 km.
    sph = Sphere(0, 0, 0, 1 * co.kilo)

    # We deform it to an ellipsoid with axes (12 km, 10 km, 2 km) and move it to
    # (0, 0, 15 km)
    map = AffineMap(SMatrix{3, 3, Float64}(Diagonal([12, 10, 2])),
                    SVector{3, Float64}([0, 0, 15 * co.kilo]))
    ellipsoid = TransformedShape(map, sph)

    # Now the hole that we will substract from the whole, slightly shifted in the
    # y direction.
    hole = Sphere(0, 1 * co.kilo, 17 * co.kilo, 3 * co.kilo)

    # Finally we combine everything together
    cloud = shapediff(union(cyl, cone, ellipsoid), hole)

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Cylinder(7 * co.kilo, 60 * co.kilo, 0, 0, 200 * co.kilo)

    # Set a homogeneous cloud droplet density and radius.
    # See variable_droplet_radius.jl for how to deal with inhomogeneous clouds.
    # Radius of the scattering particles
    radius = 10e-6

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, radius)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    # Observers record a time curve and an image of received light.  You can
    # set up as many observers as you wish but they have a significant impact
    # on performance.

    observers = [Observer(
        # Location of the observer
        position = [200 * co.kilo, 200 * co.kilo, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    viewpoint = "{" * join(observers[1].r, ", ") * "}"
    viewdir = "{" * join([observers[1].r[1], observers[1].r[2], 0.0],  ", ") * "}"
    println("To visualize the cloud geometry with Mathematica use:")
    println()
    println("R = ", mathematica(cloud))
    println("""RegionPlot3D[R, PlotTheme -> "Classic", Axes -> True,
        PlotStyle -> Directive[Opacity[0.75], RGBColor[.6, .65, .85],
            Specularity[0.2]],
        Lighting -> "Neutral", PlotPoints -> 100,
        ViewVector -> {$viewpoint, $viewdir}]
    """)
    println()
    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                   saveto="$(basename).h5")
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()


text:
"""
function attribute_weighting(data::ClustData,attribute_weights::Dict{String,Float64})

apply the different attribute weights based on the dictionary entry for each tech or exact name
"""
function attribute_weighting(data::ClustData,
                              attribute_weights::Dict{String,Float64}
                              )
  for name in keys(data.data)
    tech=split(name,"-")[1]
    if name in keys(attribute_weights)
      attribute_weight=attribute_weights[name]
      data.data[name].*=attribute_weight
      data.sdv[name]./=attribute_weight
    elseif tech in keys(attribute_weights)
      attribute_weight=attribute_weights[tech]
      data.data[name].*=attribute_weight
      data.sdv[name]./=attribute_weight
    end
  end
  return data
end


text:
@testset "ZeroOne" begin
    @test @inferred(classify(0.2f0, LabelEnc.ZeroOne(0.3))) === 0.
    @test @inferred(classify(0.4f0, LabelEnc.ZeroOne(0.3))) === 1.
    @test classify.([0.2,0.4,0.6], LabelEnc.ZeroOne(0.3)) == [0,1,1]
    @test eltype(classify.([0.2,0.4,0.6], LabelEnc.ZeroOne(0.3))) <: Float64
    @test @inferred(classify([0.2,0.4,0.6], LabelEnc.ZeroOne(0.3))) == [0,1,1]
    @test eltype(classify([0.2,0.4,0.6], LabelEnc.ZeroOne(0.3))) <: Float64
    for T in (Float16, Float32, Float64)
        @test @inferred(classify(T(0.4), 0.5)) === zero(T)
        @test @inferred(classify(T(0.5), 0.5)) === one(T)
        @test @inferred(classify(T(0.4), 0.3)) === one(T)
        @test @inferred(classify(T(0.4), LabelEnc.ZeroOne)) === zero(T)
        @test @inferred(classify(T(0.6), LabelEnc.ZeroOne)) === one(T)
        # broadcast
        @test classify.(T[0.4,0.6], 0.5) == [0,1]
        @test classify.(T[0.4,0.6], LabelEnc.ZeroOne) == [0,1]
        @test @inferred(classify(T[0.4,0.6], 0.5)) == [0,1]
        @test @inferred(classify(T[0.4,0.6], LabelEnc.ZeroOne)) == [0,1]
        @test eltype(classify.(T[0.4,0.6], 0.5)) <: T
        @test eltype(classify.(T[0.4,0.6], LabelEnc.ZeroOne)) <: T
        @test eltype(classify(T[0.4,0.6], 0.5)) <: T
        @test eltype(classify(T[0.4,0.6], LabelEnc.ZeroOne)) <: T
        buffer = zeros(2)
        @test @inferred(classify!(buffer, [0.4,0.6], 0.5)) == [0,1]
        @test buffer == [0,1]
        buffer = zeros(2)
        @test @inferred(classify!(buffer, [0.4,0.6], LabelEnc.ZeroOne)) == [0,1]
        @test buffer == [0,1]
    end
    for T in (Float16, Float32, Float64, UInt8, Int32, Int64)
        @test @inferred(classify(0.2f0, LabelEnc.ZeroOne(T,0.3))) === zero(T)
        @test @inferred(classify(0.4f0, LabelEnc.ZeroOne(T,0.3))) === one(T)
        @test @inferred(classify(0.2f0, LabelEnc.ZeroOne(T))) === zero(T)
        @test @inferred(classify(0.4f0, LabelEnc.ZeroOne(T))) === zero(T)
        @test @inferred(classify(0.6f0, LabelEnc.ZeroOne(T))) === one(T)
        # broadcast
        @test @inferred(classify([0.4,0.6], LabelEnc.ZeroOne(T))) == [0,1]
        @test classify.([0.4,0.6], LabelEnc.ZeroOne(T)) == [0,1]
        @test classify.([0.2,0.4,0.6], LabelEnc.ZeroOne(T,0.3)) == [0,1,1]
        @test eltype(classify([0.4,0.6], LabelEnc.ZeroOne(T))) <: T
        @test eltype(classify.([0.4,0.6], LabelEnc.ZeroOne(T))) <: T
        @test eltype(classify.([0.2,0.4,0.6], LabelEnc.ZeroOne(T,0.3))) <: T
        buffer = zeros(2)
        @test @inferred(classify!(buffer, [0.4,0.6], LabelEnc.ZeroOne(T))) == [0,1]
        @test buffer == [0,1]
    end
end

@testset "MarginBased" begin
    for T in (Float16, Float32, Float64)
        @test @inferred(classify(T(0.2),  LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(0.6),  LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(0.0),  LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(-0.1), LabelEnc.MarginBased)) === -one(T)
        @test @inferred(classify(T(-5.1), LabelEnc.MarginBased)) === -one(T)
        @test classify.(T[0,-.1,0.2,3,-4], LabelEnc.MarginBased) == [1,-1,1,1,-1]
        @test @inferred(classify(T[0,-.1,0.2,3,-4], LabelEnc.MarginBased)) == [1,-1,1,1,-1]
        @test eltype(classify.(T[0,-.1,0.2,3,-4], LabelEnc.MarginBased)) <: T
        @test eltype(classify(T[0,-.1,0.2,3,-4], LabelEnc.MarginBased)) <: T
        buffer = zeros(5)
        @test @inferred(classify!(buffer, T[0,-.1,0.2,3,-4], LabelEnc.MarginBased)) == [1,-1,1,1,-1]
        @test buffer == [1,-1,1,1,-1]
    end
    for T in (Int16, Int32, Int64)
        @test @inferred(classify(T(12), LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(2),  LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(0),  LabelEnc.MarginBased)) === one(T)
        @test @inferred(classify(T(-1), LabelEnc.MarginBased)) === -one(T)
        @test @inferred(classify(T(-5), LabelEnc.MarginBased)) === -one(T)
        @test classify.(T[0,-1,2,3,-4], LabelEnc.MarginBased) == [1,-1,1,1,-1]
        @test @inferred(classify(T[0,-1,2,3,-4], LabelEnc.MarginBased)) == [1,-1,1,1,-1]
        @test eltype(classify.(T[0,-1,2,3,-4], LabelEnc.MarginBased)) <: T
        @test eltype(classify(T[0,-1,2,3,-4], LabelEnc.MarginBased)) <: T
    end
    for T in (Float16, Float32, Float64, Int16, Int32, Int64)
        @test @inferred(classify(12,  LabelEnc.MarginBased(T))) === one(T)
        @test @inferred(classify(2f0, LabelEnc.MarginBased(T))) === one(T)
        @test @inferred(classify(0,   LabelEnc.MarginBased(T))) === one(T)
        @test @inferred(classify(-1., LabelEnc.MarginBased(T))) === -one(T)
        @test @inferred(classify(-5,  LabelEnc.MarginBased(T))) === -one(T)
        @test classify.([0,-.1,0.2,3,-4], LabelEnc.MarginBased(T)) == [1,-1,1,1,-1]
        @test @inferred(classify([0,-.1,0.2,3,-4], LabelEnc.MarginBased(T))) == [1,-1,1,1,-1]
        @test eltype(classify.([0,-.1,0.2,3,-4], LabelEnc.MarginBased(T))) <: T
        @test eltype(classify([0,-.1,0.2,3,-4], LabelEnc.MarginBased(T))) <: T
    end
end

@testset "OneOfK" begin
    for T in (Float16, Float32, Float64)
        @test @inferred(classify(T[0.4,0.2,0.9,.1], LabelEnc.OneOfK)) === 3
        @test @inferred(classify(T[0.4,0.2,0.9,.1], LabelEnc.OneOfK{T})) === 3
        @test @inferred(classify(T[0.4,0.2,0.9,.1], LabelEnc.OneOfK(Val{4}))) === 3
    end
    for T in (Int16, Int32, Int64)
        @test @inferred(classify(T[4,2,9,10,2], LabelEnc.OneOfK)) === 4
        @test @inferred(classify(T[4,2,9,10,2], LabelEnc.OneOfK{T})) === 4
        @test @inferred(classify(T[4,2,9,10,2], LabelEnc.OneOfK(Val{4}))) === 4
    end
    for T in (Float16, Float32, Float64, Int16, Int32, Int64)
        A = T[2 1 0 2 1 0 1;
              1 0 2 1 2 1 2;
              0 2 1 0 0 2 0]
        @test @inferred(classify(A, LabelEnc.OneOfK)) == [1,3,2,1,2,3,2]
        @test @inferred(classify(A, LabelEnc.OneOfK, ObsDim.Last())) == [1,3,2,1,2,3,2]
        @test @inferred(classify(A, LabelEnc.OneOfK, obsdim=2)) == [1,3,2,1,2,3,2]
        @test @inferred(classify(A, LabelEnc.OneOfK(Val{3}))) == [1,3,2,1,2,3,2]
        @test @inferred(classify(A, LabelEnc.OneOfK(Val{3}), ObsDim.Last())) == [1,3,2,1,2,3,2]
        @test @inferred(classify(A, LabelEnc.OneOfK(Val{3}), obsdim=2)) == [1,3,2,1,2,3,2]
        @test eltype(classify(A, LabelEnc.OneOfK)) <: Int
        @test eltype(classify(A, LabelEnc.OneOfK, ObsDim.Last())) <: Int
        @test eltype(classify(A, LabelEnc.OneOfK, obsdim=2)) <: Int
        @test eltype(classify(A, LabelEnc.OneOfK(Val{3}))) <: Int
        @test eltype(classify(A, LabelEnc.OneOfK(Val{3}), ObsDim.Last())) <: Int
        @test eltype(classify(A, LabelEnc.OneOfK(Val{3}), obsdim=2)) <: Int
        buffer = zeros(7)
        @test @inferred(classify!(buffer, A, LabelEnc.OneOfK(3), ObsDim.Last())) == [1,3,2,1,2,3,2]
        @test buffer == [1,3,2,1,2,3,2]
        buffer = zeros(7)
        @test @inferred(classify!(buffer, A, LabelEnc.OneOfK(3); obsdim=:last)) == [1,3,2,1,2,3,2]
        @test buffer == [1,3,2,1,2,3,2]
        At = A'
        @test @inferred(classify(At, LabelEnc.OneOfK, ObsDim.First())) == [1,3,2,1,2,3,2]
        @test @inferred(classify(At, LabelEnc.OneOfK, obsdim=1)) == [1,3,2,1,2,3,2]
        @test @inferred(classify(At, LabelEnc.OneOfK(Val{3}), ObsDim.First())) == [1,3,2,1,2,3,2]
        @test @inferred(classify(At, LabelEnc.OneOfK(Val{3}), obsdim=1)) == [1,3,2,1,2,3,2]
        @test eltype(classify(At, LabelEnc.OneOfK, ObsDim.First())) <: Int
        @test eltype(classify(At, LabelEnc.OneOfK, obsdim=1)) <: Int
        @test eltype(classify(At, LabelEnc.OneOfK(Val{3}))) <: Int
        @test eltype(classify(At, LabelEnc.OneOfK(Val{3}), ObsDim.First())) <: Int
        @test eltype(classify(At, LabelEnc.OneOfK(Val{3}), obsdim=1)) <: Int
        buffer = zeros(7)
        @test @inferred(classify!(buffer, At, LabelEnc.OneOfK(3), ObsDim.First())) == [1,3,2,1,2,3,2]
        @test buffer == [1,3,2,1,2,3,2]
        buffer = zeros(7)
        @test @inferred(classify!(buffer, At, LabelEnc.OneOfK, obsdim=1)) == [1,3,2,1,2,3,2]
        @test buffer == [1,3,2,1,2,3,2]
    end
end


text:
import LocalCoverage
import Pkg
using TestReports

testdir = "test-reports"
cd(mkpath(testdir))

Pkg.test(p::Symbol; kwargs...) = TestReports.test(String(p); kwargs...)
LocalCoverage.pkgdir(p::Symbol) = LocalCoverage.pkgdir(String(p))
LocalCoverage.generate_coverage(:AdaStress)


text:
using Test
using Stripeline

@testset "Instrument database" begin
    include("instrumentdb_tests.jl")
end

@testset "Scanning strategy" begin
    include("scanning.jl")
end

@testset "Map maker" begin
    include("map_maker_tests.jl")
end

@testset "TOD splitter" begin
    include("tod_splitter_tests.jl")
end

@testset "Noise generation" begin
    include("noisegeneration_tests.jl")
end

@testset "ADC simulation" begin
    include("adc_tests.jl")
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Revise
using OrdinaryDiffEq, Flux, DiffEqFlux, LinearAlgebra, Plots
using DiffEqSensitivity, JLD2
const EIGEN_EST = Ref(0.0f0)
const USE_GPU = Ref(false)
#USE_GPU[] = false # the network is small enough such that CPU works just great
USE_GPU[] = true

USE_GPU[] && (using CuArrays)

_gpu(arg) = USE_GPU[] ? gpu(arg) : cpu(arg)
_cu(arg) = USE_GPU[] ? cu(arg) : identity(arg)

function getops(grid, T=Float32)
    N, dz = length(grid), step(grid)
    d = ones(N-2)
    dl = ones(N-3) # super/lower diagonal
    zv = zeros(N-2) # zero diagonal used to extend D* for boundary conditions

    # D1 first order discretization of ∂_z
    D1 = diagm(-1 => -dl, 0 => d)
    D1[1, :] .= 0
    D1[end, :] .= 0

    # D2 discretization of ∂_zz
    D2 = diagm(-1 => dl, 0 => -2*d, 1 => dl)
    κ = 0.05
    D2[1, 1] = D2[end, end] = -1
    D2 = (κ/(dz^2)).*D2 #add the constant κ as the equation requires and finish the discretization
    display(D2), display(D1)

    EIGEN_EST[] = maximum(abs, eigvals(D2))
    D1 = _cu(D1)
    D2 = _cu(D2)
    return (D1=T.(D1), D2=T.(D2))
end

file = jldopen("../DataGeneration/rayleigh_taylor_instability_3d_horizontal_averages.jld2")

Is = keys(file["timeseries/t"])

N = file["grid/Nz"]
Lz = file["grid/Lz"]
Nt = length(Is)

t = Float32.(zeros(Nt))
soldata = Float32.(zeros(Nt, N))

for (i, I) in enumerate(Is)
    t[i] = file["timeseries/t/$I"]
    soldata[i, :] .= file["timeseries/b/$I"]
end

grid = range(0, 1, length = N)
tspan = (t[1], t[end])
u0 = _cu(soldata[1,2:end-1])
ops = getops(grid)

ann = Chain(Dense(N-2,N-2,tanh), Dense(N-2,N-2,tanh), Dense(N-2,N-2,tanh),
            Dense(N-2,N-2,tanh), Dense(N-2,N-2,tanh)) |> _gpu
pp = param(Flux.data(DiffEqFlux.destructure(ann)))
lyrs = Flux.params(pp)

function dudt(u::TrackedArray,p,t)
    Φ = DiffEqFlux.restructure(ann, p)
    return ops.D1*Φ(u) + ops.D2*u
end
function dudt(u::AbstractArray,p,t)
    Φ = DiffEqFlux.restructure(ann, p)
    return ops.D1*Tracker.data(Φ(u)) + ops.D2*u
end

predict_adjoint() = diffeq_adjoint(pp,
                              prob,
                              ROCK4(eigen_est = (integ)->integ.eigen_est = EIGEN_EST[]),
                              u0=u0, saveat = saveat,
                              reltol=1e-5, abstol=1e-6,
                              # no back solve
                              sensealg=SensitivityAlg(quad=false, backsolve=false))

function loss_adjoint()
    pre = predict_adjoint()
    sum(abs2, training_data - pre)
end

cb = function ()
    arr = Array(training_data)
    cur_pred = collect(Flux.data(predict_adjoint()))
    n = size(arr, 1)
    pl = scatter(1:n,arr[:,10],label="data", legend =:bottomright, title = "10th time over space")
    scatter!(pl,1:n,cur_pred[:,10],label="prediction")
    pl2 = scatter(saveat,arr[end÷2,:],label="data", legend =:bottomright, title = "middle point over time")
    scatter!(pl2,saveat,cur_pred[end÷2,:],label="prediction")
    #display(plot(pl, pl2, size=(600, 300)))
    display(loss_adjoint())
end

prob = ODEProblem{false}(dudt,u0,tspan,pp)
nn = 1
saveat = t
soldata = soldata'
soldata = soldata[2:end-1, :]
training_data = _cu(soldata)
epochs = Iterators.repeated((), 20)
learning_rate = ADAM(0.01)
Flux.train!(loss_adjoint, lyrs, epochs, learning_rate, cb=cb)
#=
learning_rate = ADAM(0.001)
epochs = Iterators.repeated((), 300)
Flux.train!(loss_adjoint, lyrs, epochs, learning_rate, cb=cb)
@time loss_adjoint()
=#