text:
using OrdinaryDiffEq, DiffEqDevTools
probArr = Vector{ODETestProblem}(2)
bigprobArr = Vector{ODETestProblem}(2)

const linear_bigα = parse(BigFloat,"1.01")
f = (t,u) -> (linear_bigα*u)
analytic = (t,u0) -> u0*exp(linear_bigα*t)
"""Linear ODE on Float64"""
prob_ode_bigfloatlinear = ODETestProblem(f,parse(BigFloat,"0.5"),analytic)

f = (t,u,du) -> begin
  for i in 1:length(u)
    du[i] = linear_bigα*u[i]
  end
end
"""2D Linear ODE, bigfloats"""
prob_ode_bigfloat2Dlinear = ODETestProblem(f,map(BigFloat,rand(4,2)).*ones(4,2)/2,analytic)

probArr[1] = prob_ode_linear
probArr[2] = prob_ode_2Dlinear
bigprobArr[1] = prob_ode_bigfloatlinear
bigprobArr[2] = prob_ode_bigfloat2Dlinear
setprecision(400)
srand(100)
## Convergence Testing
println("Convergence Test on Linear")
dts = 1.//2.^(8:-1:4)
testTol = 0.3
superduperbool = Vector{Bool}(2)
alg = ExplicitRK()

for i = 1:2 # 1 = num, 2 = ExplicitRK
  if i>1
    prob = probArr[2]
    bigprob = bigprobArr[2]
  else
    prob = probArr[1]
    bigprob = bigprobArr[1]
  end

  # Order 2

  tab = constructHeun()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  tab = constructRalston()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  # Order 3

  tab = constructBogakiShampine3()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-3) < testTol

  # Order 4

  tab = constructRKF4()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-4) < testTol

  # Order 5

  dts = 1.//2.^(7:-1:4)
  tab = constructRKF5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1.//2.^(7:-1:4)
  tab = constructDormandPrince()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  sim2 = test_convergence(dts,prob,DP5())
  @test (abs(sim.𝒪est[:l∞]-5) < testTol && (maximum(sim[end][end]-sim2[end][end]) < 1e-10))

  tab = constructCashKarp()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1.//2.^(7:-1:4)
  tab = constructRungeFirst5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructCassity5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructLawson5()
  sim = test_convergence(dts,prob,alg,tableau=tab) #10
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructLutherKonen5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructLutherKonen52()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructLutherKonen53()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructPapakostasPapaGeorgiou5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructPapakostasPapaGeorgiou52()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1.//2.^(6:-1:4)
  tab = constructTsitouras5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1.//2.^(6:-1:4)
  tab = constructBogakiShampine5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tab = constructSharpSmart5()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  # Order 6

  dts = 1.//2.^(6:-1:4)
  tab = constructButcher6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructButcher62()
  sim = test_convergence(dts,prob,alg,tableau=tab) #20
  @test abs(sim.𝒪est[:l∞]-6) < testTol # Less stringent

  dts = 1.//2.^(6:-1:4)
  tab = constructButcher63()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1.//2.^(5:-1:1)
  tab = constructDormandPrince6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol # Better on linear

  tab = constructSharpVerner6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructVerner916()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructVerner9162()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructVernerRobust6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructVernerEfficient6(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6.6) < testTol

  tab = constructPapakostas6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructLawson6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1.//2.^(3:-1:1)
  tab = constructTsitourasPapakostas6()
  sim = test_convergence(dts,prob,alg,tableau=tab) #30
  @test abs(sim.𝒪est[:l∞]-6.7) < testTol # Better on linear

  dts = 1.//2.^(5:-1:1)
  tab = constructDormandLockyerMcCorriganPrince6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructTanakaKasugaYamashitaYazaki6D()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructTanakaKasugaYamashitaYazaki6C()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructTanakaKasugaYamashitaYazaki6B()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructTanakaKasugaYamashitaYazaki6A()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructMikkawyEisa()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6.53) < testTol # Odd behavior

  tab = constructChummund6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructChummund62()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructHuta6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-5.5) < testTol # Low convergence, error noted in Stone notes

  dts = 1.//2.^(5:-1:1)
  tab = constructHuta62()
  sim = test_convergence(dts,prob,alg,tableau=tab)#40
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructVerner6()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6.7) < testTol # Better on linear

  dts = 1.//2.^(4:-1:1)
  tab = constructDverk()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tab = constructClassicVerner6()
  sim = test_convergence(dts,prob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  # Order 7

  dts = 1.//2.^(5:-1:1)
  tab = constructButcher7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1.//2.^(5:-1:2)
  tab = constructClassicVerner7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  tab = constructVernerRobust7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1.//2.^(5:-1:1)
  tab = constructEnrightVerner7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7.15) < testTol # Better on linear

  tab = constructTanakaYamashitaStable7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7.3) < testTol

  tab = constructTanakaYamashitaEfficient7(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1.//2.^(8:-1:3)
  tab = constructSharpSmart7(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab) #50
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1.//2.^(3:-1:1)
  tab = constructSharpVerner7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-6.5) < testTol # Coefficients aren't accurate enough, drop off error

  tab = constructVernerEfficient7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  tab = constructVerner7()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  # Order 8
  dts = 1.//2.^(4:-1:1)
  tab = constructClassicVerner8()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(5:-1:1)
  tab = constructVerner8(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructCooperVerner8()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol #Coefficients not accurate enough

  dts = 1.//2.^(4:-1:1)
  tab = constructCooperVerner82()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol #Coefficients not accurate enough

  tab = constructTsitourasPapakostas8(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructdverk78(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructEnrightVerner8(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructCurtis8()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1.//2.^(4:-1:1)
  tab = constructRKF8(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  tab = constructDormandPrince8(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8.4) < testTol

  dts = 1.//2.^(3:-1:1)
  tab = constructDormandPrince8_64bit(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-8.4) < testTol

  # Order 9

  tab = constructVernerRobust9(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-9) < testTol

  tab = constructVernerEfficient9(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-9) < testTol

  dts = 1.//2.^(3:-1:1)
  tab = constructSharp9(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-9) < testTol #Only works to Float64 precision

  dts = 1.//2.^(2:-1:1)
  tab = constructTsitouras9(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10.5) < testTol #Only works to Float64

  dts = 1.//2.^(3:-1:1)
  tab = constructTsitouras92(BigFloat)
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-9) < testTol  #Only works to Float64

  ## Order 10

  dts = 1.//2.^(5:-1:1)
  tab = constructCurtis10()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  tab = constructOno10()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  dts = 1.//2.^(5:-1:1)
  tab = constructFeagin10Tableau()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  tab = constructCurtis10()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  dts = 1.//2.^(6:-1:1)
  tab = constructBaker10()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10.8) < testTol

  tab = constructHairer10()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-10.8) < testTol

  ## Order 12

  dts = 1.//2.^(6:-1:1)
  tab = constructFeagin12Tableau()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-12.6) < testTol

  tab = constructOno12()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-11.6) < testTol

  ## Order 14

  dts = 1.//2.^(6:-1:1)
  tab = constructFeagin14Tableau()
  sim = test_convergence(dts,bigprob,alg,tableau=tab)
  @test abs(sim.𝒪est[:l∞]-15.5) < testTol
end

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
function modalparam_test()
	# State space model for testing modalparam
	# and modalparami
	#
	# javier.cara@upm.es

	M = [35.0 0.0;0.0 17.5] # mass matrix
	K = [12250.0 -3500.0;-3500.0 3500.0] # stiffness matrix
	g = [0.02;0.02] # damping ratios

	# teorethical eigenvalues and eigenvectors
	W2,V = eig(K,M)
	W = sqrt(W2)

	# damping matrix
	Mm = V'*M*V
	Gm = 2*Mm*diagm(W)*diagm(g)
	G = inv(V')*Gm*inv(V) # damping matrix

	# state-space matrices
	Minv = [ 1/M[1,1] 0.0;0.0 1/M[2,2] ]
	Ac = [zeros(2,2) eye(2,2);-Minv*K -Minv*G] # continuous A matrix
	dt = 0.02 # time step
	A = expm(Ac*dt) # discrete A matrix
	C = [-Minv*K -Minv*G] # discrete C matrix
	Bc = [zeros(2,2);Minv]
	B = (A - eye(4))*inv(Ac)*Bc

	# eigenvectors with max. component = 1
	V1 = zeros(2,2)
	for j in 1:2
		maxval = V[1,j]
		if abs(V[2,j]) > abs(maxval)
			maxval = V[2,j]
		end
		V1[:,j] = V[:,j]/maxval
	end
	# updated mass masses
	Mm1 = V1'*M*V1
	mm1 = diag(Mm1)

	ssm = Dict("M"=>M,"K"=>K,"G"=>G,"wm"=>W,"zm"=>g,"Vm"=>V1,"mm"=>mm1,
							"dt"=>dt,"A"=>A,"B"=>B,"C"=>C)

	return ssm

end

ext: jl


text:
import DataFrames
import PredictMD
import Random
import Test

Test.@testset "linearly independent" begin
    num_rows = 1_000
    x = randn(Random.MersenneTwister(1), num_rows)
    x1 = randn(Random.MersenneTwister(2), num_rows)
    x11 = randn(Random.MersenneTwister(3), num_rows)
    y = randn(Random.MersenneTwister(4), num_rows)
    y1 = randn(Random.MersenneTwister(5), num_rows)
    y11 = randn(Random.MersenneTwister(6), num_rows)
    z = randn(Random.MersenneTwister(7), num_rows)
    z1 = randn(Random.MersenneTwister(8), num_rows)
    z11 = randn(Random.MersenneTwister(9), num_rows)
    df = DataFrames.DataFrame(
        :x => x,
        :x1 => x1,
        :x11 => x11,
        :y => y,
        :y1 => y1,
        :y11 => y11,
        :z => z,
        :z1 => z1,
        :z11 => z11,
        )
    Test.@test(PredictMD.columns_are_linearly_independent(df))
    Test.@test(
        PredictMD.columns_are_linearly_independent(
            df,
            [:x, :x1, :x11, :y, :y1, :y11, :z, :z1, :z11],
            )
        )
    Test.@test(length(PredictMD.linearly_dependent_columns(df)) == 0)
    Test.@test(
        length(
            PredictMD.linearly_dependent_columns(
                df,
                [:x, :x1, :x11, :y, :y1, :y11, :z, :z1, :z11],
                ),
            ) == 0
        )
end

Test.@testset "linearly dependent" begin
    num_rows = 1_000
    x = randn(Random.MersenneTwister(10), num_rows)
    x1 = randn(Random.MersenneTwister(11), num_rows)
    x11 = randn(Random.MersenneTwister(12), num_rows)
    y = randn(Random.MersenneTwister(13), num_rows)
    y1 = randn(Random.MersenneTwister(14), num_rows)
    y11 = randn(Random.MersenneTwister(15), num_rows)
    z = 2*x .+ 3*y
    z1 = randn(Random.MersenneTwister(17), num_rows)
    z11 = randn(Random.MersenneTwister(18), num_rows)
    df = DataFrames.DataFrame(
        :x => x,
        :x1 => x1,
        :x11 => x11,
        :y => y,
        :y1 => y1,
        :y11 => y11,
        :z => z,
        :z1 => z1,
        :z11 => z11,
        )
    Test.@test(!PredictMD.columns_are_linearly_independent(df))
    Test.@test(
        !PredictMD.columns_are_linearly_independent(
            df,
            [:x, :x1, :x11, :y, :y1, :y11, :z, :z1, :z11],
            )
        )
    Test.@test(length(PredictMD.linearly_dependent_columns(df)) == 1)
    Test.@test(
        length(
            PredictMD.linearly_dependent_columns(
                df,
                [:x, :x1, :x11, :y, :y1, :y11, :z, :z1, :z11],
                )
            ) == 1
        )
    result1 = PredictMD.linearly_dependent_columns(df)
    result2 = PredictMD.linearly_dependent_columns(
        df,
        [:x, :x1, :x11, :y, :y1, :y11, :z, :z1, :z11],
        )
    Test.@test(
        all(result1.==[:x]) || all(result1.==[:y]) || all(result1.==[:z])
        )
    Test.@test(
        all(result2.==[:x]) || all(result2.==[:y]) || all(result2.==[:z])
        )

    DataFrames.deletecols!(df, [:z])
    Test.@test(PredictMD.columns_are_linearly_independent(df))
    Test.@test(
        PredictMD.columns_are_linearly_independent(
            df,
            [:x, :x1, :x11, :y, :y1, :y11, :z1, :z11],
            )
        )
    Test.@test(length(PredictMD.linearly_dependent_columns(df)) == 0)
    Test.@test(
        length(
            PredictMD.linearly_dependent_columns(
                df,
                [:x, :x1, :x11, :y, :y1, :y11, :z1, :z11],
                )
            ) == 0
        )
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.linegraph" begin
    p5L = linegraph(pathgraph(5))
    @test nodecount(p5L) == 4
    @test edgecount(p5L) == 3
    c5L = linegraph(cyclegraph(5))
    @test nodecount(c5L) == 5
    @test edgecount(c5L) == 5
    lad5L = linegraph(laddergraph(5))
    @test nodecount(lad5L) == 13
    @test edgecount(lad5L) == 22
end # graph.linegraph


text:
"""
    _inceptionblock(inplanes, out_1x1, red_3x3, out_3x3, red_5x5, out_3x3, pool_proj)

Create an inception module for use in GoogLeNet
([reference](https://arxiv.org/abs/1409.4842v1)).

# Arguments
- `inplanes`: the number of input feature maps
- `out_1x1`: the number of output feature maps for the 1x1 convolution (branch 1)
- `red_3x3`: the number of output feature maps for the 3x3 reduction convolution (branch 2)
- `out_3x3`: the number of output feature maps for the 3x3 convolution (branch 2)
- `red_5x5`: the number of output feature maps for the 5x5 reduction convolution (branch 3)
- `out_5x5`: the number of output feature maps for the 5x5 convolution (branch 3)
- `pool_proj`: the number of output feature maps for the pooling projection (branch 4)
"""
function _inceptionblock(inplanes, out_1x1, red_3x3, out_3x3, red_5x5, out_5x5, pool_proj)
    branch1 = Chain(Conv((1, 1), inplanes => out_1x1))

    branch2 = Chain(Conv((1, 1), inplanes => red_3x3),
                    Conv((3, 3), red_3x3 => out_3x3; pad = 1))

    branch3 = Chain(Conv((1, 1), inplanes => red_5x5),
                    Conv((5, 5), red_5x5 => out_5x5; pad = 2))

    branch4 = Chain(MaxPool((3, 3), stride=1, pad = 1),
                    Conv((1, 1), inplanes => pool_proj))

    return Parallel(cat_channels,
                    branch1, branch2, branch3, branch4)
end

"""
    googlenet(; nclasses = 1000)

Create an Inception-v1 model (commonly referred to as GoogLeNet)
([reference](https://arxiv.org/abs/1409.4842v1)).

# Arguments
- `nclasses`: the number of output classes
"""
function googlenet(; nclasses = 1000)
  layers = Chain(Chain(Conv((7, 7), 3 => 64; stride = 2, pad = 3),
                       MaxPool((3, 3), stride = 2, pad = 1),
                       Conv((1, 1), 64 => 64),
                       Conv((3, 3), 64 => 192; pad = 1),
                       MaxPool((3, 3), stride = 2, pad = 1),
                       _inceptionblock(192, 64, 96, 128, 16, 32, 32),
                       _inceptionblock(256, 128, 128, 192, 32, 96, 64),
                       MaxPool((3, 3), stride = 2, pad = 1),
                       _inceptionblock(480, 192, 96, 208, 16, 48, 64),
                       _inceptionblock(512, 160, 112, 224, 24, 64, 64),
                       _inceptionblock(512, 128, 128, 256, 24, 64, 64),
                       _inceptionblock(512, 112, 144, 288, 32, 64, 64),
                       _inceptionblock(528, 256, 160, 320, 32, 128, 128),
                       MaxPool((3, 3), stride = 2, pad = 1),
                       _inceptionblock(832, 256, 160, 320, 32, 128, 128),
                       _inceptionblock(832, 384, 192, 384, 48, 128, 128),
                       AdaptiveMeanPool((1, 1))),
                 Chain(flatten,
                       Dropout(0.4),
                       Dense(1024, nclasses)))

  return layers
end

"""
    GoogLeNet(; pretrain = false,  nclasses = 1000)

Create an Inception-v1 model (commonly referred to as `GoogLeNet`)
([reference](https://arxiv.org/abs/1409.4842v1)).

# Arguments
- `pretrain`: set to `true` to load the model with pre-trained weights for ImageNet
- `nclasses`: the number of output classes

See also [`googlenet`](#).
"""
struct GoogLeNet{T}
  layers::T
end

function GoogLeNet(; pretrain = false, nclasses = 1000)
  layers = googlenet(nclasses = nclasses)
  pretrain && Flux.loadparams!(layers, weights("googlenet"))

  GoogLeNet(layers)
end

@functor GoogLeNet

(m::GoogLeNet)(x) = m.layers(x)


text:
using Documenter
using Phylo

makedocs(modules = [Phylo, Phylo.API],
         clean   = false)

deploydocs(deps = Deps.pip("pygments",
                           "mkdocs",
                           "mkdocs-material",
                           "python-markdown-math"),
           repo = "github.com/richardreeve/Phylo.jl.git",
           julia="0.6",
           osname="linux")


text:
using LinearAlgebra, SparseArrays, Random, Statistics
using Plots, Plots.PlotMeasures

function gomp(y0::AbstractVector, Ψ::AbstractMatrix,
    K::Int = 2, S::Int = 1;
    max_iter::Int = 100, ϵ::Real = eps())

	# Get the dimensions
	m = length(y0)
	m_psi, n = size(Ψ)
	# Assert the dimensionality
	@assert m == m_psi "Please provide consistent input sizes"
	# Assert the selector
	@assert S <=  min(K, m/K) "S <= min(K, m/K)"

	# Normalize
	ψ = deepcopy(Ψ)
    # Inplace
	normalize!.(eachcol(ψ), 2)

	# Preconditioning
	P = ψ'pinv(ψ*ψ')
	# New matrix
	ψ = P*ψ
	y = P*y0
	# Iteration
	iters = 0

	# Support
	Λ = zeros(Bool, n)
	u = zeros(eltype(y), n)
	r = y
	amps = zeros(eltype(y), n)

	# Find the magnitudes
	for i in 1:max_iter
		# Compute the similarities via magnitude
		amps .= ψ'r
		# Get the largest entry
		idx = sortperm(abs.(amps), rev = true)[1:S]
		# Update the support
		Λ[idx] .= true
		# Update the coefficients
		u[Λ] .= (y' / ψ')[1, Λ]
		# Update r
		r .= y - ψ*u
		# Convergence
		(norm(r,2) < ϵ || sum(Λ) >= K) && break
	end

    # Last time to get the right coefficients
	u[Λ] .= (y0' / Ψ')[1, Λ]

	return u
end

## Test data
Random.seed!(1111) # Lucky
x̂ = sprandn(Float64, 100, 0.1)
k_opt = Int(norm(x̂, 0))
A = randn(Float64, 300, 100)
y = A*x̂

# Evaluate
x = gomp(y, A, k_opt)

# Check for different sparsties levels
sp_n = map(0.01:0.01:0.2) do i
	x̃ = sprandn(Float64, 100, i)
	ỹ = A*x̃
	norm(gomp(ỹ, A, 50, 1, ϵ = 0.01), 0), norm(x̃, 0)
end

gr()
pl_1 = scatter(x̂, label = "Ground Truth", markeralpha = 0.5, xlabel = "Index", ylabel = "x", legend = :bottomright)
scatter!(x, label = "Estimated", markershape = :cross, markersize = 12, markerstrokewidth = 5)
savefig(joinpath(pwd(), "figures", "gomp.pdf"))

xs, ys = first.(sp_n), last.(sp_n)
pl_2 = scatter(ys, xlabel = "Run", ylabel = "Sparsity", label = "Ground Truth", legend = :topleft)
scatter!(xs, label = "Estimated", markershape = :cross, markersize = 12, markerstrokewidth = 5)
savefig(joinpath(pwd(), "figures", "gomp_vary.pdf"))

plot(pl_1, pl_2, layout = (1,2), size = (900, 300), margins = 5mm)
savefig(joinpath(pwd(), "figures", "merged.pdf"))


text:
using Test, FinancialMonteCarlo;
@show "MertonProcess"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
lam = 5.0;
mu1 = 0.03;
sigma1 = 0.02;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianData = AsianFloatingStrikeOption(T)
Model = MertonProcess(sigma, lam, mu1, sigma1, Underlying(S0, d));

@show FwdPrice = pricer(Model, rfCurve, mc, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc, EUData);
@show AmPrice = pricer(Model, rfCurve, mc, AMData);
@show BarrierPrice = pricer(Model, rfCurve, mc, BarrierData);
@show AsianPrice = pricer(Model, rfCurve, mc, AsianData);

@test abs(FwdPrice - 99.1188767166039) < toll
@test abs(EuPrice - 9.084327245917533) < toll
@test abs(BarrierPrice - 7.880881290426765) < toll
@test abs(AsianPrice - 5.129020349580892) < toll

@show FwdPrice = pricer(Model, rfCurve, mc1, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc1, EUData);
@show AmPrice = pricer(Model, rfCurve, mc1, AMData);
@show BarrierPrice = pricer(Model, rfCurve, mc1, BarrierData);
@show AsianPrice = pricer(Model, rfCurve, mc1, AsianData);
tollanti = 0.6;
@test abs(FwdPrice - 99.1188767166039) < tollanti
@test abs(EuPrice - 9.084327245917533) < tollanti
@test abs(BarrierPrice - 7.880881290426765) < tollanti
@test abs(AsianPrice - 5.129020349580892) < tollanti

@show "Test Merton Parameters"

@test_throws(ErrorException, MertonProcess(-sigma, lam, mu1, sigma1, Underlying(S0, d)))
@test_throws(ErrorException, MertonProcess(sigma, lam, mu1, -sigma1, Underlying(S0, d)))
@test_throws(ErrorException, MertonProcess(sigma, -lam, mu1, sigma1, Underlying(S0, d)))


text:
# Plotting macros
macro makeLivePlot(plotExpression)
    # Turn plotting function into smooth version
    plotFunction = plotExpression.args[1]
    plotArgs = plotExpression.args[2:end]

    arrayArgs = map(x -> :($(esc(x))), plotArgs)
    splatArgs = :([$(arrayArgs...)])

    outExpression = :(smoothLivePlotGeneral($(esc(plotFunction)), $splatArgs))
    return outExpression
end


text:
module Persistence

using SQLite
using Niancat
using Niancat.Users
using Niancat.Gameface

struct GameInstanceDescription
    instanceinfo::InstanceInfo
    gamename::String
    gameinstanceid::Int
    state::String
end

function GameInstanceDescription(row)
    instanceinfo = InstanceInfo(row[:instance_id], row[:instance_name])
    GameInstanceDescription(instanceinfo, row[:game_name], row[:game_instance_id], row[:game_state])
end

function initializedatabase!(db::SQLite.DB)
    sqlscript = read(joinpath(pkgdir(Niancat), "extra", "niancat.sql"), String)
    stms = split(sqlscript, "\n\n")
    for stm in stms
        DBInterface.execute(db, stm)
    end
end

struct GamePersistence <: GameEventPersistence
    db::SQLite.DB

    function GamePersistence()
        db = SQLite.DB()
        initializedatabase!(db)
        new(db)
    end

    GamePersistence(db::SQLite.DB) = new(db)
end

function getgameinstances(persistence::GamePersistence) :: Vector{GameInstanceDescription}
    sql = """
        SELECT
            instances.instance_id, instances.instance_name,
            games.game_name,
            gameinstances.game_state,
            gameinstances.game_instance_id
        FROM gameinstances
            JOIN games ON gameinstances.game_id = games.game_id
            JOIN instances ON gameinstances.instance_id = instances.instance_id;
    """
    results = DBInterface.execute(persistence.db, sql)

    [GameInstanceDescription(row) for row in results]
end

function newgameinstance!(persistence::GamePersistence, instancename::String, gamename::String, initialstate::String)
    sql = """
        INSERT INTO gameinstances
            (game_id, instance_id, game_state)
        VALUES (
            (SELECT game_id FROM games WHERE game_name = ?),
            (SELECT instance_id FROM instances WHERE instance_name = ?),
            ?
        )
    """
    DBInterface.execute(persistence.db, sql, (gamename, instancename, initialstate))
    game_instance_id = SQLite.last_insert_rowid(persistence.db)

    readback_sql = """
        SELECT
            instances.instance_id, instances.instance_name,
            games.game_name,
            gameinstances.game_state,
            gameinstances.game_instance_id
        FROM gameinstances
            JOIN games ON gameinstances.game_id = games.game_id
            JOIN instances ON gameinstances.instance_id = instances.instance_id
        WHERE gameinstances.game_instance_id = ?
    """
    results = DBInterface.execute(persistence.db, readback_sql, (game_instance_id, ))
    row = first(results)

    GameInstanceDescription(row)
end

function getuser(persistence::GamePersistence, userid::String, teamname::String) :: User
    sql = """
    INSERT OR IGNORE INTO users
        (team_user_id, team_id, display_name)
    VALUES
        (?,
         (SELECT team_id FROM teams WHERE team_name = ?),
         ?);
    """

    DBInterface.execute(persistence.db, sql, (userid, teamname, userid))

    getsql = """
    SELECT users.user_id, users.team_id, users.display_name, teams.icon
    FROM users JOIN teams ON users.team_id = teams.team_id
    WHERE
        users.team_user_id = ? AND
        teams.team_name = ?;
    """

    results = DBInterface.execute(persistence.db, getsql, (userid, teamname))
    row = first(results)

    userdatabaseid = row[1]
    teamdatabaseid = row[2]
    displayname = row[3]
    teamicon = row[4]
    team = Team(teamdatabaseid, teamname, teamicon)

    User(userdatabaseid, userid, displayname, team)
end

function updatedisplayname!(persistence::GamePersistence, user::User, displayname::String)
    sql = """
        UPDATE users
        SET display_name = ?
        WHERE user_id = ?
    """
    DBInterface.execute(persistence.db, sql, (displayname, user.databaseid))
end

function listinstancenames(persistence::GamePersistence) :: Vector{String}
    sql = "SELECT instance_name FROM instances;"

    results = DBInterface.execute(persistence.db, sql)

    [row[:instance_name] for row in results]
end

function declareinstance!(persistence::GamePersistence, instancename::String)
    sql = """
        INSERT OR IGNORE INTO instances (instance_name)
        VALUES (?);
    """
    DBInterface.execute(persistence.db, sql, (instancename,))
end

function updatenotificationendpoint!(persistence::GamePersistence, teamname::AbstractString, uri::AbstractString)
    sql = """
        REPLACE INTO teamnotifications (team_id, uri)
        VALUES
            (
                (SELECT team_id FROM teams WHERE team_name = ?),
                ?
            );
    """
    DBInterface.execute(persistence.db, sql, (teamname, uri))
end

function getnotificationendpoints(persistence::GamePersistence, instancedatabaseid::Int) :: Vector{String}
    sql = """
            SELECT uri
            FROM teamnotifications
            JOIN teams ON teams.team_id = teamnotifications.team_id
            WHERE teams.instance_id = ?
            ;"""

    results = DBInterface.execute(persistence.db, sql, (instancedatabaseid,))

    [row[:uri] for row in results]
end

function recordevent!(persistence::GamePersistence, gameinstanceid::Int, event::GameUserEvent)
    sql = """
        INSERT INTO userevents
            (game_instance_id, user_id, round, event_type, event_data)
        VALUES
            (?,                ?,       ?,     ?,          ?);
    """
    params = (gameinstanceid, event.user.databaseid, event.round, event.eventtype, event.data)
    DBInterface.execute(persistence.db, sql, params)
end

function recordevent!(persistence::GamePersistence, gameinstanceid::Int, event::GameEvent)
    sql = """
        INSERT INTO gameevents
            (game_instance_id, round, event_type, event_data)
        VALUES
            (?,                ?,     ?,          ?);
    """
    params = (gameinstanceid, event.round, event.eventtype, event.data)
    DBInterface.execute(persistence.db, sql, params)
end

export GamePersistence, GameInstanceDescription, getgameinstances, getuser, initializedatabase!
export updatenotificationendpoint!, getnotificationendpoints, recordevent!
export updatedisplayname!, newgameinstance!

end


text:
function Base.Matrix(F::BunchKaufman)
	U = F.U[:, F.ipiv]
	*(U, F.D, U')
end
Base.AbstractMatrix(F::BunchKaufman) = Matrix(F)
function LinearAlgebra.rdiv!(X::AbstractMatrix, F::BunchKaufman)
	# ldiv!(F, X')' # doesn't work since X' does not count as strided vector (requirement for ldiv!)
	X .= adjoint(F \ X')
end
LinearAlgebra.logdet(B::Bidiagonal) = logabsdet(B)[1]
function LinearAlgebra.logabsdet(B::Bidiagonal)
	d = @view B[diagind(B)]
	D = Diagonal(d)
	logabsdet(D)
end


text:
using CUDNN
using CuArrays
using NNlib
using Base.Test

cu_isapprox(a, ca; kw...) = isapprox(CuArray(a), ca; kw...)

@testset "conv" begin

    for T in [Float32, Float64]
        x = rand(T, 7, 7, 3, 2)
        w = rand(T, 3, 3, 3, 2)
        cx = CuArray(x)
        cw = CuArray(w)

        @test cu_isapprox(conv2d(x, w), conv2d(cx, cw))
        @test cu_isapprox(conv2d(x, w; stride=2), conv2d(cx, cw; stride=2))
        @test cu_isapprox(conv2d(x, w; padding=1), conv2d(cx, cw; padding=1))
        @test cu_isapprox(conv2d(x, w; mode=1), conv2d(cx, cw; mode=1))

        y = conv2d(x, w)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(conv2d_grad_x(x, w, dy), conv2d_grad_x(cx, cw, cdy); atol=1e-5)
        @test cu_isapprox(conv2d_grad_w(x, w, dy), conv2d_grad_w(cx, cw, cdy); atol=1e-5)

        # with stride
        y = conv2d(x, w; stride=2)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(conv2d_grad_x(x, w, dy; stride=2),
                       conv2d_grad_x(cx, cw, cdy; stride=2); atol=1e-5)
        @test cu_isapprox(conv2d_grad_w(x, w, dy; stride=2),
                       conv2d_grad_w(cx, cw, cdy; stride=2); atol=1e-5)

        # with padding
        y = conv2d(x, w; padding=1)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(conv2d_grad_x(x, w, dy; padding=1),
                       conv2d_grad_x(cx, cw, cdy; padding=1); atol=1e-5)
        @test cu_isapprox(conv2d_grad_w(x, w, dy; padding=1),
                       conv2d_grad_w(cx, cw, cdy; padding=1); atol=1e-5)
    end

end

@testset "pool" begin

    for T in [Float32, Float64]
        x = rand(T, 8, 8, 3, 2)
        cx = CuArray(x)

        @test cu_isapprox(pool(x), pool(cx))
        @test cu_isapprox(pool(x; stride=2), pool(cx; stride=2))
        @test cu_isapprox(pool(x; padding=1), pool(cx; padding=1))

        y = pool(x)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(pool_grad(x, y, dy), pool_grad(cx, cy, cdy); atol=1e-5)

        # with stride
        y = pool(x; stride=2)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(pool_grad(x, y, dy; stride=2), pool_grad(cx, cy, cdy; stride=2))

        # with padding
        y = pool(x; padding=1)
        dy = randn(T, size(y))
        cy = CuArray(y)
        cdy = CuArray(dy)
        @test cu_isapprox(pool_grad(x, y, dy; padding=1), pool_grad(cx, cy, cdy; padding=1))
        @test cu_isapprox(pool_grad(x, y, dy; padding=1), pool_grad(cx, cy, cdy; padding=1))
    end

end

@testset "batchnorm" begin

    # smoke tests
    for T in [Float32, Float64]
        x = CuArray(randn(T, 5, 4, 3, 2))
        s = BatchNormState(x)

        y = batchnorm_train(x, s)
        y = batchnorm_infer(x, s)
        dy = similar(y)
        dx = batchnorm_grad(x, dy, s)
    end

end

@testset "softmax" begin

    # softmax() in NNlib acts on 2D arrays, while cuDNN seems to only accept 4D ones
    # so it's unclear how to test it and where we need it at all
    # I keep a smoke test here to make sure functions are at least callable, but
    # their correctness is in question
    for T in [Float32, Float64]
        x = CuArray(randn(T, 5, 4, 3, 2))
        y = softmax4d(x)
        dy = similar(y)

        dx = softmax4d_grad(y, dy)
    end

end


text:
# ∫_{SO3} D^{ll}_{μμmm} D_{L1}^*(Q) e^{ab} D_{L2}(Q) dQ
# -> (2L_1+1)×(2L_2+1) matrix
# Has been checked to consist with `vec_cou_coe` when L2 = 0; b = 1;
function mat_cou_coe(rotc::Rot3DCoeffs{T},
					   ll::StaticVector{N},
	                   mm::StaticVector{N},
					   μμ::StaticVector{N},
					   L1::Integer, L2::Integer,
					   a::Integer, b::Integer) where {T,N}
	if a > 2L1 + 1 || a <= 0 || b > 2L2 +1 || b <= 0
		error("Rotation D matrices has no such element!")
	end
	Z = zeros(2 * L1 + 1, 2 * L2 + 1)
	Dp = rotation_D_matrix_ast(L1)
	Dq = rotation_D_matrix(L2)
	Dpa = Dp[:,a]
	Dqb = Dq[b,:]
	μa = [Dpa[i].μ for i in 1:2L1+1]
	ma = [Dpa[i].m for i in 1:2L1+1]
	μb = [Dqb[i].μ for i in 1:2L2+1]
	mb = [Dqb[i].m for i in 1:2L2+1]
	LL = [ll; L1; L2]
	for i = 1:(2 * L1 + 1)
		for j = 1:(2 * L2 + 1)
			MM = [mm; ma[i]; mb[j]]
			KK = [μμ; μa[i]; μb[j]]
			Z[i,j] = Dpa[i].sign * Dqb[j].sign * rotc(LL, MM, KK)
		end
	end
	return Z
end

# Has been checked to consist with vector case by choosing L2 = 0
function gramian(A::Rot3DCoeffs{T}, ll::StaticVector{N},
	                 L1::Integer, L2::Integer) where {T,N}
	LenM = 0
	Tempm = zeros(2 * L1 + 1,2 * L2 + 1)
	Dp = rotation_D_matrix_ast(L1)
	Dq = rotation_D_matrix(L2)
	Tempμ = SVector([i for i in -(L1+L2):(L1+L2)]...)
	μ_list = collect_m(ll,Tempμ)

	## TODO: the upper bound of #m_list(`(2L1 + 1)*length(μ_list)`) seems to be
	## too large and can be further reduced; just keep it for now...
	Z = fill(zeros(2L1 + 1, 2L2 + 1), (length(μ_list), (2L1 + 1)*length(μ_list)))

	for a = 1:2 * L1+1
		for b = 1:2 * L2+1
			for i = 1:2 * L1 + 1
				for j = 1: 2 * L2 +1
					Dpa = Dp[:,a]
					Dqb = Dq[b,:]
					ma = [Dpa[k].m for k in 1:2L1+1]
					mb = [Dqb[k].m for k in 1:2L2+1]
					Tempm[i,j] = ma[i]+mb[j]
				end
			end
			m_list = collect_m(ll,Tempm)
			for (im, mm) in enumerate(m_list), (iμ, μμ) in enumerate(μ_list)
				Z[iμ, im + LenM] = mat_cou_coe(A, ll, mm, μμ, L1, L2, a, b)
			end
			LenM += length(m_list)
		end
	end
	Z = Z[:, 1:LenM]

	gram = zeros(LenM, LenM)
	for i = 1:LenM
		for j = 1:LenM
			for t = 1:length(μ_list)
				gram[i,j] += tr(Z[t,i]' * Z[t,j])
			end
		end
	end

	return gram, Z, μ_list
end

# Has been checked to consist with vector case by choosing L2 = 0
function rc_basis_all(A::Rot3DCoeffs{T}, ll::StaticVector{N},
	                  L1::Integer, L2::Integer) where {N,T}
	G, C, μ_list = gramian(A, ll, L1, L2)
	S = svd(G)
	rk = rank(G; rtol =  1e-8)
	Urcpi = fill(zeros(2L1 + 1, 2L2 + 1), (rk, length(μ_list)))
	U = S.U[:, 1:rk]
	Sigma = S.S[1:rk]
	Urcpi = C * U * Diagonal(sqrt.(Sigma))^(-1)
	return Urcpi', μ_list
end

# Gramian over nn
function Gramian(A::Rot3DCoeffs,
			     nn::StaticVector{N},
			     ll::StaticVector{N},
				 L1::Integer, L2::Integer) where {N}
	Uri, μ_list = rc_basis_all(A, ll, L1, L2)
	G = zeros(size(Uri)[1], size(Uri)[1])
	for σ in permutations(1:N)
       if (nn[σ] != nn) || (ll[σ] != ll); continue; end
       for (iU1, mm1) in enumerate(μ_list), (iU2, mm2) in enumerate(μ_list)
          if mm1[σ] == mm2
             for i1 = 1:size(Uri)[1]
				 for i2 = 1:size(Uri)[1]
                 	G[i1, i2] += tr(Uri[i1, iU1] * Uri[i2, iU2]')
				end
             end
          end
       end
    end
    return G, Uri, μ_list
end

## LI matrix-valued coupling coefficients(& corresponding μ) over nn, ll
function mat_symm_basis(A::Rot3DCoeffs,
				        nn::StaticVector{N},
					    ll::StaticVector{N},
						L1::Integer, L2::Integer) where {N}
	G, C, μ_list= Gramian(A, nn, ll, L1, L2)
	S = svd(G)
	rk = rank(G; rtol =  1e-8)
	Urcpi = fill(zeros(2L1 + 1, 2L2 + 1), (rk, length(μ_list)))
	U = S.U[:, 1:rk]
	Sigma = S.S[1:rk]
	Urcpi = C' * U * Diagonal(sqrt.(Sigma))
	Z = fill(zeros(2L1 + 1, 2L2 + 1), (rk, length(μ_list)))
	for i = 1:rk
		for j = 1:length(μ_list)
			Z[i,j] = Urcpi[j,i];
		end
	end
	return Z, μ_list
end


text:
using AnalyticalSolutions
using Documenter

DocMeta.setdocmeta!(AnalyticalSolutions, :DocTestSetup, :(using AnalyticalSolutions); recursive=true)

makedocs(;
    modules=[AnalyticalSolutions],
    authors="Adrian Egger <ae.ethz@gmail.com> and contributors",
    repo="https://github.com/awegger/AnalyticalSolutions.jl/blob/{commit}{path}#{line}",
    sitename="AnalyticalSolutions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://awegger.github.io/AnalyticalSolutions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/awegger/AnalyticalSolutions.jl",
    devbranch="main",
)


text:
using ExecutableSpecifications.Gherkin

struct StepDefinitionLocation
    filename::String
    lineno::Int
end

struct NoMatchingStepDefinition <: Exception end
struct NonUniqueStepDefinition <: Exception
    locations::Vector{StepDefinitionLocation}
end

abstract type StepDefinitionMatcher end

# A StepDefinitionMatcher should define a method
# findstepdefinition(::StepDefinitionMatcher, ::Gherkin.ScenarioStep)

struct StepDefinition
    description::String
    definition::Function
    location::StepDefinitionLocation
end

#
# Step definition context
#
struct StepDefinitionContext
    variables::Dict{Symbol, Any}

    StepDefinitionContext() = new(Dict{Symbol, Any}())
end

Base.getindex(context::StepDefinitionContext, sym::Symbol) = context.variables[sym]
Base.setindex!(context::StepDefinitionContext, value::Any, sym::Symbol) = context.variables[sym] = value

#
# Global state
#

currentdefinitions = Vector{StepDefinition}()
currentfilename = ""

#
# Step definition macros
#

function step_definition_(description::String, definition::Expr)
    definitionfunction = :(context -> $definition)
    quote
        push!(currentdefinitions, StepDefinition($description, (context) -> begin
            try
                $(esc(definitionfunction))(context)
                SuccessfulStepExecution()
            catch ex
                if ex isa StepAssertFailure
                    StepFailed(ex.assertion)
                else
                    rethrow()
                end
            end
        end,
        StepDefinitionLocation(currentfilename, 0)))
    end
end

macro given(description, definition)
    step_definition_(description, definition)
end

macro when(description, definition)
    step_definition_(description, definition)
end

macro then(description, definition)
    step_definition_(description, definition)
end

#
# Step Definition Matcher implementation
#

struct FromMacroStepDefinitionMatcher <: StepDefinitionMatcher
    stepdefinitions::Vector{StepDefinition}
    filename::String

    function FromMacroStepDefinitionMatcher(source::AbstractString; filename::String="<no filename>")
        global currentdefinitions
        global currentfilename
        currentfilename = filename
        include_string(source)
        mydefinitions = currentdefinitions
        this = new(mydefinitions, filename)
        currentdefinitions = Vector{StepDefinition}()
        this
    end
end

function findstepdefinition(matcher::FromMacroStepDefinitionMatcher, step::Gherkin.ScenarioStep)
    matchingindexes = find(x -> x.description == step.text, matcher.stepdefinitions)
    matchingstepdefinitions = [matcher.stepdefinitions[i] for i in matchingindexes]
    if isempty(matchingstepdefinitions)
        throw(NoMatchingStepDefinition())
    end
    if length(matchingstepdefinitions) > 1
        locations = map(stepdefinition -> StepDefinitionLocation(matcher.filename, 0),
                        matchingstepdefinitions)
        throw(NonUniqueStepDefinition(locations))
    end
    matchingstepdefinitions[1]
end

#
# Composite matcher
#
struct CompositeStepDefinitionMatcher <: StepDefinitionMatcher
    matchers::Vector{StepDefinitionMatcher}

    CompositeStepDefinitionMatcher(matchers...) = new([matchers...])
end

function findstepdefinition(composite::CompositeStepDefinitionMatcher, step::Gherkin.ScenarioStep)
    matching = StepDefinition[]
    nonuniquesfound = StepDefinitionLocation[]
    for m in composite.matchers
        try
            stepdefinition = findstepdefinition(m, step)
            push!(matching, stepdefinition)
        catch ex
            if ex isa NonUniqueStepDefinition
                append!(nonuniquesfound, ex.locations)
            end
        end
    end
    if length(matching) > 1 || !isempty(nonuniquesfound)
        locations = vcat(nonuniquesfound, [d.location for d in matching])
        throw(NonUniqueStepDefinition(locations))
    end
    if isempty(matching)
        throw(NoMatchingStepDefinition())
    end
    matching[1]
end


text:
module species
using utility, core
importall representation, fitness, common
export Species, convert!, genotype, SpeciesDefinition, empty_individual_collection

"""
The base type used by all species definitions.
"""
abstract SpeciesDefinition

# Load species stages.
include("species/stage.jl")

"""
Used to hold information about a given species.
"""
type Species
  stages::Dict{AbstractString, SpeciesStage}
  fitness::FitnessScheme
  genotype::SpeciesStage

  Species(st::Dict{AbstractString, SpeciesStage}, f::FitnessScheme) =
    new(st, f, root(st))
end

"""
Constructs an empty individual collection for this species.
"""
function empty_individual_collection(s::Species)
  F = uses(s.fitness)
  d = Dict{AbstractString, Any}()
  for (label, st) in s.stages
    d[label] = chromosome(st.representation)[]
  end
  return IndividualCollection{F}(F[], d)
end

"""
Returns the canonical genotype of a given species.
"""
genotype(s::Species) = s.genotype

"""
Returns the representation used by a given stage of a provided species.
"""
rep(species::Species, stage::AbstractString) =
  species.stages[stage].representation

"""

"""
__wallace__ = Base
startUp(m::Module) =
  __wallace__ = m

# Include all other components of the species module.
include("species/graph.jl")
include("species/simple.jl")
include("species/complex.jl")
end


text:
@testset "euler number" begin

    euler_numbers = ([0], [1], [1], [1,1,1,1], [-3], [-4], [0],
                        [1], [1,1], [0], [0], [-1], [0],
                        [0], [-1], [0,0,0], [0], [0,1],
                        [0,0], [0,0])

    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((BasicTopology(),), labels)
            i == 1 ? (@test isempty(t)) : (@test all(select(t,:euler₈) .== euler_numbers[i]))
        end
    end
end

@testset "holes" begin

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

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        labels = Images.label_components(img,trues(3,3))
        t = measure_components((BasicTopology(),), labels)
        @test length(t) == 1
        row₁ = t[1]
        @test row₁.holes == 3
    end
end


text:
#Julia Unit Testing

using FactCheck
f(x) = x^3
facts("cubes") do
    @fact f(2) --> 8
    @fact f(2) --> 7
end

using Base.Test
f(x) = x^3
@test f(2) == 8
@test f(2) == 7


text:
using Base64
using Random

# Will convert all relative links in th E-Mail to cid links
# and return a list of HTML-multipart sections containing the images
function encodeImages(html::String, imgPath::String)::Tuple{String, Array{String,1}}

    local pos = 1
    local plotId = 1
    local parts = []

    while true
        local range = findnext(r"src=\"[^\"]+\"", html, pos)
        range === nothing && break
        pos = range.start + 5

        # Base64-encode the image file
        local io = IOBuffer();
        local iob64_encode = Base64EncodePipe(io);
        open( joinpath(imgPath, html[range][6:end-1]) ) do file
            local str =  read(file)
            write(iob64_encode, str)
        end
        close(iob64_encode);
        local base64 = String(take!(io))

        # Generate a beautiful cid
        local cid = "plot$(plotId)"
        plotId += 1

        push!(parts, """Content-Type: image/png;
 name="file.png"
Content-Transfer-Encoding: base64
Content-ID: <$cid>
Content-Disposition: inline;
 filename="file.png"

$base64
""")

        html = replace(html, html[range]=>"src=\"cid:$cid\"", count=1)
    end

    return html, parts
end

# given a html and a path relative to which images can be found on the disk,
# generate a e-mail with given subject and recipient
function makeEmail(html::String, imgPath::String, to::String, subject::String)::String
    logger("subject=$subject, to=$to", "Generating Email")
    local html2, parts = encodeImages(html, imgPath)

    # we have two kinds of boundaries - those of the "multipart/alternative" and
    # those of the embedded "multipart/related"
    local boundary = "----------" * randstring(14)
    local boundaryA = "----------" * randstring(14)

    # insert boundary between parts
    local moreParts = ""
    for part in parts
        moreParts *= "\n--$boundary\n"
        moreParts *= part
    end

    # return the E-Mail
    return """To: $to
Subject: $subject
MIME-Version: 1.0
Content-Type: multipart/alternative;
    boundary="$boundaryA"
Content-Language: en-US

This is a multi-part message in MIME format.
--$boundaryA
Content-Type: text/plain; charset=utf-8; format=flowed
Content-Transfer-Encoding: 7bit

Message is in the HTML Section.

--$boundaryA
Content-Type: multipart/related;
    boundary="$boundary"

--$boundary
Content-Type: text/html; charset=utf-8
Content-Transfer-Encoding: 7bit

$html2

$moreParts
--$boundary--

--$boundaryA--
"""

end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
using Random: randperm
import Flux

export make_minibatches

function make_minibatches(x::AbstractMatrix, y::AbstractVecOrMat, z::Union{AbstractVecOrMat,Nothing} = nothing; n::Integer = 128)
  N = size(x,2)
  batch_idxs = Iterators.partition(randperm(N), n)
  if z isa Nothing
    [(x[:,i], y isa Matrix ? y[:,i] : y[i]) for i in batch_idxs]
  else
    [(x[:,i], y isa Matrix ? y[:,i] : y[i], z isa Matrix ? z[:,i] : z[i]) for i in batch_idxs]
  end
end


text:
export VERefValue

# Base.RefValue isn't VE compatible, so provide a compatible alternative
mutable struct VERefValue{T} <: Ref{T}
  x::T
  function VERefValue{T}(v::T) where {T}
    return new{T}(v)
  end
end

function VERefValue(v::T) where {T}
  return VERefValue{T}(v)
end

@inline Base.getindex(r::VERefValue) = r.x
@inline function Base.setindex!(r::VERefValue, x)
    r.x = x
    x
end


text:
using DiscreteMarkovChains
using Test
using LinearAlgebra

@testset "Custom eltype" begin
    struct MyType
        value
    end
    value(x::MyType) = x.value

    Base.one(::Type{MyType}) = MyType(1)
    Base.one(::MyType) = one(MyType)
    Base.zero(::Type{MyType}) = MyType(0)
    Base.zero(::MyType) = zero(MyType)

    Base.:+(x::MyType, y::MyType) = MyType(value(x)+value(y))
    Base.:*(x::MyType, y::MyType) = MyType(value(x)*value(y))
    Base.:-(x::MyType, y::MyType) = MyType(value(x)-value(y))
    Base.:/(x::MyType, y::MyType) = MyType(value(x)/value(y))
    Base.:+(x::MyType, y) = MyType(value(x)+y)
    Base.:*(x::MyType, y) = MyType(value(x)*y)
    Base.:-(x::MyType, y) = MyType(value(x)-y)
    Base.:/(x::MyType, y) = MyType(value(x)/y)

    Base.adjoint(x::MyType) = MyType(adjoint(value(x)))
    Base.abs(x::MyType) = MyType(abs(value(x)))
    Base.:<(x::MyType, y::MyType) = value(x) < value(y)
    Base.:>(x::MyType, y::MyType) = value(x) > value(y)
    Base.inv(x::MyType) = MyType(inv(value(x)))

    function Base.isapprox(x::MyType, y::MyType; atol, rtol)
        return isapprox(value(x), value(y), atol=atol, rtol=rtol)
    end
    function Base.isapprox(x::MyType, y::MyType)
        return isapprox(value(x), value(y))
    end
    function Base.rtoldefault(::Type{MyType}, ::Type{MyType}, z)
        return Base.rtoldefault(Float16, Float16, z)
    end
    Base.promote_rule(::Type{Any}, ::Type{MyType}) = MyType
    Base.convert(::Type{MyType}, x::MyType) = x
    Base.convert(::Type{MyType}, x) = MyType(x)
    Base.iterate(x::MyType) = x

    T = first.([
        ([0.5 0; 0 0.2],) ([0.5 0; 0 0.8],);
        ([0.7 0; 0 0.2],) ([0.3 0; 0 0.8],);
    ])

    T = [
        MyType(0.5) MyType(0.5);
        MyType(0.2) MyType(0.8);
    ]

    X = DiscreteMarkovChain(T)

    @test communication_classes(X) == ([[1, 2]], Any[true])
    @test is_ergodic(X)
    s = stationary_distribution(X)
    for (i, v) in enumerate(s)
        @test v ≈ [MyType(2/7), MyType(5/7)][i]
    end
    #@test stationary_distribution(X) ≈ [MyType(2/7), MyType(5/7)]

    T = first.([
        ([1 0; 0 1],) ([0 0; 0 0],);
        ([0 0; 0 0],) ([1 0; 0 1],);
    ])

    X = DiscreteMarkovChain(T)

    @test communication_classes(X) == ([[1], [2]], Any[true, true])
    @test !is_ergodic(X)
    # @test stationary_distribution(X) == T

end


text:
module TestStateSelection

using TinyModia
using Test

function checkStateSelection(model, x_names, linearEquations=[])
    # Check names of the states
    for (i, state_i) in enumerate(model.equationInfo.x_info)
        @test state_i.x_name == x_names[i]
    end

    # Check linear equation systems
    if length(linearEquations) > 0
        @test model.equationInfo.linearEquations == linearEquations
    end
end

@testset "\nTest ModiaBase/StateSelection.jl" begin

    @testset "... Test FirstOrder" begin
        FirstOrder = Model(
            T = 0.2,
            x = Var(init=0.3),
            equations = :[
                    u = 1.0
                der(x) = (u - x)/T
                    y = 2*x]
        )
        firstOrder = @instantiateModel(FirstOrder)
        checkStateSelection(firstOrder, ["x"])
    end

    @testset "... Test TwoCoupledInertias" begin
        TwoCoupledInertias = Model(
            J1_J = 1,
            J2_J = 1,
            J1_phi = Var(start=1.0),
            J1_w   = Var(start=0.0),
            J2_phi = Var(start=0.0),
            J2_w   = Var(start=0.0),
            equations = :[
                J1_w = der(J1_phi)
                J1_J * der(J1_w) = J1_tau

                J2_w = der(J2_phi)
                J2_J * der(J2_w) = J2_tau

                J2_phi = J1_phi
                J2_tau + J1_tau = 0]
        )
        twoCoupledInertias = @instantiateModel(TwoCoupledInertias)
        checkStateSelection(twoCoupledInertias, ["J1_w", "J1_phi"], [(["J1_tau"], [1], 1, 1)])
    end

    @testset "... Test ODE with linear equations 1" begin
        ODEwithLinearEquations1 = Model(
            p1=4,
            p2=2,
            x6 = Var(start=1.0),
            equations = :[
                         x1 = sin(time)
                      p2*x2 = p1*x1
                2*x3 + 4*x4 = x1
                    x3 + x5 = 2*x2
                4*x4 + 5*x5 = x1 + x2
                    der(x6) = -x5*x6]
        )
        oDEwithLinearEquations1 = @instantiateModel(ODEwithLinearEquations1, unitless=true)
        checkStateSelection(oDEwithLinearEquations1, ["x6"], [(["x5"], [1], 1, 1)])
    end

    @testset "... Test ODE with linear equations 2" begin
        ODEwithLinearEquations2 = Model(
            p1=4,
            p2=2,
            x6 = Var(start=1.0),
            equations = :[
                        x1 = sin(time)
                     p2*x2 = p1*x1
                p2*x3 + p1*x4 + 5*der(x6) = x1
                    x3 + x5 = 2*x2
                4*x4 + 5*x5 = x1 + x2
                  3*der(x6) = -5*x3 - 2*x4]
        )
        oDEwithLinearEquations2 = @instantiateModel(ODEwithLinearEquations2, unitless=true)
        checkStateSelection(oDEwithLinearEquations2, ["x6"], [(["x5"], [1], 1, 1)])
    end

    @testset "... Test Multi-index DAE" begin
        MultiIndexDAE = Model(
            x2 = Var(init = 0.0),
            x2d = Var(init = 0.0),
            equations = :[
                u1 = sin(1*time)
                u2 = sin(2*time)
                u3 = sin(3*time)
                u4 = sin(4*time)
                u5 = sin(5*time)
                u6 = sin(6*time)
                u7 = sin(7*time)
                0 = u1 + x1 - x2
                0 = u2 + x1 + x2 - x3 + der(x6)
                0 = u3 + x1 + der(x3) - x4
                x1d  = der(x1)
                x2d  = der(x2)
                x3d  = der(x3)
                x6d  = der(x6)
                x6dd = der(x6d)
                0 = u4 + 2*der(x1d) + der(x2d) + der(x3d) + der(x4) + der(x6dd)
                0 = u5 + 3*der(x1d) + 2*der(x2d) + x5
                0 = u6 + 2*x6 + x7
                0 = u7 + 3*x6 + 4*x7]
        )
        multiIndexDAE = @instantiateModel(MultiIndexDAE, unitless=true)
        #@show multiIndexDAE.equationInfo.linearEquations
        checkStateSelection(multiIndexDAE, ["x2", "x2d"],
                            [(["x7"], [1], 1, 1),
                             (["der(x7)"], [1], 1, 1),
                             (["der(der(x7))"], [1], 1, 1),
                             (["der(der(der(x7)))"], [1], 1, 1),
                             (["der(x2d)"], [1], 1, 1)])
    end

    @testset "... Test free flying mass" begin
        FreeFlyingMass = Model(
            m =1.0,
            f=[1.0,2.0,3.0],
            r = Var(init=[0.1, 0.2, 0.3]),
            v = Var(init=[-0.1, -0.2, -0.3]),
            equations = :[
                v = der(r)
                m*der(v) = f]
        )
        freeFlyingMass = @instantiateModel(FreeFlyingMass)
        checkStateSelection(freeFlyingMass, ["v", "r"], [])
    end

#=

    @testset "... Test sliding mass" begin
        SlidingMass = Model(
            m = 1.0,
            c = 1e4,
            d = 10.0,
            g = 9.81,
            n = [1.0, 1.0, 1.0],
            s = Var(init=1.0),
            equations = :[
                r = n*s
                v = der(r)
                m*der(v) = f + m*g*[0,-1,0] + u
                0 = n*f
                u = -(c*s + d*der(s))*n]
        )
        slidingMass = @instantiateModel(SlidingMass, logStateSelection=true)
        @show slidingMass.equationInfo.linearEquations
        #checkStateSelection(slidingMass, ["x6"], [(["x5"], [1], 1, 1)])
    end
=#

end

using Unitful
include("../models/AllModels.jl")

#=
Gear = Model(
    flange_a = Flange,
    flange_b = Flange,
    gear     = IdealGear_withSupport | Map(ratio = 105.0),
    fixed    = Fixed,
    spring   = Spring | Map(c=5.84e5u"N*m/rad"),
    damper1  = Damper | Map(d=500.0u"N*m*s/rad"),
    damper2  = Damper | Map(d=100.0u"N*m*s/rad"),

    connect = :[
        (flange_a     , gear.flange_a)
        (fixed.flange , gear.support, damper2.flange_b)
        (gear.flange_b, damper2.flange_a, spring.flange_a, damper1.flange_a)
        (flange_b     , spring.flange_b, damper1.flange_b)]
)

Drive = Model(
    inertia1 = Inertia,
    inertia2 = Inertia,
    gear     = Gear,
    equations = :[
        inertia1.flange_a.tau = 0u"N*m"
        inertia2.flange_b.tau = 0u"N*m"],
    connect = :[
        (inertia1.flange_b, gear.flange_a),
        (gear.flange_b    , inertia2.flange_a)
    ]
)
=#

# Linear 1D rotational damper
AbsoluteDamper = Model(
    flange = Flange,
    d = 1.0u"N*m*s/rad", # (info = "Damping constant"),
    phi = Var(init=0.0u"rad"),
    equation = :[
        phi = flange.phi
        w   = der(phi)
        flange.tau = d * w ]
)

Drive = Model(
    inertia = Inertia,
    damper  = AbsoluteDamper,
    equations = :[
        inertia.flange_a.tau = 0u"N*m"],
    connect = :[
        (inertia.flange_b, damper.flange)]
)
Drive2 = Drive | Map(damper = Map(phi=Var(init=1.0u"rad")))

drive1 = @instantiateModel(Drive)
drive2 = @instantiateModel(Drive2)

using DifferentialEquations
simulate!(drive1, Tsit5(), stopTime = 4.0)

println("Next simulate! should result in an error:\n")
simulate!(drive2, Tsit5(), stopTime = 4.0)

end


text:
prior = ones(6,6);

uniform_prior = prior / 36;

likelyhood = zeros(6,6);

#p(t=9|sa,sb)
likelyhood[3,6] = 1;
likelyhood[4,5] = 1;
likelyhood[5,4] = 1;
likelyhood[6,3] = 1;

model = likelyhood .* uniform_prior

#p(sa,sb|t=9) = p(t=9|sa,sb)p(sa)p(sb)/s(t=9)
#p(p=t) = 4/36 = 1/9

model = model / (1/9)


text:
# ImageHDU methods

# Display the image datatype and dimensions
function show(io::IO, hdu::ImageHDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    bitpix = fits_get_img_type(hdu.fitsfile)
    equivbitpix = fits_get_img_equivtype(hdu.fitsfile)
    sz = fits_get_img_size(hdu.fitsfile)

    if bitpix == equivbitpix
        datainfo = string(TYPE_FROM_BITPIX[equivbitpix])
    else
        datainfo = @sprintf("%s (physical: %s)",
                            TYPE_FROM_BITPIX[equivbitpix],
                            TYPE_FROM_BITPIX[bitpix])
    end

    print(io, """
    File: $(fits_file_name(hdu.fitsfile))
    HDU: $(hdu.ext)$(fits_get_ext_info_string(hdu.fitsfile))
    Type: Image
    Datatype: $datainfo
    Datasize: $(tuple(sz...))""")
end

# Get image dimensions
"""
    ndims(hdu::ImageHDU)

Get number of image dimensions, without reading the image into memory.
"""
function ndims(hdu::ImageHDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    fits_get_img_dim(hdu.fitsfile)
end

"""
    size(hdu::ImageHDU)
    size(hdu::ImageHDU, i)

Get image dimensions (or `i`th dimension), without reading the image
into memory.
"""
function size(hdu::ImageHDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    sz = fits_get_img_size(hdu.fitsfile)
    tuple(sz...)
end

function size(hdu::ImageHDU, i::Integer)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    sz = fits_get_img_size(hdu.fitsfile)
    sz[i]
end

"""
    length(hdu::ImageHDU)

Get total number of pixels in image (product of ``size(hdu)``).
"""
length(hdu::ImageHDU) = prod(size(hdu))

# `lastindex` is needed so that hdu[:] can throw DimensionMismatch
# when ndim != 1, rather than no method.
lastindex(hdu::ImageHDU) = length(hdu::ImageHDU)

# Read a full image from an HDU
"""
    read(hdu::ImageHDU)
    read(hdu::ImageHDU, range...)

Read the data array or a subset thereof from disk. The first form
reads the entire data array. The second form reads a slice of the array
given by the specified ranges or integers.
"""
function read(hdu::ImageHDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    sz = fits_get_img_size(hdu.fitsfile)
    bitpix = fits_get_img_equivtype(hdu.fitsfile)
    data = Array{TYPE_FROM_BITPIX[bitpix]}(undef, sz...)
    fits_read_pix(hdu.fitsfile, data)
    data
end

# _checkbounds methods copied from Julia v0.4 Base.
_checkbounds(sz, i::Integer) = 1 <= i <= sz
_checkbounds(sz, i::Colon) = true
_checkbounds(sz, r::AbstractRange{Int}) =
    (isempty(r) || (minimum(r) >= 1 && maximum(r) <= sz))

# helper functions for constructing cfitsio indexing vectors in read(hdu, ...)
_first(i::Union{Integer, AbstractRange}) = first(i)
_first(::Colon) = 1
_last(sz, i::Union{Integer, AbstractRange}) = last(i)
_last(sz, ::Colon) = sz
_step(r::AbstractRange) = step(r)
_step(::Union{Integer, Colon}) = 1

# Shape of array to create for read(hdu, ...), dropping trailing
# scalars. This is simpler than in Base because we are guaranteed that
# length(I) == length(sz).
@inline _index_shape(sz, I...) = _index_shape_dim(sz, 1, I...)
_index_shape_dim(sz, dim, I::Integer...) = ()
_index_shape_dim(sz, dim, ::Colon) = (sz[dim],)
_index_shape_dim(sz, dim, r::AbstractRange) = (length(r),)
@inline _index_shape_dim(sz, dim, ::Colon, I...) =
    tuple(sz[dim], _index_shape_dim(sz, dim+1, I...)...)
@inline _index_shape_dim(sz, dim, ::Integer, I...) =
    _index_shape_dim(sz, dim+1, I...)
@inline _index_shape_dim(sz, dim, r::AbstractRange, I...) =
    tuple(length(r), _index_shape_dim(sz, dim+1, I...)...)

# Read a subset of an ImageHDU
function read_internal(hdu::ImageHDU, I::Union{AbstractRange{Int}, Integer, Colon}...)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    sz = fits_get_img_size(hdu.fitsfile)

    # check number of indices and bounds. Note that number of indices and
    # array dimension must match, unlike in Arrays. Array-like behavior could
    # be supported in the future with care taken in constructing first, last,
    if length(I) != length(sz)
        throw(DimensionMismatch("number of indices must match dimensions"))
    end
    for i=1:length(sz)
        _checkbounds(sz[i], I[i]) || throw(BoundsError())
    end

    # construct first, last and step vectors
    firsts = Clong[_first(idx) for idx in I]
    lasts = Clong[_last(sz[i], I[i]) for i=1:length(sz)]
    steps = Clong[_step(idx) for idx in I]

    # construct output array
    bitpix = fits_get_img_equivtype(hdu.fitsfile)
    data = Array{TYPE_FROM_BITPIX[bitpix]}(undef, _index_shape(sz, I...))

    fits_read_subset(hdu.fitsfile, firsts, lasts, steps, data)
    data
end

# general method and version that returns a single value rather than 0-d array
read(hdu::ImageHDU, I::Union{AbstractRange{Int}, Int, Colon}...) =
    read_internal(hdu, I...)
read(hdu::ImageHDU, I::Int...) = read_internal(hdu, I...)[1]

"""
    write(f::FITS, data::Array; header=nothing, name=nothing, ver=nothing)

Add a new image HDU to FITS file `f` with contents `data`. The
following array element types are supported: `UInt8`, `Int8`,
`UInt16`, `Int16`, `UInt32`, `Int32`, `Int64`, `Float32`,
`Float64`. If a `FITSHeader` object is passed as the `header` keyword
argument, the header will also be added to the new HDU.
"""
function write(f::FITS, data::Array{T};
               header::Union{Nothing, FITSHeader}=nothing,
               name::Union{Nothing, String}=nothing,
               ver::Union{Nothing, Integer}=nothing) where T
    fits_assert_open(f.fitsfile)
    s = size(data)
    fits_create_img(f.fitsfile, T, [s...])
    if isa(header, FITSHeader)
        fits_write_header(f.fitsfile, header, true)
    end
    if isa(name, String)
        fits_update_key(f.fitsfile, "EXTNAME", name)
    end
    if isa(ver, Integer)
        fits_update_key(f.fitsfile, "EXTVER", ver)
    end
    fits_write_pix(f.fitsfile, ones(Int, length(s)), length(data), data)
    nothing
end

# Copy a rectangular section of an image and write it to a new FITS
# primary image or image extension. The new image HDU is appended to
# the end of the output file; all the keywords in the input image will
# be copied to the output image. The common WCS keywords will be
# updated if necessary to correspond to the coordinates of the section.

# convert a range to a string that cfitsio understands
cfitsio_range_string(r::UnitRange) = @sprintf "%d:%d" first(r) last(r)
cfitsio_range_string(r::StepRange) =
    @sprintf "%d:%d:%d" first(r) last(r) step(r)

"""
    copy_section(hdu, dest, r...)

Copy a rectangular section of an image and write it to a new FITS
primary image or image extension in `FITS` object `dest`. The new
image HDU is appended to the end of `dest`. All the keywords in the
input image will be copied to the output image. The common WCS
keywords will be updated if necessary to correspond to the coordinates
of the section.

# Examples

Copy the lower-left 200 x 200 pixel section of the image in `hdu`
to an open file, `f`

```
copy_section(hdu, f, 1:200, 1:200)
```

Same as above but only copy odd columns in y:

```
copy_section(hdu, f, 1:200, 1:2:200)
```
"""
function copy_section(hdu::ImageHDU, dest::FITS, r::AbstractRange{Int}...)
    fits_assert_open(hdu.fitsfile)
    fits_assert_open(dest.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    fits_copy_image_section(hdu.fitsfile, dest.fitsfile,
                            join([cfitsio_range_string(ri) for ri in r], ','))
end


text:
const PG_T = 0.64
const π²_8 = π^2 / 8

@doc raw"""
    PolyaGamma(b::Real, c::Real) <: ContinuousUnivariateDistribution

## Arguments
- `b::Real`
- `c::Real` exponential tilting

Create a PolyaGamma sampler with parameters `b` and `c`.
Note that sampling will differ if `b` is a `Int` or a `Real`.
"""
struct PolyaGamma{Tb,Tc} <: Distributions.ContinuousUnivariateDistribution
    b::Tb
    c::Tc
end

Base.eltype(::PolyaGamma{T,Tc}) where {T,Tc} = Tc

Distributions.params(d::PolyaGamma) = (d.b, d.c)

function Statistics.mean(d::PolyaGamma)
    if iszero(d.c)
        return d.b / 4
    else
        return d.b / (2 * d.c) * tanh(d.c / 2)
    end
end

Base.minimum(d::PolyaGamma) = zero(eltype(d))
Base.maximum(::PolyaGamma) = Inf
Distributions.insupport(::PolyaGamma, x::Real) = zero(x) <= x < Inf

function Distributions.logpdf(d::PolyaGamma, x::Real)
    b, c = Distributions.params(d)
    if iszero(b)
        return iszero(x) ? zero(x) : -Inf # The limit of PG when b->0
    # is the delta dirac at 0.
    else
        iszero(x) && return -Inf # The limit to p(x) for x-> 0 is 0.
        # valₘₐₓ = Γb - b^2 / (8x)
        ext = logtilt(x, b, c) + (b - 1) * logtwo - loggamma(b) - (log2π + 3 * log(x)) / 2
        xmax = loggamma(b) - abs2(b) / (8x)
        sumval = sum(1:201) do n
            (iseven(n) ? 1 : -1) * exp(_pdf_val_log_series(n, b, x) - xmax) * (2n + b) / b
        end
        return ext + xmax + log(b) + log(1 + sumval)
        # series = sum(1:200) do n
        #     v = 2 * n + b
        #     val = loggamma(n + b) - loggamma(n + 1) - abs2(v) / (8x)
        #     return (iseven(n) ? 1 : -1) *
        #     exp(val - valₘₐₓ)
        # end
        # return ext + log(b) + valₘₐₓ + log(series)
    end
end

function _pdf_val_log_series(n::Integer, b::Real, x)
    return loggamma(n + b) - loggamma(n + 1) - abs2(2n + b) / (8x)
end

Distributions.logpdf(d::PolyaGamma, x::NamedTuple{(:ω,),<:Tuple{<:Real}}) = logpdf(d, x.ω)

# Shortcut for computing KL(PG(ω|b,c) || PG(b, 0))
function Distributions.kldivergence(q::PolyaGamma, p::PolyaGamma)
    (q.b == p.b && iszero(p.c)) || error(
        "cannot compute the KL divergence for Polya-Gamma distributions",
        " with different parameter b, (q.b = $(q.b), p.b = $(p.b))",
        " or if p.c ≠ 0 (p.c = $(p.c))",
    )
    return logtilt(mean(q), q.b, q.c)
end

function logtilt(ω, b, c)
    return b * logcosh(c / 2) - abs2(c) * ω / 2
end

function ntrand(rng::AbstractRNG, d::PolyaGamma)
    return (; ω=rand(rng, d))
end

# Dispatch for MeasureTheory
function Base.rand(rng::AbstractRNG, ::Type{T}, d::PolyaGamma) where {T}
    return ntrand(rng, d)
end

function Distributions.rand(rng::AbstractRNG, d::PolyaGamma)
    if iszero(d.b)
        return zero(eltype(d))
    end
    return draw_sum(rng, d)
end

## Sampling when `b` is an integer
function draw_sum(rng::AbstractRNG, d::PolyaGamma{<:Integer})
    b, c = Distributions.params(d)
    return sum(1:b) do _
        sample_pg1(rng, c)
    end
end

## Sampling when `b` is a Real (and might need to be truncated)
function draw_sum(rng::AbstractRNG, d::PolyaGamma{<:Real})
    b, c = Distributions.params(d)
    if b < 1
        return rand_gamma_sum(rng, d, b)
    end
    trunc_b = floor(Int, b)
    trunc_term = sum(1:trunc_b) do _
        sample_pg1(rng, c)
    end

    res_b = b - trunc_b
    if iszero(res_b)
        return trunc_term
    else
        res_term = rand_gamma_sum(rng, d, res_b)
        return trunc_term + res_term
    end
end

# Sample ω as the series of Gamma variables (truncated at 200)
function rand_gamma_sum(rng::AbstractRNG, d::PolyaGamma, e::Real)
    inv2π² = inv2π * invπ
    w = (d.c * inv2π)^2
    ga = Gamma(e, 1)
    return inv2π² * sum(1:200) do k
        rand(rng, ga) / ((k - 0.5)^2 + w)
    end
end

## Utility functions
function a(n::Int, x::Real)
    k = (n + 1//2) * π
    if x > PG_T
        return k * exp(-k^2 * x / 2)
    elseif x > 0
        expnt = -3 / 2 * (log(halfπ) + log(x)) - 2 * (n + 1//2)^2 / x
        return k * exp(expnt)
    else
        throw(DomainError(x, "x should be a positive real"))
    end
end

function mass_texpon(z::Real, K::Real)
    t = PG_T

    b = sqrt(inv(t)) * (t * z - 1)
    a = -sqrt(inv(t)) * (t * z + 1)

    x0 = log(K) + K * t
    xb = x0 - z + logcdf(Distributions.Normal(), b)
    xa = x0 + z + logcdf(Distributions.Normal(), a)

    qdivp = fourinvπ * (exp(xb) + exp(xa))

    return 1 / (1 + qdivp)
end

# Sample from a truncated inverse gaussian
function rand_truncated_inverse_gaussian(rng::AbstractRNG, z::Real)
    μ = inv(z)
    x = one(z) + PG_T
    if μ > PG_T
        α = zero(μ)
        while α < rand(rng)
            E = randexp(rng)
            E′ = randexp(rng)
            while E^2 > (2E′ / PG_T)
                E = randexp(rng)
                E′ = randexp(rng)
            end
            x = PG_T / abs2(1 + E * PG_T)
            α = exp(-z^2 * x / 2)
        end
    else
        while (x > PG_T)
            y = randn(rng)^2
            μy = μ * y
            x = μ + μ * μy / 2 - μ * sqrt(4 * μy + abs2(μy)) / 2
            if μ / (μ + x) < rand(rng)
                x = μ^2 / x
            end
        end
    end
    return x
end

# Sample from PG(1, z)
# Algorithm 1 from "Bayesian Inference for logistic models..." p. 26
function sample_pg1(rng::AbstractRNG, c::Real)
    # Change the parameter.
    z = abs(c) / 2

    # Now sample 0.25 * J^*(1, Z := Z/2).
    if iszero(z) # We specialize on c = 0
        r = 0.5776972428360435
        K = π²_8
    else
        K = π²_8 + z^2 / 2
        r = mass_texpon(z, K)
    end
    while true
        if r > rand(rng) # sample from truncated exponential
            x = PG_T + randexp(rng) / K
        else # sample from truncated inverse Gaussian
            x = rand_truncated_inverse_gaussian(rng, z)
        end
        s = a(0, x)
        y = rand(rng) * s
        n = 0
        while true
            n += 1
            if isodd(n)
                s -= a(n, x)
                y <= s && return x / 4
            else
                s += a(n, x)
                y > s && break
            end
        end
    end
end # Sample PG(1, c)


text:
using AlgebraicPetri
using AlgebraicPetri.Epidemiology

using Catlab
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Programs.RelationalPrograms

using LabelledArrays
using OrdinaryDiffEq
using Random
using Plots

# helper function to visualize categorical representation
display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

# population x spontaneously moves to population y
spontaneous_petri(x::Symbol, y::Symbol, transition::Symbol) =
    Open(LabelledPetriNet(unique([x,y]), transition=>(x, y)))
# population y causes population x to move to population z
exposure_petri(x::Symbol, y::Symbol, z::Symbol, transition::Symbol) =
    Open(LabelledPetriNet(unique([x,y,z]), transition=>((x,y)=>(z,y))))

infection = exposure_petri(:S, :I, :I, :inf)
exposure = exposure_petri(:S, :I, :E, :exp)
illness = spontaneous_petri(:E,:I,:ill)
recovery = spontaneous_petri(:I,:R,:rec)
death = spontaneous_petri(:I,:D,:death)

epi_dict = Dict(:infection=>infection,
                :exposure=>exposure,
                :illness=>illness,
                :recovery=>recovery,
                :death=>death)

oapply_epi(ex, args...) = oapply(ex, epi_dict, args...)

Graph(infection)

Graph(recovery)

sir_wiring_diagram = @relation (s, i, r) begin
    infection(s, i)
    recovery(i, r)
end
display_uwd(sir_wiring_diagram)

sir_model = apex(oapply_epi(sir_wiring_diagram));
Graph(sir_model)

tmax = 40.0
tspan = (0.0,tmax);

u0 = LVector(S=990.0, I=10.0, R=0.0)

p = LVector(inf=0.05*10.0/sum(u0), rec=0.25); # β*c/N,γ

Random.seed!(1234);

prob_ode = ODEProblem(vectorfield(sir_model),u0,tspan,p)
sol_ode = solve(prob_ode, Tsit5());
plot(sol_ode)


text:
function main_timing(args)
    println("Début de l'action timing")

    # @error "main_timing: désolé la méthode main_timing n'est pas implantée"
    # Résolution de l'action
    println("="^70)
    println("Résolution du timing pour $(args[:infile])\n")
    inst = Instance(args[:infile])
    # args.names=[1,6,8,4,12,9,11,3,10,2,19,24,20,7,5,50,15,23,18,14,13,25,17,26,
    #          43,16,44,27,32,28,22,33,29,47,34,48,49,46,38,35,45,39,31,40,36,30,42,41,37,21]
    # args.names = %w(3 4 5 6 8 9 7 1 10 2)
    # puts inst.to_s
    if inst.nb_planes != length(args[:planes])
        println("\nERREUR taille (=$(inst.nb_planes) de l'instance différente ")
        println(" de la taille (=$(length(args[:planes])) de la solution")
        println()
        exit(1)
    end
    sol = Solution(inst, update=false)
    set_from_names!(sol, args[:planes])

    solve!(sol, do_update_cost=true)
    print_sol(sol)

    lg1("Création du fichier \"$(guess_solname(sol))\"... ")
    write(sol)

    println("Fin de l'action timing")
end

main_timing(Args.args)


text:
import Distributions: logpdf, Poisson, Normal, MvNormal

################################
# poisson convenience wrappers #
################################

function elementwise_poisson(lam)
    C, H, W = size(lam)
    samp    = Array(Float64, size(lam))
    for c in 1:C, h in 1:H, w in 1:W
        samp[c, h, w] = float(rand(Distributions.Poisson(lam[c, h, w])))
    end
    return samp
end

# output a single poisson log likelihood
function poisson_lnpdf(data, lam)
    return float(logpdf(Distributions.Poisson(lam), data))
end

#######################################################
# save images + plotting stuff for numpy/matplotlib   #
#######################################################

function save_image_array(img_stack, fname)
    npzwrite(fname, img_stack)
end

function save_blob_images(images, fname)
    imgs = permutedims(cat(3, [i.pixels for i in mod_images]...), (3, 1, 2))
    npzwrite(fname, imgs)
end

function save_plot(domain, range, fname; true_val=nothing)
    domain_name = @sprintf "%s-%s" "fin" fname
    range_name  = @sprintf "%s-%s" "fval" fname
    npzwrite(domain_name, domain)
    npzwrite(range_name, range)

    if true_val!=nothing
      tname = @sprintf "%s-%s" "ftrue" fname
      npzwrite(tname, [true_val])
    end
end

"""
Given a function (lnpdf), that takes in a NUM_COLOR_COMPONENTS dimensional vector (th),
this function fixes NUM_COLOR_COMPONENTS-1 of the dimensions of th, and returns a function
handle that only varies with respect to dimension d.
"""
function curry_dim(lnpdf::Function, th0::Vector{Float64}, dim::Int)
    thd0 = th0[dim]

    function clnpdf(thd::Float64)
        th0[dim] = thd
        return lnpdf(th0)
    end

    return clnpdf
end

"""
Function that mimics the AccuracyBenchmark module function
variational_parameters_to_data_frame_row, for direct comparison of ground
truth parameters to VB fit parameters.
"""
function catalog_to_data_frame_row(catalog_entry; objid="truth")
    ce = catalog_entry
    df = DataFrame()
    fluxes = ce.is_star ? ce.star_fluxes : ce.gal_fluxes
    df[:objid]                          = [objid]
    df[:right_ascension_deg]            = [ce.pos[1]]
    df[:declination_deg]                = [ce.pos[2]]
    df[:is_saturated]                   = [false]
    df[:is_star]                        = [ce.is_star]
    df[:de_vaucouleurs_mixture_weight]  = [ce.gal_frac_dev]
    df[:minor_major_axis_ratio]         = [ce.gal_ab]
    df[:half_light_radius_px]           = [ce.gal_scale]
    df[:angle_deg]                      = [ce.gal_angle]
    df[:reference_band_flux_nmgy]       = [fluxes[3]]
    df[:log_reference_band_flux_stderr] = [NaN]
    df[:color_log_ratio_ug]             = [log(fluxes[2]) - log(fluxes[1])]
    df[:color_log_ratio_gr]             = [log(fluxes[3]) - log(fluxes[2])]
    df[:color_log_ratio_ri]             = [log(fluxes[4]) - log(fluxes[3])]
    df[:color_log_ratio_iz]             = [log(fluxes[5]) - log(fluxes[4])]
    df[:color_log_ratio_ug_stderr]      = [NaN]
    df[:color_log_ratio_gr_stderr]      = [NaN]
    df[:color_log_ratio_ri_stderr]      = [NaN]
    df[:color_log_ratio_iz_stderr]      = [NaN]
    return df
end

function samples_to_data_frame_rows(chain; is_star=true)
    df = DataFrame()

    # reference band flux (+ log flux) and colors
    df[:log_reference_band_flux]  = chain[:,3]
    df[:reference_band_flux_nmgy] = exp.(chain[:,3])
    df[:color_log_ratio_ug] = chain[:,2] .- chain[:,1]
    df[:color_log_ratio_gr] = chain[:,3] .- chain[:,2]
    df[:color_log_ratio_ri] = chain[:,4] .- chain[:,3]
    df[:color_log_ratio_iz] = chain[:,5] .- chain[:,4]

    u_samps = [constrain_pos(chain[p,6:7]) for p in 1:size(chain, 1)]
    u_samps = transpose(hcat(u_samps...))

    df[:right_ascension_deg] = u_samps[:,1]
    df[:declination_deg]     = u_samps[:,2]

    return df
end

function samples_to_data_frame_row(sampdf; objid="mcmc")
    """ only for stars right now """
    df = DataFrame()
    df[:objid]                          = [objid]
    df[:right_ascension_deg]            = [mean(sampdf[:right_ascension_deg])]
    df[:declination_deg]                = [mean(sampdf[:declination_deg])]
    df[:is_saturated]                   = [false]
    df[:is_star]                        = [true]
    df[:de_vaucouleurs_mixture_weight]  = [NaN]
    df[:minor_major_axis_ratio]         = [NaN]
    df[:half_light_radius_px]           = [NaN]
    df[:angle_deg]                      = [NaN]
    df[:reference_band_flux_nmgy]       = [mean(sampdf[:reference_band_flux_nmgy])]
    df[:log_reference_band_flux_stderr] = [std(sampdf[:log_reference_band_flux])]
    df[:color_log_ratio_ug]             = [mean(sampdf[:color_log_ratio_ug])]
    df[:color_log_ratio_gr]             = [mean(sampdf[:color_log_ratio_gr])]
    df[:color_log_ratio_ri]             = [mean(sampdf[:color_log_ratio_ri])]
    df[:color_log_ratio_iz]             = [mean(sampdf[:color_log_ratio_iz])]
    df[:color_log_ratio_ug_stderr]      = [std(sampdf[:color_log_ratio_ug])]
    df[:color_log_ratio_gr_stderr]      = [std(sampdf[:color_log_ratio_gr])]
    df[:color_log_ratio_ri_stderr]      = [std(sampdf[:color_log_ratio_ri])]
    df[:color_log_ratio_iz_stderr]      = [std(sampdf[:color_log_ratio_iz])]

    return df
end

# run multiple chains
function run_multi_chain_mcmc(lnpdf::Function, th0s)
  chains, logprobs = [], []
  for c in 1:num_chains
    # initialize a star state (around the existing catalog entry)
    println("---- chain ", c, " ----")
    th0 = init_star_params(star_state)
    samples, lls =
      run_mh_sampler_with_warmup(star_logpdf, th0,
                                 num_samples,
                                 num_warmup;
                                 print_skip=print_skip,
                                 keep_warmup=true,
                                 warmup_prop_scale=warmup_prop_scale,
                                 chain_prop_scale=chain_prop_scale)

    push!(chains, samples[(num_warmup+1):end, :])
    push!(logprobs, lls[(num_warmup+1):end])

    # report PSRF after c = 3 chains
    if c > 2
      psrfs = MCMC.potential_scale_reduction_factor(chains)
      println(" potential scale red factor", psrfs)
    end
  end
end

######################################################################
# unit flux image utils --- writes unit flux images to pixel array   #
# Very similar to Synthetic.jl functions                             #
######################################################################

function write_star_unit_flux(img0::Image,
                              pos::Array{Float64, 1},
                              pixels::Matrix{Float64})
    iota = median(img0.nelec_per_nmgy)
    for k in 1:length(img0.psf)
        the_mean = SVector{2}(WCS.world_to_pix(img0.wcs, pos)) + img0.psf[k].xiBar
        the_cov = img0.psf[k].tauBar
        intensity = iota * img0.psf[k].alphaBar
        write_gaussian(the_mean, the_cov, intensity, pixels)
    end
end

function write_galaxy_unit_flux(img0::Image,
                                pos::Array{Float64,1},
                                gal_frac_dev::Float64,
                                gal_ab::Float64,
                                gal_angle::Float64,
                                gal_scale::Float64,
                                pixels::Matrix{Float64})
    iota = median(img0.nelec_per_nmgy)
    gal_fracdevs = [gal_frac_dev, 1 - gal_frac_dev]
    #XiXi = DeterministicVI.get_bvn_cov(ce.gal_ab, ce.gal_angle, ce.gal_scale)
    XiXi = Model.get_bvn_cov(gal_ab, gal_angle, gal_scale)

    for i in 1:2
        for gproto in galaxy_prototypes[i]
            for k in 1:length(img0.psf)
                the_mean = SVector{2}(WCS.world_to_pix(img0.wcs, pos)) +
                           img0.psf[k].xiBar
                the_cov = img0.psf[k].tauBar + gproto.nuBar * XiXi
                intensity = iota * img0.psf[k].alphaBar * gal_fracdevs[i] *
                    gproto.etaBar
                write_gaussian(the_mean, the_cov, intensity, pixels)
            end
        end
    end
end

function write_gaussian(the_mean, the_cov, intensity, pixels)
    the_precision = inv(the_cov)
    c = sqrt(det(the_precision)) / 2pi

    H, W = size(pixels)
    w_range, h_range = get_patch(the_mean, H, W)

    for w in w_range, h in h_range
        y = @SVector [the_mean[1] - h, the_mean[2] - w] # Maybe not hard code Float64
        ypy = dot(y,  the_precision * y)
        pdf_hw = c * exp(-0.5 * ypy)
        pixel_rate = intensity * pdf_hw
        pixels[h, w] += pixel_rate
    end
    pixels
end

function get_patch(the_mean::SVector{2,Float64}, H::Int, W::Int)
    radius = 50
    hm = round(Int, the_mean[1])
    wm = round(Int, the_mean[2])
    w11 = max(1, wm - radius):min(W, wm + radius)
    h11 = max(1, hm - radius):min(H, hm + radius)
    return (w11, h11)
end

"""
Potential Scale Reduction Factor --- Followed the formula from the
following website:
http://blog.stata.com/2016/05/26/gelman-rubin-convergence-diagnostic-using-multiple-chains/
"""
function potential_scale_reduction_factor(chains)
    # each chain has to be size N x NUM_COLOR_COMPONENTS, we have M chains
    N, NUM_COLOR_COMPONENTS  = size(chains[1])
    M     = length(chains)

    # mean and variance of each chain
    means = vcat([mean(s, 1) for s in chains]...)
    vars  = vcat([var(s, 1)  for s in chains]...)

    # grand mean
    gmu   = mean(means, 1)

    # between chain variance:w
    NUM_BANDS = float(N)/(float(M)-1)*sum( broadcast(-, means, gmu).^2, 1)

    # average within chain variance
    W = mean(vars, 1)

    # compute PRSF ratio
    Vhat = (float(N)-1.)/float(N) * W + (float(M)+1)/float(N*M) * NUM_BANDS
    psrf = Vhat ./ W
    return psrf
end

"""
Convert `chains` (list of sample arrays), `logprobs` (list of
log-likelihood traces) and `colnames` (list of parameter names for chains)
into a DataFrame object --- for saving and plotting
"""
function chains_to_dataframe(chains, logprobs, colnames)
    # sample, lls, chain identifier all in arrays
    samples    = vcat(chains...)
    sample_lls = vcat(logprobs...)
    chain_id   = vcat([s*ones(Integer, size(chains[s], 1))
                       for s in 1:length(chains)]...)

    # convert to data frame, rename variables
    sdf = convert(DataFrame, samples)
    rename!(sdf, names(sdf), [Symbol(s) for s in colnames])

    # add lls, chain index
    sdf[:lls]   = sample_lls
    sdf[:chain] = chain_id
    return sdf
end


text:
module problem105

using problem103

description = """
Special subset sums: testing
Problem 105

Let S(A) represent the sum of elements in set A of size n. We shall call it a special sum set if for any two non-empty disjoint subsets, B and C, the following properties are true:

    S(B) ≠ S(C); that is, sums of subsets cannot be equal.
    If B contains more elements than C then S(B) > S(C).

For example, {81, 88, 75, 42, 87, 84, 86, 65} is not a special sum set because 65 + 87 + 88 = 75 + 81 + 84, whereas {157, 150, 164, 119, 79, 159, 161, 139, 158} satisfies both rules for all possible subset pair combinations and S(A) = 1286.

Using sets.txt (right click and "Save Link/Target As..."), a 4K text file with one-hundred sets containing seven to twelve elements (the two examples given above are the first two sets in the file), identify all the special sum sets, A1, A2, ..., Ak, and find the value of S(A1) + S(A2) + ... + S(Ak).

NOTE: This problem is related to problems 103 and 106.
"""

function sum_of_special_sets()
  total = 0
  sets = read_sets()
  for s in sets
    if is_special_set(s)
      total += sum(s)
    end
  end
  return total
end

function read_sets()
  return map(s -> int(split(s, ',')), open(readlines, "sets.txt"))
end

function is_special_set(s :: Vector{Int})
  sort!(s)
  sp = SpecialSet(s[end-1:end], [s[end], s[end-1], s[end]+s[end-1]])
  for x in s[end-2:-1:1]
    if can_add(x, sp)
      sp = add(x, sp)
    else
      return false
    end
  end
  return true
end

using Base.Test

@test is_special_set([157, 150, 164, 119, 79, 159, 161, 139, 158])
@test !is_special_set([81, 88, 75, 42, 87, 84, 86, 65])

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using StatisticalRethinking, OrderedCollections

df_covid = DataFrame(
    :completed => [true, true, false, false],
    :company => [:pfizer, :moderna, :jj, :astra],
    :unv => [18261, 15000, missing, missing],
    :vac => [18262, 15000, missing, missing],
    :cov_unv => [162, 185, missing, missing],
    :cov_vac => [8, 11, missing, missing]
)

function extract_company(df::DataFrame, c::Symbol)
    cv = df[df.company .== c, :cov_vac]
    ncv = df[df.company .== c, :vac] - cv
    v = df[df.company .== c, :vac]
    cu = df[df.company .== c, :cov_unv]
    ncu = df[df.company .== c, :unv] - cu
    u = df[df.company .== c, :unv]
    NamedArray(
        [cv ncv v; cu ncu u; cv+cu ncv+ncu v+u],
    (   [:vaccinated, :unvaccinated, :total],
        [:covid, :noninfected, :total]
    ), ("Rows", "Cols"))
end

na_pfizer = extract_company(df_covid, :pfizer)
na_pfizer |> display
println()

# Risk of covid among vaccinated

@show risk_vac = na_pfizer[:vaccinated, :covid] / na_pfizer[:vaccinated, :total]

# Risk of covid among unvaccinated

@show risk_unvac = na_pfizer[:unvaccinated, :covid] / na_pfizer[:unvaccinated, :total]
@show rr_pfizer = risk_vac / risk_unvac
@show ve_pfizer = (risk_unvac - risk_vac) /risk_unvac

# Alternatively, using risk ratio

@show 1 - rr_pfizer

# So, the vaccinated group experienced 95% fewer covid cases
# than they would have if they had not been vaccinated.
println()

na_moderna = extract_company(df_covid, :moderna)
na_moderna |> display
println()

# Risk of covid among vaccinated

@show risk_vac = na_moderna[:vaccinated, :covid] / na_moderna[:vaccinated, :total]

# Risk of covid among unvaccinated

@show risk_unvac = na_moderna[:unvaccinated, :covid] / na_moderna[:unvaccinated, :total]
@show rr_moderna = risk_vac / risk_unvac
@show ve_moderna = (risk_unvac - risk_vac) /risk_unvac

# Alternatively, using risk ratio

@show 1 - rr_moderna

# So, the vaccinated group experienced 95% fewer covid cases
# than they would have if they had not been vaccinated.
println()


text: const StorageAccountType = String


text:
function expected_improvement(x::AbstractVector,
                              y_opt::Real,
                              X::Matrix,
                              K::Array{PDMats.PDMat},
                              a::Matrix,
                              k::Array{KernelFunctions.Kernel})
    N  = length(K)
    ei = zeros(N)
    @simd for i = 1:N
        μ, σ² = gp_predict(x, X, K[i], a[:,i], k[i])
        σ     = sqrt(σ²)
        Δy    = μ - y_opt
        z     = Δy / σ
        ei[i] = Δy * normcdf(z) + σ*normpdf(z)
    end
    mean(ei)
end

function expected_improvement(x::AbstractVector,
                              y_opt::Real,
                              X::Matrix,
                              K::PDMats.PDMat,
                              a::Vector,
                              k::KernelFunctions.Kernel)
    μ, σ² = gp_predict(x, X, K, a, k)
    σ     = sqrt(σ²)
    Δy    = μ - y_opt
    z     = Δy / σ
    ei    = Δy * normcdf(z) + σ*normpdf(z)
    ei
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    PlainHyperGraph,
    plainhypergraph

struct PlainHyperGraph <: OrderedHyperGraph
    incidences::Vector{Set{Int}}
    edges::Vector{Set{Int}}
    cache::Dict{Symbol,Any}
end

"""
    plainhypergraph() -> PlainHyperGraph

Generate empty `PlainGraph`.
"""
plainhypergraph() = PlainHyperGraph([], [], Dict())

"""
    plainhypergraph(size, edges) -> PlainHyperGraph

Generate hyper graph that have given size (number of nodes) and edges (Sets of
nodes).
"""
function plainhypergraph(size::Int, edges)
    incidences = [Set{Int}() for i in 1:size]
    edges = collect(edges)
    for (i, ns) in enumerate(edges)
        for n in ns
            push!(incidences[n], i)
        end
    end
    return PlainHyperGraph(incidences, edges, Dict{Symbol,Any}())
end

"""
    plainhypergraph(graph::OrderedHyperGraph) -> PlainHyperGraph

Convert an arbitrary `OrderedHyperGraph` into a `PlainHyperGraph`.
"""
function plainhypergraph(graph::OrderedHyperGraph)
    newg = plainhypergraph()
    for incset in graph.incidences
        push!(newg.incidences, copy(incset))
    end
    for eset in graph.edges
        push!(newg.edges, copy(eset))
    end
    return newg
end

"""
    plainhypergraph(graph::HyperGraph) -> PlainHyperGraph

Convert an arbitrary `HyperGraph` into a `PlainHyperGraph`.

If the given graph has non-sequential indices (ex. subgraph view), node and
edge indices are sorted in ascending order and are re-indexed.
"""
function plainhypergraph(graph::HyperGraph)
    newg = plainhypergraph()
    nkeys = sort(collect(nodeset(graph)))
    ekeys = sort(collect(edgeset(graph)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.incidences, Set())
    end
    for (i, e) in enumerate(ekeys)
        ns = Set{Int}()
        for oldn in getedge(graph, e)
            push!(ns, nmap[oldn])
            push!(newg.incidences[nmap[oldn]], i)
        end
        push!(newg.edges, ns)
    end
    return newg
end