text:
using NoteSequences.MelodyRepr
using NoteSequences.MelodyRepr: MELODY_NOTE_OFF, MELODY_NO_EVENT

@testset "MelodyOneHotEncoding" begin
    @testset "Encoding/Decoding" begin
        encoder = MelodyOneHotEncoding(43, 81)

        # (Melody event, one-hot index) pairs
        pairs = [
            (MELODY_NO_EVENT, 1),
            (MELODY_NOTE_OFF, 2),
            (43, 3),
            (80, 40),
            (79, 39),
            (50, 10),
        ]

        for (event, index) in pairs
            @test index == encode_event(event, encoder)
            @test event == decode_event(index, encoder)
        end
    end

    @testset "Out of range initialization" begin
        @test_throws ErrorException MelodyOneHotEncoding(0, 129)
        @test_throws ErrorException MelodyOneHotEncoding(-1, 100)
        @test_throws ErrorException MelodyOneHotEncoding(50, 25)
        @test_throws ErrorException MelodyOneHotEncoding(10, 10)
    end

    @testset "Test labels" begin
        @test MelodyOneHotEncoding(0, 128).labels == 1:130
        @test MelodyOneHotEncoding(43, 81).labels == 1:40
        @test MelodyOneHotEncoding(0, 1).labels == 1:3
    end
end


text:
using Revise
	using BifurcationKit, LinearAlgebra, Plots, SparseArrays, Setfield, Parameters
	const BK = BifurcationKit

f1(u, v) = u * u * v
norminf = x -> norm(x, Inf)

function plotsol(x; kwargs...)
	N = div(length(x), 2)
	plot!(x[1:N], label="u"; kwargs...)
	plot!(x[N+1:2N], label="v"; kwargs...)
end

function Fbru!(f, x, p)
	@unpack α, β, D1, D2, l = p
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h
	c1 = D1 / l^2 / h2
	c2 = D2 / l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	# Dirichlet boundary conditions
	f[1]   = c1 * (α	  - 2u[1] + u[2] ) + α - (β + 1) * u[1] + f1(u[1], v[1])
	f[end] = c2 * (v[n-1] - 2v[n] + β / α)			 + β * u[n] - f1(u[n], v[n])

	f[n]   = c1 * (u[n-1] - 2u[n] +  α   ) + α - (β + 1) * u[n] + f1(u[n], v[n])
	f[n+1] = c2 * (β / α  - 2v[1] + v[2])			 + β * u[1] - f1(u[1], v[1])

	for i=2:n-1
		  f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + α - (β + 1) * u[i] + f1(u[i], v[i])
		f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + β * u[i] - f1(u[i], v[i])
	end
	return f
end

Fbru(x, p) = Fbru!(similar(x), x, p)

function Jbru_sp(x, p)
	@unpack α, β, D1, D2, l = p
	# compute the Jacobian using a sparse representation
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h

	c1 = D1 / p.l^2 / h2
	c2 = D2 / p.l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	diag   = zeros(eltype(x), 2n)
	diagp1 = zeros(eltype(x), 2n-1)
	diagm1 = zeros(eltype(x), 2n-1)

	diagpn = zeros(eltype(x), n)
	diagmn = zeros(eltype(x), n)

	@. diagmn = β - 2 * u * v
	@. diagm1[1:n-1] = c1
	@. diagm1[n+1:end] = c2

	@. diag[1:n]    = -2c1 - (β + 1) + 2 * u * v
	@. diag[n+1:2n] = -2c2 - u * u

	@. diagp1[1:n-1]   = c1
	@. diagp1[n+1:end] = c2

	@. diagpn = u * u
	J = spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
	return J
end

function finalise_solution(z, tau, step, contResult)
	n = div(length(z.u), 2)
	printstyled(color=:red, "--> Solution constant = ", norm(diff(z.u[1:n])), " - ", norm(diff(z.u[n+1:2n])), "\n")
	return true
end

n = 500
####################################################################################################
# test for the Jacobian expression
# using ForwardDiff
# sol0 = rand(2n)
# J0 = ForwardDiff.jacobian(x-> Fbru(x, par_bru), sol0) |> sparse
# J1 = Jbru_sp(sol0, par_bru)
# J0 - J1
####################################################################################################
# different parameters to define the Brusselator model and guess for the stationary solution
par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
	sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))

# # parameters for an isola of stationary solutions
# par_bru = (α = 2., β = 4.6, D1 = 0.0016, D2 = 0.008, l = 0.061)
# 	xspace = LinRange(0, par_bru.l, n)
# 	sol0 = vcat(		par_bru.α .+ 2 .* sin.(pi*xspace/par_bru.l),
# 			par_bru.β/par_bru.α .- 0.5 .* sin.(pi*xspace/par_bru.l))

# eigls = EigArpack(1.1, :LM)
# 	opt_newton = BK.NewtonPar(eigsolver = eigls)
# 	out, hist, flag = @time BK.newton(
# 		x ->    Fbru(x, par_bru),
# 		x -> Jbru_sp(x, par_bru),
# 		sol0, opt_newton, normN = norminf)
#
# 	plot();plotsol(out);plotsol(sol0, label = "sol0",line=:dash)
####################################################################################################
eigls = EigArpack(1.1, :LM)
opts_br_eq = ContinuationPar(dsmin = 0.03, dsmax = 0.05, ds = 0.03, pMax = 1.9, detectBifurcation = 3, nev = 21, plotEveryStep = 50, newtonOptions = NewtonPar(eigsolver = eigls, tol = 1e-9), maxSteps = 1060, nInversion = 6, tolBisectionEigenvalue = 1e-20, maxBisectionSteps = 30)

	br, = @time continuation(
		Fbru, Jbru_sp, sol0, par_bru, (@lens _.l),
		opts_br_eq, verbosity = 0,
		plot = true,
		plotSolution = (x, p; kwargs...) -> (plotsol(x; label="", kwargs... )),
		printSolution = (x, p) -> x[div(n,2)], normC = norminf)
####################################################################################################
jet  = BK.get3Jet(Fbru, Jbru_sp)
hopfpt = computeNormalForm(jet..., br, 1; verbose = true)
#################################################################################################### Continuation of the Hopf Point using Jacobian expression
ind_hopf = 1
	# hopfpt = BK.HopfPoint(br, ind_hopf)
	optnew = opts_br_eq.newtonOptions
	hopfpoint, _, flag = @time newton(
		Fbru, Jbru_sp,
		br, ind_hopf;
		d2F = d2Fbru,
		options = (@set optnew.verbose=true), normN = norminf)
	flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", hopfpoint.p[1], ", ω = ", hopfpoint.p[2], ", from l = ", br.specialpoint[ind_hopf].param, "\n")

if 1==0
	br_hopf, u1_hopf = @time continuation(
		Fbru, Jbru_sp,
		br, ind_hopf, (@lens _.β),
		ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, newtonOptions = optnew); verbosity = 2, normC = norminf)

	plot(br_hopf, label="")
end

hopfpoint, hist, flag = @time newton(
	Fbru, Jbru_sp,
	br, ind_hopf;
	options = (@set optnew.verbose = true),
	d2F = d2Fbru, normN = norminf)

if 1==1
	br_hopf, u1_hopf = @time continuation(
		Fbru, Jbru_sp,
		br, ind_hopf, (@lens _.β),
		ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 10.5, pMin = 5.1, detectBifurcation = 0, newtonOptions = optnew); plot = true,
		updateMinAugEveryStep = 1,
		d2F = d2Fbru,
		verbosity = 2, normC = norminf, bothside = true)
end
plot(br_hopf)
####################################################################################################Continuation of Periodic Orbit
# number of time slices
M = 51
l_hopf, Th, orbitguess2, hopfpt, vec_hopf = guessFromHopf(br, ind_hopf, opts_br_eq.newtonOptions.eigsolver, M, 2*2.7; phase = 0.25)
#
# orbitguess_f2 = orbitguess2[1]; for ii=2:M; global orbitguess_f2 = hcat(orbitguess_f2, orbitguess2[ii]);end
orbitguess_f2 = reduce(vcat, orbitguess2)
orbitguess_f = vcat(vec(orbitguess_f2), Th) |> vec

poTrap = PeriodicOrbitTrapProblem(Fbru, Jbru_sp, real.(vec_hopf), hopfpt.u, M, 2n)

poTrap(orbitguess_f,  @set par_bru.l = l_hopf + 0.01) |> plot
BK.plotPeriodicPOTrap(orbitguess_f, n, M; ratio = 2)

ls0 = GMRESIterativeSolvers(N = 2n, reltol = 1e-9)#, Pl = lu(I + par_cgl.Δ))

poTrapMF = PeriodicOrbitTrapProblem(
			Fbru, (x, p) -> (dx -> d1Fbru(x, p, dx)),
			real.(vec_hopf),
			hopfpt.u,
			M, ls0)

deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zero(orbitguess_f)])
# deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo_f])
####################################################################################################
opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 14)
# opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 14, linsolver = BK.LSFromBLS())
outpo_f, _, flag = @time newton(poTrap,
			orbitguess_f, (@set par_bru.l = l_hopf + 0.01),
			opt_po;
			# deflationOp,
			linearPO = :BorderedSparseInplace, #  3.2 seconds (756.76 k allocations: 2.363 GiB, 6.67% gc time)
			# linearPO = :BorderedLU,
			# linearPO = :FullSparseInplace, # 3.634019 seconds (775.86 k allocations: 2.452 GiB, 6.40% gc time)
			normN = norminf,
			callback = (x, f, J, res, iteration, itl, options; kwargs...) -> (println("--> amplitude = ", BK.amplitude(x, n, M; ratio = 2));true)
			)
	flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", BK.amplitude(outpo_f, n, M; ratio = 2),"\n")
	BK.plotPeriodicPOTrap(outpo_f, n, M; ratio = 2)

opt_po = @set opt_po.eigsolver = EigKrylovKit(tol = 1e-5, x₀ = rand(2n), verbose = 2, dim = 40)
opt_po = @set opt_po.eigsolver = DefaultEig()
# opt_po = @set opt_po.eigsolver = EigArpack(; tol = 1e-5, v0 = rand(2n))
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.1, ds= 0.01, pMax = 3.0, maxSteps = 20, newtonOptions = opt_po, saveSolEveryStep = 2, plotEveryStep = 5, nev = 11, precisionStability = 1e-6, nInversion = 4, detectBifurcation = 0, dsminBisection = 1e-6, maxBisectionSteps = 15)
	br_po, = @time continuation(poTrap,
			outpo_f, (@set par_bru.l = l_hopf + 0.01), (@lens _.l),
			opts_po_cont;
			# linearPO = :FullLU,
			updateSectionEveryStep = 1,
			# linearPO = :FullSparseInplace,
			linearPO = :BorderedSparseInplace,
			# tangentAlgo = BorderedPred(),
			verbosity = 3,	plot = true,
			# callbackN = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", BK.amplitude(x, n, M));true),
			# finaliseSolution = (z, tau, step, contResult; k...) ->
				# (Base.display(contResult.eig[end].eigenvals) ;true),
			plotSolution = (x, p;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...),
			# printSolution = (x, p;kwargs...) -> BK.amplitude(x, n, M; ratio = 2),
			normC = norminf)

####################################################################################################
using IncompleteLU
Jpo = @time poTrap(Val(:JacFullSparse), orbitguess_f, @set par_bru.l = l_hopf + 0.01)

@time poTrap(Val(:JacFullSparseInplace), Jpo, orbitguess_f, @set par_bru.l = l_hopf + 0.01)
_indx = BK.getBlocks(Jpo, 2n, M)
@time poTrap(Val(:JacFullSparseInplace), Jpo, orbitguess_f, (@set par_bru.l = l_hopf + 0.01), _indx)

@time lu(Jpo)
Jpo = @time poTrap(Val(:JacCyclicSparse), orbitguess_f, @set par_bru.l = l_hopf + 0.01)
@time lu(Jpo)

Precilu = @time ilu(Jpo, τ = 0.01)
Precilu = @time lu(Jpo)
ls = GMRESIterativeSolvers(verbose = false, reltol = 1e-5, N = size(Jpo,1), restart = 30, maxiter = 150, log=true, Pl = Precilu)
	@time ls(Jpo, rand(ls.N))

opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 20)
	outpo_f, = @time newton(poTrap,
			orbitguess_f, (@set par_bru.l = l_hopf + 0.01),
			(@set opt_po.linsolver = ls); linearPO = :BorderedMatrixFree,
			normN = norminf,
			# callback = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", amplitude(x), " T = ", x[end], ", T0 = ",orbitguess_f[end]);true)
			)
	printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", BK.amplitude(outpo_f, n, M; ratio = 2),"\n")
	BK.plotPeriodicPOTrap(outpo_f, n, M; ratio = 2)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.01, pMax = 2.2, maxSteps = 3000, newtonOptions = (@set opt_po.linsolver = ls))
	br_pok, = @time continuation(poTrap,
			outpo_f, (@set par_bru.l = l_hopf + 0.01), (@lens _.l),
			opts_po_cont; linearPO = :BorderedMatrixFree,
			verbosity = 2,
			plot = true,
			# callbackN = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", BK.amplitude(x, n, M));true),
			# plotSolution = (x, p;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...)
			normC = norminf)
####################################################################################################
# automatic branch switching from Hopf point
opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.04, ds = 0.01, pMax = 2.2, maxSteps = 200, newtonOptions = opt_po, saveSolEveryStep = 2,
	plotEveryStep = 1, nev = 11, precisionStability = 1e-6,
	detectBifurcation = 0, dsminBisection = 1e-6, maxBisectionSteps = 15, tolBisectionEigenvalue = 3., nInversion = 4)

M = 51
probFD = PeriodicOrbitTrapProblem(M = M)
br_po, _ = continuation(
	# arguments for branch switching
	jet..., br, 1,
	# arguments for continuation
	opts_po_cont, probFD;
	δp = 0.01,
	verbosity = 3,	plot = true, linearPO = :FullSparseInplace,
	finaliseSolution = (z, tau, step, contResult; k...) ->
		(Base.display(contResult.eig[end].eigenvals) ;true),
	plotSolution = (x, p; kwargs...) -> heatmap!(getTrajectory(p.prob, x, par_bru).u'; ylabel="time", color=:viridis, kwargs...),
	normC = norminf)

# branches = [br_po]
push!(branches, br_po)
plot(branches..., legend = :bottomright, plotfold = false)

####################################################################################################
# semi-automatic branch switching from bifurcation BP-PO
br_po2, _ = BK.continuationPOTrapBPFromPO(
	# arguments for branch switching
	br_po, 1,
	# arguments for continuation
	opts_po_cont;
	δp = 0.01,
	verbosity = 3,	plot = true, linearPO = :FullSparseInplace,
	finaliseSolution = (z, tau, step, contResult; k...) ->
		(Base.display(contResult.eig[end].eigenvals) ;true),
	plotSolution = (x, p; kwargs...) -> (plot!(br_po,legend = :bottomright, subplot=1)),
	normC = norminf)
push!(branches, br_po2)

plot();for _b in branches; plot!(_b; label = ""); end; title!("")

plot(branches[2])

br_po.functional

ext: jl


text:
"""
# Universal Evolutionary Model for Periodical Organisms

Accompanying Julia code for the paper "Universal evolutionary model for periodical organisms"

```
Eric Goles, Ivan Slapničar and Marco A. Lardies: Universal evolutionary model for periodical organisms, arXiv:2010.00940, submitted
```

The manuscript can be found [here](https://arxiv.org/abs/2010.00940).

The code in this file (also in the Pluto notebook `Simulation_p.jl`) was used to run all simulations and produce all figures in the paper.
"""

using LinearAlgebra, Plots, Plots.PlotMeasures, Random, Primes
Random.seed!(123)
import Plots.spy
spy(A)=heatmap(A, yflip=true, legend=false, c=cgrad([:white,:gray,:red]),
	   aspectratio=1,clim=(0.0,1.0))

mkpath("figures")

# Reference images
# gcd=1
n=20
P=zeros(Float64,n,n)
for i=2:n
	for j=2:n
	    if gcd(i,j)==1
	        P[i,j]=0.5
	    end
	end
end
figG=spy(P)
xlabel!("c1")
ylabel!("c2")

# Primes
P=zeros(Float64,n,n)
for i=2:n
    for j=2:n
        if isprime(i)
            P[i,j]=0.5
        end
    end
end
figP=spy(P)
xlabel!("c1")
ylabel!("c2")

# Divisors
P=zeros(Float64,n,n)
for i=2:n
    for j=2:n
        if mod(i,j)==0
            P[i,j]=0.5
        end
    end
end
figD=spy(P)
xlabel!("c1")
ylabel!("c2")

fᵣ=plot(figG,figP,figD,layout=(1,3),titlefontsize=10,xguidefontsize=10,
    yguidefontsize=10,left_margin=[0mm 0mm],bottom_margin=[10mm 10mm],
    top_margin=[0mm 0mm],right_margin=[0mm 0mm],size=(645,200),grid=false)
savefig(fᵣ,"./figures/Reference.png")

# Local fitness functions
f₁(χ₁,χ₂,v::Tuple)= (χ₁,χ₂)==(1,0) ? v[1] : (χ₁,χ₂)==(1,1) ? v[3] : 0
f₂(χ₁,χ₂,v)= (χ₁,χ₂)==(0,1) ? v[2] : (χ₁,χ₂)==(1,1) ? v[4] : 0
# This is special case for long-time masting bamboos
α=0.8
β=1/α
f₁(χ₁,χ₂,v::Vector)= (χ₁,χ₂)==(1,0) ? v[1] : (χ₁,χ₂)==(1,1) ? α : (χ₁,χ₂)==(0,1) ? β : 0

# Global fitness functions
# v is instance, c₁ and c₂ are starting periods of species C₁ and C₂
function F₁(c₁,c₂,v)
    T=c₁*c₂
    χ₁(t)=map(Int,iszero(mod(t,c₁)))
    χ₂(t)=map(Int,iszero(mod(t,c₂)))
    fitness=0
    for t=1:T
        fitness+=f₁(χ₁(t),χ₂(t),v)
    end
    fitness/c₂
end

function F₂(c₁,c₂,v)
    T=c₁*c₂
    χ₁(t)=map(Int,iszero(mod(t,c₁)))
    χ₂(t)=map(Int,iszero(mod(t,c₂)))
    fitness=0
    for t=1:T
        fitness+=f₂(χ₁(t),χ₂(t),v)
    end
    fitness/c₁
end

# Evolution dynamics algorithm
function Evolution(c₁,c₂,v,p)
    # v is instance, c₁ and c₂ are starting periods of species C₁ and C₂,
    # p is the tuple of possible period changes
    # Initial fitnesses
    F1=F₁(c₁,c₂,v)
    F2=F₂(c₁,c₂,v)
    F1new=copy(F1)
    F2new=copy(F2)
    E=Vector{Tuple{Int32,Int32,Number,Number}}(undef,0)
    push!(E,(c₁,c₂,F1,F2))
    change=true
    steps=0
    maxsteps=30
    while change || steps<maxsteps
        # Random change of species C₁
        d₁=c₁+rand(p)
        change=false
        if d₁>1
            F1new=F₁(d₁,c₂,v)
            # Switch to new period if fitness is better
            if F1new > F1
                c₁=d₁
                F1=F1new
                F2=F₂(c₁,c₂,v)
                push!(E,(c₁,c₂,F1,F2))
                change=true
            end
        end
        # Random change of species C₂
        d₂=c₂+rand(p)
        if d₂>1
            F2new=F₂(c₁,d₂,v)
            # Switch to new period if fitness is better
            if F2new > F2
                c₂=d₂
                F2=F2new
                F1=F₁(c₁,c₂,v)
                push!(E,(c₁,c₂,F1,F2))
                change=true
            end
        end
        steps+=1
    end
    E
    return (E[end][1],E[end][2])
end

# Define fitness tuples v:
v=[(-1, 0, 1, 1), (-1, 1, 1, 1), (0, 1, 1, 1),
(0, 0, 1, 1), (-1, -1, 1, 1), (1, 1, -1, -1),
(1, 1, -1, 0), (1, 1, 1, -1), (-1, 1, 1, -1),
(-1, 1, 1, 0), (1, 0, -1, 0), (0, 1, 1, -1),
(0, 1, 1, 0), (1, 1, 1, 0), (0, -1, 1, 0),
(-1, 0, 1, 0), (1, -1, -1, 0), (1, -1, 1, 0)]

# Plots of attractors
function myPlots(n::Int, k::Int, v::Vector, quantitative::Bool=false)
    # n is size of plot
    # k defines the evolution neighbourhood p
	# v vector of tuples
    # quantitative==false means to use qualitative model (-1,0,1)

    # Define neighbourhood p
    q=collect(-k:k)
    deleteat!(q,k+1)
    p=Tuple(1q)
    # String for title
    stringp=", p=[-"*string(k)*","*string(k)*"]"

    Ev=Array{Tuple{Int32,Int32}}(undef,n,n)
    for i=1:n
        Ev[1,i]=(0,0)
        Ev[i,1]=(0,0)
    end
    # Floating-point array to plot
    FpPlot=zeros(Float64,n,n)
    # Array for figures
    fig=Array{Any}(undef,length(v))
    l=1
    for i=1:length(v)
        FpPlot=zeros(Float64,n,n)
        ν=v[i]
        if quantitative==true
			μ=Tuple([rand(1:10) for j=1:length(ν)])
            ν=μ.*v[i]
        end
        for c₁=2:n
            for c₂=2:n
                Evo=Evolution(c₁,c₂,ν,p)
                Ev[c₁,c₂]=Evo
                if Evo==(c₁,c₂)
                    if F₁(c₁,c₂,ν)>=0 && F₂(c₁,c₂,ν)>=0
                        # Both species survive
                        FpPlot[c₂,c₁]=1.0
                    else
                        # Someone dies
                        if F₁(c₁,c₂,ν)<0
                            FpPlot[c₂,c₁]=0.3
                        else
                            FpPlot[c₂,c₁]=0.3
                        end
                    end
                end
            end
        end
        # Plot
        fig[l]=spy(FpPlot)
        xlabel!("c1")
        ylabel!("c2")
        title!("v["*string(i)*"]="*string(ν)*stringp)
        l+=1
    end
    return fig
end

# Attractors for neigbourhood [-1,1]
fig1=myPlots(20,1,v)

# Attractors for neigbourhood [-4,4]
fig4=myPlots(20,4,v)

# Quantitative plots for neigbourhood [-1,1]
figQ=myPlots(20,1,v,true)

# Figures for the paper
for k=1:length(v)
    f1=fig1[k]
    f4=fig4[k]
    fQ=figQ[k]
    f=plot(f1,f4,fQ,layout=(1,3),titlefontsize=10,xguidefontsize=10,yguidefontsize=10,
        left_margin=[0mm 0mm],bottom_margin=[10mm 10mm],top_margin=[0mm 0mm],
        right_margin=[0mm 0mm],size=(645,220),grid=false)
    savefig(f,"./figures/S"*string(k)*".png")
end

# Long-period masting bamboos
# Redefine fitness vector
v₀=[[-1, 1, α, 1, β, 0, 0, 0]]
# Time period
# n=40
# k defines the evolution neighbourhood p
# k=6
# Attractors
figB=myPlots(40,6,v₀)
# Figure
fb=plot(figB[1],titlefontsize=10,xguidefontsize=10,yguidefontsize=10,
    left_margin=[0mm 0mm],bottom_margin=[0mm 0mm],top_margin=[0mm 0mm],
    right_margin=[0mm 0mm],size=(400,400),grid=:show)
savefig(fb,"./figures/Bamboos.png")


text:
using Documenter, Bitcoin

makedocs(sitename="Bitcoin",
         doctest=true)


text:
module ExponentialUtilities
using LinearAlgebra, SparseArrays, Printf, Requires
using ArrayInterface: ismutable

"""
    @diagview(A,d) -> view of the `d`th diagonal of `A`.
"""
macro diagview(A,d::Integer=0)
    s = d<=0 ? 1+abs(d) : :(m+$d)
    quote
        m = size($(esc(A)),1)
        @view($(esc(A))[($s):m+1:end])
    end
end

include("exp.jl")
include("exp_baseexp.jl")
include("exp_noalloc.jl")
include("phi.jl")
include("arnoldi.jl")
include("krylov_phiv.jl")
include("krylov_phiv_adaptive.jl")
include("kiops.jl")
include("StegrWork.jl")
include("krylov_phiv_error_estimate.jl")

export phi, phi!, KrylovSubspace, arnoldi, arnoldi!, lanczos!, ExpvCache, PhivCache,
    expv, expv!, exp_generic, phiv, phiv!, kiops, expv_timestep, expv_timestep!, phiv_timestep, phiv_timestep!,
    StegrCache, get_subspace_cache

function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        function ExponentialUtilities.expv!(w::CuArrays.CuVector{Tw},
                                    t::Real, Ks::KrylovSubspace{T, U};
                                    cache=nothing,
                                    dexpHe::CuArrays.CuVector = CuArrays.CuVector{U}(undef, Ks.m)) where {Tw, T, U}
            m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
            @assert length(w) == size(V, 1) "Dimension mismatch"
            if cache == nothing
                cache = Matrix{U}(undef, m, m)
            elseif isa(cache, ExpvCache)
                cache = get_cache(cache, m)
            else
                throw(ArgumentError("Cache must be an ExpvCache"))
            end
            copyto!(cache, @view(H[1:m, :]))
            if ishermitian(cache)
                # Optimize the case for symtridiagonal H
                F = eigen!(SymTridiagonal(cache))
                expHe = F.vectors * (exp.(lmul!(t,F.values)) .* @view(F.vectors[1, :]))
            else
                lmul!(t, cache); expH = cache
                _exp!(expH)
                expHe = @view(expH[:, 1])
            end

            copyto!(dexpHe, expHe)
            lmul!(beta, mul!(w, @view(V[:, 1:m]), dexpHe)) # exp(A) ≈ norm(b) * V * exp(H)e
        end
    end

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        function ExponentialUtilities.expv!(w::CUDA.CuVector{Tw},
                                    t::Real, Ks::KrylovSubspace{T, U};
                                    cache=nothing,
                                    dexpHe::CUDA.CuVector = CUDA.CuVector{U}(undef, Ks.m)) where {Tw, T, U}
            m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
            @assert length(w) == size(V, 1) "Dimension mismatch"
            if cache == nothing
                cache = Matrix{U}(undef, m, m)
            elseif isa(cache, ExpvCache)
                cache = get_cache(cache, m)
            else
                throw(ArgumentError("Cache must be an ExpvCache"))
            end
            copyto!(cache, @view(H[1:m, :]))
            if ishermitian(cache)
                # Optimize the case for symtridiagonal H
                F = eigen!(SymTridiagonal(cache))
                expHe = F.vectors * (exp.(lmul!(t,F.values)) .* @view(F.vectors[1, :]))
            else
                lmul!(t, cache); expH = cache
                _exp!(expH)
                expHe = @view(expH[:, 1])
            end

            copyto!(dexpHe, expHe)
            lmul!(beta, mul!(w, @view(V[:, 1:m]), dexpHe)) # exp(A) ≈ norm(b) * V * exp(H)e
        end
    end
end

end


text:
using HOODESolver
using LinearAlgebra
using Random
using Test

function fct4( du, u, p, t)
    B=[0.1 1.1 0.2 05 ; -0.9 -0.12 -0.7 0.4 ; 0.5 0.66 0.7 0.8 ; -0.34 0.8 0 0.3]
    du[:] = B*u
    missing
end

function testinterface_fct()
    @time @testset "test interface type of function" begin
        A =  [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
        B=[0.1 1.1 0.2 05 ; -0.9 -0.12 -0.7 0.4 ; 0.5 0.66 0.7 0.8 ; -0.34 0.8 0 0.3]
        u0 = 2rand(4)-ones(4)
        epsilon = 0.0001
        sol_ref = exp(1.0*(1/epsilon*A+B))*u0
        fct = (u,p,t) -> B*u
        prob = HOODEProblem(fct, u0, (0.0, 1.0), missing, A, epsilon)
        sol = solve(prob)
        @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
        fct = (u,p) -> B*u
        prob = HOODEProblem(fct, u0, (0.0, 1.0), missing, A, epsilon)
        sol = solve(prob)
        @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
        fct = (u) -> B*u
        prob = HOODEProblem(fct, u0, (0.0, 1.0), missing, A, epsilon)
        sol = solve(prob)
        @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
        prob = HOODEProblem(fct4, u0, (0.0, 1.0), missing, A, epsilon)
        sol = solve(prob)
        @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
    end
end

function testinterface_epsilon()
    @time @testset "test interface while epsilon varying" begin
        Random.seed!(199881)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        fct = (u,p,t) -> B*u
        order = 5
        for i=1:1
            epsilon = big"1.0"/big"10"^(i*3)
            eps_v = convert(Float32, epsilon)
            t_max = big"1.0"
            sol_ref = exp(t_max*(1/epsilon*A+B))*u0
            prob = HOODEProblem(fct, u0, (big"0.0", t_max), missing, A, epsilon)
            println("epsilon=$eps_v sol_ref=$sol_ref")
            nb = 100
            res_err = zeros(BigFloat,5)
            ind = 1
            while nb <= 1000
                sol = solve(prob, nb_t=nb, order=5, dense=false)
                res_err[ind] = norm(sol[end]-sol_ref,Inf)
 #               println("\nnb=$nb epsilon=$eps_v err=$(res_err[ind])")
                nb *= 10
                ind += 1
            end
            coef = Float64(log(res_err[1]) - log(res_err[2]))/log(10)
            @test isapprox( order, coef, atol=0.01)
            println("epsilon=$eps_v coef log = $coef")
        end
    end
end
function testinterface_interpolate()
    @time @testset "test interface for interpolation" begin
        seed=123871
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        order = 6
        nb = 100
        t_max = big"1.0"
        for i=1:1
            seed = seed*1341+6576123
            Random.seed!(seed)
            u0=rand(BigFloat,4)
            println("u0=$u0")
            B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
            fct = (u,p,t) -> B*u
            epsilon = big"0.4"/2.0^i
            eps_v = convert(Float32, epsilon)
            t_max = big"1.0"
            prob = HOODEProblem(fct, u0, (big"0.0",t_max), missing, A, epsilon)
   #         sol = solve(prob, getprecision=false)
            sol = solve(prob)
            m = 1/epsilon*A+B
            reftol=norm(exp(t_max*m)*u0-sol[end], Inf)*10
            for j=1:10
                t=rand(BigFloat)
                res_ex=exp(t*m)*u0
                res_ap=sol(t)
                println("t=$t norm=$(norm(res_ex-res_ap, Inf))")
                @test isapprox(res_ex, res_ap, atol=reftol, rtol=reftol*10)
            end
        end
    end
end
function testinterface_interpolate_float()
    @time @testset "test interface for interpolation" begin
        seed=124333
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        order = 6
        nb = 100
        t_max = 1.0
        for i=1:1
            seed = seed*1341+6576123
            Random.seed!(seed)
            u0=rand(4)
            println("u0=$u0")
            B = 2rand(4,4)-ones(4,4)
            fct = (u,p,t) -> B*u
            epsilon = 0.4/2.0^i
            eps_v = convert(Float32, epsilon)
            t_max = 1.0
            prob = HOODEProblem(fct, u0, (0.0,t_max), missing, A, epsilon)
   #         sol = solve(prob, getprecision=false)
            sol = solve(prob)
            m = 1/epsilon*A+B
            reftol=norm(exp(t_max*m)*u0-sol[end], Inf)*10
            @test reftol < 1e-6
            for j=1:10
                t=rand()
                res_ex=exp(t*m)*u0
                res_ap=sol(t)
#                println("type ex=$(typeof(res_ex)) type ap=$(typeof(res_ap))")
#                println("t=$t norm=$(norm(res_ex-res_ap, Inf))")
                @test isapprox(res_ex, res_ap, atol=reftol, rtol=reftol*10)
            end
        end
    end
end
function testinterface_short()
    seed=981885
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 7
    nb = 5
     for i=1:10
        seed += 10
        Random.seed!(seed)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        fct = (u,p,t) -> B*u
        epsilon = big"0.00004"/2.0^i
        eps_v = convert(Float32, epsilon)
        t_max = big"0.01"
        prob = HOODEProblem(fct, u0, (big"0.0",t_max), missing, A, epsilon)
        sol = solve(prob, getprecision=false, order=order, nb_t=nb)
        resnorm=norm(exp(t_max*(1/epsilon*A+B))*u0-sol[end], Inf)
        println("sol=$(sol[end])")
        println("solexact=$(exp(t_max*(1/epsilon*A+B))*u0)")
        println("resnorm=$resnorm")
        @test resnorm < 1e-10
    end
end
function tts_time(t_begin, t_end)
    @time @testset "test interface time from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = rand(BigFloat, 4)
        fct = (u,p,t) -> B*u
        epsilon=big"0.000000345"
        prob = HOODEProblem(fct, u0, (t_begin,t_end), missing, A, epsilon)
        sol = solve(prob, getprecision=false, nb_t=100, order=4)
        m = 1/epsilon*A+B
        solref = exp((t_end-t_begin)*m)*u0
        println("sol=$(sol[end]) solref=$solref norm=$(norm(sol[end]-solref,Inf))")
        @test isapprox(sol[end], solref,atol=1e-8, rtol=1e-7)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=exp((t-t_begin)*m)*u0
            res_ap=sol(t)
            @test isapprox(res_ex, res_ap, atol=1e-8, rtol=1e-7)
        end
    end
end
function tts_time_time(t_begin, t_end)
    @time @testset "test interface time time from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = 2rand(BigFloat, 4)-ones(BigFloat,4)
        fct = (u,p,t) -> B*u +t*p[1] +p[2]
        tuple_p = (2rand(BigFloat,4)-ones(BigFloat,4),
    2rand(BigFloat,4)-ones(BigFloat,4))
        epsilon=big"0.000003467"
        nb = 100
        order = 4
        prob = HOODEProblem(fct, u0, (t_begin,t_end), tuple_p, A, epsilon, B)
        sol = solve(prob, getprecision=false, nb_t=100, order=4)
        solref = getexactsol(sol.par_u0.parphi, u0, t_end)
        println("sol=$(sol[end]) solref=$solref norm=$(norm(sol[end]-solref,Inf))")
        @test isapprox(sol[end], solref,atol=1e-8, rtol=1e-7)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=getexactsol(sol.par_u0.parphi, u0, t)
            res_ap=sol(t)
            @test isapprox(res_ex, res_ap, atol=1e-8, rtol=1e-7)
        end
    end
end
function testinterface_time()
    seed=7887
    Random.seed!(seed)
    tts_time(big"0.0",big"1.0")
    tts_time(big"0.0",big"1.4565656565")
    tts_time(big"2.234566",big"3.45766790123")
    tts_time(-big"0.8457676",big"0.56716")
    tts_time(-big"1.8111457676",-big"0.345456716")
    tts_time(big"1.8457676",big"0.56716")
    tts_time(-big"0.8457676",-big"1.56716")
    tts_time(-big"0.845722676",-big"0.56716")
end
function testinterface_time_time()
    seed=788227
    Random.seed!(seed)
    tts_time_time(big"0.0",big"1.0")
    tts_time_time(big"0.0",big"1.4565656565")
    tts_time_time(big"11.234566",big"12.45766790123")
    tts_time_time(-big"0.8457676",big"0.56716")
    tts_time_time(-big"1.8111457676",-big"0.345456716")
    tts_time_time(big"1.8457676",big"0.56716")
    tts_time_time(-big"0.8457676",-big"1.56716")
    tts_time_time(-big"0.845722676",-big"0.56716")
end

testinterface_fct()
testinterface_interpolate_float()
testinterface_time()
testinterface_time_time()
testinterface_epsilon()
testinterface_interpolate()
testinterface_short()
# testtwoscales_pure_ab()
# testtwoscales_pure_ab2()
#testtwoscales_pure_ab3()


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

@static if VERSION >= v"0.7.0-DEV.2005"
    FLOATMAX(v) = floatmax(v)
else
    FLOATMAX(v) = realmax(v)
end

"""
    initInfo = InitInfo(...)

Data structure that is returned when calling initialize!(...)
"""
struct InitInfo
    y0::Vector{Float64}          # Initial y
    yp0::Vector{Float64}         # Initial yp0
    y_nominal::Vector{Float64}   # Nominal values of y (y_nominal[i]>0.0 required)
    y_errorControl::Vector{Bool} # = true: error control; = false, no error control
    maxTime::Float64             # Integrate at most until maxTime
    nextEventTime::Float64       # Next time event
    integrateToEvent::Bool       # = true, if exact integration to nextEventTime
    terminate::Bool              # = true, terminate simulation

    function InitInfo(y0::Vector{Float64},
                      yp0::Vector{Float64};
                      y_nominal::Vector{Float64}=ones(length(y0)),
                      y_errorControl::Vector{Bool}=fill(true, length(y0)),
                      maxTime::Float64=FLOATMAX(Float64),
                      nextEventTime::Float64=FLOATMAX(Float64),
                      integrateToEvent::Bool=true,
                      terminate::Bool=false)
        @assert(length(y0) == length(yp0))
        @assert(length(y0) == length(y_nominal))
        @static if VERSION >= v"0.7.0-DEV.2005"
            @assert(y_nominal[ argmin(y_nominal) ] > 0.0)
        else
            @assert(y_nominal[ indmin(y_nominal) ] > 0.0)
        end
        new(y0, yp0, y_nominal, y_errorControl, maxTime, nextEventTime, integrateToEvent, terminate)
    end
end

"""
    eventInfo = EventInfo(...)

Data structure in which the return values of processEvent!(...) are stored
"""
mutable struct EventInfo
    restart::EventRestart
    maxTime::Float64         # Integrate at most until maxTime
    nextEventTime::Float64   # Next time event
    integrateToEvent::Bool   # = true, if exact integration to nextEventTime
    nextModel                # If restart==FullRestart, next model

    EventInfo() = new(Restart, FLOATMAX(Float64), FLOATMAX(Float64), true, nothing)
end

function reset!(einfo::EventInfo)
    einfo.restart = Restart
    einfo.maxTime = FLOATMAX(Float64)
    einfo.nextEventTime = FLOATMAX(Float64)
    einfo.integrateToEvent = true
    einfo.nextModel = nothing
end


text:
module InfluxFlux

using HTTP
using CSV
using Dates
using DataFrames
using TimeZones

export time_spec_to_epoc_ns, influx_server, flux, flux_to_dataframe,
       measurement, aggregate_measurement, measurements, buckets

TimeSpec = Union{Int,DateTime,ZonedDateTime}

struct InfluxServer
  uri::String
  org::String
  api_token::String
end

function time_spec_to_epoc_ns(time_spec::Int)
  time_spec
end

function time_spec_to_epoc_ns(time_spec::ZonedDateTime)
  time_spec_to_epoc_ns(DateTime(time_spec, UTC))
end

function time_spec_to_epoc_ns(time_spec::DateTime)
  Int(1_000_000_000 * datetime2unix(time_spec))
end

function uri_helper(srv::InfluxServer, path::String)
  HTTP.URI(HTTP.URI("$(srv.uri)/$path"), query = Dict("org" => srv.org))
end

function token_json_headers(srv::InfluxServer)
  Dict("Authorization" => "Token $(srv.api_token)"
       , "Accept" => "application/json")
end

function influx_server(uri::String, org::String, api_token::String)::InfluxServer
  InfluxServer(uri, org, api_token)
end

function flux(srv::InfluxServer, flux_query::String)
  headers = merge(token_json_headers(srv), Dict("Content-Type" => "application/vnd.flux"))
  response = HTTP.post(uri_helper(srv, "api/v2/query"), headers, flux_query);
  if response.status == 200
    response.body
  else
    throw(response.status)
  end
end

function flux_to_dataframe(srv::InfluxServer, flux_query::String)
  CSV.File(flux(srv, flux_query), delim = ',') |> DataFrame
end

function measurement(srv::InfluxServer, bucket::String, measurement_name::String, from::TimeSpec, to::TimeSpec)
  q = """
    from(bucket: "$bucket")
    |> range(start: time(v: uint(v: $(time_spec_to_epoc_ns(from)))), stop: time(v: uint(v: $(time_spec_to_epoc_ns(to)))))
      |> filter(fn: (r) => r._measurement == "$measurement_name")
      |> pivot(rowKey:["_time"], columnKey: ["_field"], valueColumn: "_value")
      |> map(fn: (r) => ({ r with _time: uint(v: r._time) }))
  """
  #println(q)
  result = flux_to_dataframe(srv, q)
  result[:,7:end] # remove influxdb internal columns
end

function aggregate_measurement(srv::InfluxServer, bucket::String, measurement_name::String, from::TimeSpec, to::TimeSpec, window::Period; fn::String = "mean")
  q = """
    from(bucket: "$bucket")
    |> range(start: time(v: uint(v: $(time_spec_to_epoc_ns(from)))), stop: time(v: uint(v: $(time_spec_to_epoc_ns(to)))))
      |> filter(fn: (r) => r._measurement == "$measurement_name")
      |> aggregateWindow(every: $(Nanosecond(window).value)ns, fn: $fn, createEmpty: false)
      |> pivot(rowKey:["_time"], columnKey: ["_field"], valueColumn: "_value")
      |> map(fn: (r) => ({ r with _time: uint(v: r._time) }))
  """
  #println(q)
  result = flux_to_dataframe(srv, q)
  result[:,7:end] # remove influxdb internal columns
end

function measurements(srv::InfluxServer, bucket::String)
  q = """
    import "influxdata/influxdb/schema"
    schema.measurements(bucket: "$bucket")
  """
  String.(InfluxFlux.flux_to_dataframe(srv, q)[:, "_value"])
end

function buckets(srv::InfluxServer)
  String.(InfluxFlux.flux_to_dataframe(srv, "buckets()")[:, :name])
end

end # module

max_stars_repo_name: tallakt/InfluxFlux

max_issues_repo_name: tallakt/InfluxFlux

max_forks_repo_name: tallakt/InfluxFlux


text:
using Documenter, DocumenterCitations, FastDEQ

bib = CitationBibliography(joinpath(@__DIR__, "ref.bib"), sorting = :nyt)
include("pages.jl")

makedocs(
    bib,
    sitename = "Fast Deep Equilibrium Networks",
    authors = "Avik Pal et al.",
    clean = true,
    doctest = false,
    modules = [FastDEQ],
    format = Documenter.HTML(analytics = "UA-90474609-3",
                             #  assets = ["assets/favicon.ico"],
                             canonical="https://fastdeq.sciml.ai/stable/"),
    pages = pages
)

deploydocs(
   repo = "github.com/SciML/FastDEQ.jl.git";
   push_preview = true
)


text:
module QJuliaCGPCG

using QJuliaInterface
using QJuliaEnums
using QJuliaBlas
using QJuliaReduce
using QJuliaSolvers

using LinearAlgebra
using Printf

norm2    = QJuliaReduce.gnorm2
axpyZpbx = QJuliaBlas.axpyZpbx
rdot     = QJuliaReduce.reDotProduct

function solver(x::AbstractArray, b::AbstractArray, Mat::Any, MatSloppy::Any, param::QJuliaSolvers.QJuliaSolverParam_qj, K::Function)

    println("Running Chronopoulos Gear PCG solver.")

    if (param.maxiter == 0)
      if param.use_init_guess == false
        x .=@. 0.0
      end
      return
    end #if param.maxiter == 0

    mixed = (param.dtype_sloppy != param.dtype)

    if mixed == true; error("Mixed types is not supported."); end

    local r   = Vector{param.dtype_sloppy}(undef, length(x))
    # now allocate sloppy fields
    local p       = typeof(r)(undef, length(r))
    local s       = typeof(r)(undef, length(r))
    local u       = typeof(r)(undef, length(r))
    local w       = typeof(r)(undef, length(r))

    b2 = norm2(b)  #Save norm of b
    r2 = 0.0     #if zero source then we will exit immediately doing no work

    if param.use_init_guess == true
      #r = b - Ax0 <- real
      Mat(r, x)
      r .=@. b - r
      r2 = norm2(r)
    else
      r2 = b2
      r .=@. b
      x .=@. 0.0
    end

    K(u, r)

    Mat(w,u)

    # if invalid residual then convergence is set by iteration count only
    stop = b2*param.tol*param.tol

    println("CGPCG: Initial residual = ", sqrt(r2))

    converged = false; k = 0

    γ_old = 0.0
    α_old = 0.0

    case = 2
    while (k < param.maxiter && converged == false)

      γ = rdot(r,u)
      δ = rdot(u,w)

      if k > 0
        β = γ / γ_old
        α = γ / (δ-(β*γ) / α_old)
      else
        β = 0.0
        α = γ / δ
      end
      γ_old = γ
      α_old = α

      p .=@. u + β*p
      s .=@. w + β*s

      x .=@. x + α*p
      r .=@. r - α*s

      K(u, r)
      Mat(w, u)

      converged = (γ > stop) ? false : true

      @printf("CGPCG: %d iteration, iter residual: %le \n", k, sqrt(γ))

      k += 1
    end #while

    if (param.compute_true_res == true)
      Mat(r, x)

      r .=@. b - r
      r2 = norm2(r)

      param.true_res = sqrt(r2 / b2)
      println("CGPCG: converged after ", k , "  iterations, relative residual: true = ", sqrt(r2))

    end #if (param.compute_true_res == true)

end #solver

end #QJuliaCGPCG

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
mutable struct UpdateTimer
    warm_up::Int
    update_wait::Int
    t::Int
    UpdateTimer(warm_up, update_wait) = new(warm_up, update_wait, 0)
end

(ut::UpdateTimer)(replay_len) =
    replay_len >= ut.warm_up && (ut.t % ut.update_wait == 0)

next!(ut::UpdateTimer) = ut.t += 1
restart!(ut::UpdateTimer) = ut.t = 0


text:
#
# Correctness Tests
#

# require("test.jl")

using Stats

my_tests = ["test/01.jl"]
            #"test/statquiz.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
module DiffUtils

using diffutils_jll
using FIFOStreams

@noinline function check_argtype(name::Symbol, val, T::Type)
    @nospecialize
    val isa T || throw(ArgumentError(
        "expected option `$name` to be of type `$T`, but got $val::$(typeof(val)) instead.",
    ))
    nothing
end

const _BOOLEAN_ARGS = [
    "-q", "--brief",
    "-s", "--report-identical-files",
    "-e", "--ed",
    "-n", "--rcs",
    "-y", "--side-by-side",
    "--left-column",
    "--suppress-common-lines",
    "-p", "--show-c-function",
    "-t", "--expand-tabs",
    "-T", "--initial-tab",
    "--suppress-blank-empty",
    "-l", "--paginate",
    "-r", "--recursive",
    "--no-dereference",
    "-N", "--new-file",
    "--unidirectional-new-file",
    "--ignore-file-name-case",
    "--no-ignore-file-name-case",
    "-i", "--ignore-case",
    "-E", "--ignore-tab-expansion",
    "-Z", "--ignore-trailing-space",
    "-b", "--ignore-space-change",
    "-w", "--ignore-all-space",
    "-B", "--ignore-blank-lines",
    "-a", "--text",
    "--strip-trailing-cr",
    "-d", "--minimal",
    "--speed-large-files",
]
const BOOLEAN_ARGS = Dict(
    Symbol(replace(strip(arg, '-'), '-' => '_')) => arg
    for arg in _BOOLEAN_ARGS
)

const _INTEGER_ARGS = [
    "-c", "-C", "--context",
    "-u", "-U", "--unified",
    "-W", "--width",
    "--tabsize",
    "--horizon-lines",
]
const INTEGER_ARGS = Dict(
    Symbol(replace(strip(arg, '-'), '-' => '_')) => arg
    for arg in _INTEGER_ARGS
)

function parse_options(options::Dict{Symbol,Any})
    args = String[]
    for (name, val) in pairs(options)
        if haskey(BOOLEAN_ARGS, name)
            check_argtype(name, val, Bool)
            val && push!(args, BOOLEAN_ARGS[name])
        elseif haskey(INTEGER_ARGS, name)
            check_argtype(name, val, Int)
            argname = INTEGER_ARGS[name]
            if occursin(r"^\-[a-z]$", argname)
                push!(args, "$argname$val")
            elseif occursin(r"^\-[A-Z]$", argname)
                push!(args, argname)
                push!(args, string(val))
            else
                push!(args, "$argname=$val")
            end
        elseif name === :color
            check_argtype(name, val, Union{Bool, AbstractString, Symbol})
            if val isa Bool
                val && push!(args, "--color")
            else
                push!(args, "--color=$val")
            end
        elseif name === :labels
            check_argtype(name, val, NTuple{2, Union{AbstractString, Symbol}})
            push!(args, "--label=$(val[1])")
            push!(args, "--label=$(val[2])")
        end
    end
    return args
end

const DEFAULT_OPTIONS = Dict(
    :side_by_side => true,
    :color => true,
    :labels => ("stream 1", "stream 2"),
)

"""
    DiffUtils.diff(f::Function; stdout=stdout, side_by_side=true, color=true, options...)

Takes a function `f` which gets passed two output streams and may write to them. The resulting
contents are then passed to `diff`.

Optionally, `stdout` specifies the stream which the output of `diff` gets written to.

```jldoctest
julia> DiffUtils.diff() do s1, s2
           println(s1, "foo")
           println(s2, "bar")
       end
foo							      |	bar
```

## Supported Options

See the [`diff` man page](https://www.man7.org/linux/man-pages/man1/diff.1.html) for further
reference.

### Boolean Arguments

$(join(string.("`", keys(BOOLEAN_ARGS), "`"), ", "))

### Integer Arguments

$(join(string.("`", keys(INTEGER_ARGS), "`"), ", "))

### Other

$(join(string.("`", [:color, :labels], "`"), ", "))
"""
function diff(f::Function; stdout=stdout, @nospecialize(options...))
    _diff() do diff
        s = FIFOStreamCollection(2)
        local ret
        try
            args = parse_options(merge(DEFAULT_OPTIONS, Dict{Symbol, Any}(pairs(options))))
            cmd = Cmd(String[diff; path(s, 1); path(s, 2); args])
            attach(s, pipeline(ignorestatus(cmd); stdout=stdout))
            s1, s2 = s
            ret = f(s1, s2)
        finally
            close(s)
        end
        return ret
    end
end

"""
    DiffUtils.diff(x1, x2; stdout=stdout, side_by_side=true, color=true, options...)

Prints the diff of the string representations of `x1` and `x2`.

Optionally, `stdout` specifies the stream which the output of `diff` gets written to.

```jldoctest
julia> DiffUtils.diff("foo", "bar")
foo							      |	bar
```

For further details on the supported options, see [`diff(::Function)`](@ref).
"""
function diff(x1, x2; stdout=stdout, @nospecialize(options...))
    diff(; stdout=stdout, options...) do s1, s2
        print(s1, x1)
        println(s1)
        print(s2, x2)
        println(s2)
    end
    nothing
end

end


text:
@testset "Rotate90" begin
    @test typeof(@inferred(Rotate90())) <: Rotate90 <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test @inferred(Rotate90(0.7)) === Either(Rotate90(), 0.7)
        @test str_show(Rotate90()) == "Augmentor.Rotate90()"
        @test str_showconst(Rotate90()) == "Rotate90()"
        @test str_showcompact(Rotate90()) == "Rotate 90 degree"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Rotate90(), nothing)
        @test Augmentor.supports_eager(Rotate90) === true
        imgs = [
            (rect),
            (view(rect, :, :)),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for (img_in) in imgs
                res = @inferred(Augmentor.applyeager(Rotate90(), img_in))
                @test res == rotl90(rect)
                @test typeof(res) <: Array
            end
            res1, res2 = @inferred(Augmentor.applyeager(Rotate90(), (square, square2)))
            @test res1 == rotl90(square)
            @test res2 == rotl90(square2)
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Rotate90) === true
        @test_throws MethodError Augmentor.applyaffine(Rotate90(), nothing)
        @test @inferred(Augmentor.toaffinemap(Rotate90(), rect)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        @test @inferred(Augmentor.toaffinemap(Rotate90(), rect, 90)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        wv = @inferred Augmentor.applyaffine(Rotate90(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        # check that the affine map is computed for each image
        res1, res2 = @inferred(Augmentor.applyaffine(Rotate90(), (Augmentor.prepareaffine(square), Augmentor.prepareaffine(square2))))
        @test res1 == OffsetArray(rotl90(square), 0, 0)
        @test res2 == OffsetArray(rotl90(square2), 0, 0)
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Rotate90) === true
        @test_throws MethodError Augmentor.applyaffineview(Rotate90(), nothing)
        wv = @inferred Augmentor.applyaffineview(Rotate90(), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == rotl90(square)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Rotate90) === true
        v = @inferred Augmentor.applylazy(Rotate90(), rect)
        @test v === view(PermutedDimsArray(rect, (2,1)), 3:-1:1, 1:1:2)
        @test v == rotl90(rect)
        @test typeof(v) <: SubArray
        v = @inferred Augmentor.applylazy(Rotate90(), view(square,1:2,1:3))
        @test v === view(PermutedDimsArray(square, (2,1)), 3:-1:1, 1:1:2)
        @test v == rotl90(rect)
        @test typeof(v) <: SubArray
        wv = @inferred Augmentor.applylazy(Rotate90(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        res1, res2 = @inferred(Augmentor.applylazy(Rotate90(), (square, square2)))
        @test res1 == rotl90(square)
        @test res2 == rotl90(square2)
    end
    @testset "view" begin
        @test Augmentor.supports_view(Rotate90) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Rotate90) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Rotate90) === true
        v = @inferred Augmentor.applypermute(Rotate90(), rect)
        @test v === view(PermutedDimsArray(rect, (2,1)), 3:-1:1, 1:1:2)
        @test v == rotl90(rect)
        @test typeof(v) <: SubArray
        v2 = @inferred Augmentor.applypermute(Rotate90(), v)
        @test v2 === view(rect, 2:-1:1, 3:-1:1)
        @test v2 == rot180(rect)
    end
end

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

@testset "Rotate180" begin
    @test typeof(@inferred(Rotate180())) <: Rotate180 <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test @inferred(Rotate180(0.7)) === Either(Rotate180(), 0.7)
        @test str_show(Rotate180()) == "Augmentor.Rotate180()"
        @test str_showconst(Rotate180()) == "Rotate180()"
        @test str_showcompact(Rotate180()) == "Rotate 180 degree"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Rotate180(), nothing)
        @test Augmentor.supports_eager(Rotate180) === true
        imgs = [
            (rect),
            (view(rect, :, :)),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for (img_in) in imgs
                res = @inferred(Augmentor.applyeager(Rotate180(), img_in))
                @test res == rot180(rect)
                @test typeof(res) <: Array
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Rotate180) === true
        @test_throws MethodError Augmentor.applyaffine(Rotate180(), nothing)
        @test @inferred(Augmentor.toaffinemap(Rotate180(), rect)) ≈ AffineMap([-1.0 -1.22465e-16; 1.22465e-16 -1.0], [3.0,4.0])
        wv = @inferred Augmentor.applyaffine(Rotate180(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv[1:3,1:3] == rot180(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Rotate180) === true
        @test_throws MethodError Augmentor.applyaffineview(Rotate180(), nothing)
        wv = @inferred Augmentor.applyaffineview(Rotate180(), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv[1:3,1:3] == rot180(square)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Rotate180) === true
        v = @inferred Augmentor.applylazy(Rotate180(), rect)
        @test v === view(rect, 2:-1:1, 3:-1:1)
        @test v == rot180(rect)
        @test typeof(v) <: SubArray
        wv = @inferred Augmentor.applylazy(Rotate180(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv[1:3,1:3] == rot180(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(Rotate180) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Rotate180) === true
        v = @inferred Augmentor.applylazy(Rotate180(), rect)
        @test v === view(rect, 2:-1:1, 3:-1:1)
        @test v == rot180(rect)
        @test typeof(v) <: SubArray
        img = OffsetArray(rect,-2,1)
        v = @inferred Augmentor.applylazy(Rotate180(), img)
        @test v === view(img, 0:-1:-1, 4:-1:2)
        @test v == rot180(rect)
        @test typeof(v) <: SubArray
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Rotate180) === false
    end
end

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

@testset "Rotate270" begin
    @test typeof(@inferred(Rotate270())) <: Rotate270 <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test @inferred(Rotate270(0.7)) === Either(Rotate270(), 0.7)
        @test str_show(Rotate270()) == "Augmentor.Rotate270()"
        @test str_showconst(Rotate270()) == "Rotate270()"
        @test str_showcompact(Rotate270()) == "Rotate 270 degree"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Rotate270(), nothing)
        @test Augmentor.supports_eager(Rotate270) === true
        imgs = [
            (rect),
            (view(rect, :, :)),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for (img_in) in imgs
                res = @inferred(Augmentor.applyeager(Rotate270(), img_in))
                @test res == rotr90(rect)
                @test typeof(res) <: Array
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Rotate270) === true
        @test_throws MethodError Augmentor.applyaffine(Rotate270(), nothing)
        @test @inferred(Augmentor.toaffinemap(Rotate270(), rect)) ≈ AffineMap([6.12323e-17 1.0; -1.0 6.12323e-17], [-0.5,3.5])
        wv = @inferred Augmentor.applyaffine(Rotate270(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotr90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Rotate270) === true
        @test_throws MethodError Augmentor.applyaffineview(Rotate270(), nothing)
        wv = @inferred Augmentor.applyaffineview(Rotate270(), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == rotr90(square)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Rotate270) === true
        v = @inferred Augmentor.applylazy(Rotate270(), rect)
        @test v === view(PermutedDimsArray(rect, (2,1)), 1:1:3, 2:-1:1)
        @test v == rotr90(rect)
        @test typeof(v) <: SubArray
        v = @inferred Augmentor.applylazy(Rotate270(), view(square,1:2,1:3))
        @test v === view(PermutedDimsArray(square, (2,1)), 1:1:3, 2:-1:1)
        @test v == rotr90(rect)
        @test typeof(v) <: SubArray
        wv = @inferred Augmentor.applylazy(Rotate270(), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotr90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(Rotate270) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Rotate270) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Rotate270) === true
        v = @inferred Augmentor.applypermute(Rotate270(), rect)
        @test v === view(PermutedDimsArray(rect, (2,1)), 1:1:3, 2:-1:1)
        @test v == rotr90(rect)
        @test typeof(v) <: SubArray
        v2 = @inferred Augmentor.applypermute(Rotate270(), v)
        @test v2 === view(rect, 2:-1:1, 3:-1:1)
        @test v2 == rot180(rect)
    end
end

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

@testset "Rotate" begin
    @test typeof(@inferred(Rotate(1))) <: Rotate <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test_throws MethodError Rotate()
        @test_throws MethodError Rotate(:a)
        @test_throws MethodError Rotate([:a])
        @test_throws ArgumentError Rotate(3:1)
        @test @inferred(Rotate(0.7)) === Rotate(0.7:1:0.7)
        @test str_show(Rotate(0.7)) == "Augmentor.Rotate(0.7)"
        @test str_showconst(Rotate(0.7)) == "Rotate(0.7)"
        @test str_showcompact(Rotate(0.7)) == "Rotate 0.7 degree"
        @test @inferred(Rotate(10)) === Rotate(10:10)
        @test str_show(Rotate(10)) == "Augmentor.Rotate(10)"
        @test str_showconst(Rotate(10)) == "Rotate(10)"
        @test str_showcompact(Rotate(10)) == "Rotate 10 degree"
        op = @inferred(Rotate(-1:1))
        @test str_show(op) == "Augmentor.Rotate(-1:1)"
        @test str_showconst(op) == "Rotate(-1:1)"
        @test str_showcompact(op) == "Rotate by θ ∈ -1:1 degree"
        op = @inferred(Rotate([2,30]))
        @test op.degree == [2,30]
        @test str_show(op) == "Augmentor.Rotate([2, 30])"
        @test str_showconst(op) == "Rotate([2, 30])"
        @test str_showcompact(op) == "Rotate by θ ∈ [2, 30] degree"
    end
    @testset "randparam" begin
        @test @inferred(Augmentor.randparam(Rotate(45), rect)) === 45.0
        @test @inferred(Augmentor.randparam(Rotate(90), rect)) === 90.0
        @test @inferred(Augmentor.randparam(Rotate(1:10), rect)) in 1:10
        @test_throws MethodError Augmentor.toaffinemap(Rotate(45), rect)
        @test @inferred(Augmentor.toaffinemap(Rotate(45), rect, 45)) ≈ AffineMap([0.70710678118 -0.70710678118; 0.70710678118 0.70710678118], [1.85355339059,-0.47487373415])
        @test @inferred(Augmentor.toaffinemap(Rotate(1:90), rect, 45)) ≈ AffineMap([0.70710678118 -0.70710678118; 0.70710678118 0.70710678118], [1.85355339059,-0.47487373415])
        @test @inferred(Augmentor.toaffinemap(Rotate(90), rect, 90)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        @test @inferred(Augmentor.toaffinemap(Rotate(-90:1:90), rect, -90)) ≈ AffineMap([6.12323e-17 1.0; -1.0 6.12323e-17], [-0.5,3.5])
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Rotate(10), nothing)
        @test Augmentor.supports_eager(Rotate) === false
        img_out1 = rotl90(square)
        img_out2 = rotr90(square)
        imgs = [
            square,
            view(square, :, :),
            Augmentor.prepareaffine(square),
            OffsetArray(square, -1, -1),
            view(square, IdentityRange(1:3), IdentityRange(1:3))
        ]
        @testset "fixed parameter" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(Rotate(90), img_in))
                @test res == img_out1
                @test typeof(res) == typeof(img_out1)
                res = @inferred(Augmentor.applyeager(Rotate(-90), img_in))
                @test res == img_out2
                @test typeof(res) == typeof(img_out2)
                # test same with tuple of images
                res1, res2 = @inferred(Augmentor.applyeager(Rotate(90), (img_in, N0f8.(img_in))))
                @test res1 == img_out1
                @test res2 == img_out1
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: Array{N0f8}
                res1, res2 = @inferred(Augmentor.applyeager(Rotate(-90), (img_in, N0f8.(img_in))))
                @test res1 == img_out2
                @test res2 == img_out2
                @test typeof(res1) == typeof(img_out2)
                @test typeof(res2) <: Array{N0f8}
            end
            # check that the affine map is computed for each image
            res1, res2 = @inferred(Augmentor.applyeager(Rotate(90), (square, square2)))
            @test res1 == OffsetArray(rotl90(square), 0, 0)
            @test res2 == OffsetArray(rotl90(square2), 0, 0)
        end
        @testset "random parameter" begin
            for img_in in imgs
                res1, res2 = @inferred(Augmentor.applyeager(Rotate(1:90), (img_in, N0f8.(img_in))))
                # make sure same angle is used
                @test res1 == res2
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: Array{N0f8}
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Rotate) === true
        @test_throws MethodError Augmentor.applyaffine(Rotate(90), nothing)
        wv = @inferred Augmentor.applyaffine(Rotate(90), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv = @inferred Augmentor.applyaffine(Rotate(-90), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotr90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Rotate) === true
        @test_throws MethodError Augmentor.applyaffineview(Rotate(90), nothing)
        wv = @inferred Augmentor.applyaffineview(Rotate(90), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == rotl90(square)
        wv = @inferred Augmentor.applyaffineview(Rotate(-90), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == rotr90(square)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Rotate(90)) === true
        wv = @inferred Augmentor.applylazy(Rotate(90), square)
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv = @inferred Augmentor.applylazy(Rotate(-90), square)
        @test parent(wv).itp.coefs === square
        @test wv == rotr90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(Rotate) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Rotate) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Rotate) === false
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule coreutils_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("coreutils")
JLLWrappers.@generate_main_file("coreutils", UUID("5818bda4-868b-5068-b238-e370ed6eefef"))
end  # module coreutils_jll


text:
# parser

function read_instance(MyFileName::String)
  path = "./" * MyFileName
  # Si le fichier path existe
  if isfile(path)
    # L’ouvrir
    myFile = open(path)

    P = parse(Int64, split(readline(myFile), " ")[3])
    m = parse(Int64, split(readline(myFile), " ")[3])
    n = parse(Int64, split(readline(myFile), " ")[3])
    w1 = parse(Int64, split(readline(myFile), " ")[3])
    w2 = parse(Int64, split(readline(myFile), " ")[3])
    L = parse(Int64, split(readline(myFile), " ")[3])
    g = parse(Float64, split(readline(myFile), " ")[3])

    readline(myFile)

    t = Array{Int64,2}(ones(n,m)) # on met a zero pour l'instant
    for i in 1:n
    line =  split(readline(myFile), " ")
      for j in 1:m
      t[i,j] = parse(Int64,line[j + 1])
      end
    end

    return P, m, n, w1, w2, L, g, t
  end
end

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
function print_with_lines(out::IO, text::AbstractString)
    io = IOBuffer()
    for (i,line) in enumerate(split(text, "\n"))
        println(io, @sprintf("%-4d: %s", i, line))
    end
    write(out, take!(io))
end
print_with_lines(text::AbstractString) = print_with_lines(stdout, text)

"""
Style Type, which is used to choose different visualization/editing styles via multiple dispatch
Usage pattern:
visualize(::Style{:Default}, ...)           = do something
visualize(::Style{:MyAwesomeNewStyle}, ...) = do something different
"""
struct Style{StyleValue}
end
Style(x::Symbol) = Style{x}()
Style() = Style{:Default}()
mergedefault!(style::Style{S}, styles, customdata) where {S} = merge!(copy(styles[S]), Dict{Symbol, Any}(customdata))
macro style_str(string)
    Style{Symbol(string)}
end
export @style_str

"""
splats keys from a dict into variables
"""
macro materialize(dict_splat)
    keynames, dict = dict_splat.args
    keynames = isa(keynames, Symbol) ? [keynames] : keynames.args
    dict_instance = gensym()
    kd = [:($key = $dict_instance[$(Expr(:quote, key))]) for key in keynames]
    kdblock = Expr(:block, kd...)
    expr = quote
        $dict_instance = $dict # handle if dict is not a variable but an expression
        $kdblock
    end
    esc(expr)
end

"""
splats keys from a dict into variables and removes them
"""
macro materialize!(dict_splat)
    keynames, dict = dict_splat.args
    keynames = isa(keynames, Symbol) ? [keynames] : keynames.args
    dict_instance = gensym()
    kd = [:($key = pop!($dict_instance, $(Expr(:quote, key)))) for key in keynames]
    kdblock = Expr(:block, kd...)
    expr = quote
        $dict_instance = $dict # handle if dict is not a variable but an expression
        $kdblock
    end
    esc(expr)
end

"""
Needed to match the lazy gl_convert exceptions.
    `Target`: targeted OpenGL type
    `x`: the variable that gets matched
"""
matches_target(::Type{Target}, x::T) where {Target, T} = applicable(gl_convert, Target, x) || T <: Target  # it can be either converted to Target, or it's already the target
matches_target(::Type{Target}, x::Node{T}) where {Target, T} = applicable(gl_convert, Target, x)  || T <: Target
matches_target(::Function, x) = true
matches_target(::Function, x::Nothing) = false

signal_convert(T1, y::T2) where {T2<:Node} = lift(convert, Node(T1), y)

"""
Takes a dict and inserts defaults, if not already available.
The variables are made accessible in local scope, so things like this are possible:
gen_defaults! dict begin
    a = 55
    b = a * 2 # variables, like a, will get made visible in local scope
    c::JuliaType = X # `c` needs to be of type JuliaType. `c` will be made available with it's original type and then converted to JuliaType when inserted into `dict`
    d = x => GLType # OpenGL convert target. Get's only applied if `x` is convertible to GLType. Will only be converted when passed to RenderObject
    d = x => \"doc string\"
    d = x => (GLType, \"doc string and gl target\")
end
"""
macro gen_defaults!(dict, args)
    args.head == :block || error("second argument needs to be a block of form
    begin
        a = 55
        b = a * 2 # variables, like a, will get made visible in local scope
        c::JuliaType = X # c needs to be of type JuliaType. c will be made available with it's original type and then converted to JuliaType when inserted into data
        d = x => GLType # OpenGL convert target. Get's only applied if x is convertible to GLType. Will only be converted when passed to RenderObject
    end")
    tuple_list = args.args
    dictsym = gensym()
    return_expression = Expr(:block)
    push!(return_expression.args, :($dictsym = $dict)) # dict could also be an expression, so we need to asign it to a variable at the beginning
    push!(return_expression.args, :(gl_convert_targets = get!($dictsym, :gl_convert_targets, Dict{Symbol, Any}()))) # exceptions for glconvert.
    push!(return_expression.args, :(doc_strings = get!($dictsym, :doc_string, Dict{Symbol, Any}()))) # exceptions for glconvert.
    # @gen_defaults can be used multiple times, so we need to reuse gl_convert_targets if already in here
    for (i, elem) in enumerate(tuple_list)
        opengl_convert_target = :() # is optional, so first is an empty expression
        convert_target        = :() # is optional, so first is an empty expression
        doc_strings           = :()
        if Meta.isexpr(elem, :(=))
            key_name, value_expr = elem.args
            if isa(key_name, Expr) && key_name.head == :(::) # we need to convert to a julia type
                key_name, convert_target = key_name.args
                convert_target = :(GLAbstraction.signal_convert($convert_target, $key_name))
            else
                convert_target = :($key_name)
            end
            key_sym = Expr(:quote, key_name)
            if isa(value_expr, Expr) && value_expr.head == :call && value_expr.args[1] == :(=>)  # we might need to insert a convert target
                value_expr, target = value_expr.args[2:end]
                undecided = []
                if isa(target, Expr)
                    undecided = target.args
                else
                    push!(undecided, target)
                end
                for elem in undecided
                    isa(elem, Expr) && continue #
                    if isa(elem, AbstractString) # only docstring
                        doc_strings = :(doc_strings[$key_sym] = $elem)
                    elseif isa(elem, Symbol)
                        opengl_convert_target = quote
                            if GLAbstraction.matches_target($elem, $key_name)
                                gl_convert_targets[$key_sym] = $elem
                            end
                        end
                    end
                end
            end
            expr = quote
                $key_name = if haskey($dictsym, $key_sym)
                    $dictsym[$key_sym]
                else
                    $value_expr # in case that evaluating value_expr is expensive, we use a branch instead of get(dict, key, default)
                end
                $dictsym[$key_sym] = $convert_target
                $opengl_convert_target
                $doc_strings
            end
            push!(return_expression.args, expr)
        end
    end
    #push!(return_expression.args, :($dictsym[:gl_convert_targets] = gl_convert_targets)) #just pass the targets via the dict
    push!(return_expression.args, :($dictsym)) #return dict
    esc(return_expression)
end
export @gen_defaults!

makesignal(s::Node) = s
makesignal(v) = Node(v)

@inline const_lift(f::Union{DataType, Type, Function}, inputs...) = lift(f, map(makesignal, inputs)...)
export const_lift

isnotempty(x) = !isempty(x)
AND(a,b) = a&&b
OR(a,b) = a||b

#Meshtype holding native OpenGL data.
struct NativeMesh{MeshType <: GeometryBasics.Mesh}
    data::Dict{Symbol, Any}
end

export NativeMesh

NativeMesh(m::T) where {T <: GeometryBasics.Mesh} = NativeMesh{T}(m)
NativeMesh(m::Observable{T}) where {T <: GeometryBasics.Mesh} = NativeMesh{T}(m)

convert_texcoordinates(uv::AbstractVector{Vec2f0}) = uv
convert_texcoordinates(x::AbstractVector{<:Number}) = convert(Vector{Float32}, x)

function NativeMesh{T}(mesh::T) where T <: GeometryBasics.Mesh
    result = Dict{Symbol, Any}()
    attribs = GeometryBasics.attributes(mesh)
    result[:faces] = indexbuffer(faces(mesh))
    if isempty(attribs)
        # TODO position only shows up as attribute
        # when it has meta informtion
        result[:vertices] = GLBuffer(collect(metafree(coordinates(mesh))))
    end
    for (field, val) in attribs
        if val isa Makie.Sampler
            result[:image] = Texture(val.colors)
            result[:texturecoordinates] = GLBuffer(convert_texcoordinates(val.values))
            if val.scaling.range !== nothing
                result[:color_norm] = Observable(Vec2f0(val.scaling.range))
            end
        elseif field in (:position, :uv, :uvw, :normals, :attribute_id, :color)
            if field == :color
                field = :vertex_color
            elseif field in (:uv, :uvw)
                field = :texturecoordinates
            elseif field == :position
                field = :vertices
            end
            if val isa AbstractVector
                result[field] = GLBuffer(metafree(val))
            end
        else
            result[field] = Texture(val)
        end
    end
    return NativeMesh{T}(result)
end

function NativeMesh{T}(m::Node{T}) where T <: GeometryBasics.Mesh
    result = NativeMesh{T}(m[])
    on(m) do mesh

        attribs = GeometryBasics.attributes(mesh)

        update!(result.data[:faces], faces(mesh))

        if isempty(attribs)
            # TODO position only shows up as attribute
            # when it has meta informtion
            update!(result.data[:vertices], metafree(coordinates(mesh)))
        end

        for (field, val) in attribs
            if val isa Makie.Sampler
                update!(result.data[:image], val.colors)
                update!(result.data[:texturecoordinates], convert_texcoordinates(val.values))
                if val.scaling.range !== nothing
                    result.data[:color_norm][] = Vec2f0(val.scaling.range)
                end
            else
                field == :color && (field = :vertex_color)
                field == :uv && (field = :texturecoordinates)
                field == :position && (field = :vertices)
                haskey(result.data, field) && update!(result.data[field], val)
            end
        end
    end
    return result
end


text:
using DataFrames
using DataConvenience

df = DataFrame(
        a = rand(1:8, 1000),
        b = rand(1:8, 1000),
        c = rand(1:8, 1000),
    )

nest(df, :a, :meh)
unnest(nest(df, :a, :meh), :meh)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
const Mat = Matrix{Float64}
const Vec = Vector{Float64}
const Vec3 = Vector{Float64}
const Vec2 = Vector{Float64}
const SparseMatrix = Array

struct Camera
    name :: string
    focal_length :: Real
    position :: Vec3
    orientation :: Vec3
end

struct SLAM
    cameras :: Array{Camera}
    points :: Array{Vec3}
end

function loss(params :: SLAM, data :: SparseMatrix{Vec2})

end

struct dCamera
    name :: ()
    focal_length :: Real
    position :: Vec3
    orientation :: Vec3
end

struct dSLAM
    cameras :: Array{dCamera}
    points :: Array{Vec3}
end

function grad_loss(params :: SLAM, data :: SparseMatrix{Vec2}) :: Pair{dSLAM, SparseMatrix{Vec2}}

max_stars_repo_name: tomjaguarpaw/knossos-ksc

max_issues_repo_name: tomjaguarpaw/knossos-ksc

max_forks_repo_name: tomjaguarpaw/knossos-ksc


text:
#=
Copyright 2020 Ali Mert Ceylan

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

import Base: +, -, *, /, ^, convert, promote_rule, length, size, ndims, getindex, setindex!, similar, iterate, broadcasted
using LinearAlgebra

"""
    unbroadcast(x,dx)

Bring dx to x's size via unbroadcasting (reduction). This is needed
when defining gradients of multi-argument broadcasting functions where
the arguments and the result may be of different sizes.

Courtesy of AutoGrad.jl

Copyright (c) 2016: Deniz Yuret.

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""
function unbroadcast(x, dx)
    if size(x) == size(dx)
        return dx
    elseif isa(x, Number)
        return sum(dx)
    elseif isa(dx, Number)
        return fill!(similar(x), dx)
    else
        d = []
        for i = 1:ndims(dx)
            size(x, i) == size(dx, i) > 1 && continue
            size(x, i) != 1 && throw(DimensionMismatch())
            push!(d, i)
        end
        length(d) == 1 && (d = d[1])
        return reshape(sum(dx, dims = d), size(x))
    end
end

# This fixes unbroadcast when x isa Broadcasted
Base.size(x::Base.Broadcast.Broadcasted, i...) = length.(axes(x, i...))

mutable struct BD{V}
    f::Tuple{V,Function}
end

function BD{V}(v1::V) where V <: Union{Number,AbstractArray}
    return BD{V}((v1, dy->dy))
end

function BD{V}(v1::V, s::Symbol) where V <: Union{Number,AbstractArray}
    return BD{V}((v1, dy->Pair(s, dy)))
end

function BD(v1::V, s::Symbol) where V <: Union{Number,AbstractArray}
    return BD{typeof(v1)}((v1, dy->Pair(s, dy)))
end

function BD{V}(v1::V, v2) where V <: Union{Number,AbstractArray}
    return BD{V}((v1, v2))
end

function BD(v1::V, v2) where V <: Union{Number,AbstractArray}
    return BD{typeof(v1)}((v1, v2))
end

value(x::BD) = x.f[1]
func(x::BD) = x.f[2]

convert(::Type{BD}, x::U) where U <: Number = BD{U}(convert(U, x))
convert(::Type{BD{T}}, x::U) where {T <: Number,U <: Number} = BD{T}(convert(T, x))
convert(::Type{BD{T}}, x::BD{U}) where {T, U} = BD{T}(convert(T, x.f[1]), x.f[2])

convert(::Type{BD{T}}, x::T) where {T <: AbstractArray} = BD{T}(x)
convert(::Type{BD{T}}, x::BD{U}) where {T <: AbstractArray, U <: Hermitian} = BD{T}(convert(T, x.f[1]), x.f[2])

promote_rule(::Type{BD}, x::T) where {T} = BD{T}(x)
promote_rule(x::T, ::Type{BD}) where {T} = BD{T}(x)

promote_rule(::Type{BD{T}}, ::Type{U}) where {T <: Number,U <: Number} = BD{promote_type(T, U)}

promote_rule(::Type{T}, ::Type{U}) where {T <: Hermitian, U <: AbstractArray} = U
promote_rule(::Type{U}, ::Type{T}) where {T <: Hermitian, U <: AbstractArray} = U
# promote_rule(::Type{BD{T}}, ::Type{U}) where {T <: Union{Number, AbstractArray}, U <: Union{Number, AbstractArray}} = BD{promote_type(T, U)}
promote_rule(::Type{BD{T}}, ::Type{T}) where T <: AbstractArray = BD{T}

promote_rule(::Type{BD{T}}, ::Type{BD{U}}) where {T, U} = BD{promote_type(T, U)}
# promote_rule(::Type{BD{T}}, ::Type{BD{U}}) where {T<:Number, U<:Number} = BD{promote_type(T, U)}

ndims(bd::BD) = ndims(bd.f[1])
length(bd::BD) = length(bd.f[1])
size(bd::BD) = size(bd.f[1])

getindex(bd::BD{T}, i::Int) where {T <: Number} = BD{eltype(T)}(getindex(bd.f[1], i))
getindex(bd::BD{T}, i::Int) where {T <: AbstractArray} = BD{eltype(T)}(getindex(bd.f[1], i), (dy)->(dyi = zeros(size(bd)); dyi[i] = dy; bd.f[2](dyi)))
getindex(bd::BD{<:AbstractArray}, i) = BD{Array}(getindex(bd.f[1], i), (dy)->(dyi = zeros(size(bd)); dyi[i] .= dy; bd.f[2](dyi)))#

setindex!(bd::BD, v, i) = setindex!(bd.f[1], v, i)
setindex!(bd::BD, v, i::Vararg{Bool}) = setindex!(bd.f[1], i...)

length(bd::BD{T}) where T <: AbstractArray = length(bd.f[1])
similar(bd::BD{T}) where T <: Number = BD{T}(zero(T))
similar(bd::BD{T}) where T <: AbstractArray = BD{T}(zeros(size(bd)))

iterate(bd::BD{<:Union{Number, AbstractArray}}, state = 1) = state <= length(bd) ? (bd[state], state + 1) : nothing

for diffrule in DiffRules.diffrules()
    p, f, a = diffrule

    if (a == 1)

        eval(
            :($(p).$(f)(x::BD) = BD(
                $(p).$(f)(value(x)),
                (dy)->func(x)(dy.*eval(DiffRules.diffrule(Symbol($(p)), Symbol($(f)), value(x))))
                    )
                )
            )

    elseif (a == 2)

        eval(
            :($(p).$(f)(x::BD, y::BD) = BD(
                $(p).$(f)(value(x), value(y)),
                (dy)->(
                    df = DiffRules.diffrule(Symbol($(p)), Symbol($(f)), value(x), value(y));
                    (func(x)(dy*eval(df[1])'), func(y)(eval(df[2])'*dy))
                    )
                ))
            )

        eval(
            :($(f)(x::BD{T}, y::U) where {T <: Union{Number,AbstractArray}, U <: Union{Number,AbstractArray}} = $(f)(promote(x, y)...))
            )

        eval(
            :($(f)(x::T, y::BD{U}) where {T <: Union{Number,AbstractArray}, U <: Union{Number,AbstractArray}} = $(f)(promote(x, y)...))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::BD{T}, y::BD{U}) where {T <: Number,U <: AbstractArray} = BD{U}(
                o.(value(x), value(y)),
                (dy)->(func(x)(unbroadcast(value(x), dy)),  func(y)(unbroadcast(value(y), dy)))
                ))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::BD{U}, y::BD{T}) where {T <: Number,U <: AbstractArray} = BD{U}(
                o.(value(x), value(y)),
                (dy)->(func(x)(unbroadcast(x.f[1], dy)),  func(y)(unbroadcast(value(y), dy)))
                ))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::T, y::BD{U}) where {T <: Number, U <: AbstractArray} = broadcasted(o, BD{T}(x), y))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::BD{T}, y::U) where {T <: Number, U <: AbstractArray} = broadcasted(o, x, BD{U}(y)))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::T, y::BD{U}) where {T <: AbstractArray, U <: Number} = broadcasted(o, BD{T}(x), y))
            )

        eval(
            :(broadcasted(o::typeof($(p).$(f)), x::BD{T}, y::U) where {T <: AbstractArray, U <: Number} = broadcasted(o, x, BD{U}(y)))
            )

    end
end


text:
struct ProtectedException <: Exception
    msg::String
end

abstract type Permission end
struct Read <: Permission end
struct Write <: Permission end

struct Protected{T}
    data::T
    perms::NamedTuple
end

function Base.getproperty(protected::Protected, field::Symbol)
    return getfieldperm(
        getfield(protected, :data),
        field,
        get(getfield(protected, :perms), field, nothing),
    )
end

function Base.setproperty!(protected::Protected, field::Symbol, x)
    return setfieldperm!(
        getfield(protected, :data),
        field,
        x,
        get(getfield(protected, :perms), field, nothing),
    )
end

function Base.getindex(protected::Protected, idx)
    return getindexperm(
        getfield(protected, :data),
        idx,
        get(getfield(protected, :perms), idx, nothing),
    )
end

function Base.setindex!(protected::Protected, x, idx)
    return setindexperm!(
        getfield(protected, :data),
        x,
        idx,
        get(getfield(protected, :perms), idx, nothing),
    )
end

Base.haskey(d::Protected, key) = haskey(getfield(d, :data), key)

getindexperm(data::D, idx, perm::Nothing) where D =
    throw(ProtectedException("Read access to $D[$idx] disallowed."))
getindexperm(data, idx, perm::Union{Read, Write}) = getindex(data, idx)
getindexperm(data, idx, perm::NamedTuple) = protect(getindex(data, idx), perm)

setindexperm!(data::D, idx, x, perm::Union{Read, Nothing, NamedTuple}) where D =
    throw(ProtectedException("Write access to $D[$idx] disallowed."))
setindexperm!(data, x, idx, perm::Write) = setindex!(data, x, idx)

getfieldperm(data::D, field, perm::Nothing) where D =
    throw(ProtectedException("Read access to $(D).$(string(field)) disallowed."))
getfieldperm(data, field, perm::Union{Read, Write}) = getproperty(data, field)
getfieldperm(data, field, perm::NamedTuple) = protect(getproperty(data, field), perm)

setfieldperm!(data::D, field, x, perm::Union{Read, Nothing, NamedTuple}) where D =
    throw(ProtectedException("Write access to $(D).$(string(field)) disallowed."))
setfieldperm!(data, field, x, perm::Write) = setproperty!(data, field, x)

Base.fieldnames(protected::Protected) = fieldnames(typeof(getfield(protected, :data)))

protect(x, perms) = Protected(x, perms)

"""
    PropDict(dict)

Like a `Dict{Symbol}`, but attribute syntax can be used to access values.
"""
struct PropDict{V}
    d::Dict{Symbol, V}
end

Base.getproperty(d::PropDict, field::Symbol) = getfield(d, :d)[field]
Base.setproperty!(d::PropDict, field::Symbol, val) = (getfield(d, :d)[field] = val;)
Base.propertynames(d::PropDict) = Tuple(keys(getfield(d, :d)))

Base.haskey(d::PropDict, key) = haskey(getfield(d, :d), key)
Base.getindex(d::FluxTraining.PropDict{Any}, i::Symbol) = getproperty(d, i)
Base.keys(d::PropDict) = keys(getfield(d, :d))
Base.values(d::PropDict) = values(getfield(d, :d))
Base.get(d::PropDict, args...; kwargs...) = get(getfield(d, :d), args...; kwargs...)

PropDict(args...) = PropDict(Dict{Symbol, Any}(args...))

# ## Tests

@testset "`protect`" begin
    mutable struct C
        x::Int
        y::Int
    end

    mutable struct B
        c::C
    end

    mutable struct A
        b1::B
        b2::B
    end

    makea() = A(
        B(C(1, 2)),
        B(C(3, 4)),
    )

    @testset "fully protected" begin
        a_p = protect(makea(), NamedTuple())  # protect all child structs of `a`
        @test_throws ProtectedException a_p.b1 isa Protected
    end

    @testset "read direct child" begin
        a_p2 = protect(makea(), (b1=Read(),))  # protect everything but `b1`
        @test_nowarn a_p2.b1
        @test_throws ProtectedException a_p2.b1 = C(5, 5)
    end

    @testset "write direct child" begin
        a_p2 = protect(makea(), (b1=Write(),))  # protect everything but `b1`
        @test_nowarn a_p2.b1
        @test_nowarn a_p2.b1 = B(C(5, 5))
    end
    @testset "write nested" begin
        a_p3 = protect(makea(), (b1 = (c = (x=Write(),),),))  # allow mutating only a.b1.c.x
        @test_nowarn a_p3.b1.c.x = 2
        @test_nowarn a_p3.b1.c.x
        @test a_p3.b1.c isa Protected
        @test_throws ProtectedException a_p3.b1.c.y = 2
    end
end


text:
"""
Each entropy word in an entropic vector is sorted.

#### Examples
```julia-repl
julia> order_entropic1("h24-h32-h132-h2")
"h24 - h23 - h123 - h2"

julia> order_entropic1("7h32 - h243 - h13701 - h92252")
"7h23-h234-h01137-h22259"

```
"""
function order_entropic1(s::AbstractString="Hello world")

	s0=replace(s,"h"=>"h ","+"=>" + ","-"=>" - ")
	x=split(s0," ") # Split into words delim=" "

	Y=""
	for wIdx = 1:length(x)
		# w=x[wIdx] # word
		y=string(sort([a for a in x[wIdx]])...)  # sort each word
		Y= Y * " " * y
	end
	Yout=replace(Y,"h "=>"h"," " =>"")

	return Yout
end

"""
Each word in a sentence (string) is sorted alphabetically.

#### Examples
```julia-repl
julia> order_string("This is a sorted sentence; Who is 1 and two  ")
" This is a deorst ;ceeennst Who is 1 adn otw  "

```
"""
function order_string(s::AbstractString="Hello world")

	s0=replace(s,"+"=>" + ","-"=>" - ")
	x=split(s0," ") # Split into words delim=" "

	Y=""
	for wIdx = 1:length(x)
		# w=x[wIdx] # word
		y=string(sort([a for a in x[wIdx]])...)  # sort each word
		Y= Y * " " * y
	end
	Yout=replace(Y," " =>" ")

	return Yout
end

"""
	find_subset(n::Int64,p,q,RV::AbstractString="X")

Given `i` and `j` compute `Κ ⊆ \U1d4a9 \\{i,j}`; i.e., all non exclusive subsets. `i,j` can also be empty (i.e. []), in which case the non-empty superset gets listed. An optional prefix can be added (default is `X`).
#### Examples
```julia-repl
julia>find_subset(4,1,3,"X")
["" "X1" "X1,X2" "X1,X2,X5" "X1,X5" "X2" "X2,X5" "X5"]

julia>find_subset(4,1,3,"")
["" "1" "1,2" "1,2,5" "1,5" "2" "2,5" "5"]

julia> find_subset(5,4,3,"🍒")
["" "🍒3" "🍒3,🍒4" "🍒4"]

julia> find_subset(5,[],[],"🍓")
["🍓1" "🍓1,🍓2" "🍓1,🍓2,🍓3" "🍓1,🍓3" "🍓2" "🍓2,🍓3" "🍓3"]
```
"""
function find_subset(n::Int64,p,q,RV::AbstractString="X")
	s="$(RV)" .* string.(1:n)
	s1="$(RV)" .* string(p)
	s2="$(RV)" .* string(q)
	PowerSet = collect(powerset(s,1,n)) # Skip empty set
	setXk=[]
	for kk = 1:Int64(exp2(n)-1)
		subSet1 = vcat(PowerSet[kk,:][1])
		subSet2 = setdiff(subSet1,[s1,s2])
		setXk=[setXk; join(subSet2,",")]
	end

	return sort(unique(setXk))
end


text:
using JSON, Requests
reddit_url = "https://www.reddit.com/r/Julia/"
response = get("$(reddit_url)/.json")
dataReceived = JSON.parse(Requests.text(response))
nextRecord = dataReceived["data"]["after"]
counter = length(dataReceived["data"]["children"])

allPosts = []
for record in 1:counter
    url = dataReceived["data"]["children"][record]["data"]["url"]
    redditrecord_id  = dataReceived["data"]["children"][record]["data"]["id"]
    redditrecord_title  = dataReceived["data"]["children"][record]["data"]["title"]
    author  = dataReceived["data"]["children"][record]["data"]["author"]
    created = dataReceived["data"]["children"][record]["data"]["created"]
    push!(allPosts, (url, redditrecord_id, redditrecord_title, author, created))
end

for post in allPosts
    println(post[3])
end

max_stars_repo_path: Chapter 2/10.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 2/10.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 2/10.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
function perturb!(bool::BoolParameter)
    bool.value = !(bool.value)
    bool
end

function perturb!(number::NumberParameter)
    typeof(number.value) <: Integer ?
    perturb!(number, 1) : perturb!(number, 1.)
end

function perturb!(number::NumberParameter, interval::Number)
    if interval <= 0
        error("interval must be greater than zero.")
    end
    max = number.value + interval > number.max ? number.max : number.value + interval
    min = number.value - interval < number.min ? number.min : number.value - interval
    number.value = rand_in(min, max)
    number
end

function perturb!(enum::EnumParameter)
    enum.value = rand(1:length(enum.values))
    @inbounds enum.current = enum.values[enum.value]
    enum
end

function perturb!(permutation::PermutationParameter)
    shuffle!(permutation.value)
    permutation
end

function perturb!(permutation::PermutationParameter, interval::Int)
    for i = 1:interval
        a = rand(1:permutation.size)
        b = rand(1:permutation.size)
        while b == a
            b = rand(1:permutation.size)
        end
        permutation.value[a], permutation.value[b] = permutation.value[b],
                                                     permutation.value[a]
    end
    permutation
end

function perturb!(configuration::Configuration)
    for key in keys(configuration.parameters)
        if !(typeof(configuration[key]) <: StringParameter)
            perturb!(configuration[key])
        end
    end
    update!(configuration)
    configuration
end

function perturb!(configuration::Configuration, intervals::Dict{String, Any})
    for key in keys(intervals)
        @inbounds perturb!(configuration[key], intervals[key])
    end
    update!(configuration)
    configuration
end


text:
struct StochasticSearchOptimizer{P,S} <: GlobalOptimizer
    algo::Symbol
    population::P
    srand::S
    options::Dict{Symbol,Any}
end

need_gradient(::StochasticSearchOptimizer) = false

function CRSOptimizer(; population = 0, srand = nothing, options...)
    StochasticSearchOptimizer(
        :GN_CRS2_LM,
        population,
        srand,
        Dict{Symbol,Any}(options),
    )
end

function ISRESOptimizer(; population = 0, srand = nothing, options...)
    StochasticSearchOptimizer(
        :GN_ISRES,
        population,
        srand,
        Dict{Symbol,Any}(options),
    )
end

function calc_population(o::StochasticSearchOptimizer, lb, ub)
    if isa(o.population, Real)
        convert(Int, o.population)
    else
        convert(Int, o.population(lb, ub))
    end
end

function minimize(o::StochasticSearchOptimizer, f, lb, ub, x0)
    opt = create_nlopt(o.algo, lb, ub; min_objective = f, o.options...)
    pop = calc_population(o, lb, ub)
    !isnothing(pop) && NLopt.population!(opt, pop)
    !isnothing(o.srand) && NLopt.srand(o.srand)
    fx, x, ret = NLopt.optimize(opt, x0)
    NLopt.srand_time()
    ret == :FORCED_STOP && throw(InterruptException())
    return fx, x, ret
end

function maximize(o::StochasticSearchOptimizer, f, lb, ub, x0)
    opt = create_nlopt(o.algo, lb, ub; max_objective = f, o.options...)
    pop = calc_population(o, lb, ub)
    !isnothing(pop) && NLopt.population!(opt, pop)
    !isnothing(o.srand) && NLopt.srand(o.srand)
    fx, x, ret = NLopt.optimize(opt, x0)
    NLopt.srand_time()
    ret == :FORCED_STOP && throw(InterruptException())
    return fx, x, ret
end