text:
# Calibration, smoothing, then jump to the OrdinaryDiffEq._postamble!
function OrdinaryDiffEq.postamble!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})

    if isstatic(integ.cache.diffusionmodel) # Calibrate
        # @warn "sol.log_likelihood is not correct for static diffusion models!"
        integ.sol.log_likelihood = NaN
        final_diff = integ.sol.diffusions[end]
        for s in integ.sol.x
            # s.Σ .*= final_diff
            copy!(s.Σ, apply_diffusion(s.Σ, final_diff))
        end

        if isempty(size(final_diff))
            integ.sol.diffusions .= final_diff
        else
            [(d .= final_diff) for d in integ.sol.diffusions]
        end
    end

    if integ.alg.smooth
        smooth_all!(integ)
        integ.sol.pu .= [integ.cache.SolProj * x for x in integ.sol.x]
        @assert (length(integ.sol.u) == length(integ.sol.pu))
        [(su .= pu) for (su, pu) in zip(integ.sol.u, integ.sol.pu.μ)]
    end

    OrdinaryDiffEq._postamble!(integ)
end

function DiffEqBase.savevalues!(
    integrator::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    force_save=false, reduce_size=true)

    # Do whatever OrdinaryDiffEq would do
    out = OrdinaryDiffEq._savevalues!(integrator, force_save, reduce_size)

    # Save our custom stuff that we need for the posterior
    # TODO If we don't want dense output, we might not want to save these!
    # It's not completely clear how to specify that though; They are also needed for sampling.
    OrdinaryDiffEq.copyat_or_push!(integrator.sol.x, integrator.saveiter, integrator.cache.x)
    OrdinaryDiffEq.copyat_or_push!(integrator.sol.diffusions, integrator.saveiter, integrator.cache.diffusion)
    OrdinaryDiffEq.copyat_or_push!(integrator.sol.pu, integrator.saveiter, integrator.cache.SolProj*integrator.cache.x)

    return out
end


text:
using NanStats
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
#@test 1 == 2

vec1=[NaN,3,4,NaN];

mn=nanmean(vec1);
varn=nanvar(vec1);
kurn=nankurt(vec1);
sken=nanskew(vec1);

print("vec1=$(vec1)\n");
print("The mean for vec1 ignoring NaN values is $(mn)\n");
print("The variance for vec1 ignoring NaN values is $(varn)\n")
print("The kurtosis for vec1 ignoring NaN values is $(kurn)\n");
print("The skewness for vec1 ignoring NaN values is $(sken)\n");

ext: jl


text:
@otestset "Emitters" begin
    # TODO!! once emitters are finalised
end # testset Emitters


text:
"""
    space(::SiteType"S=1/2";
          conserve_qns = false,
          conserve_sz = conserve_qns,
          conserve_szparity = false,
          qnname_sz = "Sz",
          qnname_szparity = "SzParity")

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

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

ITensors.val(::ValName"Up", ::SiteType"S=1/2") = 1
ITensors.val(::ValName"Dn", ::SiteType"S=1/2") = 2

ITensors.val(::ValName"↑", st::SiteType"S=1/2") = 1
ITensors.val(::ValName"↓", st::SiteType"S=1/2") = 2

ITensors.val(::ValName"Z+", ::SiteType"S=1/2") = 1
ITensors.val(::ValName"Z-", ::SiteType"S=1/2") = 2

ITensors.state(::StateName"Up", ::SiteType"S=1/2") = [1.0, 0.0]
ITensors.state(::StateName"Dn", ::SiteType"S=1/2") = [0.0, 1.0]

ITensors.state(::StateName"↑", st::SiteType"S=1/2") = [1.0, 0.0]
ITensors.state(::StateName"↓", st::SiteType"S=1/2") = [0.0, 1.0]

ITensors.state(::StateName"Z+", st::SiteType"S=1/2") = [1.0, 0.0]
ITensors.state(::StateName"Z-", st::SiteType"S=1/2") = [0.0, 1.0]

ITensors.state(::StateName"X+", st::SiteType"S=1/2") = [1 / sqrt(2), 1 / sqrt(2)]
ITensors.state(::StateName"X-", st::SiteType"S=1/2") = [1 / sqrt(2), -1 / sqrt(2)]

ITensors.state(::StateName"Y+", st::SiteType"S=1/2") = [1 / sqrt(2), im / sqrt(2)]
ITensors.state(::StateName"Y-", st::SiteType"S=1/2") = [1 / sqrt(2), -im / sqrt(2)]

ITensors.op(::OpName"Z", ::SiteType"S=1/2") = [
  1 0
  0 -1
]

ITensors.op(::OpName"Sz", ::SiteType"S=1/2") = [
  0.5 0.0
  0.0 -0.5
]

ITensors.op(::OpName"Sᶻ", t::SiteType"S=1/2") = op(OpName("Sz"), t)

ITensors.op(::OpName"S+", ::SiteType"S=1/2") = [
  0 1
  0 0
]

ITensors.op(::OpName"S⁺", t::SiteType"S=1/2") = op(OpName("S+"), t)

ITensors.op(::OpName"Splus", t::SiteType"S=1/2") = op(OpName("S+"), t)

ITensors.op(::OpName"S-", ::SiteType"S=1/2") = [
  0 0
  1 0
]

ITensors.op(::OpName"S⁻", t::SiteType"S=1/2") = op(OpName("S-"), t)

ITensors.op(::OpName"Sminus", t::SiteType"S=1/2") = op(OpName("S-"), t)

ITensors.op(::OpName"X", ::SiteType"S=1/2") = [
  0 1
  1 0
]

ITensors.op(::OpName"Sx", ::SiteType"S=1/2") = [
  0.0 0.5
  0.5 0.0
]

ITensors.op(::OpName"Sˣ", t::SiteType"S=1/2") = op(OpName("Sx"), t)

ITensors.op(::OpName"iY", ::SiteType"S=1/2") = [
  0 1
  -1 0
]

ITensors.op(::OpName"iSy", ::SiteType"S=1/2") = [
  0.0 0.5
  -0.5 0.0
]

ITensors.op(::OpName"iSʸ", t::SiteType"S=1/2") = op(OpName("iSy"), t)

ITensors.op(::OpName"Y", ::SiteType"S=1/2") = [
  0.0 -1.0im
  1.0im 0.0
]

ITensors.op(::OpName"Sy", ::SiteType"S=1/2") = [
  0.0 -0.5im
  0.5im 0.0
]

ITensors.op(::OpName"Sʸ", t::SiteType"S=1/2") = op(OpName("Sy"), t)

ITensors.op(::OpName"S2", ::SiteType"S=1/2") = [
  0.75 0.0
  0.0 0.75
]

ITensors.op(::OpName"S²", t::SiteType"S=1/2") = op(OpName("S2"), t)

ITensors.op(::OpName"ProjUp", ::SiteType"S=1/2") = [
  1 0
  0 0
]

ITensors.op(::OpName"projUp", t::SiteType"S=1/2") = op(OpName("ProjUp"), t)

ITensors.op(::OpName"ProjDn", ::SiteType"S=1/2") = [
  0 0
  0 1
]

ITensors.op(::OpName"projDn", t::SiteType"S=1/2") = op(OpName("ProjDn"), t)

# Support the tag "SpinHalf" as equivalent to "S=1/2"

space(::SiteType"SpinHalf"; kwargs...) = space(SiteType("S=1/2"); kwargs...)

val(name::ValName, ::SiteType"SpinHalf") = val(name, SiteType("S=1/2"))

ITensors.state(name::StateName, ::SiteType"SpinHalf") = state(name, SiteType("S=1/2"))

function ITensors.op(o::OpName, ::SiteType"SpinHalf"; kwargs...)
  return op(o, SiteType("S=1/2"); kwargs...)
end

# Support the tag "S=½" as equivalent to "S=1/2"

space(::SiteType"S=½"; kwargs...) = space(SiteType("S=1/2"); kwargs...)

val(name::ValName, ::SiteType"S=½") = val(name, SiteType("S=1/2"))

ITensors.state(name::StateName, ::SiteType"S=½") = state(name, SiteType("S=1/2"))

ITensors.op(o::OpName, ::SiteType"S=½"; kwargs...) = op(o, SiteType("S=1/2"); kwargs...)


text:
using ScTenifold

using LinearAlgebra, Statistics, Distributions, Random

d=NegativeBinomial(20,0.98)
X=rand(d,100,2000)
lbszv=30

Y=copy(X)
Y[10,:]=Y[50,:]
Y[2,:]=Y[11,:]
Y[3,:]=Y[5,:]

X=X[:,vec(sum(X,dims=1).>lbszv)]
Y=Y[:,vec(sum(Y,dims=1).>lbszv)]

@show Threads.nthreads()

@time Z0,_=ScTenifold.tenrnet(X, donorm=true)
@time Z1,_=ScTenifold.tenrnet(Y, donorm=true)
Z0=0.5*(Z0+Z0')
Z1=0.5*(Z1+Z1')
@time d=ScTenifold.manialn(Z0,Z1)
fc,p,adjp=ScTenifold.drgenes(d)

@show [adjp[10] adjp[11] adjp[50] adjp[20]];
@show [adjp[70] adjp[31] adjp[55] adjp[26]];

using StatsPlots
x=rand(Chisq(1), length(fc))
qqplot(x, fc)

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
function get_celestial_body(obj::CenterName)
    return jcall(obj, "getCelestialBody", CelestialBody, ())
end

function get_celestial_body(obj::CenterName, arg0::CelestialBodies)
    return jcall(obj, "getCelestialBody", CelestialBody, (CelestialBodies,), arg0)
end

function value_of(::Type{CenterName}, arg0::JString)
    return jcall(CenterName, "valueOf", CenterName, (JString,), arg0)
end

function values(::Type{CenterName})
    return jcall(CenterName, "values", Vector{CenterName}, ())
end


text:
using Random, PGFPlotsX, KernelFunctions, Distributions, LaTeXStrings

X = reshape(collect(range(0.0,3.0,length=300)),:,1)

K = 5000

start = 25

sinactivation(x) = sqrt(2) * sin(x)
sincosactivation(x) = sin(x) + cos(x)

ffloor(x) = floor(x/π + 0.5)
_z(x) = (x - π * ffloor(x)) * (-1)^(ffloor(x))
triangeactivation(x) = π/(2*sqrt(2)) * _z(x)

_prelu(x) = π/4	* (_z(x) + _z(x + π/2))
prelu(x) = _prelu(x)

b = rand(-π:0.01:π, K)

f1 = sinactivation
f2 = sincosactivation

k1(x1, x2, w, f) = 1/K * sum( f(x1*w[k] + b[k]) * f(x2*w[k] + b[k]) for k in 1:K)
k2(x1, x2, w, f) = 1/K * sum( f(x1*w[k]) * f(x2*w[k]) for k in 1:K)

w₁ = rand(TDist(2 * 1/2), K)
k₁ = Matern12Kernel()
K₁ = kernelmatrix(k₁,X,obsdim=1)

w₂ = rand(TDist(2 * 3/2), K)
k₂ = Matern32Kernel()
K₂ = kernelmatrix(k₂,X,obsdim=1)

w₃ = rand(Normal(), K)
k₃ = SqExponentialKernel()
K₃ = kernelmatrix(k₃,X,obsdim=1)

z_sin₁ = k1.(Ref(X[start]), X, Ref(w₁), Ref(sinactivation))
z_sin₂ = k1.(Ref(X[start]), X, Ref(w₂), Ref(sinactivation))
z_sin₃ = k1.(Ref(X[start]), X, Ref(w₃), Ref(sinactivation))

z_sincos₁ = k2.(Ref(X[start]), X, Ref(w₁), Ref(sincosactivation))
z_sincos₂ = k2.(Ref(X[start]), X, Ref(w₂), Ref(sincosactivation))
z_sincos₃ = k2.(Ref(X[start]), X, Ref(w₃), Ref(sincosactivation))

z_tri₁ = k1.(Ref(X[start]), X, Ref(w₁), Ref(triangeactivation))
z_tri₂ = k1.(Ref(X[start]), X, Ref(w₂), Ref(triangeactivation))
z_tri₃ = k1.(Ref(X[start]), X, Ref(w₃), Ref(triangeactivation))

z_relu₁ = k1.(Ref(X[start]), X, Ref(w₁), Ref(prelu))
z_relu₂ = k1.(Ref(X[start]), X, Ref(w₂), Ref(prelu))
z_relu₃ = k1.(Ref(X[start]), X, Ref(w₃), Ref(prelu))

@pgf p1 = Axis(
    {
        no_markers,
        xlabel = raw"Input, $|x-x'|$",
        width = "15cm",
        height = "10cm",
        ymax = 1.2,
        no_markers,
        xticklabels={},
        yticklabels={},
        ymin = 0,
    },
    # Exponential
    PlotInc({color = "black", dashed}, Table(1:300, K₁[start,:])),
    LegendEntry("Exact"),

    PlotInc({color = "red", solid}, Table(1:300, vec(z_sin₁))),
    LegendEntry("Sinusiod"),

    PlotInc({color = "green", solid}, Table(1:300, vec(z_sincos₁))),
    LegendEntry("Sine Cosine"),

    # Matern 3/2
    PlotInc({color = "black", dashed}, Table(201:500, K₂[start,:])),
    PlotInc({color = "red", solid}, Table(201:500, vec(z_sin₂))),
    PlotInc({color = "green", solid}, Table(201:500, vec(z_sincos₂))),

    # Squared Exponential
    PlotInc({color = "black", dashed}, Table(401:700, K₃[start,:])),
    PlotInc({color = "red", solid}, Table(401:700, vec(z_sin₃))),
    PlotInc({color = "green", solid}, Table(401:700, vec(z_sincos₃))),

    [raw"\node at ",
     Coordinate(25, 1.1),
     raw"{Exponential};"],
     [raw"\node at ",
     Coordinate(225, 1.1),
     raw"{Mat\'ern-$\frac{3}{2}$};"],
     [raw"\node at ",
     Coordinate(425, 1.1),
     raw"{Squared Exponential};"]
)

@pgf p2 = Axis(
    {
        no_markers,
        xlabel = raw"Input, $|x-x'|$",
        width = "15cm",
        height = "10cm",
        ymax = 1.2,
        no_markers,
        xticklabels={},
        yticklabels={},
        ymin = 0,
    },
    # Exponential
    PlotInc({color = "black", dashed}, Table(1:300, K₁[start,:])),
    LegendEntry("Exact"),

    PlotInc({color = "blue", solid}, Table(1:300, vec(z_tri₁))),
    LegendEntry("Triangle Wave"),

    PlotInc({color = "orange", solid}, Table(1:300, vec(z_relu₁))),
    LegendEntry("Periodic ReLU"),

    # Matern 3/2
    PlotInc({color = "black", dashed}, Table(201:500, K₂[start,:])),
    PlotInc({color = "blue", solid}, Table(201:500, vec(z_tri₂))),
    PlotInc({color = "orange", solid}, Table(201:500, vec(z_relu₂))),

    # Squared Exponential
    PlotInc({color = "black", dashed}, Table(401:700, K₃[start,:])),
    PlotInc({color = "blue", solid}, Table(401:700, vec(z_tri₃))),
    PlotInc({color = "orange", solid}, Table(401:700, vec(z_relu₃))),

    [raw"\node at ",
     Coordinate(25, 1.1),
     raw"{Exponential};"],
     [raw"\node at ",
     Coordinate(225, 1.1),
     raw"{Mat\'ern-$\frac{3}{2}$};"],
     [raw"\node at ",
     Coordinate(425, 1.1),
     raw"{Squared Exponential};"]
)

@pgf gp = GroupPlot(
    {
        group_style = {
            group_size = "2 by 1",
            horizontal_sep = "0.2cm",
        },
        no_markers,
        xticklabels={},
        yticklabels={},
        ymin = 0,
    },
    p1,
    p2,
)

pgfsave("kernels_p1.tex", p1, include_preamble = false)
pgfsave("kernels_p2.tex", p2, include_preamble = false)

pgfsave("kernels.pdf", gp)
pgfsave("kernels.tex", gp)

x = range(-2*pi; stop = 2*pi, length = 500)

@pgf sin_plot = Axis(
    {
        ylabel = raw"$\sigma(x)$",
        no_markers,
        ymin = -2.5,
        ymax = 2.5,
        xticklabels={},
    },
    Plot({color = "red"},
    Table(x, sinactivation.(x))
    ), LegendEntry("Sinusiod")
    )
@pgf sincos_plot = Axis(
    {
        no_markers,
        ymin = -2.5,
        ymax = 2.5,
        xticklabels={},
    },
    PlotInc({color = "green"}, Table(x, sincosactivation.(x))), LegendEntry("Sine Cosine"),
    PlotInc({dashed, color = "red"}, Table(x, sinactivation.(x)))
    )
@pgf tri_plot = Axis(
    {
        no_markers,
        ymin = -2.5,
        ymax = 2.5,
        xticklabels={},
    },
    PlotInc({color = "blue"}, Table(x, triangeactivation.(x))), LegendEntry("Triangle Wave"))
@pgf relu_plot = Axis(
    {
        no_markers,
        ymin = -2.5,
        ymax = 2.5,
        xticklabels={},
    },
    PlotInc({color = "orange"}, Table(x, prelu.(x))), LegendEntry("Periodic ReLU"))

@pgf gp = GroupPlot(
    {
        group_style = {
            group_size = "4 by 1",
            yticklabels_at="edge left",
            xticklabels_at="edge bottom",
            horizontal_sep = "0.2cm",
        },
        no_markers,
        ymin = -2.5,
        ymax = 2.5,
        xticklabels={},
    },
    sin_plot,
    sincos_plot,
    tri_plot,
    relu_plot
)

pgfsave("sin_activation.tex", sin_plot, include_preamble = false)
pgfsave("sincosactivation.tex", sincos_plot, include_preamble = false)
pgfsave("triangle_activation.tex", tri_plot, include_preamble = false)
pgfsave("relu_activation.tex", relu_plot, include_preamble = false)

pgfsave("activations.pdf", gp)
pgfsave("activations.tex", gp)

max_stars_repo_name: AaltoML/PeriodicBNN

max_issues_repo_name: AaltoML/PeriodicBNN

max_forks_repo_name: AaltoML/PeriodicBNN


text:
"""
    DCC{p, q, VS<:UnivariateVolatilitySpec, T<:AbstractFloat, d} <: MultivariateVolatilitySpec{T, d}
"""
struct DCC{p, q, VS<:UnivariateVolatilitySpec, T<:AbstractFloat, d} <: MultivariateVolatilitySpec{T, d}
    R::Matrix{T}
    coefs::Vector{T}
    univariatespecs::Vector{VS}
    method::Symbol
    function DCC{p, q, VS, T, d}(R::Array{T}, coefs::Vector{T}, univariatespecs:: Vector{VS}, method::Symbol) where {p, q, T, VS<:UnivariateVolatilitySpec, d}
        length(coefs) == nparams(DCC{p, q})  || throw(NumParamError(nparams(DCC{p, q}), length(coefs)))
        @assert d == length(univariatespecs)
        @assert method==:twostep || method==:largescale
        new{p, q, VS, T, d}(R, coefs, univariatespecs, method)
    end
end

"""
    CCC{VS<:UnivariateVolatilitySpec, T<:AbstractFloat, d} <: MultivariateVolatilitySpec{T, d}
---
    DCC(Qbar, coefs, univariatespecs; method=:largescale)
Construct a CCC specification with the given parameters. `coefs` must be passed
as a length-zero Vector of the same element type as Qbar.
"""
const CCC = DCC{0, 0}

"""
    DCC{p, q}(Qbar, coefs, univariatespecs; method=:largescale)
Construct a DCC(p, q) specification with the given parameters.
"""
DCC{p, q}(R::Matrix{T}, coefs::Vector{T}, univariatespecs::Vector{VS}; method::Symbol=:largescale) where {p, q, T, VS<:UnivariateVolatilitySpec{T}} = DCC{p, q, VS, T, length(univariatespecs)}(R, coefs, univariatespecs, method)

nparams(::Type{DCC{p, q}}) where {p, q} = p+q

nparams(::Type{DCC{p, q, VS, T, d}}) where {p, q, VS, T, d}= p + q + d * nparams(VS)

# strange dispatch behavior. to me these methods look the same, but they aren't.

# this matches ARCHModels.presample(DCC{1,1,TGARCH{0,1,1,Float64}})
presample(::Type{DCC{p, q, VS}}) where {p, q, VS} = max(p, q, presample(VS))

# this matches ARCHModels.presample(DCC{1,1,TGARCH{0,1,1,Float64},Float64,2})
presample(::Type{DCC{p, q, VS, T, d}}) where {p, q, VS, T, d} = max(p, q, presample(VS))

fit(::Type{<:DCC}, data::Matrix{T}; meanspec=Intercept{T}, method=:largescale, algorithm=BFGS(), autodiff=:forward, kwargs...) where {T} = fit(DCC{1, 1}, data; meanspec=meanspec, method=method, algorithm=algorithm, autodiff=autodiff, kwargs...)

fit(DCCspec::Type{<:DCC{p, q}}, data::Matrix{T}; meanspec=Intercept{T},  method=:largescale, algorithm=BFGS(), autodiff=:forward, kwargs...) where {p, q, T} = fit(DCC{p, q, GARCH{1, 1}}, data; meanspec=meanspec, method=method, algorithm=algorithm, autodiff=autodiff, kwargs...)

"""
    fit(DCCspec::Type{<:DCC{p, q, VS<:UnivariateVolatilitySpec}}, data::Matrix;
        method=:largescale,  dist=MultivariateStdNormal, meanspec=Intercept,
        algorithm=BFGS(), autodiff=:forward, kwargs...)

Fit the DCC model specified by `DCCspec` to `data`. If `p` and `q` or `VS` are
unspecified, then these default to 1, 1, and `GARCH{1, 1}`.

# Keyword arguments:
- `method`: one of `:largescale` or `twostep`
- `dist`: the error distribution.
- `meanspec`: the mean specification, as a type.
- `algorithm, autodiff, kwargs, ...`: passed on to the optimizer.

# Example: DCC{1, 1, GARCH{1, 1}} model:
```jldoctest
julia> fit(DCC, DOW29)

29-dimensional DCC{1, 1} - GARCH{1, 1} - Intercept{Float64} specification, T=2785.

DCC parameters, estimated by largescale procedure:
────────────────────
       β₁         α₁
────────────────────
  0.88762  0.0568001
────────────────────

Calculating standard errors is expensive. To show them, use
`show(IOContext(stdout, :se=>true), <model>)`
```
"""
function fit(DCCspec::Type{<:DCC{p, q, VS}}, data::Matrix{T}; meanspec=Intercept{T}, method=:largescale, algorithm=BFGS(), autodiff=:forward, dist::Type{<:MultivariateStandardizedDistribution}=MultivariateStdNormal{T}) where {p, q, VS<: UnivariateVolatilitySpec, T, d}
    n, dim = size(data)
    resids = similar(data)
    if n<12 && method == :largescale
        error("largescale method requires n>11.")
    end

    m = fit(VS, data[:, 1], meanspec=meanspec)
    resids[:, 1] = residuals(m)
    univariatespecs = Vector{typeof(m)}(undef, dim)
    univariatespecs[1] = m
    Threads.@threads for i = 2:dim
        m = fit(VS, data[:, i], meanspec=meanspec)
        univariatespecs[i] = m
        resids[:, i] = residuals(m)
    end
    method == :largescale ? Σ = analytical_shrinkage(resids) : Σ = cov(resids)
    R = to_corr(Σ)
    x0 = zeros(T, p+q)
    if p+q>0
        x0[1:p] .= 0.9/p
        x0[p+1:end] .= 0.05/q
        if method == :twostep
            obj = LL2step
        elseif method==:largescale
            obj = LL2step_pairs
        else
            error("No method :$method.")
        end
        f = x -> obj(DCCspec, x, R, resids)
        x = optimize(x->-sum(f(x)), x0, algorithm, autodiff=autodiff).minimizer
    else # CCC
        x = x0
    end
    return MultivariateARCHModel(DCC{p, q}(R, x, getproperty.(univariatespecs, :spec); method=method), data; dist=MultivariateStdNormal{T, dim}(), meanspec=getproperty.(univariatespecs, :meanspec), fitted=true)
end

#LC(Θ_hat, ϕ) in Engle (2002)
@inline function LL2step!(Rt::Array{Array{T, 2}, 1}, DCCspec::Type{<:DCC{p, q}}, coef::Array{T}, R, resids::Array{T2}) where {T, T2, p, q}
    n, dims = size(resids)
    LL = zeros(T, n)
    all(0 .< coef .< 1) || (fill!(LL, T(-Inf)); return LL)
    abs(sum(coef))>1 && (fill!(LL, T(-Inf)); return LL)
    f = 1 - sum(coef)
    e = @view resids[1, :]
    R = Symmetric(R)
    Rt[1:max(p,q)] .= [R for _ in 1:max(p,q)]
    RD5 = Diagonal(zeros(T, dims))
    C = cholesky(Rt[1]).L
    u = inv(C) * e
    for t=1:n
        if t > max(p, q)
            Rt[t] .= R * f
            for i = 1:p
                Rt[t] .+=  coef[i] * Rt[t-i]
            end
            for i = 1:q
                Rt[t] .+= coef[p+i]  * resids[t-i, :]*resids[t-i, :]'
            end
            RD5 .= inv(sqrt(Diagonal(Rt[t])))
            Rt[t] .= Symmetric(RD5 * Rt[t] * RD5)
            C .= cholesky(Rt[t]).L
        end
        e = @view resids[t, :]
        u .= inv(C) * e
        L = (dot(e, e) - dot(u, u))/2-logdet(C)
        LL[t] = L
    end
    LL
end

function LL2step(DCCspec::Type{<:DCC{p, q}}, coef::Array{T}, R, resids::Array{T2}) where {T, T2, p, q}
    n, dims = size(resids)
    Rt = [zeros(T, dims, dims) for _ in 1:n]
    LL2step!(Rt, DCCspec, coef, R, resids)
end

#same as LL2step, except for init type
function LL2step2(DCCspec::Type{<:DCC{p, q}}, coef::Array{T2}, R, resids::Array{T}) where {T, T2, p, q}
    n, dims = size(resids)
    LL = zeros(T, n)
    all(0 .< coef .< 1) || (fill!(LL, T(-Inf)); return LL)
    abs(sum(coef))>1 && (fill!(LL, T(-Inf)); return LL)
    f = 1 - sum(coef)

    e = @view resids[1, :]
    Rt = [zeros(T, dims, dims) for _ in 1:n]
    R = Symmetric(R)
    Rt[1:max(p,q)] .= [R for _ in 1:max(p,q)]
    RD5 = Diagonal(zeros(T, dims))
    C = cholesky(Rt[1]).L
    u = inv(C) * e
    for t = 1:n
        if t > max(p, q)
            Rt[t] .= R * f
            for i = 1:p
                Rt[t] .+=  coef[i] * Rt[t-i]
            end
            for i = 1:q
                Rt[t] .+= coef[p+i]  * resids[t-i, :]*resids[t-i, :]'
            end
            RD5 .= inv(sqrt(Diagonal(Rt[t])))
            Rt[t] .= Symmetric(RD5 * Rt[t] * RD5)
            C .= cholesky(Rt[t]).L
        end
        e = @view resids[t, :]
        u .= inv(C) * e
        L = (dot(e, e) - dot(u, u))/2-logdet(C)
        LL[t] = L
    end
    LL
end

#doall toggles whether to return all individual likelihood contributions
function LL2step_pairs(DCCspec::Type{<:DCC{p, q}}, coef::Array{T}, R, resids::Array{T2}, doall=false) where {T, T2, p, q}
    n, dims = size(resids)
    len = doall ? n : 1
    LL = zeros(T, len, dims)
    #Threads.@threads
    for k = 1:dims-1
        thell = ll(DCCspec, coef, R[k, k+1], resids[:, k:k+1], doall)
        if doall
            LL[:, k] .= thell
        else
            LL[1, k:k] .= thell
        end
    end
    sum(LL, dims=2)
end

@inline function ll(DCCspec::Type{<:DCC{p, q}}, coef::Array{T}, rho, resids, doall=false) where {T, p, q}
    all(0 .< coef .< 1) || return T(-Inf)
    abs(sum(coef)) < 1 || return T(-Inf)
    n, dims = size(resids)
    f = 1 - sum(coef)
    len = doall ? n : 1
    LL = zeros(T, len)

    rt = zeros(T, n) # should switch this to circbuff for speed
    s1 = T(1)
    s2 = T(1)
    fill!(rt, rho)
    @inbounds for t=1:n
        if t > max(p, q)
            s1 = T(1)
            s2 = T(1)
            rt[t] = rho * f
            for i = 1:q
                s1 += coef[p+i] * (resids[t-i, 1]^2 - 1)
                s2 += coef[p+i] * (resids[t-i, 2]^2 - 1)
                rt[t] += coef[p+i] * resids[t-i, 1] * resids[t-i, 2]
            end
            for i = 1:p
                rt[t] += coef[i] * rt[t-i]
            end
            rt[t] = rt[t] / sqrt(s1 * s2)
        end
        e1 = resids[t, 1]
        e2 = resids[t, 2]
        r2 = rt[t]^2
        d = 1 - r2

        L = (((e1*e1 + e2*e2) * r2 - 2 * rt[t] *e1 * e2) / d + log(d)) / 2

        if doall
            LL[t] = -L
        else
            LL[1] -= L
        end
     end
    LL
end

function stderror(am::MultivariateARCHModel{T, d, MVS}) where {T, d, p, q, VS, MVS<:DCC{p, q, VS}}
    n, dim = size(am.data)
    r = p + q
    resids = similar(am.data)
    nunivariateparams = nparams(VS) + nparams(typeof(am.meanspec[1]))
    np = r + dim * nunivariateparams
    coefs = coef(am)
    Htt = zeros(np - r, np - r)
    dt = zeros(n, np - r)
    stderrors = zeros(np)
    Threads.@threads for i = 1:dim
        m = UnivariateARCHModel(am.spec.univariatespecs[i], am.data[:, i]; meanspec=am.meanspec[i], fitted=true)
        resids[:, i] = residuals(m)
        w=1+(i-1)*nunivariateparams:1+i*nunivariateparams-1
        Htt[w, w] .= -informationmatrix(m, expected=false)/nobs(m) # is the /nobs correct here?
        dt[:, w] = scores(m)
        stderrors[r .+ w] = stderror(m)
    end
    if p + q > 0
        if am.spec.method == :twostep
            f = x -> LL2step(MVS, x, am.spec.R, resids)
            Hpp = ForwardDiff.hessian(x->sum(f(x)), coefs[1:r])/n
            dp = ForwardDiff.jacobian(f, coefs[1:r])

            # g = x -> sum(LL2step_full(x, R, data, p, q))
            # Hpt = FiniteDiff.finite_difference_hessian(g, coefs)[1:2, 3:end]/n
            # use finite differences instead, because we don't need the whole
            # Hessian, and I couldn't figure out how to do this with ForwardDiff
            g = (x, y) -> sum(LL2step_full(MVS, VS, am.meanspec, x, y, am.spec.R, am.data))
            dg = x -> ForwardDiff.gradient(y->g(x, y), coefs[1+r:end])/n
            h = 1e-7
            Hpt = zeros(p+q, dim * nunivariateparams)
            for j=1:p+q
                dg0 = dg(coefs[1:r])
                xp = copy(coefs[1:r]); xp[j] += h
                ddg = (dg(xp)-dg0)/h
                Hpt[j, :] = ddg
            end
            A = dp-(Hpt*inv(Htt)*dt')'
            C = inv(Hpp)*A'*A*inv(Hpp)/n^2
            stderrors[1:r] = sqrt.(diag(C))
        elseif am.spec.method==:largescale
            g = x -> LL2step_pairs(MVS, x, am.spec.R, resids, true)
            sc = ForwardDiff.jacobian(g, coefs[1:r])
            I = sc'*sc/n/dim
            h = x-> LL2step_pairs_full(MVS, VS, am.meanspec, x, am.spec.R, am.data)
            H = ForwardDiff.hessian(x->sum(h(x)), coefs)/n/dim
            #J = H[1:r, 1:r] - H[1:r, r+1:end] * inv(H[1+r:end, 1+r:end]) * H[1:r, 1+r:end]'
            #std = sqrt.(diag(inv(J)*I*inv(J))/n) # from the 2014 version of the paper
            as = hcat(dt, sc) # all scores
            Sig = as'*as/n/dim
            Jnt = hcat(inv(H[1:r, 1:r])*H[1:r, 1+r:end]*inv(Htt), -inv(H[1:r, 1:r]))
            stderrors[1:r] .= sqrt.(diag(Jnt*Sig*Jnt'/n)) # from the 2018 version
        end
    end
    return stderrors
end

#LC(Θ, ϕ) in Engle (2002)
function LL2step_full(DCCspec::Type{<:DCC{p, q}}, VS, meanspec, dcccoef::Array{T}, garchcoef::Array{T2}, R, data) where {T, T2, p, q}
    n, dims = size(data)
    resids = Array{T2}(undef, size(data))
    nunivariateparams = nparams(VS) + nparams(typeof(meanspec[1]))
    for i = 1:dims
        params = garchcoef[1+(i-1)*nunivariateparams:1+i*nunivariateparams-1]
        ht = T2[]
        lht = T2[]
        zt = T2[]
        at = T2[]
        loglik!(ht, lht, zt, at, VS, StdNormal{Float64}, meanspec[i], data[:, i], params)
        resids[:, i] = zt
    end
    LL2step2(DCCspec, dcccoef, R, resids)
end

#LC(Θ, ϕ) in Engle (2002). not actually the full log-likelihood
#this method only needed for Hpt when using ForwardDiff
# function LL2step_full(coef::Array{T}, R, data, p, q) where {T}
#     n, dims = size(data)
#     resids = Array{T}(undef, size(data))
#     for i = 1:dims
#         params = coef[3+(i-1)*nparams(GARCH{1, 1}):3+i*nparams(GARCH{1, 1})-1]
#         ht = T[]
#         lht = T[]
#         zt = T[]
#         at = T[]
#         loglik!(ht, lht, zt, at, GARCH{1, 1, Float64}, StdNormal{Float64}, NoIntercept(), data[:, i], params)
#         resids[:, i] = zt
#     end
#     LL2step(coef[1:2], R, resids, p, q)
# end

function LL2step_pairs_full(DCCspec::Type{<:DCC{p, q}}, VS::Type{<:UnivariateVolatilitySpec}, meanspec, coef::Array{T}, R, data) where {T, T2, p, q}
    dcccoef = coef[1:p+q]
    garchcoef = coef[p+q+1:end]
    n, dims = size(data)
    resids = Array{T}(undef, size(data))
    nunivariateparams = nparams(VS) + nparams(typeof(meanspec[1]))
    for i = 1:dims
        params = garchcoef[1+(i-1)*nunivariateparams:1+i*nunivariateparams-1]
        ht = T[]
        lht = T[]
        zt = T[]
        at = T[]
        loglik!(ht, lht, zt, at, VS, StdNormal{Float64}, meanspec[i], data[:, i], params)
        resids[:, i] = zt
    end
    LL2step_pairs(DCCspec::Type{<:DCC{p, q}}, dcccoef, R, resids)
end

function coefnames(::Type{<:DCC{p, q}}) where {p, q}
    names = Array{String, 1}(undef, p + q)
    names[1:p] .= (i -> "β"*subscript(i)).([1:p...])
    names[p+1:p+q] .= (i -> "α"*subscript(i)).([1:q...])
    return names
end

function coef(spec::DCC{p, q, VS, T, d})  where {p, q, VS, T, d}
    vcat(spec.coefs, [spec.univariatespecs[i].coefs for i in 1:d]...)
end

function coef(am::MultivariateARCHModel{T, d, MVS}) where {T, d, MVS<:DCC}
    vcat(am.spec.coefs, [vcat(am.spec.univariatespecs[i].coefs, am.meanspec[i].coefs) for i in 1:d]...)
end

function coefnames(am::MultivariateARCHModel{T, d, MVS}) where {T, d, p, q, VS, MVS<:DCC{p, q, VS}}
    nunivariateparams = nparams(VS) + nparams(typeof(am.meanspec[1]))
    names = Array{String, 1}(undef, p + q + d * nunivariateparams)
    names[1:p+q] .= coefnames(MVS)
    for i = 1:d
            names[p + q + 1 + (i-1) * nunivariateparams : p + q +  i * nunivariateparams] = vcat(coefnames(VS) .* subscript(i), coefnames(am.meanspec[i]) .* subscript(i))
    end
    return names
end

modname(::Type{DCC{p, q, VS, T, d}})  where {p, q, VS, T, d} = "DCC{$p, $q, $(modname(VS))}"

function show(io::IO, am::MultivariateARCHModel{T, d, MVS}) where {T, d, p, q, VS, MVS<:DCC{p, q, VS}}
    r = p + q
    cc = coef(am)[1:r]
    println(io, "\n", "$d-dimensional DCC{$p, $q} - $(modname(VS)) - $(modname(typeof(am.meanspec[1]))) specification, T=", size(am.data)[1], ".\n")
    if isfitted(am) && (:se=>true) in io
        se = stderror(am)[1:r]
        z = cc ./ se
        if p + q >0
            println(io, "DCC parameters, estimated by $(am.spec.method) procedure:", "\n",
    	            CoefTable(hcat(cc, se, z, 2.0 * normccdf.(abs.(z))),
    	                      ["Estimate", "Std.Error", "z value", "Pr(>|z|)"],
    	                      coefnames(MVS), 4
    	                      )
    	            )
        end
    else
        if p + q > 0
            println(io, "DCC parameters", isfitted(am) ? ", estimated by $(am.spec.method) procedure:" : "", "\n",
    	            CoefTable(cc, coefnames(MVS), [""])
    	            )
            if isfitted(am)
                println(io, "\n","""Calculating standard errors is expensive. To show them, use
                 `show(IOContext(stdout, :se=>true), <model>)`""")
            end
        end
    end
end

"""
    correlations(am::MultivariateARCHModel)
Return the estimated conditional correlation matrices.
"""
function correlations(am::MultivariateARCHModel{T, d, MVS}) where {T, d, MVS<:DCC}
    resids = residuals(am; decorrelated=false)
    n, dims = size(resids)
    Rt = [zeros(T, dims, dims) for _ in 1:n]
    LL2step!(Rt, MVS, am.spec.coefs, am.spec.R, resids)
    return Rt
end

"""
    covariances(am::MultivariateARCHModel)
Return the estimated conditional covariance matrices.
"""
function covariances(am::MultivariateARCHModel{T, d, MVS}) where {T, d, MVS<:DCC}
    n, dims = size(am.data)
    Rt = correlations(am)
    for i = 1:d
        v = volatilities(UnivariateARCHModel(am.spec.univariatespecs[i], am.data[:, i]; meanspec=am.meanspec[i], fitted=true))
        @inbounds for t = 1:n # this is ugly, but I couldn't figure out how to do this w/ broadcasting
            Rt[t][i, :] *= v[t]
            Rt[t][:, i] *= v[t]
        end
    end
    return Rt
end

"""
    residuals(am::MultivariateARCHModel; standardized = true, decorrelated = true)
Return the residuals.
"""
function residuals(am::MultivariateARCHModel{T, d, MVS}; standardized = true, decorrelated = true) where {T, d, MVS<:DCC}
    n, dims = size(am.data)
    resids = similar(am.data)
    Threads.@threads for i = 1:dims
        m = UnivariateARCHModel(am.spec.univariatespecs[i], am.data[:, i]; meanspec=am.meanspec[i], fitted=true)
        resids[:, i] = residuals(m; standardized=standardized)
    end
    if decorrelated
        Rt = standardized ? correlations(am) : covariances(am)
        @inbounds for t = 1:n
            resids[t, :] = inv(cholesky(Rt[t]; check=false).L) * resids[t, :]
        end
    end
    return resids
end

#this assumes Ht, Rt, zt, and at are circularbuffers or vectors of arrays
Base.@propagate_inbounds @inline function update!(Ht, Rt, H, R, zt, at, MVS::Type{DCC{p, q, VS, T, d}}, coefs) where {p, q, VS, T, d}
    nvolaparams = nparams(VS)
    h5s = zeros(T, d)
    for i = 1:d
        ht = getindex.(Ht, i, i)
        lht = log.(ht)
        update!(ht, lht, getindex.(zt, i), getindex.(at, i), VS, coefs[p + q + 1 + (i-1) * nvolaparams : p + q + i * nvolaparams])
        h5s[i] = sqrt(ht[end])
    end
    Rtemp = R * (1-sum(coefs[1:p+q]))
    for i = 1:p
        Rtemp .+=  coefs[i] * Rt[end-i+1]
    end
    for i = 1:q
        Rtemp .+= coefs[p+i]  * zt[end-i+1] * zt[end-i+1]'
    end
    push!(Rt, to_corr(Rtemp))
    H5 = diagm(0 => h5s)
    push!(Ht,  H5 * Rt[end] * H5)
end

function uncond(spec::DCC{p, q, VS, T, d}) where {p, q, VS, T, d}
    h = uncond.(typeof.(spec.univariatespecs), getproperty.(spec.univariatespecs, :coefs))
    D = diagm(0 => sqrt.(h))
    return D * spec.R * D
end


text:
using Quadprog
using Base.Test

# simple test example
(x, fx, status) = quadprog([3 1; 1 4], [5, 6])
@test norm(x  - [-1.27273,-1.18182]) < 1e-5
@test norm(fx - -6.7272727272727275) < 1e-5
@test status == 1


text:
push!(LOAD_PATH, ".")
using Points

using Plots

function randomwalk(nsteps::Integer)
    ps = zeros(Point2D, nsteps+1)
    for i in 1:nsteps
        ps[i+1] = update(ps[i])
    end
    return ps
end

ps = randomwalk(1001)
xs = map(p->p.x, ps)
ys = map(p->p.y, ps)

pl = plot(xs, ys, xaxis="x", yaxis="y", color=:black, legend=false)
savefig(pl, "rw_2d.pdf")

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 83b65e34-e6d8-44da-947b-fb66ddee5872
begin
	using Revise
	using Gaston
	using PlutoUI
	include("../src/DirectionFinding.jl")
end

# ╔═╡ 9adcf5e6-751a-43e1-b2df-41f9419cd0a5
println("Gaston version: $(Gaston.GASTON_VERSION)")

# ╔═╡ 1bfb7099-2545-4be7-9c6f-af382905abe4
println("Gaston version: $(Gaston.GASTON_VERSION)")

# ╔═╡ 3d9ebfcf-2016-473b-a21a-fc10686de6e2
md"""
Rx = $(@bind Rx Slider(-10:10, show_value=true))

Ry = $(@bind Ry Slider(-10:10, show_value=true))

Ax = $(@bind Ax Slider(-10:10, show_value=true))

Ay = $(@bind Ay Slider(-10:10, show_value=true))

θ = $(@bind θ Slider(0:0.05:2pi, show_value=true))
"""

# ╔═╡ 27b94bd0-6f8a-48d5-a57b-e4531e7e4d3f
a = DirectionFinding.Antenna(DirectionFinding.Position(Ax, Ay), x->1)

# ╔═╡ fd4935cb-c229-44df-8591-356b57ca4b34
r = DirectionFinding.Position(Rx, Ry)

# ╔═╡ 773086fe-6650-40ff-9983-fedf1595d150
P = DirectionFinding.projection(a, r, θ)

# ╔═╡ 33401b43-c0e0-4d9f-977e-9e7cc93326dd
DirectionFinding.phaseshift(a, r, θ, 300e6/2e9)

# ╔═╡ 50623d25-faab-497c-9425-00576a524584
begin
	rx = r.x; ry = r.y # reference position
	ax = a.position.x; ay = a.position.y #antenna position
	# plot source direction
	m = tan(θ); y(x) = m*(x-rx)+ry;
	plot(-10:0.1:10, y.(-10:0.1:10), w="lines lc 'black'",
	     Axes(xrange=(-11,11), yrange=(-11,11), grid=:on))
	plot!([10cos(θ)],[10sin(θ)],w="points pt 'S'")
	# plot wavefront
	mw = tan(θ+π/2); yw(x) = mw*(x-rx)+ry;
	plot!(-10:0.1:10, yw.(-10:0.1:10), w="lines lc 'red' dt 3")
	# plot array reference position
	plot!([rx], [ry], w="points pt 'R'")
	# plot antenna position
	ps = DirectionFinding.phaseshift(a, r, θ, 300e6/2e9)
	if ps > 0
		plot!([ax], [ay], w="points pt 'A' lc 'blue'")
	else
		plot!([ax], [ay], w="points pt 'A' lc 'green'")
	end
	# plot projection
	px = P.x; py = P.y;
	plot!([px], [py], w="points pt 'P'")
	# plot line from source to reference position
	plot!([ax, px], [ay, py], w="lines")
end

# ╔═╡ b324892f-9c39-4dd3-994e-e393d7d3ba64
begin
	if θ ≈ 0
		if Ax > Rx
			"in front"
		else
			"behind"
		end
	elseif θ ≈ π
		if Ax < Rx
			"in front"
		else
			"behind"
		end
	else
		mm = tan(θ+π/2);
		yy(x) = mm*(x-rx)+ry
		th = yy(ax)
		if θ <= π
			if ay < th
				(ay, th, "behind")
			else
				(ay, th, "in front")
			end
		else
			if ay > th
				(ay, th, "behind")
			else
				(ay, th, "in front")
			end
		end
	end
end

# ╔═╡ a266135f-0767-47b8-951f-cc8f785ed5fd
# check orthogonality (the plot may be deceiving)
(rx-px)*(ax-px)+(ry-py)*(ay-py)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Gaston = "4b11ee91-296f-5714-9832-002c20994614"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

[compat]
Gaston = "~1.0.4"
PlutoUI = "~0.7.22"
Revise = "~3.2.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "abb72771fd8895a7ebd83d5632dc4b989b022b5b"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.2"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Gaston]]
deps = ["ColorSchemes", "DelimitedFiles", "Random"]
git-tree-sha1 = "ef62952980d19c98d00bd44d2266a98f8e9c7178"
uuid = "4b11ee91-296f-5714-9832-002c20994614"
version = "1.0.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "565564f615ba8c4e4f40f5d29784aa50a8f7bbaf"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.22"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "6f2bc8f1a444f93c52163e6f82270877224632d0"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.2.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═83b65e34-e6d8-44da-947b-fb66ddee5872
# ╠═9adcf5e6-751a-43e1-b2df-41f9419cd0a5
# ╟─1bfb7099-2545-4be7-9c6f-af382905abe4
# ╟─3d9ebfcf-2016-473b-a21a-fc10686de6e2
# ╟─27b94bd0-6f8a-48d5-a57b-e4531e7e4d3f
# ╟─fd4935cb-c229-44df-8591-356b57ca4b34
# ╟─773086fe-6650-40ff-9983-fedf1595d150
# ╠═33401b43-c0e0-4d9f-977e-9e7cc93326dd
# ╠═b324892f-9c39-4dd3-994e-e393d7d3ba64
# ╠═50623d25-faab-497c-9425-00576a524584
# ╠═a266135f-0767-47b8-951f-cc8f785ed5fd
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
### Congestion management without renewable energy sources
function CongestionManagement(time,
							  nodes::Nodes,
							  lines::Lines,
							  price,
							  P_opt,
							  powerplants::PowerPlants,
						      )

	T = time
	G = powerplants.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	lines_pmax = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		lines_pmax[j, k] = lines.pmax[id]
		lines_pmax[k, j] = lines.pmax[id]
	end

	# Create the B_prime matrix in one step
	B_prime = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		X = lines.reactance[id]				# reactance

		B_prime[j, k] = 1/X					# Fill YBUS with susceptance
		B_prime[k, j] = 1/X					# Symmetric matrix
	end

	for j in 1:size(B_prime)[1]
		B_prime[j, j] = -sum(B_prime[j, :])
	end

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	K = nodes.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_inj[T, J])
	@variable(CM_mod, P_flow[T, K, J])

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Constraints
	# Market clearing/power balance
	@constraint(CM_mod, MarketClearing[t = T],
	    sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in G) == nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g]/Sbase);

	# Ramp-up limit
	@constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
		(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
			(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
			powerplants.rup[g]/Sbase)

	# Ramp-down limit
	@constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
		-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
			(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
			powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) -	# all generation units connected to a node
			nodes.load[j][t]/Sbase == P_inj[t, j]);

	# Power injection angle
	@constraint(CM_mod, PowerInjectionAng[j = J, t = T],
		sum(B_prime[j, k] * (Θ[t, j] - Θ[t, k]) for k in K)
			== P_inj[t, j]);

	@constraint(CM_mod, LinePowerFlow[k = K, j = J, t = T; j != k],
		B_prime[j, k] * (Θ[t, j] - Θ[t, k]) == P_flow[t, j, k]);

	@constraint(CM_mod, LinePowerFlowMax[k = K, j = J, t = T],
		P_flow[t, j, k] <= lines_pmax[j, k]/Sbase);

	@constraint(CM_mod, LinePowerFlowMin[k = K, j = J, t = T],
		P_flow[t, j, k] >= -lines_pmax[j, k]/Sbase);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	P_inj_opt = JuMP.value.(P_inj)*Sbase

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
	       Θ_opt,
	       P_inj_opt,
	       P_flow_opt
	       )
end

### Congestion management with renewable energy sources
function CongestionManagement(time,
							  nodes::Nodes,
							  lines::Lines,
							  price,
							  P_opt,
						      powerplants::PowerPlants,
							  renewables::Renewables,
							  P_R_opt,
						      )

	T = time
	G = powerplants.unit
	R = renewables.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	lines_pmax = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		lines_pmax[j, k] = lines.pmax[id]
		lines_pmax[k, j] = lines.pmax[id]
	end

	# Create the B_prime matrix in one step
	B_prime = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		X = lines.reactance[id]				# reactance

		B_prime[j, k] = 1/X					# Fill YBUS with susceptance
		B_prime[k, j] = 1/X					# Symmetric matrix
	end

	for j in 1:size(B_prime)[1]
		B_prime[j, j] = -sum(B_prime[j, :])
	end

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	K = nodes.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_inj[T, J])
	@variable(CM_mod, P_flow[T, K, J])

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Constraints
	# Market clearing/power balance
	@constraint(CM_mod, MarketClearing[t = T],
	    sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in G) +
			sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in R)
				== nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.generation[r][t]/Sbase);

	# Ramp-up limit
	@constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
		(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
			(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
			powerplants.rup[g]/Sbase)

	# Ramp-down limit
	@constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
		-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
			(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
			powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) -
				nodes.load[j][t]/Sbase == P_inj[t, j]);

	# Power injection angle
	@constraint(CM_mod, PowerInjectionAng[j = J, t = T],
		sum(B_prime[j, k] * (Θ[t, j] - Θ[t, k]) for k in K)
			== P_inj[t, j]);

	@constraint(CM_mod, LinePowerFlow[k = K, j = J, t = T; j != k],
		B_prime[j, k] * (Θ[t, j] - Θ[t, k]) == P_flow[t, j, k]);

	@constraint(CM_mod, LinePowerFlowMax[k = K, j = J, t = T],
		P_flow[t, j, k] <= lines_pmax[j, k]/Sbase);

	@constraint(CM_mod, LinePowerFlowMin[k = K, j = J, t = T],
		P_flow[t, j, k] >= -lines_pmax[j, k]/Sbase);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(renewables.mc[r] * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - renewables.mc[r]) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	P_inj_opt = JuMP.value.(P_inj)*Sbase
	price = JuMP.dual.(MarketClearing)

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_inj_opt,
	       P_flow_opt,
		   price
	       )
end


text:
@concrete struct AsReal <: Real
    val
end
AsReal{T}(val) where {T} = AsReal(T(val))
# AsReal{T}(val::AsReal{T}) where {T} = AsReal(T(val.val))

const ℝ = AsReal

as_real(x::Real) = x
as_real(x) = AsReal(x)

for op in [+, -, inv, zero, sign, sqrt, deg2rad, rad2deg, abs, abs2, sin, cos, tan, sind, cosd, tand, sinh, cosh, tanh, asin, acos, atan, asind, acosd, atand, asinh, acosh, atanh, exp, exp2, exp10, log, log2, log10]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(x::ℝ) = ℝ($m.$f(x.val))
    end)
end

for op in [signbit, !]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(x::ℝ) = $m.$f(x.val)
    end)
end

for op in [+, -, *, /, //, ^, atan, atand, hypot, max, min,]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(a::ℝ, b::Number) = ℝ($m.$f(a.val, b))
        $m.$f(a::Number, b::ℝ) = ℝ($m.$f(a, b.val))
        $m.$f(a::ℝ, b::ℝ) = ℝ($m.$f(a.val, b.val))
    end)
end

for op in [+, -]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(a::ℝ, b::Complex) = Complex($m.$f(a, b.re), b.im)
        $m.$f(a::Complex, b::ℝ) = Complex($m.$f(a.re, b), a.im)
        $m.$f(a::ℝ, b::Complex{Bool}) = Complex($m.$f(a, b.re), b.im)
        $m.$f(a::Complex{Bool}, b::ℝ) = Complex($m.$f(a.re, b), a.im)
    end)
end

for op in [*, /, //]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(a::ℝ, b::Complex) = Complex($m.$f(a, b.re), $m.$f(a, b.im))
        $m.$f(a::Complex, b::ℝ) = Complex($m.$f(a.re, b), $m.$f(a.im, b))
        $m.$f(a::ℝ, b::Complex{Bool}) = Complex($m.$f(a, b.re), $m.$f(a, b.im))
        $m.$f(a::Complex{Bool}, b::ℝ) = Complex($m.$f(a.re, b), $m.$f(a.im, b))
    end)
end

for op in [==, <, >]
    f = nameof(op)
    m = Base.parentmodule(op)
    eval(quote
        $m.$f(a::ℝ, b::Number) = $m.$f(a.val, b)
        $m.$f(a::Number, b::ℝ) = $m.$f(a, b.val)
        $m.$f(a::ℝ, b::ℝ) = $m.$f(a.val, b.val)
    end)
end

Base.:(^)(a::ℝ, b::Integer) = ℝ(a.val^b)

Base.promote_rule(::Type{ℝ{T1}}, ::Type{T2}) where {T1,T2} = ℝ{promote_type(T1, T2)}
Base.promote_rule(::Type{ℝ{T1}}, ::Type{ℝ{T2}}) where {T1,T2} = ℝ{promote_type(T1, T2)}

Base.convert(::Type{ℝ{T1}}, x::T2) where {T1,T2<:Number} = promote_type(ℝ{T1}, T2)(x)
Base.convert(::Type{ℝ{T1}}, x::ℝ{T2}) where {T1,T2} = ℝ{promote_type(T1, T2)}(x.val)

Base.show(io::IO, x::ℝ) = print(io, "ℝ($(x.val))")
Base.show(io::IO, z::Complex{<:ℝ}) = print(io, "$(z.re) + $(z.im)im")


text:
using DSP
using FFTW
using IntervalSets
using SampleArrays
using Test
using Unitful
using Unitful: ns, ms, µs, s, Hz, kHz, MHz, GHz, THz, m

@testset "SampleArrays.jl" begin
    include("samplearray.jl")
    include("spectrumarray.jl")
    include("rfftspectrumarray.jl")
    include("utils.jl")
end


text:
using BinaryBuilder

name = "LAPACK"
version = v"3.9.0"

# Collection of sources required to build lapack
sources = [
    GitSource("https://github.com/Reference-LAPACK/lapack",
              "6acc99d5f39130be7cec00fb835606042101a970"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lapack*

mkdir build && cd build
cmake .. \
   -DCMAKE_INSTALL_PREFIX="$prefix" \
   -DCMAKE_FIND_ROOT_PATH="$prefix" \
   -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
   -DCMAKE_BUILD_TYPE=Release \
   -DBUILD_SHARED_LIBS=ON \
   -DBLAS_LIBRARIES="-L${libdir} -lblastrampoline"

VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("liblapack", :liblapack),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
        Dependency("CompilerSupportLibraries_jll")
        Dependency("libblastrampoline_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@info "Starting distributed tests..."

include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
disttestfile = joinpath(@__DIR__, "run_dtests.jl")
push!(test_exeflags.exec,"--color=yes")
cmd = `$test_exename $test_exeflags $disttestfile`

if !success(pipeline(cmd; stdout=stdout, stderr=stderr)) && ccall(:jl_running_on_valgrind, Cint, ()) == 0
    @error "Distributed test failed, cmd : $cmd"
end


text:
using Base.Test

import RequestsCache: Session, CachedSession
#import RequestsCache: create_query, execute
import RequestsCache: get

#session = Session()
#session = CachedSession()
session = CachedSession(cache_name="cache.jld", backend="jld", expire_after=Base.Dates.Day(1))
#println(session)

#prepared_query = create_query("GET", "http://httpbin.org/get", query = Dict("title" => "page1"))
#response = execute(prepared_query; session=session)

response = get(session, "http://httpbin.org/get"; query = Dict("title" => "page1"))

parsed_resp = Requests.json(response)
@test parsed_resp["args"]["title"] == "page1"
@test parsed_resp["headers"]["Host"] == "httpbin.org"

response = get(session, "http://httpbin.org/get"; query = Dict("title" => "page1"))

parsed_resp = Requests.json(response)
@test parsed_resp["args"]["title"] == "page1"
@test parsed_resp["headers"]["Host"] == "httpbin.org"


text:
#= src/cur.jl

References:

  S.A. Goreinov, E.E. Tyrtyshnikov, N.L. Zamarashkin. A theory of pseudoskeleton
    approximations. Linear Algebra Appl. 261: 1-21, 1997.

  S. Voronin, P.-G. Martinsson. Efficient algorithms for CUR and interpolative
    matrix decompositions. Preprint, arXiv:1412.8447 [math.NA].
=#

# CURPackedU

abstract type AbstractCURPackedU <: Factorization{Int} end

mutable struct CURPackedU <: AbstractCURPackedU
  rows::Vector{Int}
  cols::Vector{Int}
end

mutable struct HermitianCURPackedU <: AbstractCURPackedU
  cols::Vector{Int}
end
const HermCURPackedU = HermitianCURPackedU

mutable struct SymmetricCURPackedU <: AbstractCURPackedU
  cols::Vector{Int}
end
const SymCURPackedU = SymmetricCURPackedU

const HermOrSymCURPackedU = Union{HermCURPackedU, SymCURPackedU}

copy(A::CURPackedU) = CURPackedU(copy(rows), copy(cols))
copy(A::SymCURPackedU) = SymCURPackedU(copy(cols))
copy(A::HermCURPackedU) = HermCURPackedU(copy(cols))

function getindex(A::CURPackedU, d::Symbol)
  if     d == :cols  return A.cols
  elseif d == :k     return length(A.cols)
  elseif d == :rows  return A.rows
  else               throw(KeyError(d))
  end
end
function getindex(A::HermOrSymCURPackedU, d::Symbol)
  if     d in (:cols, :rows)  return A.cols
  elseif d == :k              return length(A.cols)
  else                        throw(KeyError(d))
  end
end

ndims(A::AbstractCURPackedU) = 2

size(A::AbstractCURPackedU) = (length(A[:cols]), length(A[:rows]))
size(A::AbstractCURPackedU, dim::Integer) =
  (dim == 1 ? length(A[:cols]) : (dim == 2 ? length(A[:rows]) : 1))

# CUR

abstract type AbstractCUR{T} <: Factorization{T} end

mutable struct CUR{T} <: AbstractCUR{T}
  rows::Vector{Int}
  cols::Vector{Int}
  C::Matrix{T}
  U::Factorization{T}
  R::Matrix{T}
end

mutable struct HermitianCUR{T} <: AbstractCUR{T}
  cols::Vector{Int}
  C::Matrix{T}
  U::Factorization{T}
end
const HermCUR = HermitianCUR

mutable struct SymmetricCUR{T} <: AbstractCUR{T}
  cols::Vector{Int}
  C::Matrix{T}
  U::Factorization{T}
end
const SymCUR = SymmetricCUR

const HermOrSymCUR{T} = Union{HermCUR{T}, SymCUR{T}}

function CUR(A::AbstractMatOrLinOp, U::CURPackedU)
  rows = U[:rows]
  cols = U[:cols]
  C = A[:,cols]
  R = A[rows,:]
  U, σ, V = svd!(C[rows,:])
  U2 = PartialSVD(Matrix(V), 1 ./ σ, Matrix(U'))
  CUR(rows, cols, C, U2, R)
end
CUR(A::AbstractMatOrLinOp, U::HermCURPackedU) = HermCUR(A, U)
CUR(A::AbstractMatOrLinOp, U::SymCURPackedU) = SymCUR(A, U)
CUR(A::AbstractMatOrLinOp, rows, cols) = CUR(A, CURPackedU(rows, cols))
CUR(A, args...) = CUR(LinOp(A), args...)

function HermCUR(A::AbstractMatOrLinOp, U::HermCURPackedU)
    cols = U[:cols]
    C = A[:,cols]
    F = eigen!(Hermitian(C[cols,:]))
    U = PartialHermEigen(1 ./ F.values, F.vectors)
    HermCUR(cols, C, U)
end

HermCUR(A::AbstractMatOrLinOp, cols) = HermCUR(A, HermCURPackedU(cols))
HermCUR(A, args...) = HermCUR(LinOp(A), args...)

function SymCUR(A::AbstractMatOrLinOp, U::SymCURPackedU)
  cols = U[:cols]
  C = A[:,cols]
  U, σ, V = svd!(C[cols,:])
  U2 = PartialSVD(V, 1 ./ σ, U')
  SymCUR(cols, C, U2)
end
SymCUR(A::AbstractMatOrLinOp, cols) = SymCUR(A, SymCURPackedU(cols))
SymCUR(A, args...) = SymCUR(LinOp(A), args...)

conj!(A::CUR) = CUR(A.rows, A.cols, conj!(A.C), conj!(A.U), conj!(A.R))
conj(A::CUR) = CUR(A.rows, A.cols, conj(A.C), conj(A.U), conj(A.R))
conj!(A::HermCUR) = HermCUR(A.cols, conj!(A.C), conj!(A.U))
conj(A::HermCUR) = HermCUR(A.cols, conj(A.C), conj(A.U))
conj!(A::SymCUR) = SymCUR(A.cols, conj!(A.C), conj!(A.U))
conj(A::SymCUR) = SymCUR(A.cols, conj(A.C), conj(A.U))

convert(::Type{CUR{T}}, A::CUR) where {T} =
  CUR(A.rows, A.cols, convert(Matrix{T}, A.C), convert(Factorization{T}, A.U),
      convert(Matrix{T}, A.R))
convert(::Type{HermCUR{T}}, A::HermCUR) where {T} =
  HermCUR(A.cols, convert(Matrix{T}, A.C), convert(Factorization{T}, A.U))
convert(::Type{SymCUR{T}}, A::SymCUR) where {T} =
  SymCUR(A.cols, convert(Matrix{T}, A.C), convert(Factorization{T}, A.U))
convert(::Type{AbstractCUR{T}}, A::CUR) where {T} = convert(CUR{T}, A)
convert(::Type{AbstractCUR{T}}, A::HermCUR) where {T} = convert(HermCUR{T}, A)
convert(::Type{AbstractCUR{T}}, A::SymCUR) where {T} = convert(SymCUR{T}, A)
convert(::Type{Factorization{T}}, A::CUR) where {T} = convert(CUR{T}, A)
convert(::Type{Factorization{T}}, A::HermCUR) where {T} = convert(HermCUR{T}, A)
convert(::Type{Factorization{T}}, A::SymCUR) where {T} = convert(SymCUR{T}, A)
convert(::Type{Array}, A::AbstractCUR) = Matrix(A)
convert(::Type{Array{T}}, A::AbstractCUR) where {T} = convert(Array{T}, Matrix(A))
convert(::Type{Matrix{T}}, A::AbstractCUR) where {T} = convert(Array{T}, Matrix(A))
Array(A::AbstractCUR) = convert(Array, A)
Matrix(A::AbstractCUR) = convert(Array, A)
Array{T}(A::AbstractCUR) where T = convert(Array{T}, A)
Matrix{T}(A::AbstractCUR) where T = convert(Array{T}, A)

copy(A::CUR) = CUR(copy(A.rows), copy(A.cols), copy(A.C), copy(A.U), copy(A.R))
copy(A::HermCUR) = HermCUR(copy(A.cols), copy(A.C), copy(A.U))
copy(A::SymCUR) = SymCUR(copy(A.cols), copy(A.C), copy(A.U))

Matrix(A::CUR{T}) where {T} = A[:C]*(A[:U]*A[:R])
Matrix(A::HermCUR{T}) where {T} = A[:C]*(A[:U]*A[:C]')
Matrix(A::SymCUR{T}) where {T} = A[:C]*(A[:U]*transpose(A[:C]))

  adjoint(A::AbstractCUR) = Adjoint(A)
  transpose(A::AbstractCUR) = Transpose(A)

function getindex(A::CUR{T}, d::Symbol) where T
  if     d == :C     return A.C
  elseif d == :R     return A.R
  elseif d == :U     return A.U
  elseif d == :cols  return A.cols
  elseif d == :k     return length(A.cols)
  elseif d == :rows  return A.rows
  else               throw(KeyError(d))
  end
end
function getindex(A::HermCUR{T}, d::Symbol) where T
  if     d == :C              return A.C
  elseif d == :R              return A.C'
  elseif d == :U              return A.U
  elseif d in (:cols, :rows)  return A.cols
  elseif d == :k              return length(A.cols)
  else                        throw(KeyError(d))
  end
end
function getindex(A::SymCUR{T}, d::Symbol) where T
  if     d == :C              return A.C
  elseif d == :R              return transpose(A.C)
  elseif d == :U              return A.U
  elseif d in (:cols, :rows)  return A.cols
  elseif d == :k              return length(A.cols)
  else                        throw(KeyError(d))
  end
end

ishermitian(::CUR) = false
issymmetric(::CUR) = false
ishermitian(::HermCUR) = true
issymnetric(A::HermCUR{T}) where {T} = isreal(A)
ishermitian(A::SymCUR{T}) where {T} = isreal(A)
issymmetric(::SymCUR) = true

isreal(::AbstractCUR{T}) where {T} = T <: Real

ndims(A::AbstractCUR) = 2

size(A::CUR) = (size(A.C,1), size(A.R,2))
size(A::CUR, dim::Integer) =
  (dim == 1 ? size(A.C,1) : (dim == 2 ? size(A.R,2) : 1))
size(A::HermOrSymCUR) = (size(A.C,1), size(A.C,1))
size(A::HermOrSymCUR, dim::Integer) = dim == 1 || dim == 2 ? size(A.C,1) : 1

  # BLAS/LAPACK multiplication routines

  ## CUR left-multiplication

  mul!(C::AbstractVector{T}, A::CUR{T}, B::AbstractVector{T}) where {T} =
    mul!(C, A[:C], A[:U]*(A[:R]*B))
  mul!(C::AbstractMatrix{T}, A::CUR{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:C], A[:U]*(A[:R]*B))

  ## CUR right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::CUR{T}) where {T} =
    mul!(C, (A*B[:C])*B[:U], B[:R])

  ## HermCUR left-multiplication

  mul!(C::AbstractVector{T}, A::HermCUR{T}, B::AbstractVector{T}) where {T} =
    mul!(C, A[:C], A[:U]*(A[:C]'*B))
  mul!(C::AbstractMatrix{T}, A::HermCUR{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:C], A[:U]*(A[:C]'*B))

  for Adj in (:Adjoint, :Transpose)
    @eval begin
      ## CUR left-multiplication
      function mul!(C::AbstractMatrix{T}, A::CUR{T}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        tmp = A[:R]*Bc
        mul!(C, A[:C], A[:U]*tmp)
      end

      function mul!(C::AbstractVector{T}, Ac::$Adj{T,CUR{T}}, B::AbstractVector{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * B
        tmp = $Adj(A[:U]) * tmp
        mul!(C, $Adj(A[:R]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,CUR{T}}, B::AbstractMatrix{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * B
        tmp = $Adj(A[:U]) * tmp
        mul!(C, $Adj(A[:R]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,CUR{T}}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * Bc
        tmp = $Adj(A[:U]) * tmp
        mul!(C, $Adj(A[:R]), tmp)
      end
        ## CUR right-multiplication
      function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,CUR{T}}) where T
        B = parent(Bc)
        tmp = A * $Adj(B[:R])
        tmp = tmp * $Adj(B[:U])
        mul!(C, tmp, $Adj(B[:C]))
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::CUR{T}) where T
        tmp = Ac * B[:C]
        mul!(C, tmp*B[:U], B[:R])
      end

      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, Bc::$Adj{T,CUR{T}}) where T
        B = parent(Bc)
        tmp = Ac * $Adj(B[:R])
        tmp = tmp * $Adj(B[:U])
        mul!(C, tmp, $Adj(B[:C]))
      end
    end
  end

  ## HermCUR left-multiplication

  mul!(C::AbstractMatrix{T}, A::HermCUR{T}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, A[:C], A[:U]*(A[:C]'*Bc))
  mul!(C::AbstractMatrix{T}, A::HermCUR{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T<:Real} =
    mul!(C, A[:C], A[:U]*(transpose(A[:C])*Bt))
  mul!!(C::AbstractMatrix{T}, A::HermCUR{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T<:Complex} =
    mul!(C, A, conj!(parent(Bt))')  # overwrites B
  function mul!(C::AbstractMatrix{T}, A::HermCUR{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where T<:Complex
    B = parent(Bt)
    size(B, 1) <= A[:k] && return mul!!(C, A, transpose(copy(B)))
    mul!(C, A[:C], A[:U]*((A[:C]')*Bt))
  end

  mul!(C::AbstractVector{T}, Ac::Adjoint{T,HermCUR{T}}, B::AbstractVector{T}) where {T} =
    mul!(C, parent(Ac), B)
  mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,HermCUR{T}}, B::AbstractMatrix{T}) where {T} =
    mul!(C, parent(Ac), B)
  function mul!(C::AbstractVector{T}, At::Transpose{T,HermCUR{T}}, B::AbstractVector{T}) where T
    A = parent(At)
    tmp = transpose(A[:U])*(transpose(A[:C])*B)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,HermCUR{T}}, B::AbstractMatrix{T}) where T
    A = parent(At)
    tmp = transpose(A[:U])*(transpose(A[:C])*B)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end

  mul!(C::AbstractMatrix{T}, At::Adjoint{T,HermCUR{T}}, Bt::Adjoint{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, parent(At), Bt)
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,HermCUR{T}}, Bt::Transpose{T,<:AbstractMatrix{T}}) where T
    A = parent(At)
    B = parent(Bt)
    tmp = transpose(A[:U])*(transpose(A[:C])*Bt)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end

  ## HermCUR right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::HermCUR{T}) where {T} =
    mul!(C, (A*B[:C])*B[:U], B[:C]')

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Adjoint{T,HermCUR{T}}) where {T} =
    mul!(C, A, parent(Bt))
  function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Transpose{T,HermCUR{T}}) where T
    B = parent(Bt)
    tmp = conj!(A)*B[:C]
    mul!(C, conj!(tmp)*transpose(B[:U]), transpose(B[:C]))
  end  # overwrites A
  function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Transpose{T,HermCUR{T}}) where T
    B = parent(Bt)
    size(A, 1) <= B[:k] && return mul!!(C, copy(A), Bt)
    mul!(C, (A*conj(B[:C]))*transpose(B[:U]), transpose(B[:C]))
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::HermCUR{T}) where T
        tmp = Ac * B[:C]
        mul!(C, tmp*B[:U], B[:C]')
      end
    end
  end

  mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,<:AbstractMatrix{T}}, Bc::Adjoint{T,HermCUR{T}}) where {T} =
    mul!(C, Ac, parent(Bc))
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,<:AbstractMatrix{T}}, Bt::Transpose{T,HermCUR{T}}) where {T}
    B = parent(Bt)
    A = parent(At)
    mul!(C, conj!(A'*B[:C])*transpose(B[:U]), transpose(B[:C]))
  end

  ## SymCUR left-multiplication

  mul!(C::AbstractVector{T}, A::SymCUR{T}, B::AbstractVector{T}) where {T} =
    mul!(C, A[:C], A[:U]*(transpose(A[:C])*B))
  mul!(C::AbstractMatrix{T}, A::SymCUR{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:C], A[:U]*(transpose(A[:C])*B))
  mul!(C::AbstractMatrix{T}, A::SymCUR{T}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where {T<:Real} =
    mul!(C, A[:C], A[:U]*(A[:C]'*Bc))
  mul!!(C::AbstractMatrix{T}, A::SymCUR{T}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where {T<:Complex} =
    mul!(C, A, transpose(conj!(parent(Bc))))  # overwrites B
  function mul!(C::AbstractMatrix{T}, A::SymCUR{T}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where T<:Complex
    B = parent(Bc)
    size(B, 1) <= A[:k] && return mul!!(C, A, copy(B)')
    mul!(C, A[:C], A[:U]*(transpose(A[:C])*B'))
  end
  mul!(C::AbstractMatrix{T}, A::SymCUR{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, A[:C], A[:U]*(transpose(A[:C])*Bt))

  function mul!(C::AbstractVector{T}, Ac::Adjoint{T,SymCUR{T}}, B::AbstractVector{T}) where T
    A = parent(Ac)
    tmp = A[:U]'*(A[:C]'*B)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end
  function mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,SymCUR{T}}, B::AbstractMatrix{T}) where T
    A = parent(Ac)
    tmp = A[:U]'*(A[:C]'*B)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end
  mul!(C::AbstractVector{T}, At::Transpose{T,SymCUR{T}}, B::AbstractVector{T}) where {T} =
    mul!(C, parent(At), B)
  mul!(C::AbstractMatrix{T}, At::Transpose{T,SymCUR{T}}, B::AbstractMatrix{T}) where {T} =
    mul!(C, parent(At), B)

  function mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,SymCUR{T}}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where T
    A = parent(Ac)
    tmp = A[:U]'*(A[:C]'*Bc)
    mul!(C, A[:C], conj!(tmp))
    conj!(C)
  end
  mul!(C::AbstractMatrix{T}, At::Transpose{T,SymCUR{T}}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, parent(At), Bt)

  ## SymCUR right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::SymCUR{T}) where {T} =
    mul!(C, (A*B[:C])*B[:U], transpose(B[:C]))

  function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::Adjoint{T,SymCUR{T}}) where T
    B = parent(Bc)
    tmp = conj!(A)*B[:C]
    mul!(C, conj!(tmp)*transpose(B[:U]), transpose(B[:C]))
  end  # overwrites A
  function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::Adjoint{T,SymCUR{T}}) where T
    B = parent(Bc)
    size(A, 1) <= B[:k] && return mul!!(C, copy(A), Bc)
    mul!(C, (A*conj(B[:C]))*B[:U]', B[:C]')
  end
  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Transpose{T,SymCUR{T}}) where {T} =
    mul!(C, A, parent(Bt))

  for Adj in (:Adjoint, :Transpose)
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::SymCUR{T}) where T
        tmp = Ac * B[:C]
        mul!(C, tmp*B[:U], transpose(B[:C]))
      end
    end
  end

  function mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,<:AbstractMatrix{T}}, Bc::Adjoint{T,SymCUR{T}}) where {T}
    B = parent(Bc)
    A = parent(Ac)
    mul!(C, conj!(transpose(A)*B[:C])*B[:U]', B[:C]')
  end
  mul!(C::AbstractMatrix{T}, At::Transpose{T,<:AbstractMatrix{T}}, Bt::Transpose{T,SymCUR{T}}) where {T} =
    mul!(C, At, parent(Bt))

  # standard operations

  ## left-multiplication
  function *(A::AbstractCUR{TA}, B::AbstractVector{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(AbstractCUR{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1))
    mul!(CT, AT, BT)
  end

  for Adj in (:Adjoint, :Transpose)
    @eval function *(Ac::$Adj{TA,AbstractCUR{TA}}, B::AbstractVector{TB}) where {TA,TB}
      A = parent(Ac)
      T = promote_type(TA, TB)
      AT = convert(AbstractCUR{T}, A)
      BT = (T == TB ? B : convert(Array{T}, B))
      CT = Array{T}(undef, size(A,2))
      mul!(CT, $Adj(AT), BT)
    end
  end

  function *(A::AbstractCUR{TA}, B::AbstractMatrix{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(AbstractCUR{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end
  function *(A::AbstractMatrix{TA}, B::AbstractCUR{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = (T == TA ? A : convert(Array{T}, A))
    BT = convert(AbstractCUR{T}, B)
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end

  for Adj in (:Adjoint, :Transpose)
    @eval begin
      function *(A::AbstractCUR{TA}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(AbstractCUR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:AbstractCUR{TA}}, B::AbstractMatrix{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(AbstractCUR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:AbstractCUR{TA}}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(AbstractCUR{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
      function *(A::AbstractMatrix{TA}, Bc::$Adj{TB,<:AbstractCUR{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(AbstractCUR{T}, B)
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, B::AbstractCUR{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(AbstractCUR{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, Bc::$Adj{TB,<:AbstractCUR{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(AbstractCUR{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
    end
  end

# factorization routines

for sfx in ("", "!")
  f = Symbol("curfact", sfx)
  g = Symbol("curfact_none", sfx)
  h = Symbol("cur", sfx)
  @eval begin
    function $f(A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
      opts = copy(opts; args...)
      opts.pqrfact_retval = ""
      chkopts!(opts, A)
      opts.sketch == :none && return $g(A, opts)
      if ishermitian(A)
        F = sketchfact(:left, :n, A, opts)
        cols = F[:p][1:F[:k]]
        return HermCURPackedU(cols)
      elseif issymmetric(A)
        F = sketchfact(:left, :n, A, opts)
        cols = F[:p][1:F[:k]]
        return SymCURPackedU(cols)
      end
      m, n = size(A)
      if m >= n
        F = sketchfact(:left, :c, A, opts)
        rows = F[:p][1:F[:k]]
        F = sketchfact(:left, :n, A[rows,:], opts)
        cols = F[:p][1:F[:k]]
      else
        F = sketchfact(:left, :n, A, opts)
        cols = F[:p][1:F[:k]]
        F = sketchfact(:left, :c, A[:,cols], opts)
        rows = F[:p][1:F[:k]]
      end
      kr = length(rows)
      kc = length(cols)
      k = min(kr, kc)
      rows = k < kr ? rows[1:k] : rows
      cols = k < kc ? cols[1:k] : cols
      CURPackedU(rows, cols)
    end
    $f(A, args...; kwargs...) = $f(LinOp(A), args...; kwargs...)

    function $h(A, args...; kwargs...)
      U = $f(A, args...; kwargs...)
      U[:rows], U[:cols]
    end
  end
end

function curfact_none!(A::AbstractMatrix, opts::LRAOptions)
  if ishermitian(A)
    F = pqrfact_backend!(A, opts)
    cols = F[:p][1:F[:k]]
    return HermCURPackedU(cols)
  elseif issymmetric(A)
    F = pqrfact_backend!(A, opts)
    cols = F[:p][1:F[:k]]
    return SymCURPackedU(cols)
  end
  curfact_none(A, opts)
end
function curfact_none(A::AbstractMatrix, opts::LRAOptions)
  (ishermitian(A) || issymmetric(A)) && return curfact_none!(copy(A), opts)
  m, n = size(A)
  if m >= n
    F = pqrfact_backend!(Matrix(A'), opts)
    rows = F[:p][1:F[:k]]
    F = pqrfact_backend!(A[rows,:], opts)
    cols = F[:p][1:F[:k]]
  else
    F = pqrfact_backend!(copy(A), opts)
    cols = F[:p][1:F[:k]]
    F = pqrfact_backend!(Matrix(A[:,cols]'), opts)
    rows = F[:p][1:F[:k]]
  end
  kr = length(rows)
  kc = length(cols)
  k = min(kr, kc)
  rows = k < kr ? rows[1:k] : rows
  cols = k < kc ? cols[1:k] : cols
  return CURPackedU(rows, cols)
end


text:
using Distributions

"Boostrapping for degnerate U-statistic."
bootstrap_degenerate(samples::AbstractArray, f) = begin
    # bootstrap sample
    res = 0.0

    n = size(samples, 2)

    # sample bootstrap weights
    sample_p = Multinomial(n, n)
    w = (rand(sample_p) / n) .- (1.0 / n)

    for i = 1:n
        for j = 1:n
            if i == j
                continue
            end
            r = f(samples[:, i], samples[:, j])
            res += ((w[i] * w[j]) * r)
        end
    end

    return res
end


text:
mutable struct Accumulator
    value::Float64
end

# Functor - callable type
function (acc::Accumulator)(y::Real)
    acc.value += y
    return acc.value
end

@show a = Accumulator(1.0);
@show a(0.2);
@show a(0.3);

# using a regular function
function add!(acc::Accumulator, y::Real)
    acc.value += y
    return acc.value
end

@show add!(a, 0.4);

# Additional functors - special for Int arguments
function (acc::Accumulator)(y::Int)
    acc.value += 2y
    return acc.value
end

@show a = Accumulator(0);
@show a(0.2);
@show a(0.3);
@show a(3);


text:
using Accounts
using Test

const FILE_ACCOUNTS = "./test_accounts.txt"

@testset "Accounts.jl" begin
    donald_email = create(EMAIL, "donald@duckcity.com")
    donald = create("Donald Duck", [donald_email])
    add_to_file(FILE_ACCOUNTS, [donald])
    result = read_from_file(FILE_ACCOUNTS)
    first_person = result[1]
    @test first_person.addresses[1].address == "donald@duckcity.com"
    cmd = `rm $FILE_ACCOUNTS` # linux remove file statement
    run(cmd) # remove file
end


text:
export PiecewiseUniformEnergyDistribution

"Initial piecewise uniform energy distribution."
struct PiecewiseUniformEnergyDistribution <: InitialEnergyDistribution{1}
    interior_boundary_positions::Vector{Float}
    specific_energies::Vector{Float}

    function PiecewiseUniformEnergyDistribution(
        interior_boundary_positions::AbstractVector{Float},
        specific_energies::AbstractVector{Float},
    )
        @checkeq(
            length(specific_energies),
            length(interior_boundary_positions) + 1
        )
        @check all(diff(interior_boundary_positions) .> 0)
        @check all(specific_energies .> 0)
        new(interior_boundary_positions, specific_energies)
    end
end

function (distribution::PiecewiseUniformEnergyDistribution)(
    positions::AbstractVector{SVectorF{1}},
    mass_densities::AbstractVector{Float},
    eos::EquationOfState,
)
    @checkeq(length(positions), length(mass_densities))
    [specific_energy_at(distribution, r[]) for r in positions]
end

"""
    specific_energy_at(distr::PiecewiseUniformEnergyDistribution, r) -> Float

Return the energy density of the piecewise uniform energy distribution
at the given position `r`.
"""
specific_energy_at(distribution::PiecewiseUniformEnergyDistribution, r) =
    distribution.specific_energies[1+searchsortedlast(
        distribution.interior_boundary_positions,
        r,
    )]

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
@Base.deprecate parseWKT(geom::String) readgeom(geom)
@Base.deprecate geomFromWKT(geom::String) _readgeom(geom)
@Base.deprecate geomToWKT(geom::Ptr{GEOSGeometry}) writegeom(geom)

@Base.deprecate_binding CAP_ROUND GEOSBUF_CAP_ROUND
@Base.deprecate_binding CAP_FLAT GEOSBUF_CAP_FLAT
@Base.deprecate_binding CAP_SQUARE GEOSBUF_CAP_SQUARE
@Base.deprecate_binding JOIN_ROUND GEOSBUF_JOIN_ROUND
@Base.deprecate_binding JOIN_MITRE GEOSBUF_JOIN_MITRE
@Base.deprecate_binding JOIN_BEVEL GEOSBUF_JOIN_BEVEL


text:
using Documenter, Python, Markdown

struct CustomCat{cat} end
Documenter.Utilities.doccat(::Base.Docs.Binding, ::Type{CustomCat{cat}}) where {cat} = string(cat)
struct CustomDocBlocks <: Documenter.Expanders.ExpanderPipeline end
Documenter.Expanders.Selectors.order(::Type{CustomDocBlocks}) = 20.0
Documenter.Expanders.Selectors.matcher(::Type{CustomDocBlocks}, node, page, doc) = Documenter.Expanders.iscode(node, "@customdoc")
Documenter.Expanders.Selectors.runner(::Type{CustomDocBlocks}, x, page, doc) = begin
    header, rest = split(x.code, "\n", limit=2)
    docstr = Markdown.parse(rest)
    name, cat = split(header, "-", limit=2)
    binding = Docs.Binding(Main, Symbol(strip(name)))
    object = Documenter.Utilities.Object(binding, CustomCat{Symbol(strip(cat))})
    slug = Documenter.Utilities.slugify(strip(name))
    anchor = Documenter.Anchors.add!(doc.internal.docs, object, slug, page.build)
    node = Documenter.Documents.DocsNode(docstr, anchor, object, page)
    page.mapping[x] = node
end

makedocs(
    sitename = "Python.jl",
    modules = [Python],
    pages = [
        "Home" => "index.md",
        "getting-started.md",
        "pythonjl.md",
        "juliapy.md",
        "conversion.md",
        "compat.md",
        "pycall.md",
    ]
)

deploydocs(
    repo = "github.com/cjdoris/Python.jl.git",
)


text:
# module for representing puddle

using Plots
pyplot()

mutable struct Puddle
  lower_left
  upper_right
  depth

  # init
  function Puddle(lower_left, upper_right, depth)
    self = new()
    self.lower_left = lower_left
    self.upper_right = upper_right
    self.depth = depth
    return self
  end
end

function draw!(self::Puddle)
  w = self.upper_right[1] - self.lower_left[1]
  h = self.upper_right[2] - self.lower_left[2]
  rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
  plot!(rectangle(w, h, self.lower_left[1], self.lower_left[2]),
        color=:blue, alpha=self.depth, aspect_ratio=true,
        legend=false)
end

function inside(self::Puddle, pose)
  result = true
  for i in [1, 2]
    if (self.lower_left[i] > pose[i]) ||
       (self.upper_right[i] < pose[i])
      result = false
    end
  end
  return result
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
struct Lwplsr
    X::Array{Float64}
    Y::Array{Float64}
    fm
    metric::String
    h::Real
    k::Int
    nlv::Int
    tol::Real
    verbose::Bool
end

"""
    lwplsr(X, Y; nlvdis, metric, h, k, nlv, tol = 1e-4, verbose = false)
k-Nearest-Neighbours locally weighted partial least squares regression (KNN-LWPLSR).
* `X` : X-data.
* `Y` : Y-data.
* `nlvdis` : Number of latent variables (LVs) to consider in the global PLS
    used for the dimension reduction before computing the dissimilarities.
    If `nlvdis = 0`, there is no dimension reduction.
* `metric` : Type of dissimilarity used to select the neighbors.
    Possible values are "eucl" (default; Euclidean distance)
    and "mahal" (Mahalanobis distance).
* `h` : A scalar defining the shape of the weight function. Lower is h,
    sharper is the function. See function `wdist`.
* `k` : The number of nearest neighbors to select for each observation to predict.
* `nlv` : Nb. latent variables (LVs).
* `tol` : For stabilization when very close neighbors.
* `verbose` : If true, fitting information are printed.

Function `lwplsr` fits kNN-LWPLSR models (Lesnoff et al., 2020).
The function uses functions `getknn`, `locw` and a PLSR function;
see the code for details. Many other variants of kNN-LWPLSR pipelines can be built.

The general principles of the method are as follows.

LWPLSR is a particular case of weighted PLSR (WPLSR) (e.g. Schaal et al. 2002).
In WPLSR, a priori weights, different from the usual 1/n (standard PLSR),
are given to the n training observations. These weights are used for calculating
(i) the scores and loadings of the WPLS and (ii) the regression model that fits
(by weighted least squares) the Y-response(s) to the WPLS scores.
The specificity of LWPLSR (compared to WPLSR) is that the weights are computed
from dissimilarities (e.g. distances) between the new observation to predict
and the training observations ("L" in LWPLSR comes from "localized").
Note that in LWPLSR the weights and therefore the fitted WPLSR model
change for each new observation to predict.

In the original LWPLSR, all the n training observations are used for each
observation to predict (e.g. Sicard & Sabatier 2006, Kim et al 2011).
This can be very time consuming, in particular for large n.
A faster and often more efficient strategy is to preliminary select,
in the training set, a number of `k` nearest neighbors to the observation to predict
(= "weighting 1") and then to apply LWPLSR only to this pre-selected
neighborhood (= "weighting 2"). This strategy corresponds to a kNN-LWPLSR
and is the one implemented in function `lwplsr`.

In `lwplsr`, the dissimilarities used for weightings 1 and 2 are
computed from the raw X-data or after a dimension reduction,
depending on argument `nlvdis`. In the last case, global PLS2 scores (LVs) are
computed from {`X`, `Y`} and the dissimilarities are computed over these scores.

In general, for high dimensional X-data, using the Mahalanobis distance requires
preliminary dimensionality reduction of the data.

## References

Kim, S., Kano, M., Nakagawa, H., Hasebe, S., 2011. Estimation of active
pharmaceutical ingredients content using locally weighted partial least squares
and statistical wavelength selection. Int. J. Pharm., 421, 269-274.

Lesnoff, M., Metz, M., Roger, J.-M., 2020. Comparison of locally weighted PLS
strategies for regression and discrimination on agronomic NIR data.
Journal of Chemometrics, e3209. https://doi.org/10.1002/cem.3209

Schaal, S., Atkeson, C., Vijayamakumar, S. 2002. Scalable techniques from nonparametric
statistics for the real time robot learning. Applied Intell., 17, 49-60.

Sicard, E. Sabatier, R., 2006. Theoretical framework for local PLS1 regression
and application to a rainfall data set. Comput. Stat. Data Anal., 51, 1393-1410.
"""
function lwplsr(X, Y; nlvdis, metric, h, k, nlv, tol = 1e-4, verbose = false)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    if nlvdis == 0
        fm = nothing
    else
        fm = plskern(X, Y; nlv = nlvdis)
    end
    return Lwplsr(X, Y, fm, metric, h, k, nlv, tol, verbose)
end

"""
    predict(object::Lwplsr, X)
Compute the Y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::Lwplsr, X; nlv = nothing)
    X = ensure_mat(X)
    m = size(X, 1)
    a = object.nlv
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    # Getknn
    if isnothing(object.fm)
        res = getknn(object.X, X; k = object.k, metric = object.metric)
    else
        res = getknn(object.fm.T, transform(object.fm, X); k = object.k, metric = object.metric)
    end
    listw = copy(res.d)
    for i = 1:m
        w = wdist(res.d[i]; h = object.h)
        w[w .< object.tol] .= object.tol
        listw[i] = w
    end
    # End
    pred = locwlv(object.X, object.Y, X;
        listnn = res.ind, listw = listw, fun = plskern, nlv = nlv,
        verbose = object.verbose).pred
    (pred = pred, listnn = res.ind, listd = res.d, listw = listw)
end


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2

D = 16

Q = Constant(randn(D,D))
R = Constant(randn(D,D))
Ψ = InfiniteCMPS(Q, R)

alg = ConjugateGradient(; verbosity = 2, maxiter = 10^6, gradtol = 1e-4);

#H = ∫(∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2] - 10 * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + 100 * (ψ[1]*ψ[2] + ψ[2]*ψ[1]), (-Inf,+Inf))

#H = ∫(∂ψ'*∂ψ - 1 * ψ'*ψ + 1 * (ψ')^2*ψ^2, (-Inf,+Inf))

# function finalize!(x, f, g, numiter)
#     if mod(numiter, 1000) == 0
#         @save "tempoptresults2.jld" x
#     end
#     return x, f, g
# end
# groundstate(H, Ψ; optalg = alg, linalg = GMRES(; tol = 1e-4), finalize! = finalize!)
groundstate(H, Ψ; optalg = alg, linalg = GMRES(; tol = 1e-6))


text:
function plot(
    X::Vector{Float64}, Y::Vector{Float64},
    COLOR::String, LINEWIDTH::Float64;
    dashed=false,
    dotted=false,
    dasharray=nothing
    )
    @assert length(X)==length(Y)
    svg_lines(  [([X[i],Y[i]],[X[i+1],Y[i+1]]) for i=1:length(X)-1],
                COLOR,
                LINEWIDTH;
                dashed=dashed,
                dotted=dotted,
                dasharray=dasharray )
end

max_stars_repo_name: algorithmx/PlotSVG

max_issues_repo_name: algorithmx/PlotSVG

max_forks_repo_name: algorithmx/PlotSVG


text:
using AutoHOOT, ChainRulesCore

broadcast_notangent(a) = broadcast(_ -> NoTangent(), a)

function ChainRulesCore.rrule(::typeof(ITensors.data), P::PEPS)
  return P.data, d_data -> (NoTangent(), PEPS(d_data))
end

function ChainRulesCore.rrule(::typeof(PEPS), data::Matrix{ITensor})
  return PEPS(data), dpeps -> (NoTangent(), dpeps.data)
end

function ChainRulesCore.rrule(::typeof(ITensors.prime), P::PEPS, n::Integer=1)
  return prime(P, n), dprime -> (NoTangent(), prime(dprime, -n), NoTangent())
end

function ChainRulesCore.rrule(
  ::typeof(ITensors.addtags), ::typeof(linkinds), P::PEPS, args...
)
  function pullback(dtag_peps)
    dP = ITensors.removetags(linkinds, dtag_peps, args...)
    return (NoTangent(), NoTangent(), dP, broadcast_notangent(args)...)
  end
  return ITensors.addtags(linkinds, P, args...), pullback
end

function ChainRulesCore.rrule(
  ::typeof(ITensors.removetags), ::typeof(linkinds), P::PEPS, args...
)
  function pullback(dtag_peps)
    dP = ITensors.addtags(linkinds, dtag_peps, args...)
    return (NoTangent(), NoTangent(), dP, broadcast_notangent(args)...)
  end
  return ITensors.removetags(linkinds, P, args...), pullback
end

function ChainRulesCore.rrule(
  ::typeof(ITensors.prime), ::typeof(linkinds), P::PEPS, n::Integer=1
)
  return prime(linkinds, P, n),
  dprime -> (NoTangent(), NoTangent(), prime(linkinds, dprime, -n), NoTangent())
end

function ChainRulesCore.rrule(
  ::typeof(ITensors.prime), indices::Array{<:Index,1}, P::PEPS, n::Integer=1
)
  primeinds = [prime(ind, n) for ind in indices]
  return prime(indices, P, n),
  dprime -> (NoTangent(), NoTangent(), prime(primeinds, dprime, -n), NoTangent())
end

function ChainRulesCore.rrule(::typeof(flatten), v::Array{<:PEPS})
  size_list = [size(peps.data) for peps in v]
  function adjoint_pullback(dt)
    dt = [t for t in dt]
    index = 0
    dv = []
    for (dimy, dimx) in size_list
      size = dimy * dimx
      d_peps = PEPS(reshape(dt[(index + 1):(index + size)], dimy, dimx))
      index += size
      push!(dv, d_peps)
    end
    return (NoTangent(), dv)
  end
  return flatten(v), adjoint_pullback
end

@non_differentiable inner_network(peps::PEPS, peps_prime::PEPS)

@non_differentiable inner_network(peps::PEPS, peps_prime::PEPS, projectors::Vector{ITensor})

@non_differentiable inner_network(
  peps::PEPS, peps_prime::PEPS, projectors::Vector{<:ITensor}, ::typeof(tree_w_projectors)
)

@non_differentiable inner_network(peps::PEPS, peps_prime::PEPS, ::typeof(tree))

# gradient of this function returns nothing.
@non_differentiable inner_networks(
  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hs::Array
)

@non_differentiable inner_networks(
  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, projectors::Vector{ITensor}, Hs::Array
)

@non_differentiable inner_networks(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  projectors::Vector{Vector{ITensor}},
  Hs::Vector{Tuple},
)

@non_differentiable inner_networks(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  projectors::Vector{Vector{ITensor}},
  Hs::Vector{Models.LineMPO},
  ::typeof(tree_w_projectors),
)

@non_differentiable inner_networks(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  Hs::Vector{Models.LineMPO},
  ::typeof(tree),
)

@non_differentiable ITensors.commoninds(p1::PEPS, p2::PEPS)

@non_differentiable insert_projectors(peps::PEPS, cutoff, maxdim)

@non_differentiable insert_projectors(peps::PEPS, center::Tuple, cutoff, maxdim)


text:
using DataFrames

df = DataFrame(A=[1,2], B=["hello","world"])

df.c = df.A .^ 2

max_stars_repo_name: jondea/juliab

max_issues_repo_name: jondea/juliab

max_forks_repo_name: jondea/juliab


text:
function next(x,s)
		return (2*x - floor(x/π)*2π + s/4*sin(4π*x)) % (2π)
end


text:
using Test, Random

@testset "Testing solution to Exercise 1" begin

@testset "Running ex1.jl" begin
   include("../ex1.jl")
end;

@testset "Testing that variables exist" begin
   @test @isdefined(multiply_matrix_vector_default)
   @test @isdefined(multiply_matrix_vector_rows_inner)
   @test @isdefined(test_mul_mat_vec)
   @test @isdefined(multiply_matrix_vector_cols_inner)
   @test @isdefined(response_1c)
   @test @isdefined(response_1d)
   @test @isdefined(response_1e)
   @test @isdefined(response_1f)
   @test @isdefined(response_1g)
   @test @isdefined(response_1h)
   @test @isdefined(response_1i)
end;

@testset "Testing that variables are not missing" begin
   @test !ismissing(response_1c)
   @test !ismissing(response_1d)
   @test !ismissing(response_1e)
   @test !ismissing(response_1f)
   @test !ismissing(response_1g)
   @test !ismissing(response_1h)
   @test !ismissing(response_1i)
end;

function make_A_b_for_test(nrows::Integer, ncols::Integer; seed::Integer = 42)
   @assert(1<=nrows<=8192)
   @assert(1<=ncols<=8192)
   Random.seed!(seed)
   A = rand(nrows,ncols)
   b = rand(ncols)
   return (A,b)
end

@testset "Testing multiply_matrix_vector_default" begin
   (A,b) = make_A_b_for_test(3,4)
   @test !ismissing(multiply_matrix_vector_default(A,b))
   @test !ismissing(test_mul_mat_vec(multiply_matrix_vector_default))
   @test all(multiply_matrix_vector_default(A,b) .≈ A*b )
   (A,b) = make_A_b_for_test(17,23)
   @test all(multiply_matrix_vector_default(A,b) .≈ A*b )
end;

@testset "Testing multiply_matrix_vector_cols_inner" begin
   (A,b) = make_A_b_for_test(3,4)
   @test !ismissing(multiply_matrix_vector_cols_inner(A,b))
   @test all(multiply_matrix_vector_cols_inner(A,b) .≈ A*b )
   (A,b) = make_A_b_for_test(17,23)
   @test all(multiply_matrix_vector_cols_inner(A,b) .≈ A*b )
end;

@testset "Testing multiply_matrix_vector_rows_inner" begin
   (A,b) = make_A_b_for_test(3,4)
   @test !ismissing(multiply_matrix_vector_rows_inner(A,b))
   @test all(multiply_matrix_vector_rows_inner(A,b) .≈ A*b )
   (A,b) = make_A_b_for_test(17,23)
   @test all(multiply_matrix_vector_rows_inner(A,b) .≈ A*b )
end;

end; # Exercise 1