text:
#=
In this file we define how agents are plotted and how the plots are updated while stepping
=#
struct ABMStepper{X, C, M, S, O, AC, AS, AM}
    ac::C
    am::M
    as::S
    offset::O
    scheduler::X
    pos::Observable
    colors::AC
    sizes::AS
    markers::AM
end

Base.show(io::IO, ::ABMStepper) =
println(io, "Helper structure for stepping and updating the plot of an agent based model. ",
"It is outputted by `abm_plot` and can be used in `Agents.step!`, see `abm_plot`.")

"Initialize the abmstepper and the plotted observables. return the stepper"
function abm_init_stepper_and_plot!(ax, model;
        ac = JULIADYNAMICS_COLORS[1],
        as = 10,
        am = :circle,
        scheduler = model.scheduler,
        offset = nothing,
        equalaspect = true,
        scatterkwargs = NamedTuple(),
        static_preplot! = default_static_preplot
    )

    o, e = modellims(model) # TODO: extend to 3D
    @assert length(o) == 2 "At the moment only 2D spaces can be plotted."
    # TODO: once graph plotting is possible, this will be adjusted
    @assert typeof(model.space) <: Union{Agents.ContinuousSpace, Agents.DiscreteSpace}
    # TODO: Point2f0 must be replaced by 3D version in the future

    static_preplot!(ax, model)

    ids = scheduler(model)
    colors  = ac isa Function ? Observable(to_color.([ac(model[i]) for i ∈ ids])) : to_color(ac)
    sizes   = as isa Function ? Observable([as(model[i]) for i ∈ ids]) : as
    markers = am isa Function ? Observable([am(model[i]) for i ∈ ids]) : am
    if offset == nothing
        pos = Observable(Point2f0[model[i].pos for i ∈ ids])
    else
        pos = Observable(Point2f0[model[i].pos .+ offset(model[i]) for i ∈ ids])
    end

    # Here we make the decision of whether the user has provided markers, and thus use
    # `scatter`, or polygons, and thus use `poly`:
    if user_used_polygons(am, markers)
        # For polygons we always need vector, even if all agents are same polygon
        if markers isa Observable
            markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
        else
            markers = Observable([translate(am, p) for p in pos])
        end
        poly!(ax, markers; color = colors, scatterkwargs...)
    else
        scatter!(
            ax, pos;
            color = colors, markersize = sizes, marker = markers, strokewidth = 0.0,
            scatterkwargs...
        )
    end
    # TODO: This should be expanded into 3D (and also scale and stuff)
    xlims!(ax, o[1], e[1])
    ylims!(ax, o[2], e[2])
    equalaspect && (ax.aspect = AxisAspect(1))

    return ABMStepper(ac, am, as, offset, scheduler, pos, colors, sizes, markers)
end

default_static_preplot(ax, model) = nothing

function modellims(model)
    if model.space isa Agents.ContinuousSpace
        e = model.space.extent
    elseif model.space isa Agents.DiscreteSpace
        e = size(model.space.s) .+ 1
    end
    return zero.(e), e
end

function user_used_polygons(am, markers)
    if (am isa Polygon)
        return true
    elseif (am isa Function) && (markers[][1] isa Polygon)
        return true
    else
        return false
    end
end

#=
    Agents.step!(abmstepper, model, agent_step!, model_step!, n::Int)
Step the given `model` for `n` steps while also updating the plot that corresponds to it,
which is produced with the function [`abm_plot`](@ref).

You can still call this function with `n=0` to update the plot for a new `model`,
without doing any stepping.
=#
function Agents.step!(abmstepper::ABMStepper, model, agent_step!, model_step!, n::Int)
    ac, am, as = abmstepper.ac, abmstepper.am, abmstepper.as
    offset = abmstepper.offset
    pos, colors = abmstepper.pos, abmstepper.colors
    sizes, markers =  abmstepper.sizes, abmstepper.markers

    Agents.step!(model, agent_step!, model_step!, n)
    if Agents.nagents(model) == 0
        @warn "The model has no agents, we can't plot anymore!"
        error("The model has no agents, we can't plot anymore!")
    end
    ids = abmstepper.scheduler(model)
    if offset == nothing
        pos[] = [model[i].pos for i in ids]
    else
        pos[] = [model[i].pos .+ offset(model[i]) for i in ids]
    end
    if ac isa Function; colors[] = to_color.([ac(model[i]) for i in ids]); end
    if as isa Function; sizes[] = [as(model[i]) for i in ids]; end
    if am isa Function; markers[] = [am(model[i]) for i in ids]; end
    # If we use Polygons as markers, do a final update:
    if user_used_polygons(am, markers)
        # translate all polygons according to pos
        markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
    end
    return nothing
end


text:
# Simulate an economy of F firms in J locations with varying levels of moving costs for workers to move firms
# Match the hypothetical parameter values to observed flows in the SIPP

@views @inbounds function datagenAmen(ew,sigw,siga,covwa,J=4,flag=true)
    Random.seed!(32)
    # this function generates wages and amenities for firms
    s0 = (1/J)*ones(J) # initial worker distribution
    Σ = [sigw covwa; covwa siga] # covariance between wages and amenities (should be neg. bec. compens. differentials)
    d = MvNormal(zeros(2),Σ)
    r = rand(d,J);
    η = round.(r[1,:]; digits=1)
    α = r[2,:]
    w = ew .+ η
    #w    = ew .+ ones(N)*ones(1,J)
    if flag==false
        α = 0 .* α
    end
    return s0,w,α,J
end

@views @inbounds function getθmkt(J,numMkt,th1,th2)
    #define markets
    @assert mod(J,numMkt)==0 "number of firms must be a multiple of the number of markets"
    firmpermkt = convert(Int64,J/numMkt)
    mktlist = []
    for m=1:numMkt
        append!(mktlist,m*ones(firmpermkt))
    end
    mktlist = convert(Array{Float64,1},mktlist)

    mktblkdiag = zeros(J,J)
    for m=1:numMkt
        mktblkdiag[findall(mktlist.==m),findall(mktlist.==m)] .= 1
    end

    # θ1: market moving cost
    θ1 = th1*(ones(J,J) - mktblkdiag)

    # θ2: firm switching cost
    θ2 = th2*(ones(J,J) - I)

    return θ1,θ2,mktlist,mktblkdiag
end

@views @inbounds function getdynshares(α,γ,θ1,θ2,w,J,β,T,ρ)
    #===
    α: J-vector of amenities
    γ: scalar for wage preference
    θ1: JxJ matrix of market moving costs
    θ2: JxJ matrix of firm switching costs
    w: J-vector of wages
    J: number of firms
    β: workers' discount factor
    T: time horizon
    ===#
    TransMat = zeros(J,J)
    if β==0
        for j=1:J # previous firm
            u = α .+ γ*w .+ θ1[j,:] .+ θ2[j,:]
            TransMat[j,:] .= exp.(u) ./ sum(exp.(u)) # Pr(k in t | j in t-1)
        end
    else
        v  = zeros(J)
        FV = zeros(T+1,J)
        for t=T:-1:1
            for k=1:J # previous firm
                v = α .+ γ*ρ^(t-1)*w .+ θ1[k,:] .+ θ2[k,:] .+ FV[t+1,:]
                if t==1
                    for j=1:J
                        v = α .+ γ*ρ^(t-1)*w .+ θ1[j,:] .+ θ2[j,:] .+ FV[t+1,:]
                        TransMat[j,:] .= exp.(v) ./ sum(exp.(v))
                    end
                end
                FV[t,k] = β*log(sum(exp.(v)))
            end
        end
    end
    return TransMat
end

@views @inbounds function get_εqw(step_size,TM0,α,γ,θ1,θ2,w,J,β,T,ρ)
    TM   = zeros(J,J)
    rec0 = zeros(J)
    rec1 = zeros(J)
    sep0 = zeros(J)
    sep1 = zeros(J)
    dw   = zeros(J)
    if β==0
        for k=1:J
            rec0[k] = sum(TM0[:,k])-TM0[k,k]
            sep0[k] = sum(TM0[k,:])-TM0[k,k]
            w1 = deepcopy(w)
            w1[k] += step_size
            dw[k] = w1[k]-w[k]
            for j=1:J
                u = α .+ γ*w1 .+ θ1[j,:] .+ θ2[j,:]
                TM[j,:] .= exp.(u) ./ sum(exp.(u))
            end
            rec1[k] = sum(TM[:,k])-TM[k,k]
            sep1[k] = sum(TM[k,:])-TM[k,k]
        end
    else
        v  = zeros(J)
        FV = zeros(T+1,J)
        for t=T:-1:1
            for k=1:J # previous firm
                v = α .+ γ*ρ^(t-1)*w .+ θ1[k,:] .+ θ2[k,:] .+ FV[t+1,:]
                if t==1
                    for m=1:J
                        rec0[m] = sum(TM0[:,m])-TM0[m,m]
                        sep0[m] = sum(TM0[m,:])-TM0[m,m]
                        w1 = deepcopy(w)
                        w1[m] += step_size
                        dw[m] = w1[m]-w[m]
                        for j=1:J
                            vv = α .+ γ*w1 .+ θ1[j,:] .+ θ2[j,:] .+ FV[t+1,:]
                            TM[j,:] .= exp.(vv) ./ sum(exp.(vv))
                        end
                        rec1[m] = sum(TM[:,m])-TM[m,m]
                        sep1[m] = sum(TM[m,:])-TM[m,m]
                    end
                end
                FV[t,k] = β*log(sum(exp.(v)))
            end
        end
    end
    εqw  = (sep1 .- sep0)./dw
    εrw  = (rec1 .- rec0)./dw
    return εqw,εrw
end

@views @inbounds function wrapper(exw,sgw,sga,cvwa,J,M,amenflag,wagegamma,movecost,switchcost,wagestep,discfac=0,timehorizon=40,ARcoef=.74)
    to = TimerOutput()
    @timeit to "elaps" begin
        s0,wage,amen,J = datagenAmen(exw,sgw,sga,cvwa,J,amenflag)
        θ_1,θ_2,MM,MM2 = getθmkt(J,M,movecost,switchcost)
        Tmat = getdynshares(amen,wagegamma,θ_1,θ_2,wage,J,discfac,timehorizon,ARcoef)
        quit_elast,rec_elast = get_εqw(wagestep,Tmat,amen,wagegamma,θ_1,θ_2,wage,J,discfac,timehorizon,ARcoef)
        wage_elast = rec_elast - quit_elast

        #migration rate
        migrat = sum(Tmat.*(1 .- MM2);dims=2) |> mean
        # should match .03411 in SIPP

        #job turnover rate
        switchrat = sum(Tmat.*(ones(size(Tmat)) - I);dims=2) |> mean
        # should match .2111054 in SIPP

        #get steady state of Markov Chain:
        capT = 20
        s = zeros(capT,J)
        for t=1:capT
            s[t,:]=s0'*Tmat^(t-1)
        end
    end
    result = DataFrame(NF = J, NM = M, avg_w = exw, sigma_w = sgw, sigma_amen = sga, cov_w_a = cvwa, gamma = convert(Float64,wagegamma), θ_move = convert(Float64,movecost), θ_switch = convert(Float64,switchcost), wagestep = wagestep, avgquiteps = mean(quit_elast), avgreceps = mean(rec_elast), avgelast = mean(wage_elast), medelast = median(wage_elast), stdelast = std(wage_elast), minelast = minimum(wage_elast), maxelast = maximum(wage_elast), migrate = migrat, switchrate = switchrat, avgshr = mean(s[capT,:]), minshr = minimum(s[capT,:]), maxshr = maximum(s[capT,:]), disc_fac = convert(Float64,discfac), comptime = TimerOutputs.tottime(to)/1e9, cor_w_elast = cor(wage,wage_elast), cor_a_elast = cor(amen,wage_elast))
    return result
end

function exporter(numfirm,nummkt,beta=0)
    if beta==0
        fname  = string("results_",nummkt,"_",numfirm,"_static.csv")
        fname2 = string("results_all_",nummkt,"_",numfirm,"_static.csv")
    else
        fname  = string("results_",nummkt,"_",numfirm,".csv")
        fname2 = string("results_all_",nummkt,"_",numfirm,".csv")
    end
    results = DataFrame()
    temp = wrapper(7.96,.1,.14,-.07,numfirm,nummkt,true,1,-50,-50,.0001,beta)
    append!(results,temp)
    CSV.write(fname, temp)
    for thmov in [-6:.5:0;]
        for thsw in [-6:.5:0;]
            for gam = 1:5
                temp = wrapper(7.96,.1,.14,-.07,numfirm,nummkt,true,gam,thmov,thsw,.0001,beta)
                append!(results,temp)
                CSV.write(fname, temp; append=true)
            end
        end
    end
    CSV.write(fname2, results)
    return nothing
end

ext: jl

max_stars_repo_name: tyleransom/frictions-JHR

max_issues_repo_name: tyleransom/frictions-JHR

max_forks_repo_name: tyleransom/frictions-JHR


text:
#######################################
#
# Create sweeps
#
#######################################

function create_sweeps(epochs::Array, epochsPerSweep::Int)

    epochsLen = size(epochs)[1]
    epochsNum = size(epochs)[2]
    chansNum  = size(epochs)[3]

    sweepLen = epochsLen * epochsPerSweep
    sweepNum = round(Int, floor(epochsNum / epochsPerSweep))
    sweeps = zeros(Float64, (sweepLen, sweepNum, chansNum))

    sweep = 1
    while sweep <= sweepNum

        sweepStart = (sweep-1)*(epochsPerSweep)+1
        sweepStop  = sweepStart + epochsPerSweep-1

        sweeps[:,sweep,:] = reshape(epochs[:,sweepStart:sweepStop,:], (sweepLen, 1, chansNum))

        sweep += 1
    end

    Logging.info("Generated $sweepNum sweeps of length $sweepLen for $chansNum channels")

    return sweeps
end


text:
"""
    EmptySet() :: AbstractSet{Union{}}
    EmptySet{T}() :: AbstractSet{T}
"""
EmptySet

"""
    SingletonSet(itr) :: AbstractSet
    SingletonSet{T}(itr) :: AbstractSet{T}
"""
SingletonSet

abstract type AbstractMicroSet{T,L} <: AbstractSet{T} end
abstract type AbstractEmptySet{T,L} <: AbstractMicroSet{T,L} end
abstract type AbstractSingletonSet{T,L} <: AbstractMicroSet{T,L} end

struct EmptySet{T,L} <: AbstractEmptySet{T,L} end
EmptySet() = EmptySet{Union{},Set}()
EmptySet{T}() where {T} = EmptySet{T,Set}()

struct SingletonSet{T,L} <: AbstractSingletonSet{T,L}
    value::T
end
SingletonSet{T}(x::T) where {T} = SingletonSet{T,Set}(x)

@inline getvalue(A::SingletonSet) = A.value
@inline upcast(A::AbstractMicroSet{T,L}) where {T,L} = L(A)::AbstractSet{T}

Base.length(::AbstractEmptySet) = 0
@inline Base.iterate(::AbstractEmptySet) = nothing

Base.length(::AbstractSingletonSet) = 1
@inline Base.iterate(A::AbstractSingletonSet) = (getvalue(A), nothing)
@inline Base.iterate(A::AbstractSingletonSet, ::Nothing) = nothing

emptyshim(::Type{L}, ::Type{T}) where {T,L<:AbstractSet} = EmptySet{T,L}()
singletonshim(::Type{L}, x::T) where {T,L<:AbstractSet} = SingletonSet{T,L}(x)

function EmptySet(itr)
    validate_empty_iterator(itr, :EmptySet)
    if Base.IteratorEltype(itr) isa Base.HasEltype
        return EmptySet{eltype(itr)}()
    else
        return EmptySet()
    end
end

function EmptySet{T}(itr) where {T}
    validate_empty_iterator(itr, :EmptySet)
    return EmptySet{T}()
end

function SingletonSet(itr)
    x = validate_singleton_iterator(itr, :SingletonSet)
    if Base.IteratorEltype(itr) isa Base.HasEltype
        return SingletonSet{eltype(itr)}((x,))
    else
        return SingletonSet{typeof(x)}(x)
    end
end

function SingletonSet{T}(itr) where {T}
    x = validate_singleton_iterator(itr, :SingletonSet)
    return SingletonSet{T}((x,))
end

# fast-path:
SingletonSet((x,)::Tuple{T}) where {T} = SingletonSet{T}(x)
SingletonSet{T}((x,)::Tuple{Any}) where {T} = SingletonSet{T}(convert(T, x)::T)

function Base.showarg(io::IO, x::EmptySet{<:Any,Set}, toplevel::Bool)
    @nospecialize
    print(io, "EmptySet")
    toplevel && print(io, '{', eltype(x), '}')
end

function Base.showarg(io::IO, x::SingletonSet{<:Any,Set}, toplevel::Bool)
    @nospecialize
    print(io, "SingletonSet")
    toplevel && print(io, '{', eltype(x), '}')
end

function Base.show(io::IO, ::MIME"text/plain", x::EmptySet{<:Any,Set})
    @nospecialize
    Base.showarg(io, x, true)
    print(io, "()")
end


text:
# TODO: CB this line supresses the plotlyjs outputs, but can't be included since we don't have a PlotlyJS dependency
#Base.show(io::IO, mm::MIME"text/plain", p::Plots.PlotlyJS.Plot) = show(io, mm, "PlotlyJS Plot with $(length(p.data)) traces")

function _empty_plot(backend::Plots.PlotlyJSBackend)
    return Plots.PlotlyJS.Plot()
end

function _dataframe_plots_internal(
    plot,
    variable::DataFrames.DataFrame,
    time_range::Array,
    backend::Plots.PlotlyJSBackend;
    kwargs...,
)
    names = DataFrames.names(no_datetime(variable))
    traces = plot.data
    plot_length = length(traces)
    seriescolor = permutedims(
        set_seriescolor(
            get(kwargs, :seriescolor, PLOTLY_DEFAULT),
            vcat(ones(plot_length), names),
        )[(plot_length + 1):end],
    )

    save_fig = get(kwargs, :save, nothing)
    y_label = get(kwargs, :y_label, "")
    title = get(kwargs, :title, " ")
    stack = get(kwargs, :stack, false)
    bar = get(kwargs, :bar, false)
    nofill = get(kwargs, :nofill, !bar && !stack)

    time_interval =
        IS.convert_compound_period(length(time_range) * (time_range[2] - time_range[1]))
    interval =
        Dates.Millisecond(Dates.Hour(1)) / Dates.Millisecond(time_range[2] - time_range[1])

    isnothing(plot) && _empty_plot()

    if isempty(variable)
        @warn "Plot dataframe empty: skipping plot creation"
        plot_data = Array{Float64}(undef, 0, 0)
    else
        plot_data = Matrix(no_datetime(variable))
    end

    plot_kwargs = Dict()
    if bar
        plot_data = sum(plot_data, dims = 1) ./ interval
        showtxicklabels = false
        if nofill
            plot_kwargs[:type] = "scatter"
            plot_data = [plot_data; plot_data]
            plot_kwargs[:x] = [-0.5, 0.5]
        else
            plot_kwargs[:type] = "bar"
            plot_kwargs[:fill] = "tonexty"
        end
    else
        if !nofill && stack
            plot_kwargs[:fill] = "tonexty"
        end
        plot_kwargs[:plot_type] = "scatter"
        plot_kwargs[:x] = time_range
    end

    plot_kwargs[:line_shape] = get(kwargs, :stair, false) ? "hv" : "linear"
    plot_kwargs[:mode] = "lines"
    plot_kwargs[:line_dash] = get(kwargs, :line_dash, "solid")
    plot_kwargs[:showlegend] = true

    for ix in 1:length(names)
        if bar
            plot_kwargs[:marker_color] = seriescolor[ix]
        end
        if stack
            plot_kwargs[:stackgroup] = string(plot_length + 1)
            if nofill
                plot_kwargs[:fillcolor] = "transparent"
            else
                plot_kwargs[:fillcolor] = seriescolor[ix]
            end
        elseif !nofill
            plot_kwargs[:stackgroup] = string(ix + plot_length)
        end
        plot_kwargs[:line_color] = seriescolor[ix]
        plot_kwargs[:name] = names[ix]

        trace = Plots.PlotlyJS.scatter(; y = plot_data[:, ix], plot_kwargs...)
        push!(traces, trace)
    end
    layout_kwargs = Dict{Symbol, Any}()
    layout_kwargs[:yaxis] =
        Plots.PlotlyJS.attr(; showticklabels = true, rangemode = "tozero", title = y_label)
    layout_kwargs[:xaxis] =
        Plots.PlotlyJS.attr(; showticklabels = !bar, title = "$time_interval")
    layout_kwargs[:title] = "$title"
    layout_kwargs[:barmode] = stack ? "stack" : "group"
    merge!(layout_kwargs, kwargs)
    Plots.PlotlyJS.relayout!(plot, Plots.PlotlyJS.Layout(; layout_kwargs...))

    get(kwargs, :set_display, true) && display(Plots.PlotlyJS.plot(plot))
    if !isnothing(save_fig)
        title = title == " " ? "dataframe" : title
        format = get(kwargs, :format, "png")
        save_plot(plot, joinpath(save_fig, "$title.$format"), backend; kwargs...)
    end
    return plot
end

function save_plot(plot, filename::String, backend::Plots.PlotlyJSBackend; kwargs...) # this needs to be typed but Plots.PlotlyJS.Plot doesn't exist until PlotlyJS is loaded
    save_kwargs = Dict{Symbol, Any}((
        (k, v) for (k, v) in kwargs if k in SUPPORTED_PLOTLY_SAVE_KWARGS
    ))
    @info "saving plot" filename
    if last(splitext(filename)) == ".html"
        open(filename, "w") do io
            show(io, MIME("text/html"), plot; save_kwargs...)
        end
    else
        Plots.PlotlyJS.savefig(plot, filename; save_kwargs...)
    end
    return filename
end

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
module tensioncompression1
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots
using Test
function test()
	mr = DeforModelRed3D
	E, nu = 7.0*phun("MPa"), 0.3
	m = MatDeforNeohookean(mr, E, nu)
	L= 6/2*phun("mm");
	H = 2/2*phun("mm");
	W = 2/2*phun("mm");
	umag = 1.5*phun("mm");# Magnitude of the displacement
	nincr = 48
	utol = 10e-7;
	graphics = ~true;
	maxdu_tol = W/1e7;
	maxiter = 5
	tolerance = W / 1000

	fens, fes = H8block(L, W, H, 2, 1, 1)

	l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	e1 = FDataDict("node_list"=>l1, "component"=>1)
	l2  = selectnode(fens; box = [-Inf,Inf,0,0,-Inf,Inf], inflate  =  tolerance)
	e2 = FDataDict("node_list"=>l2, "component"=>2)
	l3  = selectnode(fens; box = [-Inf,Inf,-Inf,Inf,0,0], inflate  =  tolerance)
	e3 = FDataDict("node_list"=>l3, "component"=>3)

	movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
	move(x, lambda) = table(lambda);
	e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

	# Rmout = fill(0.0, 3, 3)
	# rv = vec([-0.56 -0.1361 0.35])
	# rotmat3!(Rmout, rv)
	# mcsys = CSys(Rmout)
	# femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
	femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

	region1 = FDataDict("femm"=>femm)
	modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

	modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
	modeldata["maxdu_tol"] = maxdu_tol;
	modeldata["maxiter"] = maxiter;
	modeldata["line_search"]  = true;
	modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
		# @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
	end
	Ux = FFlt[]; Rx = FFlt[]
	modeldata["increment_observer"] = (lambda, modeldata) -> begin
		push!(Ux, mean(modeldata["un1"].values[movel1,1]));
		push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
	end

	modeldata = nonlinearstatics(modeldata);
	# @show minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)

	@test [minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)] ≈
	 [-0.0015, -6.547464398217384, 0.0014999999999999994, 2.6479612347602295]
	# pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
	# display(pl)
	return true
end
end
using .tensioncompression1
tensioncompression1.test()

module tensioncompression2
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanADModule: MatDeforNeohookeanAD
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots
using Test
function test()
	mr = DeforModelRed3D
	E, nu = 7.0*phun("MPa"), 0.3
	m = MatDeforNeohookeanAD(mr, E, nu)
	L= 6/2*phun("mm");
	H = 2/2*phun("mm");
	W = 2/2*phun("mm");
	umag = 1.5*phun("mm");# Magnitude of the displacement
	nincr = 48
	utol = 10e-7;
	graphics = ~true;
	maxdu_tol = W/1e7;
	maxiter = 5
	tolerance = W / 1000

	fens, fes = H8block(L, W, H, 2, 1, 1)

	l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	e1 = FDataDict("node_list"=>l1, "component"=>1)
	l2  = selectnode(fens; box = [-Inf,Inf,0,0,-Inf,Inf], inflate  =  tolerance)
	e2 = FDataDict("node_list"=>l2, "component"=>2)
	l3  = selectnode(fens; box = [-Inf,Inf,-Inf,Inf,0,0], inflate  =  tolerance)
	e3 = FDataDict("node_list"=>l3, "component"=>3)

	movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
	move(x, lambda) = table(lambda);
	e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

	# Rmout = fill(0.0, 3, 3)
	# rv = vec([-0.56 -0.1361 0.35])
	# rotmat3!(Rmout, rv)
	# mcsys = CSys(Rmout)
	# femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
	femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

	region1 = FDataDict("femm"=>femm)
	modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

	modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
	modeldata["maxdu_tol"] = maxdu_tol;
	modeldata["maxiter"] = maxiter;
	modeldata["line_search"]  = true;
	modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
	    # @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
	end
	Ux = FFlt[]; Rx = FFlt[]
	modeldata["increment_observer"] = (lambda, modeldata) -> begin
	    push!(Ux, mean(modeldata["un1"].values[movel1,1]));
	    push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
	end

	modeldata = nonlinearstatics(modeldata);
	# @show minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)

	@test [minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)] ≈
	 [-0.0015, -6.547464398217384, 0.0014999999999999994, 2.6479612347602295]
	# pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
	# display(pl)
	return true
end
end
using .tensioncompression2
tensioncompression2.test()

module tensioncompression3
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforStVKModule: MatDeforStVK
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots
using Test
function test()
	mr = DeforModelRed3D
	E, nu = 7.0*phun("MPa"), 0.3
	m = MatDeforStVK(mr, E, nu)
	L= 6/2*phun("mm");
	H = 2/2*phun("mm");
	W = 2/2*phun("mm");
	umag = 1.0*phun("mm");# Magnitude of the displacement
	nincr = 48
	utol = 10e-7;
	graphics = ~true;
	maxdu_tol = W/1e7;
	maxiter = 5
	tolerance = W / 1000

	fens, fes = H8block(L, W, H, 2, 1, 1)

	l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	e1 = FDataDict("node_list"=>l1, "component"=>1)
	l2  = selectnode(fens; box = [-Inf,Inf,0,0,-Inf,Inf], inflate  =  tolerance)
	e2 = FDataDict("node_list"=>l2, "component"=>2)
	l3  = selectnode(fens; box = [-Inf,Inf,-Inf,Inf,0,0], inflate  =  tolerance)
	e3 = FDataDict("node_list"=>l3, "component"=>3)

	movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
	move(x, lambda) = table(lambda);
	e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

	# Rmout = fill(0.0, 3, 3)
	# rv = vec([-0.56 -0.1361 0.35])
	# rotmat3!(Rmout, rv)
	# mcsys = CSys(Rmout)
	# femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
	femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

	region1 = FDataDict("femm"=>femm)
	modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

	modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
	modeldata["maxdu_tol"] = maxdu_tol;
	modeldata["maxiter"] = maxiter;
	modeldata["line_search"]  = true;
	modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
	    # @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
	end
	Ux = FFlt[]; Rx = FFlt[]
	modeldata["increment_observer"] = (lambda, modeldata) -> begin
	    push!(Ux, mean(modeldata["un1"].values[movel1,1]));
	    push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
	end

	modeldata = nonlinearstatics(modeldata);
	# @show minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)

	@test [minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)] ≈
	 [-0.001, -1.2962962962962954, 0.0009999999999999996, 3.629629629629627]
	# pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
	# display(pl)
	return true
end
end
using .tensioncompression3
tensioncompression3.test()

module tensioncompression4
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforStVKADModule: MatDeforStVKAD
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots
using Test
function test()
	mr = DeforModelRed3D
	E, nu = 7.0*phun("MPa"), 0.3
	m = MatDeforStVKAD(mr, E, nu)
	L= 6/2*phun("mm");
	H = 2/2*phun("mm");
	W = 2/2*phun("mm");
	umag = 1.0*phun("mm");# Magnitude of the displacement
	nincr = 48
	utol = 10e-7;
	graphics = ~true;
	maxdu_tol = W/1e7;
	maxiter = 5
	tolerance = W / 1000

	fens, fes = H8block(L, W, H, 2, 1, 1)

	l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	e1 = FDataDict("node_list"=>l1, "component"=>1)
	l2  = selectnode(fens; box = [-Inf,Inf,0,0,-Inf,Inf], inflate  =  tolerance)
	e2 = FDataDict("node_list"=>l2, "component"=>2)
	l3  = selectnode(fens; box = [-Inf,Inf,-Inf,Inf,0,0], inflate  =  tolerance)
	e3 = FDataDict("node_list"=>l3, "component"=>3)

	movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
	table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
	move(x, lambda) = table(lambda);
	e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

	# Rmout = fill(0.0, 3, 3)
	# rv = vec([-0.56 -0.1361 0.35])
	# rotmat3!(Rmout, rv)
	# mcsys = CSys(Rmout)
	# femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
	femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

	region1 = FDataDict("femm"=>femm)
	modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

	modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
	modeldata["maxdu_tol"] = maxdu_tol;
	modeldata["maxiter"] = maxiter;
	modeldata["line_search"]  = true;
	modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
	    # @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
	end
	Ux = FFlt[]; Rx = FFlt[]
	modeldata["increment_observer"] = (lambda, modeldata) -> begin
	    push!(Ux, mean(modeldata["un1"].values[movel1,1]));
	    push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
	end

	modeldata = nonlinearstatics(modeldata);
	# @show minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)

	@test [minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)] ≈
	 [-0.001, -1.2962962962962954, 0.0009999999999999996, 3.629629629629627]
	# pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
	# display(pl)
	return true
end
end
using .tensioncompression4
tensioncompression4.test()

module tensioncompression5
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforI1RivlinADModule: MatDeforI1RivlinAD
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots
using Test

function test()
	mr = DeforModelRed3D
	E, nu = 7.0*phun("MPa"), 0.3
	mu = E / 2 / (1 + nu)
	c1 = mu / 2
	c2 = c1 / 2
	K = E / 3 / (1 - 2 * nu)
	m = MatDeforI1RivlinAD(mr, c1, c2, K)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    umag = 1.5*phun("mm");# Magnitude of the displacement
    nincr = 48
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 5
    tolerance = W / 1000

    fens, fes = H8block(L, W, H, 2, 1, 1)

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    l2  = selectnode(fens; box = [-Inf,Inf,0,0,-Inf,Inf], inflate  =  tolerance)
    e2 = FDataDict("node_list"=>l2, "component"=>2)
    l3  = selectnode(fens; box = [-Inf,Inf,-Inf,Inf,0,0], inflate  =  tolerance)
    e3 = FDataDict("node_list"=>l3, "component"=>3)

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
    move(x, lambda) = table(lambda);
    e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

    # Rmout = fill(0.0, 3, 3)
    # rv = vec([-0.56 -0.1361 0.35])
    # rotmat3!(Rmout, rv)
    # mcsys = CSys(Rmout)
    # femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e3, e4])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        # @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
    end
    Ux = FFlt[]; Rx = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
        push!(Ux, mean(modeldata["un1"].values[movel1,1]));
        push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
    end

    modeldata = nonlinearstatics(modeldata);
    @test [minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx)] ≈ [-0.0015, -6.43875417356319, 0.0014999999999999994, 4.473481547305163]

    # pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
    # display(pl)
end # function stvkad_h8

end
using .tensioncompression5
tensioncompression5.test()

module m7test13bc
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using LinearAlgebra: norm
using BenchmarkTools
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3

    m1 = FinEtoolsDeforNonlinear.MatDeforNeohookeanModule.MatDeforNeohookean(mr, E, nu)
    m2 = FinEtoolsDeforNonlinear.MatDeforNeohookeanADModule.MatDeforNeohookeanAD(mr, E, nu)
    # @show m
    update! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.update!
    tangentmoduli! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.tangentmoduli!

    stress = fill(0.0, 6)
    D = fill(0.0, 6, 6)
    output = FFlt[]
    statev = FFlt[]
    tn = 0.0
    dtn = 0.0
    loc = [0.0 0.0 0.0]
    label = 0
    quantity=:nothing
    Fn1 = [1.1 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]

    update!(m1, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[863901.0097711232, 349914.99578510894, 349914.99578510894, 0.0, 0.0, 0.0]) / E < 1.0e-7
    tangentmoduli!(m1, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @btime $tangentmoduli!($m1, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    D1 = deepcopy(D)

    update!(m2, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[863901.0097711232, 349914.99578510894, 349914.99578510894, 0.0, 0.0, 0.0] ) / E < 1.0e-7
    tangentmoduli!(m2, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @btime $tangentmoduli!($m2, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    D2 = deepcopy(D)
    @test norm(D1-D2) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.01 0; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[40184.6915460777, 95648.94230769236, 40184.69154607771, 0.0, 0.0, 0.0] ) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.0 0.01; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[201.92307692305477, 740.4317307692086, 471.1538461537944, 0.0, 0.0, 26927.78846153846]) / E < 1.0e-7
end
end
using .m7test13bc
m7test13bc.test()


text:
function ueberpruefeSchuld(Schuld,f_nah,OriginalVektor)

if Schuld == f_nah
    println("Eure Antwort ist korrekt!");
elseif Schuld == f_nah-OriginalVektor
    println("Eure Antwort ist korrekt!");
else
    println("Schaut euch noch einmal den Unterschied zwischen den beiden Fehlerberechnungen an. Welcher Vektor hat sich verändert?");
end
end


text: abstract Interpolator{T}


text:
@testset "DBitstream" begin
    @testset "DBit" begin
        @test zero(DBit) == DBit(false)
        @test one(DBit) == DBit(true)
        @test float(DBit(true)) == 1
        @test float(DBit(false)) == -1
    end

    x = DBitstream()
    y = DBitstream()

    push!(x, DBit(true))
    push!(y, DBit(false))
    @test float(pop!(x)) == 1
    @test float(pop!(y)) == -1

    ops = (+, -, *, /)
    xstream = [DBit(i) for i in rand(Bool, length(ops) * 2)]
    ystream = [DBit(i) for i in rand(Bool, length(ops) * 2)]
    push!(x, xstream)
    push!(y, ystream)
    i = 1
    xstream = float.(xstream)
    ystream = float.(ystream)
    @testset for op in ops
        z = rand(Float64)
        @test op(pop!(x), pop!(y)) == op(xstream[i], ystream[i])
        @test op(pop!(x), z) == op(xstream[i + 1], z)
        @test op(z, pop!(y)) == op(z, ystream[i + 1])
        i += 2
    end
end


text:
mutable struct FGMState{VT, IT}
    iter::Integer
    x::VT
    xₙ::VT
    y::VT
    yₙ::VT
    ∇y::VT
    ∇xₙ::VT
    t::VT
    stepstate::IT
end

struct FastGradientMethodIterator{MT, VT, IT, T}
    H::MT
    b::VT
    x₀::VT
    proj::Function
    step::IT
    L::T
end

function iterate( fgm::FastGradientMethodIterator )
    VT = typeof( fgm.x₀ )
    ET = eltype( fgm.x₀ )

    len = length( fgm.x₀ )
    zer = fill( zero(ET), (len, ) )

    # Setup the step length computation iterator
    (β, stepstate) = iterate( fgm.step )

    IT = typeof( stepstate )

    state = FGMState{VT, IT}( 0, fgm.x₀, fgm.x₀, fgm.x₀, fgm.x₀, zer, zer, zer, stepstate )

    # Run the first pass
    return iterate( fgm, state )
end

function iterate( fgm::FastGradientMethodIterator, state::FGMState )
    state.iter = state.iter + 1

    # Copy the previous values into the current values
    state.x = state.xₙ
    state.y = state.yₙ

    # Compute the gradient
    state.∇y = fgm.H*state.y + fgm.b

    # Compute the new step
    state.t = state.y - (1/fgm.L) * state.∇y

    # Project into the feasible space
    state.xₙ = fgm.proj( state.t )

    # Compute the step length and apply it
    (β, state.stepstate) = iterate( fgm.step, state.stepstate )

    state.yₙ = state.xₙ + β*( state.xₙ - state.x )

    # Compute the residual of the next point for the stopping condition computation
    state.∇xₙ = fgm.H*state.xₙ + fgm.b

    return (state, state)
end

function header()
    @printf "Iteration | Residual (Conjugate)\n"
end

function disp( state::FGMState )
    @printf "   %3d    |   %.4e\n" state.iter compute_conjugate( state )
end

"""
```julia
    fastgradientmethod( H::AbstractMatrix{T}, b::AbstractVector{T}, proj::Function; x₀::AbstractVector{T} = zeros( T, size( b ) ),
                                                                                L::Union{Nothing, T} = nothing,
                                                                                μ::Union{Nothing, T} = nothing,
                                                                                step::AbstractStep = ConstantStep(),
                                                                                stopconditions::Vector{SC} = [Best(1e-4)],
                                                                                maxiter::Integer = 100,
                                                                                dispinterval::Integer = 1
                      ) where {T <: Number, SC <: AbstractStopCondition}
```

Compute the solution to the quadratic program ``min xᵀHx + bᵀx s.t. x ∈ χ`` where the set constraint
``χ`` is supplied as a projection operation that maps a value of ``x`` to a value inside the set ``χ``.
The projection operator is specified as the function `proj`, which must take a single argument that is
a vector of type `T` (the current iterate) and return a vector of type `T` (the projection of the iterate
onto the set ``Χ``).

The behavior of the algorithm can be modified using the keyword arguments to set various algorithm parameters:
* x₀ - The initial iterate
* L - User-specified value for the maximum eigenvalue of `H`. If not specified, then the initial phase of the
      algorithm will perform an eigenvalue decomposition on `H` and set `L` to its maximum eigenvalue.
* μ - User-specified value for the minimum eigenvalue of `H`. If not specified, then the initial phase of the
      algorithm will perform an eigenvalue decomposition on `H` and set `L` to its minimum eigenvalue.
* step - How the step length should be computed at each iteration.
* stopconditions - A vector containing the stopping conditions that the algorithm uses to determine convergence.
                   A boolean OR of all the stopping conditions will be used, so when the first one is met, the
                   algorithm will terminate. Set to `nothing` to only use an iteration bound.
* maxiter - The maximum number of iterations to perform
* dispinterval - Print iteration information every `dispinterval` iterations. Set to 0 to disable printing.
"""
function fastgradientmethod( H::AbstractMatrix{T}, b::AbstractVector{T}, proj::Function; x₀::AbstractVector{T} = zeros( T, size( b ) ),
                                                                                L::Union{Nothing, T} = nothing,
                                                                                μ::Union{Nothing, T} = nothing,
                                                                                step::AbstractStep = ConstantStep(),
                                                                                stopconditions::Union{Vector{SC}, Nothing} = [Best(1e-4)],
                                                                                maxiter::Integer = 100,
                                                                                dispinterval::Integer = 1
                           ) where {T <: Number, SC <: AbstractStopCondition}

    # Compute the lipschitz and convexity parameter if they aren't provided
    if isnothing( L ) || isnothing( μ )
        eig = eigvals( H )

        isnothing( L ) && ( L = maximum( eig ) )
        isnothing( μ ) && ( μ = minimum( eig ) )
    end

    n = size( b, 1 )

    init!( step, L, μ )

    iter = 0

    # Create the actual iterator that we use
    MT = typeof( H )
    VT = typeof( b )
    IT = typeof( step )
    LT = typeof( L )

    iter = FastGradientMethodIterator{MT, VT, IT, LT}( H, b, x₀, proj, step, L )

    # Setup the stopping conditions
    if stopconditions != nothing
        for cond in stopconditions
            init!( cond, n, L, μ )
        end

        iter = halt( iter, stopconditions )
    end

    iter = take( iter, maxiter )
    iter = enumerate( iter )

    # Configure the printing of iteration information
    if dispinterval > 0
        iter = apply( iter, disp, dispinterval )
        header()
    end

    # Actually loop the algorithm
    (_, finalstate) = loop( iter )

    return (finalstate.xₙ, finalstate.iter)
end


text:
using Gadfly, Distributions

function Q(x) 0.5 * erfc(x / sqrt(2)) end

sdb = [0:0.1:10] # SNR, dB per bit
sb = 10 .^ (sdb / 10) # SNR per bit
sig = sqrt(0.5 ./ sb)
sig1 = 1 ./ sig
PE = Q(sig1)

peplot = plot(x=sdb, y=PE, Geom.line, Scale.y_log10,
              Theme(default_color=color("black")),
              Guide.xticks(ticks=[0:1:10]),
              Guide.xlabel("E/N0 (dB)"), Guide.ylabel("PE"))

draw(PDF("ctKovanikov6.pdf", 4inch, 3inch), peplot)

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
# The file for defining what a one plus evolution scheme looks like

export oneplus

# Pass it your inputs, outputs, and fitness function
# The default Chromosome type is CGP, the seed is fixed, and there should be no
# expert
# String ID is optional
function oneplus(nin::Int64, nout::Int64, fitness::Function;
                 ctype::DataType=CGPChromo, seed::Int64=0, expert::Any=nothing,
                 id::String="")

    # Set up a population of CGP individuals
    population = Array{ctype}(Config.lambda)

    # For each of them, set them up with a Chromosome
    for i in eachindex(population)
        # This line actually creates random individuals
        population[i] = ctype(nin, nout)
    end

    # If we are using an expert, plug it in here
    if expert != nothing
        population[1] = expert
    end

    # We assume the best one is the first one
    best = population[1]

    # Setting some basic parameters
    max_fit = -Inf
    eval_count = 0
    fits = -Inf*ones(Float64, Config.lambda)

    log_anyways = 200
    buffer = Config.lambda

    # While we have less than the number of total evaluations we are going to
    # do
    while eval_count < Config.total_evals
        # evaluation
        log_gen = false

        # We are going to loop through this using Threads.@threads
        for p in eachindex(population)
            println("Evaluating #", p + eval_count, "(", (100 * (p + eval_count)/Config.total_evals), "%)")

            if p + eval_count > Config.total_evals
                break
            end

            if fits[p] == -Inf
                # Get fitness score for each and update the array
                fit = fitness(population[p])
                # Update the fitness
                fits[p] = fit
            end
        end
        # This is the only thing we should do in this loop

        # Afterwards, get max and argmax of population
        cur_max = maximum(fits)
        cur_max_id = indmax(fits)

        println("Finished generation.")

        # If max > max_fit, then:
        if cur_max >= max_fit
            # # Update max_fit
            max_fit = cur_max
            # # Update best
            best = clone(population[cur_max_id])
            # # Flip log_gen
            log_gen = true
            println("New Max Fit: ", max_fit)
        end

        # eval_count += len(population)
        eval_count += length(population)

        # if eval_count > Config.total_evals, then
        if eval_count > Config.total_evals
            # # Flip log_gen
            log_gen = true
        end

        if eval_count % log_anyways < buffer
            log_gen = true
        end

        # Call eval
        # Evaluate that generation
        eval(Config.log_function)(id, seed, eval_count, max_fit, best, fitness, ctype, true)

        # if eval_count > Config.total_evals, then
        if eval_count >= Config.total_evals
            # # break
            break
        end

        # selection
        fits .= -Inf
        for p in eachindex(population)
            population[p] = mutate(best)
        end

        # size limit
        for i in eachindex(population)
            # If a mutated individual exceed what we expect,
            # Generate a random one (?)
            # Side Note: Shouldn't we just generate one off the best?
            if length(population[i].nodes) > Config.node_size_cap
                population[i] = ctype(nin, nout)
                fits[i] = -Inf
            end
        end
    end

    # When done, return out the max fit candidate and its Chromosome
    max_fit, best.genes
end


text:
module LibRcl

# import Libdl

include("../src/find_ros2.jl")

# # Load in `deps.jl`, complaining if it does not exist
# const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
# if !isfile(depsjl_path)
#     error("LibXXX was not build properly. Please run Pkg.build(\"LibXXX\").")
# end
# include(depsjl_path)
# # Module initialization function
# function __init__()
#     check_deps()
# end

using CEnum

include("ctypes.jl")
# export Ctm, Ctime_t, Cclock_t
const SIZE_MAX = ~Csize_t(0)

for lib in WRAPPED_ROS2_LIBS
    @info "Loading types and API for library $lib"
    include(joinpath(@__DIR__, "lib$(lib)_common.jl"))
    include(joinpath(@__DIR__, "lib$(lib)_api.jl"))
end

# export everything
foreach(names(@__MODULE__, all=true)) do s
   if occursin(r"^(rcl|rcutils|rosidl)_", string(s))
       @eval export $s
   end
end

struct RclError <: Exception
    code :: Cint
    msg :: AbstractString
end

struct RcutilsErrorState <: Exception
    message :: String
    file :: String
    line_number :: Int
end

# Base.show(io::IO, err::RclError) = print(io, "RclError($(err.code), \"$(err.msg)\")")

@inline function fixed_size_string(s::T, field_idx, expected_field_name) where T
    U = fieldtype(T, field_idx)
    @assert fieldname(T, field_idx) == expected_field_name fieldname(T, field_idx)
    @assert U <: NTuple U
    @assert eltype(U) == Cchar eltype(U)
    ptr = Ptr{Cchar}(pointer_from_objref(s) + fieldoffset(T, field_idx))
    len = ccall(:strnlen, Csize_t, (Cstring, Csize_t), ptr, U.size)
    unsafe_string(ptr, len)
end

Base.string(err_str::rcutils_error_string_t) = fixed_size_string(err_str, 1, :str)
Base.show(io::IO, err_str::rcutils_error_string_t) = print(io, "$(typeof(err_str))(\"$(string(err_str))\")")

message(err_state::rcutils_error_state_t) = fixed_size_string(err_str, 1, :message)
file(err_state::rcutils_error_state_t) = fixed_size_string(err_str, 2, :file)
line_number(err_state::rcutils_error_state_t) = Int(err_state.line_number)

Base.convert(::Type{RcutilsErrorState}, err_state::rcutils_error_state_t) = RcutilsErrorState(
    message(err_state), file(err_state), line_number(err_state)
)

Base.show(io::IO, err_state::rcutils_error_state_t) = show(io, convert(RcutilsErrorState, err_state))

function checkcall(code :: rcl_ret_t; suppress_throw=false)
    code == RCL_RET_OK && return nothing

    msg = rcutils_get_error_string()
    str = string(msg)
    err = RclError(code, str)

    rcutils_reset_error()

    suppress_throw || throw(err)
    err
end

export RclError, checkcall

end # module


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for LabelModule

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

function label_map_sanity_test(Nv::Int, x::LabelMap, ls::Vector)
   @test nv(x) == Nv
   @test eltype(x) == eltype(ls)
   @test x == x
   @test copy(x) == x
   @test deepcopy(x) == x

   @test haslabel(x, rand(ls))

   u = rand(1 : Nv)
   v = rand(1 : Nv)

   @test decode(x, ls[v]) == v
   @test decode(x, ls) == collect(1:Nv)
   @test decode(x, Pair(ls[u], ls[v])) == EdgeID(u, v)

   @test encode(x, u) == ls[u]
   @test encode(x, 1 : 5) == ls[1 : 5]
   @test encode(x) == ls
   @test encode(x, EdgeID(u, v)) == Pair(ls[u], ls[v])
end

@testset "Label Module" begin

   ls = map(string, 1:10)

   ###
   # IDENTITYLM
   ###
   label_map_sanity_test(10, LabelMap(10), collect(1:10))

   # Add vertex without label
   @test addvertex!(LabelMap(10)) == LabelMap(11)

   # Add vertex with next label
   @test addvertex!(LabelMap(10), 11) == LabelMap(11)

   # Add vertex with discontinuous label
   label_map_sanity_test(11, addvertex!(LabelMap(10), 12), vcat(1:10, 12))

   # Remove Vertex
   @test rmvertex!(LabelMap(10), 10) == LabelMap(9)
   @test rmvertex!(LabelMap(10), [1,2,3]) == LabelMap(7)

   # Subgraph
   @test subgraph(LabelMap(10), [1,4,7]) == LabelMap(3)

   # Setlabel
   @test setlabel!(LabelMap(10), ls) == DictLM(ls)

   # Relabel
   @test isa(relabel!(LabelMap(10), 10, 15), DictLM)

   ###
   ## DICTLM
   ###

   label_map_sanity_test(10, LabelMap(ls), ls)

   # Add vertex without label
   try addvertex!(LabelMap(ls)); @test false catch @test true end

   # Add vertex with correct type
   label_map_sanity_test(11, addvertex!(LabelMap(ls), "11"), vcat(ls, "11"))

   # Add vertex with wrong type
   try addvertex!(LabelMap(ls), 11); @test false catch @test true end

   # Remove vertex
   label_map_sanity_test(9, rmvertex!(LabelMap(ls), 10), ls[1:9])
   label_map_sanity_test(8, rmvertex!(LabelMap(ls), [1, 4]), setdiff(ls, ["1", "4"]))

   # Subgraph
   label_map_sanity_test(5, subgraph(LabelMap(ls), 1:5), ls[1:5])

   # Setlabel
   label_map_sanity_test(10, setlabel!(LabelMap(ls), 1:10), collect(1:10))

   # Relabel
   lm = LabelMap(ls)
   lm = relabel!(lm, 10, "11")
   @test decode(lm, "11") == 10
   @test encode(lm , 10) == "11"
end


text:
using LogExpFunctions
using ChainRulesTestUtils
using InverseFunctions
using OffsetArrays

using Random
using Test

Random.seed!(1234)

include("basicfuns.jl")
include("chainrules.jl")
include("inverse.jl")


text:
const HDFS_OBJ_FILE = 'F'
const HDFS_OBJ_DIR = 'D'
const HDFS_OBJ_INVALID = ' '

# the following may not be entirely true. some 32 bits systems may have 64 bit time_t
if(64 == Base.WORD_SIZE)
    typealias TimeT Int64
else
    typealias TimeT Int32
end

## used to enforce types ##
type HdfsFS
    host::AbstractString
    port::Int
    user::AbstractString
    ptr::Ptr{Void}
    function HdfsFS(host::AbstractString, port::Int, user::AbstractString, pt::Ptr{Void})
        fs = new(host, port, user, pt)
        (pt != C_NULL) && finalizer(fs, finalize_hdfs_fs)
        fs
    end
end

type HdfsFile <: IO
    fs::HdfsFS
    path::AbstractString
    ptr::Ptr{Void}
    function HdfsFile(fs::HdfsFS, path::AbstractString, pt::Ptr{Void})
        f = new(fs, path, pt)
        (pt != C_NULL) && finalizer(f, hdfs_close)
        f
    end
end

type HdfsURL
    url::AbstractString
    function HdfsURL(url::AbstractString)
        uri = URI(url)
        (uri.schema != "hdfs") && error("not a HDFS URL")
        new(url)
    end
end

immutable c_hdfsfileinfo
    mKind::Cint             # file or directory. hoping the enum is int type (it is actually compiler dependent)
    mName::Ptr{UInt8}       # file name
    mLastMod::TimeT         # the last modification time for the file in seconds
    mSize::Int64            # the size of the file in bytes
    mReplication::Cshort    # the count of replicas
    mBlockSize::Int64       # the block size for the file
    mOwner::Ptr{UInt8}      # the owner of the file
    mGroup::Ptr{UInt8}      # the group associated with the file
    mPermissions::Cshort    # permissions associated with the file
    mLastAccess::TimeT      # the last access time for the file in seconds
end

type HdfsFileInfo
    kind::Int8
    name::AbstractString
    last_mod::Int64
    size::Int64
    replications::Int16
    block_sz::Int64
    owner::AbstractString
    grp::AbstractString
    permissions::Int16
    last_access::Int64

    HdfsFileInfo(pt::Ptr{c_hdfsfileinfo}) = (C_NULL != pt) ? HdfsFileInfo(unsafe_load(pt)) : new(HDFS_OBJ_INVALID, "", 0, 0, 0, 0, "", "", 0, 0)
    HdfsFileInfo(cfi::c_hdfsfileinfo) =
        new(cfi.mKind, bytestring(cfi.mName),
                        convert(Int64,cfi.mLastMod),
                        cfi.mSize, cfi.mReplication, cfi.mBlockSize,
                        bytestring(cfi.mOwner), bytestring(cfi.mGroup), cfi.mPermissions,
                        convert(Int64,cfi.mLastAccess))
end


text:
using Flux: @functor
import Flux
import Base

include("./vocab.jl")
include("./onehot.jl")
include("./etype.jl")
include("./embed.jl")
include("./position_embed.jl")


text:
# These tests checks that varying μ0, μB results in the correct scale
# transformations in various energies / fields. Note:
# does not check that the actual absolute values are
# correct.

seed = 1111
μBs = [1, 2, 1, 4]
μ0s = [1, 1, 2, 3]
crystal = Sunny.diamond_crystal()
latsize = (4, 4, 4)

function collect_energy_and_field(int, μB, μ0)
    sys = SpinSystem(crystal, [int], latsize; μB=μB, μ0=μ0)
    Random.seed!(seed)
    rand!(sys)
    return energy(sys), field(sys)
end

# All exchange interactions should be invariant under μ0, μB changes
function validate_exchanges_scaling(μBs, μ0s)
    exchange_ints = diamond_test_exchanges()

    # Check that for all values of μB, μ0, all energies
    #  and fields are unchanged compared to first set of values.
    ref_μB, ref_μ0 = μBs[1], μ0s[1]
    for int in exchange_ints
        (ref_E, ref_B) = collect_energy_and_field(int, ref_μB, ref_μ0)

        pass = true
        for (μB, μ0) in zip(μBs[2:end], μ0s[2:end])
            (E, B) = collect_energy_and_field(int, μB, μ0)
            pass &= E ≈ ref_E
            pass &= B ≈ ref_B
        end
        @test pass
    end
end

# Zeeman energies/fields should scale linearly with μB, invariant to μ0
function validate_field_scaling(μBs, μ0s)
    ref_μB, ref_μ0 = μBs[1], μ0s[1]

    ext_field = external_field(rand(3))

    (ref_E, ref_B) = collect_energy_and_field(ext_field, ref_μB, ref_μ0)

    pass = true
    for (μB, μ0) in zip(μBs[2:end], μ0s[2:end])
        (E, B) = collect_energy_and_field(ext_field, μB, μ0)
        pass &= E ≈ (μB / ref_μB) * ref_E
        pass &= B ≈ (μB / ref_μB) .* ref_B
    end
    @test pass
end

# Dipole-dipole interactions should scale linearly with μ0, and
#  quadratically with μB
function validate_dipole_scaling(μBs, μ0s)
    ref_μB, ref_μ0 = μBs[1], μ0s[1]

    dip_dip = dipole_dipole()

    (ref_E, ref_B) = collect_energy_and_field(dip_dip, ref_μB, ref_μ0)

    pass = true
    for (μB, μ0) in zip(μBs[2:end], μ0s[2:end])
        (E, B) = collect_energy_and_field(dip_dip, μB, μ0)
        pass &= E ≈ (μ0 / ref_μ0) * (μB / ref_μB)^2 * ref_E
        pass &= B ≈ (μ0 / ref_μ0) * (μB / ref_μB)^2 .* ref_B
    end
    @test pass
end

@testset verbose=true "Unit Scaling" begin
    @testset "Exchange Scaling" begin
        validate_exchanges_scaling(μBs, μ0s)
    end
    @testset "Zeeman Scaling" begin
        validate_field_scaling(μBs, μ0s)
    end
    @testset "Dipole-Dipole Scaling" begin
        validate_dipole_scaling(μBs, μ0s)
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

BASE_DIR = "."
if length(ARGS) == 1
    BASE_DIR = ARGS[1]
elseif length(ARGS) > 1
    throw("Error: create_metadata.jl accepts only 1 optional argument for file path")
end

id_to_label = Dict{AbstractString, Int}()
open("$BASE_DIR/tiny-imagenet-200/wnids.txt", "r") do ids
    for (index, line) in enumerate(eachline(ids))
        id_to_label[chomp(line)] = index - 1
    end
end

val_images = []
open("$BASE_DIR/tiny-imagenet-200/val/val_annotations.txt", "r") do val
    for line in eachline(val)
        _line = split(line)
        push!(val_images, (_line[1], id_to_label[_line[2]]))
    end
end
open("$BASE_DIR/val_metadata.txt", "w") do metadata
    for val in val_images
        write(metadata, "$BASE_DIR/tiny-imagenet-200/val/images/$(val[1]) $(val[2])\n")
    end
end

train_classes = readdir("$BASE_DIR/tiny-imagenet-200/train")
open("$BASE_DIR/train_metadata.txt", "w") do metadata
    for class in train_classes
        image_path = "$BASE_DIR/tiny-imagenet-200/train/$class/images"
        for image in readdir(image_path)
            write(metadata, "$image_path/$image $(id_to_label[class])\n")
        end
    end
end


text:
using Documenter, PiBBP

makedocs(
  # options
  modules = [PiBBP],
  doctest = false
)

deploydocs(
  deps   = Deps.pip("mkdocs", "pygments", "python-markdown-math"),
  repo   = "github.com/kalvotom/PiBBP.jl.git",
  julia  = "0.4",
  osname = "linux"
)


text:
@testset "RCF" begin
  Qx, x = PolynomialRing(FlintQQ)
  k, a = NumberField(x - 1, "a")
  Z = maximal_order(k)

  function doit(u::UnitRange, p::Int = 3)
    cnt = 0
    for i in u
      I = ideal(Z, i)
      r, mr = ray_class_group(I, n_quo=p)
      for s in index_p_subgroups(r, fmpz(p), (A,x) -> quo(A, x)[2])
        a = ray_class_field(mr, s)
        if isconductor(a, I, check=false)
          K = number_field(a)
          cnt += 1
        end
      end
    end
    return cnt
  end

  @test doit(1:100) == 16
  @test doit(10^18:10^18+100) == 18
  @test doit(10^18:10^18+1000, 11) == 2

  K, a = quadratic_field(-5)
  H = hilbert_class_field(K)
  L = number_field(H, over_subfield = true)
  @test absolute_degree(L) == 4

  f = x^3 - 36*x -1
  K, a = number_field(f, cached = false, check = false)
  H = hilbert_class_field(K)
  L1 = number_field(H)
  L2 = number_field(H, using_stark_units = true, redo = true)
  @test isisomorphic(Hecke.simplified_absolute_field(L1)[1], Hecke.simplified_absolute_field(L2)[1])[1]

  f = x^2 - x - 100
  K, a = number_field(f, cached = false, check = false)
  H = hilbert_class_field(K)
  L1 = number_field(H)
  L2 = number_field(H, using_stark_units = true, redo = true)
  @test isisomorphic(Hecke.simplified_absolute_field(L1)[1], Hecke.simplified_absolute_field(L2)[1])[1]
  @test length(closure(Hecke.absolute_automorphism_group(H), *)) == 10

  r, mr = Hecke.ray_class_groupQQ(Z, 32, true, 8);
  q, mq = quo(r, [r[1]])
  C = ray_class_field(mr, mq)
  KC = number_field(C)
  auts = Hecke.rel_auto(C)
  @test length(closure(auts, *)) == 8

  k, a = wildanger_field(3, 13)
  zk = maximal_order(k)
  r0 = hilbert_class_field(k)
  @test degree(r0) == 9
  r1 = ray_class_field(4*zk, n_quo = 2)
  r2 = ray_class_field(5*zk, n_quo = 2)
  @test isone(conductor(intersect(r1, r2))[1])
  @test conductor(r1 * r2)[1] == 20*zk
  @test Hecke.issubfield(r1, r1*r2)
  @test !Hecke.issubfield(r0, r1*r2)

  K = simple_extension(number_field(r1))[1]
  ZK = maximal_order(K)
  lp = factor(2*3*5*maximal_order(k))
  for p = keys(lp)
    t = prime_decomposition_type(r1, p)
    l = prime_decomposition(ZK, p)
    @test t[3] == length(l)
    @test valuation(norm(l[1][1]), p) == t[2]
    @test t[1] * t[2] * t[3] == degree(r1)
    @test all(x->valuation(norm(x[1]), p) == t[2], l)
  end

  ln = [(2, true), (3, false), (5, false), (13, true), (31, false)]
  for (p, b) = ln
    @test Hecke.islocal_norm(r1, zk(p)) == b
  end

  Qx, x = PolynomialRing(FlintQQ, "x");
  k, a = NumberField(x^2 - 10, "a");
  A = ray_class_field(35*maximal_order(k))
  B = Hecke.maximal_abelian_subfield(A, k)
  @test A == B
  @test conductor(A) == conductor(B)

  K, _ = compositum(k, wildanger_field(3, 13)[1])
  A = maximal_abelian_subfield(ClassField, K)
  @test degree(A) == 2
  @test degree(intersect(A, cyclotomic_field(ClassField, 10))) == 1
end

@testset "Some abelian extensions" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x - 1, "a")
  O = maximal_order(K)
  r, mr = Hecke.ray_class_groupQQ(O, 7872, true, 16)
  ls = subgroups(r, quotype = [16], fun = (x, y) -> quo(x, y, false)[2])
  @test Hecke.has_quotient(r, [16])
  class_fields = [];
  for s in ls;
    C = ray_class_field(mr, s)::Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}
    CC = number_field(C)
    if Hecke._is_conductor_minQQ(C, 16)
      push!(class_fields, CC)
    end
  end
  @test length(class_fields) == 14

  K, a = quadratic_field(2, cached = false)
  @test length(abelian_extensions(K, [2], fmpz(10)^4, absolutely_distinct = true)) == 38
end


text:
#=

  https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  Let’s consider an infinite number of dice, which we roll one after the other until we see a
  six for the first time. What is the probability of stopping after n dice? The first die is always rolled,
  those with higher numbers D are only rolled if the previous roll did not stop the process.
  """

  Cf dice_6_throws.jl
     ~/blog/rolling_dice4.blog
     ~/psi/rolling_dice4.psi
     ~/webppl/rolling_dice4.wppl

  Theoretical values (1:50)
1:0.16666666666666666
2:0.1388888888888889
3:0.11574074074074076
4:0.09645061728395063
5:0.08037551440329219
6:0.06697959533607684
7:0.05581632944673069
8:0.04651360787227558
9:0.038761339893562986
10:0.032301116577969156
11:0.02691759714830763
12:0.022431330956923026
13:0.018692775797435855
14:0.015577313164529882
15:0.012981094303774901
16:0.010817578586479085
17:0.009014648822065905
18:0.0075122073517215875
19:0.006260172793101323
20:0.005216810660917769
21:0.0043473422174314744
22:0.003622785181192896
23:0.0030189876509940797
24:0.002515823042495067
25:0.002096519202079222
26:0.0017470993350660188
27:0.0014559161125550157
28:0.0012132634271291797
29:0.0010110528559409832
30:0.0008425440466174861
31:0.0007021200388479051
32:0.0005851000323732542
33:0.0004875833603110452
34:0.00040631946692587105
35:0.00033859955577155923
36:0.0002821662964762994
37:0.00023513858039691615
38:0.00019594881699743015
39:0.0001632906808311918
40:0.0001360755673593265
41:0.00011339630613277208
42:9.449692177731006e-5
43:7.874743481442506e-5
44:6.562286234535421e-5
45:5.468571862112852e-5
46:4.5571432184273776e-5
47:3.797619348689481e-5
48:3.164682790574568e-5
49:2.63723565881214e-5
50:2.19769638234345e-5

  The values found by this model (run_model(0))
  1.00000 =>    6631  (0.165775)
2.00000 =>    5661  (0.141525)
3.00000 =>    4593  (0.114825)
4.00000 =>    3867  (0.096675)
5.00000 =>    3215  (0.080375)
6.00000 =>    2700  (0.067500)
7.00000 =>    2210  (0.055250)
8.00000 =>    1845  (0.046125)
9.00000 =>    1515  (0.037875)
10.00000 =>    1282  (0.032050)
11.00000 =>    1137  (0.028425)
12.00000 =>     881  (0.022025)
13.00000 =>     728  (0.018200)
14.00000 =>     567  (0.014175)
15.00000 =>     505  (0.012625)
16.00000 =>     459  (0.011475)
17.00000 =>     364  (0.009100)
18.00000 =>     301  (0.007525)
19.00000 =>     245  (0.006125)
20.00000 =>     225  (0.005625)
22.00000 =>     163  (0.004075)
21.00000 =>     154  (0.003850)
23.00000 =>     126  (0.003150)
24.00000 =>     117  (0.002925)
25.00000 =>      83  (0.002075)
26.00000 =>      65  (0.001625)
27.00000 =>      54  (0.001350)
30.00000 =>      48  (0.001200)
28.00000 =>      45  (0.001125)
29.00000 =>      36  (0.000900)
31.00000 =>      36  (0.000900)
32.00000 =>      27  (0.000675)
34.00000 =>      22  (0.000550)
33.00000 =>      20  (0.000500)
35.00000 =>      12  (0.000300)
37.00000 =>       9  (0.000225)
38.00000 =>       7  (0.000175)
39.00000 =>       6  (0.000150)
42.00000 =>       5  (0.000125)
40.00000 =>       5  (0.000125)
36.00000 =>       5  (0.000125)
46.00000 =>       4  (0.000100)
41.00000 =>       4  (0.000100)
44.00000 =>       3  (0.000075)
47.00000 =>       3  (0.000075)
48.00000 =>       3  (0.000075)
57.00000 =>       2  (0.000050)
49.00000 =>       1  (0.000025)
45.00000 =>       1  (0.000025)
50.00000 =>       1  (0.000025)
43.00000 =>       1  (0.000025)
72.00000 =>       1  (0.000025)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

# Closed form of the probability
function theoretical_prob(n)
    return (1/6.0) * ((5/6.0)^(n-1))
end

# For val = 0: show the full distribution of the length (len)
# otherwise show the probability of len == val
@model function dice_6_throws(val=0)

    max_len = 100
    function throws(a)
        if length(a) > max_len
           return a
        end

        # For some reason, this don't work: it just generates the same number repeatedly.
        # t ~  DiscreteUniform(1,6)
        t = rand(DiscreteUniform(1,6)) # This works.

        # If we see a 6 then we stop
        if t == 6
            return vcat(a,t)
        else
            return throws(vcat(a,t))
        end
    end

    len ~ DiscreteUniform(0,max_len)

    # a = throws(TArray{Int64}(undef,0))
    a = throws(Int64[])
    len = length(a)

    if val == 0
        return len
    else
        return len == val
    end
end

function run_model(len=0)
    model = dice_6_throws(len)

    num_chains = 4

    # chains = sample(model, Prior(), 10_000)

    # chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
    chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, MH(), 10_000)
    # chains = sample(model, MH(), 1_000)

    # chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

    # chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, SMC(1000), 10_000)
    # chains = sample(model, IS(), 10_000)
    #
    # display(chains)
    # show_var_dist_pct(chains,:len,1000)

    println("prob return value:")
    genq = generated_quantities(model, chains)
    show_var_dist_pct(genq,1000)
    if len > 0
        println("\ntheoretical_prob($len): ", theoretical_prob(len))
    end

end

for val in 0:10
    println("\nval:$val")
    run_model(val)
end

println("\nval:25")
run_model(25)

println("\nTheoretical values for 1:50")
for i in 1:50
    println("$i:", theoretical_prob(i))
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricSurface
using GeometricExamples

const Δt     = 2.
const ntime  = 5000
const nx     = 100
const ny     = 100
const nsave  = 10
const nplot  = 4
const run_id = "poincare_2nd_symmetric_vprk_dt2"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_symmetric_poincare_invariant_2nd.jl")


text:
using ManifoldsBase, Test

using ManifoldsBase: DefaultManifold, ℝ
#
# A first artificial (not real) manifold that is modelled as a submanifold
# half plane with euclidean metric is not a manifold but should test all things correctly here
#
struct HalfPlaneManifold <: AbstractDecoratorManifold{ℝ} end
struct PosQuadrantManifold <: AbstractDecoratorManifold{ℝ} end

ManifoldsBase.get_embedding(::HalfPlaneManifold) = ManifoldsBase.DefaultManifold(1, 3)
ManifoldsBase.decorated_manifold(::HalfPlaneManifold) = ManifoldsBase.DefaultManifold(2)
ManifoldsBase.representation_size(::HalfPlaneManifold) = (3,)

ManifoldsBase.get_embedding(::PosQuadrantManifold) = HalfPlaneManifold()
ManifoldsBase.representation_size(::PosQuadrantManifold) = (3,)

function ManifoldsBase.check_point(::HalfPlaneManifold, p)
    return p[1] > 0 ? nothing : DomainError(p[1], "p[1] ≤ 0")
end
function ManifoldsBase.check_vector(::HalfPlaneManifold, p, X)
    return X[1] > 0 ? nothing : DomainError(X[1], "X[1] ≤ 0")
end
function ManifoldsBase.check_point(::PosQuadrantManifold, p)
    return p[2] > 0 ? nothing : DomainError(p[1], "p[2] ≤ 0")
end
function ManifoldsBase.check_vector(::PosQuadrantManifold, p, X)
    return X[2] > 0 ? nothing : DomainError(X[1], "X[2] ≤ 0")
end

ManifoldsBase.embed(::HalfPlaneManifold, p) = reshape(p, 1, :)
ManifoldsBase.embed(::HalfPlaneManifold, p, X) = reshape(X, 1, :)

ManifoldsBase.project!(::HalfPlaneManifold, q, p) = (q .= [p[1] p[2] 0.0])
ManifoldsBase.project!(::HalfPlaneManifold, Y, p, X) = (Y .= [X[1] X[2] 0.0])

function ManifoldsBase.get_coordinates_orthonormal!(
    ::HalfPlaneManifold,
    Y,
    p,
    X,
    ::ManifoldsBase.RealNumbers,
)
    return (Y .= [X[1], X[2]])
end
function ManifoldsBase.get_vector_orthonormal!(
    ::HalfPlaneManifold,
    Y,
    p,
    c,
    ::ManifoldsBase.RealNumbers,
)
    return (Y .= [c[1] c[2] 0.0])
end

function ManifoldsBase.active_traits(f, ::HalfPlaneManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsEmbeddedSubmanifold())
end
function ManifoldsBase.active_traits(f, ::PosQuadrantManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsEmbeddedSubmanifold())
end

#
# A second manifold that is modelled as just isometrically embedded but not a submanifold
#
struct AnotherHalfPlanemanifold <: AbstractDecoratorManifold{ℝ} end

ManifoldsBase.get_embedding(::AnotherHalfPlanemanifold) = ManifoldsBase.DefaultManifold(3)
function ManifoldsBase.decorated_manifold(::AnotherHalfPlanemanifold)
    return ManifoldsBase.DefaultManifold(2)
end
ManifoldsBase.representation_size(::AnotherHalfPlanemanifold) = (2,)

function ManifoldsBase.active_traits(f, ::AnotherHalfPlanemanifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsIsometricEmbeddedManifold())
end

function ManifoldsBase.embed!(::AnotherHalfPlanemanifold, q, p)
    q[1:2] .= p
    q[3] = 0
    return q
end
function ManifoldsBase.embed!(::AnotherHalfPlanemanifold, Y, p, X)
    Y[1:2] .= X
    Y[3] = 0
    return Y
end
function ManifoldsBase.project!(::AnotherHalfPlanemanifold, q, p)
    return q .= [p[1], p[2]]
end
function ManifoldsBase.project!(::AnotherHalfPlanemanifold, Y, p, X)
    return Y .= [X[1], X[2]]
end

#
# Third example - explicitly mention an embedding.
#
function ManifoldsBase.embed!(
    ::EmbeddedManifold{𝔽,DefaultManifold{nL,𝔽},DefaultManifold{mL,𝔽2}},
    q,
    p,
) where {nL,mL,𝔽,𝔽2}
    n = size(p)
    ln = length(n)
    m = size(q)
    lm = length(m)
    (length(n) > length(m)) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) is longer than embedding dimension $(m).",
        ),
    )
    any(n .> m[1:ln]) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) has entry larger than embedding dimensions ($(m)).",
        ),
    )
    fill!(q, 0)
    q[map(ind_n -> Base.OneTo(ind_n), n)..., ntuple(_ -> 1, lm - ln)...] .= p
    return q
end

function ManifoldsBase.project!(
    ::EmbeddedManifold{𝔽,DefaultManifold{nL,𝔽},DefaultManifold{mL,𝔽2}},
    q,
    p,
) where {nL,mL,𝔽,𝔽2}
    n = size(p)
    ln = length(n)
    m = size(q)
    lm = length(m)
    (length(n) < length(m)) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) is longer than embedding dimension $(m).",
        ),
    )
    any(n .< m[1:ln]) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) has entry larger than embedding dimensions ($(m)).",
        ),
    )
    #  fill q with the „top left edge“ of p.
    q .= p[map(i -> Base.OneTo(i), m)..., ntuple(_ -> 1, lm - ln)...]
    return q
end

#
# A manifold that is a submanifold but otherwise has not implementations
#
struct NotImplementedEmbeddedSubManifold <: AbstractDecoratorManifold{ℝ} end
function ManifoldsBase.get_embedding(::NotImplementedEmbeddedSubManifold)
    return ManifoldsBase.DefaultManifold(3)
end
function ManifoldsBase.decorated_manifold(::NotImplementedEmbeddedSubManifold)
    return ManifoldsBase.DefaultManifold(2)
end
function ManifoldsBase.active_traits(f, ::NotImplementedEmbeddedSubManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsEmbeddedSubmanifold())
end

#
# A manifold that is isometrically embedded but has no implementations
#
struct NotImplementedIsometricEmbeddedManifold <: AbstractDecoratorManifold{ℝ} end
function ManifoldsBase.active_traits(f, ::NotImplementedIsometricEmbeddedManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsIsometricEmbeddedManifold())
end

#
# A manifold that is an embedded manifold but not isometric and has no other implementation
#
struct NotImplementedEmbeddedManifold <: AbstractDecoratorManifold{ℝ} end
function ManifoldsBase.active_traits(f, ::NotImplementedEmbeddedManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsEmbeddedManifold())
end

#
# A Manifold with a fallback
#
struct FallbackManifold <: AbstractDecoratorManifold{ℝ} end
function ManifoldsBase.active_traits(f, ::FallbackManifold, args...)
    return ManifoldsBase.merge_traits(ManifoldsBase.IsExplicitDecorator())
end
ManifoldsBase.decorated_manifold(::FallbackManifold) = DefaultManifold(3)

@testset "Embedded Manifolds" begin
    @testset "EmbeddedManifold basic tests" begin
        M = EmbeddedManifold(
            ManifoldsBase.DefaultManifold(2),
            ManifoldsBase.DefaultManifold(3),
        )
        @test repr(M) ==
              "EmbeddedManifold($(sprint(show, M.manifold)), $(sprint(show, M.embedding)))"
        @test base_manifold(M) == ManifoldsBase.DefaultManifold(2)
        @test base_manifold(M, Val(0)) == M
        @test base_manifold(M, Val(1)) == ManifoldsBase.DefaultManifold(2)
        @test base_manifold(M, Val(2)) == ManifoldsBase.DefaultManifold(2)
        @test get_embedding(M) == ManifoldsBase.DefaultManifold(3)
        @test get_embedding(M, [1, 2, 3]) == ManifoldsBase.DefaultManifold(3)
    end

    @testset "HalfPlanemanifold" begin
        M = HalfPlaneManifold()
        N = PosQuadrantManifold()
        @test repr(M) == "HalfPlaneManifold()"
        @test get_embedding(M) == ManifoldsBase.DefaultManifold(1, 3)
        @test representation_size(M) == (3,)
        # Check point checks using embedding
        @test is_point(M, [1 0.1 0.1], true)
        @test !is_point(M, [-1, 0, 0]) #wrong dim (3,1)
        @test_throws DomainError is_point(M, [-1, 0, 0], true)
        @test !is_point(M, [-1, 0, 0])
        @test !is_point(M, [1, 0.1]) # size (from embedding)
        @test_throws ManifoldDomainError is_point(M, [1, 0.1], true)
        @test !is_point(M, [1, 0.1])
        @test is_point(M, [1 0 0], true)
        @test !is_point(M, [-1 0 0]) # right size but <0 1st
        @test_throws DomainError is_point(M, [-1 0 0], true) # right size but <0 1st
        @test !is_vector(M, [1 0 0], [1]) # right point, wrong size vector
        @test_throws ManifoldDomainError is_vector(M, [1 0 0], [1], true)
        @test !is_vector(M, [1 0 0], [1])
        @test_throws DomainError is_vector(M, [1 0 0], [-1 0 0], true) # right point, vec 1st <0
        @test !is_vector(M, [1 0 0], [-1 0 0])
        @test is_vector(M, [1 0 0], [1 0 1], true)
        @test !is_vector(M, [-1, 0, 0], [0, 0, 0])
        @test_throws ManifoldDomainError is_vector(M, [-1, 0, 0], [0, 0, 0], true)
        @test_throws DomainError is_vector(M, [1, 0, 0], [-1, 0, 0], true)
        @test !is_vector(M, [-1, 0, 0], [0, 0, 0])
        @test !is_vector(M, [1, 0, 0], [-1, 0, 0])
        # check manifold domain error from embedding to obtain ManifoldDomainErrors
        @test !is_point(N, [0, 0, 0])
        @test_throws ManifoldDomainError is_point(N, [0, 0, 0], true)
        @test !is_vector(N, [1, 1, 0], [0, 0, 0])
        @test_throws ManifoldDomainError is_vector(N, [1, 1, 0], [0, 0, 0], true)
        p = [1.0 1.0 0.0]
        q = [1.0 0.0 0.0]
        X = q - p
        @test ManifoldsBase.check_size(M, p) === nothing
        @test ManifoldsBase.check_size(M, p, X) === nothing
        @test ManifoldsBase.check_size(M, [1, 2]) isa ManifoldDomainError
        @test ManifoldsBase.check_size(M, [1 2 3 4]) isa ManifoldDomainError
        @test ManifoldsBase.check_size(M, p, [1, 2]) isa ManifoldDomainError
        @test ManifoldsBase.check_size(M, p, [1 2 3 4]) isa ManifoldDomainError
        @test embed(M, p) == p
        pE = similar(p)
        embed!(M, pE, p)
        @test pE == p
        P = [1.0 1.0 2.0]
        Q = similar(P)
        @test project!(M, Q, P) == project!(M, Q, P)
        @test project!(M, Q, P) == [1.0 1.0 0.0]
        @test isapprox(M, p, zero_vector(M, p), [0 0 0])
        XZ = similar(X)
        zero_vector!(M, XZ, p)
        @test isapprox(M, p, XZ, [0 0 0])

        XE = similar(X)
        embed!(M, XE, p, X)
        XE2 = embed(M, p, X)
        @test X == XE
        @test XE == XE2

        @test log(M, p, q) == q - p
        Y = similar(p)
        log!(M, Y, p, q)
        @test Y == q - p
        @test exp(M, p, X) == q
        r = similar(p)
        exp!(M, r, p, X)
        @test r == q
        @test distance(M, p, r) == norm(r - p)

        @test retract(M, p, X) == q
        q2 = similar(q)
        @test retract!(M, q2, p, X) == q
        @test q2 == q
        @test inverse_retract(M, p, q) == X
        Y = similar(X)
        @test inverse_retract!(M, Y, p, q) == X
        @test Y == X

        @test vector_transport_along(M, p, X, []) == X
        @test vector_transport_along!(M, Y, p, X, []) == X
        @test parallel_transport_along(M, p, X, []) == X
        @test parallel_transport_along!(M, Y, p, X, []) == X
        @test parallel_transport_direction(M, p, X, X) == X
        @test parallel_transport_direction!(M, Y, p, X, X) == X
        @test parallel_transport_to(M, p, X, q) == X
        @test parallel_transport_to!(M, Y, p, X, q) == X

        @test get_basis(M, p, DefaultOrthonormalBasis()) isa CachedBasis
        Xc = [X[1], X[2]]
        Yc = similar(Xc)
        @test get_coordinates(M, p, X, DefaultOrthonormalBasis()) == Xc
        @test get_coordinates!(M, Yc, p, X, DefaultOrthonormalBasis()) == Xc
        @test get_vector(M, p, Xc, DefaultOrthonormalBasis()) == X
        @test get_vector!(M, Y, p, Xc, DefaultOrthonormalBasis()) == X
    end

    @testset "AnotherHalfPlanemanifold" begin
        M = AnotherHalfPlanemanifold()
        p = [1.0, 2.0]
        pe = embed(M, p)
        @test pe == [1.0, 2.0, 0.0]
        X = [2.0, 3.0]
        Xe = embed(M, pe, X)
        @test Xe == [2.0, 3.0, 0.0]
        @test project(M, pe) == p
        @test project(M, pe, Xe) == X
        # isometric passtthrough
        @test injectivity_radius(M) == Inf
        @test injectivity_radius(M, p) == Inf
        @test injectivity_radius(M, p, ExponentialRetraction()) == Inf
        @test injectivity_radius(M, ExponentialRetraction()) == Inf
    end

    @testset "Test nonimplemented fallbacks" begin
        @testset "Submanifold Embedding Fallbacks & Error Tests" begin
            M = NotImplementedEmbeddedSubManifold()
            A = zeros(2)
            # for a submanifold quite a lot of functions are passed on
            @test ManifoldsBase.check_point(M, [1, 2]) === nothing
            @test ManifoldsBase.check_vector(M, [1, 2], [3, 4]) === nothing
            @test norm(M, [1, 2], [2, 3]) ≈ sqrt(13)
            @test distance(M, [1, 2], [3, 4]) ≈ sqrt(8)
            @test inner(M, [1, 2], [2, 3], [2, 3]) == 13
            @test manifold_dimension(M) == 2 # since base is defined is defined
            @test_throws MethodError project(M, [1, 2])
            @test_throws MethodError project(M, [1, 2], [2, 3]) == [2, 3]
            @test_throws MethodError project!(M, A, [1, 2], [2, 3])
            @test vector_transport_direction(M, [1, 2], [2, 3], [3, 4]) == [2, 3]
            vector_transport_direction!(M, A, [1, 2], [2, 3], [3, 4])
            @test A == [2, 3]
            @test vector_transport_to(M, [1, 2], [2, 3], [3, 4]) == [2, 3]
            vector_transport_to!(M, A, [1, 2], [2, 3], [3, 4])
            @test A == [2, 3]
        end
        @testset "Isometric Embedding Fallbacks & Error Tests" begin
            M2 = NotImplementedIsometricEmbeddedManifold()
            @test base_manifold(M2) == M2
            A = zeros(2)
            # Check that all of these report not to be implemented, i.e.
            @test_throws MethodError exp(M2, [1, 2], [2, 3])
            @test_throws MethodError exp!(M2, A, [1, 2], [2, 3])
            @test_throws MethodError retract(M2, [1, 2], [2, 3])
            @test_throws MethodError retract!(M2, A, [1, 2], [2, 3])
            @test_throws MethodError log(M2, [1, 2], [2, 3])
            @test_throws MethodError log!(M2, A, [1, 2], [2, 3])
            @test_throws MethodError inverse_retract(M2, [1, 2], [2, 3])
            @test_throws MethodError inverse_retract!(M2, A, [1, 2], [2, 3])
            @test_throws MethodError distance(M2, [1, 2], [2, 3])
            @test_throws StackOverflowError manifold_dimension(M2)
            @test_throws MethodError project(M2, [1, 2])
            @test_throws MethodError project!(M2, A, [1, 2])
            @test_throws MethodError project(M2, [1, 2], [2, 3])
            @test_throws MethodError project!(M2, A, [1, 2], [2, 3])
            @test_throws MethodError vector_transport_along(M2, [1, 2], [2, 3], [[1, 2]])
            @test_throws MethodError vector_transport_along(
                M2,
                [1, 2],
                [2, 3],
                [[1, 2]],
                ParallelTransport(),
            )
            @test vector_transport_along!(M2, A, [1, 2], [2, 3], []) == [2, 3]
            @test A == [2, 3]
            @test_throws MethodError vector_transport_direction(M2, [1, 2], [2, 3], [3, 4])
            @test_throws MethodError vector_transport_direction!(
                M2,
                A,
                [1, 2],
                [2, 3],
                [3, 4],
            )
            @test_throws MethodError vector_transport_to(M2, [1, 2], [2, 3], [3, 4])
            @test_throws MethodError vector_transport_to!(M2, A, [1, 2], [2, 3], [3, 4])
        end
        @testset "Nonisometric Embedding Fallback Error Rests" begin
            M3 = NotImplementedEmbeddedManifold()
            @test_throws MethodError inner(M3, [1, 2], [2, 3], [2, 3])
            @test_throws StackOverflowError manifold_dimension(M3)
            @test_throws MethodError distance(M3, [1, 2], [2, 3])
            @test_throws MethodError norm(M3, [1, 2], [2, 3])
            @test_throws MethodError embed(M3, [1, 2], [2, 3])
            @test_throws MethodError embed(M3, [1, 2])
        end
    end
    @testset "Explicit Embeddings using EmbeddedManifold" begin
        M = DefaultManifold(3, 3)
        N = DefaultManifold(4, 4)
        O = EmbeddedManifold(M, N)
        # first test with same length of sizes
        p = ones(3, 3)
        q = zeros(4, 4)
        qT = zeros(4, 4)
        qT[1:3, 1:3] .= 1.0
        embed!(O, q, p)
        @test norm(qT - q) == 0
        qM = embed(O, p)
        @test norm(project(O, qM) - p) == 0
        @test norm(qT - qM) == 0
        # test with different sizes, check that it only fills first element
        q2 = zeros(4, 4, 3)
        q2T = zeros(4, 4, 3)
        q2T[1:3, 1:3, 1] .= 1.0
        embed!(O, q2, p)
        @test norm(q2T - q2) == 0
        O2 = EmbeddedManifold(M, DefaultManifold(4, 4, 3))
        q2M = embed(O2, p)
        @test norm(q2T - q2M) == 0
        # wrong size error checks
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(3, 3, 5))
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(4, 4))
        @test_throws DomainError project!(O, zeros(3, 3, 5), zeros(3, 3))
        @test_throws DomainError project!(O, zeros(4, 4), zeros(3, 3))
    end
    @testset "Explicit Fallback" begin
        M = FallbackManifold()
        # test the explicit fallback to DefaultManifold(3)
        @test inner(M, [1, 0, 0], [1, 2, 3], [0, 1, 0]) == 2
        @test is_point(M, [1, 0, 0])
    end
end


text:
# Basic algebra
#---------------
Base.copy(t::AbstractTensorMap) = Base.copyto!(similar(t), t)

Base.:-(t::AbstractTensorMap) = mul!(similar(t), t, -one(eltype(t)))
function Base.:+(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(one(T), t2, copyto!(similar(t1, T), t1))
end
function Base.:-(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(-one(T), t2, copyto!(similar(t1, T), t1))
end

Base.:*(t::AbstractTensorMap, α::Number) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), t, α)
Base.:*(α::Number, t::AbstractTensorMap) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), α, t)
Base.:/(t::AbstractTensorMap, α::Number) = *(t, one(α)/α)
Base.:\(α::Number, t::AbstractTensorMap) = *(t, one(α)/α)

LinearAlgebra.normalize!(t::AbstractTensorMap, p::Real = 2) = rmul!(t, inv(norm(t, p)))
LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) =
    mul!(similar(t), t, inv(norm(t, p)))

Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) =
    mul!(similar(t1, promote_type(eltype(t1),eltype(t2)), codomain(t1)←domain(t2)), t1, t2)
Base.exp(t::AbstractTensorMap) = exp!(copy(t))

# Special purpose constructors
#------------------------------
Base.zero(t::AbstractTensorMap) = fill!(similar(t), 0)
function Base.one(t::AbstractTensorMap)
    domain(t) == codomain(t) || throw(SectorMismatch("no identity if domain and codomain are different"))
    TensorMap(I, eltype(t), domain(t), domain(t))
end
function one!(t::AbstractTensorMap)
    domain(t) == codomain(t) || throw(SectorMismatch("no identity if domain and codomain are different"))
    for (c,b) in blocks(t)
        copyto!(b, I)
    end
    return t
end

# Equality and approximality
#----------------------------
function Base.:(==)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    (codomain(t1) == codomain(t2) && domain(t1) == domain(t2)) || return false
    for c in blocksectors(t1)
        block(t1, c) == block(t2, c) || return false
    end
    return true
end

function Base.isapprox(t1::AbstractTensorMap, t2::AbstractTensorMap; atol::Real=0, rtol::Real=Base.rtoldefault(eltype(t1), eltype(t2), atol))
    d = norm(t1 - t2)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(t1), norm(t2)))
    else
        return false
    end
end

# In-place methods
#------------------
# Wrapping the blocks in a StridedView enables multithreading if JULIA_NUM_THREADS > 1
# Copy, adjoint! and fill:
function Base.copyto!(tdest::AbstractTensorMap, tsource::AbstractTensorMap)
    codomain(tdest) == codomain(tsource) && domain(tdest) == domain(tsource) || throw(SpaceMismatch())
    for c in blocksectors(tdest)
        copyto!(StridedView(block(tdest, c)), StridedView(block(tsource, c)))
    end
    return tdest
end
function Base.fill!(t::AbstractTensorMap, value::Number)
    for (c,b) in blocks(t)
        fill!(b, value)
    end
    return t
end
function LinearAlgebra.adjoint!(tdest::AbstractTensorMap, tsource::AbstractTensorMap)
    codomain(tdest) == domain(tsource) && domain(tdest) == codomain(tsource) || throw(SpaceMismatch())
    for c in blocksectors(tdest)
        adjoint!(StridedView(block(tdest, c)), StridedView(block(tsource, c)))
    end
    return tdest
end

# Basic vector space methods: addition and scalar multiplication
LinearAlgebra.rmul!(t::AbstractTensorMap, α::Number) = mul!(t, t, α)
LinearAlgebra.lmul!(α::Number, t::AbstractTensorMap) = mul!(t, α, t)

function LinearAlgebra.mul!(t1::AbstractTensorMap, t2::AbstractTensorMap, α::Number)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), StridedView(block(t2, c)), α)
    end
    return t1
end
function LinearAlgebra.mul!(t1::AbstractTensorMap, α::Number, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), α, StridedView(block(t2, c)))
    end
    return t1
end
function LinearAlgebra.axpy!(α::Number, t1::AbstractTensorMap, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpy!(α, StridedView(block(t1, c)), StridedView(block(t2, c)))
    end
    return t2
end
function LinearAlgebra.axpby!(α::Number, t1::AbstractTensorMap, β::Number, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpby!(α, StridedView(block(t1, c)), β, StridedView(block(t2, c)))
    end
    return t2
end

# inner product and norm only valid for spaces with Euclidean inner product
function LinearAlgebra.dot(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where {S<:EuclideanSpace}
    (codomain(t1) == codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    iter = blocksectors(t1)
    if isempty(iter)
        return zero(eltype(t1))*zero(eltype(t2))
    else
        return sum(dim(c)*dot(block(t1,c), block(t2,c)) for c in blocksectors(t1))
    end
end

LinearAlgebra.norm(t::AbstractTensorMap{<:EuclideanSpace}, p::Real) =
    isempty(blocks(t)) ? zero(real(eltype(t))) : _norm(blocks(t), p)
function _norm(blockiterator, p::Real)
    if p == Inf
        maximum(norm(b, p) for (c,b) in blockiterator)
    elseif p == -Inf
        minimum(norm(b, p) for (c,b) in blockiterator)
    elseif p == 1
        sum(dim(c)*norm(b, p) for (c,b) in blockiterator)
    else
        s = sum(dim(c)*norm(b, p)^p for (c,b) in blockiterator)
        return s^inv(oftype(s, p))
    end
end

# TensorMap trace
function LinearAlgebra.tr(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("Trace of a tensor only exist when domain == codomain"))
    return sum(dim(c)*tr(b) for (c,b) in blocks(t))
end

# TensorMap multiplication:
function LinearAlgebra.mul!(tC::AbstractTensorMap, tA::AbstractTensorMap,  tB::AbstractTensorMap, α = true, β = false)
    if !(codomain(tC) == codomain(tA) && domain(tC) == domain(tB) && domain(tA) == codomain(tB))
        throw(SpaceMismatch())
    end
    for c in blocksectors(tC)
        if hasblock(tA, c) # then also tB should have such a block
            A = block(tA, c)
            B = block(tB, c)
            C = block(tC, c)
            if isbitstype(eltype(A)) && isbitstype(eltype(B)) && isbitstype(eltype(C))
                @unsafe_strided A B C mul!(C, A, B, α, β)
            else
                mul!(StridedView(C), StridedView(A), StridedView(B), α, β)
            end
        elseif β != one(β)
            rmul!(block(tC, c), β)
        end
    end
    return tC
end

# TensorMap inverse
function Base.inv(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("Inverse of a tensor only exist when domain == codomain; check pinv"))
    if sectortype(t) === Trivial
        return TensorMap(inv(block(t, Trivial())), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c,b) in blocks(t)
            data[c] = inv(b)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end
function LinearAlgebra.pinv(t::AbstractTensorMap; kwargs...)
    if sectortype(t) === Trivial
        return TensorMap(pinv(block(t, Trivial()); kwargs...), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c,b) in blocks(t)
            data[c] = pinv(b; kwargs...)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end
function Base.:(\)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    codomain(t1) == codomain(t2) ||
        throw(SpaceMismatch("non-matching codomains in t1 \\ t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) \ block(t2, Trivial())
        return TensorMap(data, domain(t1)←domain(t2))
    else
        cod = codomain(t1)
        data = SectorDict(c=>block(t1,c) \ block(t2,c) for c in blocksectors(codomain(t1)))
        return TensorMap(data, domain(t1)←domain(t2))
    end
end
function Base.:(/)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    domain(t1) == domain(t2) ||
        throw(SpaceMismatch("non-matching domains in t1 / t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) / block(t2, Trivial())
        return TensorMap(data, codomain(t1)←codomain(t2))
    else
        data = SectorDict(c=>block(t1,c) / block(t2,c) for c in blocksectors(domain(t1)))
        return TensorMap(data, codomain(t1)←codomain(t2))
    end
end

# TensorMap exponentation:
function exp!(t::TensorMap)
    domain(t) == codomain(t) ||
        error("Exponentional of a tensor only exist when domain == codomain.")
    for (c,b) in blocks(t)
        copyto!(b, LinearAlgebra.exp!(b))
    end
    return t
end

# # concatenate tensors
# function concatenate(t1::AbstractTensorMap{S}, ts::AbstractTensorMap{S}...; direction::Symbol) where {S}
#     if direction = :domain
#         numin(t1) == 1 || throw(SpaceMismatch("concatenation along domain"))
#         for t2 in ts
#             domain(t1) == domain(t2) || throw(SpaceMismatch())
#         end
#         t = t1
#         for t2 in ts
#             cat(t1)

function catdomain(t1::AbstractTensorMap{S,N₁,1}, t2::AbstractTensorMap{S,N₁,1}) where {S,N₁}
    codomain(t1) == codomain(t2) || throw(SpaceMismatch())

    V1, = domain(t1)
    V2, = domain(t2)
    isdual(V1) == isdual(V2) || throw(SpaceMismatch("cannot hcat tensors whose domain has non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1), eltype(t2)), codomain(t1), V)
    for c in sectors(V)
        block(t, c)[:,1:dim(V1, c)] .= block(t1, c)
        block(t, c)[:,dim(V1, c) .+ (1:dim(V2, c))] .= block(t2, c)
    end
    return t
end
function catcodomain(t1::AbstractTensorMap{S,1,N₂}, t2::AbstractTensorMap{S,1,N₂}) where {S,N₂}
    domain(t1) == domain(t2) || throw(SpaceMismatch())

    V1, = codomain(t1)
    V2, = codomain(t2)
    isdual(V1) == isdual(V2) || throw(SpaceMismatch("cannot hcat tensors whose domain has non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1),eltype(t2)), V, domain(t1))
    for c in sectors(V)
        block(t, c)[1:dim(V1, c),:] .= block(t1, c)
        block(t, c)[dim(V1,c) .+ (1:dim(V2, c)),:] .= block(t2, c)
    end
    return t
end

# tensor product of tensors
function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where S
    cod1, cod2 = codomain(t1), codomain(t2)
    dom1, dom2 = domain(t1), domain(t2)
    cod = cod1 ⊗ cod2
    dom = dom1 ⊗ dom2
    t = TensorMap(zeros, promote_type(eltype(t1), eltype(t2)), cod, dom)
    if sectortype(S) === Trivial
        d1 = dim(cod1)
        d2 = dim(cod2)
        d3 = dim(dom1)
        d4 = dim(dom2)
        m1 = reshape(t1[], (d1, 1, d3, 1))
        m2 = reshape(t2[], (1, d2, 1, d4))
        m = reshape(t[], (d1, d2, d3, d4))
        m .= m1 .* m2
    else
        for (f1l, f1r) in fusiontrees(t1)
            for (f2l, f2r) in fusiontrees(t2)
                for (fl, coeff1) in merge(f1l, f2l)
                    for (fr, coeff2) in merge(f1r, f2r)
                        fl.coupled == fr.coupled || continue
                        d1 = dim(cod1, f1l.uncoupled)
                        d2 = dim(cod2, f2l.uncoupled)
                        d3 = dim(dom1, f1r.uncoupled)
                        d4 = dim(dom2, f2r.uncoupled)
                        m1 = reshape(t1[f1l,f1r], (d1, 1, d3, 1))
                        m2 = reshape(t2[f2l,f2r], (1, d2, 1, d4))
                        m = reshape(t[fl, fr], (d1, d2, d3, d4))
                        m .+= coeff1 .* coeff2 .* m1 .* m2
                    end
                end
            end
        end
    end
    return t
end