text:
# quadratic constraints

MOI.supports_constraint(o::Optimizer, ::Type{SQF}, ::Type{<:BOUNDS}) = true

function MOI.add_constraint(o::Optimizer, func::SQF, set::S) where {S <: BOUNDS}
    if func.constant != 0.0
        error("SCIP does not support quadratic constraints with a constant offset.")
    end

    allow_modification(o)

    # affine terms
    linrefs = [VarRef(t.variable.value) for t in func.affine_terms]
    lincoefs = [t.coefficient for t in func.affine_terms]

    # quadratic terms
    quadrefs1 = [VarRef(t.variable_1.value) for t in func.quadratic_terms]
    quadrefs2 = [VarRef(t.variable_2.value) for t in func.quadratic_terms]
    # Divide coefficients by 2 iff they come from the diagonal:
    # Take coef * x * y as-is, but turn coef * x^2 into coef/2 * x^2.
    factor = 1.0 .- 0.5 * (quadrefs1 .== quadrefs2)
    quadcoefs = factor .* [t.coefficient for t in func.quadratic_terms]

    # range
    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    cr = add_quadratic_constraint(o.inner, linrefs, lincoefs,
                                  quadrefs1, quadrefs2, quadcoefs, lhs, rhs)
    ci = CI{SQF, S}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.delete(o::Optimizer, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    allow_modification(o)
    delete!(o.constypes[SQF, S], ConsRef(ci.value))
    delete!(o.reference, cons(o, ci))
    delete(o.inner, ConsRef(ci.value))
    return nothing
end

function MOI.set(o::SCIP.Optimizer, ::MOI.ConstraintSet, ci::CI{SQF,S}, set::S) where {S <: BOUNDS}
    allow_modification(o)

    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    @SCIP_CALL SCIPchgLhsQuadratic(o, cons(o, ci), lhs)
    @SCIP_CALL SCIPchgRhsQuadratic(o, cons(o, ci), rhs)

    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    c = cons(o, ci)

    affterms = AFF_TERM[]
    quadterms = QUAD_TERM[]

    # variables that appear only linearly
    nlin = SCIPgetNLinearVarsQuadratic(o, c)
    linvars = unsafe_wrap(Vector{Ptr{SCIP_VAR}}, SCIPgetLinearVarsQuadratic(o, c), nlin)
    lincoefs = unsafe_wrap(Vector{Float64}, SCIPgetCoefsLinearVarsQuadratic(o, c), nlin)
    for i=1:nlin
        push!(affterms, AFF_TERM(lincoefs[i], VI(ref(o, linvars[i]).val)))
    end

    # variables that appear squared, and linearly
    nquadvarterms = SCIPgetNQuadVarTermsQuadratic(o, c)
    quadvarterms = unsafe_wrap(Vector{SCIP_QUADVARTERM}, SCIPgetQuadVarTermsQuadratic(o, c), nquadvarterms)
    for term in quadvarterms
        vi = VI(ref(o, term.var).val)
        push!(affterms, AFF_TERM(term.lincoef, vi))
        # multiply quadratic coefficients by 2!
        push!(quadterms, QUAD_TERM(2.0 * term.sqrcoef, vi, vi))
    end

    # bilinear terms (pair of different variables)
    nbilinterms = SCIPgetNBilinTermsQuadratic(o, c)
    bilinterms = unsafe_wrap(Vector{SCIP_BILINTERM}, SCIPgetBilinTermsQuadratic(o, c), nbilinterms)
    for term in bilinterms
        # keep coefficients as they are!
        push!(quadterms, QUAD_TERM(term.coef, VI(ref(o, term.var1).val), VI(ref(o, term.var2).val)))
    end

    return SQF(quadterms, affterms, 0.0)
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    lhs = SCIPgetLhsQuadratic(o, cons(o, ci))
    rhs = SCIPgetRhsQuadratic(o, cons(o, ci))
    return from_bounds(S, lhs, rhs)
end

function MOI.get(o::Optimizer, ::MOI.ConstraintPrimal, ci::CI{SQF, S}) where {S <: BOUNDS}
    activity = Ref{Cdouble}()
    @SCIP_CALL SCIPgetActivityQuadratic(o, cons(o, ci), SCIPgetBestSol(o), activity)
    return activity[]
end


text:
using GLMakie

# Observable

x = Node(1.0)
y = Node(2.0)

z = @lift($x .+ $y)

f = Figure(resolution = (1000,600))

ax1 = f[1,1] = Axis(f,
    aspect = 1, targetlimits = BBox(-10,10,-10,10),
    title = "GLMakie & Observable test",
    titlegap = 28, titlesize=20,
    xautolimitmargin = (0,0), xgridwidth = 2, xticklabelsize=10,
    xticks = LinearTicks(20), xticksize = 18,
    yautolimitmargin = (0,0), ygridwidth = 2, yticklabelpad = 14,
    yticklabelsize=10, yticks = LinearTicks(20), yticksize=18
)

vlines!(ax1, [0], linewidth = 1)
hlines!(ax1, [0], linewidth = 1)

slope = Node(1.0)
intercept = Node(0.0)
x = -10:0.01:10
y = @lift($slope .* x .+ $intercept)
line1 = lines!(ax1, x, y, color= :red, linewidth = 2)

println("****")

xlims!(ax1, -10, 10)
ylims!(ax1, -10, 10)

println(f)
println("bleh")

f

# sleep(5)
# for i in 1:100
#     intercept[] = intercept[] + 0.02
#     sleep(0.2)
# end

ext: jl

max_stars_repo_path: GL_playground/Makie_playground 0x01.jl

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_path: GL_playground/Makie_playground 0x01.jl

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_path: GL_playground/Makie_playground 0x01.jl

max_forks_repo_name: leroy-walton/playgrounds


text:
# Descriptive Statistics

#############################
#
#   Location
#
#############################

# Geometric mean
"""
    geomean(a)

Return the geometric mean of a real-valued array.
"""
function geomean(a::RealArray)
    s = 0.0
    n = length(a)
    for i = 1 : n
        @inbounds s += log(a[i])
    end
    return exp(s / n)
end

# Harmonic mean
"""
    harmmean(a)

Return the harmonic mean of a real-valued array.
"""
function harmmean(a::RealArray)
    s = 0.0
    n = length(a)
    for i in 1 : n
        @inbounds s += inv(a[i])
    end
    return n / s
end

# Generalized mean
"""
    genmean(a, p)

Return the generalized/power mean with exponent `p` of a real-valued array,
i.e. ``\left( \frac{1}{n} \sum_{i=1}^n a_i^p \right)^{\frac{1}{p}}``, where `n = length(a)`.
It is taken to be the geometric mean when `p == 0`.
"""
function genmean(a::RealArray, p::Real)
    if p == 0
        return geomean(a)
    end
    s = 0.0
    n = length(a)
    for x in a
        #= At least one of `x` or `p` must not be an int to avoid domain errors when `p` is a negative int.
        We choose `x` in order to exploit exponentiation by squaring when `p` is an int. =#
        @inbounds s += convert(Float64, x)^p
    end
    return (s/n)^(1/p)
end

# Trimmed mean
"""
    trimmean(x, p)

Compute the trimmed mean of `x`, i.e. the mean after removing a
proportion `p` of its highest- and lowest-valued elements.
"""
function trimmean(x::RealArray, p::Real)
    n = length(x)
    n > 0 || error("x can not be empty.")
    0 <= p < 1 || error("p must be non-negative and less than 1.")
    rn = min(round(Int, n * p), n-1)

    sx = sort(x)
    nl = rn >> 1
    nh = (rn - nl)
    s = 0.0
    for i = (1+nl) : (n-nh)
        @inbounds s += sx[i]
    end
    return s / (n - rn)
end

# compute mode, given the range of integer values
"""
    mode(a, [r])

Return the first mode (most common number) of an array, optionally
over a specified range `r`.
"""
function mode{T<:Integer}(a::AbstractArray{T}, r::UnitRange{T})
    isempty(a) && error("mode: input array cannot be empty.")
    len = length(a)
    r0 = r[1]
    r1 = r[end]
    cnts = zeros(Int, length(r))
    mc = 0    # maximum count
    mv = r0   # a value corresponding to maximum count
    for i = 1:len
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
                mv = x
            end
        end
    end
    return mv
end

"""
    modes(a, [r])::Vector

Return all modes (most common numbers) of an array, optionally over a
specified range `r`.
"""
function modes{T<:Integer}(a::AbstractArray{T}, r::UnitRange{T})
    r0 = r[1]
    r1 = r[end]
    n = length(r)
    cnts = zeros(Int, n)
    # find the maximum count
    mc = 0
    for i = 1:length(a)
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
            end
        end
    end
    # find all values corresponding to maximum count
    ms = T[]
    for i = 1:n
        @inbounds if cnts[i] == mc
            push!(ms, r[i])
        end
    end
    return ms
end

# compute mode over arbitrary array
function mode{T}(a::AbstractArray{T})
    isempty(a) && error("mode: input array cannot be empty.")
    cnts = Dict{T,Int}()
    # first element
    mc = 1
    mv = a[1]
    cnts[mv] = 1
    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
                mv = x
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end
    return mv
end

function modes{T}(a::AbstractArray{T})
    isempty(a) && error("modes: input array cannot be empty.")
    cnts = Dict{T,Int}()
    # first element
    mc = 1
    cnts[a[1]] = 1
    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end
    # find values corresponding to maximum counts
    ms = T[]
    for (x, c) in cnts
        if c == mc
            push!(ms, x)
        end
    end
    return ms
end

#############################
#
#   quantile and friends
#
#############################

"""
    percentile(v, p)

Return the `p`th percentile of a real-valued array `v`.
"""
percentile{T<:Real}(v::AbstractArray{T}, p) = quantile(v, p * 0.01)

quantile{T<:Real}(v::AbstractArray{T}) = quantile(v, [.0, .25, .5, .75, 1.0])

"""
    nquantile(v, n)

Return the n-quantiles of a real-valued array, i.e. the values which
partition `v` into `n` subsets of nearly equal size.
"""
nquantile{T<:Real}(v::AbstractArray{T}, n::Integer) = quantile(v, (0:n)/n)

#############################
#
#   Dispersion
#
#############################

# span, i.e. the range minimum(x):maximum(x)
"""
    span(x)

Return the span of an integer array, i.e. the range `minimum(x):maximum(x)`.
"""
span{T<:Integer}(x::AbstractArray{T}) = ((a, b) = extrema(x); a:b)

# Variation coefficient: std / mean
"""
    variation(x, m=mean(x))

Return the coefficient of variation of an array `x`, optionally specifying
a precomputed mean `m`. The coefficient of variation is the ratio of the
standard deviation to the mean.
"""
variation{T<:Real}(x::AbstractArray{T}, m::Real) = stdm(x, m) / m
variation{T<:Real}(x::AbstractArray{T}) = variation(x, mean(x))

# Standard error of the mean: std(a) / sqrt(len)
"""
    sem(a)

Return the standard error of the mean of `a`, i.e. `sqrt(var(a) / length(a))`.
"""
sem{T<:Real}(a::AbstractArray{T}) = sqrt(var(a) / length(a))

# Median absolute deviation
"""
    mad(v)

Compute the median absolute deviation of `v`.
"""
function mad{T<:Real}(v::AbstractArray{T})
    isempty(v) && throw(ArgumentError("mad is not defined for empty arrays"))

    S = promote_type(T, typeof(middle(first(v))))

    mad!(LinAlg.copy_oftype(v, S))
end

"""
    StatsBase.mad!(v, center=median!(v); constant=k)

Compute the maximum absolute deviation (MAD) of `v` about a precomputed center
`center`, overwriting `v` in the process. Using the MAD as a consistent estimator
of the standard deviation requires a scaling factor that depends on the underlying
distribution. For normally distributed data, `k` is chosen as
`1 / quantile(Normal(), 3/4) ≈ 1.4826`, which is used as the default here.
"""
function mad!{T<:Real}(v::AbstractArray{T}, center::Real=median!(v);
                       constant::Real = 1 / (-sqrt(2 * one(T)) * erfcinv(3 * one(T) / 2)))
    for i in 1:length(v)
        @inbounds v[i] = abs(v[i]-center)
    end
    constant * median!(v)
end

# Interquartile range
"""
    iqr(v)

Compute the interquartile range (IQR) of an array, i.e. the 75th percentile
minus the 25th percentile.
"""
iqr{T<:Real}(v::AbstractArray{T}) = (q = quantile(v, [.25, .75]); q[2] - q[1])

#############################
#
#   Z-scores
#
#############################

function _zscore!(Z::AbstractArray, X::AbstractArray, μ::Real, σ::Real)
    # Z and X are assumed to have the same size
    iσ = inv(σ)
    if μ == zero(μ)
        for i = 1 : length(X)
            @inbounds Z[i] = X[i] * iσ
        end
    else
        for i = 1 : length(X)
            @inbounds Z[i] = (X[i] - μ) * iσ
        end
    end
    return Z
end

@ngenerate N typeof(Z) function _zscore!{S,T,N}(Z::AbstractArray{S,N}, X::AbstractArray{T,N}, μ::AbstractArray, σ::AbstractArray)
    # Z and X are assumed to have the same size
    # μ and σ are assumed to have the same size, that is compatible with size(X)
    siz1 = size(X, 1)
    @nextract N ud d->size(μ, d)
    if size(μ, 1) == 1 && siz1 > 1
        @nloops N i d->(d>1 ? (1:size(X,d)) : (1:1)) d->(j_d = ud_d ==1 ? 1 : i_d) begin
            v = (@nref N μ j)
            c = inv(@nref N σ j)
            for i_1 = 1:siz1
                (@nref N Z i) = ((@nref N X i) - v) * c
            end
        end
    else
        @nloops N i X d->(j_d = ud_d ==1 ? 1 : i_d) begin
            (@nref N Z i) = ((@nref N X i) - (@nref N μ j)) / (@nref N σ j)
        end
    end
    return Z
end

function _zscore_chksize(X::AbstractArray, μ::AbstractArray, σ::AbstractArray)
    size(μ) == size(σ) || throw(DimensionMismatch("μ and σ should have the same size."))
    for i=1:ndims(X)
        dμ_i = size(μ,i)
        (dμ_i == 1 || dμ_i == size(X,i)) || throw(DimensionMismatch("X and μ have incompatible sizes."))
    end
end

"""
    zscore!([Z], X, μ, σ)

Compute the z-scores of an array `X` with mean `μ` and standard deviation `σ`.
z-scores are the signed number of standard deviations above the mean that an
observation lies. If a destination array `Z` is provided, the scores are stored
in `Z`, otherwise `X` is overwritten.
"""
function zscore!{ZT<:AbstractFloat,T<:Real}(Z::AbstractArray{ZT}, X::AbstractArray{T}, μ::Real, σ::Real)
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore!(Z, X, μ, σ)
end

function zscore!{ZT<:AbstractFloat,T<:Real,U<:Real,S<:Real}(Z::AbstractArray{ZT}, X::AbstractArray{T},
                                                            μ::AbstractArray{U}, σ::AbstractArray{S})
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore_chksize(X, μ, σ)
    _zscore!(Z, X, μ, σ)
end

zscore!{T<:AbstractFloat}(X::AbstractArray{T}, μ::Real, σ::Real) = _zscore!(X, X, μ, σ)

zscore!{T<:AbstractFloat,U<:Real,S<:Real}(X::AbstractArray{T}, μ::AbstractArray{U}, σ::AbstractArray{S}) =
    (_zscore_chksize(X, μ, σ); _zscore!(X, X, μ, σ))

"""
    zscore(X, [μ, σ])

Compute the z-scores of `X`, optionally specifying a precomputed mean `μ` and
standard deviation `σ`. z-scores are the signed number of standard deviations
above the mean that an observation lies. The mean and standard deviation can be
real numbers or arrays of real numbers.
"""
function zscore{T<:Real}(X::AbstractArray{T}, μ::Real, σ::Real)
    ZT = typeof((zero(T) - zero(μ)) / one(σ))
    _zscore!(Array(ZT, size(X)), X, μ, σ)
end

function zscore{T<:Real,U<:Real,S<:Real}(X::AbstractArray{T}, μ::AbstractArray{U}, σ::AbstractArray{S})
    _zscore_chksize(X, μ, σ)
    ZT = typeof((zero(T) - zero(U)) / one(S))
    _zscore!(Array(ZT, size(X)), X, μ, σ)
end

zscore{T<:Real}(X::AbstractArray{T}) = ((μ, σ) = mean_and_std(X); zscore(X, μ, σ))
zscore{T<:Real}(X::AbstractArray{T}, dim::Int) = ((μ, σ) = mean_and_std(X, dim); zscore(X, μ, σ))

#############################
#
#   entropy and friends
#
#############################

"""
    entropy(p, [b])

Compute the entropy of an array `p`, optionally specifying a real number
`b` such that the entropy is scaled by `1/log(b)`.
"""
function entropy{T<:Real}(p::AbstractArray{T})
    s = zero(T)
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        if pi > z
            s += pi * log(pi)
        end
    end
    return -s
end

function entropy{T<:Real}(p::AbstractArray{T}, b::Real)
    return entropy(p) / log(b)
end

"""
    renyientropy(p, α)

Compute the Rényi (generalized) entropy of order `α` of an array `p`.
"""
function renyientropy{T<:Real, U<:Real}(p::AbstractArray{T}, α::U)
    α < 0 && throw(ArgumentError("Order of Rényi entropy not legal, $(α) < 0."))

    s = zero(T)
    z = zero(T)
    scale = sum(p)

    if α ≈ 0
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s += 1
            end
        end
        s = log(s / scale)
    elseif α ≈ 1
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s -= pi * log(pi)
            end
        end
        s = s / scale
    elseif (isinf(α))
        s = -log(maximum(p))
    else # a normal Rényi entropy
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s += pi ^ α
            end
        end
        s = log(s / scale) / (1 - α)
    end
    return s
end

"""
    crossentropy(p, q, [b])

Compute the cross entropy between `p` and `q`, optionally specifying a real
number `b` such that the result is scaled by `1/log(b)`.
"""
function crossentropy{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))
    s = 0.
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        @inbounds qi = q[i]
        if pi > z
            s += pi * log(qi)
        end
    end
    return -s
end

function crossentropy{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T}, b::Real)
    return crossentropy(p,q) / log(b)
end

"""
    kldivergence(p, q, [b])

Compute the Kullback-Leibler divergence of `q` from `p`, optionally specifying
a real number `b` such that the divergence is scaled by `1/log(b)`.
"""
function kldivergence{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))
    s = 0.
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        @inbounds qi = q[i]
        if pi > z
            s += pi * log(pi / qi)
        end
    end
    return s
end

function kldivergence{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T}, b::Real)
    return kldivergence(p,q) / log(b)
end

#############################
#
#   summary
#
#############################

immutable SummaryStats{T<:AbstractFloat}
    mean::T
    min::T
    q25::T
    median::T
    q75::T
    max::T
end

"""
    summarystats(a)

Compute summary statistics for a real-valued array `a`. Returns a
`SummaryStats` object containing the mean, minimum, 25th percentile,
median, 75th percentile, and maxmimum.
"""
function summarystats{T<:Real}(a::AbstractArray{T})
    m = mean(a)
    qs = quantile(a, [0.00, 0.25, 0.50, 0.75, 1.00])
    R = typeof(convert(AbstractFloat, zero(T)))
    SummaryStats{R}(
        convert(R, m),
        convert(R, qs[1]),
        convert(R, qs[2]),
        convert(R, qs[3]),
        convert(R, qs[4]),
        convert(R, qs[5]))
end

function Base.show(io::IO, ss::SummaryStats)
    println(io, "Summary Stats:")
    @printf(io, "Mean:         %.6f\n", ss.mean)
    @printf(io, "Minimum:      %.6f\n", ss.min)
    @printf(io, "1st Quartile: %.6f\n", ss.q25)
    @printf(io, "Median:       %.6f\n", ss.median)
    @printf(io, "3rd Quartile: %.6f\n", ss.q75)
    @printf(io, "Maximum:      %.6f\n", ss.max)
end

"""
    describe(a)

Pretty-print the summary statistics provided by `summarystats`:
the mean, minimum, 25th percentile, median, 75th percentile, and
maximum.
"""
describe{T<:Real}(a::AbstractArray{T}) = show(summarystats(a))


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

try
    AzureClusterlessHPC == Main.TestCore.AzureClusterlessHPC
catch
    ENV["PARAMETERS"] = joinpath(pwd()[1:end-4], "params.json")
    using AzureClusterlessHPC, PyCall, Test, SyntaxTree, Random, Serialization
end

###################################################################################################
# Concatenate expressions

# Case 1 (no current expressions)
AzureClusterlessHPC.batch_clear()
@test isnothing(AzureClusterlessHPC.__packages__)

expr = :(using Test)
AzureClusterlessHPC.append_batchdef_expression(expr)

@test AzureClusterlessHPC.__packages__.head == :block
@test AzureClusterlessHPC.__packages__.args[1].args[3].args[1].args[1] == :Test

# Case 2 (prior expression)
AzureClusterlessHPC.batch_clear()
@batchdef a = 1     # add some prior expression to AST
@test isnothing(AzureClusterlessHPC.__packages__)

expr = :(using Test)
AzureClusterlessHPC.append_batchdef_expression(expr)

@test AzureClusterlessHPC.__packages__.head == :block
@test AzureClusterlessHPC.__packages__.args[1].args[3].args[1].args[1] == :Test

###################################################################################################
# Expressions

AzureClusterlessHPC.batch_clear()

# Set index
expr = :(1 + 1)
@test expr.args[2] == 1
expr[2] = 2
@test expr.args[2] == 2

# Set index for subexpressions
expr = :((1 + 1), (2 + 2))
@test expr.args[1].args[2] == 1
expr[[1,2]] = 2
@test expr.args[1].args[2] == 2

# Locate symbol
level = []
locations = []
symbols = []
index_symbol = :i

expr = :(
    pmap(i -> hello_world(i, arg1, arg2, arg3; kwargs...), 1:10)
)

num_symbols = 4
AzureClusterlessHPC.locate_symbols!(expr, level, locations, symbols; symbol_exception=index_symbol)
@test length(locations) == num_symbols
@test locations[1] == [2,2,2,2,1,1]
@test locations[2] == [2,2,2,4]
@test locations[3] == [2,2,2,5]
@test locations[4] == [2,2,2,6]
@test symbols == [:kwargs, :arg1, :arg2, :arg3]

# Locate symbol w/o index symbol
level = []
locations = []
symbols = []

num_symbols = 6
AzureClusterlessHPC.locate_symbols!(expr, level, locations, symbols; symbol_exception=nothing)
@test length(locations) == num_symbols
@test locations[1] == [2,1]
@test locations[2] == [2,2,2,2,1,1]
@test locations[3] == [2,2,2,3]
@test locations[4] == [2,2,2,4]
@test locations[5] == [2,2,2,5]
@test locations[6] == [2,2,2,6]
@test symbols == [:i, :kwargs, :i, :arg1, :arg2, :arg3]

# Replace given symbol in AST
expr = :(
    pmap(i -> hello_world(i, arg1, arg2, arg3; kwargs...), 1:10)
)
symbol = :arg1
val = 1

@test expr.args[2].args[2].args[2].args[4] == :arg1
replace_symbol_in_expression!(expr, symbol, val)
@test expr.args[2].args[2].args[2].args[4] == val

# If array in exression in indexed, extract indexed symbol_exception
A = randn(2,2)
expr = :(
    hello_world(A[:, 1])
)
replace_symbol_in_expression!(expr, :A, A)  # replace symbol for A with array
replace_symbol_in_expression!(expr, :(:), :)  # replace symbol for indexing

@test size(expr.args[2].args[1]) == size(A)
eval_symbol_get_index!(expr)
@test size(expr.args[2]) == size(A[:, 1])

# Create macro expansion for single-task expression
AzureClusterlessHPC.batch_clear()
expr_single = :(
    hello_world(arg1, arg2, arg3; kwargs...)
)
task_expr_single = AzureClusterlessHPC.create_single_task_expression_list(expr_single)
@test typeof(task_expr_single) == Expr
@test length(task_expr_single.args) == 7

# Create macro expansion for multi-task expression
AzureClusterlessHPC.batch_clear()
expr_multi = :(
    pmap(i -> hello_world(i, arg1, arg2, arg3; kwargs...), 1:10)
)
task_expr_multi = AzureClusterlessHPC.create_multi_task_expression_list_pmap(expr_multi)
@test typeof(task_expr_multi) == Expr
@test length(task_expr_multi.args) == 6

# Macro expansion for generic expression
task_expr = AzureClusterlessHPC.create_expression_to_submit_batch_job(expr_single)
@test typeof(task_expr) == Expr
@test length(task_expr.args) == 7

task_expr = AzureClusterlessHPC.create_expression_to_submit_batch_job(expr_multi)
@test typeof(task_expr) == Expr
@test length(task_expr.args) == 6

# Submit batch job
AzureClusterlessHPC.batch_clear()
push!(AzureClusterlessHPC.__active_pools__, Dict("pool_id" => 1, "clients" => AzureClusterlessHPC.__clients__[1], "credentials" =>  nothing, "resources" => [nothing]))
function hello_world(name1, name2; kwargs...)
    print("Hello ", name1, " and ", name2)
    return name1
end

name1 = "Bob"
name2 = "Jane"
kwargs = (kw1 = 1, kw2 = 2)

expr = :(pmap(i -> hello_world(name1, name2; kwargs...), 1:10))
task_expr = AzureClusterlessHPC.create_expression_to_submit_batch_job(expr)
bctrl = eval(task_expr)

@test typeof(bctrl) == BatchController
@test typeof(bctrl.output[1]) == BlobFuture
@test length(bctrl.output) == 10
@test bctrl.num_tasks == 10

#######################################################################################################################
# Output

# Write
expr = :(
    function hello_world(arg1, arg2; kwargs...)
        out = arg1 + arg2
        write(IOBuffer(), out)
    end
)

filenames = []
AzureClusterlessHPC.find_output_files!(expr, filenames)
@test length(filenames) == 1
@test typeof(filenames[1]) == String

# Replace "return"
expr = :(
    function hello_world(arg1, arg2; kwargs...)
        out1 = arg1 + arg2
        out2 = arg1 - arg2
        return out1, out2
    end
)
@test expr.args[2].args[7].head == :return

filelist = []
AzureClusterlessHPC.replace_return_with_serialization!(expr, filelist)
linefilter!(expr)
@test expr.args[2].head == :block
@test expr.args[2].args[4].args[1] == :serialize

# Replace "return" for given function name
expr = :(
    function hello_world(arg1, arg2; kwargs...)
        out1 = arg1 + arg2
        out2 = arg1 - arg2
        return out1, out2
    end
)
@test expr.args[2].args[7].head == :return

fname = :hello_world
filelist = []
AzureClusterlessHPC.find_function_in_ast_and_replace_return!(expr, fname, filelist)
linefilter!(expr)

@test length(filelist) == 2
@test typeof(filelist[1]) == String
@test expr.args[2].head == :block
@test expr.args[2].args[4].args[1] == :serialize

#######################################################################################################################
# Bcast

# Broadcasting
AzureClusterlessHPC.batch_clear()

expr = :(A)
bcast_expr = AzureClusterlessHPC.bcast_expression(expr)
linefilter!(bcast_expr)

@test bcast_expr.head == :escape
@test bcast_expr.args[1].head == :block
@test length(bcast_expr.args[1].args) == 5


text:
	include("JUL.auxiliar.jl") ;
	include("JUL.routines.jl") ;


text:
function correlate_eigenmodes(U, M, V; tracked_modes=1:size(U,1), rtol=1e-3, atol=1e-3)

    # get number of (not necessarily unique) eigenvalues
    nev = size(U,1)

    # construct the correlation matrix
    C = U*M*V

    # get the row permutation that puts maximum values on the diagonals
    rval = Vector{real(eltype(C))}(undef, nev) # magnitude of values in chosen row
    cval = Vector{real(eltype(C))}(undef, nev) # magnitude of values in chosen column
    perm = zeros(Int, nev) # permutation
    corruption = Vector{real(eltype(C))}(undef, nev) # corruption index for each correlation

    # tracked modes
    for i in tracked_modes

        # get correlation matrix entry magnitudes for this row
        for j = 1:nev
            rval[j] = abs(C[i,j])
        end

        # rank modes correlations from best to worst for this row
        ranked_modes = sortperm(rval, rev=true)

        # choose the best unassigned mode correlation for this row
        j1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]

        # get correlation matrix entry magnitudes for the chosen column
        for k = 1:nev
            cval[k] = abs(C[k,j1])
        end

        # find approximate multiplicity of left and/or right eigenvector
        multiplicity = max(
            count(x -> isapprox(x, rval[j1]; rtol, atol), rval),
            count(x -> isapprox(x, cval[i]; rtol, atol), cval)
        )

        if multiplicity < nev
            # determine n+1 best unassigned fit, where n is the approximate multiplicity
            if j1 in ranked_modes[1:multiplicity]
                j2 = ranked_modes[multiplicity+1]
            else
                j2 = ranked_modes[1]
            end

            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = rval[j2]/rval[j1]
        else
            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = 0.0
        end
    end

    # remaining modes
    for i = 1:nev

        if !(i in tracked_modes)

            # get correlation matrix entry magnitudes for this row
            for j = 1:nev
                rval[j] = abs(C[i,j])
            end

            # rank modes correlations from best to worst for this row
            ranked_modes = sortperm(rval, rev=true)

            # choose the best unassigned mode correlation for this row
            j1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]

            # get correlation matrix entry magnitudes for the chosen column
            for k = 1:nev
                cval[k] = abs(C[k,j1])
            end

            # find approximate multiplicity of left and/or right eigenvector
            multiplicity = max(
                count(x -> isapprox(x, rval[j1]; rtol, atol), rval),
                count(x -> isapprox(x, cval[i]; rtol, atol), cval)
            )

            # determine n+1 best unassigned fit, where n is the approximate multiplicity
            if j1 in ranked_modes[1:multiplicity]
                j2 = ranked_modes[multiplicity+1]
            else
                j2 = ranked_modes[1]
            end

            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = rval[j2]/rval[j1]
        end

    end

    return perm, corruption
end


text:
function APMFile()
    return APMFile(())
end

function add_maneuver(obj::APMFile, arg0::APMFile_APMManeuver)
    return jcall(obj, "addManeuver", void, (APMFile_APMManeuver,), arg0)
end

function get_attitude_quaternion_direction(obj::APMFile)
    return jcall(obj, "getAttitudeQuaternionDirection", JString, ())
end

function get_epoch(obj::APMFile)
    return jcall(obj, "getEpoch", AbsoluteDate, ())
end

function get_epoch_comment(obj::APMFile)
    return jcall(obj, "getEpochComment", List, ())
end

function get_euler_comment(obj::APMFile)
    return jcall(obj, "getEulerComment", List, ())
end

function get_euler_direction(obj::APMFile)
    return jcall(obj, "getEulerDirection", JString, ())
end

function get_euler_frame_a_string(obj::APMFile)
    return jcall(obj, "getEulerFrameAString", JString, ())
end

function get_euler_frame_b_string(obj::APMFile)
    return jcall(obj, "getEulerFrameBString", JString, ())
end

function get_euler_rot_seq(obj::APMFile)
    return jcall(obj, "getEulerRotSeq", JString, ())
end

function get_has_maneuver(obj::APMFile)
    return jcall(obj, "getHasManeuver", jboolean, ())
end

function get_i11(obj::APMFile)
    return jcall(obj, "getI11", jdouble, ())
end

function get_i12(obj::APMFile)
    return jcall(obj, "getI12", jdouble, ())
end

function get_i13(obj::APMFile)
    return jcall(obj, "getI13", jdouble, ())
end

function get_i22(obj::APMFile)
    return jcall(obj, "getI22", jdouble, ())
end

function get_i23(obj::APMFile)
    return jcall(obj, "getI23", jdouble, ())
end

function get_i33(obj::APMFile)
    return jcall(obj, "getI33", jdouble, ())
end

function get_inertia_ref_frame_string(obj::APMFile)
    return jcall(obj, "getInertiaRefFrameString", JString, ())
end

function get_maneuver(obj::APMFile, arg0::jint)
    return jcall(obj, "getManeuver", APMFile_APMManeuver, (jint,), arg0)
end

function get_maneuvers(obj::APMFile)
    return jcall(obj, "getManeuvers", List, ())
end

function get_meta_data(obj::APMFile)
    return jcall(obj, "getMetaData", ADMMetaData, ())
end

function get_meta_data_comment(obj::APMFile)
    return jcall(obj, "getMetaDataComment", List, ())
end

function get_nb_maneuvers(obj::APMFile)
    return jcall(obj, "getNbManeuvers", jint, ())
end

function get_nutation(obj::APMFile)
    return jcall(obj, "getNutation", jdouble, ())
end

function get_nutation_period(obj::APMFile)
    return jcall(obj, "getNutationPeriod", jdouble, ())
end

function get_nutation_phase(obj::APMFile)
    return jcall(obj, "getNutationPhase", jdouble, ())
end

function get_quaternion(obj::APMFile)
    return jcall(obj, "getQuaternion", Quaternion, ())
end

function get_quaternion_dot(obj::APMFile)
    return jcall(obj, "getQuaternionDot", Quaternion, ())
end

function get_quaternion_frame_a_string(obj::APMFile)
    return jcall(obj, "getQuaternionFrameAString", JString, ())
end

function get_quaternion_frame_b_string(obj::APMFile)
    return jcall(obj, "getQuaternionFrameBString", JString, ())
end

function get_rate_frame_string(obj::APMFile)
    return jcall(obj, "getRateFrameString", JString, ())
end

function get_rotation_angles(obj::APMFile)
    return jcall(obj, "getRotationAngles", Vector3D, ())
end

function get_rotation_rates(obj::APMFile)
    return jcall(obj, "getRotationRates", Vector3D, ())
end

function get_spacecraft_comment(obj::APMFile)
    return jcall(obj, "getSpacecraftComment", List, ())
end

function get_spin_alpha(obj::APMFile)
    return jcall(obj, "getSpinAlpha", jdouble, ())
end

function get_spin_angle(obj::APMFile)
    return jcall(obj, "getSpinAngle", jdouble, ())
end

function get_spin_angle_vel(obj::APMFile)
    return jcall(obj, "getSpinAngleVel", jdouble, ())
end

function get_spin_comment(obj::APMFile)
    return jcall(obj, "getSpinComment", List, ())
end

function get_spin_delta(obj::APMFile)
    return jcall(obj, "getSpinDelta", jdouble, ())
end

function get_spin_direction(obj::APMFile)
    return jcall(obj, "getSpinDirection", JString, ())
end

function get_spin_frame_a_string(obj::APMFile)
    return jcall(obj, "getSpinFrameAString", JString, ())
end

function get_spin_frame_b_string(obj::APMFile)
    return jcall(obj, "getSpinFrameBString", JString, ())
end

function set_attitude_quaternion_direction(obj::APMFile, arg0::JString)
    return jcall(obj, "setAttitudeQuaternionDirection", void, (JString,), arg0)
end

function set_epoch(obj::APMFile, arg0::AbsoluteDate)
    return jcall(obj, "setEpoch", void, (AbsoluteDate,), arg0)
end

function set_epoch_comment(obj::APMFile, arg0::List)
    return jcall(obj, "setEpochComment", void, (List,), arg0)
end

function set_euler_comment(obj::APMFile, arg0::List)
    return jcall(obj, "setEulerComment", void, (List,), arg0)
end

function set_euler_direction(obj::APMFile, arg0::JString)
    return jcall(obj, "setEulerDirection", void, (JString,), arg0)
end

function set_euler_frame_a_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setEulerFrameAString", void, (JString,), arg0)
end

function set_euler_frame_b_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setEulerFrameBString", void, (JString,), arg0)
end

function set_euler_rot_seq(obj::APMFile, arg0::JString)
    return jcall(obj, "setEulerRotSeq", void, (JString,), arg0)
end

function set_i11(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI11", void, (jdouble,), arg0)
end

function set_i12(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI12", void, (jdouble,), arg0)
end

function set_i13(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI13", void, (jdouble,), arg0)
end

function set_i22(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI22", void, (jdouble,), arg0)
end

function set_i23(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI23", void, (jdouble,), arg0)
end

function set_i33(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setI33", void, (jdouble,), arg0)
end

function set_inertia_ref_frame_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setInertiaRefFrameString", void, (JString,), arg0)
end

function set_nutation(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setNutation", void, (jdouble,), arg0)
end

function set_nutation_period(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setNutationPeriod", void, (jdouble,), arg0)
end

function set_nutation_phase(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setNutationPhase", void, (jdouble,), arg0)
end

function set_quaternion(obj::APMFile, arg0::Quaternion)
    return jcall(obj, "setQuaternion", void, (Quaternion,), arg0)
end

function set_quaternion_dot(obj::APMFile, arg0::Quaternion)
    return jcall(obj, "setQuaternionDot", void, (Quaternion,), arg0)
end

function set_quaternion_frame_a_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setQuaternionFrameAString", void, (JString,), arg0)
end

function set_quaternion_frame_b_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setQuaternionFrameBString", void, (JString,), arg0)
end

function set_rate_frame_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setRateFrameString", void, (JString,), arg0)
end

function set_rotation_angles(obj::APMFile, arg0::Vector3D)
    return jcall(obj, "setRotationAngles", void, (Vector3D,), arg0)
end

function set_rotation_rates(obj::APMFile, arg0::Vector3D)
    return jcall(obj, "setRotationRates", void, (Vector3D,), arg0)
end

function set_spacecraft_comment(obj::APMFile, arg0::List)
    return jcall(obj, "setSpacecraftComment", void, (List,), arg0)
end

function set_spin_alpha(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setSpinAlpha", void, (jdouble,), arg0)
end

function set_spin_angle(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setSpinAngle", void, (jdouble,), arg0)
end

function set_spin_angle_vel(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setSpinAngleVel", void, (jdouble,), arg0)
end

function set_spin_comment(obj::APMFile, arg0::List)
    return jcall(obj, "setSpinComment", void, (List,), arg0)
end

function set_spin_delta(obj::APMFile, arg0::jdouble)
    return jcall(obj, "setSpinDelta", void, (jdouble,), arg0)
end

function set_spin_direction(obj::APMFile, arg0::JString)
    return jcall(obj, "setSpinDirection", void, (JString,), arg0)
end

function set_spin_frame_a_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setSpinFrameAString", void, (JString,), arg0)
end

function set_spin_frame_b_string(obj::APMFile, arg0::JString)
    return jcall(obj, "setSpinFrameBString", void, (JString,), arg0)
end


text:
export L2Loss, Regularization, LogLikeLoss, prior_loss, l2lossgradient!,
        colloc_grad

struct Regularization{L,P} <: DiffEqBase.DECostFunction
  λ::L
  penalty::P
end
Regularization(λ) = Regularization{typeof(λ),typeof(L2Penalty())}(λ,L2Penalty())

function (f::Regularization)(p)
  f.λ*value(f.penalty, p)
end

function prior_loss(prior,p)
  ll = 0.0
  if eltype(prior) <: UnivariateDistribution
    for i in 1:length(prior)
      ll -= logpdf(prior[i],p[i])
    end
  else
    ll -= logpdf(prior,p)
  end
  ll
end

struct L2Loss{T,D,U,W,G} <: DiffEqBase.DECostFunction
  t::T
  data::D
  differ_weight::U
  data_weight::W
  colloc_grad::G
  dudt::G
end

function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution)
  data = f.data
  weight = f.data_weight
  diff_weight = f.differ_weight
  colloc_grad = f.colloc_grad
  dudt = f.dudt

  @show sol.u

  if sol isa DiffEqBase.AbstractEnsembleSolution
    failure = any(s.retcode !== :Success && s.retcode !== :Terminated for s in sol)
  else
    failure = sol.retcode !== :Success && sol.retcode !== :Terminated
  end
  failure && return Inf

  sumsq = 0.0

  if weight == nothing
    @inbounds for i in 1:length(sol)
      sumsq +=(data[i] - sol[i])^2
    end
  else
    @inbounds for i in 1:length(sol)
      if typeof(weight) <: Real
        sumsq = sumsq + ((data[i] - sol[i])^2)*weight
      else
        sumsq = sumsq + ((data[i] - sol[i])^2)*weight[i]
      end
    end
  end
  sumsq
end

function (f::L2Loss)(sol::DiffEqBase.DESolution)
  data = f.data
  weight = f.data_weight
  diff_weight = f.differ_weight
  colloc_grad = f.colloc_grad
  dudt = f.dudt

  if sol isa DiffEqBase.AbstractEnsembleSolution
    failure = any(s.retcode !== :Success && s.retcode !== :Terminated for s in sol)
  else
    failure = sol.retcode !== :Success && sol.retcode !== :Terminated
  end
  failure && return Inf

  sumsq = 0.0

  if weight == nothing
    @inbounds for i in 1:length(sol)
      for j in 1:length(sol[i])
        sumsq +=(data[j,i] - sol[j,i])^2
      end
      if diff_weight != nothing && i != 1
          for j in 1:length(sol[i])
            if typeof(diff_weight) <: Real
              sumsq += diff_weight*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
            else
             sumsq += diff_weight[j,i]*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
            end
          end
      end
    end
  else
    @inbounds for i in 1:length(sol)
      if typeof(weight) <: Real
        for j in 1:length(sol[i])
          sumsq = sumsq + ((data[j,i] - sol[j,i])^2)*weight
        end
      else
        for j in 1:length(sol[i])
          sumsq = sumsq + ((data[j,i] - sol[j,i])^2)*weight[j,i]
        end
      end
      if diff_weight != nothing && i != 1
        for j in 1:length(sol[i])
          if typeof(diff_weight) <: Real
            sumsq += diff_weight*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
          else
            sumsq += diff_weight[j,i]*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
          end
        end
      end
    end
  end
  if colloc_grad != nothing
    for i = 1:size(colloc_grad)[2]
      sol.prob.f.f(@view(dudt[:,i]), sol.u[i], sol.prob.p, sol.t[i])
    end
    sumsq += sum(abs2, x - y for (x,y) in zip(dudt, colloc_grad))
  end
  sumsq
end

# Cost functions are written assuming a data matrix
# Turn vectors into a 1xN matrix
matrixize(x) = typeof(x) <: Vector ? reshape(x,1,length(x)) : x

L2Loss(t,data;differ_weight=nothing,data_weight=nothing,colloc_grad=nothing,
  dudt=nothing) = L2Loss(t,matrixize(data),matrixize(differ_weight),
    matrixize(data_weight),matrixize(colloc_grad),
      colloc_grad == nothing ? nothing : zeros(size(colloc_grad)))

function (f::L2Loss)(sol::DiffEqBase.AbstractEnsembleSolution)
  mean(f.(sol.u))
end

#t - 1xN array, data - mxN array, returns mxN array
function colloc_grad(t::T, data::D) where {T, D}
  splines = [Dierckx.Spline1D(t, data[i,:]) for i = 1:size(data)[1]]
  grad = [Dierckx.derivative(spline, t[1:end]) for spline in splines]
  grad = [[grad[1][i],grad[2][i]] for i = 1:length(grad[1])]
  grad = convert(Array, VectorOfArray(grad))
  return grad
end

struct LogLikeLoss{T,D} <: DiffEqBase.DECostFunction
  t::T
  data_distributions::D
  diff_distributions::L where L<:Union{Nothing,D}
  weight
end
LogLikeLoss(t,data_distributions) = LogLikeLoss(t,matrixize(data_distributions),nothing,nothing)
LogLikeLoss(t,data_distributions,diff_distributions) = LogLikeLoss(t,matrixize(data_distributions),matrixize(diff_distributions),1)

function (f::LogLikeLoss)(sol::DESolution)
  distributions = f.data_distributions
  if sol isa DiffEqBase.AbstractEnsembleSolution
    failure = any(s.retcode !== :Success && s.retcode !== :Terminated for s in sol)
  else
    failure = sol.retcode !== :Success && sol.retcode !== :Terminated
  end
  failure && return Inf
  ll = 0.0

  if eltype(distributions) <: UnivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      ll -= logpdf(distributions[i,j],sol[i,j])
    end
  else # MultivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      ll -= logpdf(distributions[i],sol[i])
    end
  end

  if f.diff_distributions != nothing
    distributions = f.diff_distributions
    diff_data = sol.u[2:end] - sol.u[1:end-1]
    fill_length = length(f.t)-length(diff_data)
    for i in 1:fill_length
      push!(diff_data,fill(Inf,size(sol[1])))
    end
    fdll = 0
    if eltype(distributions) <: UnivariateDistribution
      for j in 1:length(f.t)-1, i in 1:length(diff_data[1])
        fdll -= logpdf(distributions[j,i],diff_data[j][i])
      end
    else
      for j in 1:length(f.t)-1
        fdll -= logpdf(distributions[j],diff_data[j])
      end
    end
    ll += f.weight*fdll
  end

  ll
end

function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution)
  distributions = f.data_distributions
  if sol_tmp isa DiffEqBase.AbstractEnsembleSolution
    failure = any(s.retcode !== :Success && s.retcode !== :Terminated for s in sol_tmp)
  else
    failure = sol_tmp.retcode !== :Success && sol_tmp.retcode !== :Terminated
  end
  failure && return Inf
  ll = 0.0
  if eltype(distributions) <: UnivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      vals = [s[i,j] for s in sol]
      ll -= loglikelihood(distributions[i,j],vals)
    end
  else
    for j in 1:length(f.t)
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      vals = [s[i,j] for i in 1:length(sol[1]), s in sol]
      ll -= loglikelihood(distributions[j],vals)
    end
  end

  if f.diff_distributions != nothing
    distributions = f.diff_distributions
    fdll = 0
    if eltype(distributions) <: UnivariateDistribution
      for j in 2:length(f.t), i in 1:length(sol[1][1])
        vals = [s[i,j] - s[i,j-1] for s in sol]
        fdll -= logpdf(distributions[j-1,i],vals)[1]
      end
    else
      for j in 2:length(f.t)
        vals = [s[i,j] - s[i,j-1] for i in 1:length(sol[1]), s in sol]
        fdll -= logpdf(distributions[j-1],vals)[1]
      end
    end
    ll += f.weight*fdll
  end

  ll
end

function l2lossgradient!(grad,sol,data,sensitivities,num_p)
  fill!(grad,0.0)
  data_x_size = size(data,1)
  my_grad = @. 2 * (data - sol)
  u0len = length(data[1])
  K = size(my_grad,2)
  for k in 1:K, i in 1:num_p, j in 1:data_x_size
    grad[i] -= my_grad[j,k]*sensitivities[i][j,k]
  end
end


text:
using MyQuantum
using LinearAlgebra
include("base.jl")

function test_gates()
    display(MyQuantum.is_hermitian.([X, Y, Z, H]))
    display(MyQuantum.is_unitary.([X, Y, Z, H]))
end

function test_conj_trans(s)
    conj_transpose(s) == s'
end

function hermitian_test(s)
    # just reaffirms
    @assert conj.(s') == conj.(s)'
end

function test_outer(a, b)
    # not true
    outer(a, b) != kron(a, b)
end

function test_inner_outer(a, b)
    # the trace of the outer product is equal to the inner product
    tr(outer(a, b)) == dot(a, b)
end

"""
    test_Cgates()

verifies generated C-ops from Cgate() are correct using X, Y, Z, H.
returns true if equal
"""
function test_Cgates()
    CX = Complex[1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]
    CY = Complex[1 0 0 0; 0 1 0 0; 0 0 0 -1im; 0 0 1im 0]
    CZ = Complex[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 -1]
    CH = Complex[1 0 0 0; 0 1 0 0; 0 0 1/sqrt(2) 1/sqrt(2); 0 0 1/sqrt(2) -1/sqrt(2)]
    Cgate.([X, Y, Z, H]) == [CX, CY, CZ, CH]
end

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    struct ColumnTable{T}

This structure helps to access elements that comply with the column access
specification of Tables.jl.

"""
struct ColumnTable{T,V<:AbstractVector{Symbol}}
    table::T
    column_names::V
    size::Tuple{Int,Int}
end

"""
    struct RowTable{T}

This structure helps to access elements that comply with the row access
specification of Tables.jl.

"""
struct RowTable{T,V<:AbstractVector{Symbol}}
    table::T
    column_names::V
    size::Tuple{Int,Int}
end

"""
    PrintInfo{Td,Th,Trn}

This structure stores the information required so that the backends can print
the tables.

"""
@with_kw struct PrintInfo{Td,Th,Trn,Tc <: Tuple, Tf <: Tuple}
    data::Td
    header::Th
    id_cols::Vector{Int}
    id_rows::Vector{Int}
    num_rows::Int
    num_cols::Int
    num_printed_cols::Int
    num_printed_rows::Int
    header_num_rows::Int
    header_num_cols::Int
    show_row_names::Bool
    row_names::Trn
    row_name_alignment::Symbol
    row_name_column_title::String
    alignment::Vector{Symbol}
    cell_alignment::Tc
    formatters::Tf
end


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

"""
    module Modia3D.Graphics

Visual elements used for animation. The visual elements are passed to an
external renderer. Currently, the (free) community edition and the
(commercial) professional editions of the DLR Visualization library are supported
(not all visualization elements of this library are yet interfaced to Modia3D).
Modia3D is designed so that other renderers can be supported as well.

Visualization elements with geometry and visualization material:

![VisuElements](../../resources/images/VisuElementsWith.jpg)

Visualization elements that do not have a visualization material:

![VisuElementsWithout](../../resources/images/VisuElementsWithout.jpg)

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
```
"""
module Graphics

# Colors and material
export rgb, defaultColor, Material, colorPalette

# Visual elements
# (note: Pipe cannot be exported, due to a conflict with Base.Pipe)
export Sphere, Ellipsoid, Box, Cylinder, Capsule, Beam, Cone
export Spring, GearWheel, CoordinateSystem, Grid, FileMesh
export changeDimensionsBox
export Font, TextShape, Screen, XY_Plane, XZ_Plane, YZ_Plane, Left, Right, Center

using  StaticArrays
import Modia3D
import Modia3D.Basics

include("color.jl")
include("material.jl")
include("geometry.jl")
include("text.jl")

end


text:
"""
    example(Model, size)

generates an example graph using model `Model` with `hypergraphsize` of approximately `size`.

# Examples

    example(DCHSBM_sampler, 100)

    example(hyper_pa, 1_000_00)

    example(Kronecker_sampler, 10)

    example(Typing_sampler, 10_00)

    example(ER_sampler, 100_00)
"""
function example end

function example(::Type{DCHSBM_sampler}, size)
    Z = vcat([fill(i, ceil(Integer, ∛size*proportion)) for (i,proportion) in enumerate([.4,.1,.03,.2,.2,.07])]...)
    γ = 1.63
    θ = ((γ+1) .* rand(length(Z))) .^ (-γ-1)
    intensity = inverse_power_intensity_function(8)

    sampler = DCHSBM_sampler(intensity, Z, θ, 4)

    rand(sampler, size ÷ 4)
end

function example(::typeof(hyper_pa), size)
    degree_distribution = AliasTable([1,415,267,152,140,101,93,66,61,47,53,36,44,30,36,32,29,25,31,21,
    16,16,23,21,16,15,18,19,12,20,12,5,12,19,13,11,13,12,6,9,10,11,4,15,5,6,6,7,4,7,7,6,2,6,8,6,5,1,8,
    10,8,3,1,7,3,4,4,3,5,5,4,3,4,4,5,3,2,3,7,5,4,3,2,3,3,2,5,4,1,2,4,4,3,4,1,2,2,2,4,3,1,5,6,1,4,1,4,2,
    1,5,4,1,2,2,2,2,2,2,7,1,1,2,2,2,5,1,1,2,4,3,2,3,2,1,2,1,2,1,2,1,2,2,1,1,3,2,2,1,1,2,2,1,3,3,3,1,2,
    4,1,1,1,1,1,2,1,2,1,1,4,1,1,1,1,1,1,1,1,2,4,1,2,2,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,2,1,2,
    1,2,1,2,1,2,1,1,1,1,2,1,2,1,1,1,1,1,1,1,1,1,1,3,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,2,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1], cumsum([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,1,1,1,1,2,1,1,2,2,2,4,
    1,1,1,1,2,1,1,1,1,1,1,1,2,1,2,1,2,3,1,2,2,1,1,1,1,1,1,4,2,3,1,2,1,2,1,1,1,1,1,1,3,1,1,4,4,1,2,3,1,
    1,1,5,1,1,1,3,1,2,2,2,1,1,1,1,1,2,1,1,2,1,1,1,2,5,1,3,4,6,2,2,3,1,1,9,1,1,4,4,6,1,2,3,4,3,2,11,2,3,
    6,3,5,3,2,1,3,13,19,3,2,2,5,13,5,5,4,4,12,12,7,3,2,4,8,2,1,1,1,2,14,1,21,1,11,4,3,4,4,7,10,18,27,
    14,10,38,10,1,4,2,4,15,9,9,1,1,2,15,20,13,10,46,25,23,53,9,12,97,4,51,53,49,41,7,45,11,7,84,16,13,
    1,26,104,71,31,159]))

    edgesize_distribution = AliasTable([2345, 30991, 41226, 29829, 15690, 8247, 3983, 2279,
    1669, 1158, 927, 1114, 384, 324, 292, 629])

    max_edgesize = length(edgesize_distribution.accept)

    nodes = size ÷ 222 + 18

    hyper_pa(degree_distribution, edgesize_distribution, max_edgesize, nodes)
end

function example(::Type{Kronecker_sampler}, size)
    # parameters copied from https://www.cs.purdue.edu/homes/dgleich/codes/hyperkron/
    initializer = [0.99; 0.2;; 0.2; 0.3;;; 0.2; 0.3;; 0.3; 0.05]
    @static if VERSION < v"1.7"
        initializer = reshape(initializer, (2,2,2))
    end

    edges = size ÷ 3

    sampler = Kronecker_sampler(initializer, max(1, floor(Integer, log(size))), space=edges÷100)

    rand(sampler, edges)
end

function example(::Type{Typing_sampler}, size)

    sampler = Typing_sampler(3, .4, [.1, .2, .3], ['a', 'b', 'c'])

    rand(sampler, size ÷ 3)

end

function example(::Type{ER_sampler}, size)

    s = round(Integer, ∛size)

    sampler = ER_sampler(size, s)

    rand(sampler, s^2)

end

function example(::typeof(er), size)

    s = round(Integer, ∛size)

    er(size, s^2, s)

end

hypergraphsize(graph) = sum(length.(graph))

function MEPS(;generators = [DCHSBM_sampler, Kronecker_sampler, hyper_pa, Typing_sampler, ER_sampler, er], size=100_000, trials=5)
    [begin
        speed = median(begin
            time = @elapsed graph = example(gen, size)
            sum(length.(graph))/time
        end for i in 1:trials)
        gen => round(Integer, speed/10^6)
    end for gen in generators]
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]

############################### starting point
wi = collect(range(0.1, length=2, stop=3))
Ai,Ci,Qi,Ri,Si,m1i,P1i = ACQR_em_starting(y,wi,1,true)

#######################
max_iter=10
tol=1e-10
txo=false
print("Testing ACQR_em_s: ")
A,C,Q,R,m1,P1,loglik,aic = ACQR_em_s(y,Ai,Ci,Qi,Ri,m1i,max_iter,tol,txo)
if sum(A.^2) > 1e-6
    println("OK")
end


text:
using POMDPModelTools
using POMDPs
using POMDPModels
using Random
using Test
using Pkg
using POMDPSimulators
using BeliefUpdaters
using POMDPPolicies
import Distributions.Categorical
using SparseArrays

@testset "POMDPModelTools" begin
    @testset "ordered" begin
        include("test_ordered_spaces.jl")
    end

    # require POMDPModels
    @testset "genbeliefmdp" begin
        include("test_generative_belief_mdp.jl")
    end
    @testset "implement" begin
        include("test_implementations.jl")
    end
    @testset "weightediter" begin
        include("test_weighted_iteration.jl")
    end
    @testset "sparsecat" begin
        include("test_sparse_cat.jl")
    end
    @testset "bool" begin
        include("test_bool.jl")
    end
    @testset "deterministic" begin
        include("test_deterministic.jl")
    end
    @testset "uniform" begin
        include("test_uniform.jl")
    end
    @testset "terminalstate" begin
        include("test_terminal_state.jl")
    end

    # require POMDPModels
    @testset "info" begin
        include("test_info.jl")
    end
    @testset "obsweight" begin
        include("test_obs_weight.jl")
    end

    # require DiscreteValueIteration
    @testset "visolve" begin
        POMDPs.add_registry()
        Pkg.add("DiscreteValueIteration")
        using DiscreteValueIteration
        include("test_fully_observable_pomdp.jl")
        include("test_underlying_mdp.jl")
    end

    @testset "vis" begin
        include("test_visualization.jl")
    end

    @testset "evaluation" begin
        include("test_evaluation.jl")
    end

    @testset "pretty printing" begin
        include("test_pretty_printing.jl")
    end

    @testset "sparse tabular" begin
        include("test_tabular.jl")
    end

end


text:
#############################################################################################

print("   running BazCtx test...")

before_time = time()

@context BazCtx

struct Baz
    x::Int
    y::Float64
    z::String
end

baz_identity(x::Int) = Baz(x, float(x), "$x").x

Cassette.metadatatype(::Type{<:BazCtx}, ::Type{<:Number}) = Float64
x, n = rand(Int), rand()
ctx = enabletagging(BazCtx(), baz_identity)
result = overdub(ctx, baz_identity, tag(x, ctx, n))
@test untag(result, ctx) === x
@test untagtype(typeof(result), typeof(ctx)) === typeof(x)
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === n
@test metameta(result, ctx) === Cassette.NoMetaMeta()
@test hasmetadata(result, ctx)
@test !hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running FooBarCtx test...")

before_time = time()

struct Bar{X,Y,Z}
    x::X
    y::Y
    z::Z
end

mutable struct Foo
    a::Bar{Int}
    b
end

function foo_bar_identity(x)
    bar = Bar(x, x + 1, x + 2)
    foo = Foo(bar, "ha")
    foo.b = bar
    foo.a = Bar(4,5,6)
    foo2 = Foo(foo.a, foo.b)
    foo2.a = foo2.b
    array = Float64[]
    push!(array, foo2.a.x)
    return [array[1]][1]
end

@context FooBarCtx
Cassette.metadatatype(::Type{<:FooBarCtx}, ::Type{T}) where T<:Number = T
x, n = 1, 2
ctx = enabletagging(FooBarCtx(), foo_bar_identity)
result = overdub(ctx, foo_bar_identity, tag(x, ctx, n))
@test untag(result, ctx) === float(x)
@test untagtype(typeof(result), typeof(ctx)) === Float64
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === float(n)
@test metameta(result, ctx) === Cassette.NoMetaMeta()
@test hasmetadata(result, ctx)
@test !hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running TaggedTupleCtx test...")

before_time = time()

@context TaggedTupleCtx
Cassette.metadatatype(::Type{<:TaggedTupleCtx}, ::DataType) = Float64
x = rand()
ctx = enabletagging(TaggedTupleCtx(), 1)
result = overdub(ctx, x -> (x, [x], 1), x)
@test untag(result, ctx) == (x, [x], 1)
@test untagtype(typeof(result), typeof(ctx)) === typeof((x, [x], 1))
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
@test isa(metameta(result, ctx), Tuple{
    Cassette.Immutable{Cassette.Meta{Float64,Cassette.NoMetaMeta}},
    Cassette.Immutable{Cassette.Meta{Float64,Array{Cassette.Meta{Float64,Cassette.NoMetaMeta},1}}},
    Cassette.Immutable{Cassette.Meta{Float64,Cassette.NoMetaMeta}}
})
@test !hasmetadata(result, ctx)
@test hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running ApplyCtx test...")

before_time = time()

@context ApplyCtx
x = rand()
applytest(x) = Core._apply(hypot, (x,), (1, x), 1, x, (1, 2))
ctx = enabletagging(ApplyCtx(), 1)
@test overdub(ctx, applytest, tag(x, ctx)) === applytest(x)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running VATupleCtx test...")

before_time = time()

@context VATupleCtx
x = rand(5)
ctx = enabletagging(VATupleCtx(), 1)
result = overdub(ctx, broadcast, sin, x)
@test untag(result, ctx) == sin.(x)
@test untagtype(typeof(result), typeof(ctx)) === typeof(sin.(x))
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
@test metameta(result, ctx) === Cassette.NoMetaMeta()
@test !hasmetadata(result, ctx)
@test !hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running BroadcastCtx test...")

before_time = time()

@context BroadcastCtx

v, m = rand(5), rand(5)
ctx = enabletagging(BroadcastCtx(), 1)
Cassette.metadatatype(::Type{<:BroadcastCtx}, ::Type{T}) where {T<:Number} = T

result = overdub(ctx, broadcast, (v, m) -> tag(v, ctx, m), v, m)
@test untag(result, ctx) == v
@test untagtype(typeof(result), typeof(ctx)) === typeof(v)
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
@test m == map(metameta(result, ctx)) do x
    @test x.meta === Cassette.NoMetaMeta()
    return x.data
end
@test !hasmetadata(result, ctx)
@test hasmetameta(result, ctx)

result = @overdub(ctx, ((v, m) -> tag(v, ctx, m)).(v, m[1]))
@test untag(result, ctx) == v
@test untagtype(typeof(result), typeof(ctx)) === typeof(v)
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
foreach(metameta(result, ctx)) do x
    @test x.meta === Cassette.NoMetaMeta()
    @test x.data === m[1]
end
@test !hasmetadata(result, ctx)
@test hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running BroadcastCtx2 test...")

before_time = time()

@context BroadcastCtx2

Cassette.metadatatype(::Type{<:BroadcastCtx2{Int}}, ::Type{T}) where {T<:Number} = Vector{T}
Cassette.metadatatype(::Type{<:BroadcastCtx2{Val{N}}}, ::Type{T}) where {N,T<:Number} = NTuple{N,T}

v, m = rand(5), rand(5)
ctx = enabletagging(BroadcastCtx2(metadata=10), 1)
result = overdub(ctx, broadcast, (v, m) -> tag(v, ctx, m), v, [m])
@test untag(result, ctx) == v
@test untagtype(typeof(result), typeof(ctx)) === typeof(v)
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
@test all(e -> e == m, map(metameta(result, ctx)) do x
    @test x.meta === Cassette.NoMetaMeta()
    return x.data
end)
@test !hasmetadata(result, ctx)
@test hasmetameta(result, ctx)

v, m = rand(5), (1.0,2.0,3.0)
ctx = enabletagging(BroadcastCtx2(metadata=Val(3)), 1)
result = overdub(ctx, broadcast, (v, m) -> tag(v, ctx, m), v, [m])
@test untag(result, ctx) == v
@test untagtype(typeof(result), typeof(ctx)) === typeof(v)
@test istagged(result, ctx)
@test istaggedtype(typeof(result), typeof(ctx))
@test metadata(result, ctx) === Cassette.NoMetaData()
foreach(metameta(result, ctx)) do x
    @test x.meta === Cassette.NoMetaMeta()
    @test x.data === m
end
@test !hasmetadata(result, ctx)
@test hasmetameta(result, ctx)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running MetaTypeCtx test...")

before_time = time()

@context MetaTypeCtx
@test Cassette.metatype(typeof(MetaTypeCtx()), DataType) === Cassette.Meta{Cassette.NoMetaData,Cassette.NoMetaMeta}
ctx = enabletagging(MetaTypeCtx(), 1)
@test overdub(ctx, T -> (T,T), Float64) === (Float64, Float64)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running TagConditionalCtx test...")

before_time = time()

@context TagConditionalCtx
ctx = enabletagging(TagConditionalCtx(), 1)
@test overdub(ctx, x -> x ? 1 : 2, tag(true, ctx)) === 1
function condtest(b)
    i = 1
    while i > b
        i -= 1
    end
    if b end
    return i
end
@test overdub(ctx, condtest, tag(false, ctx)) === 0

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running KwargCtx test...")

before_time = time()

@context KwargCtx
ctx = enabletagging(KwargCtx(), 1)
kwargtest(x; y = 1) = x + y
@test overdub(ctx, _y -> kwargtest(3; y = _y), tag(2, ctx)) === 5

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running RecurTagCtx test...")

before_time = time()

@context RecurTagCtx
mutable struct RecurType
    r::RecurType
    RecurType() = new()
end
ctx = enabletagging(RecurTagCtx(), 1)
x = tag(RecurType(), ctx)
@overdub(ctx, x.r = RecurType())
@test istagged(x, ctx)
@test !hasmetadata(x, ctx)
@test hasmetameta(x, ctx)
@test isa(metameta(x, ctx), NamedTuple{(:r,),Tuple{Cassette.Mutable{Cassette.Meta}}})
@test isdefined(untag(x, ctx), :r)
@test !(isdefined(untag(x, ctx).r, :r))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running CrazyPropCtx test...")

before_time = time()

@context CrazyPropCtx

module CrazyPropModule
    const CONST_BINDING = Float64[]

    global GLOBAL_BINDING = 0.0

    struct Foo
        vector::Vector{Float64}
    end

    mutable struct FooContainer
        foo::Foo
    end

    mutable struct PlusFunc
        x::Float64
    end

    (f::PlusFunc)(x) = f.x + x

    const PLUSFUNC = PlusFunc(0.0)

    # implements a very convoluted `sum(x) * sum(y)`
    function crazy_sum_mul(x::Vector{Float64}, y::Vector{Float64})
        @assert length(x) === length(y)
        fooc = FooContainer(Foo(x))
        tmp = y

        # this loop sets:
        # `const_binding == x`
        # `global_binding == prod(y)`
        for i in 1:length(y)
            if iseven(i) # `fooc.foo.vector === y && tmp === x`
                v = fooc.foo.vector[i]
                push!(CONST_BINDING, tmp[i])
                global GLOBAL_BINDING = PLUSFUNC(v)
                PLUSFUNC.x = GLOBAL_BINDING
                fooc.foo = Foo(x)
                tmp = y
            else # `fooc.foo.vector === x && tmp === y`
                v = fooc.foo.vector[i]
                push!(CONST_BINDING, v)
                global GLOBAL_BINDING = PLUSFUNC(tmp[i])
                PLUSFUNC.x = GLOBAL_BINDING
                fooc.foo = Foo(y)
                tmp = x
            end
        end

        # accumulate result
        z = sum(CONST_BINDING) * GLOBAL_BINDING

        # reset global state
        empty!(CONST_BINDING)
        PLUSFUNC.x = 0.0
        global GLOBAL_BINDING = 0.0
        return z
    end
end

x, y = rand(100), rand(100)
primal_result = CrazyPropModule.crazy_sum_mul(x, y)
@test isapprox(primal_result, sum(x) * sum(y))
Cassette.metadatatype(::Type{<:CrazyPropCtx}, ::Type{T}) where T<:Number = T
function Cassette.overdub(ctx::CrazyPropCtx, ::typeof(*), x, y)
    z = untag(x, ctx) * untag(y, ctx)
    if hasmetadata(x, ctx) && hasmetadata(y, ctx)
        return tag(z, ctx, metadata(x, ctx) * metadata(y, ctx))
    elseif hasmetadata(x, ctx)
        return tag(z, ctx, metadata(x, ctx))
    elseif hasmetadata(y, ctx)
        return tag(z, ctx, metadata(y, ctx))
    else
        return z
    end
end
function Cassette.overdub(ctx::CrazyPropCtx, ::typeof(+), x, y)
    z = untag(x, ctx) + untag(y, ctx)
    if hasmetadata(x, ctx) && hasmetadata(y, ctx)
        return tag(z, ctx, metadata(x, ctx) + metadata(y, ctx))
    elseif hasmetadata(x, ctx)
        return tag(z, ctx, metadata(x, ctx))
    elseif hasmetadata(y, ctx)
        return tag(z, ctx, metadata(y, ctx))
    else
        return z
    end
end
ctx = enabletagging(CrazyPropCtx(), CrazyPropModule.crazy_sum_mul)
xm, ym = rand(100), rand(100)
tx = overdub(ctx, broadcast, (v, m) -> tag(v, ctx, m), x, xm)
ty = overdub(ctx, broadcast, (v, m) -> tag(v, ctx, m), y, ym)
tagged_result = overdub(ctx, CrazyPropModule.crazy_sum_mul, tx, ty)
@test isapprox(untag(tagged_result, ctx), primal_result)
@test isapprox(metadata(tagged_result, ctx), CrazyPropModule.crazy_sum_mul(xm, ym))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running DiffCtx test...")

before_time = time()

@context DiffCtx

const DiffCtxWithTag{T} = DiffCtx{Nothing,T}

Cassette.metadatatype(::Type{<:DiffCtx}, ::Type{T}) where {T<:Real} = T

tangent(x, context) = hasmetadata(x, context) ? metadata(x, context) : zero(untag(x, context))

function D(f, x)
    ctx = enabletagging(DiffCtx(), f)
    result = overdub(ctx, f, tag(x, ctx, oftype(x, 1.0)))
    return tangent(result, ctx)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(sin), x::Tagged{T,<:Real}) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    return tag(sin(vx), ctx, cos(vx) * dx)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(cos), x::Tagged{T,<:Real}) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    return tag(cos(vx), ctx, -sin(vx) * dx)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(*), x::Tagged{T,<:Real}, y::Tagged{T,<:Real}) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    vy, dy = untag(y, ctx), tangent(y, ctx)
    return tag(vx * vy, ctx, vy * dx + vx * dy)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(*), x::Tagged{T,<:Real}, y::Real) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    return tag(vx * y, ctx, y * dx)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(*), x::Real, y::Tagged{T,<:Real}) where {T}
    vy, dy = untag(y, ctx), tangent(y, ctx)
    return tag(x * vy, ctx, x * dy)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(+), x::Tagged{T,<:Real}, y::Tagged{T,<:Real}) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    vy, dy = untag(y, ctx), tangent(y, ctx)
    return tag(vx + vy, ctx, dx + dy)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(+), x::Tagged{T,<:Real}, y::Real) where {T}
    vx, dx = untag(x, ctx), tangent(x, ctx)
    return tag(vx + y, ctx, dx)
end

function Cassette.overdub(ctx::DiffCtxWithTag{T}, ::typeof(+), x::Real, y::Tagged{T,<:Real}) where {T}
    vy, dy = untag(y, ctx), tangent(y, ctx)
    return tag(x + vy, ctx, dy)
end

Cassette.overdub(ctx::DiffCtx, ::typeof(sin), x::Real) = sin(x)
Cassette.overdub(ctx::DiffCtx, ::typeof(cos), x::Real) = cos(x)
Cassette.overdub(ctx::DiffCtx, ::typeof(*), x::Real, y::Real) = x * y
Cassette.overdub(ctx::DiffCtx, ::typeof(+), x::Real, y::Real) = x + y
Cassette.overdub(ctx::DiffCtx, ::typeof(*), x, y, z) = Cassette.overdub(ctx, *, Cassette.overdub(ctx, *, x, y), z)
Cassette.overdub(ctx::DiffCtx, ::typeof(+), x, y, z) = Cassette.overdub(ctx, +, Cassette.overdub(ctx, +, x, y), z)

@test D(sin, 1) === cos(1)
@test D(x -> D(sin, x), 1) === -sin(1)
@test D(x -> sin(x) * cos(x), 1) === cos(1)^2 - sin(1)^2
@test D(x -> x * D(y -> x * y, 1), 2) === 4
@test D(x -> x * D(y -> x * y, 2), 1) === 2
@test D(x -> x  * D(y -> 5*x*y, 3), 2) === 20
@test D(x -> x * foo_bar_identity(x), 1) === 2.0

x = rand()
@test D(x -> (x + 2) * (3 + x), x) === 2x + 5
@test D(x -> CrazyPropModule.crazy_sum_mul([x], [x]), x) === (x + x)
@test D(x -> CrazyPropModule.crazy_sum_mul([x, 2], [3, x]), x) === 2x + 5

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running ArrayIndexCtx test...")

before_time = time()

@context ArrayIndexCtx

matrixliteral(x) = Int[x x; x x]

Cassette.metadatatype(::Type{<:ArrayIndexCtx}, ::Type{Int}) = String

ctx = enabletagging(ArrayIndexCtx(), matrixliteral)
result = overdub(ctx, matrixliteral, tag(1, ctx, "hi"))

@test untag(result, ctx) == matrixliteral(1)
@test metameta(result, ctx) == fill(Cassette.Meta("hi", Cassette.NoMetaMeta()), 2, 2)

println("done (took ", time() - before_time, " seconds)")


text:
module Lens

using Compat
import Base: convert, get, getindex

export
  Result,
  Listener,
  Capture,
  lens,
  capture,
  @capture,
  register_benchmarks!,
  disable_benchmarks!,
  disable_benchmarks,

  register!,
  getlisteners,
  clear_all_listeners!,
  enable_all_listeners!,
  enable_listener!,
  disable_listener!,
  disable_all_listeners!,
  nlisteners,

  getindex

include("listener.jl")
include("result.jl")
include("capture.jl")
include("std.jl")

end


text:
include("../examples/solenoid.jl")
using LinearAlgebra
using JLD
using Zygote
function sens(s,nSteps)
    d = 3
    s_ind = 3
    u = rand(d)

    n_runup = 500
    for i = 1:n_runup
        u = next(u,s)
    end

    dui = rand(d,d)
    ppi = rand(d)
    dppi = rand(d,d)

    dpert(x) = jacobian(x -> pert(x,s,s_ind), x)[1]
    d2next(x) = jacobian(x -> dstep(x,s), x)[1]

    vs = zeros(d)
    q = rand(d)
    q /= norm(q)

    pp = rand(d)
    q1 = rand(d)
    vs1 = zeros(d)

    dJds_st = 0.
    dJds_ust = 0.
    nSteps = nSteps + 1
    g = zeros(nSteps)

    # Df represents the derivative of f
    # on the
    # unstable manifold.
    # Xu = a q
    # d = 2
    Xu = zeros(d, nSteps)
    a = zeros(nSteps)
    J = zeros(nSteps)
    Da = zeros(nSteps)
    Dq = zeros(d)
    Dvs = zeros(d)
    Dvs1 = zeros(d)
    # nSteps large number.
	N = 24
    for i = 1:nSteps-1
        J[i] = u[1]*u[1] + u[2]*u[2]
        dJdu = [2*u[1], 2*u[2], 0.]
        ppi .= pert(u,s,s_ind)
        dui .= dstep(u,s) # for large systems, can't store Jacobian.
        dppi .= dpert(u)*inv(dui) # profile against using solve (backslash) here.
        q1 .= dui*q
        z = norm(q1)
        q1 ./= z

        z2 = z*z

        vs1 .= dui*vs + ppi
        a[i+1] = dot(vs1, q1)
        vs1 .-= a[i+1]*q1

        ddui = d2next(u)
        d2q = reshape(ddui*q,d,d)
        Dq = d2q*q/z2 + dui*Dq/z2
        dzdx = dot(z*Dq, q1)
        Dq .= Dq .- dot(Dq,q1)*q1

        Dvs1 .= d2q*vs/z + dui*Dvs/z + dppi*q1
        Da[i+1] = dot(vs1, Dq) + dot(Dvs1, q1)
        Dvs1 .= Dvs1 - Da[i+1]*q1 - a[i+1]*Dq

        g[i+1] = g[i]/z - dzdx/z

        dJds_st += dot(dJdu, vs)/nSteps

        vs .= vs1
        q .= q1
        Dvs .= Dvs1
        u .= next(u,s)

    end
	dJds_ust_K = zeros(N)
	for n = 1:N
        J_shift = J[n+1:end]
        nJ = length(J_shift)
        g_shift = g[2:nJ+1]
        Da_shift = Da[2:nJ+1]
        a_shift = a[2:nJ+1]
		dJds_ust_K[n] = -dot(J_shift,Da_shift .+ g_shift.*a_shift)/nJ
        dJds_ust += dJds_ust_K[n]
    end
    save("../data/sens/solenoid/metrics_$(nSteps)_K.jld", "unstable_contribution_k", dJds_ust_K, "dJds_ust", dJds_ust, "dJds_st", dJds_st)

    return dJds_st + dJds_ust
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
###################################
#####                         #####
#####     KOLMOGOROV FLOW     #####
#####                         #####
###################################

mutable struct KolmogorovFlow{T} <: NS2DProblem{T}

    # --> Unforced dynamics.
    unforced_dyn::UnforcedProblem{T}

    # --> Steady forcing.
    f::Array{Complex{T}, 2}

end

function KolmogorovFlow(p::SimParams, ω::Array{Complex{T}, 2}, n::Int) where T <: Real

    # --> Extract the parameters of the simulation.
    @unpack Lx, Ly, nx, ny, ν = p

    # --> Unforced part of the dynamics.
    unforced_dyn = UnforcedProblem(p, ω)

    # --> Kolmogorov flow forcing.
    x, y = mesh(p)
    x, y = collect(x), collect(y)

    f = -n .* cos.(n .* y * one.(x)')
    f = fft(f)

    # --> Return the container.
    prob = KolmogorovFlow(unforced_dyn, f)

    return prob
end

function rhs!(dΩ::Array{Complex{T},2}, Ω::Array{Complex{T},2}, p::KolmogorovFlow{T}, t::T) where T <: Real

    # --> Unforced dynamics.
    rhs!(dΩ, Ω, p.unforced_dyn, t)

    # --> Add forcing.
    dΩ .+= p.f

    return
end

function simulate(nsprob::KolmogorovFlow, T ; alg=Tsit5(), kwargs...)

    # --> Type of numbers used.
    N = findparams(nsprob)

    # --> Set the ODE Problem.
    tspan = (zero(N), convert(N, T))
    prob = ODEProblem(rhs!, copy(nsprob.unforced_dyn.ω₀), tspan, nsprob)

    # --> Integrate forward in time.
    sol = solve(prob, alg ; kwargs...)

    return sol
end


text:
using Flux
using Statistics
using CalibrationErrors
using ReliabilityDiagrams

# TODO: Work on GPU friendly implementations of these metrics

function get_mean_std_dev(x::AbstractArray; dims = ndims(x))
    μ = mean(x, dims = dims)
    σ = std(x, dims = dims, corrected = false)
    return μ, σ
end

function calculate_entropy(probs, eps = 1e-6; from_logits = true)
    if from_logits
        probs = softmax(probs, dims = 1)
    end
    return -sum(probs .* log.(probs .+ eps), dims = 1)
end

function mean_loglikelihood(preds, labels)
    return mean(log(p[s]) for (s, p) in zip(labels, preds))
end

function brier_score(preds, labels)
    return mean(
        sum(abs2(pi - (i == s)) for (i, pi) in enumerate(p)) for
        (s, p) in zip(labels, preds)
    )
end

function expected_calibration_error(preds, labels, num_bins = 10; from_logits = true)
    if from_logits
        preds = softmax(preds, dims = 1)
    end
    preds = cpu(preds)
    labels = cpu(labels)
    preds = [x for x in eachcol(preds)]
    ece_estimator = ECE(UniformBinning(num_bins), SqEuclidean())
    ece = ece_estimator(preds, labels)
    return ece
end

function prediction_metrics(preds, labels)
    metrics = Dict()
    preds = [x for x in eachcol(preds)]
    # Calulate mean loglikelihood
    metrics["mean_log_likelihood"] = mean_loglikelihood(preds, labels)
    # Brier Score
    metrics["brier_score"] = brier_score(preds, labels)
    return metrics
end


text:
##### Diffusion approximation for semi-infinite media #####

#=
@article{contini1997photon,
  title={Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory},
  author={Contini, Daniele and Martelli, Fabrizio and Zaccanti, Giovanni},
  journal={Applied optics},
  volume={36},
  number={19},
  pages={4587--4599},
  year={1997},
  publisher={Optical Society of America}
}

@article{kienle1997improved,
  title={Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium},
  author={Kienle, Alwin and Patterson, Michael S},
  journal={JOSA A},
  volume={14},
  number={1},
  pages={246--254},
  year={1997},
  publisher={Optical Society of America}
}
=#

### TD Reflectance ###
"""
    refl_DA_semiinf_TD(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

Compute the time-domain reflectance from a semi-infinite medium from Eqn. 36 Contini 97 (m = 0).

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction

# Examples
```jldoctest
julia> refl_DA_semiinf_TD(0:1:5, [0.1,10.0], 1.0, 1.0, 1.0)
6-element Array{Float64,1}:
 0.0
 0.0001440103022493725
 1.446739954231315e-6
 2.7354735244571076e-8
 6.794070985483474e-10
 1.9657536202689858e-11
```
"""
function refl_DA_semiinf_TD(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)
    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
  	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		  A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		  5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		  A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		  A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    Threads.@threads for n in eachindex(t)

        Rt[n] = -(ρ^2/(4D*ν*t[n]))
        Rt[n] = Rt[n] - μa*ν*t[n]
        Rt[n] = -exp(Rt[n])/(2*(4π*D*ν)^(3/2)*sqrt(t[n]*t[n]*t[n]*t[n]*t[n]))
        Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

        if isnan(Rt[n])
            Rt[n] = 0
        end

    end

     return Rt
end

### TD Fluence ###
"""
    fluence_DA_semiinf_TD(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64, z::Float64)

Compute the time-domain fluence in a semi-infinite medium (Eqn. 33 Contini).

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `z::Float64`: the z-depth in medium

# Examples
julia> fluence_DA_semiinf_TD(0:1:5, [0.1,10.0], 1.0, 1.0, 1.0, 1.0)
"""
function fluence_DA_semiinf_TD(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64, z::Float64)
    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
  	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		  A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		  5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		  A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		  A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

    Threads.@threads for n in eachindex(t)

        Rt[n] = -(ρ^2/(4D*ν*t[n]))
        Rt[n] = Rt[n] - μa*ν*t[n]
        Rt[n] = exp(Rt[n])/(2*(4π*D*ν*t[n])^(3/2))
        Rt[n] = Rt[n]*(exp(-((z - zs)^2/(4D*ν*t[n]))) - exp(-((z + 2*ze + zs)^2/(4D*ν*t[n]))))

        if isnan(Rt[n])
            Rt[n] = 0
        end

    end

     return Rt
end

### CW Fluence ###
"""
    fluence_DA_semiinf_CW(ρ::Float64, β::Array{Float64,1}, ndet::Float64, nmed::Float64, z::Float64)

Compute the steady-state fluence in a semi-infinite geometry according to Eqn. 3 of Kienle 1997.

# Arguments
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `z::Float64`: the z-depth within slab (cm)

# Examples
julia> fluence_DA_semiinf_CW(1.0, [0.1,10.0], 1.0,1.0, 0.0)
"""
function fluence_DA_semiinf_CW(ρ::Float64, β::Array{Float64,1}, ndet::Float64, nmed::Float64, z::Float64)
      n = nmed/ndet
      μa = β[1]
      μsp = β[2]
      D = 1/3μsp
      μeff = sqrt(3*μa*μsp)

	  ϕ = 0.0

		if n == 1.0
			A= 1.0
		elseif n > 1.0
			A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
			5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
		else
			A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
		end

		zs = 1/μsp
		ze = 2A*D

        ϕ += exp(-μeff*sqrt(ρ^2 + (z - zs)^2))/(sqrt(ρ^2 + (z - zs)^2))
        ϕ -= exp(-μeff*sqrt(ρ^2 + (z + 2*ze + zs)^2))/(sqrt(ρ^2 + (z + 2*ze + zs)^2))

        if isnan(ϕ)
            ϕ = 0
        end

	return ϕ/(4*π*D)
end

function fluence_DA_semiinf_FD(ρ::Float64, β::Array{Float64,1}, ndet::Float64, nmed::Float64, z::Float64)
    n = nmed/ndet
    ν = 29.9792345/nmed
    μa = β[1] + ω*im/ν
    μsp = β[2]
    D = 1/3μsp
    μeff = sqrt(3*μa*μsp)

    ϕ = 0.0

      if n == 1.0
          A= 1.0
      elseif n > 1.0
          A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
          5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
      else
          A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
      end

      zs = 1/μsp
      ze = 2A*D

      ϕ += exp(-μeff*sqrt(ρ^2 + (z - zs)^2))/(sqrt(ρ^2 + (z - zs)^2))
      ϕ -= exp(-μeff*sqrt(ρ^2 + (z + 2*ze + zs)^2))/(sqrt(ρ^2 + (z + 2*ze + zs)^2))

      if isnan(ϕ)
          ϕ = 0
      end

  return ϕ/(4*π*D)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using Base.Test
using JuLIP, ForwardDiff
using JuLIPMaterials
using StaticArrays

MST = JuLIPMaterials
CB = MST.CauchyBorn
CLE = MST.CLE
FD = ForwardDiff

Id = one(SMatrix{3,3,Float64})

println("------------------------------------------------------------")
println(" Testing Simple Lattice Cauchy--Born Implementation")
println("------------------------------------------------------------")

at = bulk(:Fe)
r0 = rnn(:Fe)
calc = LennardJones(σ = r0) * C2Shift(2.7*r0)
set_calculator!(at, calc)
variablecell!((at))
# minimise!(at)
println("Constructing Wcb . . .")
W = CB.Wcb(at, calc)

println("generate a cauchy Born potential . . . ")
println("check W, grad(W..) evaluate ok . . .")
@test W(Id) == energy(at) / det(defm(at))
@test CB.grad(W, Id) isa AbstractMatrix

println("Finite-difference consistency test")
F = Id + (rand(3,3) - 0.5) * 0.01
W0 = W(F)
dW0 = CB.grad(W, F)[:]
errors = []
for p = 2:10
   h = .1^p
   dWh = zeros(9)
   for i = 1:9
      F[i] += h
      dWh[i] = (W(F) - W0) / h
      F[i] -= h
   end
   push!(errors, norm(dWh-dW0, Inf))
   @printf(" %2d  |  %.3e \n", p, errors[end])
end
passed = minimum(errors) <= 1e-3 * maximum(errors)
if passed
   println("passed")
else
   warn("""It seems the finite-difference test has failed, which indicates
   that there is an inconsistency between the function and gradient
   evaluation. Please double-check this manually / visually. (It is
   also possible that the function being tested is poorly scaled.)""")
end

# ------ some more random tests left-over from debugging -------

# setup a Cauchy-Born model for bulk Fe
atu = bulk(:Fe)
r0 = rnn(:Fe)
calc = LennardJones(σ = r0) * C2Shift(2.7*r0)
set_calculator!(atu, calc)
variablecell!(atu)
println("Constructing Wcb . . .")
W = CB.Wcb(atu, calc, normalise = :atoms)

energy(atu)
W(Id)

@test CB.grad(W, Id) ≈ -virial(atu)
@test stress(atu) == CB.grad(W, Id) / det(cell(atu))

at1 = deepcopy(atu)
F0 = defm(at1)
for n = 1:5
   F = Id + 0.1 * rand(3,3)
   set_defm!(at1, F * F0)
   @test energy(calc, at1) == W(F)
end


text:
using Test
using FunctionalStateMachine

include("testStateMachine.jl")


text:
function velocity_index(mechanism::Mechanism{T,Nn,Ne}) where {T,Nn,Ne}
    ind = []
    off = 0
	for id in mechanism.root_to_leaves
        (id > Ne) && continue # only treat joints
        joint = mechanism.joints[id]
        nu = input_dimension(joint)
        push!(ind, Vector(off + nu .+ (1:nu)))
        off += 2nu
    end
    return vcat(ind...)
end

# find all the joints parents of a body
function parent_joints(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, body::Body) where {T,Nn,Ne,Nb,Ni}
	ids = parents(mechanism.system, body.id)
	ids = intersect(ids, 1:Ne) # filter out the bodies
	return [get_node(mechanism, id) for id in ids]
end

# dimensions
"""
    minimal_dimension(mechanism)

    dimension of a mechanism's minimal representation

    mechanism: Mechanism
"""
function minimal_dimension(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}) where {T,Nn,Ne,Nb,Ni}
    nx = 0
    free_rot_base = false # we are going to check if the link attached to the base has free orientation
    nx = 2 * input_dimension(mechanism,
        ignore_floating_base=false)
    free_rot_base && (nx += 1)
    return nx
end

"""
    maximal_dimension(mechanism)

    dimension of a mechanism's maximal representation

    mechanism: Mechanism
"""
maximal_dimension(mechanism::Mechanism{T,Nn,Ne,Nb}; attjac::Bool=false) where {T,Nn,Ne,Nb} = attjac ? 12Nb : 13Nb

"""
    input_dimension(mechanism)

    return the number of inputs for mechanism

    mechanism: Mechanism
"""
function input_dimension(mechanism::Mechanism{T,Nn,Ne,Nb,Ni};
    ignore_floating_base::Bool=false) where {T,Nn,Ne,Nb,Ni}
    nu = 0
    for joint in mechanism.joints
        nu += input_dimension(joint, ignore_floating_base=ignore_floating_base)
    end
    return nu
end

"""
    set_floating_base(mechanism, name)

    returns a mechanism with modified joints having body identified with 'name' as the floating base

    mechanism: Mechanism
    name: Symbol, identifier for floating-base body
"""
function set_floating_base(mechanism::Mechanism, name::Symbol)
    # copy
    mechanism = deepcopy(mechanism)

    # original floating base
    body1 = mechanism.origin

    # new floating base
    body2 = get_body(mechanism, name)

    # find joint for new floating base
    body = body2
    joint = mechanism.joints[findall(j -> j.child_id == body.id, mechanism.joints)[1]]
    joint_list = JointConstraint[]
    push!(joint_list, joint)

    while joint_list[end].parent_id != body1.id
        @show length(joint_list)
        push!(joint_list, mechanism.joints[findall(j -> j.child_id == joint_list[end].parent_id, mechanism.joints)[1]])
    end

    @show [j.name for j in joint_list]
    @show [(get_body(mechanism, j.parent_id).name, get_body(mechanism, j.child_id).name) for j in joint_list]

    # reverse JointConstraint (p_old -> c_old) to (c_old -> p_old)
    new_joint_list = JointConstraint[]
    for j in joint_list
        j = deepcopy(j)
        if j.parent_id == body1.id
            j.child_id = body2.id
            push!(new_joint_list, j)
            continue
        else
            parent_id = j.parent_id
            child_id = j.child_id
            pbody = get_body(mechanism, parent_id)
            cbody = get_body(mechanism, child_id)

            # translational conversion
            t = j.translational
            t_reverse = Translational{typeof(t).parameters[1], typeof(t).parameters[2]}(cbody, pbody;
                    parent_vertex=t.vertices[2],
                    child_vertex=t.vertices[1],
                    axis=-t.axis,
                    spring=t.spring,
                    damper=t.damper,
                    spring_offset=t.spring_offset,
                    joint_limits=t.joint_limits,
                    spring_type=t.spring_type)[1]

            r = j.rotational
            r_reverse = Rotational{typeof(r).parameters[1], typeof(r).parameters[2]}(cbody, pbody;
                    axis=-r.axis,
                    axis_offset=inv(r.axis_offset),
                    spring=r.spring,
                    damper=r.damper,
                    spring_offset=r.spring_offset,
                    joint_limits=r.joint_limits,
                    spring_type=r.spring_type)[1]

            j.parent_id = child_id
            j.child_id = parent_id
            j.translational = t_reverse
            j.rotational = r_reverse
        end
        push!(new_joint_list, j)
    end

    # set joints
    for j in new_joint_list
        mechanism.joints[j.id] = j
    end

    return Mechanism(mechanism.origin, mechanism.bodies, mechanism.joints, mechanism.contacts,
        gravity=mechanism.gravity,
        timestep=mechanism.timestep,
        # constructor will skip setting springs and dampers, using previously defined mech.joints
        spring=-1.0,
        damper=-1.0)
end


text:
using Test
using Distributions

import GlobalSensitivityAnalysis: scale_samples!

# scale_samples!
seq = rand(100, 8)
original_seq = copy(seq)
dists = [Normal(1, 0.2), Uniform(0.75, 1.25), LogNormal(0, 0.5), TriangularDist(0, 4, 1)]
scale_samples!(seq, dists)
@test size(seq) == size(original_seq)
@test seq != original_seq

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
"""
    AbstractProgressiveHedgingAttribute

Abstract supertype for attribute objects specific to the progressive-hedging algorithm.
"""
abstract type AbstractProgressiveHedgingAttribute <: AbstractStructuredOptimizerAttribute end
"""
    PrimalTolerance

An optimizer attribute for specifying the primal tolerance in the progressive-hedging algorithm.
"""
struct PrimalTolerance <: AbstractProgressiveHedgingAttribute end
"""
    DualTolerance

An optimizer attribute for specifying the dual tolerance in the progressive-hedging algorithm.
"""
struct DualTolerance <: AbstractProgressiveHedgingAttribute end
"""
    Regularizer

An optimizer attribute for specifying a penalization procedure to be used in the progressive-hedging algorithm. Options are:

- [`Fixed`](@ref):  Fixed penalty (default) ?Fixed for parameter descriptions.
- [`Adaptive`](@ref): Adaptive penalty update ?Adaptive for parameter descriptions.
"""
struct Penalizer <: AbstractProgressiveHedgingAttribute end
"""
    Penaltyterm

An optimizer attribute for specifying what proximal term should be used in subproblemsof the progressive-hedging algorithm. Options are:

- [`Quadratic`](@ref) (default)
- [`Linearized`](@ref)
- [`InfNorm`](@ref)
- [`ManhattanNorm`](@ref)
"""
struct Penaltyterm <: AbstractProgressiveHedgingAttribute end
"""
    PenalizationParameter

Abstract supertype for penalization-specific attributes.
"""
abstract type PenalizationParameter <: AbstractProgressiveHedgingAttribute end