text:
#!/usr/bin/env julia

# Project 1: Vanilla options in a Black-Scholes World
#
# Page 437+, Joshi's Mathematical Finance

# Best to run this with JULIA_NUM_THREADS=$(nproc) and julia 1.3

# Notation: if f is a function of x, fv is one realisation of that function

import ThreadsX
using SpecialFunctions: erf

export C, P, B, present_value, rand_price, mc_pricer, binary, mc_pricer_pathdep

N(x) = 0.5 * (1 + erf(x/sqrt(2)))

d_1(S,t,T,K,r,σ) = (log(S/K) + (r + (σ^2)/2)*(T-t)) / (σ*√(T-t))
d_2(S,t,T,K,r,σ) = (log(S/K) + (r - (σ^2)/2)*(T-t)) / (σ*√(T-t))

# Price of a call option (i.e. the holder has the right to buy stock at maturity)
# S = current stock (underlying) price
# t = current time (time since option created)
# T = time of maturity
# K = the price at which the stock may be bought at maturity
# r = the continuously compounded risk-free rate
# σ = volatility of the underlying (standard deviation per root unit of time)
# TODO: d = dividend rate
function C(S,t; T=1, K=100, r=0.02, σ=0.5)
    d_1v = d_1(S,t,T,K,r,σ)
    d_2v = d_2(S,t,T,K,r,σ)
    N(d_1v)*S - N(d_2v)*present_value(K,r,t,T)
end

# Price of a put option (i.e. holder may sell at maturity)
function P(S,t; T=1, K=100, r=0.02, σ=0.5)
    present_value(K,r,t,T) - S + C(S,t, T=T, K=K, r=r, σ=σ)
end

# Price of a binary option:
# Payoff of 1 if S(T) > E, 0 otherwise
#
# Source: http://www.iam.fmph.uniba.sk/institute/stehlikova/fd14en/lectures/05_black_scholes_1.pdf
#
# Can construct a less than binary thusly:
# less_than_binary = present_value(1) - greater_than_binary
# (cf. put call parity)
function binary(S,t; T=1, E=100, r=0.02, σ=0.5)
    present_value(1,r,t,T) * N(d_2(S,t,T,E,r,σ))
end

##############################
#                            #
# Validation via Monte Carlo #
#                            #
##############################

# Brownian motion
B(T;B0=100,r=0.02,d=0.00,σ=0.5) = B0*exp((r-d)*T-0.5*σ^2*T+σ*√T*randn())

# Sketch of pricer: Brownian motion generates a final stock price; payoff of option is calculated and then discounted. Price of option is average of these.

present_value(v,r,t,T) = v*exp(-r*(T-t))

function rand_price(payoff;S_0=100,t=0,T=1,K=S_0,r=0.02,σ=0.05)
    S = B(T-t;B0=S_0,r=r,d=0,σ=σ)
    present_value(payoff(S),r,t,T)
end

# Strictly we should use a thread-safe RNG here
function mc_pricer(
    payoff;S_0=100,t=0,T=1,K=S_0,r=0.02,σ=0.05, trials=100_000_000
)
    linear, squared = ThreadsX.mapreduce(
        x->begin
            p = rand_price(payoff;S_0=S_0,t=t,T=T,K=K,r=r,σ=σ)
            (p,p^2) # Second term is to keep track of standard error
        end,
        (a,b)->a.+b,
        Base.OneTo(trials),
        init=(0.0, 0.0)
    )

    mean = linear/trials

    (mean = mean, sem = sqrt(squared/trials - mean^2)/sqrt(trials))
end

function mc_pricer_pathdep(
    payoff;S_0=100,t=[0],T=1,K=S_0,r=0.02,σ=0.05, trials=100_000_000
)
    linear, squared = ThreadsX.mapreduce(
        x->begin
            p = B.(T.-t;B0=S_0,r=r,d=0,σ=σ)
            price = present_value(payoff(p),r,t[1],T) # Payoff must operate on time series
            (price,price^2)
        end,
        (a,b)->a.+b,
        Base.OneTo(trials),
        init=(0.0, 0.0)
    )

    mean = linear/trials

    (mean = mean, sem = sqrt(squared/trials - mean^2)/sqrt(trials))
end


text:
using .TestUtils

@testgroup "builtin_fixtures" begin

@testcase "run_options" begin
    tcs = Jute.collect_testobjs() do
        @testcase "tc" for opts in run_options
            @test opts[:test_file_postfix] == ".test2.jl"
        end
    end

    # FIXME: when custom options are available, we will need to use one of those instead.
    # Meanwhile, we will use an option that will not affect the nested test run.
    exitcode = nested_run(tcs, Dict(:test_file_postfix => ".test2.jl"))
    @test exitcode == 0
end

@testcase "temporary_dir" begin

    tcs = Jute.collect_testobjs() do
        tdir = "---"

        @testcase "tc1" for dir in temporary_dir
            @test isdir(dir)
            tdir = dir
            open(joinpath(dir, "testfile"), "w") do f
                write(f, "test")
            end
        end

        @testcase "tc2" begin
            @test tdir != "---"
            @test !isdir(tdir) # check that the directory was deleted
        end
    end

    exitcode = nested_run(tcs)
    @test exitcode == 0
end

@testcase "fixed_rng" begin

    rand_ints = [-1, -1, -1]

    tcs = Jute.collect_testobjs() do
        @testcase "tc1" for rng in fixed_rng("abc")
            rand_ints[1] = rand(rng, 1:1000000)
        end
        @testcase "tc2" for rng in fixed_rng("abc")
            rand_ints[2] = rand(rng, 1:1000000)
        end
        @testcase "tc3" for rng in fixed_rng("def")
            rand_ints[3] = rand(rng, 1:1000000)
        end
    end

    exitcode = nested_run(tcs)
    @test exitcode == 0

    @test rand_ints[1] == rand_ints[2]
    @test rand_ints[1] != rand_ints[3]
end

end

ext: jl


text:
@testset "Maps and Reductions" begin
    @testset "map" begin
        rng = MersenneTwister(42)
        n = 10
        x = randn(rng, n)
        vx = randn(rng, n)
        ȳ = randn(rng, n)
        rrule_test(map, ȳ, (sin, nothing), (x, vx))
        rrule_test(map, ȳ, (+, nothing), (x, vx), (randn(rng, n), randn(rng, n)))
    end
    @testset "mapreduce" begin
        rng = MersenneTwister(6)
        n = 10
        x = randn(rng, n)
        vx = randn(rng, n)
        ȳ = randn(rng)
        rrule_test(mapreduce, ȳ, (sin, nothing), (+, nothing), (x, vx))

        # With keyword arguments (not yet supported in rrule_test)
        X = randn(rng, n, n)
        y, pullback = rrule(mapreduce, abs2, +, X; dims=2)
        ȳ = randn(rng, size(y))
        (_, _, _, x̄_ad) = pullback(ȳ)
        x̄_fd = j′vp(central_fdm(5, 1), x->mapreduce(abs2, +, x; dims=2), ȳ, X)
        @test x̄_ad ≈ x̄_fd atol=1e-9 rtol=1e-9
    end
    @testset "$f" for f in (mapfoldl, mapfoldr)
        rng = MersenneTwister(10)
        n = 7
        x = randn(rng, n)
        vx = randn(rng, n)
        ȳ = randn(rng)
        rrule_test(f, ȳ, (cos, nothing), (+, nothing), (x, vx))
    end
    @testset "sum" begin
        @testset "Vector" begin
            rng, M = MersenneTwister(123456), 3
            frule_test(sum, (randn(rng, M), randn(rng, M)))
            rrule_test(sum, randn(rng), (randn(rng, M), randn(rng, M)))
        end
        @testset "Matrix" begin
            rng, M, N = MersenneTwister(123456), 3, 4
            frule_test(sum, (randn(rng, M, N), randn(rng, M, N)))
            rrule_test(sum, randn(rng), (randn(rng, M, N), randn(rng, M, N)))
        end
        @testset "Array{T, 3}" begin
            rng, M, N, P = MersenneTwister(123456), 3, 7, 11
            frule_test(sum, (randn(rng, M, N, P), randn(rng, M, N, P)))
            rrule_test(sum, randn(rng), (randn(rng, M, N, P), randn(rng, M, N, P)))
        end
        @testset "function argument" begin
            rng = MersenneTwister(1)
            n = 8
            rrule_test(sum, randn(rng), (cos, nothing), (randn(rng, n), randn(rng, n)))
            rrule_test(sum, randn(rng), (abs2, nothing), (randn(rng, n), randn(rng, n)))
        end
        @testset "keyword arguments" begin
            rng = MersenneTwister(33)
            n = 4
            X = randn(rng, n, n+1)
            y, pullback = rrule(sum, X; dims=2)
            ȳ = randn(rng, size(y))
            _, x̄_ad = pullback(ȳ)
            x̄_fd = j′vp(central_fdm(5, 1), x->sum(x, dims=2), ȳ, X)
            @test x̄_ad ≈ x̄_fd atol=1e-9 rtol=1e-9
        end
    end  # sum
end


text:
module FESetShellQ4Module

using FinEtools
using FinEtools.MeshQuadrilateralModule
import FinEtools.FESetModule: cat, subset, nodesperelem
using FinEtools.MatrixUtilityModule: complete_lt!
using LinearAlgebra: norm, Transpose, mul!

struct FESetShellQ4 <: AbstractFESet2Manifold{4}
end

function local_frame!(fes::FESetShellQ4, F0, J0)
    # This is the tangent to the coordinate curve 1
    a = @view J0[:, 1]
    L0 = norm(a);
    F0[:,1] = a/L0;
    # This is the tangent to the coordinate curve 2
    b = @view J0[:, 2]
    #     F0(:,3)=skewmat(F0(:,1))*b;
    # Now compute the normal
    F0[:, 3] .= (-F0[3,1]*b[2]+F0[2,1]*b[3],
                  F0[3,1]*b[1]-F0[1,1]*b[3],
                 -F0[2,1]*b[1]+F0[1,1]*b[2]);
    F0[:, 3] /= norm(@view F0[:, 3]);
    #     F0(:,2)=skewmat(F0(:,3))*F0(:,1);
    F0[:, 2] .= (-F0[3,3]*F0[2,1]+F0[2,3]*F0[3,1],
                  F0[3,3]*F0[1,1]-F0[1,3]*F0[3,1],
                 -F0[2,3]*F0[1,1]+F0[1,3]*F0[2,1]);
    return F0
end

end # module


text:
# Optimize a network derived from a cicular ladder shaped connectome.

using Plots
using SpikeNetOpt
SNO = SpikeNetOpt
using SpikingNeuralNetworks
using Evolutionary
SNN = SpikingNeuralNetworks
using Metaheuristics
#using DrWatson
#using Pkg
SNN.@load_units

##
# Ground truths
##
const Ne = 200;
const Ni = 50
const σee = 1.0
const pee = 0.5
const σei = 1.0
const pei = 0.5
MU = 10

const E
const spkd_ground
const GT = 26

g, Cg = SpikeNetOpt.make_net_from_graph_structure(GT)
P, C = SpikeNetOpt.make_net_SNN(Ne, Ni, σee = 0.5, pee = 0.8, σei = 0.5, pei = 0.8)
E, I = P
EE, EI, IE, II = C
SNN.monitor([E, I], [:fire])
sim_length = 1000
@inbounds for t = 1:sim_length*ms
    E.I = vec([11.5 for i = 1:sim_length])
    SNN.sim!(P, C, 1ms)

end
spkd_ground = SpikeNetOpt.get_trains(P[1])
display(SNN.raster(P[1]))
SNN.raster(P[1]) |> display
META_HEUR_OPT = true

function loss(model)
    """
    A loss function for calculating errors in machine learning.
    """
    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = SNO.make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 1000
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display
    error = SNO.spike_train_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)
    error
end

D = 10
bounds = [3ones(D) 40ones(D)]'
a = view(bounds, 1, 1)
b = view(bounds, 1, 2)
information = Information(f_optimum = 0.0)
options = Options( seed = 1, iterations=10, f_calls_limit =10)

D = size(bounds, 2)
nobjectives=1

options = Options( seed = 1, iterations=10000, f_calls_limit = 25000)
npartitions = nobjectives == 2 ? 100 : 12

methods = [
        SMS_EMOA(N = 5, n_samples=5, options=options),
        NSGA2(options=options)
        ]

for method in methods
    f_calls = 0
    result = optimize(SpikeNetOpt.loss, bounds, method)
    @show(result)

end


text:
"""
Compatibility layer for transferring from FixedSizeArrays. This provides
alternative definitions of `Vec`, `Mat`, `Point`, `FixedVectorNoTuple`, `@fsa`,
etc, using StaticArrays as a backend.

The type definitions are not "perfect" matches because the type parameters are
different. However, it should cover common method signatures and constructors.
"""
module FixedSizeArrays

using ..StaticArrays

export FixedArray
export FixedVector
export FixedMatrix
export Mat, Vec, Point
export @fsa
export FixedVectorNoTuple

const FixedArray = StaticArray
const FixedVector = StaticVector
const FixedMatrix = StaticMatrix
const Vec = SVector
const Mat = SMatrix
const FixedVectorNoTuple = FieldVector

macro fsa(ex)
    @assert isa(ex, Expr)
    if ex.head == :vect # Vector
        return Expr(:call, SVector{length(ex.args)}, Expr(:tuple, ex.args...))
    elseif ex.head == :hcat # 1 x n
        s1 = 1
        s2 = length(ex.args)
        return Expr(:call, SMatrix{s1, s2}, Expr(:tuple, ex.args...))
    elseif ex.head == :vcat
        if isa(ex.args[1], Expr) && ex.args[1].head == :row # n x m
            # Validate
            s1 = length(ex.args)
            s2s = map(i -> ((isa(ex.args[i], Expr) && ex.args[i].head == :row) ? length(ex.args[i].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i].args[j] for i = 1:s1, j = 1:s2]
            return Expr(:call, SMatrix{s1, s2}, Expr(:tuple, exprs...))
        else # n x 1
            return Expr(:call, SMatrix{length(ex.args), 1}, Expr(:tuple, ex.args...))
        end
    end
end

###########
## Point ##
###########

immutable Point{S, T} <: StaticVector{T}
    data::NTuple{S, T}
end

@inline (::Type{Point}){S}(x::NTuple{S}) = Point{S}(x)
@inline (::Type{Point{S}}){S, T}(x::NTuple{S,T}) = Point{S,T}(x)
@inline (::Type{Point{S}}){S, T <: Tuple}(x::T) = Point{S,promote_tuple_eltype(T)}(x)

Base.@pure Base.size{S}(::Union{Point{S},Type{Point{S}}}) = (S, )
Base.@pure Base.size{S,T}(::Type{Point{S,T}}) = (S,)

Base.@propagate_inbounds function Base.getindex(v::Point, i::Integer)
    v.data[i]
end

@inline Base.Tuple(v::Point) = v.data

end


text:
##RayTrace_lat.jl
# Todd Anderson
# 5 October 2020
#
# Ray-tracing of electromagnetic waves from upper ionosphere to plasmasphere using Haselgrove equations, following ray-tracing equations in Rice 1997 combined with plasmasphere model adapted from Carpenter and Anderson 1992.  Includes main ODE solver function calls.
#

using DifferentialEquations
using BenchmarkTools
using LSODA
using Sundials

c = 2.99792458e8
re = 6.3712e6      	# radius of Earth in meters
B0 = 3.0696381e-5   # magnitude of B field at Earth's equator surface in Tesla
e = 1.602e-19 		# electric charge in Coulombs
me = 9.1093e-31		# electron rest mass
mp = 1.6726219e-27	# proton rest mass
eps = 8.854e-12 	# permittivity of free space in mks

# pre-solve plasmasphere
r_range = [re:1000:(10*re);]
λ_range = 0
L = r_range./(re*cos(λ_range)^2)

Kp_max = 3                             # for example
Lppi = 5.6 - 0.46*Kp_max                # plasmapause inner limit in L
#Lppo = 5.75                             # placeholder
d = 0                                   # day number
R̄ = 90                                  # 13-month average sunspot number
mlt = 2                                   # magnetic local time

ne_Lppi = 10^((-0.3145*Lppi + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-Lppi)/1.5))

# pre-solve individual components in order to find Lppo
ne_plasma_1 = @. 10 .^((-0.3145.*L + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2 .-L)./1.5))
ne_plasma_2 = ne_Lppi*10 .^(-1.0.*(L .- Lppi)./0.1)
ne_plasma_3 = (5800 + 300*mlt).*L.^(-4.5) + (1 .- exp.((2 .-L)./10))

min_index = findmin(abs.(ne_plasma_2 - ne_plasma_3))[2]
Lppo = L[min_index]

# calculate phase refractive index μ
u = function phase_refractive_index(r,λ,χ,freq)
	# u = [μ, dμdψ]
    # convert from radial angle to wave normal angle
    dip = atan(2.0*tan(λ)) # dip angle: angle between horizontal and B field;

	#	ψ = 2*π - (3π/2 - dip - χ) 		# this should be ≡ to the below, but the plotted ray paths are not all that similar!
	ψ = π/2.0 + dip + χ
	#	ψ = χ - (3.0/2.0)*pi + dip		# wave normal angle: angle between wave direction and B; 11/9: fixed error (χ - ϕ) ➡ (ϕ - χ)
	# NOTE 11/24/2020: Rice 1997 uses ψ = χ - 3π/2 + dip, as well as a strange method of calculating χ (see FORTRAN code on p. 112)

    # convert from frequency to angular frequency
    ω = 2.0*pi*freq

    # find magnetic field at (r,θ) from dipole field model
    # Bmag = B0*((re^3)/(r^3))*sqrt(4.0-3.0*cos((pi/2.)-θ)*cos((pi/2.0)-θ))
	# 	   = B0*((re^3)/(r^3))*sqrt(4 - 3sin²θ)
	# 	   ≡ B0*((re^3)/(r^3))*sqrt(1 + 3cos²θ)
	# dipole field magnitude from Parks p. 61:
	# B(r,λ) = μ₀M/(4πr³)*(1+3sin²λ)^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3sin²(π/2 - θ))^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3cos²θ)^(1/2)
	# here B0*re³ = μ₀M/(4π)
	Bmag = B0*(re^3/(r^3))*sqrt(1+3*sin(λ)*sin(λ))

    # calculate electron and proton density profiles
	Lshell = r./(re*cos(λ)^2)
	ne_iono = (1.8e5.*exp.(-4.183119.*((r./re).-1.0471))) # cm^-3

	if Lshell <= Lppi
	     log_ne = (-0.3145*Lshell + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-Lshell)/1.5)
	     ne_plasma = 10^(log_ne)
	elseif Lppi < Lshell <= Lppo
	     ne_plasma = ne_Lppi*10^((Lppi-Lshell)/0.1)
	elseif Lppo < Lshell
	     ne_plasma = (5800 + 300*mlt)*Lshell^(-4.5) + (1 - exp((2-Lshell)/10))
	else
	     ne_plasma = 0.0
	end

	n_e = (ne_iono + ne_plasma)*1e6
	n_p = n_e

	# electron plasma angular frequency squared
    ω_e2 = (n_e*(e^2))/(eps*me)
	# proton plasma angular frequency squared
    ω_p2 = (n_p*(e^2))/(eps*mp)

	# electron angular gyrofrequency
    Ω_e = (e*Bmag)/(me)
	# proton angular gyrofrequency
    Ω_p = (e*Bmag)/(mp)

	# define R, L, P, D, S,
	# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
	R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
	#println("R = ",R)

	# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
	L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
	#println("L = ",L)

	# P ≡ 1 - Σ(ωₖ²/ω²)
	P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
	#println("P = ",P)

	# D ≡ ½(R-L); S ≡ ½(R+L)
	D = (R - L)/2.0
	S = (R + L)/2.0

	# define parts of dispersion relation
	# Aμ⁴ - Bμ² + C = 0
	# A = Ssin²ψ + Pcos²ψ
	A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0

	# B = RLsin²ψ + PS(1 + cos²ψ)
	B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0)

	# C  = PRL
	C = P*R*L

	# solve the dispersion relation for μ:
	# μ² = (B +- F)/2A
	# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
	F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
	F = sqrt(F2)

	# Rice 1997: typical solution to dispersion relation with quadratic formula
	μ2_minus = (B - F)/(2.0*A)
	μ2_plus = (B + F)/(2.0*A)

	# Bortnik 2004:
	# if B > 0
	# 	μ2_minus = (B - F)/(2.0*A)
	# else
	# 	μ2_plus = (2.0*C)/(B + F)
	# end

	μ_minus = try
		sqrt(abs(μ2_minus))	# abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			println("DomainError in μ_minus!")
			println("B=", B)
			println("F=", F)
			println("A=", A)
		else
			# println(err)
		end
	end

	μ_plus = try
		sqrt(abs(μ2_plus)) # abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			# println("DomainError in μ_plus!")
			# println("B=", B)
			# println("F=", F)
			# println("A=", A)
		else
			# println(err)
		end
	end

	# Electron whistler case: ψ = 0, μ² = R; this is the μ_plus case
	#μ = μ_minus		# EMIC case
	μ = μ_plus		# electron whistler case

	# Find dA/dψ, dB/dψ, dC/dψ, dμ/dψ
	dAdψ = 2.0*(S-P)*sin(ψ)*cos(ψ)
	dBdψ = 2.0*(R*L-P*S)*sin(ψ)*cos(ψ)
    dCdψ = 0.0
    #dμdψ = ((μ^4.0)*dAdψ-(μ^2.0)*dBdψ+dCdψ)/(4.0*A*(μ^3.0)-2.0*B*μ)

	dFdψ = 1/(2.0*F)*((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ))
	#dFdψ = sqrt(abs((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)))
	dμdψ = 1/(2.0*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B + F)/(2*A^2))
	# NOTE 11/24/2020: from Rice 1997; choose '+' for (B +- F) and (dBdψ +- dFdψ)

	#DEBUG check values
	#println("μ = ",μ)
	#println("dμdψ = ", dμdψ)

	# non-mutating output
	u = [μ,dμdψ]

end

# Derivatives of phase refractive index w.r.t. r, θ, χ,
# calculate derivative of μ w.r.t. r
dμdr = function ddr(r,λ,χ,freq)
	#μ,dμdψ = u			# unpack
	dr = 1.0e-11		# r step size

	μ_l = phase_refractive_index(r-dr/2.0,λ,χ,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r+dr/2.0,λ,χ,freq)
	#μ_r = μ

	dμdr = (μ_r[1] - μ_l[1])/dr
end

dμdλ = function ddλ(r,λ,χ,freq)
	#μ,dμdψ = u			# unpack
	dλ = 1.0e-11		# θ step size

	μ_l = phase_refractive_index(r,λ-dλ/2.0,χ,freq)
	#μ_l = μ
	#println("μ_l = ", μ)
	μ_r = phase_refractive_index(r,λ+dλ/2.0,χ,freq)
	#μ_r = μ
	#println("μ_r = ", μ)

	dμdλ = (μ_r[1] - μ_l[1])/dλ
end

dμdχ = function ddχ(r,λ,χ,freq)
	#μ,dμdψ = u			# unpack
	dχ = 1.0e-11		# χ step size

	μ_l = phase_refractive_index(r,λ,χ-dχ/2.0,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r,λ,χ+dχ/2.0,freq) # NOTE Nov 3: corrected error 'χ-dχ/2.0'➡'χ+dχ/2.0'
	#μ_r = μ

	dμdχ = (μ_r[1] - μ_l[1])/dχ
end

dμdf = function ddf(r,λ,χ,freq)
	#μ,dμdψ = u			# unpack
	df = 1.0e-5		# f step size

	μ_l = phase_refractive_index(r,λ,χ,freq-df/2.0)
	#μ_l = μ
	μ_r = phase_refractive_index(r,λ,χ,freq+df/2.0)
	#μ_r = μ

	dμdf = (μ_r[1] - μ_l[1])/df
end

# calculate derivatives w.r.t. time

# calculate derivatives w.r.t. time
function haselgrove!(du,u,p,t)
# 	du = drdt, dλdt, dχdt, dfdt
# 	u = r, λ, χ, f
#  	p = []
# 	t = time

	r, λ, χ, freq = u 						# unpack u
	#freq = p[1] 	# unpack p

	dμdr = ddr(r,λ,χ,freq)			# mutating form: ddr!([μ,dμdψ],r,λ,χ,freq,dμdr)
	dμdλ = ddλ(r,λ,χ,freq)			# calls phase_refractive_index(r,λ,χ,freq)
	dμdχ = ddχ(r,λ,χ,freq)
	dμdf = ddf(r,λ,χ,freq)

	v = phase_refractive_index(r,λ,χ,freq)
	μ = v[1]
	dμdψ = v[2]

	# 11/20: try dμdχ --> -1*dμdψ; these should be equal but are calculated differently
    du[1] = 1/(μ^2)*(μ*cos(χ)+dμdψ*sin(χ))
    du[2] = 1/(r*μ^2)*(μ*sin(χ)-dμdψ*cos(χ))
    du[3] = 1/(r*μ^2)*(dμdλ*cos(χ)-(r*dμdr + μ)*sin(χ))
	du[4] = 1/c*(1+(freq/μ)*dμdf)

	#du = [drdt, dλdt, dχdt, dfdt]

	#println("drdt = ", drdt)
	#println("dλdt = ", dθdt)
	#println("dχdt = ", dχdt)

	#DEBUG check values
	#println("drdt = ",drdt)
	#println("dλdt = ",dθdt)
	#println("dχdt = ",dχdt)

    # dμdψ = 1/(2*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B+F)/(2*A^2))
    # dμdω = 1/(2*μ)*((dBdω + dFdω)/(2*A) - 2*dAdω*(B+F)/(2*A^2))
	#
    # dμdr = (μ*(r+δr/2) - μ*(r-δr/2))/δr
    # dμdθ = (μ*(θ+δθ/2) - μ*(θ-δθ/2))/δθ

end

#DEBUG check values
#haselgrove!(1.0,1.0,10000.0,pi/4.0,0.0,5000.0,0.0,[1.0,1.0,1.0,1.0],[1.0,1.0,1.0])

## callbacks!
# re_cb: terminate when ray intersects Earth surface
function re_term_condition(u,t,integrator)
	u[1] - re
end

# function μ2_term_condition(u,t,integrator)
# 	body
# end

function terminate_affect!(integrator)
	terminate!(integrator)
end

re_cb = ContinuousCallback(re_term_condition,terminate_affect!)

# saveμ_cb: save μ, dμ values at each timestep
function save_func(u,t,integrator)
	r, λ, χ, f = u
	v = phase_refractive_index(r, λ, χ, f)
	dip = atan(2.0*tan(λ))
	ψ = π/2.0 + dip + χ
	vcat(v,dip,ψ)#,r,λ)
end

saved_μ = SavedValues(Float64, Array{Float64,1})
saveμ_cb = SavingCallback(save_func, saved_μ)

cb = CallbackSet(re_cb, saveμ_cb)
#cb = CallbackSet(saveμ_cb)

## Plot results
u0 = [re+1.0e+6, 1.0*pi/4, 0.0*pi/4,  1000.0]					# r0, λ0, χ0, f0
p = []	# f0, dμdψ, dμdr, dμdθ, dμdχ, dμdf
tspan = (0.0,5.0e+9)

hasel_prob = ODEProblem(haselgrove!,u0,tspan,p)
hasel_soln = solve(hasel_prob, CVODE_BDF(), reltol=1e-7, callback=cb, dtmax = 1e6, dtmin = 1e-8)

using Plots
plotly()
plot(hasel_soln)
plot(hasel_soln, vars=(1,2,3))

t = hasel_soln.t
r = hasel_soln[1,:]
λ = hasel_soln[2,:]
χ = hasel_soln[3,:]
f = hasel_soln[4,:]

x = r.*cos.(λ)
y = r.*sin.(λ)

plot(re.*cos.([0:0.01:2*π;]),re.*sin.([0:0.01:2*π;]), aspect_ratio = 1, label="Earth")
plot!(x,y, aspect_ratio=1, xlim = [-2e+7,2e+7], ylim = [-2e+7,2e+7], label=string(u0[end],"Hz"))

## saved values
t_s = saved_μ.t
val_s = saved_μ.saveval
val_s_r = reduce(hcat, val_s)
val_s_r = val_s_r'
#r_s = val_s_r[:,5]
#λ_s = val_s_r[:,6]
μ_s = val_s_r[:,1]
dμdψ_s = val_s_r[:,2]
dip_s = val_s_r[:,3]
ψ_s = val_s_r[:,4]

dt_s = zeros(length(t_s))
for i = 1:length(t_s)-1
	dt_s[i] = t_s[i+1]-t_s[i]
end

plot(t_s, μ_s, title="refractive index v. time")
plot(t_s, dμdψ_s, title="d/dpsi refractive index v. time")
plot(t_s, dip_s.*(180/pi), title="dip angle v. time")
plot(t_s, ψ_s.*(180/pi), title="wave normal angle v. time")
plot(1:1:length(dt_s), dt_s, title="time step size v. step number")

## refractive index surface
# initial conditions
x_test = 1*9.031E+6
y_test = 1*1.537E+6
λ_test = atan(abs(y_test/x_test))
r_test = sqrt(x_test^2 + y_test^2) #2.0*re*(sin(θ_test))^2 # on plotted B field line
#θ_test = θ_pw
#r_test = r_pw
ψ_test = [0:0.001:2π;]
dip_test = atan(2.0*tan(λ_test))
#χ_test = ψ_test .+ 3π/2 .- dip_test
χ_test = -1.0.*ψ_test .+ 3.0*π/2 .- dip_test
f_test = 5000.0

# calculate μ, dμdψ
u_test = phase_refractive_index.(r_test, λ_test, χ_test, f_test)
u_test_r = reduce(hcat,u_test)
u_test_r = u_test_r'
μ_test =  u_test_r[:,1]
dμdψ_test = u_test_r[:,2]

# xμ_test = μ_test.*sin.(ψ_test .+ π/2)
# yμ_test = μ_test.*cos.(ψ_test .+ π/2)
# # plot surface rotated into B|| frame
xμ_test_B = μ_test.*sin.(ψ_test)
yμ_test_B = μ_test.*cos.(ψ_test)

plot(xμ_test_B, yμ_test_B, aspect_ratio = 1)

# plot surface in x, y frame
xμ_test_xy = μ_test.*sin.(χ_test .- (π/2 - λ_test))
yμ_test_xy = μ_test.*cos.(χ_test .- (π/2 - λ_test))

plot(xμ_test_xy, yμ_test_xy, aspect_ratio = 1)

xdμ_test = dμdψ_test.*sin.(ψ_test)
ydμ_test = dμdψ_test.*cos.(ψ_test)

max_stars_repo_name: andersontodds/ray-trace

max_issues_repo_name: andersontodds/ray-trace

max_forks_repo_name: andersontodds/ray-trace


text:
__precompile__()

"""
Independent Component Analysis for Convolutive Mixtures
"""
module CICA

#using MultivariateStats
using Distributions
import JuMIT.Grid:M1D
import JuMIT.Misfits
import JuMIT.Data
import JuMIT.Grid
import JuMIT.DSP
import JuMIT.Acquisition

type ICA{T}
	"blended data at receivers(nr, nt)"
	x::Array{T,2}
	"blended data after whitening (ns, nt)"
	xhat::Array{T,2}
	"deblended data (ns, nt)"
	s::Array{T,2}
	"number of sources or independent components"
	ns::Int64
	"number of receivers"
	nr::Int64
	"number of samples"
	nt::Int64
	"store G(Wx)"
	GWx::Matrix{T}
	"store g(Wx)"
	gWx::Matrix{T}
	"store dg(Wx)"
	dgWx::Matrix{T}
	"whitening matrix --- this is updated after preprocessing"
	Q::Array{T,3}
	"unmixing matrix --- this is updated after ICA"
	W::Array{T,3}
	"previous unmixing matrix"
	Wp::Array{T,3}
	"store mean vector"
	mv::Array{T,2}
	"real or complex ica"
	attrib::Symbol
	"maximum iterations"
	maxiter::Int64
	"receiver at which deblending happens"
	magic_recv::Int64
	"tolerance"
	tol::Float64
	"bins"
	bins::Vector{UnitRange{Int64}}
end

function ICA(dat::Data.TD;
	     ica_func::Function=(x,ns)->ICA(x,ns),)
	Acquisition.Geom_isfixed(dat.acqgeom) || error("implemented only for fixed spread geom")
	acqgeom=dat.acqgeom
	ns=acqgeom.nss
	nr=acqgeom.nr[1] # fixed spread
	nt=dat.tgrid.nx
	xx=zeros(nt, nr, ns)
	for iss in 1:ns
		for ir in 1:nr
			for it in 1:nt
				xx[it,ir,iss] = dat.d[iss,1][it,ir]
			end
		end
	end
	# real fft gves only +ve freqs
	xx=rfft(xx, [1])

	xstack=complex.(zeros(size(xx,2),size(xx,1)))
	A=zeros(nr, ns)
	for iss in 1:ns
		for ir in 1:nr
			for i in 1:size(xx,1)
				# stack over supersources
				xstack[ir,i] += xx[i,ir,iss]
			end
			# mixing matrix (wrong!!!)
			A[ir,iss]=var(xx[:,ir,iss])
		end
	end

	ica = ica_func(xstack, ns), complex.(A)
	ica=ica[1]

	fastica!(ica)

	s=transpose(ica.s)

	xx=irfft(s, nt, [1])

	return xx

end
"""
`ns`	:
`dat`	: dat

"""
function ICA(x,
	     ns;
	     maxiter::Int=1000,
	     tol::Real=1e-10,                 # convergence tolerance
	     Win=nothing,  # init guess of W, size (m, k)
	     nbins=1, # number of bins
	     )

	T = eltype(x)
	nr, nt = size(x)

	# check input arguments
	(nt < 1) && error("there must be more than one samples, i.e. nt > 1.")
	(ns > nr) && error("ns must not exceed nr")

	(maxiter < 1) && error("maxiter must be greater than 1.")
	(tol < 0.) && error("tol must be positive.")

	# bins
	bindices=[round(Int, s) for s in linspace(0,nt,nbins+1)]
	bins=Array{UnitRange{Int64}}(nbins)
	for ib in 1:nbins
		bins[ib]=bindices[ib]+1:bindices[ib+1]
	end

	println("FastICA\t")
	println("===================")
	println("total samples\t",nt)
	println("total number of bins\t",nbins)
	println("average samples in each bin\t",round(Int,mean(length.(bins))))
	println("")

	# allocating vectors, that store G, g, and dg
	GWx = Matrix{T}(ns, nt)
	gWx = Matrix{T}(ns, nt)
	dgWx = Matrix{T}(ns, nt)

	# random initialization of unmixing matrix
	W = Array{T}(ns,ns, nbins)
	Wp = Array{T}(ns,ns, nbins)

	Q = Array{T}(nr, ns, nbins)
	xhat = Matrix{T}(ns, nt)

	mv = Matrix{T}(nr, nbins)
	# deblended data
	s = Matrix{T}(ns, nt)

	(T == Float64) ? attrib=:real : attrib=:complex

	ica=ICA(x, xhat, s, ns, nr, nt, GWx, gWx, dgWx, Q, W, Wp, mv, attrib,
	 maxiter, 1, tol, bins)

	unmixing_initialize!(ica, Win)

	return ica
end

function unmixing_initialize!(ica::ICA, Win=nothing)
	T = eltype(ica.x)
	ns=ica.ns
	W=ica.W
	Wp=ica.Wp
	bins=ica.bins

	if(!(Win===nothing))
		copy!(W, Win)
	else
		if(ica.attrib==:real)
			for ib in 1:length(bins), is1=1:ns, is2=1:ns,
				W[is1,is2,ib]=randn()
			end
		elseif(ica.attrib==:complex)
			for ib in 1:length(bins), is1=1:ns, is2=1:ns,
				W[is1,is2,ib]=complex.(randn()+randn())
			end
		else
			error("unknown attrib")
		end
	end

	# normalize each column of W necessary if arbitary Win is given
	for ib in 1:length(bins)
		for is = 1:ns
			w = view(W,:,is,ib)
			scale!(w, 1.0 / sqrt(sum(abs2, w)))
		end
	end
	return ica
end

"""
* update mv by storing mean of each component of x
* remove mean from each component of x
"""
function remove_mean!(ica::ICA)
	for ib in 1:length(ica.bins)
		for ir=1:ica.nr
			ica.mv[ir, ib]=zero(eltype(ica.mv))
			# compute mean
			for it=ica.bins[ib]
				ica.mv[ir,ib] += ica.x[ir,it]
			end
			ica.mv[ir, ib] *= inv(length(ica.bins[ib]))
			# remove mean
			for it=ica.bins[ib]
				ica.x[ir,it] -= ica.mv[ir,ib]

			end
		end
	end
end

"""
add stored mean to x
"""
function add_mean!(ica::ICA)
	for ib in 1:length(ica.bins)
		for ir=1:ica.nr
			# add mean
			for it=ica.bins[ib]
				ica.x[ir, it] += ica.mv[ir, ib]
			end
			ica.mv[ir, ib] = zero(eltype(ica.mv))
		end
	end
end

"""
Apply deblending matrix W to preprocessed data xhat and return deblended data ica.s
"""
function deblend!(ica::ICA)
	for ib in 1:length(ica.bins)
		s=view(ica.s,:,ica.bins[ib])
		W=view(ica.W,:,:,ib)
		xhat=view(ica.xhat,:,ica.bins[ib])
		Ac_mul_B!(s, W, xhat)
	end
	return ica.s
end

"""
update xhat
"""
function preprocess!(ica::ICA)
	# remove and store mean
	remove_mean!(ica)

	do_whiten!(ica)

	# adding mean back to x; avoided making a copy of x
	add_mean!(ica)
end

# Fast ICA
#
# Reference:
#
#   Aapo Hyvarinen and Erkki Oja
#   Independent Component Analysis: Algorithms and Applications.
#   Neural Network 13(4-5), 2000.
#
function fastica!(ica::ICA; verbose::Bool=false, A=nothing)

	preprocess!(ica)

	ns=ica.ns
	nt=ica.nt
	bins=ica.bins
	nbins=length(ica.bins)

	# aliases for x, W, Wp
	x=ica.xhat
	s=ica.s
	W=ica.W
	Q=ica.Q
	Wp=ica.Wp

	# some other aliases
	maxiter = ica.maxiter

	# aliases GWx, gWx, dgWx
	GWx = ica.GWx;	gWx = ica.gWx;	dgWx = ica.dgWx

	# main loop
	t=0
	converged=false
	while !converged && t < maxiter
		t += 1
		copy!(Wp, W)
		deblend!(ica)

		ICAfunc!(s, GWx, gWx, dgWx, ica.attrib)

		updateW!(W,gWx,dgWx,s,x,bins,ica.attrib)

		symmetric_decorrelation!(W)

		converged=compareWs(W,Wp,ica.tol)

		# compute the distance between A and W'Q
		# this gives the measure of the success of the algorithm
		# note that W'Q and A can still be different by a
		# permutation and scaling matrix, so...
		# if(!(A === nothing))
		#	SSE[t] = error_scale_perm(W'*Q', A)
		# end
	end

	# store expectation of G -- it measures Gaussianity
	EG = zeros(ns, nbins)

	# storing the expectation of G(Wx) --- measure of Gaussianity
	for ib=1:nbins
		for it in bins[ib]
			for ic=1:ns
				EG[ic,ib] += (GWx[ic, it])
			end
		end
		for ic=1:ns
			EG[ic,ib] /= length(bins[ib])
		end
	end

	fix_permutation!(W, EG)

	fix_scaling!(ica)

	# store error in the unmixing matrix -- only for testing (commented)
	SSE = (!(A === nothing)) ? error_unmixing(ica, A) : zeros(nbins)

	# perform final deblending after fixing permutation and scaling
	deblend!(ica)

	# test the performance of the algorithm
	return ica, SSE, EG
end

function error_unmixing(ica, A)
	nbins=length(ica.bins)
	Q=ica.Q
	W=ica.W
	SSE=zeros(nbins)
	for ib=1:nbins
		QQ=view(Q,:,:,ib)
		WW=view(W,:,:,ib)
		SSE[ib] = error_scale_perm(WW'*QQ', A)
	end
	return SSE
end

"""
compute the distance between s and
* `s` actual source matrix [nt, ns]
* `ica` ica output s is used

There is still global permutation
and scaling to be fixed.
"""
function error_after_ica(ica, s)
	(size(s) ≠ size(ica.s)) && error("dimension")

	errr=zeros(size(s,1))

	for is in 1:2
		ss0=view(s,is,:)
		err=zeros(size(s,1))
		for is1 in 1:2
			ss=view(ica.s,is1,:)
			err[is1] = Misfits.error_after_scaling(ss,ss0)[1]
		end
		errr[is] = minimum(err)
	end
	return errr
end

"""
Minimum Distortion Principle.
This method is not yet memory optimized.
Performs ICA for convolutive mixtures.
# Arguments
# Output
* deblended data at `magic_recv`
"""
function fix_scaling!(ica)
	nbins=length(ica.bins)
	ns=ica.ns
	W=ica.W
	Q=ica.Q
	WQ=zeros(eltype(W),size(W,1),size(W,2))
	for ib in 1:nbins
		WW=view(W,:,:,ib)
		QQ=view(Q,:,:,ib)

		Ac_mul_Bc!(WQ, WW, QQ)
		# A is the mixing matrix
		A = inv(WQ);

		# normalize mixing matrix, w.r.t. magic_recv
		A[:,1] /= A[ica.magic_recv,1];
		A[:,2] /= A[ica.magic_recv,2];

	 	copy!(WQ, (inv(A)));

		# scalings are only put in WW, so using inv(QQ)
		copy!(WW, inv(QQ)*WQ')
	end
end

"""
Order the components according to the EG values.
EG measures Gaussianity and hence, we are ordering sources
in each bin according to their distance from Gaussian.
"""
function fix_permutation!(W, EG)
	nbins=size(W,3)
	for ib=1:nbins
		WW=view(W,:,:,ib)
		# permutation order
		ord = sortperm(EG[:,ib]; rev=true)
		copy!(WW, WW[:,ord])
	end
end

function updateW!(W,gWx,dgWx,s,x,bins,attrib)
	nbins=length(bins)
	ns=size(W,2)
	if(attrib == :complex)

		# loop over bins; W for each bin is updated
		for ib=1:nbins
			for ic=1:ns
				b=zero(eltype(W))
				@simd for it in bins[ib]
					b+=gWx[ic,it]+(abs.(s[ic,it]).^2*dgWx[ic,it])
				end
				b /= length(bins[ib])

				for ic1=1:ns
					a=zero(eltype(x))
					@simd for it in bins[ib]
						a+=x[ic1,it]*conj(s[ic,it])*gWx[ic,it]
					end
					a /= length(bins[ib])
					W[ic1, ic, ib]=a-b*W[ic1,ic, ib]
				end
			end
		end
	elseif(attrib == :real)
		# loop over bins; W for each bin is updated
		for ib=1:nbins
			for ic = 1:ns
				b=zero(eltype(x))
				@simd for it=bins[ib]
					b+=dgWx[ic,it]
				end
				b /= length(bins[ib])
				for ic1=1:ns
					a=zero(eltype(x))
					@simd for it=bins[ib]
						a+=x[ic1,it]*gWx[ic,it]
					end
					a /= length(bins[ib])
					W[ic1, ic, ib]=a-b*W[ic1,ic, ib]
				end
			end
		end
	else
		error("invalid attrib")
	end
end

function ICAfunc!(s, GWx, gWx, dgWx, attrib)

	if(attrib==:complex)
		# for complex ICA (they used this constant in the paper)
		eps = 0.1
		@simd for it=1:size(s,2)
			for is=1:size(s,1)
				ss=eps+abs(s[is, it])^2
				GWx[is, it] = log(ss)
				gWx[is, it] = inv(ss)
				dgWx[is, it] = -1. * inv(ss*ss)
			end
		end
	elseif(attrib==:real)
		@simd for it=1:size(s,2)
			for is=1:size(s,1)
				GWx[is,it]=exp(-0.5*s[is,it]*s[is,it])
				gWx[is,it]=s[is,it]*GWx[is,it]
				dgWx[is,it]=GWx[is,it]*(1.-(s[is,it]*s[is,it]))
			end
		end
	end
end

# compare W with Wp as a convergence criteria
function compareWs(W,Wp,tol)
	nbins=size(W,3)
	ns=size(W,2)
	chgall=zeros(nbins)
	for ib in 1:nbins
		for ic in 1:ns
			chg=0.0
			for ic1 in 1:ns
				chg += abs(W[ic1, ic, ib]-Wp[ic1, ic, ib])
			end
			chgall[ib] = max(chg, chgall[ib])
		end
	end
	return all(chgall .< tol)
end

function symmetric_decorrelation!(W)
	# Symmetric decorrelation
	for ib=1:size(W,3)
		WW=view(W,:,:,ib)
		copy!(WW, WW * sqrtm(inv(WW'*WW)))
	end
end

"""
Compute the error in the unmixing matrix.

* x is the unmixing matrix
* y is the original mixing matrix

Ideally, x should be inverse of y upto scaling and permutation.
}
"""
function error_scale_perm{T}(x::Matrix{T}, y::Matrix{T})
	K=abs.(x*y);
	nc=size(x,2)
	err=0.0
	for ic in 1:nc
		kk = K[ic,:]
		ii = indmax(abs.(kk))
		kk ./= kk[ii]
		err += sum(abs.(kk[1:ii-1])) + sum(abs.(kk[ii+1:end]))
	end
	return err
end

"""
Whiten the data that are input to ICA. Whietning matrix is stored in `Q`.
"""
function do_whiten!(ica::ICA)
    nr = ica.nr
    for ib in 1:length(ica.bins)
	    x=view(ica.x,:,ica.bins[ib])
	    xhat=view(ica.xhat,:,ica.bins[ib])
	    Q=view(ica.Q,:,:,ib)
	    C=Matrix{eltype(x)}(nr, nr)
	    ntI = 1/length(ica.bins[ib])

	    # construct covariance matrix
	    A_mul_Bc!(C, x, x)
	    scale!(C, ntI)

	    # eigenvalue decomposition of the covariance matrix
	    Efac = eigfact(C)
	    ord = sortperm(Efac.values; rev=true)
	    v, P = Efac.values[ord], Efac.vectors[:, ord]
	    # whitening matrix
	    copy!(Q, scale!(P, sqrt.(inv.((v)))))

	    # for testing --  this should be a diagonal matrix
	    #println(Q' * C * Q)
	    Ac_mul_B!(xhat, Q, x)
    end
    return ica.Q
end

end # module


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: digits_sum

compute(n::Int, p::Int)::Int = digits_sum(BigInteger(n) ^ p)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
            # partition
            aab = AABox2D(Point2D{T}(0,0), Point2D{T}(3, 2))
            xdiv = SVector{2, T}(2,1)
            ydiv = SVector{1, T}(1)
            aabs = partition(aab, xdiv, ydiv)
            @test size(aabs) == (3,2)
            @test aabs[1].minima ≈ Point2D{T}(0, 0)
            @test aabs[1].maxima ≈ Point2D{T}(1, 1)
            @test aabs[2].minima ≈ Point2D{T}(1, 0)
            @test aabs[2].maxima ≈ Point2D{T}(2, 1)
            @test aabs[3].minima ≈ Point2D{T}(2, 0)
            @test aabs[3].maxima ≈ Point2D{T}(3, 1)
            @test aabs[4].minima ≈ Point2D{T}(0, 1)
            @test aabs[4].maxima ≈ Point2D{T}(1, 2)
            @test aabs[5].minima ≈ Point2D{T}(1, 1)
            @test aabs[5].maxima ≈ Point2D{T}(2, 2)
            @test aabs[6].minima ≈ Point2D{T}(2, 1)
            @test aabs[6].maxima ≈ Point2D{T}(3, 2)

            xdiv = SVector{2, T}(2,1)
            ydiv = SVector{0, T}()
            aabs = partition(aab, xdiv, ydiv)
            @test size(aabs) == (3,1)
            @test aabs[1].minima ≈ Point2D{T}(0, 0)
            @test aabs[1].maxima ≈ Point2D{T}(1, 2)
            @test aabs[2].minima ≈ Point2D{T}(1, 0)
            @test aabs[2].maxima ≈ Point2D{T}(2, 2)
            @test aabs[3].minima ≈ Point2D{T}(2, 0)
            @test aabs[3].maxima ≈ Point2D{T}(3, 2)

            xdiv = SVector{0, T}()
            ydiv = SVector{1, T}(1)
            aabs = partition(aab, xdiv, ydiv)
            @test size(aabs) == (1,2)
            @test aabs[1].minima ≈ Point2D{T}(0, 0)
            @test aabs[1].maxima ≈ Point2D{T}(3, 1)
            @test aabs[2].minima ≈ Point2D{T}(0, 1)
            @test aabs[2].maxima ≈ Point2D{T}(3, 2)

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
export @logical_calculus, InferenceRule, name, premises, conclusion, rule_by_name, PropositionalCalculus, ⊢

# TODO: implement set theoretic abstractions so a proper representation of the "universe of discourse" can be explicitly
#       defined in certain rules of inference (like the Principle of Explosion)

# the reason for this redeclaration is because a custom implementation of AbstractSet
# will be necessary for StatementSet to perform certain optimizations
Statement = AbstractExpression
StatementSet = Set{Statement}

"""
    InferenceRule(name::String, premise::StatementSet, conclusion::Statement)

Define a rule which maps a set of statements to a logical conclusion. The name is cited each time
a rule is applied in a proof step.
"""
struct InferenceRule
    name::String
    premise::StatementSet
    conclusion::Statement
end
name(ir::InferenceRule) = ir.name
premises(ir::InferenceRule) = ir.premise
conclusion(ir::InferenceRule) = ir.conclusion

InferenceRule(name::String, premise::Tuple{Vararg{T} where T}, conclusion::Statement) = InferenceRule(name, StatementSet(premise), conclusion)
InferenceRule(name::String, premise::T, conclusion::Statement) where {T <: Statement} = InferenceRule(name, (premise,), conclusion)

"""
    rule_by_name(calculus::Set{InferenceRule}, name::String)

Searches for an inference rule in a logical calculus based upon a given name.

## Examples
```julia
rule_by_name(PropositionalCalculus, "Modus Ponens")
rule_by_name(PropositionalCalculus, "Double Negation Introduction")
```
"""
rule_by_name(calculus::Set{InferenceRule}, rule_name::String) = collect(calculus)[findfirst(x->name(x)==rule_name, collect(calculus))]

_leaf_symbols(::Symbol) = Set{Symbol}()
function _leaf_symbols(expr::Expr)
    args = expr.args[2:end]
    if all((x -> x isa Symbol).(args))
        return Set{Symbol}(args)
    end

    Set(Iterators.flatten(_leaf_symbols.(args)))
end

"""
Defines a logical calculus based on a set of inference rules which can be applied repeatedly to
a given set of statements (premises). See `PropositionalCalculus` for an example definition using
this macro.
"""
macro logical_calculus(var, defs)
    rules_list = defs.args[2].args

    if length(rules_list) % 3 != 0
        throw(ErrorException("Each inference rule requires a name, premises, and a conclusion"))
    end

    variables_list = Set{Symbol}()
    inference_expressions = Expr[]

    for i ∈ 1:3:length(rules_list)
        push!(inference_expressions, :(
            InferenceRule($(rules_list[i]), $(esc(rules_list[i+1])), $(esc(rules_list[i+2])))
        ))
        union!(variables_list, _leaf_symbols(rules_list[i+1]))
        union!(variables_list, _leaf_symbols(rules_list[i+2]))
    end

    variable_declarations = [:(
        $(esc(var_name)) = LogicalSymbol($(esc(Symbol))($(string(var_name))))
    ) for var_name ∈ variables_list]

    quote
        $(variable_declarations...)
        $(esc(var)) = Set{InferenceRule}(
            [$(inference_expressions...)]
        )
    end
end

"""
    A ⊢ B

Provability operator used exclusively in calculus definitions. A ⊢ B means that B can be proven using A as the premise. Typed with \\vdash.
"""
⊢ = LogicalOperation((a, b) -> true, :⊢, 2, false, false)

@logical_calculus PropositionalCalculus begin
    "Modus Ponens", (
        p → q,
        p
    ), q,
    "Modus Tollens", (
        p → q,
        ¬q
    ), ¬p,
    "Associative Conjunction", (
        (p ∧ q) ∧ r
    ), p ∧ (q ∧ r),
    "Associative Disjunction", (
        (p ∨ q) ∨ r
    ), p ∨ (q ∨ r),
    "Commutative Conjunction", (
        p ∧ q
    ), q ∧ p,
    "Commutative Disjunction", (
        p ∨ q
    ), q ∨ p,
    "Law of Biconditional Propositions", (
        p → q,
        q → p
    ), p ⟷ q,
    "Exportation", (
        (p ∧ q) → r
    ), p → (q → r),
    "Transposition", (
        p → q
    ), ¬q → ¬p,
    "Hypothetical Syllogism", (
        p → q,
        q → r
    ), p → r,
    "Material Implication", (
        p → q
    ), ¬p ∨ q,
    "Distributive Conjunction", (
        (p ∨ q) ∧ r
    ), (p ∧ r) ∨ (q ∧ r),
    "Distributive Disjunction", (
        (p ∧ q) ∨ r
    ), (p ∨ r) ∧ (q ∨ r),
    "Absorption", (
        p → q
    ), p → (p ∧ q),
    "Disjunctive Syllogism", (
        p ∨ q,
        ¬p
    ), q,
    "Addition", (
        p
    ), p ∨ q,
    "Simplification", (
        p ∧ q
    ), p,
    "Conjunction", (
        p,
        q
    ), p ∧ q,
    "Double Negation Introduction", (
        p
    ), ¬¬p,
    "Double Negation Elimination", (
        ¬¬p
    ), p,
    "Disjunctive Simplification", (
        p ∨ p
    ), p,
    "Resolution", (
        p ∨ q,
        ¬p ∨ r
    ), q ∨ r,
    "Disjunction Elimination", (
        p → q,
        r → q,
        p ∨ r
    ), q
end
"""
    PropositionalCalculus

Default calculus definition for zeroth-order (propositional) logic.
"""
PropositionalCalculus

@logical_calculus ExtendedPropositionalCalculus begin
    # Rules for negations
    "Negation Introduction", (
        ϕ ⊢ ψ,
        ϕ ⊢ ¬ψ
    ), ¬ϕ,
    "Negation Elimination", (
        ¬ϕ ⊢ ψ,
        ¬ϕ ⊢ ¬ψ
    ), ϕ,
    "Principle of Explosion", (
        ϕ,
       ¬ϕ,
    ),  ψ,
    "Double Negation Elimination", (
        ¬¬ϕ
    ), ϕ,
    "Double Negation Introduction", (
        ϕ
    ), ¬¬ϕ,

    # Rules for conditionals
    "Deduction Theorem", (
        ϕ ⊢ ψ
    ), ϕ → ψ,
    "Modus Ponens", (
        ϕ → ψ,
        ϕ
    ), ψ,
    "Modus Tollens", (
        ϕ → ψ,
        ¬ψ
    ), ¬ϕ,

    # Rules for conjunctions
    "Adjunction", (
        ϕ,
        ψ
    ), ϕ ∧ ψ,
    "Simplification (Left)", (
        ϕ ∧ ψ
    ), ϕ,
    "Simplification (Right)", (
        ϕ ∧ ψ
    ), ψ,

    # Rules for disjunctions
    "Addition", (
        ϕ
    ), ϕ ∨ ψ,
    "Case Analysis", (
        ϕ → χ,
        ψ → χ,
        ϕ ∨ ψ
    ), χ,
    "Disjunctive Syllogism (Right)", (
        ϕ ∨ ψ,
        ¬ϕ
    ), ψ,
    "Disjunctive Syllogism (Left)", (
        ϕ ∨ ψ,
        ¬ψ
    ), ϕ,
    "Constructive Dilemma", (
        ϕ → χ,
        ψ → ξ,
        ϕ ∨ ψ
    ), χ ∨ ξ,

    # Rules for biconditionals
    "Biconditional Introduction", (
        ϕ → ψ,
        ψ → ϕ
    ), ϕ ⟷ ψ,
    "Biconditional Elimination (Asserted Left)", (
        ϕ ⟷ ψ,
        ϕ
    ), ψ,
    "Biconditional Elimination (Asserted Right)", (
        ϕ ⟷ ψ,
        ψ
    ), ϕ,
    "Biconditional Elimination (Negated Left)", (
        ϕ ⟷ ψ,
        ¬ϕ
    ), ¬ψ,
    "Biconditional Elimination (Negated Left)", (
        ϕ ⟷ ψ,
        ¬ψ
    ), ¬ϕ,
    "Biconditional Elimination (Asserted Disjunction)", (
        ϕ ⟷ ψ,
        ϕ ∨ ψ
    ), ϕ ∧ ψ,
    "Biconditional Elimination (Asserted Disjunction)", (
        ϕ ⟷ ψ,
        ¬ϕ ∨ ¬ψ
    ), ¬ϕ ∧ ¬ψ
end


text:
"""
Sort the original vector as well as an auxillary index vector in a least
significant digit (lsd) fashion. Good for sortperm implemenations
"""
function sorttwo_lsd16!(vs::AbstractVector{T}, index::AbstractVector{S}) where {T,S}
    l = length(vs)
    mask = mask16bit(T)

    ts = similar(vs)
    index1 = similar(index)

    # Init
    iters = sizeof(T) >> 1
    bin = zeros(UInt32, 65536, iters)

    # a function to compute numeric value of two bytes to increment the count with
    function getidx(vsi, j)
        Int(Base.bswap(UInt16((vsi >> (j-1)*16) & mask))) + 1
    end

    # Histogram for each element, radix
    for i = 1:l
        for j = 1:iters
            # idx = Int((vs[i] >> (j-1)*16) & 0xffff) + 1
            idx = getidx(vs[i], j)
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    swaps = 0
    for j = iters:-1:1
        # Unroll first data iteration, check for degenerate case
        idx = getidx(vs[l], j)

        # are all values the same at this radix?
        if bin[idx,j] == l;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]
        ts[ci] = vs[l]
        index1[ci] = index[l]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in l-1:-1:1
            idx = getidx(vs[i], j)
            ci = cbin[idx]
            ts[ci] = vs[i]
            index1[ci] = index[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        index, index1 = index1, index
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        index, index1 = index1, index
        for i = 1:l
            @inbounds vs[i] = ts[i]
            @inbounds index[i] = index1[i]
        end
    end
    (vs, index)
end


text:
using Primes, Printf

function addmsdigit(p::Integer, b::Integer, s::Integer)
    a = Vector{typeof(p)}()
    q = p
    for i in 1:(b-1)
        q += s
        isprime(q) || continue
        push!(a, q)
    end
    return a
end

function lefttruncprime(pbase::Integer)
    a = Vector{BigInt}()
    append!(a, primes(pbase - 1))
    mlt = zero(BigInt)
    s = one(BigInt)
    while !isempty(a)
        mlt = maximum(a)
        s *= pbase
        for i in 1:length(a)
            p = popfirst!(a)
            append!(a, addmsdigit(p, pbase, s))
        end
    end
    return mlt
end

lo, hi = 3, 17
println("The largest left truncatable primes for bases", @sprintf(" %d to %d.", lo, hi))
for i in lo:hi
    mlt = lefttruncprime(i)
    @printf("%10d %-30d (%s)\n", i, mlt, string(mlt, base=i))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
This module provides:
- An extension of the `AbstractArray` interface in order to properly deal with mutable caches.
- A mechanism to generate lazy arrays resulting from operations between arrays.
- A collection of concrete implementations of `AbstractArray`.

The exported names in this module are:

$(EXPORTS)
"""
module Arrays

using Gridap.Helpers
using Gridap.Algebra

using DocStringExtensions
using Test
using FillArrays
using LinearAlgebra
using Base: @propagate_inbounds
using ForwardDiff

import Base: size
import Base: getindex, setindex!
import Base: similar
import Base: IndexStyle
import Base: unaliascopy

# CachedArray

export CachedArray
export CachedMatrix
export CachedVector
export setsize!

# Map

export Map
export return_cache
export evaluate!
export evaluate
export test_map
export return_type
export return_value
export testargs
export inverse_map

export Broadcasting

export Operation

# LazyArray

export LazyArray
export array_cache
export getindex!
export lazy_map

export array_cache
export getindex!
export testitem
export testvalue
# export uses_hash
export test_array
export get_array
# # export add_to_array!

export CompressedArray

export Reindex
#export reindex

export PosNegReindex
export PosNegPartition
#export posneg_reindex

export FilterMap

export MulAddMap

export Table
export identity_table
export empty_table
export rewind_ptrs!
export length_to_ptrs!
export append_ptrs
export append_ptrs!
export get_ptrs_eltype
export get_data_eltype
export generate_data_and_ptrs
export find_inverse_index_map
export find_inverse_index_map!
export append_tables_globally
export append_tables_locally
export flatten_partition
export collect1d
export UNSET
export get_local_item
export find_local_index

export IdentityVector

export SubVector
export pair_arrays
export unpair_arrays

export AppendedArray
export lazy_append
export lazy_split

export autodiff_array_gradient
export autodiff_array_jacobian
export autodiff_array_hessian

export VectorWithEntryRemoved
export VectorWithEntryInserted

import AbstractTrees
export TreeNode
export print_op_tree
export similar_tree_node

import Gridap.Io: to_dict
import Gridap.Io: from_dict

export AddEntriesMap
export TouchEntriesMap

export ∑

include("Interface.jl")

include("Maps.jl")

include("AlgebraMaps.jl")

include("CachedArrays.jl")

include("LazyArrays.jl")

include("CompressedArrays.jl")

include("Tables.jl")

include("IdentityVectors.jl")

include("PosNegReindex.jl")

include("Reindex.jl")

include("FilteredArrays.jl")

include("SubVectors.jl")

include("ArrayPairs.jl")

include("AppendedArrays.jl")

include("Autodiff.jl")

include("VectorsWithEntryRemoved.jl")

include("VectorsWithEntryInserted.jl")

include("PrintOpTrees.jl")

const ∑ = sum

end # module


text:
export mutate

# This version added for compatibility, the mutation rate and the
# functions are taken from the parameters which are passed as a
# component of the chromosome.
function mutate(old_c::Chromosome)
    return mutate(old_c, old_c.params.funcs)
end

# Mutate Chromosome old_c.
function mutate(old_c::Chromosome, funcs::Vector{Func} )
    p = old_c.params
    mutrate = p.mutrate
    numinputs = p.numinputs
    numoutputs = p.numoutputs

    # Count number of genes
    num_genes = numoutputs # one connection per output gene
    for level = 1:p.numlevels
        for index = 1:p.numperlevel
            num_genes += 1+old_c.interiors[level,index].func.arity
        end
    end

    # Determine number of mutations, but make sure that it is always at least 1
    num_mutations = max(1,convert(Int,floor(num_genes*mutrate))) # doesn't seem right, but what CGP version 1.1 does)
    if num_mutations >= num_genes
        error("more mutations than genes in function mutate")
    end

    new_c = Chromosome(p)
    new_c.active_set = false

    # Choose the genes that will be modified, and store their numbers in the array genes_to_mutate
    genes_to_mutate = Array(Int,num_mutations)
    chosen_genes = Set()
    for i in 1:num_mutations
        gene = rand(1:num_genes)
        while gene in chosen_genes
            gene = rand(1:num_genes)
        end
        push!(chosen_genes,gene)
        genes_to_mutate[i] = gene
    end
    sort!(genes_to_mutate)

    # Input nodes
    for i = 1:numinputs
        new_c.inputs[i] = deepcopy(old_c.inputs[i])
        new_c.inputs[i].active = false
    end

    gene_index = 1
    genes_to_mutate_index = 1
    # Interior nodes
    new_c.interiors =  Array(InteriorNode, p.numlevels, p.numperlevel) # Always use a new array
    for level = 1:p.numlevels
        minlevel = max(level - p.numlevelsback, 0)
        maxlevel = level - 1
        for index = 1:p.numperlevel
            old_func = old_c[level, index].func
            # check if none of the genes in node will be mutated
            if genes_to_mutate_index > length(genes_to_mutate) || gene_index + old_func.arity < genes_to_mutate[genes_to_mutate_index]
                # no mutations in this node
                new_c.interiors[level, index] = deepcopy(old_c.interiors[level, index])
                new_c.interiors[level, index].active = false
                gene_index += 1+old_func.arity
                continue
            end
            new_inputs = Array(NodePosition, old_func.arity)
            new_func = old_func
            for i in 0:old_func.arity # iterate through the genes in this node with 0 corresponding to the function
                if genes_to_mutate_index <= length(genes_to_mutate) && gene_index == genes_to_mutate[genes_to_mutate_index]  # mutate
                    genes_to_mutate_index += 1
                    if i == 0    # mutate the function
                        new_func = funcs[rand(1:end)]
                        if new_func.arity != old_func.arity # change in arity
                            new_inputs = Array(NodePosition, new_func.arity) # we must generate a new node with the new arity
                            for j = 1:new_func.arity
                                if j <= old_func.arity
                                    new_inputs[j] = old_c[level, index].inputs[j]
                                else
                                    new_inputs[j] = random_node_position(p, minlevel, maxlevel)
                                end
                            end
                            gene_index += 1+old_func.arity # skip over mutating inputs
                            # skip over mutating inputs
                            while genes_to_mutate_index <= length(genes_to_mutate) &&  genes_to_mutate[genes_to_mutate_index] < gene_index
                                genes_to_mutate_index += 1
                            end
                            break # break out of for i loop
                        end
                    else  # mutate the connection
                        new_inputs[i] = random_node_position(p, minlevel, maxlevel)
                    end
                else # don't mutate this gene
                    if i == 0
                        new_func = old_func
                    else
                        new_inputs[i] = old_c.interiors[level,index].inputs[i]
                    end
                end # if genes_to_mutate_index >= . . .
                gene_index += 1
            end # for i in 0:old_func.arity
            new_c.interiors[level, index] = InteriorNode(new_func, new_inputs)
        end # for index =
    end # for level =

    # Output nodes
    minlevel = max(p.numlevels - p.numlevelsback + 1, 0)
    maxlevel = p.numlevels
    for i = 1:numoutputs
        if  genes_to_mutate_index <= length(genes_to_mutate) && gene_index == genes_to_mutate[genes_to_mutate_index]  # mutate
            genes_to_mutate_index += 1
            (level, index) = random_node_position(p, minlevel, maxlevel)
            new_c.outputs[i] = OutputNode((level, index))
            new_c[level, index].active = false
        else
            new_c.outputs[i] = deepcopy(old_c.outputs[i])
            new_c.outputs[i].active = false
        end
        gene_index += 1
    end # for i =
    return new_c
end


text:
"""
    getsolution(x::BnBModel)

Returns the solution stored in the BnBModel.
"""
getsolution(x::BnBModel) = x.soln

"""
    getobjval(x::BnBModel)

Returns the objective value stored in BnBModel (global upper bound).
"""
getobjval(x::BnBModel) = x.UBDg

"""
    getobjbound(x::BnBModel)

Returns the objective value stored in BnBModel (global upper bound).
"""
getobjbound(x::BnBModel) = x.UBDg

"""
    getfeasibility(x::BnBModel)

Returns feasibility of problem (feasible point found?).
"""
getfeasibility(x::BnBModel) = x.feas_fnd

"""
    LBDtime(x::BnBModel)

Returns time spent solving lower bounding problem.
"""
LBDtime(x::BnBModel) = x.LBDgtime[end]

"""
    UBDtime(x::BnBModel)

Returns time spent solving upper bounding problem.
"""
UBDtime(x::BnBModel) = x.UBDgtime[end]


text:
module UnivariateFunctions

using Dates
using SchumakerSpline
using GLM

# This includes timing functions for use of UnivariateFunctions with dates.
include("date_conversions.jl")
export years_between, years_from_global_base, period_length
# The abstrat type UnivariateFunction and all structs are implemented here.
# In addition we have operator reversals and some supporting functions.
include("0_structs_and_generic_reversals.jl")
export UnivariateFunction, Undefined_Function, PE_Function, Sum_Of_Functions, Piecewise_Function
export change_base_of_PE_Function, trim_piecewise_function, convert_to_linearly_rescale_inputs

# These all implement evaluation, calculus and operators for the main structs.
include("1_undefined_function.jl")
include("2_pe_functions.jl")
include("3_sum_of_functions.jl")
include("4_piecewise_functions.jl")
export evaluate, derivative, indefinite_integral

include("5_calculus.jl")
export  evaluate_integral, right_integral, left_integral

include("chebyshevs.jl")
export get_chevyshevs_up_to

include("6_splines_and_interpolation.jl")
export create_quadratic_spline, create_constant_interpolation_to_right
export create_constant_interpolation_to_left, create_linear_interpolation

include("7_regressions_and_approximation.jl")
export create_ols_approximation, create_chebyshev_approximation

# The following operators also have overloads
export +, -, *, /, ^
end

# In terms of UnivariateFunction types we have the following seniority:
# Undefined_Function
# PE_Function
# Sum_Of_Functions
# Piecewise_Function

# In particular a Sum_Of_Functions contains the first two types of functions but never
# another Sum_Of_Functions (although one can be used to construct a Sum_Of_Functions). If a sum
# of functions contains an undefined function then it contains nothing else.
# A Piecewise_Function can have pieces composed of the first three categories but not another
# Piecewise_Function

# In terms of multiple dispatch rules this means that we:
# * define operators of each type of UnivariateFunction with another of its own type.
# * For each we define operators (f + x) with scalars. The reverse operations (x + f) can be done
#    all at once except for powers, subtraction and division where it matters.
#  All operations between these types should be defined for every pair. The seniority rule means that
#  Where a * b (where a and b are types) is the same as b * a then we define b*a = *(b,a) and put the
# relevent code in b * a place. Undefined_Function should take prevedence over PE_Function and etc
# in the seniority list.


text:
#=
Utilities used by UMAP.jl
=#

@inline fit_ab(_, __, a, b) = a, b

"""
    fit_ab(min_dist, spread, _a, _b) -> a, b

Find a smooth approximation to the membership function of points embedded in ℜᵈ.
This fits a smooth curve that approximates an exponential decay offset by `min_dist`,
returning the parameters `(a, b)`.
"""
function fit_ab(min_dist, spread, ::Nothing, ::Nothing)
    ψ(d) = d >= min_dist ? exp(-(d - min_dist)/spread) : 1.
    xs = LinRange(0., spread*3, 300)
    ys = map(ψ, xs)
    @. curve(x, p) = (1. + p[1]*x^(2*p[2]))^(-1)
    result = curve_fit(curve, xs, ys, [1., 1.], lower=[0., -Inf])
    a, b = result.param
    return a, b
end

knn_search(dist_mat, k, metric::Symbol) = knn_search(dist_mat, k, Val(metric))

"""
    knn_search(dist_mat, k, :precomputed) -> knns, dists

Find the `k` nearest neighbors of each point in a precomputed distance
matrix.
"""
knn_search(dist_mat, k, ::Val{:precomputed}) = _knn_from_dists(dist_mat, k)

"""
    knn_search(X, k, metric) -> knns, dists

Find the `k` nearest neighbors of each point in `X` by `metric`.
"""
function knn_search(X,
                    k,
                    metric::SemiMetric)
    if size(X, 2) < 4096
        return knn_search(X, k, metric, Val(:pairwise))
    else
        return knn_search(X, k, metric, Val(:approximate))
    end
end

# compute all pairwise distances
# return the nearest k to each point v, other than v itself
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:pairwise}) where {S <: Real}
    num_points = size(X, 2)
    dist_mat = Array{S}(undef, num_points, num_points)
    pairwise!(dist_mat, metric, X, dims=2)
    # all_dists is symmetric distance matrix
    return _knn_from_dists(dist_mat, k)
end

# find the approximate k nearest neighbors using NNDescent
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:approximate}) where {S <: Real}
    knngraph = DescentGraph(X, k, metric)
    return knngraph.indices, knngraph.distances
end

function _knn_from_dists(dist_mat::AbstractMatrix{S}, k) where {S <: Real}
    knns_ = [partialsortperm(view(dist_mat, :, i), 2:k+1) for i in 1:size(dist_mat, 1)]
    dists_ = [dist_mat[:, i][knns_[i]] for i in eachindex(knns_)]
    knns = hcat(knns_...)::Matrix{Int}
    dists = hcat(dists_...)::Matrix{S}
    return knns, dists
end

# combine local fuzzy simplicial sets
@inline function combine_fuzzy_sets(fs_set::AbstractMatrix{T},
                                    set_op_ratio) where {T}
    return set_op_ratio .* fuzzy_set_union(fs_set) .+
           (one(T) - set_op_ratio) .* fuzzy_set_intersection(fs_set)
end

@inline function fuzzy_set_union(fs_set::AbstractMatrix)
    return fs_set .+ fs_set' .- (fs_set .* fs_set')
end

@inline function fuzzy_set_intersection(fs_set::AbstractMatrix)
    return fs_set .* fs_set'
end


text:
import Base: push!, eltype, close
export Signal, push!, value, preserve, unpreserve, close

##### Signal #####

const debug_memory = false # Set this to true to debug gc of nodes

const nodes = WeakKeyDict()
const io_lock = ReentrantLock()

if !debug_memory
    type Signal{T}
        value::T
        parents::Tuple
        actions::Vector
        alive::Bool
        preservers::Dict
    end
else
    type Signal{T}
        value::T
        parents::Tuple
        actions::Vector
        alive::Bool
        preservers::Dict
        bt
        function Signal(v, parents, actions, alive, pres)
            n=new(v,parents,actions,alive,pres,backtrace())
            nodes[n] = nothing
            finalizer(n, log_gc)
            n
        end
    end
end

log_gc(n) =
    @async begin
        lock(io_lock)
        print(STDERR, "Signal got gc'd. Creation backtrace:")
        Base.show_backtrace(STDERR, n.bt)
        println(STDOUT)
        unlock(io_lock)
    end

immutable Action
    recipient::WeakRef
    f::Function
end
isrequired(a::Action) = a.recipient.value != nothing && a.recipient.value.alive

Signal{T}(x::T, parents=()) = Signal{T}(x, parents, Action[], true, Dict{Signal, Int}())
Signal{T}(::Type{T}, x, parents=()) = Signal{T}(x, parents, Action[], true, Dict{Signal, Int}())
# A signal of types
Signal(t::Type) = Signal(Type, t)

# preserve/unpreserve nodes from gc
"""
    preserve(signal::Signal)

prevents `signal` from being garbage collected as long as any of its parents are around. Useful for when you want to do some side effects in a signal.
e.g. `preserve(map(println, x))` - this will continue to print updates to x, until x goes out of scope. `foreach` is a shorthand for `map` with `preserve`.
"""
function preserve(x::Signal)
    for p in x.parents
        p.preservers[x] = get(p.preservers, x, 0)+1
        preserve(p)
    end
    x
end

"""
    unpreserve(signal::Signal)

allow `signal` to be garbage collected. See also `preserve`.
"""
function unpreserve(x::Signal)
    for p in x.parents
        n = get(p.preservers, x, 0)-1
        if n <= 0
            delete!(p.preservers, x)
        else
            p.preservers[x] = n
        end
        unpreserve(p)
    end
    x
end

Base.show(io::IO, n::Signal) =
    write(io, "Signal{$(eltype(n))}($(n.value), nactions=$(length(n.actions))$(n.alive ? "" : ", closed"))")

value(n::Signal) = n.value
value(::Void) = false
eltype{T}(::Signal{T}) = T
eltype{T}(::Type{Signal{T}}) = T

##### Connections #####

function add_action!(f, node, recipient)
    a = Action(WeakRef(recipient), f)
    push!(node.actions, a)
    a
end

function remove_action!(f, node, recipient)
    node.actions = filter(a -> a.f != f, node.actions)
end

function close(n::Signal, warn_nonleaf=true)
    finalize(n) # stop timer etc.
    n.alive = false
    if !isempty(n.actions)
        any(map(isrequired, n.actions)) && warn_nonleaf &&
            warn("closing a non-leaf node is not a good idea")
        empty!(n.actions)
    end

    for p in n.parents
        delete!(p.preservers, n)
    end
end

function send_value!(node::Signal, x, timestep)
    # Dead node?
    !node.alive && return

    # Set the value and do actions
    node.value = x
    for action in node.actions
        do_action(action, timestep)
    end
end
send_value!(wr::WeakRef, x, timestep) = wr.value != nothing && send_value!(wr.value, x, timestep)

do_action(a::Action, timestep) =
    isrequired(a) && a.f(a.recipient.value, timestep)

# If any actions have been gc'd, remove them
cleanup_actions(node::Signal) =
    node.actions = filter(isrequired, node.actions)

##### Messaging #####

const CHANNEL_SIZE = 1024
abstract AbstractMessage
immutable StopMessage <: AbstractMessage
end
immutable EmptyMessage <: AbstractMessage
end
immutable Message <: AbstractMessage
    node
    value
    onerror::Function
end
# Global channel for signal updates
const _messages = Channel{AbstractMessage}(CHANNEL_SIZE)

"""
`push!(signal, value, onerror=Reactive.print_error)`

Queue an update to a signal. The update will be propagated when all currently
queued updates are done processing.

The third optional argument is a callback to be called in case the update
ends in an error. The callback receives 3 arguments: the signal, the value,
and a `CapturedException` with the fields `ex` which is the original exception
object, and `processed_bt` which is the backtrace of the exception.

The default error callback will print the error and backtrace to STDERR.
"""
Base.push!(n::Signal, x, onerror=print_error) = _push!(n, x, onerror)

function _push!(n, x, onerror=print_error)
    taken = Base.n_avail(_messages)
    if taken >= CHANNEL_SIZE
        warn("Message queue is full. Ordering may be incorrect.")
        @async put!(_messages, Message(n, x, onerror))
    else
        put!(_messages, Message(WeakRef(n), x, onerror))
    end
    nothing
end
_push!(::Void, x, onerror=print_error) = nothing

"""
Convenience function to notify the Reactive event loop without pushing to
a Signal.
"""
post_empty() = put!(_messages, EmptyMessage())

# remove messages from the channel and propagate them
global run, stop

let timestep::Int = 0, stop_runner::Bool = false

    function stop()
        run_till_now() # process all remaining events
        stop_runner = true
        post_empty() # make sure it's not waiting on an empty message list
    end

    """
    Processes `n` messages from the Reactive event queue.
    """
    function run(n::Int)
        for i=1:n
            timestep += 1
            let message = take!(_messages)
                isa(message, EmptyMessage) && continue # ignore emtpy messages
                try
                    send_value!(message.node, message.value, timestep)
                catch err
                    if isa(err, InterruptException)
                        println("Reactive event loop was inturrupted.")
                        rethrow()
                    else
                        bt = catch_backtrace()
                        message.onerror(message.node, message.value, CapturedException(err, bt))
                    end
                end
            end
        end
        return
    end
    """
    Processes messages as long as `stop` isn't called.
    """
    function run()
        stop_runner = false
        while !stop_runner
            run(1)
        end
    end
end

# Default error handler function
function print_error(node, value, ex)
    lock(io_lock)
    io = STDERR
    println(io, "Failed to push!")
    print(io, "    ")
    show(io, value)
    println(io)
    println(io, "to node")
    print(io, "    ")
    show(io, node)
    println(io)
    showerror(io, ex)
    println(io)
    unlock(io_lock)
end

# Run everything queued up till the instant of calling this function
run_till_now() = run(Base.n_avail(_messages))

# A decent default runner task
function __init__()
    global runner_task = @async begin
        Reactive.run()
    end
end


text:
@testset "859.buddy-strings.jl" begin
    A = "ab"
    B = "ba"
    @test buddy_strings(A, B) == true

    A = "ab"
    B = "ab"
    @test buddy_strings(A, B) == false

    A = "aa"
    B = "aa"
    @test buddy_strings(A, B) == true

    A = "aaaaaaabc"
    B = "aaaaaaacb"
    @test buddy_strings(A, B) == true

    A = ""
    B = "aa"
    @test buddy_strings(A, B) == false

    @test buddy_strings("abc", "bca") == false
end


text:
using GitHub
using Base.Test

event_request = open(joinpath(Pkg.dir("GitHub"), "test", "commit_comment.jls"), "r") do file
    return Base.deserialize(file)
end
event_json = Requests.json(event_request)
event = GitHub.event_from_payload!("commit_comment", event_json)

################
# WebhookEvent #
################

@test get(event.repository.name) == "BenchmarkTrackers.jl"
@test get(event.sender.login) == "jrevels"

#################
# EventListener #
#################

@test !(GitHub.has_valid_secret(event_request, "wrong"))
@test GitHub.has_valid_secret(event_request, "secret")
@test !(GitHub.is_valid_event(event_request, ["wrong"]))
@test GitHub.is_valid_event(event_request, ["commit_comment"])
@test !(GitHub.from_valid_repo(event, ["JuliaWeb/GitHub.jl"]))
@test GitHub.from_valid_repo(event, ["JuliaCI/BenchmarkTrackers.jl"])
@test GitHub.handle_event_request(event_request, x -> true,
                                  secret = "secret",
                                  events = ["commit_comment"],
                                  repos = ["JuliaCI/BenchmarkTrackers.jl"])

@test begin
    listener = EventListener(x -> true;
                             secret = "secret",
                             repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                             events = ["commit_comment"],
                             forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")])

    r = listener.server.http.handle(HttpCommon.Request(),HttpCommon.Response())
    r.status == 400
end

###################
# CommentListener #
###################

trigger_results = GitHub.extract_trigger_string(event, GitHub.AnonymousAuth(), "RunBenchmarks", false)

@test trigger_results == (true,"(\"binary\", \"unary\")")

@test begin
    listener = CommentListener((x, y) -> true, "trigger";
                               secret = "secret",
                               repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                               forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")],
                               check_collab = false)

    r = listener.listener.server.http.handle(HttpCommon.Request(),HttpCommon.Response())
    r.status == 400
end


text:
export load

using PyCall

# https://github.com/dhoegh/Hawk.jl/blob/master/src/Hawk.jl

filename = abspath(joinpath(dirname(@__FILE__), "kinematics.py"))

@pyimport imp
(path, name) = dirname(filename), basename(filename)
(name, ext) = rsplit(name, '.')

(file, filename, data) = imp.find_module(name, [path])
conv = imp.load_module(name, file, filename, data)

function load(action="walking_0"; sample=Pkg.dir("HumanMotionPrediction", "samples", "walking_10000.h5"))
    conv[:load](action, sample)
end


text:
using Test
using Distributions, Random

rng = MersenneTwister(123)

@testset "Truncated normal, mean and variance" begin
    @test mean(truncated(Normal(0,1),100,115)) ≈ 100.00999800099926070518490239457545847490332879043
    @test mean(truncated(Normal(-2,3),50,70)) ≈ 50.171943499898757645751683644632860837133138152489
    @test mean(truncated(Normal(0,2),-100,0)) ≈ -1.59576912160573071175978423973752747390343452465973
    @test mean(truncated(Normal(0,1))) == 0
    @test mean(truncated(Normal(0,1); lower=-Inf, upper=Inf)) == 0
    @test mean(truncated(Normal(0,1); lower=0)) ≈ +√(2/π)
    @test mean(truncated(Normal(0,1); lower=0, upper=Inf)) ≈ +√(2/π)
    @test mean(truncated(Normal(0,1); upper=0)) ≈ -√(2/π)
    @test mean(truncated(Normal(0,1); lower=-Inf, upper=0)) ≈ -√(2/π)
    @test var(truncated(Normal(0,1),50,70)) ≈ 0.00039904318680389954790992722653605933053648912703600
    @test var(truncated(Normal(-2,3); lower=50, upper=70)) ≈ 0.029373438107168350377591231295634273607812172191712
    @test var(truncated(Normal(0,1))) == 1
    @test var(truncated(Normal(0,1); lower=-Inf, upper=Inf)) == 1
    @test var(truncated(Normal(0,1); lower=0)) ≈ 1 - 2/π
    @test var(truncated(Normal(0,1); lower=0, upper=Inf)) ≈ 1 - 2/π
    @test var(truncated(Normal(0,1); upper=0)) ≈ 1 - 2/π
    @test var(truncated(Normal(0,1); lower=-Inf, upper=0)) ≈ 1 - 2/π
    # https://github.com/JuliaStats/Distributions.jl/issues/827
    @test mean(truncated(Normal(1000000,1),0,1000)) ≈ 999.999998998998999001005011019018990904720462367106
    @test var(truncated(Normal(),999000,1e6)) ≥ 0
    @test var(truncated(Normal(1000000,1),0,1000)) ≥ 0
    # https://github.com/JuliaStats/Distributions.jl/issues/624
    @test rand(truncated(Normal(+Inf, 1), 0, 1)) ≈ 1
    @test rand(truncated(Normal(-Inf, 1), 0, 1)) ≈ 0
end
@testset "Truncated normal $trunc" begin
    trunc = truncated(Normal(0, 1), -2, 2)
    @testset "Truncated normal $trunc with $func" for func in
        [(r = rand, r! = rand!),
         (r = ((d, n) -> rand(rng, d, n)), r! = ((d, X) -> rand!(rng, d, X)))]
        repeats = 1000000

        @test abs(mean(func.r(trunc, repeats))) < 0.01
        @test abs(median(func.r(trunc, repeats))) < 0.01
        @test abs(var(func.r(trunc, repeats)) - var(trunc)) < 0.01

        X = Matrix{Float64}(undef, 1000, 1000)
        func.r!(trunc, X)
        @test abs(mean(X)) < 0.01
        @test abs(median(X)) < 0.01
        @test abs(var(X) - var(trunc)) < 0.01
    end
end

@testset "Truncated normal should be numerically stable at low probability regions" begin
    original = Normal(-5.0, 0.2)
    trunc = truncated(original, 0.0, 5.0)
    for x in LinRange(0.0, 5.0, 100)
        @test isfinite(logpdf(original, x))
        @test isfinite(logpdf(trunc, x))
        @test isfinite(pdf(trunc, x))
    end
end


text:
@testset "area measurement" begin

    img = [0 0 0 0 0;
           0 1 0 1 0;
           0 1 1 0 0;
           0 0 0 0 0]

    labels  = Images.label_components(img,trues(3,3))
    t = measure_components((Measurement(),), labels)
    @test length(t) == 1

    row₁ = t[1]
    @test row₁.area == 4.375

    test_image_areas = ([0], [25], [1], [1, 1, 1, 1], [17], [16], [8.5],
                        [9.25], [4, 4], [21.25], [19.25], [20], [19.875],
                        [21.5], [17], [20.75, 8.5, 12.5], [38.5], [19.625, 8.75],
                        [25.5, 11], [31.5, 11.5])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((Measurement(),), labels)
            i == 1 ? (@test isempty(t)) : (@test all(select(t,:area) .== test_image_areas[i]))
        end
    end
end

@testset "perimeter" begin

    test_image_perimeters = ([0], [20], [4], [4, 4, 4, 4], [52], [48], [16],
                        [20], [10, 10], [28], [32], [34], [32],
                        [28], [30], [40, 16, 24], [76], [40, 24],
                        [72, 28], [64, 30])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((Measurement(),), labels)
            i == 1 ? (@test isempty(t)) : (@test all(select(t,:perimeter₁) .== test_image_perimeters[i]))
        end
    end
end


text:
"""
    get_string(x::VPtr, fn_sym)::Maybe{String}

C-call the SBML function `fn_sym` with a single parameter `x`, interpret the
result as a string and return it, or throw exception in case the pointer is
NULL.
"""
function get_string(x::VPtr, fn_sym)::String
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        throw(DomainError(x, "Calling $fn_sym returned NULL, valid string expected."))
    end
end

"""
    get_optional_string(x::VPtr, fn_sym)::Maybe{String}

Like [`get_string`](@ref), but returns `nothing` instead of throwing an
exception.

This is used to get notes and annotations and several other things (see
`get_notes`, `get_annotations`)
"""
function get_optional_string(x::VPtr, fn_sym)::Maybe{String}
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        return nothing
    end
end

"""
    get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}

Helper for getting out boolean flags.
"""
function get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cint, (VPtr,), x) != 0
    else
        return nothing
    end
end

"""
    get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{UInt}

Helper for getting out unsigned integers.
"""
function get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{Int}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cint, (VPtr,), x)
    else
        return nothing
    end
end

"""
    get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}

Helper for getting out C doubles aka Float64s.
"""
function get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cdouble, (VPtr,), x)
    else
        return nothing
    end
end

function get_string_from_xmlnode(xmlnode::VPtr)::String
    if ccall(sbml(:XMLNode_isText), Bool, (VPtr,), xmlnode)
        str_ptr = ccall(sbml(:XMLNode_getCharacters), Cstring, (VPtr,), xmlnode)
        str_ptr == C_NULL ? "" : unsafe_string(str_ptr)
    else
        children_num = ccall(sbml(:XMLNode_getNumChildren), Cuint, (VPtr,), xmlnode)
        join(
            (
                get_string_from_xmlnode(
                    ccall(sbml(:XMLNode_getChild), VPtr, (VPtr, Cint), xmlnode, n),
                ) for n = 0:(children_num-1)
            ),
            "\n",
        )
    end
end

function _readSBML(
    symbol::Symbol,
    fn::String,
    sbml_conversion,
    report_severities,
)::SBML.Model
    doc = ccall(sbml(symbol), VPtr, (Cstring,), fn)
    try
        get_error_messages(
            doc,
            AssertionError("Opening SBML document has reported errors"),
            report_severities,
        )

        sbml_conversion(doc)

        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
            throw(AssertionError("SBML document contains no model"))
        end

        model = ccall(sbml(:SBMLDocument_getModel), VPtr, (VPtr,), doc)

        return extract_model(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
    end
end

"""
    readSBML(
        fn::String,
        sbml_conversion = document -> nothing;
        report_severities = ["Fatal", "Error"],
    )::SBML.Model

Read the SBML from a XML file in `fn` and return the contained `SBML.Model`.

The `sbml_conversion` is a function that does an in-place modification of the
single parameter, which is the C pointer to the loaded SBML document (C type
`SBMLDocument*`). Several functions for doing that are prepared, including
[`set_level_and_version`](@ref), [`libsbml_convert`](@ref), and
[`convert_simplify_math`](@ref).

`report_severities` switches on and off reporting of certain errors; see the
documentation of [`get_error_messages`](@ref) for details.

To read from a string instead of a file, use [`readSBMLFromString`](@ref).

# Example
```
m = readSBML("my_model.xml", doc -> begin
    set_level_and_version(3, 1)(doc)
    convert_simplify_math(doc)
end)
```
"""
function readSBML(
    fn::String,
    sbml_conversion = document -> nothing;
    report_severities = ["Fatal", "Error"],
)::SBML.Model
    isfile(fn) || throw(AssertionError("$(fn) is not a file"))
    _readSBML(:readSBML, fn, sbml_conversion, report_severities)
end
"""
    readSBML(
        str::String,
        sbml_conversion = document -> nothing;
        report_severities = ["Fatal", "Error"],
    )::SBML.Model

Read the SBML from the string `str` and return the contained `SBML.Model`.

For the other arguments see the docstring of [`readSBML`](@ref), which can be
used to read from a file instead of a string.
"""
readSBMLFromString(
    str::AbstractString,
    sbml_conversion = document -> nothing;
    report_severities = ["Fatal", "Error"],
)::SBML.Model =
    _readSBML(:readSBMLFromString, String(str), sbml_conversion, report_severities)

get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)

"""
    function get_association(x::VPtr)::GeneProductAssociation

Convert a pointer to SBML `FbcAssociation_t` to the `GeneProductAssociation`
tree structure.
"""
function get_association(x::VPtr)::GeneProductAssociation
    # libsbml C API is currently missing functions to check this in a normal
    # way, so we use a bit of a hack.
    typecode = ccall(sbml(:SBase_getTypeCode), Cint, (VPtr,), x)
    if typecode == 808 # SBML_FBC_GENEPRODUCTREF
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
    elseif typecode == 809 # SBML_FBC_AND
        return GPAAnd([
            get_association(
                ccall(sbml(:FbcAnd_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcAnd_getNumAssociations), Cuint, (VPtr,), x)
        ])
    elseif typecode == 810 # SBML_FBC_OR
        return GPAOr([
            get_association(
                ccall(sbml(:FbcOr_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcOr_getNumAssociations), Cuint, (VPtr,), x)
        ])
    else
        throw(ErrorException("Unsupported FbcAssociation type"))
    end
end

extract_parameter(p::VPtr)::Pair{String,Tuple{Float64,String}} =
    get_string(p, :Parameter_getId) => (
        ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p),
        mayfirst(get_optional_string(p, :Parameter_getUnits), ""),
    )

""""
    function extract_model(mdl::VPtr)::SBML.Model

Take the `SBMLModel_t` pointer and extract all information required to make a
valid [`SBML.Model`](@ref) structure.
"""
function extract_model(mdl::VPtr)::SBML.Model
    # get the FBC plugin pointer (FbcModelPlugin_t)
    mdl_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), mdl, "fbc")

    # get the parameters
    parameters = Dict{String,Tuple{Float64,String}}()
    for i = 1:ccall(sbml(:Model_getNumParameters), Cuint, (VPtr,), mdl)
        p = ccall(sbml(:Model_getParameter), VPtr, (VPtr, Cuint), mdl, i - 1)
        id, v = extract_parameter(p)
        parameters[id] = v
    end

    # parse out the unit definitions
    units = Dict{String,Vector{SBML.UnitPart}}()
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
        id = get_string(ud, :UnitDefinition_getId)
        units[id] = [
            begin
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
                SBML.UnitPart(
                    unsafe_string(
                        ccall(
                            sbml(:UnitKind_toString),
                            Cstring,
                            (Cint,),
                            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
                        ),
                    ),
                    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getScale), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u),
                )
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
        ]
    end

    # parse out compartment names
    compartments = Dict{String,Compartment}()
    for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
        co = ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1)

        compartments[get_string(co, :Compartment_getId)] = Compartment(
            get_optional_string(co, :Compartment_getName),
            get_optional_bool(co, :Compartment_isSetConstant, :Compartment_getConstant),
            get_optional_int(
                co,
                :Compartment_isSetSpatialDimensions,
                :Compartment_getSpatialDimensions,
            ),
            get_optional_double(co, :Compartment_isSetSize, :Compartment_getSize),
            get_optional_string(co, :Compartment_getUnits),
            get_notes(co),
            get_annotation(co),
        )
    end

    # parse out species
    species = Dict{String,Species}()
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t

        formula = nothing
        charge = nothing
        if sp_fbc != C_NULL
            # if the FBC plugin is present, try to get the chemical formula and charge
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
                formula = get_string(sp_fbc, :FbcSpeciesPlugin_getChemicalFormula)
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
        end

        ia = nothing
        if ccall(sbml(:Species_isSetInitialAmount), Cint, (VPtr,), sp) != 0
            ia = (
                ccall(sbml(:Species_getInitialAmount), Cdouble, (VPtr,), sp),
                get_optional_string(sp, :Species_getSubstanceUnits),
            )
        end

        ic = nothing
        if ccall(sbml(:Species_isSetInitialConcentration), Cint, (VPtr,), sp) != 0
            ic = (
                ccall(sbml(:Species_getInitialConcentration), Cdouble, (VPtr,), sp),
                get_optional_string(sp, :Species_getSubstanceUnits),
            )
        end

        species[get_string(sp, :Species_getId)] = Species(
            get_optional_string(sp, :Species_getName),
            get_string(sp, :Species_getCompartment),
            get_optional_bool(
                sp,
                :Species_isSetBoundaryCondition,
                :Species_getBoundaryCondition,
            ),
            formula,
            charge,
            ia,
            ic,
            get_optional_bool(
                sp,
                :Species_isSetHasOnlySubstanceUnits,
                :Species_getHasOnlySubstanceUnits,
            ),
            get_notes(sp),
            get_annotation(sp),
        )
    end

    # parse out the flux objectives (these are complementary to the objectives
    # that appear in the reactions, see comments lower)
    objective = Dict{String,Float64}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumObjectives), Cuint, (VPtr,), mdl_fbc)
            o = ccall(
                sbml(:FbcModelPlugin_getObjective),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )
            for j = 1:ccall(sbml(:Objective_getNumFluxObjectives), Cuint, (VPtr,), o)
                fo = ccall(sbml(:Objective_getFluxObjective), VPtr, (VPtr, Cuint), o, j - 1)
                objective[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
            end
        end
    end

    # reactions!
    reactions = Dict{String,Reaction}()
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
        kinetic_parameters = Dict{String,Tuple{Float64,String}}()
        lower_bound = nothing
        upper_bound = nothing
        math = nothing

        # kinetic laws store a second version of the bounds and objectives
        kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
        if kl != C_NULL
            for j = 1:ccall(sbml(:KineticLaw_getNumParameters), Cuint, (VPtr,), kl)
                p = ccall(sbml(:KineticLaw_getParameter), VPtr, (VPtr, Cuint), kl, j - 1)
                id, v = extract_parameter(p)
                parameters[id] = v
                kinetic_parameters[id] = v
            end

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
            end
        end

        re_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), re, "fbc")
        if re_fbc != C_NULL
            lower_bound = get_optional_string(re_fbc, :FbcReactionPlugin_getLowerFluxBound)
            upper_bound = get_optional_string(re_fbc, :FbcReactionPlugin_getUpperFluxBound)
        end

        # extract stoichiometry
        reactants = Dict{String,Float64}()
        products = Dict{String,Float64}()

        add_stoi(sr, coll) = begin
            s = get_string(sr, :SpeciesReference_getSpecies)
            coll[s] =
                ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr)
        end

        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
            add_stoi(sr, reactants)
        end

        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
            add_stoi(sr, products)
        end

        # gene product associations
        association = nothing
        if re_fbc != C_NULL
            gpa = ccall(
                sbml(:FbcReactionPlugin_getGeneProductAssociation),
                VPtr,
                (VPtr,),
                re_fbc,
            )
            if gpa != C_NULL
                a = ccall(sbml(:GeneProductAssociation_getAssociation), VPtr, (VPtr,), gpa)
                a != C_NULL
                association = get_association(a)
            end
        end

        # explicit reversible flag (defaults to true in SBML)
        reversible = Bool(ccall(sbml(:Reaction_getReversible), Cint, (VPtr,), re))

        reid = get_string(re, :Reaction_getId)
        reactions[reid] = Reaction(;
            reactants,
            products,
            kinetic_parameters,
            lower_bound,
            upper_bound,
            gene_product_association = association,
            kinetic_math = math,
            reversible,
            notes = get_notes(re),
            annotation = get_annotation(re),
        )
    end

    # extract gene products
    gene_products = Dict{String,GeneProduct}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumGeneProducts), Cuint, (VPtr,), mdl_fbc)
            gp = ccall(
                sbml(:FbcModelPlugin_getGeneProduct),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )

            id = get_optional_string(gp, :GeneProduct_getId) # IDs don't need to be set

            if id != nothing
                gene_products[id] = GeneProduct(
                    get_optional_string(gp, :GeneProduct_getName),
                    get_optional_string(gp, :GeneProduct_getLabel),
                    get_notes(gp),
                    get_annotation(gp),
                )
            end
        end
    end

    function_definitions = Dict{String,FunctionDefinition}()
    for i = 1:ccall(sbml(:Model_getNumFunctionDefinitions), Cuint, (VPtr,), mdl)
        fd = ccall(sbml(:Model_getFunctionDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
        def = nothing
        if ccall(sbml(:FunctionDefinition_isSetMath), Cint, (VPtr,), fd) != 0
            def = parse_math(ccall(sbml(:FunctionDefinition_getMath), VPtr, (VPtr,), fd))
        end

        function_definitions[get_string(fd, :FunctionDefinition_getId)] =
            FunctionDefinition(
                get_optional_string(fd, :FunctionDefinition_getName),
                def,
                get_notes(fd),
                get_annotation(fd),
            )
    end

    initial_assignments = Dict{String,Math}()
    num_ias = ccall(sbml(:Model_getNumInitialAssignments), Cuint, (VPtr,), mdl)
    for n = 0:(num_ias-1)
        ia = ccall(sbml(:Model_getInitialAssignment), VPtr, (VPtr, Cuint), mdl, n)
        sym = ccall(sbml(:InitialAssignment_getSymbol), Cstring, (VPtr,), ia)
        math_ptr = ccall(sbml(:InitialAssignment_getMath), VPtr, (VPtr,), ia)
        if math_ptr != C_NULL
            initial_assignments[unsafe_string(sym)] = parse_math(math_ptr)
        end
    end

    # events
    events = Dict{String,Event}()
    num_events = ccall(sbml(:Model_getNumEvents), Cuint, (VPtr,), mdl)
    for n = 0:(num_events-1)
        ev = ccall(sbml(:Model_getEvent), VPtr, (VPtr, Cuint), mdl, n)

        event_assignments = EventAssignment[]
        for j = 0:(ccall(sbml(:Event_getNumEventAssignments), Cuint, (VPtr,), ev)-1)
            eva = ccall(sbml(:Event_getEventAssignment), VPtr, (VPtr, Cuint), ev, j)
            eva_math_ptr = ccall(sbml(:EventAssignment_getMath), VPtr, (VPtr,), eva)
            push!(
                event_assignments,
                EventAssignment(
                    unsafe_string(
                        ccall(sbml(:EventAssignment_getVariable), Cstring, (VPtr,), eva),
                    ),
                    eva_math_ptr == C_NULL ? nothing : parse_math(eva_math_ptr),
                ),
            )
        end

        trig_math_ptr = ccall(
            sbml(:Trigger_getMath),
            VPtr,
            (VPtr,),
            ccall(sbml(:Event_getTrigger), VPtr, (VPtr,), ev),
        )

        events[unsafe_string(ccall(sbml(:Event_getId), Cstring, (VPtr,), ev))] = SBML.Event(
            get_optional_string(ev, :Event_getName),
            trig_math_ptr == C_NULL ? nothing : parse_math(trig_math_ptr),
            event_assignments,
        )
    end

    # Rules
    rules = Rule[]
    num_rules = ccall(sbml(:Model_getNumRules), Cuint, (VPtr,), mdl)
    for n = 0:(num_rules-1)
        rule_ptr = ccall(sbml(:Model_getRule), VPtr, (VPtr, Cuint), mdl, n)
        type = if ccall(sbml(:Rule_isAlgebraic), Bool, (VPtr,), rule_ptr)
            AlgebraicRule
        elseif ccall(sbml(:Rule_isAssignment), Bool, (VPtr,), rule_ptr)
            AssignmentRule
        elseif ccall(sbml(:Rule_isRate), Bool, (VPtr,), rule_ptr)
            RateRule
        end
        if type in (AssignmentRule, RateRule)
            var = ccall(sbml(:Rule_getVariable), Cstring, (VPtr,), rule_ptr)
        end
        math_ptr = ccall(sbml(:Rule_getMath), VPtr, (VPtr,), rule_ptr)
        if math_ptr != C_NULL
            math = parse_math(math_ptr)
            rule = if type in (AssignmentRule, RateRule)
                type(unsafe_string(var), math)
            else
                type(math)
            end
            push!(rules, rule)
        end
    end

    # Constraints
    constraints = Constraint[]
    num_constraints = ccall(sbml(:Model_getNumConstraints), Cuint, (VPtr,), mdl)
    for n = 0:(num_constraints-1)
        constraint_ptr = ccall(sbml(:Model_getConstraint), VPtr, (VPtr, Cuint), mdl, n)
        xml_ptr = ccall(sbml(:Constraint_getMessage), VPtr, (VPtr,), constraint_ptr)
        message = get_string_from_xmlnode(xml_ptr)
        math_ptr = ccall(sbml(:Constraint_getMath), VPtr, (VPtr,), constraint_ptr)
        if math_ptr != C_NULL
            math = parse_math(math_ptr)
            constraint = Constraint(math, message)
            push!(constraints, constraint)
        end
    end

    return Model(;
        parameters,
        units,
        compartments,
        species,
        initial_assignments,
        rules,
        constraints,
        reactions,
        objective,
        gene_products,
        function_definitions,
        events,
        area_units = get_optional_string(mdl, :Model_getAreaUnits),
        extent_units = get_optional_string(mdl, :Model_getExtentUnits),
        length_units = get_optional_string(mdl, :Model_getLengthUnits),
        substance_units = get_optional_string(mdl, :Model_getSubstanceUnits),
        time_units = get_optional_string(mdl, :Model_getTimeUnits),
        volume_units = get_optional_string(mdl, :Model_getVolumeUnits),
        notes = get_notes(mdl),
        annotation = get_annotation(mdl),
    )
end