text:
module IrregularInterpolate

using LinearAlgebra

export PolyharmonicSpline, interpolate

struct PolyharmonicSpline{T}
    dim::Int
    order::Int
    coeff::Vector{T}
    centers::Matrix{T}
    error::T
end

function polyharmonicK(r, K)
    if iseven(K)
        return r < 1 ? r^(K-1) * log(r^r) : (r^K) * log(r)
    else
        return r^K
    end
end

function PolyharmonicSpline(K::Int, centers::Matrix{T}, values::Array{T}; s = zero(T)) where {T}
    N, dim = size(centers)
    N != length(values) && throw(DimensionMismatch())

    phK = Matrix{T}(undef, N, N)
    for i in 1:N
        for j in 1:N
            dist = 0.0
            for d in 1:dim
                dist += (centers[i, d] - centers[j, d])^2
            end
            phK[i, j] = polyharmonicK(sqrt(dist), K)
        end
    end

    A = copy(phK)
    B = zeros(T, N, dim+1)
    for i in 1:N
        B[i, 1] = 1
        B[i, 2:end] = centers[i, :]
    end

    A .+= s .* Diagonal(I, N)

    L = [A B; B' zeros(T, dim+1, dim+1)]
    w = pinv(L) * [values; zeros(T, dim+1)]

    ivalues = zeros(T, N)
    for i in 1:N
        tmp = 0.0
        for j in 1:N
            tmp += w[j] * phK[i, j]
        end
        tmp += w[N+1]
        for j in 2:dim+1
            tmp += w[N+j] * centers[i, j-1]
        end
        ivalues[i] = tmp
    end
    error = norm(values .- ivalues)

    return PolyharmonicSpline(dim, K, w, centers, error)
end

function PolyharmonicSpline(K::Int, centers::Vector{T}, values::Vector{T}; s = zero(T)) where {T}
    # PolyharmonicSpline(K, centers'', values, s = s)
    PolyharmonicSpline(K, centers, values, s = s)
end

function interpolate(S::PolyharmonicSpline{T}, x::Matrix{T}) where {T}
    m, n = size(x)
    n != S.dim && throw(DimensionMismatch("$m != $(S.dim)"))

    l = length(S.coeff) - (n+1)

    interpolates = zeros(T, m)
    for i in 1:m
        tmp = 0.0
        for j in 1:l
            dist = 0.0
            for d in 1:n
                dist += (x[i, d] - S.centers[j, d])^2
            end
            tmp += S.coeff[j] * polyharmonicK(sqrt(dist), S.order)
        end
        tmp += S.coeff[l+1]
        for j in 2:n+1
            tmp += S.coeff[l+j] * x[i, j-1]
        end
        interpolates[i] = tmp
    end
    return interpolates
end

function interpolate(S::PolyharmonicSpline{T}, x::Vector{T}) where {T}
    return interpolate(S, x)
end

function interpolate(S::PolyharmonicSpline{T}, x::Vector{T}, y::Vector{T}) where {T}
    return interpolate(S, [x y])
end

function interpolate(S::PolyharmonicSpline{T}, x::Vector{T}, y::Vector{T}, z::Vector{T}) where {T}
    return interpolate(S, [x y z])
end

end  # module


text:
import ProjectAssigner
ProjectAssigner.match(students="students.csv", projects="projects.csv", output="output.csv")

ext: jl


text:
"""
A HMC implementation where the kinetic energy is motiviated by relativistic kinetic energy following

X. Lu, V. Perrone, L. Hasenclever, Y. W. Teh, Sebastian J. Vollmer relativistic Monte Carlo https://arxiv.org/abs/1609.04388

Fields:
    - x: position
    - p: momentum
    - xi: scalar factor for the thermostat
    - ntiers: number of iteration
    - mass: diagonal mass matrix
    - c: speed of light
    - D: diagonal preconditioner
    - independent_momenta: true means kinetic energy is applied to each component of momenta
                           false mean kinetic energy is applied to euclidean length of momenta
Constructors:
    - construct using `RelSGHMCState(x::Array{Float64,1},p::Array{Float64,1},stepsize::Float64; mass=ones(length(x)),c=[1.0], D=[1.0], independent_momenta=true)`

"""

type RelSGHMCState <: SamplerState
    x::Array{Float64,1}
    p::Array{Float64,1}

    t::Int64
    stepsize::Function
    mass::Array{Float64,1}
    c::Array{Float64,1}
    D::Array{Float64,1}
    independent_momenta::Bool
    """
    Constructs a RelSGHMCState given an initial state, momentum and a fixed stepsize
    """
    function RelSGHMCState(x::Array{Float64,1},p::Array{Float64,1},stepsize::Float64; mass=ones(length(x)),c=[1.0], D=[1.0], independent_momenta=true)
        # momentum, mass and state need to have the same size.
        @assert size(x) == size(p)
        @assert size(mass) == size(x)
        # D matrix should be scalar or have the same size as x
        @assert size(x) == size(D) || length(D) == 1
        @assert size(x) == size(c) || length(c) == 1
        # stepsize is used as a function
        s(niters) = stepsize

    new(x,p,0,s,mass,c,D,independent_momenta)
    end
    """
    Constructs a RelSGHMCState given an initial state, momentum and a stepsize schedule
    """
    function RelSGHMCState(x::Array{Float64,1},p::Array{Float64,1},stepsize::Function; mass=ones(length(x)),c=[1.0], D=[1.0],  independent_momenta=true)
        # momentum, mass and state need to have the same size.
        @assert size(x) == size(p)
        @assert size(mass) == size(x)
        # D matrix should be scalar or have the same size as x
        @assert size(x) == size(D) || length(D) == 1
        @assert size(x) == size(c) || length(c) == 1
        new(x,p,0,stepsize,mass,c,D,independent_momenta)
    end
    """
    Constructs a RelSGHMCState given an initial state and a fixed stepsize
    """
    function RelSGHMCState(x::Array{Float64,1},stepsize::Float64; mass=ones(length(x)),c=[1.0], D=[1.0], independent_momenta=true)
        # momentum, mass and state need to have the same size.
        @assert size(mass) == size(x)
        # D matrix should be scalar or have the same size as x
        @assert size(x) == size(D) || length(D) == 1
        @assert size(x) == size(c) || length(c) == 1
        if independent_momenta
            p=zeros(length(x))
            for i=1:length(x) p[i]=sample_rel_p(mass,c,1)[1]  end
        else
            p = sample_rel_p(mass, c, length(x))
        end
        # stepsize is used as a function
        s(niters) = stepsize
    new(x,p,0,s,mass,c,D,independent_momenta)
    end
    """
    Constructs a RelSGHMCState given an initial state and a stepsize schedule
    """
    function RelSGHMCState(x::Array{Float64,1},stepsize::Function; mass=ones(length(x)),c=[1.0], D=[1.0], independent_momenta=true)
        @assert size(mass) == size(x)
        # D matrix should be scalar or have the same size as x
        @assert size(x) == size(D) || length(D) == 1
        @assert size(x) == size(c) || length(c) == 1
        if independent_momenta
            p=zeros(length(x))
            for i=1:length(x) p[i]=sample_rel_p(mass,c,1)[1]  end
        else
            p = sample_rel_p(mass, c, length(x))
        end
        new(x,p,0,stepsize,mass,c,D,independent_momenta)
    end
end

function sample!(s::RelSGHMCState, grad)
    #increase iteration counter
    s.t += 1

    # improve readability
    x = s.x
    p = s.p
    D = s.D
    t = s.t
    stepsize = s.stepsize(t)
    m = s.mass
    c = s.c
    independent_momenta=s.independent_momenta

    tmp = zeros(length(x))
    if independent_momenta
        tmp = m .* sqrt(p.^2 ./ (m.^2 .* c.^2) + 1)
    else
        tmp = m .* sqrt(p'p ./ (m.^2 .* c.^2) + 1)
    end

    n = sqrt(stepsize.*(2*D)).*randn(length(x)) # noise term

    #update momentum
    p[:] += stepsize.*(grad(x)-D.*p./tmp) + n
    #update tmp

    if independent_momenta
        tmp = m .* sqrt(p.^2 ./ (m.^2 .* c.^2) + 1)
    else
        tmp = m .* sqrt(p'p ./ (m.^2 .* c.^2) + 1)
    end
    #update state
    x[:] += stepsize.*p./tmp

    s
end


text:
using Elemental, TSVD, LinearAlgebra, Test
using Random
A = Elemental.DistMatrix(Float64)
Elemental.gaussian!(A, 500, 200)

# Warning! If seed is not set, each worker will end up with a different initial vector
# making the Lanczos processes inconsistent
Random.seed!(123)
vals_Elemental = tsvd(A, 5)[2]
localA = Array(A)
vals_LAPACK    = svdvals(localA)[1:5]
@test vals_Elemental ≈ vals_LAPACK


text:
function compute(obj::DOPComputer, arg0::AbsoluteDate, arg1::List)
    return jcall(obj, "compute", DOP, (AbsoluteDate, List), arg0, arg1)
end

function create(::Type{DOPComputer}, arg0::OneAxisEllipsoid, arg1::GeodeticPoint)
    return jcall(DOPComputer, "create", DOPComputer, (OneAxisEllipsoid, GeodeticPoint), arg0, arg1)
end

function get_elevation_mask(obj::DOPComputer)
    return jcall(obj, "getElevationMask", ElevationMask, ())
end

function get_min_elevation(obj::DOPComputer)
    return jcall(obj, "getMinElevation", jdouble, ())
end

function with_elevation_mask(obj::DOPComputer, arg0::ElevationMask)
    return jcall(obj, "withElevationMask", DOPComputer, (ElevationMask,), arg0)
end

function with_min_elevation(obj::DOPComputer, arg0::jdouble)
    return jcall(obj, "withMinElevation", DOPComputer, (jdouble,), arg0)
end


text:
const base_unary_real = (
    (Base.acos, (-1, 1)),
    (Base.asin, (-1, 1)),
    (Base.atan, (-50, 50)),
    (Base.cos, (-1000, 1000)),
    (Base.sin, (-1000, 1000)),
    (Base.tan, (-1000, 1000)),
    (Base.acosh, (1, 1000)),
    (Base.asinh, (-1000, 1000)),
    (Base.atanh, (-1, 1)),
    (Base.cosh, (0, 89.415985f0)),
    (Base.sinh, (-89.415985f0, 89.415985f0)),
    (Base.tanh, (-8.66434f0, 8.66434f0)),
    (Base.cbrt, (-1000, 1000)),
    (Base.sqrt, (0, 1000)),
    (Base.exp, (-88.72284f0, 88.72284f0)),
    (Base.expm1, (-88.72284f0, 88.72284f0)),
    (Base.log, (0, 1000)),
    (Base.log10, (0, 1000)),
    (Base.log1p, (-1, 1000)),
    (Base.abs, (-1000, 1000)),
    (Base.abs2, (-1000, 1000)),
    (Base.ceil, (-1000, 1000)),
    (Base.floor, (-1000, 1000)),
    (Base.round, (-1000, 1000)),
    (Base.trunc, (-1000, 1000)),
    (Base.erf, (-3.8325067f0, 3.8325067f0)),
    (Base.erfc, (-3.7439213f0, 10.019834f0)),
    (Base.erfinv, (-1, 1)),
    (Base.erfcinv, (0, 2)),
    (Base.lgamma, (0, 1000)),
    (Base.gamma, (0, 36))
)

const base_binary_real = (
    (Base.atan2, (-1, 1), (-1, 1)),
    (Base.hypot, (-1000, 1000), (-1000, 1000)),
    (Base.(:./), (-1000, 1000), (-1000, 1000)),
    (Base.(:.^), (0, 100), (-5, 20))
)

const base_unary_complex = (
    (Base.acos, (-1, 1)),
    (Base.asin, (-1, 1)),
    # (Base.atan, (-50, 50)),
    # (Base.cos, (-10, 10)),
    # (Base.sin, (-10, 10)),
    # (Base.tan, (-10, 10)),
    (Base.acosh, (1, 1000)),
    (Base.asinh, (-1000, 1000)),
    # (Base.atanh, (-1, 1)),
    # (Base.cosh, (0, 89.415985f0)),
    # (Base.sinh, (-89.415985f0, 89.415985f0)),
    # (Base.tanh, (-8.66434f0, 8.66434f0)),
    (Base.sqrt, (0, 1000)),
    (Base.exp, (-88.72284f0, 88.72284f0)),
    (Base.log, (0, 1000)),
    # (Base.log10, (0, 1000)),
    (Base.abs, (-1000, 1000)),
    (Base.angle, (-1000, 1000))
    # (Base.conj, (-1000, 1000))
)

const base_binary_complex = (
    (Base.(:./), (-1000, 1000), (-1000, 1000)),
    # (Base.(:.^), (0, 100), (-2, 10))
)

function randindomain{T<:Real}(t::Type{T}, n, domain)
    d1 = oftype(t, domain[1])
    d2 = oftype(t, domain[2])
    ddiff = d2 - d1
    @assert isfinite(ddiff)
    v = rand(t, n)
    for i = 1:length(v)
        v[i] = v[i]*ddiff+d1
    end
    v
end

function randindomain{T<:Complex}(t::Type{T}, n, domain)
    d1 = oftype(t, domain[1])
    d2 = oftype(t, domain[2])
    ddiff = d2 - d1
    @assert isfinite(ddiff)
    v = rand(t, 2*n)
    for i = 1:length(v)
        v[i] = v[i]*ddiff+d1
    end
    reinterpret(t, v)
end


text:
module ConsistencyResampling

using StatsBase: StatsBase
using Random: Random

export Consistent

include("alias.jl")
include("consistent.jl")
include("many.jl")
include("single.jl")

end # module


text:
######################################################################
# ergodic_manager.jl
######################################################################
abstract type ErgodicManager end

decompose!(em::ErgodicManager) = decompose!(em, em.phi)

"""
`phi = reconstruct(em)`

Reconstructs distribution from Fourier coefficients `em.phik`.

`phi = reconstruct(em, ck::Matrix{Float64})`

Reconstructs distribution from Fourier coefficients `ck`.

Both options return a distribution `phi`, a `em.bins` by `em.bins` matrix.
The value at `phi[xi,yi]` gives the value at the cell index `xi,yi`.

`phi = reconstruct(em, xd::VVF)`

Reconstructs from trajectory `xd`.
"""
reconstruct(em::ErgodicManager) = reconstruct(em, em.phik)
function reconstruct(em::ErgodicManager, xd::VVF)
	reconstruct(em, decompose(em, xd))
end

x_min(em::ErgodicManager) = em.domain.mins[1]
y_min(em::ErgodicManager) = em.domain.mins[2]
z_min(em::ErgodicManager) = em.domain.mins[3]

x_max(em::ErgodicManager) = em.domain.maxes[1]
y_max(em::ErgodicManager) = em.domain.maxes[2]
z_max(em::ErgodicManager) = em.domain.maxes[3]

x_size(em::ErgodicManager) = em.domain.cell_lengths[1]
y_size(em::ErgodicManager) = em.domain.cell_lengths[2]
z_size(em::ErgodicManager) = em.domain.cell_lengths[3]

x_cells(em::ErgodicManager) = em.domain.cells[1]
y_cells(em::ErgodicManager) = em.domain.cells[2]
z_cells(em::ErgodicManager) = em.domain.cells[3]

# decomposition on a group of trajectories
function decompose(em::ErgodicManager, xds::VVVF)
	num_agents = length(xds)

	ck = decompose(em, xds[1])
	for i = 2:num_agents
		ck += decompose(em, xds[i])
	end

	# divide by the number of agents (that is, length(xds))
	return ck / num_agents
end

"""
`ergodic_score(em, traj::VVF)`

First breaks down the trajectory into components ck.
"""
function ergodic_score(em::ErgodicManager, traj::VVF)
	ck = decompose(em, traj)
	return ergodic_score(em, ck)
end

function ergodic_score(em::ErgodicManager, ck)
	val = 0.0
	for (i, L) in enumerate(em.Lambda)
		d = em.phik[i] - ck[i]
		val += L * d * d
	end
	return val
end

normalize!(em::ErgodicManager) = normalize!(em.phi, em.domain.cell_size)


text:
function render(list::Tuple)
    for elem in list
        render(elem)
    end
    return
end
"""
When rendering a specialised list of Renderables, we can do some optimizations
"""
function render(list::Vector{RenderObject{Pre}}) where Pre
    isempty(list) && return nothing
    first(list).prerenderfunction()
    vertexarray = first(list).vertexarray
    program = vertexarray.program
    glUseProgram(program.id)
    glBindVertexArray(vertexarray.id)
    for renderobject in list
        Bool(to_value(renderobject.uniforms[:visible])) || continue # skip invisible
        # make sure we only bind new programs and vertexarray when it is actually
        # different from the previous one
        if renderobject.vertexarray != vertexarray
            vertexarray = renderobject.vertexarray
            if vertexarray.program != program
                program = renderobject.vertexarray.program
                glUseProgram(program.id)
            end
            glBindVertexArray(vertexarray.id)
        end
        for (key, value) in program.uniformloc
            if haskey(renderobject.uniforms, key)
                if length(value) == 1
                    gluniform(value[1], renderobject.uniforms[key])
                elseif length(value) == 2
                    gluniform(value[1], value[2], renderobject.uniforms[key])
                else
                    error("Uniform tuple too long: $(length(value))")
                end
            end
        end
        renderobject.postrenderfunction()
    end
    # we need to assume, that we're done here, which is why
    # we need to bind VertexArray to 0.
    # Otherwise, every glBind(::GLBuffer) operation will be recorded into the state
    # of the currently bound vertexarray
    glBindVertexArray(0)
    return
end

"""
Renders a RenderObject
Note, that this function is not optimized at all!
It uses dictionaries and doesn't care about OpenGL call optimizations.
So rewriting this function could get us a lot of performance for scenes with
a lot of objects.
"""
function render(renderobject::RenderObject, vertexarray=renderobject.vertexarray)
    if Bool(to_value(renderobject.uniforms[:visible]))
        renderobject.prerenderfunction()
        program = vertexarray.program
        glUseProgram(program.id)
        for (key, value) in program.uniformloc
            if haskey(renderobject.uniforms, key)
                # uniform_name_type(program, value[1])
                if length(value) == 1
                    gluniform(value[1], renderobject.uniforms[key])
                elseif length(value) == 2
                    gluniform(value[1], value[2], renderobject.uniforms[key])
                else
                    error("Uniform tuple too long: $(length(value))")
                end
            end
        end
        glBindVertexArray(vertexarray.id)
        renderobject.postrenderfunction()
        glBindVertexArray(0)
    end
    return
end

"""
Renders a vertexarray, which consists of the usual buffers plus a vector of
unitranges which defines the segments of the buffers to be rendered
"""
function render(vao::GLVertexArray{T}, mode::GLenum=GL_TRIANGLES) where T <: VecOrSignal{UnitRange{Int}}
    for elem in to_value(vao.indices)
        glDrawArrays(mode, max(first(elem) - 1, 0), length(elem) + 1)
    end
    return nothing
end

function render(vao::GLVertexArray{T}, mode::GLenum=GL_TRIANGLES) where T <: TOrSignal{UnitRange{Int}}
    r = to_value(vao.indices)
    offset = first(r) - 1 # 1 based -> 0 based
    ndraw = length(r)
    nverts = length(vao)
    if (offset < 0 || offset + ndraw > nverts)
        error("Bounds error for drawrange. Offset $(offset) and length $(ndraw) aren't a valid range for vertexarray with length $(nverts)")
    end
    glDrawArrays(mode, offset, ndraw)
    return nothing
end

function render(vao::GLVertexArray{T}, mode::GLenum=GL_TRIANGLES) where T <: TOrSignal{Int}
    r = to_value(vao.indices)
    glDrawArrays(mode, 0, r)
    return nothing
end

"""
Renders a vertex array which supplies an indexbuffer
"""
function render(vao::GLVertexArray{GLBuffer{T}}, mode::GLenum=GL_TRIANGLES) where T <: Union{Integer,AbstractFace}
    glDrawElements(
        mode,
        length(vao.indices) * cardinality(vao.indices),
        julia2glenum(T), C_NULL
    )
    return
end

"""
Renders a normal vertex array only containing the usual buffers buffers.
"""
function render(vao::GLVertexArray, mode::GLenum=GL_TRIANGLES)
    glDrawArrays(mode, 0, length(vao))
    return
end

"""
Render instanced geometry
"""
renderinstanced(vao::GLVertexArray, a, primitive=GL_TRIANGLES) = renderinstanced(vao, length(a), primitive)

"""
Renders `amount` instances of an indexed geometry
"""
function renderinstanced(vao::GLVertexArray{GLBuffer{T}}, amount::Integer, primitive=GL_TRIANGLES) where T <: Union{Integer,AbstractFace}
    glDrawElementsInstanced(primitive, length(vao.indices) * cardinality(vao.indices), julia2glenum(T), C_NULL, amount)
    return
end

"""
Renders `amount` instances of an not indexed geoemtry geometry
"""
function renderinstanced(vao::GLVertexArray, amount::Integer, primitive=GL_TRIANGLES)
    glDrawElementsInstanced(primitive, length(vao), GL_UNSIGNED_INT, C_NULL, amount)
    return
end
# handle all uniform objects

##############################################################################################
#  Generic render functions
#####
function enabletransparency()
    glEnablei(GL_BLEND, 0)
    glDisablei(GL_BLEND, 1)
    # This does:
    # target.rgb = source.a * source.rgb + (1 - source.a) * target.rgb
    # target.a = 0 * source.a + 1 * target.a
    # the latter is required to keep target.a = 1 for the OIT pass
    glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ZERO, GL_ONE)
    return
end


text:
import .Basic: tofloat

# old utility

tofloat(::Type{F}, o::OneHotArray) where {F<:AbstractFloat} = CuArray{F}(o)

# avoid scalar operation by broadcasting assigment
function CUDA.CuArray{F}(o::OneHotArray{K, N, var"N+1", A}) where {F <: AbstractFloat, K, N, var"N+1", A <: CuArray}
  dest = similar(o, F)
  dest .= o
  return dest
end

Base.BroadcastStyle(::Type{<: OneHotArray{K, N, var"N+1", A}}) where {K, N, var"N+1", A <: CuArray} = CUDA.CuArrayStyle{var"N+1"}()


text:
include("constructed_images.jl")
include("readremote.jl")

FactCheck.exitstatus()


text:
import Base: view

export MatFacModel, save_model, load_model

mutable struct MatFacModel
    X::AbstractMatrix
    Y::AbstractMatrix
    row_transform::Any  # These transforms may be functions
    col_transform::Any  # or callable structs (with trainable parameters)
    noise_model::CompositeNoise

    X_reg::Any                   # These regularizers may be functions
    Y_reg::Any                   # or callable structs (with trainable
    row_transform_reg::Any       # parameters)
    col_transform_reg::Any
    noise_model_reg::Any
end

@functor MatFacModel

function MatFacModel(X::AbstractMatrix, Y::AbstractMatrix,
                     noise_model::CompositeNoise;
                     row_transform=x->x,
                     col_transform=x->x,
                     X_reg=x->0.0, Y_reg=x->0.0,
                     row_transform_reg=x->0.0,
                     col_transform_reg=x->0.0,
                     noise_model_reg=x->0.0)

    row_transform = make_viewable(row_transform)
    col_transform = make_viewable(col_transform)

    return MatFacModel(X, Y,
                       row_transform,
                       col_transform,
                       noise_model,
                       X_reg, Y_reg,
                       row_transform_reg,
                       col_transform_reg,
                       noise_model_reg)

end

"""
    MatFacModel(M::Int, N::Int, K::Int, col_losses::Vector{String};
                row_transform=identity,
                col_transform=identity,
                X_reg=x->0, Y_reg=y->0,
                row_transform_reg=x->0,
                col_transform_reg=x->0,
                noise_model_reg=x->0)

Construct a matrix factorization model for an
M x N dataset. Assume K latent factors.

Specify the loss for each column via `col_losses`,
a length-N vector of strings. Permissible string values are:

* "normal"
* "bernoulli"
* "poisson"
* "ordinalN", where "N" is replaced by some integer.
  E.g., "ordinal3", "ordinal4".

All of the keyword arguments are either (a) functions or
(b) callable structs (i.e., functors). We allow for the possibility
for each of these to contain trainable parameters.

* `row_transform`: a (trainable) transformation on the _rows_ of the model.
* `col_transform`: a (trainable) transformation on the _columns_ of the model.
* `*_reg`: a (trainable) regularizer for the specified object.

"""
function MatFacModel(M::Integer, N::Integer, K::Integer,
                     col_losses::Vector{String}; kwargs...)

    X = randn(K,M) ./ (sqrt(K)*100)
    Y = randn(K,N) ./ 100

    noise_model = CompositeNoise(col_losses)

    return MatFacModel(X, Y, noise_model; kwargs...)
end

"""
    MatFacModel(M, N, K, loss::String; kwargs...)

Identical to the other constructor, but applies the same
`loss` function to all columns of data.
"""
function MatFacModel(M::Integer, N::Integer, K::Integer,
                     loss::String; kwargs...)
    losses = fill(loss, N)
    return MatFacModel(M, N, K, losses; kwargs)
end

function (bm::MatFacModel)()
    return invlink(bm.noise_model,
            bm.col_transform(
             bm.row_transform(
               transpose(bm.X)*bm.Y
             )
            )
           )
end

function view(bm::MatFacModel, idx1, idx2)
    return MatFacModel(view(bm.X, :, idx1),
                            view(bm.Y, :, idx2),
                            view(bm.row_transform, idx1, idx2),
                            view(bm.col_transform, idx1, idx2),
                            view(bm.noise_model, idx2),
                            nothing, nothing, nothing,
                            nothing, nothing,
                           )

end

function Base.size(bm::MatFacModel)
    return (size(bm.X, 2), size(bm.Y, 2))
end

##############################################
# Equality operation
EqTypes = Union{CompositeNoise,NormalNoise,PoissonNoise,BernoulliNoise,OrdinalNoise,
                MatFacModel,ViewableFunction}

NoEqTypes = Function

function Base.:(==)(a::T, b::T) where T <: EqTypes
    for fn in fieldnames(T)
        af = getfield(a, fn)
        bf = getfield(b, fn)
        if !(af == bf)
            if !((typeof(af) <: NoEqTypes) & (typeof(bf) <: NoEqTypes))
                return false
            end
        end
    end
    return true
end

################################################
# Model file I/O

"""
    save_model(filename, model)

Save `model` to a BSON file located at `filename`.
"""
function save_model(filename, model)
    BSON.@save filename model
end

"""
    load_model(filename)

load a model from the BSON located at `filename`.
"""
function load_model(filename)
    d = BSON.load(filename, @__MODULE__)
    return d[:model]
end


text:
function filter_noise(N::Int, N_mask::Float64, Map::Array{Float64, 2})

    inds = ((0:(N-1)) .+ 0.5 .- (N/2))
    X = reshape(repeat(inds, N), N ,N)
    Y = X'
    R = ((X .^2 .+ Y .^2).^0.5)

    mask = ones(N, N)
    mask[abs.(Y) .< N_mask] .= 0

    FMap = fftshift(ifft(fftshift(Map)))
    FMap_filtered = FMap .* mask
    Map_filtered = real.(fftshift(fft(fftshift(FMap_filtered))))

    return Map_filtered

end


text:
#
using diffMesh, Plots

ng = 10;
np = 10*ng;

xp = Array(range(-1,1,length=np))'

pou = intervalPOU(ng,x -> @. x^2)

plt = plot()
plt = plot!(xp',pou(xp)',width=3)


text:
using ExampleNeuralNetworks
using Base.Test

import ExampleNeuralNetworks.RectifiedLinearUnit
import ExampleNeuralNetworks.forward

include("feedforwardunit.jl")

max_stars_repo_name: helgith/ExampleNeuralNetworks

max_issues_repo_name: helgith/ExampleNeuralNetworks

max_forks_repo_name: helgith/ExampleNeuralNetworks


text:
using FFTW

export QFTCircuit, QFTBlock, invorder_firstdim

export openbox

"""
    openbox(block::AbstractBlock) -> AbstractBlock

For a black box, like QFTBlock, you can get its white box (loyal simulation) using this function.
"""
function openbox end

CRk(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<<k)))
CRot(n::Int, i::Int) = chain(n, i==j ? kron(i=>H) : CRk(j, i, j-i+1) for j = i:n)
QFTCircuit(n::Int) = chain(n, CRot(n, i) for i = 1:n)

struct QFTBlock{N} <: PrimitiveBlock{N} end
YaoBlocks.mat(::Type{T}, q::QFTBlock{N}) where {T, N} = T.(applymatrix(q))

YaoBlocks.apply!(reg::DefaultRegister{B}, ::QFTBlock) where B = (reg.state = ifft!(invorder_firstdim(reg |> state), 1)*sqrt(1<<nactive(reg)); reg)
YaoBlocks.apply!(reg::DefaultRegister{B}, ::Daggered{N, <:QFTBlock}) where {B,N} = (reg.state = invorder_firstdim(fft!(reg|>state, 1)/sqrt(1<<nactive(reg))); reg)

# traits
YaoBlocks.ishermitian(q::QFTBlock{N}) where N = N==1
YaoBlocks.isreflexive(q::QFTBlock{N}) where N = N==1
YaoBlocks.isunitary(q::QFTBlock{N}) where N = true

openbox(q::QFTBlock{N}) where N = QFTCircuit(N)
openbox(q::Daggered{<:QFTBlock, N}) where {N} = QFTCircuit(N)'

function print_block(io::IO, pb::QFTBlock{N}) where N
    printstyled(io, "QFT(1-$N)"; bold=true, color=:blue)
end

function print_block(io::IO, pb::Daggered{N,<:QFTBlock}) where {N, T}
    printstyled(io, "IQFT(1-$N)"; bold=true, color=:blue)
end

function invorder_firstdim(v::Matrix)
    w = similar(v)
    n = size(v, 1) |> log2i
    n_2 = n ÷ 2
    mask = [bmask(i, n-i+1) for i in 1:n_2]
    @simd for b in basis(n)
        @inbounds w[breflect(b, mask; nbits=n)+1,:] = v[b+1,:]
    end
    w
end

function invorder_firstdim(v::Vector)
    n = length(v) |> log2i
    n_2 = n ÷ 2
    w = similar(v)
    #mask = SVector{n_2, Int}([bmask(i, n-i+1)::Int for i in 1:n_2])
    mask = [bmask(i, n-i+1)::Int for i in 1:n_2]
    @simd for b in basis(n)
        @inbounds w[breflect(b, mask; nbits=n)+1] = v[b+1]
    end
    w
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot, ColorSchemes, Colors
x = -2π:0.001:2π
function circShape(h,k,r)
    θ = LinRange(0,2π,500)
    h .+ r*sin.(θ), k .+ r*cos.(θ)
end
bgcp1 = "set object rectangle from screen 0,0 to screen 1,1"
bgcp2 = " behind fillcolor rgb 'black' fillstyle solid noborder"
bgcolor = bgcp1*bgcp2
cmap = get(colorschemes[:viridis], LinRange(0,1,15))
@gp circShape(0,0,1)... "w l notit lw 2 lc '#$(hex(cmap[3]))'"
@gp :- "set auto fix" "set size square"
@gp :- "set offsets graph .05, graph .05, graph .05, graph .05"
@gp :- "set border lw 1 lc rgb 'white'"
@gp :- "set ylabel 'y' textcolor rgb 'white'"
@gp :- "set xlabel 'x' textcolor rgb 'white'"
@gp :- "set xzeroaxis linetype 3 linewidth 1"
@gp :- "set yzeroaxis linetype 3 linewidth 1"
@gp :- bgcolor
for (indx, r) in enumerate(0.9:-0.1:0.1)
    @gp :- circShape(0,0,r)... "w l notit lw 2 lc '#$(hex(cmap[indx+3]))'" :-
end
@gp
save(term="pngcairo size 400,400 ", output="./code/plt2_ex7.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#!/usr/bin/env julia
using Printf
using Distributed
using TimerOutputs
using CUDAdrv
using CUDAnative
using CuArrays

function ζ_gen_thr_range(mstart,mend,e)
    s = 0
    cls  = 8  # cache line size in units of Float64 storage size
              # so we can avoid false sharing
    Nthr = Threads.nthreads()
    psum = zeros(Float64,Nthr*cls)
    ms   = zeros(Float64,Nthr*cls)
    me   = zeros(Float64,Nthr*cls)

    # preparing each thread, start, end, and partial sums
    for i=1:Nthr
        ms[i*cls]   = (N/Nthr)*i
        me[i*cls]   = 1 + (N/Nthr)*(i-1)
        psum[i*cls] = 0.0
    end

    # actual function we will execute per thread
    # psums will be updated per thread id.  This is why
    # the psums are spaced 1 cache line apart, versus adjacent indices

    # Threaded inner loop, each thread has no dependence upon others
    Threads.@threads for t=1:Nthr
        innercpu_1!(psum,ms,me,2,cls)
    end

    # now sum up the partial sums and return s
    s = sum(psum)
    s
end

function ζ_2_sqr_thr_range(mstart,mend)
    s = 0
    cls  = 8  # cache line size in units of Float64 storage size
              # so we can avoid false sharing
    Nthr = Threads.nthreads()
    psum = zeros(Float64,Nthr*cls)
    ms   = zeros(Float64,Nthr*cls)
    me   = zeros(Float64,Nthr*cls)

    # preparing each thread, start, end, and partial sums
    for i=1:Nthr
        ms[i*cls]   = (N/Nthr)*i
        me[i*cls]   = 1 + (N/Nthr)*(i-1)
        psum[i*cls] = 0.0
    end

    # actual function we will execute per thread
    # psums will be updated per thread id.  This is why
    # the psums are spaced 1 cache line apart, versus adjacent indices

    # Threaded inner loop, each thread has no dependence upon others
    Threads.@threads for t=1:Nthr
        innercpu_2!(psum,ms,me,2,cls)
    end

    # now sum up the partial sums and return s
    s = sum(psum)
    s
end

function ζ_2_sqr_gpu_range(N)
    #dev = CuDevice(0)
    #ctx = CuContext(dev)
    #MAX_THREADS_PER_BLOCK = CUDAdrv.attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
	MAX_THREADS_PER_BLOCK = 1024
    Nthr    = min(N, MAX_THREADS_PER_BLOCK)
    Nblocks = ceil.(Int, N / Nthr)

    psum_d = CuArrays.fill(0.0::Float64, Nthr)
	p  = zeros(Float64,Nthr)
	@printf("Nblocks = %i\nNthr = %i\n",Nblocks,Nthr)

    @cuda blocks=Nblocks threads=Nthr  innergpu_1!(psum_d)
    copyto!(p,psum_d)
    s = sum(p)
    s
end

function innercpu_1!(psum,ms,me,x,cls)
	ndx     = Threads.threadid()
	@simd for i=ms[ndx*cls]:-1:me[ndx*cls]
		r = 1.0/i
		@inbounds psum[ndx*cls] += r^x
	end
nothing
end

function innercpu_2!(psum,ms,me,x,cls)
	ndx     = Threads.threadid()
	@simd for i=ms[ndx*cls]:-1:me[ndx*cls]
		@inbounds psum[ndx*cls] += (1.0/i)*(1.0/i)
	end
	return nothing
end

function innergpu_1!(ps)
	ndx     = threadIdx().x
	i = (blockIdx().x-1) * blockDim().x + threadIdx().x
	r = 1.0/i
	@atomic ps[ndx] += r*r
	return nothing
end

function ζ_2_sqr_simd_range(mstart,mend)
    s = 0.0
    @simd for i=mstart:-1:mend
      s += (1.0/i)*(1.0/i)
    end
    s
end

to = TimerOutput()

N = 16000000000
#N = 16000000
#Nthr = 1024
st = 1
s =  zeros(Float64,32)
pd = zeros(Float64,32)
ndx = 1
#ENV["JULIA_DEBUG"] = "CUDAnative"

# first bit of code to handle the compilation stage, otherwise you
# see this in the time output

s[1] = ζ_2_sqr_simd_range(N,1)       #
s[2] = ζ_gen_thr_range(N,1,2)        #
s[3] = ζ_2_sqr_thr_range(N,1)        #
s[4] = ζ_2_sqr_gpu_range(N)  #

for i=1:4
  @printf("s[%i] = %f\n",i,s[i])
end

@timeit to "ζ_gen_thr_range"       pd[ndx] = ζ_gen_thr_range(N,1,2)
ndx+=1

@timeit to "ζ_2_sqr_thr_range"     pd[ndx] = ζ_2_sqr_thr_range(N,1)
ndx+=1

@timeit to "ζ_2_sqr_simd_range"    pd[ndx] = ζ_2_sqr_simd_range(N,1)
ndx+=1

@timeit to "ζ_2_sqr_gpu_range"     pd[ndx] = ζ_2_sqr_gpu_range(N)
ndx+=1

for i=1:4
  @printf("pd[%i] = %f\n",i,pd[i])
end

show(to)
@printf("\n");

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
module Connections

using Interpolations
using CSV, DataFrames

export cfs_rot_screwfastened_k, cfs_trans_screwfastened_k, cfs_pull_through_plate_stiffness, cfs_load_deformation_interpolation_model

#Cold-formed steel
#Rotational stiffness for a screw-fastened connection
#Distributed spring!
#https://doi.org/10.1016/j.tws.2012.06.005
function cfs_rot_screwfastened_k(b, c, S, t, kp, E, CorZ)

    I = 1/12*S*t^3

    if CorZ == 0   #C

        kϕ = (S/(c^2*kp)+(((b^2*c/2)+c^2*b+(c^3/3))*S/(c^2*E*I)))^-1

    elseif CorZ == 1  #Z

        kϕ = (S/(c^2*kp) + c*S/(3*E*I))^-1

    end

    return kϕ

end

#Cold-formed steel
#Translational stiffness for a screw-fastened connection
#Discrete spring!
#https://www.buildusingsteel.org/aisi-design-resources/-/media/doc/buildusingsteel/research-reports/CFSD%20-%20Report%20-%20RP17-2.pdf
function cfs_trans_screwfastened_k(t1, t2, E1, E2, Fss, Fu1, Fu2, D)

    Ka = (1/(E1*t1)+1/(E2*t2))^-1

    ψ=(Fss/(t1*D*Fu1))*(Fss/(t2*D*Fu2))

    #monotonic elastic, steel to steel
    α = 0.27
    β = - 0.69

    Ke = α*ψ^β*Ka

    return Ka, ψ, α, β, Ke

end

"""
        cfs_pull_through_plate_stiffness(x, c, tw)

`x` is the fastener distance from an R-panel major deck rib, in mm.
`c` is the fastener distance from the purlin rotation point, in mm.
'tw' is the R-panel deck base metal thickness, in mm.

The fastener pull-through elastic stiffness 'kp' is interpolated and output in N/mm.

These results come from elastic finite element studies in
https://www.buildusingsteel.org/aisi-design-resources/-/media/doc/buildusingsteel/research-reports/CFSD%20-%20Report%20-%20RP17-2.pdf

"""

function cfs_pull_through_plate_stiffness(x, c, tw)

    #Define source file for pullover stiffness values.
    filename = string(@__DIR__, "/assets/kp.csv")

    #Import the parameters and stiffness values.

    data = CSV.File(filename)
    data = DataFrame(x = Vector(data.x), c = Vector(data.c), tw = Vector(data.tw), kp=Vector(data.kp))
    # data = CSV.read(filename, DataFrame, header=true)

    #Define the range of screw locations relative to an R-panel deck rib.
    x_range = sort(unique(data.x))

    #Define the range of fastener locations from the purlin pivot point.
    c_range = sort(unique(data.c))

    #Define the range of deck panel thicknesses considered.
    t_range = sort(unique(data.tw))

    #Check if the screw location in the R-panel is within the tested range.

    if (tw > maximum(t_range)) | (tw < minimum(t_range))
        print("Beware!  The deck base metal thickness, tw, is outside the tested range.")
    end

    if (c > maximum(c_range)) | (c < minimum(c_range))
        print("Beware!  The fastener distance from the purlin pivot point, c, is outside the tested range.")
    end

    if (x > maximum(x_range)) | (x < minimum(x_range))
        print("Beware!  The fastener distance from the R-panel major rib, x, is outside the tested range.")
    end

    #Define the number of parameters that will be used to predict kp.

    x_dim = length(x_range)
    c_dim = length(c_range)
    t_dim = length(t_range)

    #Initialize the parameter array.
    A = zeros(Float64, (x_dim,c_dim, t_dim))

    #Fill the parameter array with kp.
    for i = 1:x_dim
        for j = 1:c_dim
            for k = 1:t_dim
                x_entry = x_range[i]
                c_entry = c_range[j]
                t_entry = t_range[k]

                index_x = findall(==(x_entry), data.x)
                index_c = findall(==(c_entry), data.c)
                index_t = findall(==(t_entry), data.tw)

                index = intersect(index_x, index_c)
                index = intersect(index, index_t)

                A[i, j, k] = data.kp[index[1]]

            end
        end
    end

    #Define the interpolation operator.  Only linear interpolation is available for multi-dimensional parameter sets in Interpolations.jl.
    #Allow extrapolation with warnings above.
    stiffness_model = LinearInterpolation((x_range, c_range, t_range), A, extrapolation_bc = Line())

    kp = stiffness_model(x, c, tw)

    return kp

end

function cfs_load_deformation_interpolation_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)

    #Define source file for experimental data used to define the interpolation model.
    filename = string(@__DIR__, "/assets/Connections/cfs_fastener_test_data.csv")

    #Load experimental data from the CSV file.
    data = CSV.File(filename)

    #Define sorted, unique ranges of each parameter.
    t1_range = sort(unique(data.t1))
    t2_range = sort(unique(data.t2))
    fy1_range = sort(unique(data.fy1))
    fy2_range = sort(unique(data.fy2))
    fu1_range = sort(unique(data.fu1))
    fu2_range = sort(unique(data.fu2))
    screw_diameter_range = sort(unique(data.screw_diameter))

    #Define interpolation operators for all the deformation points (D1, D2, D3, D4) and the peak load (P3).
    D1_operator = define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, data.D1)
    D2_operator = define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, data.D2)
    D3_operator = define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, data.D3)
    D4_operator = define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, data.D4)
    P3_operator = define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, data.P3)

    #Define interpolation models for all the deformation points (D1, D2, D3, D4) and the peak load (P3).
    D1_model =  LinearInterpolation((t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range), D1_operator, extrapolation_bc = Line())
    D2_model =  LinearInterpolation((t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range), D2_operator, extrapolation_bc = Line())
    D3_model =  LinearInterpolation((t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range), D3_operator, extrapolation_bc = Line())
    D4_model =  LinearInterpolation((t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range), D4_operator, extrapolation_bc = Line())
    P3_model =  LinearInterpolation((t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range), P3_operator, extrapolation_bc = Line())

    #Perform interpolation predictions based on user-defined values.
    D1 = D1_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)
    D2 = D2_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)
    D3 = D3_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)
    D4 = D4_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)
    P3 = P3_model(t1, t2, fy1, fy2, fu1, fu2, screw_diameter)

    #Calculate P1, P2, and P4.
    P1 = 0.4 * P3
    P2 = 0.8 * P3
    P4 = 0.1 * P3

    #Warn user if the interpolation is out of range.
    interpolation_warning(t1, t2, fy1, fy2, fu1, fu2, screw_diameter, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range)

    return D1, D2, D3, D4, P1, P2, P3, P4

end

function define_interpolation_operator(data, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range, D)

    t1_dim = length(t1_range)
    t2_dim = length(t2_range)
    fy1_dim = length(fy1_range)
    fy2_dim = length(fy2_range)
    fu1_dim = length(fu1_range)
    fu2_dim = length(fu2_range)
    screw_diameter_dim = length(screw_diameter_range)

    #Initialize the parameter array.
    A = zeros(Float64, (t1_dim, t2_dim, fy1_dim, fy2_dim, fu1_dim, fu2_dim, screw_diameter_dim))

    #Fill the parameter array with kp.
    for i = 1:t1_dim
        for j = 1:t2_dim
            for k = 1:fy1_dim
                for n = 1:fy2_dim
                    for m = 1:fu1_dim
                        for h = 1:fu2_dim
                            for s  = 1:screw_diameter_dim
                                t1_entry = t1_range[i]
                                t2_entry = t2_range[j]
                                fy1_entry = fy1_range[k]
                                fy2_entry = fy2_range[n]
                                fu1_entry = fu1_range[m]
                                fu2_entry = fu2_range[h]
                                screw_diameter_entry = screw_diameter_range[s]

                                index_t1 = findall(==(t1_entry), data.t1)
                                index_t2 = findall(==(t2_entry), data.t2)
                                index_fy1 = findall(==(fy1_entry), data.fy1)
                                index_fy2 = findall(==(fy2_entry), data.fy2)
                                index_fu1 = findall(==(fu1_entry), data.fu1)
                                index_fu2 = findall(==(fu2_entry), data.fu2)
                                index_screw_diameter = findall(==(screw_diameter_entry), data.screw_diameter)

                                index = intersect(index_t1, index_t2)
                                index = intersect(index, index_fy1)
                                index = intersect(index, index_fy2)
                                index = intersect(index, index_fu1)
                                index = intersect(index, index_fu2)
                                index = intersect(index, index_screw_diameter)

                                if size(index,1) == 0
                                    A[i, j, k, n, m, h, s] = 0
                                else
                                    A[i, j, k, n, m, h, s] = D[index[1]]
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return A
end

function interpolation_warning(t1, t2, fy1, fy2, fu1, fu2, screw_diameter, t1_range, t2_range, fy1_range, fy2_range, fu1_range, fu2_range, screw_diameter_range)

    if (t1 > maximum(t1_range))
        println("Beware!  The thickness of ply 1, t1, is outside the tested range, too large.")
    elseif (t1 < minimum(t1_range))
        println("Beware!  The thickness of ply 1, t1, is outside the tested range, too small.")
    end

    if (t2 > maximum(t2_range))
        println("Beware!  The thickness of ply 2, t2, is outside the tested range, too large.")
    elseif (t2 < minimum(t2_range))
        println("Beware!  The thickness of ply 2, t2, is outside the tested range. too small.")
    end

    if (fy1 > maximum(fy1_range))
        println("Beware!  The yielding strength of ply 1, fy1, is outside the tested range, too large.")
    elseif (fy1 < minimum(fy1_range))
        println("Beware!  The yielding strength of ply 1, fy1, is outside the tested range, too small.")
    end

    if (fy2 > maximum(fy2_range))
        println("Beware!  The yielding strength of ply 2, fy2, is outside the tested range, too large.")
    elseif (fy2 < minimum(fy2_range))
        println("Beware!  The yielding strength of ply 2, fy2, is outside the tested range, too small.")
    end

    if (fu1 > maximum(fu1_range))
        println("Beware!  The ultimate strength of ply 1, fu1, is outside the tested range, too large.")
    elseif (fu1 < minimum(fu1_range))
        println("Beware!  The ultimate strength of ply 1, fu1, is outside the tested range, too small.")
    end

    if (fu2 > maximum(fu2_range))
        println("Beware!  The ultimate strength of ply 2, fu2, is outside the tested range, too large.")
    elseif (fu2 < minimum(fu2_range))
        println("Beware!  The ultimate strength of ply 2, fu2, is outside the tested range, too small.")
    end

    if (screw_diameter > maximum(screw_diameter_range))
        println("Beware!  The screw diameter, screw_diameter, is outside the tested range, too large.")
    elseif (screw_diameter < minimum(screw_diameter_range))
        println("Beware!  The screw diameter, screw_diameter, is outside the tested range, too small.")
    end

end

end #module


text:
using DataFrames, FixedEffectModels
N = 10000000
K = 100
id1 = rand(1:(N/K), N)
id2 = rand(1:K, N)
x1 =  randn(N)
x2 =  randn(N)
w = cos(id1)
y= 3 .* x1 .+ 5 .* x2 .+ cos(id1) .+ cos(id2).^2 .+ randn(N)
df = DataFrame(id1 = pool(id1), id2 = pool(id2), x1 = x1, x2 = x2, w = w, y = y)
@time reg(y ~ x1 + x2, df)
# 1.258554 seconds (723 allocations: 1.205 GB, 18.70% gc time)
#852MB to obtain the matrix etc, and then 1.205 for Regressions part
@time reg(y ~ x1 + x2, df, VcovCluster(:id2))
#  1.569679 seconds (843 allocations: 1.293 GB, 25.21% gc time)
@time reg(y ~ x1 + x2 |> id1, df)
# 1.476390 seconds (890 allocations: 1.175 GB, 20.15% gc time)
@time reg(y ~ x1 + x2 |> id1, df, VcovCluster(:id1))
# 1.974738 seconds (1.04 k allocations: 1.255 GB, 15.85% gc time)
@time reg(y ~ x1 + x2 |> id1 + id2, df)
# 4.554836 seconds (1.01 k allocations: 1.188 GB, 9.56% gc time)
@time reg(y ~ x1 + x2 |> id1, df, weight = :w)
# 2.000850 seconds (20.00 M allocations: 1.010 GB, 18.82% gc time)
@time reg(y ~ x1 + x2 |> id1 + id2, df, weight = :w)
# 4.118974 seconds (20.00 M allocations: 1.018 GB, 9.60% gc time)


text:
export TparamsST

mutable struct TparamsST{T} <: Tparams{T}  # specific parameters for this Krylov variant
    ζ :: T

    function TparamsST{T}() where T
        ζ = 0.5
        return new{T}(ζ)
    end
end


text:
using Markdown
using Mustache
using Dates
using SQLite
using Morsel
using DataFrames

include("setup.jl")
include("models.jl")

app = Morsel.app()
include("views.jl")
include("routes.jl")

global_context = { "staticroot" => staticroot }

function render(tpl, context)
	Mustache.render(tpl, merge(global_context, context))
end


text:
import SSCant
using Base.Test

function test_amplitude()
    x = [0, 1, 0, -1.]
    v = [1, 0, -1, 0.]
    omega_c = 1
    exp_A = [1, 1, 1, 1.]
    A = SSCant.amplitude(x, v, omega_c=omega_c)
    @test_approx_eq(A, exp_A)
end

function test_phase()
    x = [0, 1, 0, -1.]
    v = [1, 0, -1, 0.]
    omega_c = 1
    exp_phi = [-pi/2, 0, pi/2, pi]
    phi = SSCant.phase(x, v, omega_c=omega_c)
    @test_approx_eq(phi, exp_phi)
end

function test_fit_phase()
    t = 0:3
    phi = [-pi/2, 0, pi/2, pi]
    dphi, df, f0 = SSCant.fit_phase(t, phi)
    @test_approx_eq_eps(dphi, [0., 0, 0, 0], 1e-14)
    @test_approx_eq_eps(df, [0., 0, 0, 0], 1e-14)
    @test_approx_eq(f0, 0.25)
end

function test_get_specified()
    x = 0:0.5:5
    y = map(x -> x^2, x)
    x_specified = 0:5
    y_specified_exp =  map(x -> x^2, x_specified)
    y_specified = SSCant.get_specified(x, y, x_specified)
    @test_approx_eq(y_specified, y_specified_exp)
end

function test_chunk_range()
    t0 = -10.
    tf = 14.
    @test_approx_eq(SSCant.chunk_range(t0, tf, 10.), [-10, 0.0, 10, 14])
end

test_amplitude()
test_phase()
test_fit_phase()
test_get_specified()
test_chunk_range()


text:
# Structs and functions related to the codes_nearest type

export Nearest, find, findmultiple, destroy

mutable struct codes_nearest
end

struct Nearest
    ptr::Ptr{codes_nearest}
end

"""
    Nearest(handle::Message)

Create a Nearest from a Message.
"""
function Nearest(handle::Message)::Nearest
    eref = Ref(Int32(0))
    ptr = ccall((:codes_grib_nearest_new, eccodes), Ptr{codes_nearest},
                (Ptr{codes_handle}, Ref{Cint}), handle.ptr, eref)
    errorcheck(eref[])
    return Nearest(ptr)
end

"""
    Nearest(f::Function, handle::Message)

Create a Nearest from a Message and automatically destroy when finished.

# Example
```Julia
Nearest(handle) do near
    find(near, handle, inlat, inlon)
end
```
"""
function Nearest(f::Function, handle::Message)
    near = Nearest(handle)
    try
        f(near)
    finally
        destroy(near)
    end
end

"""
    find(near::Nearest, handle::Message, inlon::Float64, inlat::Float64; samepoint=true, samegrid=true)

Find the nearest 4 points to the given point sorted by distance.

Setting samepoint and samegrid to true (default) speeds up the calculation but restricts the point
and grid, respectively, from changing from call to call. Distance is in kilometers.

# Example
```Julia
Nearest(handle) do near
    lons, lats, values, distances = find(near, handle, inlon, inlat)
end
```
"""
function find(near::Nearest, handle::Message, inlon, inlat; samepoint=true, samegrid=true)
    flags = samepoint && samegrid ? CODES_NEAREST_SAME_POINT | CODES_NEAREST_SAME_GRID :
            samepoint ? CODES_NEAREST_SAME_POINT :
            samegrid ? CODES_NEAREST_SAME_GRID :
            UInt64(0)
    outlats = Vector{Float64}(undef, 4)
    outlons = similar(outlats)
    values = similar(outlats)
    distances = similar(outlats)
    indexes = Vector{Int32}(undef, 4)
    len = Csize_t(4)
    err = ccall((:codes_grib_nearest_find, eccodes), Cint,
                (Ptr{codes_nearest}, Ptr{codes_handle}, Cdouble, Cdouble, Culong,
                 Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ref{Csize_t}),
                near.ptr, handle.ptr, inlat, inlon, flags, outlats, outlons, values, distances,
                indexes, Ref(len))
    errorcheck(err)
    sorting = sortperm(distances)
    return outlons[sorting], outlats[sorting], values[sorting], distances[sorting]
end

"""
    findmultiple(handle::Message, inlons::Vector{Float64}, inlats::Vector{Float64}; islsm=false)

Find the nearest point to each point given.

Setting lsm to true finds the nearest land point and only works when `handle`
represents a land-sea-mask. Distance is in kilometers.

# Example
```Julia
lons, lats, values, distances = findmultiple(handle, inlons, inlats)
```
"""
function findmultiple(handle::Message, inlons::Vector{Float64}, inlats::Vector{Float64}; islsm=false)
    if length(inlats) != length(inlons)
        throw(DomainError("Length of inlats must be equal to length of inlons"))
    end
    npoints = Clong(length(inlats))
    outlats = similar(inlats)
    outlons = similar(inlats)
    outvals = similar(inlats)
    outdists = similar(inlats)
    outidx = Vector{Int32}(undef, npoints)
    intlsm = Int32(islsm)
    err = ccall((:codes_grib_nearest_find_multiple, eccodes), Cint,
                (Ptr{codes_handle}, Cint, Ref{Cdouble}, Ref{Cdouble}, Clong, Ref{Cdouble},
                 Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cint}),
                handle.ptr, intlsm, inlats, inlons, npoints, outlats, outlons,
                outvals, outdists, outidx)
    errorcheck(err)
    return outlons, outlats, outvals, outdists
end

""" Safely destroy a nearest object. """
function destroy(near::Nearest)::Nothing
    err = ccall((:codes_grib_nearest_delete, eccodes), Cint, (Ptr{codes_nearest},), near.ptr)
    errorcheck(err)
end


text:
using .StringMatch

export align

const AlignCount = Dict()

function align(pdftxtfile::String, xmlfile::String)
    pdcontents = readpdftxt(pdftxtfile)
    pdchars = filter(c -> isa(c,PDText), pdcontents)
    iddict = Dict{String,Int}()
    pdids = map(x -> get!(iddict,x.c,length(iddict)+1), pdchars)

    xmltree = readjats(xmlfile)
    xmltree = xmltree[findfirst(c -> c.name == "body", xmltree.children)]
    tokenize!(xmltree)
    xmlchars = findall(isempty, xmltree)
    xmlids = map(c -> get!(iddict,string(c.name),length(iddict)+1), xmlchars)
    pairs = lcsmatch(pdids, xmlids)

    xml2pdf = Dict()
    for (k,v) in pairs
        xml2pdf[v] = k
        push!(pdchars[k].tags, xmlchars[v].name)
    end
    ranges = Dict()
    i = 1
    bottomup(xmltree) do node
        if isempty(node)
            ranges[node] = i:i
            i += 1
        else
            ranges[node] = start(ranges[node[1]]):last(ranges[node[end]])
        end
    end
    for (node,r) in ranges
        i = start(r)
        j = last(r)
        i == j && continue
        haskey(AlignCount,node.name) || (AlignCount[node.name] = [0,0])
        countd = AlignCount[node.name]
        if haskey(xml2pdf,i) && haskey(xml2pdf,j)
            countd[1] += 1
        else
            countd[2] += 1
        end
    end
end

function tokenize2!(tree::Tree)
    nodes = findall(isempty, tree)
    dict = Dict(n.parent => n.parent for n in nodes)
    for node in keys(dict)
        children = Tree[]
        for c in node.children
            if isempty(c)
                for char in Vector{Char}(c.name)
                    char == ' ' && continue
                    push!(children, Tree(string(char)))
                end
            else
                push!(children, Tree("[$(c.name)]"))
            end
        end
        setchildren!(node, children)
    end
end


text:
#=== memory command queue apis ===#
@ocl_func(clCreateCommandQueueWithProperties, CL_command_queue,
              (CL_context, CL_device_id, CL_queue_properties, Ptr{CL_int}))

#=== memory object apis ===#
@ocl_func(clCreatePipe, CL_mem,
              (CL_context, CL_mem_flags, CL_uint, CL_uint, Ptr{CL_pipe_properties}, CL_int))

@ocl_func(clGetPipeInfo, CL_int,
              (CL_mem, CL_pipe_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== SVM Allocation API ===#
@ocl_func(clSVMAlloc, Ptr{Void},
              (CL_context, CL_svm_mem_flags, Csize_t, CL_uint))

@ocl_func(clSVMFree, Void,
              (CL_context, Ptr{Void}))

#=== sampler apis ===#

@ocl_func(clCreateSamplerWithProperties, CL_sampler,
              (CL_context, Ptr{CL_sampler_properties}, Ptr{CL_int}))

#=== kernel object apis ===#
@ocl_func(clSetKernelArgSVMPointer, CL_int,
              (CL_kernel, CL_uint, Ptr{Void}))

@ocl_func(clSetKernelExecInfo, CL_int,
              (CL_kernel, CL_kernel_exec_info, Csize_t, Ptr{Void}))

#=== Enqueued Commands APIs ===#
@ocl_func(clEnqueueSVMFree, CL_int,
               (CL_command_queue, CL_uint, Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void},
                CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueSVMMemcpy, CL_int,
               (CL_command_queue, CL_bool, Ptr{Void}, Ptr{Void}, Csize_t,
                CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueSVMMemFill, CL_int,
               (CL_command_queue, Ptr{Void}, Ptr{Void}, Csize_t, Csize_t,
                CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueSVMMap, CL_int,
               (CL_command_queue, CL_bool, CL_map_flags, Ptr{Void}, Csize_t,
                CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueSVMUnmap, CL_int,
              (CL_command_queue, Ptr{Void}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

#=== deprecation ===#

# @deprecate clCreateCommandQueue clCreateCommandQueueWithProperties
# @deprecate clCreateSampler clCreateSamplerWithProperties
# @deprecate clEnqueueTask


text:
import QMTK.TN
using Compat.Test


text:
module UOT_plot

using LinearAlgebra, Statistics, PyPlot
using UOT, UOT_toolbox

function plot_points_sources(bounds::Array{Array{Float64,1},1}, x::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}, a::Array{Float64, 1},
                                            b::Array{Float64, 1})
    n, m = length(x), length(y);
    figure(figsize = (5, 5))
    for k in 1:n
        ms = 10.0*a[k]/(mean(a)+mean(b));
        if k < n
            plot([x[k][1]], [x[k][2]], linestyle = "none", marker = ".", ms = ms, color = "red")
        else
            plot([x[k][1]], [x[k][2]], linestyle = "none", marker = ".", ms = ms, color = "red", label = L"$\mu_{GT}$")
        end
    end
    for k in 1:m
        ms = 10.0*b[k]/(mean(a)+mean(b));
        if k < m
            plot([y[k][1]], [y[k][2]], linestyle = "none", marker = ".", ms = ms, color = "green")
        else
            plot([y[k][1]], [y[k][2]], linestyle = "none", marker = ".", ms = ms, color = "green", label = L"$\mu_{est}$")
        end
    end

    ax = gca();
    ax.set_xlim([bounds[1][1], bounds[2][1]]);
    ax.set_ylim([bounds[1][2], bounds[2][2]]);
    legend()
    tight_layout()
    # savefig("figures/test.pdf", dpi = 100)
    show()
end

function plot_points_sources_uot(x::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}, a::Array{Float64,1},
                                            b::Array{Float64, 1}, bounds::Array{Array{Float64,1},1}, uot_out::UOT.UOT_output)
    n, m = length(x), length(y);
    figure(figsize = (5, 5))
    for k in 1:n
        ms = 25.0*a[k]/(mean(a)+mean(b));
        if k < n
            plot([x[k][1]], [x[k][2]], linestyle = "none", marker = ".", ms = ms, color = "red")
        else
            plot([x[k][1]], [x[k][2]], linestyle = "none", marker = ".", ms = ms, color = "red", label = L"$\mu_{GT}$")
        end
    end
    for k in 1:m
        ms = 25.0*b[k]/(mean(a)+mean(b));
        if k < m
            plot([y[k][1]], [y[k][2]], linestyle = "none", marker = ".", ms = ms, color = "green")
        else
            plot([y[k][1]], [y[k][2]], linestyle = "none", marker = ".", ms = ms, color = "green", label = L"$\mu_{est}$")
        end
    end
    #
    for i in 1:n
        for j in 1:m
            if uot_out.P[i, j] > 1e-8
                lx, ly = UOT_toolbox.line(x[i], y[j]);
                lw = 2*sum(uot_out.P.>0.0)*uot_out.P[i,j]*uot_out.C[i,j]/uot_out.cost;
                if lw > 0.0
                    plot(lx, ly, linestyle = "-", color = "black", lw = lw)
                else
                    println("### error ###")
                end
            end
        end
    end
    #
    for i in 1:n
        if abs(uot_out.mass_variation[i]) > 1e-8
            ms = 25.0*abs(uot_out.mass_variation[i])/(mean(a)+mean(b));
            if uot_out.mass_variation[i] > 0.0
                plot([x[i][1]], [x[i][2]], linestyle = "none", marker = "+", ms = ms, color = "black")
            else
                plot([x[i][1]], [x[i][2]], linestyle = "none", marker = "x", ms = ms, color = "black")
            end
        end
    end
    for j in 1:m
        if abs(uot_out.mass_variation[n + j]) > 1e-8
            ms = 25.0*abs(uot_out.mass_variation[n + j])/(mean(a)+mean(b));
            if uot_out.mass_variation[n + j] > 0.0
                plot([y[j][1]], [y[j][2]], linestyle = "none", marker = "+", ms = ms, color = "black")
            else
                plot([y[j][1]], [y[j][2]], linestyle = "none", marker = "x", ms = ms, color = "black")
            end
        end
    end

    ax = gca();
    ax.set_xlim([bounds[1][1], bounds[2][1]]);
    ax.set_ylim([bounds[1][2], bounds[2][2]]);
    title(string(L"Flat metric between $\mu_{est}$ to $\mu_{GT}$, cost=", round(uot_out.cost; digits = 5)))
    legend()
    tight_layout()
    # savefig("figures/test.pdf", dpi = 100)
    show()
end

end

max_stars_repo_name: qdenoyelle/OTMSL

max_issues_repo_name: qdenoyelle/OTMSL

max_forks_repo_name: qdenoyelle/OTMSL


text:
using SimplexSplitting, InvariantDistribution

@testset "InvDist" begin
    e_ex = Embedding(invariant_gaussian_embedding(npts = 25))
    t_ex = triang_from_embedding(e_ex)
    mm = mm_dd(t_ex)
    @test all(sum(mm, 2) .≈ 1)

    invdist = estimate_invdist(mm)
    @test sum(invdist.dist) ≈ 1
end


text:
module Optimizations

using ITensors

export gradient_descent, gd_error_tracker, loss_grad_wrap

include("peps/optimize.jl")
include("peps/accurate_peps.jl")
include("peps/approximate_peps.jl")
include("peps/peps_w_projectors.jl")
include("optimizers.jl")

end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function calc_error_norms(func, u, t, analyzer,
                          mesh::ParallelTreeMesh{2}, equations, initial_condition,
                          dg::DGSEM, cache, cache_analysis)
  l2_errors, linf_errors = calc_error_norms_per_element(func, u, t, analyzer,
                                                        mesh, equations, initial_condition,
                                                        dg, cache, cache_analysis)

  # Collect local error norms for each element on root process. That way, when aggregating the L2
  # errors, the order of summation is the same as in the serial case to ensure exact equality.
  # This facilitates easier parallel development and debugging (see
  # https://github.com/trixi-framework/Trixi.jl/pull/850#pullrequestreview-757463943 for details).
  # Note that this approach does not scale.
  if mpi_isroot()
    global_l2_errors = zeros(eltype(l2_errors), cache.mpi_cache.n_elements_global)
    global_linf_errors = similar(global_l2_errors)

    n_elements_by_rank = parent(cache.mpi_cache.n_elements_by_rank) # convert OffsetArray to Array
    l2_buf = MPI.VBuffer(global_l2_errors, n_elements_by_rank)
    linf_buf = MPI.VBuffer(global_linf_errors, n_elements_by_rank)
    MPI.Gatherv!(l2_errors, l2_buf, mpi_root(), mpi_comm())
    MPI.Gatherv!(linf_errors, linf_buf, mpi_root(), mpi_comm())
  else
    MPI.Gatherv!(l2_errors, nothing, mpi_root(), mpi_comm())
    MPI.Gatherv!(linf_errors, nothing, mpi_root(), mpi_comm())
  end

  # Aggregate element error norms on root process
  if mpi_isroot()
    # sum(global_l2_errors) does not produce the same result as in the serial case, thus a
    # hand-written loop is used
    l2_error = zero(eltype(global_l2_errors))
    for error in global_l2_errors
      l2_error += error
    end
    linf_error = reduce((x, y) -> max.(x, y), global_linf_errors)

    # For L2 error, divide by total volume
    total_volume_ = total_volume(mesh)
    l2_error = @. sqrt(l2_error / total_volume_)
  else
    l2_error = convert(eltype(l2_errors), NaN * zero(eltype(l2_errors)))
    linf_error = convert(eltype(linf_errors), NaN * zero(eltype(linf_errors)))
  end

  return l2_error, linf_error
end

function calc_error_norms_per_element(func, u, t, analyzer,
                                      mesh::ParallelTreeMesh{2}, equations, initial_condition,
                                      dg::DGSEM, cache, cache_analysis)
  @unpack vandermonde, weights = analyzer
  @unpack node_coordinates = cache.elements
  @unpack u_local, u_tmp1, x_local, x_tmp1 = cache_analysis

  # Set up data structures
  T = typeof(zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations)))
  l2_errors = zeros(T, nelements(dg, cache))
  linf_errors = copy(l2_errors)

  # Iterate over all elements for error calculations
  for element in eachelement(dg, cache)
    # Interpolate solution and node locations to analysis nodes
    multiply_dimensionwise!(u_local, vandermonde, view(u,                :, :, :, element), u_tmp1)
    multiply_dimensionwise!(x_local, vandermonde, view(node_coordinates, :, :, :, element), x_tmp1)

    # Calculate errors at each analysis node
    volume_jacobian_ = volume_jacobian(element, mesh, cache)

    for j in eachnode(analyzer), i in eachnode(analyzer)
      u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations)
      diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations)
      l2_errors[element] += diff.^2 * (weights[i] * weights[j] * volume_jacobian_)
      linf_errors[element] = @. max(linf_errors[element], abs(diff))
    end
  end

  return l2_errors, linf_errors
end

function calc_error_norms(func, u, t, analyzer,
                          mesh::ParallelP4estMesh{2}, equations,
                          initial_condition, dg::DGSEM, cache, cache_analysis)
  @unpack vandermonde, weights = analyzer
  @unpack node_coordinates, inverse_jacobian = cache.elements
  @unpack u_local, u_tmp1, x_local, x_tmp1, jacobian_local, jacobian_tmp1 = cache_analysis

  # Set up data structures
  l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations))
  linf_error = copy(l2_error)
  volume = zero(real(mesh))

  # Iterate over all elements for error calculations
  for element in eachelement(dg, cache)
    # Interpolate solution and node locations to analysis nodes
    multiply_dimensionwise!(u_local, vandermonde, view(u,                :, :, :, element), u_tmp1)
    multiply_dimensionwise!(x_local, vandermonde, view(node_coordinates, :, :, :, element), x_tmp1)
    multiply_scalar_dimensionwise!(jacobian_local, vandermonde, inv.(view(inverse_jacobian, :, :, element)), jacobian_tmp1)

    # Calculate errors at each analysis node
    @. jacobian_local = abs(jacobian_local)

    for j in eachnode(analyzer), i in eachnode(analyzer)
      u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations)
      diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations)
      l2_error += diff.^2 * (weights[i] * weights[j] * jacobian_local[i, j])
      linf_error = @. max(linf_error, abs(diff))
      volume += weights[i] * weights[j] * jacobian_local[i, j]
    end
  end

  # Accumulate local results on root process
  global_l2_error = Vector(l2_error)
  global_linf_error = Vector(linf_error)
  MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm())
  MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm())
  total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
  if mpi_isroot()
    l2_error   = convert(typeof(l2_error),   global_l2_error)
    linf_error = convert(typeof(linf_error), global_linf_error)
    # For L2 error, divide by total volume
    l2_error = @. sqrt(l2_error / total_volume)
  else
    l2_error   = convert(typeof(l2_error),   NaN * global_l2_error)
    linf_error = convert(typeof(linf_error), NaN * global_linf_error)
  end

  return l2_error, linf_error
end

function integrate_via_indices(func::Func, u,
                               mesh::ParallelTreeMesh{2}, equations, dg::DGSEM, cache,
                               args...; normalize=true) where {Func}
  # call the method accepting a general `mesh::TreeMesh{2}`
  # TODO: MPI, we should improve this; maybe we should dispatch on `u`
  #       and create some MPI array type, overloading broadcasting and mapreduce etc.
  #       Then, this specific array type should also work well with DiffEq etc.
  local_integral = invoke(integrate_via_indices,
    Tuple{typeof(func), typeof(u), TreeMesh{2}, typeof(equations),
          typeof(dg), typeof(cache), map(typeof, args)...},
    func, u, mesh, equations, dg, cache, args..., normalize=normalize)

  # OBS! Global results are only calculated on MPI root, all other domains receive `nothing`
  global_integral = MPI.Reduce!(Ref(local_integral), +, mpi_root(), mpi_comm())
  if mpi_isroot()
    integral = convert(typeof(local_integral), global_integral[])
  else
    integral = convert(typeof(local_integral), NaN * local_integral)
  end

  return integral
end

function integrate_via_indices(func::Func, u,
                               mesh::ParallelP4estMesh{2}, equations,
                               dg::DGSEM, cache, args...; normalize=true) where {Func}
  @unpack weights = dg.basis

  # Initialize integral with zeros of the right shape
  integral = zero(func(u, 1, 1, 1, equations, dg, args...))
  volume = zero(real(mesh))

  # Use quadrature to numerically integrate over entire domain
  for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
      integral += volume_jacobian * weights[i] * weights[j] * func(u, i, j, element, equations, dg, args...)
      volume += volume_jacobian * weights[i] * weights[j]
    end
  end

  global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm())
  total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
  if mpi_isroot()
    integral = convert(typeof(integral), global_integral[])
  else
    integral = convert(typeof(integral), NaN * integral)
    total_volume = volume # non-root processes receive nothing from reduce -> overwrite
  end

  # Normalize with total volume
  if normalize
    integral = integral / total_volume
  end

  return integral
end

end # @muladd


text:
#===========================java.util.Properties===============================#
#"The Properties class represents a persistent set of properties...Each key and
#its corresponding value in the property list is a string."
#https://docs.oracle.com/javase/7/docs/api/java/util/Properties.html
#
#WRAPPER
type Properties
    jproperties::JProperties
end

#CONSTRUCTORS
##key, value argument
function Properties(key::AbstractString, value::AbstractString)
    jproperties = JProperties(())
    jcall(jproperties, "setProperty", JObject, (JString, JString), key, value)
    return Properties(jproperties)
end

##Dictionary argument
function Properties(d::Dict{String, String})
    jproperties = JProperties(())
    for (k, v) in d
        jcall(jproperties, "setProperty", JObject, (JString, JString), k, v)
    end
    return Properties(jproperties)
end

##Value-only argument. Key defaults to "annotators"
function Properties(value::AbstractString)
    return Properties("annotators", value)
end

#CONVENIENCE METHODS
##Returns a string 'val' corresponding to property 'key'
function get_property(jprops::JProperties, key::AbstractString="annotators")
    return jcall(jprops, "getProperty", JString, (JString,), key)
end
#------------------------------------------------------------------------------#

#===============edu.stanford.nlp.pipeline.StanfordCoreNLP======================#
#"This is a pipeline that takes in a string and returns various analyzed linguistic
#forms. The String is tokenized via a tokenizer (using a TokenizerAnnotator),
#and then other sequence model style annotation can be used to add things like
#lemmas, POS tags, and named entities. These are returned as a list of CoreLabels.
#Other analysis components build and store parse trees, dependency graphs, etc."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/pipeline/StanfordCoreNLP.html
#
#WRAPPER
##StanfordCoreNLP wrapper
type StanfordCoreNLP
    jstanfordcorenlp::JStanfordCoreNLP
end

#CONSTRUCTORS
##Usage:
##      StanfordCoreNLP(Dict("annotations" => "tokenize,ssplit,pos")
##
##Annotations must be to the right of their dependencies. See--
##https://stanfordnlp.github.io/CoreNLP/dependencies.html
function StanfordCoreNLP(dict::Dict{String, String})
    properties = Properties(dict)
    return StanfordCoreNLP(JStanfordCoreNLP((JProperties,), properties))
end

##Usage:
##      StanfordCoreNLP("pos", "parse")
##
##Annotation arguments may be unordered. Dependencies are added automatically.
function StanfordCoreNLP(annotators::String...)
    if check_annotators(annotators...)
        annList = order!(dependencies(annotators...))
        annCSVString = join(annList, ",")
        jprops = JProperties(())
        jcall(jprops, "setProperty", JObject, (JString, JString), "annotators", annCSVString)
        return StanfordCoreNLP(JStanfordCoreNLP((JProperties,), jprops))
    else
        println("Invalid annotator in `StanfordCoreNLP` call.")
        return 1
    end
end

#CONVENIENCE METHODS
##Returns true if all arguments are valid annotators, false otherwise.
function check_annotators(annotators...)
    valid = true
        for ann in annotators
            valid = valid && in(ann, keys(ANNOTATOR_DEPENDENCIES))
        end
    return valid
end

##Accepts an unordered list of annotators and adds all dependencies. For example,
##annotator_list("ssplit") returns ["ssplit", "tokenize"].
function dependencies(annotators...)
    list = String[]
    for a in annotators
        for b in ANNOTATOR_DEPENDENCIES[a]
            push!(list, b)
        end
        push!(list, a)
    end
    return unique(list)
end

##Accepts a list of annotators and orders it in place according to ANNOTATORS.
function order!(annotators::Array{String,1})
    list = String[]
    for ann in ORDERED_ANNOTATOR_LIST
        if in(ann, annotators)
            push!(list, ann)
        end
    end
    annotators = list
end

#Returns a JProperties of a StanfordCoreNLP pipeline corresponding to `key`.
#Defaults to key = "annotators".
function get_properties(pipeline::StanfordCoreNLP)
    return jcall(pipeline.jstanfordcorenlp, "getProperties", JProperties, ())
end

#Return a String that gives detailed human-readable information about how much
#time was spent by each annotator and by the entire annotation pipeline.
function timinginformation(pipeline::StanfordCoreNLP)
    return jcall(pipeline.jstanfordcorenlp, "timingInformation", JString, ())
end
#------------------------------------------------------------------------------#

#====================edu.stanford.nlp.pipeline.Annotation======================#
#"An annotation representing a span of text in a document."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/pipeline/Annotation.html
#
#WRAPPER
type Annotation
    jannotation::JAnnotation
end

#CONSTRUCTORS
##function Annotation(text::AbstractString)
##    jannotation = JAnnotation((JString,), text)
##    return Annotation(jannotation)
##end
function Annotation(text::AbstractString)
    return Annotation(JAnnotation((JString,), text))
end

#CONVENIENCE METHODS
##Returns iterable list of sentences in an annotated text.
function sentences(doc::Annotation)
    return narrow(jcall(doc.jannotation, "get", JObject, (JClass,), JSentencesAnnotationClass))
end

##Returns a CorefChainAnnotation HashMap. Call 'values' to get an iterator.
function corefchainannotation(doc::Annotation)
    return narrow(jcall(doc.jannotation, "get", JObject, (JClass,), JCorefChainAnnotationClass))
end

##Returns iterable list of tokens in an annotated sentence.
function tokens(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JTokensAnnotationClass))
end

##Returns a collection of relation triples from an annotated sentence.
function relationtriples(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JRelationTriplesAnnotationClass))
end

##Returns a LabeledScoredTreeNode object, the syntactic parse tree of a sentence.
function labeledscoredtreenode(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JTreeAnnotationClass))
end

##Returns a SemanticGraph.
function basicdependencies(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JBasicDependenciesAnnotationClass))
end

##Returns a SemanticGraph.
function collapseddependencies(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JCollapsedDependenciesAnnotationClass))
end

##Returns a SemanticGraph.
function collapsedccprocesseddependencies(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JCollapsedCCProcessedDependenciesAnnotationClass))
end

##Returns an iterable set of edu.stanford.nlp.coref.data.Mention.
function corefmentionsannotation(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JCorefMentionsAnnotationClass))
end

##Returns Sentiment
function sentimentannotation(sentence::JAnnotation)
    return narrow(jcall(sentence, "get", JObject, (JClass,), JSentimentAnnotatedTreeClass))
end

##Performs the annotation in place.
function annotate!(pipeline::StanfordCoreNLP, doc::Annotation)
    jcall(pipeline.jstanfordcorenlp, "annotate", Void, (JAnnotation,), doc.jannotation)
end
#------------------------------------------------------------------------------#

#=============================java.util.HashMap================================#
#https://docs.oracle.com/javase/7/docs/api/java/util/HashMap.html
#
function values(hashmap::JHashMap)
    return jcall(hashmap, "values", JCollection, ())
end
#------------------------------------------------------------------------------#

#==========================edu.stanford.nlp.coref.data.Mention=================#
function to_string(mention::JMention)
    return jcall(mention, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#======edu.stanford.nlp.dcoref.CorefCoreAnnotations.CorefChainAnnotation=======#
function to_string(corefchain::JCorefChain)
    return jcall(corefchain, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#================edu.stanford.nlp.trees.LabeledScoredTreeNode==================#
function children(jlstn::JLabeledScoredTreeNode)
    return jcall(jlstn, "children", JArray, ())  #Unsure of return type: JArray? JList? JCollection?
end

##Returns a Set of Constituents generated by the parse tree.
function constituents(jlstn::JLabeledScoredTreeNode)
    return jcall(jlstn, "constituents", JSet, ())
end
#------------------------------------------------------------------------------#

#====================edu.stanford.nlp.trees.Constituent========================#
#"A Constituent object defines a generic edge in a graph."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/trees/Constituent.html
#
function to_string(constituent::JConstituent)
    return jcall(constituent, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#===================edu.stanford.nlp.ie.util.RelationTriple====================#
#"A (subject, relation, object) triple."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/ie/util/RelationTriple.html
#
##The subject of this relation triple, as a String
function subjectGloss(triple::JRelationTriple)
    return jcall(triple, "subjectGloss", JString, ())
end

##The object of this relation triple, as a String
function objectGloss(triple::JRelationTriple)
    return jcall(triple, "objectGloss", JString, ())
end

##The relation of this relation triple, as a String
function relationGloss(triple::JRelationTriple)
    return jcall(triple, "relationGloss", JString, ())
end

##The subject lemma of this relation triple, as a String
function subjectLemmaGloss(triple::JRelationTriple)
    return jcall(triple, "subjectLemmaGloss", JString, ())
end

##The object lemma of this relation triple, as a String
function objectLemmaGloss(triple::JRelationTriple)
    return jcall(triple, "objectLemmaGloss", JString, ())
end

##The relation lemma of this relation triple, as a String
function relationLemmaGloss(triple::JRelationTriple)
    return jcall(triple, "relationLemmaGloss", JString, ())
end

##Returns a JList of JCoreLabels that represent the given relation triple as a flat sentence.
function asSentence(triple::JRelationTriple)
    return jcall(triple, "asSentence", JList, ())
end

##If true, this relation expresses a "to be" relation.
function isPrefixBe(triple::JRelationTriple)
    return jcall(triple, "isPrefixBe", jboolean, ())
end

##If true, this relation expresses a "to be" relation (with the be at the end of the sentence).
function isSuffixBe(triple::JRelationTriple)
    return jcall(triple, "	isSuffixBe", jboolean, ())
end

##If true, this relation has an ungrounded "of" at the end of the relation.
function isSuffixOf(triple::JRelationTriple)
    return jcall(triple, "	isSuffixOf", jboolean, ())
end

##Returns human-readable description of this relation triple, as a tab-separated line.
function to_string(triple::JRelationTriple)
    return jcall(triple, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#====================edu.stanford.nlp.ling.CoreLabel===========================#
##"A CoreLabel represents a single word with ancillary information attached using
##CoreAnnotations. A CoreLabel also provides convenient methods to access tags,
##lemmas, etc. (if the proper annotations are set). A CoreLabel is a Map from keys
##(which are Class objects) to values, whose type is determined by the key."
##https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/ling/CoreLabel.html
##
#WRAPPER
type CoreLabel
    token::JCoreLabel
    index::Int
    lemma::String
    pos::String
    ner::String
    originaltext::String
    word::String
    value::String
end

#CONSTRUCTORS
function CoreLabel(token::JCoreLabel)
    index = jcall(token, "index", jint, ())
    lemma = jcall(token, "lemma", JString, ())
    pos = jcall(token, "tag", JString, ())
    ner = jcall(token, "ner", JString, ())
    originaltext = jcall(token, "originalText", JString, ())
    word = jcall(token, "word", JString, ())
    value = jcall(token, "value", JString, ())
    return CoreLabel(token, index, lemma, pos, ner, originaltext, word, value)
end

#CONVENIENCE METHODS
##Converts a JCoreLabel object to String.
function to_string(token::JCoreLabel)
    return jcall(token, "toString", JString, ())
end

##Converts the .token field of a CoreLabel object to String.
function to_string(corelabel::CoreLabel)
    return jcall(corelabel.jcorelabel, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#=======================edu.stanford.nlp.trees.Tree============================#
#WRAPPER
type Tree
    jtree::JTree
    size::jint
    depth::jint
    constituents::JSet
    iterator::JIterator
end

#CONSTRUCTORS
function Tree(jtree::JTree)
    size = jcall(tree, "size", jint, ())
    depth = jcall(tree, "depth", jint, ())
    constituents = jcall(tree, "constituents", JString, ())
    iterator = jcall(tree, "iterator", JString, ())
    return(jtree, size, depth, constituents, iterator)
end

#CONVENIENCE METHODS
##Converts parse tree to string in Penn Treebank format.
function to_string(jtree::JTree)
    return jcall(jtree, "toString", JString, ())
end

##Converts the .jtree field of the Tree wrapper to string in Penn Treebank format.
function to_string(tree::Tree)
    return jcall(tree.jtree, "toString", JString, ())
end

function to_string(jlstn::JLabeledScoredTreeNode)
    return jcall(jlstn, "toString", JString, ())
end

##Converts parse tree to string in Penn Treebank format.
function pennstring(jtree::JTree)
    return jcall(jtree, "pennString", JString, ())
end

##Converts parse tree to string in Penn Treebank format.
function pennstring(tree::Tree)
    return jcall(tree.jtree, "pennString", JString, ())
end

#Returns a list of children (trees) for the current node.
function getChildrenAsList(jlstn::JLabeledScoredTreeNode) #JTree?
    return jcall(jlstn, "getChildrenAsList", JList, ())
end
#------------------------------------------------------------------------------#

#=================edu.stanford.nlp.tagger.maxent.MaxentTagger==================#
##The main class for users to run, train, and test the part of speech tagger.
##https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/tagger/maxent/MaxentTagger.html
##
##WRAPPER
#=
type MaxentTagger
    jmet::JMaxentTagger
end

##CONSTRUCTORS
##Path is the location of parameter files for a trained tagger.
function MaxentTagger(path::AbstractString)
    jmet = JMaxentTagger((JString,), path)
    return MaxentTagger(jmet)
end

##List<TaggedWord> tagged = tagger.tagSentence(sentence)
function tagsentence(tagger::MaxentTagger, sentence::JList)
    jcall(tagger.jmet, "tagSentence", JList, (JList,), sentence)
end
=#
#------------------------------------------------------------------------------#

#edu.stanford.nlp.semgraph.SemanticGraph=======================================#
#"Represents a semantic graph of a sentence or document, with IndexedWord objects for nodes."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/semgraph/SemanticGraph.html
#SemanticGraph's used mainly to represent dependencies.

#Convert SemanticGraph to String
function to_string(semgraph::JSemanticGraph)
    jcall(semgraph, "toString", JString, ())
end

#Convert SemanticGraph to list.
function to_list(semgraph::JSemanticGraph)
    return jcall(semgraph, "toList", JString, ())
end

function size(semgraph::JSemanticGraph)
    return jcall(semgraph, "size", jint, ())
end

function typedDependencies(semgraph::JSemanticGraph)
    return jcall(semgraph, "typedDependencies", JCollection, ())
end
#------------------------------------------------------------------------------#

#===============edu.stanford.nlp.parser.nndep.DependencyParser=================#
#Unused in `function test' as of Oct 2, 2017
#
#DependencyParser parser = DependencyParser.loadFromModelFile(modelPath);

function dependency_parser(modelPath::AbstractString)
    return narrow(jcall(JDependencyParser, "loadFromModelFile",
                                JDependencyParser, (JString,), modelPath))
end

#GrammaticalStructure gs = parser.predict(tagged);
function predict(parser::JDependencyParser , tagged::JList)
    jcall(parser, "predict", JGrammaticalStructure, (JList,), tagged)
end

function to_string(gs::JGrammaticalStructure)
    return jcall(gs, "toString", JString, ())
end

#------------------------------------------------------------------------------#

#===================edu.stanford.nlp.trees.TypedDependency=====================#
#"A TypedDependency is a relation between two words in a GrammaticalStructure.
#Each TypedDependency consists of a governor word, a dependent word, and a
#relation, which is normally an instance of GrammaticalRelation."
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/trees/TypedDependency.html
#
function to_string(typeddependency::JTypedDependency)
    return jcall(typeddependency, "toString", JString, ())
end
#------------------------------------------------------------------------------#

#============================java.io.StringReader==============================#
#Unused in `function test' as of Oct 2, 2017
#=
type StringReader
    jsr::JStringReader
end

function StringReader(text::AbstractString)
    jsr = JStringReader((JString,), text)
    return StringReader(jsr)
end
=#
#------------------------------------------------------------------------------#

#=================edu.stanford.nlp.process.DocumentPreprocessor================#
#Unused in `function test' as of Oct 2, 2017
#
#=
type DocumentPreprocessor
    jdp::JDocumentPreprocessor
end

function DocumentPreprocessor(stringreader::StringReader)
    jdp = JDocumentPreprocessor((JReader,), stringreader.jsr)
    return DocumentPreprocessor(jdp)
end
=#
#------------------------------------------------------------------------------#

#==============edu.stanford.nlp.parser.lexparser.LexicalizedParser=============#
#https://nlp.stanford.edu/nlp/javadoc/javanlp/edu/stanford/nlp/parser/lexparser/LexicalizedParser.html
#
#WRAPPER
type LexicalizedParser
    jlp::JLexicalizedParser
end

#CONSTRUCTORS
function LexicalizedParser(parserFile::String)
    jlp = JLexicalizedParser((JString,), parserFile)
    return LexicalizedParser(jlp)
end

#CONVENIENCE METHODS
function load_lexalized_parser_model()
    return jcall(JLexicalizedParser, "loadModel", JLexicalizedParser, ())
end

#Parse sentence::JList to JTree using lp::JLexicalizedParser
function lexparse(lp::JLexicalizedParser, sentence::JAnnotation)
    return jcall(lp, "apply", JTree, (JAnnotation,), sentence)
end

#Convert Tree to String
function to_string(parse_tree::JTree)
    jcall(parse_tree, "toString", JString, ())
end
#------------------------------------------------------------------------------#


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

@doc raw"""ScaleStatus represents the current status of a scale subresource.

    IoK8sApiAutoscalingV1ScaleStatus(;
        replicas=nothing,
        selector=nothing,
    )

    - replicas::Int32 : actual number of observed instances of the scaled object.
    - selector::String : label query over pods that should match the replicas count. This is same as the label selector but in the string format to avoid introspection by clients. The string will be in the same format as the query-param syntax. More info about label selectors: http://kubernetes.io/docs/user-guide/labels#label-selectors
"""
mutable struct IoK8sApiAutoscalingV1ScaleStatus <: SwaggerModel
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    selector::Any # spec type: Union{ Nothing, String } # spec name: selector

    function IoK8sApiAutoscalingV1ScaleStatus(;replicas=nothing, selector=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV1ScaleStatus, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAutoscalingV1ScaleStatus, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        o
    end
end # type IoK8sApiAutoscalingV1ScaleStatus

const _property_map_IoK8sApiAutoscalingV1ScaleStatus = Dict{Symbol,Symbol}(Symbol("replicas")=>Symbol("replicas"), Symbol("selector")=>Symbol("selector"))
const _property_types_IoK8sApiAutoscalingV1ScaleStatus = Dict{Symbol,String}(Symbol("replicas")=>"Int32", Symbol("selector")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV1ScaleStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV1ScaleStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV1ScaleStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV1ScaleStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV1ScaleStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV1ScaleStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV1ScaleStatus)
    (getproperty(o, Symbol("replicas")) === nothing) && (return false)
    true
end

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


text:
int_rules_4_1_8 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.8*trig^m*(a+b*cos^p+c*sin^q)^n =#
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * cos(~(d') + ~(e') * ~x) ^ ~p + ~(c') * sin(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-f / ~e) * Subst(Antiderivative(ExpandToSum(~c + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2 - ~p / 2) + ~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~q) / 2 + 1), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && (GtQ(~p, 0) && LeQ(~p, ~q))))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * sin(~(d') + ~(e') * ~x) ^ ~p + ~(c') * cos(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f / ~e) * Subst(Antiderivative(ExpandToSum(~c + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2 - ~p / 2) + ~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~q) / 2 + 1), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && (GtQ(~p, 0) && LeQ(~p, ~q))))))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * cos(~(d') + ~(e') * ~x) ^ ~p + ~(c') * sin(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-f / ~e) * Subst(Antiderivative(ExpandToSum(~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2) + ~b * f ^ ~p * (~x) ^ ~p + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2 - ~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && LtQ(0, ~q, ~p)))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * sin(~(d') + ~(e') * ~x) ^ ~p + ~(c') * cos(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f / ~e) * Subst(Antiderivative(ExpandToSum(~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2) + ~b * f ^ ~p * (~x) ^ ~p + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2 - ~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && LtQ(0, ~q, ~p)))))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * cos(~(d') + ~(e') * ~x) ^ ~p + ~(c') * sin(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-f / ~e) * Subst(Antiderivative(ExpandToSum(~c + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2 - ~p / 2) + ~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~q) / 2 + 1), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && (GtQ(~p, 0) && LeQ(~p, ~q))))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * sin(~(d') + ~(e') * ~x) ^ ~p + ~(c') * cos(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f / ~e) * Subst(Antiderivative(ExpandToSum(~c + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2 - ~p / 2) + ~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~q) / 2 + 1), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && (GtQ(~p, 0) && LeQ(~p, ~q))))))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * cos(~(d') + ~(e') * ~x) ^ ~p + ~(c') * sin(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-f / ~e) * Subst(Antiderivative(ExpandToSum(~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2) + ~b * f ^ ~p * (~x) ^ ~p + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2 - ~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && LtQ(0, ~q, ~p)))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~a + ~(b') * sin(~(d') + ~(e') * ~x) ^ ~p + ~(c') * cos(~(d') + ~(e') * ~x) ^ ~q) ^ ~n, ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f / ~e) * Subst(Antiderivative(ExpandToSum(~a * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2) + ~b * f ^ ~p * (~x) ^ ~p + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ (~p / 2 - ~q / 2), ~x) ^ ~n / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~p / 2) && (IntegerQ(~q / 2) && (IntegerQ(~n) && LtQ(0, ~q, ~p)))))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.8 trig^m (a+b cos^p+c sin^q)^n.jl


text:
function cross_map(r::MersenneTwister)
    map = YAML.load_file("maps/cross.yaml")
    exits = [(4, 1), (1, 4), (4, 7), (7, 4)]
    start = rand(r, exits)
    exit = rand(r, setdiff(exits, [start]))
    map["starts"] = [Dict("x"=>start[2], "y"=>start[1])]
    map["exits"] = [Dict("x"=>exit[2], "y"=>exit[1])]
    push!(map["objects"], Dict("type"=>"food", "color"=>"green",
                               "x"=>exit[2], "y"=>exit[1]))
    map
end

function cross_map(seed::Int64)
    r = MersenneTwister(seed)
    cross_map(r)
end

function cross_meta()
    Dict{String, Any}("max_step" => 100)
end

function cross_search_fitness(cont_f::Function; seed::Int64=0)
    e = Episode(Grid(cross_map(seed)); reward=food_reward, meta=cross_meta())
    run!(e, cont_f)
    steps = terminate(e) ? e.meta["max_step"] : e.meta["step"]
    1 - steps / e.meta["max_step"]
end

function cross_memorize_fitness(cont_f::Function; seed::Int64=0)
    map_function() = cross_map(seed)
    learn_fitness(cont_f, map_function, cross_meta)
end

function cross_map(r::MersenneTwister, strategy::Int64)
    map = YAML.load_file("maps/cross.yaml")
    exits = [(4, 1), (1, 4), (4, 7), (7, 4)]
    si = rand(r, 1:4)
    start = exits[si]
    exit = exits[mod((si-1)+strategy, 4)+1]
    map["starts"] = [Dict("x"=>start[2], "y"=>start[1])]
    map["exits"] = [Dict("x"=>exit[2], "y"=>exit[1])]
    push!(map["objects"], Dict("type"=>"food", "color"=>"green",
                               "x"=>exit[2], "y"=>exit[1]))
    map
end

function cross_strategy_fitness(cont_f::Function; seed::Int64=0)
    r = MersenneTwister(seed)
    strategy = rand(r, 1:3)
    map_function() = cross_map(r, strategy)
    learn_fitness(cont_f, map_function, cross_meta)
end


text:
# ---
# title: 44. Wildcard Matching
# id: problem44
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming, Backtracking, Greedy
# link: <https://leetcode.com/problems/wildcard-matching/description/>
# hidden: true
# ---
#
# Given an input string (`s`) and a pattern (`p`), implement wildcard pattern
# matching with support for `'?'` and `'*'` where:
#
#   * `'?'` Matches any single character.
#   * `'*'` Matches any sequence of characters (including the empty sequence).
#
# The matching should cover the **entire** input string (not partial).
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aa", p = "a"
#     Output: false
#     Explanation: "a" does not match the entire string "aa".
#
#
# **Example 2:**
#
#
#
#     Input: s = "aa", p = "*"
#     Output: true
#     Explanation:  '*' matches any sequence.
#
#
# **Example 3:**
#
#
#
#     Input: s = "cb", p = "?a"
#     Output: false
#     Explanation:  '?' matches 'c', but the second letter is 'a', which does not match 'b'.
#
#
# **Example 4:**
#
#
#
#     Input: s = "adceb", p = "*a*b"
#     Output: true
#     Explanation:  The first '*' matches the empty sequence, while the second '*' matches the substring "dce".
#
#
# **Example 5:**
#
#
#
#     Input: s = "acdcb", p = "a*c?b"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length, p.length <= 2000`
#   * `s` contains only lowercase English letters.
#   * `p` contains only lowercase English letters, `'?'` or `'*'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "rendering.jl" begin
    img = WordCloud.rendertext("test", 88.3, color="blue", angle=20, border=2)
    mat, svg = WordCloud.rendertext(Random.randstring(rand(1:8)), rand(5:50), angle=rand(0:180), type=:both)
    @test all(WordCloud.svg2bitmap(svg) .≈ mat)
end

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
# Autogenerated wrapper script for libmicrohttpd_jll for x86_64-linux-gnu
export libmicrohttpd

using GnuTLS_jll
JLLWrappers.@generate_wrapper_header("libmicrohttpd")
JLLWrappers.@declare_library_product(libmicrohttpd, "libmicrohttpd.so.12")
function __init__()
    JLLWrappers.@generate_init_header(GnuTLS_jll)
    JLLWrappers.@init_library_product(
        libmicrohttpd,
        "lib/libmicrohttpd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()