text:
_maybepromote(::Type{T}) where {T} = T
_maybepromote(::Type{T}) where {T<:AbstractFloat} = Float32

function innerproduct(
    _src::Memory{T},
    _weights::Memory{T},
    bias::Memory{T};
    kind = Inference(),
    attributes::Attributes = noattributes(),
    # Allow for callers to maybe convert these descriptions to "any", allowing callers
    # to later convert their data to a better format.
    src_md = memorydesc_ptr(_src),
    weights_md = memorydesc_ptr(_weights),
    callback = (x...) -> nothing,
) where {T}
    #@show (size(_src), size(_weights), size(bias))
    dst_dims = (size(bias, 1), size(_src, 2))
    dst_md = memorydesc(T, dst_dims, dnnl_format_any())
    opdesc = Ref{Lib.dnnl_inner_product_desc_t}()

    @apicall dnnl_inner_product_forward_desc_init(
        opdesc, kind, src_md, weights_md, bias, dst_md
    )

    return temp_primitive(opdesc, attributes, global_engine(), noforward()) do p, pd
        # Maybe the caller wants to know what the optimal weights format is.
        callback(p, pd)

        # Convert if needed.
        src = reorder_if_any(src_md, pd, _src, @query(src))
        weights = reorder_if_any(weights_md, pd, _weights, @query(weights))

        # Materialize destinations
        dst_md = query_md(pd, @query(dst))
        dst = similar(src, T, dst_dims, dst_md)

        # Execute primitive
        execute!(p, @dnnl_args src weights bias dst; descriptor = pd)
        return dst
    end
end

function innerproduct_backward_data(
    diff_src_dims::NTuple, _weights::Memory, _diff_dst::Memory{T}
) where {T}
    diff_src_md = memorydesc(T, diff_src_dims, dnnl_format_any())
    opdesc = Ref{Lib.dnnl_inner_product_desc_t}()
    @apicall dnnl_inner_product_backward_data_desc_init(
        opdesc, diff_src_md, toany(_weights), toany(_diff_dst)
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, pd
        weights = maybe_reorder(pd, _weights, @query(weights))
        diff_dst = maybe_reorder(pd, _diff_dst, @query(diff_dst))

        diff_src_md = query_md(pd, @query(diff_src))
        diff_src = similar(diff_dst, T, diff_src_dims, diff_src_md)

        execute!(p, @dnnl_args diff_src weights diff_dst)
        return diff_src
    end
end

function innerproduct_backward_weights(
    diff_weights_dims::NTuple, _src::Memory{T}, _diff_dst::Memory
) where {T}
    diff_bias_dims = (diff_weights_dims[2],)
    diff_weights_md = memorydesc(_maybepromote(T), diff_weights_dims, dnnl_format_any())
    diff_bias_md = memorydesc(_maybepromote(T), diff_bias_dims)
    opdesc = Ref{Lib.dnnl_inner_product_desc_t}()
    @apicall dnnl_inner_product_backward_weights_desc_init(
        opdesc, toany(_src), diff_weights_md, diff_bias_md, toany(_diff_dst)
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, pd
        #@show (diff_weights_dims, size(_src), size(_diff_dst))
        # Maybe convert arguments
        src = maybe_reorder(pd, _src, @query(src))
        diff_dst = maybe_reorder(pd, _diff_dst, @query(diff_dst))

        # Materialize destinations
        diff_weights_md = query_md(pd, @query(diff_weights))
        diff_weights = similar(
            diff_dst, _maybepromote(T), diff_weights_dims, diff_weights_md
        )
        diff_bias = similar(diff_dst, _maybepromote(T), diff_bias_dims, diff_bias_md)

        # Execute and return
        execute!(p, @dnnl_args diff_weights diff_bias src diff_dst)
        return (diff_weights, diff_bias)
    end
end

#####
##### Dense Layer
#####

mutable struct Dense{W<:Memory,B<:Memory,F}
    weights::W
    bias::B
    activation::F
    attributes::Attributes
    # Memory format optimization
    optimized_weights::Bool
end

function Dense(
    weights::AbstractMatrix{T},
    bias::AbstractVector{T},
    activation::F = identity,
) where {T,F}
    weights_memory = Memory(weights)
    attributes = Attributes()
    postops = PostOps()
    eltwise!(postops, activation)
    append!(attributes, postops)

    return Dense(
        Memory(weights), Memory(bias), activation, attributes, false
    )
end
Dense(m::Flux.Dense) = Dense(OneDNN.Memory(transpose(m.weight)), OneDNN.Memory(m.bias), m.σ)

Flux.@functor Dense (weights, bias)
function ZygoteRules._pullback(
    cx::AContext, ::typeof(literal_getproperty), x::Dense, ::Val{f}
) where {f}
    return pullback_for_default_literal_getproperty(cx, x, Val{f}())
end

function (dense::Dense{<:OneDNN.Memory{T}})(_src, fuse_activation = true) where {T}
    src = toeltype(T, Memory(_src))
    attributes = fuse_activation ? dense.attributes : noattributes()

    if dense.optimized_weights == false
        optimized_weights_desc = Ref{MemoryDesc}()
        function weight_opt_callback(_, pd)
            optimized_weights_desc[] = query_md(pd, Lib.dnnl_query_weights_md)
            return nothing
        end

        dst = innerproduct(
            src,
            dense.weights,
            dense.bias;
            attributes = attributes,
            src_md = toany(src),
            weights_md = toany(dense.weights),
            callback = weight_opt_callback,
        )

        dense.weights = maybe_reorder(optimized_weights_desc[], dense.weights)
        dense.optimized_weights = true
        return dst
    else
        return innerproduct(
            src,
            dense.weights,
            dense.bias;
            attributes = attributes,
            src_md = toany(src),
        )
    end
end

function ChainRulesCore.rrule(
    dense::Dense{<:OneDNN.Memory{T}}, _src::AbstractMatrix{U}, fuse_activation::Bool
) where {T,U}
    src = toeltype(T, Memory(_src))
    # The result of `canfuse` is known at compile time, so Julia can optimize out the branch.
    return if canfuse(dense.activation)
        rrule_fused(dense, src, true)
    else
        rrule_unfused(dense, src, false)
    end
end

function rrule_fused(dense::Dense{<:OneDNN.Memory{T}}, _src, _fuse_activation) where {T}
    src = Memory(_src)
    src_size = size(src)
    dst = dense(src, true)
    function dense_fused_pullback(_diff_dst)
        # Maybe convert argument
        diff_dst = toeltype(eltype(dst), Memory(_diff_dst))
        # Reverse activation function.
        diff_dst_pre = toeltype(T, eltwise_backward(dense.activation, diff_dst, dst))

        # Backprop innerproduct kernel
        diff_src = innerproduct_backward_data(src_size, dense.weights, diff_dst_pre)
        (diff_weights, diff_bias) = innerproduct_backward_weights(
            size(dense.weights), src, diff_dst_pre
        )
        return (
            ChainRulesCore.Tangent{typeof(dense)}(;
                weights = diff_weights, bias = diff_bias
            ),
            diff_src,
            ChainRulesCore.NoTangent(),
        )
    end
    return dst, dense_fused_pullback
end

function rrule_unfused(dense::T, _src, _fuse_activation) where {T<:Dense}
    src = Memory(_src)
    dst_pre = dense(src, false)
    dst = eltwise(dense.activation, dst_pre)
    src_size = size(src)
    function dense_pullback(_diff_dst)
        diff_dst = Memory(_diff_dst)
        diff_dst_pre = eltwise_backward(dense.activation, diff_dst, dst_pre)
        diff_src = innerproduct_backward_data(src_size, dense.weights, diff_dst_pre)
        (diff_weights, diff_bias) = innerproduct_backward_weights(
            size(dense.weights), src, diff_dst_pre
        )

        return (
            ChainRulesCore.Tangent{T}(; weights = diff_weights, bias = diff_bias),
            diff_src,
            ChainRulesCore.NoTangent(),
        )
    end
    return dst, dense_pullback
end


text:
import  Base: zero, one, zeros, start, next, done, length, eltype,
        size, setindex!, getindex, similar, fill!, count #, print # , copy, deepcopy, fill!, getindex, setindex!

"""
`SequenceWeights` is an alias for `Union{ClusteringResult, AbstractVector}`.

This is the type of the keyword argment `weight` of `count`, `probabilities` and derived functions.
The type should define `getweight` to be useful in those functions.
"""
typealias SequenceWeights Union{ClusteringResult, AbstractVector} # AbstractVector because getweight(cl) returns a Vector

# Counts and Pseudocounts

abstract Pseudocount{T<:Real}

"""
**Additive Smoothing** or fixed pseudocount  `λ`  for `ResidueCount` (in order to estimate probabilities when the number of samples is low).

Common values of `λ` are:

- `0` :  No cell frequency prior, gives you the maximum likelihood estimator.
- `0.05` is the optimum value for `λ` found in Buslje et. al. 2009, similar results was obtained for `λ` in the range [0.025, 0.075].
- `1 / p` : Perks prior (Perks, 1947) where `p` the number of parameters (i.e. residues, pairs of residues) to estimate. If `p` is the number of residues (`20` without counting gaps), this gives you `0.05`.
- `sqrt(n) / p` : Minimax prior (Trybula, 1958) where `n` is the number of samples and `p` the number of parameters to estimate.  If the number of samples `n` is 400 (minimum number of sequence clusters for achieve good performance in Buslje et. al. 2009) for estimating 400 parameters (pairs of residues without counting gaps) this gives you `0.05`.
- `0.5` : Jeffreys prior (Jeffreys, 1946).
- `1` : Bayes-Laplace uniform prior, aka. Laplace smoothing.
"""
immutable AdditiveSmoothing{T} <: Pseudocount{T}
  λ::T
end

zero{T}(::Type{AdditiveSmoothing{T}}) = AdditiveSmoothing{T}(zero(T))
one{T}(::Type{AdditiveSmoothing{T}}) = AdditiveSmoothing{T}(one(T))

## ResidueContingencyTables

abstract ResidueContingencyTables{T, N, UseGap} <: AbstractArray{T, N}

### Abstract Array Interface

Base.linearindexing(::ResidueContingencyTables)	= Base.LinearFast()

"""
`getindex(n::ResidueContingencyTables, i::Int)`

Gives access to the table. Use `Int()` for indexing with `Residue`s.
"""
getindex(n::ResidueContingencyTables, i::Int)	= getindex(n.table, i)

"""
`setindex!(n::ResidueCount, v, i::Int)` set a value into the table, but doesn't update the marginals (and total).
Use `Int()` for indexing with `Residues`. Use `update!` in order to calculate again the marginals (and total).
"""
setindex!(n::ResidueContingencyTables, v, i::Int) = setindex!(n.table, v, i)

#### Length & Size

length{T}(n::ResidueContingencyTables{T, 1, true}) = 21
length{T}(n::ResidueContingencyTables{T, 1,false}) = 20

length{T}(n::ResidueContingencyTables{T, 2, true}) = 441
length{T}(n::ResidueContingencyTables{T, 2,false}) = 400

length{T,N,UseGap}(n::ResidueContingencyTables{T, N, UseGap}) = length(n.table)

size{T}(n::ResidueContingencyTables{T, 1, true}) = (21,)
size{T}(n::ResidueContingencyTables{T, 1,false}) = (20,)

size{T}(n::ResidueContingencyTables{T, 2, true}) = (21, 21)
size{T}(n::ResidueContingencyTables{T, 2,false}) = (20, 20)

size{T,N,UseGap}(n::ResidueContingencyTables{T, N, UseGap}) = size(n.table)

"""
Number of residues used in the `ResidueContingencyTables`.
A 20x20 table returns 20 (ungapped alphabet).
"""
@inline nresidues{T, N}(n::ResidueContingencyTables{T, N, true})  = 21
@inline nresidues{T, N}(n::ResidueContingencyTables{T, N, false}) = 20

#### Iteration Interface

start(n::ResidueContingencyTables) = 1

@inbounds next(n::ResidueContingencyTables, state::Int) = (n.table[state], state + 1)

done{T,N,UseGap}(n::ResidueContingencyTables{T, N, UseGap}, state) = state > length(n)

eltype{T,N,UseGap}(::Type{ResidueContingencyTables{T, N, UseGap}}) = T

## Counts

"""
`ResidueCount{T, N, UseGap}` is used for counting residues in columns (or sequences) of an MSA.
`N` is the dimensionality and should be an `Int`, i.e. 2 if 2 columns are used for counting pairs.
`UseGap` is a `Bool`, `true` means that **ResidueCount** counts gaps in the position 21.

- The field marginal is used for pre allocation of marginal sums.
- The field total is used for storing the table sum.
"""
type ResidueCount{T, N, UseGap} <: ResidueContingencyTables{T, N, UseGap}
  table::Array{T, N}
  marginals::Array{T, 2}
  total::T
end

call{T}(::Type{ResidueCount{T, 1, true}}) = ResidueCount{T, 1, true}(Array(T, 21), Array(T, (21,1)), zero(T))
call{T}(::Type{ResidueCount{T, 2, true}}) = ResidueCount{T, 2, true}(Array(T, (21, 21)), Array(T, (21,2)), zero(T))

call{T}(::Type{ResidueCount{T, 1, false}}) = ResidueCount{T, 1, false}(Array(T, 20), Array(T, (20,1)), zero(T))
call{T}(::Type{ResidueCount{T, 2, false}}) = ResidueCount{T, 2, false}(Array(T, (20, 20)), Array(T, (20,2)), zero(T))

function call{T, N, UseGap}(::Type{ResidueCount{T, N, UseGap}})
  nres = UseGap ? 21 : 20
  ResidueCount{T, N, UseGap}(Array(T, (Int[ nres for d in 1:N]...)), Array(T, (nres, N)), zero(T))
end

zeros{T}(::Type{ResidueCount{T, 1, true}}) = ResidueCount{T, 1, true}(zeros(T, 21), zeros(T, (21, 1)), zero(T))
zeros{T}(::Type{ResidueCount{T, 2, true}}) = ResidueCount{T, 2, true}(zeros(T, (21, 21)), zeros(T, (21, 2)), zero(T))

zeros{T}(::Type{ResidueCount{T, 1, false}}) = ResidueCount{T, 1, false}(zeros(T, 20), zeros(T, (20, 1)), zero(T))
zeros{T}(::Type{ResidueCount{T, 2, false}}) = ResidueCount{T, 2, false}(zeros(T, (20, 20)), zeros(T, (20, 2)), zero(T))

function zeros{T, N, UseGap}(::Type{ResidueCount{T, N, UseGap}})
  nres = UseGap ? 21 : 20
  ResidueCount{T, N, UseGap}(zeros(T, (Int[ nres for d in 1:N]...)), zeros(T, (nres, N)), zero(T))
end

#### Similar

similar{T,N,UseGap}(n::ResidueCount{T,N,UseGap}) = ResidueCount{T,N,UseGap}()
similar{T,N,S,UseGap}(n::ResidueCount{T,N,UseGap}, ::Type{S}) = ResidueCount{S,N,UseGap}()
similar{T,N,UseGap}(n::ResidueCount{T,N,UseGap}, D::Int) = ResidueCount{T,D,UseGap}()
similar{T,N,S,UseGap}(n::ResidueCount{T,N,UseGap}, ::Type{S}, D::Int) = ResidueCount{S,D,UseGap}()

### Update!

function _tuple_without_index(index::Int, N::Int)
	used = Array(Int, N - 1)
	j = 1
	for i in 1:N
		if i != index
			used[j] = i
			j += 1
		end
	end
	(used...)
end

"Updates the marginals values (and total value for `ResidueCount`) of a `ResidueContingencyTables`"
function update!{T,UseGap}(n::ResidueCount{T,1,UseGap})
	n.marginals[:] = n.table
	n.total = sum(n.table)
	n
end

function update!{T,UseGap}(n::ResidueCount{T,2,UseGap})
	n.marginals[:,1] = sum(n.table, 2) # this is faster than sum(n,2)
	n.marginals[:,2] = sum(n.table, 1)
	n.total = sum(n.marginals[:,1])
	n
end

function update!{T, N, UseGap}(n::ResidueCount{T, N, UseGap})
	@inbounds for i in 1:N
		n.marginals[:,i] = sum(n.table, _tuple_without_index(i,N))
	end
	n.total = sum(n.marginals[:,1])
	n
end

### Apply Pseudocount

# This is faster than array[:] += value
function _sum!(array, value)
  @inbounds for i in eachindex(array)
    array[i] += value
  end
  array
end

for (dim, gap, margi_exp, total_exp) in [ (:1, :true,  :(pse.λ), :(pse.λ * 21)),
																					(:1, :false, :(pse.λ), :(pse.λ * 20)),
																					(:2, :true,  :(pse.λ * 21), :(pse.λ * 441)),
																					(:2, :false, :(pse.λ * 20), :(pse.λ * 400)) ]
	@eval begin

		function apply_pseudocount!{T}(n::ResidueCount{T, $(dim), $(gap)}, pse::AdditiveSmoothing{T})
			margi_sum = $(margi_exp)
			total_sum = $(total_exp)
			_sum!(n.table, pse.λ)
			_sum!(n.marginals, margi_sum)
			n.total += total_sum
			n
		end

		function fill!{T}(n::ResidueCount{T, $(dim), $(gap)}, pse::AdditiveSmoothing{T})
			margi_sum = $(margi_exp)
			total_sum = $(total_exp)
			fill!(n.table, pse.λ)
			fill!(n.marginals, margi_sum)
			n.total = total_sum
			n
		end

	end

end

"""
`apply_pseudocount!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, pse::AdditiveSmoothing{T})`

Uses an instance of `AdditiveSmoothing` to efficiently fill with a constant value each element of the table.
"""
function apply_pseudocount!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, pse::AdditiveSmoothing{T})
	nres = nresidues(n)
	margi_sum = pse.λ * (nres^(N-1))
	total_sum = pse.λ * (nres^N)
	_sum!(n.table, pse.λ)
	_sum!(n.marginals, margi_sum)
	n.total += total_sum
	n
end

"""
`fill!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, pse::AdditiveSmoothing{T})`
fills a preallocated `ResidueCount` (`p`) with the pseudocount (`pse`).
"""
function fill!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, pse::AdditiveSmoothing{T})
	nres = nresidues(n)
	margi_sum = pse.λ * (nres^(N-1))
	total_sum = pse.λ * (nres^N)
	fill!(n.table, pse.λ)
	fill!(n.marginals, margi_sum)
	n.total = total_sum
	n
end

### Counting

function count!{T}(n::ResidueCount{T, 1, true}, weights, res::AbstractVector{Residue})
    @simd for i in 1:length(res)
        @inbounds n.table[Int(res[i])] += getweight(weights,i)
    end
    update!(n)
end

function count!{T}(n::ResidueCount{T, 1, false}, weights, res::AbstractVector{Residue})
    @inbounds for i in 1:length(res)
        aa = Int(res[i])
        aa == Int(GAP) && continue
        n.table[aa] += getweight(weights,i)
    end
    update!(n)
end

count!{T, UseGap}(n::ResidueCount{T, 1, UseGap}, res::AbstractVector{Residue}) = count!(n, NoClustering(), res)

function count!{T}(n::ResidueCount{T, 2, true}, weights, res1::AbstractVector{Residue}, res2::AbstractVector{Residue})
  @assert length(res1) == length(res2)
  @simd for i in 1:length(res1)
    @inbounds n.table[Int(res1[i]), Int(res2[i])] += getweight(weights,i)
  end
  update!(n)
end

function count!{T}(n::ResidueCount{T, 2, false}, weights, res1::AbstractVector{Residue}, res2::AbstractVector{Residue})
  for i in 1:length(res1)
    aa1 = res1[i]
    aa2 = res2[i]
    if (aa1 != GAP) && (aa2 != GAP)
      n.table[Int(aa1), Int(aa2)] += getweight(weights,i)
    end
  end
  update!(n)
end

count!{T}(n::ResidueCount{T, 2, true}, res1::AbstractVector{Residue}, res2::AbstractVector{Residue}) = count!(n, NoClustering(), res1, res2)
count!{T}(n::ResidueCount{T, 2, false},res1::AbstractVector{Residue}, res2::AbstractVector{Residue}) = count!(n, NoClustering(), res1, res2)

function count!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, weights, res::AbstractVector{Residue}...)
  if length(res) == N
    for i in 1:length(res[1])
      aa_list = [ Int(aa[i]) for aa in res ]
      if UseGap || (findfirst(aa_list, Int(GAP)) == 0)
        n.table[aa_list...] += getweight(weights,i)
      end
    end
    update!(n)
  else
    throw("Number of arrays ($(length(res))) != $N")
  end
end

"""
`count!` adds counts from vector of residues to a `ResidueCount` object.
It can take a SequenceWeights object as second argument.

```julia

julia> using MIToS.Information

julia> using MIToS.MSA

julia> seq = Residue[ i for i in 1:20];

julia> Ni = count(seq)
20-element MIToS.Information.ResidueCount{Float64,1,false}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

julia> count!(Ni, seq)
20-element MIToS.Information.ResidueCount{Float64,1,false}:
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0

```
"""
count!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, res::AbstractVector{Residue}...) = count!(n, NoClustering(), res...)

#### Default counters

"""
```
count(res::AbstractVector{Residue}...; usegap=false, weight=NoClustering())
count(pseudocount::Pseudocount, res::AbstractVector{Residue}...; usegap=false, weight=NoClustering())
```

`count` creates a new ResidueCount counting the number of residues, pairs of residues, etc. in the sequences/columns.
"""
count(res::AbstractVector{Residue}...; usegap::Bool=false,
	weight::SequenceWeights=NoClustering()) = count!(zeros(ResidueCount{Float64, length(res), usegap}), weight, res...)

count(pseudocount::Pseudocount, res::AbstractVector{Residue}...; usegap::Bool=false,
	weight::SequenceWeights=NoClustering()) = apply_pseudocount!(count!(zeros(ResidueCount{Float64, length(res), usegap}), weight, res...), pseudocount)

count{T}(pseudocount::AdditiveSmoothing{T}, res::AbstractVector{Residue}...; usegap::Bool=false,
	weight::SequenceWeights=NoClustering()) = apply_pseudocount!(count!(zeros(ResidueCount{T, length(res), usegap}), weight, res...), pseudocount)

## Probabilities

"""
`ResidueProbability{T, N, UseGap}` is used to store residue probabilities.
`N` is the dimensionality and should be an `Int`, e.g. 2 to store probabilities of residue pairs.
`UseGap` is a `Bool`, `true` means that gap probabilities are stored in the position 21 of each dimension.

- The field marginal is used for pre allocation of marginal sums.
"""
type ResidueProbability{T, N, UseGap} <: ResidueContingencyTables{T, N, UseGap}
  table::Array{T, N}
  marginals::Array{T, 2}
end

call{T}(::Type{ResidueProbability{T, 1, true}}) = ResidueProbability{T, 1, true}(Array(T, 21), Array(T, (21,1)))
call{T}(::Type{ResidueProbability{T, 2, true}}) = ResidueProbability{T, 2, true}(Array(T, (21, 21)), Array(T, (21,2)))

call{T}(::Type{ResidueProbability{T, 1, false}}) = ResidueProbability{T, 1, false}(Array(T, 20), Array(T, (20,1)))
call{T}(::Type{ResidueProbability{T, 2, false}}) = ResidueProbability{T, 2, false}(Array(T, (20, 20)), Array(T, (20,2)))

function call{T, N, UseGap}(::Type{ResidueProbability{T, N, UseGap}})
  nres = UseGap ? 21 : 20
  ResidueProbability{T, N, UseGap}(Array(T, (Int[ nres for d in 1:N]...)), Array(T, (nres, N)))
end

zeros{T}(::Type{ResidueProbability{T, 1, true}}) = ResidueProbability{T, 1, true}(zeros(T, 21), zeros(T, (21, 1)))
zeros{T}(::Type{ResidueProbability{T, 2, true}}) = ResidueProbability{T, 2, true}(zeros(T, (21, 21)), zeros(T, (21, 2)))

zeros{T}(::Type{ResidueProbability{T, 1, false}}) = ResidueProbability{T, 1, false}(zeros(T, 20), zeros(T, (20, 1)))
zeros{T}(::Type{ResidueProbability{T, 2, false}}) = ResidueProbability{T, 2, false}(zeros(T, (20, 20)), zeros(T, (20, 2)))

function zeros{T, N, UseGap}(::Type{ResidueProbability{T, N, UseGap}})
  nres = UseGap ? 21 : 20
  ResidueProbability{T, N, UseGap}(zeros(T, (Int[ nres for d in 1:N]...)), zeros(T, (nres, N)))
end

#### Similar

similar{T, N,UseGap}(p::ResidueProbability{T, N,UseGap}) = ResidueProbability{T, N,UseGap}()
similar{T, N,UseGap}(p::ResidueProbability{T, N,UseGap}, D::Int) = ResidueProbability{T, D,UseGap}()

### Update!

function update!{T, UseGap}(p::ResidueProbability{T, 1,UseGap})
	p.marginals[:] = p.table
	p
end

function update!{T, UseGap}(p::ResidueProbability{T, 2,UseGap})
	p.marginals[:,1] = sum(p.table, 2) # this is faster than sum(p,2)
	p.marginals[:,2] = sum(p.table, 1)
	p
end

function update!{T, N, UseGap}(p::ResidueProbability{T, N, UseGap})
	for i in 1:N
		p.marginals[:,i] = sum(p.table, _tuple_without_index(i,N))
	end
	p
end

### Normalize

# This is faster than array[:] /= value
function _div!(array, value)
  @inbounds for i in eachindex(array)
    array[i] /= value
  end
  array
end

"""
`normalize!(p::ResidueProbability; updated::Bool=false)`

This function makes the sum of the probabilities to be one.
The sum is calculated using the `probabilities` field by default (It is assumed that the marginal are not updated).
The marginals are updated in the normalization.

If the marginals are updated, you can use `updated=true` for a faster normalization.
"""
function normalize!{T, N, UseGap}(p::ResidueProbability{T, N, UseGap}; updated::Bool=false)
  if !updated
    update!(p)
  end
  total = sum(p.marginals[:,1])
  if total != T(1.0)
    _div!(p.table, total)
    _div!(p.marginals, total)
  end
  p
end

### ResidueProbability calculation from ResidueCount and others

# This is faster than p[:] = ( n ./ total )
function _fill_probabilities!{TP, TN, N}(p::Array{TP,N}, n::Array{TN,N}, total::TP)
  @inbounds for i in 1:length(p) # p and n should have the same length
    p[i] = ( n[i] / total)
  end
  p
end

"""
`fill!{T, N, UseGap}(p::ResidueProbability{T, N, UseGap}, n::ResidueCount{T, N, UseGap}; updated::Bool=false)`

This function fills a preallocated `ResidueProbability` (`p`) with the probabilities calculated from `n` (`ResidueCount`).
This function updates `n` unless `updated=true`.

If `n` is updated, you can use `updated=true` for a faster calculation.
"""
function fill!{TP, TC, N, UseGap}(p::ResidueProbability{TP, N, UseGap}, n::ResidueCount{TC, N, UseGap}; updated::Bool=false)
  if !updated
    update!(n)
  end
  total = TP(n.total)
  _fill_probabilities!(p.table, n.table, total)
  _fill_probabilities!(p.marginals, n.marginals, total)
	p
end

### BLOSUM based pseudofrequencies

"""
`blosum_pseudofrequencies!(Gab::ResidueProbability{T, 2,false}, Pab::ResidueProbability{T, 2,false})`

This function uses the conditional probability matrix `BLOSUM62_Pij` to fill a preallocated `Gab` with pseudo frequencies.
`blosum_pseudofrequencies!` also needs the real frequencies/probabilities `Pab`.
This observed probabilities are then used to estimate the pseudo frequencies.

`Gab = Σcd  Pcd ⋅ BLOSUM62( a | c ) ⋅ BLOSUM62( b | d )`
"""
function blosum_pseudofrequencies!{T}(Gab::ResidueProbability{T, 2,false}, Pab::ResidueProbability{T, 2,false})
  @inbounds for a in 1:20, b in 1:20
    Gab[a,b] = zero(T)
    for i in 1:20, j in 1:20
      P = Pab[i,j]
      if P != 0
	      # BLOSUM62_P_i_j[i,a] is p(a | i)
	      Gab[a,b] += ( P * BLOSUM62_Pij[i,a] * BLOSUM62_Pij[j,b] )
      end
    end
  end
  normalize!(Gab)
end

### Apply pseudofrequencies

"""
`apply_pseudofrequencies!{T}(Pab::ResidueProbability{T, 2,false}, Gab::ResidueProbability{T, 2,false}, α, β)`

Apply pseudofrequencies `Gab` over `Pab`, as weighted mean.
Where α is the weight of the real frequencies `Pab` and β the weight of the pseudofrequencies.

`Pab = (α ⋅ Pab + β ⋅ Gab )/(α + β)`
"""
function apply_pseudofrequencies!{T}(Pab::ResidueProbability{T, 2,false}, Gab::ResidueProbability{T, 2,false}, α, β)
	if β == 0.0
		return(Pab)
	end
  frac = one(T) / ( α + β )
  total = zero(T)
  @inbounds for i in 1:20, j in 1:20
    value = ( α * Pab[i,j] + β * Gab[i,j] ) * frac
    Pab[i,j] = value
    total += value
  end
  if total != one(T)
    _div!(Pab.table, total)
    update!(Pab)
  else
    update!(Pab)
  end
end

## Default probabilities

"""
```
probabilities(res::AbstractVector{Residue}...; usegap=false, weight=NoClustering())
probabilities(pseudocount::Pseudocount, res::AbstractVector{Residue}...; usegap=false, weight=NoClustering())
```

`probabilities` creates a new ResidueProbability with the probabilities of residues, pairs of residues, etc. in the sequences/columns.
"""
probabilities{T}(::Type{T}, res::AbstractVector{Residue}...; usegap::Bool=false,
	weight::SequenceWeights=NoClustering()) = fill!(ResidueProbability{T, length(res), usegap}(), count(res..., usegap=usegap, weight=weight))

probabilities{T}(::Type{T}, pseudocount::Pseudocount, res::AbstractVector{Residue}...; usegap::Bool=false,
	weight::SequenceWeights=NoClustering()) = fill!(ResidueProbability{T, length(res), usegap}(), count(pseudocount, res..., usegap=usegap, weight=weight))

"""
`probabilities(T, α, β, res1, res2, [weight])` use BLOSUM62 based pseudofrequencies.
α is the weight of the evidence, and β the weight of the pseudofrequencies.
"""
function probabilities{T}(::Type{T}, α, β, res1::AbstractVector{Residue}, res2::AbstractVector{Residue}; weight::SequenceWeights=NoClustering())
	Pab = fill!(ResidueProbability{T, 2, false}(), count(res1, res2, usegap=false, weight=weight))
	Gab = blosum_pseudofrequencies!(ResidueProbability{T, 2,false}(), Pab)
	apply_pseudofrequencies!(Pab, Gab, α, β)
end

#### Delete dimensions
#### Useful for get Nxy from Nxyz

function _list_without_dimensions(len::Int, out_len::Int, dimensions::Int...)
  ndim = length(dimensions)
  (len - ndim) != out_len && throw("$out_len should be $len minus the number of dimensions = $(len - ndim)")
  index_list = Array(Int, out_len)
  j = 1
  @inbounds for i in 1:len
    if ! (i in dimensions)
      index_list[j] = i
      j += 1
    end
  end
  index_list
end

"""
`delete_dimensions!(out::ResidueContingencyTables, in::ResidueContingencyTables, dimensions::Int...)`

This function fills a ResidueContingencyTables with the counts/probabilities on `in` after the deletion of `dimensions`.
i.e. This is useful for getting Pxy from Pxyz.
"""
function delete_dimensions!{T,N,S,UseGap}(out::ResidueCount{T, S, UseGap}, in::ResidueCount{T, N, UseGap}, dimensions::Int...)
  out.total = in.total
  out.marginals[:] = in.marginals[:, _list_without_dimensions(N, S, dimensions...)]
  out.table[:] = sum(in.table, dimensions)
  out
end

function delete_dimensions!{T,N,S,UseGap}(out::ResidueProbability{T, S, UseGap}, in::ResidueProbability{T, N, UseGap}, dimensions::Int...)
  out.marginals[:] = in.marginals[:, _list_without_dimensions(N, S, dimensions...)]
  out.table[:] = sum(in.table, dimensions)
  out
end

"""
`delete_dimensions(in::ResidueContingencyTables, dimensions::Int...)`

This function creates a ResidueContingencyTables with the counts/probabilities on `in` after the deletion of `dimensions`.
i.e. This is useful for getting Pxy from Pxyz.
"""
delete_dimensions{T,N,UseGap}(in::ResidueCount{T, N, UseGap}, dimensions::Int...) = delete_dimensions!(ResidueCount{T, N-length(dimensions), UseGap}(), in, dimensions...)
delete_dimensions{T,N,UseGap}(in::ResidueProbability{T, N, UseGap}, dimensions::Int...) = delete_dimensions!(ResidueProbability{T, N-length(dimensions), UseGap}(), in, dimensions...)

ext: jl


text:
## Sensor model that returns the exact state of the other vehicles

"""
Sensor that returns the exact state of the other vehicles
"""
struct PerfectSensor <: AbstractSensor end

function measure(sensor::PerfectSensor, ego::Vehicle, scene::Scene, roadway::Roadway, obstacles::Vector{ConvexPolygon})
    obs = Vector{Vehicle}()
    sizehint!(obs, 10) #XXX
    for veh in scene
        if veh.id == ego.id
            continue
        end
        push!(obs, veh)
    end
    return obs
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    modul Simulate_BouncingBall

Simulate bouncing ball (model with state events and discontinuous change of DAE states).
"""
module Simulate_BouncingBall

include(joinpath("models", "BouncingBall.jl"))
import .BouncingBall
import ModiaMath

# Simulate
model  = BouncingBall.Model()
result = ModiaMath.simulate!(model, stopTime=3.0, log=true)

ModiaMath.plot(result, [:h, (:v, :flying)], heading="Simulate_BouncingBall.jl")

end


text:
# Check that serializer hasn't gone out-of-frame
@test Base.ser_tag[Symbol] == 2
@test Base.ser_tag[()] == 47
@test Base.ser_tag[false] == 123

# issue #1770
let
    a = ['T', 'e', 's', 't']
    f = IOBuffer()
    serialize(f, a)
    seek(f, 0)
    @test deserialize(f) == a
    f = IOBuffer()
    serialize(f, a)
    seek(f, 0)
    @test deserialize(f) == a

    # issue #4414
    seek(f,0)
    serialize(f, :β)
    seek(f,0)
    @test deserialize(f) === :β
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""
The abstract base for types that manage the objective function evaluation.
`P` is the optimization problem it is used for.
"""
abstract type Evaluator{P <: OptimizationProblem} end

fitness_scheme(e::Evaluator) = fitness_scheme(problem(e))
fitness_type(e::Evaluator) = fitness_type(fitness_scheme(e))
worst_fitness(e::Evaluator) = worst_fitness(fitness_scheme(e))
numdims(e::Evaluator) = numdims(problem(e))
search_space(e::Evaluator) = search_space(problem(e))
describe(e::Evaluator) = "Problem: $(name(e.problem)) (dimensions = $(numdims(e)))"
problem_summary(e::Evaluator) = "$(name(e.problem))_$(numdims(e))d"

shutdown!(e::Evaluator) = e # do nothing

"""
    update_fitness!([f], eval::Evaluator, candidate::Candidate; force::Bool=false) -> Candidate

Calculate fitness of `candidate` and optionally apply `f`.
`force` specifies whether to re-evaluate fitness, if the candidate already has non-NA one.
"""
function update_fitness!(f::Any, e::Evaluator, candidate::Candidate; force::Bool=false)
    # evaluate fitness if not known yet
    if force || isnafitness(candidate.fitness, fitness_scheme(e.archive))
        candidate.fitness = fitness(candidate.params, e, candidate.tag)
        (f !== nothing) && f(candidate)
    end
    return candidate
end

update_fitness!(e::Evaluator, candidate::Candidate; force::Bool=false) =
    update_fitness!(nothing, e, candidate, force=force)

"""
    update_fitness!([f], eval::Evaluator, candidates; force::Bool=false)

Calculate fitness of `candidates` and optionally apply `f` to each processed one.
`force` specifies if already existing non-NA fitnesses should be re-evaluated.
"""
function update_fitness!(f::Any, e::Evaluator, candidates::Any; force::Bool=false)
    fs = fitness_scheme(e.archive)
    for candi in candidates
        if force || isnafitness(fitness(candi), fs)
            update_fitness!(f, e, candi)
        end
    end
    return candidates
end

update_fitness!(e::Evaluator, candidates::Any; force::Bool=false) =
    update_fitness!(nothing, e, candidates, force=force)

"""
The abstract base type for asynchronous evaluators.

Asynchronous evaluator allows submitting fitness calculation requests via
[`async_update_fitness!`](@ref) method that returns
[`AbstractFitnessEvaluationJob`](@ref) object without waiting for
the calculation completion. Job object could then be passed to
[`sync_update_fitness`](@ref) that waits until fitnesses for all specified jobs
is calculated.
"""
abstract type AbstractAsyncEvaluator{P <: OptimizationProblem} <: Evaluator{P} end

"""
The abstract base type for the fitness evaluation jobs returned by
[`async_update_fitness!`](@ref) and accepted by [`sync_update_fitness`](@ref).

One can also check whether the fitness completed with [`Base.isready`](@ref) method.
"""
abstract type AbstractFitnessEvaluationJob{F} end

"""
    async_update_fitness!(eval::AbstractAsyncEvaluator, candidates::Any;
                          force::Bool=false) ->
                          Union{AbstractFitnessEvaluationJob, Nothing}

Asynchronously calculate the fitnesses.
*candidates* should support the `iterate()` interface and return
[`Candidate`](@ref) elements.

Returns [`AbstractFitnessEvaluationJob`](@ref) or `nothing` if *candidates*
does not contain any candidates.

# Arguments
 - `force::Bool`: force fitness calculation even if the candidate already has
   non-NA fitness

See also: [`sync_update_fitness`](@ref)
"""
async_update_fitness!

"""
    sync_fitness_update([f], job::AbstractFitnessEvaluationJob,
                        eval::AbstractAsyncEvaluator)

Waits until fitnesses for `job` are calculated and optionally applies `f`
function to each processed candidate.

See also: [`async_update_fitness!`](@ref), [`update_fitness!`](@ref).
"""
sync_update_fitness

"""
    is_stopping(eval::AbstractAsyncEvaluator) -> Bool

Check if the evaluator is in the shutdown sequence.

When the evaluator is being shut down, it waits for the workers to finish
their jobs and stop, no new jobs could be submitted.
"""
is_stopping

"""
Default implementation of the `Evaluator`.

`FP` is the original problem's fitness type
`FA` is the fitness type actually stored by the archive.
"""
mutable struct ProblemEvaluator{FP, FA, A<:Archive, P<:OptimizationProblem} <: Evaluator{P}
    # FIXME F is the fitness type of the problem, but with current Julia it's
    # not possible to get and use it at declaration type
    problem::P
    archive::A
    num_evals::Int
    last_fitness::FP

    ProblemEvaluator(problem::P, archive::A) where {P<:OptimizationProblem, A<:Archive} =
        new{fitness_type(fitness_scheme(problem)),fitness_type(archive),A,P}(
                problem, archive,
                0, nafitness(fitness_scheme(problem)))
end

ProblemEvaluator(problem::OptimizationProblem; archiveCapacity::Int = 10) =
    ProblemEvaluator(problem, TopListArchive(fitness_scheme(problem), numdims(problem), archiveCapacity))

problem(e::Evaluator) = e.problem
num_evals(e::ProblemEvaluator) = e.num_evals

"""
    fitness(params::AbstractIndividual, e::ProblemEvaluator, tag::Int=0)

Evaluate the fitness and implicitly update the archive with the provided
parameters and calculated fitness.

Returns the fitness in the archived format.
"""
function fitness(params::AbstractIndividual, e::ProblemEvaluator, tag::Int=0)
    e.last_fitness = fit = fitness(params, e.problem)
    e.num_evals += 1
    fita = archived_fitness(fit, e.archive)
    candi = add_candidate!(e.archive, fita, params, tag, e.num_evals)
    return fita
end

"""
    last_fitness(e::Evaluator)

Get the fitness of the last evaluated candidate.

Leads to nicer code if we can first test if it is better or worse than existing candidates
and only want the fitness itself if the criteria fulfilled.
"""
last_fitness(e::Evaluator) = e.last_fitness

is_better(f1::F, f2::F, e::Evaluator) where {F} = is_better(f1, f2, fitness_scheme(e))

is_better(candidate, f, e::Evaluator) = is_better(fitness(candidate, e), f, fitness_scheme(e))

function best_of(candidate1::AbstractIndividual, candidate2::AbstractIndividual, e::Evaluator)
    f1 = fitness(candidate1, e)
    f2 = fitness(candidate2, e)
    if is_better(f1, f2, e)
        return candidate1, f1
    else
        return candidate2, f2
    end
end

function rank_by_fitness!(e::Evaluator, candidates::AbstractVector{<:Candidate{F}}) where F
    fs = fitness_scheme(e)
    sort!(update_fitness!(e, candidates);
          # FIXME use lt=fitness_scheme(a) when v0.5 #14919 would be fixed
          by=fitness, lt=(x, y) -> is_better(x, y, fs))
end

# called by check_stop_condition(OptRunController)
check_stop_condition(e::Evaluator, ctrl) = check_stop_condition(e.archive, e.problem, ctrl)


text:
export step!, develop!

function update_cell_states!(m::Model)
    for cell in m.cells
        inputs = vcat(cell.params, cell.state, cell.input)
        outputs = m.cont.cell_state_update(inputs)
        cell.state = outputs[1:m.cfg["n_cell_state"]]
        cell.output = outputs[end]
    end
end

function update_conn_states!(m::Model)
    for cell in m.cells
        for conn in cell.conns
            inputs = vcat(conn.params, conn.state, cell.output, conn.input)
            outputs = m.cont.conn_state_update(inputs)
            conn.state = outputs[1:m.cfg["n_conn_state"]]
            conn.output = outputs[end]
        end
    end
end

function conn_communication!(m::Model)
    for cell in m.cells
        cell.input = 0
        for conn in cell.conns
            conn.input = 0
        end
    end
    for cell in m.cells
        for conn in cell.conns
            conn.dest[].input += conn.output
        end
    end
    for cell in m.cells
        cell.input = min(1.0, max(-1.0, cell.input))
        for conn in cell.conns
            conn.input = min(1.0, max(-1.0, conn.input))
        end
    end
end

function update_params!(m::Model)
    for cell in m.cells
        inputs = vcat(cell.params, cell.state)
        cell.params = m.cont.cell_param_update(inputs)
        for conn in cell.conns
            inputs = vcat(cell.params, get_params(conn.dest[], m.cfg["n_cell_params"]),
                          conn.params, conn.state)
            conn.params = m.cont.conn_param_update(inputs)
        end
    end
end

function add_cells!(m::Model)
    new_cells = Array{Cell}(undef, 0)
    for cell in m.cells
        div = m.cont.cell_division(cell.params)
        if div && ((length(m.cells) + length(new_cells)) < m.cfg["cells_max"])
            params = m.cont.new_cell_params(cell.params)
            cell = Cell(m.cfg, params)
            push!(new_cells, cell)
        end
    end
    append!(m.cells, new_cells)
    nothing
end

function remove_cells!(m::Model)
    deleted = false
    apop_inds = findall(map(c->~c.interface && m.cont.cell_death(c.params),
                            m.cells))
    if length(apop_inds) > 0
        deleted = true
        deleteat!(m.cells, apop_inds)
    end
    deleted
end

function connect!(m::Model, source::Cell, dest::Obj, new_conns::Array{Conn})
    inputs = vcat(source.params, get_params(dest, m.cfg["n_cell_params"]))
    do_connect = m.cont.connect(inputs)
    if do_connect
        params = m.cont.new_conn_params(inputs)
        conn = Conn(m.cfg, source, dest, params)
        push!(new_conns, conn)
    end
    Int64(do_connect)
end

function connect_cells!(m::Model)
    nconns = 0
    for c in m.cells
        nconns += length(c.conns)
    end
    if nconns == m.cfg["conns_max"]; return nothing; end
    all_new_conns = Array{Array{Conn, 1}}(undef, 0)
    for source in m.cells
        new_conns = Array{Conn}(undef, 0)
        dests = [c.dest[] for c in source.conns]
        for dest in setdiff(m.cells, dests)
            if nconns < m.cfg["conns_max"]
                nconns += connect!(m, source, dest, new_conns)
            end
            for conn in setdiff(dest.conns, dests)
                if nconns < m.cfg["conns_max"]
                    nconns += connect!(m, source, conn, new_conns)
                end
            end
        end
        push!(all_new_conns, new_conns)
    end
    for i in eachindex(m.cells)
        append!(m.cells[i].conns, all_new_conns[i])
    end
end

function disconnect_cells!(m::Model)
    deleted = false
    for source in m.cells
        dinds = Array{Int64}(undef, 0)
        for i in eachindex(source.conns)
            conn = source.conns[i]
            inputs = vcat(source.params,
                          get_params(conn.dest[], m.cfg["n_cell_params"]),
                          conn.params, conn.state)
            do_disconnect = m.cont.disconnect(inputs)
            if do_disconnect
                push!(dinds, i)
            end
        end
        if length(dinds) > 0
            deleted = true
            deleteat!(source.conns, dinds)
        end
    end
    deleted
end

function clean_connections!(m::Model)
    while true
        deleted = false
        for cell in m.cells
            dinds = Array{Int64}(undef, 0)
            for ci in eachindex(cell.conns)
                conn = cell.conns[ci]
                if typeof(conn.dest[]) == Cell
                    if !(conn.dest[] in m.cells)
                        push!(dinds, ci)
                    end
                else
                    in_conns = false
                    for c2 in m.cells
                        if conn.dest[] in c2.conns
                            in_conns = true
                            break
                        end
                    end
                    if !in_conns
                        push!(dinds, ci)
                    end
                end
            end
            if length(dinds) > 0
                deleted = true
                deleteat!(cell.conns, dinds)
            end
        end
        if !deleted
            break
        end
    end
end

function step!(m::Model)
    t = m.state[1] + 1
    update_cell_states!(m)
    update_conn_states!(m)
    conn_communication!(m)
    if mod(t, m.cfg["T_learn"]) == 0
        update_params!(m)
    end
    if mod(t, m.cfg["T_devo"]) == 0
        add_cells!(m)
        connect_cells!(m)
        d1 = remove_cells!(m)
        d2 = disconnect_cells!(m)
        if (d1 || d2)
            clean_connections!(m)
        end
    end
    m.state[1] = t
    nothing
end

function step!(m::Model, inputs::Array{Float64})
    set_input!(m, inputs)
    step!(m)
    get_output(m)
end

function develop!(m::Model)
    for i in 1:m.cfg["init_devo"]
        update_cell_states!(m)
        update_conn_states!(m)
        conn_communication!(m)
        update_params!(m)
        add_cells!(m)
        connect_cells!(m)
        d1 = remove_cells!(m)
        d2 = disconnect_cells!(m)
        if (d1 || d2)
            clean_connections!(m)
        end
    end
end


text:
module BasicTests

using Photon, Test

function KorATest()
    a = randn(10,10)
    b = KorA(a)
    @test size(a) == size(b)
    @test isapprox(sum(a), sum(b))
end

function contextTest()
    setContext(device=:cpu)
    @test getContext().device == :cpu

    setContext( dtype=Float32)
    @test getContext().dtype == Float32

    setContext(deviceId = 2)
    @test getContext().deviceId == 2

    resetContext()
end

function cpuTest()
    setContext(device=:cpu)

    a = KorA(randn(10,10))
    @test a isa Array
end

@testset "Basic" begin
    KorATest()
    contextTest()
    cpuTest()
end

end


text:
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

@testset "CompactSizeUInt" begin
    @testset "Read" begin
        tests = [([0x01], 1),
                 ([0xfd, 0xd0, 0x24], 9424),
                 ([0xff, 0x70, 0x9a, 0xeb, 0xb4, 0xbb, 0x7f, 0x00, 0x00], 140444170951280)]
        for t in tests
            n = CompactSizeUInt(IOBuffer(t[1]))
            @test n.value == t[2]
        end
    end
    @testset "Serialize" begin
        want = [0xfd, 0xfe, 0x00]
        @test serialize(CompactSizeUInt(254)) == want
        want = [0xfe, 0xff, 0xff, 0xff, 0xff]
        @test serialize(CompactSizeUInt(0x100000000-1)) == want
        want = [0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff]
        @test serialize(CompactSizeUInt(0x10000000000000000-1)) == want
        @test_throws ErrorException CompactSizeUInt(0x10000000000000000)
    end
end


text:
"""
Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis (Farquhar et al., 1980;
von Caemmerer and Farquhar, 1981).

The definition:

- `Tᵣ`: the reference temperature (°C) at which other parameters were measured
- `VcMaxRef`: maximum rate of Rubisco activity (``μmol\\ m^{-2}\\ s^{-1}``)
- `JMaxRef`: potential rate of electron transport (``μmol\\ m^{-2}\\ s^{-1}``)
- `RdRef`: mitochondrial respiration in the light at reference temperature (``μmol\\ m^{-2}\\ s^{-1}``)
- `TPURef`: triose phosphate utilization-limited photosynthesis rate (``μmol\\ m^{-2}\\ s^{-1}``)
- `Eₐᵣ`: activation energy (``J\\ mol^{-1}``), or the exponential rate of rise for Rd.
- `O₂`: intercellular dioxygen concentration (``ppm``)
- `Eₐⱼ`: activation energy (``J\\ mol^{-1}``), or the exponential rate of rise for JMax.
- `Hdⱼ`: rate of decrease of the function above the optimum (also called EDVJ) for JMax.
- `Δₛⱼ`: entropy factor for JMax.
- `Eₐᵥ`: activation energy (``J\\ mol^{-1}``), or the exponential rate of rise for VcMax.
- `Hdᵥ`: rate of decrease of the function above the optimum (also called EDVC) for VcMax.
- `Δₛᵥ`: entropy factor for VcMax.
- `α`: quantum yield of electron transport (``mol_e\\ mol^{-1}_{quanta}``). See also eq. 4 of
Medlyn et al. (2002) and its implementation in [`get_J`](@ref)
- `θ`: determines the curvature of the light response curve for `J~PPFD`. See also eq. 4 of
Medlyn et al. (2002) and its implementation in [`get_J`](@ref)

The default values of the temperature correction parameters are taken from
[plantecophys](https://remkoduursma.github.io/plantecophys/). If there is no negative effect
of high temperatures on the reaction (Jmax or VcMax), then Δₛ can be set to 0.0.

# Note

Medlyn et al. (2002) found relatively low influence ("a slight effect") of α and θ. They also
say that Kc, Ko and Γ* "are thought to be intrinsic properties of the Rubisco enzyme
and are generally assumed constant among species".

# See also

- [`get_J`](@ref)
- [`photosynthesis`](@ref)

# References

Caemmerer, S. von, et G. D. Farquhar. 1981. « Some Relationships between the Biochemistry of
Photosynthesis and the Gas Exchange of Leaves ». Planta 153 (4): 376‑87.
https://doi.org/10.1007/BF00384257.

Farquhar, G. D., S. von von Caemmerer, et J. A. Berry. 1980. « A biochemical model of
photosynthetic CO2 assimilation in leaves of C3 species ». Planta 149 (1): 78‑90.

Medlyn, B. E., E. Dreyer, D. Ellsworth, M. Forstreuter, P. C. Harley, M. U. F. Kirschbaum,
X. Le Roux, et al. 2002. « Temperature response of parameters of a biochemically based model
of photosynthesis. II. A review of experimental data ». Plant, Cell & Environment 25 (9): 1167‑79.
https://doi.org/10.1046/j.1365-3040.2002.00891.x.

# Examples

```julia
Get the fieldnames:
fieldnames(Fvcb)
# Using default values for the model:
A = Fvcb()

A.Eₐᵥ
```
"""
struct Fvcb{T} <: AbstractAModel
    Tᵣ::T
    VcMaxRef::T
    JMaxRef::T
    RdRef::T
    TPURef::T
    Eₐᵣ::T
    O₂::T
    Eₐⱼ::T
    Hdⱼ::T
    Δₛⱼ::T
    Eₐᵥ::T
    Hdᵥ::T
    Δₛᵥ::T
    α::T
    θ::T
end

function Fvcb(;Tᵣ = 25.0, VcMaxRef = 200.0, JMaxRef = 250.0, RdRef = 0.6, TPURef = 9999., Eₐᵣ = 46390.0,
    O₂= 210.0, Eₐⱼ = 29680.0, Hdⱼ = 200000.0, Δₛⱼ = 631.88, Eₐᵥ = 58550.0, Hdᵥ = 200000.0,
    Δₛᵥ = 629.26, α = 0.425, θ = 0.90)

    Fvcb(promote(Tᵣ, VcMaxRef, JMaxRef, RdRef, TPURef, Eₐᵣ, O₂, Eₐⱼ, Hdⱼ, Δₛⱼ, Eₐᵥ, Hdᵥ, Δₛᵥ, α, θ)...)
end

function inputs(::Fvcb)
    (:PPFD,:Tₗ,:Cₛ)
end

function outputs(::Fvcb)
    (:A,:Gₛ,:Cᵢ)
end

Base.eltype(x::Fvcb) = typeof(x).parameters[1]

"""
    assimilation!(leaf::LeafModels{I,E,<:Fvcb,<:AbstractGsModel,S},constants = Constants())

Photosynthesis using the Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis
 (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981).
Computation is made following Farquhar & Wong (1984), Leuning et al. (1995), and the
MAESPA model (Duursma et al., 2012).
The resolution is analytical as first presented in Baldocchi (1994), and needs Cₛ as input.
Triose phosphate utilization (TPU) limitation is taken into account as proposed in Lombardozzi (2018) (i.e.
Aₚ = 3*TPURef, making the assumption that glycolate recycling is set to 0). TPURef is by default at 9999., meaning that if you
don't have the TPU value, TPU limitation won't have impact.

If you prefer to use Gbc, you can use the iterative implementation of the Fvcb model
[`FvcbIter`](@ref)

# Returns

Modify the first argument in place for A, Gₛ and Cᵢ:

- A: carbon assimilation (μmol[CO₂] m-2 s-1)
- Gₛ: stomatal conductance for CO₂ (mol[CO₂] m-2 s-1)
- Cᵢ: intercellular CO₂ concentration (ppm)

# Arguments

- `leaf::LeafModels{.,.,<:Fvcb,<:AbstractGsModel,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model with initialisations for:
    - `Tₗ` (°C): leaf temperature
    - `PPFD` (μmol m-2 s-1): absorbed Photosynthetic Photon Flux Density
    - `Cₛ` (ppm): surface CO₂ concentration.
    - `Dₗ` (kPa): vapour pressure difference between the surface and the air saturation
    vapour pressure in case you're using the stomatal conductance model of [`Medlyn`](@ref).
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

`Tₗ`, `PPFD`, `Cₛ` (and `Dₗ` if you use [`Medlyn`](@ref)) must be initialised by providing
them as keyword arguments (see examples). If in doubt, it is simpler to compute the energy
balance of the leaf with the photosynthesis to get those variables. See
[`energy_balance`](@ref) for more details.

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

leaf = LeafModels(photosynthesis = Fvcb(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Tₗ = 25.0, PPFD = 1000.0, Cₛ = 400.0, Dₗ = meteo.VPD)
# NB: we need  to initalise Tₗ, PPFD and Cₛ

assimilation!(leaf,meteo,Constants())
leaf.status.A
leaf.status.Cᵢ
```

# References

Baldocchi, Dennis. 1994. « An analytical solution for coupled leaf photosynthesis and
stomatal conductance models ». Tree Physiology 14 (7-8‑9): 1069‑79.
https://doi.org/10.1093/treephys/14.7-8-9.1069.

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5
(4): 919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Farquhar, G. D., S. von von Caemmerer, et J. A. Berry. 1980. « A biochemical model of
photosynthetic CO2 assimilation in leaves of C3 species ». Planta 149 (1): 78‑90.

Leuning, R., F. M. Kelliher, DGG de Pury, et E.D. Schulze. 1995. « Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.

Lombardozzi, L. D. et al. 2018.« Triose phosphate limitation in photosynthesis models
reduces leaf photosynthesis and global terrestrial carbon storage ». Environmental Research
Letters 13.7: 1748-9326. https://doi.org/10.1088/1748-9326/aacf68.
"""
function assimilation!(leaf::LeafModels{I,E,<:Fvcb,<:AbstractGsModel,S}, meteo,
    constants = Constants()) where {I,E,S}

    # Tranform Celsius temperatures in Kelvin:
    Tₖ = leaf.status.Tₗ - constants.K₀
    Tᵣₖ = leaf.photosynthesis.Tᵣ - constants.K₀

    # Temperature dependence of the parameters:
    Γˢ = Γ_star(Tₖ,Tᵣₖ,constants.R) # Gamma star (CO2 compensation point) in μmol mol-1
    Km = get_km(Tₖ,Tᵣₖ,leaf.photosynthesis.O₂,constants.R) # effective Michaelis–Menten coefficient for CO2

    # Maximum electron transport rate at the given leaf temperature (μmol m-2 s-1):
    JMax = arrhenius(leaf.photosynthesis.JMaxRef,leaf.photosynthesis.Eₐⱼ,Tₖ,Tᵣₖ,
                        leaf.photosynthesis.Hdⱼ,leaf.photosynthesis.Δₛⱼ,constants.R)
    # Maximum rate of Rubisco activity at the given leaf temperature (μmol m-2 s-1):
    VcMax = arrhenius(leaf.photosynthesis.VcMaxRef,leaf.photosynthesis.Eₐᵥ,Tₖ,Tᵣₖ,
                        leaf.photosynthesis.Hdᵥ,leaf.photosynthesis.Δₛᵥ,constants.R)
    # Rate of mitochondrial respiration at the given leaf temperature (μmol m-2 s-1):
    Rd = arrhenius(leaf.photosynthesis.RdRef,leaf.photosynthesis.Eₐᵣ,Tₖ,Tᵣₖ,constants.R)
    # Rd is also described as the CO2 release in the light by processes other than the PCO
    # cycle, and termed "day" respiration, or "light respiration" (Harley et al., 1986).

    # Actual electron transport rate (considering intercepted PAR and leaf temperature):
    J = get_J(leaf.status.PPFD, JMax, leaf.photosynthesis.α, leaf.photosynthesis.θ) # in μmol m-2 s-1
    # RuBP regeneration
    Vⱼ = J / 4

    # Stomatal conductance (mol[CO₂] m-2 s-1), dispatched on type of first argument (Gs_mod):
    gs_mod = gs_closure(leaf,meteo)

    Cᵢⱼ = get_Cᵢⱼ(Vⱼ,Γˢ,leaf.status.Cₛ,Rd,leaf.stomatal_conductance.g0,gs_mod)

    # Electron-transport-limited rate of CO2 assimilation (RuBP regeneration-limited):
    Wⱼ = Vⱼ * (Cᵢⱼ - Γˢ) / (Cᵢⱼ + 2.0 * Γˢ) # also called Aⱼ
    # See Von Caemmerer, Susanna. 2000. Biochemical models of leaf photosynthesis.
    # Csiro publishing, eq. 2.23.
    # NB: here the equation is modified because we use Vⱼ instead of J, but it is the same.

    # If Rd is larger than Wⱼ, no assimilation:
    if Wⱼ - Rd < 1.0e-6
        Cᵢⱼ = Γˢ
        Wⱼ = Vⱼ * (Cᵢⱼ - Γˢ) / (Cᵢⱼ + 2.0 * Γˢ)
    end

    Cᵢᵥ = get_Cᵢᵥ(VcMax,Γˢ,leaf.status.Cₛ,Rd,leaf.stomatal_conductance.g0,gs_mod,Km)

    # Rubisco-carboxylation-limited rate of CO₂ assimilation (RuBP activity-limited):
    if Cᵢᵥ <= 0.0 || Cᵢᵥ > leaf.status.Cₛ
        Wᵥ = 0.0
    else
        Wᵥ = VcMax * (Cᵢᵥ - Γˢ) / (Cᵢᵥ + Km)
    end

    # Net assimilation (μmol m-2 s-1)
    leaf.status.A = min(Wᵥ,Wⱼ,3*leaf.photosynthesis.TPURef) - Rd

    # Stomatal conductance (mol[CO₂] m-2 s-1)
    leaf.status.Gₛ = gs(leaf,gs_mod)
    # replace by ifelse directly ? Should be faster as `max()` add some tests.\

    # Intercellular CO₂ concentration (Cᵢ, μmol mol)
    leaf.status.Cᵢ = min(leaf.status.Cₛ, leaf.status.Cₛ - leaf.status.A / leaf.status.Gₛ)
    nothing
end

# With keyword arguments (better for users)
# function assimilation(A_mod::Fvcb, Gs_mod::AbstractGsModel; Tₗ, PPFD, Cₛ, Rh = missing,
#                         VPD = missing, ψₗ = missing)

#     environment = MutableNamedTuple(Tₗ = Tₗ, PPFD = PPFD, Rh = Rh, Cₛ= Cₛ, VPD = VPD, ψₗ = ψₗ)

#     assimilation(A_mod,Gs_mod,environment,Constants())
# end

"""
Rate of electron transport J (``μmol\\ m^{-2}\\ s^{-1}``), computed using the smaller root
of the quadratic equation (eq. 4 from Medlyn et al., 2002):

    θ * J² - (α * PPFD + JMax) * J + α * PPFD * JMax

NB: we use the smaller root because considering the range of values for θ and α (quite stable),
and PPFD and JMax, the function always tends to JMax with high PPFD with the smaller root (behavior we
are searching), and the opposite with the larger root.

# Returns

A tuple with (A, Gₛ, Cᵢ):

- A: carbon assimilation (μmol m-2 s-1)
- Gₛ: stomatal conductance (mol m-2 s-1)
- Cᵢ: intercellular CO₂ concentration (ppm)
# Arguments

- `PPFD`: absorbed photon irradiance (``μmol_{quanta}\\ m^{-2}\\ s^{-1}``)
- `α`: quantum yield of electron transport (``mol_e\\ mol^{-1}_{quanta}``)
- `JMax`: maximum rate of electron transport (``μmol\\ m^{-2}\\ s^{-1}``)
- `θ`: determines the shape of the non-rectangular hyperbola (-)

# References

Medlyn, B. E., E. Dreyer, D. Ellsworth, M. Forstreuter, P. C. Harley, M. U. F. Kirschbaum,
X. Le Roux, et al. 2002. « Temperature response of parameters of a biochemically based model
of photosynthesis. II. A review of experimental data ». Plant, Cell & Environment 25 (9): 1167‑79.
https://doi.org/10.1046/j.1365-3040.2002.00891.x.

Von Caemmerer, Susanna. 2000. Biochemical models of leaf photosynthesis. Csiro publishing.

# Examples

```jldoctest; setup = :(using PlantBiophysics)
# Using default values for the model:
julia> A = Fvcb()
Fvcb{Float64}(25.0, 200.0, 250.0, 0.6, 46390.0, 210.0, 29680.0, 200000.0, 631.88, 58550.0, 200000.0, 629.26, 0.425, 0.9)

julia> PlantBiophysics.get_J(1500, A.JMaxRef, A.α, A.θ)
236.11111111111111
```
"""
function get_J(PPFD, JMax, α, θ)
  (α * PPFD + JMax - sqrt((α * PPFD + JMax)^2 - 4 * α * θ * PPFD * JMax)) / (2 * θ)
end

"""
Analytic resolution of Cᵢ when the rate of electron transport is limiting (``μmol\\ mol^{-1}``)

# Arguments

- `Vⱼ`: RuBP regeneration (J/4.0, ``μmol\\ m^{-2}\\ s^{-1}``)
- `Γˢ`: CO2 compensation point ``Γ^⋆`` (``μmol\\ mol^{-1}``)
- `Cₛ`: stomatal CO₂ concentration (``μmol\\ mol^{-1}``)
- `Rd`: day respiration (``μmol\\ m^{-2}\\ s^{-1}``)
- `g0`: residual stomatal conductance (``μmol\\ m^{-2}\\ s^{-1}``)
- `gs_mod`: stomatal conductance term computed from a given implementation of a Gs model,
e.g. [`Medlyn`](@ref).

# References

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5
(4): 919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Wang and Leuning, 1998
"""
function get_Cᵢⱼ(Vⱼ,Γˢ,Cₛ,Rd,g0,gs_mod)
    a = g0 + gs_mod * (Vⱼ - Rd)
    b = (1.0 - Cₛ * gs_mod) * (Vⱼ - Rd) + g0 * (2.0 * Γˢ - Cₛ) -
        gs_mod * (Vⱼ * Γˢ + 2.0 * Γˢ * Rd)
    c = -(1.0 - Cₛ * gs_mod) * Γˢ * (Vⱼ + 2.0 * Rd) -
        g0 * 2.0 * Γˢ * Cₛ

    return max_root(a,b,c)
end

"""
Analytic resolution of Cᵢ when the Rubisco activity is limiting (``μmol\\ mol^{-1}``)

# Arguments

- `VcMAX`: maximum rate of Rubisco activity(``μmol\\ m^{-2}\\ s^{-1}``)
- `Γˢ`: CO2 compensation point ``Γ^⋆`` (``μmol\\ mol^{-1}``)
- `Cₛ`: stomatal CO₂ concentration (``μmol\\ mol^{-1}``)
- `Rd`: day respiration (``μmol\\ m^{-2}\\ s^{-1}``)
- `g0`: residual stomatal conductance (``μmol\\ m^{-2}\\ s^{-1}``)
- `gs_mod`: stomatal conductance term computed from a given implementation of a Gs model,
e.g. [`Medlyn`](@ref).
- `Km`: effective Michaelis–Menten coefficient for CO2 (``μ mol\\ mol^{-1}``)
"""
function get_Cᵢᵥ(VcMAX,Γˢ,Cₛ,Rd,g0,gs_mod,Km)
    a = g0 + gs_mod * (VcMAX - Rd)
    b = (1.0 - Cₛ * gs_mod) * (VcMAX - Rd) + g0 * (Km - Cₛ) - gs_mod * (VcMAX * Γˢ + Km * Rd)
    c = -(1.0 - Cₛ * gs_mod) * (VcMAX * Γˢ + Km * Rd) - g0 * Km * Cₛ

    return max_root(a,b,c)
end

"""
Maximum value between two roots of a quadratic equation.
"""
function max_root(a,b,c)
    Δ = b^2.0 - 4.0 * a * c
    x1 = (-b + sqrt(Δ)) / (2.0 * a)
    x2 = (-b - sqrt(Δ)) / (2.0 * a)
    return max(x1,x2)
end


text:
using MCMCChain
using Test

## CHAIN TESTS
# Define the experiment
n_iter = 500
n_name = 3
n_chain = 2

# some sample experiment results
val = randn(n_iter, n_name, n_chain) .+ [1, 2, 3]'
val = hcat(val, rand(1:2, n_iter, 1, n_chain))

# construct a Chains object
chn = Chains(val, start = 1, thin = 2)

@testset "basic chains functions" begin
    @test first(chn) == 1
    @test step(chn) == 2
    @test last(chn) == 999
    @test size(chn) == (999, 4, 2)
    @test keys(chn) == ["Param#1", "Param#2", "Param#3", "Param#4"]
    @test isa(chn[:,1,:], MCMCChain.AbstractChains)
    @test length(vec(chn[:,1,:].value)) == n_chain * n_iter
    @test all(collect(skipmissing(chn[:,1,1].value)) .== val[:,1,1])
    @test all(chn[:,1,2].value .== val[:,1,2])
    @test all(MCMCChain.indiscretesupport(chn) .== [false, false, false, true])
end

@testset "function tests" begin
    # do not test the following functions
    # - wrtsp
    # - window2inds (tested above, see getindex)
    # -

    # the following tests only check if the function calls work!
    @test MCMCChain.diag_all(rand(100, 2), :weiss, 1, 1, 1) != nothing
    @test MCMCChain.diag_all(rand(100, 2), :hangartner, 1, 1, 1) != nothing
    @test MCMCChain.diag_all(rand(100, 2), :billingsley, 1, 1, 1) != nothing

    @test isa(discretediag(chn[:,4,:]), MCMCChain.ChainSummary)
    @test isa(gelmandiag(chn[:,1,:]), MCMCChain.ChainSummary)
    @test isa(gewekediag(chn[:,1,:]), MCMCChain.ChainSummary)
    @test isa(heideldiag(chn[:,1,:]), MCMCChain.ChainSummary)
    @test isa(rafterydiag(chn[:,1,:]), MCMCChain.ChainSummary)
end


text:
# magma_dgesvd wrapper
function magmaDgesvd(jobu, jobvt, m, n, A, lda, s, U,
	 ldu, VT, ldvt, work, lwork, info)
    ccall((:magma_dgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{Cdouble},
			    Cint, PtrOrCuPtr{Cdouble}, PtrOrCuPtr{Cdouble},
				Cint, PtrOrCuPtr{Cdouble}, Cint,
				PtrOrCuPtr{Cdouble}, Cint, PtrOrCuPtr{Cint}),

                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, info)
end


text:
function itr2bin(x)
    res = 0
    for v in x
        res += v
        res <<= 1
    end
    return res >> 1
end

function itr2ind(x, sizes)
    n = length(sizes)
    res = 0
    _mul = 1
    for i in n : -1 : 1
        res += x[i] * _mul
        _mul *= sizes[i]
    end
    return res
end

# fill array v with bits of number x
function fill_bits!(v, x)
    for i in length(v) : -1 : 1
        v[i] = x & 1
        x >>= 1
    end
end

# fill array v with digits of number x in base-system
function fill_numbers!(v, x; base :: Int)
    for i in length(v) : -1 : 1
        v[i] = x % base
        x = div(x, base)
    end
end

# set bits with numbers vs[1],
# ..., vs[k] to 1
# return 2^(vs[1] - 1) + ... + 2^(vs[k] - 1)
function set_bits(vs)
    res = 0
    for num in vs
        res += 2^(num - 1)
    end
    return res
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
using Documenter, DiffImages

makedocs(;
    modules = [DiffImages],
    doctest = true,
    authors = "Som V Tambe <somvt@iitk.ac.in>",
    repo = "https://github.com/SomTambe/DiffImages.jl/blob/{commit}{path}#L{line}",
    sitename = "DiffImages.jl",
    # format = Documenter.HTML(;
    #     prettyurls = get(ENV, "CI", "false") == "true",
    #     canonical = "https://fluxml.ai/Flux3D.jl",
    #     assets = String["assets/favicon.ico"],
    #     analytics = "UA-154580699-2",
    # ),
    pages = [
        "Home" => "index.md",
        # "Tutorials" => [
        # ],
        "Colorspace Transforms" => "colors/index.md",
        "Geometry modules" =>
            ["Warp Modules" => "geometry/warp.md"],
        "Tutorials" =>
            ["Rotation transforms" => "tutorials/rotate/rotate.md",
            "Homography warps" => "tutorials/homog/warp.md"]
    ],
)

deploydocs(devbranch = "main", repo = "github.com/SomTambe/DiffImages.jl.git")