text:
# Model
include_model("hopper")

# Horizon
T = 31

# Time step
tf = 3.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
z_h = 0.5
q1 = [0.0, 0.5, 0.0, 0.5]
v1 = [0.0, 0.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q1]
qT = [3.0, 0.5, 0.0, 0.5]

xl, xu = state_bounds(model, T, [model.qL; model.qL], [model.qU; model.qU],
    x1 = x1, xT = [Inf*ones(model.nq); qT])

# Objective
include_objective("velocity")
obj_velocity = velocity_objective([Diagonal(ones(model.nq)) for t = 1:T],
	model.nq, idx_angle = (3:3), h = h)
obj_tracking = quadratic_tracking_objective(
    [Diagonal(zeros(model.n)) for t = 1:T],
    [Diagonal([1.0e-1, 1.0e-1, zeros(model.m - model.nu)...]) for t = 1:T-1],
    [zeros(model.n) for t = 1:T],
    [zeros(model.m) for t = 1:T-1])
obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity])

# Constraints
include_constraints("contact")
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = configuration_to_state(linear_interpolation(q1, qT, T)) # linear interpolation on state
u0 = [1.0e-3 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()
@time z̄ , info = solve(prob, copy(z0),
 	nlp = :SNOPT7,
	tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)


text:
using LinearAlgebra
using SparseArrays
using Plots

"This routine solves the heat equation using Forward Euler in time and finite differences
in space. Note that this is unstable unless dt ~ O(h^2)."

m = 100 # number of points
T = .01
dt = .0005

u0(x) = 0.0 # initial condition
f(x) = 5*(Float64((x > -.5) && (x <= 0.0)) - Float64((x < .5) && (x >= 0.0)))
f(x,t) = f(x)*exp(-t)
α(t) = 1.0
β(t) = pi

# define spatial grid
x = LinRange(-1,1,m+2)
xint = x[2:end-1]
h = x[2]-x[1]
A = (1/h^2) * spdiagm(0=>2*ones(m),-1=>-ones(m-1),1=>-ones(m-1))

function F(t)
    b = f.(xint,t) # look for f(x,t)
    b[1] += α(t)/h^2
    b[end] += β(t)/h^2
    return b
end

u = u0.(xint)
Nsteps = ceil(Int,T/dt)
dt = T / Nsteps

interval = 1
@gif for k = 1:Nsteps
    tk = k*dt
    u .= u + dt * (F(tk) - A*u)
    if k % interval==0
        plot(xint,u,linewidth=2,label="Solution",ylims=(-1.0,3.0))
        println("on timestep $k out of $Nsteps.")
    end
end every interval

ext: jl


text:
using Random, Test
using CommonLH.MatrixLH

@testset "MatrixLH" begin
	@testset "Accumarray" begin
		subsM = [1 1; 2 1; 1 1; 2 2];
		valV = [1, 2, 3, 4];
		A = MatrixLH.accumarray(subsM, valV, sum);
		@test A == [4.0 0.0; 2.0 4.0]

		# Vector
		A = MatrixLH.accumarray(subsM[:,1], valV, sum);
		@test A == [4.0, 6.0]

		# Tim Holy version that is faster but only works on vectors
		A2 = MatrixLH.accumarrayTH(subsM[:,1], valV)
		@test A2 == A
	end

	@testset "RoundToGrid" begin
		gridV = collect(range(1.0, 9.5, length = 7));
		Random.seed!(1234);
		xM = rand(4,3,2);
		idxM = round_to_grid(xM, gridV);

		isValid = true;
		for i1 in 1 : size(xM)[3]
			for i2 in 1 : size(xM)[2]
				for i3 in 1 : size(xM)[1]
					d1V = abs.(xM[i3,i2,i1] .- gridV);
					d2 = abs(xM[i3,i2,i1] - gridV[idxM[i3,i2,i1]]);
					isValid = isValid && all(d2 .< d1V .+ 1e-5)
				end
			end
		end
		@test isValid
	end
end

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


text:
using Distributed

include("$(@__DIR__)/../config/config.jl")
addprocs(num_workers; exeflags="--project=$(abspath(joinpath(@__DIR__, "..")))")

@everywhere begin
    modulepath = abspath(joinpath(@__DIR__, "../src/KernelLearning"));
    push!(LOAD_PATH, modulepath);
    using KernelLearning
end

@everywhere using OneClassActiveLearning, SVDD
@everywhere using MLKernels
@everywhere using JSON, Serialization
@everywhere using Logging
@everywhere using Random

@everywhere function evaluate_parameters(gamma, C, data, labels, solver)
    model = VanillaSVDD(data)
    SVDD.initialize!(model, SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(gamma)), FixedCStrategy(C)))
    set_adjust_K!(model, true)
    SVDD.fit!(model, solver)
    return ConfusionMatrix(SVDD.classify.(SVDD.predict(model, data)), labels)
end

@everywhere function run_experiment(e::Dict{Symbol, Any})
    Random.seed!(e[:seed])

    io = open(joinpath(e[:log_dir], "$(e[:hash]).log"), "a")
    logger = SimpleLogger(io)
    global_logger(logger)

    if isfile(e[:output_file])
        @info("Skipping $(e[:hash])")
        return nothing
    end

    @info("Running $(e[:hash]), $(e[:data_set_name]), $(string(e[:method]))")
    @show e[:hash], e[:data_set_name], e[:method]

    result = deepcopy(e)
    result[:method] = string(e[:method])
    data, labels = load_data(e[:data_file])

    result[:datastats] = Dict(:n_observations => size(data, 2),
                               :n_dims => size(data, 1),
                               :outlier_ratio => sum(labels .== :outlier) / size(data, 2))

    gamma, C = KernelLearning.find_parameters(e[:method], data, labels, e[:oracle], e[:solver], e[:quality_metrics][e[:selection_metric]], e[:seed])
    cm = evaluate_parameters(gamma, C, data, labels, e[:solver])
    result[:result] = Dict(:gamma => gamma, :C => C, :metrics => Dict{Symbol, Float64}())
    for (metric_name, metric) in e[:quality_metrics]
        m = metric(cm)
        result[:result][:metrics][metric_name] = m
        @info("$metric_name: $m")
    end

    open(e[:output_file], "w") do f
        JSON.print(f, result)
    end

    close(io)
    return nothing
end

exp_dir = joinpath(data_output_root, exp_name)
@info("Running experiment in $exp_dir")

# store experiments
experiments = deserialize(open(joinpath(exp_dir, "experiments.jser")))
open(joinpath(exp_dir, "experiments.json"), "w") do f
   JSON.print(f, experiments)
end

# run experiments
pmap(run_experiment, experiments, on_error=ex->print(ex))

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


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

using Markdown
using InteractiveUtils

# ╔═╡ 6316f37d-24d7-488a-9916-47a75512ea54
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 9fa2be99-1c4a-4e53-8a01-546a3670f19f
begin
    using FFTW
    using Flux
    using Plots
    plotly()
end

# ╔═╡ b682a76e-75c6-4e0f-b542-210725501e5b
html"""<style>
main {
    max-width: 700pt;
}
"""

# ╔═╡ b5ba5d3a-eb07-4d5f-a802-6d85075ed528
md"
# Fourier Net

*“Can neural networks learn the Fourier transform?”*

This notebook is a implementation of [The Fourier transform is a neural network](https://sidsite.com/posts/fourier-nets/).
"

# ╔═╡ 21b015ea-aa4b-11eb-182b-a7590312f5ba
md"
**By the definition of DFT:**

A single amplitude for a singal duration of frequency can be formulated as follow:

$x̂_k = \sum_{n=0}^{N-1} x_n \cdot exp(\frac{-i2\pi k}{N}n)$

$x̂_k =
\begin{bmatrix} x_0 & x_1 & ... & x_{N-1} \end{bmatrix}
\begin{bmatrix}
    exp(0) \\
    exp(\frac{-i 2 \pi k}{N}1) \\
    exp(\frac{-i 2 \pi k}{N}2) \\
    exp(\frac{-i 2 \pi k}{N}3) \\
    exp(\frac{-i 2 \pi k}{N}4) \\
    ... \\
    exp(\frac{-i 2 \pi k}{N}(N-1))
\end{bmatrix}$

Therefore, N of amplitide for N durations of frequency can be illustrated like:

$\begin{bmatrix} x̂_0 \\ ... \\ x̂_{N-1} \end{bmatrix} =
\begin{bmatrix}
    exp(0) & exp(0) & exp(0) & ... & exp(0) \\
    exp(0) & exp(\frac{-i2\pi}{N}1) & exp(\frac{-i2\pi 2}{N}1) & ... & exp(\frac{-i2\pi (N-1)}{N}1) \\
    exp(0) & exp(\frac{-i2\pi}{N}2) & exp(\frac{-i2\pi 2}{N}2) & ... & exp(\frac{-i2\pi (N-1)}{N}2) \\
    exp(0) & exp(\frac{-i2\pi}{N}3) & exp(\frac{-i2\pi 2}{N}3) & ... & exp(\frac{-i2\pi (N-1)}{N}3) \\
    exp(0) & exp(\frac{-i2\pi}{N}4) & exp(\frac{-i2\pi 2}{N}4) & ... & exp(\frac{-i2\pi (N-1)}{N}4) \\
    ... \\
    exp(0) & exp(\frac{-i2\pi}{N}(N-1)) & exp(\frac{-i2\pi 2}{N}(N-1)) & ... & exp(\frac{-i2\pi (N-1)}{N}(N-1))
\end{bmatrix}_{N \times N}^T
\begin{bmatrix} x_0 \\ ... \\ x_{N-1} \end{bmatrix}$

Or in short:

$x̂⃗ = 𝐖_{ComplexFourier}^T x⃗$
"

# ╔═╡ 212ac582-1270-49da-9cd4-078ab695be66
md"
## Generate random signal

$x⃗ = \begin{bmatrix} x_0 \\ ... \\ x_{N-1} \end{bmatrix}$
"

# ╔═╡ 7efa004c-57c4-4b9c-bc74-90252de25218
begin
    len = 128
    x = rand(len)
end;

# ╔═╡ 0d6607c2-0881-40f7-81fd-494d49faa697
md"
## Computing the Fourier transform using Fourier weights
"

# ╔═╡ eeded1aa-3edd-4aa7-b147-9c69e2293267
md"
**Create Fourier weighrts**

$𝐖_{ComplexFourier}^T = \begin{bmatrix}
    exp(0) & exp(0) & exp(0) & ... & exp(0) \\
    exp(0) & exp(\frac{-i2\pi}{N}1) & exp(\frac{-i2\pi 2}{N}1) & ... & exp(\frac{-i2\pi (N-1)}{N}1) \\
    exp(0) & exp(\frac{-i2\pi}{N}2) & exp(\frac{-i2\pi 2}{N}2) & ... & exp(\frac{-i2\pi (N-1)}{N}2) \\
    exp(0) & exp(\frac{-i2\pi}{N}3) & exp(\frac{-i2\pi 2}{N}3) & ... & exp(\frac{-i2\pi (N-1)}{N}3) \\
    exp(0) & exp(\frac{-i2\pi}{N}4) & exp(\frac{-i2\pi 2}{N}4) & ... & exp(\frac{-i2\pi (N-1)}{N}4) \\
    ... \\
    exp(0) & exp(\frac{-i2\pi}{N}(N-1)) & exp(\frac{-i2\pi 2}{N}(N-1)) & ... & exp(\frac{-i2\pi (N-1)}{N}(N-1))
\end{bmatrix}_{N \times N}^T$
"

# ╔═╡ 1f66b59e-bdcd-4cdb-9f28-3bbf0f7dd08e
function create_𝐰_fourier(len::Integer)
    𝐰 = Matrix{ComplexF64}(undef, len, len)
    for n in 0:len-1, k in 0:len-1
        𝐰[n+1, k+1] = exp(-im*2π*k*n/len)
    end

    return 𝐰'
end;

# ╔═╡ 036016cd-e4f5-40f7-96d9-cac8fdbebd01
md"
**Do the Fourier transform by Fourier weights and FFT**

$x̂⃗ = 𝐖_{ComplexFourier}^T x⃗$

$x̂⃗_{fft} = 𝐹x⃗$
"

# ╔═╡ d5933fdb-b9ed-4475-af43-58dac24a47d5
begin
    # Fourier transform via Fourier weights
    𝐰 = create_𝐰_fourier(len)
    x̂ = 𝐰 * x

    # Fourier transform via fft
    x̂_fft = fft(x)

    # root-mean-sqrt error between those two
    ϵ_x̂ = sqrt(sum(abs.(x̂_fft - x̂).^2)/len)
end;

# ╔═╡ 9d937bef-f325-4299-8c4b-7ae358aeb1c9
md"
Spectrum's RMSE: $ϵ_x̂
"

# ╔═╡ f65c7059-9c2e-4993-a7e6-32cb9a848f23
begin
    plot(title="Spectrum", xlabel="Freq (arb. unit)", ylabel="Amp (arb. unit)")
    plot!(real(x̂_fft), label="FFT")
    plot!(real(x̂), label="Fourier weights")
end

# ╔═╡ 6e200b44-97fb-43e0-8026-aabef11f4f17
md"
**Reconstruct signal by inverse Fourier weights and iFFT**

$x̂_k = \sum_{n=0}^{N-1} x_n \cdot exp(\frac{-i2\pi k}{N}n)$

$x_t = \frac{1}{N} \sum_{k=0}^{N-1} x̂_k \cdot exp(\frac{i2\pi k}{N}n)$
"

# ╔═╡ 45844ffc-9f5c-4952-abde-4ce770ab77f6
function create_𝐰¯¹_fourier(len::Integer)
    𝐰¯¹ = Matrix{ComplexF64}(undef, len, len)
    for n in 0:len-1, k in 0:len-1
        𝐰¯¹[n+1, k+1] = exp(im*2π*k*n/len)
    end

    return (𝐰¯¹/len)'
end;

# ╔═╡ 856e62d3-53c8-49a9-9061-d223dbb92ef1
md"
**Do the inverse Fourier transform by inverse Fourier weights and iFFT**

$x⃗ = 𝐖¯¹_{ComplexFourier}^T x̂⃗$

$x⃗_{fft} = 𝐹¯¹x⃗$
"

# ╔═╡ ac36e1c5-dade-4b0f-974e-2d261ea260f8
begin
    # reconstruct signal via Fourier weights
    𝐰¯¹ = create_𝐰¯¹_fourier(len)
    𝑓¯¹x̂ = real(𝐰¯¹ * x̂)

    # reconstruct signal via ifft
    𝑓¯¹x̂_fft = real(ifft(x̂_fft))

    # root-mean-sqrt error between those two
    ϵ_𝑓¯¹x̂_fft = sqrt(sum(abs.(𝑓¯¹x̂_fft - x).^2)/len)
    ϵ_𝑓¯¹x̂ = sqrt(sum(abs.(𝑓¯¹x̂ - x).^2)/len)
end;

# ╔═╡ 7339df20-651e-42ac-8286-e9696fb72458
md"
Inverse Fourier transform's RMSE:

* Via IFFT: $ϵ_𝑓¯¹x̂_fft
* Via Fourier waights: $ϵ_𝑓¯¹x̂
"

# ╔═╡ e9c246d8-ae77-4236-8ca0-15de7cced221
begin
    plot(title="Signal", xlabel="Time (arb. unit)", ylabel="Amp (arb. unit)", ticks=[], ylims=(0, 3.5))
    plot!(x .+2 , label="Original signal")
    plot!(𝑓¯¹x̂ .+ 1, label="Reconstructed signal (Fourier weights)")
    plot!(𝑓¯¹x̂_fft, label="Reconstructed signal (ifft)")
end

# ╔═╡ 3c03ec6d-8994-4886-9006-3cddbf0ac027
md"
## Learning the Fourier transform via gradient-descent

We can consider the discrete Fourier transform (DFT) to be an artificial neural network: it is a single layer network, with no bias, no activation function, and particular values for the weights. The number of output nodes is equal to the number of frequencies we evaluate.
"

# ╔═╡ 0b4853ce-8dd0-4fe5-881d-cdf5d0289cc8
md"
**Build Fourier Net**
"

# ╔═╡ 1cbd900c-b74b-41c8-a57c-6a6ae644b88c
begin
    mutable struct FourierNet
        𝐰
    end

    Flux.@functor FourierNet

    (m::FourierNet)(x) = m.𝐰 * x

    loss(m::FourierNet, x, x̂) = sum(abs2, x̂ .- m(x)) / len
end;

# ╔═╡ 0abb581e-6835-4dfc-bd01-5b2b6306d242
md"
**Initiating the training process**
"

# ╔═╡ 895ee250-e3ac-4879-b8ee-f1ba8729adaa
begin
    # define the net and loss function
    m = FourierNet(rand(ComplexF64, len, len))
    loss(x, x̂) = loss(m, x, x̂)

    # prepare data
    xs = rand(ComplexF64, len, 50000)
    x̂s = 𝐰 * xs
    data = Flux.Data.DataLoader((xs, x̂s), batchsize=15, shuffle=true)

    # define call-back function for plotting
    training_loss = []
    function cb()
        lossₜ = loss(xs, x̂s)
        push!(training_loss, lossₜ)
    end

    # training process
    train_𝐰!(n, η) = Flux.@epochs n Flux.Optimise.train!(
        loss,
        params(m),
        data,
        Descent(η),
        # cb=Flux.throttle(@show(loss(xs, x̂s)), 1, leading=false, trailing=true)
        cb=Flux.throttle(cb, 1, leading=false, trailing=true)
    )

    # start training
    train_𝐰!(15, 1e-1)
    train_𝐰!(30, 1e-2)
    train_𝐰!(20, 1e-3)
end;

# ╔═╡ 3c1157bd-9780-429a-9cbf-cd5a18d676da
plot(training_loss, title="Loss", xlabel="time (sec)", ylabel="loss", label="training")

# ╔═╡ 9ef928f0-5577-4507-8fbe-193ce6925716
md"
Weights' residual: $(sum(abs, 𝐰 - m.𝐰)/length(𝐰))

Weights' RMSE: $(sqrt(sum(abs2, 𝐰 - m.𝐰)/length(𝐰)))
"

# ╔═╡ 04189afa-6cdd-45f4-8eaa-12140bde3ada
begin
    x̂_grad = m(x)

    plot(title="Spectrum", xlabel="Freq (arb. unit)", ylabel="Amp (arb. unit)")
    plot!(real(x̂_grad), label="Fourier weights by GD")
    plot!(real(x̂_fft), label="FFT")
end

# ╔═╡ Cell order:
# ╟─b682a76e-75c6-4e0f-b542-210725501e5b
# ╟─b5ba5d3a-eb07-4d5f-a802-6d85075ed528
# ╠═6316f37d-24d7-488a-9916-47a75512ea54
# ╠═9fa2be99-1c4a-4e53-8a01-546a3670f19f
# ╟─21b015ea-aa4b-11eb-182b-a7590312f5ba
# ╟─212ac582-1270-49da-9cd4-078ab695be66
# ╠═7efa004c-57c4-4b9c-bc74-90252de25218
# ╟─0d6607c2-0881-40f7-81fd-494d49faa697
# ╟─eeded1aa-3edd-4aa7-b147-9c69e2293267
# ╠═1f66b59e-bdcd-4cdb-9f28-3bbf0f7dd08e
# ╟─036016cd-e4f5-40f7-96d9-cac8fdbebd01
# ╠═d5933fdb-b9ed-4475-af43-58dac24a47d5
# ╟─9d937bef-f325-4299-8c4b-7ae358aeb1c9
# ╠═f65c7059-9c2e-4993-a7e6-32cb9a848f23
# ╟─6e200b44-97fb-43e0-8026-aabef11f4f17
# ╠═45844ffc-9f5c-4952-abde-4ce770ab77f6
# ╟─856e62d3-53c8-49a9-9061-d223dbb92ef1
# ╠═ac36e1c5-dade-4b0f-974e-2d261ea260f8
# ╟─7339df20-651e-42ac-8286-e9696fb72458
# ╠═e9c246d8-ae77-4236-8ca0-15de7cced221
# ╟─3c03ec6d-8994-4886-9006-3cddbf0ac027
# ╟─0b4853ce-8dd0-4fe5-881d-cdf5d0289cc8
# ╠═1cbd900c-b74b-41c8-a57c-6a6ae644b88c
# ╟─0abb581e-6835-4dfc-bd01-5b2b6306d242
# ╠═895ee250-e3ac-4879-b8ee-f1ba8729adaa
# ╠═3c1157bd-9780-429a-9cbf-cd5a18d676da
# ╟─9ef928f0-5577-4507-8fbe-193ce6925716
# ╠═04189afa-6cdd-45f4-8eaa-12140bde3ada


text:
## --- Dealing with systematic uncertainty

    function add_systematic_uncert_UPb(agedistmyr::Array{Float64,1})
        λ238 = log(2)/44683 # Jaffey decay constant, 1/Myr
        λ238_ = log(2)./(44683 .+ 24 .* randn(size(agedistmyr))) # Jaffey decay with uncertainty 1/Myr
        σtracer = 1 .+ (0.03/2/100) .* randn(size(agedistmyr)) # Earthtime tracer uncertainty
        r = exp.(agedistmyr .* λ238) .- 1
        return log.(r .* σtracer .+ 1) ./ λ238_
    end
    function add_systematic_uncert_UPb(agedistmyr::Array{Float64,2})
        λ238 = log(2)/44683 # Jaffey decay constant, 1/Myr
        λ238_ = log(2)./(44683 .+ 24 .* randn(1, size(agedistmyr,2))) # Jaffey decay with uncertainty 1/Myr
        σtracer = 1 .+ (0.03/2/100) .* randn(1, size(agedistmyr,2)) # Earthtime tracer uncertainty
        r = exp.(agedistmyr .* λ238) .- 1
        return log.(r .* σtracer .+ 1) ./ λ238_
    end

    export add_systematic_uncert_UPb

    function add_systematic_uncert_ArAr(agedistmyr::Array{Float64,1})
        # Optimization intercalibration values from Renne 2011 reply
        # Values
        κ = 1.6417E-03
        λϵ = 0.5757E-10 # 1/years
        λβ = 4.9548E-10 # 1/years
        # Uncertainties
        σκ = 0.0045E-03
        σλϵ = 0.0017E-10 # 1/years
        σλβ = 0.0134E-10 # 1/years
        # Covariances
        σκσλϵ =  7.1903E-19
        σκσλβ = -6.5839E-19
        σλϵσλβ = -3.4711E-26

        # Create a multivariate normal distribution from the decay constant
        # means and covariance matrix
        μ = [κ, λϵ, λβ]
        Σ = [ σκ^2   σκσλϵ  σκσλβ;
              σκσλϵ  σλϵ^2 σλϵσλβ;
              σκσλβ σλϵσλβ  σλβ^2;]
        κλλ_ = rand(MvNormal(μ,Σ),length(agedistmyr)) # Draw from MvNormal
        # Exctract individual variables
        κ_ = κλλ_[1,:]
        λϵ_ = κλλ_[2,:]
        λβ_ = κλλ_[3,:]

        # Convert age distribution to R distribution. Assume Myr input
        λ = λβ + λϵ
        Rdist = (exp.(agedistmyr.*λ.*1000000) .- 1) .* λϵ./(κ.*λ)

        # Now convert back to age distribution
        λ_ = λβ_ .+ λϵ_
        return log.(Rdist.*κ_.*λ_./λϵ_ .+ 1)./(λ_.*1000000)
    end
    function add_systematic_uncert_ArAr(agedistmyr::Array{Float64,2})
        # Optimization intercalibration values from Renne 2011 reply
        # Values
        κ = 1.6417E-03
        λϵ = 0.5757E-10
        λβ = 4.9548E-10
        # Uncertainties
        σκ = 0.0045E-03
        σλϵ = 0.0017E-10
        σλβ = 0.0134E-10
        # Covariances
        σκσλϵ =  7.1903E-19
        σκσλβ = -6.5839E-19
        σλϵσλβ = -3.4711E-26

        # Create a multivariate normal distribution from the decay constant
        # means and covariance matrix
        μ = [κ, λϵ, λβ]
        Σ = [ σκ^2   σκσλϵ  σκσλβ;
              σκσλϵ  σλϵ^2 σλϵσλβ;
              σκσλβ σλϵσλβ  σλβ^2;]
        κλλ_ = rand(MvNormal(μ,Σ),size(agedistmyr,2)) # Draw from MvNormal
        # Exctract individual variables
        κ_ = κλλ_[1:1,:]
        λϵ_ = κλλ_[2:2,:]
        λβ_ = κλλ_[3:3,:]

        # Convert age distribution to R distribution
        λ = λβ + λϵ
        Rdist = (exp.(agedistmyr.*λ.*1000000) .- 1) .* λϵ./(κ.*λ)

        # Now convert back to age distribution
        λ_ = λβ_ .+ λϵ_
        return log.(Rdist.*κ_.*λ_./λϵ_ .+ 1)./(λ_.*1000000)
    end

    export add_systematic_uncert_ArAr

## --- End of File


text:
# These are copied from LinarAlgebra because they are not exported,
# and are perhaps not part of the API. This should be checked.
# Julia's API definition is semi-formal and not easily verified.

"""
    QMatrices.checksquare(A)

Check that a matrix is square, then return its common dimension.
For multiple arguments, return a vector.

# Examples
```jldoctest
julia> A = fill(1, (4,4)); B = fill(1, (5,5));

julia> LinearAlgebra.checksquare(A, B)
2-element Array{Int64,1}:
 4
 5
```
"""
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

function checksquare(A...)
    sizes = Int[]
    for a in A
        size(a,1)==size(a,2) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(a))"))
        push!(sizes, size(a,1))
    end
    return sizes
end


text:
using MultipleViewGeometry, Test
using StaticArrays, LinearAlgebra

# Tests for a set of two-dimensional Cartesian points represented by homogeneous
# coordinates.
ℳ = map(Point2D,
        [(-10.0, -10.0),
         (-10.0,  10.0),
         ( 10.0, -10.0),
         ( 10.0,  10.0)])

ℳʹ, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
@test ℳʹ == map(Point2D,
                                        [(-1.0,-1.0),
                                         (-1.0, 1.0),
                                         (1.0, -1.0),
                                         (1.0,  1.0)])

Λ =  [MMatrix{3,3}(Matrix(Diagonal([1.0,1.0,0.0]))) for i = 1:length(ℳ)]
Λʹ = transform(CovarianceMatrices(), CanonicalToHartley(), Λ , 𝐓)
for 𝚲 ∈ Λ
    @test 𝚲 == Matrix(Diagonal([1.0, 1.0, 0.0]))
end

for 𝚲 ∈ Λʹ
    @test 𝚲 ≈ Matrix(Diagonal([1.0, 1.0, 0.0])/100)
end

# Tests for a set of three-dimensional Cartesian points represented by homogeneous
# coordinates.
ℳ = map(Point3D,
           [(-10.0, -10.0, -10.0),
            (-10.0, -10.0,  10.0),
            (-10.0,  10.0, -10.0),
            (-10.0,  10.0,  10.0),
            ( 10.0, -10.0, -10.0),
            ( 10.0, -10.0,  10.0),
            ( 10.0,  10.0, -10.0),
            ( 10.0,  10.0,  10.0)])

ℳʹ, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
@test ℳʹ == map(Point3D,
                                         [(-1.0,-1.0, -1.0),
                                          (-1.0,-1.0,  1.0),
                                          (-1.0, 1.0, -1.0),
                                          (-1.0, 1.0,  1.0),
                                          (1.0, -1.0, -1.0),
                                          (1.0, -1.0,  1.0),
                                          (1.0,  1.0, -1.0),
                                          (1.0,  1.0,  1.0)])
@test 𝐓 == [0.1 0.0 0.0 -0.0;
            0.0 0.1 0.0 -0.0;
            0.0 0.0 0.1 -0.0;
            0.0 0.0 0.0 1.0]

Λ = [MMatrix{4,4}(Matrix(Diagonal([1.0, 1.0, 1.0, 0.0]))) for i = 1:length(ℳ)]
Λʹ = transform(CovarianceMatrices(), CanonicalToHartley(), Λ , 𝐓)
for 𝚲 ∈ Λ
    @test 𝚲 == Matrix(Diagonal([1.0, 1.0, 1.0, 0.0]))
end
for 𝚲 ∈ Λʹ
    @test 𝚲 ≈ Matrix(Diagonal([1.0, 1.0, 1.0, 0.0])/100)
end


text:
type Lattice{T <: Real, I <: Int} <: AbstractLattice
    celerities::Matrix{I}
    weights::Vector{T}
    inversion::Vector{I}

    # Creates inversion list automatically
    function Lattice(celerities::Matrix{I}, weights::Vector{T})
        @assert(size(celerities, 2) == length(weights))
        inversion = Int64[-1 for i in 1:length(weights)]
        for i in 1:length(weights)
            if inversion[i] != -1
                continue
            end
            for j in i:length(weights)
                if all(abs(celerities[:, i] + celerities[:, j]) .< 1e-8)
                   inversion[i], inversion[j] = j, i
                   break
                end
            end
        end
        @assert all(inversion .!= -1)

        new(celerities, weights, inversion)
    end
end

const D2Q9 = Lattice{Float64, Int64}(
  Int64[0 1 0 -1  0 1  1 -1 -1;
          0 0 1  0 -1 1 -1  1 -1],
  vcat(4./9., [1./9. for i=1:4], [1./36. for i=1:4])
)

const D3Q19 = Lattice{Float64, Int64}(
  Int64[0 1 0 0 -1  0  0 1 -1  1 -1 1 -1  1 -1 0  0  0  0
          0 0 1 0  0 -1  0 1 -1 -1  1 0  0  0  0 1 -1  1 -1
          0 0 0 1  0  0 -1 0  0  0  0 1 -1 -1  1 1 -1 -1  1],
  vcat(1./3., [1./18. for i=1:6], [1./36. for i=7:18])
)

const speed_of_sound_squared = 1.0 / 3.0
const speed_of_sound = sqrt(speed_of_sound_squared)

# Expand usage of neighbor_index
function neighbor_index(
        indexing::Indexing, site::GridCoords, lattice::AbstractLattice, direction::Integer)
    neighbor_index(indexing, site, lattice.celerities[:, direction])
end
function neighbor_index(
        indexing::Indexing, site::Integer, lattice::AbstractLattice, direction::Integer)
    neighbor_index(indexing, gridcoords(indexing, site), lattice.celerities[:, direction])
end
function neighbor_index(sim::Simulation, site::Union{Integer, GridCoords}, direction::Integer)
    neighbor_index(sim.indexing, site, sim.lattice, direction)
end

length(lattice::AbstractLattice) = length(lattice.weights)
ndims(lattice::AbstractLattice) = size(lattice.celerities, 1)


text:
using PyCall

@pyimport __builtin__
@pyimport Bio.SeqIO as SeqIO
@pyimport Bio.Seq as Seq

function extract_seqs(seq_IDs, filenames, format)
	output_array = []
	for filename in filenames
		print("Checking "*filename*"\n")
		for record in SeqIO.parse(filename, format)
			if record[:id] in seq_IDs
				push!(output_array, record)
			end
		end
	end
	return output_array
end

function extract_nt_seqs(models,project,filenames,output;format="fastq")
	if size(project["seq_info_table"]) == (0,0)
		error("This project contains no sequence information table. Rerun new_project() with seqinfo=true to be able to extract sequences")
	end

	seqs_to_retrieve = []
	for model in models
		append!(seqs_to_retrieve, project["seq_info_table"][project["seq_info_table"][:model] .== model, :seqID])
	end
	print(seqs_to_retrieve)
	output_sequences = extract_seqs(seqs_to_retrieve, filenames, format)
	output_handle = __builtin__.open(output, "w")
	SeqIO.write(output_sequences,output_handle,format)
end

function extract_aa_seqs(models,project,filenames,output;format="fasta")
	if size(project["seq_info_table"]) == (0,0)
		error("This project contains no sequence information table. Rerun new_project() with seqinfo=true to be able to extract sequences")
	end

	nt_seqs_to_retrieve = []
	translations = []
	aa_seqs_to_retrieve = []
	for model in models
		append!(nt_seqs_to_retrieve, project["seq_info_table"][project["seq_info_table"][:model] .== model, :seqID])
		append!(translations, project["seq_info_table"][project["seq_info_table"][:model] .== model, :translation])
		for (i, seqID) in enumerate(nt_seqs_to_retrieve)
			push!(aa_seqs_to_retrieve, join([seqID,translations[i]],'.'))
		end
	end
	output_sequences = extract_seqs(aa_seqs_to_retrieve, filenames, format)
	output_handle = __builtin__.open(output, "w")
	SeqIO.write(output_sequences,output_handle,format)
end

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
using MolecularGraph: bondsymbol!, bond!

@testset "smarts.bond" begin

@testset "bondsymbol" begin
    state = SmartsParser("", false)
    implicit1 = bondsymbol!(state)
    @test implicit1 === nothing

    state = SmartsParser("-", false)
    explicit1 = bondsymbol!(state)
    @test explicit1 == QueryFormula(:bondorder, 1)

    state = SmartsParser("\\?", false)
    stereo4 = bondsymbol!(state)
    @test stereo4 == QueryFormula(:not, QueryFormula(:stereo, :up))
    @test state.pos == 3
end

@testset "bond" begin
    state = SmilesParser("#", true)
    triple = bond!(state)
    @test triple.order == 3

    state = SmilesParser(":", true)
    arom = bond!(state)
    @test arom.isaromatic == true
end

@testset "smartsbond" begin
    state = SmartsParser("~", false)
    anyb = bond!(state)
    @test anyb.query == QueryFormula(:any, true)

    state = SmartsParser("-!@", false)
    notring = bond!(state)
    @test notring.query == QueryFormula(:and, Set([
        QueryFormula(:bondorder, 1),
        QueryFormula(:isringbond, false)
    ]))
end

end # smiles.bond


text:
# TODO: Make a parallel simulation file.
# - Call with
# julia  parsim.jl  <path-to-config>
# - Put all config.jl in sim-config/

using Distributed

# NOTE: These libs have to be imported in
#       main process and worker processes.
println("Loading packages. This may take a few minutes ..."); flush(stdout)
using Cytof5  # This is in dev mode
import Pkg; Pkg.activate("../../")  # sims
using Random, Distributions, BSON

# NOTE: Change this.
if length(ARGS) == 3
  SIMDIR = ARGS[1]
  RESULTS_DIR_PREFIX = ARGS[2]
  AWS_BUCKET_PREFIX = ARGS[3]
else
  println("Usage: julia parsim.jl <simdir> <results_dir_prefix> <aws_bucket_prefix>")
end

function setnumcores(n::Int)
  children_procs = filter(w -> w > 1, workers())
  rmprocs(children_procs)
  addprocs(n)
end

# NOTE: read configs
include("$(SIMDIR)/settings.jl")

for setting in settings
  setting[:results_dir] = "$(RESULTS_DIR_PREFIX)/$(setting[:outdir_suffix])"
  setting[:aws_bucket] = "$(AWS_BUCKET_PREFIX)/$(setting[:outdir_suffix])"
end

ncores = min(20, length(settings))  # use at most 20 cores on server
setnumcores(ncores)

println("Sourcing files ..."); flush(stdout)
@everywhere SIMDIR = $SIMDIR
@everywhere include("../Util/Util.jl")
@everywhere include("$(SIMDIR)/simfn.jl")

@everywhere function sim(setting)
  results_dir = setting[:results_dir]
  aws_bucket = setting[:aws_bucket]
  mkpath(results_dir)

  println("Running $(results_dir)")
  flush(stdout)

  Util.redirect_all("$(results_dir)/log") do
    simfn(setting)
  end

  # Send to S3.
  Util.s3sync(results_dir, aws_bucket, `--exclude '*.nfs'`)

  # Remove results to save space.
  # rm(results_dir, recursive=true)
  return
end

# NOTE:
# - f::Function: A function which takes one argument of type Dict{Any}
# - settings::Vector{Dict}): A vector of settings
println("Starting jobs ..."); flush(stdout)
result = pmap(sim, settings, on_error=identity);
println(result)
println("DONE with all runs!")


text:
function theme_black()
    Theme(
        backgroundcolor = :black,
        textcolor = :white,
        linecolor = :white,
        Axis = (
            backgroundcolor = :transparent,
            bottomspinecolor = :white,
            topspinecolor = :white,
            leftspinecolor = :white,
            rightspinecolor = :white,
            xgridcolor = RGBAf0(1, 1, 1, 0.16),
            ygridcolor = RGBAf0(1, 1, 1, 0.16),
            xtickcolor = :white,
            ytickcolor = :white,
        ),
        Legend = (
            framecolor = :white,
            bgcolor = :black,
        ),
        Axis3 = (
            xgridcolor = RGBAf0(1, 1, 1, 0.16),
            ygridcolor = RGBAf0(1, 1, 1, 0.16),
            zgridcolor = RGBAf0(1, 1, 1, 0.16),
            xspinecolor_1 = :white,
            yspinecolor_1 = :white,
            zspinecolor_1 = :white,
            xspinecolor_2 = :white,
            yspinecolor_2 = :white,
            zspinecolor_2 = :white,
            xspinecolor_3 = :white,
            yspinecolor_3 = :white,
            zspinecolor_3 = :white,
            xticklabelpad = 3,
            yticklabelpad = 3,
            zticklabelpad = 6,
            xtickcolor = :white,
            ytickcolor = :white,
            ztickcolor = :white,
        ),
        Colorbar = (
            tickcolor = :white,
            spinecolor = :white,
            topspinecolor = :white,
            bottomspinecolor = :white,
            leftspinecolor = :white,
            rightspinecolor = :white,
        )
    )
end


text:
using PkgDev
using Compat, Compat.Test, Compat.Pkg, Compat.Random, Compat.LibGit2

function temp_pkg_dir(fn::Function, remove_tmp_dir::Bool=true)
    # Used in tests below to set up and tear down a sandboxed package directory
    tmpdir = joinpath(tempdir(),Random.randstring())
    withenv("JULIA_PKGDIR" => tmpdir) do
        @test !isdir(Pkg.dir())
        try
            Pkg.init()
            @test isdir(Pkg.dir())
            Pkg.resolve()

            fn(Pkg.Dir.path())
        finally
            remove_tmp_dir && Base.rm(tmpdir, recursive=true)
        end
    end
end

temp_pkg_dir() do pkgdir

    @testset "testing a package with test dependencies causes them to be installed for the duration of the test" begin
        PkgDev.generate("PackageWithTestDependencies", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]
        @test read(Pkg.dir("PackageWithTestDependencies","REQUIRE"), String) == "julia $(PkgDev.Generate.versionfloor(VERSION))\n"

        isdir(Pkg.dir("PackageWithTestDependencies","test")) || mkdir(Pkg.dir("PackageWithTestDependencies","test"))
        open(Pkg.dir("PackageWithTestDependencies","test","REQUIRE"),"w") do f
            println(f,"Example")
        end

        open(Pkg.dir("PackageWithTestDependencies","test","runtests.jl"),"w") do f
            println(f, "using " * (VERSION < v"0.7.0-DEV.2005" ? "Base." : "") * "Test")
            if VERSION >= v"0.7.0-DEV.3656"
                println(f, "using Pkg")
            end
            println(f, "@test haskey(Pkg.installed(), \"Example\")")
        end

        Pkg.resolve()
        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

        Pkg.test("PackageWithTestDependencies")

        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

        # trying to pin an unregistered package errors
        try
            Pkg.pin("PackageWithTestDependencies", v"1.0.0")
            error("unexpected")
        catch err
            @test isa(err, Pkg.PkgError)
            @test err.msg == "PackageWithTestDependencies cannot be pinned – not a registered package"
        end
    end

    @testset "generating a package with .jl extension" begin
        PkgDev.generate("PackageWithExtension.jl", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        @show keys(Pkg.installed())
        @test "PackageWithExtension" in keys(Pkg.installed())
    end

    @testset "testing a package with no runtests.jl errors" begin
        PkgDev.generate("PackageWithNoTests", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        if isfile(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
            Base.rm(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
        end

        try
            Pkg.test("PackageWithNoTests")
            error("unexpected")
        catch err
            @test err.msg == "PackageWithNoTests did not provide a test/runtests.jl file"
        end
    end

    @testset "testing a package with failing tests errors" begin
        PkgDev.generate("PackageWithFailingTests", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        isdir(Pkg.dir("PackageWithFailingTests","test")) || mkdir(Pkg.dir("PackageWithFailingTests","test"))
        open(Pkg.dir("PackageWithFailingTests", "test", "runtests.jl"),"w") do f
            println(f, "using " * (VERSION < v"0.7.0-DEV.2005" ? "Base." : "") * "Test")
            println(f, "@test false")
        end

        try
            Pkg.test("PackageWithFailingTests")
            error("unexpected")
        catch err
            @test err.msg == "PackageWithFailingTests had test errors"
        end
    end

    # FIXME coverage is currently broken on windows?
    Compat.Sys.isunix() && @testset "testing with code-coverage" begin
        PkgDev.generate("PackageWithCodeCoverage", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        src = """
module PackageWithCodeCoverage

export f1, f2, f3, untested

f1(x) = 2x
f2(x) = f1(x)
function f3(x)
    3x
end
untested(x) = 7

end"""
        linetested = [false, false, false, false, true, true, false, true, false, false]
        open(Pkg.dir("PackageWithCodeCoverage", "src", "PackageWithCodeCoverage.jl"), "w") do f
            println(f, src)
        end
        isdir(Pkg.dir("PackageWithCodeCoverage","test")) || mkdir(Pkg.dir("PackageWithCodeCoverage","test"))
        open(Pkg.dir("PackageWithCodeCoverage", "test", "runtests.jl"),"w") do f
            println(f, "using PackageWithCodeCoverage")
            println(f, "using " * (VERSION < v"0.7.0-DEV.2005" ? "Base." : "") * "Test")
            println(f, "@test f2(2) == 4")
            println(f, "@test f3(5) == 15")
        end

        let codecov_yml = Pkg.dir("PackageWithCodeCoverage", ".codecov.yml")
            @test isfile(codecov_yml)
            @test readchomp(codecov_yml) == "comment: false"
        end

        Pkg.test("PackageWithCodeCoverage")
        covdir = Pkg.dir("PackageWithCodeCoverage","src")
        covfiles = filter!(x -> occursin("PackageWithCodeCoverage.jl", x) && occursin(".cov", x), readdir(covdir))
        @test isempty(covfiles)
        Pkg.test("PackageWithCodeCoverage", coverage=true)
        covfiles = filter!(x -> occursin("PackageWithCodeCoverage.jl", x) && occursin(".cov", x), readdir(covdir))
        @test !isempty(covfiles)
        for file in covfiles
            @test isfile(joinpath(covdir,file))
            covstr = read(joinpath(covdir,file), String)
            srclines = split(src, '\n')
            covlines = split(covstr, '\n')
            for i = 1:length(linetested)
                covline = (linetested[i] ? "        1 " : "        - ")*srclines[i]
                @test covlines[i] == covline
            end
        end
    end

    if haskey(ENV, "CI") && lowercase(ENV["CI"]) == "true"
        Compat.@info("setting git global configuration")
        run(`git config --global user.name "Julia Test"`)
        run(`git config --global user.email test@julialang.org`)
        run(`git config --global github.user JuliaTest`)
    end

    @testset "testing package tags" begin
        PkgDev.generate("PackageWithTags", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        PkgDev.register("PackageWithTags")
        PkgDev.tag("PackageWithTags")
        PkgDev.tag("PackageWithTags")

        # check tags
        repo = LibGit2.GitRepo(joinpath(pkgdir, "PackageWithTags"))
        meta_repo = LibGit2.GitRepo(joinpath(pkgdir, "METADATA"))
        try
            tags = LibGit2.tag_list(repo)
            @test ("v0.0.1" in tags)
            @test ("v0.0.2" in tags)
            # test that those version files exist
            @test ispath(joinpath(pkgdir, "METADATA", "PackageWithTags",
                                  "versions", "0.0.1", "sha1"))
            # check that we actually commited those files to METADATA
            # (this test always succeeds for some reason -- we should be using
            #  git_diff_index_to_workdir here)
            @test LibGit2.isdirty(meta_repo) == false
        finally
            finalize(repo)
            finalize(meta_repo)
        end
    end

    @testset "testing freeable" begin
        Pkg.add("Example")
        io = IOBuffer()
        f = PkgDev.freeable(io)
        @test !(any(f .== "Example") || occursin("Example", String(take!(io))))
        Pkg.checkout("Example")
        f = PkgDev.freeable(io)
        @test any(f .== "Example") || occursin("Example", String(take!(io)))
    end

    @testset "testing package registration" begin
        PkgDev.generate("GreatNewPackage", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        PkgDev.register("GreatNewPackage")
        @test !isempty(read(joinpath(pkgdir, "METADATA", "GreatNewPackage", "url"), String))
    end
end

@testset "Testing package utils" begin
    @test PkgDev.getrepohttpurl("https://github.com/JuliaLang/PkgDev.jl.git") == "https://github.com/JuliaLang/PkgDev.jl"
    @test PkgDev.getrepohttpurl("git://github.com/JuliaLang/PkgDev.jl.git")  == "https://github.com/JuliaLang/PkgDev.jl"
end


text:
# FFSCoStateInitializer : (F)orward (S)ingle (S)hooting CoState Initializer
# This struct consaints information and has associated methods for
# initializing the initial time co-state variables for solving
# indirect trajectory optimization problems using a forward single
# shooting based approach.
#
# Note: Currently, implementation assumes the trajectory optimization
# problem is for a single spacecraft with a 6-element state vector + mass,
# and therefore, 7 co-state variables. Other assumptions also made in cost
# function(s) that should likely be addressed

mutable struct FSSCoStateInitializer{HOT,HOOT} <: HeuristicsCoStateInitializer{HOT}
    # Heuristic optimizer
    ho::HOT

    # Heuristic optimizer options
    hoOpts::HOOT

    # Initialized co-state vector
    λh::Vector{Float64}
end

function FSSCoStateInitializer(bvpFunc, ICS, FCS;
                               costFunc = :WSS, optimizer = :PSO, numParticles = 100, numSwarms = 4,
                               initMethod = :Uniform,
                               UBs = [100, 100, 100, 50, 50, 50, 50],
                               LBs = [-100, -100, -100, -50, -50, -50, -50],
                               iUBs = nothing,
                               iLBs = nothing,
                               weights = nothing,
                               display = true,
                               displayInterval = 1,
                               maxIters = 1000,
                               funcTol = 1e-6,
                               maxStallIters = 25,
                               maxStallTime = 500,
                               maxTime = 1800,
                               useParallel = true
                               )

    # Check size of ICS and FCS
    if length(ICS) != 7 && length(FCS) != 7
        throw(ArgumentError("Initial and final boundary condition vectors must be of length 7."))
    end

    # Generate cost function and heuristic optimization problem
    if costFunc == :WSS # Weighted Sum of Squares
        if weights === nothing
            weights = [10, 10, 10, 1, 1, 1, 1]
        elseif length(weights) != 7
            throw(ArgumentError("Weights vector must be of length 7 to use weighted sum of squares cost function."))
        end
        prob = Problem(x -> cfFSSWSS(x, bvpFunc, ICS, FCS, weights), LBs, UBs)
    elseif costFunc == :WSSWM
        if weights === nothing
            weights = [10, 10, 10, 1, 1, 1, 1, 10]
        elseif length(weights) != 8
            throw(ArgumentError("Weights vector must be of length 8 to use weighted sum of squares w/ mass cost function."))
        end
        prob = Problem(x -> cfFSSWSSWM(x, bvpFunc, ICS, FCS, weights), LBs, UBs)
    else
        throw(ArgumentError("Cost function type not implemented."))
    end

    # Initialize optimizer and options
    opts = Options(;display = display, displayInterval = displayInterval,
                    maxIters = maxIters, useParallel = useParallel,
                    maxStallIters = maxStallIters, maxStallTime = maxStallTime,
                    maxTime = maxTime, iUB= iUBs, iLB = iLBs, funcTol = funcTol)

    if optimizer == :PSO
        ho   = PSO(prob; numParticles = numParticles, initMethod = initMethod)
    elseif optimizer == :MS_PSO
        ho   = MS_PSO(prob; numParticlesPerSwarm = numParticles, numSwarms = numSwarms, initMethod = initMethod)
    end

    # Instantiate FFS initializer
    FSSCoStateInitializer{typeof(ho), typeof(opts)}(ho, opts, Vector{Float64}(undef, 7))
end

# Cost functions
function cfFSSWSS(x, bvpFunc, ICS, FCS, ws)
    # Initialize state/co-state vector
    y0 = @SVector [ICS[1], ICS[2], ICS[3], ICS[4], ICS[5], ICS[6], ICS[7],
                   x[1], x[2], x[3], x[4], x[5], x[6], x[7]]

    # Evaluate bvp function
    yf, timeToFinalTime = bvpFunc(y0)

    # Compute cost
    cost = ws[1]*(yf[1] - FCS[1])^2 +
           ws[2]*(yf[2] - FCS[2])^2 +
           ws[3]*(yf[3] - FCS[3])^2 +
           ws[4]*(yf[4]  - FCS[4])^2 +
           ws[5]*(yf[5]  - FCS[5])^2 +
           ws[6]*(yf[6]  - FCS[6])^2 +
           ws[7]*(yf[14] - FCS[7])^2 +
           timeToFinalTime

    # Shift cost if terminated early
    return cost
end

function cfFSSWSSWM(x, bvpFunc, ICS, FCS, ws)
    # Initialize state/co-state vector
    y0 = @SVector [ICS[1], ICS[2], ICS[3], ICS[4], ICS[5], ICS[6], ICS[7],
                   x[1], x[2], x[3], x[4], x[5], x[6], x[7]]

    # Evaluate bvp function
    yf, timeToFinalTime = bvpFunc(y0)

    # Compute cost
    cost = ws[1]*(yf[1]  - FCS[1])^2 +
           ws[2]*(yf[2]  - FCS[2])^2 +
           ws[3]*(yf[3]  - FCS[3])^2 +
           ws[4]*(yf[4]  - FCS[4])^2 +
           ws[5]*(yf[5]  - FCS[5])^2 +
           ws[6]*(yf[6]  - FCS[6])^2 +
           ws[7]*(yf[14] - FCS[7])^2 +
           ws[8]*(yf[7]  - y0[7])^2 +
           timeToFinalTime

    # Shift cost if terminated early
    return cost
end


text:
listMyRepos(c::Client) = getParsedResponse(Vector{Repository},c,Requests.get,"/user/repos")

listUserRepos(c::Client, user::String) = getParsedResponse(Vector{Repository},c,Requests.get,"/user/$(user)/repos")

listOrgRepos(c::Client, org::String) = getParsedResponse(Vector{Repository},c,Requests.get,"/orgs/$(org)/repos")

FieldTags.@tag immutable CreateRepoOption
	name::String
	description::Nullable{String}
	private::Nullable{Bool}
	autoInit::Nullable{Bool} => json:"auto_init"
	gitignores::Nullable{String}
	license::Nullable{String}
	readme::Nullable{String}
end

function CreateRepoOption(name; kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[name]

	for field in fieldnames(CreateRepoOption)[2:end]
		push!(args,get(kwargsDict,field, fieldtype(CreateRepoOption,field)()))
	end

	CreateRepoOption(args...)
end

createRepo(c::Client,opt::CreateRepoOption) = getParsedResponse(Repository,c,Requests.post,"/user/repos"; json = marshalJSON(opt))
createOrgRepo(c::Client, org::String, opt::CreateRepoOption) = getParsedResponse(Repository,c,Requests.post,"/org/$(org)/repos"; json = marshalJSON(opt))

getRepo(c::Client, owner::String, repo::String) = getParsedResponse(Repository,c,Requests.get,"/repos/$owner/$repo")

getRepo(c::Client, owner::User, repo::String) = getRepo(c,owner.userName,repo)

getRepo(c::Client, index::Int) = getParsedResponse(Repository,c,Requests.get,"/repositories/$(index)")

deleteRepo(c::Client, owner::String, repo::String) = (getResponse(c,Requests.delete,"/repos/$owner/$repo"); nothing)

deleteRepo(c::Client, owner::User, repo::String) = deleteRepo(c,owner.userName,repo)

deleteRepo(c::Client, repo::Repository) = deleteRepo(c,repo.owner,repo.name)

FieldTags.@tag immutable MigrateRepoOption
	cloneAddr::String => json:"clone_addr"
	uid::Int
	repoName::String
	authUsername::Nullable{String} => json:"auth_username"
	authPassword::Nullable{String} => json:"auth_password"
	mirror::Nullable{Bool}
	private::Nullable{Bool}
	description::Nullable{String}
end
 function MigrateRepoOption(cloneAddr,uid,name; kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[cloneAddr,uid,name]

	for field in fieldnames(MigrateRepoOption)[4:end]
		push!(args,get(kwargsDict,field, fieldtype(MigrateRepoOption,field)()))
	end

	MigrateRepoOption(args...)
end

migrateRepo(c::Client, opt::MigrateRepoOption) = getParsedResponse(Repository,c,Requests.post,"/repos/migrate"; json = marshalJSON(opt))

################
### Branches ###
################

listRepoBranches(c::Client, user::String, repo::String) = getParsedResponse(Vector{Branch},c,Requests.get,"/repos/$(user)/$(repo)/branches")

listRepoBranches(c::Client, repo::Repository) = listRepoBranches(c,repo.owner.userName,repo.name)

getRepoBranch(c::Client, user::String, repo::String, branch::String) = getParsedResponse(Branch,c,Requests.get,"/repos/$(user)/$(repo)/branches/$(branch)")

getRepoBranch(c::Client, repo::Repository, branch::String) = listRepoBranches(c,repo.owner.userName,repo.name,branch)

####################
### Collaborator ###
####################

listCollaborators(c::Client, user::String, repo::String) = getParsedResponse(Vector{User},c,Requests.get,"/repos/$(user)/$(repo)/collaborators")

function isCollaborator(c::Client, user::String, repo::String, collaborator::String)
	code = getStatusCode(c,Requests.get,"/repos/$(user)/$(repo)/collaborators/$(collaborator)")
	return code == 204
end

function addCollaborator(c::Client, user::String, repo::String, collaborator::String, permission::String)
	data = Dict("permission" => permission)
	getResponse(c,Requests.post,"/repos/$(user)/$(repo)/collaborators/$(collaborator)";json = data)
	return nothing
end

deleteCollaborator(c::Client, user::String, repo::String, collaborator::String) = (getResponse(Repository,c,Requests.delete,"/repos/$(user)/$(repo)/collaborators/$(collaborator)"); nothing)

############
### File ###
############

function getFile(c::Client, user::String, repo::String, ref::String, tree::String)
	f = getResponse(c,Requests.get,"/repos/$(user)/$(repo)/raw/$(ref)/$(tree)")
	IOBuffer(Requests.bytes(f))
end

getFile(c::Client, repo::Repository, ref::String, tree::String) = listRepoBranches(c,repo.owner.userName,repo.name,ref,tree)

##################
### Deploy Key ###
##################

listDeployKeys(c::Client, user::String, repo::String) = getParsedResponse(Vector{DeployKey},c,Requests.get,"/repos/$(user)/$(repo)/keys")

getDeployKey(c::Client, user::String, repo::String, keyID::Int64) = getParsedResponse(DeployKey,c,Requests.get,"/repos/$(user)/$(repo)/keys/$(keyID)")

function createDeployKey(c::Client, user::String, repo::String, title::String, key::String)
	data = Dict("title" => title, "key" => key)
	getParsedResponse(DeployKey,c,Requests.put,"/repos/$(user)/$(repo)/keys"; json = data)
end

deleteDeployKey(c::Client, user::String, repo::String, keyID::Int64) = (getResponse(c,Requests.delete,"/repos/$(user)/$(repo)/keys/$(keyID)"); nothing)


text:
using SimpleBirthDeathProcess;
using DelimitedFiles;
using Plots;

N = 200;
F = Float64;
I = 200;
J = 100;
T = 1;
L = F(10);
M = F(10);

# values obtained with Maple 2018.2
correct_values = readdlm("log_trans_prob_values_bivariate.csv", ',', BigFloat);
correct_values = reshape(correct_values, N, :);

λset = [L * F(n) / N for n = 1:N];
μset = [M * F(n) / N for n = 1:N];
logp = map((λ, μ) -> trans_prob(I, J, T, [λ; μ]),
           repeat(reshape(λset, 1, :), N, 1),
           repeat(μset, 1, N));

relative_error = abs.(1 .- logp ./ correct_values);

colgr = ColorGradient([colorant"#ffffff", colorant"#525252"], [0.0, 1.0]);
pl = plot(λset, μset, log10.(relative_error), seriestype=:contour,
          levels=[-20; -16; -14; -13; -12], fill=true, contour_labels=false,
          legend=true, linecolour=:black, aspect_ratio=:equal, fillcolour=colgr,
          xlim=(0, 10), ylim=(0, 10), xlab="\\lambda", ylab="\\mu");

savefig(pl, "Figure_4.svg")


text:
# This file provides useful tools for meta-programming (e.g. code generation)

##########################################################################
#
#   code construction tools
#
##########################################################################

function code_block(stmts...)
    stmts_ = Any[]
    for s in stmts
        if !(s == nothing)
            push!(stmts_, s)
        end
    end
    isempty(stmts_) ? nothing : Expr(:block, stmts_...)
end

function flatten_code_block_to(q, stmt_lst)
    for s in stmt_lst
        if s == nothing || s.head == (:line)
            continue
        end
        if s.head == (:block)
            flatten_code_block_to(q, s.args)
        else
            push!(q, s)
        end
    end
end

function flatten_code_block(stmts...)
    q = Any[]
    flatten_code_block_to(q, stmts)
    code_block(q...)
end

function fun_call(funsym, args...)
    Expr(:call, funsym, args...)
end

function qname(m::Symbol, x::Symbol)
    Expr(:(.), m, Meta.quot(x))
end

qname(x::Symbol) = qname(:Devectorize, x)

function assignment(lhs, rhs)
    Expr(:(=), lhs, rhs)
end

function if_statement(con, tblock)
    Expr(:if, con, tblock)
end

function if_statement(con, tblock, fblock)
    Expr(:if, con, tblock, fblock)
end

function for_statement(head, body)
    Expr(:for, head, body)
end

function for_statement(i::Symbol, si::TIndex, ei::TIndex, body)
    Expr(:for, :( ($i) = ($si) : ($ei) ), body)
end

##########################################################################
#
#   shape inference
#
##########################################################################

# helper for returning length

to_length(s::Int) = s
to_length(s::@compat Tuple{Int}) = s[1]
to_length(s::@compat Tuple{Int, Int}) = s[1] * s[2]

# helper for returning 2D size

to_size2d(s::Int) = (s, 1)
to_size2d(s::@compat Tuple{Int}) = (s[1], 1)
to_size2d(s::@compat Tuple{Int, Int}) = s

# length inference (from LHS)

length_inference(s::Symbol) = fun_call(:length, s)

length_inference(s::Symbol, rgn::TColon) = fun_call(:length, s)

length_inference_(base_len::Expr, rgn::TColon) = base_len

function length_inference_(base_len::Expr, rgn::TInterval)
    first = rgn.first
    last = rgn.last

    blen = last == nothing ? base_len : last

    if isa(first, Int)
        if first == 1
            blen
        else
            b = first - 1
            if isa(blen, Int)
                max(blen - b, 0)
            else
                :( max($blen - $b, 0) )
            end
        end
    else
        if isa(blen, Int)
            a = blen + 1
            :( max($a - $first, 0) )
        else
            :( max($blen + 1 - $first, 0) )
        end
    end
end

function length_inference(s::Symbol, rgn::TInterval)
    length_inference_(fun_call(:length, s), rgn)
end

# the hyper-function to generate codes for size inference

function vec_size_inference(s::Symbol, rgn::Any)
    :( ($(length_inference(rgn)),) )
end

function vec_size_inference(s::Symbol, rgn::TColon)
    :( (length($s),) )
end

function vec_size_inference(s::Symbol, rgn::TInterval)
    :( ($(length_inference(s, rgn)),) )
end

function col_size_inference(s::Symbol, rgn::TColon)
    :( (size($s, 1),) )
end

function col_size_inference(s::Symbol, rgn::TInterval)
    baselen_ = :( size($s, 1) )
    :( ($(length_inference_(baselen_, rgn)),) )
end

function row_size_inference(s::Symbol, rgn::TColon)
    :( (1, size($s, 2)) )
end

function row_size_inference(s::Symbol, rgn::TInterval)
    baselen_ = :( size($s, 2) )
    :( (1, $(length_inference_(baselen_, rgn))) )
end

function mat_size_inference(s::Symbol, rrgn::TColon, crgn::TColon)
    fun_call(qname(:to_size2d), fun_call(:size, s))
end

function mat_size_inference(s::Symbol, rrgn::TColon, crgn::TInterval)
    bn = :( size($s, 2) )
    :( (size($s, 1), $(length_inference_(bn, crgn))) )
end

function mat_size_inference(s::Symbol, rrgn::TInterval, crgn::TColon)
    bm = :( size($s, 1) )
    :( ($(length_inference_(bm, rrgn)), size($s, 2)) )
end

function mat_size_inference(s::Symbol, rrgn::TInterval, crgn::TInterval)
    bm = :( size($s, 1) )
    bn = :( size($s, 2) )
    :( ($(length_inference_(bm, rrgn)), $(length_inference_(bn, crgn))) )
end

# ewise_shape

ewise_shape(s) = s

ewise_shape(::(@compat Tuple{}), ::@compat Tuple{}) = ()
ewise_shape(::(@compat Tuple{}), s) = s
ewise_shape(s, ::@compat Tuple{}) = s
ewise_shape(s1, s2) = promote_shape(s1, s2)

ewise_shape(::(@compat Tuple{}), ::(@compat Tuple{}), ::@compat Tuple{}) = ()
ewise_shape(s1, ::(@compat Tuple{}), ::@compat Tuple{}) = s1
ewise_shape(::(@compat Tuple{}), s2, ::@compat Tuple{}) = s2
ewise_shape(::(@compat Tuple{}), ::(@compat Tuple{}), s3) = s3
ewise_shape(s1, s2, ::@compat Tuple{}) = promote_shape(s1, s2)
ewise_shape(s1, ::(@compat Tuple{}), s3) = promote_shape(s1, s3)
ewise_shape(::(@compat Tuple{}), s2, s3) = promote_shape(s2, s3)
ewise_shape(s1, s2, s3) = promote_shape(promote_shape(s1, s2), s3)

ewise_shape(s1, s2, s3, s4...) = promote_shape(ewise_shape(s1, s2), ewise_shape(s3, s4...))

function ewise_size_inference(a1_siz::Expr)
    a1_siz
end

function ewise_size_inference(a1_siz::Expr, other_sizs::Expr...)
    fun_call(qname(:ewise_shape), a1_siz, other_sizs...)
end

##########################################################################
#
#   type inference
#
##########################################################################

type_inference(ex::TNum) = :( typeof($(ex.val)) )
type_inference(ex::TVar) = :( eltype($(ex.name)) )
type_inference(ex::TScalarVar) = :( typeof($(ex.name)) )
type_inference(ex::TRef) = :( eltype($(ex.host)) )
type_inference(ex::TAssign) = :( $(type_inference(ex.rhs)) )

function result_type_inference(f::Symbol, arg_tys::Expr...)
    tf = TFun{f}()
    fun_call(qname(:result_type), tf, arg_tys...)
end

##########################################################################
#
#   indexer
#
##########################################################################

indexer(rgn::TColon, i::TIndex) = i

function indexer(rgn::TInterval, i::TIndex)
    f = rgn.first
    if isa(f, Int)
        b = f - 1
        if b == 0
            i
        else
            :($i + $b)
        end
    else
        :($i + ($f - 1))
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, SHA-256

clear = "Rosetta code"
standard = "764faf5c61ac315f1497f9dfa542713965b785e5cc2f707d6468d7d1124cdfcf"

using SHA

crypt = sha256(clear)

println("Testing Julia's SHA-256:")
if crypt == standard
    println("  OK, \"", clear, "\" => ", crypt)
else
    println("The hash does not match the standard value.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
export KroneckerOperator

##########
# KroneckerOperator gives the kronecker product of two 1D operators
#########

immutable KroneckerOperator{S,V,DS,RS,DI,RI,T} <: Operator{T}
    ops::Tuple{S,V}
    domainspace::DS
    rangespace::RS
    domaintensorizer::DI
    rangetensorizer::RI
end

KroneckerOperator(A,B,ds::Space,rs::Space,di,ri) =
    KroneckerOperator{typeof(A),typeof(B),typeof(ds),typeof(rs),typeof(di),typeof(ri),
                        promote_type(eltype(A),eltype(B))}((A,B),ds,rs,di,ri)

KroneckerOperator(A,B,ds::Space,rs::Space) = KroneckerOperator(A,B,ds,rs,
                    CachedIterator(tensorizer(ds)),CachedIterator(tensorizer(rs)))
function KroneckerOperator(A,B)
    ds=domainspace(A)⊗domainspace(B)
    rs=rangespace(A)⊗rangespace(B)
    KroneckerOperator(A,B,ds,rs)
end
KroneckerOperator(A::UniformScaling,B::UniformScaling) =
    KroneckerOperator(ConstantOperator(A.λ),ConstantOperator(B.λ))
KroneckerOperator(A,B::UniformScaling) = KroneckerOperator(A,ConstantOperator(B.λ))
KroneckerOperator(A::UniformScaling,B) = KroneckerOperator(ConstantOperator(A.λ),B)
KroneckerOperator(A::Fun,B::Fun) = KroneckerOperator(Multiplication(A),Multiplication(B))
KroneckerOperator(A::UniformScaling,B::Fun) = KroneckerOperator(ConstantOperator(A.λ),Multiplication(B))
KroneckerOperator(A::Fun,B::UniformScaling) = KroneckerOperator(Multiplication(A),ConstantOperator(B.λ))
KroneckerOperator(A,B::Fun) = KroneckerOperator(A,Multiplication(B))
KroneckerOperator(A::Fun,B) = KroneckerOperator(Multiplication(A),B)

Base.eye{T,D}(S::Space{T,D,2}) = KroneckerOperator(IdentityOperator(),IdentityOperator(),S,S)

function promotedomainspace(K::KroneckerOperator,ds::TensorSpace)
    A=promotedomainspace(K.ops[1],ds[1])
    B=promotedomainspace(K.ops[2],ds[2])
    KroneckerOperator(A,B,ds,rangespace(A)⊗rangespace(B))
end

function promoterangespace(K::KroneckerOperator,rs::TensorSpace)
    A=promoterangespace(K.ops[1],rs[1])
    B=promoterangespace(K.ops[2],rs[2])
    KroneckerOperator(A,B,domainspace(K),rs)
end

function Base.convert{T<:Number}(::Type{Operator{T}},K::KroneckerOperator)
    if T == eltype(K)
        K
    else
        ops=Operator{T}(K.ops[1]),Operator{T}(K.ops[2])
        KroneckerOperator{typeof(ops[1]),typeof(ops[2]),typeof(K.domainspace),typeof(K.rangespace),
                            typeof(K.domaintensorizer),typeof(K.rangetensorizer),T}(ops,
              K.domainspace,K.rangespace,
              K.domaintensorizer,K.rangetensorizer)
    end
end

function colstart(A::KroneckerOperator,k::Integer)
    K=block(A.domaintensorizer,k)
    blockstart(A.rangetensorizer,max(1,K-blockbandwidth(A,2)))
end

function colstop(A::KroneckerOperator,k::Integer)
    K=block(A.domaintensorizer,k)
    st=blockstop(A.rangetensorizer,K+blockbandwidth(A,1))
    # zero indicates above dimension
    st==0?size(A,1):min(size(A,1),st)
end

function rowstart(A::KroneckerOperator,k::Integer)
    K=block(rangespace(A),k)
    blockstart(domainspace(A),max(1,K-blockbandwidth(A,1)))
end

function rowstop(A::KroneckerOperator,k::Integer)
    K=block(rangespace(A),k)
    st=blockstop(domainspace(A),K+blockbandwidth(A,2))
    # zero indicates above dimension
    st==0?size(A,2):min(size(A,2),st)
end

bandinds(K::KroneckerOperator) = all(isdiag,K.ops) ? (0,0) : (-∞,∞)

isbandedblockbanded(K::KroneckerOperator) = all(isbanded,K.ops)

blockbandinds(K::KroneckerOperator) =
    bandinds(K.ops[1],1)+bandinds(K.ops[2],1),bandinds(K.ops[1],2)+bandinds(K.ops[2],2)
blockbandinds(K::Operator,k::Integer) = blockbandinds(K)[k]
blockbandwidth(K::Operator,k::Integer) = k==1?-blockbandinds(K,k):blockbandinds(K,k)
subblockbandinds(K::KroneckerOperator,k::Integer) =
    k==1?min(bandinds(K.ops[1],1),-bandinds(K.ops[2],2)):max(bandinds(K.ops[1],2),-bandinds(K.ops[2],1))
subblockbandinds(::Union{ConstantOperator,ZeroOperator},::Integer) = 0

typealias Wrappers Union{ConversionWrapper,MultiplicationWrapper,DerivativeWrapper,LaplacianWrapper,
                       SpaceOperator,ConstantTimesOperator}

isbandedblockbanded(P::Union{PlusOperator,TimesOperator}) = all(isbandedblockbanded,P.ops)

blockbandinds(P::PlusOperator,k::Int) =
    mapreduce(op->blockbandinds(op,k),k==1?min:max,P.ops)
blockbandinds(P::PlusOperator) = blockbandinds(P,1),blockbandinds(P,2)
subblockbandinds(K::PlusOperator,k::Integer) =
    mapreduce(v->subblockbandinds(v,k),k==1?min:max,K.ops)

blockbandinds(P::TimesOperator,k::Int) = mapreduce(op->blockbandinds(op,k),+,P.ops)
subblockbandinds(P::TimesOperator,k::Int) = mapreduce(op->subblockbandinds(op,k),+,P.ops)
blockbandinds(P::TimesOperator) = blockbandinds(P,1),blockbandinds(P,2)

domaintensorizer(P::PlusOperator) = domaintensorizer(P.ops[1])
rangetensorizer(P::PlusOperator) = rangetensorizer(P.ops[1])

domaintensorizer(P::TimesOperator) = domaintensorizer(P.ops[end])
rangetensorizer(P::TimesOperator) = rangetensorizer(P.ops[1])

subblockbandinds(K::Wrappers,k::Integer) = subblockbandinds(K.op,k)
for FUNC in (:blockbandinds,:isbandedblockbanded,:domaintensorizer,:rangetensorizer)
    @eval $FUNC(K::Wrappers) = $FUNC(K.op)
end

function subblockbandindssum(P,k)
    ret=0
    for op in P
        ret+=subblockbandinds(op,k)::Int
    end
    ret
end

subblockbandinds(P::TimesOperator,k)=subblockbandindssum(P.ops,1)

subblockbandinds(K::Operator)=subblockbandinds(K,1),subblockbandinds(K,2)

domainspace(K::KroneckerOperator) = K.domainspace
rangespace(K::KroneckerOperator) = K.rangespace

domaintensorizer(K::KroneckerOperator) = K.domaintensorizer
rangetensorizer(K::KroneckerOperator) = K.rangetensorizer

# we suport 4-indexing with KroneckerOperator
# If A is K x J and B is N x M, then w
# index to match KO=reshape(kron(B,A),N,K,M,J)
# that is
# KO[k,n,j,m] = A[k,j]*B[n,m]
# TODO: arbitrary number of ops

getindex(KO::KroneckerOperator,k::Integer,n::Integer,j::Integer,m::Integer) =
    KO.ops[1][k,j]*KO.ops[2][n,m]

function getindex(KO::KroneckerOperator,kin::Integer,jin::Integer)
    j,m=KO.domaintensorizer[jin]
    k,n=KO.rangetensorizer[kin]
    KO[k,n,j,m]
end

function getindex(KO::KroneckerOperator,k::Integer)
    if size(KO,1) == 1
        KO[1,k]
    elseif size(KO,2) == 1
        KO[k,1]
    else
        throw(ArgumentError("[k] only defined for 1 x ∞ and ∞ x 1 operators"))
    end
end

*(A::KroneckerOperator,B::KroneckerOperator) =
    KroneckerOperator(A.ops[1]*B.ops[1],A.ops[2]*B.ops[2])

## Shorthand

⊗(A,B) = kron(A,B)

Base.kron(A::Operator,B::Operator) = KroneckerOperator(A,B)
Base.kron(A::Operator,B) = KroneckerOperator(A,B)
Base.kron(A,B::Operator) = KroneckerOperator(A,B)
Base.kron{T<:Operator}(A::Vector{T},B::Operator) =
    Operator{promote_type(eltype(T),eltype(B))}[kron(a,B) for a in A]
Base.kron{T<:Operator}(A::Operator,B::Vector{T}) =
    Operator{promote_type(eltype(T),eltype(A))}[kron(A,b) for b in B]
Base.kron{T<:Operator}(A::Vector{T},B::UniformScaling) =
    Operator{promote_type(eltype(T),eltype(B))}[kron(a,1.0B) for a in A]
Base.kron{T<:Operator}(A::UniformScaling,B::Vector{T}) =
    Operator{promote_type(eltype(T),eltype(A))}[kron(1.0A,b) for b in B]

## Conversion

conversion_rule(a::TensorSpace,b::TensorSpace) = conversion_type(a[1],b[1])⊗conversion_type(a[2],b[2])
maxspace(a::TensorSpace,b::TensorSpace) = maxspace(a[1],b[1])⊗maxspace(a[2],b[2])

# TODO: we explicetly state type to avoid type inference bug in 0.4

ConcreteConversion(a::BivariateSpace,b::BivariateSpace) =
    ConcreteConversion{typeof(a),typeof(b),
                        promote_type(eltype(a),eltype(b),real(eltype(domain(a))),real(eltype(domain(b))))}(a,b)

Conversion(a::TensorSpace,b::TensorSpace) = ConversionWrapper(promote_type(eltype(a),eltype(b)),
                KroneckerOperator(Conversion(a[1],b[1]),Conversion(a[2],b[2])))

Multiplication{D,T}(f::Fun{D,T},sp::BivariateSpace) =
    ConcreteMultiplication{D,typeof(sp),T,T}(chop(f,maxabs(f.coefficients)*40*eps(eltype(f))),sp)
function Multiplication{CS<:ConstantSpace,T,V}(f::Fun{TensorSpace{Tuple{CS,V},T,2}},sp::TensorSpace)
    a=Fun(f.coefficients,space(f)[2]) # coefficients are the same
    #Hack to avoid auto-typing bug.  TODO: incorporate basis
    MultiplicationWrapper(eltype(f),f,eye(sp[1])⊗Multiplication(a,sp[2]))
end
function Multiplication{CS<:ConstantSpace,T,V}(f::Fun{TensorSpace{Tuple{V,CS},T,2}},sp::TensorSpace)
    if isempty(f.coefficients)
        a=Fun(zeros(eltype(f),1),space(f)[1])
    else
        a=Fun(f.coefficients,space(f)[1])
    end
    MultiplicationWrapper(eltype(f),f,Multiplication(a,sp[1])⊗eye(sp[2]))
end

Multiplication{D<:UnivariateSpace,SV,TT,T}(f::Fun{D,T},sp::SumSpace{SV,TT,AnyDomain,2})=Multiplication(f⊗1,sp)
Multiplication{D<:UnivariateSpace,T}(f::Fun{D,T},sp::BivariateSpace)=Multiplication(f⊗1,sp)

## PDE Factorization

isfiniterange(::,k) = false
isfiniterange(B::KroneckerOperator,k::Integer) = isfinite(size(B.ops[k],1))
isfiniterange(B::PlusOperator,k::Integer) = isfiniterange(first(B.ops),k)

function findfunctionals(A::Vector,k::Integer)
    T=eltype(eltype(eltype(A)))
    indsBx=find(f->isfiniterange(f,k),A)
    if k==1
        indsBx,Operator{T}[(@assert dekron(Ai,2)==ConstantOperator(Float64,1.0); dekron(Ai,1)) for Ai in A[indsBx]]
    else
        @assert k==2
        indsBx,Operator{T}[(@assert dekron(Ai,1)==ConstantOperator(Float64,1.0); dekron(Ai,2)) for Ai in A[indsBx]]
    end
end

## transpose

Base.transpose(K::KroneckerOperator)=KroneckerOperator(K.ops[2],K.ops[1])

for TYP in (:ConversionWrapper,:MultiplicationWrapper,:DerivativeWrapper,:IntegralWrapper,:LaplacianWrapper),
    FUNC in (:domaintensorizer,:rangetensorizer,:blockbandinds,:subblockbandinds)
    @eval $FUNC(S::$TYP) = $FUNC(S.op)
end

Base.transpose(S::SpaceOperator) =
    SpaceOperator(transpose(S.op),domainspace(S).',rangespace(S).')
Base.transpose(S::ConstantTimesOperator) = sp.c*S.op.'

### Calculus

#TODO: general dimension
function Derivative{SV,TT}(S::TensorSpace{SV,TT,2},order::Vector{Int})
    @assert length(order)==2
    if order[1]==0
        Dy=Derivative(S[2],order[2])
        K=eye(S[1])⊗Dy
        T=eltype(Dy)
    elseif order[2]==0
        Dx=Derivative(S[1],order[1])
        K=Dx⊗eye(S[2])
        T=eltype(Dx)
    else
        Dx=Derivative(S[1],order[1])
        Dy=Derivative(S[2],order[2])
        K=Dx⊗Dy
        T=promote_type(eltype(Dx),eltype(Dy))
    end
    # try to work around type inference
    DerivativeWrapper{typeof(K),typeof(domainspace(K)),Vector{Int},T}(K,order)
end

## Functionals
Evaluation(sp::TensorSpace,x::Vec) = EvaluationWrapper(sp,x,zeros(Int,length(x)),⊗(map(Evaluation,sp.spaces,x)...))
Evaluation(sp::TensorSpace,x::Tuple) = Evaluation(sp,Vec(x...))

### Copy

# finds block lengths for a subrange
function blocklengthrange(rt,kr)
    KR=block(rt,first(kr)):block(rt,last(kr))
    Klengths=Array(Int,length(KR))
    for ν in eachindex(KR)
        Klengths[ν]=blocklength(rt,KR[ν])
    end
    Klengths[1]+=blockstart(rt,KR[1])-kr[1]
    Klengths[end]+=kr[end]-blockstop(rt,KR[end])
    Klengths
end

function Base.convert(::Type{BandedBlockBandedMatrix},S::SubOperator)
    kr,jr=parentindexes(S)
    KO=parent(S)
    l,u=blockbandinds(KO)
    λ,μ=subblockbandinds(KO)

    rt=rangetensorizer(KO)
    dt=domaintensorizer(KO)
    ret=bbbzeros(S)

    for J=1:blocksize(ret,2)
        jshft=jr[1]+blockstart(dt,J)-2
        for K=blockcolrange(ret,J)
            Bs=viewblock(ret,K,J)
            kshft=kr[1]+blockstart(rt,K)-2
            for j=1:size(Bs,2),k=colrange(Bs,j)
                Bs[k,j]=KO[k+kshft,j+jshft]
            end
        end
    end

    ret
end

#TODO: there's a bug here
function Base.convert{KKO<:KroneckerOperator,T}(::Type{BandedBlockBandedMatrix},S::SubOperator{T,KKO})
    kr,jr=parentindexes(S)
    KO=parent(S)
    l,u=blockbandinds(KO)
    λ,μ=subblockbandinds(KO)

    rt=rangetensorizer(KO)
    dt=domaintensorizer(KO)
    ret=bbbzeros(eltype(KO),-l,u,-λ,μ,
                blocklengthrange(rt,kr),
                blocklengthrange(dt,jr))

    A,B=KO.ops
    K=block(rt,kr[end]);J=block(dt,jr[end])
    AA=A[1:K,1:J]
    BB=B[1:K,1:J]

    for J=1:blocksize(ret,2)
        jshft=jr[1]+blockstart(dt,J)-2
        for K=blockcolrange(ret,J)
            Bs=viewblock(ret,K,J)
            kshft=kr[1]+blockstart(rt,K)-2
            for j=1:size(Bs,2),k=colrange(Bs,j)
                κ,ν=subblock2tensor(rt,K,k)
                ξ,μ=subblock2tensor(dt,J,j)
                Bs[k,j]=AA[κ,ξ]*BB[ν,μ]
            end
        end
    end

    ret
end


text:
using Test

using Distributions

@test convert(Binomial, Bernoulli(0.75)) == Binomial(1, 0.75)

@test convert(Gamma, Exponential(3.0)) == Gamma(1.0, 3.0)
@test convert(Gamma, Erlang(5, 2.0)) == Gamma(5.0, 2.0)
@test convert(GeneralizedGamma,  Weibull(2,3)) == GeneralizedGamma(3.0, 2.0, 2.0)


text:
module GLUE
using DataDeps

using ..Datasets: Dataset, get_channels
import ..Datasets: testfile, devfile, trainfile, get_labels

export CoLA, SST, MRPC, QQP, STS, MNLI, SNLI, QNLI, RTE, WNLI, Diagnostic

function __init__()
    cola_init()
    sst_init()
    mrpc_init()
    qqp_init()
    sts_init()
    mnli_init()
    snli_init()
    qnli_init()
    rte_init()
    wnli_init()
    diagnostic_init()
end

include("./cola.jl")
include("./sst.jl")
include("./mrpc.jl")
include("./qqp.jl")
include("./sts.jl")
include("./mnli.jl")
include("./snli.jl")
include("./qnli.jl")
include("./rte.jl")
include("./wnli.jl")
include("./diagnostic.jl")

end


text:
############################################################
# joAbstractSAparallelLinearOperator subtypes ##############
############################################################

############################################################
# joSALinearOperator #######################################

export joSALinearOperator, joSALinearOperatorException
export joSAdistributedLinearOperator, joSAdistributedLinearOperatorException
export joSAdistributingLinearOperator, joSAdistributingLinearOperatorException
export joSAgatheringLinearOperator, joSAgatheringLinearOperatorException

# type definition
"""
joSALinearOperator is glueing type & constructor

!!! Do not use it to create the operators

"""
struct joSALinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelLinearOperator{DDT,RDT,N}
    name::String
    m::Integer
    n::Integer
    nvc::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    PAs_in::joPAsetup    # input distributor
    PAs_out::joPAsetup   # output distributor
    fclean::Bool         # clean input vector post forward
    rclean::Bool         # clean input vector post reverse
end
"""
joSAdistributedLinearOperator is glueing type & constructor

!!! Do not use it to create the operators

"""
struct joSAdistributedLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelLinearOperator{DDT,RDT,N}
    name::String
    m::Integer
    n::Integer
    nvc::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    PAs_in::joPAsetup   # input distributor
    PAs_out::joPAsetup  # output distributor
    fclean::Bool        # clean input vector post forward
    rclean::Bool        # clean input vector post reverse
end
"""
joSAdistributingLinearOperator is glueing type & constructor

!!! Do not use it to create the operators

"""
struct joSAdistributingLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelLinearOperator{DDT,RDT,N}
    name::String
    m::Integer
    n::Integer
    nvc::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    PAs_out::joPAsetup  # output distributor
    gclean::Bool        # clean input vector post gathering
end
"""
joSAgatheringLinearOperator is glueing type & constructor

!!! Do not use it to create the operators

"""
struct joSAgatheringLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelLinearOperator{DDT,RDT,N}
    name::String
    m::Integer
    n::Integer
    nvc::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    PAs_in::joPAsetup  # input distributor
    gclean::Bool       # clean input vector post gathering
end

# type exception
struct joSALinearOperatorException <: Exception
    msg :: String
end
struct joSAdistributedLinearOperatorException <: Exception
    msg :: String
end
struct joSAdistributingLinearOperatorException <: Exception
    msg :: String
end
struct joSAgatheringLinearOperatorException <: Exception
    msg :: String
end

# unions
const joSALinOpsUnion{D,R,N}=Union{joSALinearOperator{D,R,N},joSAdistributedLinearOperator{D,R,N}} where {D,R,N}
const joSAdistributeLinOpsUnion{D,R,N}=Union{joSAdistribute{D,R,N},joSAdistributingLinearOperator{D,R,N}} where {D,R,N}
const joSAgatherLinOpsUnion{D,R,N}=Union{joSAgather{D,R,N},joSAgatheringLinearOperator{D,R,N}} where {D,R,N}


text:
using EDA
using Test

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


text:
using Pkg
Pkg.activate(".")

using Distributed

function s_rand()
    n = 10^4
    x = 0.0
    for i in 1:n
        x += sum(rand(10^4))
    end
    x / n
end

function p_rand()
    n = 10^4
    x = @distributed (+) for i in 1:n
        sum(rand(10^4))
    end
    x / n
end

@time s_rand()
@time s_rand()
@time p_rand()
@time p_rand()


text:
export isefficient

"""
An `AbstractOperator` is the supertype of all objects that map between function
spaces.
"""
abstract type AbstractOperator
end

dest(op::AbstractOperator) = _dest(op, dest_space(op))
_dest(op::AbstractOperator, span::Span) = dictionary(span)
_dest(op::AbstractOperator, space) = error("Generic operator does not map to the span of a dictionary.")

hasspan_dest(op::AbstractOperator) = typeof(dest_space(op)) <: Span

function apply(op::AbstractOperator, f; options...)
	if hasspan_dest(op)
		result = zeros(dest(op))
		apply!(result, op, f; options...)
	else
		error("Don't know how to apply operator $(string(op)). Please implement.")
	end
end

size(op::AbstractOperator) = _size(op, dest_space(op), src_space(op))
_size(op::AbstractOperator, dest, src) = (length(dest), length(src))

size(op::AbstractOperator, j::Int) = j <= 2 ? size(op)[j] : 1

"""
`DictionaryOperator` represents any linear operator that maps coefficients of
a source set to coefficients of a destination set. Typically, source and
destination are of type `Dictionary`.
The action of the operator is defined by providing a method for apply!.

The dimension of an operator are like a matrix: `(length(dest),length(src))`.

Source and destination should at least implement the following:
- `length`
- `size`
- `eltype`

The element type should be equal for src and dest.
"""
abstract type DictionaryOperator{T} <: AbstractOperator
end

eltype(::Type{<:DictionaryOperator{T}}) where {T} = T

# Default implementation of src and dest: assume they are fields
src(op::DictionaryOperator) = op.src
dest(op::DictionaryOperator) = op.dest

src_space(op::DictionaryOperator) = Span(src(op))
dest_space(op::DictionaryOperator) = Span(dest(op))

isreal(op::DictionaryOperator) = isreal(eltype(op)) && isreal(src(op)) && isreal(dest(op))

"""
True if the action of the operator has a better computational complexity than
the corresponding matrix-vector product.
"""
isefficient(op::DictionaryOperator) = false

size(op::DictionaryOperator) = (length(dest(op)), length(src(op)))

"Is the action of the operator in-place?"
isinplace(op::DictionaryOperator) = false

"Is the operator diagonal?"
isdiag(op::DictionaryOperator) = false

"Is the operator (knowingly) an identity operator?"
isidentity(op::DictionaryOperator) = false

function apply(op::DictionaryOperator, coef_src)
	coef_dest = zeros(eltype(op), dest(op))
	apply!(op, coef_dest, coef_src)
	coef_dest
end

# Catch applications of an operator, and do:
# - call inline implementation of operator if available
# - call apply!(op, dest(op), src(op), coef_dest, coef_src), which can be
#   implemented by operators whose action depends on src and/or dest.
function apply!(op::DictionaryOperator, coef_dest, coef_src)
	if isinplace(op)
		copyto!(coef_dest, coef_src)
		apply_inplace!(op, coef_dest)
	else
		apply_not_inplace!(op, coef_dest, coef_src)
	end
	# We expect each operator to return coef_dest, but we repeat here to make
	# sure our method is type-stable.
	coef_dest
end

# Provide a general dispatchable definition for in-place operators also
function apply!(op::DictionaryOperator, coef_srcdest)
		apply_inplace!(op, coef_srcdest)
		coef_srcdest
end

# Catch-all for missing implementations
apply_not_inplace!(op::DictionaryOperator, coef_dest, coef_src) =
	error("Operation of ", typeof(op), " on ", typeof(dest(op)), " and ", typeof(src(op)), " not implemented.")

# Catch-all for missing implementations
apply_inplace!(op::DictionaryOperator, coef_srcdest) =
	error("In-place operation of ", typeof(op), " on ", typeof(dest(op)), " and ", typeof(src(op)), " not implemented.")

"""
Apply an operator multiple times, to each column of the given argument.
"""
function apply_multiple(op::DictionaryOperator, matrix_src)
	T = eltype(op)
	matrix_dest = zeros(T, size(op,1), size(matrix_src)[2:end]...)
	coef_src = zeros(T, size(src(op)))
	coef_dest = zeros(T, size(dest(op)))
	apply_multiple!(op, matrix_dest, matrix_src, coef_dest, coef_src)
end

function apply_multiple!(op::DictionaryOperator, matrix_dest, matrix_src,
	coef_dest = zeros(eltype(matrix_dest), size(dest(op))),
	coef_src  = zeros(eltype(matrix_src), size(src(op))))

	# Make sure the first dimensions of the matrices agree with the dimensions
	# of the operator
	@assert size(matrix_dest, 1) == size(op,1)
	@assert size(matrix_src, 1) == size(op,2)
	# and that the remaining dimensions agree with each other
	for i in 2:ndims(matrix_dest)
		@assert size(matrix_src,i) == size(matrix_dest,i)
	end
	for i in 1:size(matrix_src,2)
		for j in eachindex(coef_src)
			coef_src[j] = matrix_src[j,i]
		end
		apply!(op, coef_dest, coef_src)
		for j in eachindex(coef_dest)
			matrix_dest[j,i] = coef_dest[j]
		end
	end
	matrix_dest
end

collect(op::DictionaryOperator) = matrix(op)

function sparse(op::DictionaryOperator; sparse_tol = 1e-14, options...)
	T = eltype(op)
	coef_src  = zeros(T, src(op))
    coef_dest = zeros(T, dest(op))
    R = spzeros(T, size(op,1), 0)
    for (i,si) in enumerate(eachindex(coef_src))
        coef_src[si] = 1
        apply!(op, coef_dest, coef_src)
        coef_src[si] = 0
        coef_dest[abs.(coef_dest).<sparse_tol] .= 0
        R = hcat(R,reshape(sparse(coef_dest),length(coef_dest)))
    end
    R
end

Base.Matrix(op::DictionaryOperator) = matrix(op)

matrix(op::DictionaryOperator) =
	default_matrix(op)

function default_matrix(op::DictionaryOperator)
    a = Array{eltype(op)}(undef, size(op))
    matrix!(op, a)
end

function matrix!(op::DictionaryOperator, a)
	T = eltype(op)
    n = length(src(op))
    m = length(dest(op))

    @assert (m,n) == size(a)

    coef_src  = zeros(T, src(op))
    coef_dest = zeros(T, dest(op))

    matrix_fill!(op, a, coef_src, coef_dest)
end

function matrix_fill!(op::DictionaryOperator, a, coef_src, coef_dest)
	T = eltype(op)
    for (i,si) in enumerate(eachindex(coef_src))
		coef_src[si] = one(T)
        apply!(op, coef_dest, coef_src)
		coef_src[si] = zero(T)
        for (j,dj) in enumerate(eachindex(coef_dest))
            a[j,i] = coef_dest[dj]
        end
    end
    a
end

function checkbounds(op::DictionaryOperator, i::Int, j::Int)
	1 <= i <= size(op,1) || throw(BoundsError())
	1 <= j <= size(op,2) || throw(BoundsError())
end

function getindex(op::DictionaryOperator, i, j)
	checkbounds(op, i, j)
	unsafe_op_getindex(op, i, j)
end

function unsafe_op_getindex(op::DictionaryOperator, i, j)
	T = eltype(op)
	coef_src = zeros(T, src(op))
	coef_dest = zeros(T, dest(op))
	coef_src[j] = one(T)
	apply!(op, coef_dest, coef_src)
	coef_dest[i]
end

"Return the diagonal of the operator."
function diag(op::DictionaryOperator)
    if isdiag(op)
        # Make data of all ones in the native representation of the operator
        all_ones = ones(src(op))
        # Apply the operator: this extracts the diagonal because the operator is diagonal
        diagonal_native = apply(op, all_ones)
        # Convert to vector
        linearize_coefficients(dest(op), diagonal_native)
    else
		# Compute the diagonal by calling unsafe_op_diag for each index
        [unsafe_op_diag(op, i) for i in 1:min(length(src(op)),length(dest(op)))]
    end
end

"Return the diagonal element op[i,i] of the operator."
function diag(op::DictionaryOperator, i)
	# Perform bounds checking and call unsafe_op_diag
	checkbounds(op, i, i)
	unsafe_op_diag(op, i)
end

# Default behaviour: call unsafe_op_getindex
unsafe_op_diag(op::DictionaryOperator, i) = unsafe_op_getindex(op, i, i)

# We provide a default implementation for diagonal operators
function pinv(op::DictionaryOperator, tolerance=eps(real(eltype(op))))
    @assert isdiag(op)
    newdiag = copy(diag(op))
    for i = 1:length(newdiag)
        newdiag[i] = abs(newdiag[i])>tolerance ? newdiag[i].^(-1) : 0
    end
    DiagonalOperator(dest(op),src(op), newdiag)
end

for f in (:eigvals, :svdvals, :norm, :rank)
    @eval $f(op::DictionaryOperator) = $f(Matrix(op))
end#matrix related features

function ≈(op1::DictionaryOperator,op2::DictionaryOperator; options...)
    r = rand(src(op1))
    if isapprox(op1*r,op2*r; options...)
		return true
	else
		@debug "Approx gives difference of $(norm(op1*r-op2*r))"
		return false
	end
end

"""
The function wrap_operator returns an operator with the given source and destination,
and with the action of the given operator.
"""
function wrap_operator(w_src::Dictionary, w_dest::Dictionary, op::DictionaryOperator)
    # We do some consistency checks
    @assert size(w_src) == size(src(op))
    @assert size(w_dest) == size(dest(op))
    @assert promote_type(eltype(op),operatoreltype(w_src,w_dest)) == eltype(op)
    unsafe_wrap_operator(w_src, w_dest, op)
end

show(io::IO, mime::MIME"text/plain", op::DictionaryOperator) = composite_show(io, mime, op)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function precess_xyz(x::T, y::T, z::T, equinox1::T, equinox2::T) where {T<:AbstractFloat}
    ra  = atan(y, x)
    del = norm((x, y, z)) #  Magnitude of distance to Sun
    dec = asin(z/del)
    # precess the ra and dec
    ra, dec = precess(ra, dec, equinox1, equinox2, radians=true)
    # convert back to x, y, z
    zunit, cos_dec = sincos(dec)
    yunit, xunit = cos_dec .* sincos(ra)
    return xunit*del, yunit*del, zunit*del
end

"""
    precess_xyz(x, y, z, equinox1, equinox2) -> prec_x, prec_y, prec_z

### Purpose ###

Precess equatorial geocentric rectangular coordinates.

### Arguments ###

* `x`, `y`, `z`: scalars or vectors giving heliocentric rectangular coordinates.
* `equinox1`: original equinox of coordinates, numeric scalar.
* `equinox2`: equinox of precessed coordinates, numeric scalar.

Input coordinates can be given also a 3-tuple `(x, y, z)`.

### Output ###

The 3-tuple `(x, y, z)` of coordinates modified by precession.

### Example ###

Precess 2000 equinox coordinates `(1, 1, 1)` to 2050.

```jldoctest
julia> using AstroLib

julia> precess_xyz(1, 1, 1, 2000, 2050)
(0.9838854500981734, 1.0110925876508692, 1.0048189888146941)
```

### Method ###

The equatorial geocentric rectangular coordinates are converted to right
ascension and declination, precessed in the normal way, then changed back to
`x`, `y` and `z` using unit vectors.

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
precess_xyz(x::Real, y::Real, z::Real,
            equinox1::Real, equinox2::Real) =
                precess_xyz(promote(float(x), float(y), float(z),
                                    float(equinox1), float(equinox2))...)

precess_xyz(xyz::Tuple{Real, Real, Real}, equinox1::Real, equinox2::Real) =
    precess_xyz(xyz..., equinox1, equinox2)

function precess_xyz(x::AbstractArray{X}, y::AbstractArray{<:Real},
                     z::AbstractArray{<:Real}, equinox1::Real,
                     equinox2::Real) where {X<:Real}
    @assert length(x) == length(y) == length(z) "x, y, z arrays should be of the same length"
    typex = float(X)
    x_out = similar(x, typex)
    y_out = similar(x, typex)
    z_out = similar(x, typex)
    for i in eachindex(x)
        x_out[i], y_out[i], z_out[i] = precess_xyz(x[i], y[i], z[i],
                                                   equinox1, equinox2)
    end
    return x_out, y_out, z_out
end


text:
# Formulations of relaxed AC-OPF based on PowerModels.jl

"""
InfrastructureModels build function to formulate AC OPF with variable loads
"""
function build_opf_var_load(pm)
    PM.variable_bus_voltage(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)
	variable_load(pm)

    PM.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        PM.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_power_balance_var_loads(pm, i)
    end

    for i in ids(pm, :branch)
        PM.constraint_ohms_yt_from(pm, i)
        PM.constraint_ohms_yt_to(pm, i)

        PM.constraint_voltage_angle_difference(pm, i)

        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

"""
InfrastructureModels build function to formulate AC OPF with variable loads constrained
by the minimum powerfactor, and minimum and maximum real power if given in the
PowerModels network data dictionary used to create the PowerModels model.
"""
function build_opf_constr_load(pm)
    PM.variable_bus_voltage(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)
	variable_load(pm)
	constraint_load(pm)

    PM.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        PM.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_power_balance_var_loads(pm, i)
    end

    for i in ids(pm, :branch)
        PM.constraint_ohms_yt_from(pm, i)
        PM.constraint_ohms_yt_to(pm, i)

        PM.constraint_voltage_angle_difference(pm, i)

        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

"""
Add variables for the active and reactive components of loads in the network to
the given PowerModels model
"""
function variable_load(pm)
	model = pm.model

	bus_loads = ref(pm, :bus_loads)
	num_loads = length(PM.ref(pm, :load))

	Pd = JuMP.@variable(model, pd[1:num_loads], container=DenseAxisArray)
	Qd = JuMP.@variable(model, qd[1:num_loads], container=DenseAxisArray)
	pm.var[:it][:pm][:nw][0][:pd] = Pd
	pm.var[:it][:pm][:nw][0][:qd] = Qd
end

"""
Constrains the load variables in a PowerModels AC OPF formulation model.
Constrains the reactive power with the given minimum power factor, pf_min,
maximum active power, pd_max, and minimum active power, pd_min. These constraints
are pulled from the PowerModels network data dictionary used to create the
PowerModels network model if the dictionary has these constraint keys (These keys
are to be added to the dictionary after loading it in using PowerModels).
Additionally, the power factor is constrained to lagging values by default but
can be specified to include leading values by including the constraint option
pf_lagging => false.
"""
function constraint_load(pm; nw::Int=NW_DEFAULT, pd_min=nothing)
	model = pm.model
	load_ids = ids(pm, :load)  # Load indeces, 1 to num_loads
	num_loads = length(load_ids)

	if haskey(pm.data, "pd_max")
		pd_max = pm.data["pd_max"]  # Pre-calculated & stored max active powers
		@assert length(pd_max) == num_loads
	else
		pd_max = nothing
	end

	if haskey(pm.data, "pd_min")
		pd_min = pm.data["pd_min"]  # Specified min active powers
		@assert length(pd_min) == num_loads
	else
		pd_min = nothing
	end

	if haskey(pm.data, "pf_min")
		pf_min = pm.data["pf_min"]  # Selected minimum power factor
		@assert ((length(pf_min) == num_loads) | (length(pf_min) == 1))
	else
		pf_min = nothing
	end

	if haskey(pm.data, "pf_lagging")
		pf_lagging = pm.data["pf_lagging"]  # Contrain load pf to lagging or allow leading
	else
		pf_lagging = true
	end

	pd = get(PM.var(pm, nw), :pd, Dict())
	qd = get(PM.var(pm, nw), :qd, Dict())

	# Load Bus Constraints
	for load in load_ids  #TASK: Determine how load variables align with bus order
		if !isnothing(pd_min)
			JuMP.set_lower_bound(pd[load], pd_min[load])
		else
			JuMP.set_lower_bound(pd[load], 0)
		end
		if !isnothing(pd_max)  # Contrain max real load if data is available
			JuMP.set_upper_bound(pd[load], pd_max[load])
		end
	end

	# Reactive load powers contrained to a ratio of Real load power if pf_min is not nothing
	if !isnothing(pf_min)
		for load_idx in load_ids
			d = tan(acos(pf_min[load_idx]))
			JuMP.@constraint(model, qd[load_idx] <= (pd[load_idx] * d))
			if pf_lagging
				JuMP.@constraint(model, qd[load_idx] >= 0)
			else
				JuMP.@constraint(model, qd[load_idx] >= (pd[load_idx] * -d))
			end
		end
	end
end

"""
Adds the AC OPF power balance constraint to the formulation of a PowerModels model with variable loads.
Modified from the PowerModels function 'constraint_power_balance' which can be found here,
https://github.com/lanl-ansi/PowerModels.jl/blob/460b310288787d4196f9e50b1a81127ee4677a97/src/core/constraint_template.jl#L174-L191

Copyright (c) 2016, Los Alamos National Security, LLC All rights reserved. Copyright 2016. Los Alamos National Security, LLC.
"""
function constraint_power_balance_var_loads(pm::PM.AbstractPowerModel, i::Int; nw::Int=NW_DEFAULT)
    bus = ref(pm, nw, :bus, i)
    bus_arcs = ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
    bus_arcs_sw = ref(pm, nw, :bus_arcs_sw, i)
    bus_gens = ref(pm, nw, :bus_gens, i)
    bus_loads = ref(pm, nw, :bus_loads, i)
    bus_shunts = ref(pm, nw, :bus_shunts, i)
    bus_storage = ref(pm, nw, :bus_storage, i)

	# PowerModels gathers float values for loads in the system
    #bus_pd = Dict(k => ref(pm, nw, :load, k, "pd") for k in bus_loads)
    #bus_qd = Dict(k => ref(pm, nw, :load, k, "qd") for k in bus_loads)

	# Replace with load variables
	pd = get(PM.var(pm, nw), :pd, Dict())
	qd = get(PM.var(pm, nw), :qd, Dict())

	bus_pd = Dict(k => pd[k] for k in bus_loads)
    bus_qd = Dict(k => qd[k] for k in bus_loads)

    bus_gs = Dict(k => ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bus_bs = Dict(k => ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)

    PM.constraint_power_balance(pm, nw, i, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
end

"""
Determines the maximum servable active load for each load bus with the given system information.
Formulates the problem in JuMP to maximize the load at each load bus individually using
the given solver.

min_load: The minimum active load constraint as a percentage of the nominal load. Given as either
		  a single number or a list of number of length num_loads
print_level: An integer from 0-2 that determines how much information is printed to console.
"""
function find_max_loads(pm; print_level=0, nw=NW_DEFAULT, min_load=0.0,
				     solver=optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	model = pm.model

	JuMP.set_optimizer(model, solver)
	if print_level < 2
		JuMP.set_silent(pm.model)
	end

	pd = get(PM.var(pm, nw), :pd, Dict())
	qd = get(PM.var(pm, nw), :qd, Dict())
	loads = ref(pm, :load)
	num_loads = length(PM.ref(pm, :load))
	num_buses = length(ref(pm, :bus))

	# Constrain the load pf and pd max/min
	nominal_load = pm.data["base_load"]
	min_load isa Number && (min_load = ones(num_loads) .* min_load)
	pd_min = min_load .* nominal_load[1:num_loads]

	constraint_load(pm, pd_min=pd_min)

	solved_statuses = []
	pd_max = zeros(num_loads, 1)
	for load_idx in sort(collect(keys(loads)))
		# Maximize Load Objective
		JuMP.@objective(model, Max, pd[load_idx]) # Set objective to maximize the load at a load bus

		JuMP.optimize!(model)
		load_bus = loads[load_idx]["load_bus"]
		pd_max[load_idx] = JuMP.objective_value(model)

		if print_level > 0
			println("")
			println("Load Bus: ", load_bus)
			println(JuMP.solution_summary(model))
		end

		# Check if solution was found & no errors occured
		status = JuMP.termination_status(model)
		solved = status in SOLVED_STATUSES
		append!(solved_statuses, solved)
		!solved && @warn "Load $(load) max load solution is $(status)"
	end
	return pd_max, solved_statuses
end

"""
Initializes a PowerModels model, pm, with variable loads, to be used for finding
the nearest feasible load profile.
"""
function find_nearest_feasible_model(pm; nw=NW_DEFAULT, print_level=0,
							       solver=optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	!haskey(pm.data, "pd_max") && (@warn "Max load constraint values not provided.")

	model = pm.model

	JuMP.set_optimizer(model, solver)
	if print_level < 2
		JuMP.set_silent(pm.model)
	end

	constraint_load(pm)
	return pm
end

"""
Finds the nearest feasible load profile to given load profile, with active and reactive loads, px and qx,
for the PowerModels network, pm.
"""
function find_nearest_feasible(pm, px, qx; nw=NW_DEFAULT, print_level=0)
	model = pm.model
	num_loads = length(PM.ref(pm, :load))

	# Objective, RMS of distance between sample load points in Px & Qx and potentail feasible loads Pd and Qd
	#quad_form(Pd - Px, eye(numbus))+ quad_form(Qd - Qx, eye(numbus)
	pd = get(PM.var(pm, nw), :pd, Dict())
	qd = get(PM.var(pm, nw), :qd, Dict())

	JuMP.unregister(model, :R)

	JuMP.@expression(model, R, sum((pd[i] - px[i])^2 + (qd[i] - qx[i])^2 for i in 1:num_loads))
	JuMP.@objective(model, MOI.FEASIBILITY_SENSE, 0)  # Set sense to delete the bridge, Needed for SCS & Hypatia
	JuMP.@objective(model, Min, R)  # Add the new objective

	JuMP.optimize!(model)

	r = JuMP.objective_value(model)

	PD_vals = JuMP.value.(model[:pd]).data
	QD_vals = JuMP.value.(model[:qd]).data

	if print_level > 1
		println(JuMP.solution_summary(model))
	end

	# Check if solution was found & no errors occured
	status = JuMP.termination_status(model)
	solved = status in SOLVED_STATUSES
	!solved && @warn "Nearest feasible load solution is $(status)"

	return r, PD_vals, QD_vals, solved
end


text:
using OffsetArrays

mutable struct BCTypeDirichlet end

mutable struct BCTypeNeumann end

mutable struct BCTypeRobbins end

mutable struct WaterLevelBC end

mutable struct PressureBC end

mutable struct VelocityBC end

mutable struct FluxBC end

mutable struct BCLocation{T}
    N::Int64
    dim::Int64
    cells::T
end

mutable struct BCValues{T}
N::Int64
    dim::Int64
    values::T
 end

mutable struct BoundaryCondition{string,BCType,BCParameter,BCLocation,BCValues}
    name::string
    type::BCType
    parameter::BCParameter
    location::BCLocation
    values::BCValues
end

ndims = 3
nx  = 10
ny  = 5
nz = 20

dirichlet_bc = BCTypeDirichlet()
neumann_bc   = BCTypeNeumann()
robbins_bc   = BCTypeRobbins()

pressure_bc    = PressureBC()
water_level_bc = WaterLevelBC()
velocity_bc    = VelocityBC()
flux_bc        = FluxBC()

bc_cells  = OffsetArray{Int64}(undef, 0:nx + 1, 0:ny + 1, 0:nz + 1)
bc_values = OffsetArray{Float64}(undef, 0:nx + 1, 0:ny + 1, 0:nz + 1)

bc_cells .= 0
bc_cells .= 0.0

for k = 1:nz
    for j = 1:ny
        bc_cells[0, j,k] = 1
        bc_values[0, j,k] = 1.0
    end
end

location_bc = BCLocation(ny, ndims, bc_cells)
values_bc  = BCValues(ny, ndims, bc_values)

bc_name = "west_bc"

bc_west = BoundaryCondition(bc_name, dirichlet_bc, pressure_bc, location_bc, values_bc)

bc_west.type
bc_west.parameter