text:
using LinearAlgebra
function could_be_sperm_tail(tail_length, centerline_points; smoothness=3.0)
    #Make a smoothed spline so we don't get changes from 'noise'
    L = size(centerline_points,1)
    t = curve_lengths(centerline_points)
    spl = ParametricSpline(t,centerline_points',k=2, s=smoothness)
    tspl = range(0, t[end], length=L)
    smoothed = spl(tspl)'

    #plot!(smoothed[:,1], smoothed[:,2])

    prime = derivative(spl,tspl, nu=1)'
    primeprime = derivative(spl,tspl, nu=2)'

    k = (prime[:,1].*primeprime[:,2].-prime[:,2].*primeprime[:,1])./((prime[:,1].^2+prime[:,2].^2).^(3/2))
    cuvature_changes = count(x-> x!=0,(sign.(k)[2:end])-(circshift(sign.(k),1)[2:end]))

    #length check
    length_ok = t[end] <= tail_length+0.1*tail_length && t[end] >= tail_length-0.1*tail_length

    #tail is within +/- 10% estimated length, and curvature changes sign at most once.
    return cuvature_changes, k, length_ok, prime
end

function findMinAvgSubarray(arr, k)
    if length(arr) < k
        #TODO error here!
        return
    end

    result_index = 1;

    current_sum = sum(arr[1:k])

    minimum_sum = current_sum;

    for i = k+1:(length(arr)-k)
        current_sum = current_sum+arr[i] - arr[i-k]
        if current_sum < minimum_sum
            result_index = (i-k+1)
            minimum_sum = current_sum
        end
    end

    return result_index
end

function best_parts(residual, centerline_points, ang, bins, k)
    first_third = Int64(round(size(centerline_points,1)/3))
    second_third = Int64(round(2*size(centerline_points,1)/3))

    head = centerline_points[1:first_third,:]
    mid = centerline_points[first_third+1:second_third,:]
    tail = centerline_points[second_third+1:end,:]

    magnitude = abs.(residual)

    F = Spline1D(bins, magnitude[:,1]);

    projection = [cos(ang) sin(ang)]'
    projected_head_part = (head*projection)[:,1]
    interp_head_residuals = F(projected_head_part)

    projected_mid_part = (mid*projection)[:,1]
    interp_mid_residuals = F(projected_mid_part)

    projected_tail_part = (tail*projection)[:,1]
    interp_tail_residuals = F(projected_tail_part)

    idx_best_head_part_start = findMinAvgSubarray(interp_head_residuals, k)
    idx_best_mid_part_start = findMinAvgSubarray(interp_mid_residuals, k)
    idx_best_tail_part_start = findMinAvgSubarray(interp_tail_residuals, k)
    head_indexes = unique(prepend!(collect(idx_best_head_part_start:idx_best_head_part_start+k-1),1))
    return cat(head[head_indexes,:], mid[idx_best_mid_part_start:idx_best_mid_part_start+k-1,:], tail[idx_best_tail_part_start:idx_best_tail_part_start+k-1,:], dims=1)
end

function estimate_projected_length(projection, ang, r,bins)
    projector = [cos(ang) sin(ang)]'

    #get the 'end points' of the projection, by getting the first and last value where value is greater than tail diameter, which is minimum
    projection_end1 = findfirst(p -> p > 2*r(0.0), projection)
    projection_end2 = findlast(p -> p > 2*r(0.0), projection)

    if projection_end1 == nothing || projection_end2 == nothing
        @warn "less than minimimum projection values"
        #TODO understand why!
        #just return some default value
        return 30
    else
        return abs(bins[projection_end1]-bins[projection_end2])
    end
end
#(detector)   (a)
#|             /|
#|            / |
#|           /  |
#v         (b)_(c)
function keep_best_parts(residual, centerline_points, ang, bins, k, num_points, length, projection, r)
    best = best_parts(residual, centerline_points, ang, bins, k)

    estimated_length = curve_lengths(best)[end]
    ab = length - estimated_length

    fp = parallel_forward(get_outline(best, r)[1], [ang], bins)
    ac = estimate_projected_length(projection, ang, r, bins)-estimate_projected_length(fp[:,1], ang, r, bins)
    if ab > 0.0 && ac > 0.0
        a = best[end,:]
        a_prev = best[end-1,:]
        detector = [cos(ang), sin(ang)]
        ray = [sin(ang), -cos(ang)]
        cplus = a+detector*ac
        cminus = a-detector*ac
        c = norm(cplus-a_prev) > norm(cminus-a_prev) ? cplus : cminus
        bplus = c+ray*ab
        bminus = c-ray*ab
        b = norm(bplus-a_prev) > norm(bminus-a_prev) ? bplus : bminus
        best = cat(best,b', dims=1)
    end
    #create spline
    t = curve_lengths(best)
    spl = ParametricSpline(t,best', k=1, s=1.0)
    tspl = range(0, t[end], length=num_points)
    good_parts_of_tail = spl(tspl)'
    return good_parts_of_tail
end

function translate_points(list_points::Array{T,2}, translation::AbstractArray{T,2}) where T<:AbstractFloat
    return list_points.+translation
end

function get_line(start_point, direction, n)
    return start_point.+direction.*(collect(0:n))
end

function get_projection_ends(projection, r)
    #get the 'end points' of the projection, by getting the first and last value where value is greater than tail diameter, which is minimum
    projection_end1 = findfirst(p -> p > 2*r(0.0), projection)
    projection_end2 = findlast(p -> p > 2*r(0.0), projection)
    return projection_end1, projection_end2
end

function get_straight_template(projection, r, head, ang, num_points,bins)
    projector = [cos(ang) sin(ang)]'
    projected_head = (head*projector)[1,1]

    projection_end1, projection_end2 = get_projection_ends(projection, r)

    #determine the distance from head to each end point (in projection)
    dist_2_head1 = bins[projection_end1]-projected_head
    dist_2_head2 = bins[projection_end2]-projected_head

    #pick the largest distance
    projected_distance = abs(dist_2_head1) > abs(dist_2_head2) ? dist_2_head1 : dist_2_head2

    #Create line segment starting at head, going in direction of detector, and ending where the projection ends.
    return head.+projector'.*collect(range(0.0,projected_distance,length=num_points))
end

function reflection(v,l)
    return 2*(v*l./dot(l,l)).*repeat(l',size(v,1))-v
end

function count_parallel(projection, r; limit=5)
    #get the 'end points' of the projection, by getting the first and last value where value is greater than tail diameter, which is minimum
    projection_end1 = findfirst(p -> p > 2*r(0.0), projection)
    projection_end2 = findlast(p -> p > 2*r(0.0), projection)

    first = min(projection_end1, projection_end2)
    last = max(projection_end1, projection_end2)
    all_parallel = findall(p -> p <= 2.0*r(0.0), projection[first:last])
    neighbor_diff = all_parallel[2:end]-all_parallel[1:end-1]
    changes = count(c -> c > limit, neighbor_diff)+1
    return changes
end

#TODO could check if they are parallel in same place
function try_improvement(best_residual, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)
    residual1 = norm(parallel_forward(get_outline(recon1, r)[1], [ang], bins) - projection)
    residual2 = norm(parallel_forward(get_outline(recon2, r)[1], [ang], bins) - projection)
    #ok1, k1, length_ok1, prime1 = could_be_sperm_tail(tail_length, recon1)
    #ok2, k2, length_ok2, prime2 = could_be_sperm_tail(tail_length, recon2)

    #k = count_parallel(projection[:,1], r, limit = 30)

    # c1 = norm(residual1)+max(0,abs(ok1-1))#+abs(maximum(abs.(prime1))-1.0)
    if residual1 < best_residual#  && length_ok1
        best_residual = residual1#norm(residual1)
        best_recon = recon1
    end
    #
    # c2 = norm(residual2)+max(0,abs(ok2-1))#+abs(maximum(abs.(prime2))-1.0)
    if residual2 < best_residual#  && length_ok2
        best_residual = residual2#norm(residual2)
        best_recon = recon2
    end

    return best_residual, best_recon, residual1, residual2
end

function flip(centerline_points, flip_point, ang)
    reflected = reflection(centerline_points[flip_point:end,:], [cos(ang), sin(ang)])
    first_part = centerline_points[1:flip_point-1,:]
    needed_translation = centerline_points[flip_point,:]-reflected[1,:]
    translated_and_reflected = reflected.+needed_translation'
    return cat(first_part, translated_and_reflected, dims=1)
end

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
# Transformed distributions
struct TransformedDistribution{D, B, V} <: Distribution{V, Continuous} where {D<:Distribution{V, Continuous}, B<:Bijector}
    dist::D
    transform::B

    TransformedDistribution(d::UnivariateDistribution, b::Bijector{0}) = new{typeof(d), typeof(b), Univariate}(d, b)
    TransformedDistribution(d::MultivariateDistribution, b::Bijector{1}) = new{typeof(d), typeof(b), Multivariate}(d, b)
    TransformedDistribution(d::MatrixDistribution, b::Bijector{2}) = new{typeof(d), typeof(b), Matrixvariate}(d, b)
end

const UnivariateTransformed = TransformedDistribution{<:Distribution, <:Bijector, Univariate}
const MultivariateTransformed = TransformedDistribution{<:Distribution, <:Bijector, Multivariate}
const MvTransformed = MultivariateTransformed
const MatrixTransformed = TransformedDistribution{<:Distribution, <:Bijector, Matrixvariate}
const Transformed = TransformedDistribution

"""
    transformed(d::Distribution)
    transformed(d::Distribution, b::Bijector)

Couples distribution `d` with the bijector `b` by returning a `TransformedDistribution`.

If no bijector is provided, i.e. `transformed(d)` is called, then
`transformed(d, bijector(d))` is returned.
"""
transformed(d::Distribution, b::Bijector) = TransformedDistribution(d, b)
transformed(d) = transformed(d, bijector(d))

"""
    bijector(d::Distribution)

Returns the constrained-to-unconstrained bijector for distribution `d`.
"""
bijector(d::Distribution) = DistributionBijector(d)
bijector(d::UnivariateDistribution) = TruncatedBijector(minimum(d), maximum(d))
bijector(d::Normal) = Identity{0}()
bijector(d::MvNormal) = Identity{1}()
bijector(d::PositiveDistribution) = Log{0}()
bijector(d::MvLogNormal) = Log{1}()
bijector(d::SimplexDistribution) = SimplexBijector()
bijector(d::KSOneSided) = Logit(zero(eltype(d)), one(eltype(d)))

bijector_bounded(d, a=minimum(d), b=maximum(d)) = Logit(a, b)
bijector_lowerbounded(d, a=minimum(d)) = Log() ∘ Shift(-a)
bijector_upperbounded(d, b=maximum(d)) = Log() ∘ Shift(b) ∘ Scale(- one(typeof(b)))

const BoundedDistribution = Union{
    Arcsine, Biweight, Cosine, Epanechnikov, Beta, NoncentralBeta
}
bijector(d::BoundedDistribution) = bijector_bounded(d)

const LowerboundedDistribution = Union{Pareto, Levy}
bijector(d::LowerboundedDistribution) = bijector_lowerbounded(d)

##############################
# Distributions.jl interface #
##############################

# size
Base.length(td::Transformed) = length(td.dist)
Base.size(td::Transformed) = size(td.dist)

function logpdf(td::UnivariateTransformed, y::Real)
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

# TODO: implement more efficiently for flows in the case of `Matrix`
function logpdf(td::MvTransformed, y::AbstractMatrix{<:Real})
    # batch-implementation for multivariate
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

function logpdf(td::MvTransformed{<:Dirichlet}, y::AbstractMatrix{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    return logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
end

function _logpdf(td::MvTransformed, y::AbstractVector{<:Real})
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

function _logpdf(td::MvTransformed{<:Dirichlet}, y::AbstractVector{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    return logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
end

# TODO: should eventually drop using `logpdf_with_trans` and replace with
# res = forward(inv(td.transform), y)
# logpdf(td.dist, res.rv) .- res.logabsdetjac
function _logpdf(td::MatrixTransformed, y::AbstractMatrix{<:Real})
    return logpdf_with_trans(td.dist, inv(td.transform)(y), true)
end

# rand
rand(td::UnivariateTransformed) = td.transform(rand(td.dist))
rand(rng::AbstractRNG, td::UnivariateTransformed) = td.transform(rand(rng, td.dist))

# These ovarloadings are useful for differentiating sampling wrt. params of `td.dist`
# or params of `Bijector`, as they are not inplace like the default `rand`
rand(td::MvTransformed) = td.transform(rand(td.dist))
rand(rng::AbstractRNG, td::MvTransformed) = td.transform(rand(rng, td.dist))
# TODO: implement more efficiently for flows
function rand(rng::AbstractRNG, td::MvTransformed, num_samples::Int)
    res = hcat([td.transform(rand(td.dist)) for i = 1:num_samples]...)
    return res
end

function _rand!(rng::AbstractRNG, td::MvTransformed, x::AbstractVector{<:Real})
    rand!(rng, td.dist, x)
    x .= td.transform(x)
end

function _rand!(rng::AbstractRNG, td::MatrixTransformed, x::DenseMatrix{<:Real})
    rand!(rng, td.dist, x)
    x .= td.transform(x)
end

#############################################################
# Additional useful functions for `TransformedDistribution` #
#############################################################
"""
    logpdf_with_jac(td::UnivariateTransformed, y::Real)
    logpdf_with_jac(td::MvTransformed, y::AbstractVector{<:Real})
    logpdf_with_jac(td::MatrixTransformed, y::AbstractMatrix{<:Real})

Makes use of the `forward` method to potentially re-use computation
and returns a tuple `(logpdf, logabsdetjac)`.
"""
function logpdf_with_jac(td::UnivariateTransformed, y::Real)
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

# TODO: implement more efficiently for flows in the case of `Matrix`
function logpdf_with_jac(td::MvTransformed, y::AbstractVector{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

function logpdf_with_jac(td::MvTransformed, y::AbstractMatrix{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

function logpdf_with_jac(td::MvTransformed{<:Dirichlet}, y::AbstractVector{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    lp = logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
    return (lp, res.logabsdetjac)
end

# TODO: should eventually drop using `logpdf_with_trans`
function logpdf_with_jac(td::MatrixTransformed, y::AbstractMatrix{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf_with_trans(td.dist, res.rv, true), res.logabsdetjac)
end

"""
    logpdf_forward(td::Transformed, x)
    logpdf_forward(td::Transformed, x, logjac)

Computes the `logpdf` using the forward pass of the bijector rather than using
the inverse transform to compute the necessary `logabsdetjac`.

This is similar to `logpdf_with_trans`.
"""
# TODO: implement more efficiently for flows in the case of `Matrix`
logpdf_forward(td::Transformed, x, logjac) = logpdf(td.dist, x) - logjac
logpdf_forward(td::Transformed, x) = logpdf_forward(td, x, logabsdetjac(td.transform, x))

function logpdf_forward(td::MvTransformed{<:Dirichlet}, x, logjac)
    T = eltype(x)
    ϵ = _eps(T)

    return logpdf(td.dist, mappedarray(z->z+ϵ, x)) - logjac
end

# forward function
const GLOBAL_RNG = Distributions.GLOBAL_RNG

function _forward(d::UnivariateDistribution, x)
    y, logjac = forward(Identity{0}(), x)
    return (x = x, y = y, logabsdetjac = logjac, logpdf = logpdf.(d, x))
end

forward(rng::AbstractRNG, d::Distribution) = _forward(d, rand(rng, d))
function forward(rng::AbstractRNG, d::Distribution, num_samples::Int)
    return _forward(d, rand(rng, d, num_samples))
end
function _forward(d::Distribution, x)
    y, logjac = forward(Identity{length(size(d))}(), x)
    return (x = x, y = y, logabsdetjac = logjac, logpdf = logpdf(d, x))
end

function _forward(td::Transformed, x)
    y, logjac = forward(td.transform, x)
    return (
        x = x,
        y = y,
        logabsdetjac = logjac,
        logpdf = logpdf_forward(td, x, logjac)
    )
end
function forward(rng::AbstractRNG, td::Transformed)
    return _forward(td, rand(rng, td.dist))
end
function forward(rng::AbstractRNG, td::Transformed, num_samples::Int)
    return _forward(td, rand(rng, td.dist, num_samples))
end

"""
    forward(d::Distribution)
    forward(d::Distribution, num_samples::Int)

Returns a `NamedTuple` with fields `x`, `y`, `logabsdetjac` and `logpdf`.

In the case where `d isa TransformedDistribution`, this means
- `x = rand(d.dist)`
- `y = d.transform(x)`
- `logabsdetjac` is the logabsdetjac of the "forward" transform.
- `logpdf` is the logpdf of `y`, not `x`

In the case where `d isa Distribution`, this means
- `x = rand(d)`
- `y = x`
- `logabsdetjac = 0.0`
- `logpdf` is logpdf of `x`
"""
forward(d::Distribution) = forward(GLOBAL_RNG, d)
forward(d::Distribution, num_samples::Int) = forward(GLOBAL_RNG, d, num_samples)

# utility stuff
params(td::Transformed) = params(td.dist)

# logabsdetjac for distributions
logabsdetjacinv(d::UnivariateDistribution, x::T) where T <: Real = zero(T)
logabsdetjacinv(d::MultivariateDistribution, x::AbstractVector{T}) where {T<:Real} = zero(T)

"""
    logabsdetjacinv(td::UnivariateTransformed, y::Real)
    logabsdetjacinv(td::MultivariateTransformed, y::AbstractVector{<:Real})

Computes the `logabsdetjac` of the _inverse_ transformation, since `rand(td)` returns
the _transformed_ random variable.
"""
logabsdetjacinv(td::UnivariateTransformed, y::Real) = logabsdetjac(inv(td.transform), y)
function logabsdetjacinv(td::MvTransformed, y::AbstractVector{<:Real})
    return logabsdetjac(inv(td.transform), y)
end

ext: jl


text:
squaredigits(num) =
    foldr( (x,accu) -> (s=x^2; muladd(s<10 ? 10 : 100,accu,s)), digits(num); init=0 )

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
export
    Noise,
    GaussianNoise,
    pointsPerOrientation,
    tripletFactor,
    doubletFactor,
    histogramCorrelationsInPicture_alltoall,
    get_noise_volume,
    compton_scattering,
    compton_scattering_q,
    get_compton_noise_volume,
    generateHistogram,
    loadHistograms,
    countTriplets,
    countDoublets,
    detector_to_Ewald_sphere,
    alpha_star_inverse,
    alpha_star,
    ImageVolume,
    getVolumeInterpolated

"""Abstract noise type, to be extended in the future"""
abstract type Noise end

"""Radially symmetric Gaussian noise with sigma(width) and gamma(height)"""
type GaussianNoise <: Noise
    gamma::Float64
    sigma::Float64
    photons::Int64
end

function detector_to_Ewald_sphere(vec::Vector{Float64}, lambda::Float64=1.0)
    q = norm(vec[1:2])
    theta = acos(q*lambda/(4*pi))
    c = cos(theta) #Note that theta should be pi/2 for lambda -> 0
    s = sin(theta)
    return Float64[s*vec[1:2]; c*q]
end

"""Corrects for Ewald sphere curvature"""
function alpha_star_inverse(alpha_star::Float64, k1::Int64, k2::Int64, dq::Float64, lambda::Float64)
    theta1 = acos(k1*dq*lambda/(4*pi))
    theta2 = acos(k2*dq*lambda/(4*pi))
    return acos((cos(alpha_star) - cos(theta1)*cos(theta2))/(sin(theta1)*sin(theta2)))
end

"""Corrects for Ewald sphere curvature"""
function alpha_star(alpha::Float64, k1::Int64, k2::Int64, dq::Float64, lambda::Float64)
    theta1 = acos(k1*dq*lambda/(4*pi))
    theta2 = acos(k2*dq*lambda/(4*pi))
    return acos(sin(theta1)*sin(theta2)*cos(alpha) + cos(theta1)*cos(theta2))
end

"""Interface that allows to treat images as dissections through the 3D Fourier space"""
type ImageVolume <: Volume
    img::Array{Float64,2}
    size::Tuple{Int64,Int64}
    rmax::Float64

    function ImageVolume(img::Array{Float64,2}, rmax::Float64)
        new(img, Base.size(img), rmax)
    end

end

"""Accesing ImageVolume per position"""
function getVolumeInterpolated(img::ImageVolume, pos::Vector{Float64})
    center = Float64[img.size[1]/2.0, img.size[2]/2.0] .+ 0.5
    half_width = center .- 1.0
    shifted_pixel_pos = clamp.(round.(Int64, (pos[1:2]./img.rmax).*half_width + center), 1, minimum(img.size))
    return img.img[shifted_pixel_pos[1], shifted_pixel_pos[2]]
end

function maximum(img::ImageVolume)
    return Base.maximum(img.img)
end

"""
Returns a scattering image for one random orientation given a cubic Intensity volume
@param volume:
@params qcut:
@params number_incident_photons
@params rot:
@params gamma: contribution of noise as ratio to intensity
@params sigma: width of noise
"""
function pointsPerOrientation(volume::Volume, qcut::Float64, envelope_sigma::Float64, number_incident_photons::Int64; incident_photon_variance::Int64 = 0, rot::Matrix{Float64}=random_rotation(3), lambda::Float64=1.0, beamstop_width::Float64=0.0, print_warning::Bool=true)
    maxIntensity = real(maximum(volume))
    if typeof(volume) == SurfaceVolume
        maxIntensity *= 1.3
    end

    photons_per_image = number_incident_photons
    #Introduce a variable number of photons per image
    if incident_photon_variance > 0
        photon_dist = Distributions.Normal(number_incident_photons, incident_photon_variance)
        photons_per_image = round(Int64,Distributions.rand(photon_dist, 1)[1])
    end

    #Generate gaussian distributed incident photons
    dist = Distributions.Normal(0.0, envelope_sigma)
    incidents = Distributions.rand(dist, 3, photons_per_image)

    #Delete z dimension
    incidents[3,:] = 0.0

    #Make array of points
    incidents = Vector{Float64}[ incidents[:,i] for i = 1:photons_per_image]

    #Filter out photons that are in the beamstop
    filter!((p)-> abs(p[1]) > beamstop_width && abs(p[2]) > beamstop_width, incidents)
    #e.g. beamstop_width=qcut/20.0

    #Define envolpe function g(x)
    dist_3D = Distributions.MvNormal(zeros(3), [envelope_sigma,envelope_sigma,envelope_sigma])
    gx = (x) -> pdf(dist_3D, x)
    max_gauss = gx([0.0, 0.0, 0.0])

    #counts how many photons were above the envelope Gaussian
    overflow_count = 0

    accepted = Vector{Float64}[]
    for i = 1:length(incidents)
        p = incidents[i]

        if norm(p) <= qcut
            rp = rot*detector_to_Ewald_sphere(p, lambda)
            intensity_value = real(getVolumeInterpolated(volume, rp))
            ratio = intensity_value/(maxIntensity/max_gauss*gx(rp))
            if ratio > 1.0
                overflow_count += 1
                if print_warning
                    println(STDERR, "Probability ratio is larger than 1.0: ratio:$ratio, norm:$(norm(rp)), incident_photons:$(number_incident_photons), sigma:$(envelope_sigma), qcut:$(qcut).")
                end
            end
            if rand() < ratio
                push!(accepted, p)
            end
        end
    end
    println("Overflow percentage: $(overflow_count/number_incident_photons*100.0)%, scattered photons: $(length(accepted))")
    return accepted,rot
end

"""Uniform random points except beamstop"""
function beamstop_shot(qcut::Float64, photons_per_image::Int64, beamstop_width::Float64=0.0)

    incidents = 2*qcut*(rand(3, photons_per_image) - 0.5)

    #Delete z dimension
    incidents[3,:] = 0.0

    #Make array of points
    incidents = Vector{Float64}[ incidents[:,i] for i = 1:photons_per_image]

    #Filter out photons that are in the beamstop
    filter!((p)-> abs(p[1]) > beamstop_width && abs(p[2]) > beamstop_width, incidents)
    filter!((p)-> norm(p) <= qcut, incidents)

    return incidents
end

"""Calculates, how often a triple of k1,k2,k3 appears in the three photon correlation"""
function tripletFactor(k1::Int64,k2::Int64,k3::Int64)
    if k1 == k2 == k3
        return 1.0
    elseif (k1 == k2 != k3) || (k1 == k3 != k2) || (k2 == k3 != k1)
        return 3.0
    elseif (k1 != k2) && (k1 != k3) &&  (k2 != k3)
        return 6.0
    end
end

"""Calculates, how often a pair of k1,k2 appears in the two photon correlation"""
function doubletFactor(k1::Int64,k2::Int64)
    if k1 != k2
        return 2.0
    else
        return 1.0
    end
end

# """This calculates the 3 photon correlation for a sparse image"""
# function histogramCorrelationsInPicture_random(picture::Array{Array{Float64,1},1}, c2::C2, c3::C3, dq::Float64, N::Int64)
#   da = pi/N
#   l = length(picture)
#   if l < 3 return end
#   indices = collect(1:l)
#
#   for i = 1:Int64(l*(l+1)/2)
#     #Choose random triplet and sort according to k1>k2
#     shuffle!(indices)
#     doublet = Array{Float64}[picture[indices[j]] for j = 1:2 ]
#     sort!(doublet,lt=function(a,b) norm(a)>norm(b) end)
#     p1, p2 = doublet[1], doublet[2]
#
#     #Calculate k1, k2, alpha
#     k1, k2 = round(Int64,norm(p1)/dq), round(Int64,norm(p2)/dq)
#     alpha = angle_between_simple(p1,p2)
#     ai = Int64(mod(floor(Int64, alpha/da),N)+1)
#
#     @inbounds c2[ai,k2,k1] += 1.0 / (doubletFactor(k1,k2)*k1*k2)
#     @inbounds c2[N-ai+1,k2,k1] += 1.0 / (doubletFactor(k1,k2)*k1*k2)
#   end
#
#   for i = 1:Int64(l*(l+1)*(l+2)/6)
#     #Sort them according to k1 > k2 > k3
#     shuffle!(indices)
#     triplet = Array{Float64,1}[picture[indices[j]] for j = 1:3 ]
#     sort!(triplet,lt=function(a,b) norm(a)>norm(b) end)
#     p1,p2,p3 = triplet[1],triplet[2],triplet[3]
#
#     #Calculate k1, k2, k3, alpha, beta
#     k1,k2,k3 = round(Int64,norm(p1)/dq),round(Int64,norm(p2)/dq),round(Int64,norm(p3)/dq)
#     alpha,beta = mod(angle_between(p1,p2), pi),mod(angle_between(p1,p3), pi)
#     ai,bi = Int64(mod(floor(Int64, alpha/da),N)+1), Int64(mod(floor(Int64, beta/da),N)+1)
#
#     #histogram
#     @inbounds c3[ai,bi,k3,k2,k1] += 1.0 / (tripletFactor(k1,k2,k3)*k1*k2*k3)
#     @inbounds c3[N-ai+1,N-bi+1,k3,k2,k1] += 1.0 / (tripletFactor(k1,k2,k3)*k1*k2*k3)
#   end
# end

"""This calculates the 3 photon correlation for a sparse image"""
function histogramCorrelationsInPicture_alltoall(picture::Vector{Vector{Float64}}, c1::C1, c2::C2, c3::C3, dq::Float64, N::Int64, K2::Int64, K3::Int64, lambda::Float64)
    da = pi/N

    #Define picture filter function
    picturefilter = function(p::Vector{Float64})
        k = round(Int64,norm(p)/dq)
        return k >= 1 && k <= K2
    end

    filter!(picturefilter,picture)
    l = length(picture)
    if l < 3 return end
    picture = sort(picture,lt=function(a,b) norm(a)<norm(b) end)

    for i = 1:l
        p1 = picture[i]
        k1 = round(Int64,norm(p1)/dq)

        c1[k1] += Float64(1.0 / k1)

        for j = 1:(i-1)

            p2 = picture[j]
            k2 = round(Int64,norm(p2)/dq)

            alpha = angle_between(p1,p2)
            alpha_s = alpha
            if alpha > pi
                alpha_s = 2*pi - alpha
            end
            ai = Int64(mod(floor(Int64, alpha_s/da),N)+1)

            @fastmath val2 = Float64(1.0 / (doubletFactor(k1,k2)*k1*k2))
            @inbounds c2[ai,k2,k1] += val2

            for k = 1:(j-1) #this implies k2 >= k3

                p3 = picture[k]
                k3 = round(Int64,norm(p3)/dq)

                if k1 <= K3 && k2 <= K3 && k3 <= K3

                    beta = angle_between(p1,p3)
                    if alpha > pi
                        beta = 2*pi - beta
                    end
                    bi = Int64(mod(floor(Int64, beta/da),2*N)+1)

                    @fastmath val3 = Float64(1.0 / (tripletFactor(k1,k2,k3)*k1*k2*k3))
                    @inbounds c3[ai,bi,k3,k2,k1] += val3
                end
            end
        end
    end
end

"""Noise volume in CubeVolume"""
function get_noise_volume(intensity::CubeVolume, sigma::Float64)
    volume = deepcopy(intensity)
    #Let's create a normalized noise volume
    r = linspace(-volume.rmax, volume.rmax, volume.cubesize)
    noise_volume = Float64[ gaussian_distribution([x,y,z], [0.0, 0.0, 0.0], [sigma,sigma,sigma]) for x=r, y=r, z=r]
    noise_volume /= sum(abs, noise_volume)
    volume.cube = noise_volume
    return volume
end

"""Noise volume in SphericalHarmonicsVolume"""
function get_noise_volume(intensity::SurfaceVolume, sigma::Float64)
    intensity_sh = getSphericalHarmonicsVolume(intensity)
    intensity_sh = deleteTerms(intensity_sh, 0, 0)
    for k = 1:intensity.KMAX
        q = k*dr(intensity)
        setc(intensity_sh,k,0,0, gaussian_distribution(q, 0.0, sigma))
    end
    return getSurfaceVolume(intensity_sh)
end

#http://rcwww.kek.jp/research/shield/photon_r.pdf
#http://ac.els-cdn.com/S0969806X15300621/1-s2.0-S0969806X15300621-main.pdf?_tid=1dcbacf4-3a16-11e7-87fd-00000aab0f6c&acdnat=1494925412_b9ac6c901a449780b30c1076c1a57ccc
# """Calculates the compton scattering distribution"""
function compton_scattering(theta::Float64)
    r0_in_m = 2.817940e-15 #classical electron radius [m]
    r0 = r0_in_m * 1.0e10 #classical electron radius [A]
    mec2 = 0.5109989461#electron rest mass energy [MeV] : m_e * c^2
    hnu0 = 0.005#energy of incident photon [MeV] : 5keV
    k0 = hnu0 / mec2
    hnu1 = hnu0 / (1 + k0 * (1 - cos(theta)))#energy of outbound photon
    k1 = hnu1 / mec2
    #Klein and Nishina differential cross section formula
    #From "Compton scattering I: Angular distribution and polarization degree" paper
    #     return 0.25 * r0^2 * (k1/k0)^2 * (k1/k0 + k0/k1 - 2 + 4*sin(theta)^2)
    #From wikipedia
    return 0.5 * r0^2 * (k1/k0)^2 * (k1/k0 + k0/k1 - sin(theta)^2)
end

function compton_scattering_q(q::Float64)
    h = 4.135667516e-15# [eV*s]
    c = 299792458 #[m/s]
    E = 5000 #[eV]
    lambda_in_m = h*c / E #wavelength [m]
    lambda = lambda_in_m * 1.0e10 #wavelength [A]
    val = q * lambda / (4*pi)
    if abs(val) > 1.0
        return 0.0
    else
        theta = asin(val)
        return compton_scattering(theta) * (4*pi/lambda * cos(theta))
    end
end

function get_compton_noise_volume(intensity::CubeVolume)
    volume = deepcopy(intensity)
    #Let's create a normalized noise volume
    r = linspace(-volume.rmax, volume.rmax, volume.cubesize)
    noise_volume = Float64[ compton_scattering_q(norm([x,y,z])) for x=r, y=r, z=r]
    noise_volume /= sum(abs, noise_volume)
    volume.cube = noise_volume
    return volume
end

"""
Generates pictures and histograms the triplets in these pictures without reusing photons
@params intensity:
@params qcut:
@params K:
@params N:
@params maxpictures:
@params number_incident_photons:
@params numprocesses:
@params three:
@params two:
@params gamma:
@params sigma:
"""
function generateHistogram(intensity::Volume; qcut::Float64=1.0, K2::Int64=38, K3::Int64=26, N::Int64=32, max_pictures::Int64=Int64(0), number_incident_photons::Int64=1000, incident_photon_variance::Int64 = 0, numprocesses::Int64=1, file::String="histo.dat", noise::Noise=GaussianNoise(0.0, 1.0, false), batchsize::Int64 = 1000, histogramMethod=histogramCorrelationsInPicture_alltoall, number_particles::Int64=1, lambda::Float64=1.0, beamstop_width::Float64=0.0, beamstop_only::Bool=false)

    #cutoff parameters
    dq = qcut/K2

    c1_full = zeros(Float64, K2)
    c2_full = zeros(Float64,(N, K2,K2))
    c3_full = zeros(Float64,(N,2*N,K3,K3,K3))
    params = Dict("num_pictures"=>0, "num_incident_photons"=>number_incident_photons, "noise"=>noise, "qcut"=>qcut, "K2"=>K2, "K3"=>K3, "N"=>N, "qcut"=>qcut, "dq"=>dq)

    #Load preexisting histogram file so we can continue from where we left off
    if isfile(file)
        println("Continue from preexisting correlation in $file.")
        (params, c2_full, c3_full, c1_full) = deserializeFromFile(file)
    end

    #Calculate noise volume if necessary
    noise_volume = get_noise_volume(intensity, noise.sigma)

    #Count current triplets
    current_triplets = countTriplets(c3_full)

    println("Starting to collect images")
    flush(STDOUT)

    #Make batches of pictures so we can save often
    while max_pictures > 0 && params["num_pictures"] < max_pictures

        c1_part,c2_part,c3_part = @sync @parallel ( (a,b) -> (a[1]+b[1], a[2]+b[2], a[3]+b[3])) for i = 1:numprocesses
            c1 = zeros(Float64,K2)
            c2 = zeros(Float64,N,K2,K2)
            c3 = zeros(Float64,N,2*N,K3,K3,K3)

            for j=1:batchsize
                photon_list = Vector{Float64}[]

                if beamstop_only == false
                    for n=1:number_particles
                        single_molecule,rot = pointsPerOrientation(intensity, qcut, qcut/3.0, number_incident_photons, incident_photon_variance=incident_photon_variance, lambda=lambda, beamstop_width=beamstop_width)

                        if noise.gamma > 0.0
                            noise_photons,_ = pointsPerOrientation(noise_volume,qcut, noise.sigma*1.05, noise.photons, incident_photon_variance=0, rot=rot, lambda=lambda, beamstop_width=beamstop_width)
                            append!(single_molecule, noise_photons)
                        end
                        append!(photon_list, single_molecule)
                    end
                else
                    #beamstop only
                    append!(photon_list, beamstop_shot(qcut, number_incident_photons, beamstop_width))
                end

                histogramMethod(photon_list, c1, c2, c3, dq, N, K2, K3, lambda)
            end
            (c1, c2, c3)
        end
        c1_full += c1_part
        c2_full += c2_part
        c3_full += c3_part

        params["num_pictures"] += numprocesses*batchsize #In numprocesses, we created 'batchsize' pictures
        current_triplets = countTriplets(c3_full)
        serializeToFile(file, (params, c2_full, c3_full, c1_full))

        println("triplets: $current_triplets\tpictures: $(params["num_pictures"]) / $(max_pictures) ($(params["num_pictures"]/max_pictures*100.0) perc.)")
        flush(STDOUT)
    end
end

"Loading a histogram from a file and store two/three photon histograms in global variables."
function loadHistograms(K2::Int64, K3::Int64, file::String, load_c1::Bool=true)

    params, c2_full, c3_full, c1_full = 0,0,0,0
    #Try loading two different structure types
    if load_c1
        params, c2_full, c3_full, c1_full = deserializeFromFile(file)
    else
        params, c2_full, c3_full = deserializeFromFile(file)
    end

    println("Loaded $(countDoublets(c2_full)) doublets and $(countTriplets(c3_full)) triplets from $file generated from $(params["num_pictures"]) pictures.")

    #Doubling old histograms
    NA,NB,K3,_,_ = Base.size(c3_full)
    if NA == NB
        println("Doubling old histograms.")
        new_c3 = zeros(NA, 2*NB,K3,K3,K3)
        for k1=1:K3
            for k2=1:k1
                for k3=1:k2
                    new_c3[:,1:NB,k3,k2,k1] = c3_full[:,:,k3,k2,k1]
                    new_c3[:,NB+1:2*NB,k3,k2,k1] = c3_full[:,:,k3,k2,k1]
                end
            end
        end
        c3_full = new_c3
    end

    c1_full = max.(c1_full, 1e-30)
    c2_full = max.(c2_full, 1e-30)
    c3_full = max.(c3_full, 1e-30)

    #Cut down to designated K
    c2 = c2_full[:,1:K2,1:K2]
    c2 = c2 / sum(abs, c2)
    c3 = c3_full[:,:,1:K3,1:K3,1:K3]
    c3 = c3 / sum(abs, c3)

    if load_c1
        return c2_full, c2, c3_full, c3, c1_full
    else
        return c2_full, c2, c3_full, c3
    end
end

"""Calculates the total number of triplets in a histogram"""
function countTriplets(c3::C3)
    s = size(c3)
    replace_nan!(c3)
    replace_inf!(c3)
    return sum(abs, c3[a,b,k3,k2,k1]*Float64(tripletFactor(k1,k2,k3)*k1*k2*k3) for a=1:s[1],b=1:s[2],k3=1:s[3],k2=1:s[4],k1=1:s[5])
end

"""Calculates the total number of doublets in a histogram"""
function countDoublets(c2::C2)
    s = size(c2)
    replace_nan!(c2)
    replace_inf!(c2)
    return sum(abs, c2[a,k2,k1]*doubletFactor(k1,k2)*k1*k2 for a=1:s[1],k2=1:s[2],k1=1:s[3])
end


text:
goodSymbols = ["Circle",
				"TriangleUp",
				"TriangleDown",
				"Star",
				"Box"]


text:
# Market_risk.jl
module Market_risk

# Market_risk related code here

end # module


text:
"""
simulator(r::Vector,n::Int,nhist::Int,nalleles::Int,total::Int=10000000,tol::Float64=1e-6,count=false)

Modified Gibson and Bruch next reaction algorithm to simulate G state telegraph model (classic telegraph is 2 state)

r = vector of rates in order of state transitions (forward, backward alternating), creation rate, decay rate
n  = number of states - 1, labeled 0,1,...,n
nhist = number of mRNA bins in histogram
total = maximum number of simulator steps
tolerance in histogram change

tau = matrix that keeps track of the putative reaction times of each transition, rows are transitions and columns are alleles
t = time of next reaction and reaction[1] is which reaction and reaction[2] is which allele

"""
function simulatorGM(r::Vector,n::Int,nhist::Int,nalleles::Int,total::Int=10000000,tol::Float64=1e-6,count=false)

	tau = initialize_times(r,n,nalleles)
	mhist,mhist0,m,steps,t,ts,t0,tsample,err = initialize_sim(r,nhist,tol)

	while err > tol && steps < total
		steps += 1
		t,reaction = findmin(tau)
		dt = t-t0
		t0 = t

		update_mhist!(mhist,m,dt,nhist)

		if t-ts > tsample
			err,mhist0 = update_error(mhist,mhist0)
			ts = t
		end
		if reaction[1] <= 2*n
			if isodd(reaction[1])
				gforward!(tau,reaction[1],t,r,reaction[2])
			else
				greverse!(tau,reaction[1],t,r,reaction[2])
			end
		else
			if reaction[1] == 2*n + 1
				m = eject!(tau,reaction[1],t,m,r,reaction[2])
			else
				m = decay!(tau,reaction[1],t,m,r)
			end
		end
	end  # while

	counts = max(sum(mhist),1)
	mhist /= counts

	if count
		return mhist[1:nhist],counts,steps,err
	else
		return mhist[1:nhist]
	end

end

function simulatorGRM(r::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int,total::Int=10000000,tol::Float64=1e-6,count=false)
	tau,Rstep = initialize_times(r,n,nr,nalleles)
	mhist,mhist0,m,steps,t,ts,t0,tsample,err = initialize_sim(r,nhist,tol)
	while err > tol && steps < total
		steps += 1
		t,reaction = findmin(tau)
		dt = t-t0
		t0 = t
		update_mhist!(mhist,m,dt,nhist)
		if t-ts > tsample
			err,mhist0 = update_error(mhist,mhist0)
			ts = t
		end
		if reaction[1] <= 2*n
			if isodd(reaction[1])
				gforward!(tau,Rstep,reaction[1],t,r,n,reaction[2])
			else
				greverse!(tau,reaction[1],t,r,reaction[2])
			end
		else
			if reaction[1] <= 2*n + nr
				if reaction[1] == 2*n + 1
					initiate!(tau,Rstep,reaction[1],t,r,nr,reaction[2])
				else
					rstep!(tau,Rstep,reaction[1],t,r,n,nr,reaction[2])
				end
			else
				if reaction[1] == 2*n + nr + 1
					m = eject!(tau,Rstep,reaction[1],t,m,r,n,nr,reaction[2])
				else
					m = decay!(tau,reaction[1],t,m,r)
				end
			end
		end
	end  # while
	counts = max(sum(mhist),1)
	mhist /= counts
	if count
		return mhist[1:nhist],counts,steps,err
	else
		return mhist[1:nhist]
	end
end

function simulator1(r::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int,total::Int=10000000,tol::Float64=1e-6,count=false)
	tau,taumin,reactmin,Rstep,taudecay = initialize(r,n,nr,nalleles)
	mhist,mhist0,m,steps,t,ts,t0,tsample,err = initialize_sim(r,nhist,tol)
	while err > tol && steps < total
		steps += 1
		t,allele = findmin(taumin)
		reaction = reactmin[allele]
		if taudecay < t
			t = taudecay
			dt = t-t0
			update_mhist!(mhist,m,dt,nhist)
			m,taudecay = decay(t,m,r)
		else
			dt = t-t0
			update_mhist!(mhist,m,dt,nhist)
			if reaction <= 2*n
				if isodd(reaction)
					gforward!(tau,Rstep,reaction,t,r,n,allele)
				else
					greverse!(tau,reaction,t,r,allele)
				end
			else
				if reaction <= 2*n + nr
					if reaction == 2*n + 1
						initiate!(tau,Rstep,reaction,t,r,nr,allele)
					else
						rstep!(tau,Rstep,reaction,t,r,n,nr,allele)
					end
				else
					m,taudecay = eject1!(tau,Rstep,reaction,taudecay,t,m,r,n,nr,allele)
				end
			end
			taumin[allele],reactmin[allele] = findmin(tau[:,allele])
		end
		t0 = t
		if t-ts > tsample
			err,mhist0 = update_error(mhist,mhist0)
			ts = t
		end
	end  # while
	counts = max(sum(mhist),1)
	mhist /= counts
	if count
		return mhist[1:nhist],counts,steps,err
	else
		return mhist[1:nhist]
	end
end

function simulator2(r::Vector,n::Int,nr::Int,nhist::Int,nalleles::Int,total::Int=10000000,tol::Float64=1e-6,count=false)
	tau,taumin,reactmin,Rstep,taudecay = initialize(r,n,nr,nalleles)
	mhist,mhist0,m,steps,t,ts,t0,tsample,err = initialize_sim(r,nhist,tol)
	while err > tol && steps < total
		steps += 1
		t,allele = findmin(taumin)
		reaction = reactmin[allele]
		if taudecay < t
			t = taudecay
			dt = t-t0
			update_mhist!(mhist,m,dt,nhist)
			m,taudecay = decay(t,m,r)
		else
			dt = t-t0
			update_mhist!(mhist,m,dt,nhist)
			if reaction <= 2*n
				greaction!(tau,Rstep,reaction,t,r,n,allele)
			else
				if reaction <= 2*n + nr
					rreaction!(tau,Rstep,reaction,t,r,n,nr,allele)
				else
					m,taudecay = eject1!(tau,Rstep,reaction,taudecay,t,m,r,n,nr,allele)
				end
			end
			taumin[allele],reactmin[allele] = findmin(tau[:,allele])
		end
		t0 = t
		if t-ts > tsample
			err,mhist0 = update_error(mhist,mhist0)
			ts = t
		end
	end  # while
	counts = max(sum(mhist),1)
	mhist /= counts
	if count
		return mhist[1:nhist],counts,steps,err
	else
		return mhist[1:nhist]
	end
end
function greaction!(tau,Rstep,reaction,t,r,n,allele)
	if isodd(reaction)
		gforward!(tau,Rstep,reaction,t,r,n,allele)
	else
		greverse!(tau,reaction,t,r,allele)
	end
end

function rreaction!(tau,Rstep,reaction,t,r,n,nr,allele)
	if reaction == 2*n + 1
		initiate!(tau,Rstep,reaction,t,r,nr,allele)
	else
		rstep!(tau,Rstep,reaction,t,r,n,nr,allele)
	end
end
"""
technical_loss(m,yieldfactor)

Accounts for loss by sampling from a Binomial distribution with probability yieldfactor
"""
function technical_loss(m,yield)
	d = Binomial(m,clamp(yield,0.,1.))
	rand(d)
end

initialize_sim(r,nhist,tol,samplefactor=20.,errfactor=10.) = zeros(nhist+1),ones(nhist+1),0,0,0.,0.,0.,samplefactor/minimum(r),errfactor*tol
update_error(mhist,mhist0) = (norm(mhist/sum(mhist)-mhist0/sum(mhist0),Inf),copy(mhist))
function update_mhist!(mhist,m,dt,nhist)
	if m + 1 <= nhist
		mhist[m+1] += dt
	else
		mhist[nhist+1] += dt
	end
end
"""
initialize_times(r,n,nalleles)

All alleles are initialized to state 0
"""
function initialize_times(r,n,nalleles,m=0)
	tau = fill(Inf,2*n+2,nalleles)
	for n in 1:nalleles
		tau[1,n] = -log(rand())/r[1]
		tau[end,n] = -log(rand())/(m*r[end])
	end
	return tau
end
function initialize_times(r,n,nr,nalleles)
	tau = fill(Inf,2*n+nr+2,nalleles)
	for n in 1:nalleles
		tau[1,n] = -log(rand())/r[1]
	end
	return tau,zeros(Bool,nr,nalleles)
end
function initialize(r,n,nr,nalleles)
tau = fill(Inf,2*n+nr+1,nalleles)
for n in 1:nalleles
	tau[1,n] = -log(rand())/r[1]
end
tau,tau[1,:],ones(Int,nalleles),zeros(Bool,nr,nalleles),Inf
end

"""
gforward!(tau,reaction,t,r,nallele)

update tau matrix for forward transition
"""
function gforward!(tau,reaction,t,r,nallele)
	tau[reaction,nallele] = Inf
	tau[reaction+1,nallele] = -log(rand())/r[reaction+1] + t
	if reaction > 1
		tau[reaction-1,nallele] = Inf
	end
	tau[reaction+2,nallele] = -log(rand())/r[reaction+2] + t
	nothing
end

function gforward!(tau,state,reaction,t,r,n,nallele)
	tau[reaction,nallele] = Inf
	tau[reaction+1,nallele] = -log(rand())/r[reaction+1] + t
	if reaction > 1
		tau[reaction-1,nallele] = Inf
	end
	if reaction == 2*n-1 && state[1,nallele]
			tau[2*n+1,nallele] = Inf
	else
		tau[reaction+2,nallele] = -log(rand())/r[reaction+2] + t
	end
	nothing
end

"""
greverse!(tau,reaction,t,r,nallele)

update tau matrix for reverse state transition
"""
function greverse!(tau,reaction,t,r,nallele)
	tau[reaction,nallele] = Inf
	tau[reaction+1,nallele] = Inf
	tau[reaction-1,nallele]= -log(rand())/r[reaction-1] + t
	if reaction > 2
		tau[reaction-2,nallele] = -log(rand())/r[reaction-2] + t
	end
	nothing
end

"""
eject!(tau,reaction,t,m,r,nallele)

update tau matrix for mRNA ejection

"""
function eject!(tau,reaction,t,m,r,nallele)
	m += 1
	tau[reaction,nallele] = -log(rand())/r[reaction] + t
	set_decay!(tau,reaction,t,m,r)
end

function eject!(tau,state,reaction,t,m,r,n,nr,nallele)
	m += 1
	tau[reaction,nallele] = Inf
	state[nr,nallele] = 0
	if nr == 1
		set_initiate!(tau,t,r,n,nallele)
	elseif state[nr-1,nallele]
		tau[reaction-1,nallele] = -log(rand())/(r[reaction-1])+ t
	end
	set_decay!(tau,reaction,t,m,r)
end
function eject1!(tau,state,reaction,taudecay,t,m,r,n,nr,nallele)
	m += 1
	tau[reaction,nallele] = Inf
	state[nr,nallele] = 0
	if nr == 1
		set_initiate!(tau,t,r,n,nallele)
	elseif state[nr-1,nallele]
		tau[reaction-1,nallele] = -log(rand())/(r[reaction-1])+ t
	end
	m,set_decay(taudecay,t,m)
end

"""
set_decay!(tau,reaction,t,m)

update tau matrix for decay rate

"""

function set_decay!(tau,reaction,t,m,r)
	if m == 1
		tau[reaction+1,1] = -log(rand())/(r[reaction+1])+ t
	else
		tau[reaction+1,1] = (m-1)/m*(tau[reaction+1] - t) + t
	end
	m
end
function set_decay(taudecay,t,m)
	if m == 1
		return -log(rand())/(taudecay)+ t
	else
		return (m-1)/m*(taudecay - t) + t
	end
end
"""
set_initiate!(tau,t,r,n,nallele)

update tau matrix for initiation rate

"""
function set_initiate!(tau,t,r,n,nallele)
	if isfinite(tau[2*n,nallele])
		tau[2*n+1,nallele] = -log(rand())/(r[2*n+1])+ t
	end
end
"""
decay!(tau,reaction,t,m,r)

update tau matrix for decay transition

"""
function decay!(tau,reaction,t,m,r)
	m -= 1
	tau[reaction,1] = -log(rand())/(m*r[reaction]) + t
	m
end

function decay(t,m,r)
	m -= 1
	# tau[reaction,1] = -log(rand())/(m*r[reaction]) + t
	m,-log(rand())/(m*r[end]) + t
end

"""
initiate!(tau,state,reaction,t,r,nr,nallele)

update tau matrix for initiation reaction from G to R

"""
function initiate!(tau,state,reaction,t,r,nr,nallele)
	tau[reaction,nallele] = Inf
	state[1,nallele] = 1
	if nr < 2 || ~state[2,nallele]
		tau[reaction+1,nallele] =  -log(rand())/(r[reaction+1])+ t
	end
end
"""
rstep!(tau,state,reaction,t,r,n,nr,nallele)

update tau matrix for R step transition

"""
function rstep!(tau,state,reaction,t,r,n,nr,nallele)
	step = reaction - 2*n - 1   #current Rstep
	tau[reaction,nallele] = Inf
	state[step,nallele] = 0
	state[step+1,nallele] = 1
	if  step == nr-1 || ~state[step+2,nallele]
		tau[reaction+1,nallele] =  -log(rand())/(r[reaction+1])+ t
	end
	if step > 1 && state[step-1,nallele]
		tau[reaction-1,nallele] =  -log(rand())/(r[reaction-1])+ t
	elseif step == 1
		set_initiate!(tau,t,r,n,nallele)
	end
end
"""
splice!(tau,state,reaction,t,r,n,nr,nallele)

update tau matrix for splice transition

"""
function splice!(tau,state,reaction,t,r,n,nr,nallele)
	tau[reaction,nallele] = Inf
end

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
"""
```
resample(weights::AbstractArray; method::Symbol = :systematic)
```

Reindexing and reweighting samples from a degenerate distribution

### Arguments:
- `weight`: wtsim[:,i]
        the weights of a degenerate distribution.
- `method`: :systematic, :multinomial, or :polyalgo
        the method for resampling

### Output:

- `indx`: the newly assigned indices of parameter draws.
"""
function resample(weights::Vector{Float64}; method::Symbol = :systematic)
    if method == :multinomial
        n_parts = length(weights)
        indx = Vector{Int64}(undef, n_parts)

        # Stores cumulative weights until given index
        cumulative_weights = cumsum(weights/sum(weights))
        offset = rand(n_parts)

        # This could be parallelized
        for i in 1:n_parts
            indx[i] = findfirst(x -> offset[i] < x, cumulative_weights)
        end

        return indx
    elseif method == :systematic
        n_parts = length(weights)

        # Stores cumulative weights until given index
        cumulative_weights = cumsum(weights/sum(weights))
        offset = rand()

        # Function solves where an individual "spoke" lands
        function subsys(i::Int, offset::Float64, n_parts::Int64, start_ind::Int64,
                        cumulative_weights::Vector{Float64})
            threshold = (i - 1 + offset)/n_parts
            range = start_ind:n_parts
            for j in range
                if cumulative_weights[j] > threshold
                    return j
                end
            end
            return 0
        end

        indx = Vector{Int64}(undef, n_parts)
        for i in 1:n_parts
            if i == 1
                indx[i] = subsys(i, offset, n_parts, 1, cumulative_weights)
            else
                indx[i] = subsys(i, offset, n_parts, indx[i-1], cumulative_weights)
            end
        end

        return indx
    elseif method == :polyalgo
        n_parts = length(weights)
        weights = Weights(weights./sum(weights))
        return sample(1:n_parts, weights, n_parts, replace = true)
    else
        throw("Invalid resampler in SMC. Set model setting :resampler_smc to either :systematic, :multinomial, or :polyalgo")
    end
end


text:
module RRTMGP

include("Device.jl")

include(joinpath("optics", "Vmrs.jl"))
include(joinpath("optics", "LookUpTables.jl"))
include(joinpath("optics", "AngularDiscretizations.jl"))
include(joinpath("optics", "AtmosphericStates.jl"))
include(joinpath("optics", "Sources.jl"))
include(joinpath("optics", "Optics.jl"))
include(joinpath("optics", "Fluxes.jl"))
include(joinpath("optics", "GrayUtils.jl"))
include(joinpath("optics", "BCs.jl"))
include(joinpath("optics", "RTE.jl"))

include(joinpath("rte", "RTESolver.jl"))

end # module


text:
export move

function move(x::Array{Float64,2}; by=[0.0, 0.0, 0.0])
    x .+= vec(by)
    return x
end

function move(obj::T; by=[0.0, 0.0, 0.0]) where T <: Union{Shape, PostOpObj}
    # x, v, y, vol, type
    function func(out)
        x = out[:x]
        x = move(x, by=by)
        return repack!(out, [:x], [x])
    end
    if isa(obj, Shape)
        return PostOpObj(obj, func)
    elseif isa(obj, PostOpObj)
        push!(obj.operations, func)
        return obj
    else
        error("Not allowed.")
    end
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
# types.jl

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space

# Structure to store data of polaron solution + other parameters, for each temperature
struct Polaron
    T
    # Mobilities
    Kμ; Hμ; FHIPμ
    # Spring constant and renormalised (phonon-drag) mass
    k; M
    # Osaka free energy components (A,B,C) and total (F). See Hellwarth et al. 1999 PRB Part IV
    A; B; C; F
    # Relaxation time from Kadanoff Boltzmann transport equation
    Tau
    # Raw variational parameters
    v; w
    # Reduced thermodynamic beta
    βred
    # Feynman polaron radius (Schultz), in SI units. Then also the small-alpha asymptotic approx
    rfsi; rfsmallalpha
    # Setup of simulation. These parameters are sent to the function.
    # Alpha = Frohlich alpha
    α
    # Band effective mass
    mb
    # Effective dielectric frequency
    ω
end
Polaron()=Polaron([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]) # structure initialisation


text:
# Configurations and forces
struct NewtonJacobianConfigurationForce{T,Vq,Vu,Vγ,Vb,VI,VIT,Vq0,Vq0T,Vq1,Vq1T,Vu1,Vu1T,Vpr,Vdu} <: NewtonJacobian
    R::SparseMatrixCSC{T,Int}                 # jacobian
    obj_q2::Vector{Vq}                          # obj views
    obj_u1::Vector{Vu}                          # obj views
    obj_γ1::Vector{Vγ}                          # obj views
    obj_b1::Vector{Vb}                          # obj views

    obj_q1q2::Vector{Vq}
    obj_q2q1::Vector{Vq}

    IV::Vector{VI}                          # dynamics -I views [q2, γ1, b1]
    ITV::Vector{VIT}                        # dynamics -I views [q2, γ1, b1] transposed
    q0::Vector{Vq0}                         # dynamics q0 views
    q0T::Vector{Vq0T}                       # dynamics q0 views transposed
    q1::Vector{Vq1}                         # dynamics q1 views
    q1T::Vector{Vq1T}                       # dynamics q1 views transposed
    u1::Vector{Vu1}                         # dynamics u1 views
    u1T::Vector{Vu1T}                       # dynamics u1 views transposed
    reg_pr::Vpr                             # primal regularization views
    reg_du::Vdu                             # dual regularization views
end

function NewtonJacobianConfigurationForce(model::Model, env::Environment, H::Int)
    nq = model.nq # configuration
    nu = model.nu # control
    nw = model.nw # disturbance
    nc = model.nc # contact
    nb = nc * friction_dim(env) # linear friction
    nd = nq + nc + nb # implicit dynamics constraint
    nr = nq + nu + nc + nb# + nd # size of a one-time-step block

    off = 0
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iγ = SizedVector{nc}(off .+ (1:nc)); off += nc # index of the impact γ1
    ib = SizedVector{nb}(off .+ (1:nb)); off += nb # index of the linear friction b1
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2

    iz = vcat(iq, iγ, ib) # index of the IP solver solution [q2, γ1, b1]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    # iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iν = SizedVector{nd}(1:nd) # index of the dynamics lagrange multiplier ν1

    R = spzeros(H * (nr + nd), H * (nr + nd))
    # R = zeros(H * (nr + nd), H * (nr + nd))

    obj_u1  = [view(R, (t - 1) * nr .+ iu, (t - 1) * nr .+ iu) for t = 1:H]
    obj_γ1  = [view(R, (t - 1) * nr .+ iγ, (t - 1) * nr .+ iγ) for t = 1:H]
    obj_b1  = [view(R, (t - 1) * nr .+ ib, (t - 1) * nr .+ ib) for t = 1:H]
    obj_q2  = [view(R, (t - 1) * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H]
    obj_q1q2  = [view(R, (t - 1) * nr .+ iq, t * nr .+ iq) for t = 1:H-1]
    obj_q2q1  = [view(R, t * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H-1]

    IV  = [view(R, CartesianIndex.((t - 1) * nr .+ iz, H * nr + (t - 1) * nd .+ iν)) for t = 1:H]
    ITV = [view(R, CartesianIndex.(H * nr + (t - 1) * nd .+ iν, (t - 1) * nr .+ iz)) for t = 1:H]
    q0  = [view(R, H * nr + (t - 1) * nd .+ iν, (t - 3) * nr .+ iq) for t = 3:H]
    q0T = [view(R, (t - 3) * nr .+ iq, H * nr + (t - 1) * nd .+ iν) for t = 3:H]
    q1  = [view(R, H * nr + (t - 1) * nd .+ iν, (t - 2) * nr .+ iq) for t = 2:H]
    q1T = [view(R, (t - 2) * nr .+ iq, H * nr + (t - 1) * nd .+ iν) for t = 2:H]
    u1  = [view(R, collect(H * nr + (t - 1) * nd .+ iν), collect((t - 1) * nr .+ iu)) for t = 1:H]
    u1T = [view(R, collect((t - 1) * nr .+ iu), collect(H * nr + (t - 1) * nd .+ iν)) for t = 1:H]
    reg_pr = view(R, CartesianIndex.(1:H*nr, 1:H*nr))
    reg_du = view(R, CartesianIndex.(H * nr .+ (1:H * nd), H * nr .+ (1:H * nd)))

    return NewtonJacobianConfigurationForce(
        R, obj_q2, obj_u1, obj_γ1, obj_b1,
        obj_q1q2, obj_q2q1,
        IV, ITV, q0, q0T, q1, q1T, u1, u1T, reg_pr, reg_du)
end

# Configurations
struct NewtonJacobianConfiguration{T,Vq,Vu,VI,VIT,Vq0,Vq0T,Vq1,Vq1T,Vu1,Vu1T,Vpr,Vdu} <: NewtonJacobian
    R::SparseMatrixCSC{T,Int}                 # jacobian
    obj_q2::Vector{Vq}                          # obj views
    obj_u1::Vector{Vu}                          # obj views

    obj_q1q2::Vector{Vq}
    obj_q2q1::Vector{Vq}

    IV::Vector{VI}                          # dynamics -I views [q2]
    ITV::Vector{VIT}                        # dynamics -I views [q2] transposed
    q0::Vector{Vq0}                         # dynamics q0 views
    q0T::Vector{Vq0T}                       # dynamics q0 views transposed
    q1::Vector{Vq1}                         # dynamics q1 views
    q1T::Vector{Vq1T}                       # dynamics q1 views transposed
    u1::Vector{Vu1}                         # dynamics u1 views
    u1T::Vector{Vu1T}                       # dynamics u1 views transposed
    reg_pr::Vpr                             # primal regularization views
    reg_du::Vdu                             # dual regularization views
end

function NewtonJacobianConfiguration(model::Model, env::Environment, H::Int)

    nq = model.nq # configuration
    nu = model.nu # control
    nw = model.nw # disturbance
    nc = model.nc # contact
    nd = nq # implicit dynamics constraint
    nr = nq + nu # size of a one-time-step block

    off = 0
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iz = copy(iq)                                  # index of the IP solver solution [q2]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    iν = SizedVector{nd}(1:nd) # index of the dynamics lagrange multiplier ν1

    R = spzeros(H * (nr + nd), H * (nr + nd))
    # R = zeros(H * (nr + nd), H * (nr + nd))

    obj_u1  = [view(R, (t - 1) * nr .+ iu, (t - 1) * nr .+ iu) for t = 1:H]
    obj_q2  = [view(R, (t - 1) * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H]

    obj_q1q2  = [view(R, (t - 1) * nr .+ iq, t * nr .+ iq) for t = 1:H-1]
    obj_q2q1  = [view(R, t * nr .+ iq, (t - 1) * nr .+ iq) for t = 1:H-1]

    IV  = [view(R, CartesianIndex.((t - 1) * nr .+ iz, H * nr + (t - 1) * nd .+ iν)) for t = 1:H]
    ITV = [view(R, CartesianIndex.(H * nr + (t - 1) * nd .+ iν, (t - 1) * nr .+ iz)) for t = 1:H]
    q0  = [view(R, H * nr + (t - 1) * nd .+ iν, (t - 3) * nr .+ iq) for t = 3:H]
    q0T = [view(R, (t - 3) * nr .+ iq, H * nr + (t - 1) * nd .+ iν) for t = 3:H]
    q1  = [view(R, H * nr + (t - 1) * nd .+ iν, (t - 2) * nr .+ iq) for t = 2:H]
    q1T = [view(R, (t - 2) * nr .+ iq, H * nr + (t - 1) * nd .+ iν) for t = 2:H]
    u1  = [view(R, collect(H * nr + (t - 1) * nd .+ iν), collect((t - 1) * nr .+ iu)) for t = 1:H]
    u1T = [view(R, collect((t - 1) * nr .+ iu), collect(H * nr + (t - 1) * nd .+ iν)) for t = 1:H]
    reg_pr = view(R, CartesianIndex.(1:H*nr, 1:H*nr))
    reg_du = view(R, CartesianIndex.(H * nr .+ (1:H * nd), H * nr .+ (1:H * nd)))

    return NewtonJacobianConfiguration(
        R, obj_q2, obj_u1,
        obj_q1q2, obj_q2q1,
        IV, ITV, q0, q0T, q1, q1T, u1, u1T, reg_pr, reg_du)
end

function NewtonJacobian(model::Model, env::Environment, H::Int;
    mode = :configurationforce)

    if mode == :configurationforce
        return NewtonJacobianConfigurationForce(model, env, H)
    elseif mode == :configuration
        NewtonJacobianConfiguration(model, env, H)
    else
        @error "mode not implemented"
    end
end

function initialize_jacobian!(jac::NewtonJacobian, obj::Objective, H::Int; update_hessian::Bool=true)

    fill!(jac.R, 0.0)
    # fill!(jac.R[301:480,:], 0.0)
    # fill!(jac.R[:,301:480], 0.0)

    # Objective
    update_hessian && hessian!(jac, obj)

    for t = 1:H
        # Implicit dynamics
        jac.IV[t] .-= 1.0
        jac.ITV[t] .-= 1.0
    end

    return nothing
end

function update_jacobian!(jac::NewtonJacobian, im_traj::ImplicitTrajectory, obj::Objective,
    H::Int, β::T, window::Vector{Int}) where T

    for (i, t) in enumerate(window[1:end-2])
        if i >= 3
            jac.q0[i-2]  .+= im_traj.δq0[t]
            jac.q0T[i-2] .+= transpose(im_traj.δq0[t])
        end

        if i >= 2
            jac.q1[i-1]  .+= im_traj.δq1[t]
            jac.q1T[i-1] .+= transpose(im_traj.δq1[t])
        end

        jac.u1[i]  .+= im_traj.δu1[t]
        jac.u1T[i] .+= transpose(im_traj.δu1[t])

        # Dual regularization
        # jac.reg_pr .+= 1.0 * β * im_traj.ip[t].κ[1]
        jac.reg_du .-= β * im_traj.ip[t].κ[1]
    end

    return nothing
end

function jacobian!(jac::NewtonJacobian, im_traj::ImplicitTrajectory, obj::Objective,
        H::Int, β::T, window::Vector{Int}; update_hessian::Bool=true) where T

    initialize_jacobian!(jac, obj, H; update_hessian=update_hessian)
    update_jacobian!(jac, im_traj, obj, H, β, window)

    return nothing
end

function hessian!(jac::NewtonJacobianConfigurationForce, obj::TrackingObjective)
    for t = 1:length(obj.u)
        # Cost function
        jac.obj_q2[t] .+= obj.q[t]
        jac.obj_u1[t] .+= obj.u[t]
        jac.obj_γ1[t] .+= obj.γ[t]
        jac.obj_b1[t] .+= obj.b[t]
    end
end

function hessian!(jac::NewtonJacobianConfiguration, obj::TrackingObjective)
    for t = 1:length(obj.u)
        # Cost function
        jac.obj_q2[t] .+= obj.q[t]
        jac.obj_u1[t] .+= obj.u[t]
    end
end

function hessian!(jac::NewtonJacobianConfigurationForce, obj::TrackingVelocityObjective)
    for t = 1:length(obj.u)
        # Cost function
        jac.obj_q2[t] .+= obj.q[t]
        jac.obj_u1[t] .+= obj.u[t]
        jac.obj_γ1[t] .+= obj.γ[t]
        jac.obj_b1[t] .+= obj.b[t]

        # velocity
        jac.obj_q2[t] .+= obj.v[t]
        t == 1 && continue
        jac.obj_q2[t-1] .+= obj.v[t]
        jac.obj_q1q2[t-1] .-= 1.0*obj.v[t]
        jac.obj_q2q1[t-1] .-= 1.0*obj.v[t]
    end
end

function hessian!(jac::NewtonJacobianConfiguration, obj::TrackingVelocityObjective)
    for t = 1:length(obj.u)
        # Cost function
        jac.obj_q2[t] .+= obj.q[t]
        jac.obj_u1[t] .+= obj.u[t]

        # velocity
        jac.obj_q2[t] .+= obj.v[t]
        t == 1 && continue
        jac.obj_q2[t-1] .+= obj.v[t]
        jac.obj_q1q2[t-1] .-= 1.0*obj.v[t]
        jac.obj_q2q1[t-1] .-= 1.0*obj.v[t]
    end
end


text:
using ADFPEPS
using ArgParse
using CUDA
using Random
using Test
using OMEinsum
using Optim
using Zygote
CUDA.allowscalar(false)

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--tol"
            help = "tol error for vumps"
            arg_type = Float64
            default = 1e-10
        "--maxiter"
            help = "max iterition for vumps"
            arg_type = Int
            default = 10
        "--opiter"
            help = "iterition for optimise"
            arg_type = Int
            default = 1000
        "--f_tol"
            help = "tol error for optimise"
            arg_type = Float64
            default = 1e-10
        "--D"
            help = "ipeps virtual bond dimension"
            arg_type = Int
            required = true
        "--chi"
            help = "vumps virtual bond dimension"
            arg_type = Int
            required = true
        "--folder"
            help = "folder for output"
            arg_type = String
            default = "./data/"
        "--t"
            help = "t"
            arg_type = Float64
            required = true
        "--U"
            help = "U"
            arg_type = Float64
            required = true
        "--mu"
            help = "μ"
            arg_type = Float64
            required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    Random.seed!(100)
    D = parsed_args["D"]
    χ = parsed_args["chi"]
    tol = parsed_args["tol"]
    maxiter = parsed_args["maxiter"]
    opiter = parsed_args["opiter"]
    f_tol = parsed_args["f_tol"]
    folder = parsed_args["folder"]
    t = parsed_args["t"]
    U = parsed_args["U"]
    μ = parsed_args["mu"]
    ipeps, key = init_ipeps(Hubbard(t,U,μ); Ni = 2, Nj = 2, atype = CuArray, folder = folder, D=D, χ=χ, tol=tol, maxiter=maxiter)
    optimiseipeps(ipeps, key; f_tol = f_tol, opiter = opiter, verbose = true)
end

main()


text:
using DynamicAxisWarping, SlidingDistancesBase
using DelimitedFiles, ReverseDiff, DiffResults, Optim
cd(@__DIR__)

data = readdlm("../data/CBF_TRAIN.txt")
c1,c2,c3 = ntuple(i->data[data[:,1] .== i, 2:end], 3)

c1v = copy.(collect(eachrow(c1)))

plot(c1', layout=3, sp=1, legend=false)
plot!(c2', sp=2)
plot!(c3', sp=3)

input = mean(c1v) # Our initial guess will be the Euclidean mean

function soft_dtw_barycenter(input, data, args...; kwargs...)
    function soft_dtw_mean_cost(data)
        function (x)
            mean(data) do d
                soft_dtw_cost(x, d, args...; kwargs...) / lastlength(d)
            end
        end
    end

    costfun = soft_dtw_mean_cost(data)
    # costfun(input)

    cfg = ReverseDiff.GradientConfig(input)
    tape = ReverseDiff.GradientTape(costfun, input)
    ctape = ReverseDiff.compile(tape)
    results = DiffResults.GradientResult(similar(input))

    function fg!(F, G, x)
        if G != nothing
            ReverseDiff.gradient!(results, ctape, x)
            G .= results.derivs[1]
            return results.value
        end
        if F != nothing
            return costfun(x)
        end
    end

    ##
    res = Optim.optimize(
        Optim.only_fg!(fg!),
        input,
        length(input) < 200 ? BFGS() : LBFGS(),
        Optim.Options(
            store_trace = true,
            show_trace = true,
            show_every = 1,
            iterations = 100,
            allow_f_increases = false,
            time_limit = 150,
            x_tol = 1e-3,
            f_tol = 1e-8,
            g_tol = 1e-4,
        ),
    )
    res.minimizer, res
end

##

bc,res = soft_dtw_barycenter(input, c1v)

using Plots, Plots.PlotMeasures
f1 = plot(c1', lab="", axis=false, legend=:bottom)
plot!(input, l=(4, :red), lab="Euclidean mean")
plot!(bc, l=(4, :green), lab="Soft-DTW mean")
f2 = plot(c1[1:3,:]', layout=(1,3), legend=false, axis=false, margin=-5mm)
plot(f1,f2, layout=(2,1))


text:
using Distributions, Random, Plots, LaTeXStrings; pyplot()

seed = 1
N = 100
lamGrid = 0.01:0.01:0.99

theorM(lam) = mean(Uniform(0,2*lam*(1-lam)))
estM(lam) = mean(rand(Uniform(0,2*lam*(1-lam)),N))

function estM(lam,seed)
    Random.seed!(seed)
    estM(lam)
end

trueM = theorM.(lamGrid)
estM0 = estM.(lamGrid)
estMCRN = estM.(lamGrid,seed)

plot(lamGrid,trueM,
	c=:black, label="Expected curve")
plot!(lamGrid,estM0,
	c=:blue, label="No CRN estiamte")
plot!(lamGrid,estMCRN,
	c=:red, label="CRN estimate",
	xlims=(0,1), ylims=(0,0.4), xlabel=L"\lambda", ylabel = "Mean")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
function cc_unique(xtags::Array{tag})
  ccstack = tag[];
  for xtag in xtags
    if( !cyclic_check(xtag, ccstack))
      push!(ccstack, xtag);
      end;
    end;
  return ccstack;
  end;

function cc_bunch(xtags::Array{tag})
  ztags = tag[];
  for x in xtags
    append!(ztags, cc_unique(x.attach));
    end;
  return ztags;
  end;

function cc_forward(atag::tag, depth=4)
  n=0;
  xtags = cc_unique(atag.attach);
  while(n<=depth)
    if cyclic_check(atag, xtags)
      return (term = :false, depth = n);
      end;
    if size(xtags,1) == 0
      return (term = :true, depth = n);
      end;
    xtags = cc_bunch(xtags);
    n += 1;
    end;
  return (term = :unknown, depth = n);
  end;


text:
## test the initialisation of the optimisation

rootdir = "/home/stephane/Science/ALMA/ArrayConfig/GASS"

push!(LOAD_PATH,"$rootdir/master/src")
using GASS

## directory
wdir    = "$rootdir/products"

cd(wdir)

function main(inpfile)
    res= read_cfg(inpfile)
    println(res)
end

main("../master/data/GA_Parameters_O-3.txt.julia")

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
function add_junk_file(scfc, size, path_prefix="")
    new_file = path_prefix * randstring(12)
    filepath = add!(scfc, new_file, UInt64(size))
    mkpath(dirname(filepath))
    open(filepath, "w") do io
        write(io, rand(UInt8, size))
    end
    return filepath
end

@testset "Basics" begin
    mktempdir() do scfc_root
        # Create an cache that keeps a constant size of 10KB,
        # and discards the least-recently-used entries first.
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeConstant(10*1024), DiscardLRU())

        filepath = add_junk_file(scfc, 1024)
        @test isfile(filepath)
        @test hit!(scfc, basename(filepath))
        @test filesize(filepath) == 1024
        @test scfc.total_size == 1024
        first_filepath = filepath

        filepath = add_junk_file(scfc, 512)
        @test isfile(filepath)
        @test filesize(filepath) == 512
        @test scfc.total_size == 1024 + 512

        # Delete that new one and see if it goes down:
        @test delete!(scfc, basename(filepath))
        @test !isfile(filepath)
        @test !hit!(scfc, basename(filepath))
        @test scfc.total_size == 1024

        # Test that deleting somethign twice doesn't do anything:
        @test !delete!(scfc, basename(filepath))
        @test scfc.total_size == 1024

        # Create a bunch more to fill up the cache:
        for idx in 1:9
            filepath = add_junk_file(scfc, 1024)
            @test isfile(filepath)
            @test scfc.total_size == 1024*(idx + 1)
        end

        # Add one more, and ensure that the first file gets removed
        add_junk_file(scfc, 16)
        @test scfc.total_size == 1024*9 + 16
        @test !isfile(first_filepath)

        # Add lots more small ones, ensuring they do not get removed until much later:
        for idx in 1:62
            filepath = add_junk_file(scfc, 16)
            @test isfile(filepath)
            @test scfc.total_size == 1024*9 + 16*(idx + 1)
        end
        add_junk_file(scfc, 16)
        @test scfc.total_size == 1024*10

        # Now, add a big, honking new object that is not perfectly sized;
        # this will clear most of the cache, but a few objects will remain
        filepath = add_junk_file(scfc, 9*1024 + 17)
        @test isfile(filepath)
        @test scfc.total_size == 9*1024 + 17 + 62*16

        # Ensure that the file sizes within the cache are one big guy and 62 small dudes
        @test count([v.size == 9*1024 + 17 for (k, v) in scfc.entries]) == 1
        @test count([v.size == 16 for (k, v) in scfc.entries]) == 62

        # Next, delete the big guy and rebuild the cache:
        rm(filepath; force=true)
        rebuild!(scfc)
        @test scfc.total_size == 62*16

        # Next, try adding files within folders:
        for idx in 1:3
            filepath = add_junk_file(scfc, 1024, "subdir/")
            @test scfc.total_size == idx*1024 + 62*16
        end

        # Now instead of just rebuilding, we'll create an entirely new scfc object and make sure
        # it correctly picks up all files
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeConstant(10*1024), DiscardLRU())
        @test scfc.total_size == 3*1024 + 62*16
        @test count([startswith(k, "subdir/") for (k, v) in scfc.entries]) == 3
        @test count([!startswith(k, "subdir/") for (k, v) in scfc.entries]) == 62

        # This time with a smaller target size, to ensure it runs a shrink immediately.
        # We can't tell exactly what would be deleted here, if the platform we're running on
        # doesn't track `mtime`'s properly; so we'll just assert that the total size satisfies
        # the constraint and call it good:
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeConstant(2*1024), DiscardLRU())
        @test scfc.total_size <= 2*1024

        # Also ensure that trying to cache something too big fails:
        @test_throws ArgumentError add_junk_file(scfc, 100*1024)

        # Test that add!()'ing something that already exists replaces it:
        add_junk_file(scfc, 1*1024)
        add_junk_file(scfc, 1*1024)
        k, e = first(scfc.entries)
        add!(scfc, k, UInt64(100))
        @test scfc.total_size == 1024 + 100
        @test length(scfc.entries) == 2
    end
end

@testset "LRU semantics" begin
    mktempdir() do scfc_root
        # Create an cache that keeps a constant size of 10KB,
        # and discards the least-recently-used entries first.
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeConstant(10*1024), DiscardLRU())

        # Add ten entries
        filepaths = String[]
        for idx in 1:10
            filepath = add_junk_file(scfc, 1024)
            push!(filepaths, filepath)
            @test isfile(filepaths[end])
            @test scfc.total_size == 1024*idx
        end

        # Add one more thing, ensure it was the first-created thing that gets removed:
        filepath = add_junk_file(scfc, 1024)
        push!(filepaths, filepath)
        @test isfile(filepaths[end])
        @test !isfile(popfirst!(filepaths))
        @test scfc.total_size == 10*1024

        # `hit!()` the oldest key now, and ensure that adding another file doesn't kill it,
        # but instead kills the second-to-oldest key:
        @test hit!(scfc, basename(filepaths[1]))
        filepath = add_junk_file(scfc, 1024)
        push!(filepaths, filepath)
        @test isfile(filepaths[end])
        @test isfile(filepaths[1])
        @test !isfile(filepaths[2])
        @test scfc.total_size == 10*1024

        # Ensure that `hit!()`'ing filepaths[2] fails
        @test hit!(scfc, basename(filepaths[2])) == false
    end
end

@testset "LFU semantics" begin
    mktempdir() do scfc_root
        # Create an cache that keeps a constant size of 10KB,
        # and discards the least-frequently-used entries first.
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeConstant(10*1024), DiscardLFU())

        # Add ten entries, and hit them such that the 2nd and 7th entries will get removed first:
        filepaths = String[]
        for idx in 1:10
            filepath = add_junk_file(scfc, 1024)
            push!(filepaths, filepath)
            @test isfile(filepaths[end])
            @test scfc.total_size == 1024*idx
            if idx != 2 && idx != 7
                @test hit!(scfc, basename(filepath))
            end
        end

        # Add two more files, ensure that the 2nd and 7th entries get removed first.
        # Note that here we rely upon the fallback behavior of the LFU discarding LRU
        # entries when the frequency is a tie; otherwise, we might add a single entry
        # and drop the 2nd element, then when we try to add another and drop the 7th
        # we would instead accidentally drop the entry we had just added.  The LFU
        # will instead always drop the 7th element, as it will be less-recently-used.
        add_junk_file(scfc, 1024)
        add_junk_file(scfc, 1024)
        @test !isfile(filepaths[2])
        @test !isfile(filepaths[7])
        @test scfc.total_size == 10*1024
    end
end

@testset "TargetSizeKeepFree" begin
    mktempdir() do scfc_root
        # Keep (current_freespace - 100MB) free, essentiall allowing this scfc to grow up to 100MB
        keep_free = get_disk_freespace(scfc_root) - 100*1024*1024
        scfc = SizeConstrainedFileCache(scfc_root, TargetSizeKeepFree(keep_free), DiscardLRU())

        # Add five entries, checking each time that the disk space is going down roughly as we expect.
        # Note that other activity on the machine can cause this to be a little flaky, we might have
        # to disable this test.
        filepaths = String[]
        starting_freespace = get_disk_freespace(scfc_root)
        for idx in 1:5
            filepath = add_junk_file(scfc, 10*1024*1024)
            push!(filepaths, filepath)
            @test isfile(filepaths[end])
            @test scfc.total_size == 10*1024*1024*idx
            @test abs(Int64(starting_freespace - (get_disk_freespace(scfc_root) + idx*10*1024*1024))) < 5*1024*1024
        end

        # If we now try to add an 80MB file, test that at least the first two chunks we added get removed:
        add_junk_file(scfc, 80*1024*1024)
        @test !isfile(filepaths[1])
        @test !isfile(filepaths[2])
    end
end


text:
@rule typeof(*)(:out, Marginalisation) (m_A::PointMass, m_in::PointMass) = PointMass(mean(m_A) * mean(m_in))

@rule typeof(*)(:out, Marginalisation) (m_A::PointMass{ <: Real }, m_in::GammaDistributionsFamily) = begin
    return GammaShapeRate(shape(m_in), rate(m_in) / mean(m_A))
end

@rule typeof(*)(:out, Marginalisation) (m_A::PointMass{ <: AbstractMatrix }, m_in::F) where { F <: NormalDistributionsFamily } = begin
    A = mean(m_A)
    μ_in, Σ_in = mean_cov(m_in)
    return convert(promote_variate_type(F, NormalMeanVariance), A * μ_in, A * Σ_in * A')
end

#------------------------
# AbstractVector * UnivariateNormalDistributions
#------------------------
# We consider the following updates as a special case of the MatrixVariate * Multivariate updates.
# Namely, Ax = y, where A ∈ R^{nx1}, x ∈ R^1, and y ∈ R^n. In this case, the matrix A
# can be represented by a n-dimensional vector, and x by a scalar. Before computation,
# quantities are converted to their proper dimensions (see situational sketch below).
#
#     | a ~ AbstractVector -> R^{nx1}
#     v  out ~ Multivariate -> R^n
# -->[x]-->
# in1 ~ Univariate -> R^1
@rule typeof(*)(:out, Marginalisation) (m_A::PointMass{ <: AbstractVector }, m_in::UnivariateNormalDistributionsFamily) = begin
    a = mean(m_A)

    μ_in, v_in = mean_var(m_in)

    # TODO: check, do we need correction! here? (ForneyLab does not have any correction in this case)
    # TODO: Σ in this rule is guaranteed to be ill-defined, has rank equal to one and has determinant equal to zero
    μ = μ_in * a
    Σ = mul_inplace!(v_in, a * a')

    return MvNormalMeanCovariance(μ, Σ)
end

#------------------------
# Real * UnivariateNormalDistributions
#------------------------
@rule typeof(*)(:out, Marginalisation) (m_A::PointMass{ <: Real }, m_in::UnivariateNormalDistributionsFamily) = begin
    a = mean(m_A)
    μ_in, v_in = mean_var(m_in)

    return NormalMeanVariance(a*μ_in, a^2*v_in)
end


text:
#=
Replicating Fig 1A in Gjorgjieva et al 2011
See also examples/plastictyFFtest.jl
=#

push!(LOAD_PATH, abspath(@__DIR__,".."))

using Test
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using ProgressMeter
using BenchmarkTools
function oneDSparse(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

##
dt = 0.1E-3
Ttot = 60.0
function ps_trains(rate::R,Δt_ro::R,Ttot::R;
    tstart::R = 0.05,popτ::R=1E6) where R
  post = collect(range(tstart,Ttot; step=inv(rate)))
  pre = post .- Δt_ro
  return S.PSFixedSpiketrain([pre,post],popτ)
end
function ps_trains(rate::R,Δt_ro::R,Δt_oo::R,Ttot::R;
    tstart::R = 0.05,popτ::R=1E6) where R
  _post1 = collect(range(tstart,Ttot-2Δt_oo; step=inv(rate)))
  post = sort(vcat(_post1,_post1.+Δt_oo))
  pre = _post1 .- Δt_ro
  return S.PSFixedSpiketrain([pre,post],popτ)
end
## generic connection pre to post (2<-1)
wstart = sparse([ 0.0 0.0 ; 1000.0 0.0])

function get_weight_change(Δt_ro::R,plasticity::S.PlasticityRule;
   Δt_oo::Union{R,Nothing}=nothing,Ttot::R=80.0,dt::R=0.1E-3,
   rate::R=10.0,tstart::R=0.05,wstart::R=1E3) where R
  ps_all = if isnothing(Δt_oo)
      ps_trains(rate,Δt_ro,Ttot;tstart=tstart)
    else
      ps_trains(rate,Δt_ro,Δt_oo,Ttot;tstart=tstart)
    end
  conn_all = S.ConnectionPlasticityTest(sparse([0.0 0.0 ; wstart 0.0]),
     plasticity)
  pop_all = S.Population(ps_all,(conn_all,ps_all))
  myntw = S.RecurrentNetwork(dt,pop_all)
  rec_weights = S.RecWeights(conn_all,krec,dt,Ttot)
  times = (0:myntw.dt:Ttot)
  S.reset!(rec_weights)
  S.reset!(ps_all)
  S.reset!(conn_all)
  for (k,t) in enumerate(times)
    rec_weights(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
  nsessions = length(ps_all.neurontype.trains[1])
  _wdiff = diff(rec_weights.weights_now[:])
  return sum(_wdiff)/nsessions
end

myplasticity = let τplus = 20E-3,
  τminus = 20E-3
  τx = 20E-3
  τy = 20E-3
  A2plus = 0.1
  A3plus = 0.1
  A2minus = 0.2
  A3minus = 0.2
  (n_post,n_pre) = size(wstart)
  S.PlasticityTriplets(τplus,τminus,τx,τy,A2plus,A3plus,
    A2minus,A3minus,n_post,n_pre)
end

## negative part without second post spike
Δt_ros_minus = range(-80E-3,-1E-3;length=60)
wchange_minus =@showprogress map(
    Δt_ro-> get_weight_change(Δt_ro,myplasticity;rate=1.0),Δt_ros_minus)
plot(Δt_ros_minus,wchange_minus; marker=:circle)

## positive part with different post spikes

Δt_ros_plus = range(1E-3,80E-3;length=60)
wchange_plus1 = let Δt_oo = 100E-3
   @showprogress map(
    Δt_ro-> get_weight_change(Δt_ro,myplasticity;rate=1.0,Δt_oo=Δt_oo),Δt_ros_plus)
end
wchange_plus2 = let Δt_oo = 20E-3
   @showprogress map(
    Δt_ro-> get_weight_change(Δt_ro,myplasticity;rate=1.0,Δt_oo=Δt_oo),Δt_ros_plus)
end
wchange_plus3 = let Δt_oo = 10E-3
   @showprogress map(
    Δt_ro-> get_weight_change(Δt_ro,myplasticity;rate=1.0,Δt_oo=Δt_oo),Δt_ros_plus)
end

plt=plot(Δt_ros_plus,[wchange_plus1 wchange_plus2 wchange_plus3];
   marker=:circle,label=["100ms" "20ms" "10ms"])
plot!(plt, Δt_ros_minus,wchange_minus; marker=:circle,label="")

##

# to visualize a triplet
#=
_ = let _xlims = (0.,0.3)
  _ylims = (0.0 , 0.3)
  plt=plot(leg=false)
  spkdict = S.get_spiketimes_dictionary(rec_spikes)
  yvals = [0.1,0.2]
  for (k,v) in pairs(spkdict)
    scatter!(plt,_->yvals[k],v;marker=:vline,
      color=:white,markersize=11,
      ylims=_ylims,xlims=_xlims)
  end
  plt
end
=#


text:
"""
    age_eff_pop(J, nty, nn)

Generates the age efficieny and population vectors
"""
function age_eff_pop(J, nty, nn)

    # Age-Efficiency Units from Hansen (1993)
    ephansen = zeros(J)

    ephansen[1]=1.0000
    ephansen[2]=1.0719
    ephansen[3]=1.1438
    ephansen[4]=1.2158
    ephansen[5]=1.2842
    ephansen[6]=1.3527
    ephansen[7]=1.4212
    ephansen[8]=1.4897
    ephansen[9]=1.5582
    ephansen[10]=1.6267
    ephansen[11]=1.6952
    ephansen[12]=1.7217
    ephansen[13]=1.7438
    ephansen[14]=1.7748
    ephansen[15]=1.8014
    ephansen[16]=1.8279
    ephansen[17]=1.8545
    ephansen[18]=1.8810
    ephansen[19]=1.9075
    ephansen[20]=1.9341
    ephansen[21]=1.9606
    ephansen[22]=1.9623
    ephansen[23]=1.9640
    ephansen[24]=1.9658
    ephansen[25]=1.9675
    ephansen[26]=1.9692
    ephansen[27]=1.9709
    ephansen[28]=1.9726
    ephansen[29]=1.9743
    ephansen[30]=1.9760
    ephansen[31]=1.9777
    ephansen[32]=1.9700
    ephansen[33]=1.9623
    ephansen[34]=1.9546
    ephansen[35]=1.9469
    ephansen[36]=1.9392
    ephansen[37]=1.9315
    ephansen[38]=1.9238
    ephansen[39]=1.9161
    ephansen[40]=1.9084
    ephansen[41]=1.9007
    ephansen[42]=1.8354
    ephansen[43]=1.7701
    ephansen[44]=1.7048
    ephansen[45]=1.6396

    # Population Numbers from Bell and Miller (2002)
    pop = zeros(J)
    pop[1]=	197316
    pop[2]=	197141
    pop[3]=	196959
    pop[4]=	196770
    pop[5]=	196580
    pop[6]=	196392
    pop[7]=	196205
    pop[8]=	196019
    pop[9]=	195830
    pop[10]=195634
    pop[11]=195429
    pop[12]=195211
    pop[13]=194982
    pop[14]=194739
    pop[15]=194482
    pop[16]=194211
    pop[17]=193924
    pop[18]=193619
    pop[19]=193294
    pop[20]=192945
    pop[21]=192571
    pop[22]=192169
    pop[23]=191736
    pop[24]=191271
    pop[25]=190774
    pop[26]=190243
    pop[27]=189673
    pop[28]=189060
    pop[29]=188402
    pop[30]=187699
    pop[31]=186944
    pop[32]=186133
    pop[33]=185258
    pop[34]=184313
    pop[35]=183290
    pop[36]=182181
    pop[37]=180976
    pop[38]=179665
    pop[39]=178238
    pop[40]=176689
    pop[41]=175009
    pop[42]=173187
    pop[43]=171214
    pop[44]=169064
    pop[45]=166714
    pop[46]=164147
    pop[47]=161343
    pop[48]=158304
    pop[49]=155048
    pop[50]=151604
    pop[51]=147990
    pop[52]=144189
    pop[53]=140180
    pop[54]=135960
    pop[55]=131532
    pop[56]=126888
    pop[57]=122012
    pop[58]=116888
    pop[59]=111506
    pop[60]=105861
    pop[61]=99957
    pop[62]=93806
    pop[63]=87434
    pop[64]=80882
    pop[65]=74204
    pop[66]=67462
    pop[67]=60721
    pop[68]=54053
    pop[69]=47533
    pop[70]=41241
    pop[71]=35259
    pop[72]=29663
    pop[73]=24522
    pop[74]=19890
    pop[75]=15805
    pop[76]=12284
    pop[77]=9331
    pop[78]=6924
    pop[79]=5016
    pop[80]=3550
    pop[81]=2454

    # Labor efficiency
    ep = zeros(nty,J)
    for i in 1:nty
        ep[i,1:J] = ephansen[1:J]
    end

    # WHat is measty? Probability of high and low ability
    measty = zeros(nty)
    measty[1:nty] .= 1.0/nty

    # Survival probabilities: surv(i)=prob(alive in i+1|alive in i)
    surv = zeros(J)
    for i in  1:(J-1)
    surv[i] = pop[i+1]/pop[i]
    end
    surv[J] = 0.0

    # Number of Agents in population
    Nu = zeros(J)
    Nu[1] = 1.0
    for i in 2:J
        Nu[i] = surv[i-1]*Nu[i-1]/(1.0+nn)
    end

    # Fraction of agents in population
    mu = zeros(J)
    for i in 1:J
        mu[i] = Nu[i]/sum(Nu)
    end

    # Total population
    topop=sum(Nu)

    return (ephansen = ephansen, pop = pop, surv = surv, Nu = Nu, mu = mu, ep = ep, measty = measty, topop = topop)
end


text:
#
# This example assumes below packages are available in your environment
#

using StatsModelComparisons, StanSample
using StatsFuns, CSV, Random
using ParetoSmooth
#using StatisticalRethinking

ProjDir = @__DIR__

df = CSV.read(joinpath(ProjDir, "cars_logprob.csv"), DataFrame)

log_lik = Matrix(df)';
n_sam, n_obs = size(log_lik)
lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs);

pwaic = [var(log_lik[:, i]) for i in 1:n_obs];
@show -2(sum(lppd) - sum(pwaic))
println()

@show waic(log_lik)
println()

#Random.seed!(123)
loo, loos, pk = psisloo(log_lik);
@show -2loo

#Random.seed!(123)
ll = reshape(Matrix(df), 50, 1000, 4);
psis_ll = psis(ll);

lwp = deepcopy(ll);
lwp += psis_ll.weights;
lwpt = Matrix(reshape(lwp, 50, 4000)');
loos = reshape(logsumexp(lwpt; dims=1), size(lwpt, 2));

@show loo = sum(loos)
@show 2loo

if isdefined(Main, :StatisticalRethinking)
    pk_plot(pk)
    savefig(joinpath(ProjDir, "pk.png"))

    pk_plot(psis_ll.pareto_k)
    savefig(joinpath(ProjDir, "pareto_k.png"))
end


text:
#int64o-math.jl
doc"""
  `Unums.i64add!(carry, a::ArrayNum, b::ArrayNum)`
  adds the array b into the array a, and increments the carry value if there's
  been a carry event.
"""
@gen_code function i64add!{FSS}(carry::UInt64, b::ArrayNum{FSS}, a::ArrayNum{FSS})
  #this algorithm follows exactly from the elementary school addition algorithm.
  #keep a "carry" variable around and go from least significant to most significant
  #bits.
  l = __cell_length(FSS)
  #initialize the carry variable.
  @code quote
    carryevent::Bool = false  #tracks if a carry event occurred in the last cells
    overflow::Bool = false    #tracks if the addition of the two cells causes overflow.
  end

  for idx = l:-1:1
    @code quote
      @inbounds begin
        (b.a[$idx] += a.a[$idx])     #perform the addition
        overflow = b.a[$idx] < a.a[$idx]
        b.a[$idx] += carryevent * o64
        carryevent = overflow | ((carryevent) & (b.a[$idx] == z64))
      end
    end
  end
  #add in the carry from the most significant segment to the entire carry.
  @code :(carry + carryevent * o64)
end
doc"""
  `Unums.i64add(carry, a::ArrayNum, b::ArrayNum)`
  adds the array b into the array a, and increments the carry value if there's
  been a carry event.
"""
function i64add(carry::UInt64, a::UInt64, b::UInt64)
  result::UInt64 = a + b
  augment::UInt64 = (result < a) * o64
  (carry + augment, result)
end

doc"""
  Unums.add_ubit(value, bit) adds the (zero-indexed) bit to the value.  Returns
  (result, carried)
"""
function add_ubit(value::UInt64, bit::UInt16)
  result_value = value + (t64 >> bit)
  (result_value, result_value < value)
end
@gen_code function add_ubit!{FSS}(value::ArrayNum{FSS}, bit::UInt16)
  @code quote
    cell_index = div(bit, 0x0040) + o16
    added_cell = t64 >> (bit % 0x0040)
    carried = false
  end

  for idx = __cell_length(FSS):-1:1
    @code quote
      @inbounds begin
        oldvalue = value.a[$idx]
        value.a[$idx] += ($idx == cell_index) * added_cell + carried * o64
        carried = (value.a[$idx] < oldvalue)
      end
    end
  end
  @code :(carried)
end

@universal function frac_add_ubit!(u::Unum, bit::UInt16)
  if (FSS < 7)
    (u.fraction, carried) = add_ubit(u.fraction, bit)
  else
    carried = add_ubit!(u.fraction, bit)
  end
  return carried
end

function sub_ubit(value::UInt64, bit::UInt16)
  result_value = value - (t64 >> bit)
  (result_value, result_value > value)
end
@gen_code function sub_ubit!{FSS}(value::ArrayNum{FSS}, bit::UInt16)
  @code quote
    cell_index = div(bit, 0x0040) + o16
    subbed_cell = t64 >> (bit % 0x0040)
    borrowed = false
  end

  for idx = __cell_length(FSS):-1:1
    @code quote
      @inbounds begin
        oldvalue = value.a[$idx]
        value.a[$idx] -= ($idx == cell_index) * subbed_cell + borrowed * o64
        borrowed = (value.a[$idx] > oldvalue)
      end
    end
  end
  @code :(borrowed)
end
@universal function frac_sub_ubit!(u::Unum, bit::UInt16)
  if (FSS < 7)
    (u.fraction, borrowed) = sub_ubit(u.fraction, bit)
  else
    borrowed = sub_ubit!(u.fraction, bit)
  end
  return borrowed
end

doc"""
  `Unums.i64sub(carry, a, b)` performs the calculation that is the result of
  b - a (the reverse order is very important!)
"""
function i64sub(carry::UInt64, a::UInt64, b::UInt64)
  #first do a direct substitution.
  result = b - a
  #then check to see if we lost a bit.
  (carry - ((result > b) * o64), result)
end

doc"""
  `Unums.i64sub!(carry, subtrahend, minuend)` performs the calculation that is
  the result of minuend - subtrahend (the reverse passing order is very important!).
  The reason why the order is reversed is because the minend is likely
  to have to been bitshifted prior to subtraction and it's more sensible to do
  that the value that's going to be copied.
"""
@gen_code function i64sub!{FSS}(carry::UInt64, a::ArrayNum{FSS}, b::ArrayNum{FSS})
  @code quote
    borrowed::Bool = false
    underflowed::Bool = false
  end

  for idx = __cell_length(FSS):-1:1
    @code quote
      @inbounds begin
        #first do the subtraction of the two values.
        a.a[$idx] = b.a[$idx] - a.a[$idx]
        underflowed = a.a[$idx] > b.a[$idx]
        a.a[$idx] -= borrowed * o64
        borrowed = underflowed | ((a.a[$idx] == f64) & borrowed)
      end
    end
  end

  @code :(return carry - borrowed * o64)
end

@gen_code function invert!{FSS}(a::ArrayNum{FSS})
  @code :(borrow = z64)
  for idx = __cell_length(FSS):-1:1
    @code quote
      @inbounds begin
        a.a[$idx] = - a.a[$idx] - borrow
        borrow = ((a.a[$idx] | borrow) != z64) * o64
      end
    end
  end
  @code :(a)
end

#these two functions get inlined by handling a quad load as the individual register.
top_part(n::UInt128) = UInt64((n & 0xFFFF_FFFF_FFFF_FFFF_0000_0000_0000_0000) >> 64)
bottom_part(n::UInt128) = UInt64(n & 0x0000_0000_0000_0000_FFFF_FFFF_FFFF_FFFF)

doc"""
  `Unums.i64mul_simple(a::UInt64, b::UInt64)` returns the 64-bit value that
  is the product of the top 32 bits of both numbers.
"""
function i64mul_simple(a::UInt64, b::UInt64)
  top_part(UInt128(a) * UInt128(b))
end

doc"""
  `Unums.i64mul_extended(a::UInt64, b::UInt64)` returns the 128-bit value
  that is the product of both 64-bit numbers.  This is useful for the FSS = 6 case,
  as well as for breaking apart vaules in the higher FSS case.
"""
function i64mul_extended(a::UInt64, b::UInt64)
  UInt128(a) * UInt128(b)
end

doc"""
  `Unums.i64mul(a::UInt64, b::UInt64, ::Type{Val{FSS}})` performs a fraction-
  multiplication, incorporating addition of both a and b as part of the process.
  This requires knowledge of FSS, as the reporting requiremnets for FSS == 6 are
  distinct from the reporting requirements of FSS < 6.

  Outputs a triple:

  `(carry, fraction, ubit)`
"""
@generated function i64mul{FSS}(a::UInt64, b::UInt64, ::Type{Val{FSS}})
  if FSS < 6
    tmask = mask_top(FSS)
    bmask = mask_bot(FSS)
    quote
      carry = z64
      old_res = i64mul_simple(a, b)
      new_res = old_res + a
      carry += (new_res < old_res) * o64
      old_res = new_res
      new_res = old_res + b
      carry += (new_res < old_res) * o64
      #NB: We don't filter out extraneous bits in the new result here because
      #this function will be used by the division algorithm to store extra bits.
      #filtering must take place on calling i64mul.
      return (carry, new_res, (new_res & $bmask != 0) * UNUM_UBIT_MASK)
    end
  else
    quote
      carry = z64
      res = i64mul_extended(a, b)
      old_res = top_part(res)
      ubit = (bottom_part(res) != z64) * UNUM_UBIT_MASK
      new_res = old_res + a
      carry += (new_res < old_res) * o64
      old_res = new_res
      new_res = old_res + b
      carry += (new_res < old_res) * o64
      return (carry, new_res, ubit)
    end
  end
end

#=
how extended multiplication works.

1) chunk into 64-bit words (already done)
2) multiply as 128-bit numbers
3) stack 128-bit numbers

EG.  2-word x 2-word multiply

A1 A2
B1 B2
      [A2 B2] * CALCULATE TOP (IDX == 4 == N + 2)
   [A1 B2] * CALCULATE FULL, USE TOP (IDX == 5 == N + 1)
   [A2 B1] * CALCULATE FULL, USE TOP
[A1 B1]
[R1 R2]

4-word x 4-word multiply

[A1 A2 A3 A4]
[B1 B2 B3 B4]

                  [A4 B4]
               [A3 B4]
               [A4 B3]
            [A2 B4] * CALCULATE TOP (IDX == 6 == N+2)
            [A3 B3] * CALCULATE TOP
            [A4 B2] * CALCULATE TOP
         [A1 B4] * CALCULATE FULL, USE TOP (IDX == 5 == N + 2)
         [A2 B3] * CALCULATE FULL, USE TOP
         [A3 B2] * CALCULATE FULL, USE TOP
         [A4 B1] * CALCULATE FULL, USE TOP
      [A1 B3]
      [A2 B2]
      [A3 B1]
         [B4] * ADD ADDITIVE PART
         [A4] * ADD ADDITIVE PART
   [A1 B2]
   [A2 B1]
      [B3] * ADD ADDITIVE PART
      [A3] * ADD ADDITIVE PART
[A1 B1]
   [B2] * ADD ADDITIVE PART
   [A2] * ADD ADDITIVE PART
[B1] * ADD ADDITIVE PART
[A1] * ADD ADDITIVE PART
[R1 R2 R3 R4]
=#

doc"""
 `Unums.i64mul!(::ArrayNum, ::ArrayNum)` multiplies two arraynums of equal FSS.
 The result will be the most significant digits of the resulting multiplication.
 Algorithm used is the most basic triangular multiplication, omitting as many
 64-bit integer multiplications as is reasonable to achieve the result (~half)

 This function unrolls itself by being a generated function.
"""
i64mul!{FSS}(a::ArrayNum{FSS}, b::ArrayNum{FSS}) = i64mul!(a, b, Val{__cell_length(FSS)}, Val{true})
#generic version.
@gen_code function i64mul!{FSS, cells, addin}(a::ArrayNum{FSS}, b::ArrayNum{FSS}, ::Type{Val{cells}}, ::Type{Val{addin}})
  (cells >= __cell_length(FSS)) || throw(ArgumentError("multiplication should be calculated on at least as many cells as the arraynum."))
  #set up some internal function variables.
  @code quote
    prev_sum::UInt128 = zero(UInt128)
    accum_sum::UInt128 = zero(UInt128)
    carry::UInt64 = zero(UInt64)

    ubit = (ctz(a) + ctz(b) <= max_fsize(FSS)) * UNUM_UBIT_MASK
  end

  #PHASE 1.  CALCULATE top only...  Note the A indices always start at 2
  idx_sum = cells + 2
  for idx = 2:cells
    jdx = idx_sum - idx #calculate the index for the b arraynum.
    @code quote
      @inbounds accum_sum += i64mul_simple(a.a[$idx], b.a[$jdx])
    end
  end

  #PHASE 2.  CALCULATE FULL, USE TOP
  idx_sum = cells + 1
  for idx = 1:cells
    jdx = idx_sum - idx
    @code quote
      prev_sum = accum_sum  #cache the previous value
      @inbounds accum_sum += i64mul_extended(a.a[$idx], b.a[$jdx])
    end
  end

  #PHASE 3.  CALCULATE FULL, USE ALL.
  for idx_sum = cells:-1:2
    #at the beginning of this group, we need to shift everything over 64 bits,
    #then add in the carry value, then clear it.
    @code quote
      accum_sum >>= 64
      #accum_sum += carry
    end
    for idx = 1:(idx_sum - 1)
      jdx = idx_sum - idx
      @code quote
        prev_sum = accum_sum  #cache the previous value
        @inbounds accum_sum += i64mul_extended(a.a[$idx], b.a[$jdx])
      end
    end
    #PHASE 3.5 ADDITIVE PARTS
    if (addin)
      @code quote
        prev_sum = accum_sum
        @inbounds accum_sum += a.a[$idx_sum]
        prev_sum = accum_sum
        @inbounds accum_sum += b.a[$idx_sum]
      end
    end
    #at the end of each loop, be sure to copy over the lower 64 bits to our
    #destination array.  Magically, we won't need A[idx_sum] anymore, so this is a
    #good place to stash the result.
    @code :(@inbounds a.a[$idx_sum] = bottom_part(accum_sum))
  end

  #PHASE 4.  Transfer the last part of the accumulated sum to the first index in
  #in the arraynum.
  @code quote
    carry = z64
    old_top = new_top = top_part(accum_sum)
  end

  if (addin)
    @code quote
      #additive part for last segment.
      @inbounds new_top += a.a[1]
      carry += (old_top > new_top) * o64
      old_top = new_top
      @inbounds new_top += b.a[1]
      carry += (old_top > new_top) * o64
    end
  end

  @code quote
    #store the last word.
    @inbounds a.a[1] = new_top

    return (carry, ubit)
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
## Axis transformations
#==========================================================================================#

"""
    body_to_stability_axes(coords, α)

Convert coordinates into stability axes with angle ``\\alpha``.
"""
body_to_stability_axes(coords, α :: T) where T <: Real = RotY{T}(α) * coords

"""
    body_to_stability_axes(coords, α)

Convert coordinates into stability axes with angle ``\\alpha``.
"""
stability_to_body_axes(coords, α :: T) where T <: Real = RotY{T}(-α) * coords

"""
    body_to_wind_axes(coords, α, β)

Convert coordinates from body axes to wind axes with angles ``\\alpha,~ \\beta``.
"""
body_to_wind_axes(coords, α :: T, β :: T) where T <: Real = RotZY{T}(β, α) * coords

"""
    body_to_wind_axes(coords, α, β)

Convert coordinates from wind axes to body axes with angles ``\\alpha,~ \\beta``.
"""
wind_to_body_axes(coords, α :: T, β :: T) where T <: Real = RotZY{T}(-α, -β) * coords

## Line methods
#==========================================================================================#

"""
    body_to_wind_axes(line, freestream)

Transform a Line from body to wind axes in a given `Freestream`.
"""
body_to_wind_axes(line :: Line, α :: T, β :: T) where T <: Real = Line(body_to_wind_axes(line.r1, α, β), body_to_wind_axes(line.r2, α, β))

## Reflections and projections
#==========================================================================================#

"""
    reflect_xz(vector)

Reflect the y-coordinate of a given 3-dimensional vector about the x-z plane.
"""
reflect_xz(vector) = SVector(vector[1], -vector[2], vector[3])

"""
    project_yz(vector)

Project a given 3-dimensional vector into the y-z plane.
"""
project_yz(vector) = SVector(0, vector[2], vector[3])

"""
    reflect_yz(line)

Reflect a Line onto the ``x``-``z`` plane of its reference coordinate system.
"""
reflect_xz(line :: Line) = Line((reflect_xz ∘ p2)(line), (reflect_xz ∘ p1)(line))

"""
    reflect_yz(horseshoe)

Reflect a Horseshoe onto the ``x``-``z`` plane of its reference coordinate system.
"""
reflect_xz(horseshoe :: Horseshoe) = (Horseshoe ∘ reflect_xz ∘ bound_leg)(horseshoe)

"""
    reflect_yz(panel :: Panel3D)

Reflect a Panel3D onto the ``x``-``z`` plane of its reference coordinate system.
"""
reflect_xz(panel :: Panel3D) = Panel3D((reflect_xz ∘ p1)(panel), (reflect_xz ∘ p2)(panel), (reflect_xz ∘ p3)(panel), (reflect_xz ∘ p4)(panel))

"""
    project_yz(line)

Project a Line onto the ``y``-``z`` plane of its reference coordinate system.
"""
project_yz(line :: Line) = Line((project_yz ∘ p1)(line), (project_yz ∘ p2)(line))

"""
    stability_flip(vector)

Reflect the x- and z- coordinates of a given 3-dimensional vector about the y-z and x-y planes respectively.
"""
stability_flip(vector) = SVector(-vector[1], vector[2], -vector[3])

max_forks_repo_name: GodotMisogi/AeroMDAO