text:
struct InceptionBlock
  path_1
  path_2
  path_3
  path_4
end

@functor InceptionBlock

function InceptionBlock(in_chs, chs_1x1, chs_3x3_reduce, chs_3x3, chs_5x5_reduce, chs_5x5, pool_proj)
  path_1 = Conv((1, 1), in_chs=>chs_1x1, relu)

  path_2 = (Conv((1, 1), in_chs=>chs_3x3_reduce, relu),
            Conv((3, 3), chs_3x3_reduce=>chs_3x3, relu, pad = (1, 1)))

  path_3 = (Conv((1, 1), in_chs=>chs_5x5_reduce, relu),
            Conv((5, 5), chs_5x5_reduce=>chs_5x5, relu, pad = (2, 2)))

  path_4 = (MaxPool((3,3), stride = (1, 1), pad = (1, 1)),
            Conv((1, 1), in_chs=>pool_proj, relu))

  InceptionBlock(path_1, path_2, path_3, path_4)
end

function (m::InceptionBlock)(x)
  cat(m.path_1(x), m.path_2[2](m.path_2[1](x)), m.path_3[2](m.path_3[1](x)), m.path_4[2](m.path_4[1](x)), dims = 3)
end

_googlenet() = Chain(Conv((7, 7), 3=>64, stride = (2, 2), relu, pad = (3, 3)),
      MaxPool((3, 3), stride = (2, 2), pad = (1, 1)),
      Conv((1, 1), 64=>64, relu),
      Conv((3, 3), 64=>192, relu, pad = (1, 1)),
      MaxPool((3, 3), stride = (2, 2), pad = (1, 1)),
      InceptionBlock(192, 64, 96, 128, 16, 32, 32),
      InceptionBlock(256, 128, 128, 192, 32, 96, 64),
      MaxPool((3, 3), stride = (2, 2), pad = (1, 1)),
      InceptionBlock(480, 192, 96, 208, 16, 48, 64),
      InceptionBlock(512, 160, 112, 224, 24, 64, 64),
      InceptionBlock(512, 128, 128, 256, 24, 64, 64),
      InceptionBlock(512, 112, 144, 288, 32, 64, 64),
      InceptionBlock(528, 256, 160, 320, 32, 128, 128),
      MaxPool((3, 3), stride = (2, 2), pad = (1, 1)),
      InceptionBlock(832, 256, 160, 320, 32, 128, 128),
      InceptionBlock(832, 384, 192, 384, 48, 128, 128),
      MeanPool((7, 7), stride = (1, 1), pad = (0, 0)),
      x -> reshape(x, :, size(x, 4)),
      Dropout(0.4),
      Dense(1024, 1000), softmax)

function googlenet_layers()
  weight = Metalhead.weights("googlenet.bson")
  weights = Dict{Any, Any}()
  for ele in keys(weight)
    weights[string(ele)] = convert(Array{Float64, N} where N, weight[ele])
  end
  ls = _googlenet()
  ls[1].weight .= weights["conv1/7x7_s2_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[1].bias .= weights["conv1/7x7_s2_b_0"]
  ls[3].weight .= weights["conv2/3x3_reduce_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[3].bias .= weights["conv2/3x3_reduce_b_0"]
  ls[4].weight .= weights["conv2/3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[4].bias .= weights["conv2/3x3_b_0"]
  for (a, b) in [(6, "3a"), (7, "3b"), (9, "4a"), (10, "4b"), (11, "4c"), (12, "4d"), (13, "4e"), (15, "5a"), (16, "5b")]
    ls[a].path_1.weight .= weights["inception_$b/1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_1.bias .= weights["inception_$b/1x1_b_0"]
    ls[a].path_2[1].weight .= weights["inception_$b/3x3_reduce_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_2[1].bias .= weights["inception_$b/3x3_reduce_b_0"]
    ls[a].path_2[2].weight .= weights["inception_$b/3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_2[2].bias .= weights["inception_$b/3x3_b_0"]
    ls[a].path_3[1].weight .= weights["inception_$b/5x5_reduce_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_3[1].bias .= weights["inception_$b/5x5_reduce_b_0"]
    ls[a].path_3[2].weight .= weights["inception_$b/5x5_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_3[2].bias .= weights["inception_$b/5x5_b_0"]
    ls[a].path_4[2].weight .= weights["inception_$b/pool_proj_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]; ls[a].path_4[2].bias .= weights["inception_$b/pool_proj_b_0"]
  end
  ls[20].W .= transpose(weights["loss3/classifier_w_0"]); ls[20].b .= weights["loss3/classifier_b_0"]
  return ls
end

struct GoogleNet <: ClassificationModel{ImageNet.ImageNet1k}
  layers::Chain
end

GoogleNet() = GoogleNet(googlenet_layers())

Base.show(io::IO, ::GoogleNet) = print(io, "GoogleNet()")

@functor GoogleNet

(m::GoogleNet)(x) = m.layers(x)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenBLAS32_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenBLAS32")
JLLWrappers.@generate_main_file("OpenBLAS32", UUID("656ef2d0-ae68-5445-9ca0-591084a874a2"))
end  # module OpenBLAS32_jll

ext: jl


text:
using BayesianDataFusion
using Distributed

if nprocs() <= 1
  addprocs(8)
  @everywhere using BayesianDataFusion
end

W = sprand(1_500_000, 1000, 0.01)
#W = sprand(150_000, 1000, 0.01)
rd = RelationData(W, class_cut = 0.5)
assignToTest!(rd.relations[1], 50)

function set_threads(n)
   if VERSION >= v"0.5.0"
     BLAS.set_num_threads(n)
   else
     blas_set_num_threads(n)
   end
end

############ num_latent = 10 ############
## single threaded
set_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

set_threads(8)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

set_threads(16)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

## multi-threaded (1)
set_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2], latent_blas_threads=1);

## multi-threaded (2)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=1);
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=2);

############ num_latent = 30 ############
set_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

set_threads(2)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

set_threads(4)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

## multi-threaded (1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2], latent_blas_threads=1);

## multi-threaded (2)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=1);
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=2);

## multi-threaded (4)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3,4,5], latent_blas_threads=1);

## multi-threaded (8)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2:9;], latent_blas_threads=1);


text:
"""

Definitions for Twists and Wrenches for use in 3D rigid body dynamics

    $(EXPORTS)

    $(SIGNATURES)

    $(METHODLIST)

"""
module Screws
using LinearAlgebra
using StaticArrays
using Rotations
using DocStringExtensions
using Cartesian

using Base: @propagate_inbounds

import Base.convert, Base.*
import LinearAlgebra.cross, LinearAlgebra.dot
import Cartesian.outer

export
    Twist,
    Wrench,
    Matrix6,
    SpatialInertia,
    SpatialMobility,
    spi,
    spm,
    dot

const Matrix6 = SMatrix{2,2, Matrix3}

for SpatialVector in (:Twist, :Wrench)
    @eval begin
        # Data Structures
        struct $SpatialVector
            linear::Vector3
            angular::Vector3

            @inline function $SpatialVector(linear::AbstractVector, angular::AbstractVector) where T
                @boundscheck( size(linear,1)==3 || throw(DimensionMismatch()))
                @boundscheck( size(angular,1)==3 || throw(DimensionMismatch()))
                return new(linear, angular)
            end
        end

        # Basics
        Base.eltype(::Type{$SpatialVector})  = Float64
        """
        $(SIGNATURES)
        Linear (first) part of the screw
        """
        linear(v::$SpatialVector) = v.linear
        """
        $(SIGNATURES)
        Angular (second) part of the screw
        """
        angular(v::$SpatialVector) = v.angular

        # Construct/convert given another spatial vector
        $SpatialVector(v::$SpatialVector) = v
        Base.convert(::Type{V}, v::$SpatialVector) where {V<:$SpatialVector} = V(v)

        # Construct/convert to SVector
        StaticArrays.SVector{6}(v::$SpatialVector) = convert(SVector{6}, [linear(v); angular(v)])
        StaticArrays.SVector(v::$SpatialVector) = SVector{6}(v)
        StaticArrays.SArray(v::$SpatialVector) = SVector(v)
        Base.convert(::Type{SA}, v::$SpatialVector) where {SA <: SArray} = SA(v)

        function Base.Array(v::$SpatialVector)
            Base.depwarn("Array(v) has been deprecated. Please use SVector(v) instead.", :Array)
            SVector(v)
        end
    end
end

"""
$(TYPEDEF)

A Twist represents a motion screw. A screw groups together a linear and an angular vector
into a single 6-vector. For twists the linear part transforms with location, and the angular
part remains constant.

"""
Twist

"""
$(TYPEDEF)

A Wrench represents a force screw. A screw groups together a linear and an angular vector
into a single 6-vector. For wrenches  the linear remains constant, while  the angular
part transforms with location.

"""
Wrench

"""
$(SIGNATURES)

Construct a [`Twist`](@ref) from a value vector, a position vector and a pitch

```math
\\mathrm{Twist}(\\boldsymbol{v}, \\boldsymbol{r}, h) =
\\pmatrix{ \\boldsymbol{r}\\times \\boldsymbol{v} + h\\,\\boldsymbol{v} \\\\ \\boldsymbol{v} }
```
"""
Twist(value::Vector3, position::Vector3, pitch::Float64 ) = Twist( cross(position, value) + pitch*value, value)
#Twist(value::Vector3, position::Vector3) = Twist(value, position, 0.0)
"""
$(SIGNATURES)

Construct a pure Twist
```math
\\mathrm{Twist}(\\boldsymbol{v}) =
\\pmatrix{ \\boldsymbol{v} \\\\ \\boldsymbol{0} }
```
"""
Twist(pure::Vector3 ) = Twist(pure, ô)
Twist() = Twist(ô)

function Base.:/(v::Twist, g::AbstractFloat)::Twist
    Wrench(v.linear/g, v.angular/g)
end

"""
$(SIGNATURES)
Construct a [`Wrench`](@ref) from a value vector, a position vector and a pitch

```math
\\mathrm{Wrench}(\\boldsymbol{v}, \\boldsymbol{r}, h) =
\\pmatrix{ \\boldsymbol{v} \\\\ \\boldsymbol{r}\\times \\boldsymbol{v} + h\\,\\boldsymbol{v} }
```
"""
Wrench(value::Vector3, position::Vector3, pitch::Float64 ) = Wrench(value, cross(position, value) + pitch*value)
#Wrench(value::Vector3, position::Vector3) = Wrench(value, position, 0.0)

"""
$(SIGNATURES)

Construct a pure Wrench

```math
\\mathrm{Wrench}(\\boldsymbol{\tau}) =
\\pmatrix{ \\boldsymbol{0} \\\\ \\boldsymbol{\tau}  }
```
"""
Wrench(pure::Vector3) = Wrench(ô, pure)
Wrench() = Wrench(ô)

value(v::Twist) = v.angular
axis(v::Twist) = v.linear
value(v::Wrench) = v.linear
axis(v::Wrench) = v.angular

function Base.:/(f::Wrench, g::AbstractFloat)::Wrench
    Wrench(f.linear/g, f.angular/g)
end

"""
$(SIGNATURES)

Dot products between Twists and Wrenches, used in power calculations.
"""
dot(v::Twist, f::Wrench) = dot(v.linear, f.linear) + dot(v.angular, f.angular)
dot(f::Wrench, v::Twist) = dot(v,f)

function outer(a::Twist, b::Twist)::Matrix6
    Matrix6(a.linear .* b.linear', a.linear .* b.angular',
    a.angular .* b.linear', a.angular .* b.angular')
end
function outer(a::Wrench, b::Wrench)::Matrix6
    Matrix6(a.linear .* b.linear', a.linear .* b.angular',
    a.angular .* b.linear', a.angular .* b.angular')
end

"""
$(SIGNATURES)

Cross products between Twist and Wrencehs, used in derivatives of rotating frames.
"""
cross(v::Twist, w::Twist) = Twist( cross(v.angular, w.linear) + cross(v.linear, w.angular), cross(v.angular, w.angular) )
cross(v::Twist, f::Wrench) = Wrench( cross(v.angular, w.linear), cross(v.linear, w.linear) + cross(v.angular, w.angular) )
cross(f::Wrench, v::Twist) = Twist( cross(f.linear, v.linear) + cross(f.angular, v.angular) + cross(f.linear, v.angular) )
cross(f::Wrench, g::Wrench) = Wrench( cross(f.linear, g.linear), cross(f.angular, g.linear) + cross(f.linear, g.angular) )

for SpatialMatrix in (:SpatialInertia, :SpatialMobility)
    @eval begin
        # Data Structures
        struct $SpatialMatrix
            ll::Matrix3
            la::Matrix3
            al::Matrix3
            aa::Matrix3

            @inline function $SpatialMatrix(ll::AbstractMatrix, la::AbstractMatrix, al::AbstractMatrix, aa::AbstractMatrix) where {T}
                @boundscheck( (size(ll,1)==3 && size(ll,2)==3) || throw(DimensionMismatch("$SpatialMatrix needs $ll to be a 3×3 matrix.")))
                @boundscheck( (size(la,1)==3 && size(la,2)==3) || throw(DimensionMismatch("$SpatialMatrix needs $la to be a 3×3 matrix.")))
                @boundscheck( (size(al,1)==3 && size(al,2)==3) || throw(DimensionMismatch("$SpatialMatrix needs $al to be a 3×3 matrix.")))
                @boundscheck( (size(aa,1)==3 && size(aa,2)==3) || throw(DimensionMismatch("$SpatialMatrix needs $aa to be a 3×3 matrix.")))
                # if !(size(ll,1) == size(ll,2) == 3) || !(size(la,1) == size(la,2) == 3) || !(size(al,1) == size(al,2) == 3) || !(size(aa,1) == size(aa,2) == 3)
                #     throw(DimensionMismatch("$SpatialMatrix needs four matrices of size 3×3"))
                # end
                return new(ll,la,al,aa)
            end
        end
        # Basics
        Base.eltype(::Type{$SpatialMatrix})  = Float64
        """
        $(SIGNATURES)
        LinearLinear (top-left) part of the matrix
        """
        linearlinear(m::$SpatialMatrix) = m.ll
        """
        $(SIGNATURES)
        LinearAngular (top-right) part of the matrix
        """
        linearangular(m::$SpatialMatrix) = m.la
        """
        $(SIGNATURES)
        AngularLinear (bottome-left) part of the matrix
        """
        angularlinear(m::$SpatialMatrix) = m.al
        """
        $(SIGNATURES)
        AngularAngular (bottome-right) part of the matrix
        """
        angularangular(m::$SpatialMatrix) = m.aa

        # Construct/convert given another spatial vector
        $SpatialMatrix(m::$SpatialMatrix) = m
        function $SpatialMatrix(m::Matrix6)
            $SpatialMatrix(m[1,1], m[1,2], m[2,1], m[2,2])
        end
        Base.convert(::Type{MType}, m::$SpatialMatrix) where {MType<:$SpatialMatrix} = MType(m)

        # Construct/convert to SVector
        StaticArrays.SMatrix{6,6}(m::$SpatialMatrix) = convert(SMatrix{6,6},
            [m.ll m.la; m.al m.aa])
        StaticArrays.SMatrix(m::$SpatialMatrix) = SMatrix{6,6}(m)
        StaticArrays.SArray(m::$SpatialMatrix) = SMatrix(m)
        Base.convert(::Type{SA}, m::$SpatialMatrix) where {SA <: SArray} = SA(m)

        function Base.Array(m::$SpatialMatrix)
            Base.depwarn("Array(m) has been deprecated. Please use SMatrix(m) instead.", :Array)
            SMatrix(m)
        end
    end #eval
end #for

"""
$(TYPEDEF)

"""
SpatialInertia

"""
$(TYPEDEF)

"""
SpatialMobility

spi(m::AbstractFloat, Ic::Matrix3, cg::Vector3) = SpatialInertia( m*Î, -m*cross(cg), m*cross(cg), Ic-m*cross2(cg))
spm(m::AbstractFloat, Ic::Matrix3, cg::Vector3) = SpatialMobility( (1/m)*Î-cross(cg)*inv(Ic)*cross(cg), cross(cg)*inv(Ic), -inv(Ic)*cross(cg), inv(Ic))

function Base.:*(I::SpatialInertia, v::Twist)::Wrench
    Wrench(I.ll*v.linear + I.la*v.angular, I.al*v.linear + I.aa*v.angular)
end

function Base.:*(M::SpatialMobility, f::Wrench)::Twist
    Twist(M.ll*f.linear + M.la*f.angular, M.al*f.linear + M.aa*f.angular)
end

end # module


text:
## still needed now that we can support formula?
function bootstrap(data, response, statistic::Function, model, sampling::BootstrapSampling; transform::Function = identity)
    f0 = fit(model, data, response)
    t0 = statistic(transform(f0))
    #r0 = predict(f0) - response
    r0 = resid(f0, response) # -> resid
    t1 = zeros(sampling.nrun)
    r1 = copy(r0)
    for i in 1:sampling.nrun
        sample!(r0, r1)
        f1 = fit(model, data, response + r1) # -> combine
        t1[i] = statistic(transform(f1))
    end
    return t0, t1
end


text:
module TestIteration

  using GBIF
  using Test

  set = occurrences()

  @test iterate(set) == (set[1], 2)
  @test iterate(set, 2) == (set[2], 3)

  # Example that has some issues
  plotor = taxon("Procyon lotor")
  plotor_occ = occurrences(plotor)
  next!(plotor_occ)
  for o in plotor_occ
      @test typeof(o) <: GBIFRecord
      @test o.taxon.species.second == plotor.species.second
  end
  function is_in_us(o::GBIFRecord)
      o.countryCode == "US"
  end
  filter!(is_in_us, plotor_occ)
  for o in plotor_occ
      @test o.countryCode == "US"
  end

end


text:
"""
    _compute_endpoints(a::Arf, b::Arf, x::Arb)

Given `a`, `b` and `x` with `x = Arb((a, b))` compute
```
c = a - midpoint(x)
d = b - midpoint(x)
```
rounded towards zero. It returns `c, c_exact, d, d_exact` where `c`
and `d` are as above and `c_exact` is true if `c` was computed exactly
and false if rounding was performed, similarly for `d_exact`.

This is used in [`extrema_series`](@ref), [`minimum_series`](@ref) and
[`maximum_series`](@ref) to compute the interval on which the extrema
of the polynomial should be enclosed.

If `c_exact` is false then an enclosure of the endpoint can be
computed with [`_enclose_inexact_endpoint`](@ref).
"""
function _compute_endpoints(a::Arf, b::Arf, x::Arb)
    c, d = zero(a), zero(b)
    xmid = Arblib.midref(x)
    c_exact = iszero(Arblib.sub!(c, a, xmid, rnd = Arblib.ArbRoundToZero))
    d_exact = iszero(Arblib.sub!(d, b, xmid, rnd = Arblib.ArbRoundToZero))
    return c, c_exact, d, d_exact
end

"""
    _enclose_inexact_endpoint(endpoint::Arf)

Compute a ball centered at `endpoint` with a radius of 1 ulp of
`endpoint`.

This is used to get an enclosure of the endpoint in case
[`_compute_endpoints`](@ref) gives an inexact result. The semantics of
Arb guarantees that the error of the endpoint is at most 1 ulp when
using directed rounding and a ball with `endpoint` as midpoint and a
radius set to the value of this 1 ulp will therefore always enclose
the exact result.
"""
function _enclose_inexact_endpoint(endpoint::Arf)
    res = Arb(endpoint)
    Arblib.set_ulp!(Arblib.radref(res), endpoint, prec = precision(endpoint))
    return res
end

"""
    extrema_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute both the minimum and maximum of the function `f` on the
interval `[a, b]`. It returns a three tuple `(fmin, fmax, fmid)` where
`fmin` and `fmax` are enclosures of the minimum and maximum
respectively and `fmid` is an enclosure of `f` evaluated on the
midpoint of the interval.

The extrema are computed by computing the extrema of the Taylor series
of `f` on the midpoint of the interval and then adding an error term.
The degree of the Taylor series used is set with the `degree`
argument.

The reason that `fmid` is returned is that as part of computing the
Taylor series we get this value for free. In some cases the Taylor
series for at the midpoint is however not computed, in this case
`fmid` will just be `NaN` instead.

If `abs_value = true` then compute the extrema of `abs(f(x))` on the
interval `[a, b]`. For the computation of the maximum this is mostly
the same, only difference is that we have to compute the maximum of
the absolute value of the Taylor series. For the minimum we have to
take into account that the polynomial might cross zero, in which case
the minimum is zero.

If `verbose = true` then output information about the process.

This method is mainly intended for internal use in the function
[`extrema_enclosure`](@ref) but can be used directly as well.

TODO: Handle crossing zero when computing absolute minimum better.
Currently we always add the rest term in the end.
"""
function extrema_series(
    f,
    a::Arf,
    b::Arf;
    degree::Integer = 8,
    abs_value = false,
    verbose = false,
)
    check_interval(a, b)

    maybe_abs = abs_value ? abs : identity

    if a == b
        fa = maybe_abs(f(Arb(a)))
        return fa, fa, fa
    end

    x = Arb((a, b))

    # Compute rest term of the Taylor series
    p = f(ArbSeries((x, 1), degree = degree + 1))

    # Check if f happens to be monotonic on the interval, in that case
    # evaluate on the endpoints only
    if !Arblib.contains_zero(Arblib.ref(p, 1))
        verbose && @info "monotonic on interval - evaluate on endpoints"
        fa, fb = f(Arb(a)), f(Arb(b))
        if abs_value
            sgn = _check_signs(fa, fb)
            if sgn == -1
                # The sign differs return zero
                verbose && @info "sign of endpoints differ - minimum is zero"

                return zero(fa), max(abs(fa), abs(fb)), Arb(NaN, prec = precision(a))
            else
                res = minmax(abs(fa), abs(fb))
                Arblib.nonnegative_part!(res[1], res[1])

                return res..., Arb(NaN, prec = precision(a))
            end
        else
            return min(fa, fb), max(fa, fb), Arb(NaN, prec = precision(a))
        end
    end

    remainder = taylor_remainder(p, x)

    # If the remainder is finite the result will never be finite,
    # return the zeroth order enclosure
    if !isfinite(remainder)
        verbose && @info "non-finite remainder"
        res = maybe_abs(p[0])
        return res, copy(res), copy(res)
    end

    # Compute the Taylor series at the midpoint of x
    q = f(ArbSeries((Arblib.midref(x), 1); degree))

    # Compute interval to enclose extrema of Taylor series on
    c, c_exact, d, d_exact = _compute_endpoints(a, b, x)

    # Enclose the extrema of the Taylor series
    res = extrema_polynomial(ArbPoly(q), c, d; abs_value)

    if !c_exact
        # Compute value on an enclosure of the left endpoint
        y = maybe_abs(q(_enclose_inexact_endpoint(c)))
        Arblib.min!(res[1], res[1], y)
        Arblib.max!(res[2], res[2], y)
    end
    if !d_exact
        # Compute value on an enclosure of the right endpoint
        y = maybe_abs(q(_enclose_inexact_endpoint(d)))
        Arblib.min!(res[1], res[1], y)
        Arblib.max!(res[2], res[2], y)
    end

    return res[1] + remainder, res[2] + remainder, maybe_abs(q[0])
end
"""
    minimum_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute the minimum of the function `f` on the interval `[a, b]`.

Takes the same arguments as [`extrema_series`](@ref). The algorithm is
also the same except that it only looks for the minimum.

This method is mainly intended for internal use in the function
[`minimum_enclosure`](@ref) but can be used directly as well.
"""
function minimum_series(
    f,
    a::Arf,
    b::Arf;
    degree::Integer = 8,
    abs_value = false,
    verbose = false,
)
    check_interval(a, b)

    maybe_abs = abs_value ? abs : identity
    if a == b
        fa = maybe_abs(f(Arb(a)))
        return fa, fa
    end

    x = Arb((a, b))

    # Compute rest term of the Taylor series
    p = f(ArbSeries((x, 1), degree = degree + 1))

    # Check if f happens to be monotonic on the interval, in that case
    # evaluate on the endpoints where the maximum could be attained
    if !Arblib.contains_zero(Arblib.ref(p, 1))
        verbose && @info "monotonic on interval - evaluate on endpoints"
        if abs_value
            fa, fb = f(Arb(a)), f(Arb(b))
            sgn = _check_signs(fa, fb)
            if sgn == -1
                # The sign differs return zero
                verbose && @info "sign of endpoints differ - minimum is zero"

                res = zero(fa)
            else
                res = min(abs(fa), abs(fb))
                Arblib.nonnegative_part!(res, res)
            end
        elseif Arblib.ispositive(Arblib.ref(p, 1))
            # Minimum is attained at the left endpoint
            res = f(Arb(a))
        else
            # Minimum is attained at the right endpoint
            res = f(Arb(b))
        end
        return res, Arb(NaN, prec = precision(a))
    end

    remainder = taylor_remainder(p, x)

    # If the remainder is finite the result will never be finite,
    # return the zeroth order enclosure
    if !isfinite(remainder)
        verbose && @info "non-finite remainder"
        res = maybe_abs(p[0])
        return res, copy(res)
    end

    # Compute the Taylor series at the midpoint of x
    q = f(ArbSeries((Arblib.midref(x), 1); degree))

    # Compute interval to enclose minimum of Taylor series on
    c, c_exact, d, d_exact = _compute_endpoints(a, b, x)

    # Enclose the minimum of the Taylor series
    res = minimum_polynomial(ArbPoly(q), c, d; abs_value)

    if !c_exact
        # Compute value on an enclosure of the left endpoint
        Arblib.min!(res, res, maybe_abs(q(_enclose_inexact_endpoint(c))))
    end
    if !d_exact
        # Compute value on an enclosure of the right endpoint
        Arblib.min!(res, res, maybe_abs(q(_enclose_inexact_endpoint(d))))
    end

    return res + remainder, maybe_abs(q[0])
end

"""
    maximum_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute the maximum of the function `f` on the interval `[a, b]`.

Takes the same arguments as [`extrema_series`](@ref). The algorithm is
also the same except that it only looks for the maximum.

This method is mainly intended for internal use in the function
[`maximum_enclosure`](@ref) but can be used directly as well.
"""
function maximum_series(
    f,
    a::Arf,
    b::Arf;
    degree::Integer = 8,
    abs_value = false,
    verbose = false,
)
    check_interval(a, b)

    maybe_abs = abs_value ? abs : identity

    if a == b
        fa = maybe_abs(f(Arb(a)))
        return fa, fa
    end

    x = Arb((a, b))

    # Compute rest term of the Taylor series
    p = f(ArbSeries((x, 1), degree = degree + 1))

    # Check if f happens to be monotonic on the interval, in that case
    # evaluate on the endpoints where the maximum could be attained
    if !Arblib.contains_zero(Arblib.ref(p, 1))
        verbose && @info "monotonic on interval - evaluate on endpoints"
        if abs_value
            # Maximum could be attained at either endpoint
            res = max(abs(f(Arb(a))), abs(f(Arb(b))))
        elseif Arblib.ispositive(Arblib.ref(p, 1))
            # Maximum is attained at the right endpoint
            res = f(Arb(b))
        else
            # Maximum is attained at the left endpoint
            res = f(Arb(a))
        end
        return res, Arb(NaN, prec = precision(a))
    end

    remainder = taylor_remainder(p, x)

    # If the remainder is finite the result will never be finite,
    # return the zeroth order enclosure
    if !isfinite(remainder)
        verbose && @info "non-finite remainder"
        res = maybe_abs(p[0])
        return res, copy(res)
    end

    # Compute the Taylor series at the midpoint of x
    q = f(ArbSeries((Arblib.midref(x), 1); degree))

    # Compute interval to enclose maximum of Taylor series on
    c, c_exact, d, d_exact = _compute_endpoints(a, b, x)

    # Enclose the maximum of the Taylor series
    res = maximum_polynomial(ArbPoly(q), c, d; abs_value)

    if !c_exact
        # Compute value on an enclosure of the left endpoint
        Arblib.max!(res, res, maybe_abs(q(_enclose_inexact_endpoint(c))))
    end
    if !d_exact
        # Compute value on an enclosure of the right endpoint
        Arblib.max!(res, res, maybe_abs(q(_enclose_inexact_endpoint(d))))
    end

    return res + remainder, maybe_abs(q[0])
end


text:
abstract type Service <: Component end

"""
All subtypes of Service define contributingdevices::Vector{Device}. The values get populated
with references to existing devices. The following functions override JSON encoding to
replace the devices with their UUIDs. This solves two problems:

1. The information is redundant, consuming unnecessary space. Also, a user could modify
   information in the JSON file in one place but not the other, which would be very
   problematic.
2. The contributingdevices are stored in an vector of abstract types. JSON doesn't provide a
   way to encode the actual concrete type, so deserialization would have to infer the type.
   There is no guessing if the UUIDs are stored instead. The deserialization process can
   replace the references with actual devices.

These functions could be re-defined to accept subtypes of Component or PowerSystemType.
This is the minimum amount needed for now.
"""

function JSON2.write(io::IO, service::T) where {T <: Service}
    return JSON2.write(io, encode_for_json(service))
end

function JSON2.write(service::T) where {T <: Service}
    return JSON2.write(encode_for_json(service))
end

function encode_for_json(service::T) where {T <: Service}
    fields = fieldnames(T)
    vals = []

    for name in fields
        val = getfield(service, name)
        if val isa Vector{<:Device}
            push!(vals, IS.get_uuid.(val))
        else
            push!(vals, val)
        end
    end

    return NamedTuple{fields}(vals)
end

"""Creates a Service object by decoding the data that was in JSON. This data stores the
values for the field contributingdevices as UUIDs, so this will lookup each device in
devices.
"""
function IS.convert_type(::Type{T}, data::NamedTuple, devices::Dict) where {T <: Service}
    @debug T data
    values = []
    for (fieldname, fieldtype) in zip(fieldnames(T), fieldtypes(T))
        val = getfield(data, fieldname)
        if fieldtype <: Vector{<:Device}
            real_devices = []
            for item in val
                uuid = Base.UUID(item.value)
                service = devices[uuid]
                push!(real_devices, service)
            end
            push!(values, real_devices)
        else
            obj = IS.convert_type(fieldtype, val)
            push!(values, obj)
        end
    end

    return T(values...)
end

function IS.convert_type(::Type{T}, data::Any) where {T <: Service}
    error("This form of convert_type is not supported for Services")
end


text:
using IndexedVectors
using Test

@testset "IndexedVectors.jl" begin
    data = [("a", 100), ("b", 200), ("c", 300), ("1", 400)]
    iv = IndexedVector(data, x -> x[1])

    @test iv[1] == ("a", 100)
    @test iv["b"] == ("b", 200)
    @test iv[1:2] == [("a", 100), ("b", 200)]
    @test iv[[2, 1, 3]] == [("b", 200), ("a", 100), ("c", 300)]
    @test sort(iv) == [("1", 400), ("a", 100), ("b", 200), ("c", 300)]
end


text:
function residuals(model::AModel, calibration)
    m = calibration[:exogenous]
    s = calibration[:states]
    x = calibration[:controls]
    p = calibration[:parameters]
    res = Dolo.arbitrage(model, m, s, x, m, s, x, p)
    S = Dolo.transition(model, m, s, x, m, p)
    return Dict(:arbitrage=>res, :transition=>S-s)
end

function residuals(model::AModel)
    return residuals(model, model.calibration)
end

function find_deterministic_equilibrium(model::AModel, calibration::ModelCalibration)
    m, p, s0, x0 = calibration[:exogenous, :parameters, :states, :controls]
    ns = length(s0)
    nx = length(x0)

    function obj!(sx)
        s = sx[1:ns]
        x = sx[ns+1:ns+nx]

        # update state part of residual
        s_new = transition(model, m, s, x, m, p)

        # now update control part
        resids = arbitrage(model, m, s, x, m, s, x, p)
        vcat(s_new - s, resids)
    end

    sol = nlsolve(obj!, vcat(s0, x0),  autodiff=:central, inplace=false)
    NLsolve.converged(sol) || error("Nonlinear solver failed to find steady state")

    # otherwise set controls
    out = deepcopy(calibration)
    out[:states] = sol.zero[1:ns]
    out[:controls] = sol.zero[ns+1:end]
    out
end

function find_deterministic_equilibrium(model::AModel)
    return find_deterministic_equilibrium(model, model.calibration)
end

# ---------- #
# Docstrings #
# ---------- #

"""
    find_deterministic_equilibrium(model::AModel, [calib::ModelCalibration])

Solve for the steady state equilibrium of `model`, data in `cailb` to fill
in parameter values and provide an initial guess for the states and controls.
When no calibration is passed `model.calibration` is used

The `exogenous` variables at time t-1 (`m`) and t (`M`) are set to
`calib[:exogenous]`.

The deterministic equilibrium is found by solving for vectors `s` and `x`, such
that

1. `s = transition(m, s, x, m, p)`
2. `0 = arbitrage(m, s, x, m, s, x, p)`

"""
find_deterministic_equilibrium

"""
    residuals(model::AModel, [calib::ModelCalibration])::Dict

Compute the steady state residuals for the aribtrage and transition equations
of `model`, when these functions are evaluated at the data in `calib`. If no
`calib` is provided, `model.calibration` will be used.

See the docstring for `find_deterministic_equilibrium` for more information
"""
residuals


text:
module KUDNN
using Compat
using CUDArt,CUDNN

const libkudnn = Libdl.find_library(["libkudnn"], [Pkg.dir("KUDNN/src")])
isempty(libkudnn) && error("Cannot find libkudnn")

export kudnnConvolutionForward, kudnnConvolutionBackwardFilter, kudnnConvolutionBackwardData, kudnnGetConvolutionNdForwardOutputDim
export kudnnPoolingForward, kudnnPoolingBackward, kudnnGetPoolingNdForwardOutputDim

# high level C interface (following CUDNN.jl)

using CUDNN: cudnnHandle, ConvolutionDescriptor, CUDNN_CROSS_CORRELATION, CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM, TD, FD, cptr

# TODO: do not construct this every time:
kudnnConvolutionDescriptor(n)=ConvolutionDescriptor(; padding=ntuple(i->0,n), stride=ntuple(i->1,n), upscale=ntuple(i->1,n), mode=CUDNN_CROSS_CORRELATION)

function kudnnGetConvolutionNdForwardOutputDim(src, filter; convDesc=kudnnConvolutionDescriptor(ndims(src)-2))
    nd = ndims(src)
    nd == ndims(filter) || error("Dimension mismatch")
    size(src,nd-1) == size(filter,nd-1) || error("Dimension mismatch")
    od = Array(Int,nd)
    for i=1:nd-2
        od[i] = 1 + div(size(src,i) - size(filter,i) + 2*convDesc.padding[i], convDesc.stride[i])
    end
    od[nd-1] = size(filter,nd)
    od[nd] = size(src,nd)
    tuple(od...)
end

function kudnnConvolutionForward(src::AbstractCudaArray, filter::AbstractCudaArray, dest=nothing;
                                 handle=cudnnHandle, alpha=1.0, beta=0.0,
                                 convDesc=kudnnConvolutionDescriptor(ndims(src)-2),
                                 algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM,
                                 workSpace=Int8[], workSpaceSizeInBytes=0)
    @assert ndims(filter) == ndims(src)
    @assert eltype(filter) == eltype(src)
    osize = kudnnGetConvolutionNdForwardOutputDim(src,filter;convDesc=convDesc)
    (dest == nothing) && (dest = CudaArray(eltype(src), osize))
    @assert osize == size(dest)
    @assert eltype(dest) == eltype(src)
    kudnnConvolutionForward(handle,
                            cptr(alpha,src),TD(src),src,
                            FD(filter),filter,
                            convDesc,algorithm,workSpace,workSpaceSizeInBytes,
                            cptr(beta,dest),TD(dest),dest)
    return dest
end

function kudnnConvolutionBackwardFilter(src::AbstractCudaArray, diff::AbstractCudaArray, grad::AbstractCudaArray;
                                        handle=cudnnHandle, alpha=1.0, beta=0.0,
                                        convDesc=kudnnConvolutionDescriptor(ndims(src)-2))
    kudnnConvolutionBackwardFilter(handle,
                                   cptr(alpha,src),TD(src),src,
                                   TD(diff),diff,convDesc,
                                   cptr(beta,grad),FD(grad),grad)
    return grad
end

function kudnnConvolutionBackwardData(filter::AbstractCudaArray, diff::AbstractCudaArray, grad::AbstractCudaArray;
                                      handle=cudnnHandle, alpha=1.0, beta=0.0,
                                      convDesc=kudnnConvolutionDescriptor(ndims(filter)-2))
    kudnnConvolutionBackwardData(handle,cptr(alpha,diff),
                                 FD(filter),filter,
                                 TD(diff),diff,convDesc,
                                 cptr(beta,grad),TD(grad),grad)
    return grad
end

function kudnnGetPoolingNdForwardOutputDim(pd::PoolingDescriptor, input::AbstractCudaArray)
    n = ndims(input)
    ntuple(n) do i
        i >= n-1 ? size(input, i) :
        1 + ceil(Int, (size(input, i) + 2*pd.padding[i] - pd.dims[i]) / pd.stride[i])
    end
end

function kudnnPoolingForward(pd::PoolingDescriptor, src::AbstractCudaArray, dest::AbstractCudaArray;
                             handle=cudnnHandle, alpha=1.0, beta=0.0)
    kudnnPoolingForward(handle, pd,
                        cptr(alpha,src), TD(src), src,
                        cptr(beta,dest), TD(dest), dest)
    return dest
end

function kudnnPoolingBackward(pd::PoolingDescriptor, src::AbstractCudaArray, srcDiff::AbstractCudaArray, dest::AbstractCudaArray, destDiff::AbstractCudaArray;
                              handle=cudnnHandle, alpha=1.0, beta=0.0)
    cudnnSetTensor(destDiff, 0)  # TODO: make sure kudnnPoolingBackward zeros out its output
    kudnnPoolingBackward(handle, pd,
                         cptr(alpha,src), TD(src), src,
                         TD(srcDiff), srcDiff,
                         TD(dest), dest,
                         cptr(beta,destDiff), TD(destDiff), destDiff)
    return destDiff
end

# low level C interface (similar to libcudnn.jl)

using CUDNN: cudnnCheck,cudnnStatus_t,cudnnHandle_t,cudnnTensorDescriptor_t,cudnnFilterDescriptor_t,cudnnConvolutionDescriptor_t,cudnnConvolutionFwdAlgo_t,cudnnPoolingDescriptor_t

function kudnnConvolutionForward(handle,alpha,srcDesc,srcData,filterDesc,filterData,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,destDesc,destData)
    cudnnCheck(ccall((:kudnnConvolutionForward,libkudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,cudnnConvolutionFwdAlgo_t,Ptr{Void},Csize_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,srcDesc,srcData,filterDesc,filterData,convDesc,algo,workSpace,workSpaceSizeInBytes,beta,destDesc,destData))
end

function kudnnConvolutionBackwardFilter(handle,alpha,srcDesc,srcData,diffDesc,diffData,convDesc,beta,gradDesc,gradData)
    cudnnCheck(ccall((:kudnnConvolutionBackwardFilter,libkudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void}),handle,alpha,srcDesc,srcData,diffDesc,diffData,convDesc,beta,gradDesc,gradData))
end

function kudnnConvolutionBackwardData(handle,alpha,filterDesc,filterData,diffDesc,diffData,convDesc,beta,gradDesc,gradData)
    cudnnCheck(ccall((:kudnnConvolutionBackwardData,libkudnn),cudnnStatus_t,(cudnnHandle_t,Ptr{Void},cudnnFilterDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnConvolutionDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,alpha,filterDesc,filterData,diffDesc,diffData,convDesc,beta,gradDesc,gradData))
end

function kudnnPoolingForward(handle,poolingDesc,alpha,srcDesc,srcData,beta,destDesc,destData)
    cudnnCheck(ccall((:kudnnPoolingForward,libkudnn),cudnnStatus_t,(cudnnHandle_t,cudnnPoolingDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,poolingDesc,alpha,srcDesc,srcData,beta,destDesc,destData))
end

function kudnnPoolingBackward(handle,poolingDesc,alpha,srcDesc,srcData,srcDiffDesc,srcDiffData,destDesc,destData,beta,destDiffDesc,destDiffData)
    cudnnCheck(ccall((:kudnnPoolingBackward,libkudnn),cudnnStatus_t,(cudnnHandle_t,cudnnPoolingDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void},Ptr{Void},cudnnTensorDescriptor_t,Ptr{Void}),handle,poolingDesc,alpha,srcDesc,srcData,srcDiffDesc,srcDiffData,destDesc,destData,beta,destDiffDesc,destDiffData))
end

end # module KUDNN

max_stars_repo_name: denizyuret/KUDNN

max_issues_repo_name: denizyuret/KUDNN

max_forks_repo_name: denizyuret/KUDNN


text:
using DistributionVectors
using Test
using Distributions
#using MissingStrategies
using LinearAlgebra

if !@isdefined(MyNewDistribution)
  struct MyNewDistribution <: ContinuousUnivariateDistribution; end
end

@testset "sumDistributionVector: not defined for D" begin
  dv = SimpleDistributionVector(MyNewDistribution(), MyNewDistribution());
  acf = AutoCorrelationFunction([1,0.4,0.1]);
  nsum = 5
  corrmat = cormatrix_for_acf(nsum, acf);
  #@test_throws ErrorException sum(dv, acf, PassMissing())
  @test_throws ErrorException sum(dv, Symmetric(corrmat))
  @test_throws ErrorException sum(dv, acf)
  @test_throws ErrorException sum(dv)
  #@test_throws ErrorException mean(dv, acf, PassMissing())
  @test_throws ErrorException mean(dv, Symmetric(corrmat))
  @test_throws ErrorException mean(dv, acf)
  @test_throws ErrorException mean(dv)
end;


text:
using Test, Bonsai, StructTypes
using StructTypes: @Struct
using Bonsai: json_schema, JSONSchema, json_schema_format
using NamedTupleTools
using Dates

@enum Fruit apple orange

"""
Fruit basket order
"""
struct FruitBasket
	total::Int
	prices::Vector{Int}
	fruit::Vector{Fruit}
end

@Struct FruitBasket

@testset "json_schema_format" begin
	@test json_schema_format(DateTime) == "date-time"
end

@testset "json_schema" begin

	@test !isnothing(json_schema(Fruit).enum)
	@test !isnothing(json_schema(Vector{Int}).items)

	@test !isnothing(json_schema(Tuple{Int, String}).prefixItems)
	@test Bonsai.json_schema(typeof(Body(x=10))) isa JSONSchema
	@test Bonsai.json_schema(typeof(Headers(x_next=12))) isa JSONSchema
end

@testset "doc_str" begin
	@test !isnothing(Bonsai.doc_str(FruitBasket))
end

@testset "JSONSchema" begin
	@test json_schema(FruitBasket) isa JSONSchema
	@test length(Bonsai.union_types(Union{Missing, Int, Nothing})) == 3
    @test length(Bonsai.json_schema(Union{Missing, Int, Nothing}).oneOf) == 2
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using SystemSnoop
using Test
using Dates: Dates
using StructArrays

@testset "Testing SmartSample" begin
    times = Dates.DateTime[]
    timestep = Dates.Millisecond(100)
    sampler = SystemSnoop.SmartSample(timestep)
    for i in 1:10
        sleep(sampler)
        push!(times, Dates.now())
    end

    # The first sample is unreliable for some reason - just trash it.
    gaps = Dates.value.(Dates.Millisecond.(diff(times)))
    popfirst!(gaps)

    for g in gaps
        @test Dates.value(timestep) - 2 <= g
        @test g <= Dates.value(timestep) + 2
    end

    # Now, make sure the logic for skipping time steps works. The following should work
    start_iter = sampler.iteration
    sleep(0.1)
    sleep(sampler)
    stop_iter = sampler.iteration

    # The log here is that the smart sampler should skip iterations.
    @test stop_iter > start_iter + 1
end

#####
##### Measurement Routine
#####

# Require that `kw` are forwarded
mutable struct Incrementer
    prepare_count::Int
    measurement_count::Int
    clean_count::Int
end

function SystemSnoop.prepare(x::Incrementer)
    x.prepare_count += 1
    return nothing
end

function SystemSnoop.measure(x::Incrementer)
    x.measurement_count += 1
    return x.measurement_count
end

function SystemSnoop.clean(x::Incrementer)
    x.clean_count += 1
    return nothing
end

@testset "snooped" begin
    # Use the `TestA` code from base.jl
    incrementer = Incrementer(0, 0, 0)
    measurements = (
        timestamp = SystemSnoop.Timestamp(),
        incrementer = incrementer,
    )

    # Call the outer `snoop` function with our dummy process. Make sure everything runs
    # smoothly.
    iters = 3
    increment = 20
    trace = @snooped measurements 1000 sleep(2)

    # Test that inference works correctly
    expected = Base.promote_op(
        SystemSnoop.snooploop,
        typeof(measurements),
        typeof(1000),
        NamedTuple{},
        Threads.Atomic{Int},
        Threads.Atomic{Int},
    )
    @test typeof(trace) == expected

    # Were the correct methods called?
    @test incrementer.prepare_count == 1
    @test incrementer.clean_count == 1

    num_measurements = incrementer.measurement_count
    @test num_measurements > 0
    @test trace.incrementer == 1:num_measurements
end


text:
function density(obj::MultivariateRealDistribution, arg0::Vector{jdouble})
    return jcall(obj, "density", jdouble, (Vector{jdouble},), arg0)
end

function get_dimension(obj::MultivariateRealDistribution)
    return jcall(obj, "getDimension", jint, ())
end

function reseed_random_generator(obj::MultivariateRealDistribution, arg0::jlong)
    return jcall(obj, "reseedRandomGenerator", void, (jlong,), arg0)
end

function sample(obj::MultivariateRealDistribution)
    return jcall(obj, "sample", Vector{jdouble}, ())
end

function sample(obj::MultivariateRealDistribution, arg0::jint)
    return jcall(obj, "sample", Vector{Vector{jdouble}}, (jint,), arg0)
end


text:
digits_bw_labels() = CSV.read(datadir("digits_bw/labels.csv"), DataFrame)
digits_rgb_labels() = CSV.read(datadir("digits_rgb/labels.csv"), DataFrame)
morpho_mnist_labels() = CSV.read(datadir("morpho_mnist/labels.csv"), DataFrame)

"""
	filter_data(full_labels::DataFrame, filter_dict::Dict)

Filter the data given the full labels and a dict of filters.

Example:
filter_dict_a = Dict(
    "shape" => [1],
    "normalized_orientation" => ["<=", 0.1],
    "scale" => ["==", 0.8],
    "posX" => ["<", 0.5],
    "posY" => ["<", 0.5],
    )
"""
function filter_data(full_labels, filter_dict; categorical_key=nothing)
    ks = keys(filter_dict)
    bin_inds = !isnothing(categorical_key)  ?
        map(x-> x in filter_dict[categorical_key], full_labels[!,Symbol(categorical_key)]) :
        Bool.(ones(size(full_labels,1)))
    for k in filter(k->k != categorical_key, ks)
        kp = filter_dict[k]
        #ff(x) = eval(Meta.parse("$x $(kp[1]) $(kp[2])")) # this is slow
        val = Float64(Meta.parse("$(kp[2])"))
        # this is ugly but fast
        ff(x) = if kp[1] == "=="
            x == val
        elseif kp[1] == "<="
            x <= val
        elseif kp[1] == "<"
            x < val
        elseif kp[1] == ">="
            x >= val
        elseif kp[1] == ">"
            x > val
        elseif kp[1] == "!="
            x != val
        else
            error("Unknown operator, please support one of [==, <, <=, >, >=, !=]")
        end
        bin_inds = bin_inds .& map(ff, full_labels[!,Symbol(k)])
    end
    return bin_inds
end

hyphen_join(x,y) = "$(x)-$(y)"

"""
	safe_savename(params, args...; kwargs...)

Equivalent to DrWatson.savename but preserves arrays and tuples.
"""
function safe_savename(params, args...; kwargs...)
    outparams = Dict()
    for (k,v) in pairs(params)
        if typeof(v) <: Vector
            outparams[k] = reduce(hyphen_join,v)
        else
            outparams[k] = v
        end
    end
    return replace(savename(outparams, args...; kwargs...), " "=>"")
end

function _load_mnist(type::String)
	(type in ["train", "test"]) ? nothing : error("Suggested MNIST data type $type not available.")
    x = eval(Meta.parse("MLDatasets.MNIST.$(type)tensor(Float32)"))
    x = permutedims(x, (2,1,3))
    s = size(x)
    reshape(x, s[1], s[2], 1, s[3])
end

"""
	load_mnist(type="")

Loads either all/test/train MNIST data.
"""
function load_mnist(type="")
	if type == ""
		return cat(_load_mnist("train"), _load_mnist("test"), dims=4)
	else
		return _load_mnist(type)
	end
end

"""
	load_train_val_test_data(dataset, filter_dict=Dict(); ratios=(0.6,0.2,0.2), seed=nothing,
        categorical_key=nothing)

Dataset is one of "morhpo_mnist"/"digits_bw"/"digits_rgb".
"""
function load_train_val_test_data(dataset, filter_dict=Dict(); ratios=(0.6,0.2,0.2),
    seed=nothing, categorical_key=nothing)
    if dataset == "morpho_mnist"
	    full_data = load_mnist()
    elseif dataset == "shapes2D"
        full_data,_ = load_shapes2D()
	else
		raw_data = load(datadir("$(dataset)/digits.bson"))[:digits]
		N = length(raw_data)
		full_data = zeros(Float32, size(raw_data[1])..., N);
		for i in 1:N
		    full_data[:,:,:,i] .= raw_data[i]
		end
	end
    full_labels = CSV.read(datadir("$(dataset)/labels.csv"), DataFrame)
    included_inds = filter_data(full_labels, filter_dict; categorical_key=categorical_key)

    # now split the data
    normal_data = full_data[:,:,:,included_inds]
    anomalous_data = full_data[:,:,:,.!included_inds]

    # further split the normal into train and val
    trinds, valinds, tstinds = train_val_test_inds(1:size(normal_data,4),
        ratios; seed=seed)
    tr_x = normal_data[:,:,:,trinds]
    val_x = normal_data[:,:,:,valinds]
    tst_x = normal_data[:,:,:,tstinds]
    a_x = anomalous_data

    tr_y = full_labels[included_inds,:][trinds,:]
	val_y = full_labels[included_inds,:][valinds,:]
	tst_y = full_labels[included_inds,:][tstinds,:]
	a_y = full_labels[.!included_inds,:]

    (tr_x, tr_y), (val_x, val_y), (tst_x, tst_y), (a_x, a_y)
end

# stuff for run scripts
function get_filter_info(experiment_argnames, args, arg_table)
    # construct the Dictionary with filters and their values
    filter_keys = filter(k->!(k in experiment_argnames),keys(args))
    filter_dict = Dict(zip(filter_keys, [args[k] for k in filter_keys]))

    # also, set which arguments are non-default
    non_default_filters = []
    for k in filter_keys
        argind = findfirst(map(f->f.dest_name == k,arg_table.args_table.fields))
        (arg_table.args_table.fields[argind].default == filter_dict[k]) ? nothing :
            push!(non_default_filters, k)
    end

    filter_dict, non_default_filters
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
using MeshViz
using Meshes
using GeoTables
using CategoricalArrays
using ReferenceTests
using ImageIO
using Random
using Test

import CairoMakie as Mke

@testset "MeshViz.jl" begin
  # 1D point set
  p = PointSet(Point1[(0.,), (1.,), (2.,)])
  # TODO

  # 2D point set
  p = PointSet(Point2[(0.,0.), (1.,0.), (1.,1.), (0.,1.)])
  @test_reference "data/pset2D-1.png" viz(p)
  @test_reference "data/pset2D-2.png" viz(p, color = :red)
  @test_reference "data/pset2D-3.png" viz(p, color = 1:4)
  @test_reference "data/pset2D-4.png" viz(p, color = 1:4, colorscheme = :inferno)
  @test_reference "data/pset2D-5.png" viz(p, color = :red, alpha = 0.5)
  @test_reference "data/pset2D-6.png" viz(p, color = 1:4, alpha = 0.5)

  # 3D point set
  p = PointSet(Point3[(0.,0.,0.), (1.,0.,0.), (1.,1.,0.), (0.,1.,0.),
                      (0.,0.,1.), (1.,0.,1.), (1.,1.,1.), (0.,1.,1.)])
  @test_reference "data/pset3D-1.png" viz(p)
  @test_reference "data/pset3D-2.png" viz(p, color = :red)
  @test_reference "data/pset3D-3.png" viz(p, color = 1:8)
  @test_reference "data/pset3D-4.png" viz(p, color = 1:8, colorscheme = :inferno)
  @test_reference "data/pset3D-5.png" viz(p, color = :red, alpha = 0.5)
  @test_reference "data/pset3D-6.png" viz(p, color = 1:8, alpha = 0.5)

  # 1D Cartesian grid
  d = CartesianGrid(10)
  # TODO

  # 2D Cartesian grid
  d = CartesianGrid(10,10)
  @test_reference "data/grid2D-1.png" viz(d)
  @test_reference "data/grid2D-2.png" viz(d, showfacets = true)
  @test_reference "data/grid2D-3.png" viz(d, showfacets = true, facetcolor = :red)
  @test_reference "data/grid2D-4.png" viz(d, color = 1:100)
  @test_reference "data/grid2D-5.png" viz(d, color = 1:100, colorscheme = :inferno)
  @test_reference "data/grid2D-6.png" viz(d, color = :red)
  @test_reference "data/grid2D-7.png" viz(d, color = :red, alpha = 0.5)
  @test_reference "data/grid2D-8.png" viz(d, color = 1:100, alpha = 0.5)

  # 3D Cartesian grid
  d = CartesianGrid(10,10,10)
  @test_reference "data/grid3D-1.png" viz(d)
  @test_reference "data/grid3D-2.png" viz(d, showfacets = true)
  @test_reference "data/grid3D-3.png" viz(d, showfacets = true, facetcolor = :red)
  @test_reference "data/grid3D-4.png" viz(d, color = 1:1000)
  @test_reference "data/grid3D-5.png" viz(d, color = 1:1000, colorscheme = :inferno)
  @test_reference "data/grid3D-6.png" viz(d, color = :red)
  @test_reference "data/grid3D-7.png" viz(d, color = :red, alpha = 0.5)
  @test_reference "data/grid3D-8.png" viz(d, color = 1:1000, alpha = 0.5)

  # 2D N-gons
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  @test_reference "data/tri2D-1.png" viz(t)
  @test_reference "data/tri2D-2.png" viz(t, showboundary = false)
  @test_reference "data/tri2D-3.png" viz(t, color = :orange)
  @test_reference "data/tri2D-4.png" viz(t, color = :cyan, boundarycolor = :red)
  @test_reference "data/tri2D-5.png" viz(t, color = :orange, alpha = 0.5)
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  @test_reference "data/quad2D-1.png" viz(q)
  @test_reference "data/quad2D-2.png" viz(q, showboundary = false)
  @test_reference "data/quad2D-3.png" viz(q, color = :orange)
  @test_reference "data/quad2D-4.png" viz(q, color = :cyan, boundarycolor = :red)
  @test_reference "data/quad2D-5.png" viz(q, color = :orange, alpha = 0.5)

  # 3D N-gons
  o = Octagon([(0.0,0.0,1.0), (0.5,-0.5,0.0), (1.0,0.0,0.0), (1.5,0.5,-0.5),
               (1.0,1.0,0.0), (0.5,1.5,0.0), (0.0,1.0,0.0), (-0.5,0.5,0.0)])
  @test_reference "data/oct3D-1.png" viz(o)
  @test_reference "data/oct3D-2.png" viz(o, showboundary = false)
  @test_reference "data/oct3D-3.png" viz(o, color = :orange)
  @test_reference "data/oct3D-4.png" viz(o, color = :cyan, boundarycolor = :red)
  @test_reference "data/oct3D-5.png" viz(o, color = :orange, alpha = 0.5)

  # Polygonal areas
  Random.seed!(2020)
  p = PolyArea((0.,0.), (0.5,-1.5), (1.,0.), (1.5,0.5),
               (1.,1.), (0.5,1.5), (-0.5,0.5), (0.,0.))
  @test_reference "data/poly2D-1.png" viz(p)
  @test_reference "data/poly2D-2.png" viz(p, showboundary = false)
  @test_reference "data/poly2D-3.png" viz(p, color = :orange)
  @test_reference "data/poly2D-4.png" viz(p, color = :cyan, boundarycolor = :red)
  @test_reference "data/poly2D-5.png" viz(p, color = :orange, alpha = 0.5)

  # Multi-geometries
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  m = Multi([t, q])
  @test_reference "data/multi2D-1.png" viz(m)
  @test_reference "data/multi2D-2.png" viz(m, showboundary = false)
  @test_reference "data/multi2D-3.png" viz(m, color = :orange)
  @test_reference "data/multi2D-4.png" viz(m, color = :cyan, boundarycolor = :red)
  @test_reference "data/multi2D-5.png" viz(m, color = :orange, alpha = 0.5)

  # 2D boxes
  b = Box((0.,0.), (1.,1.))
  @test_reference "data/box2D-1.png" viz(b)
  @test_reference "data/box2D-2.png" viz(b, showboundary = false)
  @test_reference "data/box2D-3.png" viz(b, color = :orange)
  @test_reference "data/box2D-4.png" viz(b, color = :cyan, boundarycolor = :red)
  @test_reference "data/box2D-5.png" viz(b, color = :orange, alpha = 0.5)

  # 3D boxes
  b = Box((0.,0.,0.), (1.,1.,1.))
  @test_reference "data/box3D-1.png" viz(b)
  @test_reference "data/box3D-2.png" viz(b, color = :orange)
  @test_reference "data/box3D-3.png" viz(b, color = :orange, alpha = 0.5)

  # 3D spheres
  s = Sphere((0.,0.,0.), 1.)
  @test_reference "data/sphere3D-1.png" viz(s)
  @test_reference "data/sphere3D-2.png" viz(s, color = :orange)
  @test_reference "data/sphere3D-3.png" viz(s, color = :orange, alpha = 0.5)

  ## Polyhedron
  # Triangle
  o = Triangle([(0.0,0.0,1.0), (0.0,1.0,0.0), (1.0,0.0,0.0)])
  viz(o)
  viz(o, showboundary = false)
  viz(o, color = :orange)
  viz(o, color = :cyan, boundarycolor = :red)

  # Tetrahedron
  t = Tetrahedron([(0.0,0.0,0.0), (0.0,0.0,1.0), (0.0,1.0,0.0), (1.0,0.0,0.0)])
  viz(t)
  viz(t, showboundary = false)
  viz(t, color = :orange)
  viz(t, color = :cyan, boundarycolor = :red)

  # Hexahedron
  # o = Hexahedron([(0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
  #                 (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0)])
  h = Hexahedron([(5.0,0.0,0.0), (6.0,0.0,0.0), (6.0,1.0,0.0), (5.0,1.0,0.0),
                  (5.0,0.0,1.0), (6.0,0.0,1.0), (6.0,1.0,1.0), (5.0,1.0,1.0)])

  viz([t, h], color = [:blue, :red])
  viz([t, h], showboundary = false)
  viz([t, h], color = 1:2)
  viz([t, h], color = :cyan, boundarycolor = :red)

  # Collections of geometries
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  d = Collection([t, q])
  @test_reference "data/collec2D-1.png" viz(d)
  @test_reference "data/collec2D-2.png" viz(d, showfacets = true)
  @test_reference "data/collec2D-3.png" viz(d, showfacets = true, facetcolor = :red)
  @test_reference "data/collec2D-4.png" viz(d, color = 1:2)
  @test_reference "data/collec2D-5.png" viz(d, color = 1:2, colorscheme = :inferno)
  @test_reference "data/collec2D-6.png" viz(d, color = [:red,:green], alpha = 0.5)
  @test_reference "data/collec2D-7.png" viz(d, color = 1:2, alpha = 0.5)

  # Simple meshes
  s = Sphere((0.,0.,0.), 1.)
  m = discretize(s, RegularDiscretization(10))
  nv = nvertices(m)
  ne = nelements(m)
  @test_reference "data/mesh3D-1.png" viz(m)
  @test_reference "data/mesh3D-2.png" viz(m, color = 1:nv)
  @test_reference "data/mesh3D-3.png" viz(m, color = 1:ne)
  @test_reference "data/mesh3D-4.png" viz(m, color = 1:ne, showfacets = true)
  @test_reference "data/mesh3D-5.png" viz(m, color = :orange, showfacets = true, facetcolor = :cyan)
  @test_reference "data/mesh3D-6.png" viz(m, color = :orange, alpha = 0.5)
  @test_reference "data/mesh3D-7.png" viz(m, color = 1:ne, alpha = 0.5)
  @test_reference "data/mesh3D-8.png" viz(m, color = 1:ne, alpha = range(0.1, 1.0, length=ne))

  # Simple meshes 2
  ps  = Meshes.Point.([(0.0,0.0,0.0), (0.0,0.0,1.0), (0.0,1.0,0.0), (1.0,0.0,0.0)])
  c   = connect.([(1, 2, 3, 4), ], Tetrahedron)
  m   = SimpleMesh(ps, c)
  viz(m, shading = false)
  viz(m, showboundary = false)
  viz(m, color = :orange)
  viz(m, color = :cyan, boundarycolor = :red)

  # Simple meshes 3
  ps  = Meshes.Point.([ (0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
                        (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0)])
  c   = connect.([(1, 2, 3, 4, 5, 6, 7, 8), ], Hexahedron)
  m   = SimpleMesh(ps, c)
  viz(m, shading = false)
  viz(m, showboundary = false)
  viz(m, color = :orange)
  viz(m, color = :cyan, boundarycolor = :red)

  # 2D partitions
  Random.seed!(2020)
  g = CartesianGrid(10, 10)
  p = partition(g, PlanePartition((1.,1.)))
  @test_reference "data/part2D-1.png" viz(p)
  @test_reference "data/part2D-2.png" viz(p, alpha = 0.2)

  # 3D partitions
  Random.seed!(2020)
  g = CartesianGrid(10, 10, 10)
  p = partition(g, PlanePartition((1.,1.,1.)))
  @test_reference "data/part3D-1.png" viz(p)
  @test_reference "data/part3D-2.png" viz(p, alpha = 0.5)

  # Vector of points
  p = Point2[(0.,0.), (1.,0.), (1.,1.), (0.,1.)]
  @test_reference "data/points2D-1.png" viz(p)
  @test_reference "data/points2D-2.png" viz(p, color = 1:4)
  @test_reference "data/points2D-3.png" viz(p, color = 1:4, colorscheme = :inferno)
  @test_reference "data/points2D-4.png" viz(p, size = 20)
  @test_reference "data/points2D-5.png" viz(p, color = :red)
  @test_reference "data/points2D-6.png" viz(p, color = :red, alpha = 0.5)
  @test_reference "data/points2D-7.png" viz(p, color = 1:4, alpha = 0.5)

  # Vector of geometries
  t = Triangle((1.,0.), (2.,0.), (2.,1.))
  q = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  g = [t, q]
  @test_reference "data/geoms2D-1.png" viz(g)
  @test_reference "data/geoms2D-2.png" viz(g, color = 1:2)
  @test_reference "data/geoms2D-3.png" viz(g, color = 1:2, colorscheme = :inferno)
  @test_reference "data/geoms2D-4.png" viz(g, color = :red)
  @test_reference "data/geoms2D-5.png" viz(g, color = :red, alpha = 0.5)
  @test_reference "data/geoms2D-6.png" viz(g, color = 1:2, alpha = 0.5)

  # Data over grid
  d = meshdata(CartesianGrid(10,10), etable = (z=1:100,w=1:100))
  @test_reference "data/griddata2D-1.png" viz(d, variable = :z)
  @test_reference "data/griddata2D-2.png" viz(d, variable = :w)
  @test_reference "data/griddata2D-3.png" viz(d, colorscheme = :inferno)
  @test_reference "data/griddata2D-4.png" viz(d, variable = :z, alpha = 0.5)
  @test_reference "data/griddata2D-5.png" viz(d, variable = :z, colorscheme = :inferno, alpha = 0.5)

  # Data over point set
  p = PointSet(centroid.(CartesianGrid(10,10)))
  d = meshdata(p, etable = (z=1:100,))
  @test_reference "data/psetdata2D-1.png" viz(d)
  @test_reference "data/psetdata2D-2.png" viz(d, colorscheme = :inferno)
  @test_reference "data/psetdata2D-3.png" viz(d, variable = :z, alpha = 0.5)
  @test_reference "data/psetdata2D-4.png" viz(d, variable = :z, colorscheme = :inferno, alpha = 0.5)

  # custom values as colors
  d = CartesianGrid(10,10)
  c = categorical(rand(1:4, 100))
  @test_reference "data/values-1.png" viz(d, color = c, colorscheme = :Accent_4)
  c = categorical(rand(1:10, 100))
  @test_reference "data/values-2.png" viz(d, color = c, colorscheme = :BrBG_10)
  c = [fill(missing, 50); categorical(rand(1:4, 50))]
  @test_reference "data/values-3.png" viz(d, color = c, colorscheme = :BrBG_4)
  d = CartesianGrid(2,2)
  c = [1,missing,3,NaN]
  @test_reference "data/values-4.png" viz(d, color = c, colorscheme = :Accent_4)
end


text:
module MDL

setprecision(BigFloat, 64)

"""
scalable binomial coefficient
"""
function choose(n::Int64, k::Int64)
    if 0 < k <= n
        p = BigInt(1)
        for t=0:(min(k,n-k)-1)
            p = floor((p * (n - t)) / (t + 1))
        end
        return p
    else
        return 0.
    end
end

"""
binary logarithm
"""
function log2(n::Real)
    return log(2, n)
end

"""
binary logarithm returning 0. at zero
"""
function log2_zero(n::Real)
    if n == 0
        return 0.
    else
        return log2(n)
    end
end

"""
binary logarithm of binomial coefficient
"""
function log2_choose(n::Int64, k::Int64)
    if n == 0 || k == 0
        return 0.
    else
        return log2_zero(choose(n,k))
    end
end

"""
iterated logarithm (log-star)
"""
function log2_star(n::Union{Int64,Float64,BigInt,BigFloat})
    if n <= 1
        return 0.
    else
        return 1. + log2_star(log2(n))
    end
end

"""
:param z: an int >= 1
:return: the encoded size of the int according to Rissanen 1983's universal code for integers
"""
function universal_integer(z::Int64)
    NORMALIZATION_CONSTANT = 2.865064
    if z <= 0
        return 0.
    end
    c = log2(NORMALIZATION_CONSTANT)
    logstar_z = log2_star(z)
    return logstar_z + c
end

end

max_stars_repo_name: dataspider/momo

max_issues_repo_name: dataspider/momo

max_forks_repo_name: dataspider/momo


text:
const VTK_CELL_TYPE = Dict(
    0=>WriteVTK.VTKCellTypes.VTK_EMPTY_CELL,
    1=>WriteVTK.VTKCellTypes.VTK_VERTEX,
    2=>WriteVTK.VTKCellTypes.VTK_POLY_VERTEX,
    3=>WriteVTK.VTKCellTypes.VTK_LINE,
    4=>WriteVTK.VTKCellTypes.VTK_POLY_LINE,
    5=>WriteVTK.VTKCellTypes.VTK_TRIANGLE,
    6=>WriteVTK.VTKCellTypes.VTK_TRIANGLE_STRIP,
    7=>WriteVTK.VTKCellTypes.VTK_POLYGON,
    8=>WriteVTK.VTKCellTypes.VTK_PIXEL,
    9=>WriteVTK.VTKCellTypes.VTK_QUAD,
    10=>WriteVTK.VTKCellTypes.VTK_TETRA,
    11=>WriteVTK.VTKCellTypes.VTK_VOXEL,
    12=>WriteVTK.VTKCellTypes.VTK_HEXAHEDRON,
    13=>WriteVTK.VTKCellTypes.VTK_WEDGE,
    14=>WriteVTK.VTKCellTypes.VTK_PYRAMID,
    15=>WriteVTK.VTKCellTypes.VTK_PENTAGONAL_PRISM,
    16=>WriteVTK.VTKCellTypes.VTK_HEXAGONAL_PRISM,
    21=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_EDGE,
    22=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_TRIANGLE,
    23=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_QUAD,
    24=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_TETRA,
    25=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON,
    26=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_WEDGE,
    27=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_PYRAMID,
    28=>WriteVTK.VTKCellTypes.VTK_BIQUADRATIC_QUAD,
    29=>WriteVTK.VTKCellTypes.VTK_TRIQUADRATIC_HEXAHEDRON,
    30=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_LINEAR_QUAD,
    31=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_LINEAR_WEDGE,
    32=>WriteVTK.VTKCellTypes.VTK_BIQUADRATIC_QUADRATIC_WEDGE,
    33=>WriteVTK.VTKCellTypes.VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON,
    34=>WriteVTK.VTKCellTypes.VTK_BIQUADRATIC_TRIANGLE,
    35=>WriteVTK.VTKCellTypes.VTK_CUBIC_LINE,
    36=>WriteVTK.VTKCellTypes.VTK_QUADRATIC_POLYGON,
    41=>WriteVTK.VTKCellTypes.VTK_CONVEX_POINT_SET,
    42=>WriteVTK.VTKCellTypes.VTK_POLYHEDRON,
    51=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_CURVE,
    52=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_SURFACE,
    53=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_TRI_SURFACE,
    54=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_QUAD_SURFACE,
    55=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_TETRA_REGION,
    56=>WriteVTK.VTKCellTypes.VTK_PARAMETRIC_HEX_REGION,
    60=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_EDGE,
    61=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_TRIANGLE,
    62=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_QUAD,
    63=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_POLYGON,
    64=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_TETRAHEDRON,
    65=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_WEDGE,
    66=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_PYRAMID,
    67=>WriteVTK.VTKCellTypes.VTK_HIGHER_ORDER_HEXAHEDRON)

const POINT_CELLS = Int[1, 2, 41]
const LINE_CELLS = Int[3, 4, 21, 35, 51, 60]
const FACE_CELLS = Int[5, 6, 7, 8, 9, 22, 23, 28, 30, 34, 36, 52, 53, 54, 61, 62, 63]
const VOLUME_CELLS = Int[10, 11, 12, 13, 14, 15, 16, 24, 25, 26, 27, 29, 31, 32, 33, 42, 55, 56, 64, 65, 66, 67]
const POLY_CELLS = Int[1, 2, 3, 4, 5, 6, 7, 8, 9]

const LINEAR_CELLS = [0:16...]
const QUADRATIC_CELLS = [21:34...; 36]
const CUBIC_CELLS = [35]
const HIGHER_ORDER_CELLS = [60:67...]
const HIGHER_ORDER_PARAMETRIC_CELLS = [51:56...]
const SPECIAL_CELLS = [41, 42]


text:
using PyPlot
using OffsetArrays

function bary_grid(n, a)
    x = OffsetArray{Float64}(undef, 0:n, 0:n)
    y = OffsetArray{Float64}(undef, 0:n, 0:n)
    pts = OffsetArray(range(0, 1, length=n+1), 0:n)
    for i = 0:n
        ξ₁ = pts[i]
        for j = 0:n-i
            ξ₂ = pts[j]
            ξ₃ = 1 - ξ₁ - ξ₂
            x[i,j] = ξ₁ * a[1,1] + ξ₂ * a[1,2] + ξ₃ * a[1,3]
            y[i,j] = ξ₁ * a[2,1] + ξ₂ * a[2,2] + ξ₃ * a[2,3]
        end
        for j = n-i+1:n
            x[i,j] = NaN
            y[i,j] = NaN
        end
    end
    return x, y
end

function bary_func(f, n, a)
    z = OffsetArray{Float64}(undef, 0:n, 0:n)
    pts = OffsetArray(range(0, 1, length=n+1), 0:n)
    ξ = Vector{Float64}(undef, 3)
    for i = 0:n
        ξ[1] = pts[i]
        for j = 0:n-i
            ξ[2] = pts[j]
            ξ[3] = 1 - ξ[1] - ξ[2]
            z[i,j] = f(ξ)
        end
    end
    return z
end

function bary_plot(x, y, z)
    n = size(x, 1) - 1
    xd = [ x[i,n-i] for i = 0:n ]
    yd = [ y[i,n-i] for i = 0:n ]
    zd = [ z[i,n-i] for i = 0:n ]
    idx = [1, 2, 3, 1]
    plot3D(a[1,idx], a[2,idx], zeros(4), "k")
    for i = 0:n
        plot3D(x[i,0:n-i], y[i,0:n-i], z[i,0:n-i], "c")
    end
    for j = 0:n
        plot3D(x[0:n-j,j], y[0:n-j,j], z[0:n-j,j], "c")
    end
    for i = 0:n
        plot3D(xd, yd, zd, "c")
    end
    axis("off")
#    xlabel(L"$x$")
#    ylabel(L"$y$")
end

function midpoints(a)
    prev = [ 3, 1, 2 ]
    succ = [ 2, 3, 1 ]
    m = Matrix{Float64}(undef, 2, 3)
    for i = 1:3
        m[:,i] = ( a[:,prev[i]] + a[:,succ[i]] ) / 2
    end
    return m
end

ψ₁(ξ) = 2 * ξ[1] * ( ξ[1] - 0.5 )
ψ₂(ξ) = 2 * ξ[2] * ( ξ[2] - 0.5 )
ψ₃(ξ) = 2 * ξ[3] * ( ξ[3] - 0.5 )
ψ₄(ξ) = 4 * ξ[2] * ξ[3]
ψ₅(ξ) = 4 * ξ[3] * ξ[1]
ψ₆(ξ) = 4 * ξ[1] * ξ[2]

a = Float64[ 0  1  2
             2  0  3 ]
m = midpoints(a)
n = 20
x, y = bary_grid(n, a)

figure(1)
z = bary_func(ψ₁, n, a)
bary_plot(x, y, z)
plot3D(a[1,1], a[2,1], "ro")
savefig("psi1.pdf")

figure(2)
z = bary_func(ψ₂, n, a)
bary_plot(x, y, z)
plot3D(a[1,2], a[2,2], "ro")
savefig("psi2.pdf")

figure(3)
z = bary_func(ψ₃, n, a)
bary_plot(x, y, z)
plot3D(a[1,3], a[2,3], "ro")
savefig("psi3.pdf")

figure(4)
z = bary_func(ψ₄, n, a)
bary_plot(x, y, z)
plot3D(m[1,1], m[2,1], "ro")
savefig("psi4.pdf")

figure(5)
z = bary_func(ψ₅, n, a)
bary_plot(x, y, z)
plot3D(m[1,2], m[2,2], "ro")
savefig("psi5.pdf")

figure(6)
z = bary_func(ψ₆, n, a)
bary_plot(x, y, z)
plot3D(m[1,3], m[2,3], "ro")
savefig("psi6.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
# Test the branching of the SoilWaterParameterizations functions,
# test the instantiation of them works as expected, and that they
# return what we expect.
using MPI
using OrderedCollections
using StaticArrays
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations

@testset "Land water parameterizations" begin
    FT = Float64
    test_array = [0.5, 1.0]
    vg_model = vanGenuchten{FT}()
    mm = MoistureDependent{FT}()
    bc_model = BrooksCorey{FT}()
    hk_model = Haverkamp{FT}()

    #Use an array to confirm that extra arguments are unused.
    @test viscosity_factor.(Ref(ConstantViscosity{FT}()), test_array) ≈
          [1.0, 1.0]
    @test impedance_factor.(Ref(NoImpedance{FT}()), test_array, test_array) ≈
          [1.0, 1.0]
    @test moisture_factor.(
        Ref(MoistureIndependent{FT}()),
        Ref(vg_model),
        test_array,
    ) ≈ [1.0, 1.0]

    viscosity_model = TemperatureDependentViscosity{FT}(; T_ref = FT(1.0))
    @test viscosity_factor(viscosity_model, FT(1.0)) == 1
    impedance_model = IceImpedance{FT}(; Ω = 2.0)
    @test impedance_factor(impedance_model, 0.2, 0.2) == FT(0.1)

    @test moisture_factor(mm, vg_model, FT(1)) == 1
    @test moisture_factor(mm, vg_model, FT(0)) == 0

    @test moisture_factor(mm, bc_model, FT(1)) == 1
    @test moisture_factor(mm, bc_model, FT(0)) == 0

    @test moisture_factor(mm, hk_model, FT(1)) == 1
    @test moisture_factor(mm, hk_model, FT(0)) == 0

    @test hydraulic_conductivity(
        impedance_model,
        viscosity_model,
        MoistureDependent{FT}(),
        vanGenuchten{FT}(),
        0.5,
        0.5,
        1.0,
        1.0,
    ) == FT(0.1)

    @test_throws Exception effective_saturation(0.5, -1.0, 0.0)
    @test effective_saturation(0.5, 0.25, 0.0) == 0.5

    test_array = [0.5, 1.0]
    n = FT(1.43)
    m = 1.0 - 1.0 / n
    α = FT(2.6)
    ps = SoilParamFunctions{FT}(
        porosity = 1.0,
        Ksat = 0.0,
        S_s = 0.001,
        θ_r = 0.0,
    )
    @test pressure_head.(Ref(vg_model), Ref(ps), test_array, Ref(0.0)) ≈
          .-((-1 .+ test_array .^ (-1 / m)) .* α^(-n)) .^ (1 / n)
    #test branching in pressure head
    @test pressure_head(vg_model, ps, 1.5, 0.0) == 500

    @test pressure_head.(Ref(hk_model), Ref(ps), test_array, Ref(0.0)) ≈
          .-((-1 .+ test_array .^ (-1 / m)) .* α^(-n)) .^ (1 / n)

    m = FT(0.5)
    ψb = FT(0.1656)
    @test pressure_head(bc_model, ps, 0.5, 0.0) ≈ -ψb * 0.5^(-1 / m)
    @test volumetric_liquid_fraction.([0.5, 1.5], Ref(0.75)) ≈ [0.5, 0.75]

    @test inverse_matric_potential(
        bc_model,
        matric_potential(bc_model, FT(0.5)),
    ) == FT(0.5)
    @test inverse_matric_potential(
        vg_model,
        matric_potential(vg_model, FT(0.5)),
    ) == FT(0.5)
    @test inverse_matric_potential(
        hk_model,
        matric_potential(hk_model, FT(0.5)),
    ) == FT(0.5)

    @test_throws Exception inverse_matric_potential(bc_model, 1)
    @test_throws Exception inverse_matric_potential(hk_model, 1)
    @test_throws Exception inverse_matric_potential(vg_model, 1)

end


text:
# Problematic problems from Convex.jl test files
using Convex

### Test 1
solver = DR(eps=1e-6, checki=100, max_iters=100000)
solver = GAPA(eps=1e-8, max_iters=10000)
solver = GAPA(eps=1e-6, max_iters=100000, checki=100, direct=true)
solver = LongstepWrapper(GAPA(eps=1e-6, max_iters=100000, checki=100, direct=true), longinterval=900, nsave=4)
solver = AP(eps=1e-8, max_iters=10000, checki=10, direct=true)

x = Variable(Positive())
y = Variable((3, 3))
p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1)

x = Variable(Positive())

solve!(p, solver)
p.optval

solver = GAPA(eps=1e-8, max_iters=100000, checki=100, direct=true)
solver = DR(eps=1e-8, max_iters=100000, checki=10, direct=true)
solver = LongstepWrapper(GAPA(eps=1e-8, max_iters=10000, checki=10, direct=true), longinterval=50, nsave=2)
solver = LongstepWrapper(GAPA(eps=1e-8, max_iters=10000, checki=10, direct=true), longinterval=50, nsave=5)

using Random
x = Variable(200, 1)
Random.seed!(1)
A = randn(500,200)
b = randn(500)
p = minimize(norm2(A * x + b))
solve!(p, solver)
p.optval

v0 = p.optval


text:
using Knet: relu

function build(m::Model, dtrn::Data)
    chain = iscategorical(m) ? CategoricalChain : LinearChain
    inputsize = size(dtrn.x, 1)
    outputsize = iscategorical(m) ? length(unique(dtrn.y)) : size(dtrn.y, 1)

    hidden = m.params["hidden"]
    pdrop = m.params["pdrop"]

    layer = LinearLayer(inputsize, hidden, relu; pdrop=pdrop)
    layer2 = LinearLayer(hidden, outputsize, identity; pdrop=pdrop)

    m.model = chain(layer, layer2)
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
using MacroTools: @capture, postwalk

using ..Basics: ⊕, projections

# If no upper or lower bound for a summation index can be determined, this cutoff is used,
# i.e., the range of summation indices is considered a subset of -CUTOFF:1/2:CUTOFF.
const CUTOFF = 30

"""
    @racahsum expression var1[=range1] var2[=range2] …

Evaluate the sum over `expression` with the summation indices `var1`, `var2`, etc. The
summand `expression` must contain at least one call to one of the functions
[`clebschgordan`](@ref), [`wigner3j`](@ref), [`wigner6j`](@ref), and [`wigner9j`](@ref). The
bounds for each of the specified summation indices are determined from the selection rules
for the arguments of these functions.

The summation index may be restricted to a range of values (which may be reduced further by
applying the selection rules). The type of such a range should be a subtype of
`OrdinalRange{<:Union{HalfInt,Integer}}`, e.g. `UnitRange{Int}` or
`StepRange{HalfInt,HalfInt}`.

# Examples

A typical call to this macro might look like this:
```jldoctest
julia> @racahsum (2j₃+1)*wigner6j(1,3/2,j₃,0,5/2,3/2)*wigner6j(1,3/2,j₃,0,5/2,3/2) j₃
0.24999999999999986
```
Here, we have verified the orthogonality relation for the Wigner 6-j symbol by summing over
the index `j₃`, which according to the selection rules can only take the value `5/2`, so the
sum consists of only one summand.

To verify the orthogonality relation of Wigner 9-j symbols, two summation indices are needed:
```jldoctest
julia> @racahsum (2j₇+1)*(2j₈+1)*wigner9j(1/2,1/2,1,3,3,2,j₇,j₈,1)*wigner9j(1/2,1/2,1,3,3,2,j₇,j₈,1) j₇ j₈
0.06666666666666662
```

If no upper bound can be placed on a summation index, it is set to $CUTOFF (unless a range is
specified explicitly):
```jldoctest
julia> @racahsum wigner6j(j,0,j,0,j,0) j       # j ∈ 0:1/2:$CUTOFF
0.6686547347577827

julia> @racahsum wigner6j(j,0,j,0,j,0) j=1:10  # j ∈ 0:10
3.2630912526125706
```
"""
macro racahsum(ex::Expr, sumindices...)
    # The following arrays hold all information about the summation indices:
    # * sumvars: the names of the variables
    # * sumtypes: the types of the variables, determined from an explicit type annotation or
    #   from their position in the Wigner symbols or Clebsch–Gordan coefficients. Possible
    #   values are :AngularL, :AngularJ, :AngularM, and :undef (if no type has been
    #   determined yet). The types are only used to restrict the possible range of each
    #   summation index (i.e., an AngularM may be an arbitrary half-integer, an AngularJ
    #   must be a non-negative half-integer, and an AngularL must be a non-negative
    #   integer), the actual type of a summation index is always HalfInt. If a symbol
    #   appears in places for different types (e.g., in an AngularM field of a 3-j symbol
    #   and in a 6-j symbol), the most restrictive type is chosen (AngularL ⊂ AngularJ ⊂
    #   AngularM).
    # * sumranges: expressions that evaluate to ranges which restrict the summation indices,
    #   e.g. `:(projections(j))` or `:(j₁ ⊕ j₂)`.
    # * sumvalues: expressions that evaluate to values that the summation index must be
    #   equal to, e.g. `:(-m₁-m₂)` for `m₃` in a Wigner 3-j symbol.
    # * sumlbounds: expressions that evaluate lower bounds on summation indices, i.e.,
    #   `:(abs(m₁))` implies that the summation index may only take values abs(m₁):1:∞.
    # Any expression that restricts a summation index (in sumranges, sumvalues, or
    # sumlbounds) may only use the values of outermore summation indices (i.e., the ones
    # that were declared earlier in the list).
    sumvars  = Vector{Symbol}(undef, length(sumindices))
    sumtypes = Vector{Symbol}(undef, length(sumindices))
    sumranges  = [[] for _ = sumindices]
    sumvalues  = [[] for _ = sumindices]
    sumlbounds = [[] for _ = sumindices]
    # Parse list of summation indices (in reverse order, so innermost index comes first)
    # with optional type annotation and range
    for (i,sumvar) in enumerate(Iterators.reverse(sumindices))
        var, typ, range = _parse_sumvar(sumvar)
        var in sumvars[1:i-1] && error("summation index $var declared multiple times.")
        sumvars[i] = var
        sumtypes[i] = something(typ, :undef)
        if range !== nothing
            _in_sumvars(range, sumvars[1:i]) &&
                error("summation range contains innermore summation index: $sumvar.")
            push!(sumranges[i], esc(range))
        end
    end
    # Analyse “Racah expression” to determine bounds for each summation index
    postwalk(ex) do x
        if @capture(x, clebschgordan(args__))
            _clebschgordanranges!(sumranges, sumtypes, sumvalues, sumlbounds, sumvars, args)
        elseif @capture(x, wigner3j(args__))
            _wigner3jranges!(sumranges, sumtypes, sumvalues, sumlbounds, sumvars, args)
        elseif @capture(x, wigner6j(args__))
            _wigner6jranges!(sumranges, sumtypes, sumvars, args)
        elseif @capture(x, wigner9j(args__))
            _wigner9jranges!(sumranges, sumtypes, sumvars, args)
        end
        x
    end
    # Assemble loops and return expression. The summation index that was declared last
    # (i.e., the *first* variable in sumvars) defines the innermost for loop or let block.
    loops = :(sum += $(esc(ex)))
    for (var,typ,ranges,values,lbounds) in zip(sumvars,sumtypes,sumranges,sumvalues,sumlbounds)
        typ === :undef && error("$var does not appear in any Wigner symbol.")
        loops = _addloop(loops, var, typ, ranges, values, lbounds)
    end
    :(local sum = 0.0; $loops; sum)
end

# Return a for loop with the loop body `body`, looping over the variable `var` which is
# restricted by the type `typ` (:AngularJ, :AngularL, or :AngularM), the ranges `ranges` and
# the lower bounds `lbounds`. If `values` is not empty, a let block is created instead of a
# loop, setting `var` to the first value in `values` and executing `body` only if `var` is
# equal to all other elements of `values` and is contained in each range specified by
# `ranges` and `lbounds`.
function _addloop(body, var, typ, ranges, values, lbounds)
    allranges = copy(ranges)
    _append_lbounds!(allranges, lbounds, isempty(ranges) & isempty(values))
    if isempty(values)
        _append_typerange!(allranges, typ, !isempty(lbounds) || _hasoplus(ranges))
        if isempty(allranges)
            error("failed to determine summation bounds for $var.")
        elseif length(allranges) == 1
            ex = :(for $(esc(var)) = StepRange{HalfInt,HalfInt}($(allranges[1])); $body; end)
        else
            ex = :(for $(esc(var)) = StepRange{HalfInt,HalfInt}(∩($(allranges...))); $body; end)
        end
    else
        conds = [:($(esc(var)) == $val) for val = values[2:end]]
        if isempty(lbounds) && !_hasoplus(ranges) && typ ∈ (:AngularJ,:AngularL)
            push!(conds, :($(esc(var)) ≥ 0))
        end
        if typ === :AngularL
            push!(conds, :(isinteger($(esc(var)))))
        end
        append!(conds, [:($(esc(var)) ∈ $range) for range = allranges])
        condition = :(true)
        for cond in Iterators.reverse(conds)
            condition = :($cond && $condition)
        end
        ex = :(let $(esc(var)) = HalfInt($(values[1])); if $condition; $body; end; end)
    end
    ex
end

# Append lbounds to ranges, using lbound:CUTOFF or lbound:typemax(HalfInt)
# depending on the value of usecutoff
function _append_lbounds!(ranges, lbounds, usecutoff)
    for lbound in lbounds
        if usecutoff
            push!(ranges, :($lbound:$CUTOFF))
        else
            push!(ranges, :($lbound:$(typemax(HalfInt))))
        end
    end
end

# Append boundaries based on the type (:AngularJ, :AngularL, or :AngularM) to ranges.
# If `isnonnegative` is true, the ranges are assumed to already exclude negative numbers.
function _append_typerange!(ranges, typ, isnonnegative)
    if typ === :AngularJ
        if isempty(ranges)
            push!(ranges, 0:HalfInt(1,2):CUTOFF)
        elseif !isnonnegative
            push!(ranges, 0:HalfInt(1,2):typemax(HalfInt))
        end
    elseif typ === :AngularL
        isempty(ranges) ? push!(ranges, 0:CUTOFF) : push!(ranges, 0:typemax(HalfInt))
    elseif typ === :AngularM
        isempty(ranges) && push!(ranges, -CUTOFF:onehalf(HalfInt):CUTOFF)
    else
        error("invalid type: $typ.")
    end
end

# Check whether any of the expressions in `ranges` is an angular momentum addition `a ⊕ b`,
# where `a` and `b` may be any other object.
_hasoplus(ranges) = any(x -> @capture(x, a_ ⊕ b_), ranges)

# For an expression `var::Type=value` where `::Type` and `=value` are optional, return the
# tuple (var, Type, value). If `Type` and/or `value` are not given, `nothing` is returned in
# their place.
function _parse_sumvar(expr)
    if @capture(expr, var_=range_)
        _extract_sumvartype(var)..., range
    else
        _extract_sumvartype(expr)..., nothing
    end
end
_parse_sumvar(x::Symbol) = x, nothing, nothing

# For an expression `var::Type` where `::Type` is optional, return the tuple (var, Type). If
# `Type` is not given, `nothing` is returned in its place.
function _extract_sumvartype(expr)
    if @capture(expr, var_::T_) && var isa Symbol && T isa Symbol
        T in (:AngularL, :AngularJ, :AngularM) ||
            error("invalid type annotation in summation index: $expr.")
        return var, T
    end
    error("invalid summation index expression: $expr.")
end
_extract_sumvartype(var::Symbol) = var, nothing

# Return true iff the expression `a` contains any of the symbols in `otherjs`
_in_sumvars(a, otherjs) = a ∈ otherjs
_in_sumvars(a::Expr, otherjs) =
    _in_sumvars(a.head, otherjs) || any(x -> _in_sumvars(x, otherjs), a.args)

# For a clebschgordan call with the arguments `args`, extract the restrictions that can be
# placed on each variable in `sumvars` due to selection rules and append the resulting
# ranges/types/values/lbounds to the appropriate arrays. Only restrictions that do not
# depend on innermore summation indices (i.e., symbols that appear earlier in `sumvars`) or
# the summation index itself are considered.
function _clebschgordanranges!(sumranges, sumtypes, sumvalues, sumlbounds, sumvars, args)
    _check_argcount(:clebschgordan, 6, args)
    for ind = 1:2:5
        j = args[ind]
        m = args[ind+1]
        jaindex, jbindex = mod1(ind+2, 6), mod1(ind+4, 6)
        jindex = findfirst(isequal(j), sumvars)
        if jindex !== nothing
            _setvartype!(sumtypes, jindex, :AngularJ)
            _addjrange!(sumranges[jindex], sumvars[1:jindex], args[jaindex], args[jbindex])
            _addlbound!(sumlbounds[jindex], sumvars[1:jindex], m)
        end
    end
    j₁, m₁, j₂, m₂, J, M = args
    if @capture(m₁, -mm₁_)
        if (mindex=findfirst(isequal(mm₁),sumvars); mindex !== nothing)
            _setvartype!(sumtypes, mindex, :AngularM)
            _addmrange!(sumranges[mindex], sumvars[1:mindex], j₁)
            _addvalue_pm!(sumvalues[mindex], sumvars[1:mindex], m₂, M)
        end
    elseif (mindex=findfirst(isequal(m₁),sumvars); mindex !== nothing)
        _setvartype!(sumtypes, mindex, :AngularM)
        _addmrange!(sumranges[mindex], sumvars[1:mindex], j₁)
        _addvalue_pm!(sumvalues[mindex], sumvars[1:mindex], M, m₂)
    end
    if @capture(m₂, -mm₂_)
        if (mindex=findfirst(isequal(mm₂),sumvars); mindex !== nothing)
            _setvartype!(sumtypes, mindex, :AngularM)
            _addmrange!(sumranges[mindex], sumvars[1:mindex], j₂)
            _addvalue_pm!(sumvalues[mindex], sumvars[1:mindex], m₁, M)
        end
    elseif (mindex=findfirst(isequal(m₂),sumvars); mindex !== nothing)
        _setvartype!(sumtypes, mindex, :AngularM)
        _addmrange!(sumranges[mindex], sumvars[1:mindex], j₂)
        _addvalue_pm!(sumvalues[mindex], sumvars[1:mindex], M, m₁)
    end
    if @capture(M, -mM_)
        if (mindex=findfirst(isequal(mM),sumvars); mindex !== nothing)
            _setvartype!(sumtypes, mindex, :AngularM)
            _addmrange!(sumranges[mindex], sumvars[1:mindex], J)
            _addvalue_mm!(sumvalues[mindex], sumvars[1:mindex], m₁, m₂)
        end
    elseif (mindex=findfirst(isequal(M),sumvars); mindex !== nothing)
        _setvartype!(sumtypes, mindex, :AngularM)
        _addmrange!(sumranges[mindex], sumvars[1:mindex], J)
        _addvalue_pp!(sumvalues[mindex], sumvars[1:mindex], m₁, m₂)
    end
end

# Like _clebschgordanranges!, but for wigner3j
function _wigner3jranges!(sumranges, sumtypes, sumvalues, sumlbounds, sumvars, args)
    _check_argcount(:wigner3j, 6, args)
    for ind = 1:3
        j = args[ind]
        m = args[ind+3]
        jaindex, jbindex = mod1(ind+1, 3), mod1(ind+2, 3)
        maindex, mbindex = jaindex+3, jbindex+3
        if (jindex=findfirst(isequal(j),sumvars); jindex !== nothing)
            _setvartype!(sumtypes, jindex, :AngularJ)
            _addjrange!(sumranges[jindex], sumvars[1:jindex], args[jaindex], args[jbindex])
            _addlbound!(sumlbounds[jindex], sumvars[1:jindex], m)
        end
        if @capture(m, -mm_)
            if (mindex=findfirst(isequal(mm),sumvars); mindex !== nothing)
                _setvartype!(sumtypes, mindex, :AngularM)
                _addmrange!(sumranges[mindex], sumvars[1:mindex], j)
                _addvalue_pp!(sumvalues[mindex], sumvars[1:mindex], args[maindex], args[mbindex])
            end
        elseif (mindex=findfirst(isequal(m),sumvars); mindex !== nothing)
            _setvartype!(sumtypes, mindex, :AngularM)
            _addmrange!(sumranges[mindex], sumvars[1:mindex], j)
            _addvalue_mm!(sumvalues[mindex], sumvars[1:mindex], args[maindex], args[mbindex])
        end
    end
end

# For a wigner6j call with the arguments `args`, extract the restrictions that can be placed
# on each variable in `sumvars` due to selection rules and append the resulting ranges/types
# to the appropriate arrays. Only restrictions that do not depend on innermore summation
# indices (i.e., symbols that appear earlier in `sumvars`) or the summation index itself are
# considered.
function _wigner6jranges!(sumranges, sumtypes, sumvars, js)
    _check_argcount(:wigner6j, 6, js)
    for (ind,j) = enumerate(js)
        jindex = findfirst(isequal(j), sumvars)
        jindex === nothing && continue
        _setvartype!(sumtypes, jindex, :AngularJ)
        div, rem = divrem(ind-1, 3)
        aindex, bindex = mod1(rem+2, 3), mod1(rem, 3)
        cindex, dindex = aindex+3, bindex+3
        if iszero(div) # j is in the top row of the 6-j symbol
            _addjrange!(sumranges[jindex], sumvars[1:jindex], js[aindex], js[bindex])
            _addjrange!(sumranges[jindex], sumvars[1:jindex], js[cindex], js[dindex])
        else # j is in the bottom row of the 6-j symbol
            _addjrange!(sumranges[jindex], sumvars[1:jindex], js[aindex], js[dindex])
            _addjrange!(sumranges[jindex], sumvars[1:jindex], js[cindex], js[bindex])
        end
    end
end

# Same as _wigner6jranges! but for wigner9j
function _wigner9jranges!(sumranges, sumtypes, sumvars, js)
    _check_argcount(:wigner9j, 9, js)
    for (ind,j) = enumerate(js)
        jindex = findfirst(isequal(j), sumvars)
        jindex === nothing && continue
        _setvartype!(sumtypes, jindex, :AngularJ)
        div, rem = divrem(ind-1, 3)
        aindex = 3div + mod1(rem+2, 3)
        bindex = 3div + mod1(rem, 3)
        cindex = 3mod(div+1, 3) + rem + 1
        dindex = 3mod(div+2, 3) + rem + 1
        _addjrange!(sumranges[jindex], sumvars[1:jindex], js[aindex], js[bindex])
        _addjrange!(sumranges[jindex], sumvars[1:jindex], js[cindex], js[dindex])
    end
end

function _check_argcount(func::Symbol, expected, args)
    length(args) == expected || error("wrong number of arguments for $func: " *
                                      "expected $expected, got $(length(args)).")
end

# Set `types[index] = newtype`, but only if `newtype` is more restrictive than the type that
# is currently stored in `types[index]`.
function _setvartype!(types, index, newtype)
    order = (:undef, :AngularM, :AngularJ, :AngularL)
    oldindex = findfirst(isequal(types[index]), order)
    newindex = findfirst(isequal(newtype), order)
    oldindex === nothing && error("invalid type annotation: $(types[index]).")
    newindex === nothing && error("invalid type annotation: $newtype.")
    types[index] = order[max(oldindex,newindex)]
end

# Add the range `j₁ ⊕ j₂` to `jranges`, unless it is already in there or contains
# symbols from `otherjs`.
function _addjrange!(jranges, otherjs, j₁, j₂)
    if !_in_sumvars(j₁, otherjs) && !_in_sumvars(j₂, otherjs) &&
        (:($(esc(j₁)) ⊕ $(esc(j₂))) ∉ jranges) && (:($(esc(j₂)) ⊕ $(esc(j₁))) ∉ jranges)
        push!(jranges, :($(esc(j₁)) ⊕ $(esc(j₂))))
    end
end

# Add the lower bound `abs(m)` to `jlbounds`, unless it is already in there or contains
# symbols from `otherjs`.
function _addlbound!(jlbounds, otherjs, m)
    if !_in_sumvars(m, otherjs) && (:(abs(HalfInt($(esc(m))))) ∉ jlbounds)
        push!(jlbounds, :(abs(HalfInt($(esc(m))))))
    end
end

# Add the range `projections(j)` (i.e., `-j:j`) to `mranges`, unless it is already in there
# or contains symbols from `otherjs`.
function _addmrange!(mranges, otherjs, j)
    if !_in_sumvars(j, otherjs) && (:(projections($(esc(j)))) ∉ mranges)
        push!(mranges, :(projections($(esc(j)))))
    end
end

# Add the value `-m₁-m₂` to `mvalues`, unless it is already in there or contains symbols
# from `otherjs`.
function _addvalue_mm!(mvalues, otherjs, m₁, m₂)
    if !_in_sumvars(m₁, otherjs) && !_in_sumvars(m₂, otherjs) &&
        (:(-$(esc(m₁))-$(esc(m₂))) ∉ mvalues) && (:(-$(esc(m₂))-$(esc(m₁))) ∉ mvalues)
        push!(mvalues, :(-$(esc(m₁))-$(esc(m₂))))
    end
end

# Add the value `M-m` to `mvalues`, unless it is already in there or contains symbols
# from `otherjs`.
function _addvalue_pm!(mvalues, otherjs, M, m)
    if !_in_sumvars(M, otherjs) && !_in_sumvars(m, otherjs) && (:($(esc(M))-$(esc(m))) ∉ mvalues)
        push!(mvalues, :($(esc(M))-$(esc(m))))
    end
end

# Add the value `m₁+m₂` to `mvalues`, unless it is already in there or contains symbols
# from `otherjs`.
function _addvalue_pp!(mvalues, otherjs, m₁, m₂)
    if !_in_sumvars(m₁, otherjs) && !_in_sumvars(m₂, otherjs) &&
        (:($(esc(m₁))+$(esc(m₂))) ∉ mvalues) && (:($(esc(m₂))+$(esc(m₁))) ∉ mvalues)
        push!(mvalues, :($(esc(m₁))+$(esc(m₂))))
    end
end


text:
function diffeq!(du, u, p, t)
    du[V.M] = p[C.kM] * p[C.S] / (p[C.KS] + p[C.S]) - p[C.dM] * u[V.M]
    du[V.E] = p[C.kE] * (u[V.M] / (p[C.KM] + u[V.M])) * (u[V.E] / (p[C.KE] + u[V.E])) + p[C.kb] * u[V.M] / (p[C.KM] + u[V.M]) +
                p[C.kP1] * u[V.CD] * u[V.RE] / (p[C.KCD] + u[V.RE]) + p[C.kP2] * u[V.CE] * u[V.RE] / (p[C.KCE] + u[V.RE]) - p[C.dE] * u[V.E] - p[C.kRE] * u[V.R] * u[V.E]
    du[V.CD] = p[C.kCD] * u[V.M] / (p[C.KM] + u[V.M]) + p[C.kCDS] * p[C.S] / (p[C.KS] + p[C.S]) - p[C.dCD] * u[V.CD]
    du[V.CE] = p[C.kCE] * u[V.E] / (p[C.KE] + u[V.E]) - p[C.dCE] * u[V.CE]
    du[V.R] = p[C.kR] + p[C.kDP] * u[V.RP] / (p[C.KRP] + u[V.RP]) - p[C.kRE] * u[V.R] * u[V.E] - p[C.kP1] * u[V.CD] * u[V.R] / (p[C.KCD] + u[V.R]) -
                p[C.kP2] * u[V.CE] * u[V.R] / (p[C.KCE] + u[V.R]) - p[C.dR] * u[V.R]
    du[V.RP] = p[C.kP1] * u[V.CD] * u[V.R] / (p[C.KCD] + u[V.R]) + p[C.kP2] * u[V.CE] * u[V.R] / (p[C.KCE] + u[V.R]) + p[C.kP1] * u[V.CD] * u[V.RE] / (p[C.KCD] + u[V.RE]) +
                p[C.kP2] * u[V.CE] * u[V.RE] / (p[C.KCE] + u[V.RE]) - p[C.kDP] * u[V.RP] / (p[C.KRP] + u[V.RP]) - p[C.dRP] * u[V.RP]
    du[V.RE] = p[C.kRE] * u[V.R] * u[V.E] - p[C.kP1] * u[V.CD] * u[V.RE] / (p[C.KCD] + u[V.RE]) - p[C.kP2] * u[V.CE] * u[V.RE] / (p[C.KCE] + u[V.RE]) - p[C.dRE] * u[V.RE]

end

function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.S] = 2.0
    p[C.kE] = 0.4
    p[C.kM] = 1.0
    p[C.kCD] = 0.03
    p[C.kCDS] = 0.45
    p[C.kR] = 0.18
    p[C.kRE] = 180
    p[C.kb] = 0.003
    p[C.KS] = 0.5
    p[C.kCE] = 0.35
    p[C.dM] = 0.7
    p[C.dE] = 0.25
    p[C.dCD] = 1.5
    p[C.dCE] = 1.5
    p[C.dR] = 0.06
    p[C.dRP] = 0.06
    p[C.dRE] = 0.03
    p[C.kP1] = 18.0
    p[C.kP2] = 18.0
    p[C.kDP] = 3.6
    p[C.KM] = 0.15
    p[C.KE] = 0.15
    p[C.KCD] = 0.92
    p[C.KCE] = 0.92
    p[C.KRP] = 0.01

    return p
end

function get_derivatives(u::Vector{Float64}, p::Vector{Float64})
    # derivatives: dF/d[bifurcation_param]
    dFdp::Vector{Float64} = zeros(V.NUM)

    dFdp[V.M] = p[C.kM] * p[C.KS] / (p[C.KS] + p[C.S])^2
    dFdp[V.CD] = p[C.kCDS] * p[C.KS] / (p[C.KS] + p[C.S])^2

    return dFdp
end

function get_steady_state(p::Vector{Float64})
    tspan::Tuple{Float64,Float64} = (0.0, Inf)
    u0::Vector{Float64} = zeros(V.NUM)

    prob = ODEProblem(diffeq!, u0, tspan, p)
    prob = SteadyStateProblem(prob)
    sol = solve(prob, DynamicSS(CVODE_BDF()), dt=1.0)

    return sol.u
end


text:
function main()
    nums::Array{Int} = [1,2,3,4,5] # 배열 생성
    println("전부: $nums")

    first = nums[1] # 첫번째 요소 접근
    println("첫번째요소: $first")
end

main()

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example