text:
function Base.show(io::IO, dsr::DatasetRow;
                   allcols::Bool = !get(io, :limit, false),
                   rowlabel::Symbol = :Row,
                   eltypes::Bool = true,
                   truncate::Int = 32,
                   kwargs...)
    r, c = parentindices(dsr)
    _show(io, view(parent(dsr), [r], c); allcols=allcols, rowlabel=rowlabel,
          summary=false, rowid=r, eltypes=eltypes, truncate=truncate,
          title="DatasetRow", kwargs...)
end

Base.show(io::IO, mime::MIME"text/plain", dsr::DatasetRow;
          allcols::Bool = !get(io, :limit, false),
          rowlabel::Symbol = :Row,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    show(io, dsr; allcols=allcols, rowlabel=rowlabel, eltypes=eltypes,
         truncate=truncate, kwargs...)

Base.show(dsr::DatasetRow;
          allcols::Bool = !get(stdout, :limit, true),
          rowlabel::Symbol = :Row,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    show(stdout, dsr; allcols=allcols, rowlabel=rowlabel, eltypes=eltypes,
         truncate=truncate, kwargs...)


text:
# A very complicated example to use `BEE.jl`
#
# We will *prove* the following fact

using BEE

@beeint x  0 5
@beeint y -4 9
@beeint z -5 10

x + y == z

@beeint w 0 10

xl = [beebool("x$i") for i=1:4]

xl[1] == -xl[2]
xl[2] == true

sum([-xl[1], xl[2], -xl[3], xl[4]]) == w

BEE.render()

ext: jl

max_stars_repo_name: newptcai/BeeEncoder

max_issues_repo_name: newptcai/BeeEncoder

max_forks_repo_name: newptcai/BeeEncoder


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "iso_codes"
version = v"4.3"

# Collection of sources required to build iso-codes
sources = [
    "https://salsa.debian.org/iso-codes-team/iso-codes/-/archive/iso-codes-$(version.major).$(version.minor)/iso-codes-iso-codes-$(version.major).$(version.minor).tar.bz2" =>
    "6b539f915d02c957c45fce8133670811f1c36a1f1535d5af3dd95dc519d3c386"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/iso-codes-*/
apk add gettext
./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    FileProduct("share/iso-codes", :iso_codes_dir),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
## CachedOperator

mutable struct CachedArray{T,N,DM<:AbstractArray{T,N},M<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::DM
    array::M
    datasize::NTuple{N,Int}
end

const CachedVector{T,DM<:AbstractVector{T},M<:AbstractVector{T}} = CachedArray{T,1,DM,M}
const CachedMatrix{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} = CachedArray{T,2,DM,M}

# CachedArray(data::AbstractArray{T,N}, array::AbstractArray{T,N}, sz::NTuple{N,Int}) where {T,N} =
#     CachedArray{T,N,typeof(data),typeof(array)}(data, array, sz)
CachedArray(data::AbstractArray, array::AbstractArray) = CachedArray(data, array, size(data))

CachedArray(::Type{Array}, array::AbstractArray{T,N}) where {T,N} =
    CachedArray(Array{T,N}(undef, ntuple(zero,N)), array)

CachedArray(array::AbstractArray{T,N}) where {T,N} =
    CachedArray(similar(array, ntuple(zero,N)), array)

"""
    cache(array::AbstractArray)

Caches the entries of an array.
"""
cache(O::AbstractArray) = CachedArray(O)
cache(::Type{MT}, O::AbstractArray) where {MT<:AbstractArray} = CachedArray(MT,O;kwds...)

convert(::Type{AbstractArray{T}}, S::CachedArray{T}) where {T} = S
convert(::Type{AbstractArray{T}}, S::CachedArray) where {T} =
    CachedArray(convert(AbstractArray{T}, S.data), convert(AbstractArray{T}, S.array), S.datasize)

size(A::CachedArray) = size(A.array)

@propagate_inbounds function Base.getindex(B::CachedArray{T,N}, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(Bool, B, kj...)
    resizedata!(B, kj...)
    B.data[kj...]
end

@propagate_inbounds function Base.setindex!(B::CachedArray{T,N}, v, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(Bool, B, kj...)
    resizedata!(B,kj...)
    @inbounds B.data[kj...] = v
    v
end

## Array caching

function resizedata!(B::CachedArray{T,N,Array{T,N}},nm::Vararg{Integer,N}) where {T<:Number,N}
    @boundscheck checkbounds(Bool, B, nm...) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = B.data
    νμ = size(olddata)
    nm = max.(νμ,nm)
    if νμ ≠ nm
        B.data = Array{T}(undef, nm)
        B.data[axes(olddata)...] = olddata
    end

    for k in 1:N-1
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k],Base.OneTo.(B.datasize[k+1:end])...)
        B.data[inds...] .= view(B.array,inds...)
    end
    let k = N
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k])
        B.data[inds...] .= view(B.array,inds...)
    end
    B.datasize = nm

    B
end


text:
# # Example 2 on MibS

# Model of the problem
# First level
# ```math
# \min_{x} 2x -4y + 10z,\\
# \notag s.t.\\
# -3x + 2y + 2z \leq 12,\\
# x + 2y \leq 20,\\
# x \leq 10,\\
# x \in \mathbb{Z}, z \in \mathbb{B},\\
# ```
# Second level
# ```math
# \min_{y} y,\\
# \notag s.t.\\
# 2x - y + 3z<= 7,\\
# -2x + 4y <= 16,\\
# y <= 5\\
# y \in \mathbb{Z}\\
# ```

using BilevelJuMP
using JuMP
using Test

model = BilevelModel()

# First we need to create all of the variables in the upper and lower problems:

# Upper level variables
@variable(Upper(model), x, Int)
@variable(Upper(model), z, Bin)

#Lower level variables
@variable(Lower(model), y, Int)

# Then we can add the objective and constraints of the upper problem:

# Upper level objecive function
@objective(Upper(model), Min, 2x - 4y + 10z)

# Upper constraints
@constraints(Upper(model), begin
    u1, -3x + 2y + 5z  <= 12
    u2, x + 2y <= 20
    u3, x <= 10
end)

# Followed by the objective and constraints of the lower problem:

# Lower objective function
@objective(Lower(model), Min, y)

# Lower constraints
@constraint(Lower(model), l1,  2x -  y + 3z <= 7)
@constraint(Lower(model), l2, -2x +  4y <= 16)
@constraint(Lower(model), l3, y <= 5)

# Using MibS Solver
solution = BilevelJuMP.solve_with_MibS(model, silent = true)

# Auto testing
@test solution.status == true
@test solution.objective ≈ -8.0
@test solution.nonzero_upper == Dict(0 => 6.0)
@test solution.nonzero_lower == Dict(0 => 5.0)
@test solution.all_upper["x"] == 6.0
@test solution.all_upper["z"] == 0
@test solution.all_lower["y"] == 5.0


text:
# Unit tests for implementations/molly/molly_base.jl

@testset "implementations/molly/molly_base.jl" begin
    particles = [
        AtomsBase.Atom(:Ar, [-2.0, -1.0, 3.0]u"bohr", [3.0, 5.0, 21.0]u"bohr * hartree / ħ_au"; a=:b),
        AtomsBase.Atom(:Ar, [-7.0, 2.0, 13.0]u"bohr", [-3.0, 7.0, 1.0]u"bohr * hartree / ħ_au"; c=:d)
    ]
    box = (@SVector [(@SVector [1.5, 0.0, 0.0]), (@SVector [0.0, 1.5, 0.0]), (@SVector [0.0, 0.0, 1.5])])u"bohr"
    boundary_conditions = @SVector [Periodic(), Periodic(), Periodic()]
    flexible = FlexibleSystem(particles, box, boundary_conditions)
    fast = FastSystem(particles, box, boundary_conditions)

    system1 = System(flexible)
    system2 = System(fast)
    system3 = System(system1; loggers=Dict("c" => CoordinateLogger(1)))

    @test system1.atoms == [
        Molly.Atom(; index=1, mass=39.9481u"u")
        Molly.Atom(; index=2, mass=39.9481u"u")
    ]
    @test system1.atoms_data == [
        AugmentedAtomData(:Ar, Dict{Symbol,Any}(:a => :b)),
        AugmentedAtomData(:Ar, Dict{Symbol,Any}(:c => :d))
    ]
    @test system1.coords == [
        (@SVector [-2.0, -1.0, 3.0])u"bohr",
        (@SVector [-7.0, 2.0, 13.0])u"bohr"
    ]
    @test system1.velocities == [
        (@SVector [3.0, 5.0, 21.0])u"bohr * hartree / ħ_au",
        (@SVector [-3.0, 7.0, 1.0])u"bohr * hartree / ħ_au"
    ]
    @test system1.box_size == (@SVector [1.5, 1.5, 1.5])u"bohr"
    @test system1.force_units == u"hartree / bohr"
    @test system1.energy_units == u"hartree"

    @test system2.atoms == [
        Molly.Atom(; index=1, mass=39.9481u"u")
        Molly.Atom(; index=2, mass=39.9481u"u")
    ]
    @test system2.atoms_data == [
        AugmentedAtomData(:Ar, Dict{Symbol,Any}()),
        AugmentedAtomData(:Ar, Dict{Symbol,Any}())
    ]
    @test system2.coords == [
        (@SVector [-2.0, -1.0, 3.0])u"bohr",
        (@SVector [-7.0, 2.0, 13.0])u"bohr"
    ]
    @test system2.velocities == [
        (@SVector [0.0, 0.0, 0.0])u"bohr * hartree / ħ_au",
        (@SVector [0.0, 0.0, 0.0])u"bohr * hartree / ħ_au"
    ]
    @test system2.box_size == (@SVector [1.5, 1.5, 1.5])u"bohr"
    @test system2.force_units == u"hartree / bohr"
    @test system2.energy_units == u"hartree"

    @test system3.atoms == [
        Molly.Atom(; index=1, mass=39.9481u"u")
        Molly.Atom(; index=2, mass=39.9481u"u")
    ]
    @test system3.atoms_data == [
        AugmentedAtomData(:Ar, Dict{Symbol,Any}(:a => :b)),
        AugmentedAtomData(:Ar, Dict{Symbol,Any}(:c => :d))
    ]
    @test system3.coords == [
        (@SVector [-2.0, -1.0, 3.0])u"bohr",
        (@SVector [-7.0, 2.0, 13.0])u"bohr"
    ]
    @test system3.velocities == [
        (@SVector [3.0, 5.0, 21.0])u"bohr * hartree / ħ_au",
        (@SVector [-3.0, 7.0, 1.0])u"bohr * hartree / ħ_au"
    ]
    @test system3.box_size == (@SVector [1.5, 1.5, 1.5])u"bohr"
    @test system3.force_units == u"hartree / bohr"
    @test system3.energy_units == u"hartree"
    @test system3.loggers["c"] isa typeof(CoordinateLogger(1))

    data = AugmentedAtomData(:Ar, a=:b)

    @test hasproperty(data, :element)
    @test hasproperty(data, :a)
    @test !hasproperty(data, :b)
    @test data.element == :Ar
    @test data.a == :b
    @test_throws KeyError data.b
    @test propertynames(data) == (:element, :a)
    @test propertynames(data, true) == (:element, :data, :a)

    system = System(;
        atoms=[Molly.Atom(; index=1, mass=auconvert(39.9481u"u"))],
        atoms_data=[data],
        coords=[(@SVector [-2.0, -1.0, 3.0])u"bohr"],
        box_size=1.5u"bohr"
    )
    view = system[1]

    @test hasproperty(view, :a)
    @test !hasproperty(view, :b)
    @test view.a == :b
    @test_throws KeyError view.b
    @test propertynames(view) == (:system, :index, :a)
    @test propertynames(view, true) == (:system, :index, :a)

    atom = AtomsBase.Atom(data, (@SVector [-2.0, -1.0, 3.0])u"bohr", (@SVector [3.0, 5.0, 21.0])u"bohr * hartree / ħ_au")

    @test atom == particles[1]

    @test atomic_symbol(system1, 1) == :Ar
end


text:
using CUDArt, PyPlot

CUDArt.init([0])

md=CuModule("vadd.ptx",false)
kernel=CuFunction(md,"vadd")

function vadd(A::CudaArray,B::CudaArray,C::CudaArray)
    nblocks=round(Int,ceil(length(B)/1024))
    launch(kernel,nblocks,1024,(length(A),A,B,C))
end

N=2000
A=CudaArray(rand(N,N))
B=CudaArray(rand(N,N))
C=CudaArray(zeros(N,N))

M=5000
tm=zeros(M)
hA=to_host(A);
hB=to_host(B);
hC=to_host(C);

GPU=true
tic()
for i in 1:M
    if GPU vadd(A,B,C)
    else
        hC=hA+hB
    end
end
if GPU CUDArt.device_synchronize() end
#to_host(C)
avtime=toc()/M;
println("average time per vadd = $avtime")

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
    kets_to_dm(kets; parallel_type, traceout)

Parallely computes a density matrix from an array of kets denoted `kets`.
If an array of indices is provided through the key argument `traceout`, then the
reduced density matrix is reconstructed only on the complementary subspace.

# Arguments
* `kets`: 1-dimensional array of kets.
* `parallel_type=:none`: The type of parallelism to employ. The types of
parallelism included are: `:none`, `:threads`, `:split_threads`, `:parfor` and
`:pmap`.
* `traceout`: Indices of some subspace to be traced out.

See also: [`kets_to_obs`](@ref)
"""
function kets_to_dm(kets::Array{T,1}; parallel_type::Symbol = :none,
        traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    @assert all([ket.basis == kets[1].basis for ket in kets]) "All kets must share a common basis"

    if parallel_type == :none
        𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
        ρ = ismissing(traceout) ? DenseOperator(first(kets).basis) : DenseOperator(𝒫(first(kets).basis));
        for i in 1:length(kets)
            @inbounds ρ.data += 𝒫(kets[i]).data;
        end
        return ρ / length(kets);
    elseif parallel_type == :threads
        return dm_threads(kets; traceout=traceout);
    elseif parallel_type == :pmap
        return dm_pmap(kets; traceout=traceout);
    elseif parallel_type == :parfor
        return dm_parfor(kets; traceout=traceout);
    elseif parallel_type == :split_threads
        return dm_split_threads(kets; traceout=traceout);
    else
        valptypes = [:none, :threads, :pmap, :parfor, :split_threads];
        error("Invalid parallel type. Type :$parallel_type not available.\n"*
              "Available types are: "*reduce(*,[":$t " for t in valptypes]))
    end
end;

function dm_threads(kets::Array{T,1}; traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    N = length(kets);
    𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
    b = ismissing(traceout) ? first(kets).basis : (𝒫(first(kets).basis)).basis_l;
    ρs = [DenseOperator(b) for i in 1:Threads.nthreads()];
    Threads.@threads for i in 1:length(kets)
        @inbounds ρs[Threads.threadid()].data .+= 𝒫(kets[i]).data;
    end
    @inbounds for i in 2:length(ρs)
        ρs[1].data .+= ρs[i].data;
    end
    return ρs[1]/N;
end

function dm_pmap(kets::Array{T,1}; traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    N = length(kets);
    𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
    wp = CachingPool(workers())
    r = RemoteChannel(1)
    b = ismissing(traceout) ? first(kets).basis : (𝒫(first(kets).basis)).basis_l;
    put!(r, b)
    @sync for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :ρ, DenseOperator(fetch(r)))))
    end
    pmap(wp,1:N;batch_size=cld(N,length(wp.workers))) do i
        @inbounds ρ.data .+= 𝒫(kets[i]).data;
    end
    R = DenseOperator(b);
    for w in wp.workers
        @inbounds R.data .+= @fetchfrom w ρ.data
    end
    @sync for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :ρ, nothing)))
    end
    return R/N;
end

function dm_parfor(kets::Array{T,1}; traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    N = length(kets);
    𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
    b = ismissing(traceout) ? first(kets).basis : (𝒫(first(kets).basis)).basis_l;
    wp = CachingPool(workers())
    r = RemoteChannel(1)
    put!(r, b)
    @sync for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :ρ, DenseOperator(fetch(r)))))
    end
    # @distributed (+) is here much slower
    @sync @distributed for i in 1:N
        @inbounds ρ.data .+= 𝒫(kets[i]).data;
    end
    R = DenseOperator(b);
    for w in wp.workers
        @inbounds R.data .+= @fetchfrom w ρ.data
    end
    @sync for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :ρ, nothing)))
    end
    return R/N;
end

function dm_split_threads(kets::Array{T,1}; traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    N = length(kets);
    𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
    wp = CachingPool(workers())
    r = RemoteChannel(1)
    b = ismissing(traceout) ? first(kets).basis : (𝒫(first(kets).basis)).basis_l;
    put!(r, b)
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :hb, fetch(r))));
    end
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, :(ρs = [DenseOperator(hb) for i in 1:Threads.nthreads()])));
    end
    batches = nfolds(1:length(kets),length(wp.workers))

    pmap(wp,1:length(wp.workers)) do i
        Threads.@threads for j in 1:length(batches[i])
            @inbounds ρs[Threads.threadid()].data .+= dm(kets[batches[i][j]]).data;
        end
        @inbounds for i in 2:length(ρs)
            ρs[1].data .+= ρs[i].data;
        end
    end
    R = DenseOperator(b);
    for w in wp.workers
        @inbounds R.data .+= @fetchfrom w ρs[1].data
    end
    @sync for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :ρs, nothing)))
    end
    return R/N;
end

# Surprisingly, slightly faster
function dm_split_threads_v2(kets::Array{T,1}; traceout::Union{Vector{Integer},Missing}=missing) where {T<:StateVector}
    N = length(kets);
    ρ = ismissing(traceout) ? DenseOperator(first(kets).basis) : DenseOperator(𝒫(first(kets).basis));
    𝒫(x) = ismissing(traceout) ? dm(x) : ptrace(x, traceout);
    wp = CachingPool(workers());
    remch = RemoteChannel(()->Channel{Any}(Inf)); # TO DO: add some finite buffer size
    acc = @async begin
        @inbounds for i in 1:Threads.nthreads()*length(wp.workers)
            ρ.data .+= take!(remch).data;
        end
        nothing
    end
    batches = nfolds(1:length(kets),length(wp.workers))
    pmap(wp,1:length(wp.workers)) do i
        ρs = ismissing(traceout) ? [DenseOperator(first(kets).basis) for i in 1:Threads.nthreads()] : [DenseOperator(𝒫(first(kets).basis)) for i in 1:Threads.nthreads()];
        Threads.@threads for ket in kets[batches[i]]
            @inbounds ρs[Threads.threadid()].data .+= 𝒫(ket).data;
        end
        for j in 1:length(ρs)
            put!(remch,ρs[j])
        end
    end
    fetch(acc);
    clear!(wp);
    return ρ/N;
end

"""
    kets_to_obs(op, kets; parallel_type, index)

Computes the quantum average of an operator over an array of kets
denoted `kets`. If an index is provided through the key argument `index`, then
the operator is only evaluated on the corresponding subspace.

# Arguments
* `op`: Arbitrary Operator.
* `parallel_type=:none`: The type of parallelism to employ. The types of
parallelism included are: `:none`, `:threads`, `:split_threads`, `:parfor` and
`:pmap`. In practice, use only `:threads` or `:none`.
* `index`: Indices of the subspace where one wants to evaluate `op`.

See also: [`kets_to_dm`](@ref)
"""
function kets_to_obs(op::AbstractOperator, kets::Array{T,1}; parallel_type::Symbol = :none,
        index=missing) where {T<:StateVector}
    if parallel_type == :none
        obs = zero(ComplexF64);
        for i in 1:length(kets)
            @inbounds obs += ismissing(index) ? expect(op,kets[i]) : expect(index,op,kets[i]);
        end
        return obs / length(kets);
    elseif parallel_type == :threads
        expect_threads(op, kets; index=missing);
    elseif parallel_type == :pmap
        expect_pmap(op, kets; index=missing)
    elseif parallel_type == :parfor
        expect_parfor(op, kets; index=missing)
    elseif parallel_type == :split_threads
        expect_split_threads(op, kets; index=missing)
    else
        valptypes = [:none, :threads, :pmap, :parfor, :split_threads];
        error("Invalid parallel type. Type :$parallel_type not available.\n"*
              "Available types are: "*reduce(*,[":$t " for t in valptypes]))
    end
end;

function expect_threads(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    obss = zeros(ComplexF64,Threads.nthreads());
    Threads.@threads for i in 1:length(kets)
        @inbounds obss[Threads.threadid()] += ismissing(index) ? expect(op,kets[i]) : expect(index,op,kets[i]);
    end
    for i in 2:Threads.nthreads()
        @inbounds obss[1] += obss[i];
    end
    return obss[1] / length(kets);
end

function expect_pmap(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    N = length(kets);
    wp = CachingPool(workers())
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :obs_p, zeros(ComplexF64,1))));
    end
    pmap(wp,1:N;batch_size=cld(N,length(wp.workers))) do i
        @inbounds obs_p[] += ismissing(index) ? expect(op,kets[i]) : expect(index,op,kets[i]);
    end
    obs = zero(ComplexF64);
    @inbounds for w in wp.workers
        obs += @fetchfrom w obs_p[];
    end
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :obs_p, nothing)))
    end
    return obs / N;
end

function expect_parfor(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    obs = @distributed (+) for i in 1:length(kets)
        @inbounds ismissing(index) ? expect(op,kets[i]) : expect(index,op,kets[i]);
    end
    return obs / length(kets);
end

function expect_parfor_v2(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    N = length(kets);
    wp = CachingPool(workers())
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :obs_p, zeros(ComplexF64,1))));
    end
    @sync @distributed for i in 1:N
        @inbounds obs_p[] += ismissing(index) ? expect(op,kets[i]) : expect(index,op,kets[i]);
    end
    obs = zero(ComplexF64);
    @inbounds for w in wp.workers
        obs += @fetchfrom w obs_p[];
    end
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :obs_p, nothing)))
    end
    return obs / N;
end

function expect_split_threads(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    N = length(kets);
    wp = CachingPool(workers())
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, :(obss = zeros(ComplexF64,Threads.nthreads()))));
    end
    batches = nfolds(1:length(kets),length(wp.workers))
    pmap(wp,1:length(wp.workers)) do i
        Threads.@threads for j in 1:length(batches[i])
            @inbounds obss[Threads.threadid()] += ismissing(index) ? expect(op,kets[batches[i][j]]) : expect(index,op,kets[batches[i][j]]);
        end
        for i in 2:Threads.nthreads()
            @inbounds obss[1] += obss[i];
        end
    end
    obs = zero(ComplexF64);
    for w in wp.workers
        obs += @fetchfrom w obss[1];
    end
    @sync @async for w in wp.workers
        @spawnat(w, Core.eval(@__MODULE__, Expr(:(=), :obs_p, nothing)))
    end
    return obs / N;
end

# Slightly slower
function expect_split_threads_v2(op::AbstractOperator, kets::Array{T,1}; index=missing) where {T<:StateVector}
    N = length(kets);
    wp = CachingPool(workers());
    remch = RemoteChannel(()->Channel{Any}(Inf)); # TO DO: add some finite buffer size
    obs = zero(ComplexF64);
    acc = @async begin
        for i in 1:Threads.nthreads()*length(wp.workers)
            obs += take!(remch);
        end
        nothing
    end
    batches = nfolds(1:length(kets),length(wp.workers))
    pmap(wp,1:length(wp.workers)) do i
        obss = zeros(ComplexF64,Threads.nthreads());
        Threads.@threads for j in 1:length(batches[i])
            @inbounds obss[Threads.threadid()] += ismissing(index) ? expect(op,kets[batches[i][j]]) : expect(index,op,kets[batches[i][j]]);
        end
        for j in 1:length(obss)
            put!(remch,obss[j])
        end
    end
    fetch(acc);
    clear!(wp);
    return obs / N;
end


text:
module DashCoreComponents
using DashBase

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "1.15.0"

include("dcc_checklist.jl")
include("dcc_confirmdialog.jl")
include("dcc_confirmdialogprovider.jl")
include("dcc_datepickerrange.jl")
include("dcc_datepickersingle.jl")
include("dcc_dropdown.jl")
include("dcc_graph.jl")
include("dcc_input.jl")
include("dcc_interval.jl")
include("dcc_link.jl")
include("dcc_loading.jl")
include("dcc_location.jl")
include("dcc_logoutbutton.jl")
include("dcc_markdown.jl")
include("dcc_radioitems.jl")
include("dcc_rangeslider.jl")
include("dcc_slider.jl")
include("dcc_store.jl")
include("dcc_tab.jl")
include("dcc_tabs.jl")
include("dcc_textarea.jl")
include("dcc_upload.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_core_components",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "async-datepicker.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-datepicker.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-dropdown.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-dropdown.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-graph.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-graph.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-highlight.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-highlight.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-markdown.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-markdown.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-slider.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-slider.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-upload.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-upload.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-datepicker.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-datepicker.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-dropdown.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-dropdown.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-graph.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-graph.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-highlight.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-highlight.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-markdown.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-markdown.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-slider.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-slider.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-upload.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-upload.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_core_components.min.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/dash_core_components.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_core_components.min.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/dash_core_components.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_core_components-shared.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/dash_core_components-shared.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_core_components-shared.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/dash_core_components-shared.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "plotly.min.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/plotly.min.js",
    dynamic = nothing,
    async = :eager,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-plotlyjs.js",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-plotlyjs.js",
    dynamic = nothing,
    async = :lazy,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-plotlyjs.js.map",
    external_url = "https://unpkg.com/dash-core-components@1.15.0/dash_core_components/async-plotlyjs.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
module Model

using StaticArrays

# parameter types
export Image, SkyPatch, PsfComponent,
       GalaxyComponent, GalaxyPrototype,
       PriorParams, UnconstrainedParams,
       CanonicalParams, BrightnessParams, StarPosParams,
       GalaxyPosParams, GalaxyShapeParams,
       PsfParams, RawPSF, CatalogEntry,
       init_source, ParamSet

# functions
export align

# constants
export band_letters, D, Ia, B,
       galaxy_prototypes, prior,
       shape_standard_alignment,
       brightness_standard_alignment,
       gal_shape_alignment,
       ids_names, ids_free_names,
       ids, ids_free, star_ids, gal_ids, gal_shape_ids, psf_ids, bids

import Base.convert
import Base.+
import Distributions
import FITSIO, WCS
import WCS.WCSTransform
import ForwardDiff
import ..Log

import Base.length

const cfgdir = joinpath(Pkg.dir("Celeste"), "cfg")

const band_letters = ['u', 'g', 'r', 'i', 'z']

# The number of bands (colors).
const B = length(band_letters)

include("model/light_source_model.jl")
include("model/psf_model.jl")
include("model/image_model.jl")
include("model/param_set.jl")
include("model/imaged_sources.jl")
include("model/wcs_utils.jl")

import ..SensitiveFloats: SensitiveFloat, clear!
include("bivariate_normals.jl")
include("model/fsm_util.jl")
include("model/log_prob.jl")

end  # module


text:
# ====================================================================
#                         PRIVATE FUNCTIONS
# ====================================================================

# ____________________________________________________________________
# Model
#
function getparams(model::Model)
    out = OrderedDict{Symbol, Parameter}()
    for (cname, wcomp) in model.comp
        for (pname, par) in getparams(wcomp)
            out[Symbol(cname, compsep, pname)] = par
        end
    end
    return out
end

getparamvalues(v::Union{Model, AbstractComponent}) = [par.val for (pname, par) in getparams(v)]
function setparamvalues!(v::Union{Model, AbstractComponent}, pval::Vector{Float64})
    i = 0
    for (pname, par) in getparams(v)
        i += 1
        par.val = pval[i]
    end
end

function prepareindex1D!(model::Model)
    out = Vector{Int}()
    push!(out, 0)
    for id in 1:length(model.instruments)
        for jj in 1:length(model.instruments[id].exprnames)
            (model.instruments[id].exprcmp[jj])  ||  (continue)
            push!(out, out[end] + length(model.instruments[id].exprevals[jj]))
        end
    end
    model.index1d = out
    model.buffer1d = Vector{Float64}(undef, model.index1d[end])
end

function data1D(model::Model, data::Vector{T}) where T<:AbstractMeasures
    @assert length(data) >= (length(model.index1d)-1) "Not enough dataset(s) for model"
    @assert length(data) <= (length(model.index1d)-1) "Too many dataset(s) for model"

    out = Vector{Measures_1D}()
    ii = 1
    for id in 1:length(model.instruments)
        for jj in 1:length(model.instruments[id].exprnames)
            (model.instruments[id].exprcmp[jj])  ||  (continue)
            tmp = length(model.instruments[id].exprevals[jj])
            @assert(length(data[ii]) == tmp,
                    "Length of dataset $ii do not match corresponding model: " *
                    string(length(data[ii])) * " != " * string(tmp))
            d1 = flatten(data[ii], model.instruments[id].domain)
            push!(out, d1)
            ii += 1
        end
    end
    return out
end

# ____________________________________________________________________
# Component
#
function getparams(wcomp::WComponent)
    out = OrderedDict{Symbol, Parameter}()
    for pname in fieldnames(typeof(wcomp.comp))
        par = getfield(wcomp.comp, pname)
        if typeof(par) == Parameter
            par.cfixed = (par.fixed  ||  wcomp.fixed)
            if par.log
                par.val  = log10(par.val)
                par.low  = log10(par.low)
                par.high = log10(par.high)
            end
            out[pname] = par
        elseif typeof(par) == Vector{Parameter}
            for i in 1:length(par)
                par[i].cfixed = (par[i].fixed  ||  wcomp.fixed)
                out[Symbol(pname, i)] = par[i]
            end
        end
    end
    return out
end

# ____________________________________________________________________
# Parameter
#
isequal(a::Parameter, b::Parameter) = ((a._private.cname == b._private.cname)  &&
                                       (a._private.pname == b._private.pname)  &&
                                       (a._private.index == b._private.index))

# ____________________________________________________________________
# Instruments
#
# Instrument constructors
Instrument(dom::AbstractDomain) =
    Instrument(false, flatten(dom), dom, "", ()->nothing, 0,
               Vector{Symbol}(), Vector{CompEvaluation}(),
               Vector{Symbol}(), Vector{Expr}(), Vector{Bool}(), Vector{Vector{Float64}}())

function Instrument(model::Model, domain::AbstractDomain,
                    labels::Vector{Symbol}, exprs::Vector{Expr}, cmp::Vector{Bool})
    (length(cmp) == 0)  && (cmp = fill(true, length(exprs)))
    @assert length(labels) == length(exprs)
    @assert length(cmp) == length(exprs)

    function parse_model_expr(expr::Union{Symbol, Expr}, cnames, accum=Vector{Symbol}())
        if typeof(expr) == Expr
            # Parse the expression to check which components are involved
            for i in 1:length(expr.args)
                arg = expr.args[i]

                if typeof(arg) == Symbol
                    if arg in cnames # if it is one of the model components...
                        push!(accum, arg)
                    end
                elseif typeof(arg) == Expr
                    parse_model_expr(arg, cnames, accum)
                end
            end
        else
            push!(accum, expr)
        end
        return accum
    end

    # Check which components are involved
    compinvolved = Vector{Symbol}()
    for expr in exprs
        compinvolved = unique([compinvolved; parse_model_expr(expr, keys(model.comp))])
    end

    # Sort involved components according to insertion order
    kk = keys(model.comp)
    sort!(compinvolved, lt=(a, b) -> findall(a .== kk)[1] < findall(b .== kk)[1])

    # Prepare the code for model evaluation
    code = Vector{String}()
    tmp = ""
    for (cname, wcomp) in model.comp
        for (pname, par) in getparams(wcomp)
            tmp *= ", $(cname)$(compsep)$(pname)::Float64"
        end
    end
    push!(code, "(_instr::Instrument $tmp, _unused_...) -> begin")
    push!(code, "  domain = _instr.ldomain")
    # The last argument, _unused_, allows to push! further
    # components in the model after an expression has already been
    # prepared
    for (cname, wcomp) in model.comp
        if cname in compinvolved
            i = findall(cname .== compinvolved)[1]
            tmp = ""
            for (pname, par) in getparams(wcomp)
                (par.expr != "")  &&  (push!(code, "  $(cname)$(compsep)$(pname) = " *
                                             replace(par.expr, "this$(compsep)" => "$(cname)$(compsep)")))
                tmp *= ", $(cname)$(compsep)$(pname)"
            end
            if isfunction(wcomp.comp)
                push!(code, "  $cname(args...) = _evaluate!(_instr.compevals[$i], domain $tmp, args...)")
            else
                push!(code, "  $cname = _evaluate!(_instr.compevals[$i], domain $tmp)")
            end
        end
    end

    exprevals = Vector{Vector{Float64}}()
    if length(exprs) > 1
        for i in 1:length(exprs)
            push!(code, "  " * string(labels[i]) * " = " * string(exprs[i]))
        end
        for i in 1:length(exprs)
            push!(code, "  if length(_instr.exprevals[$i]) == 0")
            push!(code, "    _instr.exprevals[$i]  = " * string(labels[i]))
            push!(code, "  else")
            push!(code, "    _instr.exprevals[$i] .= " * string(labels[i]))
            push!(code, "  end")
            push!(exprevals, Vector{Float64}(undef, 0))
        end
    else
        i = 1
        push!(code, "  if length(_instr.exprevals[$i]) == 0")
        push!(code, "    _instr.exprevals[$i]  = " * string(exprs[i]))
        push!(code, "  else")
        push!(code, "    _instr.exprevals[$i] .= " * string(exprs[i]))
        push!(code, "  end")
        push!(exprevals, Vector{Float64}(undef, 0))
    end
    push!(code, "  _instr.counter += 1")
    push!(code, "  return nothing")
    push!(code, "end")

    funct = eval(Meta.parse(join(code, "\n")))

    compevals = CompEvals(model, domain, compinvolved)
    instr = Instrument(false, flatten(domain), domain, join(code, "\n"), funct, 0,
                       compinvolved, compevals,
                       deepcopy(labels), deepcopy(exprs), cmp, exprevals)
    return instr
end

function CompEvals(model::Model, domain::AbstractDomain, compinvolved::Vector{Symbol})
    compevals = Vector{CompEvaluation}()
    for (cname, wcomp) in model.comp
        if cname in compinvolved
            cd = cdata(wcomp.comp, domain)
            npar = length(getparams(wcomp))
            tmp = CompEvaluation(true, npar, fill(false, npar), 0, cd,
                                 Vector{Float64}(undef, npar),
                                 Vector{Float64}(undef, outsize(cd, domain)), [NaN])
            tmp.lastParams .= NaN
            push!(compevals, tmp)
        end
    end
    return compevals
end

# ____________________________________________________________________
# Expressions
#
newexprlabel(model::Model, id::Int) = Symbol(:expr, length(model.instruments[id].exprs)+1)
newexprlabel(model::Model, id::Int, n::Int) = Symbol.(Ref(:expr), length(model.instruments[id].exprs) .+ collect(1:n))

addexpr!(model::Model                                 , expr::Expr         ; cmp=true) = addexpr!(model,  1, [newexprlabel(model, 1)]              , [expr]       ; cmp=[cmp])
addexpr!(model::Model                                 , symbol::Symbol     ; cmp=true) = addexpr!(model,  1, [newexprlabel(model, 1)]              , [:(+$symbol)]; cmp=[cmp])
addexpr!(model::Model                                 , exprs::Vector{Expr}; cmp=true) = addexpr!(model,  1,  newexprlabel(model, 1, length(exprs)), exprs        ; cmp=fill(cmp, length(exprs)))

addexpr!(model::Model, id::Int                        , expr::Expr         ; cmp=true) = addexpr!(model, id, [newexprlabel(model, 1)]              , [expr]       ; cmp=[cmp])
addexpr!(model::Model, id::Int                        , symbol::Symbol     ; cmp=true) = addexpr!(model, id, [newexprlabel(model, 1)]              , [:(+$symbol)]; cmp=[cmp])
addexpr!(model::Model, id::Int                        , exprs::Vector{Expr}; cmp=true) = addexpr!(model, id,  newexprlabel(model, 1, length(exprs)), exprs        ; cmp=fill(cmp, length(exprs)))

addexpr!(model::Model,          label::Symbol         , expr::Expr         ; cmp=true) = addexpr!(model,  1, [label]                               , [expr]       ; cmp=[cmp])
addexpr!(model::Model,          label::Symbol         , symbol::Symbol     ; cmp=true) = addexpr!(model,  1, [label]                               , [:(+$symbol)]; cmp=[cmp])
addexpr!(model::Model,          labels::Vector{Symbol}, exprs::Vector{Expr}; cmp=true) = addexpr!(model,  1,  labels                               , exprs        ; cmp=fill(cmp, length(exprs)))

function addexpr!(model::Model, id::Int, labels::Vector{Symbol}, exprs::Vector{Expr}; cmp=Vector{Bool}())
    instr = model.instruments[id]
    append!(instr.exprnames, labels)
    append!(instr.exprs, exprs)
    append!(instr.exprcmp, cmp)
    instr.compile = true
    return labels[end]
end

# ____________________________________________________________________
function _evaluate!(c::CompEvaluation, d::AbstractDomain, args...)
    if c.fixed  &&  (c.counter > 0)
        isfinite(c.value[1])  &&  (return c.value)
    else
        @assert c.npar <= length(args) # "less than" to account for compound components
        #(length(c.result) == 1)  &&  empty!(c.result) # reset after using `c.value`
        if c.npar > 0
            neweval = false
            for i in 1:c.npar
                (c.lastParams[i] != args[i])  &&  (neweval = true)
            end
            if neweval  ||  (c.counter == 0)
                c.lastParams .= args[1:c.npar]
                for i in 1:c.npar;  (c.log[i])  &&  (args[i] = 10. ^args[i]);  end
                c.counter += 1
                evaluate!(c.cdata, c.result, d, args...)
                for i in 1:c.npar;  (c.log[i])  &&  (args[i] = log10(args[i]));  end
            end
        elseif c.counter == 0
            c.counter += 1
            evaluate!(c.cdata, c.result, d, args...)
        end
    end
    return c.result
end

_evaluate!(instr::Instrument, pvalues::Vector{Float64}) =
    Base.invokelatest(instr.funct, instr, pvalues...)

function _evaluate!(model::Model, pvalues::Vector{Float64})
    for ii in 1:length(model.instruments)
        _evaluate!(model.instruments[ii], pvalues)
    end
end

function _evaluate!(model::Model)
    for id in 1:length(model.instruments)
        instr = model.instruments[id]
        if instr.compile
            tmp = Instrument(model, instr.domain, instr.exprnames, instr.exprs, instr.exprcmp)
            deleteat!(model.instruments, id)
            insert!(model.instruments, id, tmp)
            _evaluate!(model)
            prepareindex1D!(model)
        end
    end

    for instr in model.instruments
        for ii in 1:length(instr.compnames)
            wcomp = model.comp[instr.compnames[ii]]
            instr.compevals[ii].fixed = wcomp.fixed
            jj = 1
            for (pname, par) in getparams(wcomp)
                instr.compevals[ii].log[jj] = par.log
            end
        end
    end

    _evaluate!(model, getparamvalues(model))

    for instr in model.instruments
        for ceval in instr.compevals
            ceval.value[1] = NaN
            if ceval.fixed
                mm = extrema(ceval.result)
                (mm[1] == mm[2])  &&  (ceval.value[1] = mm[1])
            end
        end
    end
    return nothing
end

# ____________________________________________________________________
# Minimizer
#
function pvalues2FitComp(model::Model, pvalues::Vector{Float64}, uncert::Vector{Float64})
    bestfit = OrderedDict{Symbol, FitComp}()
    ii = 0
    for (cname, comp) in model.comp
        fitcomp = OrderedDict{Symbol, Union{FitParam, Vector{FitParam}}}()
        accum = Vector{FitParam}()
        lastpname = :-
        for (pname, par) in getparams(comp)
            if (par._private.index == 0)  &&  (length(accum) > 0)
                fitcomp[lastpname] = deepcopy(accum)
                empty!(accum)
            end

            ii += 1
            par.val = pvalues[ii] # Update model parameter values
            par._private.fitval = pvalues[ii]
            par._private.fitunc = uncert[ii]
            if par._private.index == 0
                fitcomp[pname] = FitParam(pvalues[ii], uncert[ii], par.cfixed || par.fixed)
            else
                push!(accum, FitParam(pvalues[ii], uncert[ii], par.fixed))
                lastpname = par._private.pname
            end
        end
        (length(accum) > 0)  &&  (fitcomp[lastpname] = deepcopy(accum))
        bestfit[cname] = FitComp(fitcomp)
    end
    return bestfit
end

function residuals1d(model::Model, data1d::Vector{Measures_1D}) where T<:AbstractMeasures
    ii = 1
    for id in 1:length(model.instruments)
        for jj in 1:length(model.instruments[id].exprnames)
            (model.instruments[id].exprcmp[jj])  ||  (continue)
            i1 = model.index1d[ii]+1
            i2 = model.index1d[ii+1]
            model.buffer1d[i1:i2] .=
                ((model.instruments[id].exprevals[jj] .- data1d[ii].val) ./ data1d[ii].unc)
            ii += 1
        end
    end
    return model.buffer1d
end

function _fit(model::Model, data::Vector{T}; kw...) where T<:AbstractMeasures
    pval = getparamvalues(model)
    res = _fit!(model, data; kw...)
    setparamvalues!(model, pval)
    _evaluate!(model)
    return res
end

function _fit!(model::Model, data::Vector{T}; dry=false, minimizer=lsqfit()) where T<:AbstractMeasures
    elapsedTime = Base.time_ns()

    @assert typeof(minimizer) <: AbstractMinimizer
    @assert length(model.instruments) >= 1
    _evaluate!(model)
    params = collect(values(getparams(model)))
    pvalues = getfield.(params, :val)
    uncert = fill(NaN, length(pvalues))
    ifree = findall(.! getfield.(params, :cfixed))
    @assert length(ifree) > 0 "No free parameter in the model"

    # Inner function to evaluate all the models and store the result in a 1D array
    data1d = data1D(model, data)
    function eval_residuals1d(freepvalues::Vector{Float64},
                              model=model, pvalues=pvalues, ifree=ifree, data1d=data1d)
        pvalues[ifree] .= freepvalues
        _evaluate!(model, pvalues)
        return residuals1d(model, data1d)
    end
    eval_residuals1d(getfield.(params[ifree], :val))
    #Main.code_warntype(eval_residuals1d, (Vector{Float64},))

    status = :NonOptimal
    if !dry
        (status, bestfit_val, bestfit_unc) = minimize(minimizer, eval_residuals1d, params[ifree])
        @assert length(bestfit_val) == length(ifree)
        pvalues[ifree] .= bestfit_val
        uncert[ifree] .= bestfit_unc
        eval_residuals1d(pvalues[ifree]) # Re-evaluate model to ensure best fit values are used
    end

    cost = sum(abs2, model.buffer1d)
    dof = length(model.buffer1d) - length(ifree)
    result = FitResult(pvalues2FitComp(model, pvalues, uncert),
                       length(model.buffer1d), dof,
                       cost, status, logccdf(Chisq(dof), cost) * log10(exp(1)),
                       float(Base.time_ns() - elapsedTime) / 1.e9)
    return UI{FitResult}(result)
end

# ====================================================================
using LsqFit
mutable struct lsqfit <: AbstractMinimizer
end

function minimize(minimizer::lsqfit, func::Function, params::Vector{Parameter})
    ndata = length(func(getfield.(params, :val)))
    bestfit = LsqFit.curve_fit((dummy, pvalues) -> func(pvalues),
                               1.:ndata, fill(0., ndata),
                               getfield.(params, :val),
                               lower=getfield.(params, :low),
                               upper=getfield.(params, :high))
    status = :NonOptimal
    (bestfit.converged)  &&  (status = :Optimal)
    error = LsqFit.margin_error(bestfit, 0.6827)
    return (status, getfield.(Ref(bestfit), :param), error)
end

macro enable_CMPFit()
    return esc(:(
        import DataFitting.minimize;

        mutable struct cmpfit <: DataFitting.AbstractMinimizer;
        config::CMPFit.Config;
        cmpfit() = new(CMPFit.Config());
        end;

        function minimize(minimizer::cmpfit, func::Function, params::Vector{DataFitting.Parameter});
        guess = getfield.(params, :val);
        low   = getfield.(params, :low);
        high  = getfield.(params, :high);
        parinfo = CMPFit.Parinfo(length(guess));
        for i in 1:length(guess);
        llow  = isfinite(low[i])   ?  1  :  0;
        lhigh = isfinite(high[i])  ?  1  :  0;
        parinfo[i].limited = (llow, lhigh);
        parinfo[i].limits  = (low[i], high[i]);
        end;
        bestfit = CMPFit.cmpfit((pvalues) -> func(pvalues),
                                guess, parinfo=parinfo, config=minimizer.config);
        return (:Optimal, getfield.(Ref(bestfit), :param), getfield.(Ref(bestfit), :perror));
        end;
    ))
end


text:
using LinearAlgebra
using Test
using Random

using ProximalOperators
using ProximalAlgorithms

@testset "Lasso small (strongly convex, $T)" for T in [Float32, Float64]

    Random.seed!(777)
    dim = 5
    μf = T(1)
    Lf = T(10)

    x_star = convert(Vector{T}, 1.5 * rand(T, dim) .- 0.5)

    lam = (μf + Lf) / 2
    @test typeof(lam) == T

    D = Diagonal(sqrt(μf) .+ (sqrt(Lf) - sqrt(μf)) * rand(T, dim))
    D[1] = sqrt(μf)
    D[end] = sqrt(Lf)
    Q = qr(rand(T, (dim, dim))).Q
    A = Q * D * Q'
    b = A * x_star + lam * inv(A') * sign.(x_star)

    f = LeastSquares(A, b)
    h = NormL1(lam)

    TOL = T(1e-4)

    x0 = A \ b
    x0_backup = copy(x0)

    @testset "SFISTA" begin
        solver = ProximalAlgorithms.SFISTA(tol = TOL)
        y, it = solver(x0, f = f, h = h, Lf = Lf, μf = μf)
        @test eltype(y) == T
        @test norm(y - x_star) <= TOL
        @test it < 45
        @test x0 == x0_backup
    end

    @testset "ForwardBackward" begin
        solver = ProximalAlgorithms.ForwardBackward(tol = TOL)
        y, it = solver(x0, f = f, g = h, Lf = Lf)
        @test eltype(y) == T
        @test norm(y - x_star, Inf) <= TOL
        @test it < 110
        @test x0 == x0_backup
    end

    @testset "FastForwardBackward" begin
        solver = ProximalAlgorithms.FastForwardBackward(tol = TOL)
        y, it = solver(x0, f = f, g = h, Lf = Lf, m = μf)
        @test eltype(y) == T
        @test norm(y - x_star, Inf) <= TOL
        @test it < 45
        @test x0 == x0_backup
    end

    @testset "DRLS" begin
        solver = ProximalAlgorithms.DRLS(tol = TOL)
        u, v, it = solver(x0, f = f, g = h, muf = μf)
        @test eltype(u) == T
        @test eltype(v) == T
        @test norm(v - x_star, Inf) <= TOL
        @test it < 14
        @test x0 == x0_backup
    end

    @testset "PANOC" begin
        solver = ProximalAlgorithms.PANOC(tol = TOL)
        y, it = solver(x0, f = f, g = h, Lf = Lf)
        @test eltype(y) == T
        @test norm(y - x_star, Inf) <= TOL
        @test it < 45
        @test x0 == x0_backup
    end

    @testset "PANOCplus" begin
        solver = ProximalAlgorithms.PANOCplus(tol = TOL)
        y, it = solver(x0, f = f, g = h, Lf = Lf)
        @test eltype(y) == T
        @test norm(y - x_star, Inf) <= TOL
        @test it < 45
        @test x0 == x0_backup
    end

end


text:
using PSSFSS
using LinearAlgebra: norm, diagind
using PSSFSS.GSMs: GSM, cascade, cascade!, initialize_gsm_file,
                  append_gsm_data, read_gsm_file, choose_gblocks,
                  Gblock
using PSSFSS.Layers: Layer, TEorTM, TE, TM
using PSSFSS.PSSFSSLen
using PSSFSS.Elements: rectstrip
using PSSFSS.Sheets: RWGSheet
using PSSFSS.Constants: c₀, twopi
using Test

atol = 1e-16 # For comparison to zero

@testset "cascade1" begin
    na1 = 10
    na2 = nb1 = 20
    nb2 = 15

    a = GSM(na1, na2)
    b = GSM(nb1, nb2)
    c = cascade(a, b)

    @test norm(c.s11) == norm(c.s22) == 0
    @test all(c.s12[diagind(c.s12)] .== 1)
    @test all(c.s12[setdiff(LinearIndices(c.s12), diagind(c.s12))] .== 0)
    @test all(c.s21[diagind(c.s21)] .== 1)
    @test all(c.s21[setdiff(LinearIndices(c.s21), diagind(c.s21))] .== 0)
end

@testset "cascade2" begin
    na1 = 1
    na2 = nb1 = 1
    nb2 = 1
    α = 0.35
    β = sqrt(1 - α^2)
    a = GSM(na1, na2)
    a.s11 .= a.s22 .= α
    a.s12 .= a.s21 .= β

    b = GSM(nb1, nb2)
    γ = 0.23
    δ = sqrt(1 - γ^2)
    a = GSM(na1, na2)
    b.s11 .= b.s22 .= γ
    b.s12 .= b.s21 .= δ
    c = cascade(a, b)
    Δ = 1 / (1 - a.s22[1,1] * b.s11[1,1])
    s11 = a.s11 + a.s12 .* a.s21 .* b.s11 ./ Δ
    s12 = a.s12 .* b.s12 ./ Δ
    s21 = a.s21 .* b.s21 ./ Δ
    s22 = b.s22 + b.s12 .* b.s21 .* a.s22 ./ Δ
    @test norm(c.s11 - s11) ≤ atol
    @test norm(c.s12 - s12) ≤ atol
    @test norm(c.s21 - s21) ≤ atol
    @test norm(c.s22 - s22) ≤ atol
end

@testset "cascade!" begin
    na1 = 1
    na2 = nb1 = 1
    nb2 = 1
    α = 0.35
    β = sqrt(1 - α^2)
    a = GSM(na1, na2)
    a.s11 .= a.s22 .= α
    a.s12 .= a.s21 .= β
    b = deepcopy(a)

    width = 25mm
    layer = Layer(width=width)
    layer.γ = [1.0im]

    prop = exp(-layer.γ[1] * ustrip(Float64, u"m", width))
    s11 = a.s11
    s12 = a.s12 * prop
    s21 = a.s21 * prop
    s22 = b.s22 * prop^2

    cascade!(a, layer, ustrip(Float64, mm, width))
    cascade!(b, layer)

    @test norm(a.s11 - b.s11) ≤ atol
    @test norm(a.s12 - b.s12) ≤ atol
    @test norm(a.s21 - b.s21) ≤ atol
    @test norm(a.s22 - b.s22) ≤ atol

    @test norm(a.s11 - s11) ≤ atol
    @test norm(a.s12 - s12) ≤ atol
    @test norm(a.s21 - s21) ≤ atol
    @test norm(a.s22 - s22) ≤ atol
end

@testset "GSM file" begin
    na1 = 1
    na2 = nb1 = 1
    nb2 = 1
    α = 0.35
    β = sqrt(1 - α^2)
    a = GSM(na1, na2)
    a.s11 .= a.s22 .= α
    a.s12 .= a.s21 .= β

    l1 = Layer(ϵᵣ = 2.0)
    l1.P = [TE, TM]
    l1.M = [0, 0]
    l1.N = [0, 0]

    l2 = deepcopy(l1)
    l2.ϵᵣ = 1

    sh = rectstrip(Lx=1, Ly=1, Px=1, Py=1, Nx=2, Ny=1, units=mm, fufp=true)

    fname = joinpath(tempdir(), "gsmfile.gsm")

    initialize_gsm_file(fname, l1, l2, sh, sh)
    case = Dict("θ"=>0.0, "ϕ"=>0.0, "FGHz"=>2.0)
    append_gsm_data(fname, "1", a, l1, l2, case)
    dat = read_gsm_file(fname)
    @test dat["ModelistIn"] == (TEorTM[TE, TM], [0, 0], [0, 0])
    @test dat["ModelistOut"] == (TEorTM[TE, TM], [0, 0], [0, 0])
    @test dat["s1in"] == [1.0, 0.0]
    @test dat["s1out"] == [1.0, 0.0]
    @test dat["s2in"] == [0.0, 1.0]
    @test dat["s2out"] == [0.0, 1.0]
    @test dat["1/gsm"].s11 == a.s11
    @test dat["1/gsm"].s12 == a.s12
    @test dat["1/gsm"].s21 == a.s21
    @test dat["1/gsm"].s22 == a.s22
    @test dat["1/epsmu1"] == (2.0 - 0.0im, 1.0 - 0.0im)
    @test dat["1/epsmu2"] == (1.0 - 0.0im, 1.0 - 0.0im)
    @test dat["1/case"] == case
    @test dat["unitsin"] == mm
    @test dat["unitsout"] == mm
end

@testset "choose_gblocks1" begin
    strata =
    [Layer()
     rectstrip(units=cm, Px = 1, Py = 1, Lx=0.5, Ly=0.5, Nx=10, Ny=10)
     Layer(width=1mil, ϵᵣ=2.2)
     Layer(width=1mil)
     Layer(width=3mm)
     rectstrip(units=cm, Px = 1, Py = 1, Lx=0.5, Ly=0.5, Nx=10, Ny=10)
     Layer(width=1mil, ϵᵣ=2.2)
     Layer()
    ]
    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there

    FGHz = 10.0; k0 = twopi * FGHz*1e9 / c₀
    gbl = choose_gblocks(layers, sheets, junc, k0)
    @test length(gbl) == 2
    @test gbl[1] == Gblock(1:3,1)
    @test gbl[2] == Gblock(4:5,4)
end

@testset "choose_gblocks2" begin
    duroid = Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
    strata =
    [
        Layer(name="Vacuum")
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum", width=0.05inch)
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum")
    ]

    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀
    gbl = choose_gblocks(layers, sheets, junc, k0min)
    @test gbl[1] == Gblock(1:3,2)
    @test gbl[2] == Gblock(4:6,5)
end

@testset "choose_gblocks3" begin
    duroid = Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
    strata =
    [
        Layer(name="Vacuum")
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum")
    ]

    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀
    @test choose_gblocks(layers, sheets, junc, k0min) == [Gblock(1:2,2), Gblock(3:5,4)]
end

@testset "choose_gblocks4" begin
    duroid = Layer(name="Duroid 6035", ϵᵣ=3.6, tanδ=0.0013, width=20mil)
    strata =
    [
        Layer(name="Vacuum")
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum", width=5mil)
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum", width=5mil)
        duroid
        polyring(sides=4, orient=45, a=[0.1902], b=[0.2169], units=inch,
                    s1=[0.5117, 0], s2=[0, 0.5117], ntri=100)
        duroid
        Layer(name="Vacuum")
    ]

    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    layers = convert(Vector{Layer}, strata[islayer])
    sheets = convert(Vector{RWGSheet}, strata[issheet])
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    fmin = 2.0 * 1e9  # minimum frequency [Hz]
    k0min = twopi * fmin / c₀

    @test choose_gblocks(layers, sheets, junc, k0min) == [Gblock(1:2,2), Gblock(3:5,5),Gblock(6:9,8)]

end


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

module NetworkManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

# export models
export convert, AzureAsyncOperationResult
export convert, Error
export convert, ErrorDetails
export convert, Resource
export convert, SubResource

# export operations
export convert

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using RowMajorArrays
using Test

@testset "RowMajorArrays.jl" begin
    a = [1 2 3; 4 5 6]
    a_r = RowMajorArray(a)

    @test a[1, 3] == a_r[3, 1]
    @test size(a) == reverse(size(a_r)) == (2, 3)
    @test size(a_r, 2) == 2

    @test RowMajorArray{Float64}(a) == RowMajorArray(convert.(Float64, a))

    setindex!(a_r, 42, 3, 2)
    @test getindex(a_r, 3, 2) == 42
    setindex!(a_r, 43, 2,1)
    @test getindex(a_r, 2, 1) == 43

    @test a_r[4] == a[4] # linear indices are the same for standard and RowMajorArrays

    z = RowMajorArrays.zeros(Int, 5, 3)
    @test size(z) == (5, 3)
    @test zeros(Int, 3, 5) == z.data
    @test eltype(z) == Int
    @test all(z .== 0)

    z2 = RowMajorArrays.zeros(5, 3)
    @test size(z2) == (5, 3)
    @test zeros(3, 5) == z.data
    @test eltype(z2) == Float64
    @test all(z2 .== 0)

    o = RowMajorArrays.ones(Int, 5, 3)
    @test size(o) == (5, 3)
    @test ones(Int, 3, 5) == o.data
    @test all(o .== 1)

    r1 = RowMajorArrays.rand(7, 2)
    @test size(r1.data) == (2, 7)
    @test eltype(r1) == Float64
    r1b = RowMajorArrays.rand(Float32, 7, 2)
    @test size(r1b.data) == (2, 7)
    @test eltype(r1b) == Float32

    r2 = RowMajorArrays.randn(7, 2)
    @test size(r2.data) == (2, 7)
    @test eltype(r2) == Float64
    r2b = RowMajorArrays.randn(Float32, 7, 2)
    @test size(r2b.data) == (2, 7)
    @test eltype(r2b) == Float32

    a42 = RowMajorArrays.fill(42, 5, 3)
    @test size(a42) == (5, 3)
    @test a42[2, 1] == 42

    a42b = RowMajorArrays.fill(42, (5, 3))
    @test size(a42b) == (5, 3)
    @test a42 == a42b

    undef_array = RowMajorArray{Float64, 2, Array{Float64, 2}}(undef, 3, 2)
    @test size(undef_array) == (3, 2)
    undef_array2 = RowMajorArray{Float64, 3, Array{Float64, 3}}(undef, 5, 4, 3)
    @test size(undef_array2) == (5, 4, 3)
    @test size(undef_array2.data) == (3, 4, 5)

    a_r = RowMajorArray([1 2 3; 4 5 6])
    one_r = RowMajorArrays.ones(3, 2)
    @test a_r + one_r == RowMajorArray([2. 3. 4.; 5. 6. 7.])
    @test a_r - one_r == RowMajorArray([0. 1. 2.; 3. 4. 5.])

    b_r = RowMajorArray([1 2; 3 4; 4 6])
    @test (a_r * b_r).data == transpose(a_r.data) * transpose(b_r.data)

    # inconsistent array types
    a = [1 2 3; 4 5 6]
    @test RowMajorArray(a) == RowMajorArray{Int64, 2, Array{Int64, 2}}(a) # consistent types
    @test_throws TypeError RowMajorArray(a) == RowMajorArray{Int64, 3, Array{Int64, 2}}(a)
    @test_throws TypeError RowMajorArray(a) == RowMajorArray{Int64, 2, Array{Float64, 2}}(a)

    # similar
    similar_1 = similar(a_r)
    @test size(similar_1) == (3,2)
    @test typeof(similar_1) == typeof(a_r) == RowMajorArray{Int64,2,Array{Int64,2}}
    similar_2 = similar(a_r, (10,4))
    @test size(similar_2) == (10,4)
    @test typeof(similar_2) == typeof(a_r) == RowMajorArray{Int64,2,Array{Int64,2}}
    similar_3 = similar(a_r, Int64, (10,4)) # same eltype as original array
    @test size(similar_3) == (10,4)
    @test typeof(similar_3) == typeof(a_r) == RowMajorArray{Int64,2,Array{Int64,2}}
    similar_4 = similar(a_r, Float64, (10,4)) # other eltype than original array
    @test size(similar_4) == (10,4)
    @test typeof(similar_4) == RowMajorArray{Float64,2,Array{Float64,2}}

    # broadcasting
    @test typeof(one_r .* a_r) == RowMajorArray{Float64,2,Array{Float64,2}}
    @test typeof(one_r .+ a_r) == RowMajorArray{Float64,2,Array{Float64,2}}
    @test typeof(a_r .* one_r) == RowMajorArray{Float64,2,Array{Float64,2}}
    @test typeof(a_r .+ one_r) == RowMajorArray{Float64,2,Array{Float64,2}}

    # check that standard array broadcasting is not broken
    @test typeof(a .+ ones(2, 3)) == Array{Float64,2}
    @test typeof(ones(2, 3) .+ a) == Array{Float64,2}

    # test pretty printing - 2d RowMajorArrays are printed with reversed dimensions of the underlying array
    io = IOBuffer(append=true)
    print(io, a_r)
    @test read(io, String) == "[1 4; 2 5; 3 6]"

end


text:
export debian_file

function debian_file(name, args...)
  mkpath("debian")
  open(joinpath("debian", name), "w") do stream
    for u in args
      println(stream, u)
    end
  end
end

type Debianizer <: AbstractDebianizer end

# Creates a control file for a source package declaring one or more binary packages
function control(package::Package, debianize::AbstractDebianizer; kwargs...)
  function as_string(value, joiner=", ")
    if isa(value, Array) || isa(value, Tuple)
      join(value, joiner)
    else
      string(value)
    end
  end
  # Create list of lines
  args = String[
    "Source: $(package.name)",
    "Section: $(package.section)",
    "Priority: $(package.priority)"
  ]
  if length(package.email) > 0 && length(package.maintainer) > 0
    push!(args, "Maintainer: $(package.maintainer) <$(package.email)>")
  elseif length(package.maintainer) > 0
    push!(args, "Maintainer: $(package.maintainer)")
  end
  push!(args, "Standards-Version: 3.9.5")
  if length(package.homepage) > 0; push!(args, "Homepage: $(package.homepage)"); end

  build_depends = join(all_packages(package), ", ")
  if length(build_depends) > 0
    push!(args, "Build-Depends: debhelper (>=9), $build_depends")
  else
    push!(args, "Build-Depends: debhelper (>=9)")
  end

  for binary in deepcopy(package.binaries)
    name = pop!(binary, :package, package.name) |> as_string
    desc = pop!(binary, :description, "") |> x -> as_string(x, "\n")
    arch = pop!(binary, :architecture, "any") |> as_string
    push!(args, "")
    push!(args, "Package: $name")
    push!(args, "Architecture: $arch")
    for (key, value) in filter((k, v) -> k != :install, binary)
      key = "$(uppercase(string(key)[1]))$(lowercase(string(key)[2:end]))"
      push!(args, "$key: $(as_string(value))")
    end
    if length(desc) > 0; push!(args, "Description: $desc"); end
  end

  # now call function
  debian_file("control", args...)
end

# Adds installation files
function install(package::Package, debianize::AbstractDebianizer; kwargs...)
  for binary in package.binaries
    if !haskey(binary, :install); continue; end
    name = get(binary, :package, package.name)
    debian_file("$name.install", binary[:install]...)
  end
end

function change(package::Package, args...)
  if length(args) == 0
    args = ["  * No reported changes"]
  else
    args = ["  * $u" for u in args]
  end
  debian_file("changelog",
    "$(package.name) ($(package.version)) $(package.distribution); urgency=$(package.urgency)",
    "",
    args...,
    "",
    " -- $(package.maintainer) <$(package.email)>  $(strftime("%a, %d %b %Y %T %z", time()))"
  )
end

# Adds the different debian files
function debianize(pack::AbstractPackage, ::AbstractDebianizer, workdir="workspace"; kwargs...)
  const name = package_name(pack)
  const build = builddir(pack, workdir)
  cd(joinpath(build, name)) do
    control(pack)
    rules(pack)
    change(pack, pack.changes...)
    install(pack)
  end
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
@testset "frames" begin
    Random.seed!(7)
    f1name = "1"
    f1 = CartesianFrame3D(f1name)
    f2 = CartesianFrame3D()
    f3 = CartesianFrame3D()
    @test string(f1) == f1name
    string(f2) # just to make sure it doesn't crash
    @test f2 != f3
    @boundscheck begin
        # only throws when bounds checks are enabled:
        @test_throws ArgumentError @framecheck(f1, f2)
        @test_throws ArgumentError @framecheck(f2, f3)
        @test_throws ArgumentError @framecheck(f2, (f1, f3))
    end
    @framecheck(f1, f1)
    @framecheck(f2, (f2, f3))

    t1 = rand(Transform3D, f2, f1)
    @test isapprox(t1 * inv(t1), one(Transform3D, f1))
    @test isapprox(inv(t1) * t1, one(Transform3D, f2))

    @test isapprox(t1 * Point3D(Float64, f2), Point3D(f1, translation(t1)))

    p = rand(Point3D, Float64, f2)
    v = FreeVector3D(f2, p.v)
    @test isapprox(inv(t1) * (t1 * p), p)
    @test isapprox(inv(t1) * (t1 * v), v)
    @test isapprox(t1 * p - t1 * v, Point3D(f1, translation(t1)))

    @test_throws DimensionMismatch Point3D(f2, rand(2))
    @test_throws DimensionMismatch Point3D(f2, rand(4))

    @test isapprox(transform(p, t1), t1 * p)
    @test isapprox(t1 \ (t1 * p), p)

    @test isapprox(transform(v, t1), t1 * v)
    @test isapprox(t1 \ (t1 * v), v)

    @test isapprox(p, Point3D(p.frame, p.v[1], p.v[2], p.v[3]))
    @test isapprox(-v, zero(v) - v)

    @test isapprox(norm(v), sqrt(dot(v, v)))

    show(devnull, t1)
    show(devnull, p)
    show(devnull, v)
end


text:
module Causality

import Base: union
import CausalInference
import Combinatorics

using Base
using Reexport
@reexport using LightGraphs
@reexport using MetaGraphs
@reexport using SymbolicUtils

function _update_module_doc()
    path = joinpath(@__DIR__, "..", "README.md")
    text = read(path, String)
    text = replace(text, "```julia" => "```jldoctest")
    @doc text Causality
end
_update_module_doc()

SU = SymbolicUtils
export SU

include("syms.jl")
export P, Do, ¦

include("undirected_paths.jl")
include("multivar_predicate.jl")
include("graphs.jl")
export graph, name2node

include("rules.jl")
export rule2

include("simplify_rules.jl")

end # module


text:
export FitInput
export SimulateInput
export FitResult
export has_exogenous

"""
FitInput(parameters::Dict{String,Any},
            dependent::Vector{TimeSeries{T}},
            exogenous::Vector{TimeSeries{T}}) where T

FitInput is the only information needed by Fit Functions
Fit Functions will just use dependent, exogenous and parameters
All extra information passed to the Fit Function that is not the dataset must be inside the parameters Dict.
You can also include keys 'args' ands 'kwargs' inside the dictionary and they will be passed directly to the FitFunction
"""
mutable struct FitInput{T}
    parameters::Dict{String,Any}
    dependent::Vector{TimeSeries{T}}
    exogenous::Vector{TimeSeries{T}}

    function FitInput(parameters::Dict{String,Any},
                           dependent::Vector{TimeSeries{T}},
                           exogenous::Vector{TimeSeries{T}}) where T
        # Test if the dependent vector is empty.
        if isempty(dependent)
            throw(ErrorException("Must have at least one dependent time series."))
        end
        # Warn if the exogenous vector is empty.
        has_exogenous_variable = !isempty(exogenous)
        # Test if all the dependent have the same timestamps
        if !assert_vector_time_series_timestamps(dependent)
            throw(DimensionMismatch("dependent timestamps must be the same."))
        end
        # Test if all the exogenous have the same timestamps
        if has_exogenous_variable && (!assert_vector_time_series_timestamps(exogenous))
            throw(DimensionMismatch("exogenous timestamps must be the same."))
        end
        # Test if dependent and exogenous have the same timestaamps
        if has_exogenous_variable && (!assert_two_vectors_time_series_timestamps(dependent[1], exogenous[1]))
            throw(DimensionMismatch("exogenous and dependent timestamps must be the same."))
        end

        return new{T}(parameters,
                      dependent,
                      exogenous)
    end
end

"""
    has_exogenous(fit_input::FitInput{T}) where T

Returns true if fit_input has one or more exogenous time series and false otherwise.
"""
function has_exogenous(fit_input::FitInput{T}) where T
    return length(fit_input.exogenous) >= 1
end

"""
    FitResult(hyperparameters, other)

May be developed specifically for each model.
Ideally all the essential information that flows from fit to simulate should be in the hyperparameters.
"""
mutable struct FitResult
    hyperparameters
    other

    function FitResult(hyperparameters,
                       other)
        return new(hyperparameters, other)
    end
end

# Additional constructor for FitResult
function FitResult(hyperparameters)
    return FitResult(hyperparameters, nothing)
end

"""
    SimulateInput(fit_input::FitInput{T}
                  timestamps_forecast::Vector{DateTime}
                  exogenous_forecast::Vector{TimeSeries{T}}
                  fit_result::FitResult)

SimulateInput is the only information needed by Simulate Functions
FitResult may be developed specifically for each FitFunction or you can use the generic version.
Ideally all the essential information that flows from fit to simulate should be in the hyperparameters.
"""
mutable struct SimulateInput{T}
    fit_input::FitInput{T}
    timestamps_forecast::Vector{DateTime}
    exogenous_forecast::Vector{TimeSeries{T}}
    fit_result::FitResult

    function SimulateInput(fit_input::FitInput{T},
                           timestamps_forecast::Vector{DateTime},
                           exogenous_forecast::Vector{TimeSeries{T}},
                           fit_result::FitResult) where T

        # Having parameters, dependent and exogenous here makes possible to perform all the tests
        parameters = fit_input.parameters
        dependent = fit_input.dependent
        exogenous = fit_input.exogenous
        has_exogenous_variable = !isempty(exogenous)

        # Test if exogenous and exogenous_forecast have the same number of time series
        if !(length(exogenous) == length(exogenous_forecast))
            throw(ErrorException("exogenous and exogenous_forecast must have the same numbers of time_series."))
        end

        # Test if all the exogenous_forecast have the same timestamps
        if has_exogenous_variable && (!assert_vector_time_series_timestamps(exogenous_forecast))
            throw(DimensionMismatch("exogenous_forecast timestamps must be the same."))
        end
        # Warn if exogenous_forecast timestamps are not greater than dependent timestamps
        if has_exogenous_variable && (dependent[1].timestamps[end] >= exogenous_forecast[1].timestamps[1])
            @warn("timestamps of exogenous forecast are not greater than dependent timestamps.")
        end
        # Test if each exogenous_forecast timestamps is equal timestamps_forecast
        if has_exogenous_variable && (!all([isequal(timestamps_forecast, exogenous_forecast[i].timestamps) for i = 1:length(exogenous_forecast)]))
            throw(DimensionMismatch("timestamps_forecast must be equal to each exogenous_forecast.timestamps"))
        end
        return new{T}(fit_input,
                    timestamps_forecast,
                    exogenous_forecast,
                    fit_result)
    end
end

"""
    has_exogenous(simulate_input::SimulateInput{T}) where T

Returns true if simulate_input has one or more exogenous time series and false otherwise.
"""
function has_exogenous(simulate_input::SimulateInput{T}) where T
    return length(simulate_input.exogenous_forecast) >= 1
end

# Assert Functions
function assert_vector_time_series_timestamps(vector_ts::Vector{TimeSeries{T}}) where T
    reference = vector_ts[1].timestamps
    for ts in vector_ts
        if !(isequal(ts.timestamps, reference))
            return false
        end
    end
    return true
end

function assert_two_vectors_time_series_timestamps(reference_ts::TimeSeries{T},
                                                   comparision_ts::TimeSeries{T}) where T
    return isequal(reference_ts.timestamps, comparision_ts.timestamps)
end

function assert_length_time_series_timestamp(ts::TimeSeries{T}, len_timestamps::Int) where T
    return isequal(ts.timestamps, len_timestamps)
end


text:
module WeightedOnlineStats

export WeightedSum, WeightedMean,
    WeightedVariance, WeightedCovMatrix,
    WeightedHist, WeightedAdaptiveHist, WeightedAdaptiveBins,
    fit!, merge!, weightsum, value,
    mean, std, cov, cor, median, quantile

import OnlineStats
import OnlineStats: Tup, VectorOb,
    TwoThings,
    Algorithm, Extrema,
    smooth, smooth!, smooth_syr!
import OnlineStatsBase
import OnlineStatsBase:
    OnlineStat, name,
    fit!, merge!,
    _fit!, _merge!,
    eachrow, eachcol,
    nobs, value
import Statistics
import Statistics: mean, var, std, cov, cor, median, quantile
import LinearAlgebra
import LinearAlgebra: Hermitian, lmul!, rmul!, Diagonal, diag
import StatsBase: midpoints

include("interface.jl")
include("sum.jl")
include("mean.jl")
include("var.jl")
include("covmatrix.jl")
include("pca.jl")
include("histogram.jl")

end # module WeightedOnlineStats


text:
"""
    SumToZero(x<:AbstractIntVar, v::Int)

Summing constraint, states that `x[1] + x[2] + ... + x[length(x)] == 0`
"""
struct SumToZero <: Constraint
    x                   ::Array{<:AbstractIntVar}
    active              ::StateObject{Bool}
    numberOfFreeVars    ::StateObject{Int}
    sumOfFixedVars      ::StateObject{Int}
    freeIds             ::Array{Int}
    function SumToZero(x::Array{<:AbstractIntVar}, trailer)
        @assert !isempty(x)

        freeIds = zeros(length(x))
        for i in 1:length(x)
            freeIds[i] = i
        end

        constraint = new(x, StateObject{Bool}(true, trailer), StateObject{Int}(length(x), trailer), StateObject{Int}(0, trailer), freeIds)
        for xi in x
            addOnDomainChange!(xi, constraint)
        end
        return constraint
    end
end

"""
    propagate!(constraint::SumToZero, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`SumToZero` propagation function. The pruning is quite superficial.
"""
function propagate!(constraint::SumToZero, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    # Computing maxSum, minSum, and refreshing other variables
    newNumberOfFreeVars = constraint.numberOfFreeVars.value
    sumOfMax, sumOfMin = constraint.sumOfFixedVars.value, constraint.sumOfFixedVars.value
    for i in newNumberOfFreeVars:-1:1
        currentId = constraint.freeIds[i]
        sumOfMin += minimum(constraint.x[currentId].domain)
        sumOfMax += maximum(constraint.x[currentId].domain)

        if isbound(constraint.x[currentId])
            setValue!(constraint.sumOfFixedVars, constraint.sumOfFixedVars.value + assignedValue(constraint.x[currentId]))
            constraint.freeIds[i] = constraint.freeIds[newNumberOfFreeVars]
            constraint.freeIds[newNumberOfFreeVars] = currentId
            newNumberOfFreeVars -= 1
        end
    end
    setValue!(constraint.numberOfFreeVars, newNumberOfFreeVars)

    # Checking feasibility
    if sumOfMin > 0 || sumOfMax < 0
        return false
    end

    ### Filtering ###
    # Here we must have: x_i = - sum(x_j for j != i)
    # But we know that: - sum(x_j for j != i) <= - sum(min(x_j) for j != i)
    # And: - sum(min(x_j) for j != i) = min(x_i) - sumOfMin
    # Hence we remove everything above that last value from the domain of x_i
    # The reasoning is equivalent for the minimization
    for i in newNumberOfFreeVars:-1:1
        currentId = constraint.freeIds[i]
        currentMin = minimum(constraint.x[currentId].domain)
        currentMax = maximum(constraint.x[currentId].domain)
        pruned = vcat(removeAbove!(constraint.x[currentId].domain, currentMin - sumOfMin), removeBelow!(constraint.x[currentId].domain, currentMax - sumOfMax))
        if !isempty(pruned)
            addToPrunedDomains!(prunedDomains, constraint.x[currentId], pruned)
            triggerDomainChange!(toPropagate, constraint.x[currentId])

        end
    end

    if newNumberOfFreeVars == 0
        setValue!(constraint.active, false)
    end

    return true
end

variablesArray(constraint::SumToZero) = constraint.x

function Base.show(io::IO, ::MIME"text/plain", con::SumToZero)
    ids = [var.id for var in con.x]
    print(io, typeof(con), ": ", join(ids, " + "), " == 0, active = ", con.active)
    for var in con.x
        print(io, "\n   ", var)
    end
end

function Base.show(io::IO, con::SumToZero)
    ids = [var.id for var in con.x]
    print(io, typeof(con), ": ", join(ids, " + "), " == 0")
end


text:
using AbstractPlotting:
    NoConversion,
    convert_arguments,
    conversion_trait

using IntervalSets

@testset "Conversions" begin

    # NoConversion
    struct NoConversionTestType end
    conversion_trait(::NoConversionTestType) = NoConversion()

    let nctt = NoConversionTestType(),
        ncttt = conversion_trait(nctt)
        @test convert_arguments(ncttt, 1, 2, 3) == (1, 2, 3)
    end

end

@testset "functions" begin
    x = -pi..pi
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end

    x = range(-pi, stop=pi, length=100)
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end
end

@testset "Categorical values" begin
    # AbstractPlotting.jl#345
    a = Any[Int64(1), Int32(1), Int128(2)] # vector of categorical values of different types
    ilabels = AbstractPlotting.categoric_labels(a)
    @test ilabels == [1, 2]
    @test AbstractPlotting.categoric_position.(a, Ref(ilabels)) == [1, 1, 2]
end


text:
x = rand(AFArray{Float32}, 10, 5)
@test length(x) == 50
@test size(x) == (10,5)
@test size(x,1) == 10
@test size(x,2) == 5
@test eltype(x) == Float32
@test ndims(x) == 5


text:
@doc raw"""
    adjoint_differential_bezier_control(
        M::Manifold,
        b::BezierSegment,
        t::Float64,
        η::Q)

evaluate the adjoint of the differential of a Bézier curve on the manifold `M`
with respect to its control points `b` based on a point `t` $\in[0,1]$ on the
curve and a tangent vector $\eta\in T_{\beta(t)}\mathcal M$.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(
    M::Manifold, b::BezierSegment, t::Float64, η::Q
) where {Q}
    n = length(b.pts)
    if n == 2
        return BezierSegment([
            adjoint_differential_geodesic_startpoint(M, b.pts[1], b.pts[2], t, η),
            adjoint_differential_geodesic_endpoint(M, b.pts[1], b.pts[2], t, η),
        ])
    else
        c = [b.pts, [similar.(b.pts[1:l]) for l in (n - 1):-1:2]...]
        for i in 2:(n - 1) # casteljau on the tree -> forward with interims storage
            c[i] .=
                shortest_geodesic.(Ref(M), c[i - 1][1:(end - 1)], c[i - 1][2:end], Ref(t))
        end
        Y = [η, [similar(η) for i in 1:(n - 1)]...]
        for i in (n - 1):-1:1 # propagate adjoints -> backward without interims storage
            Y[1:(n - i + 1)] .=
                [ # take previous results and add start&end point effects
                    adjoint_differential_geodesic_startpoint.(
                        Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y[1:(n - i)]
                    )...,
                    zero_tangent_vector(M, c[i][end]),
                ] .+ [
                    zero_tangent_vector(M, c[i][1]),
                    adjoint_differential_geodesic_endpoint.(
                        Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y[1:(n - i)]
                    )...,
                ]
        end
    end
    return BezierSegment(Y)
end
@doc raw"""
    adjoint_differential_bezier_control(
        M::Manifold,
        b::BezierSegment,
        t::Array{Float64,1},
        X::Array{Q,1}
    )

evaluate the adjoint of the differential of a Bézier curve on the manifold `M`
with respect to its control points `b` based on a points `T`$=(t_i)_{i=1}^n that
are pointwise in $ t_i\in[0,1]$ on the curve and given corresponding tangential
vectors $X = (\eta_i)_{i=1}^n$, $\eta_i\in T_{\beta(t_i)}\mathcal M$

See [`de_casteljau`](@ref) for more details on the curve and[^BergmannGousenbourger2018].

[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
    > by minimizing the acceleration of a Bézier curve_.
    > Frontiers in Applied Mathematics and Statistics, 2018.
    > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
"""
function adjoint_differential_bezier_control(
    M::Manifold, b::BezierSegment, t::AbstractVector{Float64}, X::AbstractVector{Q}
) where {Q}
    effects = [bt.pts for bt in adjoint_differential_bezier_control.(Ref(M), Ref(b), t, X)]
    return BezierSegment(sum(effects))
end

@doc raw"""
    adjoint_differential_bezier_control(
        M::MAnifold,
        B::AbstractVector{<:BezierSegment},
        t::Float64,
        X
    )

evaluate the adjoint of the differential of a composite Bézier curve on the
manifold `M` with respect to its control points `b` based on a points `T`$=(t_i)_{i=1}^n$
that are pointwise in $t_i\in[0,1]$ on the curve and given corresponding tangential
vectors $X = (\eta_i)_{i=1}^n$, $\eta_i\in T_{\beta(t_i)}\mathcal M$

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(
    M::Manifold, B::AbstractVector{<:BezierSegment}, t::Float64, X::Q
) where {Q}
    # doubly nested broadbast on the Array(Array) of CPs (note broadcast _and_ .)
    if (0 > t) || (t > length(B))
        error(
            "The parameter ",
            t,
            " to evaluate the composite Bézier curve at is outside the interval [0,",
            length(B),
            "].",
        )
    end
    Y = broadcast(b -> BezierSegment(zero_tangent_vector.(Ref(M), b.pts)), B) # Double broadcast
    seg = max(ceil(Int, t), 1)
    localT = ceil(Int, t) == 0 ? 0.0 : t - seg + 1
    Y[seg].pts .= adjoint_differential_bezier_control(M, B[seg], localT, X).pts
    return Y
end
@doc raw"""
    adjoint_differential_bezier_control(M,B,t,η)
evaluate the adjoint of the differential with respect to the controlpoints.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(
    M::Manifold,
    B::AbstractVector{<:BezierSegment},
    T::AbstractVector{Float64},
    X::AbstractVector{Q},
) where {P,Q}
    BT = adjoint_differential_bezier_control.(Ref(M), Ref(B), T, X)
    Y = broadcast(b -> BezierSegment(zero_tangent_vector.(Ref(M), b.pts)), B) # Double broadcast
    for Bn in BT # for all times
        for i in 1:length(Bn)
            Y[i].pts .+= Bn[i].pts
        end
    end
    return Y
end

@doc raw"""
    adjoint_differential_geodesic_startpoint(M,p, q, t, X)

Compute the adjoint of $D_p γ(t; p, q)[X]$.

# See also

[`differential_geodesic_startpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_geodesic_startpoint(M::MT, p, q, t, X) where {MT<:Manifold}
    return adjoint_Jacobi_field(M, p, q, t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    adjoint_differential_geodesic_endpoint(M, p, q, t, X)

Compute the adjoint of $D_q γ(t; p, q)[X]$.

# See also

[`differential_geodesic_endpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_geodesic_endpoint(M::MT, p, q, t, X) where {MT<:Manifold}
    return adjoint_Jacobi_field(M, q, p, 1 - t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    adjoint_differential_exp_basepoint(M, p, X, Y)

Computes the adjoint of $D_p \exp_p X[Y]$.

# See also

[`differential_exp_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_basepoint(M::MT, p, X, Y) where {MT<:Manifold}
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_basepoint)
end

@doc raw"""
    adjoint_differential_exp_argument(M, p, X, Y)

Compute the adjoint of $D_X\exp_p X[Y]$.
Note that $X ∈  T_p(T_p\mathcal M) = T_p\mathcal M$ is still a tangent vector.

# See also

[`differential_exp_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_argument(M::mT, p, X, Y) where {mT<:Manifold}
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end

@doc raw"""
    adjoint_differential_log_basepoint(M, p, q, X)

computes the adjoint of $D_p log_p q[X]$.

# See also
[`differential_log_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_basepoint(M::Manifold, p, q, X)
    return adjoint_Jacobi_field(M, p, q, 0.0, X, βdifferential_log_basepoint)
end

@doc raw"""
    adjoint_differential_log_argument(M, p, q, X)

Compute the adjoint of $D_q log_p q[X]$.

# See also
[`differential_log_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_argument(M::Manifold, p, q, X)
    # order of p and q has to be reversed in this call, cf. Persch, 2018 Lemma 2.3
    return adjoint_Jacobi_field(M, q, p, 1.0, X, βdifferential_log_argument)
end

@doc raw"""
    Y = adjoint_differential_forward_logs(M, p, X)

Compute the adjoint differential of [`forward_logs`](@ref) $F$ orrucirng,
in the power manifold array `p`, the differential of the function

$F_i(p) = \sum_{j ∈ \mathcal I_i} \log_{p_i} p_j$

where $i$ runs over all indices of the `PowerManifold` manifold `M` and $\mathcal I_i$
denotes the forward neighbors of $i$
Let $n$ be the number dimensions of the `PowerManifold` manifold (i.e. `length(size(x)`)).
Then the input tangent vector lies on the manifold $\mathcal M' = \mathcal M^n$.

# Input

* `M`     – a `PowerManifold` manifold
* `p`     – an array of points on a manifold
* `X`     – a tangent vector to from the n-fold power of `p`, where n is the `ndims` of `p`

# Ouput

`Y` – resulting tangent vector in $T_p\mathcal M$ representing the adjoint
  differentials of the logs.
"""
function adjoint_differential_forward_logs(
    M::PowerManifold{𝔽,TM,TSize,TPR}, p, X
) where {𝔽,TM,TSize,TPR}
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R).I
    N = PowerManifold(M.manifold, TPR(), power_size..., d)
    Y = zero_tangent_vector(M, p)
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = [i.I...] # array of index
            J = I .+ 1 .* (1:d .== k) #i + e_k is j
            if all(J .<= maxInd) # is this neighbor in range?
                j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
                Y[M, I...] =
                    Y[M, I...] + adjoint_differential_log_basepoint(
                        M.manifold, p[M, I...], p[M, J...], X[N, I..., k]
                    )
                Y[M, J...] =
                    Y[M, J...] + adjoint_differential_log_argument(
                        M.manifold, p[M, J...], p[M, I...], X[N, I..., k]
                    )
            end
        end # directions
    end # i in R
    return Y
end