text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DistributionTransform{
    VF <: VariateForm,
    ST <: VariateSpace,
    SF <: VariateSpace,
    DT <: ContinuousDistribution,
    DF <: Distribution{VF,Continuous}
} <: VariateTransform{VF,ST,SF}
    target_dist::DT
    source_dist::DF
    target_space::ST
    source_space::SF
end

# ToDo: Add field to cache dist-specific pre-calculated values? May be useful
# for truncated dists and others.

function _distrafo_ctor_impl(target_dist::Distribution, source_dist::Distribution)
    @argcheck totalndof(varshape(target_dist)) == totalndof(varshape(source_dist))
    target_space = getspace(target_dist)
    source_space = getspace(source_dist)
    DistributionTransform(target_dist, source_dist, target_space, source_space)
end

DistributionTransform(target_dist::Distribution{VF,Continuous}, source_dist::Distribution{VF,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{Multivariate,Continuous}, source_dist::Distribution{VF,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{VF,Continuous}, source_dist::Distribution{Multivariate,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{Multivariate,Continuous}, source_dist::Distribution{Multivariate,Continuous}) =
    _distrafo_ctor_impl(target_dist, source_dist)

# apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)
function apply_dist_trafo end

target_space(trafo::DistributionTransform) = trafo.target_space
source_space(trafo::DistributionTransform) = trafo.source_space

import Base.∘
function ∘(a::DistributionTransform, b::DistributionTransform)
    @argcheck a.source_dist == b.target_dist
    DistributionTransform(a.target_dist, b.source_dist)
end

Base.inv(trafo::DistributionTransform) = DistributionTransform(trafo.source_dist, trafo.target_dist)

ValueShapes.varshape(trafo::DistributionTransform) = varshape(trafo.target_dist)

function apply_vartrafo_impl(trafo::DistributionTransform, v::Any, prev_ladj::Real)
    apply_dist_trafo(trafo.target_dist, trafo.source_dist, v, prev_ladj)
end

function apply_vartrafo_impl(trafo::InvVT{<:DistributionTransform}, v::Any, prev_ladj::Real)
    apply_vartrafo_impl(inv(trafo.orig), v, prev_ladj)
end

const _StdDistType = Union{Uniform, Normal}

_trg_disttype(::Type{Uniform}, ::Type{Univariate}) = StandardUvUniform
_trg_disttype(::Type{Uniform}, ::Type{Multivariate}) = StandardMvUniform
_trg_disttype(::Type{Normal}, ::Type{Univariate}) = StandardUvNormal
_trg_disttype(::Type{Normal}, ::Type{Multivariate}) = StandardMvNormal

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::Distribution{Univariate,Continuous})
    trg_dt = _trg_disttype(disttype, Univariate)
    trg_dt()
end

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::Distribution{Multivariate,Continuous})
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_dt(length(source_dist))
end

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::ContinuousDistribution)
    trg_vs = varshape(source_dist)
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_dt(totalndof(trg_vs))
end

function DistributionTransform(disttype::Type{<:_StdDistType}, source_dist::Distribution{Univariate,Continuous})
    trg_dt = _trg_disttype(disttype, Univariate)
    DistributionTransform(trg_dt(), source_dist)
end

function DistributionTransform(disttype::Type{<:_StdDistType}, source_dist::Distribution{Multivariate,Continuous})
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_d = trg_dt(length(source_dist))
    DistributionTransform(trg_d, source_dist)
end

function DistributionTransform(disttype::Type{<:_StdDistType}, source_dist::ContinuousDistribution)
    trg_vs = varshape(source_dist)
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_d = trg_dt(totalndof(trg_vs))
    DistributionTransform(trg_d, source_dist)
end

function apply_dist_trafo(trg_d::Distribution{Multivariate}, src_d::ReshapedDist, src_v::Any, prev_ladj::Real)
    src_vs = varshape(src_d)
    @argcheck length(trg_d) == totalndof(src_vs)
    apply_dist_trafo(trg_d, unshaped(src_d), unshaped(src_v, src_vs), prev_ladj)
end

function apply_dist_trafo(trg_d::ReshapedDist, src_d::Distribution{Multivariate}, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    trg_vs = varshape(trg_d)
    @argcheck totalndof(trg_vs) == length(src_d)
    r = apply_dist_trafo(unshaped(trg_d), src_d, src_v, prev_ladj)
    (v = trg_vs(r.v), ladj = r.ladj)
end

function apply_dist_trafo(trg_d::ReshapedDist, src_d::ReshapedDist, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    trg_vs = varshape(trg_d)
    src_vs = varshape(src_d)
    @argcheck totalndof(trg_vs) == totalndof(src_vs)
    r = apply_dist_trafo(unshaped(trg_d), unshaped(src_d), unshaped(src_v, src_vs), prev_ladj)
    (v = trg_vs(r.v), ladj = r.ladj)
end

function apply_dist_trafo(trg_d::Distribution{VF,VS}, src_d::Distribution{VF,VS}, src_v::Any, prev_ladj::Real) where {VF,VS}
    uniform_dist = _trg_dist(Uniform, src_d)
    @assert uniform_dist == _trg_dist(Uniform, trg_d)
    uniform_v, uniform_ladj = apply_dist_trafo(uniform_dist, src_d, src_v, prev_ladj)
    apply_dist_trafo(trg_d, uniform_dist, uniform_v, uniform_ladj)
end

function apply_dist_trafo(trg_d::DT, src_d::DT, src_v::Real, prev_ladj::Real) where {DT <: StdUvDist}
    (v = src_v, ladj = prev_ladj)
end

function apply_dist_trafo(trg_d::DT, src_d::DT, src_v::AbstractVector{<:Real}, prev_ladj::Real) where {DT <: StdMvDist}
    @argcheck length(trg_d) == length(src_d) == length(eachindex(src_v))
    (v = src_v, ladj = prev_ladj)
end

function apply_dist_trafo(trg_d::Distribution{Univariate}, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    @argcheck length(src_d) == length(eachindex(src_v)) == 1
    apply_dist_trafo(trg_d, view(src_d, 1), first(src_v), prev_ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::Distribution{Univariate}, src_v::Real, prev_ladj::Real)
    @argcheck length(trg_d) == 1
    r = apply_dist_trafo(view(trg_d, 1), src_d, first(src_v), prev_ladj)
    (v = unshaped(r.v), ladj = r.ladj)
end

@inline function _value_and_ladj(::typeof(cdf), d::Distribution{Univariate,Continuous}, x::Real)
    u = cdf(d, x)
    ladj = + logpdf(d, x)
    (u, ladj)
end

@inline function _value_and_ladj(::typeof(quantile), d::Distribution{Univariate,Continuous}, u::Real)
    x = quantile(d, u)
    ladj = - logpdf(d, x)
    (x, ladj)
end

function _eval_distr_trafo_error(f::Function, d::Distribution, x::Any)
    ErrorException("Evaluating transformation via $f for distribution of type $(typeof(d)) failed at $x")
end

function _eval_dist_trafo_func(f::Function, d::Distribution{Univariate,Continuous}, src_v::Real, prev_ladj::Real)
    try
        if isnan(prev_ladj)
            trg_v = f(d, src_v)
            var_trafo_result(trg_v, src_v)
        else
            trg_v, trafo_ladj = _value_and_ladj(f, d, src_v)
            var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
        end
    catch err
        if err isa DomainError
            var_trafo_result(NaN, src_v)
        else
            rethrow(_eval_distr_trafo_error(f, d, src_v))
        end
    end
end

function apply_dist_trafo(::StandardUvUniform, src_d::Distribution{Univariate,Continuous}, src_v::Real, prev_ladj::Real)
    _eval_dist_trafo_func(cdf, src_d, src_v, prev_ladj)
end

function apply_dist_trafo(trg_d::Distribution{Univariate,Continuous}, ::StandardUvUniform, src_v::Real, prev_ladj::Real)
    TV = float(typeof(src_v))
    # Avoid src_v ≈ 0 and src_v ≈ 1 to avoid infinite variate values for target distributions with infinite support:
    mod_src_v = ifelse(src_v == 0, zero(TV) + eps(TV), ifelse(src_v == 1, one(TV) - eps(TV), convert(TV, src_v)))
    trg_v, ladj = _eval_dist_trafo_func(quantile, trg_d, mod_src_v, prev_ladj)
    (v = trg_v, ladj = ladj)
end

function _dist_trafo_rescale_impl(trg_d, src_d, src_v::Real, prev_ladj::Real)
    R = float(typeof(src_v))
    trg_offs, trg_scale = location(trg_d), scale(trg_d)
    src_offs, src_scale = location(src_d), scale(src_d)
    rescale_factor = trg_scale / src_scale
    trg_v = (src_v - src_offs) * rescale_factor + trg_offs
    trafo_ladj = log(rescale_factor)
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

apply_dist_trafo(trg_d::Uniform, src_d::Uniform, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
apply_dist_trafo(trg_d::StandardUvUniform, src_d::Uniform, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
apply_dist_trafo(trg_d::Uniform, src_d::StandardUvUniform, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
apply_dist_trafo(trg_d::Normal, src_d::Normal, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
apply_dist_trafo(trg_d::StandardUvNormal, src_d::Normal, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
apply_dist_trafo(trg_d::Normal, src_d::StandardUvNormal, src_v::Real, prev_ladj::Real) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)

# ToDo: Optimized implementation for Distributions.Truncated <-> StandardUvUniform

function apply_dist_trafo(trg_d::StandardMvUniform, src_d::StandardMvNormal, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    @argcheck length(trg_d) == length(src_d)
    apply_dist_trafo(convert(Product, trg_d), convert(Product, src_d), src_v, prev_ladj)
end

function apply_dist_trafo(trg_d::StandardMvNormal, src_d::StandardMvUniform, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    @argcheck length(trg_d) == length(src_d)
    apply_dist_trafo(convert(Product, trg_d), convert(Product, src_d), src_v, prev_ladj)
end

function apply_dist_trafo(trg_d::Distributions.Product, src_d::Distributions.Product, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    rs = apply_dist_trafo.(trg_d.v, src_d.v, src_v, zero(Float32))
    trg_v = broadcast(r -> r.v, rs)
    trafo_ladj = sum(map(r -> r.ladj, rs))
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::Distributions.Product, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    apply_dist_trafo(convert(Distributions.Product, trg_d), src_d, src_v, prev_ladj)
end

function apply_dist_trafo(trg_d::Distributions.Product, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    apply_dist_trafo(trg_d, convert(Distributions.Product, src_d), src_v, prev_ladj)
end

function _ntdistelem_to_stdmv(trg_d::StdMvDist, sd::Distribution, src_v_unshaped::AbstractVector{<:Real}, acc::ValueAccessor)
    td = view(trg_d, ValueShapes.view_range(Base.OneTo(length(trg_d)), acc))
    sv = stripscalar(view(src_v_unshaped, acc))
    apply_dist_trafo(td, sd, sv, 0)
end

function _ntdistelem_to_stdmv(trg_d::StdMvDist, sd::ConstValueDist, src_v_unshaped::AbstractVector{<:Real}, acc::ValueAccessor)
    (v = Bool[], ladj = zero(Float32))
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::NamedTupleDist, src_v::Union{NamedTuple,ShapedAsNT}, prev_ladj::Real)
    src_vs = varshape(src_d)
    @argcheck length(trg_d) == totalndof(src_vs)
    src_v_unshaped = unshaped(src_v, src_vs)
    rs = map((acc, sd) -> _ntdistelem_to_stdmv(trg_d, sd, src_v_unshaped, acc), values(src_vs), values(src_d))
    trg_v = vcat(map(r -> r.v, rs)...)
    trafo_ladj = sum(map(r -> r.ladj, rs))
    var_trafo_result(trg_v, src_v_unshaped, trafo_ladj, prev_ladj)
end

function _stdmv_to_ntdistelem(td::Distribution, src_d::StdMvDist, src_v::AbstractVector{<:Real}, acc::ValueAccessor)
    sd = view(src_d, ValueShapes.view_range(Base.OneTo(length(src_d)), acc))
    sv = view(src_v, ValueShapes.view_range(axes(src_v, 1), acc))
    apply_dist_trafo(td, sd, sv, 0)
end

function _stdmv_to_ntdistelem(td::ConstValueDist, src_d::StdMvDist, src_v::AbstractVector{<:Real}, acc::ValueAccessor)
    (v = Bool[], ladj = zero(Float32))
end

function apply_dist_trafo(trg_d::NamedTupleDist, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::Real)
    trg_vs = varshape(trg_d)
    @argcheck totalndof(trg_vs) == length(src_d)
    rs = map((acc, td) -> _stdmv_to_ntdistelem(td, src_d, src_v, acc), values(trg_vs), values(trg_d))
    trg_v_unshaped = vcat(map(r -> unshaped(r.v), rs)...)
    trafo_ladj = sum(map(r -> r.ladj, rs))
    unshaped_result = var_trafo_result(trg_v_unshaped, src_v, trafo_ladj, prev_ladj)
    (v = trg_vs(unshaped_result.v), ladj = unshaped_result.ladj)
end

#=

# Optimized transformations
# ToDo: Test and compare performance with generic version.

function apply_dist_trafo(::StandardUvUniform, src_d::Exponential, src_v::Real, prev_ladj::Real)
    R = typeof(src_v)
    theta = scale(src_d)
    trg_v = exponential_cdf(src_v, theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = exponential_cdf_ladj(src_v, theta))
        var_trafo_result(trg_v, src_v, convert(R, trafo_ladj, prev_ladj)
    end
end

function apply_dist_trafo(trg_d::Exponential, ::StandardUvUniform, src_v::Real, prev_ladj::Real)
    R = typeof(src_v)
    theta = scale(trg_d)
    trg_v = exponential_invcdf(src_v, theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - exponential_cdf_ladj(trg_v, theta))
        var_trafo_result(trg_v, src_v, convert(R, trafo_ladj, prev_ladj)
    end
end

function apply_dist_trafo(::StandardUvUniform, src_d::Logistic, src_v::Real, prev_ladj::Real)
    R = typeof(src_v)
    mu, theta = location(src_d), scale(src_d)
    trg_v = logistic_cdf(src_v, mu, theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - logistic_invcdf_ladj(trg_v, mu, theta))
        var_trafo_result(trg_v, src_v, convert(R, trafo_ladj, prev_ladj)
    end
end

function apply_dist_trafo(trg_d::Logistic, ::StandardUvUniform, src_v::Real, prev_ladj::Real)
    R = typeof(src_v)
    mu, theta = location(trg_d), scale(trg_d)
    trg_v = logistic_invcdf(src_v, mu, theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = logistic_invcdf_ladj(src_v, mu, theta))
        var_trafo_result(trg_v, src_v, convert(R, trafo_ladj, prev_ladj)
    end
end

# cdf(Normal(mu, sigma), x) may be just as fast:
normal_cdf(x::Real, mu::Real, sigma::Real) = erfc((mu - x) / sigma * invsqrt2) / 2

# quantile(Normal(mu, sigma), x) may be just as fast:
normal_invcdf(u::Real, mu::Real, sigma::Real) = mu - erfcinv(2 * u) * sigma / invsqrt2

# lefttruncexp_cdf(x::Real, x0::Real, theta::Real) = exp((x - x0) / -theta)
# lefttruncexp_invcdf(u::Real, x0::Real, theta::Real) = -theta * log(u) + x0

weibull_cdf(x::Real, alpha::Real, theta::Real) = 1 - exp(- (x / theta)^alpha)
weibull_invcdf(x::Real, alpha::Real, theta::Real) = (-log(1 - x))^(1/alpha) * theta

=#


text:
@mlj_model mutable struct S1 <: Supervised
end

@mlj_model mutable struct U1 <: Unsupervised
end

@mlj_model mutable struct D1 <: Deterministic
    a::Int = 1::(_ > 0)
end

@mlj_model mutable struct P1 <: Probabilistic
    a::Int = 1::(_ > 0)
end

@mlj_model mutable struct I1 <: Interval
end

@mlj_model mutable struct SA <: SupervisedAnnotator
end

@mlj_model mutable struct UA <: UnsupervisedAnnotator
end

foo(::P1) = 0
bar(::P1) = nothing

@testset "traits" begin
    ms = S1()
    mu = U1()
    md = D1()
    mp = P1()
    mi = I1()
    sa = SA()
    ua = UA()

    @test input_scitype(ms)  == Unknown
    @test output_scitype(ms) == Unknown
    @test target_scitype(ms) == Unknown
    @test is_pure_julia(ms)  == false

    @test package_name(ms)    == "unknown"
    @test package_license(ms) == "unknown"
    @test load_path(ms)       == "unknown"
    @test package_uuid(ms)    == "unknown"
    @test package_url(ms)     == "unknown"

    @test is_wrapper(ms)       == false
    @test supports_online(ms)  == false
    @test supports_weights(ms) == false
    @test iteration_parameter(ms) === nothing

    @test hyperparameter_ranges(md) == (nothing,)

    @test docstring(ms) == "S1 from unknown.jl.\n[Documentation](unknown)."
    @test name(ms)      == "S1"

    @test is_supervised(ms)
    @test is_supervised(sa)
    @test !is_supervised(mu)
    @test !is_supervised(ua)
    @test prediction_type(ms) == :unknown
    @test prediction_type(md) == :deterministic
    @test prediction_type(mp) == :probabilistic
    @test prediction_type(mi) == :interval

    @test hyperparameters(md) == (:a,)
    @test hyperparameter_types(md) == ("Int64",)

    # implemented methods is deferred
    setlight()
    @test_throws M.InterfaceError implemented_methods(mp)
    setfull()
    M.implemented_methods(::FI, M::Type{<:MLJType}) =
        getfield.(methodswith(M), :name)
    @test Set(implemented_methods(mp)) == Set([:clean!,:bar,:foo])
end

module Fruit

import MLJModelInterface.MLJType

struct Banana <: MLJType end

end

import .Fruit

@testset "extras" begin
    @test docstring(Float64) == "Float64"
    @test docstring(Fruit.Banana) == "Banana"
end

@testset "`_density` - helper for predict_scitype fallback" begin
    for T in [Continuous, Count, Textual]
        @test M._density(AbstractArray{T,3}) ==
            AbstractArray{Density{T},3}
    end

    for T in [Finite,
              Multiclass,
              OrderedFactor,
              Infinite,
              Continuous,
              Count,
              Textual]
        @test M._density(AbstractVector{<:T}) ==
            AbstractVector{Density{<:T}}
        @test M._density(Table(T)) == Table(Density{T})
    end

    for T in [Finite, Multiclass, OrderedFactor]
        @test M._density(AbstractArray{<:T{2},3}) ==
            AbstractArray{Density{<:T{2}},3}
        @test M._density(AbstractArray{T{2},3}) ==
            AbstractArray{Density{T{2}},3}
        @test M._density(Table(T{2})) == Table(Density{T{2}})
    end
end

@mlj_model mutable struct P2 <: Probabilistic end
M.target_scitype(::Type{<:P2}) = AbstractVector{<:Multiclass}
M.input_scitype(::Type{<:P2}) = Table(Continuous)

@mlj_model mutable struct U2 <: Unsupervised end
M.output_scitype(::Type{<:U2}) = AbstractVector{<:Multiclass}
M.input_scitype(::Type{<:U2}) = Table(Continuous)

@mlj_model mutable struct S2 <: Static end
M.output_scitype(::Type{<:S2}) = AbstractVector{<:Multiclass}
M.input_scitype(::Type{<:S2}) = Table(Continuous)

@testset "operation scitypes" begin
    @test predict_scitype(P2()) == AbstractVector{Density{<:Multiclass}}
    @test transform_scitype(P2()) == Unknown
    @test transform_scitype(U2()) == AbstractVector{<:Multiclass}
    @test inverse_transform_scitype(U2()) == Table(Continuous)
    @test predict_scitype(U2()) == Unknown
    @test transform_scitype(S2()) == AbstractVector{<:Multiclass}
    @test inverse_transform_scitype(S2()) == Table(Continuous)
end

@testset "abstract_type, fit_data_scitype" begin
    @test abstract_type(P2()) == Probabilistic
    @test abstract_type(S1()) == Supervised
    @test abstract_type(U1()) == Unsupervised
    @test abstract_type(D1()) == Deterministic
    @test abstract_type(P1()) == Probabilistic

    @test fit_data_scitype(P2()) ==
        Tuple{Table(Continuous),AbstractVector{<:Multiclass}}
    @test fit_data_scitype(U2()) == Tuple{Table(Continuous)}
    @test fit_data_scitype(S2()) == Tuple{}
end

true

ext: jl


text:
## Load the packages.
using LoopVectorization
using Interpolations
using OffsetArrays
using Statistics
using BenchmarkTools
using PyPlot

## Create the refine functions.
function refine_field_int!(hi, hi_tmp, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for k in 0:size(lo, 3)-3, j in 0:size(lo, 2)-3, i in 0:size(lo, 1)-3
            i_lo = i + 2
            j_lo = j + 2
            k_lo = k + 2
            i_hi = i*2 + 2
            j_hi = j*2 + 2
            k_hi = k*2 + 2
            for kk in 0:1, jj in 0:1
                hi_tmp[i_hi, j_hi+jj, k_hi+kk] = lo[i_lo, j_lo, k_lo]
                hi_tmp[i_hi+1, j_hi+jj, k_hi+kk] = 1//2*(lo[i_lo, j_lo, k_lo] + lo[i_lo+1, j_lo, k_lo])
            end
        end

        @views begin
            @tturbo a_hi_tmp[1, :, :] .= a_hi_tmp[end-1, :, :]
            @tturbo a_hi_tmp[end, :, :] .= a_hi_tmp[2, :, :]
            @tturbo a_hi_tmp[:, 1, :] .= a_hi_tmp[:, end-1, :]
            @tturbo a_hi_tmp[:, end, :] .= a_hi_tmp[:, 2, :]
            @tturbo a_hi_tmp[:, :, 1] .= a_hi_tmp[:, :, end-1]
            @tturbo a_hi_tmp[:, :, end] .= a_hi_tmp[:, :, 2]
        end

        coef = OffsetArray(zeros(3, 3), -1:1, -1:1)
        for kk in -1:1, jj in -1:1
            coef[jj, kk] = (1//2)^(2 + abs(jj) + abs(kk))
        end

        @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
            hi[i, j, k] = 0
            for kk in -1:1, jj in -1:1
                hi[i, j, k] += coef[jj, kk] * hi_tmp[i, j+jj, k+kk]
            end
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for k in 0:size(lo, 3)-3, j in 0:size(lo, 2)-3, i in 0:size(lo, 1)-3
            i_lo = i + 2
            j_lo = j + 2
            k_lo = k + 2
            i_hi = i*3 + 2
            j_hi = j*3 + 2
            k_hi = k*3 + 2
            for kk in 0:2, jj in 0:2
                hi_tmp[i_hi  , j_hi+jj, k_hi+kk] = lo[i_lo  , j_lo, k_lo]
                hi_tmp[i_hi+1, j_hi+jj, k_hi+kk] = lo[i_lo  , j_lo, k_lo]
                hi_tmp[i_hi+2, j_hi+jj, k_hi+kk] = lo[i_lo+1, j_lo, k_lo]
            end
        end

        @views begin
            @tturbo a_hi_tmp[1, :, :] .= a_hi_tmp[end-1, :, :]
            @tturbo a_hi_tmp[end, :, :] .= a_hi_tmp[2, :, :]
            @tturbo a_hi_tmp[:, 1, :] .= a_hi_tmp[:, end-1, :]
            @tturbo a_hi_tmp[:, end, :] .= a_hi_tmp[:, 2, :]
            @tturbo a_hi_tmp[:, :, 1] .= a_hi_tmp[:, :, end-1]
            @tturbo a_hi_tmp[:, :, end] .= a_hi_tmp[:, :, 2]
        end

        @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
            hi[i, j, k] = 0
            for kk in -1:1, jj in -1:1, ii in -1:1
                hi[i, j, k] += 1//27*hi_tmp[i+ii, j+jj, k+kk]
            end
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

function refine_field_int_ref!(hi, hi_tmp, lo, x_hi, xh_hi, x_lo, xh_lo, n_hi, n_lo)
    interp = interpolate((xh_lo, x_lo, x_lo), lo, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    # interp = interpolate((xh_lo, x_lo, x_lo), lo, (Gridded(Constant()), Gridded(Constant()), Gridded(Constant())))
    hi[2:end-1, 2:end-1, 2:end-1] = interp(xh_hi[2:end-1], x_hi[2:end-1], x_hi[2:end-1])
end

## Set up the grids.
n_hi = 384
n_lo = n_hi ÷ 2

a_lo = rand(n_lo, n_lo, n_lo)
a_hi = zeros(n_hi, n_hi, n_hi)

a_hi_gc = zeros(n_hi+2, n_hi+2, n_hi+2)
a_lo_gc = zeros(n_lo+2, n_lo+2, n_lo+2)

a_hi_gc[2:n_hi+1, 2:n_hi+1, 2:n_hi+1] = a_hi[:, :, :]
a_lo_gc[2:n_lo+1, 2:n_lo+1, 2:n_lo+1] = a_lo[:, :, :]

## Cyclic BCs in all 3 directions.
a_hi_gc[1, :, :] = a_hi_gc[end-1, :, :]
a_hi_gc[end, :, :] = a_hi_gc[2, :, :]
a_hi_gc[:, 1, :] = a_hi_gc[:, end-1, :]
a_hi_gc[:, end, :] = a_hi_gc[:, 2, :]
a_hi_gc[:, :, 1] = a_hi_gc[:, :, end-1]
a_hi_gc[:, :, end] = a_hi_gc[:, :, 2]

a_lo_gc[1, :, :] = a_lo_gc[end-1, :, :]
a_lo_gc[end, :, :] = a_lo_gc[2, :, :]
a_lo_gc[:, 1, :] = a_lo_gc[:, end-1, :]
a_lo_gc[:, end, :] = a_lo_gc[:, 2, :]
a_lo_gc[:, :, 1] = a_lo_gc[:, :, end-1]
a_lo_gc[:, :, end] = a_lo_gc[:, :, 2]

dx_hi = 1/n_hi; dx_lo = 1/n_lo
x_hi = dx_hi/2:dx_hi:1 |> collect
x_lo = dx_lo/2:dx_lo:1 |> collect
xh_hi = 0:dx_hi:1-dx_hi/2 |> collect
xh_lo = 0:dx_lo:1-dx_lo/2 |> collect

x_hi_gc = -dx_hi/2:dx_hi:1+dx_hi |> collect
x_lo_gc = -dx_lo/2:dx_lo:1+dx_lo |> collect
xh_hi_gc = -dx_hi:dx_hi:1+dx_hi/2 |> collect
xh_lo_gc = -dx_lo:dx_lo:1+dx_lo/2 |> collect

## Compute.
a_hi_tmp = zeros(size(a_hi_gc))
@btime refine_field_int!(a_hi_gc, a_hi_tmp, a_lo_gc, n_hi, n_lo)
a_hi_int = a_hi_gc[2:end-1, 2:end-1, 2:end-1]

## Compute a reference using Interpolations.jl
a_hi_tmp = zeros(size(a_hi_gc))
a_hi_gc[:, :, :] .= 0.
@btime refine_field_int_ref!(a_hi_gc, a_hi_tmp, a_lo_gc, x_hi_gc, xh_hi_gc, x_lo_gc, xh_lo_gc, n_hi, n_lo)
a_hi_int_ref = a_hi_gc[2:end-1, 2:end-1, 2:end-1]

println("Mean equal to lo: ", mean(a_lo) ≈ mean(a_hi_int))
println("Mean equal to ref: ", mean(a_hi_int_ref) ≈ mean(a_hi_int))
println("Interpolations equal: ", a_hi_int ≈ a_hi_int_ref)

x_hi_p = copy(x_hi)
x_lo_p = copy(x_lo)
pushfirst!(x_hi_p, 0)
pushfirst!(x_lo_p, 0)
xh_hi_p = copy(xh_hi)
xh_lo_p = copy(xh_lo)
push!(xh_hi_p, 1)
push!(xh_lo_p, 1)

## Plot the output.
figure(figsize=(10, 5))
subplot(121)
pcolormesh(x_lo_p, xh_lo_p, a_lo[:, :, 1]', vmin=0, vmax=1)
xlim(0, 1)
ylim(0, 1)
subplot(122)
pcolormesh(x_hi_p, xh_hi_p, a_hi_int[:, :, 1]', vmin=0, vmax=1)
xlim(0, 1)
ylim(0, 1)
tight_layout()
display(gcf())
show()


text:
for i in [Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,Bool]
    println(typemin(i),' ',typemax(i))
end

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
using Random
Random.seed!(76543210)

module PureVecs
    using Test, TestExtras
    using LinearAlgebra
    using Random
    using KrylovKit

    precision(T::Type{<:Number}) = eps(real(T))^(2/3)
    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    wrapvec(v) = v
    unwrapvec(v) = v
    wrapvec2(v) = v
    unwrapvec2(v) = v
    wrapop(A::AbstractMatrix) = A

    t = time()
    include("factorize.jl")
    include("gklfactorize.jl")
    include("linsolve.jl")
    include("eigsolve.jl")
    include("schursolve.jl")
    include("geneigsolve.jl")
    include("svdsolve.jl")
    include("expintegrator.jl")
    t = time() - t
    println("Julia Vector type: tests finisthed in $t seconds")
end

module MinimalVecs
    using Test, TestExtras
    using LinearAlgebra
    using Random
    using KrylovKit

    precision(T::Type{<:Number}) = eps(real(T))^(2/3)
    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    include("minimalvec.jl")

    wrapvec(v) = MinimalVec(v)
    unwrapvec(v::MinimalVec) = getindex(v)
    wrapvec2(v) = MinimalVec(v)
    unwrapvec2(v::MinimalVec) = getindex(v)
    wrapop(A::AbstractMatrix) = function (v, flag = Val(false))
        if flag === Val(true)
            return wrapvec(A'*unwrapvec2(v))
        else
            return wrapvec2(A*unwrapvec(v))
        end
    end

    t = time()
    include("factorize.jl")
    include("gklfactorize.jl")
    include("linsolve.jl")
    include("eigsolve.jl")
    include("schursolve.jl")
    include("geneigsolve.jl")
    include("svdsolve.jl")
    include("expintegrator.jl")
    t = time() - t
    println("Minimal vector type: tests finisthed in $t seconds")
end

module MixedSVD
    using Test, TestExtras
    using LinearAlgebra
    using Random
    using KrylovKit

    precision(T::Type{<:Number}) = eps(real(T))^(2/3)
    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    include("minimalvec.jl")

    wrapvec(v) = MinimalVec(v)
    unwrapvec(v::MinimalVec) = getindex(v)
    wrapvec2(v) = reshape(v, (length(v), 1)) # vector type 2 is a n x 1 Matrix
    unwrapvec2(v) = reshape(v, (length(v),))
    wrapop(A::AbstractMatrix) =
        (x->wrapvec2(A * unwrapvec(x)), y->wrapvec(A' * unwrapvec2(y)))

    t = time()
    include("gklfactorize.jl")
    include("svdsolve.jl")
    t = time() - t
    println("Mixed vector type for GKL/SVD: tests finisthed in $t seconds")
end

module ExtrasTest
    using Test, TestExtras
    using LinearAlgebra
    using Random
    using KrylovKit

    precision(T::Type{<:Number}) = eps(real(T))^(2/3)
    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    include("linalg.jl")
    include("recursivevec.jl")
end

if VERSION >= v"1.6"
    include("ad.jl")
end

module AquaTests
    using KrylovKit
    using Aqua
    Aqua.test_all(KrylovKit; ambiguities = false)
    # treat ambiguities special because of ambiguities between ChainRulesCore and Base
    if VERSION >= v"1.6" # ChainRulesCore leads to more ambiguities on julia < v1.6
        Aqua.test_ambiguities([KrylovKit,Base,Core]; exclude=[Base.:(==)])
    end
end


text:
using FromFile
@from "../types.jl" import MyInt
@from "utils.jl" import add_same_rand

add_same_rand2(x::MyInt)::MyInt = add_same_rand(x)


text:
#
# Example script to wrap PETSc
#

using Clang.cindex
using Clang.wrap_c

PETSC_INCLUDE = abspath("/usr/include/petsc")
MPI_INCLUDE = "/usr/include/openmpi"
JULIA_ROOT=abspath(JULIA_HOME, "../../")

LLVM_VER = "3.4"
LLVM_BUILD_TYPE = "Release+Asserts"

LLVM_PATH = joinpath(JULIA_ROOT, "deps/llvm-$LLVM_VER")
clanginc_path = joinpath(LLVM_PATH, "build_$LLVM_BUILD_TYPE/tools/clang/include/clang")

petsc_hdrs = [joinpath(PETSC_INCLUDE, "petsc.h")]

clang_includes = map(x::ASCIIString->joinpath(LLVM_PATH, x), [
    "build_$LLVM_BUILD_TYPE/tools/clang/include/clang",
    "include/llvm",
    "include"
    ])
push!(clang_includes, PETSC_INCLUDE)
push!(clang_includes, MPI_INCLUDE)
clang_extraargs = ["-v"]
# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

# Callback to test if a header should actually be wrapped (for exclusion)
function should_wrap(hdr::ASCIIString, name::ASCIIString)
    return beginswith(dirname(hdr), PETSC_INCLUDE)
end

function lib_file(hdr::ASCIIString)
    return "petsc"
end

function output_file(hdr::ASCIIString)
    return "PETSc.jl"
end

function should_wrap_cu(name::ASCIIString, cursor)
    exc = false
    exc |= contains(name, "MPI")
    return !exc
end

const wc = wrap_c.init(;
                        output_file = "libPETSc_h.jl",
                        common_file = "libPETSc_common.jl",
                        clang_includes = clang_includes,
                        clang_args = clang_extraargs,
                        header_wrapped = should_wrap,
                        header_library = lib_file,
                        header_outputfile = output_file,
                        cursor_wrapped = should_wrap_cu)
push!(wc.cache_wrapped, "ompi_file_errhandler_fn")
function wrap_libPETSc(wc::WrapContext, wrap_hdrs)
    wrap_c.wrap_c_headers(wc, wrap_hdrs)
end


text:
using QuasiGeostrophy, Test
include(pwd() * "/test/test_utils.jl")
@wrapper u=1 σ=1

@macroexpand  equ = @to_equation σ=u
equ = @to_equation σ=u
@pde_system pde_system = [
    σ = ∂x(u),
    ∂t(u) = -∂x(u * u - ∂x(σ)),
]

@testset "Basic Equation Test" begin
    @pde_system pde_system = [
        σ = ∂x(u),
        ∂t(u)= -∂x(u * u - ∂x(σ)),
    ]
    @test pde_system[1].rhs == ∂x(u)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module ArrayViews

using Compat

import Base: eltype, ndims, size, length, stride, strides
import Base: to_index, getindex, setindex!, parent, similar
import Base: Ptr, pointer

if VERSION >= v"0.4.0-dev+2085"
    import Base: iscontiguous
end

if VERSION < v"0.4.0-dev+3768"
    import Base: convert
else
    import Base: unsafe_convert
end

export
    StridedArrayView,
    ArrayView,
    ContiguousView,
    StridedView,
    UnsafeArrayView,
    UnsafeContiguousView,
    UnsafeStridedView,

    ContiguousArray,
    ContiguousVector,
    ContiguousMatrix,

    contiguous_view,
    strided_view,
    view,
    unsafe_view,
    ellipview,
    diagview,
    rowvec_view,
    flatten_view,
    reshape_view,
    iscontiguous,
    contiguousrank

## source files

include("common.jl")
include("arrviews.jl")
include("layout.jl")
include("indexing.jl")
include("contrank.jl")
include("subviews.jl")
include("convenience.jl")

include("deprecates.jl")

end  # module ArrayViews


text:
## Run this file using LIKWID:
# ```bash
# likwid-perfctr -C 0 -g MEM_DP -m julia --check-bounds=no --threads=1 measure_volume_terms.jl
# ```
# You can also inspect other groups, e.g., `-g FLOPS_DP`.

## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## load packages
using LIKWID
using Trixi

## set up benchmark code
"""
    initial_condition_isentropic_vortex(x, t, equations)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws.
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1])
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations3D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0, 0.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1], 0.0)
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

vortex_center(x, L) = mod(x + L/2, L) - L/2

function many_volume_terms!(du, u, semi, t, n)
  nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)
  for _ in 1:n
    Trixi.calc_volume_integral!(du, u, semi.mesh,
      nonconservative_terms, semi.equations,
      semi.solver.volume_integral, semi.solver, semi.cache)
  end
end

function run_measurements(polydeg=3, n=5*10^3)

  equations = CompressibleEulerEquations3D(1.4)

  Marker.init()

  for numerical_flux in (flux_shima_etal_turbo, flux_ranocha_turbo)

    solver = DGSEM(polydeg=polydeg, surface_flux=numerical_flux,
                   volume_integral=VolumeIntegralFluxDifferencing(numerical_flux))

    coordinates_min = (-5.0, -5.0, -5.0)
    coordinates_max = ( 5.0,  5.0,  5.0)
    initial_refinement_level = 3
    cells_per_dimension = 2^initial_refinement_level .* (1, 1, 1)
    t = 0.0

    # TreeMesh
    let
      mesh = TreeMesh(coordinates_min, coordinates_max,
                      initial_refinement_level=initial_refinement_level,
                      n_cells_max=100_000, periodicity=true)

      semi = SemidiscretizationHyperbolic(mesh, equations,
                                          initial_condition_isentropic_vortex, solver)

      u_ode = Trixi.compute_coefficients(t, semi)
      du_ode = zero(u_ode)

      GC.@preserve u_ode du_ode begin
        u = Trixi.wrap_array(u_ode, semi)
        du = Trixi.wrap_array(du_ode, semi)

        # compile and cool down
        many_volume_terms!(du, u, semi, t, 1)
        sleep(1.0)

        @region "TreeMesh-$(numerical_flux)" begin
          many_volume_terms!(du, u, semi, t, n)
        end
      end
    end

    # StructuredMesh
    let
      mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
                            periodicity=true)

      semi = SemidiscretizationHyperbolic(mesh, equations,
                                          initial_condition_isentropic_vortex, solver)

      u_ode = Trixi.compute_coefficients(t, semi)
      du_ode = zero(u_ode)

      GC.@preserve u_ode du_ode begin
        u = Trixi.wrap_array(u_ode, semi)
        du = Trixi.wrap_array(du_ode, semi)

        # compile and cool down
        many_volume_terms!(du, u, semi, t, 1)
        sleep(1.0)

        @region "StructuredMesh-$(numerical_flux)" begin
          many_volume_terms!(du, u, semi, t, n)
        end
      end
    end

    # P4estMesh
    let
      mesh = P4estMesh(cells_per_dimension; coordinates_min, coordinates_max,
                       polydeg=1, periodicity=true)

      semi = SemidiscretizationHyperbolic(mesh, equations,
                                          initial_condition_isentropic_vortex, solver)

      u_ode = Trixi.compute_coefficients(t, semi)
      du_ode = zero(u_ode)

      GC.@preserve u_ode du_ode begin
        u = Trixi.wrap_array(u_ode, semi)
        du = Trixi.wrap_array(du_ode, semi)

        # compile and cool down
        many_volume_terms!(du, u, semi, t, 1)
        sleep(1.0)

        @region "P4estMesh-$(numerical_flux)" begin
          many_volume_terms!(du, u, semi, t, n)
        end
      end
    end

  end

  Marker.close()

  return nothing
end

run_measurements()


text:
# calculate the least squares line of the errors

dirs = Array(1:3)  # do directories 1 through 3.  This range does not have
                   # to start at 1, but it must be contiguous
npts = length(dirs)
err_vals = zeros(maximum(dirs))
h_vals = zeros(maximum(dirs))
for i in dirs
  println("i = ", i)
  fname = "m$i/error_calc.dat"  # this file is written by the solver
  val_i = readdlm(fname)
  err_vals[i] = val_i[1]  # get error value
  h_vals[i] = val_i[2]  # get nominal mesh spacing h

end

min_idx = minimum(dirs)
max_idx = maximum(dirs)
err_vals = err_vals[min_idx:max_idx]
h_vals = h_vals[min_idx:max_idx]
npts = max_idx - min_idx + 1

# print the raw data
println("h  |   error")
println(hcat(h_vals, err_vals))

# tag the log of the data points
err_vals_log = zeros(npts)
h_vals_log = zeros(npts)
for i=1:npts
  err_vals_log[i] = log(err_vals[i])
  h_vals_log[i] = log(h_vals[i])
end

# now form system
A = ones(npts, 2)
A[:, 2] = h_vals_log

# solve system
AAt = A.'*A
rhs = A.'*err_vals_log

x = AAt\rhs

# compute error of each point compared to best fit line
for i=1:npts
  h_i = h_vals_log[i]
  err_exp = x[1] + h_i*x[2]
  err_act = err_vals_log[i]
  err_pct = 100*(err_exp - err_act)/err_act
  println("point $i percent error of best fit line = ", err_pct, "%")
end

# calculate slope from last two points
slope2 = (err_vals_log[npts - 1] - err_vals_log[npts])/(h_vals_log[npts - 1] - h_vals_log[npts])

println("2 point slope: ", slope2)
println("least squares slope: ", x[2])

# write data to file for plotting
data_arr = hcat(h_vals, err_vals)
writedlm("err_data.dat", data_arr, ' ')


text:
@testset "Discretization" begin
  @testset "FIST" begin
    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4,5]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == []

    𝒫 = Chain(P2[(0.443339268495331, 0.283757618605357),
                 (0.497822414616971, 0.398142813114205),
                 (0.770343126156527, 0.201815462842808),
                 (0.761236456732531, 0.330085709922366),
                 (0.985658085510286, 0.221530395507904),
                 (0.877899962498139, 0.325516131702896),
                 (0.561404274882782, 0.540334008885703),
                 (0.949459768187313, 0.396227653478068),
                 (0.594962560615951, 0.584927547374551),
                 (0.324208409133154, 0.607290684450708),
                 (0.424085089823892, 0.493532112641353),
                 (0.209843417261654, 0.590030658255966),
                 (0.27993878548962, 0.525162463476181),
                 (0.385557753911967, 0.322338556632868),
                 (0.443339268495331, 0.283757618605357)])
    @test Meshes.ears(𝒫) == [1,3,5,6,8,10,12,14]

    points = P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)]
    connec = connect.([(4,5,6),(3,4,6),(3,6,1),(1,2,3)], Triangle)
    target = UnstructuredMesh(points[1:end-1], connec)
    𝒫 = PolyArea(points)
    mesh = discretize(𝒫, FIST())
    @test mesh == target

    poly = readpoly(joinpath(datadir, "poly1.line"))
    points = P2.(coordinates.(vertices(first(chains(poly)))))
    connec = connect.([(26,27,28), (26,28,29), (23,24,25),
                       (23,25,26), (23,26,29), (23,29,30),
                       (21,22,23), (20,21,23), (16,17,18),
                       (16,18,19), (19,20,23), (19,23,30),
                       (13,14,15), (15,16,19), (11,12,13),
                       (13,15,19), (13,19,30), (8,9,10),
                       (7,8,10), (6,7,10), (4,5,6),
                       (4,6,10), (4,10,11), (11,13,30),
                       (11,30,1), (3,4,11), (3,11,1),
                       (1,2,3)], Triangle)
    target = UnstructuredMesh(points, connec)
    𝒫 = PolyArea([points; points[1]])
    mesh = discretize(𝒫, FIST())
    @test mesh == target

    poly = readpoly(joinpath(datadir, "poly2.line"))
    points = P2.(coordinates.(vertices(first(chains(poly)))))
    connec = connect.([(28, 29, 30), (24, 25, 26), (23, 24, 26),
                       (22, 23, 26), (21, 22, 26), (21, 26, 27),
                       (21, 27, 28), (19, 20, 21), (19, 21, 28),
                       (19, 28, 30), (18, 19, 30), (15, 16, 17),
                       (17, 18, 30), (17, 30, 1), (17, 1, 2),
                       (14, 15, 17), (14, 17, 2), (11, 12, 13),
                       (13, 14, 2), (10, 11, 13), (9, 10, 13),
                       (8, 9, 13), (8, 13, 2), (7, 8, 2),
                       (6, 7, 2), (3, 4, 5), (5, 6, 2),
                       (2, 3, 5)], Triangle)
    target = UnstructuredMesh(points, connec)
    𝒫 = PolyArea([points; points[1]])
    mesh = discretize(𝒫, FIST())
    @test mesh == target

    poly = readpoly(joinpath(datadir, "poly3.line"))
    points = P2.(coordinates.(vertices(first(chains(poly)))))
    connec = connect.([(28, 29, 30), (26, 27, 28), (24, 25, 26),
                       (23, 24, 26), (23, 26, 28), (21, 22, 23),
                       (21, 23, 28), (21, 28, 30), (19, 20, 21),
                       (19, 21, 30), (19, 30, 1), (15, 16, 17),
                       (12, 13, 14), (14, 15, 17), (12, 14, 17),
                       (10, 11, 12), (9, 10, 12), (8, 9, 12),
                       (5, 6, 7), (4, 5, 7), (3, 4, 7),
                       (18, 19, 1), (18, 1, 2), (18, 2, 3),
                       (17, 18, 3), (17, 3, 7), (17, 7, 8),
                       (8, 12, 17)], Triangle)
    target = UnstructuredMesh(points, connec)
    𝒫 = PolyArea([points; points[1]])
    mesh = discretize(𝒫, FIST())
    @test mesh == target

    poly = readpoly(joinpath(datadir, "poly4.line"))
    points = P2.(coordinates.(vertices(first(chains(poly)))))
    connec = connect.([(24, 25, 26), (23, 24, 26), (23, 26, 27),
                       (23, 27, 28), (20, 21, 22), (22, 23, 28),
                       (22, 28, 29), (22, 29, 30), (17, 18, 19),
                       (16, 17, 19), (15, 16, 19), (13, 14, 15),
                       (13, 15, 19), (13, 19, 20), (20, 22, 30),
                       (11, 12, 13), (11, 13, 20), (11, 20, 30),
                       (11, 30, 1), (6, 7, 8), (6, 8, 9),
                       (6, 9, 10), (10, 11, 1), (10, 1, 2),
                       (4, 5, 6), (4, 6, 10), (4, 10, 2),
                       (2, 3, 4)], Triangle)
    target = UnstructuredMesh(points, connec)
    𝒫 = PolyArea([points; points[1]])
    mesh = discretize(𝒫, FIST())
    @test mesh == target

    poly = readpoly(joinpath(datadir, "poly5.line"))
    points = P2.(coordinates.(vertices(first(chains(poly)))))
    connec = connect.([(25, 26, 27), (24, 25, 27), (23, 24, 27),
                       (23, 27, 28), (23, 28, 29), (19, 20, 21),
                       (18, 19, 21), (17, 18, 21), (15, 16, 17),
                       (12, 13, 14), (7, 8, 9), (6, 7, 9),
                       (6, 9, 10), (4, 5, 6), (2, 3, 4),
                       (2, 4, 6), (1, 2, 6), (30, 1, 6),
                       (22, 23, 29), (22, 29, 30), (11, 12, 14),
                       (10, 11, 14), (10, 14, 15), (10, 15, 17),
                       (6, 10, 17), (6, 17, 21), (6, 21, 22),
                       (6, 22, 30)], Triangle)
    target = UnstructuredMesh(points, connec)
    𝒫 = PolyArea([points; points[1]])
    mesh = discretize(𝒫, FIST())
    @test mesh == target
  end
end


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

"""
    Line(a, b)

A line passing through points `a` and `b`.

See also [`Segment`](@ref).
"""
struct Line{Dim,T} <: Primitive{Dim,T}
  a::Point{Dim,T}
  b::Point{Dim,T}
end

Line(a::Tuple, b::Tuple) = Line(Point(a), Point(b))

paramdim(::Type{<:Line}) = 1

(l::Line)(t) = l.a + t * (l.b - l.a)

function Base.in(p::Point, l::Line)
  w = norm(l.b-l.a)
  d = evaluate(Euclidean(), p, l)
  # d ≈ 0.0 will be too precise, and d < atol{T} can't scale.
  d+w ≈ w
end

==(l1::Line, l2::Line) = l1.a ∈ l2 && l1.b ∈ l2 && l2.a ∈ l1 && l2.b ∈ l1


text:
function loadBrukerfolderautophase(full_path::String)

    pdata_path = "$full_path/pdata/1"

    PyCall.py"""
    import numpy
    import nmrglue as ng

    s = $full_path
    procdata_dir = $pdata_path

    dic, data = ng.bruker.read(s)

    # ---time correction before FT
    precorr_time =  ng.bruker.remove_digital_filter(dic, data)
    precorr_frq0 = ng.proc_base.fft(precorr_time)
    precorr_frq = ng.proc_autophase.autops(precorr_frq0, 'acme')
    # precorr_frq = ng.proc_autophase.autops(precorr_frq0, 'peak_minima')

    #---time correction after FT
    postcorr_frq0 = ng.proc_base.fft(data)
    postcorr_frq1 = ng.bruker.remove_digital_filter(dic, postcorr_frq0, post_proc=True)
    postcorr_frq = ng.proc_autophase.autops(postcorr_frq1, 'acme')

    #---data processed using TopSpin (3.5pl7)
    bruker_frq = ng.bruker.read_pdata(procdata_dir)[1][::-1]

    guess_udic = ng.fileio.bruker.guess_udic(dic, data)
    """

    dic = PyCall.py"dic"
    data = PyCall.py"data"

    precorr_time = PyCall.py"precorr_time"
    precorr_frq0 = PyCall.py"precorr_frq0"
    precorr_frq = PyCall.py"precorr_frq"

    postcorr_frq0 = PyCall.py"postcorr_frq0"
    postcorr_frq1 = PyCall.py"postcorr_frq1"
    postcorr_frq = PyCall.py"postcorr_frq"

    bruker_frq = PyCall.py"bruker_frq"
    guess_udic = PyCall.py"guess_udic"

    return dic, data, precorr_time, precorr_frq0, precorr_frq,
            postcorr_frq0, postcorr_frq1, postcorr_frq, bruker_frq,
            guess_udic
end

"""
    loadspectrum(full_path::String;
        solvent_ppm = 4.8,
        solvent_window_ppm = 0.3,
        N_Ω::Int = 20000,
        N_optim::Int = 1000,
        n_particles::Int = 3,
        max_iters::Int = 100,
        α_lower = 1e-6,
        α_upper_factor = 50.0,
        λ_lower = 1e-10,
        λ_upper = 20.0,
        λ_initial = 3.0,
        st_ind_default = 100,
        verbose_flag::Bool = false)

Loads a Bruker 1D 1H NMR experiment. Uses the NMRGlue Python library. Estimates the 0 ppm and solvent resonance components. The solvent resonance component is estimated in the frequency interval `solvent_ppm` +/- `solvent_window_ppm`, in units ppm.

On a Linux terminal, run the following command if you don't have NMRGlue installed.
pip install nmrglue --user

Notes:
Either `st_ind_default` samples or `dead_time_st` number of samples are removed from the raw time-domain data to create the output time-domain data, s_t. This is to remove the silence at the beginning of the acquired data due to the dead-time property of practical NMR spectrometers. 'dead_time_st' is calculated automatically from the metadata of the experiment files.

Therefore, the returned sampling frequency `fs` from loadspectrum() and the sampling frequency in the returned metadata ditionary `dic` are slightly different.

...
# Outputs, in the order returned.

- `s_t::Vector{Complex{Float64}}`: time-domain experiment data.
- `S::Function`: Discrete-time Fourier transform of s_t.
- `hz2ppmfunc::Function`: conversion function from Hz to ppm for this experiment. The formula used is `hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs`.
- `ppm2hzfunc::Function`: conversion function from ppm to Hz for this experiment. The formula used is `ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)`.

- `ν_0ppm::Float64`: where the estimated frequency of the 0 ppm resonance component, in Hz.
- `fs::Float64`: the estimated sampling frequency of `s_t`.
- `SW::Float64`: the spectral window of the spectrometer, as provided by the metadata of the experiment. In units of ppm.
- `α_0ppm`: The amplitude parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `β_0ppm`: The phase parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `λ_0ppm`: The T2 parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `Ω_0ppm`: The radial frequency parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model. In units Hz/(2*π).
- `results_0ppm`: The NLopt.jl optimization results for the estimation of the 0 ppm reference. Inspect this to see if the estimated parameters for the 0 ppm reference are reliable.

- `dic`; The dictionary that contains all metadata from NMRGlue for the NMR experiment.

- `α_solvent`: The amplitude parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `β_solvent`: The phase parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `λ_solvent`: The T2 parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model.
- `Ω_solvent`: The radial frequency parameter of the estimate resonance component for 0 ppm reference, in the free-induction decay model. In units Hz/(2*π).
- `results_solvent`: The NLopt.jl optimization results for the estimation of the 0 ppm reference. Inspect this to see if the estimated parameters for the 0 ppm reference are reliable.
, , , , ,

# Examples
See load_experiment.jl in the /examples folder for an example.
"""
function loadspectrum(full_path::String;
    solvent_ppm = 4.8,
    solvent_window_ppm = 0.3,
    N_Ω::Int = 20000,
    N_optim::Int = 1000,
    n_particles::Int = 3,
    max_iters::Int = 100,
    α_lower = 1e-6,
    α_upper_factor = 50.0,
    λ_lower = 1e-10,
    λ_upper = 20.0,
    λ_initial = 3.0,
    st_ind_default = 100, # future: see if we can derive this from DE dead time later.
    verbose_flag::Bool = false)

    dic, data, precorr_time, precorr_frq0, precorr_frq, postcorr_frq0, postcorr_frq1,
    postcorr_frq, Bruker_spectrums_nmrglue,
    guess_udic = loadBrukerfolderautophase(full_path)

    DE = dic["acqus"]["DE"] # in microseconds.

    BF1 = dic["acqus"]["BF1"]
    ## see Bruker TopSpin acquistion commands and parameters v 003.
    # TD - Time Domain; Number Of Raw Data Points. i.e., length(data)*2 - TD = 0.
    # total samples from both time-series.
    TD = dic["acqus"]["TD"]

    # SW - Spectral Width in ppm
    SW = dic["acqus"]["SW"]

    # SFO1 - SFO8 - Irradiation (carrier) Frequencies For Channels f1 to f8
    SFO1 = dic["acqus"]["SFO1"]

    O1 = dic["acqus"]["O1"]
    dic["acqus"]["NUC1"]

    CAR = O1*1.0

    fs_dic = dic["acqus"]["SW_h"]

    N_diff = length(data)*2 - TD

    t_data = gettimerange(length(data), fs_dic)

    AQ = TD/(2*SW*SFO1)
    dead_time = t_data[end] - AQ

    dead_time_st = round(Int, dead_time/t_data[end] * length(data))
    st_ind = max(st_ind_default, dead_time_st)

    s_t = data[st_ind:end]
    t = gettimerange(length(s_t), fs_dic)

    fs = (length(s_t)-1)/t[end] # this is really fs_dic.
    @assert abs(fs-fs_dic) < 1e-9

    α_0ppm, β_0ppm, λ_0ppm, Ω_0ppm,
    results_0ppm = estimatereferencecompoundfreq(CAR, fs, SW, s_t, t;
    N_Ω = N_Ω,
    N_optim = N_optim,
    n_particles = n_particles,
    max_iters = max_iters,
    α_lower = α_lower,
    α_upper_factor = α_upper_factor,
    λ_lower = λ_lower,
    λ_upper = λ_upper,
    λ_initial = λ_initial,
    verbose_flag = verbose_flag)

    ν_0ppm = Ω_0ppm ./ (2*π)

    hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
    ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

    # S is scaled DTFT of data.
    S = vv->(computeDTFTch3eq29(s_t, vv, t)/fs)

    # solvent.
    u_solvent_lower = ppm2hzfunc(solvent_ppm-solvent_window_ppm)
    u_solvent_upper = ppm2hzfunc(solvent_ppm+solvent_window_ppm)

    α_solvent, β_solvent, λ_solvent, Ω_solvent,
    results_solvent = estimatesinglet(fs, s_t, t,
                u_solvent_lower, u_solvent_upper;
                N_Ω = N_Ω,
                N_optim = N_optim,
                n_particles = n_particles,
                max_iters = max_iters,
                α_lower = α_lower,
                α_upper_factor = α_upper_factor,
                λ_lower = λ_lower,
                λ_upper = λ_upper,
                λ_initial = λ_initial,
                verbose_flag = verbose_flag)

    return s_t, S, hz2ppmfunc, ppm2hzfunc, ν_0ppm, fs, SW, α_0ppm, β_0ppm, λ_0ppm, Ω_0ppm,
    results_0ppm, dic, α_solvent, β_solvent, λ_solvent, Ω_solvent, results_solvent
end

function estimatereferencecompoundfreq( CAR::T,
    fs::T,
    SW::T,
    s_t::Vector{Complex{T}},
    t;
    N_Ω::Int = 20000,
    N_optim::Int = 1000,
    n_particles::Int = 3,
    max_iters::Int = 100,
    window_ppm::T = 0.3,
    α_lower::T = 1e-6,
    α_upper_factor::T = 50.0,
    λ_lower::T = 1e-10,
    λ_upper::T = 20.0,
    λ_initial::T = 3.0,
    verbose_flag::Bool = false) where T <: Real

    # prepare frequency range.
    ν0_initial, hz2ppmfunc0, ppm2hzfunc0 = getinitalguessreferencecompoundfreq(CAR, fs, SW)

    bp_a = ppm2hzfunc0(-window_ppm)
    bp_b = ppm2hzfunc0(window_ppm)

    return estimatesinglet(fs, s_t, t, bp_a, bp_b;
    N_Ω = N_Ω,
    N_optim = N_optim,
    n_particles = n_particles,
    max_iters = max_iters,
    α_lower = α_lower,
    α_upper_factor = α_upper_factor,
    λ_lower = λ_lower,
    λ_upper = λ_upper,
    λ_initial = λ_initial,
    verbose_flag = verbose_flag)

end

# λ_initial: m a value of 3 seems decent for the λ of DSS under 600 MHz.
# no error checking on whether λ_initial ∈ [λ_lower, λ_upper].
function estimatesinglet( fs::T,
    s_t::Vector{Complex{T}},
    t,
    u_lower,
    u_upper;
    N_Ω::Int = 20000,
    N_optim::Int = 1000,
    n_particles::Int = 3,
    max_iters::Int = 100,
    α_lower::T = 1e-6,
    α_upper_factor::T = 100.0,
    λ_lower::T = 1e-10,
    λ_upper::T = 20.0,
    λ_initial::T = 3.0,
    verbose_flag::Bool = false,
    xtol_rel::T = 1e-12,
    nlopt_seed::Int = 25) where T <: Real

    # prepare data.
    DTFT_s_scaled = vv->(computeDTFTch3eq29(s_t, vv, t)/fs)

    u_range_optim = LinRange(u_lower, u_upper, N_optim)
    S_U = DTFT_s_scaled.(u_range_optim)

    # set up optimization.
    max_val, max_ind = findmax(abs.(S_U))
    Ω = u_range_optim[max_ind] *2*π
    #println("Ω = ", Ω)
    costfunc = pp->evalcostreferencecompound(pp,
    S_U, u_range_optim, Ω)

    optim_lower_limit = [α_lower; -π; λ_lower]
    optim_upper_limit = [max_val*α_upper_factor; π; λ_upper]

    p0 = [max_val; 0.0; λ_initial]

    # optimize.
    op = Optim.Options( iterations = max_iters,
    store_trace = false,
    show_trace = verbose_flag)

    results = Optim.optimize( costfunc,
    p0,
    Optim.ParticleSwarm(; lower = optim_lower_limit,
        upper = optim_upper_limit,
        n_particles = n_particles),
    op)

    #α, β, λ = results.minimizer
    # println("α = ", α)
    # Q = [1221.2329692741823;
    # 1.5072497616902667;
    # 4.266537181116996]
    # println("cost Q = ", costfunc(Q))
    # println("cost p = ", costfunc(results.minimizer))

    # use NLopt.
    grad_func = xx->FiniteDiff.finite_difference_gradient(costfunc, xx)

    NLopt.srand(nlopt_seed)

    p_initial = results.minimizer
    #p_initial = p0

    opt = NLopt.Opt(:GN_ESCH, length(p_initial))
    #opt = NLopt.Opt(:LD_MMA, length(p_initial))

    min_cost_p, min_p, ret_p, numevals_p = runNLopt!(  opt,
        p_initial,
        costfunc,
        grad_func,
        optim_lower_limit,
        optim_upper_limit;
        max_iters = 100*max_iters,
        xtol_rel = xtol_rel)
    #
    opt = NLopt.Opt(:LD_MMA, length(min_p))
    min_cost_p, min_p, ret_p, numevals_p = runNLopt!(  opt,
        min_p,
        costfunc,
        grad_func,
        optim_lower_limit,
        optim_upper_limit;
        max_iters = 100*max_iters,
        xtol_rel = xtol_rel)

    α, β, λ = min_p

    # #
    # println("p_initial = ", p_initial)
    # println("α         = ", α)

    return α, β, λ, Ω, results
end

function evalcostreferencecompound( p::Vector{T},
    S_U::Vector{Complex{T}},
    u_range,
    Ω) where T <: Real

    # parse.
    α, β, λ = p

    cost = zero(T)
    for i = 1:length(u_range)
        q_u = evalcomplexLorentzian(u_range[i], α, β, λ, Ω)
        cost += abs2(q_u - S_U[i])
        #cost += (abs(q_u)-abs(S_U[i]))^2
    end

    return cost
end

function getinitalguessreferencecompoundfreq(CAR::T, fs::T, SW::T) where T <: Real

    ν0_initial = fs - CAR

    hz2ppmfunc0 = uu->(uu - ν0_initial)*SW/fs
    ppm2hzfunc0 = pp->(ν0_initial + pp*fs/SW)

    return ν0_initial, hz2ppmfunc0, ppm2hzfunc0
end


text:
@testset "kDPP" begin
    seed!(42)
    N = 20
    D = 3
    nInd = 10
    kernel = SqExponentialKernel()
    X = ColVecs(rand(D, N))
    alg = kDPP(nInd, kernel)
    @test repr(alg) == "k-DPP selection of inducing points"
    Z = inducingpoints(alg, X)
    @test length(Z) == nInd
    @test_throws ArgumentError kDPP(-1, kernel)
    @test_throws ErrorException inducingpoints(kDPP(100, kernel), X)
end


text:
config = MOI.Test.TestConfig()
_optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 1.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1.0]])
mock = bridged_mock(_optimize!)
Tests.zero_polynomial_in_fixed_variables_set_test(mock, config)


text:
include("1.jl")
function mark_allzic(r::Robot)
    angleofsquare1111(r::Robot)
    putmarker!(r)
    snake11(r)
    if (isborder(r,Ost))
        putmarkhole(r,West)
    else
        putmarkhole(r,inverse(West))
    end
end

function snake11(r::Robot)
    while !isborder(r,Nord)
        if (isborder(r,Ost))
            putmarkhole(r,West)
        else
            putmarkhole(r,inverse(West))
        end
        if (isborder(r,Ost) || isborder(r,West))
            move!(r,Nord)
            putmarker!(r)
        end
    end

end
function movehole11111(r::Robot,side::HorizonSide)
    while (!isborder(r,side))
        move!(r,side)
        if (isborder(r,Nord))
            break
        end
    end

end
function angleofsquare1111(r::Robot)
    while ((!isborder(r,Sud)) || (!isborder(r,West)))
        if (!isborder(r,Sud))
            move!(r,Sud)
        end
        if (!isborder(r,West))
            move!(r,West)
        end
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
############################################################################################
"""
$(TYPEDEF)
Holds information for array structure of data and reference particle.

# Fields
$(TYPEDFIELDS)
"""
mutable struct ParticleFilterConfig{A,B,C,D}
    "Data size and sorting information."
    data::BaytesCore.ArrayConfig{A,B}
    "Particle size and sorting information."
    particle::BaytesCore.ArrayConfig{C,D}
    function ParticleFilterConfig(
        dataconfig::ArrayConfig{A,B}, particleconfig::ArrayConfig{C,D}
    ) where {A,B,C,D}
        return new{A,B,C,D}(dataconfig, particleconfig)
    end
end
function ParticleFilterConfig(
    data::AbstractArray{T}, particle::AbstractArray{F}
) where {T,F}
    return ParticleFilterConfig(
        BaytesCore.ArrayConfig(data), BaytesCore.ArrayConfig(particle)
    )
end

function update!(
    config::ParticleFilterConfig, data::AbstractArray{T}, particle::AbstractArray{F}
) where {T,F}
    config.data = BaytesCore.ArrayConfig(data)
    config.particle = BaytesCore.ArrayConfig(particle)
    return nothing
end

############################################################################################
"""
$(TYPEDEF)
Holds tuning information for Particle Filter.

# Fields
$(TYPEDFIELDS)
"""
struct ParticleFilterTune{
    T<:Tagged,
    A<:BaytesCore.ParameterWeighting,B<:BaytesCore.ResamplingMethod,C<:ParticleReferencing,
    D,E,F,G,
    U<:BaytesCore.UpdateBool
} <: AbstractTune
    "Tagged Model parameter."
    tagged::T
    "Weighting Methods for particles."
    weighting::A
    "Resampling methods for particle trajectories."
    resampling::B
    "Referencing type for last particle at each iteration - either Conditional, Ancestral or Marginal Implementation."
    referencing::C
    "Contains data and reference size and sorting."
    config::ParticleFilterConfig{D,E,F,G}
    "Number of particle chains and tuning information."
    chains::BaytesCore.ChainsTune
    "Memory for latent and observed data."
    memory::ParticleFilterMemory
    "Boolean if generated quantities should be generated while sampling"
    generated::U
    "Current iteration."
    iter::Iterator
    function ParticleFilterTune(
        objective::O,
        weighting::A,
        resampling::B,
        referencing::C,
        config::ParticleFilterConfig{D,E,F,G},
        chains::BaytesCore.ChainsTune,
        memory::ParticleFilterMemory,
        generated::U,
        iter::Iterator,
    ) where {
        O<:Objective,
        A<:BaytesCore.ParameterWeighting,
        B<:BaytesCore.ResamplingMethod,
        C<:ParticleReferencing,
        D,
        E,
        F,
        G,
        U<:BaytesCore.UpdateBool
    }
        return new{typeof(objective.tagged),A,B,C,D,E,F,G,U}(
            objective.tagged,
            weighting,
            resampling,
            referencing,
            config,
            chains,
            memory,
            generated,
            iter,
        )
    end
end

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

function update!(tune::ParticleFilterTune, data::AbstractArray, reference::AbstractArray)
    ## Update Iteration count
    init!(tune.iter, 1)
    ## Update Data/Particle Configuration
    update!(tune.config, data, reference)
    Ndata = maximum(tune.config.data.size)
    ## Update Number of Particles
    update_Nparticles, update_Ndata = update!(tune.chains, Ndata)
    ## Return if adjustments in either dimension are needed
    return update_Nparticles, update_Ndata, Ndata
end

############################################################################################
#export
export ParticleFilterTune


text:
type Sink <: Data.Sink
    schema::Data.Schema
    dsn::DSN
    table::String
    columns::Vector{Any}
    indcols::Vector{Any}
end

# DataStreams interface
function Sink{T}(source, ::Type{T}, append::Bool, dsn::DSN, table::AbstractString)
    cols = size(source, 2)
    sink = Sink(Data.schema(source), dsn, table, Vector{Any}(cols), Vector{Any}(cols))
    !append && ODBC.execute!(dsn, "delete from $table")
    stmt = sink.dsn.stmt_ptr2
    ODBC.execute!(sink.dsn, "select * from $(sink.table)", stmt)
    return sink
end
function Sink{T}(sink, source, ::Type{T}, append::Bool)
    !append && ODBC.execute!(sink.dsn, "delete from $(sink.table)")
    stmt = sink.dsn.stmt_ptr2
    ODBC.execute!(sink.dsn, "select * from $(sink.table)", stmt)
    return sink
end

Data.streamtypes{T<:ODBC.Sink}(::Type{T}) = [Data.Column]

prep!{T}(::Type{T}, A) = A, 0
prep!{T}(::Type{Nullable{T}}, A) = A.values, 0
prep!(::Union{Type{Date},Type{Nullable{Date}}}, A) = ODBC.API.SQLDate[isnull(x) ? ODBC.API.SQLDate() : ODBC.API.SQLDate(x) for x in A], 0
prep!(::Union{Type{DateTime},Type{Nullable{DateTime}}}, A) = ODBC.API.SQLTimestamp[isnull(x) ? ODBC.API.SQLTimestamp() : ODBC.API.SQLTimestamp(x) for x in A], 0
if is_unix()
prep!(::Union{Type{Dec64},Type{Nullable{Dec64}}}, A) = Float64[isnull(x) ? 0.0 : Float64(get(x)) for x in A], 0
end

getptrlen(x::AbstractString) = pointer(x.data), length(x), UInt8[]
getptrlen{T}(x::WeakRefString{T}) = convert(Ptr{UInt8}, x.ptr), codeunits2bytes(T, x.len), UInt8[]
getptrlen{T}(x::Nullable{T}) = isnull(x) ? (convert(Ptr{UInt8}, C_NULL), 0, UInt8[]) : getptrlen(get(x))
function getptrlen(x::CategoricalArrays.CategoricalValue)
    ref = String(x).data
    return pointer(ref), length(ref), ref
end

prep!{T<:AbstractString}(::Type{T}, A) = _prep!(T, A)
prep!{T<:AbstractString}(::Type{Nullable{T}}, A) = _prep!(T, A)
prep!{T<:CategoricalArrays.CategoricalValue}(::Type{T}, A) = _prep!(T, A)
prep!{T<:CategoricalArrays.CategoricalValue}(::Type{Nullable{T}}, A) = _prep!(T, A)

function _prep!{T}(::Type{T}, column)
    maxlen = maximum(clength, column)
    maxlen = typeof(maxlen) <: Nullable ? get(maxlen) : maxlen
    data = zeros(UInt8, maxlen * length(column))
    ind = 1
    for i = 1:length(column)
        ptr, len, ref = getptrlen(column[i])
        unsafe_copy!(pointer(data, ind), ptr, len)
        ind += maxlen
    end
    return data, maxlen
end

function prep!{T}(source, ::Type{T}, col, columns, indcols)
    column = Data.getcolumn(source, T, col)
    columns[col], maxlen = prep!(T, column)
    indcols[col] = ODBC.API.SQLLEN[clength(x) for x in column]
    return length(column), maxlen
end

getCtype{T}(::Type{T}) = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)
getCtype{T}(::Type{Nullable{T}}) = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)

function Data.streamcolumn!{T}(sink::ODBC.Sink, source, ::Type{T}, col, row)
    stmt = sink.dsn.stmt_ptr2
    rows, len = prep!(source, T, col, sink.columns, sink.indcols)
    ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLBindCols(stmt, col, getCtype(T), sink.columns[col], len, sink.indcols[col])
    if col == size(source, 2)
        ODBC.API.SQLSetStmtAttr(stmt, ODBC.API.SQL_ATTR_ROW_ARRAY_SIZE, rows, ODBC.API.SQL_IS_UINTEGER)
        ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLBulkOperations(stmt, ODBC.API.SQL_ADD)
    end
    return rows
end

function load{T}(dsn::DSN, table::AbstractString, ::Type{T}, args...; append::Bool=false)
    source = T(args...)
    cols = size(source, 2)
    sink = Sink(Data.schema(source), dsn, table, Vector{Any}(cols), Vector{Any}(cols))
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end
function load(dsn::DSN, table::AbstractString, source; append::Bool=false)
    cols = size(source, 2)
    sink = Sink(Data.schema(source), dsn, table, Vector{Any}(cols), Vector{Any}(cols))
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end

load{T}(sink::Sink, ::Type{T}, args...; append::Bool=false) = (sink = Data.stream!(T(args...), sink, append); Data.close!(sink); return sink)
load(sink::Sink, source; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)

# function Data.stream!(source, ::Type{Data.Column}, sink::ODBC.Sink, append::Bool=false)
#     Data.types(source) == Data.types(sink) || throw(ArgumentError("schema mismatch: \n$(Data.schema(source))\nvs.\n$(Data.schema(sink))"))
#     rows, cols = size(source)
#     Data.isdone(source, 1, 1) && return sink
#     ODBC.execute!(sink.dsn, "select * from $(sink.table)")
#     stmt = sink.dsn.stmt_ptr
#     types = Data.types(source)
#     columns = Vector{Any}(cols)
#     indcols = Array{Vector{ODBC.API.SQLLEN}}(cols)
#     row = 0
#     # get the column names for a table from the DB to generate the insert into sql statement
#     # might have to try quoting
#     # SQLPrepare (hdlStmt, (SQLTCHAR*)"INSERT INTO customers (CustID, CustName,  Phone_Number) VALUES(?,?,?)", SQL_NTS) ;
#     try
#         # SQLSetConnectAttr(hdlDbc, SQL_ATTR_AUTOCOMMIT, SQL_AUTOCOMMIT_OFF, SQL_NTS)
#         while !Data.isdone(source, row+1, cols+1)
#
#             for col = 1:cols
#                 T = types[col]
#                 # SQLBindParameter(hdlStmt, 1, SQL_PARAM_INPUT, SQL_C_LONG, SQL_INTEGER, 0, 0, (SQLPOINTER)custIDs, sizeof(SQLINTEGER) , NULL);
#                 rows, cT = ODBC.bindcolumn!(source, T, col, columns, indcols)
#                 ret = ODBC.API.SQLBindCols(stmt, col, cT, pointer(columns[col]), sizeof(eltype(columns[col])), indcols[col])
#                 println("$col: $ret")
#             end
#             ODBC.API.SQLSetStmtAttr(stmt, ODBC.API.SQL_ATTR_ROW_ARRAY_SIZE, rows, ODBC.API.SQL_IS_UINTEGER)
#             # SQLSetStmtAttr( hdlStmt, SQL_ATTR_PARAMSET_SIZE, (SQLPOINTER)NUM_ENTRIES, 0 );
#             # ret = SQLExecute(hdlStmt);
#             ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLBulkOperations(stmt, ODBC.API.SQL_ADD)
#             row += rows
#         end
#         # SQLEndTran(SQL_HANDLE_DBC, hdlDbc, SQL_COMMIT);
#     # finally
#         # SQLSetConnectAttr(hdlDbc, SQL_ATTR_AUTOCOMMIT, SQL_AUTOCOMMIT_ON, SQL_NTS);
#     end
#     Data.setrows!(source, row)
#     return sink
# end


text:
using QuadGK: quadgk

export fourier_coefficient, FourierSeries, eval_series, rotating_frame_series

"""
    fourier_coefficient(f::Function, frequency::Real, harmonic::Int)

Calculate the fourier coefficient `1/T ∫_0^T f(τ) e^(-2πi * frequency * harmonic * τ) dτ`
of a periodic function `f`.

# args
* `f`: the function of one variable.
* `frequency`: the frequency of periodicity of `f`.
* `harmonic`: the desired harmonic.

# returns
The fourier coefficient.
"""
function fourier_coefficient(f::Function, frequency::Real, harmonic::Int)
    T = 1/frequency
    g(t) = f(t) * exp(-2π * 1im * frequency * harmonic * t)
    integral, _ = quadgk(g, 0.0, T, atol=sqrt(eps(T)))
    return integral/T
end

"""
Representation of a Fourier series `f(t) = ∑ Aₙ e^(iωnt)`
"""
struct FourierSeries
    frequency::Real
    terms::Dict{Int, <:Number} # harmonic => coefficient
end

"""
    FourierSeries(f::Function, frequency::Real, harmonics::AbstractVector{Int})

Construct a FourierSeries from a periodic function `f` on given harmonics.

# args
* `f`: the periodic function.
* `frequency`: the frequency of periodicity of f.
* `harmonics`: a vector of harmonics to compute.

# returns
A FourierSeries.
"""
function FourierSeries(f::Function, frequency::Real, harmonics::AbstractVector{Int})
    terms = Dict(h => fourier_coefficient(f, frequency, h) for h in harmonics)
    return FourierSeries(frequency, terms)
end

"""
    eval_series(fs::FourierSeries, times::AbstractVector{<:Real})

Evaluate a FourierSeries on given times.

# args
* `fs`: a FourierSeries.
* `times`: a vector of times.

# returns
A vector of values.
"""
function eval_series(fs::FourierSeries, times::AbstractVector{<:Real})
    return reduce(+, coeff * exp.(2π * 1im * fs.frequency * harmonic * times) for (harmonic, coeff) in fs.terms)
end

"""
    rotating_frame_series(fs::FourierSeries, harmonics::AbstractVector{Int})

Compute the fourier series for `exp(i∫_0^t f(τ) dτ)` where `f` is the function given by
the FourierSeries `fs`. This function ignores a constant term in `fs` since that maps onto
an overall exponential factor which is not part of a FourierSeries.

# args
* `fs`: a FourierSeries.
* `harmonics`: the desired harmonics.

# returns
The FourierSeries.
"""
function rotating_frame_series(fs::FourierSeries, harmonics::AbstractVector{Int})
    iω = 2π * 1im * fs.frequency
    function exp_of_integral(t::Real)
        int = 0.0
        for (h, c) in fs.terms
            if h != 0
                int += c * (exp(iω * h * t) - 1)/(iω * h)
            end
        end
        return exp(1im * int)
    end
    return FourierSeries(exp_of_integral, fs.frequency, harmonics)
end


text:
# equations 11-13 for spatial frequencies converted as cosh = e^x - e^-x /2
function fluence_spatial_freq(z, s, μa, μsp, l)

    D = @. 1/3μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    if z < l

        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

        if z < z0
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))
        end

    else

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    if isnan(ϕ) || isinf(ϕ)
        ϕ = 0
    end

    return ϕ
end

# compute equation 14 integrand
function compute_integral(ρ, z, s, μa, μsp, l)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

# numerically integrate equation 14
function fluence_steadystate(ρ, z, μa, μsp, l)

    ϕ = quadgk((s) -> compute_integral(ρ, z, s, μa, μsp, l), 0, Inf, rtol=1e-8)
    return ϕ[1]/2π

 end

function ss_compare(ρ, z, μa, μsp)

    D = 1/3μsp
    A = 1
    z0 = 1/(μsp + μa)
    zb = 2*A*D
    μe = sqrt(3*μa*(μa + μsp))

    ϕ = exp(-μe*sqrt((z - z0)^2 + ρ^2))/sqrt((z - z0)^2 + ρ^2)
    ϕ = ϕ - exp(-μe*sqrt((z + z0 + 2*zb)^2 + ρ^2))/sqrt((z + z0 + 2*zb)^2 + ρ^2)
    ϕ = ϕ/(4*π*D)
end

function fluence_spatial_freq(z, s, μa, μsp, l)

    D = @. 1/3μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])
    tolr = 1e-50
    ϕ = []

    ## equation 11
    if z < z0

        # limit as neg exp gets large
        if exp(-2*α[1]*z0) < tolr
            ϕ = 0.0
        elseif exp(-2*l) < tolr
            ϕ = exp(α[1]*z0 - α[1]*z)/(2*D[1]*α[1])
            ϕ = ϕ*(exp(-2*α[1]*(z0 - z)) - exp(-2*α[1]*(zb + z0)))
        elseif exp(-2*l*α[1]) < tolr
            ϕ = 0.0
        else
            ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
            ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
            ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
            ϕ = ϕ*exp(α[1]*z0 - α[1]*z)
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))
        end

    ## equation 12
    elseif z > z0 && z < l
        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

    #equation 13
    elseif z > l

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

function fluence_spatial_freq(z, s, μa, μsp, l)
    #get numerical instabillity at large s

    D = @. 1/(3(μsp + μa))
    α = @. sqrt((D*s^2 + μa)/D)
    A = 0.493
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])
    tolr = 1e-50
    ϕ = zeros(Float64, 1)

    ## equation 11
    if z < z0

        # have to account for exp blowing up to inf
        if exp(α[1]*z0 - α[1]*z) > 1e8
            ϕ = 0.0
        else
            ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
            ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
            ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
            ϕ = ϕ*exp(α[1]*z0 - α[1]*z)
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))
        end

    ## equation 12
    elseif z > z0 && z < l
        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

    #equation 13
    elseif z >= l

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

function manual_integral(ρ, z, μa, μsp, l)

    xs = 0:0.001:300
    a = map((s) -> fluence_spatial_freq(z, s, μa, μsp, l), xs)
    b = a.*xs
    c  = b.*[besselj0(ρ*x) for x in xs]
    return c
end

# compute equation 14 integrand
function compute_integral(ρ, z, s, μa, μsp, l)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

# numerically integrate equation 14
function fluence_steadystate(ρ, z, μa, μsp, l)

    ϕ = quadgk((s) -> compute_integral(ρ, z, s, μa, μsp, l), 0, Inf, rtol=1e-8)
    return ϕ[1]/2π

 end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using LightNLP
using LightNLP.DepParser
# using JLD2, FileIO

include("config.jl")

if CONFIG["training"]
    parser = DepParser.Decoder(CONFIG)
    # save("parser.jld2", "parser", parser)
else
    parser = load("parser.jld2", "parser")
    decode(parser, testfile)
end


text:
@everywhere using FWI
using Test
# using PyPlot

include("setupFWItests.jl")
# include("../src/plotResults.jl")

Dobs, = getData(m,pFor)
# setup inversion parameters
Wd    = 1./(abs(Dobs) + 0.1)
alpha = 1e-3
mref  = ones(nx*nz)
pInv  = InversionParam(alpha,gamma,mref,Mesh)

opti = (GNinv, NLCGinv, SDinv)

println("Test optimizers\n")
for k=1:length(opti)
    println("\ttesting $(string(opti[k]))")

    mc,his = opti[k](copy(mref),Dobs,Wd,pFor,pInv,maxIter=4,doPlot=false)
    mp,hisp = opti[k](copy(mref),Dobs,Wd,pForp,pInv,maxIter=4)

    @test norm(mc-mp)/norm(mc) < 1e-10
    @test norm(vec(his-hisp))/norm(vec(his)) < 1e-10
    @test size(his,1)==5
    @test all(his[5,:] .!= 0)
    @test all(diff(his[:,3]).<0) # descent method
    @test all(diff(his[:,4]).<0) # reduction in gradient

    println("\t\tOK!")
end


text:
module TestJost
using Test

# Checking Jost's diversities
using Diversity
using Diversity.ShortNames
using Diversity.Jost

qs = [0, 1, 2, 3, 4, 5, 6, Inf];

numspecies = 100;
numcommunities = 8;
communities = rand(numspecies, numcommunities);
communities /= sum(communities);
probs = reshape(mapslices(sum, communities, dims=2), (size(communities, 1)));
colweights = rand(numcommunities);
colweights /= sum(colweights);
allthesame = probs * colweights';

@testset "Jost" begin
    @test jostbeta(communities, 1)[!,:diversity] ≈
        (1.0 ./ metadiv(ρ̄(Metacommunity(communities)), 1)[!,:diversity])
    @test jostbeta(allthesame, qs)[!,:diversity] ≈ fill!(similar(qs), 1)
    @test jostbeta(Metacommunity(allthesame), 1)[!,:diversity] ==
        jostbeta(allthesame, 1)[!,:diversity]
    @test_throws ErrorException jostbeta(Metacommunity(allthesame,
        GeneralTypes(rand(numspecies, numspecies))), 1)
    ## Check Jost's alpha diversity works for all the same subcommunity
    @test jostalpha(allthesame, qs)[!,:diversity] ≈
          metadiv(ᾱ(Metacommunity(allthesame)), qs)[!,:diversity]
    @test jostalpha(Metacommunity(allthesame), 1)[!,:diversity] ==
        jostalpha(allthesame, 1)[!,:diversity]
    @test_throws ErrorException jostalpha(Metacommunity(allthesame,
        GeneralTypes(rand(numspecies, numspecies))), 1)

    ## And for all different subcommunities and any subcommunities with the same sizes
    weights = rand(numspecies);
    weights /= sum(weights);
    communitylist = rand(1:numcommunities, numspecies)
    distinct = zeros(Float64, (numspecies, numcommunities))
    for i in 1:numspecies
        distinct[i, communitylist[i]] = weights[i]
    end
    evendistinct = mapslices((x) -> x / (sum(x) * numcommunities), distinct, dims=1)

    @test jostalpha(evendistinct, qs)[!,:diversity] ≈
          metadiv(ᾱ(Metacommunity(evendistinct)), qs)[!,:diversity]

    # Now some even communities, should see that raw and normalised
    # diversities are the same
    smoothed = communities ./ mapslices(sum, communities, dims=1);
    smoothed /= numcommunities;
    @test jostalpha(smoothed, qs)[!,:diversity] ≈
          metadiv(ᾱ(Metacommunity(smoothed)), qs)[!,:diversity]
end

end


text:
using ProbNumODE
using ProbNumODE: remake_prob_with_jac, multitime, make_savable, compute_errors
using BSON: @save
using DifferentialEquations
using Plots
using ProgressMeter

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_lotkavoltera, prob_ode_rigidbody
lotkavolterra = remake(prob_ode_lotkavoltera, tspan=(0.0, 10.0))

failed_runs = []

problist = [
    (logistic_equation(), "logistic"),
    (lotkavolterra, "lotkavolterra"),
    (fitzhugh_nagumo_iip(), "fitzhughnagumo"),
    (prob_ode_rigidbody, "rigidbody"),
    (brusselator(), "brusselator"),
]

for (prob, probname) in problist
    @info "New Problem: $probname"
    prob = remake_prob_with_jac(prob)
    appxsol = DiffEqBase.has_analytic(prob.f) ? nothing :
        solve(remake(prob, u0=big.(prob.u0)), abstol=1e-18, reltol=1e-18)

    result_list = []

    algs = [EKF0(), EKF1()]
    qs = 1:5
    sigmas = [:fixedMLE, :fixedMLEMV, :dynamic, :dynamicMV,]
    steprules = [:standard]
    smooths = [false, true]

    abstols, reltols = Dict(
        "logistic" => (4:13, 1:10),
        "lotkavolterra" => (5:13, 2:10),
        "fitzhughnagumo" => (4:13, 1:10),
        "rigidbody" => (4:13, 1:10),
        "brusselator" => (4:13, 1:10),
    )[probname]
    abstols, reltols = 10.0 .^ -abstols, 10.0 .^ -reltols

    n_combinations = length(algs) * length(qs) * length(steprules) * length(sigmas) * length(smooths) * length(abstols)
    p = Progress(n_combinations, 1, "Working on $probname: ")

    for alg in algs,
        q in qs,
        sigmarule in sigmas,
        steprule in steprules,
        smooth in smooths,
        (abstol, reltol) in zip(abstols, reltols)

        next!(p, showvalues=[(:alg, alg), (:q, q), (:sigmarule, sigmarule),
                             (:smooth, smooth),
                             (:abstol, abstol), (:reltol, reltol)])

        if alg isa EKF1 && sigmarule in (:dynamicMV, :fixedMLEMV) continue end

        config = (
            alg=string(nameof(typeof(alg))),
            q=q,
            steprule=steprule,
            sigmarule=sigmarule,
            abstol=abstol,
            reltol=reltol,
            smooth=smooth,
        )

        @debug "Run started:" probname config

        try

            solve_call = () -> solve(prob, alg; config..., maxiters=1e6)
            sol = solve_call()
            if sol.retcode != :Success
                if sol.retcode == :MaxIters
                    # @info "MaxIters" sol.retcode probname config
                else
                    @warn "Run not successful!" sol.retcode probname config
                    push!(failed_runs, (sol.retcode, probname, config))
                end
                continue
            end

            runtime = multitime(solve_call; numruns=10)
            summaries = Dict(
                :runtime=>runtime,
                :nf => sol.destats.nf,
                :num_evals => sol.destats.nf + sol.destats.njacs,
                :naccept => sol.destats.naccept,
                :nreject => sol.destats.nreject,
            )
            savable_sol = make_savable(sol; appxsol=appxsol)
            final_vals = (
                x=savable_sol.x[end],
                t=savable_sol.t[end],
                pu=savable_sol.pu[end],
                u_ref=savable_sol.u_ref[end],
            )

            errors = compute_errors(savable_sol)
            run = (
                problem=probname,
                config=config,
                # solution=savable_sol,
                final_vals=final_vals,
                summaries=merge(summaries, errors),
                destats=sol.destats,
            )

            push!(result_list, run)
            @debug "Run finished:" run
        catch err
            if err isa InterruptException
                throw(err)
            end
            @warn "Encountered error!" err probname config
            push!(failed_runs, (err, probname, config))
            # throw(err)
        end
    end

    varname = "results"
    @save "data/$(probname)_fixedMLEMV.bson" result_list
end

@info "All unsuccessful runs again:" failed_runs

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
using SafeTestsets
using Metatheory
using Test

function test(file::String)
  @info file
  @eval @time @safetestset $file begin
    include(joinpath(@__DIR__, $file))
  end
end

const TEST_FILES = ["reductions.jl", "EGraphs/egraphs.jl", "EGraphs/ematch.jl", "EGraphs/analysis.jl"]
const INTEGRATION_TEST_FILES = map(
  x -> joinpath("integration", x),
  [
    "custom_types.jl",
    "fibonacci.jl",
    "kb_benchmark.jl",
    "logic.jl",
    "mu.jl",
    "stream_fusion.jl",
    "taylor.jl",
    "while_interpreter.jl",
    "while_superinterpreter.jl",
  ],
)

@timev begin
  @timev map(test, TEST_FILES)
  @timev map(test, INTEGRATION_TEST_FILES)
end

# TODO: use Aqua.jl
# exported consistency test
for m in [Metatheory, Metatheory.EGraphs, Metatheory.EGraphs.Schedulers]
  for i in propertynames(m)
    xxx = getproperty(m, i)
  end
end


text:
# Grigoriadis Khachiyan Matrix Games.
#
# Example run: gk(10000, [0.1])
#
# Code from Dilys Thomas <dilys@cs.stanford.edu>

##### generate skew symmetric matrix  #######
function myunifskew(n)
    A = zeros(n, n)

    #print("A[i,j] initialized with zeros \n");

    for i=1:n
    	for j=1:i-1
            temp=rand()
	    #print(temp)
	    if (temp < 0.5)
		temp = rand()
         	A[i,j]= temp
         	A[j,i]= -A[i,j]
		#print("welcome");
            else
		temp = rand()
         	A[j,i]= temp
         	A[i,j]= -A[j,i]
		#print("welcome");
            end

    	end
    	if rem(i,1000) == 0
	    #print(i)
	    #print("\n")
    	end
    end

    return A
end

############ GK Algorithm starts ##################
#@profile begin
function gk(n, myeps)

    A = myunifskew(n)

    g = length(myeps)
    iteration = zeros(g)
    times = zeros(g)

    for KK=1:g

	eps = myeps[KK]

	e = ones(n)
	X = zeros(n)
	U = zeros(n)
	p = e./n

	t = 0
	#tm=U/t
	stop = 0
	iter = 0
	epse = eps .* e

	csum = zeros(n)

#	tic()
	while(stop != 1)
     	    t=t+1
     	    iter=t
     	    #iteration number

     	    if rem(iter, 100) == 0
         	#print(iter)
		#print("\n")
     	    end

	    for i=1:n
        	csum[i] = sum(p[1:i])
    	    end

	    marker=rand()
	    k=1

  	    for i=2:n
		if csum[i-1] <= marker && marker <= csum[i]
                    k=i
    		    break
            	end
            end

    	    X[k] += 1
    	    for i=1:n
        	U[i] += A[i,k]
    	    end

    	    s = sum(p[1:n] .* exp((eps/2)*A[1:n,k]))
    	    for i=1:n
       		p[i]=(p[i]*exp((eps/2)*A[i,k])) / s
    	    end

	    u = U ./ t

	    True=0
	    for i=1:n
    		if u[i] <= epse[i]
		    True = True+1
		end
    	    end

	    if True == n
		stop=1
	    end

	end

	times[KK] = 0#toc()
	iteration[KK] = iter

	x = X/t
	etx=sum(x)

	AX=A*X
	Ax = A*x
	error=abs(AX)-abs(U)
	#print(Ax)

	Axepse=0
	for i=1:n
	    if Ax[i]<=epse[i]
		Axepse = Axepse+1
	    end
    	end

	if Axepse==n
	    #print(" \n Ax <= eps*e  \n")
	end

	#if A*x <= eps*e
    	#    print(" Ax <= eps*e \n ")
	#end

	errorlmt = 0
	for i=1:n
	    if error[i]<1e-8
    		errorlmt = errorlmt+1
	    end
	end

        if errorlmt ==n
    	    #print("Assertion condition is satisfied i.e. AX-U<10^-8 \n")
	else
    	    print("Error:  AX-U<10^-8 not satisfied \n")
	end

	#print("welcome")

	# print("Time for \n")
	# print(eps)
	# print("\n")
	# print("is  \n")
	# print(times[KK])
	# print("\n")
	# print("Number of iteration is \n")
	# print(iteration[KK])
	# print("\n")
    end

    # print("Epsilon vector is \n")
    # print(myeps)
    # print("\n")
    # print("time vector is \n ")
    # print(times)
    # print("\n")
    # print("Iteration Vector is \n")
    # print(iteration)
    # print("\n \n")

    #out = [myeps, time, iteration]
    #print("Epsilon-Time-Iteration tradeoff \n")
    #print(out)
    #print("\n \n")

end
#end # @profile begin

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function calc_phi_deriv(x,xi,n,nn,c,vars,var_seq)
    phi_der=zeros(nn,n)
    for i=1:nn
        for j=1:n
            dr=0
            for k=1:vars
                dr = dr + (xi[i,k] - x[j,k])^2
            end
            r = (dr)

            # 2 r^4
            # phi_der[i,j] = -(r^3)*(8*(xi[i,var_seq]-x[j,var_seq]))/4

            # 1 Gaussian
            phi_der[i,j]=(-2*xi[i,var_seq] + 2*x[j,var_seq])*exp(-r/c)/c

        end
    end
    return phi_der
end

# phi[i,j]=exp(-dr^2/c^2)*(-2*dr/c)
# phi(i,j) = exp(-r^2)/2 + (r*pi^(1/2)*erf(r))/2;
# phi[i,j] = 2*(1./(1+exp.(-r)))-1;
# phi[i,j] = (-2*r/c[i,j])*exp(-r^2/c[i,j])
# phi[i,j]=tanh.(r/c[i,j])
# phi_der[i,j] = -(r^2)*(2*(xi[i,var_seq]-x[j,var_seq]))
# phi_der[i,j] = (-2*xi[i,var_seq] + 2*x[j,var_seq])*exp(-r/c)/c
# phi_der[i,j] = (c*sqrt(pi)*erf(sqrt(r)/c))/2


text:
import HTTP
import JSON

import SQLite

include("constants.jl")
include("config.jl")
include("database.jl")

function request(url::String, requestType::RequestType, token::String, database::SQLite.DB)
    headers = [
        "Authorization" => "Bearer $token",
        "Accept" => "application/json",
        "Host" => HOST,
        "User-Agent" => USERAGENT,
        "Time-Zone" => TIMEZONE,
    ]

    response = HTTP.get(url, headers)
    content = JSON.parse(String(response.body))

    if requestType == TypeUserinfo
        user = User(get(content, "id", 0), get(content, "login", ""))
        SQLiteCreateUser(database, user)
    elseif requestType == TypeFollowing || requestType == TypeFollowers
        for con in content
            user = User(get(con, "id", 0), get(con, "login", ""))
            SQLiteCreateUser(database, user)
        end
        link = headerLink(response.headers)
        if link != ""
            sleep(3)
            request(link, requestType, token, database)
        end
    end
end

function headerLink(headers::Array{Pair{SubString{String},SubString{String}},1})::String
    for head in headers
        if head.first == "Link"
            for str in split(head.second, ", ")
                s = split(str, "; ")
                if s[2] == "rel=\"next\""
                    return replace(s[1], r"[<>]" => "")
                end
            end
        end
    end
    return ""
end

max_stars_repo_name: tosone/julia-spider

max_issues_repo_name: tosone/julia-spider

max_forks_repo_name: tosone/julia-spider


text:
# * Solver trace columns
# ** Energy

mutable struct EnergyColumn{T<:AbstractFloat} <: TraceColumn
    E::T
    fmt::FormatExpr
    header::String
    eV::Bool
end

function EnergyColumn(E::T, eV::Bool=true, header="Energy") where T
    fmt = if eV
        FormatExpr("{1:+10.5f} Ha = {2:+10.5f} eV")
    else
        FormatExpr("{1:+10.5f} Ha")
    end
    EnergyColumn(E, fmt,
                 rpad(header, 1+length(eV ? format(fmt, E, E) : format(fmt, E))), eV)
end

(e::EnergyColumn)(::Integer) =
    e.eV ? (e.E,27.211e.E) : (e.E,)

# ** Virial theorem

mutable struct VirialColumn{T<:AbstractFloat} <: TraceColumn
    V::T
    fmt::FormatExpr
    lc::LinearColorant{T}
    header::String
end

function VirialColumn(V::T) where T
    fmt = FormatExpr("{1:.5f} ({2:s}{3:.2f}×10{4:3s}$(crayon"reset"))")
    lc = LinearColorant(one(T),zero(T),SolverTraces.red_green_scale())
    VirialColumn(V, fmt, lc,
                 rpad("⟨V̂⟩/⟨T̂⟩ ⟨V̂⟩/⟨T̂⟩ + 2", 1+length(format(fmt, V, "", base_exp(abs(V + 2.0))...))))
end

function (v::VirialColumn{T})(::Integer) where T
    Vb,Ve = base_exp(abs(v.V + 2.0))
    (v.V, v.lc(100abs.(log10(2.0)-log10(abs(v.V)))/log10(2.0)),
     Vb,to_superscript(Ve))
end

# ** Relaxation

mutable struct RelaxationColumn{T<:AbstractFloat} <: TraceColumn
    ω::T
    fmt::FormatExpr
    header::String
end

function RelaxationColumn(ω::T) where T
    header = "1-ω"
    fmt = FormatExpr("{1:.2f}×10{2:3s}")
    RelaxationColumn(ω, fmt,
                     rpad(header, 1+length(format(fmt, base_exp(abs(1.0 - ω))...))))
end

function (r::RelaxationColumn{T})(::Integer) where T
    ωb,ωe = base_exp(abs(1.0 - r.ω))
    (ωb,to_superscript(ωe))
end

# ** Flags

mutable struct FlagsColumn <: TraceColumn
    val::Vector{String}
    fmt::FormatExpr
    header::String
end

function FlagsColumn()
    header = "Flags"
    FlagsColumn(String[], FormatExpr("{1:5s}"), "Flags")
end

(flags::FlagsColumn)(::Integer) = (join(flags.val, ""),)

# * Setup solver trace
function setup_solver_trace(verbosity, max_iter, tol, ω, num_printouts;
                            tol_header="Tolerance")
    if verbosity > 1
        tolerance = Tolerance(tol, tol_header, print_target=false)

        columns = (CurrentStep(max_iter,
                               lc=LinearColorant(max_iter,1,SolverTraces.red_green_scale()),
                               header="Iteration"),
                   tolerance)

        relaxation = if !iszero(ω)
            columns = (columns..., RelaxationColumn(ω))
            last(columns)
        else
            nothing
        end

        eng = EnergyColumn(0.0),EnergyColumn(0.0,false,"⟨T̂⟩"),EnergyColumn(0.0,false, "⟨V̂⟩")
        virial = VirialColumn(-2.0)
        flags = FlagsColumn()

        columns = (columns..., eng..., virial, flags)

        trace = SolverTrace(max_iter,
                            columns...;
                            progress_meter=false,
                            num_printouts=num_printouts)

        trace,tolerance,relaxation,eng,virial,flags
    else
        nothing,nothing,nothing,nothing,nothing,nothing
    end
end


text:
push!(LOAD_PATH, "./")

using CASIM
using PyPlot
using PyCall
@pyimport matplotlib.animation as animation

rc("font", family="Times New Roman")
rc("font", size=16)
rc("text", usetex=true)

# Settings:
sim_horizon = 1000
init_num_ac = 30
TAKEOFF_RATE = 50
spawn_rate = TAKEOFF_RATE * 100 / 3600
airspace_controller = Constant_Spawn_Rate_Controller(init_num_ac, spawn_rate)

# cas = Corrected_VICAS()
# policy.env.correction_weight = 0.075
# policy_name = "CorrectedClosest"
# global correction = true
# global factor = 1.0

# cas = VICAS()
# global correction = false
# policy_name = "VICASMulti"
# global factor = 1.0

# cas = No_CAS()
# global correction = false
# policy_name = "NOCAS"
# global factor = 1.0

cas = VICAS()
global correction = false
policy_name = "VICASClosest"
global factor = 1.0

cas_ctrl = Uniform_CAS(cas)
# sensor = Multi_Threat_Tracker()
sensor = Nearest_Intruder_Tracker()
sensor_ctrl = Uniform_Sensor(sensor)

fig = figure("SimAnim", figsize=(15, 8))

if correction_type == :closest
    if airspace_controller.spawn_rate > 4900./3600
        policy.env.correction = true
        println("Correction turned on")
    else
        policy.env.correction = false
        println("Correction turned off")
    end
elseif correction_type == :sector
    if airspace_controller.spawn_rate > 7900./3600
        policy.env.correction = true
        println("Correction turned on")
    else
        policy.env.correction = false
        println("Correction turned off")
    end
end

T = time()
global stats = Stats(T)
global t = 0
global All_AC = airspace_initialization(stats, airspace_controller.ac_number, cas_ctrl, sensor_ctrl, t)

function ft2m(ft::Float64)
    return ft * 0.3048
end

function animate(i)
    global All_AC
    global t
    global stats
    global correction
    global factor
    # Plot:
    clf()
    ax1 = subplot2grid((2, 2), (0, 0), rowspan=3)
    for ac in All_AC
        ax1 = scatter(ft2m(ac.dyn.x), ft2m(ac.dyn.y), marker="o", color="blue", s=6)
        ax1 = scatter(ft2m(ac.x_dest), ft2m(ac.y_dest), marker=",", color="magenta", s=3)
        ax1 = plot([ft2m(ac.dyn.x), ft2m(ac.x_dest)], [ft2m(ac.dyn.y), ft2m(ac.y_dest)], linestyle="--", color="black", linewidth=0.15)
    end
    ax1 = xlabel("\$x\$ (m)")
    ax1 = ylabel("\$y\$ (m)")
    ax1 = title(policy_name * ", Take-off Rate = " * string(TAKEOFF_RATE) * " flight/km\$^2\$-hour", fontsize=18)
    ax = gca()
    ax[:set_xlim]([0, 10000])
    ax[:set_ylim]([0, 10000])
    ax[:set_aspect]("equal")

    # Step:
    if t % 100 == 0
        println(string(t) * " : " * string(length(All_AC)))
    end
    update_all_dynamics(All_AC, correction=correction)
    t += 1
    run_stats(All_AC, stats, t)
    All_AC = airspace_control(stats, All_AC, airspace_controller, cas_ctrl, sensor_ctrl, t)

    ax2 = subplot2grid((2, 2), (0, 1))
    ax2 = plot(1 : length(stats.num_ac), stats.num_ac, color="red")
    ax2 = xlabel("Time (sec)")
    ax2 = ylabel("Number of Aircraft")
    ax = gca()
    ax[:set_xlim]([0, 1000])
    ax[:set_ylim]([0, 1000])

    ax3 = subplot2grid((2, 2), (1, 1))
    ax3 = plot(1 : length(stats.num_NMAC_t), stats.num_NMAC_t, color="red")
    println(stats.num_NMAC_t[end])
    ax3 = xlabel("Time (sec)")
    ax3 = ylabel("NMAC/sec")
    ax = gca()
    ax[:set_xlim]([0, 1000])
    ax[:set_ylim]([0, 25])

    frame = (ax1, ax2, ax3)

    return frame
end

sim_anim = animation.FuncAnimation(fig, animate, frames=sim_horizon, interval=200)
sim_anim[:save]("test" * policy_name * "_" * string(Int(spawn_rate * 3600)) * ".mp4", fps=24, extra_args=["-vcodec", "libx264"])

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
@testset "Reconstruct" begin
    @test prettify(reconstruct(deconstruct(f))) == prettify(f)
    @test prettify(reconstruct(deconstruct(doc_f))) == prettify(doc_f)
    @test prettify(reconstruct(deconstruct(e2))) == prettify(e2)
    @test prettify(reconstruct(deconstruct(s))) == prettify(s)
    @test prettify(reconstruct(deconstruct(doc_s))) == prettify(doc_s)
end