text:
module M

include("CodeGen.jl")
using .CodeGen: PType, Struct, Field, Method, GlobalMethod, operator
using .CodeGen: @field, @type, write_cpp_typestub, write_jl_typestub, write_cpp_methodstub

const CPP_DIR = abspath(joinpath(@__DIR__, "../ibapi-wrapper/include/generated"))
const JL_DIR = abspath(joinpath(@__DIR__, "../src/IBKR/IBKRCore/generated"))
@assert isdir(abspath(joinpath(CPP_DIR, "..")))
@assert isdir(abspath(joinpath(JL_DIR, "..")))

const MOD = :mod

rm(CPP_DIR; recursive=true, force=true)
rm(JL_DIR; recursive=true, force=true)
mkdir(CPP_DIR)
mkdir(JL_DIR)

jlfile(t::Struct) = "$(t.name).jl"
cppfile(t::Struct) = "$(t.name).h"
function write_stub(type::Struct, methods::Vector{GlobalMethod} = GlobalMethod[];  extra_cpp_typestub = [], extra_cpp_methodstub = [])
    write_jl_typestub(joinpath(JL_DIR, "$(type.name).jl"), type)
    open(joinpath(CPP_DIR, "$(type.name).h"), "w") do io
        write_cpp_typestub(io, MOD, type, extra=extra_cpp_typestub)
        write_cpp_methodstub(io, MOD, methods, extra=extra_cpp_methodstub)
    end
end

let T = Struct(
        :ComboLeg;
        fields=[
            Field(:conId, :CxxLong),
            Field(:ratio, :CxxLong),
            Field(:action, :StdString),
            Field(:exchange, :StdString),
            Field(:openClose, :CxxLong),
            Field(:shortSaleSlot, :CxxLong),
            Field(:designatedLocation, :StdString),
            Field(:exemptCode, :CxxLong),
        ],
    )
    write_stub(T, [
        operator(:(==), T),
        GlobalMethod(:ComboLegSPtr, "[]() { return ComboLegSPtr(new ComboLeg); }"),
        GlobalMethod(:showComboLeg, "&showComboLeg"),
        ],
        extra_cpp_methodstub = [
            "jlcxx::stl::apply_stl<ComboLeg>(mod)"
        ],
    )
end

let T = Struct(
        :Contract;
        fields=[
            Field(:conId, :CxxLong),
            Field(:symbol, :StdString),
            Field(:secType, :StdString),
            Field(:lastTradeDateOrContractMonth, :StdString),
            Field(:strike, :Cdouble),
            Field(:right, :StdString),
            Field(:multiplier, :StdString),
            Field(:exchange, :StdString),
            Field(:primaryExchange, :StdString),
            Field(:currency, :StdString),
            Field(:localSymbol, :StdString),
            Field(:tradingClass, :StdString),
            Field(:includeExpired, :CxxBool),
            Field(:secIdType, :StdString),
            Field(:secId, :StdString),

            Field(:comboLegsDescrip, :StdString),
            Field(:comboLegs, PType(:SharedPtr, PType(:StdVector, PType(:SharedPtr, :ComboLeg))), mutable=false),
            # Field(:comboLegs, PType(:StdVector, PType(:SharedPtr, :ComboLeg))),
        ],
    )
    write_stub(T)
end

end


text:
using DemoPackageBXW
using Test

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

ext: jl


text:
@with_kw_noshow mutable struct Penalty{T1 <: AbstractFloat}
    name :: Symbol
    isUsed :: Bool = true
    scalarValue :: T1
    # Function that computes `scalarValue`
    scalarFct :: Function
    showFct :: Union{Function, Nothing} = nothing
end

mutable struct PenaltyVector
    pv :: Vector{Penalty}
end

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


text:
module Polynomials

#  using GenericLinearAlgebra ## remove for now. cf: https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/71#issuecomment-743928205
using LinearAlgebra
import Base: evalpoly

include("abstract.jl")
include("show.jl")
include("plots.jl")
include("contrib.jl")

# Interface for all AbstractPolynomials
include("common.jl")

# Polynomials
include("polynomials/standard-basis.jl")
include("polynomials/mutable-arithmetics.jl")
include("polynomials/Polynomial.jl")
include("polynomials/ImmutablePolynomial.jl")
include("polynomials/SparsePolynomial.jl")
include("polynomials/LaurentPolynomial.jl")
include("polynomials/pi_n_polynomial.jl")
include("polynomials/factored_polynomial.jl")
include("polynomials/ngcd.jl")
include("polynomials/multroot.jl")

include("polynomials/ChebyshevT.jl")

# Rational functions
include("rational-functions/common.jl")
include("rational-functions/rational-function.jl")
include("rational-functions/fit.jl")
#include("rational-functions/rational-transfer-function.jl")
include("rational-functions/plot-recipes.jl")

# compat; opt-in with `using Polynomials.PolyCompat`
include("polynomials/Poly.jl")

end # module


text:
@defstruct FMmodel (
  w0::Float64                = 0.0,
  w::Vector{Float64}         = Vector{Float64}(),
  v::Array{Float64}          = Array{Float64}(),

  # controls
  (SOLVER::String            = "TDAP", SOLVER in ["ALS" "MCMC" "SGD" "BGD" "FTRL" "TDAP"]),
  (TASK::String              = "CLASSIFICATION", TASK in ["CLASSIFICATION" "REGRESSION"]),
  (nthreads::Int             = 1, nthreads >= 1),
  (num_attribute::Int        = 1, num_attribute > 0),
  k0::Bool                   = true,
  k1::Bool                   = true,
  (num_factor::Int           = 2, num_factor >= 0),
  (l1_regw::Float64          = 0.0, l1_regw >= 0),
  (l1_regv::Float64          = 0.0, l1_regv >= 0),
  (l2_regw::Float64          = 0.0, l2_regw >= 0),
  (l2_regv::Float64          = 0.0, l2_regv >= 0),
  (l2_reg0::Float64          = 0.0, l2_reg0 >= 0),
  init_mean::Float64         = 0.0,
  (init_stdev::Float64       = 0.0, init_stdev >= 0),

  # temporary
  m_sum::Vector{Float64}     = Vector{Float64}(),
  m_sum_sqr::Vector{Float64} = Vector{Float64}()
)

function InitModel(model::FMmodel)
  model.w = zeros(Float64, model.num_attribute)
  model.v = zeros(Float64, (model.num_factor, model.num_attribute))
  model.m_sum = zeros(Float64, model.num_factor)
  model.m_sum_sqr = zeros(Float64, model.num_factor)
end

function SaveModel(model::FMmodel, file::String)
  checkPath(file)
  open(file, "w+") do io
    write(io, "FMmodel\n")
    for field in fieldnames(model)
      writeKV(io, string(field), getfield(model, field))
    end
  end
end

function LoadModel(file::String)
  checkPath(file)
  isfile(file) || error(file * " does not exist")
  open(file, "r") do io
    ftype = readline(io)
    @assert(ftype == "FMmodel\n", "the file read in is not a FMmodel output")
    for line in eachline(io)
      key, value = readKV(line)
      set(model, convert(String, key), value)
    end
  end
  return model
end

function predict(model::FMmodel, vec::SparseVector)
  pred::Float64 = 0.0
  model.k0 && pred += model.w0
  model.m_sum = zeros(model.m_sum)
  model.m_sum_sqr = zeros(model.m_sum_sqr)

  for (_idx, _val) in zip(vec.nzind, vec.nzval):
    model.k1 &&　pred += model.w[_idx] * _val
    for i in 1:model.num_factor
      _tmp = model.v[i, _idx] * _val
      model.m_sum[i] += _tmp
      model.m_sum_sqr[i] += _tmp^2
    end
  end

  for i in 1:model.num_factor
    pred += 0.5 * (model.m_sum[i] ^ 2 - model.m_sum_sqr[i])
  end

  return pred
end


text:
using .StencilBuilder

## Kernels.
function dynamics_u_kernel!(
    ut, u, v, w,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    @fast3d begin
        @fd (ut, u, v, w) ut += (
            - gradx(interpx(u) * interpx(u)) + visc * (gradx(gradx(u)))
            - grady(interpx(v) * interpy(u)) + visc * (grady(grady(u)))
            - gradz(interpx(w) * interpz(u)) + visc * (gradz(gradz(u))) )
    end
end

function dynamics_v_kernel!(
    vt, u, v, w,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    @fast3d begin
        @fd (vt, u, v, w) vt += (
            - gradx(interpy(u) * interpx(v)) + visc * (gradx(gradx(v)))
            - grady(interpy(v) * interpy(v)) + visc * (grady(grady(v)))
            - gradz(interpy(w) * interpz(v)) + visc * (gradz(gradz(v))) )
    end
end

function dynamics_w_kernel!(
    wt, u, v, w, s,
    visc, alpha,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke) # CvH: ke represents keh here, but that is not yet supported by fast3d, needs fixing.

    @fast3d begin
        @fd (wt, u, v, w, s) wt += (
            - gradx(interpz(u) * interpx(w)) + visc * (gradx(gradx(w)))
            - grady(interpz(v) * interpy(w)) + visc * (grady(grady(w)))
            - gradz(interpz(w) * interpz(w)) + visc * (gradz(gradz(w)))
            + alpha*interpz(s) )
    end
end

function dynamics_s_kernel!(
    st, u, v, w, s,
    visc,
    dxi, dyi, dzi, dzhi,
    is, ie, js, je, ks, ke)

    @fast3d begin
        @fd (st, u, v, w, s) st += (
            - gradx(u * interpx(s)) + visc * (gradx(gradx(s)))
            - grady(v * interpy(s)) + visc * (grady(grady(s)))
            - gradz(w * interpz(s)) + visc * (gradz(gradz(s))) )
    end
end

## Dynamics kernels launcher.
function calc_dynamics_tend!(f::Fields, g::Grid)
    dynamics_u_kernel!(
        f.u_tend, f.u, f.v, f.w,
        f.visc,
        g.dxi, g.dyi, g.dzi, g.dzhi,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)

    dynamics_v_kernel!(
        f.v_tend, f.u, f.v, f.w,
        f.visc,
        g.dxi, g.dyi, g.dzi, g.dzhi,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)

    dynamics_w_kernel!(
        f.w_tend, f.u, f.v, f.w, f.s,
        f.visc, f.alpha,
        g.dxi, g.dyi, g.dzi, g.dzhi,
        # g.is, g.ie, g.js, g.je, g.ks, g.keh) # CvH temporary hack
        g.is, g.ie, g.js, g.je, g.ks+1, g.keh-1)

    dynamics_s_kernel!(
        f.s_tend, f.u, f.v, f.w, f.s,
        f.visc,
        g.dxi, g.dyi, g.dzi, g.dzhi,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)
end


text:
# Pointers with address space information

# Address spaces

export AS, addrspace

abstract type AddressSpace end

module AS

import ..AddressSpace

struct Generic  <: AddressSpace end
struct Global   <: AddressSpace end
struct Region   <: AddressSpace end
struct Local    <: AddressSpace end
struct Constant <: AddressSpace end
struct Private  <: AddressSpace end

end # module

# Device pointer

"""
    DevicePtr{T,A}
"""
DevicePtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type DevicePtr{T,A} <: Ref{T} 64 end
else
    primitive type DevicePtr{T,A} <: Ref{T} 32 end
end

# constructors

DevicePtr{T,A}(x::Union{Int,UInt,Ptr,DevicePtr}) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, x)
DevicePtr{T}(ptr::Ptr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)
DevicePtr(ptr::Ptr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

# getters

Base.eltype(::Type{<:DevicePtr{T}}) where {T} = T
addrspace(x::DevicePtr) = addrspace(typeof(x))
addrspace(::Type{DevicePtr{T,A}}) where {T,A} = A

# conversions

# to and from integers
# pointer to integer
Base.convert(::Type{T}, x::DevicePtr) where {T<:Integer} = T(UInt(x))
# integer to pointer
Base.convert(::Type{DevicePtr{T,A}}, x::Union{Int,UInt}) where {T,A<:AddressSpace} = DevicePtr{T,A}(x)
Int(x::DevicePtr)  = Base.bitcast(Int, x)
UInt(x::DevicePtr) = Base.bitcast(UInt, x)

# between host and device pointers
Base.convert(::Type{Ptr{T}},  p::DevicePtr)  where {T}                 = Base.bitcast(Ptr{T}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::Ptr) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T}}, p::Ptr)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{Ptr{T}}, x::DevicePtr{T}) where {T} = reinterpret(Ptr{T}, x)

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = Base.bitcast(DevicePtr{T,A}, p)

# identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

# convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities

# unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

# limited pointer arithmetic & comparison
isequal(x::DevicePtr, y::DevicePtr) = (x === y) && addrspace(x) == addrspace(y)
isless(x::DevicePtr{T,A}, y::DevicePtr{T,A}) where {T,A<:AddressSpace} = x < y
Base.:(==)(x::DevicePtr, y::DevicePtr) = UInt(x) == UInt(y) && addrspace(x) == addrspace(y)
Base.:(<)(x::DevicePtr,  y::DevicePtr) = UInt(x) < UInt(y)
Base.:(-)(x::DevicePtr,  y::DevicePtr) = UInt(x) - UInt(y)
Base.:(+)(x::DevicePtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::DevicePtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::DevicePtr) = y + x

# memory operations

@static if Base.libllvm_version < v"7.0"
    # Old values (LLVM 6)
    Base.convert(::Type{Int}, ::Type{AS.Private})  = 0
    Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
    Base.convert(::Type{Int}, ::Type{AS.Constant}) = 2
    Base.convert(::Type{Int}, ::Type{AS.Local})    = 3
    Base.convert(::Type{Int}, ::Type{AS.Generic})  = 4
    Base.convert(::Type{Int}, ::Type{AS.Region})   = 5
else
    # New values (LLVM 7+)
    Base.convert(::Type{Int}, ::Type{AS.Generic})  = 0
    Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
    Base.convert(::Type{Int}, ::Type{AS.Region})   = 2
    Base.convert(::Type{Int}, ::Type{AS.Local})    = 3
    Base.convert(::Type{Int}, ::Type{AS.Constant}) = 4
    Base.convert(::Type{Int}, ::Type{AS.Private})  = 5
end

function tbaa_make_child(name::String, constant::Bool=false; ctx::LLVM.Context=JuliaContext())
    tbaa_root = MDNode([MDString("roctbaa", ctx)], ctx)
    tbaa_struct_type =
        MDNode([MDString("roctbaa_$name", ctx),
                tbaa_root,
                LLVM.ConstantInt(0, ctx)], ctx)
    tbaa_access_tag =
        MDNode([tbaa_struct_type,
                tbaa_struct_type,
                LLVM.ConstantInt(0, ctx),
                LLVM.ConstantInt(constant ? 1 : 0, ctx)], ctx)

    return tbaa_access_tag
end

tbaa_addrspace(as::Type{<:AddressSpace}) = tbaa_make_child(lowercase(String(as.name.name)))

@generated function Base.unsafe_load(p::DevicePtr{T,A}, i::Integer=1,
                                     ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)
    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})
    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ld = load!(builder, ptr)

        if A != AS.Generic
            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end

        alignment!(ld, align)
        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,A}, Int}, :((p, Int(i-one(i)))))
end

@generated function Base.unsafe_store!(p::DevicePtr{T,A}, x, i::Integer=1,
                                       ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)
    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})
    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, eltyp, T_int]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[3]])
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)

        if A != AS.Generic
            metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end

        alignment!(st, align)
        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{DevicePtr{T,A}, T, Int},
                  :((p, convert(T,x), Int(i-one(i)))))
end


text:
using StorageGraphs: findnodes

@testset "Data Query" begin
    g = StorageGraph()

    dep = (A=1,)=>(D=0.4,)=>(B=0.5,)
    add_bulk!(g, dep, (E=1:3,))
    dep = (A=1,)=>(D=0.4,)=>(B=0.6,)
    add_bulk!(g, dep, (E=1:3,))

    @test get_prop(g, g[(B=0.5,)]) == (B=0.5,)
    @test get_prop(g, g[(B=0.6,)]) == (B=0.6,)

    @test isempty(setdiff(findnodes(g, :A), [(A=1,)]))
    @test isempty(setdiff(findnodes(g, :D), [(D=0.4,)]))
    @test isempty(setdiff(findnodes(g, :B), [(B=0.5,),(B=0.6,)]))
    @test isempty(setdiff(findnodes(g, :E), [(E=1,),(E=2,),(E=3,)]))
    @test findnodes(g, :x) == NamedTuple[]
    @test g[] == []

    @test isempty(setdiff(g[:A], [1]))
    @test isempty(setdiff(g[:D], [0.4]))
    @test isempty(setdiff(g[:B], [0.5, 0.6]))
    @test isempty(setdiff(g[:E], [1, 2, 3]))
    B, E = g[(:B,:E)]
    @test isempty(setdiff(B, [0.5, 0.6]))
    @test isempty(setdiff(E, [3,2,1]))

    p = paths_through(g, (A=1,)=>(B=0.5,))
    @test length(p) == 3
    @test p == Set([1, 2, 3])
    @test length(paths_through(g, (A=1,)=>(B=0.5,)=>(E=1.,))) == 0
    @test length(paths_through(g, (A=1,))) == 6

    @test g[1] == get_prop(g, 1) == (D=0.4,)
    @test g[g[(A=1,)]] == (A=1,)
    @test isempty(setdiff(get_prop.(Ref(g), g[(A=1,)=>(D=0.4,)]), [(B=0.5,),(B=0.6,)]))
    @test isempty(setdiff(g[:B, (A=1,)=>(D=0.4,)], [0.5, 0.6]))
    @test isempty(setdiff(g[:E, (A=1,)=>(D=0.4,)=>(B=0.5,)], [1, 2, 3]))
    @test isempty(setdiff(g[:B, (A=1,), (D=0.4,)], [0.5, 0.6]))

    add_bulk!(g, (A=1,)=>(D=0.5,)=>(B=0.5,), (E=[5, 6],))
    @test isempty(setdiff(g[:E, (D=0.5,), (B=0.5,)], [5, 6]))
    @test isempty(setdiff(g[:E, (A=1,), (B=0.5,)], [1, 2, 3, 5, 6]))

    @test isempty(setdiff(g[:E, (D=0.5,)], [5,6]))
    @test isempty(setdiff(g[:E, Dict(:D=>v->v.D > 0.4), (A=1,)], [5,6]))
    @test g[Dict(:D=>v->v.D>0.4), (A=1,),(B=0.5,)] == [(E=5,),(E=6,)]
    @test isempty(setdiff(g[:E, Dict(:D=>v->v.D<0.5), (A=1,),(B=0.6,)], [1,2,3]))

    add_nodes!(g, (A=1,)=>(D=0.4,)=>(B=0.5,)=>(E=1,)=>(x=1,y=2))
    @test isempty(setdiff(g[(:x,:y), (A=1,),(E=1,)], ([1],[2])))
    @test g[:x, (A=1,)] == [1]
    @test isempty(setdiff(g[(:x,:y), (A=1,),(B=0.5,)],([1],[2])))

    add_nodes!(g, (A=1,)=>(D=0.4,)=>(B=0.5,)=>(E=2,)=>(x=2,y=0))
    conditions = Dict(:D=>v->v.D < 0.5, :E=>v->v.E > 1)
    @test g[:y, conditions, (A=1,)] == [0]

    @test @inferred paths_through(g, 1) == Set{eltype(g)}(1:6)
    @test @inferred paths_through(g, (A=1,)) == Set{eltype(g)}(1:8)
    @test @inferred paths_through(g, (A=1,)=>(D=0.4,)) == Set{eltype(g)}(1:6)
end


text:
using LaserTypes, Test
using LaserTypes: immutable_cache
using Unitful
using StaticArrays
using UnitfulAtomic

@testset "Values at origin" begin
    units = [:SI_unitful, :atomic_unitful, :SI, :atomic]
    x₀ = SVector{3}(0u"μm",0u"μm",0u"μm")
    t₀ = 0u"s"
    x = (x₀, auconvert.(x₀), ustrip.(x₀), ustrip.(auconvert.(x₀)))

    @testset "$unit" for (unit, xᵢ) in zip(units, x)
        s = setup_laser(GaussLaser, unit, profile=ConstantProfile())
        Ex, Ey, Ez = E(xᵢ, s)
        @test Ex ≈ immutable_cache(s, :E₀)
        @test iszero(Ey)
        @test iszero(Ez)

        Bx, By, Bz = B(xᵢ, s)
        @test iszero(Bx)
        @test By ≈ immutable_cache(s, :E₀) * immutable_cache(s, :inv_c)
        @test iszero(Bz)
    end
end

@testset "Orientation" begin
    s = setup_laser(GaussLaser, :SI, orientation=(:y,:x))
    x₀ = SVector{3}(0,0,0)
    t₀ = 0
    Ex, Ey, Ez = E(x₀, t₀, s)
    @test iszero(Ex)
    @test Ey ≈ immutable_cache(s, :E₀)
    @test iszero(Ez)

    Bx, By, Bz = B(x₀, t₀, s)
    @test iszero(Bx)
    @test iszero(By)
    @test Bz ≈ immutable_cache(s, :E₀) * immutable_cache(s, :inv_c)
end


text:
using Kronecker
using BenchmarkTools

SUITE = BenchmarkGroup()

SUITE["mul_vec"] = BenchmarkGroup(["vector", "multiplication"])
SUITE["mul_mat"] = BenchmarkGroup(["matrix", "multiplication"])
SUITE["ldiv_vec"] = BenchmarkGroup(["vector", "solve"])
SUITE["ldiv_mat"] = BenchmarkGroup(["matrix", "solve"])

###############################################################################

for k in keys(SUITE)
    SUITE[k]["all_square"] = BenchmarkGroup()
end

KProd = kronecker([randn(2,2) for i in 1:10]...)
KPow = kronecker(randn(2,2), 10)
KSum = kroneckersum([randn(2,2) for i in 1:10]...)
x = randn(2^10)

SUITE["mul_vec"]["all_square"]["KroneckerProduct"] = @benchmarkable ($KProd * $x)
SUITE["mul_vec"]["all_square"]["KroneckerPower"] = @benchmarkable ($KPow * $x)
SUITE["mul_vec"]["all_square"]["KroneckerSum"] = @benchmarkable ($KSum * $x)
SUITE["ldiv_vec"]["all_square"]["KroneckerProduct"] = @benchmarkable ($KProd \ $x)
SUITE["ldiv_vec"]["all_square"]["KroneckerPower"] = @benchmarkable ($KPow \ $x)
SUITE["ldiv_vec"]["all_square"]["KroneckerSum"] = @benchmarkable ($KSum \ $x)

###############################################################################

for k in keys(SUITE)
    SUITE[k]["all_rectangular"] = BenchmarkGroup()
end

KProd = kronecker([randn(3,2) for i in 1:10]...)
KPow = kronecker(randn(3,2), 10)
y = randn(3^10)

SUITE["mul_vec"]["all_rectangular"]["KroneckerProduct"] = @benchmarkable ($KProd * $x)
SUITE["mul_vec"]["all_rectangular"]["KroneckerPower"] = @benchmarkable ($KPow * $x)
SUITE["ldiv_vec"]["all_rectangular"]["KroneckerProduct"] = @benchmarkable ($KProd \ $y)
SUITE["ldiv_vec"]["all_rectangular"]["KroneckerPower"] = @benchmarkable ($KPow \ $y)

###############################################################################

for k in keys(SUITE)
    SUITE[k]["mixed_squareness"] = BenchmarkGroup()
end

KP = kronecker([randn(2,2) for i in 1:5]..., [randn(3,2) for i in 1:5]...)
x = randn(2^10)
y = randn((2^5)*(3^5))

SUITE["mul_vec"]["mixed_squareness"]["KroneckerProduct"] = @benchmarkable ($KP * $x)
SUITE["ldiv_vec"]["mixed_squareness"]["KroneckerProduct"] = @benchmarkable ($KP \ $y)


text:
filename() = "day12.txt"

function parse_line(line)
    instr = line[1]
    val = parse(Int, line[2:end])
    return (; instr, val)
end

read_data() = parse_line.(readlines(filename()))

function part1()
    pos = (x = 0, y = 0)
    rad = 0
    for d in read_data()
        if d.instr == 'N'
            pos = (x = pos.x, y = pos.y - d.val)
        elseif d.instr == 'S'
            pos = (x = pos.x, y = pos.y + d.val)
        elseif d.instr == 'E'
            pos = (x = pos.x + d.val, y = pos.y)
        elseif d.instr == 'W'
            pos = (x = pos.x - d.val, y = pos.y)
        elseif d.instr == 'F'
            pos = (x = pos.x + cos(rad) * d.val, y = pos.y + sin(rad) * d.val)
        elseif d.instr == 'L'
            rad -= deg2rad(d.val)
        elseif d.instr == 'R'
            rad += deg2rad(d.val)
        end
        @info "pos=$pos rad=$rad"
    end
    return abs(pos.x) + abs(pos.y)
end

function part2()
    pos = (x = 0, y = 0)
    way = (x = 10, y = -1)
    for d in read_data()
        dx, dy = way.x - pos.x, way.y - pos.y
        if d.instr == 'N'
            way = (x = way.x, y = way.y - d.val)
        elseif d.instr == 'S'
            way = (x = way.x, y = way.y + d.val)
        elseif d.instr == 'E'
            way = (x = way.x + d.val, y = way.y)
        elseif d.instr == 'W'
            way = (x = way.x - d.val, y = way.y)
        elseif d.instr == 'F'
            pos = (x = pos.x + d.val * dx, y = pos.y + d.val * dy)
            way = (x = way.x + d.val * dx, y = way.y + d.val * dy)
        elseif d.instr == 'L'
            r = sqrt(dx^2 + dy^2)
            θ = atan(dy, dx)
            ϕ = θ - deg2rad(d.val)
            way = (x = pos.x + r * cos(ϕ), y = pos.y + r * sin(ϕ))
        elseif d.instr == 'R'
            r = sqrt(dx^2 + dy^2)
            θ = atan(dy, dx)
            ϕ = θ + deg2rad(d.val)
            way = (x = pos.x + r * cos(ϕ), y = pos.y + r * sin(ϕ))
        else
            error("wat")
        end
    end
    return abs(pos.x) + abs(pos.y)
end

# Animation

using Plots
using FileIO

function part2_animated()
    pos = (x = 0, y = 0)
    way = (x = 10, y = 1)
    anim = Animation()
    # img = FileIO.load("/Users/tomkwong/Downloads/iu.jpeg")
    # plot(img)
    # frame(anim)
    @animate for (i,d) in enumerate(read_data())
        dx, dy = way.x - pos.x, way.y - pos.y
        if d.instr == 'N'
            way = (x = way.x, y = way.y + d.val)
        elseif d.instr == 'S'
            way = (x = way.x, y = way.y - d.val)
        elseif d.instr == 'E'
            way = (x = way.x + d.val, y = way.y)
        elseif d.instr == 'W'
            way = (x = way.x - d.val, y = way.y)
        elseif d.instr == 'F'
            pos = (x = pos.x + d.val * dx, y = pos.y + d.val * dy)
            way = (x = way.x + d.val * dx, y = way.y + d.val * dy)
        elseif d.instr == 'L'
            r = sqrt(dx^2 + dy^2)
            θ = atan(dy, dx)
            ϕ = θ + deg2rad(d.val)
            way = (x = pos.x + r * cos(ϕ), y = pos.y + r * sin(ϕ))
        elseif d.instr == 'R'
            r = sqrt(dx^2 + dy^2)
            θ = atan(dy, dx)
            ϕ = θ - deg2rad(d.val)
            way = (x = pos.x + r * cos(ϕ), y = pos.y + r * sin(ϕ))
        else
            error("wat")
        end
        if i % 2 == 0
            s, adx, ady = 3, abs(dx), abs(dy)
            ms = 15
            xlims = (min(pos.x, way.x) - s * adx, max(pos.x, way.x) + s * adx)
            ylims = (min(pos.y, way.y) - s * ady, max(pos.y, way.y) + s * ady)
            scatter([pos.x], [pos.y], xlims = xlims, ylims = ylims, markersize = ms,
                legend = nothing, ticks = nothing, showaxis = false,
                background = "black", markershape = :octagon,
                title = "AoC Day 12 - Iteration $i")
            scatter!([way.x], [way.y],
                xlims = xlims, ylims = ylims, markersize = ms + 5,
                markershape = :star)
            frame(anim)
        end
    end
    gif(anim, "day12_anim.gif", fps = 15)
end


text:
module LinkedLists

using Base
import ..AbstractDataStructure

export LinkedList, ins!, del!, at, search, len, ordered_data_structure_p

struct _EOL end

EOL = _EOL()

mutable struct LinkedList{V} <: AbstractDataStructure.DataStructure{Int, V}
    data::Union{V, _EOL}
    tail::Union{LinkedList{V}, Nothing}
    function LinkedList{V}(data::Union{V, _EOL}, tail::Union{LinkedList{V}, Nothing}) where {V}
        @assert((data === EOL && tail === nothing) || (data !== EOL && tail !== nothing))
        new{V}(data, tail)
    end
end

LinkedList{V}() where {V} = LinkedList{V}(EOL, nothing)

function LinkedList{V}(vec::Array{Any, 1}) where {V}
    if length(vec) > 0
        vec::Vector{V} = convert(Array{V}, vec)
        return LinkedList{V}(vec)
    end
    LinkedList{V}(EOL, nothing)
end

function LinkedList{V}(vec::Vector{V}) where {V}
    l = LinkedList{V}()
    for e in reverse(vec)
        l = LinkedList{V}(e, l)
    end
    l
end

AbstractDataStructure.ordered_data_structure_p(::LinkedList) = false

function hopn(l::LinkedList, n::Int)::LinkedList
    @assert(n ≥ 0)
    i = 0
    while i < n
        if l.data === EOL; throw(KeyError(n)) end
        l = l.tail
        i += 1
    end
    if l === nothing; throw(KeyError(n)) end
    l
end

Base.iterate(l::LinkedList) = l.data === EOL ? nothing : (l.data, l.tail)
Base.iterate(::LinkedList, state::LinkedList) = state.data === EOL ? nothing : (state.data, state.tail)
Base.iterate(::LinkedList, ::Nothing) = nothing

"""
ins!(ds::DataStructure, val::Value) -> Ref
"""
function AbstractDataStructure.ins!(l::LinkedList{V}, val::V)::Int where {V}
    l.tail = LinkedList{V}(l.data, l.tail)
    l.data = val
    1
end

"""
del!(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function AbstractDataStructure.del!(l::LinkedList{V}, pos::Int)::V where {V}
    @assert(pos > 0)
    to_delete = hopn(l, pos-1)
    if to_delete.data === EOL; throw(KeyError(pos)) end
    removed_data = to_delete.data
    to_delete.data = to_delete.tail.data
    to_delete.tail = to_delete.tail.tail
    removed_data
end

"""
at(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function AbstractDataStructure.at(l::LinkedList{V}, pos::Int)::V where {V}
    @assert(pos > 0)
    node = hopn(l, pos-1)
    if node.data === EOL; throw(KeyError(pos)) end
    node.data
end

"""
search(ds::DataStructure, val::Value) -> Vector{Ref}
"""
function AbstractDataStructure.search(ds::LinkedList{V}, val::V)::Vector{Int} where {V}
    return [i for (i, e) in enumerate(ds) if e == val]
end

"""
len(ds::DataStructure) -> Int
"""
function AbstractDataStructure.len(ds::LinkedList)::Int
    i = 0
    while ds.data != EOL; ds = ds.tail; i += 1 end
    i
end

_iobuf = IOBuffer()
printstyled(IOContext(_iobuf, :color => true), "[ ", color=:green)
_LIST_OPEN_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " ]", color=:green)
_LIST_CLOSE_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " → ", color=:light_green)
_LIST_SEP_COLOR = String(take!(_iobuf))
_iobuf = nothing

function Base.show(io::IO, x::LinkedList)
    if get(io, :color, false)
        open = _LIST_OPEN_COLOR
        close = _LIST_CLOSE_COLOR
        sep = _LIST_SEP_COLOR
    else
        open = "[ "
        close = " ]"
        sep = " → "
    end
    print(io, "LinkedList: " * open)
    join(io, x, sep)
    print(io, close)
end

end  # module LinkedLists


text:
using DataFrames

@testset "Synthesis" begin

# get data for use in tests
data = GRASS.SolarData()
wavt = data.wav[(:c, :mu10)][:,1]
bist = data.bis[(:c, :mu10)][:,1]
dept = data.dep[(:c, :mu10)][:,1]
widt = data.wid[(:c, :mu10)][:,1]

# set trimming parameters
top = NaN
dep = 0.75

@testset "Testing function definitions" begin
    @test isdefined(GRASS, :trim_bisector_chop!)
    @test isdefined(GRASS, :line_profile!)
    @test isdefined(GRASS, :line_from_bis!)
    @test isdefined(GRASS, :synthesize_spectra)
end

@testset "Testing line trimming" begin
    # copy the data
    wavt1 = copy(wavt)
    bist1 = copy(bist)
    dept1 = copy(dept)
    widt1 = copy(widt)

    # do the trimming
    GRASS.trim_bisector_chop!(dep, wavt1, bist1, dept1, widt1, top=top)

    # test that it wasn't changed
    @test wavt1 !== wavt
    @test bist1 !== bist
    @test dept1 !== dept
    @test widt1 !== widt

    # top and bottom
    @test bist1[1] == 1.0 - dep
    @test bist1[end] == 1.0
    @test wavt1[1] == wavt[findfirst(bist .>= 1.0 - dep)]
end

@testset "Testing line synthesis" begin
    # copy the data
    wavt1 = copy(wavt)
    bist1 = copy(bist)
    dept1 = copy(dept)
    widt1 = copy(widt)

    # do the trimming
    GRASS.trim_bisector_chop!(dep, wavt1, bist1, dept1, widt1, top=top)

    # allocate memory for line
    λs = range(5434.0, 5435.0, length=1000)
    lwavgrid = zeros(100)
    rwavgrid = zeros(100)
    allwavs = zeros(200)
    allints = zeros(200)
    prof = ones(length(λs))

    # synthesize line
    mid = 5434.5
    GRASS.line_from_bis!(mid, λs, prof, wavt1, dept1, widt1, lwavgrid, rwavgrid, allwavs, allints)

    # test it
    @test !all(prof .== 1.0)
    @test prof[1] == 1.0
    @test prof[end] == 1.0
    @test minimum(prof) == 1.0 - dep
    @test maximum(prof) == 1.0
end

@testset "Testing bisector in/out" begin
    # read in data
    wavt1 = copy(wavt)
    bist1 = copy(bist)
    dept1 = copy(dept)
    widt1 = copy(widt)

    # do the trimming
    GRASS.trim_bisector_chop!(dep, wavt1, bist1, dept1, widt1, top=top)

    # allocate memory for line
    λs = range(5434.0, 5435.0, length=1000)
    lwavgrid = zeros(100)
    rwavgrid = zeros(100)
    allwavs = zeros(200)
    allints = zeros(200)
    prof = ones(length(λs))

    # synthesize line
    mid = 5434.5
    GRASS.line_from_bis!(mid, λs, prof, wavt1, dept1, widt1, lwavgrid, rwavgrid, allwavs, allints)

    # measure the bisector of synth line
    wavo, biso = GRASS.measure_bisector(λs, prof, interpolate=false)
    wavo .-= mid

    # test that output bisector matches input bisector
    @test all(isapprox.(wavo[2:99], wavt1[2:99], atol=1e-3))
end

@testset "Testing disk-integrated spectrum synthesis" begin
    # set params
    N = 132
    Nt = 5
    res = 7e5

    # simulate spectra
    spec = SpecParams(lines=[5434.5], depths=[dep], resolution=res)
    disk = DiskParams(N=N, Nt=Nt)
    wavs, flux = synthesize_spectra(spec, disk)

    # ensure spectra have correct properties
    @test size(flux,2) == 5
    @test maximum(flux[:,1]) == 1.0
    # @test isapprox(minimum(flux[:,1]), 1.0 - dep, atol=1e-1) # TODO: fix?

    # get velocities
    v_grid, ccf1 = calc_ccf(wavs, flux, spec, normalize=true)
    rvs, sigs = calc_rvs_from_ccf(v_grid, ccf1)

    # test the velocities
    @test calc_rms(rvs) < 1.0
end

end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using Test
using Assets

@testset "testall" begin

    pkgpath = joinpath("foo", "bar", "MyPkg")
    filepath = joinpath(pkgpath, "src", "code.jl")
    @test Assets.getpkgfolder(filepath) == pkgpath

    @test Assets.is_standard_julia_session()

    pkgpath = joinpath("foo", "bar", "MyPkg")
    @test Assets.assetpath(pkgpath) == joinpath(pkgpath, "assets")

end

max_issues_repo_name: nielsls/Assets

max_forks_repo_name: nielsls/Assets


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

using Markdown
using InteractiveUtils

# ╔═╡ d18893e3-89ec-46fc-9265-0daac744caea
using Manopt, Manifolds, Colors, PlutoUI

# ╔═╡ 3b55c91e-bda9-11ec-0399-b77442d8ff45
md"""
# Illustration of the Gradient of a Second Order Difference

This example explains how to compute the gradient of the second order
difference mid point model using `adjoint_Jacobi_field`

This example also illustrates the `PowerManifold` manifold as well as `ArmijoLinesearch`.
"""

# ╔═╡ b2a115cc-5d4e-4052-afc1-2972ed493958
md"We define some colors from [Paul Tol](https://personal.sron.nl/~pault/)"

# ╔═╡ c19c8ede-ea81-4b10-85c2-f73051d8090a
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
    TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
    TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
    TolVibrantTeal = RGBA{Float64}(colorant"#009988")
end;

# ╔═╡ b5ba2630-e247-46b9-b306-faf8a531f22a
begin
    T = [0:0.1:1.0...]
    #render asy yes/no. If not, images included w/ markdown are assumed to be prerendered
    render_asy = false
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # remove file to get files folder
    image_prefix = localpath * "/second_order_difference"
    @info image_prefix
end;

# ╔═╡ 2a51bd75-c47e-452c-bfea-f6d7cfcd3e45
md"""
Assume we have two points ``p,q`` on the equator of the
[Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) ``\mathcal M = \mathbb S^2``
and a point ``r`` near the north pole
"""

# ╔═╡ 3c8a6104-d380-4a32-907c-aaf205ea0bdd
begin
    M = Sphere(2)
    p = [1.0, 0.0, 0.0]
    q = [0.0, 1.0, 0.0]
    c = mid_point(M, p, q)
    r = shortest_geodesic(M, [0.0, 0.0, 1.0], c, 0.1)
end;

# ╔═╡ 9cfecd72-0c1a-4432-9628-54b71db9c22d
md"""
Now the second order absolute difference can be stated as [^BacakBergmannSteidlWeinmann2016]

```math
d_2(p_1,p:2,p_3) := \min_{c ∈ \mathcal C_{p_1,p_3}} d_{\mathcal M}(c,p_2),\qquad p_1,p_2,p_3∈\mathcal M,
```

where ``\mathcal C_{p,q}`` is the set of all mid points ``g(\frac{1}{2};p,q)``, between `p`and `q`, i.e. where ``g``
is a (not necessarily minimizing) geodesic connecting both.

For illustration we further define the point opposite of the mid point `c` defined above
"""

# ╔═╡ a9596c67-1cfb-4d94-a7c8-71bf36fe06bb
c2 = -c;

# ╔═╡ 2dfaebe0-e706-4eb6-b661-aed9d8ae19c1
md"To illustrate the second order difference let‘s look at the geodesic connecting ``r`` and the mid point ``c``"

# ╔═╡ 1c484736-bb94-4f64-bcfd-ad7ed52ade39
geoPts_rc = shortest_geodesic(M, r, c, T);

# ╔═╡ 80aa3913-7639-4752-a762-58c17a7e1c47
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/SecondOrderData.asy";
        curves=[geoPts_rc],
        points=[[p, r, q], [c, c2]],
        colors=Dict(:curves => [TolVibrantTeal], :points => [black, TolVibrantBlue]),
        dot_size=3.5,
        line_width=0.75,
        camera_position=(1.2, 1.0, 0.5), #src
    )
    render_asymptote(image_prefix * "/SecondOrderData.asy"; render=2)
end;

# ╔═╡ 717ba868-7da5-4ce1-9163-2b20f8cdc253
PlutoUI.LocalResource(image_prefix * "/SecondOrderData.png")

# ╔═╡ dcf6bbeb-5b09-4689-8bb9-e35f089cff45
md" **Figure.** _The three poins ``p, q``, and ``r`` (neart north pole, all black), which is connected by a geodesic (teal) to the mid point of the first two, ``c``, (blue). The length of this geodesic is the cost of the second order total variation._"

# ╔═╡ 07a44b4c-876c-4fff-ad72-2baedc06d090
md"""
Since we moved ``r`` 10% along the geodesic from the north pole to ``c``, the distance
to ``c`` is ``\frac{9\pi}{20}\approx 1.4137``, and this is also what the second order total variation cost, see [`costTV2`](https://manoptjl.org/stable/functions/costs.html#Manopt.costTV2), yields.
"""

# ╔═╡ ebab946c-ba27-422d-a191-56294e57177b
costTV2(M, (p, r, q))

# ╔═╡ db26c554-31ac-4d62-93a1-9440361fa821
md"""
But also its gradient can be
easily computed since it is just a distance with respect to ``r`` and a
concatenation of a geodesic, where the start or end point is the argument,
respectively, with a distance.
Hence the adjoint differentials
[`adjoint_differential_geodesic_startpoint`](https://manoptjl.org/stable/functions/adjointdifferentials.html#Manopt.adjoint_differential_geodesic_startpoint-Tuple{AbstractManifold,%20Any,%20Any,%20Any,%20Any}) and [`adjoint_differential_geodesic_endpoint`](https://manoptjl.org/stable/functions/adjointdifferentials.html#Manopt.adjoint_differential_geodesic_endpoint-Tuple{AbstractManifold,%20Any,%20Any,%20Any,%20Any}) can be employed.
The gradient is also directly implemented, see [`grad_TV2`](https://manoptjl.org/stable/functions/gradients.html#Manopt.grad_TV2).
we obtain
"""

# ╔═╡ b685cd13-8b06-4d6d-a001-576f68951bd9
(Xp, Xr, Xq) = grad_TV2(M, (p, r, q))

# ╔═╡ dddd7edf-47b7-4b6d-8b20-e49b4db211a9
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/SecondOrderGradient.asy";
        points=[[p, r, q], [c, c2]],
        tangent_vectors=[Tuple.([[p, -Xp], [r, -Xr], [q, -Xq]])],
        colors=Dict(:tvectors => [TolVibrantCyan], :points => [black, TolVibrantBlue]),
        dot_size=3.5,
        line_width=0.75,
        camera_position=(1.2, 1.0, 0.5),
    )
    render_asymptote(image_prefix * "/SecondOrderGradient.asy"; render=2)
end;

# ╔═╡ 18bc72b8-fef9-4890-bdb8-dc477e0109db
PlutoUI.LocalResource(image_prefix * "/SecondOrderGradient.png")

# ╔═╡ cdf26e32-2383-4ac3-b814-ef7c98860b55
md"""
**Figure.** _The negative gradient of the second order difference cost indicates the movement of the three points in order to reduce their cost._
"""

# ╔═╡ 77149fb5-f3cb-438f-843c-9221ae5c08a9
md"If we now perform a gradient step with constant step size 1, we obtain the three points"

# ╔═╡ e53622ca-3e08-4e16-a1a3-6536a46c8a9f
pn, rn, qn = exp.(Ref(M), [p, r, q], [-Xp, -Xr, -Xq])

# ╔═╡ 6b50fb13-45db-4cbb-a1e0-d62187c0b2d3
md"as well we the new mid point"

# ╔═╡ 0bbd56ce-9041-4974-8525-a0b5bd03e9fe
cn = mid_point(M, pn, qn)

# ╔═╡ 2eee7e24-20aa-4ae4-a336-b2eec90ce24c
md"Let‘s also again consider the geodesic connecting the new point ``r_n`` and the new mid point ``c_n`` as well as the gradient"

# ╔═╡ ce758b03-0177-45c2-b895-e64d133ca4af
begin
    geoPts_rncn = shortest_geodesic(M, rn, cn, T)
    (Xpn, Xrn, Xqn) = grad_TV2(M, (pn, rn, qn))
end;

# ╔═╡ 77e2409a-0168-4715-8770-23272e469f87
md"The new configuration of the three points looks as follows"

# ╔═╡ d90730e1-7111-42d9-8688-b4ff9cca9c6e
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/SecondOrderMin1.asy";
        points=[[p, r, q], [c, c2, cn], [pn, rn, qn]],
        curves=[geoPts_rncn],
        tangent_vectors=[
            Tuple.([[p, -Xp], [r, -Xr], [q, -Xq]]),
            Tuple.([[pn, -Xpn], [rn, -Xrn], [qn, -Xqn]]),
        ],
        colors=Dict(
            :tvectors => [TolVibrantCyan, TolVibrantMagenta],
            :points => [black, TolVibrantBlue, TolVibrantOrange],
            :curves => [TolVibrantTeal],
        ),
        dot_size=3.5,
        line_width=0.75,
        camera_position=(1.2, 1.0, 0.5),
    )
    render_asymptote(image_prefix * "/SecondOrderMin1.asy"; render=2)
end;

# ╔═╡ 533f324b-40ca-4b74-968e-e9f1ce4262eb
PlutoUI.LocalResource(image_prefix * "/SecondOrderMin1.png")

# ╔═╡ a959be34-c0fb-4dc3-b46d-1702f616a8bb
md"""
**Figure.** *The new situation of ``p_n, q_n``, and ``r_n`` (orange) and the miod point of the first two, ``c_n`` (blue), which is again connected by a geodesic to ``r_n``. Note that this geodesic is shorter, but also that ``c`` and ``r`` switched places*. The new gradient (magenta) is only slightly reduced in magnitude.*
"""

# ╔═╡ 5b4aa1ae-7ab1-4efc-a796-1fa4e771e2b2
md"""
One can see, that this step slightly “overshoots”, i.e. ``r`` is now even below ``c``.
and the cost function is still at
"""

# ╔═╡ ea3ceaad-d032-42fd-84ca-a08b2fdd6db0
costTV2(M, (pn, rn, qn))

# ╔═╡ 96600890-af72-4842-973e-438cb930f57b
md"""
But we can also search for the best step size using [`linesearch_backtrack`](https://manoptjl.org/stable/plans/index.html#Manopt.linesearch_backtrack-Union{Tuple{T},%20Tuple{TF},%20Tuple{AbstractManifold,%20TF,%20Any,%20T,%20Any,%20Any,%20Any},%20Tuple{AbstractManifold,%20TF,%20Any,%20T,%20Any,%20Any,%20Any,%20AbstractRetractionMethod},%20Tuple{AbstractManifold,%20TF,%20Any,%20T,%20Any,%20Any,%20Any,%20AbstractRetractionMethod,%20T},%20Tuple{AbstractManifold,%20TF,%20Any,%20T,%20Any,%20Any,%20Any,%20AbstractRetractionMethod,%20T,%20Any}}%20where%20{TF,%20T})
on the `PowerManifold` manifold ``\mathcal N = \mathcal M^3 = (\mathbb S^2)^3``
"""

# ╔═╡ 4b662610-0812-4ae2-809c-52f7a8979398
begin
    x = [p, r, q]
    N = PowerManifold(M, NestedPowerRepresentation(), 3)
    s = linesearch_backtrack(N, x -> costTV2(N, x), x, grad_TV2(N, x), 1.0, 0.957, 0.999)
end

# ╔═╡ 17ccdc00-35ad-4ec5-a480-6f6527e38834
md"and we obtain the new points"

# ╔═╡ db8916e6-0ee7-4e37-a879-5f35aa11d38d
begin
    pm, rm, qm = exp.(Ref(M), [p, r, q], s * [-Xp, -Xr, -Xq])
    cm = mid_point(M, pm, qm)
    geoPts_pmqm = shortest_geodesic(M, pm, qm, T)
end;

# ╔═╡ 78d44ce1-9d7b-4cb1-b9cb-2d402b1dfabb
md"we obtain"

# ╔═╡ fdf9aa5e-a25a-42f3-a74c-cef5877f2a1d
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/SecondOrderMin2.asy";
        points=[[p, r, q], [c, c2, cm], [pm, rm, qm]],
        curves=[geoPts_pmqm],
        tangent_vectors=[Tuple.([[p, -Xp], [r, -Xr], [q, -Xq]])],
        colors=Dict(
            :tvectors => [TolVibrantCyan],
            :points => [black, TolVibrantBlue, TolVibrantOrange],
            :curves => [TolVibrantTeal],
        ),
        dot_size=3.5,
        line_width=0.75,
        camera_position=(1.2, 1.0, 0.5),
    )
    render_asymptote(image_prefix * "/SecondOrderMin2.asy"; render=2)
end;

# ╔═╡ b0a75ff9-9ff7-4db6-ad82-76f183663dc0
PlutoUI.LocalResource(image_prefix * "/SecondOrderMin2.png")

# ╔═╡ 9c37f43b-dc3e-45fb-a99f-756c83f2f658
md"""
**Figure.** *For the best step size found by line search, ``r_m`` and ``c_m`` nearly agree, i.e. ``r_m`` lies on the geodesic between ``p_m`` and ``q_m`` as the geodesic drawn here indicates.*
"""

# ╔═╡ 428d1422-7521-42fd-993d-e617115edf14
md"Here, the cost function yields"

# ╔═╡ 25804825-c9f6-4cfa-b209-918ad19d2ae1
costTV2(M, (pm, rm, qm))

# ╔═╡ 718f744e-e6af-4158-9418-7d2f4b757903
md"which is much closer to zero, as one can also see, since the new center ``c_m`` and ``r_m`` are quite close."

# ╔═╡ 962936da-c154-444f-8d48-ca1176dee1a1
md"""
## Literature

[^BacakBergmannSteidlWeinmann2016]:
    > Bačák, M., Bergmann, R., Steidl, G. and Weinmann, A.: _A second order nonsmooth
    > variational model for restoring manifold-valued images_,
    > SIAM Journal on Scientific Computations, Volume 38, Number 1, pp. A567–597,
    > doi: [10.1137/15M101988X](https://doi.org/10.1137/15M101988X),
    > arXiv: [1506.02409](https://arxiv.org/abs/1506.02409)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Colors = "~0.12.8"
Manifolds = "~0.7.8"
Manopt = "~0.3.23"
PlutoUI = "~0.7.38"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

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

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

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

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "c933ce606f6535a7c7b98e1d86d5d1014f730596"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.7"

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

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

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

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

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

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

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

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

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

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

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

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

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

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

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

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

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

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

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

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

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

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

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

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

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "a970d55c2ad8084ca317a4658ba6ce99b7523571"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.12"

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

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "cc5e34ecd061bf9bedd62a8b8a28abe96c3b7c6e"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.8"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "b8518cbf0f06f1375c95479759b266bd75d89ad9"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.13"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "a8a292a422b0a6924f46d4e022fc8263721120c5"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.23"

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

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

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

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

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

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

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

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

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "bfe14f127f3e7def02a6c2b1940b39d0dabaa3ef"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.26.3"

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

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

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

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

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

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "8d7530a38dbd2c397be7ddd01a424e4f411dcc41"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.2"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

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

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

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

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

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

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

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

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

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

# ╔═╡ Cell order:
# ╟─3b55c91e-bda9-11ec-0399-b77442d8ff45
# ╠═d18893e3-89ec-46fc-9265-0daac744caea
# ╟─b2a115cc-5d4e-4052-afc1-2972ed493958
# ╠═c19c8ede-ea81-4b10-85c2-f73051d8090a
# ╠═b5ba2630-e247-46b9-b306-faf8a531f22a
# ╟─2a51bd75-c47e-452c-bfea-f6d7cfcd3e45
# ╠═3c8a6104-d380-4a32-907c-aaf205ea0bdd
# ╟─9cfecd72-0c1a-4432-9628-54b71db9c22d
# ╠═a9596c67-1cfb-4d94-a7c8-71bf36fe06bb
# ╟─2dfaebe0-e706-4eb6-b661-aed9d8ae19c1
# ╠═1c484736-bb94-4f64-bcfd-ad7ed52ade39
# ╠═80aa3913-7639-4752-a762-58c17a7e1c47
# ╟─717ba868-7da5-4ce1-9163-2b20f8cdc253
# ╟─dcf6bbeb-5b09-4689-8bb9-e35f089cff45
# ╟─07a44b4c-876c-4fff-ad72-2baedc06d090
# ╠═ebab946c-ba27-422d-a191-56294e57177b
# ╟─db26c554-31ac-4d62-93a1-9440361fa821
# ╠═b685cd13-8b06-4d6d-a001-576f68951bd9
# ╠═dddd7edf-47b7-4b6d-8b20-e49b4db211a9
# ╟─18bc72b8-fef9-4890-bdb8-dc477e0109db
# ╟─cdf26e32-2383-4ac3-b814-ef7c98860b55
# ╟─77149fb5-f3cb-438f-843c-9221ae5c08a9
# ╠═e53622ca-3e08-4e16-a1a3-6536a46c8a9f
# ╟─6b50fb13-45db-4cbb-a1e0-d62187c0b2d3
# ╟─0bbd56ce-9041-4974-8525-a0b5bd03e9fe
# ╠═2eee7e24-20aa-4ae4-a336-b2eec90ce24c
# ╠═ce758b03-0177-45c2-b895-e64d133ca4af
# ╟─77e2409a-0168-4715-8770-23272e469f87
# ╠═d90730e1-7111-42d9-8688-b4ff9cca9c6e
# ╠═533f324b-40ca-4b74-968e-e9f1ce4262eb
# ╟─a959be34-c0fb-4dc3-b46d-1702f616a8bb
# ╟─5b4aa1ae-7ab1-4efc-a796-1fa4e771e2b2
# ╠═ea3ceaad-d032-42fd-84ca-a08b2fdd6db0
# ╟─96600890-af72-4842-973e-438cb930f57b
# ╠═4b662610-0812-4ae2-809c-52f7a8979398
# ╟─17ccdc00-35ad-4ec5-a480-6f6527e38834
# ╠═db8916e6-0ee7-4e37-a879-5f35aa11d38d
# ╟─78d44ce1-9d7b-4cb1-b9cb-2d402b1dfabb
# ╠═fdf9aa5e-a25a-42f3-a74c-cef5877f2a1d
# ╠═b0a75ff9-9ff7-4db6-ad82-76f183663dc0
# ╟─9c37f43b-dc3e-45fb-a99f-756c83f2f658
# ╟─428d1422-7521-42fd-993d-e617115edf14
# ╠═25804825-c9f6-4cfa-b209-918ad19d2ae1
# ╟─718f744e-e6af-4158-9418-7d2f4b757903
# ╟─962936da-c154-444f-8d48-ca1176dee1a1
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#RSDeltaSigma: simulateMS() algorithm
#-------------------------------------------------------------------------------

#==simulateMS()
===============================================================================#
"""`simulateMS(v, M=16, mtf=zpk(1,0,1,1), d=0, dw=[1..], sx0=[0..])`

Simulate the vector-based Mismatch-Shaping for a multi-element DAC:
 - simresult = (sv, sx, sigma_se, max_sx, max_sy)

# Inputs:
 - `v`: A vector of the digital input values. v is in (-M:2:M) if dw=[1..]
   otherwise v is in [-sum(dw),sum(dw)].
 - `M`: The number of elements.
 - `mtf`: The mismatch-shaping transfer function, given in zero-pole form.
 - `d`: Dither uniformly distributed in [-d,d] is added to sy.
 - `dw`: A vector of dac element weights
 - `sx0`: A matrix whose columns are the initial states of the ESL.

# Outputs:
 - `sv`: An MxN matrix whose columns are the selection vectors.
 - `sx`: An orderxM matrix containing the final state of the ESL.
 - `sigma_se`: The rms value of the selection error.
 - `max_sx`: The maximum absolute value of the state for all modulators.
 - `max_sy`: The maximum absolute value of the input to the VQ.

# Notes:
For a brief description of the theory of mismatch-shaping DACs, see
R. Schreier and B. Zhang "Noise-shaped multibit D/A convertor employing
unit elements", Electronics Letters, vol. 31, no. 20, pp. 1712-1713,
Sept. 28 1995.
"""
function simulateMS(v; M::Int=16, mtf=[], d::Float64=0.0, dw::Vector=[], sx0=[])
	local mtf_z; local mtf_p; local mtf_k

	if isa(mtf, Array) && isempty(mtf)
		mtf_z = [1]; mtf_p = [0]
	elseif isa(mtf, ZPKData)
		(mtf_z, mtf_p, mtf_k) = _zpkdata(mtf)
	else
		_type = typeof(mtf)
		throw("mtf type not supported: $_type.")
	end

	order = length(mtf_p)
	if isempty(sx0)
		sx0 = zeros(order,M)
	end
	if isempty(dw)
		dw = ones(M)
	end

	#B/A = MTF-1
	num = poly(mtf_z)
	den = poly(mtf_p)
	A = real.(-den[2:order+1])
	B = real.(num[2:order+1]+A)
	A = A'; B = B' #Need row vectors
	N = length(v)
	sv = zeros(M,N)

	sx = sx0
	max_sx = maximum(maximum(abs.(sx)))
	max_sy = 0
	sum_se2 = 0

	for i in 1:N
		#Compute the sy vector.
		sy = B*sx
		#Normalize sy for a minimum value of zero.
		sy = sy .- minimum(sy)
		syd = (sy + d*(2*rand(1,M) .- 1))[:]
		#Pick the elements that have the largest desired_usage (sy) values.
		sv[:,i] = selectElement(v[i], syd, dw)
		#Compute the selection error.
		se = sv[:,[i]]' - sy
		#Compute the new sx matrix
		sxn = A*sx + se
		sx = [ sxn; sx[1:order-1,:]]
		#Keep track of some statistics.
		sum_se2 = sum_se2 + sum((se .- mean(se)) .^ 2)
		max_sx = maximum([max_sx abs.(sxn)])
		max_sy = maximum([max_sy abs.(sy)])
	end
	sigma_se = sqrt(sum_se2/(M*N))

	simresult = (sv=sv, sx=sx, sigma_se=sigma_se, max_sx=max_sx, max_sy=max_sy)
	return simresult
end

@warn("Implements .m version of simulateMS(). C code might be more efficient/validated.")

#Last line


text:
# # Примеры
#
# На этой странице мы делаем `include` примеров,
# которые лежат в корневом каталоге пакета в папке `"examples/"`.
#
# [Literate.jl](https://github.com/fredrikekre/Literate.jl) включает
# текст этих примеров и в md файл.
#
#src `#src` это комментарии даже для Literate.jl
#src `#md` Строки Literate.jl включит только в md файлы.
#
#src Перед всматриванием в следующий блок, посмотри `docs/make.jl`.
#src Там, скормив этот скрипт Literate.jl мы получим страницу
#src `docexample_adjoint.md`. Именно её будет обрабатывать Documenter.jl
#
#md # ```@contents
#md # Pages = [ "docexample_adjoint.md" ]
#md # Depth = 5
#md # ```

include("examples/example_adjoint.jl")

# ### Визуализация

# Мы не можем строить большие анимации на стороне `Travis`-a.
# Если мы на CI, то будем рисовать только 10 кадров.
# Если мы генерируем документацию локально, то рисуем 80 кадров.
isTravis = in("Travis", keys(ENV))
ftw = isTravis ? range(1, stop = M+1, length=9) : [1; 2:div(M+1, 80):M+1];

# На отрисовку, решение сопряженной задачи передадим в инвертированном времени.
# Передадим свою подпись к графиками с помощью keyword `label="\\psi"`,
# не забыл про экранировку спец символа
make_gif(reverse(ψ, dims=2), XX, Tₘ;
         label="\\psi",
         name="solution_adjoint_ex1.gif",
         frames_to_write = ftw)

# Результат должен быть около нулевой, ведь в качестве текущего приближения `q` мы взяли искомое,
# а при нем — градиент должен обнуляться.

include("examples/example_adjoint_nonuniform.jl")

# Результат должен быть около нулевой, ведь в качестве текущего приближения `q` мы взяли искомое,
# а при нем — градиент должен обнуляться.

# ### [Контроль параметра аппроксимации дельта-функции](@id empiric-control)
#
# Проведем контроль выбора априорного параметра в аппроксимации дельта-функции.
# На точных данных, эта неоднородность должна быть нулевой почти везде исходя
# из следующих соображений:
# * сопряженная задача — ретроспективная,
# * её начальные условие — нулевые,
# * если из уравнения вычеркнуть эту неоднородность, то решением с такими
#   ГУ и НУ будет тривиальное нулевое решение.
# * исходя из [формулировки градиента](@id function-gradient),
# * нулевое решение сопряженной задачи даст нам нулевой градиент,
#   а значит мы нашли решение обратной задачи.
#
# Выбор этого параметра **значительно** влияет на вычисление градиента, а
# значит и на решение, к которому сойдется итерационный процесс.
using  NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity;
using LaTeXStrings, Plots;
Uₙₘ = u[2:N, :];
X = XX[2:N, :];
H = [  - heterogeneity(n, m, X[:, m], N, Tₘ, M, Uₙₘ, f1_data, f2_data, w) for n in 1:N-1, m in 1:M+1]
heatmap(H', title=L"-2 δ( x - f_1(t))(u(f1(t), t) - f2(t) ")

# ### Визуализация
ftw = isTravis ? range(1, stop = M+1, length=9) : [1; 2:div(M+1, 80):M+1];

make_gif(reverse(ψ, dims=2), XX, Tₘ;
         label="\\psi",
         name="solution_adjoint_ex3.gif",
         frames_to_write = ftw)


text:
using NiLang
using Test

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

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

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

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

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

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

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

@testset "autodiff" begin
    include("autodiff/autodiff.jl")
end

@testset "stdlib" begin
    include("stdlib/stdlib.jl")
end


text:
using ReverseSearch
using FactCheck

include("../examples/complete_graph.jl")

context("reverse search graph traversal on the complete graph on 1:n") do
    for n in 1:10
        facts("for n=$n, every vertex is visited exactly once") do
            @fact reverse_search_vertices_for_complete_graph(n) --> collect(1:n)
        end
    end
end


text:
using Luxor, Random, Test

function test_tangents(fname)
    Drawing(500, 500, fname)
    background("white")
    origin()
    setline(0.5)
    sethue("black")
    for i in 1:30
        circlecenter = polar(rand(100:190), 2π * i/30)
        circleradius = rand(4:10)
        point = polar(rand(170:240), 2π * i/20)
        circle(circlecenter, 2, :fill)
        circle(point, 2, :fill)
        circle(circlecenter, circleradius, :fill)
        pt1, pt2 = pointcircletangent(point, circlecenter, circleradius)
        if pt1 != pt2
            line(point, pt1, :stroke)
            line(point, pt2, :stroke)
        end
    end
    @test finish() == true
end

fname = "point-circle-tangent.png"
test_tangents(fname)
println("...finished test: output in $(fname)")


text:
function e003(N)
    return largest_factor(N)
end

function largest_factor(N::Integer)
    if N < 2
        return 1
    else
        f::Integer = 2
        while f <= sqrt(N)
            while N % f == 0
                N ÷= f
                if N == 1
                    return f
                end
            end
            f += 1
        end
        return N
    end
end

function solution003()
    N = 600851475143
    answer = e003(N)
    println("The answer to problem 003 is $answer")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
using Random

const nb_threads = parse(Int,get(ENV, "OPENBLAS_NUM_THREADS", "-1"))
if nb_threads == -1
  error("Matmul tests require explicitly setting OPENBLAS_NUM_THREADS")
end

const nsamples = 4

function timed_mul(A,x)
  times = zeros(nsamples)
  for i in 1:length(times)
    val, times[i], bytes, gctime, memallocs = @timed b = A*x
  end
  return extrema(times[2:end])
end

function write_timings(filename, distributed_dimension, mintime, maxtime)
  filepath = filename
  if !isabspath(filename)
    filepath = joinpath(@__DIR__, filename)
  end
  open(filepath, "a") do f
    if position(f) == 0
      write(f, "# nbprocs nbthreads distributed_dimension mintime (s) maxtime (s)\n")
    end
    write(f, "$nb_procs $nb_threads $distributed_dimension $mintime $maxtime\n")
  end
end

const N = 5000

if do_serial
  const As = rand(N,N)
  const xs = rand(N)

  const ref = As*xs

  (mintime, maxtime) = timed_mul(As,xs)
  write_timings("out-serial.txt", 1, mintime, maxtime)
end


text:
using AVLTrees, BenchmarkTools
using Random
using Plots
using DataFrames

function batch_insert!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        insert!(t, i, i)
    end
end

function batch_delete!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        delete!(t, i)
    end
end

function batch_find(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        i in t
    end
end

d = DataFrame((op=[], time=[], n=[]))
x = [1_000, 10_000, 100_000, 1_000_000, 10_000_000]

function prepare_t_insert(t)
    for i in nums_test
        delete!(t, i)
    end
    t
end

function prepare_t_delete(t)
    for i in nums_test
        insert!(t, i, i)
    end
    t
end

for attempt in 1:5
    for N in x
        global t = AVLTree{Int64,Int64}()
        rng = MersenneTwister(1111)
        global nums_fill = rand(rng, Int64, N)
        nn = 100
        global nums_test = rand(rng, Int64, nn)

        for i in nums_fill
            insert!(t, i, i)
        end

        insertion = @benchmark batch_insert!(_t, nums_test) evals=1 setup=(_t=prepare_t_insert(t))
        deletion = @benchmark batch_delete!(_t, nums_test) evals=1 setup=(_t=prepare_t_delete(t))

        batch_insert!(t, nums_test)
        search = @benchmark batch_find(t, nums_test)

        push!(d, ("insert", minimum(insertion).time/nn, N))
        push!(d, ("delete", minimum(deletion).time/nn,N))
        push!(d, ("search", minimum(search).time/nn,N))
        println("done $N")
    end
end

c = combine(groupby(d, [:op,:n]), :time => minimum)

plot(
    x,
    [c[(c.op.=="insert"),:].time_minimum,c[(c.op.=="delete"),:].time_minimum, c[(c.op.=="search"),:].time_minimum],
    xscale = :log10,
    # yscale = :log10,
    ylabel = "operation time [ns]",
    xlabel = "N",
    xticks = [1e3, 1e4, 1e5, 1e6, 1e7],
    markershape =[:diamond :utriangle :dtriangle],
    labels= ["insert" "delete" "lookup"],
    legend=:topleft,
)


text:
module Movement

import Params
export Point, randomwalk

type Point{T}
    x :: T
    y :: T
end

function Point{T}(arr :: Array{T, 1})
    if length(arr) != 2
        error("Points are in 2D.")
    end
    new(arr[1], arr[2])
end

function Point(T :: DataType)
    Point(zero(T), zero(T))
end

+(a :: Point, b :: Point) = Point(a.x + b.x, a.y + b.y)

#=
Return a function that makes a random walk according to a
Wiener process with given parameters.

Parameters:

- speed: Speed of random walk (m/s)
- unit_interval: Interval of a jump (timeslots)

Returns:
- a function that returns the difference (dx, dy) of random walk

=#
function randomwalk(speed :: Float64, unit_interval :: Int64)
    unit_speed :: Float64 = speed * sqrt(Params.t_slot)
    phi :: Float64 = rand() * 2*pi
    return t -> begin
        r = unit_speed * randn()
        if t % unit_interval == 0
            phi = rand() * 2*pi
        end
        return Point(r*cos(phi), r*sin(phi))
    end
end

end


text:
# Autogenerated wrapper script for libfacedetection_jll for x86_64-w64-mingw32
export libfacedetection

JLLWrappers.@generate_wrapper_header("libfacedetection")
JLLWrappers.@declare_library_product(libfacedetection, "libfacedetection.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfacedetection,
        "bin\\libfacedetection.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function linsolve_cg( LF::LF3dGrid, b::Array{Float64,1};
                      x0 = nothing,
                      NiterMax = 1000, TOL=5.e-10,
                      convmsg=true, showprogress=false )
  #
  Npoints = size(b)[1]
  if x0 == nothing
    x = zeros(Float64, Npoints)
  else
    x = copy(x0)
  end
  #
  r = zeros( Float64, Npoints )
  p = zeros( Float64, Npoints )
  #
  L_x = apply_Laplacian( LF, x )
  for ip = 1 : Npoints
    r[ip] = b[ip] - L_x[ip]
    p[ip] = r[ip]
  end

  rsold = dot( r, r )
  for iter = 1 : NiterMax
    #
    L_x = apply_Laplacian( LF, p )
    #
    α = rsold/dot( p, L_x )
    #
    x[:] = x[:] + α * p[:]  # FIXME use x[:] to force x to be copied, not referenced
    r[:] = r[:] - α * L_x[:]
    #
    rsnew = dot( r, r )
    # deltars = rsold - rsnew
    if showprogress
      @printf("%8d %20.10f\n", iter, sqrt(rsnew))
    end
    #
    if sqrt(rsnew) < TOL
      if convmsg
        @printf("#Convergence achieved in linsolve_cg: %8d iterations.\n", iter)
      end
      break
    end
    #
    p = r + (rsnew/rsold) * p
    #
    rsold = rsnew
  end
  #
  return x
  #
end


text:
model, obj0 = Dynamics.cartpole_analytical

obj = copy(obj0)
obj.x0 = [0;0;0;0.]
obj.xf = [0.5;pi;0;0]
obj.tf = 2.
u_bnd = 50
x_bnd = [0.6,Inf,Inf,Inf]
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd, x_min=-x_bnd, x_max=x_bnd)
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd)
dt = 0.05

N,dt = TrajectoryOptimization.calc_N(obj.tf, dt)

println("Hermite Simpson")
solver = Solver(model,obj,dt=dt)
N = solver.N
U0 = ones(1,N)*1
X0 = line_trajectory(obj.x0, obj.xf, N)
@time X,U = solve_dircol(solver,X0,U0,method=:hermite_simpson,grads=:quad)

solver = Solver(model,obj,dt=dt)
println("Trapezoid (no gradients)")
@time X,U = solve_dircol(solver,X0,U0,method=:trapezoid,grads=:none)
println("Trapezoid")
@time X,U = solve_dircol(solver,X0,U0,method=:trapezoid,grads=:quad)
println("HS Separated (no gradients)")
@time X,U = solve_dircol(solver,X0,U0,method=:hermite_simpson_separated,grads=:none)
println("HS Separated")
@time X,U = solve_dircol(solver,X0,U0,method=:hermite_simpson_separated,grads=:quad)
println("Hermite Simpson (no gradients)")
@time X,U = solve_dircol(solver,X0,U0,method=:hermite_simpson,grads=:none)
println("Hermite Simpson")
@time X,U = solve_dircol(solver,X0,U0,method=:hermite_simpson,grads=:quad)

# using Juno
# @profiler solve_dircol(solver,X0,U0,method=:hermite_simpson_separated,grads=:none)


text:
#two-sided diode-type test.
@testset "optimize radius power" begin
    λ0 = 1 #doesn't matter since everything is normalized to λ0
    k0 = 2π/λ0
    kin = 4k0
    ui = [PlaneWave(); PlaneWave(π)]
    a = λ0/2
    P = 6

    centers = rect_grid(4, 8, a, a) #symmetric 4x8 grid with distance a
    M = size(centers, 1)
    ids = ones(Int, M)
    φs = zeros(M)
    centers_abs = centers[:,1] + 1im*abs.(centers[:,2])
    ids, centers_abs = ParticleScattering.uniqueind(centers_abs)
    Nrs = maximum(ids) #if some are the same, Nrs is the number of different r
    initial_rs = 0.25*a*ones(Nrs)
    r_max = 0.45*a*ones(Nrs)
    r_min = 1e-4*a*ones(Nrs)
    shapes = CircleParams.(initial_rs)
    sp = ScatteringProblem(shapes, ids, centers, φs)
    Npoints = 15
    points1 = [2λ0*ones(Npoints) range(-λ0, stop=λ0, length=Npoints)]
    points2 = [-2λ0*ones(Npoints) range(-λ0, stop=λ0, length=Npoints)]
    nhat1 = [ones(Npoints) zeros(Npoints)]
    nhat2 = [-ones(Npoints) zeros(Npoints)]
    len = norm(points1[1,:] - points1[end,:]) #same for 2
    @assert verify_min_distance(sp, [points1;points2])

    fmm_options = FMMoptions(true, acc = 6, dx = λ0)
    groups,boxSize = divideSpace(sp.centers, fmm_options)
    P2,Q = FMMtruncation(fmm_options.acc, boxSize, k0)
    optim_options = Optim.Options(f_tol = 1e-5, iterations = 10, outer_iterations = 5,
                        store_trace = false, show_trace = false)
    # Allocate buffers
    opb = [OptimProblemBuffer(k0, kin, sp.centers, ui[1], P),
            OptimProblemBuffer(k0, kin, sp.centers, ui[2], P)]
    power_buffer = [optMatrixPwr(points1, sp.centers, M, P, k0, ui[1], 1, nhat1, len),
                    optMatrixPwr(points2, sp.centers, M, P, k0, ui[2], 2, nhat2, len)]
    mFMM = [FMMbuildMatrices(k0, P, P2, Q, groups, sp.centers, boxSize, tri = true) for i=1:2]
    buf = FMMbuffer(M, P, Q, length(groups))
    scatteringMatrices = [[sparse(one(Complex{Float64})I, 2P+1, 2P+1) for i = 1:Nrs] for i = 1:2]
    dS_S = [[sparse(one(Complex{Float64})I, 2P+1, 2P+1) for i = 1:Nrs] for i = 1:2]

    function fobj_testr(sv::Array{PowerBuffer})
        if sv[1].pow < 0 || sv[2].pow < 0
            return abs(sv[2].pow/sv[1].pow)*1e9 #instead of Inf
        end
        barr = -(log(sv[1].pow) + log(sv[2].pow)) #so they're > 0
        res = sv[2].pow/sv[1].pow
        res*(1 + barr)
    end

    function gobj_testr!(sv::Array{PowerBuffer}, opb::Array{OptimProblemBuffer})
        #calculate -(∂f/∂β)ᵀ for adjoint method
        if sv[1].pow < 0 || sv[2].pow < 0
            opb[1].rhs_grad[:] .= 0
            opb[2].rhs_grad[:] .= 0
            return
        end
        barr = -(log(sv[1].pow) + log(sv[2].pow))
        res = sv[2].pow/sv[1].pow
        opb[1].rhs_grad[:] = sv[1].∂pow*((1 + barr)*sv[2].pow/sv[1].pow^2)
        opb[2].rhs_grad[:] = sv[2].∂pow*(-(1 + barr)/sv[1].pow)
        opb[1].rhs_grad .+= sv[1].∂pow*(res/sv[1].pow)
        opb[2].rhs_grad .+= sv[2].∂pow*(res/sv[2].pow)
    end

    #initial, last must be different before first iteration
    last_rs = similar(initial_rs); @assert last_rs != initial_rs #initial_rs, last_rs must be different before first iteration
    df = Optim.OnceDifferentiable(
        rs -> optimize_pwr_rs_f(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, M, buf, fmm_options, φs, mFMM, fobj_testr),
        (grad_stor, rs) -> optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, M, buf, fmm_options, φs, mFMM, gobj_testr!),
        initial_rs)

    res = Optim.optimize(df, r_min, r_max, initial_rs, Optim.Fminbox(Optim.LBFGS()), optim_options)
    sp_after = ScatteringProblem(CircleParams.(res.minimizer), ids, centers, φs)
    power_after_right = calc_power(k0, kin, P, sp_after, points1, nhat1, ui[1])*len
    power_after_left = calc_power(k0, kin, P, sp_after, points2, nhat2, ui[2])*len
    @test power_after_right > 0
    @test power_after_left > 0
    @test power_after_right/power_after_left > 1
end


text:
######## Sampler for multifractional Brownian motion (mBm) and related processes ########

"""
    val2grid(X::Vector{Float64}, δ0::Float64=2.5e-2)

Construct a regular grid from a vector of continuous values with a given step.
"""
function val2grid(X::Vector{Float64}, δ0::Float64=2.5e-2)
    xmin, xmax = minimum(X), maximum(X)
    δ0 > 0 || error("Step of grid must be > 0.")
    # (hmin>0 && hmax<1) || error("Husrt exponent must be bounded in (0,1).")
    δ = min(δ0, xmax-xmin)  # hstep smaller than 2.5e-2 is unstable
    grid::Vector{Float64} = (δ > 0) ? collect(xmin:δ:xmax) : [xmin]
    M, N = length(grid), length(X)
    idx::Vector{Integer} = (δ > 0) ? min.(M, max.(1, ceil.(Int, M*(X-xmin)/(xmax-xmin)))) : ones(Integer, N)
    return grid, idx
end

"""
    linprd_fGns(Hgrid::Vector{<:AbstractFloat}, N::Integer, epsilon=1e-8)

Compute the coefficients of linear prediction for a sequence of fGn.

# Args
- Hgrid: grid points of Hurst exponents
- N: length of stationary process
- epsilon: threshold for identical Hgrid exponents

# Returns
- J: bool index of spanning subsequence
- Ψ: coefficients of basis
- Φ: coefficients of linear prediction
- R: variance of residual
- C: covariance sequence of residual
"""
function linprd_SfGn(Hgrid::Vector{<:AbstractFloat}, N::Integer, epsilon=1e-8)
    # Initialization
    M = length(Hgrid)
    ρ2 = 1. * N   # residual
    ρ = sqrt(abs(ρ2))
#     J::Vector{Integer} = [1] # index of spanning processes
    J::Vector{Bool} = zeros(M)  # bool index of spanning processes
    J[1] = true
#     A::Matrix{<:AbstractFloat} = eye(N)  # triangular system matrix for ϕ
#     F::Matrix{<:AbstractFloat} = zeros(N,N)  # triangular system matrix for ψ
#     F[1,1] = ρ
    iF::Matrix{<:AbstractFloat} = zeros(N,N)  # triangular system matrix for ψ, inverse
    iF[1,1] = 1/ρ

    Ψ = Vector{Vector{<:AbstractFloat}}(M)  # coeff of linear prediction
    Φ = Vector{Vector{<:AbstractFloat}}(M)  # coeff of orthonormal basis
    R = zeros(M); R[1] = ρ2
    C = zeros(N, M)
    cvec = zeros(N)

    # Iteration
    for j=2:M
        Jidx = find(J)
        γ = N * [mGn_cov(0, Hgrid[j], Hgrid[i]) for i in Jidx]  #RHS vector
#         G = N * [mGn_cov(0, Hgrid[i1], Hgrid[i2]) for i1 in Jidx, i2 in Jidx]  # system matrix of process, such that ϕ = G \ γ
        k = length(Jidx)
#         Ak = view(A, 1:k, 1:k)
#         Fk = view(F, 1:k, 1:k)
        iFk = view(iF, 1:k, 1:k)

        # Update by innovation
#         ϕ = Ak' * diagm(1 ./ R[Jidx]) * Ak * γ
#         ψ = Fk' * ϕ
        # another way, which seems more stable than the previous one
        ψ = iFk * γ
        ϕ = iFk' * ψ

        # Construct covariance matrix of the residual
        for n=0:N-1
            u = ϕ' * [mGn_cov(n, Hgrid[j], Hgrid[i]) for i in Jidx]
            v = ϕ' * [mGn_cov(n, Hgrid[i1], Hgrid[i2]) for i1 in Jidx, i2 in Jidx] * ϕ
            cvec[n+1] = mGn_cov(n, Hgrid[j], Hgrid[j]) - 2*u + v
        end

        # Residual
        ρ2 = N*cvec[1]
#         ρ2 = N - 2 * ϕ' * γ + ϕ' * G * ϕ
#         ρ2 = N - ϕ' * γ  # This gives negative values!
        ρ = sqrt(abs(ρ2))

        Ψ[j] = ψ
        Φ[j] = ϕ
        C[:,j] = cvec
        R[j] = ρ2

        # current process is accepted in the list of spanning processes
        if ρ/sqrt(N) > epsilon
            J[j] = true
            k = sum(J)
#             A[k,1:(k-1)] = -ϕ'
#             F[k,1:(k-1)] = ψ'
#             F[k,k] = ρ
            iF[k,1:(k-1)] = -1/ρ * ϕ'
            iF[k,k] = 1/ρ
        end
    end
    return J, Ψ, Φ, R, C
end

"""
Conditional sampling of a sequence of fGn
"""
function condsampl_SfGn(Hgrid::Vector{<:AbstractFloat},
                        J::Vector{Bool},
                        Ψ::Vector{Vector{<:AbstractFloat}},
                        Φ::Vector{Vector{<:AbstractFloat}},
                        R::Vector{<:AbstractFloat},
                        C::Matrix{<:AbstractFloat},
                        epsilon=1e-8)
    N = size(C,1)
    M = length(Hgrid)

    # Initialization
    Z = circulant_embedding(covseq(FractionalGaussianNoise(Hgrid[1]), 1:N))
    # Z = rand(CircSampler(FractionalGaussianNoise(Hgrid[1]), DiscreteTimeGrid(1:N)))

    X = zeros(N, M)  # historical sample paths
    X[:,1] = Z
    ρ = sqrt(abs(R[1]))
#     E = zeros(N, M)  # historical sample paths of residual
#     E[:,1] = Z/ρ

    for j=2:M
        ρ = sqrt(abs(R[j]))
        ϕ, ψ = Φ[j], Ψ[j]
        Jidx = find(J[1:j-1])

        # Conditional mean
        μ = X[:,Jidx] * ϕ
        # or use psi
#         μ = E[:,Jidx] * ψ

        # Conditional error
        Z = (ρ/sqrt(N) > epsilon) ? circulant_embedding(C[:,j]) : zeros(N)
        X[:,j] = μ + Z
#         E[:,j] = Z/ρ
    end
    return X
end

"""
Conditionalized sampling of a sequence of fGn (SfGn).
"""
struct SfGnSampler <: DiscreteTimeSampler{FractionalGaussianNoise}
    Hurst::Vector{<:AbstractFloat}
    grid::DiscreteTimeGrid  # grid in use for sampling
    hgrid::Vector{<:AbstractFloat}  # grid for Hurst exponents
    hidx::Vector{Integer}  # index of Hurst exponents in `hrid`

    J::Vector{Bool}  # bool index of spanning processes
    Ψ::Vector{Vector{<:AbstractFloat}}  # coeff of linear prediction
    Φ::Vector{Vector{<:AbstractFloat}}  # coeff of orthonormal basis
    R::Vector{<:AbstractFloat}  # square of residual
    C::Matrix{<:AbstractFloat}  # sequence of covariance

    function SfGnSampler(Hurst::Vector{<:AbstractFloat}, δ::Float64=2.5e-2)
        N = length(Hurst)
        grid = 1:N
        hgrid, hidx = val2grid(Hurst, δ)
        J, Ψ, Φ, R, C = linprd_SfGn(hgrid, N)
        new(Hurst, grid, hgrid, hidx, J, Ψ, Φ, R, C)
    end
end

function rand!(x::Vector{<:AbstractFloat}, s::SfGnSampler)
    # @assert length(x) <= length(s)

    dW0 = condsampl_SfGn(s.hgrid, s.J, s.Ψ, s.Φ, s.R, s.C)
    N = length(s.grid)
    dY = ((1/N).^ s.Hurst) .* [dW0[n, s.hidx[n]] for n=1:N]
    return copyto!(x, cumsum(dY)[1:length(x)])
end

# function rand_mBm(Hurst0, )
#     J, Ψ, Φ, R, C = FracFin.linprd_fGns(Hgrid, N);

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
# State transition function
function POMDPs.transition(mdp::HCAS_MDP, s::stateType, ra::actType)
    r = s[1]; t = s[2]; p = s[3]; vown = s[4]; vint = s[5]; pra = s[6];

    # Computation is faster when using vector of static size
    nextStates = MVector{9, stateType}(undef)
    nextProbs = @MVector(zeros(9))
    next_pra = ra
    ind=1

    # Compute probabilities of next states using sigma point sampling
    ownProbs, ownTurns = mdp.turns[pra]
    intProbs, intTurns = mdp.turns[-1]
    for i = 1:3
        for j = 1:3
            next_r,next_t,next_p,next_vown,next_vint = dynamics(r,t,p,vown,vint,ownTurns[i],intTurns[j],pra)
            nextStates[ind] = (next_r,next_t,next_p,next_vown,next_vint,next_pra)
            nextProbs[ind]  = ownProbs[i]*intProbs[j]
            ind+=1
        end
    end

    return SparseCat(nextStates,nextProbs)
end

# Dynamic equations
function dynamics(r::Float64,t::Float64,p::Float64,vown::Float64,vint::Float64,ownTurn::Float64, intTurn::Float64, ra::Int)
    x_own = 0.0; y_own = 0.0; x_int = r*cos(t); y_int = r*sin(t)
    dx_own = vown; dy_own = 0.0; dx_int = vint*cos(p); dy_int = vint*sin(p)
    x_own += dx_own; y_own += dy_own; x_int += dx_int; y_int += dy_int

    x_int_new = x_int - x_own; y_int_new = y_int - y_own
    heading_own = ownTurn; heading_int = p + intTurn

    r_new = sqrt(x_int_new.^2 + y_int_new.^2)
    t_new = atan(y_int_new,x_int_new) - heading_own
    p_new = heading_int - heading_own

    t_new = t_new > pi  ? t_new-2*pi : t_new
    t_new = t_new < -pi ? t_new+2*pi : t_new
    p_new = p_new > pi  ? p_new-2*pi : p_new
    p_new = p_new < -pi ? p_new+2*pi : p_new

    return r_new, t_new, p_new, vown, vint

end

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
export concat

function concat{T,N}(xs::Vector{Array{T,N}})
    p = Ptr{Cvoid}[0]
    attr = primitive_attributes()
    #[create_layout() for x in xs]
    #dnnConcatCreate_F32(p, attr, length(xs), )
end


text:
@test all(stocks.value .>= 0)
@test all(stocks.vol .>= 0)
@test all(stocks.impact .>= 0)

@test all(funds.holdings .>= 0)
@test all(funds.stakes .>= 0)
@test all(funds.value .>= 0)

@test all(investors.assets .>= 0)
@test all(investors.horizon .>= 0)
@test all(1 .>= investors.threshold .>= -1)

@test all(buyorder.values .>= 0)
@test all(buyorder.funds .>= 0)
@test all(sellorder.values .<= 0)
@test all(sellorder.investors .>= 0)


text:
using Base.Test
using QuantumOptics

fockbasis = FockBasis(40)
spinbasis = SpinBasis(1//2)

alpha = 0.5
a = full(destroy(fockbasis))
at = full(create(fockbasis))
n = full(number(fockbasis))

sx = full(sigmax(spinbasis))
sy = full(sigmay(spinbasis))
sz = full(sigmaz(spinbasis))
sp = full(sigmap(spinbasis))
sm = full(sigmam(spinbasis))

xket = coherentstate(fockbasis, alpha)
xbra = dagger(xket)
op1 = DenseOperator(spinbasis, GenericBasis([3]), [1 1 1; 1 1 1])
op2 = DenseOperator(GenericBasis([3]), spinbasis, [1 1; 1 1; 1 1])
I = dense_identity(fockbasis)

# Test creation
@test_throws DimensionMismatch DenseOperator(spinbasis, [1 1 1; 1 1 1])
@test_throws DimensionMismatch DenseOperator(spinbasis, FockBasis(3), [1 1; 1 1; 1 1])
@test_approx_eq 0. maximum(abs((dagger(op1)-op2).data))

# Test addition
@test_approx_eq 0. tracedistance(sx, sp + sm)
@test_throws bases.IncompatibleBases op1+op2

# Test substraction
@test_approx_eq 0. tracedistance(sy, -1im*(sp - sm))
@test_throws bases.IncompatibleBases op1-op2

# Test multiplication
@test_approx_eq 0. norm(I*xket - xket)
@test_approx_eq 0. norm(xbra*I - xbra)
@test_approx_eq alpha norm(a*xket)
@test_approx_eq alpha norm(xbra*at)
@test_approx_eq_eps 0. tracedistance(n, at*a) 1e-13
@test_approx_eq_eps 0. norm((5. * at)*xket - 5 * (at*xket)) 1e-13
@test_approx_eq_eps 0. norm((at * 5.)*xket - (at*xket) * 5) 1e-13
@test_throws bases.IncompatibleBases a*op1

# Test division
@test_approx_eq_eps 0. norm((at/5.)*xket - (at*xket)/5) 1e-13

# Test trace and normalize
op = DenseOperator(GenericBasis([3]), [1 3 2;5 2 2;-1 2 5])
@test_approx_eq 8. trace(op)
op_normalized = normalize(op)
@test_approx_eq 8. trace(op)
@test_approx_eq 1. trace(op_normalized)
op_ = normalize!(op)
@test op_ === op
@test_approx_eq 1. trace(op)

# Test operator exponential
b = GenericBasis([3])
op = DenseOperator(GenericBasis([3]), [1 3 2;5 2 2;-1 2 5])
op = op + dagger(op)
op /= norm(op.data)
v = SubspaceBasis(b, eigenstates_hermitian(op))
P = projector(v, b)
op_diag = P*op*dagger(P)
op_diag_exp = DenseOperator(v, diagm(exp(diag(op_diag.data))))
@test_approx_eq_eps 0. tracedistance(expm(op), dagger(P)*op_diag_exp*P) 1e-13

# Test identity function
@test full(I) == dense_identity(a)

# Test gemv implementation
result_ket = deepcopy(xket)
operators.gemv!(complex(1.0), at, xket, complex(0.), result_ket)
@test_approx_eq 0. norm(result_ket-at*xket)

result_bra = deepcopy(xbra)
operators.gemv!(complex(1.0), xbra, at, complex(0.), result_bra)
@test_approx_eq 0. norm(result_bra-xbra*at)


text:
#
# Linear Algebra of order 2 Tensors
#
# Even though DenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray
export eigs,
       truncerror,
       entropy

function Base.:*(T1::Tensor{ElT1,2,StoreT1,IndsT1},
                 T2::Tensor{ElT2,2,StoreT2,IndsT2}) where
                                       {ElT1,StoreT1<:Dense,IndsT1,
                                        ElT2,StoreT2<:Dense,IndsT2}
  RM = matrix(T1)*matrix(T2)
  indsR = IndsT1(ind(T1,1),ind(T2,2))
  return Tensor(Dense(vec(RM)),indsR)
end

function LinearAlgebra.exp(T::DenseTensor{ElT,2}) where {ElT}
  expTM = exp(matrix(T))
  return Tensor(Dense(vec(expTM)),inds(T))
end

function expHermitian(T::DenseTensor{ElT,2}) where {ElT}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(Hermitian(matrix(T))))
  return Tensor(Dense(vec(expTM)),inds(T))
end

"""
  Spectrum
contains the (truncated) density matrix eigenvalue spectrum which is computed during a
decomposition done by `svd` or `eigenHermitian`. In addition stores the truncation error.
"""
struct Spectrum
  eigs::Vector{Float64}
  truncerr::Float64
end

eigs(s::Spectrum) = s.eigs
truncerror(s::Spectrum) = s.truncerr

function entropy(s::Spectrum)
  S = 0.0
  for p in eigs(s)
    p > 1e-13 && (S -= p*log(p))
  end
  return S
end

#svd of an order-2 tensor
function LinearAlgebra.svd(T::DenseTensor{ElT,2,IndsT};
                           kwargs...) where {ElT,IndsT}
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  fastSVD::Bool = get(kwargs,:fastSVD,false)

  if fastSVD
    MU,MS,MV = svd(matrix(T))
  else
    MU,MS,MV = svd_recursive(matrix(T))
  end
  conj!(MV)

  P = MS.^2
  truncerr,_ = truncate!(P;mindim=mindim,
                         maxdim=maxdim,
                         cutoff=cutoff,
                         absoluteCutoff=absoluteCutoff,
                         doRelCutoff=doRelCutoff)
  spec = Spectrum(P,truncerr)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    resize!(MS,dS)
    MV = MV[:,1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T,1),u))
  Sinds = IndsT((u,v))
  Vinds = IndsT((ind(T,2),v))
  U = Tensor(Dense(vec(MU)),Uinds)
  S = Tensor(Diag(MS),Sinds)
  V = Tensor(Dense(vec(MV)),Vinds)
  return U,S,V,spec
end

function eigenHermitian(T::DenseTensor{ElT,2,IndsT};
                        kwargs...) where {ElT,IndsT}
  ispossemidef::Bool = get(kwargs,:ispossemidef,false)
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)

  DM,UM = eigen(Hermitian(matrix(T)))

  # Sort by largest to smallest eigenvalues
  p = sortperm(DM; rev = true)
  DM = DM[p]
  UM = UM[:,p]

  if ispossemidef
    truncerr,_ = truncate!(DM;maxdim=maxdim,
                           cutoff=cutoff,
                           absoluteCutoff=absoluteCutoff,
                           doRelCutoff=doRelCutoff)
    dD = length(DM)
    if dD < size(UM,2)
      UM = UM[:,1:dD]
    end
  else
    dD = length(DM)
    truncerr = 0.0
  end
  spec = Spectrum(DM,truncerr)

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dD)
  v = eltype(IndsT)(dD)
  Uinds = IndsT((ind(T,1),u))
  Dinds = IndsT((u,v))
  U = Tensor(Dense(vec(UM)),Uinds)
  D = Tensor(Diag(DM),Dinds)
  return U,D,spec
end

function LinearAlgebra.qr(T::DenseTensor{ElT,2,IndsT}) where {ElT,
                                                              IndsT}
  # TODO: just call qr on T directly (make sure
  # that is fast)
  QM,RM = qr(matrix(T))
  # Make the new indices to go onto Q and R
  q,r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = Tensor(Dense(vec(Matrix(QM))),Qinds)
  R = Tensor(Dense(vec(RM)),Rinds)
  return Q,R
end

function polar(T::DenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM,RM = polar(matrix(T))
  dim = size(QM,2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = Tensor(Dense(vec(QM)),Qinds)
  R = Tensor(Dense(vec(RM)),Rinds)
  return Q,R
end