text:
"""
    noisy_distribution(;input::Input, reflectivity::Real, interf::Interferometer, exact=true, approx=true, samp=true)

Compute the exact and/or approximated and/or sampled probability distribution of
all possible output configurations of partially-distinguishable photons through a
lossy interferometer.
By default, `exact`, `approx` and `samp` are set to `true` meaning that `noisy_distribution`
returns an array containing the three distributions.

!!! note
    - The probabilities within a distribution are indexed following the same order as [`output_mode_occupation(n,m)`](@ref)
    - The approximated distribution has error and failure probability of ``1e^{-4}``.
!!! note "Reference"
    [https://arxiv.org/pdf/1809.01953.pdf](https://arxiv.org/pdf/1809.01953.pdf)
"""
function noisy_distribution(;input::Input, reflectivity::Real, interf::Interferometer, exact=true, approx=true, samp=true)

    output = []
    ϵ = 1e-4
    δ = 1e-4

    input_modes = input.r.state
    number_photons = input.n
    number_modes = input.m
    U = interf.U

    if get_parametric_type(input)[1] == OneParameterInterpolation
        distinguishability = input.distinguishability_param
    else
        get_parametric_type(input)[1] == Bosonic ? distinguishability = 1.0 : distinguishability = 0.0
    end

    number_output_photons = trunc(Int, number_photons*reflectivity)
    input_occupancy_modes = fill_arrangement(input_modes)

    if reflectivity == 1 ||distinguishability == 0
        throw(ArgumentError("invalid input parameters"))
    end

    ki = (log(ϵ*δ*(1-reflectivity*distinguishability^2)/2))/log(reflectivity*distinguishability^2)
    ki = 10^(trunc(Int, log(10, ki)))
    ki < 10 ? k = number_photons : k = min(ki, number_photons)
    kmax = k-1

    combs = collect(combinations(input_occupancy_modes, k))
    nlist = output_mode_occupation(k, number_modes)

    function compute_pi(i::Integer, ans)

        res = 0
        for j = 1:length(combs)
            perm = collect(permutations(combs[j]))
            for l = 1:length(perm)

                count = sum(collect(Int(perm[l][ll] == combs[j][ll]) for ll = 1:length(perm[l])))
                korder = k - count

                if ans == true
                    pterm = distinguishability^korder * ryser(U[combs[j], nlist[i]] .* conj(U[perm[l], nlist[i]]))
                    pterm = real(pterm / (binomial(number_photons, k) * factorial(k)))
                    res += pterm
                elseif ans == false && korder <= kmax
                    pterm = distinguishability^korder * ryser(U[combs[j], nlist[i]] .* conj(U[perm[l], nlist[i]]))
                    pterm = real(pterm / (binomial(number_photons, k) * factorial(k)))
                    res += pterm
                end

            end
        end

        return res

    end

    compute_full_distribution(ans) = map(i->compute_pi(i,ans), 1:length(nlist))

    function sampling(n_samples = 1e5)

        b = nlist[1]
        comb = input_occupancy_modes[1:k]
        proba = compute_pi(1, true)
        current_pos = 1
        lsample = []
        bprob = zeros(length(nlist))

        for i = 1:n_samples

            pos_test = rand(1:length(nlist))
            b_test = nlist[pos_test]
            comb_test = combs[rand(1:length(combs))]
            perm = collect(permutations(comb_test))
            prob_test = 0

            for j = 1:length(perm)
                count = sum(collect(Int(k - sum(perm[j]) == comb_test)))
                if count <= kmax
                    prob_test += ryser(U[comb_test, b_test] .* conj(U[perm[j], b_test]))
                end
            end

            prob_test = real(prob_test)

            if prob_test > proba
                b = b_test
                proba = prob_test
                comb = comb_test
                current_pos = pos_test
            elseif prob_test / proba > rand()
                b = b_test
                proba = prob_test
                comb = comb_test
                current_pos = pos_test
            end

            bprob[current_pos] += 1 / n_samples
            push!(lsample) = current_pos

        end

        return bprob

    end

    exact ? push!(output, compute_full_distribution(true)) : nothing
    approx ? push!(output, compute_full_distribution(false)) : nothing
    samp ? push!(output, sampling()) : nothing

    return output

end


text:
datadir = "Panda"
panda_examples_url = "https://github.com/frankaemika/franka_ros/tree/kinetic-devel/franka_description"

meshpaths = ["meshes/collision", "meshes/visuals"]
meshnames = ["finger", "hand", "link0", "link1", "link2", "link3", "link4", "link5", "link6", "link7"]

urdfnames = ["hand.urdf.xacro", "hand.xacro", "panda_arm.urdf.xacro", "panda_arm.xacro", "panda_arm_hand.urdf.xacro"]

ispath(datadir) || mkpath(datadir)

ext: jl


text:
# 1D Heat Diffusion in an Infinite Slab
# LBM
#
# Hongyang Zhou, hyzhou@umich.edu 04/30/2019

# Problem Description
# A slab initially at temperature equal to zero, T=0. For time t>0, the left
# surface of the slab is subjected to a high temperature and equal to unity,
# T=1.0: The slab length is 100 units.
# ∂T/∂t = α ∂^2 T / ∂x^2
# Calculate the temperature distribution in the slab for t = 200, α=0.25.

# Solution
# Let us divide the domain of integration into Δx=1.0 and Δt=1.0.
# We are applying D1Q2 LBM with Dirichlet and Riemann BC.
# Note that in cylindrical and spherical coordinate systems, an axis-symmetric
# radial diffusion is similar to the diffusion in 1D Cartesian coordinates with
# an extra source term.

using PyPlot

function solveHeatDiffusion1D()
# Parameters
Δx = 1.0           # Discretized space interval
Δt = 1.0           # Discretized time interval
ω  = 4/3           # Relaxation coefficient
TotalTime = 200.   # Total simulation time
nT = 200           # Total timesteps
nx = 100           # Number of total space intervals
w  = [0.5,0.5]     # Lattice weights

Tleft  = 1.0       # Left boundary temperature
Tright = 0.0       # Right boundary temperature

heatflux  = 100.   # Constant heat flux, [W/mK]
k_thermal = 20.    # Thermal conductivity, [W/mK]

BCDirichlet,BCRiemann = 1,2 # Constants for different boundary types
BCtype = 1         # Type of boundary conditions

UseSource = false  # Logical for calculating source terms
if UseSource
   qg = 1.0        # Rate of heat generation per volume
   ρ  = 1.0        # Density
   C  = 1.0        # Heat capacity
   source = qg/(ρ*C) # Source term. Note that α=k/(ρC)
end

# Initializations
f   = zeros(2,nx)
feq = zeros(2,nx)
T   = zeros(nx)
t   = 0.0     # Current simulation time

# Advance
while t < TotalTime

   # Calculate equilibrium distribution
   feq[1,:] .= w[1].*T[:]
   feq[2,:] .= w[2].*T[:]

   # Collision
   @. f[1,:] = (1.0-ω)*f[1,:] + ω*feq[1,:]
   @. f[2,:] = (1.0-ω)*f[2,:] + ω*feq[2,:]
   if UseSource
      f[1,:] .+= Δt*w[1]*source
      f[2,:] .+= Δt*w[2]*source
   end

   # Streaming
   f[1,2:nx] .= f[1,1:nx-1]
   f[2,1:nx-1] .= f[2,2:nx]

   # Remaining BCs
   if BCtype == BCDirichlet
      f[1,1] = Tleft - f[2,1]
      f[2,nx]= Tright- f[1,nx]
   elseif BCtype == BCRiemann
      f[1,1] = f[1,2] + f[2,2] - f[2,1] + heatflux*Δx/k_thermal
      f[2,nx]= Tright - f[1,nx]
   end

   # Calculate macroscopic T
   @. T = f[1,:] + f[2,:]

   t += Δt
end

# Visualization
plot(0:nx-1,T,linestyle="-",marker="o",label="LBM")
xlim(0,30)
legend()
title("1D Heat Diffusion in an Infinite Slab")
xlabel("x")
ylabel("T")
grid(true)

end


text:
using Unishox
using Test

@testset "Unishox.jl" begin
    s = "What's up, doc?"
    c = compress(s)
    @test 0 < sizeof(c) < sizeof(s)
    @test decompress(c) == s

    @test compress("") == ""
    @test decompress("") == ""

    s = "؉'s ⎨<g"
    d = decompress(s)
    @test sizeof(d) > sizeof(s)

end


text:
# test whether prox satisfies necessary conditions for optimality

using LinearAlgebra
using SparseArrays
using ProximalOperators
using Test

check_optimality(f::LeastSquares, x, gamma, y) = norm(y + gamma * f.lambda * (f.A' * (f.A * y - f.b)) - x) <= 1e-10
check_optimality(f::HuberLoss, x, gamma, y) = isapprox((x - y) / gamma, (norm(y) <= f.rho ? f.mu * y : f.rho * f.mu * y / norm(y)))
check_optimality(f::SqrHingeLoss, x, gamma, y) = isapprox((x - y) / gamma, -2 .* f.mu .* f.y .* max.(0, 1 .- f.y .* y))
check_optimality(::IndSimplex, x, _, y) = begin
    w = x - y
    tau = dot(w, y) / sum(y)
    all(w .<= tau + 10 * eps(eltype(x)))
end
check_optimality(f::IndBallL1, x, gamma, y) = begin
    if norm(x, 1) <= f.r
        return all(y .== x)
    end
    sign_is_correct = (sign.(y) .== 0) .| (sign.(x) .== sign.(y))
    return all(sign_is_correct) && check_optimality(IndSimplex(f.r), abs.(x), gamma, abs.(y))
end

check_optimality(f::TotalVariation1D, x, gamma, y) = begin
    N = length(x)
    # compute dual solution
    R = real(eltype(x))
    u = zeros(R, N+1)
    u[1] = 0
    for k in 2:N+1
        u[k] = x[k-1]-y[k-1]+u[k-1]
    end

    # check whether all duals in interval
    c1 = all(abs.(u) .<= gamma*f.lambda + 10*eps(R))
    # check whether last equals 0 (first by construction)
    c2 = isapprox(u[end], 0, atol=10*eps(R))
    # check whether equal +- gamma*lambda
    h = sign.(y[1:end-1] - y[2:end])
    c3 = all(isapprox.( u[2:end-1] .* abs.(h) , h *f.lambda*gamma))
    return c1 && c2 && c3
end

test_cases = [
    Dict(
        "f" => LeastSquares(randn(20, 10), randn(20)),
        "x" => randn(10),
        "gamma" => rand(),
    ),

    Dict(
        "f" => LeastSquares(randn(15, 40), randn(15), rand()),
        "x" => randn(40),
        "gamma" => rand(),
    ),

    Dict(
        "f" => LeastSquares(rand(Complex{Float64}, 15, 40), rand(Complex{Float64}, 15), rand()),
        "x" => rand(Complex{Float64}, 40),
        "gamma" => rand(),
    ),

    Dict(
        "f" => LeastSquares(sprandn(100,1000,0.05), randn(100), rand()),
        "x" => randn(1000),
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndSimplex(),
        "x" => randn(10),
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndSimplex(2),
        "x" => [1.5, 0.0, 0.5, 0.0, 0.0],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndSimplex(5.0),
        "x" => [0.5, 0.5, 1.0, 1.0, 2.0],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndSimplex(5.0),
        "x" => [0.5, 0.5, 1.0, 1.0, 2.0] * 3,
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndSimplex(rand()),
        "x" => randn(10),
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndBallL1(),
        "x" => [-0.39, 0.1, -0.2, 0.3],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndBallL1(),
        "x" => [0.1, -0.1, 0.2, -0.3, 0.4, 0.5],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndBallL1(1.7),
        "x" => [0.1, -0.1, 0.2, -0.3, 0.4, 0.5],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndBallL1(1.7),
        "x" => [0.4, 0.1, 0.6, 0.2, -0.1, -0.3, 0.5],
        "gamma" => rand(),
    ),

    Dict(
        "f" => IndBallL1(1.7),
        "x" => [0.5, 0.1, -0.1, 0.4, 0.6, -0.3, 0.2] * 10,
        "gamma" => rand(),
    ),

    Dict(
        "f" => HuberLoss(),
        "x" => randn(10),
        "gamma" => rand(),
    ),

    Dict(
        "f" => HuberLoss(rand()),
        "x" => randn(8, 10),
        "gamma" => rand(),
    ),

    Dict(
        "f" => HuberLoss(rand(), rand()),
        "x" => randn(20),
        "gamma" => rand(),
    ),

    Dict(
        "f" => HuberLoss(rand(), rand()),
        "x" => rand(Complex{Float64}, 12, 15),
        "gamma" => rand(),
    ),

    Dict(
        "f" => SqrHingeLoss(randn(5)),
        "x" => randn(5),
        "gamma" => 0.1 + rand(),
    ),

    Dict(
        "f" => SqrHingeLoss(randn(5), 0.1+rand()),
        "x" => randn(5),
        "gamma" => 0.1 + rand(),
    ),

    Dict(
        "f" => SqrHingeLoss(randn(3, 5), 0.1+rand()),
        "x" => randn(3, 5),
        "gamma" => 0.1 + rand(),
    ),

    Dict(
        "f" => TotalVariation1D(0.01),
        "x" => vcat(LinRange(1., -1., 10), -1*ones(3), LinRange(-1., 1., 10)),
        "gamma" => 1.,
    ),

    Dict(
        "f" => TotalVariation1D(1.0),
        "x" => [-2.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 2.4375],
        "gamma" => 1.0,
    ),
    Dict(
        "f" => TotalVariation1D(1.0),
        "x" => [0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 2.4375],
        "gamma" => 1.0,
    ),
]

@testset "Optimality conditions" begin
    @testset "$(typeof(d["f"]))" for d in test_cases
        f, x, gamma = d["f"], d["x"], d["gamma"]
        y, fy = prox(f, x, gamma)
        @test fy ≈ f(y)
        @test check_optimality(f, x, gamma, y)
    end
end


text:
export
    ssh_dir,
    ssh_key_pass,
    ssh_key_name,
    ssh_key_path,
    ssh_pub_key_path,
    ssh_known_hosts_files,
    ssh_known_hosts_file

"""
    ssh_dir() :: String

The `ssh_dir()` function returns the location of the directory where the `ssh`
program keeps/looks for configuration files. By default this is `~/.ssh` but
this can be overridden by setting the envirnoment variable `SSH_DIR`.
"""
ssh_dir() = get(ENV, "SSH_DIR", joinpath(homedir(), ".ssh"))

"""
    ssh_key_pass() :: String

The `ssh_key_pass()` function returns the value of the environment variable
`SSH_KEY_PASS` if it is set or `nothing` if it is not set. In the future, this
may be able to find a password by other means, such as secure system storage, so
packages that need a pawword to decrypt an SSH private key should use this API
instead of directly checking the environment variable so that they gain such
capabilities automatically when they are added.
"""
ssh_key_pass() = get(ENV, "SSH_KEY_PASS", nothing)

"""
    ssh_key_name() :: String

The `ssh_key_name()` function returns the base name of key files that SSH should
use for when estabilishing a connection. There is usually no reason that this
function should be called directly and libraries should generally use the
`ssh_key_path` and `ssh_pub_key_path` functions to get full paths. If the
environment variable `SSH_KEY_NAME` is set then this function returns that;
otherwise it returns `id_rsa` by default.
"""
ssh_key_name() = get(ENV, "SSH_KEY_NAME", "id_rsa")

"""
    ssh_key_path() :: String

The `ssh_key_path()` function returns the path of the SSH private key file that
should be used for SSH connections. If the `SSH_KEY_PATH` environment variable
is set then it will return that value. Otherwise it defaults to returning

    joinpath(ssh_dir(), ssh_key_name())

This default value in turn depends on the `SSH_DIR` and `SSH_KEY_NAME`
environment variables.
"""
function ssh_key_path()
    key_path = get(ENV, "SSH_KEY_PATH", "")
    !isempty(key_path) && return key_path
    return joinpath(ssh_dir(), ssh_key_name())
end

"""
    ssh_pub_key_path() :: String

The `ssh_pub_key_path()` function returns the path of the SSH public key file
that should be used for SSH connections. If the `SSH_PUB_KEY_PATH` environment
variable is set then it will return that value. If that isn't set but
`SSH_KEY_PATH` is set, it will return that path with the `.pub` suffix appended.
If neither is set, it defaults to returning

    joinpath(ssh_dir(), ssh_key_name() * ".pub")

This default value in turn depends on the `SSH_DIR` and `SSH_KEY_NAME`
environment variables.
"""
function ssh_pub_key_path()
    pub_key_path = get(ENV, "SSH_PUB_KEY_PATH", "")
    !isempty(pub_key_path) && return pub_key_path
    key_path = get(ENV, "SSH_KEY_PATH", "")
    !isempty(key_path) && return "$key_path.pub"
    return joinpath(ssh_dir(), ssh_key_name() * ".pub")
end

"""
    ssh_known_hosts_files() :: Vector{String}

The `ssh_known_hosts_files()` function returns a vector of paths of SSH known
hosts files that should be used when establishing the identities of remote
servers for SSH connections. By default this function returns

    [joinpath(ssh_dir(), "known_hosts"), bundled_known_hosts]

where `bundled_known_hosts` is the path of a copy of a known hosts file that is
bundled with this package (containing known hosts keys for `github.com` and
`gitlab.com`). If the environment variable `SSH_KNOWN_HOSTS_FILES` is set,
however, then its value is split into paths on the `:` character (or on `;` on
Windows) and this vector of paths is returned instead. If any component of this
vector is empty, it is expanded to the default known hosts paths.

Packages that use `ssh_known_hosts_files()` should ideally look for matching
entries by comparing the host name and key types, considering the first entry in
any of the files which matches to be the definitive identity of the host. If the
caller cannot compare the key type (e.g. because it has been hashes) then it
must approximate the above algorithm by looking for all matching entries for a
host in each file: if a file has any entries for a host then one of them must
match; the caller should only continue to search further known hosts files if
there are no entries for the host in question in an earlier file.
"""
function ssh_known_hosts_files()
    bundled = bundled_known_hosts()
    default = joinpath(ssh_dir(), "known_hosts")
    value = get(ENV, "SSH_KNOWN_HOSTS_FILES", nothing)
    value === nothing && return [default, bundled]
    isempty(value) && return String[]
    paths = String[]
    for path in split(value, Sys.iswindows() ? ';' : ':')
        if !isempty(path)
            path in paths || push!(paths, path)
        else
            default in paths || push!(paths, default)
            bundled in paths || push!(paths, bundled)
        end
    end
    return paths
end

"""
    ssh_known_hosts_file() :: String

The `ssh_known_hosts_file()` function returns a single path of an SSH known
hosts file that should be used when establishing the identities of remote
servers for SSH connections. It returns the first path returned by
`ssh_known_hosts_files` that actually exists. Callers who can look in more than
one known hosts file should use `ssh_known_hosts_files` instead and look for
host matches in all the files returned as described in that function's docs.
"""
function ssh_known_hosts_file()
    files = ssh_known_hosts_files()
    for file in files
        ispath(file) && return file
    end
    return !isempty(files) ? files[1] :
        isfile("/dev/null") ? "/dev/null" : tempname()
end

## helper functions

const BUNDLED_KNOWN_HOSTS_LOCK = ReentrantLock()
const BUNDLED_KNOWN_HOSTS_FILE = Ref{String}()

function bundled_known_hosts()
    lock(BUNDLED_KNOWN_HOSTS_LOCK) do
        if !isassigned(BUNDLED_KNOWN_HOSTS_FILE)
            BUNDLED_KNOWN_HOSTS_FILE[], io = mktemp()
            write(io, BUNDLED_KNOWN_HOSTS)
            close(io)
        end
    end
    return BUNDLED_KNOWN_HOSTS_FILE[]
end

const BUNDLED_KNOWN_HOSTS = """
github.com ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAq2A7hRGmdnm9tUDbO9IDSwBK6TbQa+PXYPCPy6rbTrTtw7PHkccKrpp0yVhp5HdEIcKr6pLlVDBfOLX9QUsyCOV0wzfjIJNlGEYsdlLJizHhbn2mUjvSAHQqZETYP81eFzLQNnPHt4EVVUh7VfDESU84KezmD5QlWpXLmvU31/yMf+Se8xhHTvKSCZIFImWwoG6mbUoWf9nzpIoaSjB+weqqUUmpaaasXVal72J+UX2B+2RPW3RcT0eOzQgqlJL3RKrTJvdsjE3JEAvGq3lGHSZXy28G3skua2SmVi/w4yCE6gbODqnTWlg7+wC604ydGXA8VJiS5ap43JXiUFFAaQ==
gitlab.com ecdsa-sha2-nistp256 AAAAE2VjZHNhLXNoYTItbmlzdHAyNTYAAAAIbmlzdHAyNTYAAABBBFSMqzJeV9rUzU4kWitGjeR4PWSa29SPqJ1fVkhtj3Hw9xjLVXVYrU9QlYWrOLXBpQ6KWjbjTDTdDkoohFzgbEY=
gitlab.com ssh-ed25519 AAAAC3NzaC1lZDI1NTE5AAAAIAfuCHKVTjquxvt6CM6tdG4SLp1Btn/nOeHHE5UOzRdf
gitlab.com ssh-rsa AAAAB3NzaC1yc2EAAAADAQABAAABAQCsj2bNKTBSpIYDEGk9KxsGh3mySTRgMtXL583qmBpzeQ+jqCMRgBqB98u3z++J1sKlXHWfM9dyhSevkMwSbhoR8XIq/U0tCNyokEi/ueaBMCvbcTHhO7FcwzY92WK4Yt0aGROY5qX2UKSeOvuP4D6TPqKF1onrSzH9bx9XUf2lEdWT/ia1NEKjunUqu1xOB/StKDHMoX4/OKyIzuS0q/T1zOATthvasJFoPrAjkohTyaDUz2LN5JoH839hViyEG82yB+MjcFV5MU3N1l1QL3cVUCh93xSaua1N85qivl+siMkPGbO5xR/En4iEY6K2XPASUEMaieWVNTRCtJ4S8H+9
"""


text:
using FerriteViz
using Test

@testset "FerriteViz.jl" begin
    # Write your tests here.
end


text:
struct Scale
    x::Float64
    y::Float64
end

Base.:+(s1::Scale, s2::Scale) = Scale(s1.x + s2.x, s1.y + s2.y)
Base.:-(s1::Scale, s2::Scale) = Scale(s1.x - s2.x, s1.y - s2.y)
Base.:*(s1::Scale, f::Number) = f * s1
Base.:*(f::Number, s1::Scale) = Scale(f * s1.x, f * s1.y)
Base.:*(s1::Scale, s2::Scale) = Scale(s1.x * s2.x, s1.y * s2.y)
Base.:/(s1::Scale, s2::Scale) = Scale(s1.x / s2.x, s1.y / s2.y)


text:
# Convert
# =======
#
# Conversion methods for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Conversion between sequences of different alphabet size.
for A in [DNAAlphabet, RNAAlphabet]

    # Convert from a 4 bit encoding to a 2 bit encoding.
    @eval function Base.convert(::Type{BioSequence{$(A{2})}}, seq::BioSequence{$(A{4})})
        # TODO: make it faster with bit-parallel algorithm
        newseq = BioSequence{$(A{2})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end

    # Convert from a 2 bit encoding to a 4 bit encoding.
    @eval function Base.convert(::Type{BioSequence{$(A{4})}}, seq::BioSequence{$(A{2})})
        newseq = BioSequence{$(A{4})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end
end

# Conversion between DNA and RNA sequences.
for (A1, A2) in [(DNAAlphabet, RNAAlphabet), (RNAAlphabet, DNAAlphabet)], n in (2, 4)
    # NOTE: assumes that binary representation is identical between DNA and RNA
    @eval function Base.convert(::Type{BioSequence{$(A1{n})}},
                                seq::BioSequence{$(A2{n})})
        newseq = BioSequence{$(A1{n})}(seq.data, seq.part, true)
        seq.shared = true
        return newseq
    end
end

# Convert from a DNA or RNA vector to a BioSequence.
function Base.convert{A<:DNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{DNA})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert{A<:RNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{RNA})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert(::Type{AminoAcidSequence}, seq::AbstractVector{AminoAcid})
    return AminoAcidSequence(seq, 1, endof(seq))
end

# Convert from a BioSequence to to a DNA or RNA vector
Base.convert(::Type{Vector}, seq::BioSequence) = collect(seq)
function Base.convert{A<:DNAAlphabet}(::Type{Vector{DNA}},
                                      seq::BioSequence{A})
    return collect(seq)
end
function Base.convert{A<:RNAAlphabet}(::Type{Vector{RNA}},
                                      seq::BioSequence{A})
    return collect(seq)
end
Base.convert(::Type{Vector{AminoAcid}}, seq::AminoAcidSequence) = collect(seq)

# Covert from a string to a BioSequence and _vice versa_.
function Base.convert{S<:AbstractString}(::Type{S}, seq::BioSequence)
    return convert(S, [Char(x) for x in seq])
end
Base.convert{S<:AbstractString,A}(::Type{BioSequence{A}}, seq::S) = BioSequence{A}(seq)


text:
module Everest

include("loadcode.jl")

end

max_forks_repo_name: ohinder/Everest


text:
import .Utils: UnionFind, unite!, reset!, find

"""
    kruskal_mst(n, edges, weight)

Kruskal's minimum spanning tree algorithm.
Return a vector indicating for each edge whether it is part of the
spanning tree.
"""
function kruskal_mst(n::Int, edges, weight::Vector{Float64})::Vector{Bool}
    @assert length(weight) == length(edges)
    m = length(edges)
    selected, order, uf = Vector{Bool}(undef, m), Vector{Int}(undef, m), UnionFind(n)
    return _kruskal_mst(weight, edges, selected, order, uf)
end

function _kruskal_mst(weight, edges, selected, order, unfi)
    selected .= false
    sortperm!(order, weight, rev=false)  # , alg=Base.Sort.DEFAULT_STABLE)
    reset!(unfi)
    for ei in order
        (u, v) = edges[ei]
        fu, fv = find(unfi, u), find(unfi, v)
        if fu != fv
            unite!(unfi, fu, fv)
            selected[ei] = true
        end
    end
    return selected
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
module SpatioTemporalSystems

using Parameters
const ICERROR = "Incorrect type for `ic`."

include("systems.jl")
include("barkley.jl")
include("henon.jl")
include("kuramoto-sivashinsky.jl")
include("bocf.jl")

end


text:
function use_linear_search(time_series)
    return length(time_series) < 15
end

function univariate_search_bounds(p0)
    return p0[1] - 100.0, p0[1] + 100.0
end

# create diagonal matrix given x0, t, x1 s.t. x1 = e^{At}*x0 (if possible)
function guess_dyn_linear(time_series, x0)
    n = length(time_series[1]) - 1
    t_x1 = time_series[end]
    t = t_x1[1] - time_series[1][1]
    diag = Vector{Float64}(undef, n)
    for i in 1:n
        if iszero(x0[i]) || t_x1[i+1] / x0[i] <= 1e-5
            diag[i] = 0.0  # not solvable
        else
            @assert t > 0.0
            diag[i] = log(t_x1[i+1] / x0[i]) / t
        end
    end
    A = Diagonal{Float64}(diag)
    return A
end

# use constant dynamics
function guess_dyn(time_series, x0)
    n = length(time_series[1]) - 1
    t_x1 = time_series[end]
    t = t_x1[1] - time_series[1][1]
    A = zeros(n, n)
    b = [(t_x1[i+1] - x0[i]) / t for i in 1:n]
    return AffDyn(A, b)
end

function guess_x0(time_series)
    x0 = collect(time_series[1][2:end])
    return x0
end

# convert parameter array to symbolic representation
function parametric_to_symbolic(p, n; first::Bool=false, linear::Bool=false)
    A = Matrix{Float64}(undef, n, n)
    k = 0
    for j in 1:n
        for i in 1:n
            A[i, j] = p[k + i]
        end
        k += n
    end

    if linear
        dyn = A
    else
        b = p[(n*n+1):(n*(n+1))]
        dyn = AffDyn(A, b)
    end

    if first
        x0 = p[(end-n+1):end]
        return dyn, x0
    else
        return dyn
    end
end

# solve ODE problem in continuous time
function solve_ode_continuous(dyn, x0, time)
    f(u, p, t) = dyn.A * u + dyn.b
    problem = DifferentialEquations.ODEProblem(f, x0, time)
    sol = DifferentialEquations.solve(problem)
    return sol
end

# solve ODE problem at discrete points in time
function solve_ode_discrete(dyn, x0, time)
    return (t) -> successor(dyn, t, x0)
end

# distance between ODE solution and time series
function distance(dyn, x0, time_series)
    d = -Inf
    t0 = time_series[1][1]
    for t_p_ts in time_series
        t = t_p_ts[1] - t0
        p_ts = t_p_ts[2:end]
        p_ode = successor(dyn, t, x0)
        if p_ode == nothing
            return Inf  # exact solution failed
        end
        d_raw = map(-, p_ts, p_ode)
        d_new = norm(d_raw, Inf)  # infinity norm of the distance
        d = max(d, d_new)
    end
    return d
end

# distance between parametric representation and time series
function parametric_problem(p, time_series, n, x0; linear::Bool=false)
    first = x0 == nothing

    # convert to symbolic representation
    result = parametric_to_symbolic(p, n; first=first, linear=linear)
    if first
        dyn, x0 = result
    else
        dyn = result
    end

    # compute distance between ODE solution and time series
    d = distance(dyn, x0, time_series)

    return d
end

# convert symbolic representation to parameter array
function symbolic_to_parametric(dyn, n, x0=nothing; linear::Bool=linear)
    first = x0 != nothing
    b_factor = linear ? 0 : 1
    A = linear ? dyn : dyn.A

    p = Vector{Float64}(undef,
                        first ? n * (n + 1 + b_factor) : n * (n + b_factor))
    k = 0
    for j in 1:n
        for i in 1:n
            p[k + i] = A[i, j]
        end
        k += n
    end
    if !linear
        p[(n*n+1):(n*(n+1))] = dyn.b
    end
    if first
        p[(end-n+1):end] = x0
    end
    return p
end

function time_series_to_symbolic_one_piece(time_series, dyn_init, x0_init;
                                           first::Bool=false,
                                           linear::Bool=false)
    # problem dimension
    n = length(time_series[1]) - 1

    # distance function to optimize
    x0_init_use = first ? nothing : x0_init
    dist = (p) -> parametric_problem(p, time_series, n, x0_init_use;
                                     linear=linear)

    # convert initial guess to parametric representation
    if first
        p0 = symbolic_to_parametric(dyn_init, n, x0_init; linear=linear)
    else
        p0 = symbolic_to_parametric(dyn_init, n; linear=linear)
    end

    # optimize
    if length(p0) == 1
        bounds = univariate_search_bounds(p0)
        sol = optimize(dist, bounds[1], bounds[2])
    else
        sol = optimize(dist, p0)
    end
    p = sol.minimizer
    min_d = sol.minimum

    result = parametric_to_symbolic(p, n; first=first, linear=linear)
    if first
        dyn, x0 = result
        return dyn, x0, min_d
    else
        dyn = result
        return dyn, min_d
    end
end

function time_series_to_symbolic_longest_piece(time_series, ε, dyn_init,
                                               x0_init;
                                               linear_search::Bool=true,
                                               first::Bool=false)
    guess_dyn_first = dyn_init == nothing
    solution_found = false
    i = length(time_series)
    binary_search_range = [1, i+1]  # range [l, u] where l works and u does not
    dyn_best = nothing
    x0_best = nothing
    while i > 1
        time_series_shortened = time_series[1:i]

        for linear in [true, false]
            if length(time_series_shortened) == 2 && linear
                # skipping linear fit
                continue
            end

            if guess_dyn_first
                if linear
                    dyn_init_use = guess_dyn_linear(time_series_shortened,
                                                    x0_init)
                else
                    dyn_init_use = guess_dyn(time_series_shortened, x0_init)
                end
            else
                if linear
                    dyn_init_use = dyn_init.A
                else
                    dyn_init_use = dyn_init
                end
            end

            result = time_series_to_symbolic_one_piece(time_series_shortened,
                                                       dyn_init_use, x0_init;
                                                       first=first,
                                                       linear=linear)
            if first
                dyn, x0, min_d = result
            else
                dyn, min_d = result
            end
            if min_d <= ε
                if linear
                    @assert dyn isa AbstractMatrix
                    n = length(time_series[1]) - 1
                    dyn = AffDyn(dyn, zeros(n))
                else
                    @assert dyn isa AffDyn
                end
                dyn_best = dyn
                if first
                    x0_best = x0
                end
                solution_found = true
                break
            end
        end
        guess_dyn_first = true

        if linear_search
            # linear search
            if !solution_found
                i -= 1
                solution_found = false
                continue
            end
        else
            # binary search, rounding up
            if solution_found
                binary_search_range[1] = i
            else
                binary_search_range[2] = i
            end
            if binary_search_range[1] < binary_search_range[2] - 1
                i = div(binary_search_range[1] + binary_search_range[2] + 1, 2)
                solution_found = false
                continue
            end
        end

        # found the largest solution
        if !linear_search
            i = binary_search_range[1]
        end
        @assert dyn_best != nothing
        if first
            @assert x0_best != nothing
            return dyn_best, x0_best, i
        else
            return dyn_best, i
        end
    end
    @assert false "it should always be possible to fit a function to two points"
end

function time_series_to_symbolic(time_series, ε, dyn_init, x0_init;
                                 linear_search::Bool=use_linear_search(time_series),
                                 log::Bool=false)
    @assert iszero(time_series[1][1]) "the time series should start with t = 0 but started with " *
        time_series[1][1]

    if x0_init == nothing
        x0_init = guess_x0(time_series)
    end

    dyn, x0, l = time_series_to_symbolic_longest_piece(time_series, ε, dyn_init,
        x0_init; linear_search=linear_search, first=true)
    t = time_series[l][1]
    x1 = successor(dyn, t, x0)
    dyns = [dyn]
    times = [t]
    indices = [l]
    while l < length(time_series)
        time_series_shortened = time_series[l:end]
        dyn, Δl = time_series_to_symbolic_longest_piece(time_series_shortened,
            ε, nothing, x1; linear_search=linear_search, first=false)
        l_new = l + Δl - 1
        t = time_series[l_new][1] - time_series[l][1]
        x1 = successor(dyn, t, x1)
        l = l_new
        push!(dyns, dyn)
        push!(times, time_series[l_new][1])
        push!(indices, l)
    end
    if log
        println("found a function with $(length(dyns)) pieces")
    end
    return dyns, x0, times, indices
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
# This code presents an example to illustrate how HyperbolicCrossApprox can be used

using HyperbolicCrossApprox

function test_hyperbolic_cross_approx()

  d  = 5  # Set the number of dimensions
  k = 7  # Set the level of approximation

  g1, m1 = hyperbolic_cross_grid(chebyshev_nodes,d,k)
  g2, m2 = hyperbolic_cross_grid(chebyshev_extrema,d,k)  # Construct the hyperbolic cross grid and the multi index

  # An arbitrary test function

  function test(grid)

    y_value = (grid[:,1].+1).^0.1.*exp.(grid[:,2]).*log.(grid[:,3].+2).^0.2.*(grid[:,4].+2).^0.8.*(grid[:,5].+7).^0.1

    return y_value

  end

  y1 = test(g1)  # Evaluate the test function on the hyperbolic cross grid
  y2 = test(g2)

  point = [0.75, 0.45, 0.82, -0.15, -0.95]  # Choose a point to evaluate the approximated function

  # Chebyshev nodes

  # One way of computing the weights and evaluating the approximated function

  w1       = hyperbolic_cross_weights(y1,g1,m1)    # Compute the hyperbolic cross weights
  y1_hat   = hyperbolic_cross_evaluate(w1,point,m1)  # Evaluate the approximated function

  w2       = hyperbolic_cross_weights(y2,g2,m2)    # Compute the hyperbolic cross weights
  y2_hat   = hyperbolic_cross_evaluate(w2,point,m2)  # Evaluate the approximated function

  #= A second way of computing the weights and evaluating the approximated function that
  computes the interpolation matrix just once. =#

  interp_mat = hyperbolic_cross_inverse_interpolation_matrix(g1,m1)  # Compute the interpolation matrix
  w1b = hyperbolic_cross_weights(y1,interp_mat)             # Compute the hyperbolic cross weights
  y1b_hat = hyperbolic_cross_evaluate(w1b,point,m1)  # Evaluate the approximated function

  # Evaluate the exact function at point

  y_actual = test(point')

  # Now consider the ansiotropic case

  k = 7
  n = [15, 11, 11, 11, 15]
  g3, m3 = hyperbolic_cross_grid(chebyshev_nodes,d,k,n)  # Construct the hyperbolic cross grid and the multi index
  y3 = test(g3)
  w3 = hyperbolic_cross_weights(y3,g3,m3)                   # Compute the hyperbolic cross weights
  y3_hat_ansio = hyperbolic_cross_evaluate(w3,point,m3)       # Evaluate the approximated function
  w3_th = hyperbolic_cross_weights_threaded(y3,g3,m3)                   # Compute the hyperbolic cross weights

  g4, m4 =  hyperbolic_cross_grid(chebyshev_extrema,d,k)
  y4 = test(g4)
  w4 = hyperbolic_cross_weights(y4,g4,m4)
  y4_hat_ansio = hyperbolic_cross_evaluate(w4,point,m4)
  w4_th = hyperbolic_cross_weights_threaded(y4,g4,m4)

  return y_actual, y1_hat, y2_hat, y1b_hat, y3_hat_ansio, y4_hat_ansio

end

test_hyperbolic_cross_approx()


text:
module OrderBook

    using Dates
    using Printf

    using DataFrames
    using DataStructures
    using WebSockets

    export OrderBook

    include("structs.jl")

end


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

struct DeploymentsApi <: SwaggerApi
    client::Swagger.Client
end

"""
Cancels a currently running template deployment.
You can cancel a deployment only if the provisioningState is Accepted or Running. After the deployment is canceled, the provisioningState is set to Canceled. Canceling a template deployment stops the currently running template deployment and leaves the resource group partially deployed.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function deploymentsCancel(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsCancel", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsCancel", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsCancel", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsCancel", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}/cancel", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Checks whether the deployment exists.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function deploymentsCheckExistence(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsCheckExistence", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsCheckExistence", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsCheckExistence", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsCheckExistence", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "HEAD", Nothing, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""
Deploys resources to a resource group.
You can provide the template and parameters directly in the request or link to JSON files.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: parameters::Deployment (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DeploymentExtended
"""
function deploymentsCreateOrUpdate(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsCreateOrUpdate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsCreateOrUpdate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsCreateOrUpdate", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsCreateOrUpdate", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "PUT", DeploymentExtended, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""
Deletes a deployment from the deployment history.
A template deployment that is currently running cannot be deleted. Deleting a template deployment removes the associated deployment operations. Deleting a template deployment does not affect the state of the resource group. This is an asynchronous operation that returns a status of 202 until the template deployment is successfully deleted. The Location response header contains the URI that is used to obtain the status of the process. While the process is running, a call to the URI in the Location header returns a status of 202. When the process finishes, the URI in the Location header returns a status of 204 on success. If the asynchronous request failed, the URI in the Location header returns an error-level status code.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function deploymentsDelete(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsDelete", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsDelete", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsDelete", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsDelete", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Exports the template used for specified deployment.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DeploymentExportResult
"""
function deploymentsExportTemplate(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsExportTemplate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsExportTemplate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsExportTemplate", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsExportTemplate", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "POST", DeploymentExportResult, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}/exportTemplate", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a deployment.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DeploymentExtended
"""
function deploymentsGet(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsGet", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsGet", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsGet", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsGet", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", DeploymentExtended, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Get all the deployments for a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Return: DeploymentListResult
"""
function deploymentsListByResourceGroup(_api::DeploymentsApi, resourceGroupName::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsListByResourceGroup", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsListByResourceGroup", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", DeploymentListResult, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Validates whether the specified template is syntactically correct and will be accepted by Azure Resource Manager..
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: parameters::Deployment (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DeploymentValidateResult
"""
function deploymentsValidate(_api::DeploymentsApi, resourceGroupName::String, deploymentName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentsValidate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentsValidate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentsValidate", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentsValidate", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "POST", DeploymentValidateResult, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/Microsoft.Resources/deployments/{deploymentName}/validate", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export deploymentsCancel, deploymentsCheckExistence, deploymentsCreateOrUpdate, deploymentsDelete, deploymentsExportTemplate, deploymentsGet, deploymentsListByResourceGroup, deploymentsValidate


text:
@testset "reduce" begin

# test of reduce with no map stage
sc = SparkContext(master="local")
txt = parallelize(sc, ["hello", "world"])

@test reduce(txt, *) == "helloworld"

close(sc)

end


text:
module LoadFlint

using Libdl

if VERSION > v"1.3.0-rc4"
  # this should do the dlopen for 1.3 and later
  # and imports the libxxx variables
  using GMP_jll
  using MPFR_jll
  using FLINT_jll
else
  deps_dir = joinpath(@__DIR__, "..", "deps")
  include(joinpath(deps_dir,"deps.jl"))
end

libflint_handle = C_NULL

const __isthreaded = Ref(false)

function __init__()
  if VERSION < v"1.3.0-rc4"
    # this does the dlopen for 1.0-1.2
    check_deps()
  end

  l = dllist()
  f = filter(x->occursin("libflint", x), l)
  if length(f) != 1
    error("there should be exactly one libflint, but we have: ", f)
  end

  # the [.-] at the end is important to avoid matching libgmpxx
  f = filter(x->occursin(r"libgmp[.-]", x), dllist())
  if length(f) == 0
    error("there must be at least one libgmp loaded.")
  elseif length(f) > 1
    # at the moment there doesnt seem to be a way to avoid this
    # because julia comes with libgmp
    # and GMP_jll will load another libgmp
    @debug("there should be exactly one libgmp, but we have: ", f)
  end

  # variable libflint comes from deps file or flint_jll
  global libflint_handle = dlopen(libflint,RTLD_NOLOAD)

  if !Sys.iswindows() && !__isthreaded[]
    #to match the global gmp ones
    fm = dlsym(libflint_handle, :__flint_set_memory_functions)
    ccall(fm, Nothing,
      (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}),
        cglobal(:jl_malloc),
        cglobal(:jl_calloc),
        cglobal(:jl_realloc),
        cglobal(:jl_free))
  end
end

end # module


text:
module SpgpGpu

using LinearAlgebra

export predict
include("predict.jl")

export kernel
include("kernel.jl")

end


text:
fanprojection(e::Ellipse, βᵢ, γᵢ, D) = projection(e, βᵢ + γᵢ, D * sin(γᵢ))

function fanprojection(
    ellipses::AbstractVector{<:Ellipse},
    β::AbstractVector,
    γ::AbstractVector,
    D::AbstractFloat,
)
    [sum(fanprojection(e, βᵢ, γᵢ, D) for e in ellipses) for γᵢ in γ, βᵢ in β]
end

function g(N::Integer, α::AbstractFloat)
    # p. 83 eq. 96
    g = zeros(N)
    N2 = div(N, 2)
    for (i, gi) in enumerate(g)
        n = i - N2 - 1
        if mod(n, 2) != 0
            g[i] = 0.5 * (1 / (π * sin(n * α)))^2
        elseif n == 0
            g[i] = 1 / (8 * α^2)
        end
    end
    g
end

function fanWeight!(R::AbstractMatrix, γ::AbstractVector, D::AbstractFloat)
    views = size(R, 2)
    for βᵢ = 1:views
        @inbounds @. R[:, βᵢ] *= D * cos(γ)
    end
end

function fanFilter!(Q::AbstractMatrix, α::AbstractFloat)
    samples, views = size(Q)
    N = nextpow(2, 2 * samples - 1)
    ZP(x) = vcat(x, zeros(N - samples))
    w = ZP(g(samples, α))
    # indices to slice out unpadded projection?
    i = samples ÷ 2 + 1
    j = i + samples - 1
    for βᵢ = 1:views
        # p. 83 eq. 94
        @inbounds Q[:, βᵢ] .= α * real.(ifft(fft(ZP(Q[:, βᵢ])) .* fftshift(fft(w))))[i:j]
    end
end

function fanReconstruction(
    R::AbstractMatrix,
    β::AbstractVector,
    γ::AbstractVector,
    D::AbstractFloat,
)
    Q = copy(R)
    α = step(γ)
    fanWeight!(Q, γ, D)
    fanFilter!(Q, α)

    Δβ = 2π / length(β)
    pixels = 128
    f = zeros(pixels, pixels)

    Q = LinearInterpolation((γ, β), Q, extrapolation_bc = Line())

    for (i, xᵢ) in enumerate(range(-1, 1, length = pixels)),
        (j, yⱼ) in enumerate(range(-1, 1, length = pixels))

        @inbounds f[j, i] = sum(β) do βᵢ
            r = √(xᵢ^2 + yⱼ^2)
            ϕ = atan(yⱼ, xᵢ)
            γ′ = atan((r * cos(βᵢ - ϕ)) / (D + r * sin(βᵢ - ϕ)))
            L² = (xᵢ - D * cos(βᵢ))^2 + (yⱼ - D * sin(βᵢ))^2
            Q(γ′, βᵢ) / L²
        end
    end
    Δβ .* f
end


text:
cep_prime(n::Order,q::Number,z::Number,N::Integer=matsize3(n,q)) = cep_prime(n,q,z:z,N)[1]
function cep_prime(n::Order,q::Number,z::AbstractVector,N::Integer=matsize3(n,q))
    An = cep_coeffs(n,q,N) # Calculate coefficients
    # Evaluate sine for all z values
    C = [-2k*sin(2k*z) for z in z, k in 0:N-1]
    # Calculate function at z, for all n
    return C * An
end

cea_prime(n::Order,q::Number,z::Number,N::Integer=matsize3(n,q)) = cea_prime(n,q,z:z,N)[1]
function cea_prime(n::Order,q::Number,z::AbstractVector,N::Integer=matsize3(n,q))
    An = cea_coeffs(n,q,N) # Calculate coefficients
    # Evaluate sine for all z values
    C = [-(2k+1)*sin((2k+1)*z) for z in z, k in 0:N-1]
    # Calculate function at z, for all n
    return C * An
end

sea_prime(n::Order,q::Number,z::Number,N::Integer=matsize3(n,q)) = sea_prime(n,q,z:z,N)[1]
function sea_prime(n::Order,q::Number,z::AbstractVector,N::Integer=matsize3(n,q))
    Bn = sea_coeffs(n,q,N) # Calculate coefficients
    # Evaluate cosine for all z values
    S = [(2k+1)*cos((2k+1)*z) for z in z, k in 0:N-1]
    # Calculate function at z, for all n
    return S * Bn
end

sep_prime(n::Order,q::Number,z::Number,N::Integer=matsize3(n,q)) = sep_prime(n,q,z:z,N)[1]
function sep_prime(n::Order,q::Number,z::AbstractVector,N::Integer=matsize3(n,q))
    Bn = sep_coeffs(n,q,N) # Calculate coefficients
    # Evaluate cosine for all z values
    S = [(2k+2)*cos((2k+2)*z) for z in z, k in 0:N-1]
    # Calculate function at z, for all n
    return S * Bn
end


text:
# In this file I compare three values for the model integrals

push!(LOAD_PATH,"src")
using DalitzPlotAnalysis
using amplitudes_compass
using PWAHelper

# calculate integrals with the acceptance
@time precalculate_compass_basis(joinpath("data","2300_2320_t1_mc.txt"), "mc.jld")
@time const PsiMC = read_precalc_basis("mc.jld");
@time const BmatMC = [sum(PsiMC[e,i]'*PsiMC[e,j] for e in 1:size(PsiMC,1))
    for i=1:Nwaves, j=1:Nwaves] /size(PsiMC,1);

# calculate integrals on phase space
@time precalculate_compass_basis(joinpath("data","2300_2320_t1_fu.txt"), "fu.jld")
@time const PsiFU = read_precalc_basis("fu.jld");
@time const BmatFU = [sum(PsiFU[e,i]'*PsiFU[e,j] for e in 1:size(PsiFU,1))
    for i=1:Nwaves, j=1:Nwaves] /size(PsiFU,1);

# same integral analytically
using Cuba
d2_analy = let s = 2.31^2
    function diagint(ind)
        function integrand(x,f)
            σ1 = 4mπ2+x[1]*((√s-mπ)^2-4mπ2)
            cosθ1, ϕ1  = 2*x[2]-1, π*(2*x[3]-1)
            cosθ23,ϕ23 = 2*x[4]-1, π*(2*x[5]-1)
            val = COMPASS_wave(ind,s,σ1,cosθ1,ϕ1,cosθ23,ϕ23)
            f[1] = abs2(val)*sqrt(λ(σ1,mπ2,mπ2)*λ(s,σ1,mπ2))/(2π*s*σ1)*((√s-mπ)^2-4mπ2)
        end
        cuhre(integrand, 5, 1)[1][1]
    end
    diagint(2)/diagint(1)
end
d2_mc = diag(BmatMC)[2]
d2_fu = diag(BmatFU)[2]

@show d2_analy, d2_mc, d2_fu


text:
#=
Resolution of singularities
=#

using Makie

s = Scene()
θs = -π / 2:0.05:π / 2
r = 2
lines!(zeros(length(θs)), zeros(length(θs)), θs, linewidth = 10)
for θ in θs
    lines!(s, r * [-cos(θ), cos(θ)], r * [-sin(θ), sin(θ)], [θ, θ])
end

x = 0.001:0.001:1.5
y = x.^(3 / 2)
z = y ./ x
lines!(s, x, y, fill(-π / 2, length(x)), linewidth = 5.)
lines!(s, x, y, z, linewidth = 5., color = :blue)

x = 0.001:0.001:1.5
y = -x.^(3 / 2)
z = y ./ x
lines!(s, x, y, fill(-π / 2, length(x)), linewidth = 5.)
lines!(s, x, y, z, linewidth = 5., color = :blue)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text: inversedict(dict::Dict{T,U}) where {T,U} = Dict{U,T}(v=>k for (k,v) ∈ dict)


text:
using ChainRulesCore: ChainRulesCore

unwrap_maybe(x) = x
unwrap_maybe(x::SLArray) = x.__x
unwrap_maybe(x::LArray) = x.__x
# Respect the thunk.
function unwrap_maybe(x::ChainRulesCore.Thunk)
    return ChainRulesCore.@thunk(unwrap_maybe(ChainRulesCore.unthunk(x)))
end

function ChainRulesCore.rrule(::typeof(getproperty), A::Union{SLArray,LArray}, s::Symbol)
    function getproperty_LArray_adjoint(d)
        # Hopefully this reference to `A` is optimized away.
        Δ = similar(A) .= 0
        setproperty!(Δ, s, ChainRulesCore.unthunk(d))
        return (ChainRulesCore.NoTangent(), Δ, ChainRulesCore.NoTangent())
    end
    return getproperty(A, s), getproperty_LArray_adjoint
end

function ChainRulesCore.rrule(::Type{LArray{S}}, x::AbstractArray) where {S}
    # This rule covers constructors of the form `LArray{(:a, :b)}(x)`
    # which, amongst other places, is also used in the `@LArray` macro.
    LArray_adjoint(Δ) = ChainRulesCore.NoTangent(), unwrap_maybe(Δ)
    return LArray{S}(x), LArray_adjoint
end

function ChainRulesCore.rrule(::Type{SLArray{Size,S}}, x::AbstractArray) where {Size,S}
    # This rule covers constructors of the form `SLArray{(2, ), (:a, :b)}(x)`
    # which, amongst other places, is also used in the `@LArray` macro.
    SLArray_adjoint(Δ) = ChainRulesCore.NoTangent(), unwrap_maybe(Δ)
    return SLArray{Size,S}(x), SLArray_adjoint
end


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

mutable struct IoK8sApiCoreV1PodTemplateSpec <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodSpec } # spec name: spec

    function IoK8sApiCoreV1PodTemplateSpec(;metadata=nothing, spec=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodTemplateSpec, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1PodTemplateSpec, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        o
    end
end # type IoK8sApiCoreV1PodTemplateSpec

const _property_map_IoK8sApiCoreV1PodTemplateSpec = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"))
const _property_types_IoK8sApiCoreV1PodTemplateSpec = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiCoreV1PodSpec")
Base.propertynames(::Type{ IoK8sApiCoreV1PodTemplateSpec }) = collect(keys(_property_map_IoK8sApiCoreV1PodTemplateSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodTemplateSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodTemplateSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodTemplateSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodTemplateSpec[property_name]

function check_required(o::IoK8sApiCoreV1PodTemplateSpec)
    true
end

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


text:
# run with julia --track-allocation=user

using GeometricIntegrators
using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem
using Profile
using Test

set_config(:nls_atol_break, Inf)
set_config(:nls_rtol_break, Inf)
set_config(:nls_stol_break, Inf)

Δt = 0.1
nt = 10000
ns = 100

nsave = 1
nwrte = 10000

h5file = "test.hdf5"

ode = harmonic_oscillator_ode(vcat(rand(1,ns), zeros(1,ns)))
#tab = getTableauImplicitMidpoint()
tab = getTableauERK4()

sim = Simulation(ode, tab, Δt, "Serial Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
run!(sim)
rm(h5file)

sim = Simulation(ode, tab, Δt, "Serial Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)

Profile.clear()
Profile.clear_malloc_data()

@profile run!(sim)
rm(h5file)

sim = ParallelSimulation(ode, tab, Δt, "Parallel Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
run!(sim)
rm(h5file)

sim = ParallelSimulation(ode, tab, Δt, "Parallel Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)

Profile.clear()
Profile.clear_malloc_data()

@profile run!(sim)
rm(h5file)


text:
mutable struct Lexer
    input::String
    next::Union{Tuple{Char,Int},Nothing}
end

Lexer(input::String) = Lexer(input, iterate(input))

function read_char(l::Lexer)
    if isnothing(l.next)
        return
    else
        ch, _ = l.next
        return ch
    end
end

function peek_char(l::Lexer)
    if isnothing(l.next)
        return
    else
        _, state = l.next
        nxt = iterate(l.input, state)
        if !isnothing(nxt)
            ch, _ = nxt
            return ch
        else
            return
        end
    end
end

function read_char!(l::Lexer)
    if isnothing(l.next)
        return
    else
        ch, state = l.next
        l.next = iterate(l.input, state)
        return ch
    end
end

function next_token!(l::Lexer)
    skip_whitespace!(l)
    ch = read_char(l)
    if ch == '='
        if peek_char(l) == '='
            read_char!(l)
            token = Token(EQ, "==")
        else
            token = Token(ASSIGN, "=")
        end
    elseif ch == ';'
        token = Token(SEMICOLON, ";")
    elseif ch == '('
        token = Token(LPAREN, "(")
    elseif ch == ')'
        token = Token(RPAREN, ")")
    elseif ch == ','
        token = Token(COMMA, ",")
    elseif ch == '!'
        if peek_char(l) == '='
            read_char!(l)
            token = Token(NOT_EQ, "!=")
        else
            token = Token(BANG, "!")
        end
    elseif ch == '+'
        token = Token(PLUS, "+")
    elseif ch == '-'
        token = Token(MINUS, "-")
    elseif ch == '*'
        token = Token(ASTERISK, "*")
    elseif ch == '/'
        token = Token(SLASH, "/")
    elseif ch == '<'
        token = Token(LT, "<")
    elseif ch == '>'
        token = Token(GT, ">")
    elseif ch == '{'
        token = Token(LBRACE, "{")
    elseif ch == '}'
        token = Token(RBRACE, "}")
    elseif ch == '"'
        return read_string!(l)
    elseif ch == '['
        token = Token(LBRACKET, "[")
    elseif ch == ']'
        token = Token(RBRACKET, "]")
    elseif ch == ':'
        token = Token(COLON, ":")
    elseif isnothing(ch)
        token = Token(EOF, "")
    elseif isidentletter(ch)
        return read_identifier!(l)
    elseif isdigit(ch)
        return read_number!(l)
    else
        token = Token(ILLEGAL, string(ch))
    end

    read_char!(l)
    return token
end

function read_identifier!(l::Lexer)
    chars = Char[]
    while isidentletter(read_char(l))
        push!(chars, read_char(l))
        read_char!(l)
    end
    literal = join(chars, "")
    return Token(lookup_indent(literal), literal)
end

function read_number!(l::Lexer)
    chars = Char[]
    while isdigit(read_char(l))
        push!(chars, read_char(l))
        read_char!(l)
    end
    literal = join(chars, "")
    return Token(INT, literal)
end

function read_string!(l::Lexer)
    read_char!(l)

    chars = Char[]
    while read_char(l) != '"' && !isnothing(read_char(l))
        push!(chars, read_char(l))
        read_char!(l)
    end

    read_char!(l)

    return Token(STRING, join(chars, ""))
end

function skip_whitespace!(l::Lexer)
    while isspace(read_char(l))
        read_char!(l)
    end
end

isidentletter(ch) = !isnothing(ch) && ('a' <= ch <= 'z' || 'A' <= ch <= 'Z' || ch == '_')
isspace(ch) = !isnothing(ch) && (ch == ' ' || ch == '\t' || ch == '\n' || ch == '\r')
isdigit(ch) = !isnothing(ch) && ('0' <= ch <= '9')


text:
# Decision trees as found in DecisionTree Julia package.
@reexport module DecisionTreeLearners

using DataFrames
using TSML.TSMLTypes
import TSML.TSMLTypes.fit!
import TSML.TSMLTypes.transform!
using TSML.Utils

import DecisionTree
DT = DecisionTree

export PrunedTree,
       RandomForest,
       Adaboost,
       fit!,
       transform!

# Pruned CART decision tree.

"""
    PrunedTree(
      Dict(
        :purity_threshold => 1.0,
        :max_depth => -1,
        :min_samples_leaf => 1,
        :min_samples_split => 2,
        :min_purity_increase => 0.0
      )
    )

Decision tree classifier.
See [DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparmeters:
- `:purity_threshold` => 1.0 (merge leaves having >=thresh combined purity)
- `:max_depth` => -1 (maximum depth of the decision tree)
- `:min_samples_leaf` => 1 (the minimum number of samples each leaf needs to have)
- `:min_samples_split` => 2 (the minimum number of samples in needed for a split)
- `:min_purity_increase` => 0.0 (minimum purity needed for a split)

Implements `fit!`, `transform!`
"""
mutable struct PrunedTree <: TSLearner
  model
  args
  function PrunedTree(args=Dict())
    default_args = Dict(
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Merge leaves having >= purity_threshold CombineMLd purity.
        :purity_threshold => 1.0,
        # Maximum depth of the decision tree (default: no maximum).
        :max_depth => -1,
        # Minimum number of samples each leaf needs to have.
        :min_samples_leaf => 1,
        # Minimum number of samples in needed for a split.
        :min_samples_split => 2,
        # Minimum purity needed for a split.
        :min_purity_increase => 0.0
      )
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(tree::PrunedTree, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}

Optimize the hyperparameters of `PrunedTree` instance.
"""
function fit!(tree::PrunedTree, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances=convert(Matrix,features)
  end
  impl_args = tree.args[:impl_args]
  tree.model = DT.build_tree(
    labels,
    instances,
    0, # num_subfeatures (keep all)
    impl_args[:max_depth],
    impl_args[:min_samples_leaf],
    impl_args[:min_samples_split],
    impl_args[:min_purity_increase])
  tree.model = DT.prune_tree(tree.model, impl_args[:purity_threshold])
end

"""
    transform!(tree::PrundTree, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Predict using the optimized hyperparameters of the trained `PrunedTree` instance.
"""
function transform!(tree::PrunedTree, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances=convert(Matrix,features)
  end
  return DT.apply_tree(tree.model, instances)
end

# Random forest (CART).

"""
    RandomForest(
      Dict(
        :output => :class,
        :num_subfeatures => 0,
        :num_trees => 10,
        :partial_sampling => 0.7,
        :max_depth => -1
      )
    )

Random forest classification.
See [DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparmeters:
- `:num_subfeatures` => 0  (number of features to consider at random per split)
- `:num_trees` => 10 (number of trees to train)
- `:partial_sampling` => 0.7 (fraction of samples to train each tree on)
- `:max_depth` => -1 (maximum depth of the decision trees)
- `:min_samples_leaf` => 1 (the minimum number of samples each leaf needs to have)
- `:min_samples_split` => 2 (the minimum number of samples in needed for a split)
- `:min_purity_increase` => 0.0 (minimum purity needed for a split)

Implements `fit!`, `transform!`
"""
mutable struct RandomForest <: TSLearner
  model
  args
  function RandomForest(args=Dict())
    default_args = Dict(
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Number of features to train on with trees (default: 0, keep all).
        :num_subfeatures => 0,
        # Number of trees in forest.
        :num_trees => 10,
        # Proportion of trainingset to be used for trees.
        :partial_sampling => 0.7,
        # Maximum depth of each decision tree (default: no maximum).
        :max_depth => -1
      )
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(forest::RandomForest, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}

Optimize the parameters of the `RandomForest` instance.
"""
function fit!(forest::RandomForest, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances=convert(Matrix,features)
  end
  # Set training-dependent options
  impl_args = forest.args[:impl_args]
  # Build model
  forest.model = DT.build_forest(
    labels,
    instances,
    impl_args[:num_subfeatures],
    impl_args[:num_trees],
    impl_args[:partial_sampling],
    impl_args[:max_depth]
  )
end

"""
    transform!(forest::RandomForest, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Predict using the optimized hyperparameters of the trained `RandomForest` instance.
"""
function transform!(forest::RandomForest, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances=convert(Matrix,features)
  end
  return DT.apply_forest(forest.model, instances)
end

# Adaboosted decision stumps.

"""
    Adaboost(
      Dict(
        :output => :class,
        :num_iterations => 7
      )
    )

Adaboosted decision tree stumps. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:
- `:num_iterations` => 7 (number of iterations of AdaBoost)

Implements `fit!`, `transform!`
"""
mutable struct Adaboost <: TSLearner
  model
  args
  function Adaboost(args=Dict())
    default_args = Dict(
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Number of boosting iterations.
        :num_iterations => 7
      )
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(adaboost::Adaboost, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}

Optimize the hyperparameters of `Adaboost` instance.
"""
function fit!(adaboost::Adaboost, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances = convert(Matrix,features)
  end
  # NOTE(svs14): Variable 'model' renamed to 'ensemble'.
  #              This differs to DecisionTree
  #              official documentation to avoid confusion in variable
  #              naming within CombineML.
  ensemble, coefficients = DT.build_adaboost_stumps(
    labels, instances, adaboost.args[:impl_args][:num_iterations]
  )
  adaboost.model = Dict(
    :ensemble => ensemble,
    :coefficients => coefficients
  )
end

"""
    transform!(adaboost::Adaboost, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Predict using the optimized hyperparameters of the trained `Adaboost` instance.
"""
function transform!(adaboost::Adaboost, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  instances = features
  if typeof(features) <: DataFrame
    instances = convert(Matrix,features)
  end
  return DT.apply_adaboost_stumps(
    adaboost.model[:ensemble], adaboost.model[:coefficients], instances
  )
end

end # module


text:
function xlogx(x)
    x*log(x)
end


text:
module TestDatumStorage

using Mimi
using Test

comp_first = 2003
comp_last = 2008

@defcomp foo begin
    v = Variable(index = [time])
    function run_timestep(p, v, d, ts)
        # implement "short component" via time checking
        if TimestepValue(comp_first) <= ts <= TimestepValue(comp_last)
            v.v[ts] = gettime(ts)
        end
    end
end

@defcomp bar begin
    region = Index()
    v = Variable(index = [time, region])
    function run_timestep(p, v, d, ts)
        v.v[ts, :] .= gettime(ts)
    end
end

years = 2001:2010
regions = [:A, :B]

nyears = length(years)
nregions = length(regions)

#------------------------------------------------------------------------------
# 1. Single dimension case, fixed timesteps
#------------------------------------------------------------------------------

m = Model()
set_dimension!(m, :time, years)
add_comp!(m, foo)

run(m)
v = m[:foo, :v]
@test length(v) == nyears # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years)
    if y < comp_first || y > comp_last
        @test ismissing(v[i])
    else
        @test v[i] == y
    end
end

#------------------------------------------------------------------------------
# 2. Multi-dimension case, fixed timesteps
#------------------------------------------------------------------------------

m2 = Model()

@defcomp baz begin
    region = Index()
    v = Variable(index = [time, region])
    function run_timestep(p, v, d, ts)
        # implement "short component" via time checking
        if TimestepValue(comp_first) <= ts <= TimestepValue(comp_last)
            v.v[ts, :] .= gettime(ts)
        end
    end
end

set_dimension!(m2, :time, years)
set_dimension!(m2, :region, regions)

add_comp!(m2, baz)

run(m2)
v2 = m2[:baz, :v]
@test size(v2) == (nyears, nregions) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years)
    if y < comp_first || y > comp_last
        [@test ismissing(v2[i, j]) for j in 1:nregions]
    else
        [@test v2[i, j]==y for j in 1:nregions]
    end
end

#------------------------------------------------------------------------------
# 3. Single dimension case, variable timesteps
#------------------------------------------------------------------------------

years_variable = [2000:2004..., 2005:5:2030...]
foo2_first = 2003
foo2_last = 2010

m = Model()
set_dimension!(m, :time, years_variable)

@defcomp foo2 begin
    v = Variable(index = [time])
    function run_timestep(p, v, d, ts)
        # implement "short component" via time checking
        if TimestepValue(foo2_first) <= ts <= TimestepValue(foo2_last)
            v.v[ts] = gettime(ts)
        end
    end
end

add_comp!(m, foo2)

run(m)
v = m[:foo2, :v]
@test length(v) == length(years_variable) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years_variable)
    if y < foo2_first || y > foo2_last
        @test ismissing(v[i])
    else
        @test v[i] == y
    end
end

#------------------------------------------------------------------------------
# 4. Multi-dimension case, variable timesteps
#------------------------------------------------------------------------------

m2 = Model()

buz_first = 2003
buz_last = 2010

@defcomp buz begin
    region = Index()
    v = Variable(index = [time, region])
    function run_timestep(p, v, d, ts)
        # implement "short component" via time checking
        if TimestepValue(buz_first) <= ts <= TimestepValue(buz_last)
            v.v[ts, :] .= gettime(ts)
        end
    end
end

set_dimension!(m2, :time, years_variable)
set_dimension!(m2, :region, regions)
add_comp!(m2, buz)
run(m2)
v2 = m2[:buz, :v]
@test size(v2) == (length(years_variable), nregions) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years_variable)
    if y < buz_first || y > buz_last
        [@test ismissing(v2[i, j]) for j in 1:nregions]
    else
        [@test v2[i, j]==y for j in 1:nregions]
    end
end

end # module


text:
"Add some style overrides to make things a bit prettier and more consistent with Pluto."
function _add_documenter_css(html)
    style = """

        """
    return string(style, '\n', html)
end

"""
Wrap the body text in the proper styling.
Regex matching won't work here because the body can contain arbitrary HTML.

This function is applied after the other conversions, so we know that the first and last `p` tags must denote the start and end point of the body.
"""
function _wrap_admonition_body(html::AbstractString)
    start_p = first(findfirst("
", html)::UnitRange{Int})::Int
    end_p = last(findlast("
", html)::UnitRange{Int})::Int
    body_with_p = html[start_p:end_p]

    function wrap_body(body_with_p)
        return rstrip("""

                $body_with_p

            """)
    end

    replace(html, body_with_p => wrap_body)
end

"Convert a single admonition from Pluto to Documenter."
function _convert_admonition(html::AbstractString)
    html = replace(html, """

        """<div class="admonition is-""")
    html = replace(html, r"""
([^<]*)
""" =>
        s"""\n  \n    \1\n  \n""")
    html = _wrap_admonition_body(html)
    return html
end

"""
Convert Pluto's admonitions HTML to Documenter's admonitions HTML.
This ensures that admonitions are properly styled in `documenter_output`.
"""
function _convert_admonitions(html::AbstractString)
    rx = r"""
[^\n]*\n[^\n]*"""
    return replace(html, rx => _convert_admonition)
end


text:
using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Latexify, LaTeXStrings

######### Question 1 ################

# 1-Initializing variables and practice with basic matrix operations

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

#Create the following four matrices of random numbers, setting the seed to ’1234’ .
#Name the matrices and set the dimensions as noted

# create 10x7 distruburted U[-5,10]
Random.seed!(1234)
A = rand(10,7)  .* (10 - (-5)) .- 5

# create 10x7 distruburted random numbers distributed N (−2, 15) [st dev is 15]
Random.seed!(1234)
B = (-2) .+ (randn(10,7) * 15)

# create C 5×7 - the first 5 rows and first 5 columns of A and the last two columns and first 5 rows of B
sub_A = A[1:5,1:5]; sub_B = B[1:5,6:7]
C = hcat(sub_A, sub_B)

### D 10×7 - where Di,j = Ai,j if Ai,j ≤ 0, or 0 otherwise
D = similar(A)

for i in 1:10, j in 1:7
    if A[i,j] <= 0
        D[i,j]=A[i,j]
    else
        D[i,j] = 0
    end
end

display(D)

# or D = A .* (A .≤ 0)

## Use a built-in Julia function to list the number of elements of A

for i in 1:10, j in 1:7
    println(A[i,j])
end

## Use a series of built-in Julia functions to list the number of unique elements of D

vector_D = D[:]; unique_D = unique(vector_D);
println(unique_D)
print(length(unique_D))

## Using the reshape() function, create a new matrix called E which is the ‘vec’ operator applied to B. Can you find an easier way to accomplish this?

E = reshape(B,:,1)

# easier way
# E  = B[:]

## Create a new array called F which is 3-dimensional and contains A in the first column of the third dimension and B in the second column of the third dimension

F = cat(A, B, dims=3)
# size(F)

## Use the permutedims() function to twist F so that it is now F2×10×7 instead of F10×7×2. Save this new matrix as F.

# permutedims(F, [2, 10, 7])
F = reshape(F,(2,10,7))

## Create a matrix G which is equal to B ⊗ C (the Kronecker product of B and C). What happens when you try C ⊗ F ?

G = kron(B,C);

# kron(C, F)  # Doesn't work because they are of different # dimensions

## Save the matrices A, B, C, D, E, F and G as a .jld file named matrixpractice.

matrixpractice = Dict("A"=> A, "B"=> B, "C"=> C, "D"=> D, "E"=> E, "F"=> F, "G"=> G)
@save("matrixpractice.jld", matrixpractice)

## Save only the matrices A, B, C, and D as a .jld file called firstmatrix.

firstmatrix = Dict("A"=> A, "B"=> B, "C"=> C, "D"=> D)
@save("firstmatrix.jld", firstmatrix)

## Export C as a .csv file called Cmatrix. You will first need to transform C into a DataFrame.

typeof(C)
Cmatrix = convert(DataFrame, C)
CSV.write("Cmatirx.csv", Cmatrix)

## Export D as a tab-delimited .dat file called Dmatrix. You will first need to transform D into a DataFrame.

typeof(D)
Dmatrix = convert(DataFrame, D)
@save("Dmatrix.dat", Dmatrix, delim='\t')

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

# Wrap a function definition around all of the code for question 1.
# Call the function q1(). The function should have 0 inputs and should output the arrays
# A, B, C and D. At the very bottom of your script you should add the code A,B,C,D = q1().

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

function q1()
    # 1.a
    Random.seed!(1234);

    A = rand(10,7)  .* (10 - (-5)) .- 5;

    B = (-2) .+ (randn(10,7) * 15);

    C = zeros(5, 7);

    sub_A = A[1:5,1:5]; sub_B = B[1:5,6:7];
    C = hcat(sub_A, sub_B);

    D = A .* (A .≤ 0);

    D = similar(A);
    for i in 1:10, j in 1:7
        if A[i,j] <= 0
            D[i,j]=A[i,j]
        else
            D[i,j] = 0
        end
    end

    for i in 1:10, j in 1:7
        println(A[i,j])
    end

    vector_D = D[:]; unique_D = unique(vector_D);
    print(length(unique_D))

    E = reshape(B,:,1);

    F = cat(A, B, dims=3);
    F = reshape(F,(2,10,7));

    G = kron(B,C);

    matrixpractice = Dict("A"=> A, "B"=> B, "C"=> C, "D"=> D, "E"=> E, "F"=> F, "G"=> G);
    @save("matrixpractice.jld", matrixpractice)

    Cmatrix = convert(DataFrame, C);
    CSV.write("Cmatirx.csv", Cmatrix)

    Dmatrix = convert(DataFrame, D);
    @save("Dmatrix.dat", Dmatrix, delim='\t')

    return A, B, C, D
end

A,B,C,D = q1()

############## Question 2 ###############

# Practice with loops and comprehensions

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

# Write a loop or use a comprehension that computes the element-by-element product of A and B.
# Name the new matrix AB. Create a matrix called AB2 that accomplishes this ask without a loop or comprehension.

AB = zeros(size(A)...);

for i = 1:size(A)[1], j = 1:size(A)[2]
    AB[i, j] = A[i, j] * B[i, j];
end

AB2 = A .* B;

# Write a loop that creates a column vector called Cprime which contains only
#the elements of C that are between -5 and 5 (inclusive). Create a vector called
# Cprime2 which does this calculation without a loop.

Cprime = Vector{Float64}();
for item in eachrow(C), i in item
    if   -5 ≤  i ≤ 5
        append!(Cprime,i)
    end
end

Cprime2 = C[-5 .≤ C .≤ 5];

# Using loops or comprehensions, create a 3-dimensional array called X that is of
# dimen- sionN×K×TwhereN=15,169,K=6,andT=5. Forallt,thecolumnsofX should be (in order)

N = 15_169;
K = 6;
T = 5;
X = zeros(N, K, T);
binord(;n=20, p=0.6) = sum(rand(n) .≤ p);

for t in 1:T
    X[:, 1, t] .= 1;
    X[:, 2, t] = rand(N, 1) .≤ 0.75 * (6 - t) / 5;
    X[:, 3, t] = randn(N, 1) .* 5(t - 1) .+ (15 + t - 1);
    X[:, 4, t] = randn(N, 1) .* (1 / ℯ) .+ π * (6 - t) / 3;
    X[:, 5, t] = [binord() for _ in 1:N];
    X[:, 6, t] = [binord(p=0.5) for _ in 1:N];
end

# Use comprehensions to create a matrix β which is K × T and whose elements evolve across time in the following fashion:

β = zeros(K, T);
β[1, :] = [1 + 0.25(t - 1) for t in 1:T];
β[2, :] = [log(t) for t in 1:T];
β[3, :] = [-sqrt(t) for t in 1:T];
β[4, :] = [exp(t) - exp(t + 1) for t in 1:T];
β[5, :] = [t for t in 1:T];
β[6, :] = [t / 3 for t in 1:T];

# Use comprehensions to create a matrix Y which is N × T defined by Yt = Xt βt + εt , iid

ε = Random.randn(N, T) .* 0.36;
Y = [sum(X[n, :, t] .* β[:, t]) + ε[n, t] for n in 1:N, t in 1:T];

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

# Wrap a function definition around all of the code for question 2.
# Call the function q2(). The function should have take as inputs the arrays A, B and C.
# It should return nothing. At the very bottom of your script you should add the code q2(A,B,C,D).
#Make sure q2() gets called after q1()!

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

function q2(A, B, C, D)

    AB = zeros(size(A)...);

    for i = 1:size(A)[1], j = 1:size(A)[2]
        AB[i, j] = A[i, j] * B[i, j];
    end

    AB2 = A .* B;

    Cprime = Vector{Float64}();

    for item in eachrow(C), i in item
        if   -5 ≤  i ≤ 5
            append!(Cprime,i)
        end
    end

    Cprime2 = C[-5 .≤ C .≤ 5];

    N = 15_169;
    K = 6;
    T = 5;

    X = zeros(N, K, T);
    Random.seed!(1234)
    binord(;n=20, p=0.6) = sum(rand(n) .≤ p);

    for t in 1:T
        X[:, 1, t] .= 1;
        X[:, 2, t] = rand(N, 1) .≤ 0.75 * (6 - t) / 5;
        X[:, 3, t] = randn(N, 1) .* 5(t - 1) .+ (15 + t - 1);
        X[:, 4, t] = randn(N, 1) .* (1 / ℯ) .+ π * (6 - t) / 3;
        X[:, 5, t] = [binord() for _ in 1:N];
        X[:, 6, t] = [binord(p=0.5) for _ in 1:N];
    end

    β = zeros(K, T);
    β[1, :] = [1 + 0.25(t - 1) for t in 1:T];
    β[2, :] = [log(t) for t in 1:T];
    β[3, :] = [-sqrt(t) for t in 1:T];
    β[4, :] = [exp(t) - exp(t + 1) for t in 1:T];
    β[5, :] = [t for t in 1:T];
    β[6, :] = [t / 3 for t in 1:T];

    ε = randn(N, T) .* 0.36;
    Y = [sum(X[n, :, t] .* β[:, t]) + ε[n, t] for n in 1:N, t in 1:T];

    return nothing
end

q2(A, B, C, D)

############## Question 3 ###############

# Reading in Data and calculating summary statistics

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

# Clear the workspace and import the file nlsw88.csv into Julia as a DataFrame.

data = CSV.read("/Users/firatmelihyilmaz/fall-2020/ProblemSets/PS1-julia-intro/nlsw88.csv")
describe(data)
@save("nlsw88.jld", data)

#  What percentage of the sample has never been married? What percentage are college graduates?

ptc_never_married = (sum(data["never_married"])/(sum(data["never_married"]) + sum(data["married"]))) * 100;
println("Percentage percentage of the sample has never been married: $ptc_never_married %")

# Use the tabulate command to report what percentage of the sample is in each race category

races = prop(freqtable(data["race"]));
# (prop ∘ freqtable)(data["race"])

# Use the describe() function to create a matrix called summarystats which lists the mean,
# median, standard deviation, min, max, number of unique elements, and in- terquartile range
# (75th percentile minus 25th percentile) of the data frame. How many grade observations are missing?

summarystats = describe(data);

# Show the joint distribution of industry and occupation using a cross-tabulation.

joint_table = freqtable(data, :industry, :occupation);

# Wrap a function definition around all of the code for question 3.
# Call the function q3(). The function should have no inputs and no outputs.
# At the very bottom of your script you should add the code q3().

function q3()
    data = CSV.read("/Users/firatmelihyilmaz/fall-2020/ProblemSets/PS1-julia-intro/nlsw88.csv");
    @save("nlsw88.jld", data)
    ptc_never_married = (sum(data["never_married"])/(sum(data["never_married"]) + sum(data["married"]))) * 100;
    races = prop(freqtable(data["race"]));
    summarystats = describe(data);
end

############## Question 4 ###############

# Practice with functions

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

# Load firstmatrix.jld

@load("firstmatrix.jld")

# Write a function called matrixops that takes as inputs the matrices A and B
# from question (a) of problem 1 and has three outputs: (i) the element-by-element product of the inputs,
# (ii) the product A′B, and (iii) the sum of all the elements of A + B.

function matrixops(A, B)
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end

matrixops(A,B);

function matrixops(A, B)
    """
    Takes  matrix A and B as inputs and returns following outputs:

    i = A .+ B
    ii = A' * B
    iii = sum(A + B)
    """
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end

function matrixops(A, B)
    """
    Takes  matrix A and B as inputs and returns following outputs:

    i = A .+ B
    ii = A' * B
    iii = sum(A + B)
    """
    if size(A) != size(B)
        throw(DimensionMismatch("inputs must have the same size."))
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end