text:
# # Differentiating a QP wrt a single variable

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/matrix-inversion-manual.jl)

# Consider the quadratic program

# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & \frac{1}{2} x^T Q x + q^T x \\
# \mbox{subject to} & G x \leq h, x \in \mathcal{R}^2 \\
# \end{array}
# \end{split}
# ```

# where `Q`, `q`, `G` are fixed and `h` is the single parameter.

# In this example, we'll try to differentiate the QP wrt `h`, by finding its
# jacobian by hand (using Eqn (6) of [QPTH article](https://arxiv.org/pdf/1703.00443.pdf))
# and compare the results:
# - Manual compuation
# - Using JuMP and DiffOpt

# Assuming
# ```
# Q = [[4, 1], [1, 2]]
# q = [1, 1]
# G = [1, 1]
# ```
# and begining with a starting value of `h=-1`

# few values just for reference

# | variable | optimal value | note |
# |----|------|-----|
# | x* | [-0.25; -0.75] | Primal optimal |
# | 𝜆∗ | -0.75 | Dual optimal |

# ## Finding Jacobian using matrix inversion
# Lets formulate Eqn (6) of [QPTH article](https://arxiv.org/pdf/1703.00443.pdf) for our QP. If we assume `h` as the only parameter and `Q`,`q`,`G` as fixed problem data - also note that our QP doesn't involves `Ax=b` constraint - then Eqn (6) reduces to
# ```math
# \begin{gather}
#  \begin{bmatrix}
#      Q & g^T \\
#      \lambda^* G & G x^* - h
#  \end{bmatrix}
#  \begin{bmatrix}
#      dx \\
#      d \lambda
#  \end{bmatrix}
#  =
#   \begin{bmatrix}
#    0 \\
#    \lambda^* dh
#    \end{bmatrix}
# \end{gather}
# ```

# Now to find the jacobians $$ \frac{\partial x}{\partial h}, \frac{\partial \lambda}{\partial h}$$
# we substitute `dh = I = [1]` and plug in values of `Q`,`q`,`G` to get
# ```math
# \begin{gather}
#  \begin{bmatrix}
#      4 & 1 & 1 \\
#      1 & 2 & 1 \\
#      -0.75 & -0.75 & 0
#  \end{bmatrix}
#  \begin{bmatrix}
#      \frac{\partial x_1}{\partial h} \\
#      \frac{\partial x_2}{\partial h} \\
#      \frac{\partial \lambda}{\partial h}
#  \end{bmatrix}
#  =
#   \begin{bmatrix}
#    0 \\
#    0 \\
#    -0.75
#    \end{bmatrix}
# \end{gather}
# ```

# Upon solving using matrix inversion, the jacobian is
# ```math
# \frac{\partial x_1}{\partial h} = 0.25, \frac{\partial x_2}{\partial h} = 0.75, \frac{\partial \lambda}{\partial h} = -1.75
# ```

# ## Finding Jacobian using JuMP and DiffOpt

using JuMP
import DiffOpt
import Ipopt

n = 2 # variable dimension
m = 1; # no of inequality constraints

Q = [4. 1.;1. 2.]
q = [1.; 1.]
G = [1. 1.;]
h = [-1.;]   # initial values set

# Initialize empty model

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)

# Add the variables

@variable(model, x[1:2])

# Add the constraints.

@constraint(
    model,
    cons[j in 1:1],
    sum(G[j, i] * x[i] for i in 1:2)  <= h[j]
);

@objective(
    model,
    Min,
    1/2 * sum(Q[j, i] * x[i] *x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)

# Solve problem

optimize!(model)

# primal solution

value.(x)

# dual solution

dual.(cons)

# set sentivitity

MOI.set(
    model,
    DiffOpt.ForwardInConstraint(),
    cons[1],
    0.0 * index(x[1]) - 1.0,  # to indicate the direction vector to get directional derivatives
)

# Note that `0.0 * index(x[1])` is used to make its type `typeof(0.0 * index(x[1]) - 1.0) <: MOI.AbstractScalarFunction`.
# To indicate different direction to get directional derivative, users should replace `0.0 * index(x[1]) - 1.0` as the form of `dG*x - dh`, where `dG` and `dh` correspond to the elements of direction vectors along `G` and `h` axes, respectively.

# Compute derivatives

DiffOpt.forward(model)

# Query derivative

dx = MOI.get.(
    model,
    DiffOpt.ForwardOutVariablePrimal(),
    x,
)

using Test                                  #src
@test dx ≈ [0.25 ,0.75] atol=1e-4 rtol=1e-4 #src


text:
using Clang.Generators
using libportaudio_jll

cd(@__DIR__)

include_dir = joinpath(libportaudio_jll.artifact_dir, "include") |> normpath
portaudio_h = joinpath(include_dir, "portaudio.h")

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

ctx = create_context(portaudio_h, args, options)

build!(ctx)

ext: jl


text:
"""
$(TYPEDEF)

Store both grid and density for KDE over ``ℝ²``.

Reading the fields directly is part of the API, and

```julia
sum(density) * step(x) ≈ 1
```

# Fields

$(FIELDS)
"""
mutable struct UnivariateKDE{R<:AbstractRange} <: AbstractKDE
    "Gridpoints for evaluating the density."
    x::R
    "Kernel density at corresponding gridpoints `x`."
    density::Vector{Float64}
end

# construct kernel from bandwidth
kernel_dist(::Type{Normal},w::Real) = Normal(0.0,w)
kernel_dist(::Type{Uniform},w::Real) = (s = 1.7320508075688772*w; Uniform(-s,s))

const LocationScale = Union{Laplace,Logistic,SymTriangularDist}
kernel_dist(::Type{D},w::Real) where {D} = (s = w/std(D(0.0,1.0)); D(0.0,s))

# Silverman's rule of thumb for KDE bandwidth selection
function default_bandwidth(data::RealVector, alpha::Float64 = 0.9)
    # Determine length of data
    ndata = length(data)
    ndata <= 1 && return alpha

    # Calculate width using variance and IQR
    var_width = std(data)
    q25, q75 = quantile(data, [0.25, 0.75])
    quantile_width = (q75 - q25) / 1.34

    # Deal with edge cases with 0 IQR or variance
    width = min(var_width, quantile_width)
    if width == 0.0
        if var_width == 0.0
            width = 1.0
        else
            width = var_width
        end
    end

    # Set bandwidth using Silverman's rule of thumb
    return alpha * width * ndata^(-0.2)
end

function default_weights(data::RealVector)
    UniformWeights(length(data))
end

# Roughly based on:
#   B. W. Silverman (1982) "Algorithm AS 176: Kernel Density Estimation Using
#   the Fast Fourier Transform", Journal of the Royal Statistical
#   Society. Series C (Applied Statistics) , Vol. 31, No. 1, pp. 93-99
#   URL: http://www.jstor.org/stable/2347084
# and:
#   M. C. Jones and H. W. Lotwick (1984) "Remark AS R50: A Remark on Algorithm
#   AS 176. Kernal Density Estimation Using the Fast Fourier Transform",
#   Journal of the Royal Statistical Society. Series C (Applied Statistics) ,
#   Vol. 33, No. 1, pp. 120-122
#   URL: http://www.jstor.org/stable/2347674

# default kde range
# Should extend enough beyond the data range to avoid cyclic correlation from the FFT
function kde_boundary(data::RealVector, bandwidth::Real)
    lo, hi = extrema(data)
    lo - 4.0*bandwidth, hi + 4.0*bandwidth
end

# convert boundary and npoints to Range object
function kde_range(boundary::Tuple{Real,Real}, npoints::Int)
    lo, hi = boundary
    lo < hi || error("boundary (a,b) must have a < b")

    range(lo, stop=hi, length=npoints)
end

struct UniformWeights{N} end

UniformWeights(n) = UniformWeights{n}()

Base.sum(x::UniformWeights) = 1.
Base.getindex(x::UniformWeights{N}, i) where {N} = 1/N

const Weights = Union{UniformWeights, RealVector, StatsBase.Weights}

# tabulate data for kde
function tabulate(data::RealVector, midpoints::R, weights::Weights=default_weights(data)) where R<:AbstractRange
    npoints = length(midpoints)
    s = step(midpoints)

    # Set up a grid for discretized data
    grid = zeros(Float64, npoints)
    ainc = 1.0 / (sum(weights)*s*s)

    # weighted discretization (cf. Jones and Lotwick)
    for (i,x) in enumerate(data)
        k = searchsortedfirst(midpoints,x)
        j = k-1
        if 1 <= j <= npoints-1
            grid[j] += (midpoints[k]-x)*ainc*weights[i]
            grid[k] += (x-midpoints[j])*ainc*weights[i]
        end
    end

    # returns an un-convolved KDE
    UnivariateKDE(midpoints, grid)
end

# convolve raw KDE with kernel
# TODO: use in-place fft
function conv(k::UnivariateKDE, dist::UnivariateDistribution)
    # Transform to Fourier basis
    K = length(k.density)
    ft = rfft(k.density)

    # Convolve fft with characteristic function of kernel
    # empirical cf
    #  = \sum_{n=1}^N e^{i*t*X_n} / N
    #  = \sum_{k=0}^K e^{i*t*(a+k*s)} N_k / N
    #  = e^{i*t*a} \sum_{k=0}^K e^{-2pi*i*k*(-t*s*K/2pi)/K} N_k / N
    #  = A * fft(N_k/N)[-t*s*K/2pi + 1]
    c = -twoπ/(step(k.x)*K)
    for j = 0:length(ft)-1
        ft[j+1] *= cf(dist,j*c)
    end

    dens = irfft(ft, K)
    # fix rounding error.
    for i = 1:K
        dens[i] = max(0.0,dens[i])
    end

    # Invert the Fourier transform to get the KDE
    UnivariateKDE(k.x, dens)
end

# main kde interface methods
function kde(data::RealVector, weights::Weights, midpoints::R, dist::UnivariateDistribution) where R<:AbstractRange
    k = tabulate(data, midpoints, weights)
    conv(k,dist)
end

function kde(data::RealVector, dist::UnivariateDistribution;
             boundary::Tuple{Real,Real}=kde_boundary(data,std(dist)), npoints::Int=2048, weights=default_weights(data))

    midpoints = kde_range(boundary,npoints)
    kde(data,weights,midpoints,dist)
end

function kde(data::RealVector, midpoints::R;
             bandwidth=default_bandwidth(data), kernel=Normal, weights=default_weights(data)) where R<:AbstractRange
    bandwidth > 0.0 || error("Bandwidth must be positive")
    dist = kernel_dist(kernel,bandwidth)
    kde(data,weights,midpoints,dist)
end

function kde(data::RealVector; bandwidth=default_bandwidth(data), kernel=Normal,
             npoints::Int=2048, boundary::Tuple{Real,Real}=kde_boundary(data,bandwidth), weights=default_weights(data))
    bandwidth > 0.0 || error("Bandwidth must be positive")
    dist = kernel_dist(kernel,bandwidth)
    kde(data,dist;boundary=boundary,npoints=npoints,weights=weights)
end

"""
    optimize(f, x_lower, x_upper; iterations=1000, rel_tol=nothing, abs_tol=nothing)

Minimize the function `f` in the interval `x_lower..x_upper`, using the
[golden-section search](https://en.wikipedia.org/wiki/Golden-section_search).
Return an approximate minimum `x̃` or error if such approximate minimum cannot be found.

This algorithm assumes that `-f` is unimodal on the interval `x_lower..x_upper`,
that is to say, there exists a unique `x` in `x_lower..x_upper` such that `f` is
decreasing on `x_lower..x` and increasing on `x..x_upper`.

`rel_tol` and `abs_tol` determine the relative and absolute tolerance, that is
to say, the returned value `x̃` should differ from the actual minimum `x` at most
`abs_tol + rel_tol * abs(x̃)`.
If not manually specified, `rel_tol` and `abs_tol` default to `sqrt(eps(T))` and
`eps(T)` respectively, where `T` is the floating point type of `x_lower` and `x_upper`.

`iterations` determines the maximum number of iterations allowed before convergence.

This is a private, unexported function, used internally to select the optimal bandwidth
automatically.
"""
function optimize(f, x_lower, x_upper; iterations=1000, rel_tol=nothing, abs_tol=nothing)

    if x_lower > x_upper
        error("x_lower must be less than x_upper")
    end

    T = promote_type(typeof(x_lower/1), typeof(x_upper/1))
    rtol = something(rel_tol, sqrt(eps(T)))
    atol = something(abs_tol, eps(T))

    function midpoint_and_convergence(lower, upper)
        midpoint = (lower + upper) / 2
        tol = atol + rtol * midpoint
        midpoint, (upper - lower) <= 2tol
    end

    invphi::T = 0.5 * (sqrt(5) - 1)
    invphisq::T = 0.5 * (3 - sqrt(5))

    a::T, b::T = x_lower, x_upper
    h = b - a
    c = a + invphisq * h
    d = a + invphi * h

    fc, fd = f(c), f(d)

    for _ in 1:1000
        h *= invphi
        if fc < fd
            m, converged = midpoint_and_convergence(a, d)
            converged && return m
            b = d
            d, fd = c, fc
            c = a + invphisq * h
            fc = f(c)
        else
            m, converged = midpoint_and_convergence(c, b)
            converged && return m
            a = c
            c, fc = d, fd
            d = a + invphi * h
            fd = f(d)
        end
    end

    error("Reached maximum number of iterations without convergence.")
end

# Select bandwidth using least-squares cross validation, from:
#   Density Estimation for Statistics and Data Analysis
#   B. W. Silverman (1986)
#   sections 3.4.3 (pp. 48-52) and 3.5 (pp. 61-66)

function kde_lscv(data::RealVector, midpoints::R;
                  kernel=Normal,
                  bandwidth_range::Tuple{Real,Real}=(h=default_bandwidth(data); (0.25*h,1.5*h)),
                  weights=default_weights(data)) where R<:AbstractRange

    ndata = length(data)
    k = tabulate(data, midpoints, weights)

    # the ft here is K/ba*sqrt(2pi) * u(s), it is K times the Yl in Silverman's book
    K = length(k.density)
    ft = rfft(k.density)

    ft2 = abs2.(ft)
    c = -twoπ/(step(k.x)*K)
    hlb, hub = bandwidth_range

    minimizer = optimize(hlb, hub) do h
        dist = kernel_dist(kernel, h)
        ψ = 0.0
        for j = 1:length(ft2)-1
            ks = real(cf(dist, j*c))
            ψ += ft2[j+1]*(ks-2.0)*ks
        end
        ψ*step(k.x)/K + pdf(dist,0.0)/ndata
    end

    dist = kernel_dist(kernel, minimizer)
    for j = 0:length(ft)-1
        ft[j+1] *= cf(dist, j*c)
    end

    dens = irfft(ft, K)
    # fix rounding error.
    for i = 1:K
        dens[i] = max(0.0,dens[i])
    end

    # Invert the Fourier transform to get the KDE
    UnivariateKDE(k.x, dens)
end

function kde_lscv(data::RealVector;
                  boundary::Tuple{Real,Real}=kde_boundary(data,default_bandwidth(data)),
                  npoints::Int=2048,
                  kernel=Normal,
                  bandwidth_range::Tuple{Real,Real}=(h=default_bandwidth(data); (0.25*h,1.5*h)),
                  weights::Weights = default_weights(data))

    midpoints = kde_range(boundary,npoints)
    kde_lscv(data,midpoints; kernel=kernel, bandwidth_range=bandwidth_range, weights=weights)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

make_value{T<:Integer}(::Type{T}, i::Integer) = 3*i%T
make_value{T<:AbstractFloat}(::Type{T},i::Integer) = T(3*i)

Vec{N,T} = NTuple{N,Base.VecElement{T}}

# Crash report for #15244 motivated this test.
@generated function thrice_iota{N,T}(::Type{Vec{N,T}})
    :(tuple($([:(Base.VecElement(make_value($T,$i))) for i in 1:N]...)))
end

function call_iota(n::Integer,t::DataType)
    x = thrice_iota(Vec{n,t})
    @test x[1].value === make_value(t,1)
    @test x[n].value === make_value(t,n)
end

# Try various tuple lengths and element types
for i=1:20
    for t in [Bool, Int8, Int16, Int32, Int64, Float32, Float64]
        call_iota(i,t)
    end
end

# Another crash report for #15244 motivated this test.
struct Bunch{N,T}
    elts::NTuple{N,Base.VecElement{T}}
end

unpeel(x) = x.elts[1].value

@test unpeel(Bunch{2,Float64}((Base.VecElement(5.0),
                               Base.VecElement(4.0)))) === 5.0

rewrap(x) = VecElement(x.elts[1].value+0)
b = Bunch((VecElement(1.0), VecElement(2.0)))

@test rewrap(b)===VecElement(1.0)

struct Herd{N,T}
    elts::NTuple{N,Base.VecElement{T}}
    Herd{N,T}(elts::NTuple{N,T}) where {N,T} = new(ntuple(i->Base.VecElement{T}(elts[i]), N))
end

function check{N,T}(x::Herd{N,T})
    for i=1:N
        @test x.elts[i].value === N*N+i-1
    end
end

check(Herd{1,Int}((1,)))
check(Herd{2,Int}((4,5)))
check(Herd{4,Int}((16,17,18,19)))

struct Gr{N, T}
    u::T
    v::Bunch{N,T}
    w::T
end

a = Vector{Gr{2,Float64}}(2)
a[2] = Gr(1.0, Bunch((VecElement(2.0), VecElement(3.0))), 4.0)
a[1] = Gr(5.0, Bunch((VecElement(6.0), VecElement(7.0))), 8.0)
@test a[2] == Gr(1.0, Bunch((VecElement(2.0), VecElement(3.0))), 4.0)

@test isa(VecElement((1,2)), VecElement{Tuple{Int,Int}})


text:
################################################################################
# hx_core.jl
# v 0.5 by Richard Fitzhugh
# Last edited on 10/17/2017
#
# This file contains functions for overall heat exchanger analysis using the
# methods found in Fundamentals of Heat and Mass Transfer by Incropera and
# DeWitt, published by John Wiley and Sons, Inc. The book is standard for
# teaching heat and mass trasfer to mechanical and chemical engineers in the
# United States. The chapter (11) and page numbers used in documentation
# refer to the sixth edition but with minor changes could probably be used for
# other editions as well.
#
# This module uses the ϵ-NTU method because it's much more robust numerically
# than the LMTD method. That is, the equations are formulated in terms of the '
# exponential rather than logarithmic functions, so they don't start returning
# NaNs if a numerical solver happens to guess negative temperatures or the like.
################################################################################

@enum hx_type conc_parallel=0 conc_counter=1 st_oneshell=2 st_twoshell=3 Cr_0=4

struct HeatExchanger
    oldshell::Pipe
    shell::Pipe
    tube::Pipe
    n_tubes::Int
    bafflespace::AbstractFloat
    bafflecut::AbstractFloat
    type_::hx_type
end

function HeatExchanger(shell::Pipe, tube::Pipe, ntubes::Integer,
                        bsp, bcut, type_::hx_type)
    farea = shell.flowarea - ntubes * tube.totalarea
    PW = pi*shell.diameter + ntubes*pi*tube.outer_diameter
    Dh = 4*farea/PW
    newshell = Pipe(shell.length, Dh, Dh/2, farea, shell.thickness,
        shell.outer_diameter, shell.outer_radius, shell.totalarea,
        shell.roughness, shell.elevation_change, shell.therm_cond)
    return HeatExchanger(shell, newshell, tube, ntubes, bsp, bcut, type_)
end

"""
The following are implementations of the effectiveness (ϵ) relations for various
heat exchanger geometries. These equations can be found in Incropera and DeWitt,
Table 11.3, page 689 (6th Ed.)
"""
function ϵ_Cp(NTU, Cr)
    top = 1 - exp(-NTU *(1+Cr))
    return top/(1+Cr)
end

function ϵ_Cc(NTU, Cr)

    if Cr==1
        return NTU/(1+NTU)
    end

    top = 1 - exp(-NTU *(1+Cr))
    bottom = 1 - Cr*exp(-NTU *(1+Cr))
    return top/bottom

end

function ϵ_St₁(NTU, Cr)
    term = exp(-NTU * √(1+Cr^2))
    ans = 2/(1 + Cr + √(1+Cr^2) * (1+term)/(1-term))
    return ans
end

function ϵ_Stₙ(NTU, Cr, n=2)
    term = ((1 - ϵ_St₁(NTU, Cr)*Cr) / (1 - ϵ_St₁(NTU, Cr)))^n
    return (term - 1)/(term - Cr)
end

function ϵ_Cr0(NTU, Cr=0)
    return 1 - exp(-NTU)
end

ϵ_functions = Dict(conc_parallel => ϵ_Cp,
                   conc_counter => ϵ_Cc,
                   st_oneshell => ϵ_St₁,
                   st_twoshell => ϵ_Stₙ,
                   Cr_0 => ϵ_Cr0)

function ϵ(NTU, Cr, hx::hx_type)
    return ϵ_functions[hx](NTU, Cr)
end

"""
Heat transfer correlation for liquid metals, applicable for a uniform heat flux
through the heat transfer surface, with Re between 3600 and 90500, and Pe
between 100 and 10000. Proposed by Skupinski et al. in Int. Journal of Heat and
Mass Transfer, Vol 8 p. 937, 1965.
"""
function lm_correlation(astream, Twall)
    return 4.82+0.0185*(astream.tstream.Re*astream.tstream.Pr)^0.827
end

"""
Gneilinski heat transfer correlation.
A workhorse for turbulent convective heat transfer through a tube.
"""
function gneilinski(astream, Twall)
    f = astream.f
    Re  = astream.Re
    Pr = astream.Pr
    top = (f/8)*(Re-1000)*Pr
    bottom = 1+12.7*√(f/8)*(Pr^(2/3)-1)
    return top/bottom
end

"""
Correlation for molten salt flowing through the shell side
of a segmental-baffled shell-and-tube HX. From
B.-C. Du et al. / International Journal of Heat and Mass Transfer 113 (2017) 456–465
"""
function moltensalt_shell(astream, Twall)
    return 0.0676*  astream.Re^ 0.70413*  astream.Pr^ 0.4
end

"""
Laminar, fully developed Nusselt correlation for a wall @ const heat flux
"""
function laminar_tube(astream, Twall)
    return 4.36
end

function get_correlation(r::FlowStream, descriptor="shell")

  sname = r.tstream.subst_.name
  key_ = (sname, descriptor)
  if haskey(correlations, key_)
    return correlations[key_]
  else
    return gneilinski
  end
end

"""
The purpose of this function is to provide an initial estimate of heat exchanger
performance. It calculates the heat exchanger effectiveness from the ϵ-NTU
method but doesn't take into account how the fluid proerties of either stream
change down the length of the exchanger. The outlet conditions obtained from
this analysis are fed as initial guesses to the solver function that appears
later (called hx_performance_solver)
"""
function hx_performance_calcs(hx_::HeatExchanger,
                              shellstream::Stream,
                              tubestream::Stream,
                              scorr::Function, tcorr::Function)
    # scorr and tcorr are the functions that correspond to the shell-side and
    # the tube-side of the heat exchanger, respectively. The user needs to
    # choose them correctly. Automatic assignment of heat transfer correlations
    # is something I would like someday, but will be really hard to do.
    # They must accept the FlowStream object and the wall temp as arguments,
    # and output the Nusselt number.
    # The wall temp because properties at the wall are needed for some
    # correlations.
    if shellstream.T > tubestream.T
        hot = shellstream
        cold = tubestream
    else
        hot = tubestream
        cold = shellstream
    end

    Chot = hot.ṁ * hot.Cp
    Ccold = cold.ṁ * cold.Cp

    if Chot > Ccold
        Cmin = Ccold
        Cmax = Chot
    else
        Cmin = Chot
        Cmax = Ccold
    end

    Cr = Cmin/Cmax
    Twall = (hot.T+cold.T)/2
    k = hx_.tube.therm_cond((hot.T+cold.T)/2)
    sflow = FlowStream(hx_.shell, shellstream)
    tstr_onetube = Stream(tubestream.subst_, tubestream.ṁ/hx_.n_tubes, tubestream.p, tubestream.T)
    tflow = FlowStream(hx_.tube, tstr_onetube)
    Nu_s = scorr(sflow, Twall)
    h_o = Nu_s*shellstream.k/hx_.shell.diameter
    Nu_t = tcorr(tflow, Twall)
    h_i = Nu_t*tubestream.k/hx_.tube.diameter

    Ai = π*hx_.tube.diameter*hx_.tube.length*hx_.n_tubes
    Ao = π*hx_.tube.outer_diameter*hx_.tube.length*hx_.n_tubes

    R_cond = log(hx_.tube.outer_diameter / hx_.tube.diameter) / (2π*hx_.tube.length*k*hx_.n_tubes)
    UA = 1/(1/(h_o*Ao) + 1/(h_i*Ai) + R_cond)
    NTU = UA/Cmin
    ϵ_ = ϵ(NTU, Cr, hx_.type_)
    qmax = Cmin*(hot.T - cold.T)
    q = ϵ_*qmax
    Th_out = hot.T - q/Chot
    Tc_out = cold.T+ q/Ccold
    println("Heat transferred: ",q," W")

    if hot==shellstream
        hout = Stream(hot.subst_, hot.ṁ, hot.p - sflow.ΔP, Th_out)
        cout = Stream(cold.subst_,cold.ṁ,cold.p- tflow.ΔP, Tc_out)
    else
        hout = Stream(hot.subst_, hot.ṁ, hot.p - tflow.ΔP, Th_out)
        cout = Stream(cold.subst_,cold.ṁ,cold.p- sflow.ΔP, Tc_out)
    end
    return [hout, cout]
end

"""
Function hx_performance_solver
This function takes the preliminary guesses provided by the previous function
hx_performance_calcs and uses them to find a solution
"""
function hx_performance_solver(hx_::HeatExchanger,
                               shellstream::Stream,
                               tubestream::Stream,
                               scorr::Function, tcorr::Function)
    # scorr and tcorr are the functions that correspond to the shell-side and
    # the tube-side of the heat exchanger, respectively. The user needs to
    # choose them correctly. Automatic assignment of heat transfer correlations
    # is something I would like someday, but will be really hard to do.
    # They must accept the FlowStream object and the wall temp as arguments,
    # and output the Nusselt number.
    # The wall temp because properties at the wall are needed for some
    # correlations.
    out_init = hx_performance_calcs(hx_, shellstream, tubestream, scorr, tcorr)
    hout_g = out_init[1]
    cout_g = out_init[2]
    x0 = [hout_g.T, cout_g.T, hout_g.p, cout_g.p]

    function dummysolve!(guess, fvec)

        if shellstream.T > tubestream.T
            shellhot_flag = true
        else
            shellhot_flag = false
        end

        Tho, Tco, Pho, Pco = guess

        if shellhot_flag
            hot = shellstream
            Pso = Pho
            cold = tubestream
            Pto = Pco
        else
            hot = tubestream
            Pto = Pho
            cold = shellstream
            Pso = Pco
        end

        Th_avg = (hot.T+Tho)/2
        Tc_avg = (cold.T+Tco)/2
        Ph_avg = (hot.p+Pho)/2
        Pc_avg = (cold.p+Pco)/2
        Twall = (Th_avg+Tc_avg)/2

        hot_avg = Stream(hot.subst_, hot.ṁ, Ph_avg, Th_avg)
        cold_avg = Stream(cold.subst_, cold.ṁ, Pc_avg, Tc_avg)

        Chot = hot.ṁ * hot_avg.Cp
        Ccold = cold.ṁ * cold_avg.Cp

        if Chot > Ccold
            Cmin = Ccold
            Cmax = Chot
        else
            Cmin = Chot
            Cmax = Ccold
        end

        Cr = Cmin/Cmax
        k = hx_.tube.therm_cond((hot.T+cold.T)/2)

        if shellhot_flag
            sflow = FlowStream(hx_.shell, hot_avg)
            tflow = FlowStream(hx_.tube, cold_avg)
        else
            sflow = FlowStream(hx_.shell, cold_avg)
            tflow = FlowStream(hx_.tube, hot_avg)
        end
        Nu_s = scorr(sflow, Twall)
        h_o = Nu_s*shellstream.k/hx_.shell.diameter
        Nu_t = tcorr(tflow, Twall)
        h_i = Nu_t*tubestream.k/hx_.tube.diameter

        Ai = π*hx_.tube.diameter*hx_.tube.length*hx_.n_tubes
        Ao = π*hx_.tube.outer_diameter*hx_.tube.length*hx_.n_tubes

        R_cond = log(hx_.tube.outer_diameter / hx_.tube.diameter) / (2π*hx_.tube.length*k*hx_.n_tubes)
        UA = 1/(1/(h_o*Ao) + 1/(h_i*Ai) + R_cond)
        NTU = UA/Cmin
        ϵ_ = ϵ(NTU, Cr, hx_.type_)
        qmax = Cmin*(hot.T - cold.T)
        # The WHOLE POINT of the calculations so far was to get q
        q = ϵ_*qmax

        if shellhot_flag
            hout = Stream(hot.subst_, hot.ṁ, hot.p - sflow.ΔP, Tho)
            cout = Stream(cold.subst_,cold.ṁ,cold.p- tflow.ΔP, Tco)
        else
            hout = Stream(hot.subst_, hot.ṁ, hot.p - tflow.ΔP, Tho)
            cout = Stream(cold.subst_,cold.ṁ,cold.p- sflow.ΔP, Tco)
        end

        # Calculate errors
        # I've chosen to use enthalpy and pressure balances on the two streams
        # as my four equations
        fvec[1] = q - abs(hot.H - hout.H)
        fvec[2] = q - abs(cold.H - cout.H)
        fvec[3] = (shellstream.p - Pso) - sflow.ΔP
        fvec[4] = (tubestream.p - Pto) - tflow.ΔP
    end
    answer = NLsolve.nlsolve(dummysolve!, x0, ftol=2.0)
    return answer
end

#import pipes
#import co2
#function test_hx_init()
#    tube = pipes.pipe(("40", 1.0), "Stainless Steel", 0.8, 0.0)
#    shell = pipes.pipe(("40", 10.0), "Stainless Steel", 0.8, 0.0)
#    testit = hx.HeatExchanger(shell, tube, 15, 0, 0, hx.st_oneshell)
#    hot_ = Stream("Sodium", 2.0, 1.0e5, 900.)
#    #print(hot_)
#    cold_ = Stream("CO₂", 1.0, 7.5e6, 400.)
#    #print(cold_)
#    print(hx_performance_solver(testit, hot_, cold_))
#end

#test_hx_init()


text:
using DataStructures;
using Makie
using AbstractPlotting, GLMakie

struct Coord
  x::Int32
  y::Int32
  z::Int32
end
struct LapCoord
  coord
  lap::UInt64
end
function data_pos(sizes::Coord, offset::Int32)::Vec3{Int32}
  coffset::Int32 = offset
  z = div(coffset,sizes.x*sizes.y)
  coffset -= z*sizes.x*sizes.y
  y = div(coffset,sizes.x)
  coffset -= y*sizes.x
  x = coffset
  return Vec3{Int32}(x,y,z)
end
function data_offset(sizes::Coord, pos::Coord)::Int32
  return pos.x + pos.y*sizes.x + pos.z*sizes.x*sizes.y + 1; #+1 due to Julia counting from 1 instead of 0
end
function debug_printQueue(queue::Queue{LapCoord})
  println("queue2 check_enqueue_vox:")
  next = iterate(queue) #reverse_iter
  while next !== nothing
    (el,state) = next
    println(el)
    next = iterate(queue,state)
  end
end
function check_enqueue_vox(sizes::Coord, data::Array{UInt8}, queue::Queue{LapCoord}, queue2::Queue{LapCoord}, vpos::Coord, lap::UInt64, step::Base.RefValue{Int64}, ok1::Bool)
  offset::Int32 = data_offset(sizes, vpos)
  if (data[offset] == 1)
    if ok1 == true
      if step[] == 200 #210 ### CHANGES THIS ###
        #println("length: ", length(queue2))
        #debug_printQueue(queue2)
        #isempty(queue2) == false
        if length(queue2) == 0
          #println("added ", LapCoord(vpos,lap), "to queue2")
          enqueue!(queue2, LapCoord(vpos,lap))
          #debug_printQueue(queue2)
        end
        step[] = 0
        data[offset] = 255
        return
      end
      step[] += 1
      data[offset] = 255
    else
      data[offset] = 100 #thread2 marks with 100
    end
    enqueue!(queue, LapCoord(vpos,lap))
  end
end
function processNeighbor(sizes::Coord, n::Int64, data::Array{UInt8}, queue::Queue{LapCoord}, queue2::Queue{LapCoord}, elem::LapCoord, step::Base.RefValue{Int64}, ok1::Bool)
  pos::Coord = elem.coord
  if n == 0
    if (pos.x+1 < sizes.x)
      vpos = Coord(pos.x+1,pos.y,pos.z)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  elseif n==1
    if (pos.x-1 >= 0)
      vpos = Coord(pos.x-1,pos.y,pos.z)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  elseif n==2
    if (pos.y+1 < sizes.y)
      vpos = Coord(pos.x,pos.y+1,pos.z)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  elseif n==3
    if (pos.y-1 >= 0)
      vpos = Coord(pos.x,pos.y-1,pos.z)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  elseif n==4
    if (pos.z+1 < sizes.z)
      vpos = Coord(pos.x,pos.y,pos.z+1)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  else # n==5
    if (pos.z-1 >= 0)
      vpos = Coord(pos.x,pos.y,pos.z-1)
      check_enqueue_vox(sizes, data, queue, queue2, vpos, elem.lap+1,step, ok1)
    end
  end
end
function srg(sizes::Coord, data::Array{UInt8}, queue1::Queue{LapCoord}, queue2::Queue{LapCoord}, lastround::UInt64)
  step::Base.RefValue{Int64} = Ref(1)
  ok1::Bool = true
  ok2::Bool = true
  started::Bool = false
  startedtmp::Bool = false
  _coord = Coord(0,0,0)
  elem1 = LapCoord(_coord, 0)
  elem2 = LapCoord(_coord, 0)
  while ( ok1 == true || ok2 == true )
    if (length(queue1) != 0 && first(queue1).lap <= lastround)
      elem1 = dequeue!(queue1)
      ok1 = true
    else
      ok1 = false
    end
    #println("length: ", length(queue2))
    if (length(queue2) != 0 && first(queue2).lap <= lastround)
      elem2 = dequeue!(queue2)
      #println(elem2) # elem2
      ok2 = true
    else
      ok2 = false
    end
    for n in 0:5 #n for 6 neighbors 0..5
      if ok1 == true
        processNeighbor(sizes, n, data, queue1, queue2, elem1, step, true) #called by thread1
      end
      if ok2 == true
        processNeighbor(sizes, n, data, queue2, queue2, elem2, step, false) #called by thread2
      end
    end
    if startedtmp == true
      started = true
      step[] = 0
    end
    #println("lap:",lap)
  end
end
function computePlotSizes(size::Coord)::Coord
  return Coord(2*size.x + 1, 2*size.y + 1, 2*size.z + 1)
end
function transform_into_plotcoord(pos::Vec3{Int32})::Vec3{Int32}
  #return Vec3{Int32}(pos[1]*2+1,pos[2]*2+1,pos[3]*2+1)
  return Vec3{Int32}(pos[1]+1,pos[2]+1,pos[3]+1)
end
function Cube_surface_edges(vecmin::Vec3{Int32}, vecdist::Vec3{Int32}, outx::Array{Int32,1}, outy::Array{Int32,1}, outz::Array{Int32,1})
  vecmax::Vec3 = vecmin + vecdist
  #x000-x001-x011-x010 \\
  #x000-x100-x101-x001 \\
  #x000-x010-x110-x111 \\
  #x011-x111-x101-x100 \\
  #x110
  #
  #x000
  append!(outx, vecmin[1])
  append!(outy, vecmin[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmin[1])
  append!(outy, vecmin[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmin[1])
  append!(outy, vecmax[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmin[1])
  append!(outy, vecmax[2])
  append!(outz, vecmin[3]) #
  # x000
  append!(outx, vecmin[1])
  append!(outy, vecmin[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmin[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmin[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmin[1])
  append!(outy, vecmin[2])
  append!(outz, vecmax[3]) #
  # x000
  append!(outx, vecmin[1])
  append!(outy, vecmin[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmin[1])
  append!(outy, vecmax[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmax[2])
  append!(outz, vecmin[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmax[2])
  append!(outz, vecmax[3]) #
  # x011
  append!(outx, vecmin[1])
  append!(outy, vecmax[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmax[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmin[2])
  append!(outz, vecmax[3]) #
  append!(outx, vecmax[1])
  append!(outy, vecmin[2])
  append!(outz, vecmin[3]) #
  # x110
  append!(outx, vecmax[1])
  append!(outy, vecmax[2])
  append!(outz, vecmin[3]) #
end
function plotVoxel(posmin::Vec3{Int32}, starting::Bool)
  xout::Array{Int32,1} = []
  yout::Array{Int32,1} = []
  zout::Array{Int32,1} = []
  #posmin = Vec3(0,0,0)
  posdist::Vec3{Int32} = Vec3{Int32}(1,1,1)
  Cube_surface_edges(posmin, posdist, xout,yout,zout)
  rect = Rect3D{Int32}(posmin, posdist)
  if (starting == true)
    mesh!(rect, color = (:red, 0.1), transparency = false) #true
  else
    mesh!(rect, color = (:blue, 0.1), transparency = false) #true
  end
  lines!(xout, yout, zout, color = :black)
end
function main()
  GLMakie.activate!()
  AbstractPlotting.inline!(false)
  lrs::Array{UInt64} = [11] #[0,1,2,3,4,5,6,7,8,9,10,11]
                        #,36,40,44,48]
  cacheline_size::Int32 = 8 # simplify 64 Byte cache line
  #sizehint!(s, n) #???
  for i in lrs
    scene = Scene()
    queue = Queue{LapCoord}()
    queue2 = Queue{LapCoord}()
    seed = Coord(50,50,50)
    seedvec = Vec3{Int32}(seed.x,seed.y,seed.z)
    seedPlot = transform_into_plotcoord(seedvec);
    plotVoxel(seedPlot, true)
    sizes = Coord(100,100,100)
    data = ones(UInt8, 10^6) # 1-dim arary of size 10^6
    lap = 0
    pos = Vec3{Int32}(0,0,0)
    enqueue!(queue, LapCoord(seed,lap))
    data[data_offset(sizes,seed)] = 255
    lastround::UInt64 = 0
    lastround = i
    srg(sizes,data,queue,queue2,lastround)
    for j in 1:length(data)
      if data[j] == 255
        pos = data_pos(sizes, convert(Int32,j - 1))
        posPlot = transform_into_plotcoord(pos)
        plotVoxel(posPlot,false) #plot thread0
      end
      if data[j] == 100
        pos = data_pos(sizes, convert(Int32,j - 1))
        posPlot = transform_into_plotcoord(pos)
        plotVoxel(posPlot,true) #plot thread1
      end
    end

    #while (length(queue) != 0)
    #  elem = dequeue!(queue)
    #  pos::Coord = elem.coord
    #  posvec = Vec3{Int32}(pos.x,pos.y,pos.z)
    #  posPlot = transform_into_plotcoord(posvec);
    #  plotVoxel(posPlot, false) #can not plot multiple different colors or surface not visible
    #end
    #println("length: ", length(queue2))
    #while (length(queue2) != 0)
    #  elem = dequeue!(queue2)
    #  pos::Coord = elem.coord
    #  posvec = Vec3{Int32}(pos.x,pos.y,pos.z)
    #  posPlot = transform_into_plotcoord(posvec);
    #  plotVoxel(posPlot, true)
    #end
    file = string("srg2threads",string(lastround), ".png")
    save(file, scene)
  end
end
main()

##TODO
# TESTING display(scene)
#  computation 64 * length(queue) + length(queue) * 3 * 32B
#  length(queue) * 160B


text:
struct RBFDictionary{S,T,R <: RBF{T}} <: Dictionary{S,T}
    centers ::  Vector{S}
    rbf     ::  R
    domain  ::  Domain{S}
end

RBFDictionary(centers::Vector{S}, rbf::RBF) where {S} = RBFDictionary(centers, rbf, DomainSets.FullSpace{S}())

size(dict::RBFDictionary) = size(dict.centers)

domain(dict::RBFDictionary) = dict.domain

name(dict::RBFDictionary) = "A list of radial basis functions ($(typeof(dict.rbf)))"

center(dict::RBFDictionary, idx) = dict.centers[idx]

support(dict::RBFDictionary) = dict.domain

unsafe_eval_element(dict::RBFDictionary, idx, x) =
    dict.rbf(norm(x - center(dict,idx)))

plotgrid(dict::RBFDictionary, n) = equispaced_grid(boundingbox(dict.domain), n)

import CompactTranslatesDict:
    CompactTranslationDict,
    eval_kernel,
    kernel_span

struct PeriodicRBFDictionary{T,R <: RBF{T}} <: CompactTranslatesDict.CompactTranslationDict{T}
    rbf     ::  R
    n       ::  Int
end

eval_kernel(dict::PeriodicRBFDictionary, x) = dict.rbf(x)

function kernel_span(dict::PeriodicRBFDictionary)
    D = approximate_support(dict.rbf)
    -D..D
end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
using LinearAlgebra
using RecursiveArrayTools
using Distributions
using Distances
using Random
using Profile
using ProfileView
using JLD
using MAT
using IterativeSolvers
using BlockArrays
Random.seed!(123)

# Creating uniform random data
NTrain = 1000000
d =  2
NTest = 50
N = NTest + NTrain
x = rand(d, N)
setSplit = randperm(N)
trainInds = setSplit[1:NTrain]
testInds = setSplit[(NTrain + 1) : end]
xTrain = x[ :, trainInds ]
xTest = x[:, testInds ]
x = hcat(xTrain, xTest)
y = randn(N)
y = PseudoBlockArray(y, [NTest, NTrain])

ρ = 3.0
λ = 1.5

# reordering points and forming skeletons
# @profview colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
@time colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
x = x[:, P]

LTest = SparseMatrixCSC(N, N, colptr, rowval, ones(length(rowval)))
LTest22 = LTest[(NTest + 1) : end, (NTest + 1) : end]

xTest = x[:, 1:NTest]
xTrain = x[:, (NTest + 1) : end]
skeletons = construct_skeletons(colptr, rowval, distances, λ)

l = 0.1
# σ = vcat(zeros(NTest), fill(10.0, NTrain))
σ = vcat(zeros(NTest), 1.0 * exp.(randn(NTrain)) .+ 1e-4)
σ = PseudoBlockArray(σ, [NTest, NTrain])

function covfunc!(r)
  matern32!(r,l)
end

LinearAlgebra.BLAS.set_num_threads(1)
# @profview L = assembleL(skeletons, x, covfunc!)
@time L = assembleL(skeletons, x, covfunc!)
# @time LNaive = assembleL(skeletons, x, covfunc!, σ)
# @time L22 = L[(NTest + 1) : end, (NTest + 1) : end]
# @time A22 = L22 * L22' .+ spdiagm(0 => 1 ./ max.(σ[Block(2)], eps(eltype(σ))))
# # @time U22 = sparse(triu(A22))
# @time U22 = squareSparse(L22)+ spdiagm(0 => 1 ./ max.(σ[Block(2)], 1e-30))
# @time icholU_high_level!(U22)
#
#
# # The vecchia posterior mean
# z = (L22 * (L22' * y[Block(2)]))
# z = cg(A22 , z, Pl = IChol(U22',U22), verbose=true, maxiter=5, tol=1e-16)
# z = (Diagonal(σ[Block(2)]))\ z
# μVechia = (L' \ (L \ (vcat(zero(y[Block(1)]), z))))[1:NTest]
# # display(μVechia)
# @show mean(abs.(μVechia))
# @show median(abs.(μVechia))
#
#
# # The naive application of vecchia
# μNaive = -LNaive[1:NTest, 1:NTest]' \ (LNaive[(NTest + 1) : end, 1 : NTest]' * y[Block(2)])
# # display(μNaive)
# @show mean(abs.(μNaive))
# @show median(abs.(μNaive))
#


text:
"""
    SeisAddNoise(d, snr; <keyword arguments>)

Add noise to an N-dimensional input
data `d`. Should specify the signal-to-noise ratio level `snr`.
Noise can be band limited using kewyord `L`.

# Arguments
- `d::Array{Real, N}`: N-dimensional data.
- `snr::Real`: signal-to-noise ratio.
- `db::Bool=false`: Flag is false if snr is given by amplitude. Flag is true if
snr is given in dB.
- `pdf::AbstractString="gaussian"`: random noise probability distribution:
`"gaussian"` or `"uniform"`.
- `L::Int=1`: averaging operator length to band-limit the random noise.

# Examples
```julia
julia> using PyPlot
julia> w = Ricker(); wn = SeisAddNoise(w, 2); plot(w); plot(wn);
MeasureSNR(w, wn)
julia> d = SeisHypEvents(); dn = SeisAddNoise(d, 1.0, db=true, L=9);
SeisPlotTX([d dn]); MeasureSNR(d, dn, db=true)
```
Credits: Juan I. Sabbione, 2016
"""
function SeisAddNoise(d::Array{T, N}, snr::Real; db::Bool=false,
                      pdf::AbstractString="gaussian", L::Int=1) where {T <: Real, N}

    noise = GenNoise(size(d), pdf, L=L)
    if db==false
        noise = noise/norm(noise) * norm(d)/snr
    elseif db==true
        noise = noise/norm(noise) * norm(d)/10.0^(0.05*snr)
    end
    noisy = d + noise
    @assert(abs(MeasureSNR(d, noisy, db=db) - snr) < 1e10*eps(AbstractFloat(snr)))
    return noisy

end

# Generate the noise
function GenNoise(dims::Tuple, pdf::AbstractString; L::Int=1)
    N = length(dims)
    n1 = dims[1]
    N == 1 ? nx = 1 : nx = *(dims[2:end]...)
    hamm = Hamming(L)
    L2 = floor(Int,L/2)
    if pdf=="gaussian"
        noise = randn(n1, nx)
    elseif pdf=="uniform"
        noise = -1 + 2*rand(n1, nx)
    else
        error("pdf must be gussian or uniform")
    end
    for j = 1:nx
         noise[:,j] = conv(noise[:,j], hamm)[L2+1:n1+L2]
    end
    reshape(noise, dims...)
end


text:
module DalyNotebook

using Plots, DataFrames
output = ""

macro md_str(s, t...)
    global output *= s * "\n"
end

function save_plot(plot; file::String, desc::String)
    mkpath(".fig")
    savefig(plot, ".fig/$(file)")
    plot_img_str = "\n![$desc](.fig/$(file))\n"
    global output *= plot_img_str
end

md"""
# DALY - economical measurement of the burden of disease
And some contents
"""

rectangle(x, y, w, h) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

function plotDALY(;
    skrocenie_zycia,
    czas_choroby,
    obnizenie_jakosci,
    plot_choroba = true,
    plot_pyll = false,
    plot_yld = false,
    legend = :topright,
    srednia_dlugosc_zycia = 78.58,
    labels = true
)
    wiek_smierci = srednia_dlugosc_zycia - skrocenie_zycia
    wiek_zachorowania = wiek_smierci - czas_choroby
    jakoscZyciaZdrowego = DataFrame(
        :Year => [0, srednia_dlugosc_zycia, srednia_dlugosc_zycia * 1.1],
        :Quality => [1.0, 0.0, 0.0],
    )
    jakoscZyciaChorego = DataFrame(
        :Year => [0, wiek_zachorowania, wiek_smierci],
        :Quality => [1.0, 1.0 - obnizenie_jakosci, 0.0],
    )

    p = plot(
        jakoscZyciaZdrowego.Year, jakoscZyciaZdrowego.Quality,
        label = labels ? "Życie zdrowego człowieka" : "",
        xlabel = labels ? "Lata życia" : "", ylabel = labels ? "Jakość życia [%]" : "",
        seriestype = :steppost, linewidth = 3, linecolor = :green,
        fill = (0, 0.5, :green), legend = legend,
    )
    if plot_choroba
        plot!(p,
            jakoscZyciaChorego.Year, jakoscZyciaChorego.Quality,
            label = "Życie chorego człowieka",
            seriestype = :steppost, linewidth = 3, linecolor = :red,
            fill = (0, 0.99, :darkred),
        )

        plot!(p,
            [25, jakoscZyciaChorego.Year[2]], [0.9, jakoscZyciaChorego.Quality[2]], arrow = (:closed, 2.0), color = :white, label = nothing,
            ann = [(25, 0.9, text("Zachorowanie", 8, :white))]
        )
        plot!(p,
            [62, wiek_smierci], [0.2, 0.0], arrow = (:closed, 2.0), color = :white, label = nothing,
            ann = [(62, 0.2, text("Śmierć", 8, :white))]
        )
    end

    if plot_pyll
        plot!(p,
            rectangle(wiek_smierci, 0, skrocenie_zycia, 1.0),
            color = :black, opacity = 0.9, label = "PYLL",
            ann = labels ? [(wiek_smierci + skrocenie_zycia / 2, 0.5, text("PYLL", 8, :white))] : []
        )
    end

    if plot_yld
        plot!(p,
            rectangle(wiek_zachorowania, 1.0 - obnizenie_jakosci, czas_choroby, obnizenie_jakosci),
            color = :black, opacity = 0.8, label = "YLD",
            ann = labels ? [(wiek_zachorowania + czas_choroby / 2, 1.0 - obnizenie_jakosci / 2, text("YLD", 8, :white))] : []
        )
    end

    p
end

plot_daly_choroba = plotDALY(
    skrocenie_zycia = 10,
    obnizenie_jakosci = 0.2,
    czas_choroby = 12,
)

md"""
## Subsection
"""
save_plot(plot_daly_choroba, file = "daly-choroba.svg", desc = "A figure")

open(f -> write(f, output), "index.md", "w+")
end # module

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
@testset "$TEST $G" begin

n = 10
m = n*(n-1)/2
c = ones(Int, n)
g = CompleteGraph(n, G)
@test  modularity(g, c) == 0
#
g = G(n)
@test modularity(g, c) == 0

end # testset


text:
# Filename: install.jl
# Author: Dharmesh Tarapore <dharmesh@bu.edu>

Pkg.add("Compat")
Pkg.add("Distributions")
Pkg.add("JLD")
Pkg.add("HDF5")
Pkg.add("PyPlot")
Pkg.add("ArgParse")

Pkg.clone("https://github.com/weirdindiankid/CASInterface.jl", "CASInterface")
Pkg.clone("https://github.com/weirdindiankid/RLESUtils.jl", "RLESUtils")
Pkg.clone("https://github.com/weirdindiankid/TensorFlow.jl", "TensorFlow")

Pkg.init(); run(`ln -s $(pwd()) $(Pkg.dir("SISLES"))`); Pkg.pin("SISLES"); Pkg.resolve()
using SISLES; @assert isdefined(:SISLES); @assert typeof(SISLES) === Module


text:
using Test, Random, LinearAlgebra

using MAGMA: magmaInit, magmaFinalize, gesvd!

@testset "magmaSgesvd" begin
    include("magmaSgesvd.jl")
end

@testset "magmaDgesvd" begin
    include("magmaDgesvd.jl")
end


text:
module Utils
  export posinf, neginf, hasposinf, hasneginf, hasinf, isposinf, isneginf
  """
      posinf(T)

  Positive infinity of type `T`, or `nothing` if this is not possible.
  """
  posinf(::Type{T}) where {T<:Real} = try convert(T, Inf); catch; nothing; end

  """
      neginf(T)

  Negative infinity of type `T`, or `nothing` if this is not possible.
  """
  neginf(::Type{T}) where {T<:Real} = try convert(T, -Inf); catch; nothing; end

  """
      hasposinf(T)

  True if positive infinity is representable in type `T`.
  """
  hasposinf(::Type{T}) where {T<:Real} = posinf(T) !== nothing

  """
      hasneginf(T)

  True if negative infinity is representable in type `T`.
  """
  hasneginf(::Type{T}) where {T<:Real} = neginf(T) !== nothing

  """
      hasinf(T)

  True if positive and negative infinity are representable in type `T`.
  """
  hasinf(::Type{T}) where {T<:Real} = hasposinf(T) && hasneginf(T)

  """
      isposinf(x)

  True if `x` is positive infinity.
  """
  isposinf(x::Real) = isinf(x) && !signbit(x)

  """
      isneginf(x)

  True if `x` is negative infinity.
  """
  isneginf(x::Real) = isinf(x) && signbit(x)
end

using .Utils


text:
module SurrogateCube
import Distributions: Normal, Laplace, Cauchy, Categorical
using Docile
#
#
#
# We want to model a disturbed dataset of the form X_xyt = μ + n + M
# where μ is some constant, periodic, or quasiperiodic mean, n is some kind
# of Noise and M is the disturbance

# Define the return Datastructure of genDataset, so that it can be easily saved
type DataCube
	variables::Array{Float64,4}
	components::Array{Float64,4}
	eventmap::Array{Float64,4}
	baselines::Vector
	compnoise::Vector
	events::Vector
	varnoise::Vector
	coupling::Vector
	weight::Vector
	k1::Float64
	k2::Float64
	k3::Float64
	k4::Float64
end
export DataCube
include("noise.jl")
include("baseline.jl")
include("event.jl")
include("io.jl")
include("datafarm.jl")
include("coupling.jl")
include("variableweights.jl")

"""
Function that generates a single datacube, this can serve as a basis to build an independent component of a dataset. The input is as follows:
- mt: Baseline, the process mean
- nt: Noise, the noise chosen
- dt: Event, the disturbance event
- Nlon: Number of longitudes
- Nlat: Number of Latitudes
- Ntime: Number of Time steps
- `kb` : Strength of baseline modulation by event, kb=0 means no modulation, kb=1 means double amplitude, kb=-1 means half amplitude
- `kn` : Strength of noise modulation by event, kn=0, means no modulation, kn=1 double noise, kn=-1 half noise
- `ks`: Strength of mean shift event, ks=0 no shift ks=-1 shift in negative direction by one sd(noise), ks=1 shift by one sd(noise) in positive direction
"""
function genDataStream(mt::Baseline,nt::Noise,dt::Event,Ntime,Nlat,Nlon,kb,kn,ks)
    b  = genBaseline(mt,Ntime,Nlat,Nlon)
    n  = genNoise(nt,Ntime,Nlat,Nlon)
    ev = genEvent(dt,Ntime,Nlat,Nlon)
    x = b .* 2.^(kb*ev) + n .* 2.^(kn*ev) + std(n) * ks * ev
    outevent!(dt,ev)
    return(x,ev)
end
#If only one k is given, apply it to k3
genDataStream(mt::Baseline,nt::Noise,dt::Event,Ntime,Nlat,Nlon,k)=genDataStream(mt,nt,dt,Ntime,Nlat,Nlon,0.0,0.0,k)
export genDataStream

"""
This function generates a whole multivariate datacube, the input is as follows. It needs the followings input:
- mt: Baseline or Vector{Baseline}, the process mean for the components
- nt: Noise or Vector{Noise}, the noise chosen for the independent components
- dt: Event or Vector{Event}, the distirbance events for each component
- dataNoise: Noise or Vector{Noise}, the noise added to each variable of the DataCube
- Ncomp: Number of independent components to generate
- Nvar: Number of variables to generate
- Nlon: Number of longitudes
- Nlat: Number of Latitudes
- Ntime: Number of Time steps
- `kb` : Strength of baseline modulation by event, kb=0 means no modulation, kb=1 means double amplitude, kb=-1 means half amplitude
- `kn` : Strength of noise modulation by event, kn=0, means no modulation, kn=1 double noise, kn=-1 half noise
- `ks` : Strength of mean shift event, ks=0 no shift ks=-1 shift in negative direction by one sd(noise), ks=1 shift by one sd(noise) in positive direction
"""
function genDataCube{T<:Baseline,U<:Noise,V<:Event,W<:Noise}(
	mt::Union{Vector{T},T},
	nt::Union{Vector{U},U},
	dt::Union{Vector{V},V},
	dataNoise::Union{Vector{W},W},
	coupling::Vector,
	weight::Vector,
	Ncomp,
	Nvar,
	Ntime,
	Nlat,
	Nlon,
	kb::Number,
	kn::Number,
	ks::Number,
	kw::Number)
	#Make vectors if single values were given by the user
	mt2 = isa(mt,Vector) ? mt : fill(mt,Ncomp)
	nt2 = isa(nt,Vector) ? nt : fill(nt,Ncomp)
	dt2 = isa(dt,Vector) ? dt : fill(dt,Ncomp)
  dataNoise2 = isa(dataNoise,Vector) ? dataNoise : fill(dataNoise,Nvar)

	# First generate the independent components of the datacube
	acomp = zeros(Float64,Ntime,Nlat,Nlon,Ncomp)
	dcomp = zeros(Float64,Ntime,Nlat,Nlon,Ncomp)

	for i = 1:Ncomp
		x,d = genDataStream(mt2[i],nt2[i],dt2[i],Ntime,Nlat,Nlon,kb,kn,ks)
		acomp[:,:,:,i]=x
		dcomp[:,:,:,i]=d
	end

	# Generate weights to generate variables out of the independent components
	xout = zeros(Float64,Ntime,Nlat,Nlon,Nvar)
  combineComponents!(xout,acomp,dcomp,coupling,dataNoise2,weight,kw)

	DataCube(xout,acomp,dcomp,mt2,nt2,dt2,dataNoise2,coupling,weight,kb,kn,ks,kw)
end
genDataCube{T<:Baseline,U<:Noise,V<:Event,W<:Noise}(mt::Union{Vector{T},T},nt::Union{Vector{U},U},dt::Union{Vector{V},V},dataNoise::Union{Vector{W},W},Ncomp,Nvar,Ntime,Nlat,Nlon,k)=genDataCube(mt,nt,dt,dataNoise,Ncomp,Nvar,Ntime,Nlat,Nlon,0.0,0.0,k)
export genDataCube

"""
This function generates a series of whole multivariate datacubes for a vector of k-values given, the input is as follows. It needs the followings input:
- `mt`: Baseline or Vector{Baseline}, the process mean for the components
- `nt`: Noise or Vector{Noise}, the noise chosen for the independent components
- `dt`: Event or Vector{Event}, the distirbance events for each component
- `dataNoise`: Noise or Vector{Noise}, the noise added to each variable of the DataCube
- `Ncomp`: Number of independent components to generate
- `Nvar`: Number of variables to generate
- `Nlon`: Number of longitudes
- `Nlat`: Number of Latitudes
- `Ntime`: Number of Time steps
- `kb` : Strength of baseline modulation by event, kb=0 means no modulation, kb=1 means double amplitude, kb=-1 means half amplitude, can be a vector
- `kn` : Strength of noise modulation by event, kn=0, means no modulation, kn=1 double noise, kn=-1 half noise, can be a vector
- `ks` : Strength of mean shift event, ks=0 no shift ks=-1 shift in negative direction by one sd(noise), ks=1 shift by one sd(noise) in positive direction, can be a vector
-
"""
function genDataCube{T<:Baseline,U<:Noise,V<:Event,W<:Noise}(mt::Union{Vector{T},T},nt::Union{Vector{U},U},dt::Union{Vector{V},V},dataNoise::Union{Vector{W},W},coupling,weight,Ncomp,Nvar,Ntime,Nlat,Nlon,kb,kn,ks,kw,experimentname,path = joinpath(Pkg.dir(),"SurrogateCube","data"));

    isdir(joinpath(path,experimentname)) ? error("Directory $path already exists.") : mkdir(joinpath(path,experimentname))

    for ikb = 1:length(kb), ikn=1:length(kn), iks=1:length(ks), ikw = 1:length(kw)

        dc = genDataCube(mt,nt,dt,dataNoise,coupling,weight,Ncomp,Nvar,Ntime,Nlat,Nlon,kb[ikb],kn[ikn],ks[iks],kw[ikw])

        save_cube(dc,joinpath(path,experimentname),string("cube_kb",kb[ikb],"kn_",kn[ikn],"_ks",ks[iks],"_kw",kw[ikw],".nc"))
    end
end

end  #Module

include("experiments/small2dfarm.jl")

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

using CSV

using DataFrames

using Plots

using StatsPlots

include("./functions/rectangle.jl")

#using Pyplot
#pyplot()

#Encoder's resolutio
Encoder_mm_per_Pulse = 0.02

#experiment number
experiment = "03"
speed = "3"

#read CSV of encoder logger file
Encoder = CSV.read(
    string("../measures/xp_pulse_all_seuil/encoder/20201119_gel_1cm/Logger_encoder_gel_1cm_v",speed,"_", experiment, ".txt");
    header = 1,
)

#convert the number of pulse position change into mm
Encoder.PositionChange_mm = Encoder." PositionChange" * Encoder_mm_per_Pulse

#convert the detected time of change into from ms to s
#WARNING sometime the encoder is disconnected and we lose measures
#you can track those deconnection by looking at TimeChange==0ms
#As we don't know the time duration of deconnection, I added 1s delay
time = cumsum(Encoder." TimeChange" + (Encoder." TimeChange" .== 0) .* 1000) / 1000

#plot the encoder velocity in time
plot_Encoder = plot(
    time,
    Encoder.PositionChange_mm ./ Encoder." TimeChange",
    xlabel = "time[s]",
    ylabel = "Velocity[m/s]",
    label = "Velocity",
    title = string("velocity and position of experiment ", experiment),
    xlims = (42.5,43.5)
)

#draws 1s grey rectangle during unknown deconnection duration
for k in findall(x -> x == 0, Encoder." TimeChange")
    if k == 1
        plot_Encoder = plot!(
            rectangle(time[k+1] - time[k], 2, time[k], 0),
            color = :grey,
            label = "unknown timing",
        )
    elseif k != size(time)
        if k == findfirst(x -> x == 0, Encoder." TimeChange")
            plot_Encoder = plot!(
                rectangle(time[k+1] - time[k-1], 2, time[k-1], 0),
                color = :grey,
                label = "unknown timing",
            )
        elseif k != findlast(x -> x == 0, Encoder." TimeChange")
            plot_Encoder = plot!(
                rectangle(time[k+1] - time[k-1], 2, time[k-1], 0),
                color = :grey,
                label = false,
            )
        end
    end
end

# #add driven distance in front of the velocity
# #add right y axis
# #WARNING the right axis is only visible in the pdf image, not in the Juno Plots interface
# plot_Encoder=plot!(twinx(),
#     time,
#     cumsum(Encoder.PositionChange_mm./1000),
#     color=:red,
#     ylabel = "Position[m]",
#     label = "Position")
plot_Encoder
#save the plot
#savefig(plot_Encoder, string("save_figure/velocity/velocity_", experiment, ".pdf"))


text:
abstract type AbstractTableauSplitting{T <: Real} <: AbstractTableau{T} end

"Tableau for non-symmetric splitting methods."
struct TableauSplittingNS{T} <: AbstractTableauSplitting{T}
    @HeaderTableau

    a::Vector{T}
    b::Vector{T}

    function TableauSplittingNS{T}(name, o, s, a, b) where {T}
        @assert s == length(a) == length(b)
        new(name, o, s, a, b)
    end
end

function TableauSplittingNS(name, o, a::Vector{T}, b::Vector{T}) where {T}
    TableauSplittingNS{T}(name, o, length(a), a, b)
end

"Tableau for symmetric splitting methods with general stages."
struct TableauSplittingGS{T} <: AbstractTableauSplitting{T}
    @HeaderTableau

    a::Vector{T}
    b::Vector{T}

    function TableauSplittingGS{T}(name, o, s, a, b) where {T}
        @assert s == length(a) == length(b)
        new(name, o, s, a, b)
    end
end

function TableauSplittingGS(name, o, a::Vector{T}, b::Vector{T}) where {T}
    TableauSplittingGS{T}(name, o, length(a), a, b)
end

"Tableau for symmetric splitting methods with symmetric stages."
struct TableauSplittingSS{T} <: AbstractTableauSplitting{T}
    @HeaderTableau

    a::Vector{T}

    function TableauSplittingSS{T}(name, o, s, a) where {T}
        @assert s == length(a)
        new(name, o, s, a)
    end
end

function TableauSplittingSS(name, o, a::Vector{T}) where {T}
    TableauSplittingSS{T}(name, o, length(a), a)
end

"Splitting integrator."
struct IntegratorSplitting{DT, TT, VT, ST <: AbstractTableauSplitting, FT, CT, N} <: DeterministicIntegrator{DT,TT}
    equation::SODE{DT,TT,VT,N}
    tableau::ST
    f::FT
    c::CT
    Δt::TT

    q::Array{DT,1}
    qₑᵣᵣ::Array{DT,1}
    v::Array{DT,1}

    function IntegratorSplitting{DT,TT,VT,ST,FT,CT,N}(equation::SODE, tableau::ST, f::FT, c::CT, Δt) where {DT,TT,VT,ST,FT,CT,N}
        @assert length(f) == length(c)
        D = equation.d
        new(equation, tableau, f, c, Δt, zeros(DT,D), zeros(DT,D), zeros(DT,D))
    end
end

function IntegratorSplitting(equation::SODE{DT,TT,VT,N}, tableau::ST, f::Vector{Int}, c::Vector{TT}, Δt::TT) where {DT, TT, VT, ST <: AbstractTableauSplitting, N}
    ft = Tuple(f)
    ct = Tuple(c)
    IntegratorSplitting{DT,TT,VT,ST,typeof(ft),typeof(ct),N}(equation, tableau, ft, ct, Δt)
end

function get_splitting_coefficients(r, a::Vector{TT}, b::Vector{TT}) where {TT}
    @assert length(a) == length(b)

    s = length(a)
    f = zeros(Int, 2r*s)
    c = zeros(TT,  2r*s)

    for i in 1:s
        for j in 1:r
            f[(2i-2)*r+j] = j
            c[(2i-2)*r+j] = a[i]
        end
        for j in 1:r
            f[(2i-1)*r+j] = r-j+1
            c[(2i-1)*r+j] = b[i]
        end
    end

    return f, c
end

"Construct splitting integrator for non-symmetric splitting tableau with general stages."
function IntegratorSplitting(equation::SODE{DT,TT,VT}, tableau::ST, Δt::TT) where {DT, TT, VT, ST <: TableauSplittingNS{TT}}

    # basic method: Lie composition
    # \varphi_{\tau,A} = \varphi_{\tau,v_1} \circ \varphi_{\tau,v_2} \circ \hdots \varphi_{\tau,v_{r-1}} \circ \varphi_{\tau,v_r}
    # \varphi_{\tau,B} = \varphi_{\tau,v_r} \circ \varphi_{\tau,v_{r-1}} \circ \hdots \varphi_{\tau,v_2} \circ \varphi_{\tau,v_1}

    # integrator:
    # \varphi_{NS} = \varphi_{b_s \tau, B} \circ \varphi_{a_s \tau, A} \circ \hdots \circ \varphi_{b_1 \tau, B} \circ \varphi_{a_1 \tau, A}

    f, c = get_splitting_coefficients(length(equation.v), tableau.a, tableau.b)

    # R = length(equation.v)
    # S = tableau.s
    #
    # f = zeros(Int, 2R*S)
    # c = zeros(TT,  2R*S)
    #
    # for i in 1:S
    #     for j in 1:R
    #         f[(2i-2)*R+j] = j
    #         c[(2i-2)*R+j] = tableau.a[i]
    #     end
    #     for j in 1:R
    #         f[(2i-1)*R+j] = R-j+1
    #         c[(2i-1)*R+j] = tableau.b[i]
    #     end
    # end

    IntegratorSplitting(equation, tableau, f, c, Δt)
end

"Construct splitting integrator for symmetric splitting tableau with general stages."
function IntegratorSplitting(equation::SODE{DT,TT,VT}, tableau::ST, Δt::TT) where {DT, TT, VT, ST <: TableauSplittingGS{TT}}

    # basic method: Lie composition
    # \varphi_{\tau,A} = \varphi_{\tau,v_1} \circ \varphi_{\tau,v_2} \circ \hdots \varphi_{\tau,v_{r-1}} \circ \varphi_{\tau,v_r}
    # \varphi_{\tau,B} = \varphi_{\tau,v_r} \circ \varphi_{\tau,v_{r-1}} \circ \hdots \varphi_{\tau,v_2} \circ \varphi_{\tau,v_1}

    # integrator:
    # \varphi_{GS} = \varphi_{a_1 \tau, A} \circ \varphi_{b_1 \tau, B} \circ \hdots \circ \varphi_{b_1 \tau, B} \circ \varphi_{a_1 \tau, A}

    f, c = get_splitting_coefficients(length(equation.v), tableau.a, tableau.b)

    # R = length(equation.v)
    # S = tableau.s
    #
    # f = zeros(Int, 2R*S)
    # c = zeros(TT,  2R*S)
    #
    # for i in 1:S
    #     for j in 1:R
    #         f[(2i-2)*R+j] = j
    #         c[(2i-2)*R+j] = tableau.a[i]
    #     end
    #     for j in R:-1:1
    #         f[(2i-1)*R+j] = R-j+1
    #         c[(2i-1)*R+j] = tableau.b[i]
    #     end
    # end

    IntegratorSplitting(equation, tableau, vcat(f, f[end:-1:1]), vcat(c, c[end:-1:1]), Δt)
end

"Construct splitting integrator for symmetric splitting tableau with symmetric stages."
function IntegratorSplitting(equation::SODE{DT,TT,VT}, tableau::ST, Δt::TT) where {DT, TT, VT, ST <: TableauSplittingSS{TT}}

    # basic method: symmetric Strang composition
    # \varphi_{\tau,A} = \varphi_{\tau/2,v_1} \circ \varphi_{\tau/2,v_2} \circ \hdots \varphi_{\tau/2,v_{r-1}} \circ \varphi_{\tau/2,v_r}
    #              \circ \varphi_{\tau/2,v_r} \circ \varphi_{\tau/2,v_{r-1}} \circ \hdots \varphi_{\tau/2,v_2} \circ \varphi_{\tau/2,v_1}

    # integrator:
    # \varphi_{SS} = \varphi_{a_1 \tau, A} \circ \varphi_{a_2 \tau, A} \hdots \circ \varphi_{a_s \tau, A} \circ \hdots \circ \varphi_{a_2 \tau, A} \circ \varphi_{a_1 \tau, A}

    r = length(equation.v)
    a = vcat(tableau.a, tableau.a[end-1:-1:1]) ./ 2
    s = length(a)

    f = zeros(Int, 2r*s)
    c = zeros(TT,  2r*s)

    for i in 1:s
        for j in 1:r
            f[(2i-2)*r+j] = j
            c[(2i-2)*r+j] = a[i]
            f[(2i-1)*r+j] = r-j+1
            c[(2i-1)*r+j] = a[i]
        end
    end

    IntegratorSplitting(equation, tableau, f, c, Δt)
end

equation(int::IntegratorSplitting) = int.equation
timestep(int::IntegratorSplitting) = int.Δt

"Explicit Runge-Kutta integrator cache."
mutable struct IntegratorCacheSplitting{DT,TT,D} <: ODEIntegratorCache{DT,D}
    n::Int
    t::TT
    t̅::TT
    q::Vector{TwicePrecision{DT}}
    q̅::Vector{TwicePrecision{DT}}
    v::Vector{DT}
    s̃::Vector{DT}

    function IntegratorCacheSplitting{DT,TT,D}() where {DT,TT,D}
        q = zeros(TwicePrecision{DT}, D)
        q̅ = zeros(TwicePrecision{DT}, D)
        v = zeros(DT, D)
        s̃ = zeros(DT, D)
        new(0, zero(TT), zero(TT), q, q̅, v, s̃)
    end
end

function create_integrator_cache(int::IntegratorSplitting{DT,TT}) where {DT,TT}
    IntegratorCacheSplitting{DT, TT, ndims(equation(int))}()
end

function CommonFunctions.reset!(cache::IntegratorCacheSplitting{DT,TT}, Δt::TT) where {DT,TT}
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.set_solution!(cache::IntegratorCacheSplitting, sol, n=0)
    t, q = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
end

"Integrate ODE with splitting integrator."
function integrate_step!(int::IntegratorSplitting{DT,TT,FT}, cache::IntegratorCacheSplitting{DT,TT}) where {DT,TT,FT}
    local tᵢ::TT

    # reset cache
    reset!(cache, int.Δt)

    # compute internal stages
    for i in eachindex(int.f, int.c)
        if int.c[i] ≠ zero(TT)
            tᵢ = cache.t̅ + int.Δt * int.c[i]
            int.equation.v[int.f[i]](tᵢ, cache.q, cache.v, int.c[i] * int.Δt)
            cache.q .= cache.v
        end
    end

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


text:
projectionmatrix(us...) = sum(u*u' for u in map(normalize, us)) / length(us)

struct SimpleSimulation
    N::Integer
    mol::Molecule
    μ
    Δ
    H0
    ΓL
    ΓR
    ΔΓL
    ΔΓR
    Lx
    Ly
    Lz
    Lxs
    Lys
    Lzs
end

function SimpleSimulation(N::Integer; η = 1e-15, ϕ = π/4, l=1.4, handedness=1, α=1., Γ=1/27, θ=nothing, δz=nothing)
    # Construct molecule and leads
    mol = makemol(N, ϕ, l, handedness=handedness, θ=θ, δz=δz)
    skt = loaddir(joinpath(PARAMDIR, "ob2-1-1", "base"))
    hydrogen_atoms = collect(filter(x -> x isa Hydrogen, mol.atoms))
    carbon_atoms = collect(filter(x -> x isa Carbon, mol.atoms))

    H0 = hamiltonian(Float64, skt, mol)
    energies, eigvecs = molecularorbitals(H0)
    μ, Δ = real.(chemicalpotential(mol, energies))

    # ΓL = Γ * sum((u = Vector{Float64}(mol, carbon_atoms[1], i); u*u'/3) for i in 1:4)

    # ΓR = Γ * sum((u = Vector{Float64}(mol, carbon_atoms[end], i); u*u'/3) for i in 1:4)

    Random.seed!(47)
    ΓL = 1e-2*Γ * sum((u = Vector{Float64}(mol, carbon_atoms[1+j], i); u*u') for i in 1:4, j in 0:0)

    ΓR1 = 1e-2*Γ * sum((u = Vector{Float64}(mol, carbon_atoms[end-j], i); (rand() + 0.5)*u*u') for i in 1:4, j in 0:0)
    ΓR2 = 1e-2*Γ * sum((u = Vector{Float64}(mol, carbon_atoms[end-j], i); (rand() + 0.5)*u*u') for i in 1:4, j in 0:0)
    ΓR = (ΓR1 .+ ΓR2) ./ 2
    ΔΓR = (ΓR1 .- ΓR2) ./ 2

    Lxs = [angularmomentum(:x, mol, i) for (i, atom) in enumerate(mol) if atom isa Carbon]
    Lys = [angularmomentum(:y, mol, i) for (i, atom) in enumerate(mol) if atom isa Carbon]
    Lzs = [angularmomentum(:z, mol, i) for (i, atom) in enumerate(mol) if atom isa Carbon]
    Lx = angularmomentum(:x, mol)
    Ly = angularmomentum(:y, mol)
    Lz = angularmomentum(:z, mol)

    return SimpleSimulation(N, mol, μ, Δ, H0, ΓL, ΓR, 0I, ΔΓR, Lx, Ly, Lz, Lxs, Lys, Lzs)
end

function HeliceneSimulation(N::Integer; η = 1e-15, ϕ = 0.5, l=1.4, handedness=1, α=1., Γ=1/27, oxidation=0, electric_field=(x, y, z) -> 0.0)
    # Construct molecule and leads
    mol = makehelicene2(N, ϕ, l, handedness=handedness)
    skt = loaddir(joinpath(PARAMDIR, "ob2-1-1", "base"))
    hydrogen_atoms = collect(filter(x -> x isa Hydrogen, mol.atoms))
    carbon_atoms = collect(filter(x -> x isa Carbon, mol.atoms))
    oxygen_atoms = collect(filter(x -> x isa Oxygen, mol.atoms))
    nitrogen_atoms = collect(filter(x -> x isa Nitrogen, mol.atoms))

    H0 = hamiltonian(Float64, skt, mol, electric_field=electric_field)
    energies, eigvecs = molecularorbitals(H0)
    μ, Δ = real.(chemicalpotential(mol, energies, oxidation))

    Random.seed!(47)
    first_indices = Molecules.indices(mol, oxygen_atoms[1])
    last_indices = Molecules.indices(mol, nitrogen_atoms[end])
    HLM = begin
        U = fill(0.0, size(H0))
        A = rand(4, 4)
        U[first_indices, first_indices] = (A + A') / 2
        U
    end
    HRM = begin
        U = fill(0.0, size(H0))
        A = rand(4, 4)
        U[last_indices, last_indices] = (A .+ A') ./ 2
        U ./ tr(U)
    end

    ΓL = 1e-2*Γ * sum((u = Vector{Float64}(mol, oxygen_atoms[1+j], i); u*u') for i in 1:4, j in 0:0)
    #ΓL .= HLM'ΓL*HLM

    ΓR1 = 1e-2*Γ * sum((u = Vector{Float64}(mol, oxygen_atoms[end-j], i); u*u') for i in 1:4, j in 0:0)
    ΓR2 = 1e-2*Γ * sum((u = Vector{Float64}(mol, carbon_atoms[end-j], i); u*u') for i in 1:4, j in 0:0)
    #ΓR2 .= HRM'ΓR2*HRM
    ΓR = (ΓR1 .+ ΓR2) ./ 2
    ΔΓR = (ΓR1 .- ΓR2) ./ 2

    Lxs = [angularmomentum(:x, mol, i) for (i, atom) in enumerate(mol) if typeof(atom) <: Union{Carbon, Nitrogen, Oxygen}]
    Lys = [angularmomentum(:y, mol, i) for (i, atom) in enumerate(mol) if typeof(atom) <: Union{Carbon, Nitrogen, Oxygen}]
    Lzs = [angularmomentum(:z, mol, i) for (i, atom) in enumerate(mol) if typeof(atom) <: Union{Carbon, Nitrogen, Oxygen}]
    Lx = angularmomentum(:x, mol)
    Ly = angularmomentum(:y, mol)
    Lz = angularmomentum(:z, mol)

    return SimpleSimulation(N, mol, μ, Δ, H0, ΓL, ΓR, 0I, ΔΓR, Lx, Ly, Lz, Lxs, Lys, Lzs)
end

function new_coupling(rs, mol; γ = 1/27)
    carbon_atoms = collect(filter(x -> x isa Carbon, mol.atoms))
    uL = [Vector{Float64}(mol, carbon_atoms[1], i) for i in 2:4]
    k = sum(r*u for (u, r) in zip(uL, rs))
    return γ*k*k'
end

function theorem_f(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3]), γ=1/27
)
    ΓL, ΓR, γL = sim.ΓL, sim.ΓR, new_coupling(rs, sim.mol, γ=γ)
    H0 = Matrix(sim.H0)
    ΓA = ΓL*(1 - m) + γL*m

    Λ = λ*sim.Lz
    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    A = -im.*(G .- G') ./ 2
    t = real(dot(G*ΓA*G', ΓR))
    s = 2real(dot(ΓA, G'ΓR*G*Λ*G))
    # function T(b, Ln)
    #     G = inv(E*I - H0 + im*(ΓL + ΓR)/2 - λ*b*Ln)
    #     return real(dot(G*ΓA*G', ΓR))
    # end
    # p = derivative(x -> log(T(x, sim.Lz)), 0.0)

    return t, s, A
end

function theorem_f2(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3]), γ=1/27
)
    ΓL, ΓR, γL = sim.ΓL, sim.ΓR, new_coupling(rs, sim.mol, γ=γ)
    H0 = Matrix(sim.H0)
    ΓA = ΓL*(1 - m) + γL*m

    Λ = λ*sim.Lz
    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    t = real(dot(G*ΓA*G', ΓR))
    s = 2real(dot(ΓA, G'ΓR*G*Λ*G))
    # function T(b, Ln)
    #     G = inv(E*I - H0 + im*(ΓL + ΓR)/2 - λ*b*Ln)
    #     return real(dot(G*ΓA*G', ΓR))
    # end
    # p = derivative(x -> log(T(x, sim.Lz)), 0.0)

    return t, s / t
end

function theorem_f2_each(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3])
)
    ΓL, ΓR, γL = sim.ΓL, sim.ΓR, new_coupling(rs, sim.mol)
    H0 = Matrix(sim.H0)
    ΓA = ΓL*(1 - m) + γL*m

    function T(b, Ln)
        G = inv(E*I - H0 + im*(ΓL + ΓR)/2 - λ*b*Ln)
        return real(dot(G*ΓA*G', ΓR))
    end

    t = T(0, sim.Lz)
    return t, [(
        derivative(x -> log(T(x, Lx)), 0.0),
        derivative(x -> log(T(x, Ly)), 0.0),
        derivative(x -> log(T(x, Lz)), 0.0))
            for (Lx, Ly, Lz) in zip(sim.Lxs, sim.Lys, sim.Lzs)]
end

function theorem_f2_approx(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3])
)
    ΓL, ΓR, γL = sim.ΓL, sim.ΓR, new_coupling(rs, sim.mol)
    H0 = Matrix(sim.H0)
    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    eig = eigen(E*I - H0)

    valmin =
        #(0, 36)
        findmin(inv.(eig.values))
    valmax =
        #(0, 37)
        findmax(inv.(eig.values))

    us = [
        #eig.vectors[:, valmin[2] - 1],
        eig.vectors[:, valmin[2]],
        eig.vectors[:, valmax[2]]
        #eig.vectors[:, valmax[2] + 1]
    ]
    #us = [eig.vectors[:, i] for i in 1:72]
    P = sum(u*u' for u in us)

    ΓA = ΓL*(1 - m) + γL*m

    local T, d
    let P = P, H0 = P*H0*P, ΓL = P*ΓL*P, ΓR = P*ΓR*P, ΓA = P*ΓA*P
        function T(b, Ln)
            Ln = P*Ln*P
            G = P*inv((E*I - H0 + im*(ΓL + ΓR)/2 - λ*b*Ln))*P
            return real(dot(G*ΓA*G', ΓR))
        end
        d = dot(ΓL, ΓA) / dot(ΓL, ΓL)#d = sqrt(tr((ΓL - ΓA).^2)) / sqrt(tr(ΓL.^2))
    end

    t = T(0, sim.Lz)
    p = derivative(x -> log(T(x, sim.Lz)), 0.0)
    return t, p*t, d, valmin, valmax
end

function theorem_f2_approx2(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3])
)
    ΓL, ΓR, γL = sim.ΓL, sim.ΓR, new_coupling(rs, sim.mol)
    H0 = Matrix(sim.H0)
    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    eig = eigen(E*I - H0)

    vals = collect(enumerate(abs.(inv.(eig.values))))
    sort!(vals, by=x->x[2])
    valmin =
        #(0, 36)
        #findmin(inv.(eig.values))
        vals[end]
    valmax =
        #(0, 37)
        #findmax(inv.(eig.values))
        vals[end-1]

    us = [
        #eig.vectors[:, valmin[2] - 1],
        eig.vectors[:, valmin[1]],
        eig.vectors[:, valmax[1]]
        #eig.vectors[:, valmax[2] + 1]
    ]
    P = sum(u*u' for u in us)

    ΓA = ΓL*(1 - m) + γL*m

    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    Λ = λ*sim.Lz
    A = -im.*(G .- G') ./ 2

    #t = real(tr(ΓA*P*G'P*ΓR*P*G*P))

    Λ12 = us[1]'Λ*us[2]

    A1 = us[1]'A*us[1]
    A2 = us[2]'A*us[2]
    A12 = us[1]'A*us[2]
    A = [A1 A12; A12 A2]

    γ11 = us[1]'ΓA*us[1]
    γ12 = us[1]'ΓA*us[2]
    γ22 = us[2]'ΓA*us[2]
    γ = [γ11 γ12; γ12 γ22]

    Γ11 = us[1]'ΓR*us[1]
    Γ12 = us[1]'ΓR*us[2]
    Γ22 = us[2]'ΓR*us[2]
    Γ = [Γ11 Γ12; Γ12 Γ22]

    prefactor = real(sum(A[1,i]*A[j,k]*A[3-i,2]*(γ[j,i]*Γ[3-i,k] - γ[j,3-i]*Γ[i,k]) for i in 1:2, j in 1:2, k in 1:2))

    return 2imag(Λ12), prefactor
end

function theorem_magnet(E, m=0; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3]), γ=1/27
)
    ΓL, ΓR, ΔΓR = sim.ΓL, sim.ΓR, sim.ΔΓR
    H0 = Matrix(sim.H0)
    ΓRup = (ΓR .* (2 - m) .+ m .* ΔΓR)
    ΓRdown = (ΓR .* (2 - m) .- m .* ΔΓR)

    Λ = λ*sim.Lz
    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    t = real(dot(ΓL, G'ΓR*G))
    #s = real(dot(ΓL, G'ΔΓR*G*Λ*G))
    t += 2imag(dot(ΓL, G'ΔΓR*G*ΔΓR*G))
    s = real(dot(ΓRup, G'ΓRdown*G*Λ*G))

    return t, s
end

function β(γ, a, b)
    A = similar(γ)
    N, M = size(A)
    for j in 1:M, i in 1:N
        A[j, i] = (j == a && i == b) || (j == b && i == a) ? 1.0 : 0.0
    end
    return A
end

export theorem_derivative
function theorem_derivative(E, θ, ϕ; sim, λ=1e-2 / 27, rs = normalize!([(0.5 - rand()) for _ in 1:3]), γ=1/27
)
    ΓL, ΓR = sim.ΓL, sim.ΓR
    H0 = Matrix(sim.H0)
    #βs = [-sin(ϕ)*sin(θ), sin(ϕ)*cos(θ), cos(ϕ)]
    βs = [cos(ϕ)*cos(θ), cos(ϕ)*sin(θ), -sin(ϕ)]
    γ = sum(βs[i] * β(sim.ΓL, i+2, i+2) for i in 1:3)

    G = inv(E*I - H0 + im*(ΓL + ΓR)/2)
    Λ = λ*sim.Lz

    t = real(dot(γ, G'ΓR*G))
    s = 2*real(dot(γ, G'ΓR*G*Λ*G))

    return t, s, s / t
end

## Lattice

struct LatticeSimulation
    uc::Molecules.UnitCell
    H0
    Lx
    Ly
    Lz
end

function LatticeSimulation(; η = 1e-15, ϕ = π/4, l=1.4, handedness=1, α=1., Γ=1/27, θ=nothing, δz=nothing, ξs = range(0, π, length=200))
    skt = loaddir(joinpath(PARAMDIR, "ob2-1-1", "base"))
    uc = makelattice(ϕ, l)
    Lx = Molecules.angularmomentum(:x, uc.unit)
    Ly = Molecules.angularmomentum(:y, uc.unit)
    Lz = Molecules.angularmomentum(:z, uc.unit)

    H0(ξ) = hamiltonian(Complex{Float64}, skt, uc, ξ)

    return LatticeSimulation(uc, H0, Lx, Ly, Lz)
end


text:
module CaesarM

export caesarshift, rot13

# Caesar shift function
function caesarshift(text::String, alphabet::String, shift::Int)

    # Shift nmubers
    return join([
        alphabet[i+1]
        for i in [
            ((x + shift - 1) % length(alphabet))
            for x in [
                findfirst(isequal(letter), alphabet)
                for letter in text if occursin(letter, alphabet)
            ]
        ]
    ])

end # function

# ROT13 because why not
"""
    rot13(text::String, alphabet::String)

Returns decrypted/encrypted (they're the same process) text as a string.
"""
function rot13(text::String, alphabet::String)
    return caesarshift(text, alphabet, 13)
end # function

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
module PyFOOOF

using PyCall

#####
##### init
#####

const fooof = PyNULL()

function __init__()
    # all of this is __init__() so that it plays nice with precompilation
    # see https://github.com/JuliaPy/PyCall.jl/#using-pycall-from-julia-modules
    copy!(fooof, pyimport("fooof"))
    # don't eval into the module while precompiling; this breaks precompilation
    # of downstream modules
    if ccall(:jl_generating_output, Cint, ()) == 0
        # delegate everything else to fooof
        for pn in propertynames(fooof)
            isdefined(@__MODULE__, pn) && continue
            prop = getproperty(fooof, pn)
            @eval $pn = $prop
        end
    end
    return nothing
end

end #module


text:
module Methods

export classmethods, classdeclaredmethods,
    isstatic, ispublic, MethodDescriptor

using JavaCall.JNI
using JavaCall.Signatures
using JavaCall.Conversions
using JavaCall.Core

using JavaCall.Reflection: Classes, Modifiers

#=
Struct to hold information aboud java methods
used to generate functions that call the jni
This should not be a replacement of java.lang.reflect.Method
as it should only store essential information

name:       Name of the method

rettype:    Descriptor for the return class (void methods also
            have a descriptor with Nothing as the juliatype,
            jvoid as the jni type and V signature)

paramtypes: List of descriptors with the parameter types

modifiers:  Information about the method modifiers
=#
struct MethodDescriptor
    name::String
    rettype::Classes.ClassDescriptor
    paramtypes::Vector{Classes.ClassDescriptor}
    modifiers::Modifiers.ModifiersDescriptor
end

function Base.show(io::IO, m::MethodDescriptor)
    print(
        io,
        "MethodDescriptor{name: ", m.name,
        ", ret: ", string(m.rettype),
        ", params: ", string(m.paramtypes),
        ", modifiers: ", string(m.modifiers),
        "}")
end

function Base.:(==)(x::MethodDescriptor, y::MethodDescriptor)
    x.name == y.name &&
    x.rettype == y.rettype &&
    x.paramtypes == y.paramtypes  &&
    x.modifiers == y.modifiers
end

isstatic(m::MethodDescriptor) = m.modifiers.static

ispublic(m::MethodDescriptor) = m.modifiers.public

function descriptorfrommethod(method::jobject)
    name = callinstancemethod(method, :getName, Symbol("java.lang.String"), [])
    rettype = callinstancemethod(method, :getReturnType, Symbol("java.lang.Class"), [])
    paramtypes = callinstancemethod(method, :getParameterTypes, Vector{Symbol("java.lang.Class")}, [])
    MethodDescriptor(
        convert_to_string(String, name),
        Classes.descriptorfromclass(rettype),
        map(Classes.descriptorfromclass, convert_to_vector(Vector{jclass}, paramtypes)),
        Modifiers.methodmodifiers(method))
end

classmethods(classname::Symbol) = classmethods(Classes.findclass(classname))

function classmethods(classdescriptor::Classes.ClassDescriptor)
    array = convert_to_vector(Vector{jobject}, callinstancemethod(
        classdescriptor.jniclass,
        :getMethods,
        Vector{Symbol("java.lang.reflect.Method")},
        []))
    map(descriptorfrommethod, array)
end

classdeclaredmethods(classname::Symbol) = classdeclaredmethods(Classes.findclass(classname))

function classdeclaredmethods(classdescriptor::Classes.ClassDescriptor)
    array = convert_to_vector(Vector{jobject}, callinstancemethod(
        classdescriptor.jniclass,
        :getDeclaredMethods,
        Vector{Symbol("java.lang.reflect.Method")},
        []))
    map(descriptorfrommethod, array)
end

end


text:
@testset "Convergence on the Rosenbrock function" begin
    dim = 2
    x0 = fill(-1, dim)
    @testset "By gradient norm (default)" begin
        for method in OPT_TYPES
            opt = method(x0)
            optresult = optimize!(
                rosenbrock!, opt, x0;
                maxiter=1000,
                constrain_step=(x0,d)->Inf,
            )
            @testset "$method" begin
                if opt isa FixedRateDescent
                    @test_broken optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                else
                    @test optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                end
            end
        end
    end
    @testset "By x - xpre" begin
        for method in OPT_TYPES
            opt = method(x0)
            optresult = optimize!(
                rosenbrock!, opt, x0;
                maxiter=1000,
                constrain_step=(x0,d)->Inf,
                convcond=(x, xpre, y, ypre, g) -> norm(x - xpre, 2) ≤ 1e-6,
            )
            @testset "$method" begin
                if any(T -> opt isa T, (HyperGradDescent, FixedRateDescent, MomentumDescent, NesterovMomentum))
                    @test_broken optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                else
                    @test optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                end
            end
        end
    end
    @testset "By y - ypre" begin
        for method in OPT_TYPES
            opt = method(x0)
            optresult = optimize!(
                rosenbrock!, opt, x0;
                maxiter=1000,
                constrain_step=(x0,d)->Inf,
                convcond=(x, xpre, y, ypre, g) -> isapprox(y, ypre; atol=1e-6),
            )
            @testset "$method" begin
                if any(T -> opt isa T, (HyperGradDescent, FixedRateDescent, MomentumDescent, NesterovMomentum))
                    @test_broken optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                else
                    @test optresult.converged
                    @test optresult.argument ≈ [1, 1] rtol=0.05
                end
            end
        end
    end
end


text:
using Zygote

N, W, R, B= 3, 3, 1, 1

function expand(arr::AbstractVecOrMat, dims...)
    res = rand(Float32, dims...)
    res[:, :, 1] = Float32.(arr)
    res
end

M = expand([1 2 3;
           -1 -2 -3;
           1 -2 3], N, N, B)

# Content-based read
contentread = (
    kr = expand([1, 2, 0], W, R, B),
    βr = expand([100], R, B),
    f = expand([1], R, B),
    readmode = expand([0, 1, 0], (3, R, B))
    )

# Temporal linkage based read
backwardread = (
    kr =expand([-1, -2, -5], W, R, B),
    βr = expand([100], R, B),
    f = expand([0], R, B),
    readmode = expand([1, 0, 0], 3, R, B)
    )

forwardread = (
    kr = expand([-1,-2,-5], W, R, B),
    βr = expand([100], R, B),
    f = expand([0], R, B),
    readmode = expand([0,0,1], 3, R, B)
    )

contentwrite = (
    kw = expand([1,2, 0], W, R, B),
    βw = expand([10], R, B),
    e = expand([1, 1, 1], W, B),
    v = expand([10, 20, 30], W, B),
    ga = expand([0], 1, B),
    gw = expand([1], 1, B),
    f = expand([1], R, B)
)

allocationwrite = (
    kw = expand([1,1,1], W, R, B),
    βw = expand([1], R, B),
    e = expand([1, 1, 1], W, B),
    v = expand([10, 20, 30], W, B),
    ga =expand([1], 1, B),
    gw =expand([1], 1, B),
    f = expand([1], R, B)
)

state = DNC.State(
    # Write history is 1 -> 2 -> 3
    expand([0 0 0;
            1 0 0;
            0 1 0], N, N, B),
    zeros(Float32, N, B),
    expand([1, 1, 0], N, B), # 3 is (artificially) set to unused/ free to be allocated
    expand([0, 1, 0], N, 1, B),
    expand([0, 1, 0], N, R, B) # last read was 2, so forward points to 3, back to 1
)

@testset "Sharp read/write" begin
    @testset "Readweights" begin
        L, wr = state.L, state.wr
        r1 = DNC.readweights(M, contentread, L, wr)
        @test eltype(r1) == Float32
        @test round.(r1; digits=5)[:,1] == [1.0, 0, 0]
        r2 = DNC.readweights(M, backwardread, L, wr)
        @test eltype(r2) == Float32
        @test round.(r2; digits=5)[:, 1] == [1.0, 0, 0]
        r3 = DNC.readweights(M, forwardread, L, wr)
        @test eltype(r3) == Float32
        @test round.(r3; digits=5)[:, 1] == [0.0, 0, 1]
    end
    @testset "Writeweights" begin
        ww, wr, prev_u = state.ww, state.wr, state.u
        u = DNC.usage(prev_u, ww, wr, contentwrite.f)
        ww1 = DNC.writeweights(M, contentwrite, u)
        @test eltype(ww1) == Float32
        @test round.(ww1; digits=3)[:, 1, 1] == [1.0, 0.0, 0.0]
        u = DNC.usage(prev_u, ww, wr, allocationwrite.f)
        ww2 = DNC.writeweights(M, allocationwrite, u)
        @test eltype(ww2) == Float32
        @test round.(ww2; digits=3)[:, 1, 1] == [0.0, 0.0, 1.0]
    end
end

@testset "Erase and add" begin
    B = 1
    u = DNC.usage(state.u, state.ww, state.wr, allocationwrite.f)
    new = DNC.eraseandadd(M, reshape([0.0f0, 0, 1], N, 1, B), ones(Float32, 3, B), 2*ones(Float32, W, B))
    @test eltype(new) == Float32
    @test new[3,:, 1] == [2.0, 2.0, 2.0]
    ww2 = round.(DNC.writeweights(M, allocationwrite, u); digits=3)
    new2 = DNC.eraseandadd(M, ww2, allocationwrite[:e], allocationwrite[:v])
    @test new2[3,:, 1] == allocationwrite[:v][:, 1]
end

using Random
rng = MersenneTwister(234)

@testset "MemoryAccess" begin
    insize, N, W, R, B = 20, 5, 10, 2, 1
    ma = DNC.MemoryAccess(insize, N, W, R, B)
    inputs = rand(rng, Float32, insize, B)
    @testset "Dimensions" begin
        @test size(ma.M) == (N, W, B)
        @test size(ma.state.wr) == (N, R, B)
        readvectors = ma(inputs)
        @test eltype(readvectors) == Float32
        @test size(readvectors) == (W, R, B)
    end
end

insize, N, W, R, B = 20, 5, 10, 2, 4
ma = DNC.MemoryAccess(insize, N, W, R, B)
M = rand(rng, Float32, N, W, B)
L = rand(rng, Float32, N, N, B)
prev_wr = rand(rng, Float32, N, R, B)
inputsraw = rand(rng, Float32, insize, B)
mappings = DNC.inputmappings(20, R, W)
inputs = DNC.split_ξ(inputsraw, mappings)
state = DNC.State(N, R, B)

@testset "Batch training" begin
    r = DNC.readweights(M, inputs, L, prev_wr)
    @test size(r) == (N, R, B)
    u = DNC.usage(state.u, state.ww, state.wr, inputs.f)
    w = DNC.writeweights(M, inputs, u)
    @test size(w) == (N, 1, B)
    newM = DNC.eraseandadd(M, w, inputs.e, inputs.v)
    @test size(newM) == (N, W, B)
    res = ma(inputsraw)
    @test size(res) == (W, R, B)
end


text:
# ---
# title: 539. Minimum Time Difference
# id: problem539
# author: Indigo
# date: 2021-06-23
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/minimum-time-difference/description/>
# hidden: true
# ---
#
# Given a list of 24-hour clock time points in **" HH:MM"** format, return _the
# minimum **minutes** difference between any two time-points in the list_.
#
#
#
# **Example 1:**
#
#
#
#     Input: timePoints = ["23:59","00:00"]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: timePoints = ["00:00","23:59","00:00"]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `2 <= timePoints <= 2 * 104`
#   * `timePoints[i]` is in the format **" HH:MM"**.
#
#
## @lc code=start
using LeetCode

function find_min_difference(time_points::Vector{String})
    ttp = Tuple{Int, Int}[]
    for t in time_points
        st = split(t, ':')
        push!(ttp, (parse(Int, st[1]), parse(Int, st[2])))
    end
    sort!(ttp)
    push!(ttp, (ttp[1][1] + 24, ttp[1][2]))
    return minimum((ttp[i][1] - ttp[i-1][1]) * 60 + (ttp[i][2] - ttp[i-1][2]) for i in 2:length(ttp))
end
## @lc code=end


text:
module TestLEM
	using ManifoldLearning
	using Base.Test

	k = 12
	d = 2
	X, L = swiss_roll()
	Y = transform(LEM, X, k=k, d=d)

	@test outdim(Y) == d
	@test size(projection(Y), 2) == size(X, 2)
	@test neighbors(Y) == k
	@test length(eigvals(Y)) == d
	@test length(ccomponent(Y)) > 1
end


text:
function conv(path::String)
    data = []
    lines = open(readlines, path)
    lines = map(chomp, lines)
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line)
            push!(data, "")
            continue
        end
        items = split(line, " ")
        word = String(items[1])
        tag = String(items[4])
        if startswith(tag, "I-")
            if i == 1 || isempty(lines[i-1])
                tag = replace(tag, "I-", "")
            else
                ptag = String(split(lines[i-1])[4])
                if ptag == "O"
                    tag = replace(tag, "I-", "")
                elseif startswith(ptag, "B-")
                    tag = "_"
                elseif startswith(ptag, "I-")
                    tag = "_"
                else
                    throw("$ptag")
                end
            end
        elseif startswith(tag, "B-")
        elseif tag == "O"
        else
            throw("Invalid tag: $tag.")
        end
        push!(data, "$(word)\t$(tag)")
    end
    open("a.out", "w") do f
        foreach(x -> println(f,x), data)
    end
end

function convIOBES(path::String)
    data = []
    lines = open(readlines, path)
    lines = map(chomp, lines)
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line)
            push!(data, "")
            continue
        end
        items = split(line, "\t")
        word = String(items[1])
        tag = String(items[2])
        if startswith(tag, "B-") || startswith(tag, "S-")
            tag = tag[3:end]
        elseif startswith(tag, "I-") || startswith(tag, "E-")
            tag = "_"
        end
        push!(data, "$(word)\t$(tag)")
    end
    open("a.out", "w") do f
        foreach(x -> println(f,x), data)
    end
end

workspace()
convIOBES(joinpath(dirname(@__FILE__), ".data/eng.train.IOBES"))

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
@testset "Conjugated orbitals" begin
    o = o"1s"
    co = conj(o)
    @test co == Conjugate(o)
    @test string(co) == "1s†"
end


text:
using AssigmentProblem
using Test

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


text:
module SweatVideo

using ColorTypes
using VideoIO
using ImageFiltering
using Graphs
using CSVFiles
using FileIO
# using FFTW

@assert 1==1
# const minimum_blob_size = 350 #35
include("read_video.jl")
include("find_blobs.jl")
include("changepoint.jl")
include("morphology.jl")
include("segment_blobs0.jl")
# include("grow_blobs_gpu.jl")
include("grow_blobs.jl")
include("overlay_spots.jl")
include("colorutil.jl")
#include("util/read_glands.jl")
#include("datafiles.jl") # in sweat metapaackge
include("extract_glands_table.jl")
export readpngs, readvideo, readcsv,
    colortype,
    find_split, do_split,
    catcolor, #addcolor, colortype,
    internalgradient, externalgradient, gradient,
    opening, closing,
    segment_blobs,
    # segment_blobs2,
    segment_blobs0,
    # segment_blobs_gpu,
    overlay_spots,
    align_frames,
#    readglands,
    datalist, paths,
    extractglands
# function __init__()
#     init() # GPU context
# end
end # module


text:
@testset "exceptions" begin
    @test_throws AssertionError Interval(+Inf,-Inf) == emptyinterval()
end


text:
mutable struct Model{T,Tv,Tt}
    n::Int  # Num vars
    m::Int  # Num cons
    x::Tv  # Starting and final solution
    x_L::Tv # Variables Lower Bound
    x_U::Tv # Variables Upper Bound
    g::Tv  # Final constraint values
    g_L::Tv # Constraints Lower Bound
    g_U::Tv # Constraints Upper Bound
    j_str::Tt
    h_str::Tt
    mult_g::Tv # lagrange multipliers on constraints
    mult_x_L::Tv # lagrange multipliers on lower bounds
    mult_x_U::Tv # lagrange multipliers on upper bounds
    obj_val::T  # Final objective
    status::Int  # Final status
    g_order::Vector{Int}

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Union{Function,Nothing}

    intermediate  # Can be nothing

    # For MathProgBase
    sense::Symbol
    convex_model

    parameters::Parameters
    statistics::Dict{String,Any}   # collects parameters of all iterations inside the algorithm if StatisticsFlag > 0

    Model(
        n::Int,
        m::Int,
        x_L::Tv,
        x_U::Tv,
        g_L::Tv,
        g_U::Tv,
        j_str::Tt,
        h_str::Tt,
        eval_f::Function,
        eval_g::Function,
        eval_grad_f::Function,
        eval_jac_g::Function,
        eval_h::Union{Function,Nothing},
        parameters::Parameters,
        g_order, convex_model
    ) where {T, Tv<:AbstractArray{T}, Tt<:AbstractArray{Tuple{Int64,Int64}}} = new{T,Tv,Tt}(
        n, m,
        zeros(n), x_L, x_U,
        zeros(m), g_L, g_U,
        j_str, h_str,
        zeros(m), zeros(n), zeros(n),
        0.0,
        -5,
        g_order,
        eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h,
        nothing, :Min, convex_model,
        parameters,
        Dict{String,Any}()
    )
end

function optimize!(model::Model)
    if isnothing(model.parameters.external_optimizer)
    	model.status = -12;
        @error "`external_optimizer` parameter must be set for subproblem solutions."
    else
        if model.parameters.method == "SLP"
            if model.parameters.algorithm == "Line Search"
                slp = SlpLS(model)
            elseif model.parameters.algorithm == "Trust Region"
                slp = SlpTR(model)
            end
            run!(slp)
        else
            @error "The method is not defined"
        end
    end
    return nothing
end

function add_statistic(model::Model, name::String, value)
    if model.parameters.StatisticsFlag == 0
        return
    end
    model.statistics[name] = value
end

function add_statistics(model::Model, name::String, value::T) where T
    if model.parameters.StatisticsFlag == 0
        return
    end
    if !haskey(model.statistics, name)
        model.statistics[name] = Array{T,1}()
    end
    push!(model.statistics[name], value)
end


text:
using FEASTSolver
using LinearAlgebra
using DelimitedFiles
using SparseArrays

N = 1000
A = spdiagm(-1 => fill(-1.0, N-1), 0 => fill(2.0, N), 1 => fill(-1.0, N-1))
# A = diagm(-1 => fill(-1.0, N-1), 0 => fill(2.0, N), 1 => fill(-1.0, N-1))
# A = I + sprand(100, 100, 0.1)
R = 500.0/(N^2)
C = complex(R, 0.0) # contour
# C, R = complex(0.1, 0.0), 0.1 # contour
# X = rand(ComplexF64, N, 20) # initial subspace
X = sprand(ComplexF64, N, 20, 0.05)

contour = circular_contour_trapezoidal(C, R, 16)
# contour = circular_contour_gauss(C, R, 16)
# contour = rectangular_contour_gauss(0.0 - R*im, 2*R + R*im, 32)
# contour = rectangular_contour_trapezoidal(0.0 - R*im, 2*R + R*im, 8)
# contour = rectangular_contour_gauss(0.9 - R*im, 1.1 + R*im, 128)
writedlm("test/data/contour.dat", zip(real.(contour.nodes), imag.(contour.nodes)))
writedlm("test/data/weights.dat", zip(real.(contour.weights), imag.(contour.weights)))

# z = 0.025+0.0im
# println(rational_func(z, contour))

# e, v, res = feast!(deepcopy(X), A, contour, debug=true, ϵ=10e-15, iter=0, store=false)
# @timev e, v, res = feast!(deepcopy(X), A, contour, debug=true, ϵ=10e-15, store=false, mixed_prec=false)
# @timev e, v, res = feast!(X, A, contour, debug=true, ϵ=10e-15, store=true)

println(contour_estimate_eig(A, contour, debug=true, samples=100))

# convergence_info(e,v,res,contour)
# display(e)


text:
module YaoCompiler

export @device,
    @gate,
    @ctrl,
    @measure,
    @barrier,
    compile,
    YaoInterpreter,
    YaoCompileTarget,
    JLGenericTarget,
    NativeJuliaTarget,
    MeasureResult,
    TargetHostKernel,
    HardwareFreeOptions,
    Routine,
    GenericRoutine,
    IntrinsicRoutine,
    Operation,
    AdjointOperation,
    routine_name,
    IntrinsicError,
    AnyReg,
    # reexport YaoLocations
    Locations,
    CtrlLocations,
    # reflection
    @yao_code_lowered,
    @yao_code_typed

using MLStyle
using YaoAPI
using LLVM
using Expronicon
using YaoHIR
using YaoLocations
using TimerOutputs
using LinearAlgebra
using GPUCompiler
using ZXCalculus
using Configurations
using CompilerPluginTools
using LLVM.Interop
using YaoHIR: routine_name, measure_cmp, MeasureResult, Typed
using GPUCompiler: CodeCache, CompilerJob, AbstractCompilerTarget, AbstractCompilerParams, WorldView
using YaoLocations: map_check, map_check_nothrow, map_error, plain, unsafe_mapping
using CompilerPluginTools: Argument
using Base.Meta: ParseError

const to = TimerOutput()
timings() = (TimerOutputs.print_timer(to); println())
enable_timings() = (TimerOutputs.enable_debug_timings(Compiler); return)

@as_record Locations
@as_record CtrlLocations

"""
    AnyReg <: AbstractRegister{1}

A place holder for registers when compilation is not register specific.
"""
struct AnyReg <: AbstractRegister{1} end
Base.show(io::IO, ::AnyReg) = print(io, "AnyReg()")

include("compiler/intrinsics.jl")
include("compiler/syntax.jl")
include("compiler/interp.jl")
include("compiler/reflection.jl")

include("codegen/llvmopt.jl")
include("codegen/native.jl")

# We have one global JIT and TM
const orc = Ref{LLVM.OrcJIT}()
const tm = Ref{LLVM.TargetMachine}()

function __init__()
    TimerOutputs.reset_timer!(to)
    opt_level = Base.JLOptions().opt_level
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end

    tm[] = LLVM.JITTargetMachine(; optlevel = optlevel)
    LLVM.asm_verbosity!(tm[], true)

    if LLVM.has_orc_v1()
        orc[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
    end

    atexit() do
        return LLVM.dispose(orc[])
    end
end

end # module


text:
using SQLite
using ZulipSnippetBot

include("configuration.jl")

db = SQLite.DB(DB)

ZulipSnippetBot.up(db)

# There should be `down` function and some way to control db version...

max_stars_repo_name: Arkoniak/ZulipSnippetBot

max_issues_repo_name: Arkoniak/ZulipSnippetBot

max_forks_repo_name: Arkoniak/ZulipSnippetBot


text:
# # Forest fire

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../forest.mp4" type="video/mp4">
# </video>
# ```

# The forest fire model is defined as a cellular automaton on a grid.
# A position can be empty or occupied by a tree which is ok, burning or burnt.
# We implement a slightly different ruleset to that of
# [Drossel and Schwabl (1992)](https://en.wikipedia.org/wiki/Forest-fire_model),
# so that our implementation can be compared with other ABM frameworks
#
# 1. A burning position turns into a burnt position
# 1. A tree will burn if at least one neighbor is burning

# The forest has an innate `density`, which is the proportion of trees initialized as
# `green`, however all trees that reside on the left side of the grid are `burning`.
# The model is also available from the `Models` module as [`Models.forest_fire`](@ref).

# ## Defining the core structures

# Cellular automata don't necessarily require an agent-like structure. Here we will
# demonstrate how a model focused solution is possible.
using Agents, Random
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

@agent Automata GridAgent{2} begin end
nothing # hide

# The agent type `Automata` is effectively a dummy agent, for which we will invoke
# [`dummystep`](@ref) when stepping the model.

# We then make a setup function that initializes the model.
function forest_fire(; density = 0.7, griddims = (100, 100))
    space = GridSpace(griddims; periodic = false, metric = :euclidean)
    ## The `trees` field is coded such that
    ## Empty = 0, Green = 1, Burning = 2, Burnt = 3
    forest = ABM(Automata, space; properties = (trees = zeros(Int, griddims),))
    for I in CartesianIndices(forest.trees)
        if rand(forest.rng) < density
            ## Set the trees at the left edge on fire
            forest.trees[I] = I[1] == 1 ? 2 : 1
        end
    end
    return forest
end

forest = forest_fire()

# ## Defining the step!

function tree_step!(forest)
    ## Find trees that are burning (coded as 2)
    for I in findall(isequal(2), forest.trees)
        for idx in nearby_positions(I.I, forest)
            ## If a neighbor is Green (1), set it on fire (2)
            if forest.trees[idx...] == 1
                forest.trees[idx...] = 2
            end
        end
        ## Finally, any burning tree is burnt out (2)
        forest.trees[I] = 3
    end
end
nothing # hide

# ## Running the model

Agents.step!(forest, dummystep, tree_step!, 1)
count(t == 3 for t in forest.trees) # Number of burnt trees on step 1

#

Agents.step!(forest, dummystep, tree_step!, 10)
count(t == 3 for t in forest.trees) # Number of burnt trees on step 11

# Now we can do some data collection as well using an aggregate function `percentage`:

Random.seed!(2)
forest = forest_fire(griddims = (20, 20))
burnt_percentage(f) = count(t == 3 for t in f.trees) / prod(size(f.trees))
mdata = [burnt_percentage]

_, data = run!(forest, dummystep, tree_step!, 10; mdata)
data

# Now let's plot the model. We use green for unburnt trees, red for burning and a
# dark red for burnt.
forest = forest_fire()
Agents.step!(forest, dummystep, tree_step!, 1)

plotkwargs = (
    add_colorbar = false,
    heatarray = :trees,
    heatkwargs = (
        colorrange = (0, 3),
        colormap = cgrad([:white, :green, :red, :darkred]; categorical = true),
    ),
)
fig, _ = abm_plot(forest; plotkwargs...)
fig

# or animate it
Random.seed!(10)
forest = forest_fire(density = 0.6)
add_agent!(forest) # Add one dummy agent so that abm_video will allow us to plot.
abm_video(
    "forest.mp4",
    forest,
    dummystep,
    tree_step!;
    as = 0,
    framerate = 5,
    frames = 20,
    spf = 5,
    title = "Forest Fire",
    plotkwargs...,
)
nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../forest.mp4" type="video/mp4">
# </video>
# ```


text:
function movers!(r,c,side)
    for i in 1:c
        move!(r, side)
    end
end

function convertside(side)
    return HorizonSide((Int(side)+2)%4)
end

function task_8!(r)
 n=1
 side=Ost
    while isborder(r, Nord)
        movers!(r, n, side)
        n+=1
        side=convertside(side)
    end
end

task_8!(r)