text:
function _pairwise_Hudson_lxl(data::PopData)
    !isbiallelic(data) && throw(error("Data must be biallelic to use the Hudson estimator"))
    idx_pdata = groupby(data.genodata, :population)
    pops = getindex.(keys(idx_pdata), :population)
    nloci = data.metadata.loci
    locnames = loci(data)
    allpairs = pairwisepairs(pops)
    npairs = length(collect(allpairs))
    results = Vector{Vector{Float64}}(undef, npairs)
    p1 = repeat(getindex.(allpairs,1), inner = nloci)
    p2 = repeat(getindex.(allpairs,2), inner = nloci)
    locs = repeat(locnames, outer = npairs)
    @inbounds @sync for (i,j) in enumerate(allpairs)
        Base.Threads.@spawn begin
            pop1 = reshape(idx_pdata[(population = j[1],)].genotype, :, nloci)
            pop2 = reshape(idx_pdata[(population = j[2],)].genotype, :, nloci)
            results[i] = _hudson_fst_lxl(pop1, pop2)
        end
    end
    return DataFrame(:pop1 => PooledArray(p1, compress = true), :pop2 => PooledArray(p2, compress = true),:locus => PooledArray(locs, compress = true) ,:fst => reduce(vcat, results))
end

function _hudson_fst_lxl(population_1::T, population_2::T) where T<:AbstractMatrix
    @views [_hudson_fst(population_1[:,i], population_2[:,i]) for i in 1:size(population_1,2)]
end

function _nei_fst_lxl(population_1::T, population_2::T) where T<:AbstractMatrix
    n =  hcat(map(nonmissing, eachcol(population_1)), map(nonmissing, eachcol(population_2)))
    # number of populations represented per locus
    n_pop_per_loc = map(countnonzeros, eachrow(n))
    # corrected n for population size
    corr_n_per_loc = n_pop_per_loc ./ map(reciprocalsum, eachrow(n))
    # observed heterozygosity
    het_obs_p1 = map(_hetero_obs, eachcol(population_1))
    het_obs_p2 = map(_hetero_obs, eachcol(population_2))
    # expected heterozygosity
    het_exp_p1 = map(_hetero_exp, eachcol(population_1))
    het_exp_p2 = map(_hetero_exp, eachcol(population_2))
    # genic diversity for population 1 and 2
    p1_nei = _genediversitynei87.(het_exp_p1, het_obs_p1, corr_n_per_loc)
    p2_nei = _genediversitynei87.(het_exp_p2, het_obs_p2, corr_n_per_loc)
    # mean genic diversity
    HS = map(i -> mean(skipmissing(i)), zip(p1_nei, p2_nei))
    # allele freqs for population 1 and 2
    alle_freq_p1 = map(allelefreq, eachcol(population_1))
    alle_freq_p2 = map(allelefreq, eachcol(population_2))
    # sum of the squared average allele frequencies
    avg_freq = avg_allelefreq.(zip(alle_freq_p1, alle_freq_p2),2) .|> values .|> sum
    # average observed heterozygosity per locus
    Het_obs = map(i -> mean(skipmissing(i)), zip(het_obs_p1, het_obs_p2))

    Ht = @inbounds 1.0 .- avg_freq .+ (HS ./ corr_n_per_loc ./ n_pop_per_loc) - (Het_obs ./ 2.0 ./ corr_n_per_loc ./ n_pop_per_loc)
    DST = @inbounds Ht .- HS
    DST′ = @inbounds n_pop_per_loc ./ (n_pop_per_loc .- 1) .* DST
    HT′ = @inbounds HS .+ DST′
    round.(skipinfnan(DST′) ./ skipinfnan(HT′), digits = 5)
end

function _pairwise_Nei_lxl(data::PopData)
    idx_pdata = groupby(data.genodata, :population)
    pops = getindex.(keys(idx_pdata), :population)
    nloci = data.metadata.loci
    locnames = loci(data)
    allpairs = pairwisepairs(pops)
    npairs = length(collect(allpairs))
    results = Vector{Vector{Float64}}(undef, npairs)
    p1 = repeat(getindex.(allpairs,1), inner = nloci)
    p2 = repeat(getindex.(allpairs,2), inner = nloci)
    locs = repeat(locnames, outer = npairs)
    @inbounds @sync for (i,j) in enumerate(allpairs)
        Base.Threads.@spawn begin
            pop1 = reshape(idx_pdata[(population = j[1],)].genotype, :, nloci)
            pop2 = reshape(idx_pdata[(population = j[2],)].genotype, :, nloci)
            results[i] = _nei_fst_lxl(pop1, pop2)
        end
    end
    return DataFrame(:pop1 => PooledArray(p1, compress = true), :pop2 => PooledArray(p2, compress = true),:locus => PooledArray(locs, compress = true) ,:fst => reduce(vcat, results))
end


text:
################################################################################
#
# General utilities
#
################################################################################

export extend, restrict, sz_abc, gaussian_filt, contr2abs, abs2contr, gradprec_contr2abs, gendata, compute_y, proj_bounds

function sz_abc(n::NTuple{2, Int64}, nabc::NTuple{4, Int64})

    return n.+(nabc[1]+nabc[2], nabc[3]+nabc[4])

end

function extend(m::Array{T, 2}, nabc::NTuple{4, Int64}) where T <: RuC

    m_ext = Array{T}(undef, sz_abc(size(m), nabc))
    m_ext[nabc[1]+1:end-nabc[2], nabc[3]+1:end-nabc[4]] .= m
    m_ext[1:nabc[1], :] .= reshape(m_ext[nabc[1]+1, :], 1, :) # left
    m_ext[end-nabc[2]+1:end, :] .= reshape(m_ext[end-nabc[2], :], 1, :) # right
    m_ext[:, 1:nabc[3]] .= reshape(m_ext[:, nabc[3]+1], :, 1) # top
    m_ext[:, end-nabc[4]+1:end] .= reshape(m_ext[:, end-nabc[4]], :, 1) # bottom

    return m_ext

end

function extend(m::Array{T, 3}, nabc::NTuple{4, Int64}) where T <: RuC

    n = (size(m, 1), size(m, 2))
    m_ext = Array{T}(undef, (sz_abc(size(m), nabc)..., size(m, 3)))
    m_ext[nabc[1]+1:end-nabc[2], nabc[3]+1:end-nabc[4], :] .= m
    m_ext[1:nabc[1], :, :] .= reshape(m_ext[nabc[1]+1, :, :], 1, size(m_ext, 2), :) # left
    m_ext[end-nabc[2]+1:end, :, :] .= reshape(m_ext[end-nabc[2], :, :], 1, size(m_ext, 2), :) # right
    m_ext[:, 1:nabc[3], :] .= reshape(m_ext[:, nabc[3]+1, :], :, 1, size(m_ext, 3)) # top
    m_ext[:, end-nabc[4]+1:end, :] .= reshape(m_ext[:, end-nabc[4]], :, 1, size(m_ext, 3)) # bottom

    return m_ext

end

## Restriction for 2D images

function restrict(m::Array{T, 2}, nabc::NTuple{4, Int64}) where T <: RuC

    return m[nabc[1]+1:end-nabc[2], nabc[3]+1:end-nabc[4]]

end

function restrict(m::Array{T, 3}, nabc::NTuple{4, Int64}) where T <: RuC

    return m[nabc[1]+1:end-nabc[2], nabc[3]+1:end-nabc[4], :]

end

## Utilities for model input preprocessing and gradient postprocessing: contrast to absolute properties / effective to total domain

# Contrast preconditioning
function contr2abs(xvec::Array{R, 1}, mask::BitArray{2}, mb::Array{R, 2})
# Preprocessing steps: - effective to global domain
#                      - contrast to absolute properties

    x = zeros(R, size(mb))
    x[mask] .= xvec
    return mb.*(R(1).+x)

end
function abs2contr(m::Array{R, 2}, mask::BitArray{2}, mb::Array{R, 2})

	return ((m.-mb)./mb)[mask]

end

function gradprec_contr2abs(grad::Array{R, 2}, mask::BitArray{2}, mb::Array{R, 2})
# Postprocessing steps: - preconditioning by diagonal background model

    return mb[mask].*grad[mask]

end

## Smoothing utilities

function gaussian_filt(m_raw::Array{R, 2}, σ::R)

    return imfilter(m_raw, Kernel.gaussian(σ))

end

## Generate data

function gendata(m::Array{R, 2}, h::NTuple{2, R}, freqs::Array{R, 1}, x_rcv::Array{R, 1}, z_rcv::R, x_src::Array{R, 1}, z_src::R, src_enc::Array{C, 2}; nabc::Union{Nothing, NTuple{4, Int64}} = nothing, dat::Union{Nothing, Array{Array{C, 2}, 1}} = nothing, verbose::Bool = false)

	# Initializing data output
	n = size(m)
	Nr = length(x_rcv)
	Ns = size(src_enc, 2)
	Nf = length(freqs)
	rec = Array{Array{C, 2}, 1}(undef, Nf)

    # Modeling (loop over frequencies)
    for idxf = 1:Nf

		if verbose
			println("Generating data for frequency: ", freqs[idxf], " Hz")
		end

		# Set absorbing boundary size (frequency-dependent)
        if nabc == nothing
            nabc = computeNabc(h, m, freqs[idxf])
        end

        # Setting source/receiver restriction/injection operators & source weigths
        Pr = restrOpMat(n, h, x_rcv, z_rcv, nabc)
        Ps = restrOpMat(n, h, x_src, z_src, nabc)
        Fsrc = Matrix(adjoint(Ps)*src_enc)

        # Solve for wavefield (forward)
        H = HelmholtzLinOp(h, m, freqs[idxf], nabc)
        U = C.(H\Fsrc)

        # Receiver restriction
		if dat == nothing
			rec[idxf] = Pr*U # Data recorded
		else
			rec[idxf] = dat[idxf]*src_enc.-Pr*U # Data residual
		end

    end

	return rec

end

function compute_y(m::Array{R, 2}, h::NTuple{2, R}, freqs::Array{R, 1}, x_rcv::Array{R, 1}, z_rcv::R, x_src::Array{R, 1}, z_src::R, src_enc::Array{C, 2}, dat::Array{Array{C, 2}, 1}, λ::R; nabc::Union{Nothing, NTuple{4, Int64}} = nothing, weight_array::Array{R, 2}, mode_y::String = "optimal", c::Union{Nothing, Array{R, 1}} = nothing, β::Union{R, Nothing} = nothing)

	# Initializing data output
	n = size(m)
	Nr = length(x_rcv)
	Ns = size(src_enc, 2)
	Nf = length(freqs)
	y = Array{Array{C, 2}, 1}(undef, Nf)

	# Loop over frequencies
	for idxf = 1:Nf

		# Encoding data
		dat_enc = dat[idxf]*src_enc

         # Set absorbing boundary size (frequency-dependent)
        if nabc == nothing
            nabc = computeNabc(h, m, freqs[idxf])
        end

        # Setting source/receiver restriction/injection operators & source weigths
        Pr = restrOpMat(n, h, x_rcv, z_rcv, nabc)
        Ps = restrOpMat(n, h, x_src, z_src, nabc)
        Fsrc = Matrix(adjoint(Ps)*src_enc)

        # Helmholtz operator
        H = HelmholtzLinOp(h, m, freqs[idxf], nabc)
		if mode_y != "optimal"
			H = lu(H)
		end

		# Computing y, if not provided (y = data residual)
		if c == nothing && mode_y != "residual"
			ηd = norm(dat_enc)
			ηf = norm(Fsrc)
			c = ηf./(λ.*ηd)
		end
		if mode_y == "optimal"
			W = spdiagm(0 => vec(weight_array))
			U = C.(ComplexF64.([c.*Pr; W*H])\[c.*dat_enc; W*Fsrc])
			y[idxf] = dat_enc.-Pr*U
		elseif mode_y == "residual"
			U = C.(H\Fsrc)
			y[idxf] = dat_enc.-Pr*U
		elseif mode_y == "residual2"
			Winv2 = spdiagm(0 => vec(R(1)./weight_array.^R(2)))
			Uinc = C.(H\Fsrc)
			r = dat_enc.-Pr*Uinc
			Fr = C.(adjoint(H)\(adjoint(Pr)*r))
			FWFr = Pr*(C.(H\(Winv2*Fr)))
			y[idxf] = r.-β.*c.^R(2).*FWFr
		elseif mode_y == "adjoint"
			Winv2 = spdiagm(0 => vec(R(1)./weight_array.^R(2)))
			Uinc = C.(H\Fsrc)
			r = dat_enc.-Pr*Uinc
			Fr = C.(adjoint(H)\(adjoint(Pr)*r))
			FWFr = Pr*(C.(H\(Winv2*Fr)))
			y[idxf] = r.+c.^R(2).*FWFr
		end

	end

	return y

end

## Projection on box constraints

function proj_bounds(m::Array{R, 2}, bounds::NTuple{2, Array{R, 2}})
	m_ = copy(m)
	m_[m .< bounds[1]] .= bounds[1][m .< bounds[1]]
	m_[m .> bounds[2]] .= bounds[2][m .> bounds[2]]
	return m_
end

function proj_bounds(m::Array{R, 2}, bounds::NTuple{2, R})
	m_ = copy(m)
	m_[m .< bounds[1]] .= bounds[1]
	m_[m .> bounds[2]] .= bounds[2]
	return m_
end

function proj_bounds(m::Array{R, 2}, bounds::Nothing)
	return m
end

ext: jl


text:
"""
    map(f, ds, [key]) -> KeyedDataset

Apply function `f` to each component of the [`KeyedDataset`](@ref).
Returns a new dataset with the same constraints, but new components.
The function can be applied to a subselection of components via a [`Pattern`](@ref) `key`.

# Example
```jldoctest
julia> using AxisKeys, Statistics; using AxisSets: KeyedDataset, flatten;

julia> ds = KeyedDataset(
           flatten([
               :g1 => [
                   :a => KeyedArray(zeros(3); time=1:3),
                   :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
                ],
                :g2 => [
                    :a => KeyedArray(ones(3); time=1:3),
                    :b => KeyedArray(zeros(3, 2); time=1:3, loc=[:x, :y]),
                ]
            ])...
       );

julia> r = map(a -> a .+ 100, ds, (:__, :a, :_));  # The extra `:_` is to clarify that we don't care about the dimnames.

julia> [k => mean(v) for (k, v) in r.data]  # KeyedArray printing isn't consistent in jldoctests
4-element Vector{Pair{Tuple{Symbol, Symbol}, Float64}}:
 (:g1, :a) => 100.0
 (:g1, :b) => 1.0
 (:g2, :a) => 101.0
 (:g2, :b) => 0.0
```
"""
function Base.map(f::Function, ds::KeyedDataset, keys...)
    # Ensure that our keys are patterns or converted to patterns
    patterns = isempty(keys) ? Pattern[(:__,)] : Pattern[keys...]

    # Select any components our keys match
    selected = unique(x[1:end-1] for x in dimpaths(ds) if any(p -> x in p, patterns))
    return map!(f, deepcopy(ds), ds(k -> k in selected))
end

function Base.map!(f::Function, dest::KeyedDataset, src::KeyedDataset)
    # Mutate the underlying data field to avoid calling `validate`
    # on each insertion
    data = dest.data
    for (k, v) in src.data
        data[k] = f(v)
    end

    validate(dest)
    return dest
end

"""
    mapslices(f, ds, [key]; dims) -> KeyedDataset

Apply the `mapslices` call to each of the desired components and returns a new [`KeyedDataset`](@ref).
Selection [`Pattern`](@ref)s may be provided via `key`, but components are selected by the desired `dims` by default.

# Example
```jldoctest
julia> using AxisKeys, Statistics; using AxisSets: KeyedDataset;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(zeros(3, 2); time=1:3, obj=[:a, :b]),
           :val2 => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
       );

julia> r = mapslices(sum, ds; dims=:time);  # KeyedArray printing isn't consistent in jldoctests

julia> [k => parent(parent(v)) for (k, v) in r.data]
2-element Vector{Pair{Tuple{Symbol}, Matrix{Float64}}}:
 (:val1,) => [0.0 0.0]
 (:val2,) => [3.0 3.0]
```
"""
function Base.mapslices(f::Function, ds::KeyedDataset, keys...; dims)
    patterns = if isempty(keys)
        dims isa Symbol ? Pattern[(:__, dims)] : Pattern[(:__, d) for d in dims]
    else
        Pattern[keys...]
    end
    return map(a -> mapslices(f, a; dims=dims), ds, patterns...)
end

"""
    merge(ds::KeyedDataset, others::KeyedDataset...)

Combine the constraints and data from multiple [`KeyedDataset`](@ref)s.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset;

julia> ds1 = KeyedDataset(
           :a => KeyedArray(zeros(3); time=1:3),
           :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
       );

julia> ds2 = KeyedDataset(
           :c => KeyedArray(ones(3); time=1:3),
           :d => KeyedArray(zeros(3, 2); time=1:3, loc=[:x, :y]),
       );

julia> collect(keys(merge(ds1, ds2).data))
4-element Vector{Tuple}:
 (:a,)
 (:b,)
 (:c,)
 (:d,)
```
"""
function Base.merge(ds::KeyedDataset, others::KeyedDataset...)
    return KeyedDataset(
        union(ds.constraints, getfield.(others, :constraints)...),
        merge(ds.data, getfield.(others, :data)...),
    )
end

"""
    rekey(f, ds, dim)

Apply function `f` to key values of each matching `dim` in the [`KeyedDataset`](@ref).
`dim` can either be a `Symbol` or a [`Pattern`](@ref) for the dimension paths.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset, rekey;

julia> ds = KeyedDataset(
           :a => KeyedArray(zeros(3); time=1:3),
           :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
       );

julia> r = rekey(k -> k .+ 1, ds, :time);

julia> r.time
3-element ReadOnlyArrays.ReadOnlyArray{Int64, 1, UnitRange{Int64}}:
 2
 3
 4
```
"""
rekey(f::Function, ds::KeyedDataset, args...) = rekey!(f, deepcopy(ds), args...)
rekey!(f::Function, ds::KeyedDataset, dim::Symbol) = rekey!(f, ds, Pattern(:__, dim))
rekey!(f::Function, ds::KeyedDataset, dimpath::Tuple) = rekey!(f, ds, Pattern(dimpath))
function rekey!(f::Function, ds::KeyedDataset, pattern::Pattern)
    for p in filter(in(pattern), dimpaths(ds))
        k, d = p[1:end-1], p[end]
        a = ds.data[k]

        names = dimnames(a)
        keys = map(zip(names, axiskeys(a))) do (n, vals)
            n === d ? f(vals) : vals
        end

        kw = NamedTuple{names}(keys)
        ds.data[k] = KeyedArray(parent(a); kw...)
    end
    validate(ds)
    return ds
end

# TODO: Implement cat, hcat and vcat


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    intersecttype(f, segment, triangle)

Calculate the intersection type of `segment and `triangle`
and apply function `f` to it.

## References

* Jiménez, J., Segura, R. and Feito, F. 2009. [A robust segment/triangle
  intersection algorithm for interference tests. Efficiency study]
  (https://www.sciencedirect.com/science/article/pii/S0925772109001448?via%3Dihub)
"""
function intersecttype(f::Function, s::Segment{3,T}, t::Triangle{3,T}) where {T}
  vₛ = vertices(s)
  vₜ = vertices(t)

  A = vₛ[1] - vₜ[3]
  B = vₜ[1] - vₜ[3]
  C = vₜ[2] - vₜ[3]

  W₁ = B × C
  wᵥ = A ⋅ W₁

  D = vₛ[2] - vₜ[3]
  sᵥ = D ⋅ W₁

  if wᵥ > atol(T)
    # rejection 2
    if sᵥ > atol(T)
      return NoIntersection() |> f
    end

    W₂ = A × D
    tᵥ = W₂ ⋅ C

    # rejection 3
    if tᵥ < -atol(T)
      return NoIntersection() |> f
    end

    uᵥ = -(W₂ ⋅ B)

    # rejection 4
    if uᵥ < -atol(T)
      return NoIntersection() |> f
    end

    # rejection 5
    if wᵥ < (sᵥ + tᵥ + uᵥ)
      return NoIntersection() |> f
    end
  elseif wᵥ < -atol(T)
    # rejection 2
    if sᵥ < -atol(T)
      return NoIntersection() |> f
    end

    W₂ = A × D
    tᵥ = W₂ ⋅ C

    # rejection 3
    if tᵥ > atol(T)
      return NoIntersection() |> f
    end

    uᵥ = -(W₂ ⋅ B)

    # rejection 4
    if uᵥ > atol(T)
      return NoIntersection() |> f
    end

    # rejection 5
    if wᵥ > (sᵥ + tᵥ + uᵥ)
      return NoIntersection() |> f
    end
  else
    if sᵥ > atol(T)
      W₂ = D × A
      tᵥ = W₂ ⋅ C

      # rejection 3
      if tᵥ < -atol(T)
        return NoIntersection() |> f
      end

      uᵥ = -(W₂ ⋅ B)

      # rejection 4
      if uᵥ < -atol(T)
        return NoIntersection() |> f
      end
      # rejection 5
      if -sᵥ < (tᵥ + uᵥ)
        return NoIntersection() |> f
      end
    elseif sᵥ < -atol(T)
      W₂ = D × A
      tᵥ = W₂ ⋅ C

      # rejection 3
      if tᵥ > atol(T)
        return NoIntersection() |> f
      end

      uᵥ = -(W₂ ⋅ B)

      # rejection 4
      if uᵥ > atol(T)
        return NoIntersection() |> f
      end

      # rejection 5
      if -sᵥ > (tᵥ + uᵥ)
        return NoIntersection() |> f
      end
    else
      # rejection 1, coplanar segment
      return NoIntersection() |> f
    end
  end

  λ = clamp(wᵥ / (wᵥ - sᵥ), zero(T), one(T))

  return IntersectingSegmentTriangle(s(λ)) |> f
end

"""
    intersecttype(s, t)

Calculate the intersection of a ray and triangle in 3D
and apply function `f` to it.

## References

* Möller, T. and Trumbore, B., 1997. [Fast, minimum storage ray-triangle
  intersection. Journal of graphics tools]
  (https://www.tandfonline.com/doi/abs/10.1080/10867651.1997.10487468)
"""
function intersecttype(f::Function, r::Ray{3,T}, t::Triangle{3,T}) where {T}
  vs = vertices(t)
  o = origin(r)
  d = direction(r)

  e₁ = vs[3] - vs[1]
  e₂ = vs[2] - vs[1]
  p = d × e₂
  det = e₁ ⋅ p

  # keep det > 0, modify T accordingly
  if det > atol(T)
    τ = o - vs[1]
  else
    τ = vs[1] - o
    det = -det
  end

  if det < atol(T)
    # This ray is parallel to the plane of the triangle.
    return NoIntersection() |> f
  end

  # calculate u parameter and test bounds
  u = τ ⋅ p
  if u < -atol(T) || u > det
    return NoIntersection() |> f
  end

  q = τ × e₁

  # calculate v parameter and test bounds
  v = d ⋅ q
  if v < -atol(T) || u + v > det
    return NoIntersection() |> f
  end

  λ = (e₂ ⋅ q) * (one(T) / det)

  if λ < -atol(T)
    return NoIntersection() |> f
  end

  λ = clamp(λ, zero(T), typemax(T))

  return IntersectingRayTriangle(r(λ)) |> f
end

"""
    intersecttype(s, p)

Calculate the intersection of a ray and polygon in 3D
and apply function `f` to it.
"""
function intersecttype(f::Function, r::Ray{3,T}, p::Ngon{N,3,T}) where {N,T}
  for t in triangulate(p)
    res = intersecttype(identity, r, t)
    if res != NoIntersection()
      return res |> f
    end
  end
  return NoIntersection() |> f
end


text:
#=
sino_geom.jl
sinogram geometry for 2D tomographic image reconstruction
2019-07-01, Jeff Fessler, University of Michigan
=#

export SinoGeom, sino_geom, sino_geom_par, sino_geom_fan, sino_geom_moj
export SinoPar, SinoFan, SinoMoj, SinoFanArc, SinoFanFlat
export sino_geom_help

abstract type SinoGeom end
abstract type SinoFan <: SinoGeom end
abstract type SinoParallel <: SinoGeom end

"`SinoPar` 2D parallel-beam sinogram geometry"
struct SinoPar <: SinoParallel
	units::Symbol			# :nothing | :mm | :cm etc.
	nb::Int					# # of "radial" samples, aka nr
	na::Int					# # of angular samples, aka nϕ
	d::Float32				# dr, "radial" sample spacing
	orbit::Float32			# [degrees]
	orbit_start::Float32	# [degrees]
	offset::Float32			# sample offset, cf offset_r or offset_s [unitless]
	strip_width::Float32	#
end

"`SinoMoj` 2D Mojette sinogram geometry"
struct SinoMoj <: SinoParallel
	units::Symbol			# :nothing | :mm | :cm etc.
	nb::Int					# # of "radial" samples, aka ns
	na::Int					# # of angular samples, aka nϕ
	d::Float32				# dx, pixels must be square
	orbit::Float32			# [degrees]
	orbit_start::Float32	# [degrees]
	offset::Float32			# sample offset, cf offset_r or offset_s [unitless]
	strip_width::Float32	#
end

"`SinoFanArc` 2D Fan-beam sinogram geometry for arc detector"
struct SinoFanArc <: SinoFan
	units::Symbol			# :nothing | :mm | :cm etc.
	nb::Int					# # of "radial" samples, aka ns
	na::Int					# # of angular samples, aka nβ
	d::Float32				# ds detector sample spacing
	orbit::Float32			# [degrees]
	orbit_start::Float32	# [degrees]
	offset::Float32			# sample offset, cf offset_r or offset_s [unitless]
	strip_width::Float32	#
	source_offset::Float32	# same units as d, etc., e.g., [mm]
							# use with caution!
	dsd::Float32			# dis_src_det, Inf for parallel beam
	dod::Float32			# dis_iso_det
#	dso::Float32			# dis_src_iso = dsd-dod, Inf for parallel beam
	dfs::Float32			# distance from focal spot to source
end

"`SinoFanFlat` 2D Fan-beam sinogram geometry for flat detector"
struct SinoFanFlat <: SinoFan
	units::Symbol			# :nothing | :mm | :cm etc.
	nb::Int					# # of "radial" samples, aka ns
	na::Int					# # of angular samples, aka nβ
	d::Float32				# ds detector sample spacing
	orbit::Float32			# [degrees]
	orbit_start::Float32	# [degrees]
	offset::Float32			# sample offset, cf offset_r or offset_s [unitless]
	strip_width::Float32	#
	source_offset::Float32	# same units as d, etc., e.g., [mm]
							# use with caution!
	dsd::Float32			# dis_src_det, Inf for parallel beam
	dod::Float32			# dis_iso_det
#	dso::Float32			# dis_src_iso = dsd-dod, Inf for parallel beam
	dfs::Float32			# distance from focal spot to source
end

"""
    sino_geom_help()
"""
function sino_geom_help( ; io::IO = isinteractive() ? stdout : devnull )
	print(io, "$(basename(@__FILE__)) propertynames:\n\t")
	print(io, Tuple(sort([propertynames(sino_geom(:par))...])))

	print(io,
	"\n
	Derived values

	sg.dim			dimensions: `(nb,na)`
	sg.ds|dr		radial sample spacing (NaN for :moj)
	sg.s			[nb] s sample locations
	sg.w			(nb-1)/2 + offset ('middle' sample position)
	sg.ad			source angles [degrees]
	sg.ar			source angles [radians]
	sg.ones			ones(Float32, nb,na)
	sg.zeros		zeros(Float32, nb,na)
	sg.rfov			radial fov
	sg.xds			[nb] center of detector elements (beta=0)
	sg.yds			[nb] ''
	sg.grid			(rg, phigrid) [nb na] parallel-beam coordinates
	sg.plot_grid(scatter)	plot `sg.grid` using `Plots.scatter`

	For mojette:

	sg.d_ang		[na]

	For fan beam:

	sg.gamma		[nb] gamma sample values [radians]
	sg.gamma_max	half of fan angle [radians]
	sg.dso			# dsd - dod, Inf for parallel beam

	Methods

	sg.down(down)		reduce sampling by integer factor
	sg.shape(sino)		reshape sinograms into array [nb na :]
	sg.unitv(;ib,ia)	unit 'vector' with single nonzero element
	sg.taufun(x,y)		projected s/ds for each (x,y) pair [numel(x) na]
	sg.plot!(plot!;ig)	plot system geometry (mostly for SinoFan)
	\n")
end

"""
    function sg = sino_geom(...)

Constructor for `SinoGeom`

Create the "sinogram geometry" structure that describes the sampling
characteristics of a given sinogram for a 2D parallel or fan-beam system.
Using this structure facilitates "object oriented" code.
(Use `ct_geom()` instead for 3D axial or helical cone-beam CT.)

in
- `how::Symbol`	`:fan` (fan-beam) | `:par` (parallel-beam) | `:moj` (mojette)

options for all geometries (including parallel-beam):
- `units::Symbol`	e.g. `:cm` or `:mm`; default: :none
- `orbit_start`		default: 0
- `orbit`			[degrees] default: `180` for parallel / mojette
and `360` for fan
   * can be `:short` for fan-beam short scan
- `down::Int`		down-sampling factor, for testing

- `nb`				# radial samples cf `nr` (i.e., `ns` for `:fan`)
- `na`				# angular samples (cf `nbeta` for `:fan`)
- `d`				radial sample spacing; cf `dr` or `ds`; default 1
  * for mojette this is actually `dx`
- `offset`			cf `offset_r` `channel_offset` unitless; default 0
   * (relative to centerline between two central channels).
   * Use 0.25 or 1.25 for "quarter-detector offset"
- `strip_width`		detector width; default: `d`

options for fan-beam
- `source_offset`		same units as d; use with caution! default 0
fan beam distances:
- `dsd` cf `dis_src_det` default: `Inf` (parallel beam)
- `dod` cf `dis_iso_det` default: `0`
- `dfs` cf `dis_foc_src` default: `0` (3rd generation CT arc),
   * use `Inf` for flat detector

out
- `sg::SinoGeom`	initialized structure

See also
- `sino_geom_help()` help on methods
- `sino_geom_plot_grids()` show sampling

Jeff Fessler, University of Michigan
"""
function sino_geom(how::Symbol ; kwarg...)
	if how === :par
		sg = sino_geom_par( ; kwarg...)
	elseif how === :fan
		sg = sino_geom_fan( ; kwarg...)
	elseif how === :moj
		sg = sino_geom_moj( ; kwarg...)
	elseif how === :ge1
		sg = sino_geom_ge1( ; kwarg...)
#=
	elseif how === :hd1
		sg = sino_geom_hd1( ; kwarg...)
	elseif how === :revo1fan
		tmp = ir_fan_geom_revo1(type)
		sg = sino_geom(:fan, tmp{:}, varargin{:})
=#
	else
		throw("unknown sino type $how")
	end

	return sg
end

SinoFanMaker(dfs::Real) =
	(dfs == 0) ? SinoFanArc :
	isinf(dfs) ? SinoFanFlat :
	throw("dfs $dfs") # must be 0 or Inf

# common to all
function _downsample(sg::SinoGeom, down::Int)
	nb = 2 * (sg.nb ÷ 2down) # keep it even
	na = sg.na ÷ down
	return (sg.units,
		nb, na, sg.d * down, sg.orbit, sg.orbit_start, sg.offset,
		sg.strip_width * down)
end

"""
    sg = downsample(sg, down)
down-sample (for testing with small problems)
"""
function downsample(sg::T, down::Int) where {T <: SinoParallel}
	return (down == 1) ? sg : T(_downsample(sg, down)...)
end
function downsample(sg::SinoFan, down::Int)
	down == 1 && return sg
	return SinoFanMaker(sg.dfs)(_downsample(sg, down)...,
		sg.source_offset, sg.dsd, sg.dod, sg.dfs)
end

# common to all
function _sino_geom_over(sg::SinoGeom, over::Int)
	return (sg.units,
		sg.nb * over, sg.na, sg.d / over,
		sg.orbit, sg.orbit_start, sg.offset * over,
		sg.strip_width / over)
end

"""
    sg = sino_geom_over(sg, over::Int)
over-sample in "radial" dimension
For Mojette sampling, it means that `d = dx/over`.
"""
function sino_geom_over(sg::T, over::Int) where {T <: SinoParallel}
	return (over == 1) ? sg : T(_sino_geom_over(sg, over)...)
end
function sino_geom_over(sg::SinoFan, over::Int)
	return (over == 1) ? sg :
		SinoFanMaker(sg.dfs)(_sino_geom_over(sg, over)...,
			sg.source_offset, sg.dsd, sg.dod, sg.dfs)
end

"""
    sg = sino_geom_fan()
"""
function sino_geom_fan( ;
	units::Symbol = :none,
	nb::Int = 128,
	na::Int = 2 * floor(Int, nb * π/2 / 2),
	d::Real = 1,
	orbit::Union{Symbol,Real} = 360, # [degrees]
	orbit_start::Real = 0,
	strip_width::Real = d,
	offset::Real = 0,
	source_offset::Real = 0,
	dsd::Real = 4*nb*d,	# dis_src_det
#	dso::Real = [],		# dis_src_iso
	dod::Real = nb*d,	# dis_iso_det
	dfs::Real = 0,		# dis_foc_src (3rd gen CT)
	down::Int = 1,
)

	maker = (dfs == 0) ? SinoFanArc :
			isinf(dfs) ? SinoFanFlat :
			throw("dfs $dfs") # must be 0 or Inf

	if orbit === :short # trick
		sg_tmp = maker(units,
			nb, na, d, 0, orbit_start, offset, strip_width,
			source_offset, dsd, dod, dfs)
		orbit = sg_tmp.orbit_short
	end
	isa(orbit, Symbol) && throw("orbit :orbit")

	sg = maker(units,
		nb, na, d, orbit, orbit_start, offset, strip_width,
		source_offset, dsd, dod, dfs)

	return downsample(sg, down)
end

"""
    sg = sino_geom_par( ... )
"""
function sino_geom_par( ;
	units::Symbol = :none,
	nb::Int = 128,
	na::Int = 2 * floor(Int, nb * π/2 / 2),
	down::Int = 1,
	d::Real = 1,
	orbit::Real = 180, # [degrees]
	orbit_start::Real = 0,
	strip_width::Real = d,
	offset::Real = 0,
)
	sg = SinoPar(units,
		nb, na, d, orbit, orbit_start, offset, strip_width)
	return downsample(sg, down)
end

"""
    sg = sino_geom_moj( ... )
"""
function sino_geom_moj( ;
	units::Symbol = :none,
	nb::Int = 128,
	na::Int = 2 * floor(Int, nb * π/2 / 2),
	down::Int = 1,
	d::Real = 1, # means dx for :moj
	orbit::Real = 180, # [degrees]
	orbit_start::Real = 0,
	strip_width::Real = d, # ignored ?
	offset::Real = 0,
)
	sg = SinoMoj(units,
		nb, na, d, orbit, orbit_start, offset, strip_width)
	return downsample(sg, down)
end

"gamma for general finite dfs (rarely used)"
function sino_geom_gamma_dfs(sg::SinoFan)
	dis_foc_det = sg.dfs + sg.dsd
	alf = sg.s / dis_foc_det
	atan.(dis_foc_det * sin.(alf), dis_foc_det * cos.(alf) .- sg.dfs)
	# equivalent to s/dsd when dfs=0
end

"""
    sino_geom_gamma()
gamma sample values for :fan
"""
function sino_geom_gamma(sg::SinoFan)
	return	sg.dfs == 0 ? sg.s / sg.dsd : # 3rd gen: equiangular
			isinf(sg.dfs) ? atan.(sg.s / sg.dsd) : # flat
			sino_geom_gamma_dfs(sg) # general
end

"""
    sino_geom_rfov()
radial FOV
"""
sino_geom_rfov(sg::SinoPar) = maximum(abs.(sg.r))
sino_geom_rfov(sg::SinoFan) = sg.dso * sin(sg.gamma_max)
sino_geom_rfov(sg::SinoMoj) = sg.nb/2 * minimum(sg.d_ang) # (ignores offset)

function _sino_geom_taufun(sg::SinoParallel, x, y)
	ar = sg.ar' # row vector, for outer-product
	return (x * cos.(ar) + y * sin.(ar)) / sg.dr # tau
end

function _sino_geom_taufun(sg::SinoFan, x, y)
	b = sg.ar' # row vector, for outer-product
	xb = x * cos.(b) + y * sin.(b)
	yb = -x * sin.(b) + y * cos.(b)
	tangam = (xb .- sg.source_offset) ./ (sg.dso .- yb) # e,tomo,fan,L,gam
	if sg.dfs == 0 # arc
		tau = sg.dsd / sg.ds * atan.(tangam)
	elseif isinf(sg.dfs) # flat
		tau = sg.dsd / sg.ds * tangam
#	else
#		throw("bad dfs $(sg.dfs)")
	end
	return tau
end

"""
    sino_geom_taufun()
projected `s/ds`, useful for footprint center and support
"""
function sino_geom_taufun(sg::SinoGeom, x, y)
	size(x) != size(y) && throw("bad x,y size")
	return _sino_geom_taufun(sg, vec(x), vec(y))
end

"""
    sino_geom_xds()
center positions of detectors (for beta = 0)
"""
sino_geom_xds(sg::SinoPar) = sg.s
sino_geom_xds(sg::SinoMoj) = sg.s # todo: really should be angle dependent
sino_geom_xds(sg::SinoFanArc) = sg.dsd * sin.( sg.gamma) .+ sg.source_offset
sino_geom_xds(sg::SinoFanFlat) = sg.s .+ sg.source_offset

"""
    sino_geom_yds()
center positions of detectors (for beta = 0)
"""
sino_geom_yds(sg::SinoPar) = zeros(Float32, sg.nb)
sino_geom_yds(sg::SinoMoj) = zeros(Float32, sg.nb)
sino_geom_yds(sg::SinoFanArc) = sg.dso .- sg.dsd * cos.(sg.gamma)
sino_geom_yds(sg::SinoFanFlat) = fill(-sg.dod, sg.nb)

"""
    sino_geom_unitv()
sinogram with a single ray
"""
function sino_geom_unitv(
	sg::SinoGeom ;
	ib::Int = sg.nb ÷ 2 + 1,
	ia::Int = sg.na ÷ 2 + 1,
)
	out = sg.zeros
	out[ib,ia] = 1
	return out
end

"""
    (rg, ϕg) = sino_geom_grid(sg::SinoGeom)

Return grids `rg` and `ϕg` (in radians) of size `[nb na]`
of equivalent *parallel-beam* `(r,ϕ)` (radial, angular) sampling positions,
for any sinogram geometry.
For parallel beam this is just `ndgrid(sg.r, sg.ar)`
but for fan beam and mojette this involves more complicated computations.
"""
sino_geom_grid(sg::SinoPar) = ndgrid(sg.r, sg.ar)

function sino_geom_grid(sg::SinoFan)
	gamma = sg.gamma
	rad = sg.dso * sin.(gamma) + sg.source_offset * cos.(gamma)
	rg = repeat(rad, 1, sg.na) # [nb na]
	return (rg, gamma .+ sg.ar') # [nb na] phi = gamma + beta
end

function sino_geom_grid(sg::SinoMoj)
	phi = sg.ar
	# trick: ray_spacing aka ds comes from dx which is sg.d for mojette
	wb = (sg.nb - 1)/2 + sg.offset
	dt = sg.d_ang # [na]
	pos = ((0:(sg.nb-1)) .- wb) * dt' # [nb na]
	return (pos, repeat(phi', sg.nb, 1))
end

"""
    show(io::IO, sg::SinoGeom)
    show(io::IO, ::MIME"text/plain", sg::SinoGeom)
"""
Base.show(io::IO, sg::SinoGeom) =
    print(io, "$(typeof(sg)): $(sg.dim)")
function Base.show(io::IO, ::MIME"text/plain", sg::SinoGeom)
    ir_dump(io, sg)
end

#sino_geom_how(sg::SinoGeom{G}) where {G} = G

# Extended properties

sino_geom_fun0 = Dict([
	(:help, sg -> sino_geom_help()),

	(:dim, sg -> (sg.nb, sg.na)),
	(:w, sg -> (sg.nb-1)/2 + sg.offset),
	(:ones, sg -> ones(Float32, sg.dim)),
	(:zeros, sg -> zeros(Float32, sg.dim)),

	(:dr, sg -> ((sg isa SinoMoj) ? NaN : sg.d)),
	(:ds, sg -> sg.dr),
	(:r, sg -> sg.d * ((0:sg.nb-1) .- sg.w)),
	(:s, sg -> sg.r), # sample locations ('radial')

	(:gamma, sg -> sino_geom_gamma(sg)),
	(:gamma_max, sg -> maximum(abs.(sg.gamma))),
	(:orbit_short, sg -> 180 + 2 * rad2deg(sg.gamma_max)),
	(:ad, sg -> (0:sg.na-1)/sg.na * sg.orbit .+ sg.orbit_start),
	(:ar, sg -> deg2rad.(sg.ad)),

	(:rfov, sg -> sino_geom_rfov(sg)),
	(:xds, sg -> sino_geom_xds(sg)),
	(:yds, sg -> sino_geom_yds(sg)),
	(:dso, sg -> sg.dsd - sg.dod),
	(:grid, sg -> sino_geom_grid(sg)),
	(:plot_grid, sg -> ((plot::Function) -> sino_geom_plot_grid(sg, plot))),

	(:plot!, sg ->
		((plot!::Function ; ig=nothing) -> sino_geom_plot!(sg, plot! ; ig=ig))),
	(:shape, sg -> (((x::AbstractArray) -> reshaper(x, sg.dim)))),
	(:taufun, sg -> ((x,y) -> sino_geom_taufun(sg,x,y))),
	(:unitv, sg -> ((;kwarg...) -> sino_geom_unitv(sg; kwarg...))),

	# angular dependent d for :moj
	(:d_moj, sg -> (ar -> sg.d * max(abs(cos(ar)), abs(sin(ar))))),
	(:d_ang, sg -> sg.d_moj.(sg.ar)),

	# functions that return new geometry:

	(:down, sg -> (down::Int -> downsample(sg, down))),
	(:over, sg -> (over::Int -> sino_geom_over(sg, over))),

	])

# Tricky overloading here!

Base.getproperty(sg::SinoGeom, s::Symbol) =
		haskey(sino_geom_fun0, s) ? sino_geom_fun0[s](sg) :
		getfield(sg, s)

Base.propertynames(sg::SinoGeom) =
	(fieldnames(typeof(sg))..., keys(sino_geom_fun0)...)

"""
    reshaper(x::AbstractArray, dim:Dims)
Reshape `x` to size `dim` with `:` only if needed
"""
reshaper(x::AbstractArray, dim::Dims) =
	(length(x) == prod(dim)) ? reshape(x, dim) : reshape(x, dim..., :)

"""
    sino_geom_ge1()
sinogram geometry for GE lightspeed system
These numbers are published in IEEE T-MI Oct. 2006, p.1272-1283 wang:06:pwl
"""
function sino_geom_ge1( ;
	na::Int = 984,
	nb::Int = 888,
	orbit::Union{Symbol,Real} = 360,
	units::Symbol = :mm, # default units is mm
	kwarg...,
)

	if orbit === :short
		na = 642 # trick: reduce na for short scans
		orbit = na/984*360
	end

	scale = units === :mm ? 1 :
			units === :cm ? 10 :
			throw("units $units")

	return sino_geom_fan( ; units=units,
		nb=nb, na=na,
		d = 1.0239/scale, offset = 1.25,
		dsd = 949.075/scale,
		dod = 408.075/scale,
		dfs = 0, kwarg...,
	)
end


text:
# Test collection_sets function
@testset "collection_sets" begin
    set1 = IntervalSet(0, 1)
    set2 = UniDistributionSet(Uniform())
    set = CollectionSet([set1, set2])
    @test_throws ErrorException collection_sets(set1)
    @test collection_sets(set) == [set1, set2]
end

# Test length
@testset "Base.length" begin
    sset1 = IntervalSet(0, 1)
    sset2 = UniDistributionSet(Uniform())
    set1 = CollectionSet([sset1, sset2])
    set2 = MultiDistributionSet(MvNormal(ones(4)))
    @test length(sset1) == 1
    @test length(sset2) == 1
    @test length(set1) == 2
    @test length(set2) == 4
end

# Test supports_in_set functions
@testset "supports_in_set" begin
    # supports_in_set (IntervalSet)
    @testset "IntervalSet" begin
        set = IntervalSet(0, 1)
        @test supports_in_set(0, set)
        @test !supports_in_set(-1, set)
        @test !supports_in_set(2, set)
    end
    # supports_in_set (UnivariateDistribution set)
    @testset "Univariate Distribution" begin
        set = UniDistributionSet(Uniform())
        @test supports_in_set(0, set)
        @test !supports_in_set(-1, set)
        @test !supports_in_set(2, set)
    end
    # supports_in_set (MultivariateDistribution set)
    @testset "Multivariate Distribution" begin
        set = MultiDistributionSet(MvNormal(ones(2)))
        @test supports_in_set(ones(2), set)
        @test supports_in_set(ones(2, 10), set)
        bad_supports = [1 1; 2 2; 3 3];
        @test_throws ErrorException supports_in_set(bad_supports, set)
    end
    # supports_in_set (CollectionSet)
    @testset "CollectionSet" begin
        set1 = IntervalSet(0, 1)
        set2 = UniDistributionSet(Uniform())
        set = CollectionSet([set1, set2])
        bad_supports = [1 1; 2 2; 3 3];
        @test_throws ErrorException supports_in_set(bad_supports, set)
        supports = [1 2; 3 4];
        @test !supports_in_set(supports, set)
        supports = [0 1; 0 1];
        @test supports_in_set(supports, set)
    end
    # supports_in_set (Fallback)
    @testset "Fallback" begin
        @test supports_in_set(0, BadSet())
    end
end

# Test lower bound functions
@testset "Lower Bound" begin
    set1 = IntervalSet(0, 1)
    set2 = UniDistributionSet(Uniform())
    set3 = MultiDistributionSet(MvNormal([0., 0.], [1. 0.; 0. 2.]))
    set4 = CollectionSet([set1, set2])
    set5 = BadSet()
    # JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(set1)
        @test has_lower_bound(set2)
        @test !has_lower_bound(set3)
        @test has_lower_bound(set4)
        @test !has_lower_bound(set5)
    end
    # JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(set1) == 0.0
        @test lower_bound(set2) == 0.0
        @test_throws ErrorException lower_bound(set3)
        @test lower_bound(set4) == [0.0, 0.0]
        @test_throws ErrorException lower_bound(set5)
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test set_lower_bound(set1, 0.5) == IntervalSet(0.5, 1)
        @test_throws ErrorException set_lower_bound(set2, 2)
        @test_throws ErrorException set_lower_bound(set3, 2)
        @test_throws ErrorException set_lower_bound(set4, [0.1, 0.1])
        @test_throws ErrorException set_lower_bound(set5, 2)
        set6 = CollectionSet([set1, set1])
        @test set_lower_bound(set6, [0.1, 0.1]).sets ==
                                    [IntervalSet(0.1, 1.), IntervalSet(0.1, 1.)]
    end
end

# Test upper bound functions
@testset "Upper Bound" begin
    set1 = IntervalSet(0, 1)
    set2 = UniDistributionSet(Uniform())
    set3 = MultiDistributionSet(MvNormal([0., 0.], [1. 0.; 0. 2.]))
    set4 = CollectionSet([set1, set2])
    set5 = BadSet()
    # JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(set1)
        @test has_upper_bound(set2)
        @test !has_upper_bound(set3)
        @test has_upper_bound(set4)
        @test !has_upper_bound(set5)
    end
    # JuMP.upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(set1) == 1.0
        @test upper_bound(set2) == 1.0
        @test_throws ErrorException upper_bound(set3)
        @test upper_bound(set4) == [1.0, 1.0]
        @test_throws ErrorException upper_bound(set5)
    end
    # JuMP.set_upper_bound
    @testset "JuMP.set_upper_bound" begin
        @test set_upper_bound(set1, 0.5) == IntervalSet(0, 0.5)
        @test_throws ErrorException set_upper_bound(set2, 2)
        @test_throws ErrorException set_upper_bound(set3, 2)
        @test_throws ErrorException set_upper_bound(set4, [0.1, 0.1])
        @test_throws ErrorException set_upper_bound(set5, 2)
        set6 = CollectionSet([set1, set1])
        @test set_upper_bound(set6, [0.9, 0.9]).sets ==
                                    [IntervalSet(0., 0.9), IntervalSet(0., 0.9)]
    end
end

# Test support labels
@testset "Support Labels" begin
    @testset "DataTypes" begin
        @test AbstractSupportLabel isa DataType
        @test All <: AbstractSupportLabel
        @test PublicLabel <: AbstractSupportLabel
        @test UserDefined <: PublicLabel
        @test UniformGrid <: PublicLabel
        @test SampleLabel <: PublicLabel
        @test MCSample <: SampleLabel
        @test WeightedSample <: SampleLabel
        @test Mixture <: PublicLabel
        @test UniqueMeasure isa UnionAll
        @test InternalLabel <: AbstractSupportLabel
        @test OrthogonalCollocationNode <: InternalLabel
    end
    @testset "generate_unique_label" begin
        @test InfiniteOpt.generate_unique_label() <: UniqueMeasure
    end
end

# Test support generation
@testset "generate_support_values" begin
    @testset "IntervalSet" begin
        set = IntervalSet(0., 1.)
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1] isa Vector{<:Number}
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[2] == UniformGrid
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1][2] == 0.111
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1][2] != 1/11
        @test length(generate_support_values(set, num_supports = 10, sig_digits = 3)[1]) == 10
        @test generate_support_values(set, MCSample, num_supports = 10)[1] isa Vector{<:Number}
        @test generate_support_values(set, MCSample, num_supports = 10)[2] == MCSample
        @test_throws ErrorException generate_support_values(set, Val(:a))
        @test generate_support_values(set, All)[2] == UniformGrid
    end
    @testset "Distribution Sets" begin
        dist1 = Normal(0., 1.)
        dist2 = MvNormal([0.; 0.], [1. 0.; 0. 2.])
        set1 = UniDistributionSet(dist1)
        set2 = MultiDistributionSet(dist2)
        @test generate_support_values(set1, num_supports = 10)[1] isa Vector{<:Number}
        @test generate_support_values(set2, num_supports = 10)[1] isa Array{<:Number, 2}
        @test generate_support_values(set2, num_supports = 10)[2] == WeightedSample
        @test length(generate_support_values(set1, num_supports = 10)[1]) == 10
        @test size(generate_support_values(set2, num_supports = 10)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(set1, Val(:a))
        @test_throws ErrorException generate_support_values(set2, Val(:a))
        @test generate_support_values(set1, All)[2] == WeightedSample
        @test generate_support_values(set2, All)[2] == WeightedSample
        @test generate_support_values(set1, MCSample)[2] == MCSample
    end
    @testset "Matrix Distribution Sets" begin
        dist = MatrixBeta(2, 2, 2)
        set = MultiDistributionSet(dist)
        @test generate_support_values(set, num_supports = 10)[1] isa Array{<:Number, 2}
        @test generate_support_values(set, num_supports = 10)[2] == WeightedSample
        @test size(generate_support_values(set, num_supports = 10)[1]) == (4, 10)
        @test_throws ErrorException generate_support_values(set, Val(:a))
        @test generate_support_values(set, All)[2] == WeightedSample
    end
    @testset "_generate_collection_supports" begin
        set1 = IntervalSet(0., 1.)
        set2 = IntervalSet(0., 1.)
        set = CollectionSet([set1, set2])
        @test InfiniteOpt._generate_collection_supports(set, 10, 3) isa Array{<:Number, 2}
        @test InfiniteOpt._generate_collection_supports(set, 10, 3)[2, 2] == 0.111
        @test InfiniteOpt._generate_collection_supports(set, 10, 3)[2, 2] != 1/11
        @test size(InfiniteOpt._generate_collection_supports(set, 10, 3)) == (2, 10)
    end
    @testset "CollectionSet (IntervalSets)" begin
        set1 = IntervalSet(0., 1.)
        set2 = IntervalSet(0., 1.)
        set = CollectionSet([set1, set2])
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[2] == UniformGrid
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1][2, 2] == 0.111
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1][2, 2] != 1/11
        @test size(generate_support_values(set, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(set, Val(:a))
        @test generate_support_values(set, All)[2] == UniformGrid
        @test generate_support_values(set, MCSample)[2] == MCSample
        @test generate_support_values(set, MCSample)[1] isa Array{<:Number, 2}
    end
    @testset "CollectionSet (UniDistributionSets)" begin
        set1 = UniDistributionSet(Normal())
        set2 = UniDistributionSet(Normal())
        set = CollectionSet([set1, set2])
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[2] == WeightedSample
        @test size(generate_support_values(set, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(set, Val(:a))
        @test generate_support_values(set, All)[2] == WeightedSample
    end
    @testset "CollectionSet (InfiniteScalarSets)" begin
        set1 = UniDistributionSet(Normal())
        set2 = IntervalSet(0., 1.)
        set = CollectionSet([set1, set2])
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(set, num_supports = 10, sig_digits = 3)[2] == Mixture
        @test generate_support_values(CollectionSet([set2, set2]), MCSample, num_supports = 10, sig_digits = 3)[2] == MCSample
        @test size(generate_support_values(set, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(set, Val(:a))
        @test generate_support_values(set, All)[2] == Mixture
        @test generate_support_values(set, MCSample)[2] == MCSample
        @test generate_support_values(set, MCSample)[1] isa Array{<:Number, 2}
    end
    @testset "Fallback" begin
        @test_throws ErrorException generate_support_values(BadSet())
    end
    @testset "User Interface" begin
        set = IntervalSet(0., 1.)
        @test generate_supports(set)[1] isa Vector{<:Number}
        @test generate_supports(set, MCSample)[1] isa Vector{<:Number}
    end
end


text:
"""
makes a soft analogue of the heaviside step function. useful for inputs to differential equations, as it's easier on the numerics.
"""
soft_heaviside(t, nastiness, step_time) = 1 / (1 + exp(nastiness * (step_time - t)))
soft_heaviside(nastiness, step_time) = t -> soft_heaviside(t, nastiness, step_time)

get_ids_names(opArray) = repr.(opArray)

"""
    l2_hessian(nom_sol)
gets hessian according to L2 loss under the assumption that loss(θ₀) = 0. nom_sol is the solution of the nominal ODEProblem.
The Hessian then only requires first derivatives: it is sum_ij dyi/dθ * dyj/dtheta
"""
function l2_hessian(nom_sol)
    prob = nom_sol.prob
    function pToL2(p)
        pprob = remake(prob, p=p)
        psol = solve(pprob, nom_sol.alg, saveat=nom_sol.t) |> Array
        psol = reshape(psol, 1, :)
        return psol
    end

    gr = ForwardDiff.jacobian(pToL2, prob.p)
    u, d, v = svd(gr)
    return v * diagm(d.^2) * v'
    # = gr'*gr but a bit more accurate
end


text:
function test_pi_chudnovsky()
    simple  = Brainstorm._Math._Series._Pi._Chudnovsky._Simple.pi
    inplace = Brainstorm._Math._Series._Pi._Chudnovsky._Inplace.pi

    setprecision(4_000_000) do
        expected = [trunc(BigInt, pi*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)]

        @test all([simple(d)  for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] - expected .== 0)
        @test all([inplace(d) for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] - expected .== 0)
    end
end

function test_pi_machin()
    euler           = Brainstorm._Math._Series._Pi._Machin.euler
    hermann         = Brainstorm._Math._Series._Pi._Machin.hermann
    hutton          = Brainstorm._Math._Series._Pi._Machin.hutton
    machin          = Brainstorm._Math._Series._Pi._Machin.machin

    gauss           = Brainstorm._Math._Series._Pi._Machin.gauss
    strassnitzky    = Brainstorm._Math._Series._Pi._Machin.strassnitzky
    klingenstierna  = Brainstorm._Math._Series._Pi._Machin.klingenstierna
    euler3          = Brainstorm._Math._Series._Pi._Machin.euler3
    loney           = Brainstorm._Math._Series._Pi._Machin.loney
    stoermer        = Brainstorm._Math._Series._Pi._Machin.stoermer

    takano1982      = Brainstorm._Math._Series._Pi._Machin.takano1982
    stoermer1896    = Brainstorm._Math._Series._Pi._Machin.stoermer1896

    hwang1997       = Brainstorm._Math._Series._Pi._Machin.hwang1997
    hwang2003       = Brainstorm._Math._Series._Pi._Machin.hwang2003
    wetherfield2004 = Brainstorm._Math._Series._Pi._Machin.wetherfield2004

    setprecision(40_000) do
        expected = [trunc(BigInt, pi*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000)]

        @test all([euler(d)           for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([hermann(d)         for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([hutton(d)          for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([machin(d)          for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)

        @test all([gauss(d)           for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([strassnitzky(d)    for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([klingenstierna(d)  for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([euler3(d)          for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([loney(d)           for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([stoermer(d)        for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)

        @test all([takano1982(d)      for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([stoermer1896(d)    for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)

        @test all([hwang1997(d)       for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([hwang2003(d)       for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
        @test all([wetherfield2004(d) for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
    end
end

function test_pi_gausslegendre()
    int_simple    = Brainstorm._Math._Series._Pi._Iterative._GaussLegendre._Int._Simple.pi
    int_inplace   = Brainstorm._Math._Series._Pi._Iterative._GaussLegendre._Int._Inplace.pi
    float_simple  = Brainstorm._Math._Series._Pi._Iterative._GaussLegendre._Float._Simple.pi
    float_inplace = Brainstorm._Math._Series._Pi._Iterative._GaussLegendre._Float._Inplace.pi

    setprecision(400_000) do
        expected = [trunc(BigInt, pi*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000, 100_000)]

        @test all([int_simple(d)                               for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([int_inplace(d)                              for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([trunc(BigInt, float_simple(d)*(big(10)^d))  for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([trunc(BigInt, float_inplace(d)*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
    end
end

function test_pi_sine()
    sine = Brainstorm._Math._Series._Pi._Iterative._Sine.pi

    setprecision(40_000) do
        expected = [trunc(BigInt, pi*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000)]

        @test all([trunc(BigInt, sine(d)*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000)] - expected .== 0)
    end
end

function test_pi_all()
    print(rpad("Math.Series.Pi...", 50, ' '))

    test_pi_chudnovsky()
    test_pi_machin()
    test_pi_gausslegendre()
    test_pi_sine()

    println("PASS")
end


text:
#%%
using Revise
using GdMn2O5
const Gd = GdMn2O5
using LinearAlgebra
using LaTeXStrings
#%%
l = Ledger(physical_model(J1=0.7));
function find_switching_field(hsweep)
    id = findfirst(x -> abs(hsweep[x+1][L][1].ϕ - hsweep[x][L][1].ϕ) > π/4, 1:length(hsweep)-1)
    if id === nothing
        error("Switching field could not be found.")
    else
        return hsweep[id][H][1]
    end
end

start = 0
stop = 10.0
len = 200
Hr = H_sweep_range(0, 2.5, 120)
l[E_field][1].e = -0.05
hsweep1 = generate_Hsweep_states(l, 10.0, Hr)
H1 = find_switching_field(hsweep1)
l[E_field][1].e = 0.05
hsweep2 = generate_Hsweep_states(l, 10.0, Hr)
H2 = find_switching_field(hsweep2)
l[E_field][1].e = 0.00
hsweep3 = generate_Hsweep_states(l, 10.0, Hr)
H3 = find_switching_field(hsweep3)

Hr1 = [range(0, H2.magnitude+0.02, length=120); range(H2.magnitude+0.02, 0, length=120)[2:end]; H_sweep_range(0, 2.5, 120)]

l[E_field][1].e = -0.05
hsweep1 = generate_Hsweep_states(l, 10.0, Hr1)
l[E_field][1].e = 0.05
hsweep2 = generate_Hsweep_states(l, 10.0, Hr1)
l[E_field][1].e = 0.00
hsweep3 = generate_Hsweep_states(l, 10.0, Hr1)
plot(Hr1, 5 .* Pb.(hsweep1) .+ 5 .* Pb.(hsweep2) .+ Pb.(hsweep3), color=Gd.get_colors(ColorSchemes.rainbow, Hr1))
plot!(Hr1, )
plot!(Hr1, Pb.(hsweep3))
#%%
Gd.save_movie(dio[Gd.Movie][1]);
Gd.load_movie("1");

cur = sim_snapshot(dio);

cur == l

cur[H]
@time set_sim!(dio, cur);
#%%

l = Ledger(physical_model())
Entity(l, GdMn2O5.VisualizationSettings())
l[H].data[1] = H(ϕ=deg2rad(11))
dio = Diorama(l);
expose(dio);

#%%
# l = Ledger(physical_model(χ_L = 2/20.0, K_Gd=0.40));
l = Ledger(physical_model(J1=0.7));
p = plot(Hr, Pb.(hsweep))
#%%
plot!(angles, winding_numbers, label="$v")
plot(p, legendtitle=L"K_{L}", dpi=200, yguide = "winding number", xguide = L"\phi_H", legendfontsize=10)
#%%
physical_model().J1
function calc_magnetization_diff(l)
    diffs = map(i->l[Spin].data[i].v - l[Spin].data[i-2].v, 3:2:8)
    push!(diffs, l[Spin].data[7].v - l[Spin].data[1].v)
    return diffs
end
alldiffs = map(l -> calc_magnetization_diff(l), hsweep)
sum(alldiffs)
t = Vec3[]
for i = 2:length(alldiffs)
    push!(t, alldiffs[i][2] - alldiffs[i-1][2])
end
plot(getindex.(t,2))
plot!(getindex.(alldiffs,2))
plot!(getindex.(alldiffs,3))
plot!(getindex.(alldiffs,4))
#%%
t = Ledger(physical_model())

Erange = H_sweep_range(0.2, -0.2, len)
Esweep = generate_Esweep_states(hsweep[75], Erange)

tl = copy(hsweep[75])

plot(H_sweep_range(0.2, -0.2, len), map(x->(tl[E_field][1].e = x;optimize!(tl, Spin, L, Pb); Pb(tl)), H_sweep_range(0.2, -0.2, len)))
plot(Hr, map(x->x[Pb][1].pb, hsweep))
plot(Erange, map(x->x[L][1].ϕ, Esweep), color=GdMn2O5.get_colors(ColorSchemes.rainbow, Esweep))
l[Etot][1].e
update(GdMn2O5.E_LL(), l)

#%%
#%%
Model().dipole_strength
plot(Hr, map(x->x[L][1].ϕ, hsweep), color=GdMn2O5.get_colors(ColorSchemes.rainbow, hsweep))
components(dio)
#%%
testsweep = deepcopy(hsweep)
for t in testsweep
    t[H].data[1] = GdMn2O5.H(magnitude=6.4, ϕ=deg2rad(12))
end
test = GdMn2O5.neb(testsweep, 1.0)
test = GdMn2O5.neb(test)
plot(Hr, Pb.(test))

plot(map(x->x.e, calc_etot.(test)))
findfirst(x->x.i1 == l[Bond][1].i2 && x.i2 == l[Bond][1].i1 && x.r == -l[Bond][1].r, l[Bond].data)

#%%

dio = Diorama(l)
expose(dio)
#%%J_perpendicular(best_model())
states = find_states(l)
state_ledgers = Ledger[]

for s in states
    state_ledger = deepcopy(l)
    set_angles!(state_ledger, s)
    push!(state_ledgers, state_ledger)
end
Hr = repeat([range(0.0,10.0, length=100);range(10.0, 0.0, length=100)], 2)
# h_sweep_states0 = map(state_ledgers) do sl
#     generate_Hsweep_states(sl, 6, Hr)
# end

# h_sweep_states11 = map(state_ledgers) do sl
#     generate_Hsweep_states(sl, 11, Hr)
# end
# h_sweep_states20 = map(state_ledgers) do sl
#     generate_Hsweep_states(sl, 20, Hr)
# end
h_sweep_state0 = generate_Hsweep_states(l, 6, Hr)

h_sweep_state11 = generate_Hsweep_states(l, 11, Hr)
h_sweep_state20 = generate_Hsweep_states(l, 20, Hr)

# p3 = plot(Hr, Pb(h_sweep_states0[1]), ylabel="Pb", xlabel="Hr", title="0 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
# p4 = plot(Hr, Pb(h_sweep_states11[1]), ylabel="Pb", xlabel="Hr", title="11 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
# p5 = plot(Hr, Pb(h_sweep_states20[1]), ylabel="Pb", xlabel="Hr", title="20 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
# plot(p3,p4,p5)
p3 = plot(Hr, Pb(h_sweep_state0), ylabel="Pb", xlabel="Hr", title="0 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
p4 = plot(Hr, Pb(h_sweep_state11), ylabel="Pb", xlabel="Hr", title="11 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
p5 = plot(Hr, Pb(h_sweep_state20), ylabel="Pb", xlabel="Hr", title="20 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr))
plot(p3,p4,p5)
#%%
plot(L_angles(h_sweep_state11, 1),label ="L1")
# plot!(L_angles(h_sweep_state11, 2), label="L2")
# plot!(mod.(L_angles(h_sweep_state11, 1) .- L_angles(h_sweep_state11, 2), pi), label ="L1-L2")
plot!(Gd_angles(h_sweep_state11, 2) .- L_angles(h_sweep_state11, 1), label ="Gd2-L1")
# plot!(L_angles(h_sweep_state11, 2), label="L2")
# plot!(Gd_angles(h_sweep_state11, 1), label ="Gd1")
plot!(Gd_angles(h_sweep_state11, 2), label= "Gd2")
plot!(x->deg2rad(11), label="H", legend=:bottomleft)
plot!(x->pi, label="pi", legend=:bottomleft)
#%%
h_sweep_states16 = map(state_ledgers) do sl
    generate_Hsweep_states(sl, 16, Hr)
end

#%%
l = Ledger(best_model())

Hr = repeat([range(0.0,10.0, length=100);range(10.0, 0.0, length=100)], 2)
base_sweep = generate_Hsweep_states(l, 11, Hr)

plot(Hr, [Etot(base_sweep, :E_dipole) Etot(base_sweep, :E_Gd_L) Etot(base_sweep, :E_Gd) Etot(base_sweep, :E_L) Etot(base_sweep, :E_all)], label=["dipole" "Gd L interaction" "Gd: H" "L: LL, H, easy" "total"], title="Energies", xaxis="H", yaxis="E", legend=:bottomleft, dpi=300)

savefig(papersdir("fig/energy_contributions.png"))
#%%
#%%
Hr = repeat([range(5.5,7.6, length=10);range(7.6, 5.5, length=10)], 2)
base_sweep = generate_Hsweep_states(l, 10, Hr)
Pbs = map(Pb, base_sweep)
p4 = plot(Hr, Pbs, color=Gd.get_colors(ColorSchemes.redblue, Hr), legend=false, xlabel = "H", ylabel="Pb")
ids = [5, 7, 8, 10]
scatter!(p4, [Hr[i] for i in ids], [Pbs[i] for i in ids], color=:red)
# ids = [1, 7, 10, 13, 15, 16, 20, 30, 35, 40]
# ids = [1, 10, 20, 30, 40]

base_sweep= [base_sweep[i] for i in ids]

base_sweep = map(x->set_all_angles!(x, optimize(x, Spin, L).minimizer), Gd.smart_interpolate_states(base_sweep, 10))

plot(p4, plot_fields(base_sweep, 3, 4, 1)..., plot_title="10 degrees", legend=false)

# savefig(plotsdir("fields_0degrees.png"))
#%%
      map(x->x.L[3][1], fields_on_Gds)
      map(x->x.dipole[3][1], fields_on_Gds).-map(x->x.L[3][1], fields_on_Gds)],
scatter([
Hr = repeat([range(5.5,7.6, length=10);range(7.6, 5.5, length=10)], 2)
base_sweep = generate_Hsweep_states(l, 10, Hr)
Pbs = map(Pb, base_sweep)
p4 = plot(Hr, Pbs, color=Gd.get_colors(ColorSchemes.redblue, Hr), legend=false, xlabel = "H", ylabel="Pb")
ids = [5, 7, 8, 10]
scatter!(p4, [Hr[i] for i in ids], [Pbs[i] for i in ids], color=:red)
# ids = [1, 7, 10, 13, 15, 16, 20, 30, 35, 40]
# ids = [1, 10, 20, 30, 40]
base_sweep= [base_sweep[i] for i in ids]
base_sweep = map(x->set_all_angles!(x, optimize(x, Spin, L).minimizer), Gd.smart_interpolate_states(base_sweep, 10))
plot(p4, plot_fields(base_sweep, 3, 4, 1)..., plot_title="10 degrees", legend=false)
# savefig(plotsdir("fields_0degrees.png"))
     [map(x->x.dipole[3][2], fields_on_Gds)
      map(x->x.L[3][2], fields_on_Gds)
      map(x->x.dipole[3][2], fields_on_Gds).-map(x->x.L[3][2], fields_on_Gds)],
     title="Fields",
     xaxis="x",
     yaxis="y",
     color=cs,
     legend=false)
annotate!([(0.5, 0.75, text("Field from L", 10)), (1.2, 0.75, text("field from dipole", 10)), (-1.0, -0.2, text("dipole field - L field",10))])
#%%

#%%
#L1L2 field plots
plot(plot(plot(Hr, [[atan(v[1][2],v[1][1]) for v in L_fields0] [atan(v[2][2], v[2][1]) for v in L_fields0] L_angle0[1] L_angle0[2]],
               label=["l1_field" "l2_field" "l1" "l2"],
               ylabel="angle",
               xlabel="H"),
          plot(Hr, [[norm(v[1]) for v in L_fields0] [norm(v[2]) for v in L_fields0] L_norm0[1] L_norm0[2]],
               ylabel="norm",
               xlabel="H",
               label=["l1_field" "l2_field" "l1" "l2"]),
          title="0 degrees"),
     plot(plot(Hr, [[atan(v[1][2],v[1][1]) for v in L_fields11] [atan(v[2][2], v[2][1]) for v in L_fields11] L_angle11[1] L_angle11[2]],
               label=["l1_field" "l2_field" "l1" "l2"],
               ylabel="angle",
               xlabel="H"),
          plot(Hr, [[norm(v[1]) for v in L_fields11] [norm(v[2]) for v in L_fields11] L_norm11[1] L_norm11[2]],
               ylabel="norm",
               xlabel="H",
               label=["l1_field" "l2_field" "l1" "l2"]),
          title="11 degrees"),
     plot(plot(Hr, [[atan(v[1][2],v[1][1]) for v in L_fields16] [atan(v[2][2], v[2][1]) for v in L_fields16] L_angle16[1] L_angle16[2]],
               label=["l1_field" "l2_field" "l1" "l2"],
               ylabel="angle",
               xlabel="H"),
          plot(Hr, [[norm(v[1]) for v in L_fields16] [norm(v[2]) for v in L_fields16] L_norm16[1] L_norm16[2]],
               ylabel="norm",
               xlabel="H",
               label=["l1_field" "l2_field" "l1" "l2"]),
          title="16 degrees"),
     layout=(3,1),
     size=(1920, 1080))

# savefig(plotsdir("l1_l2_field_Gd.png"))

#%%
# Pb M vs H
plot(plot(plot(Hr, map(Pb, h_sweep_state20), ylabel="Pb", title="20 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr, map(Pb, h_sweep_state11),  title="11 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr, map(Pb, h_sweep_state0),   title="0 degrees", color=Gd.get_colors(ColorSchemes.redblue, Hr)), ylims=[-3,1], layout=(1,3)),
     plot(plot(Hr, map(M, h_sweep_state20), ylabel="M", color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr, map(M, h_sweep_state11), color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr, map(M, h_sweep_state0), color=Gd.get_colors(ColorSchemes.redblue, Hr)), layout=(1,3)),
      plot(plot(Hr, map(x->calc_etot(x,:E_all).e, h_sweep_state20), xlabel="H", ylabel="Energy", color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr,  map(x->calc_etot(x,:E_all).e, h_sweep_state11), xlabel="H", color=Gd.get_colors(ColorSchemes.redblue, Hr)),
          plot(Hr,  map(x->calc_etot(x,:E_all).e, h_sweep_state0), xlabel="H", color=Gd.get_colors(ColorSchemes.redblue, Hr)), layout=(1,3)),
       layout=(3,1),legend=false)
savefig(plotsdir("Pb_M_E.png"))
#%%

plot(plot(Hr, [Etot(h_sweep_state20, :E_dipole) Etot(h_sweep_state20, :E_Gd_L) Etot(h_sweep_state20, :E_Gd) Etot(h_sweep_state20, :E_L)], label=["dipole" "Gd L interaction" "Gd: H" "L: LL, H, easy"], title="20 degrees", yaxis="E", legend=:bottomleft),
           plot(Hr, [Etot(h_sweep_state11, :E_dipole) Etot(h_sweep_state11, :E_Gd_L) Etot(h_sweep_state11, :E_Gd) Etot(h_sweep_state11, :E_L)], label=["dipole" "Gd L interaction" "Gd: H" "L: LL, H, easy"], title="11 degrees",  yaxis="E", legend=false),
           plot(Hr, [Etot(h_sweep_state0, :E_dipole) Etot(h_sweep_state0, :E_Gd_L) Etot(h_sweep_state0, :E_Gd) Etot(h_sweep_state0, :E_L)], label=["dipole" "Gd L interaction" "Gd: H" "L: LL, H, easy"], title="0 degrees", xaxis="H", yaxis="E", legend=false),
      layout=(3,1), legendfontsize=10, dpi=300)

savefig(plotsdir("energy_contributions.png"))
#%%
l = Ledger(best_model())
l[H].data[1] = H(magnitude=0.0, ϕ=deg2rad(10))

start = 2
stop = 10.5
len = 20
Hr = H_sweep_range(start, stop, len)
hsweep = generate_Hsweep_states(l, 10, Hr)
plot(Hr, Pb.(hsweep))
#%%
d = Diorama(l)
expose(d)
#%%
### MINIMAL STUFF
# Model

singleton(l, Model).J1 = 5.0
singleton(l, Model).χ_Gd = 0.4
singleton(l, Model).l2h2_strength = 0.0
singleton(l, Model).dipole_strength = 10.5
singleton(l, Model).K_Gd = 0.0
singleton(l, Model).K_L /= 2
#%%
# New energies
function GdMn2O5.Overseer.update(::GdMn2O5.E_Gd_L, m::GdMn2O5.Overseer.AbstractLedger)
    etot     = singleton(m, Etot)
    p        = singleton(m, Model)
    Gd_spins = m[Spin]
    Lc       = m[L]
    L1 = Lc[1].v
    J1 = p.J1
    etot.e += J1 * Gd_spins[1].v ⋅ L1 -
              J1 * Gd_spins[2].v ⋅ L1
end

function GdMn2O5.Overseer.update(::GdMn2O5.E_Dipole, m::GdMn2O5.Overseer.AbstractLedger)
    etot = singleton(m, Etot)
    p = singleton(m, Model)
    spin = m[Spin]
    bonds = m[Bond]
    dip_tot = 0.0
    @inbounds for b in bonds
        pr = b.r
        s1 = spin[b.i1].v
        s2 = spin[b.i2].v
        dip_tot -= 0.5 * p.dipole_strength * (3 * (s1 ⋅ pr) * (s2 ⋅ pr))*b.inv_nr5  #0.5 for double counting
    end
    etot.e += dip_tot
end

function GdMn2O5.Overseer.update(::GdMn2O5.E_Pb, m::GdMn2O5.Overseer.AbstractLedger)
    etot = singleton(m, Etot)
    p = singleton(m, Model)
    Ls = m[L]
    L1 = Ls[1].v
    Gd_spins = m[Spin]
    pb = singleton(m, Pb).pb
    etot.e -= (Gd_spins[1].v - Gd_spins[2].v) ⋅ L1 * pb * p.α
    etot.e -= 8p.g * L1 ⋅ GdMn2O5.polar2vec(deg2rad(-23.4), 0.0) * pb
    etot.e += pb^2/2 - pb * singleton(m, E_field).e
end

function GdMn2O5.set_all_angles!(l::AbstractLedger, angles)
    l[Spin].data[1] = Spin(ϕ=angles[1])
    l[Spin].data[2] = Spin(ϕ=angles[2])
    l[L].data[1] = L(ϕ=angles[3])
    return l
end

function minimal_ledger()
    model = minimal_model(dipole_cutoff=3.6, K_L = 3.95)
    l = Ledger(Stage(:E_all, [GdMn2O5.E_Gd_H(), GdMn2O5.E_L_H(), GdMn2O5.E_L_easy(), GdMn2O5.E_Gd_L(), GdMn2O5.E_Dipole(), GdMn2O5.E_Pb()]),
               Stage(:E_Gd_h, [GdMn2O5.E_Gd_H()]),
               Stage(:E_L_h, [GdMn2O5.E_L_H()]),
               Stage(:E_Gd_L, [GdMn2O5.E_Gd_L()]),
               Stage(:E_L_easy, [GdMn2O5.E_L_easy()]),
               Stage(:E_Dipole, [GdMn2O5.E_Dipole()]),
               Stage(:Visualization, [GdMn2O5.SpinUpdater(), GdMn2O5.Visualizer()]))

    #original Gd2
    Gd1_position = Position(Vec3(1.0257574253082276, 1.4649663247357028, 0.0))
    #original Gd7, shifted to -1, -1
    Gd2_position = Position(Vec3(-1.0257574253082276, -1.4649663247357028, 0.0))

    Entity(l, Gd1_position, Spin(ϕ=deg2rad(12)))
    Entity(l, Gd2_position, Spin(ϕ=deg2rad(-12)))

    Entity(l, Bond(1, 2, Gd1_position.p - Gd2_position.p))
    Entity(l, Bond(2, 1, Gd2_position.p - Gd1_position.p))

    Entity(l, L(ϕ=deg2rad(23.4)), EasyAxis(ϕ=deg2rad(23.4)))

    Entity(l, model, Etot(0.0), H(ϕ=0.0, magnitude=0.0), Pb(), E_field())
    res = optimize(l)
    set_optim_values!(l, res.minimizer)
    return l
end
#%%
using ProgressMeter
function calculate_Etots(hsweep...; N=150, sym=:E_all)
    p = Progress(length(hsweep))
    S1_ϕs = range(0,2π, length=N)
    S2_ϕs = range(0,2π, length=N)
    E_tots = [zeros(N,N) for i=1:length(hsweep)]
    for ih = 1:length(hsweep)
        thread_ls = [deepcopy(hsweep[ih]) for t=1:Threads.nthreads()]
        Threads.@threads for i1 = 1:N
            for (i2, ϕ2) in enumerate(S2_ϕs)
                ϕ1 = S1_ϕs[i1]
                tl = thread_ls[Threads.threadid()]
                tl[Spin].data[1] = GdMn2O5.Spin(ϕ = ϕ1)
                tl[Spin].data[2] = GdMn2O5.Spin(ϕ = ϕ2)
                set_optim_values!(tl,  optimize(tl, L, Pb).minimizer, (L, Pb))
                # set_angles!(tl, optimize(tl, Spin).minimizer)
                E_tots[ih][i1,i2] = calc_etot(tl, sym).e
            end
        end
        next!(p)
    end
    return E_tots
end

function plot_heatmap(h_ledger, hsweep, L1_traj_shift, L2_traj_shift, offset=1.0)
    ϕ1_r = range(0,2π, length=N)
    ϕ2_r = range(0,2π, length=N)
    E_tots = calculate_Etots(h_ledger)[1]
    cs = :rainbow
    p = heatmap(ϕ1_r, ϕ2_r, E_tots, clims=(minimum(E_tots), minimum(E_tots)+offset), colorbar_title = "E",color=cs, aspect_ratio=1,legend=false)
    plot!(p, map(x->x[L][2].ϕ+L2_traj_shift, hsweep), map(x->x[L][1].ϕ+L1_traj_shift, hsweep),color=RGB(242/255, 239/255, 26/255))
    scatter!(p, [h_ledger[L][2].ϕ+L2_traj_shift], [h_ledger[L][1].ϕ+L1_traj_shift], color=RGB(242/255, 239/255, 26/255),markersize=4)
    return p
end
#%%
l = minimal_ledger()

start = 2
stop = 10.5
len = 20
Hr = H_sweep_range(start, stop, len)
hsweep = generate_Hsweep_states(l, 12, Hr)

plot(Hr, map(x->calc_etot(x,:E_Dipole).e, hsweep))
plot(Hr, Gd_angles(hsweep,1))
GdMn2O5.Overseer.ensure_component!(l, GdMn2O5.Optimize)
Etots = calculate_Etots(hsweep...; sym=:E_Dipole)
Etots2 = calculate_Etots(hsweep...; sym=:E_Gd_L)

Etots2 = calculate_Etots(hsweep...)
#%%
t = 20
heatmap(range(0,2π, length=150), range(0,2π, length=150), Etots2[t].+Etots[t], clims=(-10,-9))
# plot(, heatmap(range(0,2π, length=150), range(0,2π, length=150), Etots2[t],clims=(-12,-11)))
#%%
plot(Hr, Pb.(hsweep))

GdMn2O5.get_colors(ColorSchemes.rainbow, Hr)
plot(H_magnitudes(hsweep), Gd_angles(hsweep,2).-Gd_angles(hsweep,1))
plot(H_magnitudes(hsweep), map(x->calc_etot(x, :E_Gd_h).e + calc_etot(x, :E_Dipole).e, hsweep),color=GdMn2O5.get_colors(ColorSchemes.rainbow, Hr))
plot(H_magnitudes(hsweep), L_angles(hsweep,1))
plot!(Hr, Gd_angles(hsweep, 2))

Entity(l, GdMn2O5.VisualizationSettings())
d = Diorama(l);
expose(d)
d.loop
#%%
l[H].data[1] = H(magnitude=3.0, ϕ=deg2rad(6))

deleteat!(l[Bond].data, findall(x->x.i1 ∈ (1,2,7, 4,5,8) || x.i2 ∈ (1,2,7, 4,5,8), l[Bond].data))
findall(x->!(x.i1 ∈ (1,2,7, 4,5,8) || x.i2 ∈ (1,2,7, 4,5,8)), l[Bond].data)
set_angles!(l, optimize(l, Spin, L).minimizer)
#%%
l = Ledger(zero_dipole_model())

d = Diorama(l)

filter(x->x.i1 == 6 && x.i2 == 3, l[Bond].data)
filter(x->x.i1 == 2, l[Bond].data)
filter(x->x.i1 == 3, l[Bond].data)
filter(x->x.i1 == 7 && x.i2 == 2, l[Bond].data)
calc_etot(hsweep[80], :E_dipole)

l[Bond][1].r
expose(d)
l[Bond][2].r
l[Spin][1].v
l[Spin][2].v
#%%
### Energy depending on L1 L2 orientation in fixed field
# l = Ledger(zero_dipole_model(γ=0.00, χ_L=0.1))
# Hr = H_sweep_range(0,10.3,20)
# hsweep = generate_Hsweep_states(l, 30, Hr)
# plot(Hr, Pb.(hsweep))
# l = Ledger(zero_dipole_model(γ=0.10, χ_L=0.3, K_L=-2.0, K_Gd=0.0, J2=05.1))
l = Ledger(zero_dipole_model(γ=0.10, χ_L=0.3, K_L=-2.0, K_Gd=0.6))
Hr = H_sweep_range(0, 5.0,20)
hsweep = generate_Hsweep_states(l, 10, Hr)
#%%
E_tots = calculate_Etots(hsweep[10:2:16]...,hsweep[26:2:30]...,hsweep[36])
#%%
Hr2 = H_sweep_range(0.0,5.0,100)
hsweep2 = generate_Hsweep_states(l, 10, Hr)
plts = []
ϕ1_r = range(0,2π, length=N)
ϕ2_r = range(0,2π, length=N)
ylab = "L1 phi(rad)"
xlab = "L2 phi(rad)"
cs = :rainbow
offset = 0.05
for (h, e) in zip((hsweep[10:2:16]...,hsweep[26:2:30]..., hsweep[36]), E_tots)
    tit = "H=$(norm(h[H][1].v))"
    p = heatmap(ϕ1_r, ϕ2_r, e, ylabel = ylab, xlabel=xlab, clims=(minimum(e), minimum(e)+offset), title=tit, colorbar_title = "E",color=cs, aspect_ratio=1,legend=false)
    plot!(p, map(x->x[L][2].ϕ+2pi, hsweep2), map(x->x[L][1].ϕ+2pi, hsweep2),color=:white)
    scatter!(p, [h[L][2].ϕ+2pi], [h[L][1].ϕ+2pi], color=:white,markersize=10)
    push!(plts, p)
end

plot(plts..., margin=0.1Plots.Measures.mm,layout=(2,4))
#%
contour(range(-π,π, length=N), range(-π,π, length=N), E_tots[28], ylabel = "L1 phi(rad)",xlabel = "L2 phi(rad)",colorbar_title = "E",color=:rainbow, fill=true)
Hr[28]
savefig(plotsdir("L1L2_energy_H_zero_dipole.png"))
#%%
l = Ledger(physical_model())

#%%


text:
@params mutable struct PCGDisplacementSolver{T, dim, TP<:AbstractPenalty{T}} <: AbstractDisplacementSolver
    problem::StiffnessTopOptProblem{dim, T}
    globalinfo::GlobalFEAInfo{T}
    elementinfo::ElementFEAInfo{dim, T}
    u::AbstractVector{T}
    lhs::AbstractVector{T}
    rhs::AbstractVector{T}
    vars::AbstractVector{T}
    penalty::TP
    prev_penalty::TP
    xmin::T
    cg_max_iter::Integer
    tol::T
    cg_statevars::CGStateVariables{T, <:AbstractVector{T}}
    preconditioner
    preconditioner_initialized::Ref{Bool}
    conv
end
function PCGDisplacementSolver(sp::StiffnessTopOptProblem{dim, T};
    conv = DefaultCriteria(),
    xmin=T(1)/1000,
    cg_max_iter=700,
    tol=xmin,
    penalty=PowerPenalty{T}(1),
    prev_penalty=copy(penalty),
    preconditioner=identity,
    quad_order=default_quad_order(sp)) where {dim, T, Tconv}

    prev_penalty = @set prev_penalty.p = T(NaN)
    elementinfo = ElementFEAInfo(sp, quad_order, Val{:Static})
    globalinfo = GlobalFEAInfo(sp)
    u = zeros(T, ndofs(sp.ch.dh))
    lhs = similar(u)
    rhs = similar(u)
    vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white))
    varind = sp.varind
    us = similar(u) .= 0
    cg_statevars = CGStateVariables{eltype(u), typeof(u)}(us, similar(u), similar(u))

    return PCGDisplacementSolver(sp, globalinfo, elementinfo, u, lhs, rhs, vars, penalty, prev_penalty, xmin, cg_max_iter, tol, cg_statevars, preconditioner, Ref(false), conv)
end

function (s::PCGDisplacementSolver{T})(::Type{Val{safe}} = Val{false}; assemble_f = true) where {T, safe}
    rhs = assemble_f ? s.globalinfo.f : s.rhs
    lhs = assemble_f ? s.u : s.lhs
    globalinfo = GlobalFEAInfo(s.globalinfo.K, rhs)
    assemble!(globalinfo, s.problem, s.elementinfo, s.vars, s.penalty, s.xmin, assemble_f = assemble_f)
    Tconv = typeof(s.conv)
    K, f = globalinfo.K, globalinfo.f
    if safe
        m = meandiag(K)
        for i in 1:size(K,1)
            if K[i,i] ≈ zero(T)
                K[i,i] = m
            end
        end
    end

    @unpack cg_max_iter, tol, cg_statevars = s
    @unpack preconditioner, preconditioner_initialized = s

    _K = K isa Symmetric ? K.data : K
    if !(preconditioner === identity)
        if !preconditioner_initialized[]
            UpdatePreconditioner!(preconditioner, _K)
            preconditioner_initialized[] = true
        end
    end
    op = MatrixOperator(_K, f, s.conv)
    if preconditioner === identity
        return cg!(lhs, op, f, tol=tol, maxiter=cg_max_iter, log=false, statevars=cg_statevars, initially_zero=false)
    else
        return cg!(lhs, op, f, tol=tol, maxiter=cg_max_iter, log=false, statevars=cg_statevars, initially_zero=false, Pl = preconditioner)
    end
end


text:
# Implementations of CGD for linear games
using LinearAlgebra
using StaticArrays
using LinearMaps
using IterativeSolvers
using BlockArrays

# The dual of the quadratic distance
# ̄ψ = sum( α .* x.^2 ) / 2
function squared_distance(α)
    # Computing the primal point associated to the dual point x̄
    function Dxψ(x̄)
      return α .\ x̄
    end

    function Dxψ_inv(x̄)
      return α .* x̄
    end

    # Hessian of the dual potential, evaluated in the primal point x
    function Dxxψ(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v ./ α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end

    # Hessian of the primal potential, evaluated in the primal point x
    function Dxxψ_inv(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v .* α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end
    return Dxψ, Dxψ_inv, Dxxψ, Dxxψ_inv
end

# The dual of the quadratic distance
# ̄ψ = sum( α .* x.^2 ) / 2
function positive_projection(α)
    # Computing the primal point associated to the dual point x̄
    function Dxψ(x̄)
      x̄ .= max.(x̄, zero(eltype(x̄)))
      return α .\ x̄
    end

    function Dxψ_inv(x̄)
      x̄ .= max.(x̄, zero(eltype(x̄)))
      return α .* x̄
    end

    # Hessian of the dual potential, evaluated in the primal point x
    function Dxxψ(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v ./ α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end

    # Hessian of the primal potential, evaluated in the primal point x
    function Dxxψ_inv(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v .* α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end
    return Dxψ, Dxψ_inv, Dxxψ, Dxxψ_inv
end

# The dual of the (negative) shannon entropy
# ̄ψ = sum( α .* (x .* log.(x) .- x) )
function shannon_entropy(α)
    # Computing the primal point associated to the dual point x̄
    function Dxψ(x̄)
      return exp.(α .\ x̄)
    end

    function Dxψ_inv(x)
      return α .* log.(x)
    end

    # Hessian of the dual potential, evaluated in the primal point x
    function Dxxψ(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v .* x ./ α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end

    # Hessian of the primal potential, evaluated in the primal point x
    function Dxxψ_inv(x)
      function Dxxψ_v!(hvp, v)
        hvp .= v ./ x .* α
      end
      return LinearMap{eltype(x)}(Dxxψ_v!, length(x); ismutating=true, issymmetric=true, ishermitian=true)
    end
    return Dxψ, Dxψ_inv, Dxxψ, Dxxψ_inv
end

# Combines a list of potentials to act on different parts of a block-vector
function combine_potentials(potentials...)
  function Dxψ(x̄::AbstractBlockVector)
    out = similar(x̄)
    @assert length(potentials) == blocklength(x̄)
    for k = 1 : length(potentials)
      out[Block(k)] .= potentials[k][1](x̄[Block(k)])
    end
    return out
  end

  function Dxψ_inv(x̄::AbstractBlockVector)
    out = similar(x̄)
    @assert length(potentials) == blocklength(x̄)
    for k = 1 : length(potentials)
      out[Block(k)] .= potentials[k][2](x̄[Block(k)])
    end
    return out
  end

  # Hessian of the dual potential, evaluated in the primal point x
  function Dxxψ(x)
    @assert length(potentials) == blocklength(x)
    block_diagonal = [potentials[k][3](x[Block(k)]) for k = 1 : length(potentials)]

    return LinearMaps.BlockDiagonalMap(block_diagonal...)
  end

  # Hessian of the primal potential, evaluated in the primal point x
  function Dxxψ_inv(x)
    @assert length(potentials) == blocklength(x)
    block_diagonal = [potentials[k][4](x[Block(k)]) for k = 1 : length(potentials)]

    return LinearMaps.BlockDiagonalMap(block_diagonal...)
  end

  return Dxψ, Dxψ_inv, Dxxψ, Dxxψ_inv
end

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
module AxisAlgorithms

using WoodburyMatrices
using LinearAlgebra, SparseArrays

export A_ldiv_B_md!,
    A_ldiv_B_md,
    A_mul_B_md!,
    A_mul_B_md,
    A_mul_B_perm!,
    A_mul_B_perm

"""
`A_ldiv_B_md(F, src, dim)` solves `F\b` for slices `b` of `src` along dimension `dim`,
storing the result along the same dimension of the output.
Currently, `F` must be an LU-factorized tridiagonal matrix or a Woodbury matrix.
"""
A_ldiv_B_md(F, src, dim::Integer) = A_ldiv_B_md!(similar(src), F, src, dim)

"""
`A_mul_B_md(M, src, dim)` computes `M*x` for slices `x` of `src` along dimension `dim`,
storing the resulting vector along the same dimension of the output.
`M` must be an `AbstractMatrix`. This uses an in-place naive algorithm.
"""
A_mul_B_md(M::AbstractMatrix, src, dim::Integer) = A_mul_B_md!(alloc_matmul(M,src,dim), M, src, dim)

"""
`A_mul_B_perm(M, src, dim)` computes `M*x` for slices `x` of `src` along dimension `dim`, storing the
resulting vector along the same dimension of the output.
`M` must be an `AbstractMatrix`. This uses `permutedims` to make dimension
`dim` into the first dimension, performs a standard matrix multiplication, and restores the original
dimension ordering. In many cases, this algorithm exhibits the best cache behavior.
"""
A_mul_B_perm(M::AbstractMatrix, src, dim::Integer) = A_mul_B_perm!(alloc_matmul(M,src,dim), M, src, dim)

function alloc_matmul(M,src::AbstractArray{S,N},dim) where {S,N}
    sz = [size(src)...]
    sz[dim] = size(M,1)
    T = Base.promote_op(*, eltype(M), S)
    Array{T,N}(undef, sz...)
end

include("tridiag.jl")
include("matmul.jl")
include("woodbury.jl")

end # module

max_stars_repo_name: timholy/FilterMD

max_issues_repo_name: timholy/FilterMD

max_forks_repo_name: timholy/FilterMD


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

mutable struct IoK8sApiCoreV1ScaleIOVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    gateway::Any # spec type: Union{ Nothing, String } # spec name: gateway
    protectionDomain::Any # spec type: Union{ Nothing, String } # spec name: protectionDomain
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalObjectReference } # spec name: secretRef
    sslEnabled::Any # spec type: Union{ Nothing, Bool } # spec name: sslEnabled
    storageMode::Any # spec type: Union{ Nothing, String } # spec name: storageMode
    storagePool::Any # spec type: Union{ Nothing, String } # spec name: storagePool
    system::Any # spec type: Union{ Nothing, String } # spec name: system
    volumeName::Any # spec type: Union{ Nothing, String } # spec name: volumeName

    function IoK8sApiCoreV1ScaleIOVolumeSource(;fsType=nothing, gateway=nothing, protectionDomain=nothing, readOnly=nothing, secretRef=nothing, sslEnabled=nothing, storageMode=nothing, storagePool=nothing, system=nothing, volumeName=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("gateway"), gateway)
        setfield!(o, Symbol("gateway"), gateway)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("protectionDomain"), protectionDomain)
        setfield!(o, Symbol("protectionDomain"), protectionDomain)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("sslEnabled"), sslEnabled)
        setfield!(o, Symbol("sslEnabled"), sslEnabled)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("storageMode"), storageMode)
        setfield!(o, Symbol("storageMode"), storageMode)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("storagePool"), storagePool)
        setfield!(o, Symbol("storagePool"), storagePool)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("system"), system)
        setfield!(o, Symbol("system"), system)
        validate_property(IoK8sApiCoreV1ScaleIOVolumeSource, Symbol("volumeName"), volumeName)
        setfield!(o, Symbol("volumeName"), volumeName)
        o
    end
end # type IoK8sApiCoreV1ScaleIOVolumeSource

const _property_map_IoK8sApiCoreV1ScaleIOVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("gateway")=>Symbol("gateway"), Symbol("protectionDomain")=>Symbol("protectionDomain"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("sslEnabled")=>Symbol("sslEnabled"), Symbol("storageMode")=>Symbol("storageMode"), Symbol("storagePool")=>Symbol("storagePool"), Symbol("system")=>Symbol("system"), Symbol("volumeName")=>Symbol("volumeName"))
const _property_types_IoK8sApiCoreV1ScaleIOVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("gateway")=>"String", Symbol("protectionDomain")=>"String", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1LocalObjectReference", Symbol("sslEnabled")=>"Bool", Symbol("storageMode")=>"String", Symbol("storagePool")=>"String", Symbol("system")=>"String", Symbol("volumeName")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ScaleIOVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ScaleIOVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ScaleIOVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ScaleIOVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ScaleIOVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ScaleIOVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ScaleIOVolumeSource)
    (getproperty(o, Symbol("gateway")) === nothing) && (return false)
    (getproperty(o, Symbol("secretRef")) === nothing) && (return false)
    (getproperty(o, Symbol("system")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ScaleIOVolumeSource }, name::Symbol, val)
end


text:
module GreekAnalyticalInv

__precompile__(true)

using Distributions

export blackScholesInvValue, blackScholesInvDelta, blackScholesInvGamma, blackScholesInvVega, blackScholesInvTheta
export blackScholesInvVanna, blackScholesInvVolga, blackScholesInvCharm
export blackScholesInvRho, blackScholesInvRho2

std_normal = Normal(0, 1)

## y = funding rate, corresponds to Espen's book, is q.
## y = r - b, b = r - y

function blackScholesInvValue(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1

    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    sign * (exp(-r * τ) * cdf(std_normal, sign * d2) - exp((y - r + σ^2) * τ) * strike / spot * cdf(std_normal, sign * d3))
end

function blackScholesInvDelta(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    (exp(-r * τ) - exp((y - r) * (1 - σ) * τ)) / spot / σ / sqrt(τ) * pdf(std_normal, d2) + sign * exp((y - r + σ*σ) * τ) * strike / spot / spot * cdf(std_normal, sign * d3)
end

function blackScholesInvGamma(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = (1.5 + (r - y) / σ + log(spot / strike) / σ / σ / τ) * exp((y - r) * (1 - σ) * τ)
    pdf_coeff2 = (0.5 + (r - y) / σ + log(spot / strike) / σ / σ / τ) * exp(-r * τ)

    (pdf_coeff1 - pdf_coeff2) / spot / spot / σ / sqrt(τ) * pdf(std_normal, d2) - sign * exp((y - r + σ * σ) * τ) / (spot^3) * strike * 2 * cdf(std_normal, sign * d3)
end

function blackScholesInvVega(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = exp((y - r) * (1 - σ) * τ) * (1.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ))
    pdf_coeff2 = exp(-r * τ) * (0.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ))

    (pdf_coeff1 - pdf_coeff2) * pdf(std_normal, d2) - sign * 2 * exp((y - r + σ * σ) * τ) * strike / spot * σ * τ * cdf(std_normal, sign * d3)
end

function blackScholesInvTheta(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = exp((y - r) * (1 - σ) * τ) * ((y - r + 1.5 * σ) / 2 / sqrt(τ) + log(spot / strike) / 2 / σ / τ / sqrt(τ))
    pdf_coeff2 = exp(-r * τ) * ((y - r + 0.5 * σ) / 2 / sqrt(τ) + log(spot / strike) / 2 / σ / τ / sqrt(τ))

    sign * exp(-r * τ) * r * cdf(std_normal, sign * d2) + sign * exp((y - r + σ * σ) * τ) * strike / spot * (y - r + σ * σ) * cdf(std_normal, sign * d3) - (pdf_coeff1 - pdf_coeff2) * pdf(std_normal, d2)
end

function blackScholesInvVanna(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = exp((y - r) * (1 - σ) * τ) * (y - r) * sqrt(τ) / spot / σ
    pdf_coeff2 = (exp(-r * τ) -  exp((y - r) * (1 - σ) * τ)) / spot / σ / σ / sqrt(τ)
    pdf_coeff3 = exp((y - r) * (1 - σ) * τ) / spot * (1.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ))
    pdf_coeff4 = σ * pdf_coeff2 * (0.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ)) * ((r - y - 0.5 * σ) * sqrt(τ) + log(spot / strike) / σ / sqrt(τ))

    sign * 2 * strike * σ * τ / spot / spot * exp((y - r + σ * σ) * τ) * cdf(std_normal, sign * d3) + (pdf_coeff1 - pdf_coeff2 - pdf_coeff3 + pdf_coeff4) * pdf(std_normal, d2)
end

function blackScholesInvVolga(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = 2 * exp((y - r) * (1 - σ) * τ) * σ * τ * (1.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ))
    pdf_coeff2 = (0.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ)) * ((r - y - 0.5 * σ) * sqrt(τ) + log(spot / strike) / σ / sqrt(τ)) * (exp((y - r) * (1 - σ) * τ) * (1.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ)) - exp(-r * τ) * (0.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ)))
    pdf_coeff3 = 2 * log(spot / strike) / σ / σ / σ / sqrt(τ) * (exp(-r * τ) -  exp((y - r) * (1 - σ) * τ)) - exp((y - r) * (1 - σ) * τ) * (y - r) * τ * (1.5 * sqrt(τ) + log(spot / strike) / σ / σ / sqrt(τ))

    -2 * sign * exp((y - r + σ * σ) * τ) * strike * τ / spot * (1 + 2 * σ * σ * τ) * cdf(std_normal, sign * d3) + (pdf_coeff1 + pdf_coeff2 + pdf_coeff3) * pdf(std_normal, d2)
end

function blackScholesInvCharm(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    pdf_coeff1 = (exp(-r * τ) -  exp((y - r) * (1 - σ) * τ)) / 2 / spot / σ / τ / sqrt(τ)
    pdf_coeff2 = (r * exp(-r * τ) + (y - r) * (1 - σ) * exp((y - r) * (1 - σ) * τ)) / spot / σ / sqrt(τ)
    pdf_coeff3 = ((y - r + 1.5 * σ) / 2 / sqrt(τ) + log(spot / strike) / 2 / σ / τ / sqrt(τ)) * exp((y - r) * (1 - σ) * τ) / spot
    pdf_coeff4 = (exp(-r * τ) -  exp((y - r) * (1 - σ) * τ)) * ((y - r + 0.5 * σ) / 2 / sqrt(τ) + log(spot / strike) / 2 / σ / τ / sqrt(τ)) * ((y - r + 0.5 * σ) * sqrt(τ) - log(spot / strike) / σ / sqrt(τ)) / spot / σ / sqrt(τ)

    (pdf_coeff1 + pdf_coeff2 + pdf_coeff3 + pdf_coeff4) * pdf(std_normal, d2) - sign * exp((y - r + σ * σ) * τ) / spot / spot * strike * (y - r + σ * σ) * cdf(std_normal, sign * d3)
end

function blackScholesInvRho2(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    (exp((y - r) * (1 - σ) * τ) - exp(-r * τ)) * sqrt(τ) * pdf(std_normal, d2) - sign * exp((y - r + σ * σ) * τ) * strike * τ / spot * cdf(std_normal, sign * d3)
end

function blackScholesInvRho(isCall::Bool, spot::T, strike::T, τ::T, r::T, y::T, σ::T) where {T}
    sign = isCall ? 1 : -1
    d1 = log(spot / strike) / σ / sqrt(τ) + (r - y + 0.5 * σ) * sqrt(τ)
    d2 = d1 - σ * sqrt(τ)
    d3 = d2 - σ * sqrt(τ)

    sign * strike / spot * τ * exp((y - r + σ * σ) * τ) * cdf(std_normal, sign * d3) - sign * exp(-r * τ) * τ * cdf(std_normal, sign * d2) + (exp(-r * τ) - exp((y - r) * (1 - σ) * τ)) * sqrt(τ) * pdf(std_normal, d2)
end

end

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
export Transition

"""
Description:

    The transition node models a transition between discrete
    random variables, with node function

    f(out, in1, a) = Cat(out | a*in1)

    Where a is a left-stochastic matrix (columns sum to one)

Interfaces:

    1. out
    2. in1
    3. a

Construction:

    Transition(out, in1, a, id=:some_id)
"""
mutable struct Transition <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Transition(out, in1, a; id=generateId(Transition))
        @ensureVariables(out, in1, a)
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:a] = self.interfaces[3] = associate!(Interface(self), a)

        return self
    end
end

slug(::Type{Transition}) = "T"

# Average energy functional
function averageEnergy(::Type{Transition}, marg_out::ProbabilityDistribution, marg_in1::ProbabilityDistribution, marg_a::ProbabilityDistribution{MatrixVariate})
    -unsafeMeanVector(marg_out)'*unsafeLogMean(marg_a)*unsafeMeanVector(marg_in1)
end

function averageEnergy(::Type{Transition}, marg_out_in1::ProbabilityDistribution{Multivariate, Contingency}, marg_a::ProbabilityDistribution{MatrixVariate})
    -tr(marg_out_in1.params[:p]'*unsafeLogMean(marg_a))
end


text:
module TestBroadcast
    using NullableArrays
    using Base.Test

    A1 = rand(10)
    M1 = rand(Bool, 10)
    n = rand(2:5)
    dims = [ rand(2:5) for i in 1:n]
    A2 = rand(10, dims...)
    M2 = rand(Bool, 10, dims...)
    C2 = Array(Float64, 10, dims...)
    i = rand(2:5)
    A3 = rand(10, [dims; i]...)
    M3 = rand(Bool, 10, [dims; i]...)
    C3 = Array(Float64, 10, [dims; i]...)

    m1 = broadcast((x,y)->x, M1, M2)
    m2 = broadcast((x,y)->x, M2, M3)
    R2 = reshape(Bool[ (x->(m1[x] | M2[x]))(i) for i in 1:length(M2) ], size(M2))
    R3 = reshape(Bool[ (x->(m2[x] | M3[x]))(i) for i in 1:length(M3) ], size(M3))
    r2 = broadcast((x,y)->x, R2, R3)
    Q3 = reshape(Bool[ (x->(r2[x] | R3[x]))(i) for i in 1:length(R3) ], size(M3))

    U1 = NullableArray(A1)
    U2 = NullableArray(A2)
    U3 = NullableArray(A3)
    V1 = NullableArray(A1, M1)
    V2 = NullableArray(A2, M2)
    V3 = NullableArray(A3, M3)
    Z2 = NullableArray(Float64, 10, dims...)
    Z3 = NullableArray(Float64, 10, [dims; i]...)

    f() = 5
    f(x) = Nullable(5) * x
    f(x, y) = x + y
    f(x, y, z) = x + y + z
    g() = 5
    g(x::Float64) = 5 * x
    g(x::Float64, y::Float64) = x * y
    g(x::Float64, y::Float64, z::Float64) = x * y * z

    i = 1
    for (dests, arrays, nullablearrays, mask) in
        ( ((C2, Z2), (A1, A2), (U1, U2), ()),
          ((C3, Z3), (A2, A3), (U2, U3), ()),
          ((C3, Z3), (A1, A2, A3), (U1, U2, U3), ()),

          ((C2, Z2), (A1, A2), (V1, V2), (R2,)),
          ((C3, Z3), (A2, A3), (V2, V3), (R3,)),
          ((C3, Z3), (A1, A2, A3), (V1, V2, V3), (Q3,)),
    )

    # Base.broadcast!(f, B::NullableArray, As::NullableArray...; lift::Bool=false)
        broadcast!(f, dests[1], arrays...)
        broadcast!(f, dests[2], nullablearrays...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))

        broadcast!(g, dests[1], arrays...)
        broadcast!(g, dests[2], nullablearrays...; lift=true)
        @test isequal(dests[2], NullableArray(dests[1], mask...))

        # Base.broadcast(f, As::NullableArray...;lift::Bool=false)
        D = broadcast(f, arrays...)
        X = broadcast(f, nullablearrays...)
        @test isequal(X, NullableArray(D, mask...))

        D = broadcast(g, arrays...)
        X = broadcast(g, nullablearrays...; lift=true)
        @test isequal(X, NullableArray(D, mask...))
    end

    # Base.broadcast!(f, X::NullableArray; lift::Bool=false)
    for (array, nullablearray, mask) in
        ( (A1, U1, ()), (A2, U2, ()), (A3, U3, ()),
          (A1, V1, (M1,)), (A2, V2, (M2,)), (A3, V3, (M3,)),
    )
        broadcast!(f, array)
        broadcast!(f, nullablearray)
        @test isequal(nullablearray, NullableArray(array, mask...))

        broadcast!(g, array)
        broadcast!(g, nullablearray; lift=true)
        @test isequal(nullablearray, NullableArray(array, mask...))
    end

    # test broadcasted arithmetic operators
    A = rand(10)
    X1 = NullableArray(A)
    n = rand(2:5)
    dims = rand(2:5, n)
    B = rand(Float64, 10, dims...)
    X2 = NullableArray(B)
    M = rand(Bool, 10, dims...)
    Y = NullableArray(B, M)

    for op in (
        (.+),
        (.-),
        (.*),
        (./),
        (.%),
        (.^),
        (.==),
        (.!=),
        (.<),
        (.>),
        (.<=),
        (.>=),
    )
        @test isequal(op(X1, X2), NullableArray(op(A, B)))
        @test isequal(op(X1, Y), NullableArray(op(A, B), M))
    end

end # module


text:
function assess(M::Model, y, ŷ)
  assess(M.obj, y, ŷ)
end

function assessgrad(M::Model, y, ŷ, x)
  assessgrad(M.obj, y, ŷ, x)
end

"""
L1 Distance Objective
"""
struct L1_Obj <: Objective end

function assess(O::L1_Obj, y, ŷ)
  value(L1DistLoss(), y, ŷ, AvgMode.Sum())
end

function assessgrad(O::L1_Obj, y, ŷ, x)
  deriv(L1DistLoss(), y, ŷ) * x
end

"""
L2 Distance Objective
"""
struct L2_Obj <: Objective end

function assess(O::L2_Obj, y, ŷ)
  value(L2DistLoss(), y, ŷ, AvgMode.Sum())
end

function assessgrad(O::L2_Obj, y, ŷ, x)
  deriv(L2DistLoss(), y, ŷ) * x
end

"""
Logistic Objective
"""
struct Logit_Obj <: Objective end

function assess(O::Logit_Obj, y, ŷ)
  value(LogitDistLoss(), y, ŷ, AvgMode.Sum())
end

function assessgrad(O::Logit_Obj, y, ŷ, x)
  deriv(LogitDistLoss(), y, ŷ) * x
end

#
# """
# L1 Distance Objective
# """
# struct L1 <: Objective end
#
# function loss(O::L1, ŷ, y)
#   return (1/length(y)) * sum(abs, y .- ŷ)
# end
#
# function lossgrad(O::L1, ŷ, y, x)
#   return (1/length(y)) * sign.(ŷ .- y) * x
# end

# """
# L2 Distance Objective
# """
# struct L2 <: Objective end
#
# function loss(O::L2, y, ŷ)
#   return (1/length(y)) * sum(abs2, y .- ŷ)
# end
#
# function lossgrad(O::L2, y, ŷ, x)
#   return (1/length(y)) * (ŷ .- y) * x
# end
#
# """
# Sigmoid Objective
# """
# struct Sigmoid <: Objective end
#
# function loss(O::Sigmoid, y, ŷ)
#   return 1 ./ (exp(1) .^ (y .- ŷ))
# end
#
# function lossgrad(O::Sigmoid, y, ŷ, x)
#   return Void()
# end


text:
"""
# set_label, set_label!

## Description
Assign a free-text label to something. Labels help to identify objects (in the computer science but not access control sense). But labels are not unique identifiers. Labels can be modified any time.
"""
function set_label end;
function set_label(thing::AbstractThing, label::String)::AbstractThing
    thing.label = label;
    return thing;
end;
function set_label!(thing::AbstractArray{<:AbstractThing}, label::AbstractArray{<:String})::Array{<:AbstractThing}
    broadcast(set_label, thing, label);
    return thing;
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
using TimerOutputs

mutable struct InvalidStoppingCriterion <: SOI.AbstractStoppingCriterion
end

@testset "Stopping Criterion" begin
    #K = 42
    z_LB = 1
    z_UB = 1
    σ = 1

    result = SOI.Result()
    result.upperbound = z_UB
    result.lowerbound = z_LB
    #result.npaths = K
    result.σ_UB = σ

    info = SOI.Info()
    @timeit info.timer "iteration 8" begin
        for i in 1:8
            push!(info.results, result)
            @timeit info.timer SOI.OCUTS_KEY begin end
        end
    end

    @test_throws ErrorException SOI.stop(InvalidStoppingCriterion(), info)
    @test SOI.stop(SOI.AndStoppingCriterion(SOI.IterLimit(8), SOI.CutLimit(8)), info)

    info = SOI.Info()
    for i in 1:7
        push!(info.results, result)
    end
    @timeit info.timer "iteration 7" begin
        for i in 1:2
            @timeit info.timer SOI.OCUTS_KEY begin end
        end
        for i in 1:5
            @timeit info.timer SOI.FCUTS_KEY begin end
        end
    end
    @test !SOI.stop(SOI.AndStoppingCriterion(SOI.IterLimit(8), SOI.CutLimit(8)), info)

    @timeit info.timer "iteration 7" begin
        @timeit info.timer SOI.FCUTS_KEY begin end
    end
    @test !SOI.stop(SOI.AndStoppingCriterion(SOI.IterLimit(8), SOI.CutLimit(8)), info)
end