text:
# %%
# Hans Martinez
# Code based off Emmanuel's Code

# cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro/Assignment8")
# mkpath("Figures")
# Pkg.add("Plots")
# Pkg.add.(["Optim" ,"Roots","Parameters","Distributions","QuadGK"])
# Pkg.add.(["Latexify","PrettyTables","StatsPlots"])
using Plots
# using LateXStrings # Pkg.add("LaTeXStrings") # https://github.com/stevengj/LaTeXStrings.jl
# Pkg.add("Dierckx")
using Dierckx # Pkg.add("Dierckx") # https://github.com/kbarbary/Dierckx.jl
# Pkg.add("Interpolations")
using Interpolations # Pkg.add("Interpolations") # https://github.com/JuliaMath/Interpolations.jl
# Pkg.add("ForwardDiff")
using ForwardDiff # Pkg.add("ForwardDiff") # https://github.com/JuliaDiff/ForwardDiff.jl
using Optim # Pkg.add("Optim") # https://julianlsolvers.github.io/Optim.jl/stable/
    using Optim: converged, maximum, maximizer, minimizer, iterations
using Roots # Pkg.add("Roots") # https://github.com/JuliaMath/Roots.jl
using Parameters # Pkg.add("Parameters") # https://github.com/mauro3/Parameters.jl
using Distributions #Pkg.add("Distributions")
using QuadGK # Pkg.add("QuadGK") # https://juliamath.github.io/QuadGK.jl/latest/
using LinearAlgebra
using Random
using Statistics
using Latexify
# Pkg.add("StatsPlots")
using StatsPlots
# Call Scaled Interpolation Functions
    include("../Scaled_Interpolation_Functions.jl")

# Pkg.add("SparseArrays")
using SparseArrays

# mkpath("Figures")
#
# using Random, Distributions
# using Statistics
# using LinearAlgebra
# using Plots
# using Interpolations
# using Dierckx
# using ForwardDiff
# using Optim
#     using Optim: converged, maximum, maximizer, minimizer, iterations
# using Roots
# using Parameters
# include("Scaled_Interpolation_Functions.jl")
# # using .VFI_Toolbox
# println(" ")
# println("------------------------")
# println("Hugget (1993) in Julia")
# println("PWD: ",pwd())
# println("Solve Hugget's model with EGM and several implementations of histogram method")
# println("------------------------")
# println(" ")

#-----------------------------------------------------------
#-----------------------------------------------------------
# Paramters and Model Structure
    # Generate structure for parameters using Parameters module
    # We can set default values for our parameters
    @with_kw struct Par
        # Model Parameters
        α::Float64 = 1/3  ; # Production function
        β::Float64 = 0.96 ; # Discount factor
        γ::Float64 = 2.0  ; # Relative risk aversion parameter
        δ::Float64 = 0.05 ; # Depreciation rate
        ρ::Float64 = 0.90 ; # Persistence of labor efficiency process
        σ::Float64 = 0.08 ; # Standard devaiation of labor efficiency innovation
        z_bar::Float64 = 1; # Reference level for productivity
        # VFI Paramters
        max_iter::Int64   = 100000; # Maximum number of iterations
        dist_tol::Float64 = 1E-6  ; # Tolerance for distance
        # Histogram iteration parameters
        Hist_max_iter     = 10000 ;
        Hist_tol          = 1E-5  ;
        # Histogram iteration parameters
        N_eq              = 1000  ;
        tol_eq            = 1E-6  ;
        # Minimum consumption for numerical optimization
        c_min::Float64    = 1E-16
    end

    # Allocate paramters to object p for future calling
    p = Par()

# Function to distretize AR(1) markov process with Rouwenhorst (1995)
function Rouwenhorst95(N,p::Par)
    @unpack ρ,σ=p
    # INPUTS:
        # ρ: persistence of unerlying AR(1) process where log(z') = ρlog(z)+η
        # σ_z: Std dev of inovation η in AR(1) process where η∼N(0,σ^2)
        # N: Size of grid for discrete process
    # OUTPUT:
        # z: All possible values of discretized AR(1) process, equally spaced grid of size N
        # Π: Matrix of transition probabilities
        # PDF_z: Stationary PDF of z
    #---------------------------------------------------------------------------
    Π = zeros(N,N)
    Π_Nm = zeros(N-1,N-1)
    P = (1+ρ)/2
    ϕ = σ*(sqrt((N-1)/(1-ρ^2)))
    z = range(-ϕ,ϕ;length=N)
    if N==2
        Π = [P 1-P;1-P P]
    else
        Π_Nm = Rouwenhorst95(N-1,p)[2]
        o = zeros(N-1)
        Π = P*[Π_Nm o; o' 0] + (1-P)*[o Π_Nm; 0 o'] + (1-P)*[o' 0; Π_Nm o] + P*[0 o';o Π_Nm]
        Π = Π./repeat(sum(Π,dims=2),1,N)
    end
    PDF_z = pdf.(Binomial(N-1,0.5),(0:N-1))
    return (z,Π,P)
end

# Function to get grid on assets
function Make_A_Grid(n_a,θ_a,a_min,a_max,p::Par)
    # Get k_grid
    if θ_a≠1
        a_grid = PolyRange(a_min,a_max;θ=θ_a,N=n_a)
    else
    a_grid = range(a_min,a_max,length=n_a)
    end
    # Return
    return a_grid
end

# %%
    # Generate structure of model objects
    @with_kw struct Model
        # Parameters
        p::Par = Par() # Model paramters in their own structure
        # Parameters for Asset Grid
        a_max::Float64  = 50                         # Default max node of a_grid
        θ_a::Float64    = 2.5                        # Default Curvature of a_grid
        n_a::Int64      = 200                        # Default Size of a_grid
        n_a_fine::Int64 = 1000                       # Default Size of fine grid for interpolation and distribution
        # Productivity process
        n_ϵ::Int64     = 11     # Default size of discretized grid for productivity as a markov process
        log_ϵ          = Rouwenhorst95(n_ϵ,p)[1]
        Π              = Rouwenhorst95(n_ϵ,p)[2]
        ϵ_grid         = exp.(log_ϵ)
        # Prices and aggregates
        r::Float64 = 0.05 # Initial guess for rate of return on assets
        # Lower bound on Asset Grid (borrowing limit)
        a_min = -ϵ_grid[1]/r
        # Capital grid
        a_grid = Make_A_Grid(n_a,θ_a,a_min,a_max,p)          # Grid on assets for model solution
        a_grid_fine = Make_A_Grid(n_a_fine,1,a_min,a_max,p)  # Fine grid on assets for interpolation
        # State matrices
        a_mat     = repeat(a_grid',n_ϵ,1)
        a_mat_fine= repeat(a_grid_fine',n_ϵ,1)
        ϵ_mat     = repeat(ϵ_grid,1,n_a)
        # Value and policy functions
        V         = Array{Float64}(undef,n_ϵ,n_a)       # Value Function
        G_ap      = Array{Float64}(undef,n_ϵ,n_a)       # Policy Function
        G_c       = Array{Float64}(undef,n_ϵ,n_a)       # Policy Function
        V_fine    = Array{Float64}(undef,n_ϵ,n_a_fine)  # Value Function on fine grid
        G_ap_fine = Array{Float64}(undef,n_ϵ,n_a_fine)  # Policy Function on fine grid
        G_c_fine  = Array{Float64}(undef,n_ϵ,n_a_fine)  # Policy Function on fine grid
        # Distribution
        Γ         = 1/(n_ϵ*n_a_fine)*ones(n_ϵ,n_a_fine) # Distribution (initiliazed to uniform)
        # Error in Euler equation
        Euler     = Array{Float64}(undef,n_ϵ,n_a_fine)  # Errors in Euler equation
    end

    # Allocate model to object M for future calling
    M = Model()

# Utility function
function utility(c,p::Par)
    if p.γ>1
    return (c).^(1-p.γ)/(1-p.γ)
    else
    return log.(c)
    end
end
function d_utility(c,p::Par)
    return (c).^(-p.γ)
end
function d_utility_inv(x,p::Par)
    return x.^(-1/p.γ)
end

# Get index for histogram method
function get_ind(a_grid,g_a)
    n = length(a_grid)
    if g_a > maximum(a_grid)
        return n
    elseif g_a < minimum(a_grid)
        return 1
    else
        return floor(((g_a-minimum(a_grid))/(maximum(a_grid)-minimum(a_grid)))*(n-1)+1)
    end
end

# %%

# PFI Fixed Point - Iterate over policy functions for a given guess of the interest rate r
function PFI_Fixed_Point(T::Function,M::Model,G_ap_old=nothing)
    # Unpack model structure
    @unpack p, n_ϵ, n_a, n_a_fine, θ_a, a_grid, a_grid_fine, r = M
    # PFI paramters
    @unpack max_iter, dist_tol = p
    # Initialize variables for loop
    if G_ap_old==nothing
        G_ap_old  = (1+r)*M.a_mat
    end
    G_dist = 1              ; # Initialize distance
    println(" ")
    println("------------------------")
    println("PFI - n_ϵ=$n_ϵ, n_a=$n_a - θ_a=$θ_a - r=$r")
    for iter=1:max_iter
        # Update value function
        G_ap_new, G_c = T(Model(M,G_ap=copy(G_ap_old)))
        # Update distance and iterations
        G_dist = sqrt(norm(G_ap_new-G_ap_old,2))
        # Update old function
        G_ap_old  = G_ap_new
        # Report progress every 250 iterations
        if mod(iter,250)==0
            println("   PFI Loop: iter=$iter, dist=",G_dist)
        end
        # Check convergence and return results
        if G_dist<=dist_tol
            println("PFI - n_ϵ=$n_ϵ, n_a=$n_a - θ_a=$θ_a - r=$r")
            println("Iterations = $iter and Distance = ",G_dist)
            println("------------------------")
            println(" ")
            # Interpolate to fine grid
            G_ap_fine = zeros(n_ϵ,n_a_fine)
            G_c_fine  = zeros(n_ϵ,n_a_fine)
            for i_ϵ=1:n_ϵ
            G_ap_ip = ScaledInterpolations(a_grid,G_ap_new[i_ϵ,:] , BSpline(Cubic(Line(OnGrid()))))
                G_ap_fine[i_ϵ,:].= G_ap_ip.(collect(a_grid_fine))
            G_c_ip  = ScaledInterpolations(a_grid,G_c[i_ϵ,:]  , BSpline(Cubic(Line(OnGrid()))))
                G_c_fine[i_ϵ,:] .= G_c_ip.(collect(a_grid_fine))
            end
            # Update model
            M = Model(M; G_ap=G_ap_new,G_c=G_c,G_ap_fine=G_ap_fine,G_c_fine=G_c_fine)
            # Return results
            return M
        end
    end
    # If loop ends there was no convergence -> Error!
    error("Error in PFI - Solution not found")
end

# Bellman operator - EGM - Iterate on Policy Functions
function T_EGM_G(M::Model)
    @unpack p, n_ϵ, n_a, G_ap, r, a_min, Π = M
    @unpack β = p
    # Define RHS of Euler equation for each (ϵ,a')
    # Rows are present ϵ and columns are tomorrow's a in fixed grid
    Euler_RHS = β*(1+r)*Π*d_utility( (1+r)*M.a_mat + M.ϵ_mat - G_ap , p )
    # Check Monotonicity
    if any( Euler_RHS.<0 )
        error("RHS must be monotone for EGM to work")
    end
    # Define consumption from Euler equation
    C_endo = max.(d_utility_inv(Euler_RHS,p),p.c_min)
    # Define endogenous grid on assets
    A_endo = (C_endo .+ M.a_mat - M.ϵ_mat)/(1+r)
    # Interpolate functions on exogenous grid
    G_c = Array{Float64}(undef,n_ϵ,n_a)
    for i_ϵ=1:n_ϵ
        # Sort A_endo for interpolation
        sort_ind = sortperm(A_endo[i_ϵ,:])
        A_aux    = A_endo[i_ϵ,:][sort_ind]
        C_aux    = C_endo[i_ϵ,:][sort_ind]
        # Check boundary condition
            # Ap(ϵ,a)=a_min for all a<min(A_aux)
            # Note that in that case C is linear between a_min and min(A_aux)
        if minimum(A_aux)>a_min
            a_vec = M.a_grid[M.a_grid.<minimum(A_aux)] # All assets today between a_min and min(A_aux)
            A_aux = [a_vec ; A_aux] # Append to existing grid of today's assets
            C_aux = [((1+r)*a_vec.+M.ϵ_grid[i_ϵ].-a_min) ; C_aux] # I can know consumption directly because all those assets imply a'=a_min
        end
        C_ip        = Spline1D(A_aux,C_aux;k=1)
        G_c[i_ϵ,:] .= C_ip.(M.a_grid) # interpolate consumption to exogenous grid of assets today
        #Ap_aux      = (1+r)*collect(M.a_grid) .+ M.ϵ_grid[i_ϵ] .- G_c[i_ϵ,:] # Find corresponding assets tomorrow
    end
    # Update policy function
    G_ap .= (1+r)*M.a_mat .+ M.ϵ_mat .- G_c
        # Adjust for numerical error
        for ind = findall(<=(1e-10),abs.(G_ap.-a_min))
            G_ap[ind] = a_min
            G_c[ind]  = (1+r)*M.a_mat[ind] + M.ϵ_mat[ind] - a_min
        end
        # Check for borrowing constraint
        if any( G_ap.<a_min )
            error("Borrowing Constraint Violated")
        end
    # Return Results
    return G_ap, G_c
end

# Value function (given policy function)
function Value_Function(M::Model)
    # Unpack model structure
    @unpack p, n_ϵ, n_a, n_a_fine, a_grid, a_grid_fine, r, G_ap, G_c, Π = M
    @unpack β, dist_tol = p

    # Compute value function with policy function iteration
    V      = zeros(n_ϵ,n_a)
    V_new  = zeros(n_ϵ,n_a)
    V_dist = 1
    U_mat  = utility(G_c,p)
    while V_dist>dist_tol
        for i_ϵ=1:n_ϵ
            Pr = Π[i_ϵ,:]'
            for i_a=1:n_a
                ap = G_ap[i_ϵ,i_a]
                Vp = zeros(n_ϵ)
                    for i_ϵp=1:n_ϵ
                        Vp_ip    = ScaledInterpolations(a_grid,V[i_ϵp,:], BSpline(Cubic(Line(OnGrid()))))
                        Vp[i_ϵp] = Vp_ip(ap)
                    end
                V_new[i_ϵ,i_a] = U_mat[i_ϵ,i_a] + β*(Pr*Vp)
            end
        end
        V_dist = maximum(abs.(V_new./V.-1))
        V      = V_new
    end
    # Interpolate to fine grid
    V_fine    = zeros(n_ϵ,n_a_fine)
    for i_ϵ=1:n_ϵ
        V_ip    = ScaledInterpolations(a_grid,V[i_ϵ,:], BSpline(Cubic(Line(OnGrid()))))
        V_fine[i_ϵ,:] .= V_ip.(collect(a_grid_fine))
    end
    # Update model
    M = Model(M; V=V,V_fine=V_fine)
    return M
end

# %%
# Compute stationary distribution with Histogram method
# Histogram method
function Histogram_Method_Loop(M::Model,N_H=nothing,Γ_0=nothing)
    @unpack a_min, p, n_ϵ, Π, n_a_fine, a_grid_fine, G_ap_fine = M
    @unpack Hist_max_iter, Hist_tol = p

    println("\n--------------------------------\nBegining Histogram Method with Loops")

    # Change max iter
    if N_H==nothing
        N_H = Hist_max_iter
    end

    # Initial distribution (uniorm in income shock and assets)
    if Γ_0==nothing
        Γ_0 = M.Γ
    end

    # Discretize distribution
    H_ind    = Array{Int64}(undef,n_ϵ,n_a_fine)
    H_weight = Array{Float64}(undef,n_ϵ,n_a_fine)
    a_max    = maximum(a_grid_fine)
    for i_ϵ=1:n_ϵ
    for i_a=1:n_a_fine
        H_ind[i_ϵ,i_a]    = get_ind(a_grid_fine,G_ap_fine[i_ϵ,i_a])
        # weight on node corresponding to H_ind
        if H_ind[i_ϵ,i_a] == n_a_fine || H_ind[i_ϵ,i_a] == 1 # Mass points above max or below min in the fine grid
            H_weight[i_ϵ,i_a] = 1
        else
            H_weight[i_ϵ,i_a] = 1-(G_ap_fine[i_ϵ,i_a]-a_grid_fine[H_ind[i_ϵ,i_a]])/(a_grid_fine[H_ind[i_ϵ,i_a]+1]-a_grid_fine[H_ind[i_ϵ,i_a]])
        end
    end
    end

    #Γ = zeros(n_ϵ,n_a_fine)
    ## Loop for updating histogram
    #H_dist = 1
    #for l = 1:N_H
    #    for i = 1:n_ϵ
    #        for j = 1:n_a_fine
    #            # all combinations of (ϵ,a) today that may lead to a_j tomorrow
    #            index = findall(x->x==j,H_ind)
    #            index_ =  findall(x->x==j,H_ind.+1)
    #            # from the above, take only ϵ for the income probability transition
    #            transindex = map(i->i[1], index)
    #            transindex_ = map(i->i[1], index_)
    #            # Update the distribution
    #            Γ[i,j] = Π[i,transindex]'*(H_weight[index].*Γ_0[index])+Π[i,transindex_]'*((-H_weight[index_].+1).*Γ_0[index_])
    #        end
    #    end
    #    # Update distance
    #    H_dist = maximum(abs.(Γ-Γ_0))
    #    # Update initial distribution
    #    Γ_0 .= Γ
    #    # Report progress
    #    println("   Histogram Loop: iter=$l, dist=$H_dist")
    #    #if mod(l,250)==0
    #    #    println("   Histogram Loop: iter=$l, dist=$H_dist")
    #    #end
    #    # Check convergence
    #    if H_dist<Hist_tol
    #        println("Histogram iteartion converged in iteration $i_H. H_dist=$H_dist\n--------------------------------\n")
    #        M = Model(M; Γ=Γ)
    #        return M
    #    end
    #end

    Γ = zeros(n_ϵ,n_a_fine)
    # Loop for updating histogram
    H_dist = 1
    for i_H=1:N_H
        # Update histogram
        Γ = zeros(n_ϵ,n_a_fine)
        for i_ϵ=1:n_ϵ # Current ϵ
        for i_a=1:n_a_fine # Current a
            i_ap = H_ind[i_ϵ,i_a]
            ω_ap = H_weight[i_ϵ,i_a]
            for i_ϵp=1:n_ϵ # Future ϵ
                Γ[i_ϵp,i_ap]   = Γ[i_ϵp,i_ap]   +    ω_ap *Π[i_ϵ,i_ϵp]*Γ_0[i_ϵ,i_a]
                if i_ap < n_a_fine
                    Γ[i_ϵp,i_ap+1] = Γ[i_ϵp,i_ap+1] + (1-ω_ap)*Π[i_ϵ,i_ϵp]*Γ_0[i_ϵ,i_a]
                end
            end
        end
        end
        # Update distance
        H_dist = maximum(abs.(Γ-Γ_0))
        # Update initial distribution
        Γ_0 .= Γ
        # Report progress
        if mod(i_H,50)==0
            println("   Histogram Loop: iter=$i_H, dist=$H_dist")
        end
        # Check convergence
        if H_dist<Hist_tol
            println("Histogram iteartion converged in iteration $i_H. H_dist=$H_dist\n--------------------------------\n")
            M = Model(M; Γ=Γ)
            return M
        end
    end

    # Return Results
    println("Histogram updated for $N_H iteartions. Current H_dist=$H_dist \n--------------------------------\n")
    M = Model(M; Γ=Γ_0)
    return M
end

# %%

# Do optimization then check market clearing
function market_clearing(r,M::Model)
    # Iterate over policy function given r
    M = PFI_Fixed_Point(T_EGM_G,Model(r=r))
    # Get value function
    M = Value_Function(M)
    # Find stationary distribution
    M = Histogram_Method_Loop(M)
    # Check for market clearing
    distance = sum(sum(M.G_ap_fine.*M.Γ,dims=2)).^2
    println("Squared distance from market clearing=$distance at intertest rate = $r")
    return distance
end

# Outer loop to find interest rate that clears market
function S_RCE(M::Model)
    @unpack p = M
    @unpack β = p
    # Maximum interest rate
    r_max = 1/β-1
    # Minimum interest rate
    r_min = 0
    # Find equilibrium interest rate
    min_result = optimize(x->market_clearing(x,M).^2,r_min,r_max)
    # Check result
    converged(min_result) || error("Failed to clear asset market in $(iterations(min_result)) iterations")
    # Upddate policy function
    r  = min_result.minimizer
    println("Equilibrium found in $(iterations(min_result)) iterations: r=$r")
    # Load equilibrium interest rate into model
    M = Model(M; r=r)
    # Compute policy function, value function and stationary distribution
    M = PFI_Fixed_Point(T_EGM_G,M)
    M = Value_Function(M)
    M = Histogram_Method_Loop(M)
    return M
end

# %%

# Function that makes 3-d plots of the value function, policy function and asset distributon
function plot_hugget(M::Model)
    # Interpolate the value function, policy function and asset distribution along the income shock dimension
    ϵ_grid_fine = range(M.ϵ_grid[1],M.ϵ_grid[end],length=M.n_a_fine)
    V_fine3d = zeros(M.n_a_fine,M.n_a_fine)
    G_ap_fine3d = zeros(M.n_a_fine,M.n_a_fine)
    Γ_fine3d = zeros(M.n_a_fine,M.n_a_fine)
    for i in 1:M.n_a_fine
        Γ_fine3d[:,i] = Spline1D(M.ϵ_grid,M.Γ[:,i];k=1).(collect(ϵ_grid_fine))
        V_fine3d[:,i] = Spline1D(M.ϵ_grid,M.V_fine[:,i];k=1).(collect(ϵ_grid_fine))
        G_ap_fine3d[:,i] = Spline1D(M.ϵ_grid,M.G_ap_fine[:,i];k=1).(collect(ϵ_grid_fine))
    end
    # Normalize distribution such that it sums to one (need because I interpolated)
    Γ_fine3d = Γ_fine3d./(sum(sum(Γ_fine3d)))
    # Plots of the value function, policy function and asset distributon
    gr()
    plot(ϵ_grid_fine,M.a_grid_fine,V_fine3d, st=:contour,xlabel="Income",ylabel="Assets",title="Value Function - Hugget") # Surface plot
    savefig("./Figures/surface_vf_hugget")
    plot(ϵ_grid_fine,M.a_grid_fine,G_ap_fine3d, st=:contour,xlabel="Income",ylabel="assets",title="Asset policy - Hugget") # Surface plot
    savefig("./Figures/surface_policy_hugget")
    plot(ϵ_grid_fine,M.a_grid_fine,Γ_fine3d, st=:surface,xlabel="Income",ylabel="Assets",title="Distribution - Hugget") # Surface plot
    savefig("./Figures/surface_dist_hugget")
end

# %%
# Call PFI
M_hugget = S_RCE(Model(n_a=100,n_ϵ=11))
# Get plots
    plot_hugget(M_hugget)

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


text:
include("header.jl")
using Knet: KnetMems, KnetPtr, gpuCount, blocksize, initKnetMems

function _testkptr(kptrs, navail, nfree)
    # sump = sum(p->p.nptr, values(mem.pools))
    # @show (mem.limit, mem.bytes, mem.avail, length(mem.pools), sump)
    # klens = kptrs == nothing ? [] : (p->p.len).(kptrs)
    # @show klens
    mem = KnetMems[gpu()+1]
    blocksizes = blocksize.(2 .^ (1:10))
    @test length(KnetMems) == gpuCount()
    @test length(mem.pools) == 10
    @test sort(collect(keys(mem.pools))) == blocksizes
    @test mem.limit >= mem.bytes
    @test mem.bytes == sum(blocksizes)
    @test mem.avail == navail
    @test all(Bool[v.nptr==1 && length(v.free)==nfree for (k,v) in mem.pools])
    if nfree == 0
        @test (p->p.len).(kptrs) == blocksizes
    end
end

_testingkptr = false

if gpu() >= 0 && KnetMems === nothing
    initKnetMems()
    @testset "kptr:alloc"   begin; _testkptr(KnetPtr.(2 .^ (1:10)), 0, 0); end
    _testingkptr = true
end

GC.gc() # gc does not work reliably inside function, module, let, if, @testset etc.

if _testingkptr
    @testset "kptr:gc"      begin; _testkptr(nothing, sum(blocksize.(2 .^ (1:10))), 1); end
    @testset "kptr:realloc" begin; _testkptr(KnetPtr.(2 .^ (1:10)), 0, 0); end
    _testingkptr = false
end

nothing

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_PolicyAssignment.jl")
include("model_PolicyAssignmentListResult.jl")
include("model_PolicyAssignmentProperties.jl")
include("model_PolicyDefinition.jl")
include("model_PolicyDefinitionListResult.jl")
include("model_PolicyDefinitionProperties.jl")


text:
using ApproxFun, Base.Test

@test_approx_eq Fun([1,1.,1.],Laurent([0,2π]))(0.1) 1+2cos(0.1+π)
@test_approx_eq Fun([1,1.,1.],Laurent([-1,1]))(0.1) 1+2cos(π*0.1)
@test_approx_eq Fun([1,1.,1.],Laurent([0,1]))(0.1) 1+2cos(2π*(0.1-1/2))

@test abs(Fun(cos,Circle())(exp(0.1im))-cos(exp(0.1im)))<100eps()
@test abs(Fun(cos,Circle())'(exp(0.1im))+sin(exp(0.1im)))<100eps()
@test abs(Fun(cos,Circle())'(exp(0.1im))+Fun(sin,Circle())(exp(0.1im)))<100eps()

@test norm(Fun(x->Fun(cos,Fourier,20)(x),20)-Fun(cos,20)) <100eps()
@test norm(Fun(x->Fun(cos,Fourier)(x))-Fun(cos)) <100eps()
@test norm(Fun(cos,Fourier)'+Fun(sin,Fourier)) < 100eps()
@test norm(Fun(x->Fun(cos,Laurent)(x))-Fun(cos)) <100eps()
@test norm(Fun(cos,Laurent)'+Fun(sin,Laurent)) < 100eps()
@test norm(Fun(cos,Circle())'+Fun(sin,Circle()))<100eps()

for f in (Fun(θ->sin(sin(θ)),SinSpace()),Fun(θ->cos(θ)+cos(3θ),CosSpace()),
            Fun(θ->sin(sin(θ)),Fourier()),Fun(θ->cos(θ)+cos(3θ),CosSpace()))
    @test norm(integrate(f)'-f)<eps()
end

let f=Fun(exp,Circle())
    @test norm(f'-f)<100eps()
    @test norm(integrate(f)+1-f)<100eps()
end

let f=Fun(x->exp(-10sin((x-.1)/2)^2),Fourier)
    @test_approx_eq real(f)(.1) f(.1)
end

@test_approx_eq((Fun(z->sin(z)*cos(1/z),Circle())*Fun(z->exp(z)*airyai(1/z),Circle()))(exp(.1im)),
                (z->sin(z)*cos(1/z)*exp(z)*airyai(1/z))(exp(.1im)))

## Calculus

let f=Fun(t->cos(t),CosSpace)
    D=Derivative(space(f))
    @test_approx_eq (D*f)(.1) -sin(.1)
    @test_approx_eq f'(.1) -sin(.1)
end

let f=Fun(t->sin(t),SinSpace)
    D=Derivative(space(f))
    @test_approx_eq (D*f)(.1) cos(.1)
    @test_approx_eq f'(.1) cos(.1)
end

let f=Fun(cos,Fourier)
    @test norm((Derivative(space(f))^2)*f+f)<10eps()
end

## Taylor

@test Fun(Taylor())  == Fun([0.,1.],Taylor())

@test Fun(Taylor())(1.0) ≈ 1.0
@test Fun(Taylor(Circle(0.1,2.2)))(1.0) ≈ 1.0
@test Fun(Taylor(Circle(0.1+0.1im,2.2)))(1.0) ≈ 1.0

@test Multiplication(Fun(Taylor()),Taylor())[1:3,1:3] == [0. 0. 0.; 1. 0. 0.; 0. 1. 0.]

for d in (Circle(),Circle(0.5),Circle(-0.1,2.))
    S=Taylor(d)
    D=Derivative(S)
    ef=Fun(exp,S)
    @test norm((D*ef-ef).coefficients)<1000eps()
    @test norm((D^2*ef-ef).coefficients)<100000eps()
    u=[Evaluation(S,0.),D-I]\[1.]
    @test norm((u-ef).coefficients)<100eps()
    @test norm((Integral(S)*Fun(exp,S)+ef.coefficients[1]-ef).coefficients)<100eps()

    f=Fun(z->exp(1/z)-1,Hardy{false}(d))
    df=Fun(z->-1/z^2*exp(1/z),Hardy{false}(d))
    @test norm((Derivative()*f-df).coefficients)<1000eps()
    @test norm((Derivative()^2*f-df').coefficients)<100000eps()
    @test norm((f'-df).coefficients)<1000eps()
end

d=Circle()
S=Taylor(d)
D=Derivative(S)
D-I
ef=Fun(exp,S)
@test norm((D*ef-ef).coefficients)<1000eps()
@test norm((D^2*ef-ef).coefficients)<100000eps()
u=[Evaluation(S,0.),D-I]\[1.]

# check's Derivative constructor works
D=Derivative(Taylor(PeriodicInterval()))

## Multiplication

s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
b=Fun(t->(sin(t)+sin(3t)),SinSpace)

@test_approx_eq (s*s)(.1) s(.1)^2
@test_approx_eq (s*b)(.1) s(.1)*b(.1)

s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

@test_approx_eq (s*s)(.1) s(.1)^2
@test_approx_eq (s*b)(.1) s(.1)*b(.1)

s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
b=Fun(t->(sin(t)+sin(3t)),SinSpace)

@test_approx_eq (s*b)(.1) s(.1)*b(.1)

s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

@test_approx_eq (s*b)(.1) s(.1)*b(.1)

a=Fun(t->exp(cos(t)+sin(t)),Fourier)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test_approx_eq (a*b)(.1) a(.1)*b(.1)

a=Fun(t->exp(cos(t)),CosSpace)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test_approx_eq (a*b)(.1) a(.1)*b(.1)

a=Fun(t->sin(sin(t)),SinSpace)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test_approx_eq (a*b)(.1) a(.1)*b(.1)

# Check bug in off centre Circle
c2=-0.1+.2im;r2=0.3;
d2=Circle(c2,r2)
z=Fun(identity,d2)

@test_approx_eq z(-0.1+.2im+0.3*exp(0.1im)) (-0.1+.2im+0.3*exp(0.1im))

# false Circle
@test_approx_eq Fun(exp,Fourier(Circle(0.,1.,false)))(exp(0.1im)) exp(exp(.1im))
@test_approx_eq Fun(exp,Laurent(Circle(0.,1.,false)))(exp(0.1im)) exp(exp(.1im))

## Reverse orientation

f=Fun(z->1/z,Taylor(1/Circle()))
@test_approx_eq f(exp(0.1im)) exp(-0.1im)

## exp(z)

z=Fun(identity,Circle())
cfs=exp(z).coefficients[1:2:end]
for k=1:length(cfs)
    @test_approx_eq_eps cfs[k] 1/factorial(1.0(k-1)) 1E-10
end

##  Norms

@test_approx_eq sum(Fun([1.],CosSpace()))/π 1.
@test_approx_eq sum(Fun([0.,1.],CosSpace())^2)/π 0.5
@test_approx_eq sum(Fun([0.,0.,1.],CosSpace())^2)/π 0.5
@test_approx_eq sum(Fun([0.,0.,0.,1.],CosSpace())^2)/π 0.5

@test_approx_eq sum(Fun([0.,1.],SinSpace())^2)/π 0.5
@test_approx_eq sum(Fun([0.,0.,1.],SinSpace())^2)/π 0.5
@test_approx_eq sum(Fun([0.,0.,0.,1.],SinSpace())^2)/π 0.5


text:
# This file was generated, do not modify it. # hide
DecisionTree = @load DecisionTreeClassifier pkg=DecisionTree # model type
model = DecisionTree(min_samples_split=5)                    # model instance

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@timedtestset "test get_matrices rrule" for ix in 1:10
    ψ_arr = rand(ComplexF64, (8, 8, 8))

    function f(arr::Array{ComplexF64, 3})
        ψ = convert_to_cmps(arr)
        Q, R = get_matrices(ψ)
        return norm(get_matrices(ψ)[1]) + norm(get_matrices(ψ)[2])
    end

    @test f'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f, ψ_arr)[1]
end

@timedtestset "test transf_mat and transf_mat_T" for ix in 1:10
    ψ = cmps(rand, 8, 12)
    ϕ = cmps(rand, 9, 12)
    lop = transf_mat(ψ, ϕ)
    lop_T = transf_mat_T(ψ, ϕ)
    v = TensorMap(rand, ComplexF64, ℂ^9, ℂ^8)

    Iψ = id(ℂ^8)
    Iϕ = id(ℂ^9)

    Qψ_arr, Rψ_arr, Iψ_arr = toarray(ψ.Q), toarray(ψ.R), toarray(Iψ)
    Qϕ_arr, Rϕ_arr, Iϕ_arr = toarray(ϕ.Q), toarray(ϕ.R), toarray(Iϕ)
    v_arr = toarray(v)

    lop_arr = ein"db,ca->dcba"(Qϕ_arr, conj(Iψ_arr)) +
              ein"db,ca->dcba"(Iϕ_arr, conj(Qψ_arr)) +
              ein"dmb,cma->dcba"(Rϕ_arr, conj(Rψ_arr))
    lop_T_arr = ein"dcba->badc"(conj(lop_arr))

    lop_v_arr = ein"dcba,ba->dc"(lop_arr, v_arr)
    @test lop_v_arr ≈ toarray(lop(v))
    lop_T_v_arr = ein"badc,dc->ba"(lop_T_arr, v_arr)
    @test lop_T_v_arr ≈ toarray(lop_T(v))
end

# todo: improve the test cases for cmps canonical
@timedtestset "test left canonical for cmps" for ix in 1:10
    ψ = cmps(rand, 8, 3)
    _, ψL = left_canonical(ψ)

    @tensor result[-1; -2] := ψL.Q'[-1, -2] + ψL.Q[-1, -2] + ψL.R[1, 2, -2] * ψL.R'[-1, 1, 2]
    @test isapprox(result, TensorMap(zeros, ComplexF64, ℂ^8, ℂ^8), atol=sqrt(eps()))
end

# todo: improve the test cases for cmps canonical
@timedtestset "test right canonical for cmps" for ix in 1:10
    ψ = cmps(rand, 8, 3)
    _, ψR = right_canonical(ψ)

    @tensor result[-1; -2] := ψR.Q'[-1, -2] + ψR.Q[-1, -2] + ψR.R[-1, 2, 1] * ψR.R'[1, -2, 2]
    @test isapprox(result, TensorMap(zeros, ComplexF64, ℂ^8, ℂ^8), atol=sqrt(eps()))
end

#@timedtestset "normalized after canonicalization" for ix in 1:10
#    ψ = cmps(rand, 9, 4)
#    _, ψR = right_canonical(ψ)
#    X, _ = left_canonical(ψR)
#    _, s, _ = tsvd(X)
#    s = diag(s.data)
#
#    @test isapprox(sum(s .^ 2), 1, rtol=eps()^(1/4))
#
#    ψ = cmps(rand, 9, 4)
#    _, ψL = left_canonical(ψ)
#    Y, _ = right_canonical(ψL)
#    _, s, _ = tsvd(Y)
#    s = diag(s.data)
#
#    @test isapprox(sum(s .^ 2), 1, rtol=eps()^(1/4))
#end
#

@timedtestset "test ovlp and act" for ix in 1:10
    ψ = cmps(rand, 9, 3)
    ϕ = cmps(rand, 5, 3)
    op = cmpo_xxz(0.321)
    idψ = id(ℂ^9)
    idϕ = id(ℂ^5)
    idop = id(ℂ^2)

    Qψ_arr, Qϕ_arr = convert_to_array(ψ.Q), convert_to_array(ϕ.Q)
    Rψ_arr, Rϕ_arr = convert_to_array(ψ.R), convert_to_array(ϕ.R)
    Qop_arr, Lop_arr, Rop_arr = convert_to_array(op.Q), convert_to_array(op.L), convert_to_array(op.R)
    idψ_arr, idϕ_arr, idop_arr = convert_to_array(idψ), convert_to_array(idϕ), convert_to_array(idop)

    T = kron(conj(Qψ_arr), idϕ_arr) +
        kron(conj(idψ_arr), Qϕ_arr)
    for ix in 1:3
        T += kron(conj(Rψ_arr[:, ix, :]), Rϕ_arr[:, ix, :])
    end

    w, _ = eigen(T)
    @test logsumexp(w*6.123) ≈ log_ovlp(ψ, ϕ, 6.123)

    T2 = kron(conj(Qψ_arr), idop_arr, idϕ_arr) +
         kron(conj(idψ_arr), Qop_arr, idϕ_arr) +
         kron(conj(idψ_arr), idop_arr, Qϕ_arr)
    for ix in 1:3
        T2 += kron(conj(Rψ_arr[:, ix, :]), Rop_arr[:, ix, :], idϕ_arr) +
              kron(conj(idψ_arr), (Lop_arr[:, ix, :])', Rϕ_arr[:, ix, :])
    end
    w2, _ = eigen(T2)
    @test logsumexp(w2*6.123) ≈ log_ovlp(ψ, act(op, ϕ), 6.123)
end

# test rrule for log_ovlp
@timedtestset "rrule for cmps log_ovlp" for ix in 1:10
    ψ_arr = rand(ComplexF64, (3, 4, 3))
    ϕ_arr = rand(ComplexF64, (4, 4, 4))

    function f1(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        return real(log_ovlp(ψ, ψ, 6))
    end
    @test f1'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f1, ψ_arr)[1]

    function f2(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        ϕ = convert_to_cmps(ϕ_arr)
        return real(log_ovlp(ψ, ϕ, 6))
    end
    @test f2'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f2, ψ_arr)[1]
end

# test rrule for acting cmpo on cmps
@timedtestset "rrule for acting cmpo on cmps" for ix in 1:10
    ψ_arr = rand(ComplexF64, (3, 4, 3))
    ϕ_arr = rand(ComplexF64, (4, 4, 4))

    op = cmpo_xxz(0.321)

    function f1(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        return real(log_ovlp(ψ, act(op, ψ), 1.234))
    end
    @test f1'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f1, ψ_arr)[1]

    function f2(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        ϕ = convert_to_cmps(ϕ_arr)

        return real(log_ovlp(ψ, act(op, ϕ), 4.321))
    end
    @test f2'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f2, ψ_arr)[1]

end

# test finite_env and its rrule
@timedtestset "test finite_env and its rrule" for ix in 1:10
    ψ = cmps(rand, 6, 6)
    K = K_mat(ψ, ψ)
    L = rand()
    @test finite_env(K, L) ≈ exp(L*K) / tr(exp(L*K))

    # test rrule
    ψ_arr = rand(ComplexF64, (3, 4, 3))
    ϕ_arr = rand(ComplexF64, (4, 4, 4))
    O = TensorMap(rand, ComplexF64, (ℂ^3)'*ℂ^4, (ℂ^3)'*ℂ^4)

    function f1(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        ϕ = convert_to_cmps(ϕ_arr)
        K = K_mat(ψ, ϕ)
        return tr(finite_env(K, L) * O) |> real
    end

    function f2(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        ϕ = convert_to_cmps(ϕ_arr)
        K = K_mat(ψ, ϕ)
        return norm(finite_env(K, L))
    end

    function f3(arr::Array{<:Complex})
        ψ = convert_to_cmps(arr)
        K = K_mat(ψ, ψ)
        return tr(finite_env(K, L) * K) |> real
    end

    @test f1'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f1, ψ_arr)[1]
    @test f2'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f2, ψ_arr)[1]
    @test f3'(ψ_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f3, ψ_arr)[1]
end

# test normalization
@timedtestset "finite periodic cmps normalization" for ix in 1:10
    ψ = cmps(rand, 8, 8)
    ψ1 = normalize(ψ, 12)

    @test log_ovlp(ψ1, ψ1, 12; sym=true) < 1e-12
end

# test gauge_fixing_proj
@timedtestset "test gauge fixing" for ix in 1:10
    L = 10*rand()
    χ, d = 4, 1
    ψ = cmps(rand, χ, d)
    K = K_mat(ψ, ψ)
    ρ = finite_env(K, L)
    P = gauge_fixing_proj(ψ, L)

    A = cmps(id(ℂ^get_chi(ψ)), copy(ψ.R))
    A = convert_to_tensormap(A)
    V = TensorMap(rand, ComplexF64, ℂ^(d*χ), ℂ^χ)
    G = P * V

    @tensor M[-1, -2; -3, -4] := A'[-3, -1, 1] * G[-2, 1, -4]
    ρM = ρ * M
    @test abs(tr(ρM)) < 1e-14

    @tensor empty[-1; -2] := ρM[1, -1; 1, -2]
    @test findmax(abs.(empty.data))[1] < 1e-14
end

# test tangent_map
@timedtestset "test tangent_map" for ix in 1:10
    L = 50
    χ, d = 4, 1
    p = rand() * 2 * pi / L
    ψ = 0.01*cmps(rand, χ, d)
    ψ = normalize(ψ, L; sym=false)
    K = K_mat(ψ, ψ)
    W, UR = eig(K)
    UL = inv(UR)

    P = gauge_fixing_proj(ψ, L)

    A = cmps(id(ℂ^get_chi(ψ)), copy(ψ.R))
    A = convert_to_tensormap(A)
    V = TensorMap(rand, ComplexF64, ℂ^(d*χ), ℂ^χ)

    lop = tangent_map(ψ, L, p)

    #manual integration
    @tensor A1[-1, -2; -3, -4] := V[1, -4] * P[-2, 2, 1] * A'[-3, -1, 2]
    normψ = tr(exp(L*K))

    function fVc2(tau::Real)
        A2 = permute(exp(K*(L-tau)) * A1 * exp(K*(tau)), (2, 3), (4, 1))
        @tensor Vtmp[-1; -2] := A2[1, 3, 2, -2] * A[2, 4, 1] * P'[-1, 3, 4]
        return exp(-p*tau*im) * Vtmp
    end
    Vc2 = quadgk(fVc2, 0, L)[1]

    ρ = permute(exp(L*K), (2, 3), (4, 1))
    D = id(ℂ^(d+1))
    D.data[1] = 0
    @tensor Vd2[-1; -2] := ρ[3, 5, 2, -2] * P[2, 4, 1] * V[1, 3] * D[6, 4] * P'[-1, 5, 6]
    V2 = L*(Vc2 + Vd2) / normψ

    # compare
    @test norm(lop(V) - V2) / L < 1e-14
end

@timedtestset "test util functions for h tangent map" for ix in 1:10
    L = rand()
    p, q = 2*pi/L, 4*pi/L
    δ=1e-6
    theta2 = (a, b) -> statmech_tm_solver.theta2(L, a, b)
    theta3 = (a, b, c) -> statmech_tm_solver.theta3(L, a, b, c)

    M = TensorMap(rand, ComplexF64, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    A = TensorMap(rand, ComplexF64, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    B = TensorMap(rand, ComplexF64, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    C = TensorMap(rand, ComplexF64, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    W, VR = eig(M)
    VL = inv(VR)
    Wvec = diag(W.data)
    At = VL * A * VR
    Bt = VL * B * VR
    Ct = VL * C * VR

    # test elem_mult_f1 and theta2
    @test isapprox(theta2(1, 2), theta2(2, 1))
    @test isapprox(theta2(1, 1+δ), theta2(1, 1); atol=sqrt(δ))

    At1 = statmech_tm_solver.elem_mult_f1(At, (ix, iy) -> theta2(Wvec[ix], Wvec[iy] + im*p))
    @test tr(At1 * Bt) ≈ quadgk(x1 -> exp(im*p*x1)tr(A * exp(x1*M) * B * exp((L-x1)*M)), 0, L)[1]

    # test elem_mult_f2 and theta3
    @test theta3(1, 2, 3) ≈ theta3(2, 3, 1) ≈ theta3(3, 1, 2)
    @test theta3(1, 1, 2) ≈ theta3(2, 1, 1) ≈ theta3(1, 2, 1)

    @test isapprox(theta3(1, 1+δ, 2), theta3(1+δ, 1, 2); atol=δ)
    @test isapprox(theta3(1, 1, 2), theta3(1+δ, 1, 2); atol=sqrt(δ))
    @test isapprox(theta3(1, 1, 1), theta3(1, 1+δ, 1); atol=sqrt(δ))
    @test isapprox(theta3(1, 1-δ, 1+δ), theta3(1, 1, 1); atol=sqrt(δ))

    ABt = statmech_tm_solver.elem_mult_f2(At, Bt, (ix, iz, iy) -> theta3(Wvec[ix] - im*p, Wvec[iz] - im*(p+q), Wvec[iy]))
    @test tr(ABt * Ct) ≈ quadgk(x1 -> quadgk(
        x2 -> exp(im*p*x1) * exp(-im*(p+q)*x2) * tr(A * exp(x2*M) * B * exp((x1 - x2)*M) * C * exp((L-x1)*M) ),
        0, x1)[1], 0, L)[1]
end

@timedtestset "test lieb liniger h tangent map; test for d=1" for ix in 1:10

    k0, μ, c = rand(3)
    χ, d, L = 2, 1, rand()
    p, q = 2*pi/L, 2*2*pi/L

    ψ = cmps(rand, χ, d)
    ψ = normalize(ψ, L; sym=false)
    P = gauge_fixing_proj(ψ, L)
    Vx2 = TensorMap(rand, ComplexF64, ℂ^(χ*d), ℂ^χ)
    Vx1= TensorMap(rand, ComplexF64, ℂ^(χ*d), ℂ^χ)

    # overlap by tangent map function
    lop = lieb_liniger_h_tangent_map(ψ, p, q, c, L, μ; k0=k0)
    aaa = tr(Vx2' * lop(Vx1))

    # manual integration + alternative construction of My
    ψtn = convert_to_tensormap(ψ)
    K = K_mat(ψ, ψ)

    # to be used in construction of M matrices (My, Mx1, Mx2, see below) with tensorkit
    A = cmps(id(ℂ^χ), copy(ψ.R))
    A = convert_to_tensormap(A)
    DR = id(ℂ^(d+1))
    DR.data[1] = 0
    DQ = zero(DR)
    DQ.data[1] = 1
    DQRRQ = zeros(d+1, d+1, d+1, d+1)
    DRQQR = zeros(d+1, d+1, d+1, d+1)
    for ix in 2:d+1
        DQRRQ[1, ix, ix, 1] = 1
        DRQQR[ix, 1, 1, ix] = 1
    end
    DQRRQ = TensorMap(DQRRQ, ℂ^(d+1)*ℂ^(d+1), ℂ^(d+1)*ℂ^(d+1))
    DRQQR = TensorMap(DRQQR, ℂ^(d+1)*ℂ^(d+1), ℂ^(d+1)*ℂ^(d+1))

    @tensor QR_comm[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] + DR[-1, -3] * DQ[-2, -4] -
                                       DQRRQ[-1, -2, -3, -4] - DRQQR[-1, -2, -3, -4]
    @tensor QR_comm2[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] - DRQQR[-1, -2, -3, -4]
    @tensor QR_comm3[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] - DQRRQ[-1, -2, -3, -4]

    # to be used in construction of M matrices by using Q, R matrices, assuming d = 2
    Q = ψ.Q
    Iψ = id(ℂ^χ)
    R = TensorMap(convert_to_array(ψ.R)[:, 1, :], ℂ^χ, ℂ^χ)
    tangent1, tangent2 = P * Vx1, P * Vx2
    V1  = TensorMap(convert_to_array(tangent1)[:, 1, :], ℂ^χ, ℂ^χ)
    W1  = TensorMap(convert_to_array(tangent1)[:, 2, :], ℂ^χ, ℂ^χ)
    V2 = TensorMap(convert_to_array(tangent2)[:, 1, :], ℂ^χ, ℂ^χ)
    W2 = TensorMap(convert_to_array(tangent2)[:, 2, :], ℂ^χ, ℂ^χ)
    tnp = (A1, A2) -> (@tensor Aout[-1, -2; -3, -4] := A1'[-3, -1] * A2[-2, -4])
    comm = (A1, A2) -> (A1 * A2 - A2 * A1)

    # Mx's
    @tensor Mx1[-1, -2; -3, -4] := A'[-3, -1, 2] * P[-2, 2, 1] * Vx1[1, -4]
    @tensor Mx2[-1, -2; -3, -4] := P'[1, -1, 2] * Vx2'[-3, 1] * A[-2, 2, -4]
    @tensor Mxx[-1, -2; -3, -4] := P'[1, -1, 3] * Vx2'[-3, 1] * DR[3, 4] * P[-2, 4, 2] * Vx1[2, -4]
    @test Mx1 ≈ tnp(Iψ, V1) + tnp(R, W1)
    @test Mx2 ≈ tnp(V2, Iψ) + tnp(W2, R)
    @test Mxx ≈ tnp(W2, W1)

    # x != x' != y
    @tensor My[-1, -2; -3, -4] := -μ * ψtn'[-3, -1, 1] * DR[1, 2] * ψtn[-2, 2, -4] +
                                  k0 * ψtn'[1, -1, 3] * ψtn'[-3, 1, 4] * QR_comm[3, 4, 5, 6] * ψtn[-2, 5, 2] * ψtn[2, 6, -4] +
                                  c * ψtn'[1, -1, 3] * ψtn'[-3, 1, 4] * DR[3, 5] * DR[4, 6] * ψtn[-2, 5, 2] * ψtn[2, 6, -4]
    @test My ≈ -μ * tnp(R, R) +
               k0 * tnp(comm(Q, R), comm(Q, R)) +
               c * tnp(R*R, R*R)
    res1 =
    quadgk(x1 -> quadgk(
        x2 -> exp(-im*(p+q)*x2) * exp(im*p*x1) * tr(My * exp(x2*K) * Mx2 * exp((x1 - x2)*K) * Mx1 * exp((L-x1)*K) ),
        0, x1)[1], 0, L)[1] +
    quadgk(x1 -> quadgk(
        x2 -> exp(-im*(p+q)*x2) * exp(im*p*x1) * tr(My * exp(x1*K) * Mx1 * exp((x2 - x1)*K) * Mx2 * exp((L-x2)*K) ),
        x1, L)[1], 0, L)[1]

    # y != x == x'
    res1_a = quadgk(x1 -> exp(-im*q*x1) * tr(My * exp(x1*K) * Mxx * exp((L-x1)*K)), 0, L)[1]

    # y == x != x'
    @tensor My[-1, -2; -3, -4] := -μ * ψtn'[-3, -1, 2] * DR[2, 3] * P[-2, 3, 1] * Vx1[1, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * QR_comm[4, 5, 6, 7] * P[-2, 6, 1] * Vx1[1, 2] * ψtn[2, 7, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * DQ[4, 6] * DR[5, 7] * (im*p)*A[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] -
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * DRQQR[4, 5, 6, 7] * (im*p)*A[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] +
                                  c * ψtn'[1, -1, 4] * ψtn'[-3, 1, 5] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * P[3, 7, 2] * Vx1[2, -4] +
                                  c * ψtn'[1, -1, 4] * ψtn'[-3, 1, 5] * DR[4, 6] * DR[5, 7] * P[-2, 6, 2] * Vx1[2, 3] * ψtn[3, 7, -4]
    @test My ≈ -μ * tnp(R, W1) +
               k0 * tnp(comm(Q, R), comm(V1, R) + comm(Q, W1) + im*p*W1) +
               c * tnp(R*R, R*W1 + W1*R)
    res2 = quadgk(x2 -> exp(-im*(p+q)*x2) * tr(My * exp(x2*K) * Mx2 * exp((L-x2)*K)), 0, L)[1]

    # y == x' != x
    @tensor My[-1, -2; -3, -4] := -μ * P'[1, -1, 2] * Vx2'[-3, 1] * DR[2, 3] * ψtn[-2, 3, -4] +
                                  k0 * P'[1, -1, 4] * Vx2'[2, 1] * ψtn'[-3, 2, 5] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  k0 * ψtn'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DQ[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] -
                                  k0 * (-im*(p+q))*A'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DQRRQ[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  c * ψtn'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  c * P'[1, -1, 4] * Vx2'[2, 1] * ψtn'[-3, 2, 5] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4]
    @test My ≈ -μ * tnp(W2, R) +
               k0 * tnp(comm(V2, R) + comm(Q, W2) + im*(p+q)*W2, comm(Q, R)) +
               c * tnp(R*W2 + W2*R, R*R)
    res3 = quadgk(x1 -> exp(im*p*x1) * tr(My * exp(x1*K) * Mx1 * exp((L-x1)*K)), 0, L)[1]

    # y == x' == x
    @tensor My[-1, -2; -3, -4] := -μ * P'[1, -1, 3] * Vx2'[-3, 1] * DR[3, 4] * P[-2, 4, 2] * Vx1[2, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm2[5, 6, 7, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm2[5, 6, 7, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm3[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm3[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DQ[5, 7] * DR[6, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DR[5, 7] * DR[6, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * DR[5, 7] * DR[6, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DR[5, 7] * DR[6, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  c * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * DR[5, 7] * DR[6, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4]
    @test My ≈ -μ * tnp(W2, W1) +
               k0 * tnp(comm(V2, R) + comm(Q, W2) + im*(p+q)*W2, comm(V1, R) + comm(Q, W1) + im*p*W1) +
               c * tnp(R*W2 + W2*R, R*W1 + W1*R)
    res4 = tr(exp(L*K) * My)

    bbb = L*(res1 + res1_a + res2 + res3 + res4)

    @test aaa ≈ bbb
end

@timedtestset "test lieb liniger h tangent map; test for d=2" for ix in 1:10

    k0, μ, c = rand(3)
    χ, d, L = 2, 2, rand()
    p, q = 2*pi/L, 2*2*pi/L

    ψ = cmps(rand, χ, d)
    ψ = normalize(ψ, L; sym=false)
    P = gauge_fixing_proj(ψ, L)
    Vx2 = TensorMap(rand, ComplexF64, ℂ^(χ*d), ℂ^χ)
    Vx1= TensorMap(rand, ComplexF64, ℂ^(χ*d), ℂ^χ)

    # overlap by tangent map function
    lop = lieb_liniger_h_tangent_map(ψ, p, q, c, L, μ; k0=k0)
    aaa = tr(Vx2' * lop(Vx1))

    # manual integration + alternative construction of My
    ψtn = convert_to_tensormap(ψ)
    K = K_mat(ψ, ψ)

    # to be used in construction of M matrices (My, Mx1, Mx2, see below) with tensorkit
    A = cmps(id(ℂ^χ), copy(ψ.R))
    A = convert_to_tensormap(A)
    DR = id(ℂ^(d+1))
    DR.data[1] = 0
    DQ = zero(DR)
    DQ.data[1] = 1
    DQRRQ = zeros(d+1, d+1, d+1, d+1)
    DRQQR = zeros(d+1, d+1, d+1, d+1)
    for ix in 2:d+1
        DQRRQ[1, ix, ix, 1] = 1
        DRQQR[ix, 1, 1, ix] = 1
    end
    DQRRQ = TensorMap(DQRRQ, ℂ^(d+1)*ℂ^(d+1), ℂ^(d+1)*ℂ^(d+1))
    DRQQR = TensorMap(DRQQR, ℂ^(d+1)*ℂ^(d+1), ℂ^(d+1)*ℂ^(d+1))

    @tensor QR_comm[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] + DR[-1, -3] * DQ[-2, -4] -
                                       DQRRQ[-1, -2, -3, -4] - DRQQR[-1, -2, -3, -4]
    @tensor QR_comm2[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] - DRQQR[-1, -2, -3, -4]
    @tensor QR_comm3[-1, -2; -3, -4] := DQ[-1, -3] * DR[-2, -4] - DQRRQ[-1, -2, -3, -4]

    # to be used in construction of M matrices by using Q, R matrices, assuming d = 2
    Q = ψ.Q
    Iψ = id(ℂ^χ)
    Ra = TensorMap(convert_to_array(ψ.R)[:, 1, :], ℂ^χ, ℂ^χ)
    Rb = TensorMap(convert_to_array(ψ.R)[:, 2, :], ℂ^χ, ℂ^χ)
    tangent1, tangent2 = P * Vx1, P * Vx2
    V1  = TensorMap(convert_to_array(tangent1)[:, 1, :], ℂ^χ, ℂ^χ)
    Wa1  = TensorMap(convert_to_array(tangent1)[:, 2, :], ℂ^χ, ℂ^χ)
    Wb1  = TensorMap(convert_to_array(tangent1)[:, 3, :], ℂ^χ, ℂ^χ)
    V2 = TensorMap(convert_to_array(tangent2)[:, 1, :], ℂ^χ, ℂ^χ)
    Wa2 = TensorMap(convert_to_array(tangent2)[:, 2, :], ℂ^χ, ℂ^χ)
    Wb2 = TensorMap(convert_to_array(tangent2)[:, 3, :], ℂ^χ, ℂ^χ)
    tnp = (A1, A2) -> (@tensor Aout[-1, -2; -3, -4] := A1'[-3, -1] * A2[-2, -4])
    comm = (A1, A2) -> (A1 * A2 - A2 * A1)

    # Mx's
    @tensor Mx1[-1, -2; -3, -4] := A'[-3, -1, 2] * P[-2, 2, 1] * Vx1[1, -4]
    @tensor Mx2[-1, -2; -3, -4] := P'[1, -1, 2] * Vx2'[-3, 1] * A[-2, 2, -4]
    @tensor Mxx[-1, -2; -3, -4] := P'[1, -1, 3] * Vx2'[-3, 1] * DR[3, 4] * P[-2, 4, 2] * Vx1[2, -4]
    @test Mx1 ≈ tnp(Iψ, V1) + tnp(Ra, Wa1) + tnp(Rb, Wb1)
    @test Mx2 ≈ tnp(V2, Iψ) + tnp(Wa2, Ra) + tnp(Wb2, Rb)
    @test Mxx ≈ tnp(Wa2, Wa1) + tnp(Wb2, Wb1)

    # x != x' != y
    @tensor My[-1, -2; -3, -4] := -μ * ψtn'[-3, -1, 1] * DR[1, 2] * ψtn[-2, 2, -4] +
                                  k0 * ψtn'[1, -1, 3] * ψtn'[-3, 1, 4] * QR_comm[3, 4, 5, 6] * ψtn[-2, 5, 2] * ψtn[2, 6, -4] +
                                  c * ψtn'[1, -1, 3] * ψtn'[-3, 1, 4] * DR[3, 5] * DR[4, 6] * ψtn[-2, 5, 2] * ψtn[2, 6, -4]
    @test My ≈ -μ * (tnp(Ra, Ra) + tnp(Rb, Rb)) +
               k0 * (tnp(comm(Q, Ra), comm(Q, Ra)) + tnp(comm(Q, Rb), comm(Q, Rb))) +
               c * (tnp(Ra*Ra, Ra*Ra) + tnp(Rb*Rb, Rb*Rb) + tnp(Ra*Rb, Ra*Rb) + tnp(Rb*Ra, Rb*Ra))
    res1 =
    quadgk(x1 -> quadgk(
        x2 -> exp(-im*(p+q)*x2) * exp(im*p*x1) * tr(My * exp(x2*K) * Mx2 * exp((x1 - x2)*K) * Mx1 * exp((L-x1)*K) ),
        0, x1)[1], 0, L)[1] +
    quadgk(x1 -> quadgk(
        x2 -> exp(-im*(p+q)*x2) * exp(im*p*x1) * tr(My * exp(x1*K) * Mx1 * exp((x2 - x1)*K) * Mx2 * exp((L-x2)*K) ),
        x1, L)[1], 0, L)[1]

    # y != x == x'
    res1_a = quadgk(x1 -> exp(-im*q*x1) * tr(My * exp(x1*K) * Mxx * exp((L-x1)*K)), 0, L)[1]

    # y == x != x'
    @tensor My[-1, -2; -3, -4] := -μ * ψtn'[-3, -1, 2] * DR[2, 3] * P[-2, 3, 1] * Vx1[1, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * QR_comm[4, 5, 6, 7] * P[-2, 6, 1] * Vx1[1, 2] * ψtn[2, 7, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] +
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * DQ[4, 6] * DR[5, 7] * (im*p)*A[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] -
                                  k0 * ψtn'[3, -1, 4] * ψtn'[-3, 3, 5] * DRQQR[4, 5, 6, 7] * (im*p)*A[-2, 6, 2] * P[2, 7, 1] * Vx1[1, -4] +
                                  c * ψtn'[1, -1, 4] * ψtn'[-3, 1, 5] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * P[3, 7, 2] * Vx1[2, -4] +
                                  c * ψtn'[1, -1, 4] * ψtn'[-3, 1, 5] * DR[4, 6] * DR[5, 7] * P[-2, 6, 2] * Vx1[2, 3] * ψtn[3, 7, -4]
    @test My ≈ -μ * (tnp(Ra, Wa1) + tnp(Rb, Wb1)) +
               k0 * (tnp(comm(Q, Ra), comm(V1, Ra) + comm(Q, Wa1) + im*p*Wa1) + tnp(comm(Q, Rb), comm(V1, Rb) + comm(Q, Wb1) + im*p*Wb1)) +
               c * (tnp(Ra*Ra, Ra*Wa1 + Wa1*Ra) + tnp(Rb*Rb, Rb*Wb1 + Wb1*Rb) + tnp(Ra*Rb, Ra*Wb1 + Wa1*Rb) + tnp(Rb*Ra, Rb*Wa1 + Wb1*Ra))
    res2 = quadgk(x2 -> exp(-im*(p+q)*x2) * tr(My * exp(x2*K) * Mx2 * exp((L-x2)*K)), 0, L)[1]

    # y == x' != x
    @tensor My[-1, -2; -3, -4] := -μ * P'[1, -1, 2] * Vx2'[-3, 1] * DR[2, 3] * ψtn[-2, 3, -4] +
                                  k0 * P'[1, -1, 4] * Vx2'[2, 1] * ψtn'[-3, 2, 5] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  k0 * ψtn'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * QR_comm[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DQ[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] -
                                  k0 * (-im*(p+q))*A'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DQRRQ[4, 5, 6, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  c * ψtn'[2, -1, 4] * P'[1, 2, 5] * Vx2'[-3, 1] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4] +
                                  c * P'[1, -1, 4] * Vx2'[2, 1] * ψtn'[-3, 2, 5] * DR[4, 6] * DR[5, 7] * ψtn[-2, 6, 3] * ψtn[3, 7, -4]
    @test My ≈ -μ * (tnp(Wa2, Ra) + tnp(Wb2, Rb)) +
               k0 * (tnp(comm(V2, Ra) + comm(Q, Wa2) + im*(p+q)*Wa2, comm(Q, Ra)) + tnp(comm(V2, Rb) + comm(Q, Wb2) + im*(p+q)*Wb2, comm(Q, Rb))) +
               c * (tnp(Ra*Wa2 + Wa2*Ra, Ra*Ra) + tnp(Rb*Wb2 + Wb2*Rb, Rb*Rb) + tnp(Ra*Wb2 + Wa2*Rb, Ra*Rb) + tnp(Rb*Wa2 + Wb2*Ra, Rb*Ra))
    res3 = quadgk(x1 -> exp(im*p*x1) * tr(My * exp(x1*K) * Mx1 * exp((L-x1)*K)), 0, L)[1]

    # y == x' == x
    @tensor My[-1, -2; -3, -4] := -μ * P'[1, -1, 3] * Vx2'[-3, 1] * DR[3, 4] * P[-2, 4, 2] * Vx1[2, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * QR_comm2[5, 6, 7, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm2[5, 6, 7, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm3[5, 6, 7, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * QR_comm3[5, 6, 7, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  k0 * (-im*(p+q))*A'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DQ[5, 7] * DR[6, 8] * (im*p)*A[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DR[5, 7] * DR[6, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * DR[5, 7] * DR[6, 8] * ψtn[-2, 7, 4] * P[4, 8, 3] * Vx1[3, -4] +
                                  c * ψtn'[2, -1, 5] * P'[1, 2, 6] * Vx2'[-3, 1] * DR[5, 7] * DR[6, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4] +
                                  c * P'[1, -1, 5] * Vx2'[2, 1] * ψtn'[-3, 2, 6] * DR[5, 7] * DR[6, 8] * P[-2, 7, 3] * Vx1[3, 4] * ψtn[4, 8, -4]
    @test My ≈ -μ * (tnp(Wa2, Wa1) + tnp(Wb2, Wb1)) +
               k0 * (tnp(comm(V2, Ra) + comm(Q, Wa2) + im*(p+q)*Wa2, comm(V1, Ra) + comm(Q, Wa1) + im*p*Wa1) +
                     tnp(comm(V2, Rb) + comm(Q, Wb2) + im*(p+q)*Wb2, comm(V1, Rb) + comm(Q, Wb1) + im*p*Wb1)) +
               c * (tnp(Ra*Wa2 + Wa2*Ra, Ra*Wa1 + Wa1*Ra) +
                    tnp(Rb*Wb2 + Wb2*Rb, Rb*Wb1 + Wb1*Rb) +
                    tnp(Ra*Wb2 + Wa2*Rb, Ra*Wb1 + Wa1*Rb) +
                    tnp(Rb*Wa2 + Wb2*Ra, Rb*Wa1 + Wb1*Ra))
    res4 = tr(exp(L*K) * My)

    bbb = L*(res1 + res1_a + res2 + res3 + res4)

    @test aaa ≈ bbb
end


text:
@testset "HLAAllele(::AbstractString)" begin
    @test HLAAllele("HLA-B*11:10") == HLAAllele("B", "11", "10", missing, missing, missing)
    @test_throws ErrorException HLAAllele("B*11:10")
    @test_throws ErrorException HLAAllele("HLA-A*1110N")
end

@testset "HLAType(::NTuple{6, HLAAllele})" begin
    @test HLAType(parse_allele("A13", "A21", "B51", "B57", "C03", "C07")) isa HLAType
    @test_throws MethodError HLAType(parse_allele("A13", "B07"))
end

@testset "hla_matrix(::Vector{HLAType})" begin
    hla_types = [HLAType(parse_allele("A11", "A03", "B7", "B7", "C7", "C7")),
                 HLAType(parse_allele("A33", "A03", "B7", "B7", "C7", "C7"))]

    @test Escape.hla_matrix(hla_types) == [1 1 0 1 1; 1 0 1 1 1]
end

@testset "hla_matrix_standardized(::Vector{HLAType})" begin
    hla_types = [HLAType(parse_allele("A11", "A03", "B7", "B7", "C7", "C7")),
                 HLAType(parse_allele("A33", "A03", "B7", "B7", "C7", "C7"))]

    @test Escape.hla_matrix_standardized(hla_types) ≈ [0 0.5 -0.5 0 0; 0 -0.5 0.5 0 0]
end

@testset "Base.in(::HLAAllele, ::HLAType)" begin
    hla_type = HLAType(parse_allele("A13", "A21:01", "B51", "B57", "C03", "C07"))
    @test parse_allele("A13") in hla_type
    @test parse_allele("A21:01") in hla_type
    @test parse_allele("B51") ∈ hla_type
    @test parse_allele("B52") ∉ hla_type
    @test parse_allele("A21") in hla_type
end

@testset "Base.:(==)(::HLAAllele, ::HLAAllele)" begin
    @test HLAAllele("HLA-B*11:10") == HLAAllele("HLA-B*11:10")
    @test HLAAllele("HLA-B*11:10") != HLAAllele("HLA-B*11:11")
end

@testset "is_valid_allele(::HLAAllele)" begin
    @test is_valid_allele("HLA-B") == false
    @test is_valid_allele("A*11:01") == false
    @test is_valid_allele("HLA-A*11:01") == true
    @test is_valid_allele("HLA-A*11:01N") == true
end

@testset "Base.isless(::HLAAllele, ::HLAAllele)" begin
    @test parse_allele("A11") > parse_allele("A03")
    @test parse_allele("A11") < parse_allele("B11")
    @test parse_allele("A11") < parse_allele("C20")
    @test parse_allele("A03") <= parse_allele("A11")
    @test parse_allele("A11") <= parse_allele("A11")
    @test ismissing(parse_allele("A11:01") <= parse_allele("A11"))
    @test parse_allele("B51:03") > parse_allele("B51:01")
end

@testset "parse_allele(::AbstractString)" begin
    @test parse_allele("A2402") == HLAAllele("A", "24", "02", missing, missing, missing)
    @test parse_allele("A2401") != parse_allele("A2402")
    @test parse_allele("Cw*0302") == HLAAllele("C", "03", "02", missing, missing, missing)
    @test parse_allele("Cw0302") == parse_allele("HLA-C*03:02")
    @test parse_allele("Cw4") == HLAAllele("C", "04", missing, missing, missing, missing)
    @test parse_allele("A*29") == HLAAllele("A", "29", missing, missing, missing, missing)
    @test parse_allele("A021010102N") == HLAAllele("A", "02", "101", "01", "02", "N")
    @test parse_allele("A02(1010102N)") == HLAAllele("A", "02", "101", "01", "02", "N")
    @test parse_allele("A02:10N") == HLAAllele("A", "02", "10", missing, missing, "N")
    @test parse_allele("HLA-B*44:02:01:02S") == HLAAllele("B", "44", "02", "01", "02", "S")
    @test ismissing(parse_allele("?"))
    @test ismissing(parse_allele("N"))
end

@testset "parse_allele(x...)" begin
    @test length(parse_allele("A*11", "A*12")) == 2
    @test parse_allele("A*11", "A*12") == (parse_allele("A*11"), parse_allele("A*12"))
end

@testset "string(::HLAAllele)" begin
    @test string(HLAAllele("HLA-A*11:01:05N")) == "HLA-A*11:01:05N"
    @test string(parse_allele("A1")) == "HLA-A*01"
end

@testset "limit_hla_accuracy(::HLAAllele)" begin
    allele = parse_allele("A*11:01:05N")
    @test limit_hla_accuracy(allele) == parse_allele("A*11")
    @test limit_hla_accuracy(allele, allele_depth = 2) == parse_allele("A*1101")
    @test limit_hla_accuracy(parse_allele("A11")) == parse_allele("A*11")
end

@testset "rand(::HLAAllele, ::Symbol)" begin
    @test typeof(rand(HLAAllele, :A)) == HLAAllele
    @test typeof(rand(HLAAllele, :B)) == HLAAllele
    @test typeof(rand(HLAAllele, :C)) == HLAAllele
    @test_throws ErrorException rand(HLAAllele, :Z)

    @test rand(HLAAllele, :A).gene == "A"
    @test rand(HLAAllele, :B).gene == "B"
    @test rand(HLAAllele, :C).gene == "C"
end

@testset "rand(::HLAType)" begin
    @test rand(HLAType) isa HLAType
    @test length(rand(HLAType, 3)) == 3
    @test rand(HLAType, 3)[2] isa HLAType
end

@testset "unique_alleles" begin
    hla_types = [HLAType(parse_allele("A11", "A03", "B03", "B51", "C03", "C03")),
                 HLAType(parse_allele("A13", "A27", "B03", "B51", "C07", "C07"))]
    @test length(Escape.unique_alleles(hla_types)) == 8
    @test allunique(Escape.unique_alleles(hla_types))
end

@testset "hla_accuracy(::HLAAllele)" begin
    @test hla_accuracy(parse_allele("A11")) == 1
    @test hla_accuracy(parse_allele("A11:01")) == 2
    @test hla_accuracy(parse_allele("B57:01:01")) == 3
    @test hla_accuracy(parse_allele("B57:01:01:01")) == 4
    @test hla_accuracy(parse_allele("B57:01:01:01N")) == 5
end


text:
# Code gen
dir = joinpath(@__DIR__)
include(joinpath(dir, "model.jl"))

# Analytical Model
model = deepcopy(rigidbody)

path_base = joinpath(dir, "test_dynamics/base.jld2")
path_dyn = joinpath(dir, "test_dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model,
	M_analytical = true,
	mapping = G_func,
	nv = model.dim.q - 1)

save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

# Lagrnagian Model
model_lag = deepcopy(rigidbody)
path_base_lag = joinpath(dir, "test_dynamics/base_lag.jld2")
path_dyn_lag = joinpath(dir, "test_dynamics/dynamics_lag.jld2")

expr_base_lag = generate_base_expressions(model_lag,
	M_analytical = false,
	mapping = G_func,
	nv = model.dim.q - 1)

save_expressions(expr_base, path_base_lag, overwrite=true)
instantiate_base!(model_lag, path_base_lag)

expr_dyn_lag = generate_dynamics_expressions(model_lag)
save_expressions(expr_dyn_lag, path_dyn_lag, overwrite=true)
instantiate_dynamics!(model_lag, path_dyn_lag)

# Data
env = environment_3D_flat()
nq = model.dim.q
nv = model.dim.q - 1
nu = model.dim.u
nw = model.dim.w
nc = model.dim.c
h = 0.01

# initial conditions
r0 = [0.0; 0.0; 1.0]
v0 = [7.5; 5.0; 0.0]

quat0 = [1.0; 0.0; 0.0; 0.0]
ω0 = [0.0; 0.0; 0.0]

q0 = SVector{s.model.dim.q}([r0; quat0])
q1 = SVector{s.model.dim.q}([r0 + v0 * h; 0.5 * h * L_multiply(quat0) * [sqrt((2.0 / h)^2.0 - ω0' * ω0); ω0]])
q2 = SVector{s.model.dim.q}([r0 + v0 * 2h; 0.5 * 2h * L_multiply(quat0) * [sqrt((2.0 / (2h))^2.0 - ω0' * ω0); ω0]])
q̇0 = zeros(SVector{nv})
u1 = zeros(SVector{nu})
w1 = zeros(SVector{nw})
λ1 = zeros(SVector{nc * dim(env)})

# Test equality
@test norm(M_fast(model, q0) - M_fast(model_lag, q0)) < 1e-10
@test norm(C_fast(model, q0, q̇0) - C_fast(model_lag, q0, q̇0)) < 1e-10
@test norm(d_fast(model, [h], q0, q1, u1, w1, λ1, q2) -
	d_fast(model_lag, [h], q0, q1, u1, w1, λ1, q2)) < 1e-10


text:
using Test
using EllZeroTrendFiltering

for problem_fcn in ["discontinuous1",
                    "discontinuous2",
                    "OEIS1",
                    "super_exponential1",
                    "super_exponential2",
                    "super_exponential3",
                    "super_exponential4",
                    "super_exponential5",
                    "square_wave",
                    "white_noise",
                    "exponential",
                    "straight_line"]

    include(joinpath(dirname(@__FILE__),"problems", problem_fcn * ".jl"))
    local g, I_sols  # TODO, can remove in julia 1.0?
    g, ζ_vec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()

    @testset "Data set: $problem_fcn, regularization with ζ=$ζ" for ζ in ζ_vec

        I_reg, _, f_reg = fit_pwl_regularized(g, ζ)
        # (the cost f_reg includes penality mζ)

        # Use costs in the solution file to find out how many segments the solution should contain
        m_expected = argmin([f_sols[m] + m*ζ for m=1:length(f_sols)])

        @test m_expected == length(I_reg) - 1
        if !isempty(I_sols[m_expected]) # Only if there is a unique solution
            @test I_reg == I_sols[m_expected]
        end
        @test f_reg ≈ f_sols[m_expected]   atol=1e-10  # + ζ*m_expected

    end

    # Iterate only over problems where solution exists
    @testset "Data set: $problem_fcn, subsampled, regularization with ζ=$ζ" for ζ in ζ_vec
        # Use costs in the solution file to find out how many segments the solution should contain
        m_expected = argmin([f_sols[m] + m*ζ for m=1:length(f_sols)])

        # Generate max(length(g)/5, m) random gridpoints in range 1:length(g)
        randI = shuffle(1:length(g))[1:max(floor(Int,length(g)/5), m_expected+1)]
        # Merge with correct solution, and [1,length(g)] if solution doesn't exist
        t_grid = sort(union(randI, I_sols[m_expected], [1, length(g)]))
        # Solve on subsampled grid
        I_reg, _, f_reg = fit_pwl_regularized(g, ζ; t = t_grid)

        # Should be same if we had a solution
        if !isempty(I_sols[m_expected])
            @test m_expected == length(I_reg) - 1
            @test f_reg ≈ f_sols[m_expected]  atol=1e-10
            @test t_grid[I_reg] == I_sols[m_expected]
        else
            @test f_reg >= f_sols[m_expected] || isapprox(f_reg, f_sols[m_expected], atol=1e-10)
        end
    end
end


text:
""" Algorithms on graphs based on C-sets.
"""
module GraphAlgorithms
export connected_components, connected_component_projection, topological_sort,
  transitive_reduction!

using Compat: isnothing
using DataStructures: Stack

using ...Theories: dom, codom
using ...CSetDataStructures, ..BasicGraphs

# Connectivity
##############

""" (Weakly) connected components of a graph.

Returns a vector of vectors, which are the components of the graph.
"""
function connected_components(g::AbstractACSet)::Vector{Vector{Int}}
  π = connected_component_projection(g)
  components = [ Int[] for c in codom(π) ]
  for v in dom(π)
    push!(components[π(v)], v)
  end
  components
end

""" Projection onto (weakly) connected components of a graph.

Returns a function in FinSet{Int} from the vertex set to the set of components.
"""
function connected_component_projection end
# Implemented elsewhere, where coequalizers are available.

# DAGs
######

abstract type TopologicalSortAlgorithm end
struct TopologicalSortByDFS <: TopologicalSortAlgorithm end

""" Topological sort of a directed acyclic graph.

The [depth-first search](https://en.wikipedia.org/wiki/Topological_sorting#Depth-first_search)
algorithm is adapted from the function
[`topological_sort_by_dfs`](https://github.com/JuliaGraphs/LightGraphs.jl/blob/1c6cf65cc0981250e430bbef39055da23bd25bd0/src/traversals/dfs.jl#L44)
in LightGraphs.jl.
"""
function topological_sort(g::AbstractACSet;
                          alg::TopologicalSortAlgorithm=TopologicalSortByDFS())
  topological_sort(g, alg)
end

function topological_sort(g::AbstractACSet, ::TopologicalSortByDFS)
  vs = Int[]
  marking = fill(Unmarked, nv(g))
  for v in reverse(vertices(g))
    marking[v] == Unmarked || continue
    marking[v] = TempMarked
    stack = Stack{Int}()
    push!(stack, v)
    while !isempty(stack)
      u = first(stack)
      u_out = outneighbors(g, u)
      i = findfirst(u_out) do w
        marking[w] != TempMarked || error("Graph is not acyclic: $g")
        marking[w] == Unmarked
      end
      if isnothing(i)
        marking[u] = Marked
        push!(vs, u)
        pop!(stack)
      else
        w = u_out[i]
        marking[w] = TempMarked
        push!(stack, w)
      end
    end
  end
  reverse!(vs)
end

@enum TopologicalSortDFSMarking Unmarked=0 TempMarked=1 Marked=2

""" Transitive reduction of a DAG.

The algorithm computes the longest paths in the DAGs and keeps only the edges
corresponding to longest paths of length 1. Requires a topological sort, which
is computed if it is not supplied.
"""
function transitive_reduction!(g::AbstractACSet; sorted=nothing)
  lengths = longest_paths(g, sorted=sorted)
  transitive_edges = filter(edges(g)) do e
    lengths[tgt(g,e)] - lengths[src(g,e)] != 1
  end
  rem_edges!(g, transitive_edges)
  return g
end

""" Longest paths in a DAG.

Returns a vector of integers whose i-th element is the length of the longest
path to vertex i. Requires a topological sort, which is computed if it is not
supplied.
"""
function longest_paths(g::AbstractACSet;
                       sorted::Union{AbstractVector{Int},Nothing}=nothing)
  if isnothing(sorted)
    sorted = topological_sort(g)
  end
  lengths = fill(0, nv(g))
  for v in sorted
    lengths[v] = mapreduce(max, inneighbors(g, v), init=0) do u
      lengths[u] + 1
    end
  end
  lengths
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

""" Heat equations.

Field equation is:

    ρc∂u/∂t = ∇⋅(k∇u) + f

Weak form is: find u∈U such that ∀v in V

    ∫k∇u∇v dx = ∫fv dx + ∫gv ds,

where

    k = temperature thermal conductivity    defined on volume elements
    f = temperature load                    defined on volume elements
    g = temperature flux                    defined on boundary elements

Parameters
----------
temperature thermal conductivity
temperature load
temperature flux
thermal conductivity
heat source
heat flux
heat transfer coefficient
external temperature

Formulations
------------
1D, 2D, 3D

References
----------
https://en.wikipedia.org/wiki/Heat_equation
https://en.wikipedia.org/wiki/Heat_capacity
https://en.wikipedia.org/wiki/Heat_flux
https://en.wikipedia.org/wiki/Thermal_conduction
https://en.wikipedia.org/wiki/Thermal_conductivity
https://en.wikipedia.org/wiki/Thermal_diffusivity
https://en.wikipedia.org/wiki/Volumetric_heat_capacity
"""
type Heat <: FieldProblem
    formulation :: AbstractString
    store_fields :: Vector{Symbol}
end

function Heat()
    return Heat("3D", [])
end

function get_unknown_field_name(problem::Problem{Heat})
    return "temperature"
end

function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time)
    formulation = Val{Symbol(problem.properties.formulation)}
    assemble!(assembly, problem, element, time, formulation)
end

# 3d heat problems

function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
    info("Unknown element type $E for 3d heat problem!")
end

typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8, Hex20, Hex27}
typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4, Quad8, Quad9}

function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
    gdofs = get_gdofs(problem, element)
    field_name = get_unknown_field_name(problem)
    nnodes = length(element)
    K = zeros(nnodes, nnodes)
    fq = zeros(nnodes)
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        if haskey(element, "$field_name thermal conductivity")
            dN = element(ip, time, Val{:Grad})
            k = element("$field_name thermal conductivity", ip, time)
            K += w*k*dN'*dN
        end
        if haskey(element, "thermal conductivity")
            dN = element(ip, time, Val{:Grad})
            k = element("thermal conductivity", ip, time)
            K += w*k*dN'*dN
        end
        if haskey(element, "$field_name load")
            f = element("$field_name load", ip, time)
            fq += w*N'*f
        end
        if haskey(element, "heat source")
            f = element("heat source", ip, time)
            fq += w*N'*f
        end
    end
    T = vec(element[field_name](time))
    fq -= K*T
    add!(assembly.K, gdofs, gdofs, K)
    add!(assembly.f, gdofs, fq)
end

function postprocess!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time)
    haskey(element, "temperature thermal conductivity") || return
    gdofs = get_gdofs(problem, element)
    field_name = get_unknown_field_name(problem)
    nnodes = length(element)
    Me = zeros(nnodes, nnodes)
    De = zeros(nnodes, nnodes)
    f = zeros(nnodes, 3)
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        De += w*diagm(vec(N))
        Me += w*N'*N
    end
    Ae = De*inv(Me)
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = vec(element(ip, time))
        Phi = transpose(Ae*N)
        k = element("temperature thermal conductivity", ip, time)
        gradT = element(field_name, ip, time, Val{:Grad})
        q = -vec(k*gradT)
        update!(ip, "heat flux", time => q)
        f += w*Phi'*q'
    end
    add!(assembly.M, gdofs, gdofs, De)
    add!(assembly.f, gdofs, [1, 2, 3], f)
end

function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
    gdofs = get_gdofs(problem, element)
    field_name = get_unknown_field_name(problem)
    nnodes = length(element)
    K = zeros(nnodes, nnodes)
    fq = zeros(nnodes)
    for ip in get_integration_points(element, 1)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        if haskey(element, "$field_name flux")
            q = element("$field_name flux", ip, time)
            fq += w*N'*q
        end
        if haskey(element, "heat flux")
            q = element("heat flux", ip, time)
            fq += w*N'*q
        end
        if haskey(element, "$field_name heat transfer coefficient") && haskey(element, "$field_name external temperature")
            h = element("$field_name heat transfer coefficient", ip, time)
            Tu = element("$field_name external temperature", ip, time)
            K += w*h*N'*N
            fq += w*N'*h*Tu
        end
        if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
            h = element("heat transfer coefficient", ip, time)
            Tu = element("external temperature", ip, time)
            K += w*h*N'*N
            fq += w*N'*h*Tu
        end
    end
    T = vec(element[field_name](time))
    fq -= K*T
    add!(assembly.K, gdofs, gdofs, K)
    add!(assembly.f, gdofs, fq)
end

# 2d heat problems

function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
    info("Unknown element type $E for 2d heat problem!")
end

typealias Heat2DVolumeElements Union{Tri3, Tri6, Quad4}
typealias Heat2DSurfaceElements Union{Seg2, Seg3}

function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
    gdofs = get_gdofs(problem, element)
    field_name = get_unknown_field_name(problem)
    nnodes = length(element)
    K = zeros(nnodes, nnodes)
    fq = zeros(nnodes)
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        if haskey(element, "$field_name thermal conductivity")
            dN = element(ip, time, Val{:Grad})
            k = element("$field_name thermal conductivity", ip, time)
            K += w*k*dN'*dN
        end
        if haskey(element, "thermal conductivity")
            dN = element(ip, time, Val{:Grad})
            k = element("thermal conductivity", ip, time)
            K += w*k*dN'*dN
        end
        if haskey(element, "$field_name load")
            f = element("$field_name load", ip, time)
            fq += w*N'*f
        end
        if haskey(element, "heat source")
            f = element("heat source", ip, time)
            fq += w*N'*f
        end
    end
    T = vec(element[field_name](time))
    fq -= K*T
    add!(assembly.K, gdofs, gdofs, K)
    add!(assembly.f, gdofs, fq)
end

function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
    gdofs = get_gdofs(problem, element)
    field_name = get_unknown_field_name(problem)
    nnodes = length(element)
    K = zeros(nnodes, nnodes)
    fq = zeros(nnodes)
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        if haskey(element, "$field_name flux")
            g = element("$field_name flux", ip, time)
            fq += w*N'*g
        end
        if haskey(element, "heat flux")
            g = element("heat flux", ip, time)
            fq += w*N'*g
        end
        if haskey(element, "$field_name heat transfer coefficient") && haskey(element, "$field_name external temperature")
            h = element("$field_name heat transfer coefficient", ip, time)
            Tu = element("$field_name external temperature", ip, time)
            K += w*h*N'*N
            fq += w*N'*h*Tu
        end
        if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
            h = element("heat transfer coefficient", ip, time)
            Tu = element("external temperature", ip, time)
            K += w*h*N'*N
            fq += w*N'*h*Tu
        end
    end
    T = vec(element[field_name](time))
    fq -= K*T
    add!(assembly.K, gdofs, gdofs, K)
    add!(assembly.f, gdofs, fq)
end


text:
println("Im Julia")

println("This is Fibonacci")
function fib(n)
    x,y = (0,1)
    for i = 1:n x,y = (y, x+y) end
    x
end

let x = 0
    while x < 100
        println(fib(x))
        x += 1
    end
end

max_stars_repo_name: millecodex/decarepo

max_issues_repo_name: millecodex/decarepo

max_forks_repo_name: millecodex/decarepo


text:
@testset "Permutation generation (example)" begin

clauses = @julog [
    permutations(T) <<=
        unifies(T, [_, _ ,_]) &
        subset([1,2,3], T),

    member(X, [X | Y]) <<= true,
    member(X, [Y | YS]) <<= member(X, YS),
    subset([], _) <<= true,
    subset([X | XS], L) <<= member(X,L) & subset(XS,L),
]

sat, subst = resolve(@julog(permutations(T)),clauses)
@test @varsub({T => [1, 2, 3]}) in subst
@test @varsub({T => [1, 3, 2]}) in subst
@test @varsub({T => [2, 1, 3]}) in subst
@test @varsub({T => [2, 3, 1]}) in subst
@test @varsub({T => [3, 1, 2]}) in subst
@test @varsub({T => [3, 2, 1]}) in subst

end


text:
using Test
using Pkg
using Pkg.Artifacts
using ArtifactHelpers

@testset "ArtifactHelpers" begin

    artifacts_toml = touch(joinpath(@__DIR__, "Artifacts.toml"))

    artifacts_dir = mkpath(joinpath(@__DIR__, "tmp_artifacts"))

    Pkg.Artifacts.ARTIFACTS_DIR_OVERRIDE[] = artifacts_dir

    function test_process(entry::Processed; force::Bool = false, verbose::Bool = false)

        tree_hash = create_artifact() do path_artifact #Note: this will create an artifact that is ready for use.
            write(joinpath(path_artifact, "test1.txt"), "test1")
            write(joinpath(path_artifact, "test2.txt"), "test2")
        end

        return tree_hash
    end

    @testset "Helpers" begin
        @test ArtifactHelpers.isurl("http") == true
        @test ArtifactHelpers.isurl("https") == true
        @test ArtifactHelpers.isurl("https://") == true
        @test ArtifactHelpers.isurl("htt") == false

        #Check setindex!.
        entry = Processed("test")
        change = setindex!(entry, "test", "test")
        @test ArtifactHelpers.metadata(entry) == change

    end #testset "Helpers"

    @testset "Building" begin
        @test_nowarn build_artifact!(artifacts_toml, File("http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/mm10.chrom.sizes"), force = true, verbose = false)#TODO: host or find reasonable download.

        @test_nowarn build_artifact!(artifacts_toml, Zip("https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip"), force = true, verbose = false) #TODO: host or find reasonable download.

        @test_nowarn build_artifact!(artifacts_toml, AutoDownloadable("http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromAgp.tar.gz"), force = true, verbose = false)#TODO: host or find reasonable download.

        @test_nowarn build_artifact!(artifacts_toml, Processed("Processed"), test_process, force = true, verbose = false)

        @test_nowarn build_artifact!(artifacts_toml, Processed("ProcessedDo"), force = true, verbose = false) do path_artifact
            write(joinpath(path_artifact, "test.txt"), "test")
        end

    end #testset "Binding"

    @testset "Initialise with artifacts present" begin

        @test initialise_artifact(artifacts_toml, "mm10.chrom.sizes") == artifact_hash("mm10.chrom.sizes", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "Trimmomatic-0.39.zip") == artifact_hash("Trimmomatic-0.39.zip", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "chromAgp.tar.gz") == artifact_hash("chromAgp.tar.gz", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "Processed", test_process) == artifact_hash("Processed", artifacts_toml)

    end #testest "Initialise with artifacts present"

    @testset "Initialise without artifacts present" begin
        artifacts = joinpath.(artifacts_dir, readdir(artifacts_dir))
        rm.(artifacts,recursive = true)

        @test initialise_artifact(artifacts_toml, "mm10.chrom.sizes") == artifact_hash("mm10.chrom.sizes", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "Trimmomatic-0.39.zip") == artifact_hash("Trimmomatic-0.39.zip", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "chromAgp.tar.gz") == artifact_hash("chromAgp.tar.gz", artifacts_toml)

        @test initialise_artifact(artifacts_toml, "Processed", test_process) == artifact_hash("Processed", artifacts_toml)

    end #testset "Initialise without artifacts present"

    @testset "Initialise without Artifacts.toml present" begin
        rm(artifacts_toml)
        rm(artifacts_dir, recursive = true) #Note: may aswell clean up.

        # Check fallthrough.

        @test_logs (:warn, "Artifacts.toml does not exist at specified path.") initialise_artifact(artifacts_toml, "something")

        touch(artifacts_toml)

        @test_logs (:warn, "An artifact entry for \"something\" does not exist in Artifacts.toml.") initialise_artifact(artifacts_toml, "something")

        rm(artifacts_toml) # Final clean up.

    end #testset "Initialise without artifacts present"

end # testset ArtifactHelpers


text:
# given a pre-configured boottest linear model with one-degree null imposed, compute boostrapped p value associated with r
function r_to_p(o::StrBootTest{T}, r::AbstractVector{T}) where T
  o.r = r
  getp(o)
end

# Chandrupatla 1997, "A new hybrid quadratic-bisection algorithm for finding the zero of a nonlinear function without using derivatives"
# x₁, x₂ must bracket the true value, with f₁=f(x₁) and f₂=f(x₂)
function search(o::StrBootTest{T}, α::T, f₁::T, x₁::T, f₂::T, x₂::T) where T<:Real
  t = half = T(.5)
  while true
		x = x₁ + t * (x₂ - x₁)
		fx = r_to_p(o, [x])
		((fx>f₁) == (fx>f₂)) && return x  # violation of monotonicity because of precision problems? That's as good as it gets.

		if (fx<α) == (f₁<α)
			x₃, x₁, f₃, f₁ = x₁, x, f₁, fx
		else
			x₃, x₂, x₁, f₃, f₂, f₁ = x₂, x₁, x, f₂, f₁, fx
		end

		((o.B>0 && abs(fx - α) < (1+(o.ptype == :equaltail)) / o.BFeas * 1.000001) || ≈(x₂, x₁, rtol=o.rtol)) &&
			return abs(f₁ - α) < abs(f₂ - α) ? x₁ : x₂

		ϕ₁ = (f₁ - f₂) / (f₃ - f₂)
		ϕ₁² = ϕ₁^2
		xi1 = (x₁ - x₂) / (x₃ - x₂)
		t = ϕ₁² > xi1 || xi1 > 2 * ϕ₁ - ϕ₁² ?
					half :
					clamp(((f₃ - α) / (f₁ - f₂) + (x₃ - x₁) / ((x₂ - x₁) * (f₃ - f₁)) * (f₂ - α)) * (f₁ - α) / (f₃ - f₂), T(0.000001), T(0.999999))
  end
end

# derive wild bootstrap-based CI, for case of linear model with one-degree null imposed
# and generate plot data
function plot!(o::StrBootTest{T}) where T
  _r = copy(o.r)
  α = one(T) - o.level

	iszero(nrows(o.gridmin   )) && (o.gridmin    = fill(T(NaN), o.q))
	iszero(nrows(o.gridmax   )) && (o.gridmax    = fill(T(NaN), o.q))
	iszero(nrows(o.gridpoints)) && (o.gridpoints = fill(Float32(NaN), o.q))

	o.gridpoints[isnan.(o.gridpoints) .| isinf.(o.gridpoints)] .= 25

  o.boottest!(o)

	Phi = quantile(Normal{T}(zero(T),one(T)), α/2)
  if o.arubin
		p = o.dist[1] * o.multiplier
		p = ccdf(Chisq{T}(T(o.dof)), o.sqrt ? p^2 : p)
		halfwidth = abs.(o.confpeak) * quantile(Normal{T}(zero(T),one(T)), p/2) / Phi
  else
		halfwidth = T(-1.5) * Phi .* sqrtNaN.(diag(getV(o)))
		any(isnan.(halfwidth)) && return
		o.confpeak = getb(o) + o.r
  end

	if isone(o.q)  # 1D plot
		α≤0 && (α = T(.05))  # if level=100, no ci constructed, but we need a reasonable α to choose graphing bounds
		p_lo = p_hi = T(NaN)
		if isnan(o.gridmin[1]) || isnan(o.gridmax[1])
			if o.B>0  # initial guess based on classical distribution
				lo = isnan(o.gridmin[1]) ? o.confpeak - halfwidth : o.gridmin
				hi = isnan(o.gridmax[1]) ? o.confpeak + halfwidth : o.gridmax
			else
				tmp = [sqrtNaN.(o.statDenom[1,1])] * (o.small ? cquantile(TDist{T}(o.dof_r), α/2) : cquantile(Normal{T}(zero(T),one(T)), α/2))
				lo = isnan(o.gridmin[1]) ? o.confpeak - tmp : o.gridmin
				hi = isnan(o.gridmax[1]) ? o.confpeak + tmp : o.gridmax
				if o.scorebs && !o.null && !o.willplot  # if doing simple Wald test with no graph, we're done
					o.ci = [lo hi]
					return
				end
			end

			if abs(lo[1] - o.r[1]) > abs(hi[1] - o.r[1])  # brute force way to ensure that first trial bound tested is the farther one from r, for better interpolation
				if isnan(o.gridmin[1]) && o.ptype ≠ :lower  # unless lower-tailed p value, try at most 10 times to bracket confidence set by symmetrically widening
					for _ ∈ 1:10
						p_lo = r_to_p(o, lo)
						p_lo < α && break
						diff = hi - lo
						lo .-= diff
						isnan(o.gridmax[1]) && o.twotailed && (hi .+= diff)  # maintain rough symmetry unless user specified upper bound
					end
				end
				if isnan(o.gridmax[1]) && o.ptype ≠ :upper  # ditto for high side
					for _ ∈ 1:10
						p_hi = r_to_p(o, hi)
						p_hi < α && break
						diff = hi - lo
						isnan(o.gridmin[1]) && o.twotailed && (lo .-= diff)
						hi .+= diff
					end
				end
			else
				if isnan(o.gridmax[1]) && o.ptype ≠ :upper  # ditto for high side
					for _ ∈ 1:10
						p_hi = r_to_p(o, hi)
						p_hi < α && break
						diff = hi - lo
						isnan(o.gridmin[1]) && o.twotailed && (lo .-= diff)
						hi .+= diff
					end
				end
				if isnan(o.gridmin[1]) && o.ptype ≠ :lower  # unless upper-tailed p value, try at most 10 times to bracket confidence set by symmetrically widening
					for _ ∈ 1:10
						p_lo = r_to_p(o, lo)
						p_lo < α && break
						diff = hi - lo
						lo .-= diff
						isnan(o.gridmax[1]) && o.twotailed && (hi .+= diff)  # maintain rough symmetry unless user specified upper bound
					end
				end
			end
		else  # both grid bounds pre-specified
			lo = [o.gridmin[1]]
			hi = [o.gridmax[1]]
		end

		_gridpoints = round.(Int32, o.gridpoints)
		o.plotX = [collect(range(lo[1], hi[1], length=_gridpoints[1]))]  # store as 1-vector since q=1
		o.plotY = fill(T(NaN), _gridpoints[1])
		o.plotY[1  ] = p_lo
		o.plotY[end] = p_hi
		p_confpeak = o.WRE ? T(NaN) : o.twotailed ? one(T) : T(.5)

		c = clamp((floor(Int, (o.confpeak[1] - lo[1]) / (hi[1] - lo[1]) * (_gridpoints[1] - 1)) + 2), 1, _gridpoints[1]+1)  # insert original point estimate into grid
		insert!(o.plotX[1], c, o.confpeak[1])
		insert!(o.plotY   , c, p_confpeak)

		isnan(o.plotY[1]) && (o.plotY[1] = r_to_p(o, [o.plotX[1][1]]))  # do in this order for widest interpolation base range
		for (i,v) ∈ Iterators.reverse(enumerate(o.plotX[1]))
			i>1 && isnan(o.plotY[i]) && (o.plotY[i] = r_to_p(o, [v]))
		end

	else  # 2D plot
		_gridpoints = round.(Int32, o.gridpoints)

		lo = Vector{T}(undef, 2)
		hi = Vector{T}(undef, 2)
		for d ∈ 1:o.dof
			lo[d] = isnan(o.gridmin[d]) ? o.confpeak[d] - halfwidth[d] : o.gridmin[d]
			hi[d] = isnan(o.gridmax[d]) ? o.confpeak[d] + halfwidth[d] : o.gridmax[d]
		end
		o.plotX = [collect(range(lo[1], hi[1], length=_gridpoints[1])),
		           collect(range(lo[2], hi[2], length=_gridpoints[2]))]
		o.plotY = fill(T(NaN), _gridpoints[1]*_gridpoints[2])

		isnan(o.plotY[1]) && (o.plotY[1] = r_to_p(o, [o.plotX[i][1] for i in 1:o.q]))  # do in this order for widest interpolation base range
		for (i,v) ∈ Iterators.reverse(enumerate(Iterators.product(o.plotX[1],o.plotX[2])))
			i>1 && isnan(o.plotY[i]) && (o.plotY[i] = r_to_p(o, [v...]))
		end
	end

	if any(isnan.(o.plotY))
		o.ci = [T(-Inf) T(Inf)]
	elseif isone(o.q) && o.level<100 # find ci bounds
		_CI = Vector{T}(undef, nrows(o.plotY))
		for i in eachindex(_CI)  # map() version hampers type inference in Julia 1.6.2
			_CI[i] = isnan(o.plotY[i]) ? o.plotY[i] : T(o.plotY[i] > α)
		end
		_CI = _CI[2:end] - _CI[1:end-1]
		lo = T.(findall(x->x== 1, _CI))
		hi = T.(findall(x->x==-1, _CI))
		if iszero(length(lo)) && iszero(length(hi))
			o.ci = [T(-Inf) T(Inf)]
		else
			if iszero(length(lo))
				lo = [T(-Inf)]
			elseif iszero(length(hi))
				hi = [T(Inf)]
			else
				lo[1  ] > hi[1  ] && (lo = [T(-Inf) ; lo    ]) # non-rejection ranges that are not within grid range
				lo[end] > hi[end] && (hi = [hi      ; T(Inf)])
			end
			o.ci = [lo hi]

			for i ∈ 1:length(lo), j ∈ 1:2
				if !isinf(o.ci[i,j])
					t = Int(o.ci[i,j])
					o.ci[i,j] = search(o, α, o.plotY[t], o.plotX[1][t], o.plotY[t+1], o.plotX[1][t+1])
				end
			end
		end
		nothing
	end

  if @isdefined c  # now that it's done helping 1-D graph look good, remove peak point from returned grid for evenness, for Bayesian sampling purposes
		o.peak = (X = [o.plotX[1][c]], p = o.plotY[c])
		deleteat!(o.plotX[1], c)
		deleteat!(o.plotY   , c)
  end

	o.r = _r  # restore backup
	nothing
end


text:
@require Mux begin

using Mux

Mux.Response(o::AbstractWidget) = Mux.Response(WebIO.render(o))

end


text:
# include("../src/LightGraphs.jl")
using LightGraphs
using Base.Test

g1 = smallgraph(:petersen)
g2 = smallgraph(:tutte)
g3 = PathGraph(5)
g4 = PathDiGraph(5)
g5 = DiGraph(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
g6 = smallgraph(:house)

h1 = Graph(5)
h2 = Graph(3)
h3 = Graph()
h4 = DiGraph(7)
h5 = DiGraph()

# self loops
s2 = DiGraph(3)
add_edge!(s2,1,2); add_edge!(s2,2,3); add_edge!(s2,3,3)
s1 = Graph(s2)

r1 = Graph(10,20)
r2 = DiGraph(5,10)

e0 = Edge(2, 3)
e1 = Edge(1, 2)
re1 = Edge(2, 1)

# polygons
triangle = random_regular_graph(3,2)
quadrangle = random_regular_graph(4,2)
pentagon = random_regular_graph(5,2)

testdir = dirname(@__FILE__)

pdict = load(joinpath(testdir,"testdata","tutte-pathdigraph.jgz"))
p1 = pdict["Tutte"]
p2 = pdict["pathdigraph"]

adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
distmx1 = [Inf 2.0 Inf; 2.0 Inf 4.2; Inf 4.2 Inf]
distmx2 = [Inf 2.0 Inf; 3.2 Inf 4.2; 5.5 6.1 Inf]

a1 = Graph(adjmx1)
a2 = DiGraph(adjmx2)

tests = [
    "core",
    "edgeiter",
    "operators",
    "graphdigraph",
    "persistence",
    "distance",
    "edit_distance",
    "spectral",
    "cliques",
    "subgraphs",
    "connectivity",
    "generators/randgraphs",
    "generators/staticgraphs",
    "generators/smallgraphs",
    "generators/euclideangraphs",
    "shortestpaths/astar",
    "shortestpaths/bellman-ford",
    "shortestpaths/dijkstra",
    "shortestpaths/floyd-warshall",
    "traversals/bfs",
    "traversals/dfs",
    "traversals/maxadjvisit",
    "traversals/graphvisit",
    "traversals/randomwalks",
    "community/core-periphery",
    "community/label_propagation",
    "community/modularity",
    "community/clustering",
    "centrality/betweenness",
    "centrality/closeness",
    "centrality/degree",
    "centrality/katz",
    "centrality/pagerank",
    "flow/edmonds_karp",
    "flow/dinic",
    "flow/boykov_kolmogorov",
    "flow/push_relabel",
    "flow/maximum_flow",
    "flow/multiroute_flow",
    "utils"
]

for t in tests
    tp = joinpath(testdir,"$(t).jl")
    println("running $(tp) ...")
    include(tp)
end


text:
using Documenter
using Yao
using Yao: YaoBlocks, YaoArrayRegister, YaoBase, YaoSym
using YaoBase: BitBasis
using YaoBlocks: AD
using YaoBlocks: Optimise
using Documenter.Writers.HTMLWriter
using Documenter.Utilities.DOM
using Documenter.Utilities.DOM: Tag, @tags
#Venerable Inventor :)

download("yaoquantum.org/assets/logo-light.png", "docs/src/assets/logo.png")

const PAGES = [
    "Home" => "index.md",
    "Manual" => Any[
        "man/array_registers.md",
        "man/symbolic.md",
        "man/blocks.md",
        "man/automatic_differentiation.md",
        "man/simplification.md",
        "man/base.md",
        "man/registers.md",
        "man/bitbasis.md",
        "man/extending_blocks.md",
    ],
    "Benchmark" => "benchmarks.md",
    "Developer Notes" => "dev/index.md",
]

makedocs(
    modules = [Yao, YaoBase, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, AD, Optimise],
    format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = ("deploy" in ARGS) ? "https://docs.yaoquantum.org/" : nothing,
        assets = [
            "assets/themes/indigo.css",
            asset("https://yaoquantum.org/assets/favicon-light.ico", class = :ico),
        ],
    ),
    doctest = ("doctest=true" in ARGS),
    clean = false,
    sitename = "Documentation | Yao",
    linkcheck = !("skiplinks" in ARGS),
    pages = PAGES,
)

x = []
for (root, dirs, files) in walkdir("docs/build")
           global x = [x; joinpath.(root, files)] # files is a Vector{String}, can be empty
end

for i in x
	if(endswith(i, ".html"))
		y = read(i, String)
		y = replace(y, """
""" => """

""")
		y = replace(y, """""" => """<link href="https://cdnjs.cloudflare.com/ajax/libs/tocbot/4.11.1/tocbot.css" rel="stylesheet" type="text/css"/>""")
		y = replace(y, """""" => """""")

		y = replace(y, """
""" => """
""")
		f = open(i, "w")
		write(f, y)
		close(f)
	end
end

deploydocs(repo = "github.com/VarLad/Yao.jl.git", target = "build")


text:
"""
    read(A::AbstractRaster)
    read(A::AbstractRasterStack)
    read(A::AbstractRasterSeries)

`read` will move a Rasters.jl object completely to memory.
"""
function Base.read(x::Union{AbstractRaster,AbstractRasterStack,AbstractRasterSeries})
    modify(x) do ds
        # Some backends don't implement `Array` properly.
        Array{eltype(ds),ndims(ds)}(undef, size(ds)) .= ds
    end
end

"""
    read!(src::Union{AbstractString,AbstractRaster}, dst::AbstractRaster)
    read!(src::Union{AbstractString,AbstractRasterStack}, dst::AbstractRasterStack)
    read!(scr::AbstractRasterSeries, dst::AbstractRasterSeries)

`read!` will copy the data from `src` to the object `dst`.

`src` can be an object or a file-path `String`.
"""
Base.read!(src::AbstractRaster, dst::AbstractArray) = dst .= src
function Base.read!(src::AbstractRasterStack, dst::AbstractRasterStack)
    map(k -> read!(src[k], dst[k]), keys(dst))
    return dst
end
function Base.read!(src::AbstractRasterSeries, dst::AbstractRasterSeries)
    map(read!, src, dst)
    return dst
end

# Filename methods
function Base.read!(filename::AbstractString, dst::AbstractRaster)
    src = Raster(filename;
        dims=dims(dst), refdims=refdims(dst), name=name(dst),
        metadata=metadata(dst), missingval=missingval(dst),
    )
    read!(src, dst)
end
function Base.read!(filenames::Union{NamedTuple,<:AbstractVector{<:AbstractString}}, dst::AbstractRasterStack)
    _readstack!(filenames, dst)
end
function Base.read!(filenames::AbstractString, dst::AbstractRasterStack)
    _readstack!(filenames, dst)
end
function Base.read!(filenames::AbstractVector{<:Union{AbstractString,NamedTuple}}, dst::AbstractRasterSeries)
    map((fn, d) -> read!(fn, d), filenames, dst)
    return dst
end

function _readstack!(filenames, dst)
    src = RasterStack(filenames;
        dims=dims(dst), refdims=refdims(dst), keys=keys(dst), metadata=metadata(dst),
        layermetadata=DD.layermetadata(dst), missingval=missingval(dst),
    )
    read!(src, dst)
end