text:
#### Interface to Clustering package

struct LMCLUSResult <: ClusteringResult
    manifolds::Vector{Manifold}
    separations::Vector{Separation}
end

nclusters(R::LMCLUSResult) = length(R.manifolds)

counts(R::LMCLUSResult) = map(size, R.manifolds)

"Return points-to-cluster assignments"
function assignments(R::LMCLUSResult)
    Ms = manifolds(R)
    A = zeros(Int, sum(map(m->size(m), Ms)))
    for (i,m) in enumerate(Ms)
        A[points(m)] .= i
    end
    return A
end

"""Get linear manifold cluster"""
manifold(R::LMCLUSResult, idx::Int) = R.manifolds[idx]
manifolds(R::LMCLUSResult) = R.manifolds

function setdim!(R::LMCLUSResult, d::Int)
    for i in 1:nclusters(R)
        R.manifolds[i].d = d
    end
end

"""Iterative refinement of LMCLUS clustering

    refine(res::LMCLUSResult, data::AbstractMatrix, dfun::Function, efun::Function; tol::Real = 10.0, maxiter::Integer = 100)

Updates assignments to linear manifold clusters `res` on dataset `data` by evaluating similarity from dataset points to each cluster by

    dfun::(AbstractMatrix, Manifold) -> Vector{<:Real}

with following evaluation of the refined clustering by

    efun::(AbstractMatrix, Manifold) -> Real

Function returns refined `LMCLUSResult` object, after algorithm converges with tolerance, `tol`, or after specified number of iterations, `maxiter`.
"""
function refine(res::LMCLUSResult, data::AbstractMatrix,
                dfun::Function, efun::Function; bounds = false,
                tol::Real = 10.0, maxiter::Integer = 100, debug = false,
                drop_last = true, min_cluster_size = 20)
    M = manifolds(res)[1:(drop_last ? end-1 : end)]
    Δ = efun(data, M)

    # main loop
    c = 0
    converged = false
    while !converged && c < maxiter
        c += 1

        D = map(m->size(m) > 0 ? dfun(data, m) : fill(Inf, size(data,2)), M)
        DM = hcat(D...)
        A = mapslices(d->last(findmin(d)), DM, dims=2)
        Alidx = LinearIndices(A)

        # update assignments
        M⁺ = Manifold[]
        for (i,m) in enumerate(M)
            Idxs = Alidx[findall(A .== i)]
            if length(Idxs) < min_cluster_size #TODO: reassign points before modifing cluster parameters
                for ii in Idxs
                    A[ii] = sortperm(DM[ii,:])[2]
                end
                continue # reassign cluster points, effectivly removing it
            end
            if length(Idxs) > 0
                X = data[:, Idxs]
                R = fit(PCA, X; pratio=1.0)
                μ = mean(R)
                B = projection(R)
                m = Manifold(outdim(m), mean(R), projection(R), Idxs)
                if bounds
                    m.θ = maximum(distance_to_manifold(X, m))
                    m.σ = maximum(distance_to_manifold(X, m, ocss=true))
                end
                push!(M⁺, m)
            end
        end
        filter!(m->size(m) != 0, M⁺)

        # evaluate new clustering
        Δ⁺ = sum(efun(data, m) for m in M⁺)
        debug && println("$c: Δ = ($Δ - $Δ⁺) = $(Δ - Δ⁺)")
        M = M⁺

        # check convergence
        converged = abs(Δ - Δ⁺) < tol
        Δ = Δ⁺
    end

    return LMCLUSResult(M, Separation[])
end

"""Reassign outliers to the closes cluster in the clustering."""
function clearoutliers(res::LMCLUSResult, data::AbstractMatrix,
                       dfun::Function; debug = false)
    M = manifolds(res)
    outdim(M[end]) != 0 && return res
    O = pop!(M)

    # evaluate distances to outliers
    D = map(m->dfun(data, m), M)
    OI = points(O)
    DM = hcat(D...)[OI, :]

    # reassing outliers to available clusters
    A = mapslices(d->last(findmin(d)), DM, dims=2)
    for (i,a) in enumerate(A)
        push!(points(M[a]), OI[i])
    end

    # refine parameters of manifolds
    for m in M
        X = data[:, points(m)]
        R = fit(PCA, X; maxoutdim = outdim(m))
        m.μ = mean(R)
        m.basis = projection(R)
        m.θ = 0.0
        m.σ = 0.0
    end

    return LMCLUSResult(M, Separation[])
end

points(cl::C, i::Int) where {C <: ClusteringResult} = findall(i .∈  assignments(cl))

function LMCLUSResult(cl::C, data::AbstractMatrix, d=0) where {C <: ClusteringResult}
    R = LMCLUSResult([Manifold(data, points(cl, i)) for i in 1:nclusters(cl)], Separation[])
    d > 0 && setdim!(R, d)
    return R
end


text:
# time ./julia day08.jl -> 255ms when using a single type to represent
# all operations:
# struct Instr
#  op :: Symbol
#  num :: Int64
# end
#
# Increases to 347ms when using separate types for
# each op (OpNop, OpJmp, OpAcc) and Union{..} for an instruction.

lines = open("../inputs/day08.txt") do f
    readlines(f)
end

# Instructions
struct OpNop
    num::Int64
end
struct OpJmp
    num::Int64
end
struct OpAcc
    num::Int64
end

isNop(i::OpNop) = true
isNop(_) = false

isJmp(i::OpJmp) = true
isJmp(_) = false

isAcc(i::OpAcc) = true
isAcc(_) = false

Instr = Union{OpNop,OpJmp,OpAcc}

# Possible terminations
struct Terminated
    acc::Int64
end

isTerminated(t::Terminated) = true
isTerminated(_) = false

struct Seen
    acc::Int64
end

# Parse

function parseLine(ln::String)::Instr
    parts = split(ln, " ")
    num = parse(Int64, parts[2])
    if parts[1] == "nop"
        return OpNop(num)
    elseif parts[1] == "jmp"
        return OpJmp(num)
    else
        @assert parts[1] == "acc"
        return OpAcc(num)
    end
end

function run(program::Vector{Instr})::Union{Terminated,Seen}
    acc = 0
    pc = 1
    seen = Set{Int64}()
    term_pc = length(program) + 1
    while pc ∉ seen && pc != term_pc
        push!(seen, pc)
        op = program[pc]
        if isNop(op)
            pc += 1
            continue
        elseif isJmp(op)
            pc += op.num
            continue
        else
            @assert isAcc(op)
            acc += op.num
            pc += 1
        end
    end
    if pc == term_pc
        return Terminated(acc)
    else
        return Seen(acc)
    end
end

part1(program) = run(program).acc

function part2(program)::Int64
    for i in 1:length(program)
        op = program[i]
        if isJmp(op)
            newOp = OpNop(op.num)
        elseif isNop(op)
            newOp = OpJmp(op.num)
        else
            continue
        end
        program[i] = newOp
        code = run(program)
        if isTerminated(code)
            return code.acc
        end
        program[i] = op
    end
end

function main()
    program::Vector{Instr} = parseLine.(lines)
    println("Part 1: ", part1(program))
    println("Part 2: ", part2(program))
end

main()

ext: jl


text:
@compat abstract type Optimizer end
function print_header(method::Optimizer)
    @printf "Iter     Function value   Gradient norm \n"
end

immutable Options{TCallback <: Union{Void, Function}}
    x_tol::Float64
    f_tol::Float64
    g_tol::Float64
    f_calls_limit::Int
    g_calls_limit::Int
    h_calls_limit::Int
    allow_f_increases::Bool
    iterations::Int
    store_trace::Bool
    show_trace::Bool
    extended_trace::Bool
    show_every::Int
    callback::TCallback
    time_limit::Float64
end
function Options(;
        x_tol::Real = 1e-32,
        f_tol::Real = 1e-32,
        g_tol::Real = 1e-8,
        f_calls_limit::Int = 0,
        g_calls_limit::Int = 0,
        h_calls_limit::Int = 0,
        allow_f_increases::Bool = false,
        iterations::Integer = 1_000,
        store_trace::Bool = false,
        show_trace::Bool = false,
        extended_trace::Bool = false,
        show_every::Integer = 1,
        callback = nothing,
        time_limit = NaN)
    show_every = show_every > 0 ? show_every : 1
    #if extended_trace && callback == nothing
    #    show_trace = true
    #end
    Options{typeof(callback)}(
        Float64(x_tol), Float64(f_tol), Float64(g_tol), f_calls_limit, g_calls_limit, h_calls_limit,
        allow_f_increases, Int(iterations), store_trace, show_trace, extended_trace,
        Int(show_every), callback, time_limit)
end
function print_header(options::Options)
    if options.show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end
end

immutable OptimizationState{T <: Optimizer}
    iteration::Int
    value::Float64
    g_norm::Float64
    metadata::Dict
end
function Base.show(io::IO, t::OptimizationState)
    @printf io "%6d   %14e   %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

@compat OptimizationTrace{T} = Vector{OptimizationState{T}}
function Base.show(io::IO, tr::OptimizationTrace)
    @printf io "Iter     Function value   Gradient norm \n"
    @printf io "------   --------------   --------------\n"
    for state in tr
        show(io, state)
    end
    return
end

@compat abstract type OptimizationResults end
type UnivariateOptimizationResults{T,O<:Optimizer} <: OptimizationResults
    method::O
    initial_lower::T
    initial_upper::T
    minimizer::T
    minimum::T
    iterations::Int
    iteration_converged::Bool
    converged::Bool
    rel_tol::T
    abs_tol::T
    trace::OptimizationTrace{O}
    f_calls::Int
end
type MultivariateOptimizationResults{O<:Optimizer,T,N} <: OptimizationResults
    method::O
    initial_x::Array{T,N}
    minimizer::Array{T,N}
    minimum::T
    iterations::Int
    iteration_converged::Bool
    x_converged::Bool
    x_tol::Float64
    x_residual::Float64
    f_converged::Bool
    f_tol::Float64
    f_residual::Float64
    g_converged::Bool
    g_tol::Float64
    g_residual::Float64
    f_increased::Bool
    trace::OptimizationTrace{O}
    f_calls::Int
    g_calls::Int
    h_calls::Int
end
function Base.append!(a::MultivariateOptimizationResults, b::MultivariateOptimizationResults)
    a.iterations += iterations(b)
    a.minimizer = minimizer(b)
    a.minimum = minimum(b)
    a.iteration_converged = iteration_limit_reached(b)
    a.x_converged = x_converged(b)
    a.f_converged = f_converged(b)
    a.g_converged = g_converged(b)
    append!(a.trace, b.trace)
    a.f_calls += f_calls(b)
    a.g_calls += g_calls(b)
end


text:
txt = (readlines("day2.txt") .|> split)
hor = 0
aim = 0
depth = 0

for i in 1:size(txt, 1)
    X = (txt[i][2] |> x -> parse(Int64,x))
    if txt[i][1] == "forward"
        hor += X
        depth += aim*X
    else aim -= (txt[i][1] == "down" ? -X : X)
    end
end

print(hor*depth)


text:
function write_dataset(x_train::Matrix, y_train::Matrix, x_test::Matrix, y_test::Matrix, windowsize::String)
    # Train dataset
    open("x_train_"*windowsize*".txt", "w") do io
        print(io, x_train)
    end

    open("y_train_"*windowsize*".txt", "w") do io
        print(io, y_train)
    end

    # Test dataset
    open("x_test_"*windowsize*".txt", "w") do io
        print(io, x_test)
    end

    open("y_test_"*windowsize*".txt", "w") do io
        print(io, y_test)
    end
end


text:
### TEST pflap

## using grid

# generate image GRID
n = 20
mgrid = ones(512, 512);
@views mgrid[1:n:end, :] .= 0;
@views mgrid[:, 1:n:end] .= 0;
img = mgrid

# genrate flow
flow = gen_tiled_flow(size(img), 50);

# generate warpped image
imgw = LAP_julia.interpolation.imWarp(img, real(flow), imag(flow))

## using mandrill image

# genrate image MANDRILL
img = testimage("mandril_gray")
imgr = reverse(img, dims=1)
img = imgr

# genrate flow
flow = gen_tiled_flow(size(img), 30);
showflow(flow)

# generate warpped image
imgw = LAP_julia.interpolation.imWarp(img, real(flow), imag(flow));

# show imgw
imgshow(Float32.(imgw))

imgshow(img)

## using chessboard image

tile_size = 5
board_size = 4 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]
chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
img = chessboard

img = gen_chess()

# genrate flow
flow = gen_tiled_flow(size(img), 20, 1000);

# generate warpped image
imgw = LAP_julia.interpolation.warp_img(img, real(flow), imag(flow))

# show imgw
imgshow(imgw)

# show flow
showflow(flow)

maximum(filter(!isnan, real(flow)))
minimum(filter(!isnan, real(flow)))

### RUN OPTIFLOW:

u_est, source_reg, figs = pflap(img, imgw);

# compare with a single lap:
u_sin_est = single_lap(img, imgw, 32, [65, 65], 3)
showflow(u_sin_est)

LAP_julia.inpaint.inpaint_nans!(u_sin_est)
single_img = LAP_julia.interpolation.imWarp_replicate(imgw, real(u_sin_est), imag(u_sin_est))
imgshow(single_img)
imgshow(imgw)

### compare results

# estimation
maximum(abs.(imag(u_est)))
showflow(u_est)
showflow(u_est .* (-1) .- flow)

# rewarp check
dewarp = LAP_julia.interpolation.imWarp(imgw, real(u_est), imag(u_est))
imgshow(Float32.(dewarp))
imgshow(Float32.(img))

imgshow(source_reg)

# truth
maximum(abs.(imag(flow)))
sum(abs.(real(flow)))/length(real(flow))
showflow(flow, mag = 1)

# show source_reg
fig = PyPlot.figure(dpi = 300, figsize = (5, 5))
PyPlot.imshow(source_reg, cmap = :gray);
gcf()

# show imgw
fig = PyPlot.figure(dpi = 300, figsize = (5, 5))
PyPlot.imshow(Float64.(imgw), cmap = :gray);
gcf()

# k = (-1:1) * ones(1,3)
# A =      [1     2     3;
#           4     5     6;
#           7     8     9];
#
# transpose(A[:]) * k[:]
#
# sum(transpose(A) * (-1:1))
#
# # l0
#
# l = transpose((-1:1) * ones(1,3))
#
# transpose(A[:]) * l[:]
#
# sum(A * (-1:1))
1


text:
# # Example (Ch. 9 Bathe)

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/example_9_1_Bathe.ipynb)
#

# The example considered next can be found in Chapter 9 of [[BATHE]](@ref).

using StructuralDynamicsODESolvers, Plots

# ## Problem formulation
M = [2 0; 0 1.]
K = [6 -2; -2 4.]
C = zeros(2, 2)
f = [0.0, 10.0]
example_9_1_Bathe = SecondOrderAffineContinuousSystem(M, C, K, f)

NSTEPS = 500
tdom = range(0, NSTEPS * 0.1, length=NSTEPS + 1)
U₀, U₀′ = zeros(2), zeros(2)
prob = InitialValueProblem(example_9_1_Bathe, (U₀, U₀′))

# ## Analytic solution
A = [1/√3  (1/2)*√(2/3);
     1/√3      -√(2/3)]
x₁(t) = (5 / √3) * (1 - cos(t*√2))
x₂(t) = (2 * √(2/3)) * (-1 + cos(t*√5))
U(t) = A * [x₁(t), x₂(t)]

# ## Central difference

sol = solve(prob, CentralDifference(Δt=0.1); NSTEPS=NSTEPS) |> displacements
ind = 150:170
fig = plot(xlab="time", ylab="x1(t)", legend=:outertopright)
fig2 = plot(xlab="time", ylab="x1(t)", legend=:outertopright)
plot!(fig, tdom, [s[1] for s in sol], lab="Central difference")
plot!(fig2, tdom[ind], [s[1] for s in sol[ind]], lab="Central difference")

# ## Houbolt

sol = solve(prob, Houbolt(Δt=0.1); NSTEPS=NSTEPS) |> displacements
plot!(fig, tdom, [s[1] for s in sol], lab="Houbolt")
plot!(fig2, tdom[ind], [s[1] for s in sol[ind]], lab="Houbolt")

# ## Newmark

sol = solve(prob, Trapezoidal(Δt=0.1); NSTEPS=NSTEPS) |> displacements
plot!(fig, tdom, [s[1] for s in sol], lab="Newmark")
plot!(fig2, tdom[ind], [s[1] for s in sol[ind]], lab="Newmark")

# ## Bathe

sol = solve(prob, Bathe(Δt=0.1); NSTEPS=NSTEPS) |> displacements
plot!(fig, tdom, [s[1] for s in sol], lab="Bathe")
plot!(fig2, tdom[ind], [s[1] for s in sol[ind]], lab="Bathe")

# ## Analytic solution

tdom = range(0, NSTEPS*0.1, length=1000)
plot!(fig, tdom, [U(t)[1] for t in tdom], lab="Analytic")
plot!(fig2, tdom[299:330], [U(t)[1] for t in tdom[299:330]], lab="Analytic")


text:
const Tup = Union{Tuple, NamedTuple}
const EmptyTup = Union{Tuple{}, NamedTuple{(), Tuple{}}}

@generated function staticschema(::Type{T}) where {T}
    name_tuple = Expr(:tuple, [QuoteNode(f) for f in fieldnames(T)]...)
    type_tuple = Expr(:curly, :Tuple, [Expr(:call, :fieldtype, :T, i) for i in 1:fieldcount(T)]...)
    Expr(:curly, :NamedTuple, name_tuple, type_tuple)
end

staticschema(::Type{T}) where {T<:Tup} = T

function createinstance(::Type{T}, args...) where {T}
    isconcretetype(T) ? bypass_constructor(T, args) : T(args...)
end

createinstance(::Type{T}, args...) where {T<:Tup} = T(args)


text:
# Copyright (c) 2019 Guido Kraemer
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

"""
    TxOut

Each output spends a certain number of satoshis, placing them under control of
anyone who can satisfy the provided pubkey script.
A `TxOut` is composed of
- `value::UInt64`, number of satoshis to spend. May be zero; the sum of all
outputs may not exceed the sum of satoshis previously spent to the outpoints
provided in the input section. (Exception: coinbase transactions spend the
block subsidy and collected transaction fees.)
- `scriptpubkey::Script` which defines the conditions which must be
satisfied to spend this output.

"""
struct TxOut
    value           :: UInt64
    scriptpubkey    :: Script
end

"""
    TxOut(io::IOBuffer)

Parse an `IOBuffer` to a `TxOut`
"""
function TxOut(io::IOBuffer)
    value = read(io, UInt64)
    scriptpubkey = Script(io)

    TxOut(value, scriptpubkey)
end

function Base.show(io::IO, output::TxOut)
    println(io, "Transaction output: ₿" * string(output.value/100000000))
    println(io, "ScriptPubKey:")
    println(io, output.scriptpubkey)
end

"""
    serialize(tx::TxOut) -> Vector{UInt8}

Returns the byte serialization of the transaction output
"""
function serialize(tx::TxOut)
    result = bytes(tx.value, len=8, little_endian=true)
    append!(result, serialize(tx.scriptpubkey))
    return result
end


text:
import ControlSystems: feedback

"""
    FRD(w,r)

Represents frequency-response data. `w` holds the frequency vector and `r` the response. Methods defined on this type include
- `+-*`
- `length, vec, sqrt`
- `plot`
- `feedback`
"""
struct FRD{WT<:AbstractVector,RT<:AbstractVector} <: LTISystem
    w::WT
    r::RT
end

FRD(w, s::LTISystem) = FRD(w, freqresp(s,w)[:,1,1])
import Base: +, -, *, length, sqrt
Base.vec(f::FRD) = f.r
*(f::FRD, f2)    = FRD(f.w, f.r .* vec(f2))
+(f::FRD, f2)    = FRD(f.w, f.r .+ vec(f2))
-(f::FRD, f2)    = FRD(f.w, f.r .- vec(f2))
-(f::FRD)        = FRD(f.w, -f.r)
length(f::FRD)   = length(f.w)
sqrt(f::FRD)     = FRD(f.w, sqrt.(f.r))

feedback(P::FRD,K) = FRD(P.w,vec(P)./(1. .+ vec(P).*vec(K)))
feedback(P,K::FRD) = FRD(K.w,vec(P)./(1. .+ vec(P).*vec(K)))
feedback(P::FRD,K::FRD) = feedback(P,vec(K))

feedback(P::FRD,K::LTISystem) = feedback(P,freqresp(K,P.w)[:,1,1])
feedback(P::LTISystem,K::FRD) = feedback(freqresp(P,K.w)[:,1,1],K)

@recipe function plot_frd(frd::FRD)
    yscale --> :log10
    xscale --> :log10
    xlabel --> "Frequency [rad/s]"
    ylabel --> "Magnitude"
    title --> "Bode Plot"
    legend --> false
    @series begin
        inds = findall(x->x==0, frd.w)
        useinds = setdiff(1:length(frd.w), inds)
        frd.w[useinds], abs.(frd.r[useinds])
    end
    nothing
end

freqvec(h,k) = LinRange(0,π/h, length(k))

"""
    H, N = tfest(h,y,u, σ = 0.05)

Estimate a transfer function model using the Correlogram approach.
    Both `H` and `N` are of type `FRD` (frequency-response data).

`σ` determines the width of the Gaussian window applied to the estimated correlation functions before FFT.
- `H` = Syu/Suu             Process transfer function
- `N` = Sy - |Syu|²/Suu     Noise PSD
"""
function tfest(h::Real,y::AbstractVector,u::AbstractVector, σ = 0.05)
    Syy,Suu,Syu = fft_corr(y,u,σ)
    w   = freqvec(h,Syu)
    H   = FRD(w, Syu./Suu)
    N   = FRD(w, @.(Syy - abs2(Syu)/Suu)./length(y))
    return H, N
end

function fft_corr(y,u,σ = 0.05)
    n = length(y)
    w = gaussian(2n-1,σ)

    Syu = rfft(ifftshift(w.*xcorr(y,u)))
    Syy = rfft(ifftshift(w.*xcorr(y,y)))
    Suu = rfft(ifftshift(w.*xcorr(u,u)))
    Syy,Suu,Syu
end

"""
    κ, N = coherence(h,y,u)

Calculates the coherence Function. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy. An estimated noise model is also returned.
κ: Coherence function (not squared)
N: Noise model
"""
function coherence(h::Real,y::AbstractVector,u::AbstractVector; n = length(y)÷10, noverlap = n÷2, window=hamming)
    Syy,Suu,Syu = wcfft(y,u,n=n,noverlap=noverlap,window=window)
    k = (abs2.(Syu)./(Suu.*Syy))#[end÷2+1:end]
    Sch = FRD(freqvec(h,k),k)
    return Sch
end

function wcfft(y,u; n = length(y)÷10, noverlap = n÷2, window=hamming)
    win, norm2 = DSP.Periodograms.compute_window(window, n)
    uw  = arraysplit(u,n,noverlap,nextfastfft(n),win)
    yw  = arraysplit(y,n,noverlap,nextfastfft(n),win)
    Syy = zeros(length(uw[1])÷2 + 1)
    Suu = zeros(length(uw[1])÷2 + 1)
    Syu = zeros(ComplexF64,length(uw[1])÷2 + 1)
    for i in eachindex(uw)
        xy      = rfft(yw[i])
        xu      = rfft(uw[i])
        Syu .+= xy.*conj.(xu)
        Syy .+= abs2.(xy)
        Suu .+= abs2.(xu)
    end
    Syy,Suu,Syu
end

function wfft_corr(y,u; n = length(y)÷10, noverlap = n÷2, window=hamming)
    win, norm2 = DSP.Periodograms.compute_window(window, n)
    Cyu = xcorr(y,u)
    Cyy = xcorr(y,y)
    Cuu = xcorr(u,u)
    uw  = arraysplit(Cuu,n,noverlap,nextfastfft(n),win)
    yw  = arraysplit(Cyy,n,noverlap,nextfastfft(n),win)
    yuw = arraysplit(Cyu,n,noverlap,nextfastfft(n),win)
    Syy = zeros(ComplexF64,length(uw[1])÷2 + 1)
    Suu = zeros(ComplexF64,length(uw[1])÷2 + 1)
    Syu = zeros(ComplexF64,length(uw[1])÷2 + 1)
    for i in eachindex(uw)
        xy    = rfft(yw[i])
        xu    = rfft(uw[i])
        xyu   = rfft(yuw[i])
        Syu .+= xyu
        Syy .+= xy
        Suu .+= xu
    end
    Syy,Suu,Syu
end

@userplot Coherenceplot
"""
coherenceplot(h,y,u)

Calculates and plots the coherence Function κ. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy.
"""
coherenceplot

@recipe function cp(p::Coherenceplot)
    ae = ArgumentError("Call like this: coherenceplot(h,y,u), where h is sample time and y/u are vectors of equal length.")
    length(p.args) == 3 || throw(ae)
    h,y,u = p.args[1:3]
    length(u) == length(y) || throw(ae)
    h isa Number || throw(ae)
    yscale --> :identity
    xscale --> :log10
    ylims --> (0,1)
    xlabel --> "Frequency [rad/s]"
    title --> "Coherence"
    legend --> false
    frd = coherence(h,y,u)
    @series begin
        inds = findall(x->x==0, frd.w)
        useinds = setdiff(1:length(frd.w), inds)
        frd.w[useinds], abs.(frd.r[useinds])
    end
    nothing
end

"""
ir, t, Σ = impulseest(h,y,u,n)

Estimates the system impulse response by fitting an `n`:th order FIR model. Returns impulse-response estimate, time vector and covariance matrix.
See also `impulseestplot`
"""
function impulseest(h,y,u,n,λ=0)
    yt,A = getARXregressor(y,u,0,n)
    ir = ls(A,yt,λ)
    t = range(h,length=n, step=h)
    Σ = parameter_covariance(yt, A, ir, λ)
    ir./h, t, Σ
end

@userplot Impulseestplot

"""
impulseestplot(h,y,u,n)

Estimates the system impulse response by fitting an `n`:th order FIR model and plots the result with a 95% confidence band.
See also `impulseestplot`
"""
impulseestplot
@recipe function impulseestplot(p::Impulseestplot)
    h,y,u = p.args[1:3]
    n = length(p.args) >= 4 ? p.args[4] : 25
    λ = length(p.args) >= 5 ? p.args[5] : 0
    ir,t,Σ = impulseest(h,y,u,n,λ)
    title --> "Estimated Impulse Response"
    xlabel --> "Time [s]"

    @series begin
        label --> ""
        t,ir
    end
    linestyle := :dash
    color := :black
    label := ""
    seriestype := :hline
    @series begin
        t, 2 .*sqrt.(diag(Σ))
    end
    @series begin
        t, -2 .*sqrt.(diag(Σ))
    end
end


text:
# ENV["JULIA_DEBUG"] = "all"

Base.push!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
import VSCodeDebugger
pop!(LOAD_PATH)

let
    local pipenames = VSCodeDebugger.DebugAdapter.clean_up_ARGS_in_launch_mode()
    try
        VSCodeDebugger.DebugAdapter.startdebug(pipenames[1])
    catch err
        VSCodeDebugger.DebugAdapter.global_err_handler(err, catch_backtrace(), pipenames[2], "Debugger")
    end
end

max_stars_repo_name: fonsp/julia-vscode

max_issues_repo_name: fonsp/julia-vscode

max_forks_repo_name: fonsp/julia-vscode


text:
module MakieLayout

using AbstractPlotting
using AbstractPlotting: Rect2D
import AbstractPlotting: IRect2D
using AbstractPlotting.Keyboard
using AbstractPlotting.Mouse
using AbstractPlotting: ispressed, is_mouseinside
using Observables: onany
import Observables
import Formatting
using Match
import Animations
import PlotUtils
using GridLayoutBase
import Showoff

const FPS = Node(30)
const COLOR_ACCENT = Ref(RGBf0(((79, 122, 214) ./ 255)...))
const COLOR_ACCENT_DIMMED = Ref(RGBf0(((174, 192, 230) ./ 255)...))

# Make GridLayoutBase default row and colgaps themeable when using MakieLayout
# This mutates module-level state so it could mess up other libraries using
# GridLayoutBase at the same time as MakieLayout, which is unlikely, though
function __init__()

    @warn """
    MakieLayout.jl has been absorbed by AbstractPlotting.jl and will not receive any more updates.
    It is enough to install AbstractPlotting to receive the newest functionality.
    Visit the Makie docs at http://makie.juliaplots.org/ for more information on how to use it.
    """

    GridLayoutBase.DEFAULT_COLGAP_GETTER[] = function()
        ct = AbstractPlotting.current_default_theme()
        if haskey(ct, :colgap)
            ct[:colgap][]
        else
            GridLayoutBase.DEFAULT_COLGAP[]
        end
    end
    GridLayoutBase.DEFAULT_ROWGAP_GETTER[] = function()
        ct = AbstractPlotting.current_default_theme()
        if haskey(ct, :rowgap)
            ct[:rowgap][]
        else
            GridLayoutBase.DEFAULT_ROWGAP[]
        end
    end
end

include("types.jl")
include("helpers.jl")
include("mousestatemachine.jl")
include("ticklocators/linear.jl")
include("ticklocators/wilkinson.jl")
include("defaultattributes.jl")
include("lineaxis.jl")
include("lobjects/laxis.jl")
include("lobjects/lcolorbar.jl")
include("lobjects/ltext.jl")
include("lobjects/lslider.jl")
include("lobjects/lbutton.jl")
include("lobjects/lrect.jl")
include("lobjects/ltoggle.jl")
include("lobjects/llegend.jl")
include("lobjects/lobject.jl")
include("lobjects/lscene.jl")
include("lobjects/lmenu.jl")

export LAxis
export LSlider
export LButton
export LColorbar
export LText
export LRect
export LToggle
export LLegend
export LegendEntry, MarkerElement, PolyElement, LineElement, LegendElement
export LScene
export LMenu
export linkxaxes!, linkyaxes!, linkaxes!
export AxisAspect, DataAspect
export autolimits!, limits!
export LinearTicks, WilkinsonTicks
export hidexdecorations!, hideydecorations!, hidedecorations!, hidespines!
export tight_xticklabel_spacing!, tight_yticklabel_spacing!, tight_ticklabel_spacing!, tightlimits!
export layoutscene
export set_close_to!
export xaxis_bottom!, xaxis_top!, yaxis_left!, yaxis_right!
export labelslider!
export addmousestate!

# from GridLayoutBase
export GridLayout, GridPosition
export GridLayoutSpec
export BBox
export LayoutObservables
export Inside, Outside, Mixed
export Fixed, Auto, Relative, Aspect
export width, height, top, bottom, left, right
export with_updates_suspended
export appendcols!, appendrows!, prependcols!, prependrows!, deletecol!, deleterow!, trim!
export gridnest!
export AxisAspect, DataAspect
export colsize!, rowsize!, colgap!, rowgap!
export Left, Right, Top, Bottom, TopLeft, BottomLeft, TopRight, BottomRight
export grid!, hbox!, vbox!
export swap!
export ncols, nrows

end # module

max_stars_repo_name: jkrumbiegel/MakieLayout


text:
testOnAllContexts("DNA") do ctx
    println("\tconstruct:")

    println("\tcreate from a number")

    dna = DNA(ctx, 42)
    @test length(dna) == 42
    values = Array(dna)
    @test typeof(values) == Vector{Float32}
    @test length(values) == 42
    for value in values
        @test value == 0f0
    end

    println("\tcreate from an context's array")

    arr = constant(ctx, 1f0, 42)
    dna = DNA(ctx, arr)
    @test length(dna) == 42
    values = Array(dna)
    @test typeof(values) == Vector{Float32}
    @test length(values) == 42
    for value in values
        @test value == 1f0
    end

    println("\tcrate from another dna")

    dna2 = DNA(dna)
    @test length(dna2) == 42
    values = Array(dna2)
    @test typeof(values) == Vector{Float32}
    @test length(values) == 42
    for value in values
        @test value == 1f0
    end

    println("\tcrate from Vector{Float32}")

    values = [1f0, 2f0, 3f0]
    dna = DNA(ctx, values)
    @test length(dna) == 3
    values2 = Array(dna)
    @test values == values2

    println("\trandomize uniformly")

    dna = DNA(ctx, 100)
    randomizeUniform!(dna, 100)
    @test length(dna) == 100

    values = Array(dna)
    sum = 0f0
    lo = false
    hi = false
    for value in values
        @test value >= 0f0 && value < 100f0
        sum += value
        if value < 50f0
            lo = true
        else
            hi = true
        end
    end
    @test sum > 0f0
    @test lo && hi

    println("\tcrossover")

    dna1 = DNA(ctx, 100)
    randomizeUniform!(dna1, 100)
    dna2 = DNA(ctx, 100)
    randomizeUniform!(dna2, 100)
    child = crossover(dna1, dna2)
    @test length(child) == 100

    values1 = Array(dna1)
    values2 = Array(dna2)
    values = Array(child)

    found1 = 0
    found2 = 0
    invalid = 0
    for i in 1:100
        if values[i] == values1[i]
            found1 += 1
        elseif values[i] == values2[i]
            found2 += 1
        else
            invalid += 1
        end
    end
    @test found1 > 10
    @test found2 > 10
    @test invalid == 0

    println("\tmutate")

    dna = DNA(ctx, 1000) # filled with zero
    mutate!(dna, 0.05f0, 0.1f0) # 5% chance of 0.1 = +/-0.05
    values = Array(dna)
    count = 0
    for value in values
        @test value >= -0.05f0 && value <= 0.5f0
        if value != 0f0
            count += 1
        end
    end
    @test count > 20 && count < 80

    println("\tnormalized")

    function verify(a)
        @test a[1] == 0.5f0
        @test a[2] == 0f0
        @test a[3] == 0.5f0
        @test a[4] == 0.75f0
        @test a[5] == 1f0
        @test a[6] == 0.5f0
        @test a[7] == 0.25f0
        @test a[8] == 0f0
        @test a[9] == 0.25f0
        @test a[10] == 0.5f0
        @test a[11] == 1f0
        @test a[12] == 0.75f0
        @test a[13] == 0.5f0
        @test a[14] == 0f0
        @test a[15] == 0.5f0
    end

    dna = DNA(ctx, [
        -2.5f0,         # 0.5
        -2f0,           # 0
        -1.5f0,         # 0.5
        -1.25f0,        # 0.75
        -1f0,           # 1
        -0.5f0,         # 0.5
        -0.25f0,        # 0.25
        0f0,            # 0
        0.25f0,         # 0.25
        0.5f0,          # 0.5
        1f0,            # 1
        1.25f0,         # 0.75
        1.5f0,          # 0.5
        2f0,            # 0
        2.5f0           # 0.5
    ])

    dna2 = normalized(dna)
    @test length(dna2) == 15
    values = Array(dna2)
    @test length(values) == 15
    verify(values)

    println("\tnormalize!")

    values = Array(dna)
    @test_throws ErrorException verify(values)
    normalize!(dna)
    values = Array(dna)
    verify(values)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Ramnet.Encoders
using Ramnet

using Random

using ..Buffers: MultiStepDynamicBuffer, add!

mutable struct MonteCarloDiscountedDiscriminatorAgent{O <: AbstractVector,A <: Real,E <: AbstractEncoder} <: AbstractAgent{O,A}
    actions::StepRange{A,Int}
    n::Int
    size::Int
    regressor_discount::Float64
    γ::Float64
    ϵ::Float64
    # episodes::Int
    encoder::E
    Q̂::Dict{A,FastRegressionDiscriminator}
    buffer::MultiStepDynamicBuffer{O,A}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,A}
    rng::MersenneTwister

    # TODO: Validate actions and encoder and ϵ
    # TODO: Do I need the size if the encoder is given? Not if a knew the observation length (static array?)
    function MonteCarloDiscountedDiscriminatorAgent{O,A,E}(actions, n, obs_size, regressor_discount, γ, ϵ, encoder::E; seed::Union{Nothing,Int}=nothing) where {O,A <: Real,E <: AbstractEncoder}
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
        rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

        n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))
        # size < n && throw(DomainError(size, "Input size may not be smaller then tuple size"))
        size = obs_size * resolution(encoder)
        # !(0 ≤ discount ≤ 1) && throw(DomainError(discount, "Discount must lie in the [0, 1] interval"))
        # episodes < 1 && throw(DomainError(episodes, "Number of episodes must be greater then one"))

        Q̂ = Dict(a => FastRegressionDiscriminator(size, n; γ=regressor_discount) for a in actions)
        buffer = MultiStepDynamicBuffer{O,A}()

        new(actions, n, size, regressor_discount, γ, ϵ, encoder, Q̂, buffer, nothing, false, nothing, rng)
    end
end

function MonteCarloDiscountedDiscriminatorAgent{O,A,E}(
    actions, obs_size, encoder::E; n, regressor_discount, gamma, epsilon, seed) where {O,A <: Real,E <: AbstractEncoder}
    MonteCarloDiscountedDiscriminatorAgent{O,A,E}(actions, n, obs_size, regressor_discount, gamma, epsilon, encoder; seed)
end

function observe!(agent::MonteCarloDiscountedDiscriminatorAgent{O,A,E}, observation::O) where {O,A <: Real,E <: AbstractEncoder}
    agent.observation = observation
    agent.done = false

    nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::MonteCarloDiscountedDiscriminatorAgent{O,A,E}, action::A, reward::Float64, observation::O, done::Bool) where {O,A <: Real,E <: AbstractEncoder}
    add!(agent.buffer, agent.observation, agent.action, reward)

    agent.observation = observation
    agent.done = done

    nothing
end

function update!(agent::MonteCarloDiscountedDiscriminatorAgent)
    if agent.done
        G = 0.0
        for transition in Iterators.reverse(agent.buffer)
            G = transition.reward + agent.γ * G
            train!(
                agent.Q̂[transition.action],
                encode(agent.encoder, transition.observation),
                G
            )
        end

        reset!(agent.buffer)
    end

    nothing
end

function Agents.reset!(agent::MonteCarloDiscountedDiscriminatorAgent)
    for action in agent.actions
        agent.Q̂[action] = FastRegressionDiscriminator(agent.size, agent.n; γ=agent.regressor_discount)
    end

    reset!(agent.buffer)

    nothing
end


text:
export SolverInfo, rownorm², nrmsd

"""
The Solver Info type is used to obtain aditional information
of linear solvers and their iteration process.

# Fields
* `convMeas::Vector{Float64}` - Vector with convergence parameters
* `nrmse::Vector{Float64}`    - NRMSD of the current iterate w.r.t some reference (if provided)
* `x_ref::Vector{T}`          - Reference for computing NRMSD of iterates
* `x_iter::Vector{Vector{T}}` - Vector for storing iterates
* `store_solutions::Bool`     - specification whether iterates are to be stored

the content of `convMeas` is specific to each solver. It can be found in the documentation of the
corresping `solve`-method.
"""
mutable struct SolverInfo{T}
  convMeas::Vector{Float64}
  nrmse::Vector{Float64}
  x_ref::Vector{T}
  x_iter::Vector{Vector{T}}
  store_solutions::Bool
end

function SolverInfo(x_ref::Vector{T}; store_solutions::Bool=false,
    kargs...) where T
  return SolverInfo(Vector{Float64}(), Vector{Float64}(), x_ref, Vector{Vector{T}}(), store_solutions)
end

function SolverInfo(T::Type=Number; store_solutions::Bool=false,
    kargs...)
  return SolverInfo(Vector{Float64}(), Vector{Float64}(), Vector{T}(), Vector{Vector{T}}(), store_solutions)
end

function storeInfo(solverInfo::SolverInfo,x::Vector{T},convMeas::Vararg{rT,N}) where {T,rT<:Real,N}
  # convergence criteria
  push!(solverInfo.convMeas,convMeas...)
  # solution
  if solverInfo.store_solutions
    push!(solverInfo.x_iter, deepcopy(x))
  end
  # nrmse
  if !isempty(solverInfo.x_ref)
    push!( solverInfo.nrmse, nrmsd(solverInfo.x_ref,x) )
  end
end

### rownorm² ###

"""
This function computes the 2-norm² of a rows of S for dense matrices.
"""
function rownorm²(B::Transpose{T,S},row::Int) where {T,S<:DenseMatrix}
  A = B.parent
  U = typeof(real(A[1]))
  res::U = BLAS.nrm2(size(A,1), pointer(A,(LinearIndices(size(A)))[1,row]), 1)^2
  return res
end

function rownorm²(A::AbstractMatrix,row::Int)
  T = typeof(real(A[1]))
  res = zero(T)
  @simd for n=1:size(A,2)
    res += abs2(A[row,n])
  end
  return res
end

"""
This function computes the 2-norm² of a rows of S for dense matrices.
"""
function rownorm²(B::Transpose{T,S},row::Int) where {T,S<:SparseMatrixCSC}
  A = B.parent
  U = typeof(real(A[1]))
  res::U = BLAS.nrm2(A.colptr[row+1]-A.colptr[row], pointer(A.nzval,A.colptr[row]), 1)^2
  return res
end

### dot_with_matrix_row ###

# Fallback implementation
#=function dot_with_matrix_row_simd{T<:Complex}(A::AbstractMatrix{T}, x::Vector{T}, k::Int64)
  res = zero(T)
  @simd for n=1:size(A,2)
    @inbounds res += conj(A[k,n])*x[n]
  end
  return res
end=#

"""
This funtion calculates ∑ᵢ Aᵢₖxᵢ for dense matrices.
"""
function dot_with_matrix_row(A::DenseMatrix{T}, x::Vector{T}, k::Int64) where {T<:Complex}
  BLAS.dotu(length(x), pointer(A, (LinearIndices(size(A)))[k,1]), size(A,1), pointer(x,1), 1)
end

function dot_with_matrix_row(B::Transpose{T,S},
                       x::Vector{T}, k::Int64) where {T<:Complex,S<:DenseMatrix}
  A = B.parent
  BLAS.dotu(length(x), pointer(A,(LinearIndices(size(A)))[1,k]), 1, pointer(x,1), 1)
end

"""
This funtion calculates ∑ᵢ Aᵢₖxᵢ for dense matrices.
"""
function dot_with_matrix_row(A::DenseMatrix{T}, x::Vector{T}, k::Int64) where {T<:Real}
  BLAS.dot(length(x), pointer(A,(LinearIndices(size(A)))[k,1]), size(A,1), pointer(x,1), 1)
end

"""
This funtion calculates ∑ᵢ Aᵢₖxᵢ for dense matrices.
"""
function dot_with_matrix_row(B::Transpose{T,S},
      x::Vector{T}, k::Int64) where {T<:Real,S<:DenseMatrix}
  A = B.parent
  BLAS.dot(length(x), pointer(A,(LinearIndices(size(A)))[1,k]), 1, pointer(x,1), 1)
end

"""
This funtion calculates ∑ᵢ Aᵢₖxᵢ for sparse matrices.
"""
function dot_with_matrix_row(B::Transpose{T,S},
                             x::Vector{T}, k::Int64) where {T,S<:SparseMatrixCSC}
  A = B.parent
  tmp = zero(T)
  N = A.colptr[k+1]-A.colptr[k]
  for n=A.colptr[k]:N-1+A.colptr[k]
    @inbounds tmp += A.nzval[n]*x[A.rowval[n]]
  end
  tmp
end

### enfReal! / enfPos! ###

"""
This function enforces the constraint of a real solution.
"""
function enfReal!(x::Vector{T}, mask=ones(Bool, length(x))) where {T<:Complex}
  #Returns x as complex vector with imaginary part set to zero
  @simd for i in 1:length(x)
    @inbounds mask[i] && (x[i] = complex(x[i].re))
  end
end

"""
This function enforces the constraint of a real solution.
"""
enfReal!(x::Vector{T}, mask=ones(Bool, length(x))) where {T<:Real} = nothing

"""
This function enforces positivity constraints on its input.
"""
function enfPos!(x::Vector{T}, mask=ones(Bool, length(x))) where {T<:Complex}
  #Return x as complex vector with negative parts projected onto 0
  @simd for i in 1:length(x)
    @inbounds (x[i].re < 0 && mask[i]) && (x[i] = im*x[i].im)
  end
end

"""
This function enforces positivity constraints on its input.
"""
function enfPos!(x::Vector{T}, mask=ones(Bool, length(x))) where {T<:Real}
  #Return x as complex vector with negative parts projected onto 0
  @simd for i in 1:length(x)
    @inbounds (x[i] < 0 && mask[i]) && (x[i] = zero(T))
  end
end

function applyConstraints(x, sparseTrafo, enforceReal, enforcePositive, constraintMask=ones(Bool, length(x)) )

  mask = (constraintMask != nothing) ? constraintMask : ones(Bool, length(x))

  if sparseTrafo != nothing
     x[:] = sparseTrafo * x
  end
  enforceReal && enfReal!(x, mask)
  enforcePositive && enfPos!(x, mask)
  if sparseTrafo != nothing
    x[:] = adjoint(sparseTrafo)*x
  end
end

### im2col / col2im ###

"""
This function rearranges distinct image blocks into columns of a matrix.
"""
function im2colDistinct(A::Array{T}, blocksize::NTuple{2,Int64}) where T

  nrows = blocksize[1]
  ncols = blocksize[2]
  nelem = nrows*ncols

  # padding for A such that patches can be formed
  row_ext = mod(size(A,1),nrows)
  col_ext = mod(size(A,2),nrows)
  pad_row = (row_ext != 0)*(nrows-row_ext)
  pad_col = (col_ext != 0)*(ncols-col_ext)

  # rearrange matrix
  A1 = zeros(T, size(A,1)+pad_row, size(A,2)+pad_col)
  A1[1:size(A,1),1:size(A,2)] = A

  t1 = reshape( A1,nrows, floor(Int,size(A1,1)/nrows), size(A,2) )
  t2 = reshape( permutedims(t1,[1 3 2]), size(t1,1)*size(t1,3), size(t1,2) )
  t3 = permutedims( reshape( t2, nelem, floor(Int,size(t2,1)/nelem), size(t2,2) ),[1 3 2] )
  res = reshape(t3,nelem,size(t3,2)*size(t3,3))

  return res
end

"""
This funtion rearrange columns of a matrix into blocks of an image.
"""
function col2imDistinct(A::Array{T}, blocksize::NTuple{2,Int64},
                 matsize::NTuple{2,Int64}) where T
  # size(A) should not be larger then (blocksize[1]*blocksize[2], matsize[1]*matsize[2]).
  # otherwise the bottom (right) lines (columns) will be cut.
  # matsize should be divisble by blocksize.

  if mod(matsize[1],blocksize[1]) != 0 || mod(matsize[2],blocksize[2]) != 0
    error("matsize should be divisible by blocksize")
  end

  blockrows = blocksize[1]
  blockcols = blocksize[2]
  matrows = matsize[1]
  matcols = matsize[2]
  mblock = floor(Int,matrows/blockrows) # number of blocks per row
  nblock = floor(Int,matcols/blockcols) # number of blocks per column

  # padding for A such that patches can be formed and arranged into a matrix of
  # adequate size
  row_ext = blockrows*blockcols-size(A,1)
  col_ext = mblock*nblock-size(A,2)
  pad_row = (row_ext > 0 )*row_ext
  pad_col = (col_ext > 0 )*col_ext

  A1 = zeros(T, size(A,1)+pad_row, size(A,2)+pad_col)
  A1[1:blockrows*blockcols, 1:mblock*nblock] = A[1:blockrows*blockcols, 1:mblock*nblock]

  # rearrange matrix
  t1 = reshape( A1, blockrows,blockcols,mblock*nblock )
  t2 = reshape( permutedims(t1,[1 3 2]), matrows,nblock,blockcols )
  res = reshape( permutedims(t2,[1 3 2]), matrows,matcols)

end

### NRMS ###

function nrmsd(I,Ireco)
  N = length(I)

  # This is a little trick. We usually are not interested in simple scalings
  # and therefore "calibrate" them away
  alpha = norm(Ireco)>0 ? (dot(vec(I),vec(Ireco))+dot(vec(Ireco),vec(I))) /
          (2*dot(vec(Ireco),vec(Ireco))) : 1.0
  I2 = Ireco.*alpha

  RMS =  1.0/sqrt(N)*norm(vec(I)-vec(I2))
  NRMS = RMS/(maximum(abs.(I))-minimum(abs.(I)) )
  return NRMS
end


text:
@with_kw struct iLQSolver{TLM, TOM, TQM}
    "The regularization term for the state cost quadraticization."
    state_regularization::Float64 = 0.0
    "The regularization term for the control cost quadraticization."
    control_regularization::Float64 = 0.0
    "The initial scaling of the feed-forward term."
    α_scale_init::Float64 = 0.5
    "The geometric scaling of the feed-forward term per scaling step in
    backtrack scaling."
    α_scale_step::Float64 = 0.5
    "Iteration is aborted if this number is exceeded."
    max_n_iter::Int = 200
    "The maximum number of backtrackings per scaling step"
    max_scale_backtrack::Int = 20
    "The maximum elementwise difference bewteen operating points for
    convergence."
    max_elwise_diff_converged::Float64 = α_scale_init/10
    "The maximum elementwise difference bewteen operating points for per
    iteration step."
    max_elwise_diff_step::Float64 = 30 * max_elwise_diff_converged
    "Preallocated memory for lq approximations."
    _lq_mem::TLM
    "Preallocated memory for quadraticization results."
    _qcache_mem::TQM
    "Preallocated memory for operting points."
    _op_mem::TOM
end

qcache(solver::iLQSolver) = solver._qcache_mem

function regularize(solver::iLQSolver, c::QuadraticPlayerCost)
    return QuadraticPlayerCost(c.l, c.Q + I * solver.state_regularization,
                               c.r, c.R + I * solver.control_regularization)
end

function iLQSolver(g, args...; kwargs...)
    return iLQSolver(args...; kwargs...,
                     _lq_mem=lqgame_preprocess_alloc(g),
                     _qcache_mem=zero(QuadCache{n_states(g),n_controls(g)}),
                     _op_mem=zero(SystemTrajectory, g))
end

function has_converged(solver::iLQSolver, last_op::SystemTrajectory{h},
                       current_op::SystemTrajectory{h}) where {h}
    return are_close(current_op, last_op, solver.max_elwise_diff_converged)
end

function are_close(op1::SystemTrajectory, op2::SystemTrajectory,
                   max_elwise_diff::Float64)
    @assert horizon(op1) == horizon(op2)
    return all(infnorm(op1.x[k] - op2.x[k]) < max_elwise_diff for k in
               eachindex(op2.x))
end

# modifies the current strategy to stabilize the update
function scale!(current_strategy::SizedVector, current_op::SystemTrajectory,
                α_scale::Float64)
    map!(current_strategy, current_strategy) do el
        return AffineStrategy(el.P, el.α * α_scale)
    end
end

function backtrack_scale!(current_strategy::SizedVector,
                          current_op::SystemTrajectory, last_op::SystemTrajectory,
                          cs::ControlSystem, solver::iLQSolver)
    for i in 1:solver.max_scale_backtrack
        # initially we do a large scaling. Afterwards, always half feed forward
        # term.
		sf = i == 1 ? solver.α_scale_init : solver.α_scale_step
        scale!(current_strategy, current_op, sf)
        # we compute the new trajectory but abort integration once we have
        # diverged more than solver.max_elwise_diff_step
        if trajectory!(current_op, cs, current_strategy, last_op,
                       first(last_op.x), solver.max_elwise_diff_step)
            return true
        end
    end
    # in this case, the result in next_op is not really meaningful because the
    # integration has not been finished, thus the `success` state needs to be
    # evaluated and handled by the caller.
    return false
end

function solve(g::AbstractGame, solver::iLQSolver, args...)
    op0 = zero(SystemTrajectory, g)
    γ0 = zero(strategytype(g))
    return solve!(op0, γ0, g, solver, args...)
end

function solve(initial_strategy::StaticVector, g, args...)
    return solve!(zero(SystemTrajectory, g), copy(initial_strategy), g, args...)
end

function solve(initial_op::SystemTrajectory, initial_strategy::StaticVector, g,
               args...)
    return solve!(copy(initial_op), copy(initial_strategy), g, args...)
end

"Copy the `initial_op` to a new operting point instance using preallocated memory."
function copyop_prealloc!(solver::iLQSolver, initial_op::SystemTrajectory)
    om = solver._op_mem
    new_op = SystemTrajectory{samplingtime(om)}(om.x, om.u, initialtime(initial_op))
    return copyto!(new_op, initial_op)
end

"""

    $(TYPEDSIGNATURES)

Computes a solution solution to a (potentially non-linear and non-quadratic)
finite horizon game g.
"""
function solve!(initial_op::SystemTrajectory, initial_strategy::StaticVector,
                g::GeneralGame, solver::iLQSolver, x0::SVector,
                verbose::Bool=false)

    converged = false
    i_iter = 0
    last_op = copyop_prealloc!(solver, initial_op)

    current_op = initial_op
    current_strategy = initial_strategy
    lqg_approx = solver._lq_mem

    # 0. compute the operating point for the first run.
    # TODO -- we could probably allow to skip this in some warm-starting scenarios
    trajectory!(current_op, dynamics(g), current_strategy, last_op, x0)

    # ... and upate the current by integrating the non-linear dynamics
    while !(converged || i_iter >= solver.max_n_iter)
        # sanity chech to make sure that we don't manipulate the wrong
        # object...
        @assert !(last_op === current_op) "current and last operating point
        refer to the *same* object."

        # 1. linearize dynamics and quadratisize costs to obtain an lq game
        lq_approximation!(lqg_approx, solver, g, current_op)

        # 2. solve the current lq version of the game
        solve_lq_game!(current_strategy, lqg_approx)

        # 3. do line search to stabilize the strategy selection and extract the
        # next operating point
        copyto!(last_op, current_op)
        success = backtrack_scale!(current_strategy, current_op, last_op,
                                   dynamics(g), solver)
        if(!success)
            verbose && @warn "Could not stabilize solution."
            # we immetiately return and state that the solution has not been
            # stabilized
            return false, current_op, current_strategy
        end

        i_iter += 1
        converged = has_converged(solver, last_op, current_op)
    end

    # NOTE: for `converged == false` the result may not be meaningful. `converged`
    # has to be handled outside this function
    return converged, current_op, current_strategy
end


text:
const block_signs = [
    0b1000 0b0010
    0b0100 0b0001
]

const n_block = 16
const block_decode = Vector{Char}(undef, typemax(UInt16))
block_decode[0b0000 + 1] = ' '
block_decode[0b0001 + 1] = '▗'
block_decode[0b0010 + 1] = '▖'
block_decode[0b0011 + 1] = '▄'
block_decode[0b0100 + 1] = '▝'
block_decode[0b0101 + 1] = '▐'
block_decode[0b0110 + 1] = '▞'
block_decode[0b0111 + 1] = '▟'
block_decode[0b1000 + 1] = '▘'
block_decode[0b1001 + 1] = '▚'
block_decode[0b1010 + 1] = '▌'
block_decode[0b1011 + 1] = '▙'
block_decode[0b1100 + 1] = '▀'
block_decode[0b1101 + 1] = '▜'
block_decode[0b1110 + 1] = '▛'
block_decode[0b1111 + 1] = '█'
block_decode[(n_block + 1):typemax(UInt16)] = unicode_table[1:(typemax(UInt16) - n_block)]

"""
The `BlockCanvas` is also Unicode-based.
It has half the resolution of the `BrailleCanvas`.
In contrast to BrailleCanvas, the pixels don't
have visible spacing between them.
This canvas effectively turns every character
into 4 pixels that can individually be manipulated
using binary operations.
"""
struct BlockCanvas <: LookupCanvas
    grid::Array{UInt16,2}
    colors::Array{ColorType,2}
    blend::Bool
    visible::Bool
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
    xscale::Union{Symbol,Function}
    yscale::Union{Symbol,Function}
end

@inline x_pixel_per_char(::Type{BlockCanvas}) = 2
@inline y_pixel_per_char(::Type{BlockCanvas}) = 2

@inline lookup_encode(c::BlockCanvas) = block_signs
@inline lookup_decode(c::BlockCanvas) = block_decode

BlockCanvas(args...; nargs...) = CreateLookupCanvas(BlockCanvas, args...; nargs...)

function char_point!(
    c::BlockCanvas,
    char_x::Int,
    char_y::Int,
    char::Char,
    color::UserColorType,
)
    if checkbounds(Bool, c.grid, char_x, char_y)
        c.grid[char_x, char_y] = n_block + char
        set_color!(c.colors, char_x, char_y, crayon_256_color(color), c.blend)
    end
    c
end


text:
function construct_device!(canonical::CanonicalModel, sys::PSY.System,
                           model::DeviceModel{St, D},
                           ::Type{S};
                           kwargs...) where {St<:PSY.Storage,
                                             D<:AbstractStorageFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(St, sys)

    if validate_available_devices(devices, St)
        return
    end

    parameters = get(kwargs, :parameters, true)

    #Variables
    active_power_variables(canonical, devices)

    reactive_power_variables(canonical, devices)

    energy_storage_variables(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, D)

    #Constraints
    active_power_constraints(canonical, devices, D, S)

    reactive_power_constraints(canonical, devices, D, S)

    energy_capacity_constraints(canonical, devices, D, S)

    feedforward!(canonical, St, model.feedforward)

    # Energy Balanace limits
    energy_balance_constraint(canonical, devices, D, S)

    return

end

function construct_device!(canonical::CanonicalModel, sys::PSY.System,
                           model::DeviceModel{St, D},
                           ::Type{S};
                           kwargs...) where {St<:PSY.Storage,
                                             D<:AbstractStorageFormulation,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(St, sys)

    if validate_available_devices(devices, St)
        return
    end

    parameters = get(kwargs, :parameters, true)

    #Variables
    active_power_variables(canonical, devices);

    energy_storage_variables(canonical, devices);

    #Initial Conditions

    initial_conditions!(canonical, devices, D)

    #Constraints
    active_power_constraints(canonical, devices, D, S)

    energy_capacity_constraints(canonical, devices, D, S)

    feedforward!(canonical, St, model.feedforward)

    # Energy Balanace limits
    energy_balance_constraint(canonical, devices, D, S)

    return

end

function construct_device!(canonical::CanonicalModel, sys::PSY.System,
                                        model::DeviceModel{St, BookKeepingwReservation},
                                        ::Type{S};
                                        kwargs...) where {St<:PSY.Storage,
                                                          S<:PM.AbstractPowerModel}

    devices = PSY.get_components(St, sys)

    if validate_available_devices(devices, St)
        return
    end

    parameters = get(kwargs, :parameters, true)

    #Variables
    active_power_variables(canonical, devices)

    reactive_power_variables(canonical, devices)

    energy_storage_variables(canonical, devices)

    storage_reservation_variables(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    active_power_constraints(canonical, devices, model.formulation, S)

    reactive_power_constraints(canonical, devices, model.formulation, S)

    energy_capacity_constraints(canonical, devices, model.formulation, S)

    feedforward!(canonical, St, model.feedforward)

    # Energy Balanace limits
    energy_balance_constraint(canonical, devices, model.formulation, S)

    return

end

function construct_device!(canonical::CanonicalModel, sys::PSY.System,
                           model::DeviceModel{St, BookKeepingwReservation},
                           ::Type{S};
                           kwargs...) where {St<:PSY.Storage,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(St, sys)

    if validate_available_devices(devices, St)
        return
    end

    parameters = get(kwargs, :parameters, true)

    #Variables
    active_power_variables(canonical, devices)

    energy_storage_variables(canonical, devices)

    storage_reservation_variables(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    active_power_constraints(canonical, devices, model.formulation, S)

    energy_capacity_constraints(canonical, devices, model.formulation, S)

    feedforward!(canonical, St, model.feedforward)

    # Energy Balanace limits
    energy_balance_constraint(canonical, devices, model.formulation, S)

    return

end


text:
using Optim
# using MPI_Tools
##This uses equations presented in "Stray capacitances of single-layer air-core inductors for high-frequency applications"
# By"Grandi, G., Kazimierczuk, M. K.,Massarini, A. Reggiani, U."

function F(p,t,r,ϵr)
    return (p/(2*r))./((1+t/r)^(1-1/ϵr))
end
function Ct(D,r,p,t,ϵr)
    #D - solenoid diameter in meters
    #r - wire radius
    #p - pitch in meters
    #t - wire insulation thickness
    #ϵr - insulation relative permittivity

    ϵ0 = 8.85e-12

    FVal = F(p,t,r,ϵr)
    CtVal = pi^2*D*ϵ0 / (log(FVal+sqrt(FVal^2 - (1+t/r)^(2/ϵr))))
end

## Using derivation in appendix
function CC(r,t,ϵr)
    return 2*π*ϵ₀*ϵr / (log(1+t/r))
end

function CG(p,r,t)
    (π * ϵ₀) ./ log( (p/(2*r))/(1+t/r) + sqrt( ( (p/(2*r))/(1+t/r) )^2 - 1) )
end

"""
This function calculates the equivalent turn-turn capacitcance of a single layer solenoid
    #D - solenoid diameter in meters
    #r - wire radius
    #p - pitch in meters
    #t - wire insulation thickness
    #ϵr - insulation relative permittivity
"""
function Ceq(D,r,p,t,ϵr)
    CC(r,t,ϵr)*CG(p,r,t) / (CC(r,t,ϵr)+2*CG(p,r,t)) * π*D
end

# Ceq(.2,.15e-3,1.33e-3,.5e-3,4)/10

## Medhurst
function CMedh(L,D)
    ϵ₀ = 8.85*1e-12
    return (4*ϵ₀*L/π + 8e-12*D + 27e-12*D*sqrt(D/L)) #Farads
end

function CMedh_Knight(L,D)
    ϵ₀ = 8.85*1e-12
    return 4*ϵ₀*L/π*(1+ 0.8249* (D/L)+ 2.329*(D/L)^(3/2))#Farads
end

## Equation 5.3, Cₗ-DAE in Knight, p.52

function CL_DAE_Knight(
    L, #Coil length (axial)
    D, #Coil diamter
    N, #N turns
    ϵᵣₓ, #external permittivity
    ϵᵣᵢ) #Internal permittivity
    kc(L,D) =0.717439(D/L) +0.933048(D/L)^(3/2) +0.106(D/L)^2

    ϵ₀ = 8.85*1e-12
    p = L/N

    Ψ = atan(p/(π*D))

    C_L = (4*ϵ₀*ϵᵣₓ)/π *L * (1+kc(L,D)*(1+ϵᵣᵢ/ϵᵣₓ)/2)*1/(cos(Ψ)^2)
    return C_L
end


text:
using Distributions
using Base.Test
import Distributions: CategoricalDirectSampler, AliasTable

function test_categoricalsampler(s, p::Vector{Float64}, n::Int, tol::Float64)
    # s: the sampler to be tested
    # p: the reference probability vector
    # n: the number of samples to be generated
    # tol: maximal allowable deviation

    K = length(p)
    @assert ncategories(s) == K

    cnts = zeros(Int, K)

    for i = 1:n
        k = rand(s)
        @assert 1 <= k <= K
        cnts[k] += 1
    end

    q = cnts ./ n
    # println(Base.maxabs(p - q))
    @test_approx_eq_eps p q tol
end

@test_throws ErrorException CategoricalDirectSampler(Float64[])
@test_throws ErrorException AliasTable(Float64[])

for S in [CategoricalDirectSampler, AliasTable]
    s = S([1.0])
    for i = 1:10
        @test rand(s) == 1
    end
    for p in {[0.3, 0.7], [0.2, 0.3, 0.4, 0.1]}
        test_categoricalsampler(S(p), p, 10^5, 0.015)
    end
end


text:
# This file was generated, do not modify it. # hide
rpt2 = report(mach_tuned_iterated_pipe);
best_booster = rpt2.best_model.model.evo_tree_classifier

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using HDF5

runtest(filename) = (println(filename); include(filename))

runtest("jldtests.jl")
runtest("require.jl")
runtest("custom_serialization.jl")
runtest("type_translation.jl")
if Pkg.installed("DataFrames") != nothing
    runtest("jld_dataframe.jl")
end


text:
#!/usr/bin/env julia

# Publish simulated data continuously over a ZMQ socket.
# Binds PUB socket to tcp://*:5556

using ZMQ
using JSON

context = Context()
socket = Socket(context, PUB)
ZMQ.bind(socket, "tcp://*:5556")

time_ms() = round(Int, time_ns()/1e6)

alpha = 0.3 # Smoothing parameter in [0,1] (smaller is stronger)
start = time_ms()
ylist = zeros(2)

println("Sending data...")
while true

    # Random walk
    y = clamp(ylist[1] + 0.1*randn(), -1, 1)

    # Filtered value is a convex sum of current and historical time series data
    # http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc431.htm
    s = alpha*y + (1 - alpha)*ylist[2]

    ylist = [y,s]

    # Print to console and send out the message string in JSON format
    msg = JSON.json(Dict("t"=>time_ms() - start, "y"=>y, "s"=>s))
    println(msg)
    ZMQ.send(socket, msg)

    sleep(0.02) # ~50 Hz
end

# Never reached
ZMQ.close(socket)
ZMQ.close(context)

max_stars_repo_name: andrewadare/realtime

max_issues_repo_name: andrewadare/realtime

max_forks_repo_name: andrewadare/realtime


text:
struct MockDofBasis{P} <: Dof
  x::Vector{P}
end

function dof_cache(b::MockDofBasis,field)
  field_cache(field,b.x)
end

@inline function evaluate_dof!(cache,b::MockDofBasis,field)
  evaluate_field!(cache,field,b.x)
end


text:
"""
    GaussLik <: Likelihood

Gaussian, a.k.a. Normal, likelihood
```math
p(y | f, σ) = 1 / √(2πσ²) \\exp(-(y - f)²/(2σ²)),
```
where standard deviation ``σ`` is a non-fixed hyperparameter and ``f`` is the latent
Gaussian process.
"""
mutable struct GaussLik <: Likelihood
    "Standard deviation"
    σ::Float64
    "Priors for likelihood parameters"
    priors::Array

    """
        GaussLik(lσ::Float64)

    Create `GaussLik` with standard deviation `exp(lσ)`.
    """
    GaussLik(lσ::Float64) = new(exp(lσ), [])
end

# log of probability density
function log_dens(gauss::GaussLik, f::AbstractVector, y::AbstractVector)
    return (-0.5 * log(2 * pi) - log(gauss.σ)) .- 0.5 * ((y - f) / gauss.σ).^2
end

# derivative of log pdf wrt latent function
function dlog_dens_df(gauss::GaussLik, f::AbstractVector, y::AbstractVector)
    return [(yi-fi)/gauss.σ^2 for (fi,yi) in zip(f,y)]
end

# derivative of log pdf wrt to parameters
function dlog_dens_dθ(gauss::GaussLik, f::AbstractVector, y::AbstractVector)
    return gauss.σ*[-1/gauss.σ + 1/gauss.σ^3*(yi-fi).^2 for (fi,yi) in zip(f,y)]
end

#mean and variance under likelihood
mean_lik(gauss::GaussLik, f::AbstractVector) = f
var_lik(gauss::GaussLik, f::AbstractVector) = ones(length(f))*gauss.σ^2

function set_params!(gauss::GaussLik, hyp::AbstractVector)
    length(hyp) == 1 || throw(ArgumentError("Gaussian/Normal likelihood has only one free parameter"))
    gauss.σ = exp(hyp[])
end

get_params(gauss::GaussLik) = Float64[log(gauss.σ)]
num_params(gauss::GaussLik) = 1

#Computes the predictive mean and variance
function predict_obs(gauss::GaussLik, fmean::AbstractVector, fvar::AbstractVector)
    return fmean, fvar + gauss.σ^2
end

function var_exp(ll::GaussLik, y::AbstractArray, m::AbstractArray, V::AbstractMatrix)
    tot = 0
    V_diag = diag(V)
    for (a, b, c) in zip(y, m, V_diag)
        tot +=  var_exp(ll, a, b, c) # convert to lgamma(y+1)
    end
    return tot
end

function var_exp(ll::GaussLik, y::AbstractArray, m::AbstractArray, V::AbstractArray)
    tot = 0
    for (a, b, c) in zip(y, m, V)
        tot +=  var_exp(ll, a, b, c)
    end
    return tot
end

function var_exp(ll::GaussLik, y::Number, m::Number, V::Number)
    return -0.5*log(2*π) - 0.5*log(ll.σ) - 0.5*((y - m)^2 + V)/ll.σ
end

function dv_var_exp(ll::GaussLik, y::Number, m::Number, V::Number)
    return gradient(x -> var_exp(ll, y, m, x), V)[1]
end

function dm_var_exp(ll::GaussLik, y::Number, m::Number, V::Number)
    return gradient(x -> var_exp(ll, y, x, V), m)[1]
end


text:
using Plots, Test, Dates

@testset "Limits" begin
    y = [1.0 * i * i for i in 1:10]
    x = [Date(2019, 11, i) for i in 1:10]

    rx = [x[3], x[5]]

    p = plot(x, y, widen = false)
    vspan!(p, rx, label = "", alpha = 0.2)

    ref_ylims = (y[1], y[end])
    ref_xlims = (x[1].instant.periods.value, x[end].instant.periods.value)
    @test Plots.ylims(p) == ref_ylims
    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset

@testset "Date xlims" begin
    y = [1.0 * i * i for i in 1:10]
    x = [Date(2019, 11, i) for i in 1:10]
    span = (Date(2019, 10, 31), Date(2019, 11, 11))

    ref_xlims = map(date -> date.instant.periods.value, span)

    p = plot(x, y, xlims = span, widen = false)

    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset

@testset "DateTime xlims" begin
    y = [1.0 * i * i for i in 1:10]
    x = [DateTime(2019, 11, i, 11) for i in 1:10]
    span = (DateTime(2019, 10, 31, 11, 59, 59), DateTime(2019, 11, 11, 12, 01, 15))

    ref_xlims = map(date -> date.instant.periods.value, span)

    p = plot(x, y, xlims = span, widen = false)
    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# %%
using Plots
x = range(-5, 5; length=300)
y = range(-3, 3; length=200)
f(x, y) = exp(-x^2/2-y^2)
z = f.(x', y)
surface(x, y, z; camera=(30, 60), color=:gist_earth)

# %%
surface(x, y, f; camera=(30, 60), color=:gist_earth)

# %%

max_stars_repo_path: 0016/3d plot example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/3d plot example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/3d plot example.jl

max_forks_repo_name: genkuroki/public