text:
function context_test_callback(arg1, arg2, arg3)
    # We're not really testing it because, nvidia doesn't seem to care about this functionality:
    # https://devtalk.nvidia.com/default/topic/497433/context-callback-never-called/
    OpenCL.cl.log_error("Callback works")
    return
end

function create_context_error(ctx)
    empty_kernel = "
    __kernel void test() {
        int c = 1 + 1;
    };"
    try
        p = cl.Program(ctx, source = empty_kernel) |> cl.build!
        k = cl.Kernel(p, "test")
        q = cl.CmdQueue(ctx)
        q(k, 1, 10000000)
    catch
    end
end

@testset "OpenCL.Context" begin
    @testset "OpenCL.Context constructor" begin
        @test_throws MethodError (cl.Context([]))
        for platform in cl.platforms()
            for device in cl.devices(platform)
                ctx = cl.Context(device)
                @test ctx != nothing
                ctx_id = ctx.id
                ctx2 = cl.Context(ctx_id)
                @test cl.is_ctx_id_alive(ctx_id)
                @test ctx.id != C_NULL
                @test ctx2.id != C_NULL
                finalize(ctx)
                @test ctx.id == C_NULL
                @test ctx2.id != C_NULL
                @test cl.is_ctx_id_alive(ctx_id)
                finalize(ctx2)
                @test ctx.id == C_NULL
                @test ctx2.id == C_NULL
                # jeez, this segfaults... WHY? I suspect a driver bug for refcount == 0?
                # NVIDIA 381.22
                #@test !cl.is_ctx_id_alive(ctx_id)
                @testset "Context callback" begin
                    ctx = cl.Context(device, callback = context_test_callback)
                    create_context_error(ctx)
                end
            end
        end
    end

    @testset "OpenCL.Context platform properties" begin
        for platform in cl.platforms()
            try
                cl.Context(cl.CL_DEVICE_TYPE_CPU)
            catch err
                @test typeof(err) == cl.CLError
                # CL_DEVICE_NOT_FOUND could be throw for GPU only drivers
                @test err.desc in (:CL_INVALID_PLATFORM,
                                         :CL_DEVICE_NOT_FOUND)
            end

            if platform[:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Context platform properties for " *
                     "Portable Computing Language Platform")
                continue
            end

            properties = [(cl.CL_CONTEXT_PLATFORM, platform)]
            for (cl_dev_type, sym_dev_type) in [(cl.CL_DEVICE_TYPE_CPU, :cpu),
                                                (cl.CL_DEVICE_TYPE_GPU, :gpu)]
                if !cl.has_device_type(platform, sym_dev_type)
                    continue
                end
                @test cl.Context(sym_dev_type, properties=properties) != nothing
                @test cl.Context(cl_dev_type, properties=properties) != nothing
                ctx = cl.Context(cl_dev_type, properties=properties)
                @test isempty(cl.properties(ctx)) == false
                test_properties = cl.properties(ctx)

                @test test_properties == properties

                platform_in_properties = false
                for (t, v) in test_properties
                    if t == cl.CL_CONTEXT_PLATFORM
                        @test v[:name] == platform[:name]
                        @test v == platform
                        platform_in_properties = true
                        break
                    end
                end
                @test platform_in_properties
            end
            try
                ctx2 = cl.Context(cl.CL_DEVICE_TYPE_ACCELERATOR,
                                  properties=properties)
            catch err
                @test typeof(err) == cl.CLError
                @test err.desc == :CL_DEVICE_NOT_FOUND
            end
        end
    end

    @testset "OpenCL.Context create_some_context" begin
        @test cl.create_some_context() != nothing
        @test typeof(cl.create_some_context()) == cl.Context
    end

   @testset "OpenCL.Context parsing" begin
        for platform in cl.platforms()
            properties = [(cl.CL_CONTEXT_PLATFORM, platform)]
            parsed_properties = cl._parse_properties(properties)

            @test isodd(length(parsed_properties))
            @test parsed_properties[end] == 0
            @test parsed_properties[1] == cl.cl_context_properties(cl.CL_CONTEXT_PLATFORM)
            @test parsed_properties[2] == cl.cl_context_properties(platform.id)
        end
    end

end


text:
import Distributions
using SpecialFunctions: lgamma, lbeta, digamma

#########################
# abstract distribution #
#########################

abstract type Distribution{T} end

function random end
function logpdf end
get_return_type(::Distribution{T}) where {T} = T

export Distribution
export random
export logpdf
export logpdf_grad

#########
# dirac #
#########

struct Dirac <: Distribution{Float64} end

"""
    dirac(val::Real)

Deterministically return a `Float64` that is equal to the given value.
`logpdf(dirac, x, y)` returns `0` if `x == y` and `-Inf` otherwise.
"""
const dirac = Dirac()

logpdf(::Dirac, x::Real, y::Real) = (x == y ? 0. : -Inf)

random(::Dirac, y::Real) = Float64(y)

(::Dirac)(y) = random(Dirac(), y)

has_output_grad(::Dirac) = false
has_argument_grads(::Dirac) = (false,)
get_static_argument_types(::Dirac) = [Float64]

export dirac

#############
# bernoulli #
#############

struct Bernoulli <: Distribution{Bool} end

"""
    bernoulli(prob_true::Real)

Samples a `Bool` value which is true with given probability
"""
const bernoulli = Bernoulli()

function logpdf(::Bernoulli, x::Bool, prob::Real)
    x ? log(prob) : log(1. - prob)
end

function logpdf_grad(::Bernoulli, x::Bool, prob::Real)
    prob_grad = x ? 1. / prob : -1. / (1-prob)
    (nothing, prob_grad)
end

random(::Bernoulli, prob::Real) = rand() < prob

(::Bernoulli)(prob) = random(Bernoulli(), prob)

has_output_grad(::Bernoulli) = false
has_argument_grads(::Bernoulli) = (true,)
get_static_argument_types(::Bernoulli) = [Float64]

export bernoulli

##########
# normal #
##########

struct Normal <: Distribution{Float64} end

"""
    normal(mu::Real, std::Real)

Samples a `Float64` value from a normal distribution.
"""
const normal = Normal()

function logpdf(::Normal, x::Real, mu::Real, std::Real)
    var = std * std
    diff = x - mu
    -(diff * diff)/ (2.0 * var) - 0.5 * log(2.0 * pi * var)
end

function logpdf_grad(::Normal, x::Real, mu::Real, std::Real)
    precision = 1. / (std * std)
    diff = mu - x
    deriv_x = diff * precision
    deriv_mu = -deriv_x
    deriv_sigma = -0.5 * precision * (1. - precision * (diff * diff))
    (deriv_x, deriv_mu, deriv_sigma)
end

random(::Normal, mu::Real, std::Real) = mu + std * randn()

(::Normal)(mu, std) = random(Normal(), mu, std)

has_output_grad(::Normal) = true
has_argument_grads(::Normal) = (true, true)
get_static_argument_types(::Normal) = [Float64, Float64]

export normal

#########
# gamma #
#########

struct Gamma <: Distribution{Float64} end

"""
    gamma(shape::Real, scale::Real)

Sample a `Float64` from a gamma distribution.
"""
const gamma = Gamma()

function logpdf(::Gamma, x::Real, shape::Real, scale::Real)
    if x > 0.
        (shape - 1.0) * log(x) - (x / scale) - shape * log(scale) - lgamma(shape)
    else
        -Inf
    end
end

function logpdf_grad(::Gamma, x::Real, shape::Real, scale::Real)
    error("Not Implemented")
    (nothing, nothing, nothing)
end

function random(::Gamma, shape::Real, scale::Real)
    rand(Distributions.Gamma(shape, scale))
end

(::Gamma)(shape, scale) = random(Gamma(), shape, scale)

has_output_grad(::Gamma) = false
has_argument_grads(::Gamma) = (false, false)
get_static_argument_types(::Gamma) = [Float64, Float64]

export gamma

#################
# inverse gamma #
#################

struct InverseGamma <: Distribution{Float64} end

"""
    inv_gamma(shape::Real, scale::Real)

Sample a `Float64` from a inverse gamma distribution.
"""
const inv_gamma = InverseGamma()

function logpdf(::InverseGamma, x::Real, shape::Real, scale::Real)
    if x > 0.
        shape * log(scale) - (shape + 1) * log(x) - lgamma(shape) - (scale / x)
    else
        -Inf
    end
end

function logpdf_grad(::InverseGamma, x::Real, shape::Real, scale::Real)
    error("Not Implemented")
    (nothing, nothing, nothing)
end

function random(::InverseGamma, shape::Real, scale::Real)
    rand(Distributions.InverseGamma(shape, scale))
end

(::InverseGamma)(shape, scale) = random(InverseGamma(), shape, scale)

has_output_grad(::InverseGamma) = false
has_argument_grads(::InverseGamma) = (false, false)
get_static_argument_types(::InverseGamma) = [Float64, Float64]

export inv_gamma

########
# beta #
########

struct Beta <: Distribution{Float64} end

"""
    beta(alpha::Real, beta::Real)

Sample a `Float64` from a beta distribution.
"""
const beta = Beta()

function logpdf(::Beta, x::Real, alpha::Real, beta::Real)
    (alpha - 1) * log(x) + (beta - 1) * log1p(-x) - lbeta(alpha, beta)
end

function logpdf_grad(::Beta, x::Real, alpha::Real, beta::Real)
    if 0 <= x <= 1
        deriv_x = (alpha - 1) / x - (beta - 1) / (1 - x)
        deriv_alpha = log(x) - (digamma(alpha) - digamma(alpha + beta))
        deriv_beta = log1p(-x) - (digamma(beta) - digamma(alpha + beta))
    else
        error("x is outside of support: $x")
    end
    (deriv_x, deriv_alpha, deriv_beta)
end

function random(::Beta, alpha::Real, beta::Real)
    rand(Distributions.Beta(alpha, beta))
end

(::Beta)(alpha, beta) = random(Beta(), alpha, beta)

has_output_grad(::Beta) = true
has_argument_grads(::Beta) = (true, true)
get_static_argument_types(::Beta) = [Float64, Float64]

export beta

###############
# categorical #
###############

struct Categorical <: Distribution{Int} end

"""
    categorical(probs::AbstractArray{U, 1}) where {U <: Real}

Given a vector of probabilities `probs` where `sum(probs) = 1`, sample an `Int` `i` from the set {1, 2, .., `length(probs)`} with probability `probs[i]`.
"""
const categorical = Categorical()

function logpdf(::Categorical, x::Int, probs::AbstractArray{U,1}) where {U <: Real}
    log(probs[x])
end

function logpdf_grad(::Beta, x::Int, probs::AbstractArray{U,1})  where {U <: Real}
    grad = zeros(length(probs))
    grad[x] = 1.0
    (nothing, grad)
end

function random(::Categorical, probs::AbstractArray{U,1}) where {U <: Real}
    rand(Distributions.Categorical(probs))
end

(::Categorical)(probs) = random(Categorical(), probs)

has_output_grad(::Categorical) = false
has_argument_grads(::Categorical) = (true,)
get_static_argument_types(::Categorical) = [Vector{Float64}]

export categorical

####################
# uniform_discrete #
####################

struct UniformDiscrete <: Distribution{Int} end

"""
    uniform_discrete(low::Integer, high::Integer)

Sample an `Int` from the uniform distribution on the set {low, low + 1, ..., high-1, high}.
"""
const uniform_discrete = UniformDiscrete()

function logpdf(::UniformDiscrete, x::Int, low::Integer, high::Integer)
    d = Distributions.DiscreteUniform(low, high)
    Distributions.logpdf(d, x)
end

function logpdf_grad(::UniformDiscrete, x::Int, lower::Integer, high::Integer)
    (nothing, nothing, nothing)
end

function random(::UniformDiscrete, low::Integer, high::Integer)
    rand(Distributions.DiscreteUniform(low, high))
end

(::UniformDiscrete)(low, high) = random(UniformDiscrete(), low, high)

has_output_grad(::UniformDiscrete) = false
has_argument_grads(::UniformDiscrete) = (false, false)
get_static_argument_types(::UniformDiscrete) = [Int, Int]

export uniform_discrete

######################
# uniform_continuous #
######################

struct UniformContinuous <: Distribution{Float64} end

const uniform_continuous = UniformContinuous()

"""
    uniform(low::Real, high::Real)

Sample a `Float64` from the uniform distribution on the interval [low, high].
"""
const uniform = uniform_continuous

function logpdf(::UniformContinuous, x::Real, low::Real, high::Real)
    (x >= low && x <= high) ? -log(high-low) : -Inf
end

function logpdf_grad(::UniformContinuous, x::Real, low::Real, high::Real)
    inv_diff = 1. / (high-low)
    (0., inv_diff, -inv_diff)
end

function random(::UniformContinuous, low::Real, high::Real)
    rand() * (high - low) + low
end

(::UniformContinuous)(low, high) = random(UniformContinuous(), low, high)

has_output_grad(::UniformContinuous) = true
has_argument_grads(::UniformContinuous) = (true, true)
get_static_argument_types(::UniformContinuous) = [Float64, Float64]

export uniform_continuous, uniform

###########
# poisson #
###########

struct Poisson <: Distribution{Int} end

"""
    poisson(lambda::Real)

Sample an `Int` from the Poisson distribution with rate `lambda`.
"""
const poisson = Poisson()

function logpdf(::Poisson, x::Integer, lambda::Real)
    x * log(lambda) - lambda - lgamma(x+1)
end

function logpdf_grad(::Poisson, x::Integer, lambda::Real)
    error("Not implemented")
    (nothing, nothing)
end

function random(::Poisson, lambda::Real)
    rand(Distributions.Poisson(lambda))
end

(::Poisson)(lambda) = random(Poisson(), lambda)

has_output_grad(::Poisson) = false
has_argument_grads(::Poisson) = (false,)
get_static_argument_types(::Poisson) = [Float64]

export poisson

ext: jl

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
#
# Superblock
#

const SUPERBLOCK_SIGNATURE = htol(0x0a1a0a0d46444889) # UInt8[0o211, 'H', 'D', 'F', '\r', '\n', 0o032, '\n']

# https://www.hdfgroup.org/HDF5/doc/H5.format.html#FileMetaData
# Superblock (Version 2)
struct Superblock
    file_consistency_flags::UInt8
    base_address::Int64
    superblock_extension_address::RelOffset
    end_of_file_address::Int64
    root_group_object_header_address::RelOffset
end

Base.sizeof(::Union{Type{Superblock},Superblock}) =
    12+sizeof(RelOffset)*4+4

function Base.read(io::IO, ::Type{Superblock})
    cio = begin_checksum_read(io)

    # Signature
    signature = read(cio, UInt64)
    signature == SUPERBLOCK_SIGNATURE || throw(InvalidDataException())

    # Version
    version = read(cio, UInt8)
    version == 2 || throw(UnsupportedVersionException())

    # Size of offsets and size of lengths
    size_of_offsets = read(cio, UInt8)
    size_of_lengths = read(cio, UInt8)
    (size_of_offsets == 8 && size_of_lengths == 8) || throw(UnsupportedFeatureException())

    # File consistency flags
    file_consistency_flags = read(cio, UInt8)

    # Addresses
    base_address = read(cio, Int64)
    superblock_extension_address = read(cio, RelOffset)
    end_of_file_address = read(cio, Int64)
    root_group_object_header_address = read(cio, RelOffset)

    # Checksum
    cs = end_checksum(cio)
    read(io, UInt32) == cs || throw(InvalidDataException())

    Superblock(file_consistency_flags, base_address, superblock_extension_address,
                end_of_file_address, root_group_object_header_address)
end

function Base.write(io::IO, s::Superblock)
    cio = begin_checksum_write(io, 8+4+4*sizeof(RelOffset))
    write(cio, SUPERBLOCK_SIGNATURE::UInt64)    # Signature
    write(cio, UInt8(2))                        # Version
    write(cio, UInt8(8))                        # Size of offsets
    write(cio, UInt8(8))                        # Size of lengths
    write(cio, s.file_consistency_flags::UInt8)
    write(cio, s.base_address::Int64)
    write(cio, s.superblock_extension_address::RelOffset)
    write(cio, s.end_of_file_address::Int64)
    write(cio, s.root_group_object_header_address::RelOffset)
    write(io, end_checksum(cio))
end


text:
using Documenter, XLSXasJSON

# copy coverage report to build folder for Documentor.jl to take along to `gh-pages` branch
function copy_coverage()
  source = joinpath(@__DIR__, "src/coverage")
  target = joinpath(@__DIR__, "build/coverage")
  @info "Copy $source to $target"

  cp(source, target; force = true)
  nothing
end

makedocs(
   modules = [XLSXasJSON],
 checkdocs = :all,
   authors = "YongHee Kim",
  sitename = "XLSXasJSON.jl",
  pages = [ "Home" => "index.md",
            "Tutorial" => "tutorial.md",
            "API Reference" => "api.md" ]
)

deploydocs(
    repo   = "github.com/YongHee-Kim/XLSXasJSON.jl.git",
    target = "build",
    deps   = copy_coverage(),
    make   = nothing
)


text:
# -*- coding: UTF-8 -*-
try:
    import biggles # Requires biggles: http://biggles.sourceforge.net/
    has_biggles = True
except ImportError:
    print "Biggles is not installed. If you wish to automatically plot some nice statistics please install it: http://biggles.sourceforge.net/"
    has_biggles = False

try:
    import pydot   # Requires PyDot: http://code.google.com/p/pydot/downloads/list
    has_pydot = True
except:
    print "PyDot is not installed. If you wish to generate a graphical representation of the resulting neural network, please install it: http://code.google.com/p/pydot/"
    has_pydot = False

import random

def draw_net(chromosome, id=''):
    ''' Receives a chromosome and draws a neural network with arbitrary topology. '''
    output = 'digraph G {\n  node [shape=circle, fontsize=9, height=0.2, width=0.2]'

    # subgraph for inputs and outputs
    output += '\n  subgraph cluster_inputs { \n  node [style=filled, shape=box] \n    color=white'
    for ng in chromosome.node_genes:
        if ng.type== 'INPUT':
            output += '\n    '+str(ng.id)
    output += '\n  }'

    output += '\n  subgraph cluster_outputs { \n    node [style=filled, color=lightblue] \n    color=white'
    for ng in chromosome.node_genes:
        if ng.type== 'OUTPUT':
            output += '\n    '+str(ng.id)
    output += '\n  }'
    # topology
    for cg in chromosome.conn_genes:
        output += '\n  '+str(cg.innodeid)+' -> '+str(cg.outnodeid)
        if cg.enabled is False:
            output += ' [style=dotted, color=cornflowerblue]'

    output += '\n }'

    if has_pydot:
        g = pydot.graph_from_dot_data(output)
        g.write('phenotype'+id+'.svg', prog='dot', format='svg')
    else:
        print 'You do not have the PyDot package.'

def draw_ff(chromosome):
    ''' Draws a feedforward neural network '''

    output = 'digraph G {\n  node [shape=circle, fontsize=9, height=0.2, width=0.2]'

    # subgraph for inputs and outputs
    output += '\n  subgraph cluster_inputs { \n  node [style=filled, shape=box] \n    color=white'
    for ng in chromosome.node_genes:
        if ng.type== 'INPUT':
            output += '\n    '+str(ng.id)
    output += '\n  }'

    output += '\n  subgraph cluster_outputs { \n    node [style=filled, color=lightblue] \n    color=white'
    for ng in chromosome.node_genes:
        if ng.type== 'OUTPUT':
            output += '\n    '+str(ng.id)
    output += '\n  }'
    # topology
    for cg in chromosome.conn_genes:
        output += '\n  '+str(cg.innodeid)+' -> '+str(cg.outnodeid)
        if cg.enabled is False:
            output += ' [style=dotted, color=cornflowerblue]'

    output += '\n }'

    if has_pydot:
        g = pydot.graph_from_dot_data(output)
        g.write('feedforward.svg', prog='dot', format='svg')
    else:
        print 'You do not have the PyDot package.'

def plot_stats(stats):
    ''' Plots the population's average and best fitness. '''
    if has_biggles:
        generation = [i for i in xrange(len(stats[0]))]

        fitness = [c.fitness for c in stats[0]]
        avg_pop = [avg for avg in stats[1]]

        plot = biggles.FramedPlot()
        plot.title = "Population's average and best fitness"
        plot.xlabel = r"Generations"
        plot.ylabel = r"Fitness"

        plot.add(biggles.Curve(generation, fitness, color="red"))
        plot.add(biggles.Curve(generation, avg_pop, color="blue"))

        #plot.show() # X11
        plot.write_img(600, 300, 'avg_fitness.svg')
        # width and height doesn't seem to affect the output!
    else:
        print 'You dot not have the Biggles package.'

def plot_spikes(spikes):
    ''' Plots the trains for a single spiking neuron. '''
    if has_biggles:
        time = [i for i in xrange(len(spikes))]

        plot = biggles.FramedPlot()
        plot.title = "Izhikevich's spiking neuron model"
        plot.ylabel = r"Membrane Potential"
        plot.xlabel = r"Time (in ms)"

        plot.add(biggles.Curve(time, spikes, color="green"))
        plot.write_img(600, 300, 'spiking_neuron.svg')
        # width and height doesn't seem to affect the output!
    else:
        print 'You dot not have the Biggles package.'

def plot_species(species_log):
    ''' Visualizes speciation throughout evolution. '''
    if has_biggles:
        plot = biggles.FramedPlot()
        plot.title = "Speciation"
        plot.ylabel = r"Size per Species"
        plot.xlabel = r"Generations"
        generation = [i for i in xrange(len(species_log))]

        species = []
        curves = []

        for gen in xrange(len(generation)):
            for j in xrange(len(species_log), 0, -1):
                try:
                    species.append(species_log[-j][gen] + sum(species_log[-j][:gen]))
                except IndexError:
                    species.append(sum(species_log[-j][:gen]))
            curves.append(species)
            species = []

        s1 = biggles.Curve(generation, curves[0])

        plot.add(s1)
        plot.add(biggles.FillBetween(generation, [0]*len(generation), generation, curves[0], color=random.randint(0,90000)))

        for i in range(1, len(curves)):
            c = biggles.Curve(generation, curves[i])
            plot.add(c)
            plot.add(biggles.FillBetween(generation, curves[i-1], generation, curves[i], color=random.randint(0,90000)))

        plot.write_img(1024, 800, 'speciation.svg')

    else:
        print 'You dot not have the Biggles package.'


text:
function vecnorm(A::Array, dim::Int)
    B = sum(x -> x^2, A, dims = dim)
    B .= sqrt.(B)
end

function vecnorm2(A::Array, dim::Int)
    A .= A.^2
    B = sum(A, dims = dim)
    B .= sqrt.(B)
end

function kernel(input::Array, dim::Int, mode::Char, mDict::Array)

    if mode == '1'
        output = [input[1] * ones(dim - length(input), 1) , input];
    elseif mode == '2'
        dist_var = 0.075;
        num_input = length(input);
        output = zeros(dim, 1);
        idx = 1;
        # TODO: Optimize this
        for k_idx = 1:num_input
           for i_idx = k_idx + 1 : num_input
               output[idx] = exp(-(input[k_idx] - input[i_idx])/(2 * dist_var));
               idx = idx + 1;
           end
        end
        output = output[1:dim];
    elseif mode == '3'
        # TODO: Implement Vecnorm - Done
        if length(mDict) == 6
            mDict = [mDict zeros(Float64, (6, 1))]
        end
        add_dim = vecnorm(input .- mDict[:, 2:end], 1)
        output = [add_dim'; input];

        output = output[1:dim, :];
    elseif mode == '4'
        dist_var = 0.75;
        c = size(mDict, 2);
        if (c < 2)
            add_dim = 0;
        else
            add_dim = vecnorm(input .- mDict[:, 1: 0 + dim - length(input)], 1);
        end
        # % covariance_matrix = eye(length(add_dim));
        new_vec = 1 / (dist_var * sqrt(2*pi)) * exp( - add_dim .^ 2 / (2 * dist_var));
        output = [new_vec'; input];
        output = output[1:dim];
    end

    return output
end


text:
using Test, SafeTestsets

@time begin
@time @safetestset "sumBetweenElementsTest.jl" begin include("./laboratorio3/sumBetweenElementsTest.jl") end
@time @safetestset "freqTwo.jl" begin include("./laboratorio3/freqTwoTest.jl") end
@time @safetestset "sumOfMult.jl" begin include("./laboratorio3/sumOfMultTest.jl") end
@time @safetestset "countInc.jl" begin include("./laboratorio3/countIncTest.jl") end
@time @safetestset "countAAA.jl" begin include("./laboratorio3/countAAATest.jl") end
end


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 2.0
x, y = .1, .1
r = l1 / 4
h = r / 10
box = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))
cyl = Cylinder(r, h, 2 * r, color = RGBA(1., 0., 0.))
box2 = Box(h, 2r, 2r, 2r, color = RGBA(1., 0., 1.))

p0 = [0;l1 / 2;0]
p1 = -p0
p2 = [0;r;0]
p3 = -p2

# Initial orientation
q1 = UnitQuaternion(RotY(pi / 2))

# Links
origin = Origin{Float64}()
link1 = Body(box)
link2 = Body(cyl)
link3 = Body(box2)

# Constraints
joint0to23 = EqualityConstraint(Revolute(origin, link2, ez; p1=p3 + p1, p2=zeros(3)), Revolute(origin, link3, ex; p1=p2, p2=zeros(3)))
joint1to23 = EqualityConstraint(Spherical(link1, link2; p1=p1, p2=p3), Spherical(link1, link3; p1=p0, p2=p3))

links = [link1;link2;link3]
constraints = [joint0to23;joint1to23]
shapes = [box,cyl,box2]

mech = Mechanism(origin, links, constraints, g = 0., shapes = shapes)
setPosition!(link1,x = -p0)
setPosition!(link2,x = p3 + p1)
setPosition!(link3,x = p2)
setVelocity!(link2,ω = [0;0;2.])
setVelocity!(link3,ω = [2.;0;0])


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 19731e2c-a09f-11ec-053c-e5c07ed14eeb
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using Revise, VOTables, DataFrames, HTTP, HTTP.URIs, PlutoUI, CSV
end

# ╔═╡ e4a3469b-e3cd-4d69-bc0e-7730037dffcf
md"""
## Simbad -- [sim-script](http://simbad.u-strasbg.fr/guide/sim-url.htx)

Simbad has its own scripting language 😮
"""

# ╔═╡ ad392537-8502-4304-b722-c79fbee6edde
script = """
output console=off script=off
votable {
    MAIN_ID
    RA(s)
    DEC(s)
    PLX(V)
    FLUX(V)
    FLUX(G)
    FLUX(H)
}
votable open
set radius 5m
query around HD 32297
votable close
""" |> escapeuri;

# ╔═╡ bf63f661-6983-4b59-b6c6-6527f165aba6
res = HTTP.get("https://simbad.u-strasbg.fr/simbad/sim-script";
	query = "script=$script"
);

# ╔═╡ fff97d31-471d-4523-8190-8ce05a30eba8
doc = VOTable(String(res.body))

# ╔═╡ 274013c1-c691-46a5-bcf0-8e3d6077fdff
tbl = DataFrame(doc[1])

# ╔═╡ 60dfbacc-42fd-40a2-a0ac-816c287cea11
names(doc[1])

# ╔═╡ aef54862-2292-4949-b3e7-333667902c9d
VOTables.units(doc[1])

# ╔═╡ caf98e0a-aa92-467d-a2eb-b646662e243b
md"""
## Exoplanet archive -- [TAP](https://exoplanetarchive.ipac.caltech.edu/docs/TAP/usingTAP.html)

Follows the [ADQL](https://www.ivoa.net/documents/latest/ADQL.html) spec used by Vizier, Simbad, others?

Let's start by downloading some data in VOTable format:
"""

# ╔═╡ fbad1e29-e582-4228-9188-292f096d9862
query_exoarchive = """
SELECT TOP 10
hostname, pl_name, tic_id, rastr, decstr, sy_jmag
FROM pscomppars
WHERE tran_flag=1
""" |> escapeuri;

# ╔═╡ 6410b817-b03e-41ff-988a-404e690f0b4c
res_exoarchive = HTTP.get("https://exoplanetarchive.ipac.caltech.edu/TAP/sync";
	query = "query=$(query_exoarchive)"
);

# ╔═╡ b39f8132-bcd4-4ef1-9615-20a2a0ef4a9a
doc_exoarchive = VOTable(String(res_exoarchive.body))

# ╔═╡ 99fb246b-2954-444a-bf62-72abe20ec550
tbl_exoarchive = DataFrame(doc_exoarchive[1])

# ╔═╡ 6c5f7bc1-8475-4e0c-b9bc-065ff3b7dad6
VOTables.units(doc_exoarchive[1])

# ╔═╡ fe26f2ab-df63-4e6a-9200-0cda0b3dbaf7
md"""
!!! question
	Should we expect unit determination to just work™ here?
"""

# ╔═╡ 56d78381-48a7-4f8a-82e6-af635738eab5
md"""
Besides VOTable, we can also download the data as a CSV and handle it as a DataFrame directly:
"""

# ╔═╡ 21f1b4df-1f10-4ca8-aee1-3e4c09b3482f
res_exoarchive_df = HTTP.get("https://exoplanetarchive.ipac.caltech.edu/TAP/sync";
	query = "query=$(query_exoarchive)&format=csv"
);

# ╔═╡ 1c30a30e-f999-44c8-a2c6-be6fdc6fbc04
CSV.read(res_exoarchive_df.body, DataFrame)

# ╔═╡ 8adb03b9-177e-4d7a-b972-b585774f1e5c
TableOfContents()

# ╔═╡ Cell order:
# ╠═19731e2c-a09f-11ec-053c-e5c07ed14eeb
# ╟─e4a3469b-e3cd-4d69-bc0e-7730037dffcf
# ╠═ad392537-8502-4304-b722-c79fbee6edde
# ╠═bf63f661-6983-4b59-b6c6-6527f165aba6
# ╠═fff97d31-471d-4523-8190-8ce05a30eba8
# ╠═274013c1-c691-46a5-bcf0-8e3d6077fdff
# ╠═60dfbacc-42fd-40a2-a0ac-816c287cea11
# ╠═aef54862-2292-4949-b3e7-333667902c9d
# ╟─caf98e0a-aa92-467d-a2eb-b646662e243b
# ╠═fbad1e29-e582-4228-9188-292f096d9862
# ╠═6410b817-b03e-41ff-988a-404e690f0b4c
# ╠═b39f8132-bcd4-4ef1-9615-20a2a0ef4a9a
# ╠═99fb246b-2954-444a-bf62-72abe20ec550
# ╠═6c5f7bc1-8475-4e0c-b9bc-065ff3b7dad6
# ╟─fe26f2ab-df63-4e6a-9200-0cda0b3dbaf7
# ╟─56d78381-48a7-4f8a-82e6-af635738eab5
# ╠═21f1b4df-1f10-4ca8-aee1-3e4c09b3482f
# ╠═1c30a30e-f999-44c8-a2c6-be6fdc6fbc04
# ╟─8adb03b9-177e-4d7a-b972-b585774f1e5c


text:
# 2D FWI on Overthrust model using minConf library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: May 2017
#

using PyCall, PyPlot, HDF5, opesciSLIM.TimeModeling, opesciSLIM.SLIM_optim, SeisIO

# Load starting model
n,d,o,m0 = h5open("../data/overthrust_2D_initial_model.h5","r") do file
	read(file, "n", "d", "o", "m0")
end

# Set up model structure
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)	# need n,d,o as tuples and m0 as array

# Bound constraints
vmin = ones(Float32,model0.n) + 0.3f0
vmax = ones(Float32,model0.n) + 5.5f0

# Slowness squared [s^2/km^2]
mmin = vec((1f0./vmax).^2)
mmax = vec((1f0./vmin).^2)

# Load data
block = segy_read("../data/overthrust_2D.segy")
dobs = joData(block)

# Set up wavelet
src_geometry = Geometry(block; key="source", segy_depth_key="SourceDepth")
wavelet = source(src_geometry.t[1],src_geometry.dt[1],0.008f0)	# 8 Hz wavelet
q = joData(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
srand(1)	# set seed of random number generator
niterations = 20
batchsize = 20
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)],2),model0.n)	# bound projection operator
fTerm = 1f-2
gradTerm = 1f-4
num_feval = 0
max_feval = 100
fval = 0
fvals = zeros(Float32,niterations)
dsub = []
qsub = []

# evaluate FWI misfit function for current model
function misfit(x)
	global num_feval; num_feval += 1
	model_new = deepcopy(model0)
	model_new.m = proj(reshape(x,model0.n))
	return f,g = fwi_objective(model_new,qsub,dsub)
end

# Backtracking line search with Armijo-Goldstein condition
function backtracking_linesearch(objective, x, grad; f_curr=[], alpha=1f0, tau=.5f0, c=.5f0, maxiter=10)

	# Parameters
	p = -grad/norm(grad)	# normalized descent direction
	m = dot(p,grad)
	t = -c*m
	isempty(f_curr) && (f_curr = objective(x)[1])	# current function value
	f_new = objective(x + alpha*p)[1]	# function value with initial step size
	iteration = 0

	# Line search
	while(f_curr - f_new < alpha*t && iteration < maxiter)
		println("	Iter LS: ", iteration, "; ",f_curr - f_new," >= ", alpha*t, "; alpha: ",alpha)
		alpha *= tau
		f_new = objective(x + alpha*p)[1]
		iteration += 1
	end
	return alpha*p
end

# Main loop
for j=1:niterations

	# select current batch
	idx = randperm(dobs.nsrc)[1:batchsize]
	dsub = subsample(dobs,idx)
	qsub = subsample(q,idx)

	# get fwi objective function value and gradient
	fval, gradient = fwi_objective(model0,qsub,dsub)
	println("FWI iteration no: ",j,"; function value: ",fval)

	# linesearch
	step = backtracking_linesearch(misfit, vec(model0.m), gradient; f_curr=fval, alpha=10.f0, tau=.1f0, c=1e-4)

	# Update model and bound projection
	model0.m = proj(model0.m + reshape(step,model0.n))

	# termination criteria
	fvals[j] = fval
	if fval <= fTerm || norm(gradient) <= gradTerm || num_feval > max_feval
		break
	end
end

# Save final velocity model, function value and history
h5open("fwi_overthrust_2D_linesearch.h5", "w") do file
	write(file, "x", sqrt.(1f0./reshape(model0.m,model0.n)), "fhistory",fvals)
end

# Plot results
figure(); imshow(sqrt.(1f0./m0)',cmap="jet",vmin=1.5,vmax=6.); title("Initial model")
figure(); imshow(sqrt.(1f0./reshape(model0.m,model0.n))',cmap="jet",vmin=1.5,vmax=6.); title("FWI")
figure(); plot(fvals); title("Function value")

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
module  TestConditionals

using PopGenCore
using Test

x = @nancycats

@testset "Conditionals.jl" begin
    @testset "biallelic" begin
        @test isbiallelic(x) == false
        @test isbiallelic(fill((1,2), 10)) == true
        @test isbiallelic(x.genodata) == false
    end

    @testset "homozygous" begin
        @test ishom(missing) == false
        @test PopGenCore._ishom(missing) === missing
        @test ishom((1,2)) == false
        @test ishom((1,1)) == true
        @test ishom((1,1,2)) == false
        @test ishom((2,2,2)) == true
        @test all(ishom(x.genodata.genotype[1:5]) .== [false, false, false, false, false])
        @test all(ishom(skipmissing(x.genodata.genotype[1:5])) .== [false, false, false])
        @test ishom((1,1), 1) == true
        @test ishom((1,2), 1) == false
        @test ishom((1,1), 2) == false
        @test ishom(missing, 2) == false
        @test PopGenCore._ishom(missing, 2) === missing
    end

    @testset "heterozygous" begin
        @test ishet(missing) == false
        @test PopGenCore._ishet(missing) === missing
        @test ishet((1,2)) == true
        @test ishet((1,1)) == false
        @test ishet((1,1,2)) == true
        @test ishet((2,2,2)) == false
        @test all(ishet(x.genodata.genotype[1:5]) .=== [false, false, true, true, true])
        @test all(ishet(skipmissing(x.genodata.genotype[1:5])) .== [true, true, true])
        @test ishet((1,1), 1) == false
        @test ishet((1,2), 1) == true
        @test ishet((1,1), 2) == false
        @test ishet(missing, 2) == false
        @test PopGenCore._ishet(missing, 2) === missing
    end
end

end # module


text:
fs2()

@testset "fig/rpath" begin
    gotd()

    mkpath(joinpath(td, "__site", "assets", "index"))
    write(joinpath(td, "__site", "assets", "index", "baz.png"), "baz")

    mkpath(joinpath(td, "__site", "assets", "blog", "kaggle"))
    write(joinpath(td, "__site", "assets", "blog", "kaggle", "foo.png"), "foo")

    s = raw"""
        @def fd_rpath = "index.md"
        ABC
        \fig{./baz.png}
        """ |> fd2html_td
    @test isapproxstr(s, """

ABC

        """)

    s = raw"""
        @def fd_rpath = "index.md"
        ABC
        \fig{./unknown.png}
        """ |> fd2html_td
    @test isapproxstr(s, """

ABC
// Image matching '/assets/index/unknown.png' not found. //

        """)

    s = raw"""
        @def fd_rpath = "blog/kaggle/index.md"
        ABC
        \fig{./foo.png}
        """ |> fd2html_td
    @test isapproxstr(s, """

ABC

        """)

    s = raw"""
        @def fd_rpath = "blog/kaggle/index.md"
        ABC
        \fig{./baz.png}
        """ |> fd2html_td
    @test isapproxstr(s, """

ABC
// Image matching '/assets/blog/kaggle/baz.png' not found. //

        """)
end


text:
using RDatasets, DataArrays, DataFrames, Gadfly

set_default_plot_size(6inch, 3inch)

df = dataset("plm", "Cigar")
df = df[df[:State] .== 1, :]
ymin = df[:Sales] .- 20*rand()
ymax = df[:Sales] .+ 20*rand()

plot(df, x="Year", y="Sales", ymin=ymin, ymax=ymax,
     Geom.bar, Geom.errorbar)


text:
module ArchGDAL

    using Dates
    using GDAL: GDAL
    using GeoFormatTypes: GeoFormatTypes
    using GeoInterface: GeoInterface
    using Tables: Tables

    const GFT = GeoFormatTypes

    include("utils.jl")
    include("types.jl")
    include("driver.jl")
    include("gcp.jl")
    include("spatialref.jl")
    include("dataset.jl")
    include("raster/rasterband.jl")
    include("raster/rasterio.jl")
    include("raster/array.jl")
    include("raster/rasterattributetable.jl")
    include("raster/colortable.jl")
    include("ogr/geometry.jl")
    include("ogr/feature.jl")
    include("ogr/featurelayer.jl")
    include("ogr/featuredefn.jl")
    include("ogr/fielddefn.jl")
    include("ogr/styletable.jl")
    include("utilities.jl")
    include("context.jl")
    include("base/iterators.jl")
    include("base/display.jl")
    include("tables.jl")
    include("geointerface.jl")
    include("convert.jl")

    mutable struct DriverManager
        function DriverManager()
            drivermanager = new()
            GDAL.gdalallregister()
            finalizer((dm,) -> GDAL.gdaldestroydrivermanager(), drivermanager)
            return drivermanager
        end
    end

    const DRIVER_MANAGER = Ref{DriverManager}()

    function __init__()
        DRIVER_MANAGER[] = DriverManager()
    end

end # module


text:
import JITrench
using JITrench: Variable, Model, Linear, sigmoid, SGD,  mean_squared_error, cleargrads!, backward!, do_optimize!
using Printf
using Random
using Plots

Random.seed!(10)

ENV["GKSwstype"] = "nul"

function generate_dataset(N)
    x = rand(N, 1)
    y = sin.(2π .* x) .+ (rand(N, 1) .* 0.5)
    return Variable(x),  Variable(y)
end

function train(model, x, y; n_iters=10000, log_interval=200, lr=1e-1)
    optimizer = SGD(layers, lr=1e-1)
    anim = Animation()
    history = []
    JITrench.plot_model(model, x)
    for iter in 1:n_iters
        y_pred = model(x)
        loss = mean_squared_error(y, y_pred)
        push!(history, loss.values)
        cleargrads!(model, layers, skip_uninit=true)
        backward!(loss)
        do_optimize!(model, optimizer)
        if (iter - 1) % log_interval == 0
            @printf "[iters] %4i [loss] %.2f\n" iter loss.values
            val = reshape(collect(0:0.01:1), :, 1)
            y_pred = model(val).values
            p1 = scatter(x.values, y.values, label="data")
            plot!(p1, 0:0.01:1, y_pred, title="[mse : $(repr(loss.values)[1:10])]", label="predict")
            p2 = plot(history, title="learning curve", xlims=(0, n_iters))
            plt = plot(p1, p2, size=(800, 1000), layout=(2, 1))
            frame(anim, plt)
        end
    end
    return history, anim
end

mutable struct MLP <: Model
    l1 :: Linear
    l2 :: Linear
    MLP(hidden_dim, out_dim) = new(Linear(hidden_dim), Linear(out_dim))
end

layers(model::MLP) = (model.l1, model.l2)

(model::MLP)(x) = x |> model.l1 .|> sigmoid |> model.l2

model = MLP(10, 1)

history, anim = train(model, generate_dataset(100)...)
gif(anim, "fitting_history.gif", fps=8)


text:
using Test
using SuffixArrays

function test_suffix(args)
    for file in args
        data = codeunits(read(file, String))
        t = @elapsed suffixes = suffixsort(data, 0)
        println("Sorting '$file' took: $t")
        @test sufcheck(data, suffixes) == 0
    end
end

function sufcheck(T, SA)
    n = length(T)
    n == 0 && (println("Done."); return 0)
    n < 0 && (println("Invalid length $n"); return -1)
    C = zeros(Int, 256)
    for i = 1:n
        if SA[i] < 0 || n <= SA[i]
            println("Out of range $n")
            println("SA[$i] = $(SA[i])")
            return -2
        end
    end
    for i = 2:n
        if T[SA[i-1]+1] > T[SA[i]+1]
            println("Suffixes in wrong order")
            println("T[SA[$(i-1)]+1] = $(T[SA[(i-1)]+1])")
            println("T[SA[$i]+1] = $(T[SA[i]+1])")
            return -3
        end
    end
    for i = 1:n
        C[Int(T[i])+1] += 1
    end
    p = 0
    for i = 1:256
        t = C[i]
        C[i] = p
        p += t
    end
    q = C[Int(T[n])+1]
    C[Int(T[n])+1] += 1
    for i = 1:n
        p = SA[i]
        if 0 < p
            p -= 1
            c = T[p+1]
            t = C[Int(c)+1]
        else
            p = n - 1
            c = T[p+1]
            t = q
        end
        if t < 0 || p != SA[t+1]
            println("Suffixes in wrong position")
            return -4
        end
        if t != q
            C[Int(c)+1] += 1
            if n <= C[Int(c)+1] || T[SA[C[Int(c)+1]+1]+1] != c
                C[Int(c)+1] = -1
            end
        end
    end
    println("Done.")
    return 0
end

function initwalk(dir, files)
    files = walkdir("$dir/src", files)
    files = walkdir("$dir/test", files)
    files
end

function walkdir(dir, files)
    t = readdir(dir)
    for f in t
        f == ".git" && continue
        j = joinpath(dir, f)
        if isdir(j)
            append!(files, walkdir(j, files))
        else
            push!(files, j)
        end
    end
    return unique(files)
end

@testset "source files" begin
    files = initwalk(dirname(dirname(@__FILE__)), [])
    test_suffix(files)
end

@testset "UTF-8 strings" begin
    s = "¡Hello, 😄 world!"
    sa = suffixsort(s)
    suffixes = [String(codeunits(s)[i:end]) for i in sa]
    @test issorted(suffixes)
end

## define a simple UTF-16 string type ##

struct UTF16 <: AbstractString
    codeunits::Vector{UInt16}
end
UTF16(s::String) = UTF16(Base.transcode(UInt16, s))

Base.codeunits(s::UTF16) = s.codeunits
Base.ncodeunits(s::UTF16) = length(s.codeunits)
Base.isvalid(s::UTF16, i::Int) = isvalid(iterate(s, i)[1])
Base.isless(s::UTF16, t::UTF16) = s.codeunits < t.codeunits

function Base.iterate(s::UTF16, i::Int=1)
    i ≤ length(s.codeunits) || return
    u = s.codeunits[i]
    0xD800 ≤ u ≤ 0xDBFF || return Char(u), i+1
    # otherwise is a high surrogate
    v = s.codeunits[i+1]
    # not followed by low surrogate
    0xDC00 ≤ v ≤ 0xDFFF || return Char(u), i+1
    # u, v are high/low surrogate pair
    Char(0x10000 + (UInt32(u & 0x03ff) << 10) | (v & 0x03ff)), i+2
end

@testset "UTF-16 strings" begin
    s = UTF16("¡Hello, 😄 world!")
    sa = suffixsort(s)
    suffixes = [UTF16(codeunits(s)[i:end]) for i in sa]
    @test issorted(suffixes)
end

function commonprefixlen(s1, s2)
    h = 0
    maxh = min(length(s1), length(s2))
    for i in 1:maxh
        if s1[i] != s2[i]
            break
        end
        h += 1
    end
    h
end

@testset "Longest common prefix" begin
    s = rand(0x00:0xff, 100)
    sa = suffixsort(s)
    suff = [s[i:end] for i in sa]
    lcparr = lcp(sa, s)
    # LCP the hard way
    lcpref = [commonprefixlen(suff[i], suff[i+1]) for i in 1:length(sa)-1]
    @test lcparr[1] == 0
    @test lcparr[2:end] == lcpref
    @test typeof(lcparr) == typeof(sa)

    # retest with base != 1
    base = 0
    sa = suffixsort(s, base)
    lcparr = lcp(sa, s, base)
    @test lcparr[1] == 0
    @test lcparr[2:end] == lcpref
    @test typeof(lcparr) == typeof(sa)

    # sequence where common prefix reaches the end
    s = [0x01, 0x02, 0x03, 0x04, 0x01, 0x02, 0x03]
    sa = suffixsort(s)
    suff = [s[i:end] for i in sa]
    lcparr = lcp(sa, s)
    # LCP the hard way
    lcpref = [commonprefixlen(suff[i], suff[i+1]) for i in 1:length(sa)-1]
    @test lcparr[1] == 0
    @test lcparr[2:end] == lcpref
    @test typeof(lcparr) == typeof(sa)
end

@testset "Issue #15 fix" begin
    s = join(rand('a':'z', 10000)) * '\$'
    sa = suffixsort(s)
    suffixes = [String(codeunits(s)[i:end]) for i in sa]
    @test issorted(suffixes)

    utext = rand(0x0001:0x0020, 4000)
    sa = suffixsort(utext)
    suffixes = [utext[i:end] for i in sa]
    @test issorted(suffixes)
end


text:
using SpecialFunctions

d = PGeneralizedGaussian() # mean zero, scale √2, shape 2.

# PDF and CDF values from R using the same default
# mean scale and shape parameters. Format is [x, pdf, cdf].
test = [
    -2.0000 0.053990966513188 0.0227501319481792;
    -1.9596 0.0584867753200998 0.0250212806912653;
    -1.9192 0.0632536242130714 0.0274795137708557;
    -1.8788 0.0682974223520507 0.0301359035293964;
    -1.8384 0.0736231464260287 0.0330017425665818;
    -1.798 0.0792347326126488 0.0360885038859181;
    -1.7576 0.0851349688324957 0.0394077966787751;
    -1.7172 0.091325388414829 0.0429713177801418;
    -1.6768 0.0978061663608371 0.0467907988768448;
    -1.6364 0.104576019449044 0.0508779495981413;
    -1.596 0.111632111473603 0.0552443966698443;
    -1.5556 0.118969964938241 0.0599016193659873;
    -1.5152 0.126583380545089 0.0648608815458657;
    -1.4748 0.134464365817196 0.0701331606184371;
    -1.4344 0.142603074175027 0.0757290738298424;
    -1.394 0.1509877557498 0.0816588023224809;
    -1.3536 0.159604721159404 0.0879320134647649;
    -1.3132 0.168438319395545 0.0945577819987581;
    -1.2728 0.177470930873585 0.101544510597386;
    -1.2324 0.186682976579551 0.108899850463094;
    -1.192 0.196052944112626 0.116630622634895;
    -1.1516 0.20555743126705 0.124742740699938;
    -1.1112 0.215171207626101 0.133241135628347;
    -1.0708 0.224867294454356 0.142129683465442;
    -1.0304 0.234617062974812 0.151411136623046;
    -0.99 0.244390350907 0.161087059510831;
    -0.9496 0.254155596923586 0.171157769239293;
    -0.9092 0.263879992459047 0.181622282107569;
    -0.8688 0.273529650077933 0.192478266561878;
    -0.8284 0.283069787385305 0.203722003273787;
    -0.788 0.292464925241515 0.215348352941949;
    -0.7476 0.301679098831105 0.227350732366586;
    -0.7072 0.31067607993471 0.239721099283345;
    -0.6668 0.319419608566778 0.252449946372615;
    -0.6264 0.327873631974065 0.265526304782692;
    -0.586 0.336002548843259 0.278937757421185;
    -0.5456 0.343771456443474 0.292670462179526;
    -0.5052 0.351146398333471 0.306709185161556;
    -0.4648 0.358094610196201 0.321037343889939;
    -0.4244 0.364584761326559 0.335637060364848;
    -0.384 0.370587189293274 0.350489223749189;
    -0.3436 0.376074125323367 0.36557356235494;
    -0.3032 0.38101990801804 0.380868724507227;
    -0.2628 0.385401183101656 0.396352367767935;
    -0.2224 0.389197087030118 0.412001255910218;
    -0.182 0.392389412439884 0.42779136295047;
    -0.1416 0.394962753602334 0.443697983466388;
    -0.1012 0.396904630257822 0.459695848359751;
    -0.0608 0.398205588436696 0.475759245161512;
    -0.0204 0.398859277127687 0.491862141925561;
    0.02 0.398862499923666 0.507978313716902;
    0.0604 0.398215241057149 0.524081470670474;
    0.1008 0.396920665528795 0.540145386578997;
    0.1412 0.394985093327439 0.556144026962162;
    0.1816 0.392417948035498 0.572051675575411;
    0.222 0.38923168040452 0.587843058334322;
    0.2624 0.385441667768147 0.603493463659977;
    0.3028 0.381066090429489 0.618978858291238;
    0.3432 0.376125786413223 0.634275997660943;
    0.3836 0.370644086205869 0.649362529994026;
    0.424 0.364646629317593 0.664217093355429;
    0.4644 0.358161164682611 0.678819404953587;
    0.5048 0.351217337070506 0.693150342089967;
    0.5452 0.343846461805699 0.707192014235536;
    0.5856 0.336081290185479 0.720927825809832;
    0.626 0.327955768047694 0.734342529336241;
    0.6664 0.319504789967085 0.747422268746731;
    0.7068 0.310763951554577 0.760154612709528;
    0.7472 0.301769302297434 0.772528577952478;
    0.7876 0.292557101311255 0.784534642652148;
    0.828 0.283163578279136 0.796164750052605;
    0.8684 0.273624701730973 0.807412302567402;
    0.9088 0.263975956669541 0.81827214670236;
    0.9492 0.254252133382257 0.828740549214512;
    0.9896 0.244487129091676 0.838815164993143;
    1.03 0.234713763897012 0.848494997211656;
    1.0704 0.224963612246627 0.857780350353393;
    1.1108 0.215266850961145 0.866672776760204;
    1.1512 0.205652124601902 0.875175017389245;
    1.1916 0.196146428753372 0.883290937490971;
    1.232 0.186775011564286 0.891025457939793;
    1.2724 0.177561293673426 0.898384482958291;
    1.3128 0.168526806435554 0.905374824976726;
    1.3532 0.159691148163201 0.912004127362079;
    1.3936 0.151071957913399 0.918280785735546;
    1.434 0.142684906177067 0.924213868574892;
    1.4744 0.13454370167416 0.929813037768908;
    1.5148 0.126660113321382 0.935088469756236;
    1.5552 0.119044006322032 0.940050777840661;
    1.5956 0.111703391230226 0.944710936230536;
    1.636 0.104644484764399 0.949080206301951;
    1.6764 0.0978717810877164 0.953170065534611;
    1.7168 0.091388132235354 0.956992139516678;
    1.7572 0.085194836349956 0.960558137361137;
    1.7976 0.0792917323859474 0.963879790822026;
    1.838 0.0736772999596982 0.966968797345075;
    1.8784 0.0683487630544158 0.969836767234444;
    1.9188 0.0633021963346127 0.97249517506594;
    1.9592 0.0585326328834392 0.97495531542798;
    1.9996 0.0540341722453821 0.977228263024933;
]
# CDF test.
for i=1:size(test, 1)
    @test cdf(d, test[i, 1]) ≈ test[i, 3]
end
# PDF test.
for i=1:size(test, 1)
    @test pdf(d, test[i, 1]) ≈ test[i, 2]
end
@test mean(d) ≈ 0
@test median(d) ≈ 0
@test mode(d) ≈ 0
@test var(d) ≈ 1 # unity variance with shape 2 and scale √2
@test std(d) ≈ 1
@test skewness(d) ≈ 0
@test kurtosis(d) ≈ 0 atol = 10e-10
@test entropy(d) ≈ 1.418938533204673

# Non-Gaussian case
@test cdf(PGeneralizedGaussian(1., 2.0, 2.5), 1.) ≈ 0.5

α, β = rand(2)
d = PGeneralizedGaussian(0.0, α, β)
@test var(d) ≈ α^2 * (gamma(3.0 * inv(β)) / gamma(inv(β)))
@test kurtosis(d) ≈ gamma(5.0 * inv(β)) * gamma(inv(β)) / (gamma(3.0 * inv(β))^2) - 3.0
@test entropy(d) ≈ inv(β) - log( β / (2.0 * α * gamma(inv(β))))


text:
export BEMQuadAdaptive,
    BEMQuadQBF_Adaptive,
    BEMQuadQBF_Graded

using CompactTranslatesDict:
    PeriodicInterval

using BasisFunctions:
    unsafe_eval_element,
    unsafe_weightfun

using DomainIntegrals:
    QuadratureStrategy,
    QuadAdaptive,
    Q_quadgk,
    Q_hcubature,
    UnitIntervalRule,
    graded_rule_left,
    graded_rule_right,
    PointSingularity

###########################################
# Definition of BEM quadrature strategies
###########################################

"""
A `BEMQuadratureStrategy` groups quadrature strategies (as used by `DomainIntegrals.jl`)
for the regular, singular and nearly singular integrals corresponding to BEM matrix entries.

A concrete type should implement `regular`, `singular` and `nearlysingular` to
return the corresponding quadrature strategies.
"""
abstract type BEMQuadratureStrategy <: DomainIntegrals.QuadratureStrategy end

# For DomainIntegral strategies, there is no difference between the three types
regular(quad::QuadratureStrategy) = quad
singular(quad::QuadratureStrategy) = quad
nearlysingular(quad::QuadratureStrategy) = quad

"Generic adaptive quadrature for all BEM matrix entries."
struct BEMQuadAdaptive{T} <: BEMQuadratureStrategy
    adaptive    ::  QuadAdaptive{T}
end

regular(quad::BEMQuadAdaptive) = quad.adaptive
singular(quad::BEMQuadAdaptive) = quad.adaptive
nearlysingular(quad::BEMQuadAdaptive) = quad.adaptive

# Accept the same arguments as the QuadAdaptive constructor of DomainIntegrals.jl
BEMQuadAdaptive(args...) = BEMQuadAdaptive(QuadAdaptive(args...))
BEMQuadAdaptive{T}(atol::Number, args...) where {T} = BEMQuadAdaptive{T}(QuadAdaptive{T}(atol, args...))

"Combination of QBF and adaptive quadrature."
struct BEMQuadQBF_Adaptive{T} <: BEMQuadratureStrategy
    regular     ::  QuadQBF{T}
    adaptive    ::  QuadAdaptive{T}
end

BEMQuadQBF_Adaptive(regular::QuadQBF{T}, args...) where {T} =
    BEMQuadQBF_Adaptive(regular, QuadAdaptive{T}(args...))

regular(quad::BEMQuadQBF_Adaptive) = quad.regular
singular(quad::BEMQuadQBF_Adaptive) = quad.adaptive
nearlysingular(quad::BEMQuadQBF_Adaptive) = quad.adaptive

splineorder(q::BEMQuadQBF_Adaptive) = splineorder(q.regular)

"Combination of QBF and graded quadrature."
struct BEMQuadQBF_Graded{T} <: BEMQuadratureStrategy
    regular         ::  QuadQBF{T}
    graded_left     ::  UnitIntervalRule{T}
    graded_right    ::  UnitIntervalRule{T}
end

BEMQuadQBF_Graded(regular::QuadQBF{T}, σ = T(3)/10, n = 1, μ = one(T)/2, args...) where {T} =
    BEMQuadQBF_Graded(regular, graded_rule_left(σ, n, μ, args...), graded_rule_right(σ, n, μ, args...))

regular(quad::BEMQuadQBF_Graded) = quad.regular
singular(quad::BEMQuadQBF_Graded) = quad
nearlysingular(quad::BEMQuadQBF_Graded) = quad

splineorder(q::BEMQuadQBF_Graded) = splineorder(q.regular)

"A generic combination of BEM quadrature strategies."
struct BEMQuad <: BEMQuadratureStrategy
    regular
    singular
    nearlysingular
end

regular(quad::BEMQuad) = quad.regular
singular(quad::BEMQuad) = quad.singular
nearlysingular(quad::BEMQuad) = quad.nearlysingular

########################################################################
# Definition of projectionintegral and doubleprojectionintegral
########################################################################

function projectionintegral(qs, f, dict, idx, measure, sing = NoSingularity(), domain = support(dict, idx))
    integrand = t -> f(t)*unsafe_eval_element(dict, idx, t)*unsafe_weightfun(measure, t)
    integral(qs, integrand, domain, sing)
end

function projectionintegral(qs::QuadQBF, f, dict, idx, measure, sing, domain)
    # For QBF we don't include the basis function
    integrand = t -> f(t)*unsafe_weightfun(measure, t)
    integral(qs, integrand, domain, sing)*sqrt(length(dict))
end

# Teach DomainIntegrals how to evaluate on a PeriodicInterval
DomainIntegrals.domain_splits(domain::PeriodicInterval) = true
DomainIntegrals.domain_split(domain::PeriodicInterval) = components(domain)

periodizer(A,B) = t -> A + mod(t-A,B-A)

# Special treatment for QBF quadrature: don't split the domain, instead periodize the integrand
function DomainIntegrals.integrate_domain(qs::QuadQBF, integrand, domain::PeriodicInterval, measure, properties...)
    if ncomponents(domain) == 1
        DomainIntegrals.integrate_domain(qs, integrand, component(domain, 1), measure, properties...)
    else
        # only periodize when crossing the boundary
        A, B = extrema(domain.periodicdomain)
        DomainIntegrals.integrate_domain(qs, integrand ∘ periodizer(A,B), domain.subdomain, measure, properties...)
    end
end

# Teach DomainIntegrals how to split a PeriodicInterval
function DomainIntegrals.splitdomain_point(x, domain::PeriodicInterval)
    if ncomponents(domain) == 1
        DomainIntegrals.splitdomain_point(x, component(domain, 1))
    else
        vcat(DomainIntegrals.splitdomain_point(x, component(domain, 1))..., DomainIntegrals.splitdomain_point(x, component(domain, 2))...)
    end
end

"Split a domain (originating from a spline) in subsets on which the spline is smooth."
piecewisedomains(domain::MappedDomain, order) = map_domain.(Ref(inverse_map(domain)), piecewisedomains(superdomain(domain), order))
function piecewisedomains(domain::AbstractInterval, order)
    a, b = extrema(domain)
    h = (b-a)/(order+1)
    [a+i*h..a+(i+1)*h for i in 0:order]
end
function piecewisedomains(domain::PeriodicInterval, order)
    realdomain = domain.subdomain
    periodicdomain = domain.periodicdomain
    a, b = extrema(realdomain)
    h = (b-a)/(order+1)
    segments = [a+i*h..a+(i+1)*h for i in 0:order]
    vcat([collect(components(PeriodicInterval(s, periodicdomain))) for s in segments]...)
end

function projectionintegral(qs::BEMQuadQBF_Graded, f, dict, idx, measure, sing, domain)
    domains = piecewisedomains(support(dict, idx), splineorder(qs))
    sum(projectionintegral(qs, f, dict, idx, measure, sing, d) for d in domains)
end

function leftofsingularity(x, domain, support)
    L = width(support)
    if x - leftendpoint(domain) > L/2
        x = x - L
    elseif rightendpoint(domain) - x > L/2
        x = x + L
    end
    leftendpoint(domain) < x
end

function projectionintegral(qs::BEMQuadQBF_Graded, f, dict, idx, measure, sing::PointSingularity, domain::AbstractInterval)
    x = DomainIntegrals.point(sing)
    domains = DomainIntegrals.splitdomain_point(x, domain)
    z = zero(typeof(f(x)))
    for d in domains
        if leftofsingularity(x, d, support(dict))
            z += projectionintegral(qs.graded_right, f, dict, idx, measure, sing, d)
        else
            z += projectionintegral(qs.graded_left, f, dict, idx, measure, sing, d)
        end
    end
    z
end

doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing = NoSingularity()) =
    doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing, support(dict1, idx1), support(dict2, idx2))

function doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
            domain1::AbstractInterval, domain2::AbstractInterval)
    integrand = t -> f(t[2],t[1])*unsafe_eval_element(dict1, idx1, t[1])*conj(unsafe_eval_element(dict2,idx2, t[2]))*
        unsafe_weightfun(measure1, t[1])*unsafe_weightfun(measure2, t[2])
    integral(qs, integrand, domain1 × domain2, sing)
end

function doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
        domain1::PeriodicInterval, domain2::PeriodicInterval)
    sum(doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing, d1, d2)
            for d1 in components(domain1), d2 in components(domain2))
end

function doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
        domain1::PeriodicInterval, domain2::AbstractInterval)
    sum(doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing, d1, domain2)
            for d1 in components(domain1))
end

function doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
        domain1::AbstractInterval, domain2::PeriodicInterval)
    sum(doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing, domain1, d2)
            for d2 in components(domain2))
end

doubleprojection(qs::QuadQBF, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
            domain1::AbstractInterval, domain2::AbstractInterval) =
        qbf_quadrature2(f, dict1, idx1, measure1, domain1, dict2, idx2, measure2, domain2,
            leftpoint(qs), rightpoint(qs), points(qs), weights(qs))

function doubleprojection(qs::QuadQBF, f, dict1, idx1, measure1, dict2, idx2, measure2, sing,
        domain1::PeriodicInterval, domain2::PeriodicInterval)

    A1, B1 = extrema(domain1.periodicdomain)
    A2, B2 = extrema(domain2.periodicdomain)
    f2 = t -> f(A1+mod(t[1]-A1,B1-A1), A2+mod(t[2]-A2,B2-A2))
    doubleprojection(qs, f, dict1, idx1, measure1, dict2, idx2, measure2, sing, domain1.subdomain, domain2.subdomain)
end


text:
export Webhook

"""
A Webhook.
More details [here](https://discordapp.com/developers/docs/resources/webhook#webhook-object).
"""
struct Webhook
    id::Snowflake
    guild_id::Optional{Snowflake}
    channel_id::Snowflake
    user::Optional{User}
    name::Nullable{String}
    avatar::Nullable{String}
    token::Optional{String}  # Missing in audit log entries.
end
@boilerplate Webhook :constructors :docs :lower :merge :mock


text:
function addToDataChain!(env::Vector{Float64}, atm::Vector{Int32}, bas::Vector{Int32},
                         bf::FloatingGTBasisFuncs)
    center = [bf.center[1](), bf.center[2](), bf.center[3]()]
    xpns = Float64[]
    cons = Float64[]
    for i in bf.gauss
        push!(xpns, i.xpn())
        push!(cons, i.con())
    end
    nGauss = bf.gauss |> length
    envEndIndex = length(env)
    gAtmIndex = length(atm) / 6 |> Int32

    append!(env, center)
    append!(atm, Int32[0, envEndIndex, 1, 0, 0, 0])
    envEndIndex += 3

    append!(env, xpns)
    norm = bf.normalizeGTO ? normOfGTOin(bf) : 1.0
    append!(env, cons.*norm)
    append!(bas, Int32[gAtmIndex, AngularMomentumList[bf.subshell], nGauss, 1, 0,
                       envEndIndex, envEndIndex+nGauss, 0])
    envEndIndex += nGauss*2
    (env, atm, bas)
end

function addToDataChain!(env::Vector{Float64}, atm::Vector{Int32},
                         nuclei::Vector{String}, nucleiCoords::Vector{<:AbstractArray})
    @compareLength nuclei nucleiCoords "nuclei" "their coordinates"
    envEndIndex = length(env)
    len = length(nuclei)
    atmsConfig = Int32[]
    for i = 1:len
        append!(env, Float64[nucleiCoords[i][1], nucleiCoords[i][2], nucleiCoords[i][3]])
        append!(atmsConfig, Int32[getCharge(nuclei[i]), envEndIndex, 1, 0, 0, 0])
        envEndIndex += 3
    end
    append!(atm, atmsConfig)
    (env, atm)
end


text:
using Strategems
using Documenter

Documenter.makedocs(
    modules = [Strategems],
    sitename = "Strategems.jl",
)


text:
"""
Tools for generating demand distributions
"""

"""
Requests for each product is generated with Poisson process (fast coin tossing) with intensity λ.
We consider the time interval of each selling period to be (0,1), so that for each product, expected number of requests is λ.
Because we discretize the time interval (0,1) into 'm' steps, we are approximate the Poisson distribution with Bernouli distribution.
As such, we have λ~mp where p is the probability of success (success is request arrival in one timestep).

Notable property of Poisson process: mix of poisson processes of different products is a poisson process with multiple possible values.
"""

"""
Create the same demand for all products, sum of which across all products will be λ.

    λ = create_λ(20., products)
"""
function create_uniform_λ_per_product(demand::Float64, products::Array{Product{n_edges}}) where n_edges
    λ = fill(demand/(length(products)-1), length(products))
    λ[1]=0. # product[1] is the empty product
    return λ
end

"""
Create same product demand for each product size. For all products of size len, the expected sum of demands (λ) is demand[len].

    λ = create_λ(Float64[10,3,3,5,4], products)
"""
function create_uniform_λ_per_resource(demand::Array{Float64}, products::Array{Product{n_edges}}) where n_edges
    λ = zeros(Float64, length(products))
    sizes = map(sum, products)
    for len in unique(sizes)
        if len==0
            λ[1]=0. # product[1] is the empty product
            continue
        end
        selector = sizes.==len
        λ[selector].= demand[len]/sum(selector)
    end
    return λ
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/rc4.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_rc4_alloc,
    av_rc4_init,
    av_rc4_crypt

function av_rc4_alloc()
    ccall((:av_rc4_alloc,libavutil),Ptr{AVRC4},())
end

function av_rc4_init(d,key,key_bits::Integer,decrypt::Integer)
    ccall((:av_rc4_init,libavutil),Cint,(Ptr{AVRC4},Ptr{UInt8},Cint,Cint),d,key,key_bits,decrypt)
end

function av_rc4_crypt(d,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_rc4_crypt,libavutil),Void,(Ptr{AVRC4},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{UInt8},Cint),d,dst,src,count,iv,decrypt)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
@testset "423.reconstruct-original-digits-from-english.jl" begin
    @test original_digits("owoztneoer") == "012"
    @test original_digits("fviefuro") == "45"
end


text:
__precompile__()

module DualNumbers
  importall Base

  import NaNMath
  import Calculus
  using Compat

  include("dual.jl")
  include("dual_n.jl")

  Base.@deprecate_binding du ɛ
  @deprecate inf{T}(::Type{Dual{T}}) convert(Dual{T}, Inf)
  @deprecate nan{T}(::Type{Dual{T}}) convert(Dual{T}, NaN)

  export
    Dual,
    Dual128,
    Dual64,
    Dual32,
    DualComplex256,
    DualComplex128,
    DualComplex64,
    dual,
    epsilon,
    realpart,
    dualpart,
    isdual,
    dual_show,
    conjdual,
    absdual,
    abs2dual,
    ɛ,
    imɛ
end


text:
using Test

import LinearAlgebra.diagm
import LinearAlgebra.triu
import LinearAlgebra.I
import LinearAlgebra.det

include("../../../lib/KTBC.jl")
import .KTBC.CreateKTBC

function createU(dim=3)
    return I(dim) - diagm(1 => ones(dim - 1))
end

function createS(dim=3)
    return triu(ones(dim, dim))
end

# 3
U₅ = createU(5)
S₅ = createS(5)

U₅ * S₅

@test S₅ == inv(U₅)

# 4
S₄ = triu(ones(4, 4))

S₄ * S₄'

# 5
KTBC₂ = CreateKTBC()
K₂ = convert(Array{Int}, KTBC₂[1])

det(K₂)
inv(K₂)

@test inv(K₂) * det(K₂) == [
  2.0 1.0
  1.0 2.0
]

KTBC₃ = CreateKTBC(3)
K₃ = convert(Array{Int}, KTBC₃[1])

det(K₃)
inv(K₃)

@test round.(Int, inv(K₃) * det(K₃)) == [
  3 2 1
  2 4 2
  1 2 3
]

# 6
function klcounter(n=4)
    j = 1
    result = zeros(n, n)
    while j <= n
        i = j
        while i <= n
            result[i, j] = (n - i + 1) * j
            i += 1
        end
        j += 1
    end
    return result
end

@test klcounter() == [
  4.0  0.0  0.0  0.0
  3.0  6.0  0.0  0.0
  2.0  4.0  6.0  0.0
  1.0  2.0  3.0  4.0
]

# 7
KTBC₃ = CreateKTBC(3)

T₃ = convert(Array{Int}, KTBC₃[2])
K₃ - T₃

u = [
  1
  0
  0
]
vᵀ = transpose(u)

@test u * vᵀ == K₃ - T₃

invT₃ = inv(T₃)
invK₃ = inv(K₃)

@test round.(Int, 4(invT₃ - invK₃)) == [
 9  6  3
 6  4  2
 3  2  1
]

u = [
  3
  2
  1
]
vᵀ = transpose(u)

@test 1 / 4 * u * vᵀ == round.(invT₃ - invK₃; digits=2)


text:
using CacheVariables, BSON, Test

## Add data directory, define data file path
dirpath = joinpath(@__DIR__, "data")
isdir(dirpath) && error("Test directory already has a data subdirectory.")
path = joinpath(dirpath, "test.bson")

## Test save functionality
@testset "Save" begin
    # 1. verify log message
    @test_logs (:info, "Saving to $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end)

    rm(path)
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end

    # 3. did variables enter workspace correctly?
    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Test load functionality
@testset "Load" begin
    # 1. set all variables to nothing
    x = nothing
    y = nothing
    z = nothing
    out = nothing

    # 2. file exists: load variables from it
    # verify log message
    @test_logs (:info, "Loading from $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end)

    # load variables
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end

    # 3. did variables enter workspace correctly?
    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Test overwrite behavior
@testset "Overwrite" begin
    # 1. change file contents
    bson(path; x = nothing, y = nothing, z = nothing, ans = nothing)

    # 2. add `true` to @cache call to overwrite
    # validate log message
    @test_logs (:info, "Overwriting $path\nx\ny\nz") (@cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end true)

    # overwrite data file
    overwrite = true
    @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end overwrite

    # 3. check that correct data was overwritten
    x = nothing
    y = nothing
    z = nothing
    out = nothing
    out = @cache path begin
        x = collect(1:3)
        y = 4
        z = "test"
        "final output"
    end

    @test x == [1, 2, 3]
    @test y == 4
    @test z == "test"
    @test out == "final output"
end

## Test no variable case (where we just save ans)
@testset "Only ans" begin
    # 0. Clean up
    rm(path)

    # 1. verify log message
    @test_logs (:info, "Saving to $path\n") (@cache path 2 + 3)

    # 2. cache ans
    rm(path)
    out = @cache path 2 + 3

    # 3. did variables enter workspace correctly?
    @test out == 5

    # 4. set variable to nothing
    out = nothing

    # 5. file exists: load variables from it
    # verify log message
    @test_logs (:info, "Loading from $path\n") (@cache path 2 + 3)

    # load variables
    out = @cache path 2 + 3

    # 6. did variables enter workspace correctly?
    @test out == 5
end

## Test functionality in module
# Motivated by Pluto and based on test case from:
# https://github.com/JuliaIO/BSON.jl/issues/25
module MyModule
using CacheVariables, Test, DataFrames

@testset "In module" begin
    # 0. module test path
    dirpath = joinpath(@__DIR__, "data")
    modpath = joinpath(dirpath, "modtest.bson")

    # 1a. save
    out = @cache modpath begin
        d = DataFrame(a = 1:10, b = 'a':'j')
        "final output"
    end

    # 1b. check: did variables enter workspace correctly?
    @test d == DataFrame(a = 1:10, b = 'a':'j')
    @test out == "final output"

    # 2. set all variables to nothing
    d = nothing
    out = nothing

    # 3a. load
    out = @cache modpath begin
        d = DataFrame(a = 1:10, b = 'a':'j')
        "final output"
    end

    # 3b. check: did variables enter workspace correctly?
    @test d == DataFrame(a = 1:10, b = 'a':'j')
    @test out == "final output"
end

end

## Test function form
@testset "Function form" begin
    funcpath = joinpath(dirpath, "functest.bson")

    # 1a. Save - verify log message
    @test_logs (:info, "Saving to $funcpath\n") cache(funcpath) do
        x = collect(1:3)
        y = 4
        z = "test"
        (; x = x, y = y, z = z)
    end

    # 1b. Save - save values to cache
    rm(funcpath)
    out = cache(funcpath) do
        x = collect(1:3)
        y = 4
        z = "test"
        (; x = x, y = y, z = z)
    end

    # 1c. Save - did variables enter workspace correctly?
    @test out == (; x = [1, 2, 3], y = 4, z = "test")

    # 2. Reset - set all variables to nothing
    out = nothing

    # 3a. Load - verify log message
    @test_logs (:info, "Loading from $funcpath\n") cache(funcpath) do
        x = collect(1:3)
        y = 4
        z = "test"
        (; x = x, y = y, z = z)
    end

    # 3b. Load - load values from cache
    out = cache(funcpath) do
        x = collect(1:3)
        y = 4
        z = "test"
        (; x = x, y = y, z = z)
    end

    # 3c. Load - did variables enter workspace correctly?
    @test out == (; x = [1, 2, 3], y = 4, z = "test")
end

module MyFuncModule
using CacheVariables, Test, DataFrames

@testset "Function form - in module" begin
    # 0. module test path
    dirpath = joinpath(@__DIR__, "data")
    modpath = joinpath(dirpath, "funcmodtest.bson")

    # 1a. save
    out = cache(modpath; mod = @__MODULE__) do
        DataFrame(a = 1:10, b = 'a':'j')
    end

    # 1b. check: did variables enter workspace correctly?
    @test out == DataFrame(a = 1:10, b = 'a':'j')

    # 2. set all variables to nothing
    out = nothing

    # 3a. load
    out = cache(modpath; mod = @__MODULE__) do
        DataFrame(a = 1:10, b = 'a':'j')
    end

    # 3b. check: did variables enter workspace correctly?
    @test out == DataFrame(a = 1:10, b = 'a':'j')
end

end

## Clean up
rm(dirpath; recursive = true)


text:
module MotifPvalue

#### Dependencies ##################
using DataStructures, DoubleFloats
####################################

#### Exported methods and types ####
export score2pvalue, pvalue2score
####################################

#### Load files ####################
include("helpers.jl")
include("score2pval.jl")
include("pval2score.jl")
####################################

end


text:
using CleanDung
using Test


text:
using Turing, MCMCChains, Distributions, StatsPlots, CSV

@model die_model(observations) = begin
    # Set prior distribution for p (probability of rolling 1)
    p ~ Uniform(0,1)#Beta(5,1)

    # Probability of rolling one of 2 - 6
    q = (1-p)/5

    # Define how samples are obtained
    for i=1:length(observations)
        observations[i] ~ Categorical([p, q, q, q, q, q])
    end

    return p
end

# Define the true probabilities
true_p = 0.25
true_q = (1-true_p)/5

## Simulate observed data
using DataFrames
using CSV
data = rand(Categorical([true_p, true_q, true_q, true_q, true_q, true_q]), 100)
df = DataFrame(X = data)
CSV.write("pp_ex_1_data.csv", df)

# Load observed data
data = CSV.read("pp_ex_1_data.csv")[:,1]

# Use HMC sampler to obtain posterior samples
num_posterior_samples = 1000
ϵ = 0.05 # Leapfrog step size
τ = 10   # Number of leapfrog iterations

# Obtain posterior samples
chain = sample(die_model(data), HMC(ϵ, τ, :p), num_posterior_samples, progress=false)

# Extract summary of p parameter and plot histogram
p_summary = chain[:p]
plot(p_summary, seriestype = :histogram, bg =  RGB(247/255, 236/255, 226/255))

savefig("../images/pp_ex_1_hist.svg")

# Extract 95% interval
println(quantile(Array(p_summary), [0.025, 0.975]))

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
"""

    Mesh

Abstract type representing a discretisation mesh for a numerical scheme.
"""
abstract type Mesh{T<:AbstractVector{Float64}} end

"""

    n_intervals(m)

Total number of cells for a Mesh or DiscretisedFluidQueue.
"""
n_intervals(mesh::Mesh) = length(mesh.nodes) - 1

"""

    Δ(m)

The width of each cell for a Mesh or DiscretisedFluidQueue.
"""
Δ(mesh::Mesh) = mesh.nodes[2:end] - mesh.nodes[1:end-1]
"""

    Δ(m,k::Int)

The width of cell k for a Mesh or DiscretisedFluidQueue.
"""
Δ(mesh::Mesh,k) = mesh.nodes[k+1] - mesh.nodes[k]
Δ(mesh::Mesh{StepRangeLen{Float64}},k) = mesh.nodes.step.hi

"""

    n_bases_per_phase(m)

Total number of basis functions used to represent the fluid queue for each phase
    of a Mesh or DiscretisedFluidQueue.
"""
n_bases_per_phase(mesh::Mesh) = n_bases_per_cell(mesh) * n_intervals(mesh)

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
export Abstract_Solver_Conf

@doc raw"""
```julia
abstract type Abstract_Solver_Conf end
```

Abstract solver configuration. These are the solvers to be used to
solve unconstrained nonlinear least squares:

```math
\min\limits_\theta \frac{1}{2}\|r(\theta)\|^2
```

Implementations:
- [`LevenbergMarquardt_Conf`](@ref)
"""
abstract type Abstract_Solver_Conf end


text:
using Test
using NeXLCore

@testset "Liljequist1989 tests" begin
    # Almost Table II
    @test isapprox(σₜ(Liljequist1989, n"Al", 100.0), 14.001 * NeXLCore.a₀^2, atol = 0.001) # 14.0
    @test isapprox(σₜ(Liljequist1989, n"Al", 500.0), 5.879 * NeXLCore.a₀^2, atol = 0.001) # 5.89
    @test isapprox(σₜ(Liljequist1989, n"Al", 1000.0), 3.523 * NeXLCore.a₀^2, atol = 0.001) # 3.53
    @test isapprox(σₜ(Liljequist1989, n"Al", 5000.0), 0.894 * NeXLCore.a₀^2, atol = 0.001) # 0.899
    @test isapprox(σₜ(Liljequist1989, n"Al", 10000.0), 0.4695 * NeXLCore.a₀^2, atol = 0.001) # 0.474
    @test isapprox(σₜ(Liljequist1989, n"Al", 50000.0), 0.1037 * NeXLCore.a₀^2, atol = 0.001) # 0.109
    @test isapprox(
        σₜ(Liljequist1989, n"Al", 100000.0),
        0.0560 * NeXLCore.a₀^2,
        atol = 0.001,
    ) # 0.0615

    @test isapprox(σₜ(Liljequist1989, n"Au", 100.0), 6.775 * NeXLCore.a₀^2, atol = 0.001) # 6.78
    @test isapprox(σₜ(Liljequist1989, n"Au", 500.0), 8.738 * NeXLCore.a₀^2, atol = 0.001) # 8.75
    @test isapprox(σₜ(Liljequist1989, n"Au", 1000.0), 9.303 * NeXLCore.a₀^2, atol = 0.001) # 9.32
    @test isapprox(σₜ(Liljequist1989, n"Au", 5000.0), 4.762 * NeXLCore.a₀^2, atol = 0.001) # 4.79
    @test isapprox(σₜ(Liljequist1989, n"Au", 10000.0), 3.095 * NeXLCore.a₀^2, atol = 0.001) # 3.13
    @test isapprox(σₜ(Liljequist1989, n"Au", 50000.0), 1.041 * NeXLCore.a₀^2, atol = 0.001) # 1.08
    @test isapprox(σₜ(Liljequist1989, n"Au", 100000.0), 0.630 * NeXLCore.a₀^2, atol = 0.001) # 0.692

    @test isapprox(Rₐ(n"Po"), 0.1208 * 1.0e-8, atol = 0.0001 * 1.0e-8)
    @test isapprox(Rₐ(n"H"), 0.5292 * 1.0e-8, atol = 0.0001 * 1.0e-8)
    @test isapprox(Rₐ(n"Fe"), 0.1786 * 1.0e-8, atol = 0.0001 * 1.0e-8)
end


text:
using Test

@testset "gRPCClient" begin
    include("runtests_routeguide.jl")
    include("runtests_errors.jl")
end


text:
mutable struct Graph{V<:AbstractNode, E<:AbstractEdge, IncList, VK} <: AbstractStaticGraph{V, E}
    nodes::Vector{V}
    edges::Vector{E}
    ilist::IncList # incidence list
    indexof::Dict{VK,Int}
end


text:
module Display_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

S = BnBSolver()
S.Verbosity = "Full"

println("test #1:")
@test_nowarn EAGOBranchBound.print_node!(S,1,3.0,[Interval(1.0,2.0)])
println("test #2:")
@test_nowarn EAGOBranchBound.print_results!(S,2.0,[1.5],false,true)
println("test #3:")
@test_nowarn EAGOBranchBound.print_results!(S,2.0,[1.5],false,false)

end


text:
module Utils

using Test

macro defmock(name::Symbol)
    esc(:($name() = $mock($(QuoteNode(__module__)), $(QuoteNode(name)))))
end

function mock(__module__::Module, name::Symbol)
    # The result of test `$__module__.$name` can be controlled by setting
    # Setting ENV["SampleMockedTests.$name"].
    @test lowercase(get(ENV, "$__module__.$name", "true")) == "true"
end

macro getbool(key, default = "true")
    m = @__MODULE__
    esc(:($lowercase($m.@getstr($key, $default)) == "true"))
end

macro getstr(key, default)
    esc(:($_getstr($(QuoteNode(__module__)), $key, $default)))
end

_getstr(__module__::Module, key::AbstractString, default) =
    get(ENV, "$__module__.$key", default)

end  # module


text:
# run CGinJulia unit tests   Jack Snoeyink Jun 2016
using CGinJulia, Test

tests = ["CG_2dTypes","CG_3dTypes", "LA_Dets"]   # the test file names are stored as strings...
for t in tests
 include("$(t).jl")						 # ... so that they can be evaluated in a loop
end

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
using Documenter, TrussMorph

makedocs(;
    modules=[TrussMorph],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/yijiangh/TrussMorph.jl/blob/{commit}{path}#L{line}",
    sitename="TrussMorph.jl",
    authors="Yijiang Huang",
    assets=[],
)

deploydocs(;
    repo="github.com/yijiangh/TrussMorph.jl",
)


text:
import Base: length, copy
using LaTeXStrings

"""
    Instance
Solution instance
## Examples
Consider the following example:
```jldoctest objective_function; setup = :(using Pidoh)
julia> x = Instance(BitArray([true,true,false,true]), OneMax(4));

julia> fitness(x)
3

julia> optimum(x)
4

```
"""
struct Instance
    individual::Any
    problem::AbstractProblem
    generator::Union{Nothing,AbstractIP}
    fitnessvalue::Number
    name::LaTeXString

    function Instance(
        individual,
        problem::AbstractProblem;
        generator::Union{Nothing,AbstractIP} = nothing,
        fitnessvalue::Number = -1,
        name::LaTeXString = L"problem",
    )
        if fitnessvalue == -1
            fitnessvalue = fitness(individual, problem)
        end
        new(individual, problem, generator, fitnessvalue, name)
    end

    function Instance(
        problem::AbstractProblem;
        generator::Union{Nothing,AbstractIP} = nothing,
        individual = nothing,
        fitnessvalue::Number = -1,
        name::LaTeXString = L"problem",
    )
        if !isnothing(individual) && fitnessvalue == -1
            fitnessvalue = fitness(individual, problem)
        end
        new(individual, problem, generator, fitnessvalue, name)
    end
end

fitness(instance::Instance) = instance.fitnessvalue
fitness(instance::Instance, flippositions) =
    fitness(instance.individual, instance.problem, instance.fitnessvalue, flippositions)
optimum(instance::Instance) = optimum(instance.problem)
isoptimum(instance::Instance) = fitness(instance) == optimum(instance)
length(instance::Instance) = length(instance.individual)

function initial(instance::Instance)
    ind = generate(instance.generator)
    if !isnothing(instance.generator)
        return Instance(instance.problem, individual = ind)
    end
end

function copy(instance::Instance; fitnessvalue::Number = fitness(instance))
    Instance(instance.individual, instance.problem, fitnessvalue = fitnessvalue)
end

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
const labels = ["A1", "A2", "E", "B1", "B2", "A"]
const label_map = Dict(key => "$num" for (num, key) ∈ enumerate(labels))
# add the extra character
label_map["-"] = ""
const column_types = [Int64, Int64, Int64, Float64, Int64, Float64, Float64]

"""Function to parse in a CSV formatted line list output from PGopher.
Kwargs specify what to skip and where to read the header in: by default,
the options expect to skip the first line, read the header from the second,
startin parsing from the third, and skips the last line.
"""
function read_pgopher_linelist(filepath; header=2, skipto=3, footerskip=1)
    csv_data = CSV.File(filepath, header=header, skipto=skipto, footerskip=footerskip)
    ν = csv_data[:Position]
    A = csv_data[:A]
    encodings = clean_linelist_encoding(csv_data)
    transitions = StructArray([Transition(ν=ν[i], I=A[i], encoding=encodings[:,i]) for i ∈ eachindex(ν, A)])
    return transitions
end

"""Parse through an energy level list output of PGopher.
"""
function read_pgopher_levels(filepath)
    open(filepath) do file
        lines = readlines(file)[3:end-3]
        return StructArray(parse_level_row.(lines))
    end
end

"""Not intended for public consumption, but encapsulating the work done
on a row actually helps the eventual `StructArray` creation infer
the correct types, which should give downstream performance improvements
(instead of just `Vector{Any}`).
"""
function parse_level_row(row)
    split_line = replace_sym_labels(row) |> split
    # for remaining columns, parse as integers
    remaining = abs(length(column_types) - length(split_line))
    expected_types = vcat(copy(column_types), [Int64 for i in 1:remaining])
    data = [parse(expec_type, value) for (expec_type, value) in zip(expected_types, split_line)]
    Level(e=data[4], g=data[5], encoding=[val for val in data[8:end]])
end

function replace_sym_labels(string_encoding)
    for (key, value) ∈ label_map
        string_encoding = replace(string_encoding, key=>value)
    end
    string_encoding
end

"""Function to clean up the linelist quantum number encoding. The
first action replaces symmetry labels with a numeric mapping, followed
by splitting each row into a vector of strings. The result is combined
to yield a `M x N` matrix for `M` quantum numbers, and `N` entries.

Finally, the matrix values are converted into integers.
"""
function clean_linelist_encoding(csv_data)
    mapreduce(x->replace_sym_labels(x) |> split, hcat, csv_data[:Label]) |> x -> parse.(Int, x)
end


text:
module StochasticBlockModel

using Distributions, Random, JuMP, GLPKMathProgInterface

include("datasets.jl")
include("sbm.jl")
include("estimate.jl")
include("results.jl")
include("exact.jl")
include("heuristic.jl")

export generate, estimate, Dataset

end # module


max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
#=
        ScalarAffineFunction -in- LessThan
        ScalarAffineFunction -in- GreaterThan
        ScalarAffineFunction -in- EqualTo
        ScalarAffineFunction -in- Interval
=#
constrdict(model::LinQuadOptimizer, ::LCI{LE})  = cmap(model).less_than
constrdict(model::LinQuadOptimizer, ::LCI{GE})  = cmap(model).greater_than
constrdict(model::LinQuadOptimizer, ::LCI{EQ})  = cmap(model).equal_to
constrdict(model::LinQuadOptimizer, ::LCI{IV})  = cmap(model).interval

function MOI.addconstraint!(model::LinQuadOptimizer, func::Linear, set::T) where T <: LinSets
    __assert_supported_constraint__(model, Linear, T)
    canonicalized_func = MOIU.canonical(func)
    add_linear_constraint(model, canonicalized_func, set)
    model.last_constraint_reference += 1
    index = LCI{T}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = get_number_linear_constraints(model)
    push!(model.constraint_primal_solution, NaN)
    push!(model.constraint_dual_solution, NaN)
    push!(model.constraint_constant, func.constant)
    return index
end

"""
    add_linear_constraint(model::LinQuadOptimizer, func::Linear, set::S)

Add a constraint of type `func`-in-`set` to `model`.
"""
function add_linear_constraint(model::LinQuadOptimizer, func::Linear, set::S) where S <: Union{LE, GE, EQ}
    add_linear_constraint(model, func, backend_type(model, set), MOIU.getconstant(set))
end

function add_linear_constraint(model::LinQuadOptimizer, func::Linear, set::IV)
    columns = [get_column(model, term.variable_index) for term in func.terms]
    coefficients = [term.coefficient for term in func.terms]
    A = CSRMatrix{Float64}([1], columns, coefficients)
    add_ranged_constraints!(model, A, [set.lower], [set.upper])
end

function add_linear_constraint(model::LinQuadOptimizer, func::Linear, sense, rhs::Float64)
    if abs(func.constant) > eps(Float64)
        Compat.@warn("Constant in scalar function moved into set.")
    end
    columns = [get_column(model, term.variable_index) for term in func.terms]
    coefficients = [term.coefficient for term in func.terms]
    A = CSRMatrix{Float64}([1], columns, coefficients)
    add_linear_constraints!(model, A, [sense], [rhs - func.constant])
end

#=
    Add linear constraints (plural)
=#

function MOI.addconstraints!(model::LinQuadOptimizer, func::Vector{Linear},
                             set::Vector{S}) where S <: LinSets
    __assert_supported_constraint__(model, Linear, S)
    @assert length(func) == length(set)
    canonicalized_functions = MOIU.canonical.(func)
    num_existing_constraints = get_number_linear_constraints(model)
    add_linear_constraints(model, canonicalized_functions, set)
    indices = Vector{LCI{S}}(undef, length(canonicalized_functions))
    for (i, foo) in enumerate(canonicalized_functions)
        model.last_constraint_reference += 1
        index = LCI{S}(model.last_constraint_reference)
        dict = constrdict(model, index)
        dict[index] = num_existing_constraints + i
        push!(model.constraint_primal_solution, NaN)
        push!(model.constraint_dual_solution, NaN)
        push!(model.constraint_constant, foo.constant)
        indices[i] = index
    end
    return indices
end

function add_linear_constraints(model::LinQuadOptimizer, functions::Vector{Linear},
                                sets::Vector{S}) where S <: LinSets
    return add_linear_constraints(model, functions,
                                  backend_type.(Ref(model), sets),
                                  MOIU.getconstant.(sets))
end

"""
    functions_to_CSRMatrix(model::LinQuadOptimizer, functions::Vector{Linear}, num_non_zeros::Int)

Convert a vector of ScalarAffineFunctions into a CSRMatrix with `num_non_zeros`
non-zero coefficients.
"""
function functions_to_CSRMatrix(model::LinQuadOptimizer, functions::Vector{Linear}, num_non_zeros::Int)
    columns = Vector{Int}(undef, num_non_zeros)
    coefficients  = Vector{Float64}(undef, num_non_zeros)
    # Compute the row pointers in the compressed sparse row matrix. The row
    # pointers are defined recursively:
    #  r[1] = 1
    #  r[i] = r[i - 1] + (number of nonzero elements in the (i - 1)th row)
    #
    # To compute this, we first count up the number of nonzero elements in each
    # row (i - 1), storing the result in r[i]. Then we perform a cumsum on r,
    # storing the result back in r.
    num_rows = length(functions)
    row_pointers = fill(0, num_rows)
    row_pointers[1] = 1
    non_zero_index = 0
    for (row, func) in enumerate(functions)
        for term in func.terms
            non_zero_index += 1
            if non_zero_index > num_non_zeros
                error("There were more non-zero entries in the function than " *
                      "indicated: >$(num_non_zeros).")
            end
            columns[non_zero_index] = get_column(model, term.variable_index)
            coefficients[non_zero_index] = term.coefficient
            if row < num_rows
                row_pointers[row + 1] += 1
            end
        end
    end
    cumsum!(row_pointers, row_pointers)
    return CSRMatrix{Float64}(row_pointers, columns, coefficients)
end

function add_linear_constraints(model::LinQuadOptimizer, functions::Vector{Linear}, sets::Vector{IV})
    # loop through once to get number of non-zeros and to move rhs across
    lower_bounds = [set.lower for set in sets]
    upper_bounds = [set.upper for set in sets]
    num_non_zeros = 0
    for (i, func) in enumerate(functions)
        if abs(func.constant) > eps(Float64)
            Compat.@warn("Constant in scalar function moved into set.")
            lower_bounds[i] -= func.constant
            upper_bounds[i] -= func.constant
        end
        num_non_zeros += length(func.terms)
    end
    A = functions_to_CSRMatrix(model, functions, num_non_zeros)
    add_ranged_constraints!(model, A, lower_bounds, upper_bounds)
end

function add_linear_constraints(model::LinQuadOptimizer, functions::Vector{Linear},
                                senses::Vector, right_hand_sides::Vector{Float64})
    # loop through once to get number of non-zeros and to move rhs across
    nnz = 0
    for (i, func) in enumerate(functions)
        if abs(func.constant) > eps(Float64)
            Compat.@warn("Constant in scalar function moved into set.")
            right_hand_sides[i] -= func.constant
        end
        nnz += length(func.terms)
    end
    A = functions_to_CSRMatrix(model, functions, nnz)
    add_linear_constraints!(model, A, senses, right_hand_sides)
end

#=
    Constraint set of Linear function
=#
MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{LCI{S}}) where S <: Union{LE, GE, EQ} = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::LCI{S}) where S <: Union{LE, GE, EQ}
    row = model[index]
    rhs = get_rhs(model, row)
    S(rhs + model.constraint_constant[row])
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{LCI{IV}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::LCI{IV})
    row = model[index]
    lowerbound, upperbound = get_range(model, row)
    return IV(lowerbound + model.constraint_constant[row],
              upperbound + model.constraint_constant[row])
end

#=
    Constraint function of Linear function
=#

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{<:LCI{<: LinSets}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::LCI{<: LinSets})
    row = model[index]
    columns, coefficients = get_linear_constraint(model, row)
    terms = map(
        (variable, coefficient)->MOI.ScalarAffineTerm{Float64}(coefficient, variable),
        model.variable_references[columns],
        coefficients
    )
    Linear(terms, -model.constraint_constant[row])
end

#=
    Scalar Coefficient Change of Linear Constraint
=#

function MOI.modify!(model::LinQuadOptimizer, index::LCI{S}, change::MOI.ScalarCoefficientChange{Float64}) where S <: LinSets
    row = model[index]
    column = get_column(model, change.variable)
    change_matrix_coefficient!(model, row, column, change.new_coefficient)
end

#=
    Change RHS of linear constraint without modifying sense
=#
MOI.supports(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{LCI{S}}) where S <: LinSets = true
function MOI.set!(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::LCI{S},
                  new_set::S) where S <: Union{LE, GE, EQ}
    change_rhs_coefficient!(model, model[index], MOIU.getconstant(new_set))
end

function MOI.set!(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::LCI{IV},
                  new_set::IV)
    modify_ranged_constraints!(model, [model[index]], [new_set.lower], [new_set.upper])
end

#=
    Delete a linear constraint
=#

function MOI.delete!(model::LinQuadOptimizer, index::LCI{<: LinSets})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    row = dict[index]
    delete_linear_constraints!(model, row, row)
    deleteat!(model.constraint_primal_solution, row)
    deleteat!(model.constraint_dual_solution, row)
    deleteat!(model.constraint_constant, row)
    # shift all the other references
    shift_references_after_delete_affine!(model, row)
    delete!(dict, index)
end

#=
    Transform scalar constraint
=#

function MOI.transform!(model::LinQuadOptimizer, index::LCI{S1}, new_set::S2) where S1 where S2 <: Union{LE, GE, EQ}
    __assert_supported_constraint__(model, Linear, S2)
    dict = constrdict(model, index)
    row = dict[index]
    change_linear_constraint_sense!(model, [row], [backend_type(model, new_set)])
    model.last_constraint_reference += 1
    index_2 = LCI{S2}(model.last_constraint_reference)
    dict_2 = constrdict(model, index_2)
    dict_2[index_2] = row
    delete!(dict, index)
    delete_constraint_name(model, index)
    return index_2
end


text:
using Distributions
using JSON

a = 3
b = 6
c = 5
d = TriangularDist( a,b,c )

x = [ -5, -2.5, 0, 2.5, 5 ]

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("a", a),
	("b", b),
	("c", c),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/partial.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

max_stars_repo_name: distributions-io/triangular-mgf

max_issues_repo_name: distributions-io/triangular-mgf

max_forks_repo_name: distributions-io/triangular-mgf


text:
function displayResults(resultados)
    numStructures =length(resultados)
    sn = resultados[1]
    sa = resultados[2]
    si = resultados[3]
    st = resultados[4]
    sm = resultados[5]
    sf = resultados[6]

    displayResults(sn, sa, si, st, sm, sf)
end

function displayResults(sn, sa, si, st, sm, sf)

    println("Características Generales de la Cabida Óptima:")
    println("----------------------------------------------")
    println("N° Deptos: ", sum(floor.(sa.numDeptosTipo)), " (Máx. Densidad = ", round(Int, sn.maxNumDeptos), ")")
    println("N° Deptos por Tipo: ", floor.(sa.numDeptosTipo)')
    println("Ocupacion: ", round(sa.ocupacion, digits = 2), " m2 (Máx. Ocupación = ", round(sn.maxOcupacion, digits = 2), " m2)")
    println("Constructibilidad: ", round(sa.constructibilidad, digits = 2), " m2 (Máx. Constructibilidad = ", round(sn.maxConstructibilidad, digits = 2), " m2)")
    println("N° Pisos: ", Int(floor(sa.numPisos)), " (Máx. Pisos = ", Int(floor(sn.maxPisos)), ")")
    println("Altura: ", round(sa.altura, digits = 2), " m2 (Máx. Altura = ", round(sn.maxAltura, digits = 2), " m2)")
    println("Superficie Útil SNT: ", round(sa.superficieUtilSNT, digits = 2), " m2")
    println("Superficie Común SNT: ", round(sa.superficieComunSNT, digits = 2), " m2")
    println("Superficie Edificada SNT: ", round(sa.superficieEdificadaSNT, digits = 2), " m2")
    println("Superficie por Piso: ", round(sa.superficiePorPiso, digits = 2))
    println("N° Estac. Vendibles: ", round(Int, sa.estacionamientosVendibles))
    println("N° Estac. Visita: ", round(Int, sa.estacionamientosVisita))
    println("N° Estac. Discapacitados: ", round(Int, sn.minEstacionamientosDiscapacitados))
    println("N° Estac. Totales: ", round(Int, sa.estacionamientosVendibles) + round(Int, sa.estacionamientosVisita) + round(Int, sn.minEstacionamientosDiscapacitados))
    println("N° Estac. Bicicletas: ", round(Int, sa.numBicicleteros))

    println("")
    println("Características del Terreno:")
    println("----------------------------")
    println("Superficie Terreno: ", round(st.superficieTerreno, digits = 2), " m2")
    println("Superficie Bruta: ", round(st.superficieBruta, digits = 2), " m2")
    println("Costo Terreno antes Corredor: UF ", round(st.costoTerrenoAntesCorr, digits = 2), " (", round(st.costoUnitTerrenoAntesCorr, digits = 2), " UF/m2)")
    println("Costo Corredor: UF ", round(st.CostoCorredor, digits = 2), (" (IVA Inc.)"))
    println("Costo Total Terreno: UF ", round(st.costoTerreno, digits = 2), " (", round(st.costoUnitTerreno, digits = 2), " UF/m2)")

    println("")
    println("Indicadores de Desempeño:")
    println("-------------------------")
    println("Ingresos por Ventas: UF ", round(si.IngresosVentas, digits = 2))
    println("Costo Total: UF ", round(si.CostoTotal, digits = 2))
    println("Margen antes de Impuesto: UF ", round(si.MargenAntesImpuesto, digits = 2))
    println("Rentabilidad Total Bruta: UF+ ", round(si.RentabilidadTotalBruta, digits = 2) * 100," %")
    println("TIR antes Impuestos: UF+ ", round(si.TirAntesImpuestos, digits = 2) * 100," %")
    println("Impuesto Renta: UF ", round(si.ImpuestoRenta, digits = 2))
    println("Utilidad después Impuesto: UF ", round(si.UtilidadDespuesImpuesto, digits = 2))
    println("Rentabilidad Total Neta: UF+ ", round(si.RentabilidadTotalNeta, digits = 2) * 100," %")
    println("Incidencia Terreno: ", round(si.IncidenciaTerreno, digits = 2))

end


text:
struct AffineParametricArray{T, N, MT<:AbstractArray{T, N}, VAR} <: AbstractArray{T, N}
    coeffs::Vector{MT}
    vars::Vector{VAR}
end

const AffineParametricMatrix{T, MT} = AffineParametricArray{T, 2, MT} where {T, MT<:AbstractMatrix{T}}

const AffineParametricVector{T, VT} = AffineParametricArray{T, 1, VT} where {T, VT <: AbstractVector{T}}

# apas are callable
function (apa::AffineParametricArray)(p::AbstractVector)
    length(p) + 1 == length(apa.coeffs) || throw(ArgumentError("dimension mismatch"))
    return sum(apa.coeffs[i] * p[i] for i in eachindex(p)) + apa.coeffs[end]
end

function (apa::AffineParametricArray)()
   return apa.coeffs[end] + sum(apa.coeffs[i] * apa.vars[i] for i in eachindex(apa.vars))
end

# Array interface
IndexStyle(::Type{<:AffineParametricArray}) = IndexLinear()
size(apa::AffineParametricArray) = size(apa.coeffs[1])

getindex(apa::AffineParametricArray, idx...) = getindex(apa(), idx...)

# operations
function ==(apa1::AffineParametricArray, apa2::AffineParametricArray)
    return apa1.vars == apa2.vars && apa1.coeffs == apa2.coeffs
end

# unary operations
+(apa::AffineParametricArray) = apa
-(apa::AffineParametricArray) = AffineParametricArray(-apa.coeffs, apa.vars)

# addition subtraction
for op in (:+, :-)
    @eval function $op(apa1::AffineParametricArray, apa2::AffineParametricArray)
        return AffineParametricArray($op(apa1.coeffs, apa2.coeffs), apa1.vars)
    end

    @eval function $op(apa::AffineParametricArray{T, N, MT, VAR}, B::MS) where {T, N, MT, VAR, S, MS<:AbstractArray{S, N}}
        MTS = promote_type(MT, MS)
        coeffs = similar(apa.coeffs, MTS)
        coeffs .= apa.coeffs
        coeffs[end] = $op(coeffs[end], B)
        return AffineParametricArray(coeffs, apa.vars)
    end

    @eval function $op(B::MS, apa::AffineParametricArray{T, N, MT, VAR}) where {T, N, MT, VAR, S, MS<:AbstractArray{S, N}}
        MTS = promote_type(MT, MS)
        coeffs = similar(apa.coeffs, MTS)
        coeffs .= $op(apa.coeffs)
        coeffs[end] = $op(B, apa.coeffs[end])
        return AffineParametricArray(coeffs, apa.vars)
    end
end

# multiplication, backslash
for op in (:*, :\)
    @eval function $op(A::AbstractMatrix, apa::AffineParametricArray)
        coeffs = [$op(A, coeff) for coeff in apa.coeffs]
        return AffineParametricArray(coeffs, apa.vars)
    end
end

function *(apa::AffineParametricMatrix, B::AbstractMatrix)
    coeffs = [coeff*B for coeff in apa.coeffs]
    return AffineParametricArray(coeffs, apa.vars)
end

function *(apa::AffineParametricMatrix, v::AbstractVector)
    coeffs = [coeff*v for coeff in apa.coeffs]
    return AffineParametricArray(coeffs, apa.vars)
end

*(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs * n, apa.vars)
*(n::Number, apa::AffineParametricArray) = AffineParametricArray(apa.coeffs * n, apa.vars)
/(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs / n, apa.vars)

# with AffineExpression
function AffineParametricArray(A::AbstractArray{<:AffineExpression}, vars::Vector{<:AffineExpression})
    ncoeffs = length(vars) + 1
    coeffs = [map(a -> a.coeffs[i], A) for i in 1:ncoeffs]
    return AffineParametricArray(coeffs, vars)
end

function AffineParametricArray(A::AbstractArray{<:Number}, vars::Vector{<:AffineExpression})
    ncoeffs = length(vars) + 1
    coeffs = [zero(A) for _ in 1:ncoeffs]
    coeffs[end] = A
    # vars = [AffineExpression(Vector(c)) for c in eachcol(Matrix{T}(I, ncoeffs, ncoeffs-1))]
    return AffineParametricArray(coeffs, vars)
end

function setindex!(apa::AffineParametricArray, ae::AffineExpression, idx...)
    @inbounds for i in eachindex(ae.coeffs)
        setindex!(apa.coeffs[i], ae.coeffs[i], idx...)
    end
end

max_stars_repo_name: lucaferranti/IntervalLinearAlgebra

max_issues_repo_name: lucaferranti/IntervalLinearAlgebra

max_forks_repo_name: lucaferranti/IntervalLinearAlgebra