text:
####
# This macro overrides mul! to call lazymul!
####

# support mul! by calling lazy mul
macro lazymul(Typ)
    ret = quote
        LinearAlgebra.mul!(dest::AbstractVector, A::$Typ, b::AbstractVector) =
            copyto!(dest, LazyArrays.Mul(A,b))

        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::AbstractMatrix) =
            copyto!(dest, LazyArrays.Mul(A,b))
        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Typ) =
            copyto!(dest, LazyArrays.Mul(A,b))

        Base.:*(A::$Typ, B::$Typ, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B,C...))
        Base.:*(A::$Typ, B::AbstractMatrix, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B,C...))
        Base.:*(A::$Typ, B::AbstractVector) = LazyArrays.materialize(LazyArrays.Mul(A,B))
        Base.:*(A::AbstractMatrix, B::$Typ, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B,C...))
    end
    for Struc in (:AbstractTriangular, :Diagonal)
        ret = quote
            $ret

            Base.:*(A::$Struc, B::$Typ, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::$Typ, B::$Struc, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
        end
    end
    for Mod in (:Adjoint, :Transpose, :Symmetric, :Hermitian)
        ret = quote
            $ret

            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Mod{<:Any,<:AbstractMatrix}) =
                copyto!(dest, LazyArrays.Mul(A,b))

            LinearAlgebra.mul!(dest::AbstractVector, A::$Mod{<:Any,<:$Typ}, b::AbstractVector) =
                copyto!(dest, LazyArrays.Mul(A,b))

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Mod{<:Any,<:$Typ}, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::AbstractMatrix, B::$Mod{<:Any,<:$Typ}, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractVector, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Typ, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::$Typ, B::$Mod{<:Any,<:$Typ}, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::Diagonal, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::Diagonal, B::$Mod{<:Any,<:$Typ}, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))

            Base.:*(A::AbstractTriangular, B::$Mod{<:Any,<:$Typ}, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractTriangular, C...) = LazyArrays.materialize(LazyArrays.Mul(A,B, C...))
        end
    end

    esc(ret)
end

macro lazylmul(Typ)
    esc(quote
        LinearAlgebra.lmul!(A::$Typ, x::AbstractVector) = copyto!(x, LazyArrays.Mul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::AbstractMatrix) = copyto!(x, LazyArrays.Mul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::StridedVector) = copyto!(x, LazyArrays.Mul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::StridedMatrix) = copyto!(x, LazyArrays.Mul(A,x))
    end)
end

@lazymul AddMatrix


text:
using StaticArrays

const Mark = SVector{2, Float64}

# disk(r::Int) = [CartesianIndex(y,x) for y in -r:r for x in -r:r if sqrt(y^2 + x^2) ≤ r]

const sz = 1024
const intensity_radius = 2 # the radius of the disk that is used for the intensity
const weight_radius = 5 # the radius of the disk that is used for the tracking
const border = max(intensity_radius, weight_radius) + 1

outside(i::Float64) = i ≤ 1 + border || i ≥ sz - border
outside(p::T) where T <: AbstractVector = any(outside(i) for i in p)

# mark2cart(p::Mark) = CartesianIndex(round.(Int, p)...)

ext: jl


text:
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["lines.linewidth"] = 2. # Change linewidth

function fill_past(m, ylims)
    domain_idx = (t(m) .> m.present_year)
    fill_between(
        t(m)[.~domain_idx],
        ones(size(t(m)[.~domain_idx])) * ylims[1] * 2.,
        ones(size(t(m)[.~domain_idx])) * ylims[2] * 2.,
        facecolor="b", alpha=0.1
    )
    ylim(ylims)
    return
end

function plot_emissions(m::ClimateModel)
    title("effective greenhouse gas emissions")
    #fill_past(m, ylims)

    fill_between(t(m), effective_emissions(m), effective_emissions(m, M=true), facecolor="C0", alpha=0.3, label="Mitigation")
    fill_between(t(m), effective_emissions(m, M=true), effective_emissions(m, M=true, R=true), facecolor="C1", alpha=0.3, label="CDR")
    plot(t(m), effective_emissions(m), "-", color="grey", lw=2.0, label=L"$rq$ (no-policy baseline)")
    plot(t(m), effective_emissions(m, M=true), "k-", lw=1, alpha=0.4)
    plot(t(m), effective_emissions(m, M=true, R=true), "k-", lw=2., label=L"$rq(1-M) - rq_{0}R$ (controlled)")
    plot(t(m), zeros(size(t(m))), dashes=(2.5, 1.75), color="grey", alpha=0.5)

    ylimit = maximum(effective_emissions(m)) * 1.1
    ylims = [-ylimit, ylimit]
    ylabel(L"effective CO$_{2e}$ emissions [ppm / yr]")
    xlim(t(m)[1],2200.)
    ylim(minimum(effective_emissions(m, M=true, R=true))-5.,1.1*maximum(effective_emissions(m)))
    xticks(t(m)[1]:40.:2200.)
    xlabel("year")
    grid(true, alpha=0.3)
    return
end

function plot_concentrations(m::ClimateModel)
    title("greenhouse gas concentrations")
    #fill_past(m, ylims)

    fill_between(t(m), c(m), c(m, M=true), facecolor="C0", alpha=0.3, label="Mitigation")
    fill_between(t(m), c(m, M=true), c(m, M=true, R=true), facecolor="C1", alpha=0.3, label="CDR")
    plot(t(m), c(m), "-", color="grey", lw=2.0, label=L"$c$ (no-policy baseline)")
    plot(t(m), c(m, M=true), "k-", lw=1, alpha=0.4)
    plot(t(m), c(m, M=true, R=true), "k-", lw=2.0, label=L"$c_{M,R}$ (controlled)")
    ylims = [0., maximum(c(m))*1.05]
    ylabel(L"CO$_{2e}$ concentration [ppm]")
    xlabel("year")
    xlim(t(m)[1],2200.)
    ylim(ylims)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.3)
    return
end

function plot_forcings(m::ClimateModel; F0=3.)
    title("forcing (greenhouse effect and SRM)")
    #fill_past(m, ylims)

    fill_between(t(m), F(m, F0=F0), F(m, M=true, F0=F0), facecolor="C0", alpha=0.3, label="Mitigation")
    fill_between(t(m), F(m, M=true, F0=F0), F(m, M=true, R=true, F0=F0), facecolor="C1", alpha=0.3, label="CDR")
    fill_between(t(m), F(m, M=true, R=true, F0=F0), F(m, M=true, R=true, G=true, F0=F0), facecolor="C3", alpha=0.3, label="SRM")
    plot(t(m), F(m, F0=F0), "-", color="grey", lw=2.0, label=L"$F$ (no-policy baseline)")
    plot(t(m), F(m, M=true, F0=F0), "k-", lw=1, alpha=0.4)
    plot(t(m), F(m, M=true, R=true, F0=F0), "k-", lw=1, alpha=0.4)
    plot(t(m), F(m, M=true, R=true, G=true, F0=F0), "k-", lw=2.0, label=L"$F_{M,R,G}$ (controlled)")
    ylims = [0., maximum(F(m, F0=F0))*1.05]
    ylabel(L"radiative forcing [W/m$^{2}$]")
    xlabel("year")
    xlim(t(m)[1],2200.)
    ylim(ylims)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.3)
    return
end

function plot_temperatures(m::ClimateModel; temp_goal=1.2)
    title("adaptive temperature change")
    #fill_past(m, ylims)

    fill_between(t(m), T_adapt(m), T_adapt(m, M=true), facecolor="C0", alpha=0.3, label="Mitigation")
    fill_between(t(m), T_adapt(m, M=true), T_adapt(m, M=true, R=true), facecolor="C1", alpha=0.3, label="CDR")
    fill_between(t(m), T_adapt(m, M=true, R=true), T_adapt(m, M=true, R=true, G=true), facecolor="C3", alpha=0.3, label="SRM")
    fill_between(t(m), T_adapt(m, M=true, R=true, G=true), T_adapt(m, M=true, R=true, G=true, A=true), facecolor="C2", alpha=0.3, label="Adaptation")
    plot(t(m), T_adapt(m), "-", color="grey", lw=2.0, label=L"$T$ (no-policy baseline)")
    plot(t(m), T_adapt(m, M=true), "k-", lw=1, alpha=0.4)
    plot(t(m), T_adapt(m, M=true, R=true), "k-", lw=1, alpha=0.4)
    plot(t(m), T_adapt(m, M=true, R=true, G=true), "k-", lw=1., alpha=0.4)
    plot(t(m), T_adapt(m, M=true, R=true, G=true, A=true), "k-", lw=2.0, label=L"$T_{M,R,G,A}$ (adaptive)")
    plot(t(m),temp_goal .* ones(size(t(m))), dashes=(2.5, 1.75), color="grey", alpha=0.75, lw=2.5)
    ylims = [0., maximum(T_adapt(m)) * 1.05]
    ylabel("temperature anomaly [°C]")
    xlabel("year")
    xlim(t(m)[1],2200.)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.3)
    legend()
    return
end

function plot_controls(m::ClimateModel)
    title("optimized control deployments")
    plot(t(m)[m.economics.baseline_emissions .> 0.], m.controls.mitigate[m.economics.baseline_emissions .> 0.],
        color="C0", lw=2, label=L"$M$ (emissions mitigation)")
    plot(t(m), m.controls.remove, color="C1", lw=2, label=L"$R$ (carbon dioxide removal)")
    plot(t(m), m.controls.adapt, color="C2", lw=2, label=L"$A$ (adaptation)")
    plot(t(m), m.controls.geoeng, color="C3", lw=2, label=L"$G$ (solar geoengineering)")
    ylims = [0., 1.075]
    yticks(0.:0.2:1.0, 0:20:100)
    ylim(ylims)
    ylabel("control deployment [%]")
    xlabel("year")
    xlim(t(m)[1],2200.)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.3)
    return
end

function plot_benefits(m::ClimateModel; discounting=true)
    domain_idx = (t(m) .> m.domain.present_year)

    fill_between(
        t(m)[domain_idx],
        0 .*ones(size(t(m)))[domain_idx],
        net_benefit(m, discounting=discounting, M=true, R=true, G=true, A=true)[domain_idx],
        facecolor="grey", alpha=0.2
    )
    plot(t(m)[domain_idx], 0 .*ones(size(t(m)))[domain_idx], "-", lw=2, color="gray", label="no-policy baseline")
    plot(t(m)[domain_idx], benefit(m, discounting=discounting, M=true, R=true, G=true, A=true)[domain_idx], color="C8", lw=2, label="benefits (of avoided damages)")
    plot(t(m)[domain_idx], cost(m, discounting=discounting, M=true, R=true, G=true, A=true)[domain_idx], color="C4", lw=2, label="costs (of climate controls)")
    plot(t(m)[domain_idx], net_benefit(m, discounting=discounting, M=true, R=true, G=true, A=true)[domain_idx], color="k", lw=2, label="net benefits (benefits - costs)")
    ylabel(L"discounted costs and benefits [10$^{12}$ \$ / year]")
    xlabel("year")
    xlim(t(m)[1],2200.)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.3)
    title("cost-benefit analysis")
    return
end

function plot_damages(m::ClimateModel; discounting=true, percent_GWP=false, temp_goal=1.2)
    Enorm = deepcopy(E(m))/100.
    if ~percent_GWP; Enorm=1.; end;

    domain_idx = (t(m) .> m.domain.present_year)
    fill_between(
        t(m)[domain_idx],
        0 .*ones(size(t(m)))[domain_idx],
        (cost(m, discounting=discounting, M=true, R=true, G=true, A=true) ./ Enorm)[domain_idx],
        facecolor="C4", alpha=0.2
    )
    damages = damage(m, discounting=discounting, M=true, R=true, G=true, A=true)
    costs = cost(m, discounting=discounting, M=true, R=true, G=true, A=true)
    plot(t(m)[domain_idx], (damage(m, discounting=discounting) ./ Enorm)[domain_idx], color="gray", lw=2.0, label="uncontrolled damages")
    plot(t(m)[domain_idx], ((damages .+ costs)./ Enorm)[domain_idx], color="k", lw=2.0, label="net costs (controlled damages + controls)")
    plot(t(m)[domain_idx], (damages ./ Enorm)[domain_idx], color="C8", lw=2.0, label="controlled damages")
    plot(t(m)[domain_idx], (costs ./ Enorm)[domain_idx], color="C4", lw=2.0, label="cost of controls")

    ylim([0, maximum((damage(m, discounting=discounting) ./ Enorm)[domain_idx]) * 0.75])

    dmg_label = string("damage threshold at ",round(temp_goal, digits=2),L"°C with $A=0$")
    plot(
        t(m)[domain_idx],
        (damage(m.economics.β, E(m), temp_goal, discount=discount(m)) ./ Enorm)[domain_idx],
        dashes=(2.5,1.75), color="grey", alpha=0.75, lw=2.0, label=dmg_label
    )

    if ~percent_GWP;
        if ~discounting;
            ylabel(L"costs [10$^{12}$ \$ / year]");
        else;
            ylabel(L"discounted costs [10$^{12}$ \$ / year]");
        end
    else
        if ~discounting
            ylabel("costs [% GWP]")
        else
            ylabel("discounted costs [% GWP]")
            print("NOT YET SUPPORTED")
        end
    end

    xlabel("year")
    xlim(t(m)[1],2200.)
    xticks(t(m)[1]:40.:2200.)
    grid(true, alpha=0.175)
    title("costs of avoiding a damage threshold")
    legend()
    return
end

function plot_state(m::ClimateModel; new_figure=true, plot_legends=true, temp_goal=1.2)
    if new_figure
        fig, axs = subplots(2, 3, figsize=(14,8))
        axs = vcat(axs...)
    end

    sca(axs[1])
    plot_emissions(m)
    title("a)", loc="left")
    sca(axs[2])
    plot_concentrations(m)
    title("b)", loc="left")
    sca(axs[3])
    plot_temperatures(m, temp_goal=temp_goal)
    title("c)", loc="left")

    sca(axs[4])
    plot_controls(m)
    title("d)", loc="left")
    sca(axs[5])
    plot_benefits(m)
    title("e)", loc="left")
    sca(axs[6])
    plot_damages(m, temp_goal=temp_goal)
    title("f)", loc="left")

    if plot_legends;
        for ii in 1:6
            sca(axs[ii])
            if ii <= 2;
                legend(loc="lower left", labelspacing=0.14, handlelength=1.75);
            else
                legend(loc="upper left", labelspacing=0.14, handlelength=1.75);
            end
        end
    end
    tight_layout()

    return fig, axs
end


text:
function from_ggm(d::AbstractString; tolowercase=false)
  str = replace(strip(d), " " => "")
  str = replace(str, "\n" => "")
  startp = findall("(", str)
  endp = findall(")", str)
  orderp = findall("order", str)
  if length(orderp) == 0
    f = split(str[startp[1].start+1:endp[end].stop-1], ",")
  else
    f = split(str[startp[1].start+1:orderp[1].start-2], ",")
  end
  dct = OrderedDict{Symbol, SymbolList}()
  for s in f
    if tolowercase
      lhs = Symbol(lowercase(split(s, "~")[1]))
      rhs = Symbol.(lowercase.(split(split(s, "~")[2], "+")))
    else
      lhs = Symbol(split(s, "~")[1])
      rhs = Symbol.(split(split(s, "~")[2], "+"))
    end
    rhs = length(rhs) == 1 ? rhs[1] : rhs
    dct[lhs] = rhs
  end
  dct
end

function to_ggm(d::OrderedDict; touppercase=false, order=false)
  str = "DAG("
  for (ind, key) in enumerate(keys(d))
    str = ind > 1 ? "$(str), " : "$(str)"
    rhs = d[key]
    if typeof(rhs) == Symbol
      rhs_str = String(rhs)
    elseif length(rhs) == 0
      # Not used yet
    elseif length(rhs) == 1
      rhs_str = String(rhs[1])
    elseif length(rhs) > 1
      rhs_str = "$(String(rhs[1]))"
      for s in rhs[2:end]
        rhs_str = "$(rhs_str) + $(String(s))"
      end
    end
    # Split multiple lhs symbols (if present)
    if typeof(key) == Symbol
      str = "$(str)$(String(key)) ~ $rhs_str"
    elseif length(key) == 1
      str = "$(str)$(String(key[1])) ~ $rhs_str"
    elseif length(key) > 1
      str = "$(str)$(String(key[1])) ~ $(rhs_str), "
      for s in key[2:end-1]
        str = "$(str)$(String(s)) ~ $(rhs_str)"
      end
      str = "$(str)$(String(key[end])) ~ $(rhs_str)"
    end
  end
  str = "$(str))"
  str = touppercase ? uppercase(str) : str
  order ? "$(str[1:end-1]), order=TRUE)" : "$(str)"
end

function to_ggm(d::DAG; touppercase=false, order=false)
  to_ggm(d.d; touppercase=touppercase, order=order)
end

export
  to_ggm,
  from_ggm


text:
# kernel execution

export AbstractGPUBackend, gpu_call, synchronize, thread_blocks_heuristic

abstract type AbstractGPUBackend end

backend(::Type{T}) where T = error("Can't choose GPU backend for $T")

"""
    gpu_call(kernel::Function, A::AbstractGPUArray, args::Tuple, configuration = length(A))

Calls function `kernel` on the GPU.
`A` must be an AbstractGPUArray and will help to dispatch to the correct GPU backend
and supplies queues and contexts.
Calls the kernel function with `kernel(state, args...)`, where state is dependant on the backend
and can be used for getting an index into `A` with `linear_index(state)`.
Optionally, a launch configuration can be supplied in the following way:

    1) A single integer, indicating how many work items (total number of threads) you want to launch.
        in this case `linear_index(state)` will be a number in the range `1:configuration`
    2) Pass a tuple of integer tuples to define blocks and threads per blocks!

"""
function gpu_call(kernel, A::AbstractArray, args::Tuple, configuration = length(A))
    ITuple = NTuple{N, Integer} where N
    # If is a single integer, we assume it to be the global size / total number of threads one wants to launch
    thread_blocks = if isa(configuration, Integer)
        thread_blocks_heuristic(configuration)
    elseif isa(configuration, ITuple)
        # if a single integer ntuple, we assume it to configure the blocks
        configuration,  ntuple(x-> x == 1 ? 256 : 1, length(configuration))
    elseif isa(configuration, Tuple{ITuple, ITuple})
        # 2 dim tuple of ints == blocks + threads per block
        if any(x-> length(x) > 3 || length(x) < 1, configuration)
            error("blocks & threads must be 1-3 dimensional. Found: $configuration")
        end
        map(x-> Int.(x), configuration) # make sure it all has the same int type
    else
        error("""Please launch a gpu kernel with a valid configuration.
            Found: $configurations
            Configuration needs to be:
            1) A single integer, indicating how many work items (total number of threads) you want to launch.
                in this case `linear_index(state)` will be a number in the range 1:configuration
            2) Pass a tuple of integer tuples to define blocks and threads per blocks!
                `linear_index` will be inbetween 1:prod((blocks..., threads...))
        """)
    end
    _gpu_call(backend(typeof(A)), kernel, A, args, thread_blocks)
end

# Internal GPU call function, that needs to be overloaded by the backends.
_gpu_call(::Any, f, A, args, thread_blocks) = error("Not implemented") # COV_EXCL_LINE

"""
    synchronize(A::AbstractArray)

Blocks until all operations are finished on `A`
"""
function synchronize(A::AbstractArray)
    # fallback is a noop, for backends not needing synchronization. This
    # makes it easier to write generic code that also works for AbstractArrays
end

function thread_blocks_heuristic(len::Integer)
    # TODO better threads default
    threads = clamp(len, 1, 256)
    blocks = max(ceil(Int, len / threads), 1)
    (blocks,), (threads,)
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x = y = -15:0.33:15
fz(x,y) = sin.(sqrt.(x.^2 + y.^2))./sqrt.(x.^2+y.^2)
fxy = [fz(x,y) for x in x, y in y]
@gsp " " :-  # reset session
@gsp(x, y, fxy, "w l lc palette", "set view map",
    "set contour base;set key off","set auto fix",
    "set cntrparam levels 15", "unset surface",
    xlab = "x", ylab = "y")
save(term="pngcairo size 800,800", output="plt3d_ex4.png")


text:
struct Rectangle
    height::Int64
    length::Int64
end

function is_square(self::Rectangle)::Bool
    return self.height == self.length
end


text: positivesum(arr::Array{Int}) = isempty(arr) ? 0 : sum( x -> max(0,x), arr )

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#####
##### 'Tupled closure' implementation: 1-tuple, 2-tuple, and then n-tuple by induction
#####

function closure_summary(closures::Tuple, padchar="│")
    Nclosures = length(closures)
    if Nclosures == 1
        return string("Tuple with 1 closure:", '\n',
                      "$padchar   └── $(dict.keys[1]) => $(typeof(dict.vals[1]).name)")
    else
        return string("Tuple with $Nclosures closures:", '\n',
         Tuple(string("$padchar   ├── ", summary(c), '\n') for c in closures[1:end-1])...,
                      "$padchar   └── ", summary(closures[end]))
    end
end

#####
##### Stress divergences
#####

for stress_div in (:∂ⱼ_τ₁ⱼ, :∂ⱼ_τ₂ⱼ, :∂ⱼ_τ₃ⱼ)
    @eval begin
        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any}, clock, U, Ks, args...) =
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any}, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
                  + $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...))

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any, <:Any}, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
                  + $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...)
                  + $stress_div(i, j, k, grid, closures[3], clock, U, Ks[3], args...))

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1:2], clock, U, Ks[1:2], args...)
                  + $stress_div(i, j, k, grid, closures[3:end], clock, U, Ks[3:end], args...))
    end
end

#####
##### Tracer flux divergences
#####

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any}, c, iᶜ, clock, Ks, args...) =
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any}, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[2], c, iᶜ, clock, Ks[2], args...))

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any, <:Any}, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[2], c, iᶜ, clock, Ks[2], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[3], c, iᶜ, clock, Ks[3], args...))

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1:2], c, iᶜ, clock, Ks[1:2], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[3:end], c, iᶜ, clock, Ks[3:end], args...))

#####
##### Utilities
#####

with_tracers(tracers, closure_tuple::Tuple) =
    Tuple(with_tracers(tracers, closure) for closure in closure_tuple)

function calculate_diffusivities!(diffusivity_fields_tuple, closure_tuple::Tuple, args...)
    for (α, closure) in enumerate(closure_tuple)
        @inbounds diffusivity_fields = diffusivity_fields_tuple[α]
        calculate_diffusivities!(diffusivity_fields, closure, args...)
    end
    return nothing
end

function add_closure_specific_boundary_conditions(closure_tuple::Tuple, bcs, args...)
    # So the last closure in the tuple has the say...
    for closure in closure_tuple
        bcs = add_closure_specific_boundary_conditions(closure, bcs, args...)
    end
    return bcs
end

#####
##### Support for VerticallyImplicit
#####

const SingleImplicitClosure = AbstractTurbulenceClosure{<:VerticallyImplicitTimeDiscretization}
const VIC = Union{AbstractTurbulenceClosure{<:VerticallyImplicitTimeDiscretization}, AbstractArray{<:SingleImplicitClosure}}

# Filter explicitly-discretized closures.
@inline z_diffusivity(closure::Tuple, iᶜ, Kstuple, args...) =
        Tuple((clo isa VIC ? z_diffusivity(clo, iᶜ, Ks, args...) : 0) for (clo, Ks) in zip(closure, Kstuple))

@inline z_viscosity(closure::Tuple, Kstuple, args...) =
        Tuple((clo isa VIC ? z_viscosity(clo, Ks, args...) : 0) for (clo, Ks) in zip(closure, Kstuple))

for coeff in (:νᶜᶜᶜ, :νᶠᶠᶜ, :νᶠᶜᶠ, :νᶜᶠᶠ, :κᶜᶜᶠ, :κᶜᶠᶜ, :κᶠᶜᶜ)
    @eval begin
        @inline $coeff(i, j, k, grid, clock, ν::Tuple{C1})     where C1       = $coeff(i, j, k, grid, clock, ν[1])
        @inline $coeff(i, j, k, grid, clock, ν::Tuple{C1, C2}) where {C1, C2} = $coeff(i, j, k, grid, clock, ν[1])   + $coeff(i, j, k, grid, clock, ν[2])
        @inline $coeff(i, j, k, grid, clock, ν::Tuple)                        = $coeff(i, j, k, grid, clock, ν[1:2]) + $coeff(i, j, k, grid, clock, ν[3:end])
    end
end

@inline combine_time_discretizations(disc) = disc

@inline combine_time_discretizations(::ExplicitTimeDiscretization, ::VerticallyImplicitTimeDiscretization)           = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VerticallyImplicitTimeDiscretization, ::ExplicitTimeDiscretization)           = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VerticallyImplicitTimeDiscretization, ::VerticallyImplicitTimeDiscretization) = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::ExplicitTimeDiscretization, ::ExplicitTimeDiscretization)                     = ExplicitTimeDiscretization()

@inline combine_time_discretizations(disc1, disc2, other_discs...) =
    combine_time_discretizations(combine_time_discretizations(disc1, disc2), other_discs...)

@inline time_discretization(closures::Tuple) = combine_time_discretizations(time_discretization.(closures)...)


text:
using LanguageServer, CSTParser, Test, SymbolServer
const LS = LanguageServer
const Range = LanguageServer.Range

@testset "LanguageServer" begin

include("test_document.jl")
include("test_communication.jl")
include("test_hover.jl")
include("text_edit.jl")

end


text:
@testset "AllDifferentExceptConstants2Reification: $(fct_type), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    base_model = if T == Int
        IntReificationEqualToModel{T}()
    elseif T == Float64
        FloatReificationEqualToModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.AllDifferentExceptConstants2Reification{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.AllDifferentExceptConstants{T},
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end
    values_set = Set([zero(T), one(T)])
    c = MOI.add_constraint(model, fct, CP.AllDifferentExceptConstants(dim, values_set))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.AllDifferentExceptConstants}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.AllDifferentExceptConstants) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
            (MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == dim * length(values_set) + dim * (dim - 1) / 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == dim * length(values_set) + dim * (dim - 1) / 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}()) == dim * length(values_set)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}}()) == dim * (dim - 1) / 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == dim * (dim - 1) / 2

        @test Set(MOI.get(bridge, MOI.ListOfVariableIndices())) == Set([vec(collect(values(bridge.vars_compare)))..., vec(collect(values(bridge.vars_different)))...])
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}())) == Set([vec(collect(values(bridge.vars_compare_bin)))..., vec(collect(values(bridge.vars_different_bin)))...])
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}())) == Set(vec(bridge.cons_compare_reif))
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}}())) == Set(collect(values(bridge.cons_different_reif)))
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}())) == Set(collect(values(bridge.cons)))
    end

    @testset "Set of variables" begin
        @test length(bridge.vars_compare) == dim * length(values_set)
        @test length(bridge.vars_compare_bin) == dim * length(values_set)

        for i in 1:dim
            for j in 1:2
                @test MOI.is_valid(model, bridge.vars_compare[i, j])
                @test MOI.is_valid(model, bridge.vars_compare_bin[i, j])
            end
            for j in 3:dim
                @test (i, j) ∉ keys(bridge.vars_compare)
                @test (i, j) ∉ keys(bridge.vars_compare_bin)
            end
        end

        @test length(bridge.vars_different) == dim * (dim - 1) / 2
        @test length(bridge.vars_different_bin) == dim * (dim - 1) / 2

        for i in 1:dim
            for j in 1:i
                @test (i, j) ∉ keys(bridge.vars_different)
                @test (i, j) ∉ keys(bridge.vars_different_bin)
            end
            for j in (i + 1):dim
                @test MOI.is_valid(model, bridge.vars_different[i, j])
                @test MOI.is_valid(model, bridge.vars_different_bin[i, j])
            end
        end
    end

    @testset "Compare array values to excluded values" begin
        @test length(bridge.cons_compare_reif) == dim * length(values_set)

        for i in 1:dim
            for j in 1:length(values_set)
                @test MOI.is_valid(model, bridge.cons_compare_reif[i, j])
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_compare_reif[i, j]) == CP.Reification(MOI.EqualTo(zero(T)))
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_compare_reif[i, j])
                @test length(f.terms) == 2
                @test f.constants == [zero(T), -T(j - 1)]

                t1 = f.terms[1]
                @test t1.output_index == 1
                @test t1.scalar_term.coefficient === one(T)
                @test t1.scalar_term.variable == bridge.vars_compare[i, j]

                t2 = f.terms[2]
                @test t2.output_index == 2
                @test t2.scalar_term.coefficient === one(T)
                @test t2.scalar_term.variable == x[i]
            end
        end
    end

    @testset "Compare array variables" begin
        @test length(bridge.cons_different_reif) == dim * (dim - 1) / 2

        for i in 1:dim
            for j in (i+1):dim
                @test MOI.is_valid(model, bridge.cons_different_reif[i, j])
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_different_reif[i, j]) == CP.Reification(CP.DifferentFrom(zero(T)))
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_different_reif[i, j])
                @test length(f.terms) == 3

                t1 = f.terms[1]
                @test t1.output_index == 1
                @test t1.scalar_term.coefficient === one(T)
                @test t1.scalar_term.variable == bridge.vars_different[i, j]

                t2 = f.terms[2]
                @test t2.output_index == 2
                @test t2.scalar_term.coefficient === one(T)
                @test t2.scalar_term.variable == x[i]

                t3 = f.terms[3]
                @test t3.output_index == 2
                @test t3.scalar_term.coefficient === -one(T)
                @test t3.scalar_term.variable == x[j]
            end
        end
    end

    @testset "Disjunction" begin
        @test length(bridge.cons) == dim * (dim - 1) / 2

        for i in 1:dim
            for j in (i+1):dim
                @test MOI.is_valid(model, bridge.cons[i, j])
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons[i, j]) == MOI.GreaterThan(one(T))
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons[i, j])
                @test length(f.terms) == 1 + 2 * length(values_set)

                for k in 1:length(values_set)
                    t1 = f.terms[k]
                    @test t1.coefficient === one(T)
                    @test t1.variable == bridge.vars_compare[i, k]

                    t2 = f.terms[length(values_set) + k]
                    @test t2.coefficient === one(T)
                    @test t2.variable == bridge.vars_compare[j, k]
                end

                t = f.terms[2 * length(values_set) + 1]
                @test t.coefficient === one(T)
                @test t.variable == bridge.vars_different[i, j]
            end
        end
    end
end


text:
using JLD2, Dierckx, FFTW, StatsBase, Infiltrator

#Construct Cell
Cell = Construct("LG M50")

#Arrhenius
Cell.Const.T = 298.15
Arr_Factor = (1/Cell.Const.T_ref-1/Cell.Const.T)/R

#Set Cell Constants
Cell.Const.SOC = 0.8
Cell.Const.κ = Cell.Const.κf(Cell.Const.ce0)*exp(Cell.Const.Ea_κ*Arr_Factor)
Cell.RA.Nfft = Cell.RA.Nfft!(Cell.RA.Fs, Cell.RA.Tlen)
Cell.RA.f = Cell.RA.f!(Cell.RA.Nfft)
Cell.RA.s = Cell.RA.s!(Cell.RA.Fs,Cell.RA.Nfft,Cell.RA.f)

tfft = @. (1/Cell.RA.Fs)*Cell.RA.f
OrgT = @. Cell.RA.SamplingT*(0:floor(tfft[end]/Cell.RA.SamplingT))
puls = Array{Float64}(undef,Cell.RA.Outs,size(OrgT,1)-1)
D = Vector{Float64}(undef,Cell.RA.Outs)
C_Aug = Vector{Float64}(undef,Cell.RA.Outs)

function tf_new!(puls,D,C_Aug,s,f)
#Initialise Loop Variables
i=Int64(1)
l=Int64(1)
u=Int64(0)

    for run in Cell.Transfer.tfs[:,1]
        tf = Array{ComplexF64}(undef,size(Cell.Transfer.tfs[i,3],1),size(s,2))
        Di = Vector{Float64}(undef,size(Cell.Transfer.tfs[i,3],1))
        res0 = Vector{Float64}(undef,size(Cell.Transfer.tfs[i,3],1))

        if Cell.Transfer.tfs[i,2] == "Pos"
            run(Cell,s,Cell.Transfer.tfs[i,3],"Pos",tf,Di,res0)
        elseif Cell.Transfer.tfs[i,2] == "Neg"
            run(Cell,s,Cell.Transfer.tfs[i,3],"Neg",tf,Di,res0)
        else
            run(Cell,Cell.RA.s,Cell.Transfer.tfs[i,3],tf,Di,res0)
        end
        jk = Cell.RA.Fs*real(ifft(tf,2)) # inverse fourier transform tranfser function response (Large Compute)
        stpsum = (cumsum(jk, dims=2).*(1/Cell.RA.Fs)) #cumulative sum of tf response * sample time
        samplingtf = Array{Float64}(undef,size(stpsum,1),length(OrgT))

        #Interpolate H(s) to obtain h_s(s) to obtain discrete-time impulse response
        for Output in 1:size(stpsum,1)
            spl1 = Spline1D(tfft,stpsum[Output,:]; k=3)
            samplingtf[Output,:] .= evaluate(spl1,OrgT)
        end
        u += size(Cell.Transfer.tfs[i,3],1)
        #puls = [puls; diff(samplingtf, dims=2)]
        puls[l:u,:] .= diff(samplingtf, dims=2)
        D[l:u,:] .= Di
        C_Aug[l:u,:] .= res0
        l = u+1
        i += 1

    end
end

function tf(Cell,s,f)
    tfft = @. (1/Cell.RA.Fs)*Cell.RA.f
    OrgT = @. Cell.RA.SamplingT*(0:floor(tfft[end]/Cell.RA.SamplingT))
    i=Int64(1)
    puls = Array{Float64}(undef,0,size(OrgT,1)-1)
    D = Array{Float64}(undef,0,1)
    Dtt = Array{String}(undef,0,1)
    C_Aug = Array{Float64}(undef,0,1)
    for run in Cell.Transfer.tfs[:,1]
        if Cell.Transfer.tfs[i,2] == "Pos"
            tf, Di, res0, Dti = run(Cell,s,Cell.Transfer.tfs[i,3],"Pos")
        elseif Cell.Transfer.tfs[i,2] == "Neg"
            tf, Di, res0, Dti = run(Cell,s,Cell.Transfer.tfs[i,3],"Neg")
        else
            tf, Di, res0, Dti = run(Cell,s,Cell.Transfer.tfs[i,3])
        end

        jk = Cell.RA.Fs*real(ifft(tf,2)') # inverse fourier transform tranfser function response (Large Compute)
        stpsum = (cumsum(jk, dims=1).*(1/Cell.RA.Fs))' # cumulative sum of tf response * sample time
        samplingtf = Array{Float64}(undef,size(stpsum,1),length(OrgT))

        # Interpolate H(s) to obtain h_s(s) to obtain discrete-time impulse response
        for Output in 1:size(stpsum,1)
            spl1 = Spline1D(tfft,stpsum[Output,:]; k=3)
            samplingtf[Output,:]= evaluate(spl1,OrgT)
        end

        puls = [puls; diff(samplingtf, dims=2)]
        D = [D; Di]
        Dtt = [Dtt; Dti]
        C_Aug = [C_Aug; res0]
        i += 1

    end
    return puls
end

#Main Branch
#puls_ref = tf(Cell,Cell.RA.s,Cell.RA.f)

#Dev Branch
puls_ref = load("puls.jld2", "puls_ref")
tf_new!(puls,D,C_Aug,Cell.RA.s,Cell.RA.f)

maxdiff = Vector{Float64}(undef,size(puls_ref,1))
rmsediff = Vector{Float64}(undef,size(puls_ref,1))
for i in 1:size(puls_ref,1)
    maxdiff[i] = maxad(puls[i,:],puls_ref[i,:])
    rmsediff[i] = rmsd(puls[i,:],puls_ref[i,:])
end
#jldsave("puls.jld2"; puls_ref)

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
struct Wind
    id::Int
    pwm::UInt8
end
Wind(id::Int, ::Missing) = Wind(id, 0x00)
Wind(id::Int, pwm::String) = Wind(id, parse(Int, pwm))

tosecond(t::T) where {T <: TimePeriod}= t/convert(T, Second(1))
sincestart(t) = tosecond(t - t₀[])

function toint(msg)
    y = zero(UInt32)
    for c in msg
        y <<= 8
        y += c
    end
    return y
end

function getrpm(t)
    t4 = 15000000
    shortest_t = t4/1.1top_rpm[]
    t < shortest_t ?  missing : t4/t
end

function set_pwm!(sp, c, pwm)
    lock(c) do
        sp_flush(sp, SP_BUF_INPUT)
        encode(sp, pwm)
    end
end

function update_rpm!(sp, c, pwm, msg, rpm)
    if pwm[] > 19
        lock(c) do
            sp_flush(sp, SP_BUF_OUTPUT)
            decode!(msg, sp, 12)
        end
        for (i, x) in enumerate(Iterators.partition(msg, 4))
            rpm[i] = getrpm(toint(x))
        end
    else
        fill!(rpm, 0.0)
    end
end

function isconnected(port)
    try
        sp = LibSerialPort.open(port, baudrate[])
        close(sp)
        true
    catch ex
        false
    end
end

mutable struct FanArduino <: AbstractArduino
    id::Int
    c::ReentrantLock
    port::String
    sp::SerialPort
    msg::Vector{UInt8}
    rpm::Vector{Union{Missing, Float64}}
    pwm::Observable{UInt8}
    function FanArduino(id::Int, port::String)
        c = ReentrantLock()
        sp = LibSerialPort.open(port, baudrate[])
        pwm = Observable(0x00)
        on(pwm) do x
            set_pwm!(sp, c, x)
        end
        msg = Vector{UInt8}(undef, 12)
        rpm = Vector{Float64}(undef, 3)
        new(id, c, port, sp, msg, rpm, pwm)
    end
end

update_rpm!(a::FanArduino) = update_rpm!(a.sp, a.c, a.pwm, a.msg, a.rpm)

get_rpm(a::FanArduino) = a.rpm

mutable struct AllWind
    arduinos::Vector{FanArduino}
    io::IOStream
    framerate::Int
    function AllWind(arduinos::Vector{FanArduino}, framerate::Int)
        io = open(tempname(), "w")
        close(io)
        new(arduinos, io, framerate)
    end
end

function get_rpms(allwind::AllWind)
    @sync for a in allwind.arduinos
        @async update_rpm!(a)
    end
    now() => get_rpm.(allwind.arduinos)
end

function record(allwind::AllWind, folder)
    !isdir(folder) && mkpath(folder)
    isopen(allwind.io) && close(allwind.io)
    allwind.io = open(folder / "fans.csv", "w")
    println(allwind.io, "time,", join([join(["fan$(a.id)_speed$j" for j in 1:3], ",") for a in allwind.arduinos], ","))
    @async while isopen(allwind.io)
        println(allwind.io, now(), ",",join(Iterators.flatten(get_rpm.(allwind.arduinos)), ","))
        sleep(1/allwind.framerate)
    end
end

Base.isopen(allwind::AllWind) = all(isopen, allwind.arduinos)
Base.close(allwind::AllWind) = close.(allwind.arduinos)


text:
module LMM_module

using Statistics
using LinearAlgebra
using Distributions
using Optim
using RCall

#####################
### Sub-functions ###
#####################

"""
# _______________________________________
# Polyalgorithmic matrix inverse function

`inverse(A::Array{Float64,2})`

# Inputs
1. *A* [Array{Float64,2}]: square or rectangular matrix

# Output
1. Inverse or pseudo inverse of A

# Examples
```
S = convert(Array{Float64,2}, reshape(collect(1:9), 3, 3))
LMM_module.inverse(S)
R = convert(Array{Float64,2}, reshape(collect(1:10), 5, 2))
LMM_module.inverse(R)
```
"""
function inverse(A::Array{Float64,2})
	out = try
			LinearAlgebra.inv(A)
		catch
			LinearAlgebra.pinv(A)
		end
	return(out)
end

"""
# ________________________________________________________________
# Linear mixed model parameter optimization and effects estimation

`LMM_optim(parameters::Array{Float64,1}; X::Array{Float64,2}, y::Array{Float64,1}, Z::Array{Float64,2}, METHOD_VAR_EST::String=["ML", "REML"][1], OPTIMIZE::Bool=true, VARFIXEF::Bool=false)`

The mixed linear model is defined as y = Xb + Zu + e, where:
	- X [n,p] is the centered matrix of allele frequencies
	- Z [n,n] is the square symmetric matrix of relatedness
	- y [n,1] is the centered vector of phenotypic values
	- no intercept is explicitly fitted but implicitly set at the mean phenotypic value as a consequence of centering y
	- u ~ N(0, σ²uI)
	- e ~ N(0, σ²eI)
	- y ~ N(0, V); V = (Z (σ²uI) Z') + (σ²eI)
	- variance component (σ²e, σ²u) are estimated via maximum likelihood (ML) or restricted maximum likelihood (REML)
	- fixed effects (b) are estimated by solving: (X' * inverse(V)) * inverse(X * X' * inverse(V)) * y
	- random effects (u) are estimated by solving: (σ²uI) * Z' * inverse(V) * (y - (X*b))

# Inputs
1. parameters [Array{Float64,1}]: the error variance (σ2e) and random effects variance (σ2u)
2. *X* [Array{Float64,2}; xij ∈ [0.00, 1.00]]: matrix of allele frequencies (Pool-seq data) set as the fixed effects
3. *y* [Array{Float64,1}]: vector of phenotypes
4. *Z* [Array{Any,2}]: matrix of relatedness (Fst or XX'/p where p is the number of columns in X) set as the random effects
5. *METHOD_VAR_EST* [String]: method of variance parameter estimation choose "ML" for maximum likelihood or "REML" for restricted maximum likelihood (default="ML")
6. *OPTIMIZE* [Bool]: logical value to perform optimization to find estimate the variance parameters if TRUE or to estimate the fixed and random effects if FALSE (default=true)
7. *VARFIXEF* [Bool]: logical value to estimate the variances of the fixed effects which can be computationally and memory intensive for large datasets (default=false)

# Outputs
1. *If OPTIMIZE==true*: Negative log-likelihood of the parameters (σ2e and σ2u) given the dataset (X, y, and Z)
2. *If OPTIMIZE==false*: Fixed effects (b_hat), random effects (u_hat), and the optional fixed effects variances (σ2b_hat)

# Examples
```
using Distributions
using LinearAlgebra
using Optim
n=5; l=200; Ve=1.00; Vu=1.00;
X = reshape(abs.(rand(n*l)), n, l)
Z = reshape(abs.(rand(n^2)), n, n)
b = zeros(l); b[[5,10,100]] = [1.00, 1.00, 1.00]
u = rand(Distributions.Normal(0.00, Vu), n)
e = rand(Distributions.Normal(0.00, Ve), n)
y = (X*b) .+ (Z*u) .+ e
par = Optim.optimize(par -> LMM_module.LMM_optim(par, X=X, y=y, Z=Z), [0.0,0.0], [Inf,Inf], [1.00,1.00]).minimizer
b_hat, u_hat, Vb_hat = LMM_module.LMM_optim(par, X=X, y=y, Z=Z, OPTIMIZE=false)
```
"""
function LMM_optim(parameters::Array{Float64,1}; X::Array{Float64,2}, y::Array{Float64,1}, Z::Array{Float64,2}, METHOD_VAR_EST::String=["ML", "REML"][1], OPTIMIZE::Bool=true, VARFIXEF::Bool=false)
	### y = Xb + Zu + e; where Xb=fixed and Za=random
	σ2e = parameters[1] # variance of the error effects (assuming homoscedasticity)
	σ2u = parameters[2] # variance of the other random effects (assuming homoscedasticity)
	n = length(y)		# number of individual samples
	l = size(X, 2)		# number of fixed effects
	q = size(Z, 2)		# number of random effects excluding error effect
	### Variance components
	R = σ2e * I					# error variance-covariance matrix
	D = σ2u * I					# random effects variance-covariance matrix
	V = (Z * D * Z') + R		# y variance-covariance matrix
	inv_V = inverse(V)
	### Mixed model equations
	###### Fixed effects:
	b_hat = (X' * inv_V) * inverse(X * X' * inv_V) * y
	###### Random effects:
	u_hat = D * Z' * inv_V * (y - (X*b_hat))
	### Calculate the log-likelihood of the parameters (variances of the error and random effects: σ2e and σ2u) given the data
	### Or estimate the fixed and random effects as well as the fixed effects variance given the optimum parameters
	if OPTIMIZE
		################
		### OPTIMIZE ###
		################
		### Calculate the -log-likelihood given the error variance (σ2e) and random effects variance (σ2u) for maximum likelihood optimization to find the best fit σ2e and σ2u
		OUT = []
		if METHOD_VAR_EST == "ML"
			mu_y = (X * b_hat)
			n = length(y)
			push!(OUT, 0.5 * ( log(abs(det(V))) + ((y - mu_y)' * inv_V * (y - mu_y)) + (n*log(2*pi)) )) ### The negative log-likelihood function y given σ2e and σ2u
		elseif METHOD_VAR_EST == "REML"
			### NOTE: Z MUST BE SQUARE!
			M = try
					inverse(LinearAlgebra.cholesky(σ2u*Z + R)) ### Cholesky decomposition
				catch
					inverse(LinearAlgebra.lu(σ2u*Z + R).L) ### LU decomposition
				end
			y_new = M' * y
			intercept_new = sum(M', dims=2) #M' * ones(n)
			V_new = M' * V * M
			n = length(y_new)
			push!(OUT, 0.5 * ( log(abs(det(V_new))) .+ ((y_new - intercept_new)' * inverse(V_new) * (y_new - intercept_new)) .+ (n*log(2*pi)) )[1,1]) ### negative log-likelihood of σ2e and σ2u given X, y, and Z
		else
			println(string("Sorry. ", METHOD_VAR_EST, " is not a valid method of estimating the variances of the random effects effects. Please pick ML or REML."))
			exit()
		end
	else
		################
		### ESTIMATE ###
		################
		### Return the fixed effects (b_hat), random effects (u_hat), and optional fixed effects variance (σ2b_hat) estimates
		OUT = []
		if VARFIXEF
			### Computationally intensive!!!
			Vb_hat = inverse(X' * inverse(V) * X)
			σ2b_hat = diag(Vb_hat)
		else
			σ2b_hat = nothing
		end
		# push!(OUT, (fixef=b_hat, ranef=u_hat, var_fixef=σ2b_hat))
		push!(OUT, [b_hat, u_hat, σ2b_hat])
	end
	return(OUT[1])
end

#####################
### main function ###
#####################
"""
# ______________________
# Linear Mixed Modelling

`LMM(;X::Array{Float64,2}, y::Array{Float64,1}, Z::Array{Float64,2}, METHOD_VAR_EST::String=["ML", "REML"][1])`

The mixed linear model is defined as y = Xb + Zu + e, where:
	- no intercept is fitted
	- X is centered
	- y is centered which means the intercept is automatically the mean of y
	- Z is a square symmetric matrix of relatedness
	- u ~ N(0, σ2uI)
	- e ~ N(0, σ2eI)
	- y ~ (Z (σ2uI) Z') + (σ2eI)

# Inputs
1. *X* [Array{Float64,2}; xij ∈ [0.00, 1.00]]: matrix of allele frequencies (Pool-seq data) set as the fixed effects
2. *y* [Array{Float64,1}]: vector of phenotypes
3. *Z* [Array{Any,2}]: matrix of relatedness (Fst or XX'/p where p is the number of columns in X) set as the random effects
4. *METHOD_VAR_EST* [String]: method of variance parameter estimation choose "ML" for maximum likelihood or "REML" for restricted maximum likelihood (default="ML")

# Outputs
1. Intercept (b0) the mean of the phenotype data
2. Fixed effects estimates of X (b_hat)
3. Random effects estimates of Z (u_hat)

# Examples
```
using Distributions
using LinearAlgebra
n=5; l=200; Ve=1.00; Vu=1.00;
X = reshape(abs.(rand(n*l)), n, l)
Z = reshape(abs.(rand(n^2)), n, n)
b = zeros(l); b[[5,10,100]] = [1.00, 1.00, 1.00]
u = rand(Distributions.Normal(0.00, Vu), n)
e = rand(Distributions.Normal(0.00, Ve), n)
y = (X*b) .+ (Z*u) .+ e
b0, b_hat, u_hat = LMM_module.LMM(X=X, y=y, Z=Z)
b0, b_hat, u_hat = LMM_module.LMM(X=X, y=y, Z=Z, METHOD_VAR_EST="REML")
```
"""
function LMM(;X::Array{Float64,2}, y::Array{Float64,1}, Z::Array{Float64,2}, METHOD_VAR_EST::String=["ML", "REML"][1])
	### MODEL: y = Xb + Zu + e, where: u~N(0,Vu) and e~(0,Ve)
	### Generate copies of the input vector and matrices so that we don't modify the original ones
	X = copy(X); y = copy(y); Z = copy(Z);
	### Center y vector, X matrix and Z matrix (This way we don't need to estimate the intercept! It will be at zero!)
	uy = mean(y); uX = mean(X); uZ = mean(Z);
	sy = std(y); sX = std(X); sZ = std(Z);
	y[:] = (y .- uy)
	X[:,:] = (X .- uX)
	Z[:,:] = (Z .- uZ)
	### Maximum likelihood estimation of the error variance (Ve) and random covariate variance (Vu)
	# lower_limit = [0.0, 0.0]
	# upper_limit = [Inf, Inf]
	lower_limit = [1.0e-50, 1.0e-50]
	upper_limit = [1.0e+50, 1.0e+50]
	initial_VeVu = [1.0, 1.0]
	# optim_out = Optim.optimize(VeVu -> LMM_optim(VeVu, X=X, y=y, Z=Z, METHOD_VAR_EST=METHOD_VAR_EST, METHOD_FIX_EST=METHOD_FIX_EST, alfa=alfa, lambda=lambda, OPTIMIZE=true), lower_limit, upper_limit, initial_VeVu)
	### Optimization in Julia 1.0 is too darn slow! Going with R then! Importing Julia functions and variables into R.
	@rput inverse;
	@rput LMM_optim;
	@rput X;
	@rput y;
	@rput Z;
	@rput METHOD_VAR_EST;
	@rput lower_limit;
	@rput upper_limit;
	@rput initial_VeVu;
	R"optim_out = optim(par=initial_VeVu, fn=LMM_optim, X=X, y=y, Z=Z, METHOD_VAR_EST=METHOD_VAR_EST, OPTIMIZE=TRUE, method='L-BFGS-B', lower=lower_limit, upper=upper_limit)"
	@rget optim_out;
	Ve, Vu = optim_out[:par]
	### Fixed effects (b), random effects (u), and fixed effects variance (Vb) estimation (NOTE: No need to back-transform these estimated effects because we simply centered the explanatory variables meaning they're still on the same space)
	b, u, Vb = LMM_optim([Ve, Vu], X=X, y=y, Z=Z, METHOD_VAR_EST=METHOD_VAR_EST, OPTIMIZE=false)
	### Output
	return (b0=uy, b=b, u=u, Ve=Ve, Vu=Vu)
end

end #end of LMM module

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
include("../benchmark.jl")
include("create_report.jl")
overlapping_results = get_CUTEst_results()

#############
############# SUMMARY TABLE
#############

problem_list = collect(keys(first(overlapping_results)[2]))
method_list = collect(keys(overlapping_results))

different_status_problems = Array{String,1}()

begin
  for problem_name in problem_list

      previous_status = :none
      same_status = false

      for method_name in method_list
        info = overlapping_results[method_name][problem_name]

        if info.status == previous_status
          same_status = true
        end
        previous_status = info.status
      end

      if !same_status
        push!(different_status_problems,problem_name)
      end
  end
end

shared_fails = shared_failures(overlapping_results)
fail_only = Dict()
for (method_name, sum_data) in overlapping_results
    fail_only[method_name] = tot_failures(sum_data) - shared_fails
    println("total failures of $method_name is ", tot_failures(sum_data))
end

different_status_results = restrict_problems(overlapping_results,different_status_problems)

## create data frame
using DataFrames
df = DataFrame(
column_name = [
"failure",
"infeasible",
"unbounded",
#"worse obj",
#"better obj",
"KKT"
],
ipopt = [
fail_only["IPOPT"],
tot(different_status_results["IPOPT"],[:primal_infeasible]),
tot(different_status_results["IPOPT"],[:dual_infeasible]),
#winners["one-phase"],
#winners["IPOPT"],
tot(different_status_results["IPOPT"],[:Optimal])
],
one_phase = [
fail_only["one-phase"],
tot(different_status_results["one-phase"],[:primal_infeasible]),
tot(different_status_results["one-phase"],[:dual_infeasible]),
#winners["IPOPT"],
#winners["one-phase"],
tot(different_status_results["one-phase"],[:Optimal])
]
)

using CSV
CSV.write("$folder/summary.csv",df)

## bar plot of data frame
#Pkg.clone("https://github.com/JuliaPlots/StatPlots.jl.git")
using StatPlots
red2 = RGBA(1.0,0.0,0.0,0.6) #[red2,red2]
green2 = RGBA(1.0,65.0/256.0,256.0/256.0,63.0/256.0)
colour = [[:red,:red] [:orange,:orange] [:blue,:blue] [:green,:green]]
label = ["failure" "infeasible" "unbounded" "KKT"] #"worse obj" "better obj"
arr = Array(df)
data = convert(Array{Float64,2},arr[:,2:3])'
space = " "^50

for i = 1:length(label)
    groupedbar(data[:,1:i], bar_position = :stack,
    bar_width=0.9,
    colour=colour[:,1:i],
    label=label[:,1:i],
    xticks=[],
    xlabel="ipopt $space one phase",
    ylabel="number of problems",
    ylims=(0,80)
    )
    savefig("$folder/bar_$i.pdf")
end


text:
using SafeTestsets
import Base: isapprox

@time @safetestset "Basic Operators Interface" begin include("basic_operators_interface.jl") end
@time @safetestset "JacVec Operators Interface" begin include("jacvec_operators.jl") end
@time @safetestset "Composite Operators Interface" begin include("composite_operators_interface.jl") end

@time @safetestset "Derivative Operators Interface" begin include("derivative_operators_interface.jl") end
@time @safetestset "Dirichlet BCs" begin include("dirichlet.jl") end
@time @safetestset "Periodic BCs" begin include("periodic.jl") end
@time @safetestset "Neumann BCs" begin include("neumann.jl") end
@time @safetestset "2nd order check" begin include("2nd_order_check.jl") end
@time @safetestset "None BCs" begin include("none.jl") end
@time @safetestset "Finite Difference Operator" begin include("generic_operator_check.jl") end
#@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails
@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end
@time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end

using StaticArrays, DiffEqOperators
function isapprox(x::DerivativeOperator{T,S,LBC,RBC},y::FiniteDifference{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC}
    der_order           = (x,y) -> x.derivative_order == y.derivative_order
    approximation_order = (x,y) -> x.approximation_order == y.approximation_order
    dx                  = (x,y) -> !any(x.dx .!= y.dx)
    dimension           = (x,y) -> x.dimension == y.dimension
    stencil_length      = (x,y) -> x.stencil_length == y.stencil_length
    stencil_coefs       = (x,y) -> !any(!isapprox.(x.stencil_coefs,y.stencil_coefs; kwargs...))
    boundary_point_count= (x,y) -> x.boundary_point_count == y.boundary_point_count
    boundary_length     = (x,y) -> x.boundary_length == y.boundary_length
    low_boundary_coefs  = (x,y) -> !any(!isapprox.(x.low_boundary_coefs,y.low_boundary_coefs; kwargs...))
    high_boundary_coefs = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...))
    boundary_condition  = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...))
    t                   = (x,y) -> x.t ≈ y.t

    #in order of fast to slow comparison
    for f in (der_order,approximation_order,dimension,t,boundary_point_count,boundary_length,dx,boundary_condition,low_boundary_coefs,high_boundary_coefs,stencil_coefs)
        if !f(x,y)
            return false
        end
    end
    return true
end

isapprox(x::FiniteDifference{T,S,LBC,RBC},y::DerivativeOperator{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC} = isapprox(y,x; kwargs...)


text:
export createnode

################################################################
# Init
################################################################
# doc or element initialize
"""
    createnode(struct_nodetype)

Function to initialize the aml
"""
function createnode(::Type{AbsHtml}, struct_name::String = "html")
    out = HTMLDocument() # no URI and external id
    htmlNode = ElementNode(struct_name)
    setroot!(out, htmlNode) # adding html node
    return out
end

function createnode(::Type{AbsXml}, struct_name::String = "xml_root")
    out = XMLDocument() # version 1
    xmlNode = ElementNode(struct_name)
    setroot!(out, xmlNode) # adding html node
    return out
end

function createnode(::Type{<:AbsNormal}, struct_name::String)
    out = ElementNode(struct_name) # element node
    return out
end

function createnode(::Type{AbsText}, struct_name::String)
    out = TextNode(struct_name) # text node
    return out
end

# no type method
function createnode(struct_name::String)
    out = ElementNode(struct_name) # element node
    return out
end


text:
InputFolder = './Images/NucleiCluster/';
OutputFolder = './Results/Images/NucleiCluster/';
Rescale = 2;
Fill = -1;

@iA = '*.tif';

@fxg_sBlobDetPatch [iA] > [S];
params.GaussianGradRad = 1.5;
params.GaussianVoteRad = 6;
params.GradMagThr = 0.2;
params.Rmin = 7;
params.Rmax = 11;
params.NAngles = 24;
params.Delta = [pi/8 pi/16];
params.DScale = [1 1];
params.Thr = 0.5;
/endf

/show iA > S;
/keep S > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
let xmin = (-1.0, -1.0, -1.0)
    xmax = ( 1.0,  1.0,  1.0)

    for n = 0:2:10
        @test nintegrate((x,y,z)->x^n, xmin, xmax)[1] ≈ 8/(n+1)
        @test nintegrate((x,y,z)->y^n, xmin, xmax)[1] ≈ 8/(n+1)
        @test nintegrate((x,y,z)->z^n, xmin, xmax)[1] ≈ 8/(n+1)
    end

    @test nintegrate((x,y,z)->x^2*y^2, xmin, xmax)[1] ≈ 8/9
    @test nintegrate((x,y,z)->z^2*y^2, xmin, xmax)[1] ≈ 8/9
    @test nintegrate((x,y,z)->x^2*z^2, xmin, xmax)[1] ≈ 8/9

    @test nintegrate((x,y,z)->x^2*y^2*z^2, xmin, xmax)[1] ≈ 8/27

    ɛ = eps()

    for n = 1:2:11
        @test nintegrate((x,y,z)->x^n, xmin, xmax)[1] < ɛ
    end

    @test nintegrate((x,y,z)->x^2*z^3, xmin, xmax)[1] < ɛ
    @test nintegrate((x,y,z)->x^3*z^2, xmin, xmax)[1] < ɛ

    @test nintegrate((x,y,z)->x^2*y*z, xmin, xmax)[1] < ɛ
    @test nintegrate((x,y,z)->x*y^2*z, xmin, xmax)[1] < ɛ
    @test nintegrate((x,y,z)->x*y*z^2, xmin, xmax)[1] < ɛ

end

let xmin = (0.0, 0.0, 0.0)
    xmax = (1.0, 1.0, 1.0)

    for n = 0:11
        @test nintegrate((x,y,z)->x^n, xmin, xmax)[1] ≈ 1/(n+1)
        @test nintegrate((x,y,z)->y^n, xmin, xmax)[1] ≈ 1/(n+1)
        @test nintegrate((x,y,z)->z^n, xmin, xmax)[1] ≈ 1/(n+1)
    end

    @test nintegrate((x,y,z)->x^2*y^2, xmin, xmax)[1] ≈ 1/9
    @test nintegrate((x,y,z)->z^2*y^2, xmin, xmax)[1] ≈ 1/9
    @test nintegrate((x,y,z)->x^2*z^2, xmin, xmax)[1] ≈ 1/9

    @test nintegrate((x,y,z)->x^2*z^3, xmin, xmax)[1] ≈ 1/12
    @test nintegrate((x,y,z)->x^3*z^2, xmin, xmax)[1] ≈ 1/12
    @test nintegrate((x,y,z)->x^2*y*z, xmin, xmax)[1] ≈ 1/12
    @test nintegrate((x,y,z)->x*y^2*z, xmin, xmax)[1] ≈ 1/12
    @test nintegrate((x,y,z)->x*y*z^2, xmin, xmax)[1] ≈ 1/12

    @test nintegrate((x,y,z)->x^3*y*z, xmin, xmax)[1] ≈ 1/16
    @test nintegrate((x,y,z)->x*y^3*z, xmin, xmax)[1] ≈ 1/16
    @test nintegrate((x,y,z)->x*y*z^3, xmin, xmax)[1] ≈ 1/16

    @test nintegrate((x,y,z)->x^2*y^2*z^2, xmin, xmax)[1] ≈ 1/27

end

let f = (x, y, z) -> x * sin(2y) * cos(3z)

    xmin = (0.0, 0.0, 0.0)
    xmax = (1.0, 1.0, 1.0)

    (I, E, n, R) = nintegrate(f, xmin, xmax)
    @test I == nintegrate(f, R)
    @test length(R) == 1

    W = weightedpoints(R)
    @test I ≈ nintegrate(f, W)

    W = weightedpoints(NIntegration.idem, R)
    @test I ≈ nintegrate(f, W)

    p = [0.5, 0.5, 0.5]
    io = IOBuffer()
    show(io, R)
    @test String(take!(io)) == "1 subregion"
    show(io, R[1])
    @test String(take!(io)) == string("(", p, ", ", p, ")")

end


text:
#
# cameras.jl --
#
# Management of cameras for the Julia interface to the C libraries of TAO, a
# Toolkit for Adaptive Optics.
#
#------------------------------------------------------------------------------
#
# This file is part of TAO software (https://git-cral.univ-lyon1.fr/tao)
# licensed under the MIT license.
#
# Copyright (C) 2018-2021, Éric Thiébaut.
#

propertynames(cam::SharedCamera) =
    (
      #tao attribute
      :owner,
      :shmid,
      :size,
      :lock,

      :listlength,
      :last,
      :next,
      :lastTS,
      :nextTS,

      :cameras,
      :attachedCam,
      :img_config,

     )
# dispatch for sharedcamera
getproperty(dev::SharedCamera, sym::Symbol) =
                            getproperty(dev,Val(sym))

#dispatch tao attribute
getproperty(dev::SharedCamera,::Val{:shmid}) = getattribute(dev,Val(:shmid))
getproperty(dev::SharedCamera,::Val{:size}) = getattribute(dev,Val(:size))
getproperty(dev::SharedCamera,::Val{:owner}) = getattribute(dev,Val(:owner))

# relay to taolib
getattribute(obj::SharedCamera, ::Val{:shmid}) =
    ccall((:tao_get_shared_data_shmid, taolib), ShmId,
          (Ptr{AbstractSharedObject},), obj)

getattribute(obj::SharedCamera, ::Val{:size}) =
    ccall((:tao_get_shared_data_size, taolib), Csize_t,
          (Ptr{AbstractSharedObject},), obj)

getattribute(obj::SharedCamera, ::Val{:owner}) =
  _pointer_to_string(ccall((:tao_get_shared_object_owner, taolib),
                             Ptr{UInt8}, (Ptr{AbstractSharedObject},), obj))

# dispatch Julia properties
getproperty(dev::SharedCamera, ::Val{sym}) where{sym} =
                                        getfield(dev, sym)

# property is read-only
setproperty!(dev::SharedCamera, sym::Symbol, val)  =
                                setattribute!(dev,Val(sym),val)

# attribute manipulation function
inc_attachedCam(dev::SharedCamera) = begin
                        val = dev.attachedCam +1
                        setfield!(dev, :attachedCam,Int8(val))
                    end

for sym in (
            :last,
            :next,
            :lastTs,
            :nextTS
                    )
    _sym = "$sym"
    @eval $(Symbol("set_",sym))(dev::SharedCamera,val::Integer) =
                                    setfield!(dev,Symbol($_sym),Int16(val))
end

set_img_config(shcam::SharedCamera, conf::ImageConfigContext) = setfield!(shcam,:img_config,conf)
#--- RemoteCamera property

propertynames(cam::RemoteCamera) =
    (
      # RemoteCamera properties
      :arrays,
      :timestamps,
      :img,
      :imgTime,
      :cmds,
      :time_origin,
      :no_cmds,

      # SharedCamera
      :owner,
      :shmid,
      :lock,
      :attachedcam,
      :listlength,
      :state,
      :counter,
      :last,
      :next,
      :lastTS,
      :nextTS

     )

#top level dispatch
getproperty(remcam::RemoteCamera,sym::Symbol) = getproperty(remcam, Val(sym))

#RemoteCamera properties (read-only)
for sym in (:timestamps,
            :img,
            :imgTime,
            :cmds,
            :shmids,
            :arrays,
            :device,
            :no_cmds,
            :time_origin,
            :worker_pid)
    _sym = "$sym"
    @eval getproperty(remoteCam::RemoteCamera,::$(Val{sym})) =
                        getfield(remoteCam,Symbol($_sym))
end

# Shared Camera properties
# getattribute(cam::RemoteCamera, ::Val{:lock}) = getfield(device(cam), :lock)
getproperty(remoteCam::RemoteCamera, key::Val)=
                                        getproperty(device(remoteCam), key)

# wall to setting properties
setproperty!(remoteCam::RemoteCamera, sym::Symbol, val) =
    throw_non_existing_or_read_only_attribute(remoteCam,sym)

# Constructors for `RemoteCamera`s.
# RemoteCamera(dev::SharedCamera) = RemoteCamera{dev.pixeltype}(dev)

# RemoteCamera(srv::ServerName) = RemoteCamera(attach(SharedCamera, srv))
# RemoteCamera{T}(srv::ServerName) where {T<:AbstractFloat} =
#     RemoteCamera{T}(attach(SharedCamera, srv))

#--- Accessors.
camera(cam::RemoteCamera) = cam
camera(cam::SharedCamera) = cam
device(cam::RemoteCamera) = getfield(cam, :device)
device(cam::SharedCamera,i::Int64) = cam.cameras[i]
eltype(::AbstractCamera{T}) where {T} = T

show(io::IO, cam::RemoteCamera{T}) where {T} =
    print(io, "RemoteCamera{$T}(owner=\"", cam.owner,")")

#--- create functions
function create(::Type{SharedCamera}; owner::AbstractString = default_owner(),
                perms::Integer = 0o600)

    length(owner) < SHARED_OWNER_SIZE || error("owner name too long")

    ptr = ccall((:tao_create_shared_object, taolib), Ptr{AbstractSharedObject},
                (Cstring, UInt32, Csize_t, Cuint),
                owner, _fix_shared_object_type(SHARED_CAMERA), 464, perms)
                # 464 bytes from sizeof(tao_shar_camera) in C program
    _check(ptr != C_NULL)

    return _wrap(SharedCamera, ptr)
end

## Camera util functions
function register(dev::SharedCamera, cam::Camera)
    inc_attachedCam(dev)
    try
        dev.cameras[dev.attachedCam] = cam
    catch
        throw(ErrorException("No more space to attach cameras"))
    end
end

function broadcast_shmids(shmids::Vector{ShmId})
    fname = "shmids.txt"
    path = "/tmp/FLICameras/"

    fname ∈ readdir(path) || touch(joinpath(path,fname))

    open(path*fname,"w") do f
        for shmid in shmids
            write(f,@sprintf("%d\n",shmid))
        end
    end
end

const img_param = fieldnames(ImageConfigContext)
const img_param_type = fieldtypes(ImageConfigContext)

#TODO: Fix this
function _read_and_update_config(shcam::SharedCamera)
    fname = "img_config.txt"
    path = "/tmp/FLICameras/"
    fname ∈ readdir(path) || throw(LoadError("image config doest not exist"))
    f = open(path*fname,"r")
    rd = readlines(f)

    new_conf = shcam.img_config

    for i in 1:length(rd)
      if !isempty(rd[i])
        setfield!(new_conf,img_param[i],parse(img_param_type[i],rd[i]))
      end
    end
    @show new_conf
    set_img_config(shcam,new_conf)
  end

#--- Command and state mapper

for (cmd, now_state, next_state) in (
        (:CMD_INIT,  :STATE_UNKNOWN   ,   STATE_INIT),
        (:CMD_WORK,  :STATE_WAIT      ,   STATE_WORK),
        (:CMD_STOP,  :STATE_WORK      ,   STATE_WAIT),
        (:CMD_QUIT,  :STATE_WAIT      ,   STATE_QUIT),
    )

    @eval sort_next_state(::Val{$cmd}, ::Val{$now_state}) = $next_state
end

# sort_next_state(Val(),::Val{SIG_ERROR})
for (sig, cmd, next_state) in (
        (:SIG_DONE,  :CMD_INIT      ,   STATE_WAIT),
        (:SIG_DONE,  :CMD_STOP      ,   STATE_WAIT),
        (:SIG_DONE,  :CMD_WORK      ,   STATE_WORK),
        (:SIG_DONE,  :CMD_ABORT     ,   STATE_WAIT),
        (:SIG_DONE,  :CMD_QUIT      ,   STATE_QUIT),
        (:SIG_ERROR,  :CMD_WORK      ,   STATE_ERROR),
    )

    @eval sort_next_state(::Val{$cmd}, ::Val{$sig}) = $next_state
end

sort_next_state(cmd::RemoteCameraCommand, current_state::RemoteCameraState) = sort_next_state(Val(cmd), Val(current_state))
sort_next_state(cmd::RemoteCameraCommand, shcam_sig::ShCamSIG) = sort_next_state(Val(cmd), Val(shcam_sig))

thread_safe_wait(r::Condition) = begin
                                lock(r)
                                try
                                  wait(r)
                                finally
                                  unlock(r)
                                end
                              end
thread_safe_notify(r::Condition) = begin
                                lock(r)
                                try
                                  wait(r)
                                finally
                                  unlock(r)
                                end
                              end
## Listening function
@inline no_new_cmds(cmds::SharedArray{Cint,1}) = begin
                        sum(cmds .== -1) == length(cmds)
                      end
pop!(cmds::SharedArray{Cint,1}) = begin
                          val = cmds[1]
                          cmds[:] = vcat(cmds[2:end],-1)
                          return val
                        end
"""
  Listening
  1. check cmds written to the shmid
  2. read the cmds
  3. start the command
"""
listening(shcam::SharedCamera, remcam::RemoteCamera) = @async _listening(shcam, remcam)

function _listening(shcam::SharedCamera, remcam::RemoteCamera)

  while true
   # check the command
   @info "wait cmds"

   @async _keep_checking_cmds(remcam.cmds, remcam.no_cmds )
   wait(remcam.no_cmds)
   #  read new cmd
   cmd = rdlock(remcam.cmds,0.5) do
       pop!(remcam.cmds)
   end
   # sent the command to the camera
   if next_camera_operation(RemoteCameraCommand(cmd),shcam, remcam)
     @info "Command successful..."
   else
     @info "Command failed..."
   end

 end
end

"""
   This function keeps checking the cmd array if there is a new command written to
   the shared array
"""
function _keep_checking_cmds(cmd_array::SharedArray{Cint,1}, no_cmds::Condition)

     cmd0 = rdlock(cmd_array,0.5) do
         cmd_array[1]
     end
     cmd_now = copy(cmd0)

   while cmd_now == cmd0
     cmd_now = rdlock(cmd_array,0.5) do
         cmd_array[1]
     end
     sleep(0.01)

   end
   notify(no_cmds)
   nothing
end

## attach extension
function attach(::Type{SharedCamera}, shmid::Integer)
    # Attach the shared object to the address space of the caller, then wrap it
    # in a Julia object.
    ptr = ccall((:tao_attach_shared_camera, taolib),
                Ptr{AbstractSharedObject}, (ShmId,), shmid)
    _check(ptr != C_NULL)
    cam = _wrap(SharedCamera, ptr)
    return cam
end

# attached image and image timestamp shared array on a worker process
function attach_remote_process()
    # read shared array shmids
    fname = "/tmp/FLICameras/shmids.txt"
    f = open(fname,"r")
    rd = readlines(f)
    shmids = []
    for line in rd
      push!(shmids,parse(Int64,line))
    end

    img_array = attach(SharedArray{UInt8},shmids[1])
    imgTime_array = attach(SharedArray{Float64},shmids[2])

    return img_array, imgTime_array
end
#--- shared camera operations
update_worker_pid(pid::Integer,remcam::RemoteCamera) = begin
  wrlock(remcam.device,0.1) do
    setfield!(remcam,:worker_pid,pid)
  end
end

# Operation - Commands table
next_camera_operation(cmd::RemoteCameraCommand,shcam::SharedCamera, remcam::RemoteCamera) = next_camera_operation(Val(cmd),shcam ,remcam)

for (cmd, cam_op) in (
    (:CMD_WORK,     :start),
    (:CMD_STOP,     :stop),
    (:CMD_CONFIG,   :config),
    (:CMD_UPDATE,   :update),
    (:CMD_RESET,    :reset)

    )
@eval begin
    function next_camera_operation(::Val{$cmd}, shcam::SharedCamera, remcam::RemoteCamera)
        try
            cameraTask = $cam_op(shcam,remcam)
            str_op = $cam_op
            @info "Executed  $str_op \n"
        catch ex
            @assert ex.code != 0
            @warn "Error occurs at $(ex.func)"

            return_val =  false
            return return_val
        end

        return_val =  true

        return return_val
    end
end
end

"""
    FLICameras.start(cam)
    Start camera acquisition loop on a worker process.
    Sending data to shared arrays
""" start

start(shcam::SharedCamera,remcam::RemoteCamera ) = begin

    shcam.attachedCam > 0 || throw(ErrorException("No attached cameras"))

    camera = device(shcam,1)
    pid =[0]
    try
      # start working
      pid[1] = working(1)
      @info "worker pid = $(pid[1])"

    catch ex
      rethrow(ex)
    end

    # broadcast worker pid
    update_worker_pid(pid[1],remcam)
    @info "start is DONE"

    nothing
end

"""
  FLICameras.update(SharedCamera, RemoteCamera)
  This function interrupts the ongoing acquisition loop to update the camera
  configuration and start a new acquisition loop
  The new camera configuration setting is read through a file.
""" update
function update(shcam::SharedCamera,remcam::RemoteCamera)
  stop(shcam,remcam)
  _read_and_update_config(shcam)
  config(shcam,remcam)
  @info "Camera configuration has been updated"
  working(1)
  nothing
end

"""
    FLICameras.stop(cam; timeout=5.0)

    stops image acquisition by shared or remote camera `cam` not waiting more than
    the limit set by `timeout` in seconds.  Nothing is done if acquisition is not
    running or about to start.
stop
"""
stop(shcam::SharedCamera,remcam::RemoteCamera) = begin

    shcam.attachedCam > 0 || throw(ErrorException("No attached cameras"))

    camera = device(shcam,1)
    pid = remcam.worker_pid
    try
      interrupt(pid)

    catch ex
      rethrow(ex)
    end

    nothing

end

"""
  FLICameras.config(SharedCamera, RemtoeCamera)
  Change camera setting according to the ImageConfigContext in SharedCamera
""" config

config(shcam::SharedCamera,remcam::RemoteCamera) = begin

    shcam.attachedCam > 0 || throw(ErrorException("No attached cameras"))

    camera = device(shcam,1)
    ROI = []    #contain ROI info
    try

      for k in 1:length(fieldnames(ImageConfigContext))
        if k <=4
          val =  getfield(shcam.img_config,param)
          push!(ROI,val)
        elseif k == 5
          set_ROI(camera,ROI)
          val =  getfield(shcam.img_config,:exposuretime)
          set_exposuretime(camera,val)
        else
          val =  getfield(shcam.img_config,:bytePerPixel)
        end

    end

    catch ex
        rethrow(ex)
    end

    nothing
end


text:
@testset "read_opf: simple_plant.opf -> read" begin
    mtg = read_opf("files/simple_plant.opf", Dict)
end

mtg = read_opf("files/simple_plant.opf", Dict)

@testset "read_opf: simple_plant.opf -> attributes" begin
    @test length(mtg) == 7
    @test haskey(mtg.attributes, :ref_meshes)
    @test names(mtg) == [:ref_meshes, :FileName, :geometry, :Length, :Width, :XEuler]
    @test descendants(mtg, :Length) == Any[nothing, nothing, 0.1f0, 0.2f0, 0.1f0, 0.2f0]
    @test descendants(mtg, :Width) == Any[nothing, nothing, 0.02f0, 0.1f0, 0.02f0, 0.1f0]
    @test descendants(mtg, :FileName) == Any["ArchiTree", nothing, nothing, nothing, nothing, nothing]
    @test descendants(mtg, :XEuler) == Any[nothing, nothing, nothing, nothing, 180.0f0, nothing]
end

@testset "read_opf: simple_plant.opf -> ref. meshes" begin
    ref_meshes = mtg[:ref_meshes].meshes
    @test length(ref_meshes) == 2
    first_mesh = ref_meshes[1]
    @test isa(first_mesh.mesh, Meshes.SimpleMesh)
    @test length(first_mesh.mesh) == 50
    @test isa(first_mesh.material, Phong)
    @test first_mesh.name == "Mesh0"
    @test isa(first_mesh.normals, SVector{50,Meshes.Point3})
    @test first_mesh.taper == true
    @test isa(first_mesh.texture_coords, SVector{50,Meshes.Point2})
end

@testset "read_opf: simple_plant.opf -> meshes" begin
    Internode = get_node(mtg, 4)

    @test collect(keys(Internode.attributes)) == [:geometry, :Length, :Width]
    @test isa(Internode[:geometry], PlantGeom.geometry)

    geom = Internode[:geometry]
    @test isnothing(geom.mesh)
    @test geom.ref_mesh === mtg[:ref_meshes].meshes[geom.ref_mesh_index]
    @test isa(geom.transformation, AffineMap{Matrix{Float64},Vector{Float64}})
    @test geom.transformation([1, 1, 1]) ≈ [-1.0, 1.0, 10.0] atol = 1.0e-6
    # NB: last one is 10 because there is some tappering
end

@testset "read_opf: read coffee.opf" begin
    mtg = read_opf("files/coffee.opf", Dict)
    @test length(mtg) == 4191
    @test get_attributes(mtg) ==
          [
        :FileName, :Variety, :ref_meshes, :Plot, :Treatment, :File,
        :Name, :geometry, :XEuler, :Length, :Width, :Stifness, :YInsertionAngle,
        :Plagiotropy, :Phyllotaxy, :StiffnessAngle, :Area, :XInsertionAngle
    ]

    @test sum(descendants(mtg, :Area, ignore_nothing = true)) ≈ 77961.414f0 atol = 1e-6
end


text:
"""
Drawings and calculations useful for research in Function Iteration Theory,
or Discrete Dynamical Systems.
"""
module SDD

using Reexport

@reexport using SDDCore, SDDGraphics
#@reexport import

include("graphicalanalysis.jl")
include("orbits.jl")
include("trappedpoints.jl")
include("mandelbrots.jl")
include("preimages.jl")
include("basins.jl")

export
    plot,
    graphicalanalysis,
    drawpointorbitR2,
    drawpointorbitC,
    drawpointssetorbitR2,
    drawpointssetorbitC,
    drawpointorbitpathR,
    drawpointorbitpathR2,
    drawpointorbitpathC,
    drawtrappedpointsR2,
    drawtrappedpointsC,
    drawmandelbrotR2,
    drawmandelbrotC,
    drawpreimageR2,
    drawpreimageC,
    drawbasinsR2,
    drawbasinsC

end # module


text:
using AdaptiveWindows
using Test

@testset verbose=true "Adaptive Mean" begin

    @testset "Mean Computation " begin
        m = AdaptiveMean(δ = 1e-9)

        r = randn(1000)

        fit!(m, r)

        m1 = sum(r) / length(r)
        m2 = value(m)

        @test m1 ≈ m2
        ad = AdaptiveMean()

        # This should not trigger a truncated window
        fit!(ad, randn(10_000))
        @test stats(ad).n == 10_000

        # Changing the distribution should trigger a truncated window
        fit!(ad, 1 .+ randn(10_000))
        @test 9_900 < stats(ad).n < 20_000

        # check truncation of shifting using the callback function
        shifted = false

        m = AdaptiveMean(onshiftdetected = ad -> shifted = true)

        for i in 1:1_000
            r = randn()
            if i > 500
                r += 1
            end
            fit!(m, r)
        end

        @test shifted
    end

    function consistent(ad)
        total = sum(nobs(v) for v in ad.window)
        total == nobs(ad.stats)
    end

    @testset "Memory Management" begin

        m = AdaptiveMean()
        fit!(m, 1)
        @test nobs(m.window[1]) == 0
        @test nobs(m.window[2]) == 1
        @test nobs(m.window[3]) == 0
        fit!(m, 1)
        @test nobs(m.window[1]) == 0
        @test nobs(m.window[2]) == 1
        @test nobs(m.window[3]) == 1
        fit!(m, 1)
        fit!(m, 1)
        fit!(m, 1)
        fit!(m, 1)
        fit!(m, 1)
        @test consistent(m)

        m = AdaptiveMean()
        n = AdaptiveWindows.M * ( 1 + 2 + 4)
        fit!(m, ones(n))

        @test length(m.window) <= AdaptiveWindows.M * log2(n)
        @test nobs(m) == n
        @test consistent(m)

        mn = AdaptiveMean()
        n = 1<<12

        # withoutdropping for speed
        fit!(withoutdropping(mn), ones(n))
        m = AdaptiveWindows.M
        expected = m * ceil(log2(n) - log2(m))
        @test length(mn.window) == expected
        @test nobs(mn) == n
        @test consistent(mn)

        # Maximum amount of memory
        mn = withmaxlength(AdaptiveMean(), 3)
        fit!(mn, rand(10000))
        @test length(mn.ad.window) == AdaptiveWindows.M * 3
        @test consistent(mn.ad)

    end
end


text:
using Base: Callable, TwicePrecision
using GeometricBase
using Test

@testset "$(rpad("Null Types",80))" begin

    invs = NullInvariants()
    @test_throws ErrorException invs[1]
    @test_throws ErrorException for i in invs end

    params = NullParameters()
    @test_throws ErrorException params[1]
    @test_throws ErrorException for p in params end

    period = NullPeriodicity()
    @test_throws ErrorException period[1]
    @test_throws ErrorException for p in period end

end

@testset "$(rpad("Optional Types",80))" begin

    @test Nothing <: OptionalArray
    @test Nothing <: OptionalArray{AbstractArray}
    @test Nothing <: OptionalArray{Matrix}

    @test Matrix <: OptionalArray
    @test Matrix <: OptionalArray{AbstractArray}
    @test Matrix <: OptionalArray{Matrix}

    @test Nothing         <: OptionalAbstractArray
    @test Matrix          <: OptionalAbstractArray
    @test Matrix{Int}     <: OptionalAbstractArray
    @test AbstractMatrix  <: OptionalAbstractArray

    @test Nothing         <: OptionalCallable
    @test Callable        <: OptionalCallable

    @test Nothing         <: OptionalFunction
    @test Function        <: OptionalFunction

    @test Nothing         <: OptionalNamedTuple
    @test NamedTuple      <: OptionalNamedTuple

    @test Nothing         <: OptionalTuple
    @test Tuple           <: OptionalTuple

    @test NullInvariants  <: OptionalInvariants
    @test NamedTuple      <: OptionalInvariants

    @test NullParameters  <: OptionalParameters
    @test NamedTuple      <: OptionalParameters

    @test NullPeriodicity <: OptionalPeriodicity
    @test AbstractArray   <: OptionalPeriodicity
    @test Array           <: OptionalPeriodicity
    @test Matrix          <: OptionalPeriodicity
    @test Matrix{Int}     <: OptionalPeriodicity

end

@testset "$(rpad("Data Types",80))" begin

    @test Number <: AbstractData
    @test Int    <: AbstractData
    @test Matrix{Number} <: AbstractData
    @test Matrix{Int}    <: AbstractData
    @test AbstractMatrix{Number} <: AbstractData
    @test AbstractMatrix{Int}    <: AbstractData

    @test Matrix{Number} <: State
    @test Matrix{Int}    <: State
    @test Matrix{Number} <: State{Number}
    @test Matrix{Int}    <: State{<:Number}
    @test Matrix{Int}    <: State{Int}
    @test AbstractMatrix{Number} <: State
    @test AbstractMatrix{Int}    <: State

    @test Vector{Matrix{Int}} <: StateVector
    @test Vector{Matrix{Number}} <: StateVector
    @test Vector{AbstractMatrix{Int}} <: StateVector
    @test Vector{AbstractMatrix{Number}} <: StateVector

    @test Vector{Int} <: SolutionVector
    @test Vector{Int} <: SolutionVector{Int}
    @test Vector{TwicePrecision{Int}} <: SolutionVector
    @test Vector{TwicePrecision{Int}} <: SolutionVector{TwicePrecision{Int}}

end


text:
@testset "constructor test" begin

    @test name(test_model) ≡ NAME

    @test final(test_model) ≡ false

    @test isdir(directory(test_model))

    @test isa(objective(test_model), Missing)

    @test isa(gradient(test_model), Missing)

    @test isa(initial_iterate(test_model), Missing)

    @test isa(dimension(test_model), Missing)

    @test isa(formula(test_model), Missing)
end

@testset "loading model test" begin

    @test objective!(test_model, f) ≡ f

    @test final(test_model) ≡ false

    @test gradient!(test_model, ∇f!) ≡ ∇f!

    @test final(test_model) ≡ true

    @test !isa(objective!(test_model, missing), Missing)

    @test !isa(gradient!(test_model, missing), Missing)

    @test initial_iterate!(test_model, x₀) ≡ x₀

    @test formula!(test_model, LaTex) ≡ LaTex
end

@testset "loaded model test" begin

    @test objective(test_model) ≡ f

    @test gradient(test_model) ≡ ∇f!

    @test initial_iterate(test_model) ≡ x₀

    @test dimension(test_model) ≡ n

    @test formula(test_model) ≡ LaTex

    @test obj(test_model, x₀) ≈ f₀

    @test grad!(test_model, similar(x₀), x₀) ≈ ∇f₀

    @test grad(test_model, x₀) ≈ ∇f₀

    @test hessAD(test_model, x₀) ≈ ∇²f

    @test hess_sample(test_model, x₀, ∇f₀) ≈ h₀
end


text:
include("../utils.jl")

@testset "Rotations" begin
    M = Manifolds.Rotations(2)
    @test repr(M) == "Rotations(2)"
    @test representation_size(M) == (2, 2)
    @test injectivity_radius(M) == π * sqrt(2.0)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0]) == π * sqrt(2.0)
    @test injectivity_radius(M, ExponentialRetraction()) == π * sqrt(2.0)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0], ExponentialRetraction()) ==
          π * sqrt(2.0)
    @test injectivity_radius(M, PolarRetraction()) ≈ π / sqrt(2)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0], PolarRetraction()) ≈ π / sqrt(2)
    @test get_embedding(M) == Euclidean(2, 2)
    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{2,2,Float64,4})
    retraction_methods = [Manifolds.PolarRetraction(), Manifolds.QRRetraction()]

    inverse_retraction_methods =
        [Manifolds.PolarInverseRetraction(), Manifolds.QRInverseRetraction()]

    basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))

    @testset "vee/hat" begin
        M = Manifolds.Rotations(2)
        Xf = [1.23]
        p = Matrix{Float64}(I, 2, 2)
        X = Manifolds.hat(M, p, Xf)
        @test isa(X, AbstractMatrix)
        @test norm(M, p, X) / sqrt(2) ≈ norm(Xf)
        @test Manifolds.vee(M, p, X) == Xf

        X = project(M, p, randn(2, 2))
        Xf = Manifolds.vee(M, p, X)
        @test isa(Xf, AbstractVector)
        @test Manifolds.hat(M, p, Xf) == X
    end

    for T in types
        angles = (0.0, π / 2, 2π / 3, π / 4)
        pts = [convert(T, [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]) for ϕ in angles]
        test_manifold(
            M,
            pts;
            test_injectivity_radius=false,
            test_project_tangent=true,
            test_musical_isomorphisms=true,
            test_default_vector_transport=true,
            retraction_methods=retraction_methods,
            inverse_retraction_methods=inverse_retraction_methods,
            point_distributions=[Manifolds.normal_rotation_distribution(M, pts[1], 1.0)],
            tvector_distributions=[Manifolds.normal_tvector_distribution(M, pts[1], 1.0)],
            basis_types_to_from=basis_types,
            test_inplace=true,
            test_rand_point=true,
        )

        @testset "log edge cases" begin
            X = Manifolds.hat(M, pts[1], [Float64(π)])
            p = exp(M, pts[1], X)
            @test isapprox(p, exp(M, pts[1], log(M, pts[1], p)))
        end

        X = log(M, pts[1], pts[2])
        @test norm(M, pts[1], X) ≈ (angles[2] - angles[1]) * sqrt(2)

        # check that exp! does not have a side effect
        q = allocate(pts[1])
        copyto!(M, q, pts[1])
        q2 = exp(M, pts[1], X)
        exp!(M, q, q, X)
        @test norm(q - q2) ≈ 0

        X14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction())
        p4_polar = retract(M, pts[1], X14_polar, Manifolds.PolarRetraction())
        @test isapprox(M, pts[4], p4_polar)

        X14_qr = inverse_retract(M, pts[1], pts[4], Manifolds.QRInverseRetraction())
        p4_qr = retract(M, pts[1], X14_qr, Manifolds.QRRetraction())
        @test isapprox(M, pts[4], p4_qr)
    end

    @testset "Distribution tests" begin
        usd_mmatrix =
            Manifolds.normal_rotation_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
        @test isa(rand(usd_mmatrix), MMatrix)

        usd1_mmatrix =
            Manifolds.normal_rotation_distribution(Rotations(1), (@MMatrix [1.0]), 1.0)
        @test isa(rand(usd1_mmatrix), MMatrix)
        @test rand(usd1_mmatrix) == @MMatrix [1.0]

        gtsd_mvector =
            Manifolds.normal_tvector_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
        @test isa(rand(gtsd_mvector), MMatrix)
    end

    Random.seed!(42)
    for n in (3, 4, 5)
        @testset "Rotations: SO($n)" begin
            SOn = Manifolds.Rotations(n)
            ptd = Manifolds.normal_rotation_distribution(SOn, Matrix(1.0I, n, n), 1.0)
            tvd = Manifolds.normal_tvector_distribution(SOn, Matrix(1.0I, n, n), 1.0)
            pts = [rand(ptd) for _ in 1:3]
            diag_basis_1 = if n == 3
                DiagonalizingOrthonormalBasis(
                    [
                        0.0 0.24800271831269094 0.30019597622794186
                        -0.24800271831269094 0.0 -0.5902347224334308
                        -0.30019597622794186 0.5902347224334308 0.0
                    ],
                )
            else
                DiagonalizingOrthonormalBasis(rand(SOn; vector_at=pts[1]))
            end
            diag_basis_2 = DiagonalizingOrthonormalBasis(
                hat(SOn, pts[1], [1.0, zeros(manifold_dimension(SOn) - 1)...]),
            )
            test_manifold(
                SOn,
                pts;
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_mutating_rand=true,
                test_default_vector_transport=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[ptd],
                tvector_distributions=[tvd],
                basis_types_to_from=basis_types,
                basis_types_vecs=(diag_basis_1, diag_basis_2),
                exp_log_atol_multiplier=250,
                retraction_atol_multiplier=12,
                test_inplace=true,
                test_rand_point=true,
                test_rand_tvector=true,
            )

            @testset "vee/hat" begin
                p = Matrix(1.0I, n, n)
                Xf = randn(manifold_dimension(SOn))
                X = Manifolds.hat(SOn, p, Xf)
                @test isa(X, AbstractMatrix)
                @test norm(SOn, p, X) / sqrt(2) ≈ norm(Xf)
                @test Manifolds.vee(SOn, p, X) == Xf

                X = project(SOn, p, randn(n, n))
                Xf = Manifolds.vee(SOn, p, X)
                @test isa(Xf, AbstractVector)
                @test Manifolds.hat(SOn, p, Xf) ≈ X
            end

            if n == 4
                @testset "exp/log edge cases" begin
                    Xs = [
                        [0, 0, π, 0, 0, π],  # θ = (π, π)
                        [0, 0, π, 0, 0, 0],  # θ = (π, 0)
                        [0, 0, π / 2, 0, 0, π],  # θ = (π, π/2)
                        [0, 0, π, 0, 0, 0] ./ 2,  # θ = (π/2, 0)
                        [0, 0, π, 0, 0, π] ./ 2,  # θ = (π/2, π/2)
                        [0, 0, 0, 0, 0, 0],  # θ = (0, 0)
                        [0, 0, 1, 0, 0, 1] .* 1e-100, # α = β ≈ 0
                        [0, 0, 1, 0, 0, 1] .* 1e-6, # α = β ⩰ 0
                        [0, 0, 10, 0, 0, 1] .* 1e-6, # α ⪆ β ⩰ 0
                        [0, 0, π / 4, 0, 0, π / 4 - 1e-6], # α ⪆ β > 0
                    ]
                    for Xf in Xs
                        @testset "rotation vector $Xf" begin
                            X = Manifolds.hat(SOn, Matrix(1.0I, n, n), Xf)
                            p = exp(X)
                            @test p ≈ exp(SOn, one(p), X)
                            p2 = exp(log(SOn, one(p), p))
                            @test isapprox(p, p2; atol=1e-6)
                        end
                    end
                end
            end

            X = Matrix(
                Manifolds.hat(SOn, pts[1], π * normalize(randn(manifold_dimension(SOn)))),
            )
            p = exp(SOn, pts[1], X)
            X2 = log(SOn, pts[1], p)
            @test p ≈ exp(SOn, pts[1], X2)
        end
    end
    @testset "Test AbstractManifold Point and Tangent Vector checks" begin
        M = Manifolds.Rotations(2)
        for p in [1, [2.0 0.0; 0.0 1.0], [1.0 0.5; 0.0 1.0]]
            @test_throws DomainError is_point(M, p, true)
            @test !is_point(M, p)
        end
        p = one(zeros(2, 2))
        @test is_point(M, p)
        @test is_point(M, p, true)
        for X in [1, [0.0 1.0; 0.0 0.0]]
            @test_throws DomainError is_vector(M, p, X, true)
            @test !is_vector(M, p, X)
        end
        X = [0.0 1.0; -1.0 0.0]
        @test is_vector(M, p, X)
        @test is_vector(M, p, X, true)
    end
    @testset "Project point" begin
        M = Manifolds.Rotations(2)
        p = Matrix{Float64}(I, 2, 2)
        p1 = project(M, p)
        @test is_point(M, p1, true)

        M = Manifolds.Rotations(3)
        p = collect(reshape(1.0:9.0, (3, 3)))
        p2 = project(M, p)
        @test is_point(M, p2, true)

        rng = MersenneTwister(44)
        x3 = project(M, randn(rng, 3, 3))
        @test is_point(M, x3, true)
    end
    @testset "Convert from Lie algebra representation of tangents to Riemannian submanifold representation" begin
        M = Manifolds.Rotations(3)
        p = project(M, collect(reshape(1.0:9.0, (3, 3))))
        x = [[0, -1, 3] [1, 0, 2] [-3, -2, 0]]
        @test is_vector(M, p, x, true)
        @test embed(M, p, x) == p * x
        Y = zeros((3, 3))
        embed!(M, Y, p, x)
        @test Y == p * x
        @test Y ≈ p * (p'Y - Y'p) / 2
    end
    @testset "Edge cases of Rotations" begin
        @test_throws OutOfInjectivityRadiusError inverse_retract(
            Rotations(2),
            [1.0 0.0; 0.0 1.0],
            [1.0 0.0; 0.0 -1.0],
            PolarInverseRetraction(),
        )
    end
end