text:
## plot.jl : functions for easy visualization of Gaussian random fields

## type aliases ##
const OneDimCov = CovarianceFunction{1}
const TwoDimCov = CovarianceFunction{2}
const ThreeDimCov = CovarianceFunction{3}
const OneDimSepCov = SeparableCovarianceFunction{1}
const TwoDimSepCov = SeparableCovarianceFunction{2}
const ThreeDimSepCov = SeparableCovarianceFunction{3}

const OneDimGRF = GaussianRandomField{C} where {C<:Union{OneDimCov,OneDimSepCov}}
const TwoDimGRF = GaussianRandomField{C} where {C<:Union{TwoDimCov,TwoDimSepCov}}
const ThreeDimGRF = GaussianRandomField{C} where {C<:Union{ThreeDimCov,ThreeDimSepCov}}
const FiniteElemGRF = GaussianRandomField{C,M,Tuple{T1,T2}} where {C<:TwoDimCov,M,T1<:AbstractMatrix,T2<:AbstractMatrix}

const OneDimSpectralGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:OneDimCov,n}
const TwoDimSpectralGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:TwoDimCov,n}
const ThreeDimSpectralGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:ThreeDimCov,n}
const OneDimSpectralSepGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:OneDimSepCov,n}
const TwoDimSpectralSepGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:TwoDimSepCov,n}
const ThreeDimSpectralSepGRF = GaussianRandomField{C,KarhunenLoeve{n}} where {C<:ThreeDimSepCov,n}
const FiniteElemSpectralGRF = GaussianRandomField{C,KarhunenLoeve{n},Tuple{T1,T2}} where {C<:TwoDimCov,n,T1<:AbstractMatrix,T2<:AbstractMatrix}

## 1D ##
function plot(grf::OneDimGRF;n=1,kwargs...)
    for i in 1:n
        plot(grf.pts[1],sample(grf);kwargs...)
    end
end

## 2D ##
plot(grf::TwoDimGRF;kwargs...) = surf(grf;kwargs...)
plot(grf::FiniteElemGRF;kwargs...) = tricontourf(grf;kwargs...)

function surf(grf::TwoDimGRF;kwargs...)
    x,y = grf.pts
    xgrid = repmat(x,1,length(y))
    ygrid = repmat(y',length(x),1)
    plot_surface(ygrid,xgrid,reshape(sample(grf),(length(x),length(y))),rstride=2,edgecolors="k",cstride=2,cmap=ColorMap("viridis");kwargs...)
end

function contourf(grf::TwoDimGRF;kwargs...)
    x,y = grf.pts
    xgrid = repmat(x,1,length(y))
    ygrid = repmat(y',length(x),1)
    contourf(xgrid,ygrid,reshape(sample(grf),(length(x),length(y)));kwargs...)
    colorbar()
end

function contour(grf::TwoDimGRF;kwargs...)
    x,y = grf.pts
    xgrid = repmat(x,1,length(y))
    ygrid = repmat(y',length(x),1)
    cp = contour(xgrid,ygrid,reshape(sample(grf),(length(x),length(y))),colors="black";kwargs...)
    clabel(cp, inline=1, fontsize=10)
end

function tricontourf(grf::FiniteElemGRF;kwargs...)
    (p,t) = grf.pts
    isempty(t) && throw(ArgumentError("cannot plot mesh when random field is computed in element centers"))
    x = p[1,:]
    y = p[2,:]
    tricontourf(x,y,sample(grf),triangles=t'-1,cmap=get_cmap("viridis");kwargs...)
    triplot(x,y,triangles=t'-1,color="k",linewidth=0.5;kwargs...)
end

function plot_trisurf(grf::FiniteElemGRF;kwargs...)
    (p,t) = grf.pts
    isempty(t) && throw(ArgumentError("cannot plot mesh when random field is computed in element centers"))
    x = p[1,:]
    y = p[2,:]
    plot_trisurf(x,y,sample(grf),triangles=t'-1,cmap=get_cmap("viridis");kwargs...)
end

## 3D ##
function plot(grf::ThreeDimGRF;kwargs...)
    kwargs = Dict(kwargs)
    x,y,z = grf.pts
    nx,ny,nz = length.(grf.pts)
    nx2 = haskey(kwargs,:ix) ? is_valid_idx(kwargs[:ix],nx) : round(Int,nx/2)
	delete!(kwargs,:ix)
    ny2 = haskey(kwargs,:iy) ? is_valid_idx(kwargs[:iy],ny) : round(Int,ny/2)
	delete!(kwargs,:iy)
    nz2 = haskey(kwargs,:iz) ? is_valid_idx(kwargs[:iz],nz) : round(Int,nz/2)
	delete!(kwargs,:iz)
    A = reshape(sample(grf),(nx,ny,nz))
    slice(x,y,z,nx2,ny2,nz2,A,:x;kwargs...)
    slice(x,y,z,nx2,ny2,nz2,A,:y;kwargs...)
    slice(x,y,z,nx2,ny2,nz2,A,:z;kwargs...)
end

is_valid_idx(idx,n) = ( typeof(idx)<:Integer && idx > 0 && idx < n+1 ) ? idx : throw(ArgumentError("invalid index $(idx)"))

function slice(x,y,z,dx,dy,dz,A,mode;kwargs...)
    for quadrant in (:I,:II,:III,:IV)
        quadrant_slice((x,y,z),(dx,dy,dz),A,quadrant,mode;kwargs...)
    end
end

function quadrant_slice(xs,dxs,A,quadrant,mode;kwargs...)
    if mode == :x
        order = (1,2,3)
    elseif mode == :y
        order = (2,1,3)
    elseif mode == :z
        order = (3,2,1)
    end

    if quadrant == :I
        idcs1 = 1:dxs[order[2]]
        idcs2 = 1:dxs[order[3]]
    elseif quadrant == :II
        idcs1 = dxs[order[2]]:length(xs[order[2]])
        idcs2 = 1:dxs[order[3]]
    elseif quadrant == :III
        idcs1 = dxs[order[2]]:length(xs[order[2]])
        idcs2 = dxs[order[3]]:length(xs[order[3]])
    elseif quadrant == :IV
        idcs1 = 1:dxs[order[2]]
        idcs2 = dxs[order[3]]:length(xs[order[3]])
    end
    mslices = (dxs[order[1]],idcs1,idcs2)
    cut = A[[mslices[i] for i in order]...]
	cut = mode == :z ? cut' : cut
    xgrid = xs[order[1]][dxs[order[1]]]*ones(length(idcs1),length(idcs2))
    ygrid = [xs[order[2]][i] for i = idcs1, j = idcs2]
    zgrid = [xs[order[3]][j] for i = idcs1, j = idcs2]
    grids = (xgrid,ygrid,zgrid)
    plot_surface([grids[i] for i in order]...,facecolors=get_cmap("viridis").o((cut-minimum(A))/(maximum(A)-minimum(A))),shade=false;kwargs...)
end

## eigenvalues and eigenfunctions ##
function plot_eigenvalues(grf::GaussianRandomField{C,KarhunenLoeve{n}} where {C<:CovarianceFunction,n})
    ev = grf.data.eigenval.^2
    loglog(1:length(ev),ev)
    xlabel("n")
    ylabel("magnitude")
end

function plot_eigenvalues(grf::GaussianRandomField{C,CirculantEmbedding} where {C<:CovarianceFunction,n})
    ev = sort(grf.data[1][:],rev=true).^2
    loglog(1:length(ev),ev)
    xlabel("n")
    ylabel("magnitude")
end

function plot_eigenvalues(grf::GaussianRandomField{S,KarhunenLoeve{n}} where {S<:SeparableCovarianceFunction}) where {n}
	(order,data) = grf.data
	ev = zeros(n)
	for i in eachindex(ev)
		ev[i] = prod([data[j].eigenval[order[i][j]] for j = 1:length(grf.cov.cov)]).^2
	end
    loglog(1:length(ev),ev)
    xlabel("n")
    ylabel("magnitude")
end

function plot_eigenfunction(grf::GaussianRandomField{C,KarhunenLoeve{n}} where {C,n}, n::Integer; kwargs...)
    ( n > 0 && n < randdim(grf)+1 ) || throw(ArgumentError("eigenfunction index n must be between 1 and $(randdim(grf))"))
    _plot_eigenfunction(grf,n;kwargs...)
end

function _plot_eigenfunction(grf::OneDimSpectralGRF, n::Integer; kwargs...)
    PyPlot.plot(grf.pts[1],grf.data.eigenfunc[:,n])
end

function _plot_eigenfunction(grf::TwoDimSpectralGRF, n::Integer; kwargs...)
    x,y = grf.pts
    xgrid = repmat(x,1,length(y))
    ygrid = repmat(y',length(x),1)
    contourf(xgrid,ygrid,reshape(grf.data.eigenfunc[:,n],(length(x),length(y)));kwargs...)
    colorbar()
end

function _plot_eigenfunction(grf::OneDimSpectralSepGRF, n::Integer; kwargs...)
	(order,data) = grf.data
	plot(grf.pts[1],data[1].eigenfunc[:,n])
end

function _plot_eigenfunction(grf::TwoDimSpectralSepGRF, n::Integer; kwargs...)
    x,y = grf.pts
    xgrid = repmat(x,1,length(y))
    ygrid = repmat(y',length(x),1)
	(order,data) = grf.data
	ef = kron([data[j].eigenfunc[:,order[n][j]] for j = 1:length(grf.cov.cov)]...)
    contourf(xgrid,ygrid,reshape(ef,(length(x),length(y)));kwargs...)
    colorbar()
end

function _plot_eigenfunction(grf::ThreeDimSpectralGRF, n::Integer; kwargs...)
    x,y,z = grf.pts
    nx = length(x); nx2 =round(Int,nx/2)
    ny = length(y); ny2 =round(Int,ny/2)
    nz = length(z); nz2 =round(Int,nz/2)
    A = reshape(grf.data.eigenfunc[:,n],(nx,ny,nz))
    slice(x,y,z,nx2,ny2,nz2,A,:x)
    slice(x,y,z,nx2,ny2,nz2,A,:y)
    slice(x,y,z,nx2,ny2,nz2,A,:z)
end

function _plot_eigenfunction(grf::ThreeDimSpectralSepGRF, n::Integer; kwargs...)
    x,y,z = grf.pts
    nx = length(x); nx2 =round(Int,nx/2)
    ny = length(y); ny2 =round(Int,ny/2)
    nz = length(z); nz2 =round(Int,nz/2)
	(order,data) = grf.data
	ef = kron([data[j].eigenfunc[:,order[n][j]] for j = 1:length(grf.cov.cov)]...)
    A = reshape(ef,(nx,ny,nz))
    slice(x,y,z,nx2,ny2,nz2,A,:x)
    slice(x,y,z,nx2,ny2,nz2,A,:y)
    slice(x,y,z,nx2,ny2,nz2,A,:z)
end

function _plot_eigenfunction(grf::FiniteElemSpectralGRF, n::Integer; kwargs...)
    (p,t) = grf.pts
    isempty(t) && throw(ArgumentError("cannot plot mesh when random field is computed in element centers"))
    x = p[1,:]
    y = p[2,:]
    tricontourf(x,y,grf.data.eigenfunc[:,n],triangles=t'-1,cmap=get_cmap("viridis");kwargs...)
    triplot(x,y,triangles=t'-1,color="k",linewidth=0.5;kwargs...)
end

## plot_covariance_matrix ##
function plot_covariance_matrix(c::CovarianceFunction,pts::V...) where {V<:Range}
    length(pts) == ndims(c) || throw(DimensionMismatch("dimension of points must agree with dimension of covariancefunction"))
    C = apply(c,pts,pts)
    imshow(C)
    colorbar()
end


text:
using LinearAlgebra

"""
RESETQR using JULIA functions. Used to test the algorithm at the beginning of the project.
"""
function resetqr_julia(A,k)
    @time begin

    min_error = 1e-16

    U = 0
    V = 0

    U,V = get_U_and_V_julia(A, k, min_error)

    end
    return U,V
end

function get_U_and_V_julia(A, k, min_error, max_steps = 20)
    (m,n) = size(A)
    U = rand(m,k)
    V = rand(k,n)
    error = min_error + 1
    num_step = 0

    f = norm(A - U*V,2)

    while(num_step < max_steps)
        #gradient = [U*V*V' - A*V', U'*U*V - U'*A]

        f_old = norm(A - U*V,2)
        Q,R = qr(V')
        U = (inv(R)*Q'*A')'
        error = norm(A - U*V,2)
        Q,R = qr(U)
        V = inv(R)*Q'*A
        f = norm(A - U*V,2)

        println(norm(A - U*V,2))

        if norm(f_old - f,2) <= min_error
            break
        end

        num_step += 1
    end

    return U,V
end

ext: jl


text:
# # A fast [Trie](https://github.com/gkappler/Tries.jl)-based parser for a collection of literal Strings.
# Matching any one of many `String`s (e.g. names) is slow in PCRE with `|` alternatives.
using Random
s = [ randstring(10) for _ in 1:1000 ];
re = Regex(join(s,"|"));

using BenchmarkTools
@benchmark match(re,s[end])

# A optimized state-machine parser will perform much better.
#
# TODO: Benchmark Automa.jl here
#
# `CombinedParsers` does not use state machines, but
# provides fast `_iterate` implementation for `Either{Trie{Char}}`,
# that is constructed with `Either(::Vector{<:AbstractString})`
using CombinedParsers
pc = Either(s);
@benchmark match(pc,s[end])

# Note that `pc` does not capture the result by default with `parse`.
parse(!pc,s[end])

# For large word collections, PCRE fails
s = [ randstring(10) for _ in 1:10000 ];
re = Regex(join(s,"|"));

# `CombinedParsers` displays good time complexity.
pc = Either(s);
@benchmark Either(s);
@benchmark match(pc,s[end])

# `Either{Trie{Char}}` can be used to perform fast text search in Julia.
# ### Next steps
# - Flesh out text search example.
# - Parallelization
# - Relate to Aho-Corasick algorithm.


text:
head = Node(1)
n = append(head, 2)
n = append(n, 3)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
__precompile__()

module EarthMoversDistance

export emd

using Compat
import Compat.String

const libemd = joinpath(dirname(@__FILE__), "..", "deps", "libemd")

function __init__()
    # If any init or destroy functions need to be called for the lib.
    # In this case not, since emd is a "pure" function (no side effects).
    #ccall((:emd_init,libemd), Void, ())
    #atexit() do
    #    ccall((:emd_destroy,libemd), Void, ())
    #end
end

# low-level functions:
# double emd(int n_x, double *weight_x,
#           int n_y, double *weight_y,
#           double *cost, double *flows);
libemd_emd(n_x, weight_x, n_y, weight_y, cost, flows) =
    ccall((:emd,libemd), Cdouble,
          (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          n_x, weight_x, n_y, weight_y, cost, flows)

# API
function emd!(Xweight::Vector{Float64}, Yweight::Vector{Float64}, cost::Matrix{Float64},
    flows::Union{Void, Matrix{Float64}} = nothing)

    n_x = length(Xweight)
    n_y = length(Yweight)

    @assert n_x > 0
    @assert n_y > 0

    @assert n_x == size(cost, 1)
    @assert n_y == size(cost, 2)

    if flows != nothing
        @assert n_y == size(flows, 2)
        @assert n_x == size(flows, 1)
        libemd_emd(n_x, Xweight, n_y, Yweight, cost, flows)
    else
        libemd_emd(n_x, Xweight, n_y, Yweight, cost, C_NULL)
    end
end

function emd(Xweight::Vector{Float64}, Yweight::Vector{Float64}, cost::Union{Void, Matrix{Float64}} = nothing)
    if cost == nothing
        cost = ones(Float64, length(Xweight), length(Yweight))
    end
    emd!(Xweight, Yweight, cost)
end

function emd_flows(Xweight::Vector{Float64}, Yweight::Vector{Float64}, cost::Union{Void, Matrix{Float64}} = nothing)
    if cost == nothing
        cost = ones(Float64, length(Xweight), length(Yweight))
    end
    flows = zeros(cost)
    distance = emd!(Xweight, Yweight, cost, flows)
    return distance, flows
end

end


text:
module GaussianBeamLens

# Set optimization level 1 for faster compilation, at potential cost of later performance
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

#include("GtkGUI.jl")

using FFTW: ifft, fft, ifftshift, fftshift

export GaussianBeamLensPropagate

const L = 1                      # wavelength of light
const n = 1                      # refractive index
const k = 2*pi*n/L               # wave number in material
const Nx = 1024                  # number of x cells
const dx = L/2                   # length of x cells
const Lx = Nx*dx                 # length of grid in x
const dkx = (2*pi)/Lx            # length of kx cell
const Nkx = Nx                   # same number of cells in real and fourier space

const x = dx.*(-(Nx/2):Nx/2-1)
const kx = dkx.*(-(Nkx/2):Nkx/2-1)
const kz = sqrt.(k^2 .- kx.^2)

const dz = 10                    # steps in z
const zmax = 6000                # end of grid in z
const z = 0:dz:zmax              # z grid

const FE = zeros(Complex{Float64}, length(z), length(x))

const cache = Dict{Tuple{Int64,Int64,Int64,Int64},Matrix{Float64}}()

"""
    GaussianBeamLensPropagate(x0,w0,f,Lens_z)
        x0:     Lateral offset of the input beam
        w0:     Width of the Gaussian beam in object space
        f:      Focal length of the lens
        Lens_z: Axial Position of the lens

Computes the intensity of the optical field at every point in x (lateral) and z (axial).
This method caches previously generated images in a dictionary at `GaussianBeamLens.cache`.
See GaussianBeamLensPropagate.GaussianBeamLensPropagateUncached for the actual implementation.
"""
function GaussianBeamLensPropagate(x0,w0,f,Lens_z)
    if haskey(cache,(x0,w0,f,Lens_z))
        return cache[(x0,w0,f,Lens_z)]
    else
        return GaussianBeamLensPropagateUncached(x0, w0, f, Lens_z)
    end
end

"""
    GaussianBeamLensPropagateUncached(x0,w0,f,Lens_z)
        x0:     Lateral offset of the input beam
        w0:     Width of the Gaussian beam in object space
        f:      Focal length of the lens
        Lens_z: Axial Position of the lens
Computes the intensity of the optical field at every point in x (lateral) and z (axial)
"""
function GaussianBeamLensPropagateUncached(x0, w0, f, Lens_z)
    E = exp.(-((x .- x0)./w0).^2)             # initial transverse electric field
    Tf = exp.(1im.*k.*(x.^2)/(2*f))           # transmission function for a thin lens

    # Propagate beam
    #FE = zeros(Complex{Float64},length(z),length(E))
    FE[1,:] = E

    for d = 1:length(z)-1
        # We could save some time if we saved the prior iteration in Fourier space
        FE[d+1,:] = ifft(ifftshift( fftshift(fft(FE[d,:])).*exp.(-1im.*kz.*dz) ))
        if abs(z[d]- Lens_z) < dz/2
            FE[d+1,:] .*= Tf
        end
    end
    I = abs2.(FE)                # intensity of the field at every point in x and z
    cache[(x0,w0,f,Lens_z)] = I
    return I
end

precompile(GaussianBeamLensPropagateUncached, (Int64, Int64, Int64, Int64))
precompile(GaussianBeamLensPropagate, (Int64, Int64, Int64, Int64))

end

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
export ARC_SecA

function ARC_SecA(h :: AbstractNLPModel;
				  kwargs...)

	(t, f, opt_res, iter, optimality, tired, status, hf, hg,hh) = ARC_generic(h; direction = "SecA", kwargs...)
	return (t, f, opt_res, iter, optimality, tired, status, hf, hg,hh)
end


text:
module BackpropNeuralNet

export
    # types
    NeuralNetwork,

    # functions
    init_network,
    train,
    net_eval

include("neural_net.jl")

end # module


text:
geom_sigdigits(::Type{Int64})::Int = 12
geom_sigdigits(::Type{Float32})::Int = 6
geom_sigdigits(::Type{Float64})::Int = 12

geom_atol_zero(::Type{Int64})::Float64 = 1e-12
geom_atol_zero(::Type{Float32})::Float32 = 1e-9
geom_atol_zero(::Type{Float64})::Float64 = 1e-12

function geom_round(v::T)::T where {T}
    return isapprox(v, 0, atol = geom_atol_zero(T)) ? T(0) : round(v, sigdigits = geom_sigdigits(T))
end
# function geom_round(cyl::CylindricalPoint{T})::CylindricalPoint{T} where T
#     return Cylindrical{T}(geom_round(cyl.r),geom_round(cyl.φ),geom_round(cyl.z))
# end

function geom_round(pt::CylindricalPoint{T})::CylindricalPoint{T} where {T <: SSDFloat}
    return CylindricalPoint{T}( geom_round(pt.r), geom_round(pt.φ), geom_round(pt.z)  )
end

function geom_round(pt::CartesianPoint{T})::CartesianPoint{T} where {T <: SSDFloat}
    return CartesianPoint{T}( geom_round(pt.x), geom_round(pt.y), geom_round(pt.z)  )
end


text:
# .title[Catwalk.jl: A profile guided dispatch optimizer]
#
# .log.blur[Krisztián Schäffer, independent researcher]
#
# .log.blur[JuliaCon 2021]
#
# ![](../../src/assets/catwalk-speeddemo.gif)
#
# ???
#
# Catwalk is an optimizing JIT compiler embedded in Julia.
# It can speed up long-running processes by minimizing the overhead of dynamic dispatch.
#
# But wait, Julia itself is already a JIT compiler, isn't?
#
# ---
# ### Ahead Of Time compiler
# ![](aot.jpg)
#
# ???
#
# A traditional compiler runs separately from the code it compiles, like
# when the machines lay the tracks first, and then execution, like a train,
# runs through the generated code.
#
# ---
# ### Julia: Just Ahead Of Time
# ![](jaot.jpg)
#
# ???
#
# Julia on the other hand can switch between execution and compilation.
# Code generation will be deferred to the latest possible point: Just before
# execution of that piece of code.
# The compiler has more information at that time, which it can use to
# optimize better.
#
# However, in Julia we typically only lay new tracks, and don't
# renew existing ones.
# Once a code was compiled, it is not really possible to recompile it,
# meaning that information gathered during execution of the code cannot be used
# to re-optimize it.
#
# ---
# ### Optimizing Just In Time
# ![](jit.png)
#
# ???
#
# And this re-optimization during runtime is what Catwalk does for you.
# It profiles dynamic dispatch in your code, measuring the distribution of
# dispatched types, and generates fast, type-stable routes for the most frequent ones.
#
# ---
# ### A Rock-Paper-Scissors "simulation"

abstract type Hand end
struct Rock <: Hand end
struct Paper <: Hand end
struct Scissors <: Hand end

play(::T, ::T) where T = 3       # 3: Tie
play(::Paper, ::Rock) = 1        # 1: First player wins
play(::Scissors, ::Paper) = 1
play(::Rock, ::Scissors) = 1
play(a, b) = 3 - play(b, a) # Reverse order. 3-1 == 2: second player wins

play(Paper(), Scissors()) # Scissors cuts paper, second player wins

# .small.right[Example was inspired by https://giordano.github.io/blog/2017-11-03-rock-paper-scissors/]
#
# ???
#
# Let's see how it works.
# This is a classic example of the rock-paper-scissors game.
#
# Hands are modeled with types and game logic is implemented using multiple dispatch.
# The play function returns one if the first player wins, two if the second,
# and three on tie.
#
# --
#

const hands = [Rock(), Paper(), Scissors()] # Vector{Hand}

function playrand(hands)
    hand1 = rand(hands)
    hand2 = rand(hands)
    return play(hand1, hand2) # Dynamic dispatch
end
nothing # hide

#
# ???
#
# The `playrand` function simulates a single play of two random-strategy players.
# By the nature of our design the call to `play` is dynamically dispatched at run-time,
# because all we know is that both hand1 and hand2 are subtypes of the
# abstract type Hand.
#
# This dynamic dispatch implements our main logic, and it is where most of the runtime is spent.
# We will try to speed it up.

# ---
# ### Play repeatedly and count results

function playmatch(hands, num_plays, result)
    for i=1:num_plays
        winner = playrand(hands)::Int
        result[winner] += 1
    end
end

result = [0, 0, 0]
playmatch(hands, 1000, result)

result

# ???
#
# Let's say that the `playrand` function is called in a hot loop.
#
# Why? Well, if it is not called in a hot loop, then there is hardly a
# reason to speed it up.
#
# So a match is a sequence of plays between two players. We count
# the results in an array.

# ---
#
# ### all-play-all tournament of num_players

function tournament(hands, num_players)
    results = [[0,0,0] for i1=1:num_players, i2=1:num_players]
    for i1 = 1:num_players
        for i2 = i1+1:num_players
            playmatch(hands, 1000, results[i1, i2])
        end
    end
    return results
end

tournament(hands, 5)

#
# ???
#
# Going further, a tournament is a matrix of matches.
# All-play-all, one thousand-long matches.
#
# --
#
# Performance:

@time tournament(hands, 300)
nothing # hide

# ???
#
# It is not very slow with more than 10 million games per second, but it can be much faster.
#
# ---
#
# ### Speed it up!

# ##### Option 1: Redesign ...

# ##### Option 2: Union spliting (does not work here) ...

# ##### Option 3: Dispatch manually:

function playrand_manualdispatch(hands)
    hand1 = rand(hands)
    hand2 = rand(hands)
    if hand1 isa Rock
        return play(hand1, hand2) # Fast, (partially) type-stable route
    elseif hand1 isa Paper
        return play(hand1, hand2)
    elseif hand1 isa Scissors
        return play(hand1, hand2)
    end
    return play(hand1, hand2) # Fallback to fully dynamic dispatch
end
nothing # hide

# (ManualDispatch.jl may help, but not in this 2-args case)
#
# ???
#
# When we do not want to redesign, and union splitting fails to work,
# then the first real solution is manual dispatch.
# The code seems strange at first, because all the branches do the same
# thing. But if you ask: "What is the semantics of performance optimization?"
# The answer is that it must be a no-op. The optimized
# version is equivalent with the original. The difference between the branches is
# not semantics, only implementation.
#
# The trick is that inside the branches the compiler knows the concrete type of hand1
# and will generate fast code based on this knowledge.
# You may have noticed that the type of hand2 is still unknown, so
# we still dispatch dynamically. But single dispatch is easier and the Julia compiler is very smart,
# it generates different code for the branches. Stabilizing `hand1` is enough for us now.
#
# The only problem with this solution is that we have to list the concrete types in the source code,
# which limits extensibility. If you want to introduce a new hand and still run fast,
# you have to change this code. That may not be feasible.
#
# So we arrive at Catwalk. It will automatically generate code that looks like this

# ---

using Catwalk

@jit play hand1 function playrand_jit(hands, jitctx)
    hand1 = rand(hands)
    hand2 = rand(hands)
    return play(hand1, hand2)
end

function playmatch_jit(hands, num_plays, result, jitctx)
    for i=1:num_plays
        winner = playrand_jit(hands, jitctx)::Int
        result[winner] += 1
    end
end

function tournament_jit(hands, num_players, jit=Catwalk.JIT())
    results = [[0,0,0] for i1=1:num_players, i2=1:num_players]
    for i1 = 1:num_players
        for i2 = i1+1:num_players
            Catwalk.step!(jit)
            playmatch_jit(hands, 1000, results[i1, i2], Catwalk.ctx(jit))
        end
    end
    return results
end
nothing # hide

# ???
#
# Let's see what we have to modify in our code.
# We mark `playrand` with the `@jit` macro and provide the name of the
# dynamically dispatched function, `play` and the argument to stabilize, `hand1`.
#
# Then we explicitly add the `jitctx` argument, which will drive the recompilation.
# When the Catwalk optimizer decides to compile a new version, it changes the type of
# the jit context, which triggers the recompilation. The type of the context
# describes the list of concrete types to stabilize.
#
# The body of the function is the same. `playmatch` just forwards the jit context
# explicitly, while at the outer level we also have to do some housekeeping.
# Catwalk works in batches: it can recompile between batches, and it profiles
# batches separately. Here a batch is a single match.
#
# ---
const jit = Catwalk.JIT()

@time tournament_jit(hands, 300, jit)
nothing # hide

# .log.blur[Debug: Found call site 'play', optimizing argument 'hand1'.]
# .log.blur[Debug: Profiled type freqs: Catwalk.Frequency[Catwalk.Frequency(Scissors, 334), Catwalk.Frequency(Paper, 330), Catwalk.Frequency(Rock, 336)\]]
# .log[Debug: Selected new fixtype list for 99.3% speedup vs last used, and 99.3% vs best: Rock, Scissors, Paper,]
# .log.blur[Debug: Profiled type freqs: Catwalk.Frequency[Catwalk.Frequency(Scissors, 364), Catwalk.Frequency(Paper, 324), Catwalk.Frequency(Rock, 312)\]]
# .log[Debug: Selected previous fixtype list: Rock, Scissors, Paper,]
# .log.blur[...]

#
# Hot start:

@time tournament_jit(hands, 300, jit)
nothing # hide

# ???
#
# And it is much faster now, we won around 50 percent. Of course, your mileage may vary.
# There is a huge initial compilation cost, so long-running processes
# will see higher wins.
#
# ---
# ### Reoptimization

function tournament_needing_reopt(hands, num_players, jit=Catwalk.JIT())
    results = [[0,0,0] for i1=1:num_players, i2=1:num_players]
    for i1 = 1:num_players
        hands = [Rock(), Paper(), (i1 > 20 ? Scissors() : Paper())] # No Scissors played by
        for i2 = i1+1:num_players                                   # the first 20 players
            Catwalk.step!(jit)
            playmatch_jit(hands, 1000, results[i1, i2], Catwalk.ctx(jit))
        end
    end
    return results
end
const reopt_jit = Catwalk.JIT()
@time tournament_needing_reopt(hands, 300, reopt_jit)
nothing # hide

# .log.blur[Debug: Profiled type freqs: Catwalk.Frequency[Catwalk.Frequency(Paper, 670), Catwalk.Frequency(Rock, 330)\]]
# .log[Debug: Selected new fixtype list for 99.4% speedup vs last used, and 99.4% vs best: **Paper, Rock,**]
# ...
# .log.blur[Debug: Profiled type freqs: Catwalk.Frequency[Catwalk.Frequency(Scissors, 298), Catwalk.Frequency(Paper, 351), Catwalk.Frequency(Rock, 351)\]]
# .log[Debug: Selected new fixtype list for 97.72% speedup vs last used, and 97.72% vs best: **Paper, Rock, Scissors,**]
#
# ???
#
# Catwalk was designed to handle situations where type distribution
# changes significantly during the run.
# This simple example shows it in action.
# The debug logs show that the first profile
# does not contain the Scissors type, but as it later appears, Catwalk recompiles.
#
# ---
# ### Catwalk.jl highlights:
#
#   - Low cost statistical profiler
#   - Tunable cost model
#   - Smart recompilation: Only if type-distribution changed significantly
#   - No world age issues, works inside a function
#   - 360 lines of readable vanilla Julia (no dependencies/ccall/llvmcall)
#
# Help, Tuning, inner workings, source of this example: https://tisztamo.github.io/Catwalk.jl/dev/
#

const configured_jit = Catwalk.JIT(;explorerfactory = Catwalk.NoExplorer)
Catwalk.add_boost!(
    configured_jit,
    Catwalk.CallBoost(
        :play,
        profilestrategy  =  Catwalk.SparseProfile(0.01),
        optimizer        =  Catwalk.TopNOptimizer(50;
                                compile_threshold = 1.04,
                                costmodel = Catwalk.DefaultDispatchCostModel(
                                    skip                = 4,
                                    static_dispatch     = 48,
                                    dynamic_dispatch    = 1000,
                                )
                            )
    ))

# .small.right[Configuration example]
#
# ???
#
# One more thing: We had no time to look inside, but there is no magic here.
# Even if you are a relative beginner to Julia, I suggest you to check the source
# code of Catwalk. It is simpler than you may think.
# But be warned: Metaprogramming is a rabbit hole!
#
# Thank you for your attention!
#


text:
using ITensors, ITensorGPU, Test
@testset "cuMPS Basics" begin
  N = 10
  sites = [Index(2, "Site") for n in 1:N]
  psi = cuMPS(sites)
  @test length(psi) == N
  @test length(cuMPS()) == 0

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteind(psi, 2) == sites[2]
  @test hasind(psi[3], linkind(psi, 2))
  @test hasind(psi[3], linkind(psi, 3))

  psi[1] = cuITensor(sites[1])
  @test hasind(psi[1], sites[1])

  L = randomMPS(sites)
  K = cuMPS(L)
  @test all(ITensors.data(cpu(K)) .== ITensors.data(cpu(L)))

  @testset "cuproductMPS" begin
    @testset "vector of string input" begin
      sites = siteinds("S=1/2", N)
      state = fill("", N)
      for j in 1:N
        state[j] = isodd(j) ? "Up" : "Dn"
      end
      psi = productCuMPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        ops = cuITensor(op(sites, "Sz", j))
        psip = prime(psi[j], "Site")
        res = psi[j] * ops * dag(psip)
        @test res[] ≈ sign / 2
      end
      @test_throws DimensionMismatch cuMPS(sites, fill("", N - 1))
      @test_throws DimensionMismatch productCuMPS(sites, fill("", N - 1))
    end

    @testset "vector of int input" begin
      sites = siteinds("S=1/2", N)
      state = fill(0, N)
      for j in 1:N
        state[j] = isodd(j) ? 1 : 2
      end
      psi = productCuMPS(sites, state)
      for j in 1:N
        sign = isodd(j) ? +1.0 : -1.0
        ops = cuITensor(op(sites, "Sz", j))
        psip = prime(psi[j], "Site")
        @test (psi[j] * ops * dag(psip))[] ≈ sign / 2
      end
    end
  end

  @testset "randomMPS" begin
    phi = randomCuMPS(sites)
    @test hasind(phi[1], sites[1])
    @test norm(phi[1]) ≈ 1.0
    @test hasind(phi[4], sites[4])
    @test norm(phi[4]) ≈ 1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1]) * psi[1]
    for j in 2:N
      phipsi *= dag(phi[j]) * psi[j]
    end
    @test phipsi[] ≈ inner(phi, psi)
    phi = randomCuMPS(sites)
    psi = randomCuMPS(sites)
    cphi = MPS([cpu(phi[i]) for i in 1:length(phi)])
    cpsi = MPS([cpu(psi[i]) for i in 1:length(psi)])
    phipsi = dag(phi[1]) * psi[1]
    cphipsi = dag(cphi[1]) * cpsi[1]
    for j in 2:N
      phipsi *= dag(phi[j]) * psi[j]
      cphipsi *= dag(cphi[j]) * cpsi[j]
    end
    @test cpu(phipsi)[] ≈ cphipsi[]
    @test cpu(phipsi)[] ≈ inner(cphi, cpsi)
    @test cpu(phipsi)[] ≈ inner(phi, psi)
    phipsi = dag(phi[1]) * psi[1]
    for j in 2:N
      phipsi = phipsi * dag(phi[j]) * psi[j]
    end
    @test cpu(phipsi)[] ≈ inner(phi, psi)

    badsites = [Index(2) for n in 1:(N + 1)]
    badpsi = randomCuMPS(badsites)
    @test_throws DimensionMismatch inner(phi, badpsi)
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(deepcopy(psi))
    ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1] * psi[1]
    for j in 2:N
      psipsi = psipsi * psidag[j] * psi[j]
    end
    @test psipsi[] ≈ inner(psi, psi)
    psi = randomCuMPS(sites)
    psidag = dag(deepcopy(psi))
    ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1] * psi[1]
    for j in 2:N
      psipsi = psipsi * psidag[j] * psi[j]
    end
    @test psipsi[] ≈ inner(psi, psi)
  end

  @testset "add MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    psi = randomCuMPS(sites)
    phi = deepcopy(psi)
    xi = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
  end

  sites = siteinds("S=1/2", N)
  psi = cuMPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomCuMPS(sites)
  orthogonalize!(psi, N - 1)
  @test ITensors.leftlim(psi) == N - 2
  @test ITensors.rightlim(psi) == N
  orthogonalize!(psi, 2)
  @test ITensors.leftlim(psi) == 1
  @test ITensors.rightlim(psi) == 3
  psi = randomCuMPS(sites)
  psi.rlim = N + 1 # do this to test qr from rightmost tensor
  orthogonalize!(psi, div(N, 2))
  @test ITensors.leftlim(psi) == div(N, 2) - 1
  @test ITensors.rightlim(psi) == div(N, 2) + 1

  #@test_throws ErrorException linkind(MPS(N, fill(cuITensor(), N), 0, N + 1), 1)

  @testset "replacebond!" begin
    # make sure factorization preserves the bond index tags
    psi = randomCuMPS(sites)
    phi = psi[1] * psi[2]
    bondindtags = tags(linkind(psi, 1))
    replacebond!(psi, 1, phi)
    @test tags(linkind(psi, 1)) == bondindtags

    # check that replaceBond! updates llim_ and rlim_ properly
    orthogonalize!(psi, 5)
    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="left")
    @test ITensors.leftlim(psi) == 5
    @test ITensors.rightlim(psi) == 7

    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="right")
    @test ITensors.leftlim(psi) == 4
    @test ITensors.rightlim(psi) == 6

    psi.llim = 3
    psi.rlim = 7
    phi = psi[5] * psi[6]
    replacebond!(psi, 5, phi; ortho="left")
    @test ITensors.leftlim(psi) == 3
    @test ITensors.rightlim(psi) == 7
  end
end

# Helper function for making MPS
function basicRandomCuMPS(N::Int; dim=4)
  sites = [Index(2, "Site") for n in 1:N]
  M = MPS(sites)
  links = [Index(dim, "n=$(n-1),Link") for n in 1:(N + 1)]
  M[1] = randomCuITensor(sites[1], links[2])
  for n in 2:(N - 1)
    M[n] = randomCuITensor(links[n], sites[n], links[n + 1])
  end
  M[N] = randomCuITensor(links[N], sites[N])
  M[1] /= sqrt(inner(M, M))
  return M
end
@testset "MPS gauging and truncation" begin
  N = 30

  @testset "orthogonalize! method" begin
    c = 12
    M = basicRandomCuMPS(N)
    orthogonalize!(M, c)

    @test ITensors.leftlim(M) == c - 1
    @test ITensors.rightlim(M) == c + 1

    # Test for left-orthogonality
    L = M[1] * prime(M[1], "Link")
    l = linkind(M, 1)
    @test cpu(L) ≈ delta(l, l') rtol = 1E-12
    for j in 2:(c - 1)
      L = L * M[j] * prime(M[j], "Link")
      l = linkind(M, j)
      @test cpu(L) ≈ delta(l, l') rtol = 1E-12
    end

    # Test for right-orthogonality
    R = M[N] * prime(M[N], "Link")
    r = linkind(M, N - 1)
    @test cpu(R) ≈ delta(r, r') rtol = 1E-12
    for j in reverse((c + 1):(N - 1))
      R = R * M[j] * prime(M[j], "Link")
      r = linkind(M, j - 1)
      @test cpu(R) ≈ delta(r, r') rtol = 1E-12
    end

    @test norm(M[c]) ≈ 1.0
  end

  @testset "truncate! method" begin
    M = basicRandomCuMPS(N; dim=10)
    M0 = copy(M)
    truncate!(M; maxdim=5)
    @test ITensors.rightlim(M) == 2
    # Test for right-orthogonality
    R = M[N] * prime(M[N], "Link")
    r = linkind(M, N - 1)
    @test cpu(R) ≈ delta(r, r') rtol = 1E-12
    for j in reverse(2:(N - 1))
      R = R * M[j] * prime(M[j], "Link")
      r = linkind(M, j - 1)
      @test cpu(R) ≈ delta(r, r') rtol = 1E-12
    end
    @test inner(M0, M) > 0.1
  end
end

#=@testset "Other MPS methods" begin

  @testset "sample! method" begin
    N = 10
    sites = [Index(3,"Site,n=$n") for n=1:N]
    psi = makeRandomCuMPS(sites,chi=3)
    nrm2 = inner(psi,psi)
    psi[1] *= (1.0/sqrt(nrm2))

    s = sample!(psi)

    @test length(s) == N
    for n=1:N
      @test 1 <= s[n] <= 3
    end

    # Throws becase not orthogonalized to site 1:
    orthogonalize!(psi,3)
    @test_throws ErrorException sample(psi)

    # Throws becase not normalized
    orthogonalize!(psi,1)
    psi[1] *= (5.0/norm(psi[1]))
    @test_throws ErrorException sample(psi)

    # Works when ortho & normalized:
    orthogonalize!(psi,1)
    psi[1] *= (1.0/norm(psi[1]))
    s = sample(psi)
    @test length(s) == N
  end

end=#


text:
## test_fem.jl : test GRF defined on FEM

@testset "finite element mesh      " begin

    ## Load node/element table ##
    p, t = star()
    m = Matern(0.2,2.0)

    ## Cholesky ##
    grf = GaussianRandomField(CovarianceFunction(2,m),Cholesky(),p,t)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.pts) == 2
    @test size(grf.pts[1],1) == 2
    @test size(grf.pts[1],2) == size(p,1)
    @test size(grf.pts[2],1) == 3
    @test size(grf.pts[2],2) == size(t,1)
    #plot_trisurf(grf); close()
    #tricontourf(grf); close()
    #plot(grf); close()

    ## Spectral ##
    grf = GaussianRandomField(CovarianceFunction(2,m),Spectral(),p,t)
    @test isa(grf,GaussianRandomField{Spectral})

    ## CirculantEmbedding ##
    @test_throws ArgumentError GaussianRandomField(CovarianceFunction(2,m),CirculantEmbedding(),p,t)

    ## KarhunenLoeve ##
    grf = GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="nodes")
    @test isa(grf,GaussianRandomField{KarhunenLoeve{100}})
    #plot_eigenvalues(grf); close()
    #plot_eigenfunction(grf,1); close()

    ## Anisotropic random field ##
    a = AnisotropicExponential([1000 0; 0 1000])
    cov = CovarianceFunction(2,a)
    grf = GaussianRandomField(cov,KarhunenLoeve(500),p,t,quad=GaussLegendre())
    @test isa(grf.cov.cov,AnisotropicExponential)

    ## Center mode ##
    grf = GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="center")
    @test isa(grf,GaussianRandomField{KarhunenLoeve{100}})
    @test_throws ArgumentError GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")
    @test_throws ArgumentError GaussianRandomField(1,CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")
    @test_throws ArgumentError GaussianRandomField(zeros(size(t,1)),CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")

    ## GaussianRandomField with mean value ##
    grf = GaussianRandomField(1,CovarianceFunction(2,m),Cholesky(),p,t,mode="center")
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.mean) == size(t,1)
    @test all(grf.mean.==1)

    grf = GaussianRandomField(ones(size(p,1)),CovarianceFunction(2,m),Spectral(),p,t)
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.mean) == size(p,1)
    @test all(grf.mean.==1)

    ## Non-bounding box method ##
    p, t = Lshape()
    grf = grf = GaussianRandomField(CovarianceFunction(2,Matern(0.2,1.0)),Spectral(),p,t,n=10)
    @test size(grf.data.eigenfunc,2) == 10

end


text:
using LennardJones
using Test

function CheckLattice(md::MD)
    nside = round( (md.n)^(1.0/3.0) );
    a = md.L/nside;
    n = md.n;
    L = md.L;
    Linv = 1.0/md.L;
    min_d2 = md.L;
    rij=zeros(3);
    for i in 1:(md.n-1)
        ri = md.r[:,i];
        for j in (i+1):md.n
            rj = md.r[:,j];
            rij = ri - rj;
            rij .-= L*round.(rij*Linv);
            d2 = sum(rij.*rij);
            min_d2 = (d2 < min_d2 ? d2 : min_d2);
        end
    end
    return sqrt(min_d2) ≈ a;
end

@testset "LennardJones.jl" begin
    @test CheckLattice( MD() );#Check the lattice spacing.
end


text:
# Supported backends
const _supportedbackends = Dict{Symbol,String}(
    :images => "Images",
    :luxor => "Luxor"
    #:cairomakie => "CairoMakie"
    #:glmakie => "GLMakie"
    #:wglmakie => "WGLMakie"
    )

# Initializaed backends
_backends = Dict{Symbol,Module}()

# Current backend
_backend_symbol = :images
_backend = ImagesBE
_backends[:images] = ImagesBE

backends() = keys(_backends)

backend() = _backend_symbol

"Set ready for use a supported graphics backend."
function backend(sym::Symbol)
    if sym in keys(_supportedbackends)
        if sym in keys(_backends) # if initialized
            global _backend_symbol = sym
            global _backend = _backends[sym]
            #_backend._init()
        else
            include(string(@__DIR__,"/backends/", sym, ".jl"))
            global _backend_symbol = sym
            namebe = Symbol(_supportedbackends[sym]*"BE")
            global _backend = @eval $namebe
            global _backends[sym] = _backend
            #_backend._init()
        end # if initialized
        #@info "Using backend :$_backend_symbol."
    else
        @warn """:$sym is not a supported backend.
            Current backend is :$_backend_symbol."""
    end # if supported
    _backend_symbol
end # function

function supported(funname)
    if !isdefined(_backend, funname)
        @error string(funname, " is not supported in ",  _supportedbackends[_backend_symbol], " backend")
        return false
    end
    true
end


text:
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: DRP_y, ∂DRP_∂y, directRP, ∂directRP_∂y;
using NonLinearReactionAdvectionDiffusionWithFrontData: ARP_y, ∂ARP_∂y, adjointRP, ∂adjointRP_∂y;
using BenchmarkTools
using LinearAlgebra

u_l(t) = -8 + 4*sin(2*π / T * t);# Прямая задача может быть с неоднородными ГУ
u_r(t) =  4 + sin(-2*π / T * t);# Но в дальнейшем, будем использовать только однородные.
q(x) = 4*sin(3 * π * x);        # Коэффициент линейного усиления, который в обратной
                                # задаче необходимо определить, но при генерации априорной
                                # информации мы задаем некоторый коэффициент, который,
                                # собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.28;                # Область по T
N, M = 50, 80;                  # Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      q.(Xₙ);               # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ);               # Начальные условия
nothing #hide

# Подготовим массивы, выбросив граничные точки
# ведь тестируемая функция — для внутреннего использования
X = Xₙ[2:N];
qq = qₙ[2:N];
y = u₀[2:N];

dl, d, du = DRP_y(y, 1, X, N, ε, ulₘ, urₘ, qq);

b = @benchmark DRP_y($y, 1, $X, $N, $ε, $ulₘ, $urₘ, $qq);
@info "Сформировать 3 вектора диагоналей якобиана прямой задачи:" b

b = @benchmark ∂DRP_∂y($y, 1, $X, $N, $ε, $ulₘ, $urₘ, $qq);
@info "Трехдиагональная матрица `∂DRP_∂y`:" b

b = @benchmark ∂directRP_∂y($y, 1, $X, $N, $ε, $ulₘ, $urₘ, $qq);
@info "Якобиан прямой задачи с помощью автоматического дифференциарования:" b

b = @benchmark solve($u₀, $Xₙ, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qₙ, $directRP, $∂DRP_∂y);
@info "Решение прямой задачи с трехдиагональным якобианом" b

b = @benchmark solve($u₀, $Xₙ, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qₙ, $directRP, $∂directRP_∂y);
@info "Решение прямой задачи с якобианом афтодифференцирования" b

# Готовимся к решению сопряженной задачи
u = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
#' ## Генерация априорной информации
ϕl = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                               # Левый вырожденный корень
ϕr = phidetermination(qₙ, urₘ, reverse(Xₙ), N, Tₘ, M);                      # Нужно подать инвертированную сетку
ϕr = reverse(ϕr, dims=1);                                                   # А после — инвертировать решение по X
ϕ = NonLinearReactionAdvectionDiffusionWithFrontData.Φ(ϕl, ϕr, N, M);       # Серединный корень
f1 = NonLinearReactionAdvectionDiffusionWithFrontData.f1(ϕ, u, Xₙ, N, M);   # Положение переходного слоя
f2 = NonLinearReactionAdvectionDiffusionWithFrontData.f2(f1, u, Xₙ, N, M);  # Значение функции на переходном слое

# Подготовим массивы, выбросив граничные точки
# ведь тестируемая функция — для внутреннего использования
U = u[2:N, :];

dl, d, du = ARP_y(y, 1, X, N, Tₘ, M, ε, ulₘ, urₘ, qq, U, f1, f2);

b = @benchmark ARP_y($y, 1, $X, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qq, $U, $f1, $f2);
@info "Сформировать 3 вектора диагоналей якобиана сопряженной задачи:" b

b = @benchmark ∂ARP_∂y($y, 1, $X, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qq, $U, $f1, $f2);
@info "Трехдиагональная матрица `∂ARP_∂y`:" b

# Здесь значения якобиана вычисляются неправильно, ведь нам нужно еще
# развернуть вектора f1, f2 и U(по второй координате).
# Помните — Функция для внутреннего использования, поэтому сейчас мы
# должны сделать это сами.
b = @benchmark ∂adjointRP_∂y($y, 1, $X, $N, $reverse(Tₘ), $M, $ε, $ulₘ, $urₘ, $qq, $U, $f1, $f2);
@info "Якобиан сопряженной задачи с помощью автоматического дифференциарования:" b

b = @benchmark solve_adjoint($u₀, $Xₙ, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qₙ, $u, $f1, $f2, $adjointRP, $∂ARP_∂y);
@info "Решение сопряженной задачи с трехдиагональным якобианом" b

b = @benchmark solve_adjoint($u₀, $Xₙ, $N, $Tₘ, $M, $ε, $ulₘ, $urₘ, $qₙ, $u, $f1, $f2, $adjointRP, $∂adjointRP_∂y);
@info "Решение сопряженной задачи с якобианом афтодифференцирования" b

############################# Benchmark minimization loop #####################

q₀ = [ 0 for x in qₙ];
ψ₀ = zeros(N+1);
ψl = zeros(M+1);
ψr = zeros(M+1);
S = 10;
β = 0.1;

qf, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, ψ₀, ψl, ψr, Xₙ, N, Tₘ, M, ε, f1, f2, S = S, β = β)

b = @benchmark minimize($q₀, $u₀, $ulₘ, $urₘ, $ψ₀, $ψl, $ψr, $Xₙ, $N, $Tₘ, $M, $ε, $f1, $f2, S = $S, β = $β)
@info "Десять итераций цикла минимизации" b


text:
using simulationQP934

data = Data()
md1  = MDinput(dt=1)
md2  = MDinput(dt=0.1)
md3  = MDinput(dt=0.01)
md4  = MDinput(dt=0.001)

println("calculos para dt = 1")
md(data,md1)
println("  ")

println("calculos para dt = 0.1")
md(data,md2)
println("  ")

println("calculos para dt = 0.01")
md(data,md3)
println("  ")

println("calculos para dt = 0.001")
md(data,md4)
println("  ")


text:
using Test
using LoggingTestSets
using AsyncLog
using UnixIO
using UnixIO: C

# FIXME use ReTest.jl
"""
# Test Cases to add

- Close while transfer is waiting.

"""

cd(@__DIR__)

@testset LoggingTestSet "UnixIO" begin

@testset LoggingTestSet "Open Type" begin

f() = UnixIO.open(UnixIO.FDType, "foobar")
info, type = code_typed(f, ())[1]
@test type == DuplexIO{UnixIO.FD{UnixIO.In,UnixIO.FDType},
                       UnixIO.FD{UnixIO.Out,UnixIO.FDType}}

f() = UnixIO.open(UnixIO.FDType, "foobar", C.O_RDWR)
info, type = code_typed(f, ())[1]
@test type == DuplexIO{UnixIO.FD{UnixIO.In,UnixIO.FDType},
                       UnixIO.FD{UnixIO.Out,UnixIO.FDType}}

f() = UnixIO.open(UnixIO.FDType, "foobar", C.O_RDONLY)
info, type = code_typed(f, ())[1]
@test type == UnixIO.FD{UnixIO.In,UnixIO.FDType}

f() = UnixIO.open(UnixIO.FDType, "foobar", C.O_WRONLY)
info, type = code_typed(f, ())[1]
@test type == UnixIO.FD{UnixIO.Out,UnixIO.FDType}

end # testset

@info "Test read(::Cmd)"
@test UnixIO.read(`uname -a`) ==
             read(`uname -a`)

@info "Test open, readline, read, readbytes!, eof, isopen (with file)"
jio = open("runtests.jl")
uio = UnixIO.open("runtests.jl")

@test readline(jio) == readline(uio)
@test read(jio, UInt8) == read(uio, UInt8)
@test read(jio, UInt32) == read(uio, UInt32)
@test readline(jio) == readline(uio)
jv = UInt8[]; readbytes!(jio, jv, 2)
uv = UInt8[]; readbytes!(uio, uv, 2)
@test jv == uv
jv = UInt8[]; readbytes!(jio, jv, 100)
uv = UInt8[]; readbytes!(uio, uv, 100)
@test jv == uv
readbytes!(jio, jv, 100)
readbytes!(uio, uv, 100)
@test jv == uv
@test read(jio) == read(uio)
@test eof(jio) == eof(uio)
@test isopen(jio) == isopen(uio)
@test close(jio) == close(uio)
@test eof(jio) == eof(uio)
@test isopen(jio) == isopen(uio)

@info "Test system() and read(::Cmd) with large data"
for _ in 1:1
    mktempdir() do d
        f = joinpath(d, "testdata")
        UnixIO.system("dd if=/dev/urandom of=$f bs=1024 count=10000 " *
                      "> /dev/null 2>&1")
        @test stat(f).size == 1024 * 10000
        UnixIO.system("hexdump $f > $f.hex")

        a = read("$f.hex")
        b = read(`hexdump $f`)
        c = UnixIO.read(`hexdump $f`)
        d = UnixIO.open(`hexdump`) do cin, cout
            @sync begin
                @async begin
                    write(cin, UnixIO.open(f))
                    close(cin)
                end
                read(cout)
            end
        end

        @test b == a
        @test c == a
        @test d == a
        if a != c
            @show length(a) length(b) length(c)
            a = String(a)
            b = String(b)
            c = String(c)
            @show a[1:100] a[end-100:end]
            @show b[1:100] b[end-100:end]
            @show c[1:100] c[end-100:end]
        end
    end
end

@info "Test eachline()"
jio = open("runtests.jl")
uio = UnixIO.open("runtests.jl")

@test [x for x in eachline(jio)] ==
      [x for x in eachline(uio)]

@info "Test open(:Cmd)"
@test UnixIO.open(`hexdump`) do cmdin, cmdout
    @async begin
        write(cmdin, read(UnixIO.open("runtests.jl")))
        close(cmdin)
    end
    read(cmdout)
end ==
open(`hexdump`, open("runtests.jl"); read = true) do io
    read(io)
end

@info "Test read(`hexdump`)"
@test UnixIO.read(`hexdump runtests.jl`) ==
             read(`hexdump runtests.jl`)

cmd = `bash -c "echo FOO; sleep 1; echo BAR"`
@info "Test open($cmd)"
@test UnixIO.open(cmd) do i, o collect(eachline(o)) end ==
    open(cmd; read = true) do io collect(eachline(io)) end

@info "Test @async read(::Cmd) with delay"
results = []
@sync for i in 1:5
    @async try
        push!(results, (i, UnixIO.read(`bash -c "sleep 2; echo $i"`, String)))
    catch err
        UnixIO.printerr("ERROR: $err")
        exception=(err, catch_backtrace())
        UnixIO.printerr(exception)
        @error "ERROR" exception
    end
end
for (i, r) in results
    @test r == "$i\n"
end

@info "Test read(::Cmd) with timeout."
@test UnixIO.read(`bash -c "echo FOO; sleep 4; echo FOO"`,
                  String; timeout=2.4) == "FOO\n"

times = []
@sync for i in 1:4
    @async try
        t0 = time()
        UnixIO.read(`bash -c "sleep 4; echo FOO"`, String; timeout=2.4)
        push!(times, time() - t0)
    catch err
        UnixIO.printerr("ERROR: $err")
        exception=(err, catch_backtrace())
        UnixIO.printerr(exception)
        @error "ERROR" exception
    end
end
@show times
for t in times
    @test abs(t - 2.4) < 0.2
end

@info "Test open(::Cmd) -> cin, cout."
sleep(1)
@test isempty(UnixIO.processes)
cmdin, cmdout = UnixIO.open(`bash -c "while true; do date ; sleep 1; done"`)
close(cmdin)
@test !isempty(UnixIO.processes)
@test readline(cmdout) != ""
@test readline(cmdout) != ""
@test readline(cmdout) != ""

@info "Test sub-process cleanup."
close(cmdout)
sleep(1)
@test isempty(UnixIO.processes)

include("pseudoterminal.jl")

end #testset


text:
using Test
include("module-2-exercise-3.jl")

@testset "$(get_module_2_excersise_3_author_name()): count_even_numbers tests" begin
    @test 0 == count_even_numbers([])
    @test 0 == count_even_numbers([1])
    @test 1 == count_even_numbers([2])
    @test 2 == count_even_numbers([2,4])
    @test 2 == count_even_numbers([2,4,1])
    @test 1 == count_even_numbers([1,3,5,2,3])
end


text:
"""
    AbstractInterpolatingPolynomial{T,X}

Abstract type for interpolating polynomials.

These are polynomial representations of `p(x)` satisfying `p(x_i) =
y_i` for a specified set of `x` values and `y` values.

For a collection of points `(x_0,y_0), ..., (x_n, y_n)` there is a
*unique* polynomial of degree `n` or less satisfying
`p(x_i)=y_i`. This fact allows the specification of `p(x)` using a
vector of coefficients relative to some set of basis vectors.

The two main types, `Lagrange` and `Newton`, store the nodes within
the instance. In particular, the type does not contain all the
information needed to describe the family. So methods like
`convert(::Type, p)` will not work. Use `fit(Type, xs, p)`, as
appropriate, instead.

"""
abstract type AbstractInterpolatingPolynomial{T,X} <: AbstractSpecialPolynomial{T,X} end

function Polynomials.showterm(
    io::IO,
    ::Type{P},
    pj::T,
    var,
    j,
    first::Bool,
    mimetype,
) where {N,T,P<:AbstractInterpolatingPolynomial}
    iszero(pj) && return false
    !first && print(io, " ")
    print(io, Polynomials.hasneg(T) && Polynomials.isneg(pj) ? "- " : (!first ? "+ " : ""))
    print(io, "$(abs(pj))⋅$(basis_symbol(P))_$(j)($var)")
    return true
end

basis_symbol(::Type{<:AbstractInterpolatingPolynomial}) = "P"

function Polynomials.variable(
    ::Type{P},
    var::Polynomials.SymbolLike=:x,
) where {P<:AbstractInterpolatingPolynomial}
    @warn "No `variable` function is defined for a type, just an instance"
    throw(MethodError())
end

# pick new nodes for product
function _new_nodes(xs::Vector{S}, ys::Vector{T}) where {S,T}
    sort!(ys)
    n = length(ys)
    out = setdiff(ys, xs)
    length(out) >= n - 1 && return out[1:(n - 1)]
    out = out / 1
    tmp = intersect(xs, ys)
    sort!(tmp)
    for i in 2:length(tmp)
        a, b = tmp[i - 1], tmp[i]
        c = (a + b) / 2
        while c in xs
            c = (a + c) / 2
        end
        push!(out, c)
    end
    out
end


text:
### Intelligent pedestrian model that stops if a crosswalk is busy

@with_kw mutable struct IntelligentPedestrian <: DriverModel{ConstantSpeedDawdling}
    a::ConstantSpeedDawdling = ConstantSpeedDawdling(0., 0.)
    v_noise::Float64 = 0.
    dt::Float64 = 0.1
    motion::DriverModel = ConstantPedestrian(v_noise = v_noise, dt = dt)
    ttc_threshold::Float64 = 4.5
    crosswalk::Lane = Lane()
    conflict_lanes::Vector{Lane{Float64}} = Lane{Float64}[]
    dtol::Float64 = 0.1
    # states
    wait::Bool = true
    crossing::Bool = false
    done::Bool = false
end

AutomotiveDrivingModels.get_name(model::IntelligentPedestrian) = "IntelligentPedestrian"
Base.rand(model::IntelligentPedestrian) = model.a

AutomotiveDrivingModels.reset_hidden_state!(model::IntelligentPedestrian) = model

function AutomotiveDrivingModels.observe!(model::IntelligentPedestrian,
                                        scene::EntityFrame{VehicleState, VehicleDef, Int64},
                                        roadway::Roadway,
                                        egoid::Int)
    ped = scene[findfirst(egoid, scene)]
    observe!(model.motion, scene, roadway, egoid)
    model.crossing = update_crossing(ped, model.crosswalk, model.conflict_lanes, roadway)
    model.done = update_done(ped, model.crossing, model.crosswalk)
    if !model.crossing
        model.wait = update_wait(ped, scene, model, roadway)
    else
        model.wait = false
    end
    if !model.wait || model.done# if engaged, keep going
        model.a = model.motion.a
    else
        if model.wait
            model.a = ConstantSpeedDawdling(0., 0.)
        else
            model.a = model.motion.a
        end
    end
end

function update_crossing(ped::Vehicle, crosswalk::Lane{T}, conflict_lanes::Vector{Lane{T}}, roadway::Roadway) where T
    # check if the pedestrian is in the conflict zone
    for lane in conflict_lanes
        ped_f = Frenet(proj(ped.state.posG, lane, roadway, move_along_curves=false), roadway)
        if abs(ped_f.t) < lane.width/2 && get_lane(roadway, ped).tag == crosswalk.tag
            # println(" position on lane : ", ped_f)
            return true
        end
    end
    return false
end

function update_done(ped::Vehicle, crossing::Bool, crosswalk::Lane)
    # check if the pedestrian is done crossing
    return !crossing && direction_from_center(ped, crosswalk) .< 0.
end

function get_distance_to_crosswalk(veh::Vehicle, crosswalk::Lane, roadway::Roadway)
    lane = get_lane(roadway, veh)
    cw_length = get_end(crosswalk)
    cw_center = posg(Frenet(crosswalk, cw_length/2), roadway)
    collision_point = VecSE2(cw_center.x, veh.state.posG.y)
    cw_proj = Frenet(collision_point, lane, roadway)
    Δs = cw_proj.s - veh.state.posF.s
    return Δs
end

function update_wait(ped::Vehicle,
                    scene::Scene,
                    model::IntelligentPedestrian,
                    roadway::Roadway)
    # check if the pedestrian should wait or not
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN
            continue
        end
        # check if vehicle is on the crosswalk
        posG = veh.state.posG
        cw_posF = Frenet(proj(posG, model.crosswalk, roadway, move_along_curves=false), roadway)
        lane = get_lane(roadway, veh)
        # if lane ∈ model.conflict_lanes
        # println("vehicle position on cw $(cw_posF)")
        # println("pedestrian position $(ped.state.posF)")
        # println(" difference: $(cw_posF.s - ped.state.posF.s)")
        # t_car = cw_posF.t
        # println( abs(cw_posF.t) < model.crosswalk.width/2 + veh.def.length/2 - model.dtol)
        # println( (cw_posF.s - ped.state.posF.s)*cos(ped.state.posF.ϕ) > 0.)
        if abs(cw_posF.t) < model.crosswalk.width/2 + veh.def.length/2 - model.dtol &&
                                (cw_posF.s - ped.state.posF.s)*cos(ped.state.posF.ϕ) > 0.
            # println("vehicle $(veh.id) is on the crosswalk")
            # vehicle is on the crosswalk
            # check if pedestrian is passed it or not
            return true
        end
        if lane ∈ model.conflict_lanes
            # check the ttc
            Δs = get_distance_to_crosswalk(veh, model.crosswalk, roadway)
            ttc = Δs/veh.state.v
            if 0. < ttc < model.ttc_threshold
                # println("vehicle $(veh.id) ttc: $ttc, Δs: $Δs")
                return true
            end
        end
    end
    return false
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs