text:
struct Rain_Gauge
    name::String
    latitude::Float64
    longitude::Float64
    precipitation::Array{Float64}
end

in_path = "Resources/WBNM2017_v001/WBNM 1998/190917_FCT_CS08a_STR02_August1998_RF11.wbn"
out_path = "data/rain/1998.csv"

file = open(in_path)
file_content = read(file, String)

precipitation_data = match(r"#####START_RECORDED_RAIN(?s).*#####END_RECORDED_RAIN", file_content).match

gauges = [x.match for x in eachmatch(r".*\n\s+\d{4,7}\.\d{2}\s+\d{4,7}\.\d{2}.*\n(?:.*\d{1,3}\.\d{2}\r\n)*", precipitation_data)]
num_gauge = length(gauges)

cells = Dict()
for gauge=1:num_gauge
    gauge_capture = gauges[gauge]
    name = match(r".*\n", gauge_capture).match
    latitude, longitude = [parse(Float64, x) for x in match(r"(?<latitude>\d{4,7}\.\d{2})\s+(?<longitude>\d{4,7}\.\d{2})", gauge_capture).captures]
    values = [parse(Float32, x.match) for x in eachmatch(r"\D+\d{1,3}\.\d{2}", gauge_capture)]
    cells[gauge] = Rain_Gauge(name, latitude, longitude, values)
end

csv_file = open(out_path, "w")

num_steps = length(cells[1].precipitation)

column_names = ["Event code", "Start time", "ID"]

append!(column_names, ["t"*string(t-1) for t in 1:num_steps])

write(csv_file, join(column_names, ",")*"\n")

for gauge=1:num_gauge
    write(csv_file, "common,common,$gauge,"*join(cells[gauge].precipitation, ",")*"\n")
end

close(csv_file)

max_stars_repo_name: smart-facility/smart-faware

max_issues_repo_name: smart-facility/smart-faware

max_forks_repo_name: smart-facility/smart-faware


text:
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/model_fit_global.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/models.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/helper_functions.jl")

using DataFrames

CD3_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/dissociations.csv",DataFrame)

delete_missing(x) = [collect(skipmissing(x[:,i])) for i in 1:ncol(x)]

ge = delete_missing(CD3_data[:,2:3])
de = delete_missing(CD3_data[:,8:9])
mix = delete_missing(CD3_data[:,14:15])
gethrml = delete_missing(CD3_data[:,5:5])

scatter(ge[1],ge[2])
scatter(de[1],de[2])
scatter(mix[1],mix[2])

# for each experiment: [forces, lifetimes]

function bin_dissociations(data, n_bins; bin_type = "force")
    F = data[1]
    t = data[2]
    # even force bins
    if bin_type == "force"
        f_min = minimum(F)
        f_max = maximum(F)

        edges = range(f_min,stop=f_max,length=n_bins+1)
        # h = fit(Histogram, F)
        bins,bin_ids = partition(F,edges)
        binned_t = [similar(t,0) for i in bins]
        for (i,bin) in enumerate(bin_ids)
            push!(binned_t[bin],t[i])
        end

    end
    return bins, bin_ids,binned_t
end

function partition(data, intervals)
    ranges = intervals[1:end-1] .=> intervals[2:end]
    bins = [similar(data, 0) for _ in 1:length(ranges)]
    ids = similar(data,0)
    for x in data
        for (i, (a, b)) in pairs(ranges)
            if a <= x < b
                push!(bins[i], x)
                push!(ids,i)
                break
            end
        end
    end
    return bins, Int64.(ids)
end

f_binned, bin_ids, binned_t = bin_dissociations(ge,8)
using StatsBase
t_avg = mean.(binned_t)
f_avg = mean.(f_binned)
plot!(f_avg,t_avg)
histogram(binned_t[5])
    # even events number bins

binned_t_sorted = sort.(binned_t)
sort_id = sortperm.(binned_t)

[f_binned[i][sort_id[i]] for i in 1:length(sort_id)]

survival_p(t) = collect((length(t) :-1:1) / length(t))

ps = (survival_p.(binned_t_sorted))

p2lnp(p) = map((x) -> log.(x), p)
# lnps = map((x) -> log.(x), ps)
plot(binned_t_sorted,lnps)

gethrml_sorted = sort(gethrml)
p_gethrml = survival_p.(gethrml_sorted)
lnp_gethrml = p2lnp(p_gethrml)
plot(gethrml_sorted,lnp_gethrml)

using StatsBase, LinearAlgebra
x = randn(1000)*10
h = fit(Histogram, x)
h.edges
h.weights

normalize(h, mode=:pdf)

plot(h)

ext: jl


text:
"Mod functions for CartesianIndex"
@inline function modc(i::CartesianIndex{D}, m) :: CartesianIndex{D} where {D}
    CartesianIndex(mod.(Tuple(i), Tuple(m)))
end
@inline function modc1(i::CartesianIndex{D}, m) :: CartesianIndex{D} where {D}
    CartesianIndex(mod1.(Tuple(i), Tuple(m)))
end
@inline function offset(i::CartesianIndex{D}, n::SVector{D,Int}, m) :: CartesianIndex{D} where {D}
    CartesianIndex(Tuple(mod1.(Tuple(i) .+ n, Tuple(m))))
end
"Splits a CartesianIndex into its first index, and the rest"
@inline function splitidx(i::CartesianIndex{D}) where {D}
    return (i[1], CartesianIndex(Tuple(i)[2:end]))
end

# Taken from:
# https://discourse.julialang.org/t/efficient-tuple-concatenation/5398/8
"Functions for joining tuples"
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

# For efficiency, may need to look into Base.unsafe_wrap
#   and pointer trickery if we want to stick with Vec3.

"Reinterprets an array of Vec3 to an equivalent array of Float64"
@inline function _reinterpret_from_spin_array(A::Array{Vec3}) :: Array{Float64}
    Ar = reinterpret(reshape, Float64, A)
end

"Reinterprets an array of Floats with leading dimension 3 to an array of Vec3"
@inline function _reinterpret_to_spin_array(A::Array{Float64}) :: Array{Vec3}
    Ar = reinterpret(reshape, Vec3, A)
end

"Reinterprets an array of Mat3 to an equivalent array of Float64"
@inline function _reinterpret_dipole_tensor(A::OffsetArray{Mat3}) :: Array{Float64}
    Ar = reinterpret(reshape, Float64, parent(A))
    return reshape(Ar, 3, 3, size(A)...)
end


text:
"infinite vector with non-trivial values in discrete indexes. For other indexes has a default value."
struct DiscreteInfiniteVector{T} <: AbstractInfiniteVector{T}
    values::Dict{Int, T}
    default::T
    DiscreteInfiniteVector{T}(values::Dict{Int, T}, default::T) where {T} = new(values, default)
end
DiscreteInfiniteVector(values::Dict{Int, T}, default::T) where {T} = DiscreteInfiniteVector{T}(values, default)
DiscreteInfiniteVector(values::Dict{Int, T}) where {T <: Number} = DiscreteInfiniteVector(values, zero(T))
DiscreteInfiniteVector(default::T) where {T} = DiscreteInfiniteVector(Dict{Int, T}(), default)

Base.getindex(x::DiscreteInfiniteVector, i::Int) = get(x.values, i, x.default)
Base.getindex(x::DiscreteInfiniteVector, I) = [x[i] for i in I]
Base.setindex!(x::DiscreteInfiniteVector, v, i) = x.values[i] = v


text:
using ApproxFun, LinearAlgebra, Test
import ApproxFun: resizedata!, CachedOperator, RaggedMatrix, testbandedblockbandedoperator,
                    testblockbandedoperator, ldiv_coefficients, mul_coefficients,
                    factor

@testset "Full PDE" begin
    @testset "Gaussian zero Dirichlet Poisson" begin
        S=JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ=Laplacian(S)

        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),rangespace(Δ))  #default is [-1,1]^2
        @time v = \(Δ,f;tolerance=1E-14)
        @test norm((Δ*v-f).coefficients) < 1E-14

    end

    @testset "Neumann" begin
        d=ChebyshevInterval()^2

        @time u=\([Neumann(d); Laplacian(d)-100.0I],[[[1,1],[1,1]],0.];tolerance=1E-12)
        @test u(.1,.9) ≈ 0.03679861429138079
    end

    @testset "Transport" begin
        dx=ChebyshevInterval(); dt=Interval(0,2.)
        d=dx × dt
        Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
        x,y=Fun(identity,d)
        @time u=\([I⊗ldirichlet(dt);Dt+x*Dx],[Fun(x->exp(-20x^2),dx);0.];tolerance=1E-12)

        @test u(0.1,0.2) ≈ 0.8745340845783758  # empirical
    end

    @testset "Bilaplacian" begin
        dx=dy=ChebyshevInterval()
        d=dx × dy
        Dx=Derivative(dx);Dy=Derivative(dy)
        L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

        testbandedblockbandedoperator(L)

        A=[ldirichlet(dx)⊗Operator(I,dy);
                rdirichlet(dx)⊗Operator(I,dy);
                Operator(I,dx)⊗ldirichlet(dy);
                Operator(I,dx)⊗rdirichlet(dy);
                lneumann(dx)⊗Operator(I,dy);
                rneumann(dx)⊗Operator(I,dy);
                Operator(I,dx)⊗lneumann(dy);
                Operator(I,dx)⊗rneumann(dy);
                 L]

        # Checks bug in constructor
        f=Fun((x,y)->real(exp(x+1.0im*y)),component(rangespace(A)[1],1),22)
        @test f(-1.,0.1) ≈ real(exp(-1+0.1im))
        f=Fun((x,y)->real(exp(x+1.0im*y)),component(rangespace(A)[1],1))
        @test f(-1.,0.1) ≈ real(exp(-1+0.1im))

        F=[Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[1]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[2]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[3]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[4]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[5]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[6]);
            Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A)[7]);
            Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A)[8]);
            0]

        @time u=\(A,F;tolerance=1E-10)

        @test u(0.1,0.2)  ≈ exp(0.1)*cos(0.2)

        dx=dy=ChebyshevInterval()
        d=dx×dy
        Dx=Derivative(dx);Dy=Derivative(dy)
        L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

        A=[(ldirichlet(dx)+lneumann(dx))⊗Operator(I,dy);
                (rdirichlet(dx)+rneumann(dx))⊗Operator(I,dy);
                Operator(I,dx)⊗(ldirichlet(dy)+lneumann(dy));
                Operator(I,dx)⊗(rdirichlet(dy)+rneumann(dy));
                (ldirichlet(dx)-lneumann(dx))⊗Operator(I,dy);
                (rdirichlet(dx)-rneumann(dx))⊗Operator(I,dy);
                Operator(I,dx)⊗(ldirichlet(dy)-lneumann(dy));
                Operator(I,dx)⊗(rdirichlet(dy)-rneumann(dy));
                 L]

        u=\(A,[fill(1.0,8);0];tolerance=1E-5)
        @test u(0.1,0.2) ≈ 1.0

        F=[2Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[1]);
            2Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[2]);
            Fun((x,y)->real(exp(x+1.0im*y))-imag(exp(x+1.0im*y)),rangespace(A)[3]);
            Fun((x,y)->real(exp(x+1.0im*y))-imag(exp(x+1.0im*y)),rangespace(A)[4]);
            0;
            0;
            Fun((x,y)->real(exp(x+1.0im*y))+imag(exp(x+1.0im*y)),rangespace(A)[7]);
            Fun((x,y)->real(exp(x+1.0im*y))+imag(exp(x+1.0im*y)),rangespace(A)[8]);
            0]

        u=\(A,F;tolerance=1E-10)

        @test u(0.1,0.2)  ≈ exp(0.1)*cos(0.2)
    end

    @testset "Operator resize" begin
        S=ChebyshevDirichlet()^2
        B=Dirichlet(S)
        f = Fun((x,y)->exp(x)*sin(y),S)
        @test norm((Fun((x,y)->exp(x)*sin(y),∂(domain(S))) - B*f).coefficients) < 100eps()

        S=JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ=Laplacian(S)

        @test cache(Δ)[1:100,1:100]  ≈ Δ[1:100,1:100]
        @test cache(Δ;padding=true)[1:100,1:100]  ≈ Δ[1:100,1:100]

        @test cache(Δ)[5:100,7:100]  ≈ Δ[5:100,7:100]
        @test cache(Δ;padding=true)[5:100,7:100]  ≈ Δ[5:100,7:100]

        # Check that QR is growing correctly
        for col in (1,2,3,10,11,40)
            QR=qr(Δ)
            resizedata!(QR.R_cache,:,col+100)
            resizedata!(QR,:,col)
            QR2=qr!(CachedOperator(RaggedMatrix,Δ;padding=true))
            resizedata!(QR2.R_cache,:,QR.ncols+100)
            resizedata!(QR2,:,QR.ncols)
            n=min(size(QR.H,1),size(QR2.H,1))
            @test QR.H[1:n,1:col] ≈ QR2.H[1:n,1:col]
            @test QR.R_cache[1:col,1:col] ≈ QR2.R_cache[1:col,1:col]
            @test QR.R_cache[1:col+10,1:col+10] ≈ QR2.R_cache[1:col+10,1:col+10]
        end

        QRR=qr(Δ)
        QR2=qr!(CachedOperator(RaggedMatrix,Δ;padding=true))
        for col in (80,200)
            resizedata!(QRR,:,col)
            resizedata!(QR2,:,QRR.ncols)
            n=min(size(QRR.H,1),size(QR2.H,1))
            @test QRR.H[1:n,1:col] ≈ QR2.H[1:n,1:col]
            @test QRR.R_cache[1:col,1:col] ≈ QR2.R_cache[1:col,1:col]
            @test QRR.R_cache[1:col+10,1:col+10] ≈ QR2.R_cache[1:col+10,1:col+10]
        end

        # this checks a bug
        QRR=qr(Δ)
        resizedata!(QRR,:,548)
        resizedata!(QRR,:,430)

        u=Fun((x,y)->sin(π*x)*sin(π*y),S)
        f=-2π^2*u

        QRR=qr(Δ)
        v=QRR\f
        @test norm((u-v).coefficients)<100eps()

        v=Δ\f
        @test norm((u-v).coefficients)<100eps()

        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),rangespace(Δ))  #default is [-1,1]^2
        @time v=\(Δ,f;tolerance=1E-14)
        @test norm((Δ*v-f).coefficients)<1E-14

        KO=Δ.op.ops[1].ops[1].op

        M=ApproxFun.BandedBlockBandedMatrix(view(KO,1:4,1:4))
        @test norm(ApproxFun.BandedBlockBandedMatrix(view(KO,1:4,2:4))-M[:,2:4]) < 10eps()
        @test norm(ApproxFun.BandedBlockBandedMatrix(view(KO,1:4,3:4))-M[:,3:4]) < 10eps()
    end

    @testset "Check resizing" begin
        d=ChebyshevInterval()^2
        A=[Dirichlet(d);Laplacian()+100I]
        QRR = qr(A)
        @time ApproxFun.resizedata!(QRR.R_cache,:,2000)
        @test norm(QRR.R_cache.data[1:200,1:200] - A[1:200,1:200]) ==0

        @time ApproxFun.resizedata!(QRR,:,200)
        j=56
        v=QRR.R_cache.op[1:100,j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=300)[j+1:end]) < 100eps()

        j=195
        v=QRR.R_cache.op[1:ApproxFun.colstop(QRR.R_cache.op,j),j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=1000)[j+1:end]) < 100eps()

        j=300
        v=QRR.R_cache.op[1:ApproxFun.colstop(QRR.R_cache.op,j),j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=1000)[j+1:end]) < j*20eps()

        @test ApproxFun.colstop(QRR.R_cache.op,195)-194 == ApproxFun.colstop(QRR.H,195)

        QR1 = qr(A)
        @time ApproxFun.resizedata!(QR1.R_cache,:,1000)
        QR2 = qr([Dirichlet(d);Laplacian()+100I])
        @time ApproxFun.resizedata!(QR2.R_cache,:,500)
        n=450;QR1.R_cache.data[1:n,1:n]-QR2.R_cache.data[1:n,1:n]|>norm
        @time ApproxFun.resizedata!(QR2.R_cache,:,1000)
        N=450;QR1.R_cache.data[1:N,1:N]-QR2.R_cache.data[1:N,1:N]|>norm
        N=1000;QR1.R_cache.data[1:N,1:N]-QR2.R_cache.data[1:N,1:N]|>norm

        QR1 = qr(A)
        @time ApproxFun.resizedata!(QR1,:,1000)
        QR2 = qr([Dirichlet(d);Laplacian()+100I])
        @time ApproxFun.resizedata!(QR2,:,500)
        @time ApproxFun.resizedata!(QR2,:,1000)

        @test norm(QR1.H[1:225,1:1000]-QR2.H[1:225,1:1000]) ≤ 10eps()

        QR1 = qr(A)
        @time ApproxFun.resizedata!(QR1,:,5000)
        @time u=\(QR1,[ones(∂(d));0.];tolerance=1E-7)

        @test norm((Dirichlet(d)*u-ones(∂(d))).coefficients) < 1E-7
        @test norm((A*u-Fun([ones(∂(d));0.])).coefficients) < 1E-7
        @test norm(((A*u)[2]-(Laplacian(space(u))+100I)*u).coefficients) < 1E-10
        @test eltype(ApproxFun.promotedomainspace(Laplacian(),space(u))) == Float64
        @test eltype(ApproxFun.promotedomainspace(Laplacian()+100I,space(u))) == Float64
        @test norm(((A*u)[2]-(Laplacian()+100I)*u).coefficients) < 1E-10
        @test norm((Laplacian()*u+100*u - (A*u)[2]).coefficients) < 10E-10
        @time v=\(A,[ones(∂(d));0.];tolerance=1E-7)
        @test norm((u-v).coefficients) < 100eps()

        @test u(0.1,1.) ≈ 1.0
        @test u(0.1,-1.) ≈ 1.0
        @test u(1.,0.1) ≈ 1.0
        @test u(-1.,0.1) ≈ 1.0
    end

    @testset "Chebyshev Dirichlet Kronecker" begin
        S=ChebyshevDirichlet()^2
        ff=(x,y)->exp(x)*cos(y)
        u=Fun(ff,S)

        for KO in [Operator(I,factor(S,1))⊗rdirichlet(factor(S,1)),
                    rdirichlet(factor(S,1))⊗Operator(I,factor(S,2))]
            testblockbandedoperator(KO)
            @test norm((KO*u-Fun(ff,rangespace(KO))).coefficients) ≤ 1E-10
        end
    end

    @testset "Periodic Poisson" begin
        d=PeriodicSegment()^2
        S=Space(d)

        f=Fun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
        A=Laplacian(d)+0.1I
        testbandedblockbandedoperator(A)
        @time u=A\f
        @test u(.1,.2) ≈ u(.2,.1)
        @test (lap(u)+.1u-f)|>coefficients|>norm < 1000000eps()
    end

    @testset "periodic x interval" begin
        dθ=PeriodicSegment(-2.,2.);dt=Interval(0,1.)
        d=dθ×dt
        Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
        u0=Fun(θ->exp(-20θ^2),dθ)

        ε = 0.1
        @time u=\([I⊗ldirichlet(dt);Dt-ε*Dθ^2-Dθ],[u0;0.];tolerance=1E-4)
        @test ≈(u(0.1,0.2),0.3103472600253807;atol=1E-2)

        A=Dt+Dθ
        testbandedblockbandedoperator(A)

        @time u=\([I⊗ldirichlet(dt);Dt+Dθ],[u0;0.0];tolerance=1E-6)
        @test ≈(u(0.2,0.1),u0(0.1);atol=1E-6)
    end

    @testset "Small diffusion" begin
        dx=ChebyshevInterval();dt=Interval(0,0.2)
        d=dx×dt
        Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
        x,t=Fun(dx×dt)

        B=0.0
        C=0.0
        V=B+C*x
        ε=0.1
        f=Fun(x->exp(-30x^2),dx)
        u=\([timedirichlet(d);Dt-ε*Dx^2-V*Dx],[f;zeros(3)];tolerance=1E-6)

        @test u(.1,.2) ≈ 0.496524222625512
        B=0.1
        C=0.2
        V=B+C*x
        u=\([timedirichlet(d);Dt-ε*Dx^2-V*Dx],[f;zeros(3)];tolerance=1E-7)
        @test u(.1,.2) ≈ 0.46810331039791464
    end

    @testset "Laplace in a strip" begin
        d=PeriodicSegment() × ChebyshevInterval()
        g=Fun((x,y)->real(cos(x+im*y)),∂(d))  # boundary data

        @test g(0.1,1.0) ≈ real(cos(0.1+im))
        @test g(0.1,-1.0) ≈ real(cos(0.1-im))
        v=[g;0]
        @test v(0.1,-1) ≈ [real(cos(0.1-im));0]

        A=[Dirichlet(d);Laplacian(d)]
        a = space(v)
        b = rangespace(A)

        @test Fun(component(v[1],1), component(b[1],1))(0.1,-1.0) ≈ v(0.1,-1.0)[1]
        @test Fun(component(v[1],2), component(b[1],2))(0.1,-1.0) ≈ v(0.1,-1.0)[1]
        @test ApproxFun.default_Fun(v[1] , b[1])(0.1,1.0) ≈ v(0.1,1.0)[1]
    end

    @testset "concatenate InterlaceOperator" begin
        a=Fun(x -> 0 ≤ x ≤ 0.5 ? 0.5 : 1, Domain(-1..1) \ [0,0.5])
        @test a(0.1) == 0.5
        @test a(0.7) == 1.0
        s=space(a)
        # Bx=[ldirichlet(s);continuity(s,0)]
        # TODO: this should concat
        dt=Interval(0,2.)
        Dx=Derivative(s);Dt=Derivative(dt)
        Bx=[ldirichlet(s);continuity(s,0)]

        @test ApproxFun.rangetype(rangespace(continuity(s,0))) == Vector{Float64}
        @test ApproxFun.rangetype(rangespace(Bx)) == Vector{Any}
        @test ApproxFun.rangetype(rangespace(Bx⊗Operator(I,Chebyshev()))) == Vector{Any}

        rhs = Fun([0,[0,[0,0]],0],rangespace([I⊗ldirichlet(dt);Bx⊗I;I⊗Dt+(a*Dx)⊗I]))
        @test rhs(-0.5,0.0) == [0,[0,[0,0]],0]

        u=\([I⊗ldirichlet(dt);Bx⊗I;I⊗Dt+(a*Dx)⊗I],
            [Fun(x->exp(-20(x+0.5)^2),s),[[0],[0,0]],0.0];tolerance=1E-2)

        @test u(-0.4,0.1) ≈ u(-0.5,0.0) atol = 0.0001
    end
end


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

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

function _swaggerinternal_usergroupsUsersList(_api::UsergroupsUsersApi, usergroup::String; include_disabled=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/usergroups.users.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "include_disabled", include_disabled)  # type Bool
    Swagger.set_param(_ctx.query, "usergroup", usergroup)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List all users in a User Group
Param: usergroup::String (required)
Param: include_disabled::Bool
Return: Dict{String, Any}
"""
function usergroupsUsersList(_api::UsergroupsUsersApi, usergroup::String; include_disabled=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usergroupsUsersList(_api, usergroup; include_disabled=include_disabled, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usergroupsUsersList(_api::UsergroupsUsersApi, response_stream::Channel, usergroup::String; include_disabled=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usergroupsUsersList(_api, usergroup; include_disabled=include_disabled, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usergroupsUsersUpdate(_api::UsergroupsUsersApi, usergroup::String, users::String; include_count=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/usergroups.users.update", ["slackAuth"])
    Swagger.set_param(_ctx.form, "include_count", include_count)  # type Bool
    Swagger.set_param(_ctx.form, "usergroup", usergroup)  # type String
    Swagger.set_param(_ctx.form, "users", users)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Update the list of users for a User Group
Param: usergroup::String (required)
Param: users::String (required)
Param: include_count::Bool
Return: Dict{String, Any}
"""
function usergroupsUsersUpdate(_api::UsergroupsUsersApi, usergroup::String, users::String; include_count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usergroupsUsersUpdate(_api, usergroup, users; include_count=include_count, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usergroupsUsersUpdate(_api::UsergroupsUsersApi, response_stream::Channel, usergroup::String, users::String; include_count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usergroupsUsersUpdate(_api, usergroup, users; include_count=include_count, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export usergroupsUsersList, usergroupsUsersUpdate


text:
"""
  LeastSquaresProblem(residuals::VectorObjective)

A LeastSquaresProblem is used to represent the
mathematical problem of finding the smallest possible sum of squared values of
the elements of the residual function. The problem is defined by `residual` which is
a `VectorObjective`.

Options are stored in `options` and are of the `NLsqOptions` type. See more information
about options using `?NLsqOptions`.

The package NLSolversAD.jl adds automatic conversion of problems to match algorithms
that require higher order derivates than provided by the user. It also adds AD
constructors for a target number of derivatives.
"""
struct LeastSquaresProblem{R, B, M, C}
  residuals::R
  bounds::B
  manifold::M
  constraints::C
end
function LeastSquaresProblem(residuals, bounds = nothing)
  manifold = Euclidean(0)
  LeastSquaresProblem{typeof(residuals), typeof(bounds), typeof(manifold), Nothing}(residuals, bounds, manifold, nothing)
end

struct LeastSquaresOptions{Tmi}
  maxiter::Tmi
end

# observations are useful if you only want to specify it as a ScalarObjective
function solve(problem::LeastSquaresProblem, x, method, options::LeastSquaresOptions)
  _error = ArgumentError("solve not implemented for $(summary(method)) and LeastSquaresProblem.")
  throw(_error)
end


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

mutable struct SubscriptionPolicies <: SwaggerModel
    locationPlacementId::Any # spec type: Union{ Nothing, String } # spec name: locationPlacementId
    quotaId::Any # spec type: Union{ Nothing, String } # spec name: quotaId
    spendingLimit::Any # spec type: Union{ Nothing, String } # spec name: spendingLimit

    function SubscriptionPolicies(;locationPlacementId=nothing, quotaId=nothing, spendingLimit=nothing)
        o = new()
        validate_property(SubscriptionPolicies, Symbol("locationPlacementId"), locationPlacementId)
        setfield!(o, Symbol("locationPlacementId"), locationPlacementId)
        validate_property(SubscriptionPolicies, Symbol("quotaId"), quotaId)
        setfield!(o, Symbol("quotaId"), quotaId)
        validate_property(SubscriptionPolicies, Symbol("spendingLimit"), spendingLimit)
        setfield!(o, Symbol("spendingLimit"), spendingLimit)
        o
    end
end # type SubscriptionPolicies

const _property_map_SubscriptionPolicies = Dict{Symbol,Symbol}(Symbol("locationPlacementId")=>Symbol("locationPlacementId"), Symbol("quotaId")=>Symbol("quotaId"), Symbol("spendingLimit")=>Symbol("spendingLimit"))
const _property_types_SubscriptionPolicies = Dict{Symbol,String}(Symbol("locationPlacementId")=>"String", Symbol("quotaId")=>"String", Symbol("spendingLimit")=>"String")
Base.propertynames(::Type{ SubscriptionPolicies }) = collect(keys(_property_map_SubscriptionPolicies))
Swagger.property_type(::Type{ SubscriptionPolicies }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SubscriptionPolicies[name]))}
Swagger.field_name(::Type{ SubscriptionPolicies }, property_name::Symbol) =  _property_map_SubscriptionPolicies[property_name]

const _allowed_SubscriptionPolicies_spendingLimit = ["On", "Off", "CurrentPeriodOff"]

function check_required(o::SubscriptionPolicies)
    true
end

function validate_property(::Type{ SubscriptionPolicies }, name::Symbol, val)
    if name === Symbol("spendingLimit")
        Swagger.validate_param(name, "SubscriptionPolicies", :enum, val, _allowed_SubscriptionPolicies_spendingLimit)
    end
end


text:
#=
    Figure 9
=#

using Inference
using Model

using JLD2
using Plots
using StatsFuns
using StatsPlots
using Random

include("../figure_defaults.jl")

##############################################################
## Load model and results
##############################################################

include("../../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

##############################################################
## "Noiseless" model: Q̄ at t = 10.0 and 120.0
##############################################################

model2 = θ -> simulate_model_noiseless(θ,[10.0,120.0],param_dist;N=1000,η)

##############################################################
## Figure 7(a)
##############################################################

Z = model2(θ)
ρ = cor([norminvcdf.(u) for u in cdf.(Z)]...)
fig7a = scatter(cdf.(Z)...,
            α      = 0.3,
            msw    = 0.0,
            c      = :black,
            legend = :none,
            box    = :on,
            xlabel = "F_20",
            ylabel = "F_120"
        )
annotate!(fig7a,0.55,0.1,"ρ = $(round(ρ,sigdigits=3))")

##############################################################
## Figure 7(b)
##############################################################

Θ = sample(C[1:100:end],1000)
ρ = [cor([norminvcdf.(u) for u in cdf.(model2(θ))]...) for θ in Θ]
fig7b = stephist(ρ,normalize=:pdf,frange=0.0,
            legend = :none,
            grid   = :x,
            axis   = :x,
            yticks = [],
            xlabel = "cor(Q̄(10),Q̄(120))",
            c      = :grey,
            widen  = false,
            xlim   = (0.49,1.01)
        )

##############################################################
## Figure 7(c)
##############################################################

idx = ρ .> 0.9
Θfilt = Θ[idx]

plts = [(   stephist(hcat(Θ...)[i,:],xlabel=param_names[i],frange=0.0,c=:grey,normalize=:pdf,label="All",α=0.8);
            stephist!(hcat(Θfilt...)[i,:],frange=0.0,c=:purple,normalize=:pdf,α=0.8,label="ρ > 0.9")
) for i = [10,11]]

fig7c = plot(plts...,xlim=(0.0,0.1),layout=grid(2,1),axis=:x,yticks=[],widen=false)

##############################################################
## Figure 7
##############################################################

fig7 = plot(fig7a,fig7b,fig7c,layout=grid(1,3),size=(700,250),bottom_margin=5Plots.mm)
add_plot_labels!(fig7)

##############################################################
## Save...
##############################################################

savefig(fig7,"$(@__DIR__)/fig7.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
    print_value(io, value)

This is the default translation of interpolated values within rawtext
tags, such as `<style>` and attribute values.

* The elements of a `Tuple` or `AbstractArray` object are printed,
  with a space between each item.

* The `Pair`, `NamedTuple`, and `Dict` objects are treated as if
  they are CSS style elements, with a colon between key and value,
  each pair delimited by a semi-colon.

* The `Nothing` object is treated as an empty string.

Otherwise, this method simply uses the standard `print` representation
for the given object.
"""
print_value(io::IO, @nospecialize value) =
    print(io, value)
print_value(io::IO, ::Nothing) =
    nothing

function print_value(io::IO, xs::Union{Tuple, AbstractArray, Base.Generator})
    prior = false
    for x in xs
        if prior
            print(io, " ")
        end
        print_value(io, x)
        prior = true
    end
end

function print_pairs(io, xs)
    prior = false
    for (key, value) in xs
        name = normalize_attribute_name(key)
        if prior
            print(io, "; ")
        end
        print(io, name)
        print(io, ": ")
        print_value(io, value)
        prior = true
    end
    print(io, ";")
end

print_value(io::IO, pair::Pair) = print_pairs(io, (pair,))
print_value(io::IO, items::Dict) = print_pairs(io, items)
print_value(io::IO, items::NamedTuple) = print_pairs(io, pairs(items))
print_value(io::IO, items::Tuple{Pair, Vararg{Pair}}) = print_pairs(io, items)

"""
    attribute_value(x)

This method may be implemented to specify a printed representation
suitable for use within a quoted attribute value.
"""
attribute_value(x::String) = x
attribute_value(x::Number) = x
attribute_value(x::Symbol) = x

mutable struct AttributeValue
    content::Any
end

Base.print(ep::EscapeProxy, x::AttributeValue) =
    print_value(ep, x.content)

attribute_value(@nospecialize x) = AttributeValue(x)

"""
    content(x)

This method may be implemented to specify a printed representation
suitable for `text/html` output. `AbstractString`, `Symbol` and `Number`
(including `Bool`) types are printed, with proper escaping.

A default implementation first looks to see if `typeof(x)` has
implemented a way to show themselves as `text/html`, if so, this is
used. Otherwise, the result is printed within a `<span>` tag, using a
`class` that includes the module and type name. Hence, `missing` is
serialized as: `<span class="Base-Missing">missing</span>`.
"""
@generated function content(x)
     if hasmethod(show, Tuple{IO, MIME{Symbol("text/html")}, x})
         return :(Render(x))
     else
         mod = parentmodule(x)
         cls = string(nameof(x))
         if mod == Core || mod == Base || pathof(mod) !== nothing
             cls = join(fullname(mod), "-") * "-" * cls
         end
         span = """"""
         return :(reprint(Bypass($span), x, Bypass("")))
     end
end

function reprint(xs...)
    # generated functions cannot contain a closure
    Reprint() do io::IO
        for x in xs
            print(io, x)
        end
    end
end

content(x::Union{AbstractString, Symbol}) = x
content(x::Nothing) = ""
content(x::Union{AbstractFloat, Bool, Integer}) = x
content(xs...) = content(xs)

function content(xs::Union{Tuple, AbstractArray, Base.Generator})
    Reprint() do io::IO
        for x in xs
            print(io, content(x))
        end
    end
end

#-------------------------------------------------------------------------
"""
    attribute_pair(name, value)

Wrap and escape attribute name and pair within a single-quoted context
so that it is `showable("text/html")`. It's assumed that the attribute
name has already been normalized.

If an attribute value is `Bool` or `Nothing`, then special treatment is
provided. If the value is `false` or `nothing` then the entire pair is
not printed.  If the value is `true` than an empty string is produced.
"""

no_content = Reprint(io::IO -> nothing)

function attribute_pair(name, value)
    Reprint() do io::IO
        print(io, " ")
        print(io, name)
        print(io, Bypass("='"))
        print(io, attribute_value(value))
        print(io, Bypass("'"))
    end
end

function attribute_pair(name, value::Bool)
    if value == false
        return no_content
    end
    Reprint() do io::IO
        print(io, " ")
        print(io, name)
        print(io, Bypass("=''"))
    end
end

attribute_pair(name, value::Nothing) = no_content

"""
    inside_tag(value)

Convert Julian object into a serialization of attribute pairs,
`showable` via `MIME"text/html"`. The default implementation of this
delegates value construction of each pair to `attribute_pair()`.
"""

function inside_tag(value::Pair)
    name = normalize_attribute_name(value.first)
    return attribute_pair(name, value.second)
end

function inside_tag(value::Union{AbstractString, Symbol})
    name = normalize_attribute_name(value)
    return attribute_pair(name, "")
end

function inside_tag(xs::AbstractDict)
    Reprint() do io::IO
        for (key, value) in xs
            name = normalize_attribute_name(key)
            print(io, attribute_pair(name, value))
        end
    end
end

inside_tag(values::NamedTuple) =
    inside_tag(pairs(values))

inside_tag(::Nothing) = no_content


text:
using ReversibleTuringMachine
using NiLang

prog = quadruples([
    1 BLANK BLANK 2;
    2 SLASH RIGHT 3;
    3 0     1     4;
    3 1     0     2;
    3 BLANK BLANK 4;
    4 SLASH LEFT  5;
    5 0     0     4;
    5 BLANK BLANK 6;
])

rtm = RTM(1,6,prog)
tape = [BLANK, 1, 1, 0, 0, 1, BLANK]

function read_tape(tape)
    sum(tape[2:end-1] .* (1 .<< (0:length(tape)-3)))
end

run!(rtm, tape, 1, 1; visualize=true)

function main(rtm, tape)
    loc = 1
    pc = 1
    println("marching forward")
    for i=1:50
        @instr run!(rtm, tape, loc, pc)
        println("number → $(read_tape(tape))")
    end

    println("marching backward")
    for i=1:50
        @instr (~run!)(rtm, tape, loc, pc)
        println("number → $(read_tape(tape))")
    end
end


text:
#julia4 --machinefile=$PBS_NODEFILE -- run100.jl 1 2 110
#run on all available cores using:
#addprocs(readcsv(ENV["PBS_NODEFILE"], ASCIIString)[:, 1])
@everywhere args = @fetchfrom 1 ARGS

include(joinpath(Pkg.dir("LatentGaussianMixtureModel"), "challenge/gamma_ifC1.jl"))
Xraw=readcsv(joinpath(Pkg.dir("LatentGaussianMixtureModel"), "challenge/X.csv"))
#groupindex_raw=convert(Vector{UInt32},readcsv(joinpath(Pkg.dir("LatentGaussianMixtureModel"), "challenge/groupindex.csv"), Int)[:,1])
@everywhere X = @fetchfrom 1 Xraw
@everywhere groupindex = @fetchfrom 1 groupindex_raw

teststat = pmap((b) -> Brun(b, X, groupindex), 1:parse(Int, args[1]))
writecsv("gamma_ifC1.csv", teststat)
run(`sed -i "s/\[//g" gamma_ifC1.csv`)
run(`sed -i "s/\]//g" gamma_ifC1.csv`)


text:
using Plots
using Statistics
using StatsBase
using LaTeXStrings
using PWS

system = GaussianSystem(delta_t=0.05, duration=2.0)
initial = generate_configuration(system)

block(arr) = 0.5 .* (arr[begin:2:end-1] .+ arr[begin+1:2:end])

function plot_block_averages!(p, values; kwargs...)
    a = values

    y_vals = Float64[]
    y_err = (Float64[], Float64[])

    while length(a) > 4
        normed_variance = var(a, corrected=false) / (length(a) - 1)
        variance_err = sqrt(2*normed_variance^2 / (length(a) - 1))

        push!(y_vals, sqrt(normed_variance))
        push!(y_err[1], -sqrt(max(normed_variance - variance_err, 0)) + sqrt(normed_variance))
        push!(y_err[2], sqrt(normed_variance + variance_err) - sqrt(normed_variance))

        a = block(a)
    end

    block_size = 2 .^ (eachindex(y_vals) .- 1)
    plot!(p, block_size, y_vals, yerror=y_err; kwargs...)
end

signal = PWS.sample(initial, system)
energies_list = []
num_samples_list = [16, 20, 22]
for num_samples ∈ num_samples_list
    chain = PWS.chain(system, θ=1.0, scale=0.1)
    sampler = PWS.MetropolisSampler(signal, chain, burn_in=2^14)
    samples = PWS.sample(x->PWS.energy(x, chain), sampler, 2^num_samples)
    push!(energies_list, samples)
end

pyplot()
p = plot()
for (e, num_samples) ∈ zip(energies_list, num_samples_list)
    plot_block_averages!(p, e, label=L"N = 2^ {%$num_samples}", xscale=:log2, size=(400, 200))
end
xlabel!(p, L"block size $= N/N_\mathrm{blocks}$")
ylabel!(p, L"\sigma^2 / N_\mathrm{blocks}")

savefig(p, projectdir("plots", "mcmc_correlations.png"))
savefig(p, projectdir("plots", "mcmc_correlations.tex"))


text:
## Real-to-complex and complex-to-real transforms.

"""
    RFFT()

Real-to-complex FFT.

See also
[`AbstractFFTs.rfft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.rfft).
"""
struct RFFT <: AbstractTransform end

"""
    BRFFT()

Unnormalised inverse of [`RFFT`](@ref).

To obtain the inverse transform, divide the output by the length of the
transformed dimension (of the real output array).

See also
[`AbstractFFTs.brfft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.brfft).
"""
struct BRFFT <: AbstractTransform end

const TransformR2C = RFFT
const TransformC2R = BRFFT

is_inplace(::Union{TransformR2C, TransformC2R}) = false

length_output(::TransformR2C, length_in::Integer) = div(length_in, 2) + 1
length_output(::TransformC2R, length_in::Integer) = 2 * length_in - 2

eltype_output(::TransformR2C, ::Type{T}) where {T <: FFTReal} = Complex{T}
eltype_output(::TransformC2R, ::Type{Complex{T}}) where {T <: FFTReal} = T

eltype_input(::TransformR2C, ::Type{T}) where {T <: FFTReal} = T
eltype_input(::TransformC2R, ::Type{T}) where {T <: FFTReal} = Complex{T}

plan(::RFFT, args...; kwargs...) = FFTW.plan_rfft(args...; kwargs...)

# NOTE: unlike most FFTW plans, this function also requires the length `d` of
# the transform output along the first transformed dimension.
plan(::BRFFT, args...; kwargs...) = FFTW.plan_brfft(args...; kwargs...)

binv(::RFFT) = BRFFT()
binv(::BRFFT) = RFFT()

# Note: the output of RFFT (BRFFT) is complex (real).
scale_factor(::BRFFT, A::RealArray, dims) = _prod_dims(A, dims)
scale_factor(::RFFT, A::ComplexArray, dims) = _scale_factor_r2c(A, dims...)

function _scale_factor_r2c(A::ComplexArray, dim1, dims...)
    # I need to normalise by the *logical* size of the real output.
    # We assume that the dimension `dim1` is the dimension with Hermitian
    # symmetry.
    s = size(A)
    2 * (s[dim1] - 1) * _intprod((s[i] for i in dims)...)
end

# r2c along the first dimension, then c2c for the other dimensions.
expand_dims(::RFFT, ::Val{N}) where N =
    N === 0 ? () : (RFFT(), expand_dims(FFT(), Val(N - 1))...)
expand_dims(::BRFFT, ::Val{N}) where N =
    N === 0 ? () : (BRFFT(), expand_dims(BFFT(), Val(N - 1))...)


text:
# from a suggestion at https://discourse.julialang.org/t/install-a-packages-from-a-list/30920/3
# These are additional dependencies found when running the code

using Pkg

dependencies = [
    "DataFrames",
    "CSV"
]

Pkg.add(dependencies)


text:
using CSV, DataFrames, StanSample, CairoMakie, Jedi, Colors, Statistics, StatsBase

CairoMakie.activate!()
# Set Style
Jedi.styles.default_makie!()

##
# Find home directory
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-4]...)

# Find date
workdir = split(dir, '/')[end]
DATE = parse(Int64, split(workdir, '_')[1])
RUN_NO = parse(Int64, split(workdir, '_')[2][end])
##

# Read output
data_df = CSV.read("$dir/output/$(DATE)_r$(RUN_NO)_growth_plate.csv", DataFrame)

# Import Stan file
gp_stan = open("/$home_dir/fit_seq/stan_code/gp_growth_rate_prior_deriv.stan") do file
    read(file, String)
end

exp_stan = open("/$home_dir/fit_seq/stan_code/exponential_model.stan") do file
    read(file, String)
end
##
# Compile stan code
println("Compiling Stan files...")
if ~isdir("$dir/stan")
    mkdir("$dir/stan")
end
sm_gp = SampleModel("gaussian_process", gp_stan, "$dir/stan")
sm_exp = SampleModel("exponential_growth", exp_stan, "$dir/stan")

function run_stan_gp(x, y, sm)
    # Set data dictionary
    data = Dict(
        "N" => length(x),
        "t"=> x,
        "y"=> y,
        "N_predict"=> length(x),
        "t_predict"=> x,
        "alpha_param"=> [0, 1],
        "sigma_param"=> [0, 1],
        "rho_param"=> [0, 500],
    )

    # Run stan
    stan_sample(sm; data, num_chains=4);
    # Import results
    df = read_samples(sm, :dataframe)
    # Run summary
    stan_summary(sm)
    # Import summary
    summary = read_summary(sm)
    return df, summary
end

function run_stan_exp(x, y, sm)
    data = Dict(
        "N" => length(x),
        "t"=> x,
        "y"=> log.(y),
        "N_ppc"=> length(x),
        "t_ppc"=> x
    )

    # Run stan
    stan_sample(sm; data, num_chains=4);
    # Import results
    df = read_samples(sm, :dataframe)
    # Run summary
    stan_summary(sm)
    # Import summary
    summary = read_summary(sm)
    return df, summary
end

function plot_samples_gp(df, data, fig, i, title=nothing)
    # Extract variable names
    y_gp = filter(x -> String(x)[1] == 'y', names(df))
    dy_gp = filter(x -> String(x)[1:2] == "dy", names(df))

    # Prepare axes
    ax1 = Axis(fig[i, 1])
    ax1.xlabel = "time [min]"
    ax1.ylabel = "OD 600"

    ax2 = Axis(fig[i, 2])
    ax2.xlabel = "time [min]"
    ax2.ylabel = "d/dt OD600"

    if ~isnothing(title)
        ax1.title = title
    end
    ax1 = Jedi.viz.predictive_regression(
        [df[!, x] for x in y_gp],
        data[1],
        ax1,
        data=data
        )

    ax2 = Jedi.viz.predictive_regression(
        [df[!, x] for x in dy_gp],
        data[1],
        ax2,
        )

    return fig
end

function plot_samples_exp(df, data, fig, i, title=nothing)
    # Extract variable names
    y_ppc = filter(x -> occursin("y_ppc", String(x)), names(df))

    # Prepare axes
    ax1 = Axis(fig[i, 1])
    ax1.xlabel = "time [min]"
    ax1.ylabel = "OD 600"

    if ~isnothing(title)
        ax1.title = title
    end
    ax1 = Jedi.viz.predictive_regression(
        [df[!, x] for x in y_ppc],
        data[1],
        ax1,
        data=data
        )

    return fig
end

function analyze_plate(data_df::DataFrames.DataFrame)
    # Get wells
    wells = data_df.well |> unique
    #wells = wells[[1, 50]]

    # Define plotting canvas
    fig_gp = Figure(resolution=(600, 300*length(wells)))
    fig_exp = Figure(resolution=(300, 300*length(wells)))

    # Define DataFrames
    return_df = DataFrames.DataFrame()
    return_sum_df = DataFrames.DataFrame()

    println("Analyzing data...")
    for (i, well) in enumerate(wells)
        println(" Running well $well...")
        # Choosing data for well
        sub_df = data_df[data_df.well .== well, :]
        x = sub_df[!, "time_min"]
        y = sub_df[!, "OD600"]

        println("  Running Gaussian Process...")
        # Run Gaussian process
        df_gp, sum_gp = run_stan_gp(x, y, sm_gp)
        if sum_gp[sum_gp.parameters .== :divergent__, "mean"][1] != 0
            println("There were divergences! $(sum_gp[sum_gp.parameters .== :divergent__, "mean"][1])")
        end

        # Find columns with predictions
        y_gp = filter(x -> String(x)[1] == 'y', names(df_gp))
        dy_gp = filter(x -> String(x)[1:2] == "dy", names(df_gp))

        # Compute mean and variance of each step
        mean_y = mean.(eachcol(df_gp[!, y_gp]))
        mean_dy = mean.(eachcol(df_gp[!, dy_gp]))
        var_y = var.(eachcol(df_gp[!, y_gp]))
        var_dy = var.(eachcol(df_gp[!, dy_gp]))

        # Compute maximum growth rate
        lambda_max = maximum(mean_dy ./ mean_y)

        # Add information to dataframe
        insertcols!(
            sub_df,
            :mean_OD_gp=>mean_y,
            :var_OD_gp=>var_y,
            :mean_lambda_gp=>mean_dy,
            :var_lambda_gp=>var_dy,
        )
        append!(return_df, sub_df)

        # Create dataframe for maximum growth rate
        _df = DataFrames.DataFrame(
            strain=sub_df.strain |> unique,
            pos_selection=sub_df.pos_selection |> unique,
            #neg_selection=sub_df.neg_selection |> unique,
            well=well,
            lambda_max=lambda_max,
        )

        # Run Exponential model
        println("  Running Exponential Growth Model...")

        ind1 = findfirst(x -> x > 0.1, y)
        ind2 = findfirst(x -> x > 0.4, y)
        x_exp = x[ind1:ind2]
        y_exp = y[ind1:ind2]
        df_exp, summary_exp = run_stan_exp(x_exp, y_exp, sm_exp)
        if summary_exp[summary_exp.parameters .== :divergent__, "mean"][1] != 0
            println("There were divergences! $(summary_exp[summary_exp.parameters .== :divergent__, "mean"][1])")
        end

        insertcols!(
            _df,
            :exp_growth_rate=>mean(df_exp[!, "lambda"]),
        )
        append!(return_sum_df, _df)
        fig_gp = plot_samples_gp(df_gp, [x, y], fig_gp, i, "Gaussian Process: $well")
        fig_exp = plot_samples_exp(df_exp, [x_exp, y_exp], fig_exp, i, "Exponential Growth: $well")
    end
    return return_df, return_sum_df, fig_gp, fig_exp

end

return_df, return_sum_df, fig_gp, fig_exp = analyze_plate(data_df)

CSV.write("$dir/output/gp_analysis.csv", return_df)
CSV.write("$dir/output/gp_analysis_summary.csv", return_sum_df)
save("$dir/output/gp_analysis.pdf", fig_gp)
save("$dir/output/exp_model_analysis.pdf", fig_exp)
println("Done!")


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

using Markdown
using InteractiveUtils

# ╔═╡ 9b06e82a-83d9-11ec-2bec-45900ccb18c2
begin
	import Pkg
	Pkg.activate(Base.current_project())
	Pkg.add("NNlib")
	Pkg.add("Distributions")
	# Pkg.add("DiffEq")
	using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC, MCMCChains
	using JLD, StatsPlots
	using BSON
	using Random
	using Glucose: prepare_all_data_meals_hypo
	using NNlib
	using Glucose
	using LinearAlgebra
	using NNlib
	using Distributions
	using Omega
	using DiffEq
	using OrdinaryDiffEq: solve
end

# ╔═╡ dc3cc130-aa7d-417e-ac66-b3b38d4a0036
# prepare data
_, ode_data = prepare_all_data_meals_hypo(10)

# ╔═╡ 0dd8222c-a031-405e-b1fb-523b4fae908f
_, ode_data_unnormal = prepare_all_data_meals_hypo(10, normal=false)

# ╔═╡ 28566c39-a010-4f61-864d-4f052b1cb784
glucose_mean, glucose_std = (mean(ode_data_unnormal[1, :]), std(ode_data_unnormal[1, :])) # (159.44942528735626, 54.68441830006641)

# ╔═╡ ab98112d-6ddc-4b51-aaed-c46711b745e1
steps_mean, steps_std = (mean(ode_data_unnormal[2, :]), std(ode_data_unnormal[2, :])) # (2.1413793103448278, 2.478409230791181)

# ╔═╡ 746fdb88-587f-417b-bbf9-6dca06110730
bolus_mean, bolus_std = (mean(ode_data_unnormal[3, :]), std(ode_data_unnormal[3, :])) # (0.030344827586206897, 0.12658940258683044)

# ╔═╡ 2cb2b71b-e327-49a4-852c-84768a379326
meal_mean, meal_std = (mean(ode_data_unnormal[4, :]), std(ode_data_unnormal[4, :])) # (0.8620689655172413, 1.8987615875887272)

# ╔═╡ 95ff457e-a43e-4429-8289-7a86fbe44b72
begin
	if isfile("new_model_hypo2.bson")
		d = BSON.load("new_model_hypo2.bson")
		theta = d[:thetas]
	else
		theta, losses, prob = optim_all_vars_exo(10, log_dir="")
	end
end
# FIXME: Generate data if not there: extract optimal weights
# d = BSON.load("new_model_hypo2.bson")

# ╔═╡ 87467485-7e2f-454b-95f6-89996703890d
"Make time t ∈ [l, r] to index 1:n"
function scaleidx(t, lb, ub, n)
  δ = ub - lb
  Int(round(t*(n - 1)/δ + (1 - (n - 1)*lb/δ)))
end

# ╔═╡ 6582679c-f481-4359-bb31-89eb23cd361b
function no_intervention(f, u, p, t, tspan, exo_data)
  datasize = size(ode_data)[2]
  scaled_t = scaleidx(t, tspan[1], tspan[2], datasize)
  f(vcat(u, exo_data[:, scaled_t]...), p)
end

# ╔═╡ 4ebd0515-b8c5-4069-83bd-a8bd17a4b90f
"Generate the network"
function gennet()
  f = FastChain(FastDense(4, 64, swish),
                FastDense(64, 32, swish),
                FastDense(32, 2))
end

# ╔═╡ 99397fa7-0cf3-447a-b03e-711514b5a505
function predict_intervene(prob_intervene, θ, tsteps)
  _prob = remake(prob_intervene, p=θ)
  res = solve(_prob, Tsit5(), saveat=tsteps, abstol = 1e-9, reltol = 1e-9)
  Array(res)
end

# ╔═╡ 5bfbf974-b7b5-4aac-ae8c-3c2b5d61eb81
function generate_series(f,
                         u0,
                         tspan,
                         tsteps,
                         weights::AbstractArray=theta;
                         exo_data,
                         intervene = (f, u, p, t) -> no_intervention(f, u, p, t, tspan, exo_data))
  # intervened model
  function f_intervene(u, p, t)
    intervene(f, u, p, t)
  end

  prob_intervene = ODEProblem(f_intervene, u0, tspan, nothing)
  Array(predict_intervene(prob_intervene, weights, tsteps))
  # Array(vcat(predict_intervene(prob_intervene, weights, tsteps), exo_data))
end

# ╔═╡ 6bc373f4-1964-4191-a283-735277652912
function unnormalize_data(data)
  Array(transpose(hcat(
    data[1,:] .* glucose_std .+ glucose_mean,
    data[2, :] .* steps_std .+ steps_mean,
    data[3, :] .* bolus_std .+ bolus_mean,
    data[4, :] .* meal_std .+ meal_mean
  )))
end

# ╔═╡ 5ef4dd73-79a3-49c3-95b0-936865aaff58

# ╔═╡ 615ca691-a2c1-423b-8199-f66b1d07cd19

# ╔═╡ ec7bdeff-b05d-4092-ac3e-bc334c1b5d4d
I = constant(no_intervention)

# ╔═╡ 8fbaec56-660d-4148-bd89-1b81f0f382a3
Data(ω) = exo_data = ode_data[3:4, :]

# ╔═╡ 79d0e2cc-fea3-42b6-b8bd-e4b36fee2bd9
"""Example"""
function examplep(ω)
  network = gennet()
  exo_data = ode_data[3:4, :]
  tspan = (0.0, 1.0)
  tsteps =  range(tspan[1], tspan[2], length = 500)
  datatsteps = range(tspan[1], tspan[2], length=size(ode_data)[2])
  datasize = size(ode_data)[2]
  u0 = ode_data[1:2, 1]

  # Generate new exogenous data
  new_exo_data = deepcopy(exo_data)
  new_exo_data[2, 8] = 3.5 # 13, 3.5
  new_exo_data[2, 9] = 3.5 # 14, 3.5
  new_exo_data[2, 10] = 3.5 # 15, 3.5
  function f_int(f, u, p, t)
    scaled_t = scaleidx(t, tspan[1], tspan[2], datasize)
      # scaled_t = Int(round(t*(datasize - 1)/(tspan[2] - tspan[1]) + (1 - (datasize - 1)*tspan[1]/(tspan[2] - tspan[1]))))
    if scaled_t in [8, 9, 10]
        f(vcat(u..., exo_data[1, scaled_t], 3.5), p) # intervening on meals!
    else
        f(vcat(u, exo_data[:, scaled_t]...), p)
    end
  end

  thetas = []
  losses1 = []
  losses2 = []

  resol1s = []
  resol2s = []
  nsamples = 100
  dist = (@~ Normal(0, 0.05))(ω)

  for i = 1:nsamples
    noise = rand(dist, length(theta))
    noisy_theta = theta .+ noise
    resol1 = generate_series(network, u0, tspan, tsteps, noisy_theta; exo_data = exo_data)
    resol2 = generate_series(network, u0, tspan, tsteps, noisy_theta; intervene =  f_int, exo_data = new_exo_data)
    push!(resol1s, resol1)
    push!(resol2s, resol2)
  end
  resol1s, resol2s, new_exo_data, tsteps, exo_data, datatsteps
end

# ╔═╡ 43e1fed4-ad83-4d6a-ba37-ed9ff68eceeb
function genplot(resol1s, resol2s, tsteps, ode_data, nsamples, exo_data, new_exo_data, datatsteps)
  pl1 = plot(datatsteps, ode_data[1,:], seriestype = :scatter, color = :red, w=1.5, label = "CGM", xlabel = "t", title = "CGM, Steps, Bolus, Meals", legend = :top)
  plot!(pl1, datatsteps, ode_data[2,:], seriestype = :scatter, color = :blue, w=1.5, label = "Steps")
  plot!(pl1, datatsteps, ode_data[3,:], seriestype = :scatter, color = :green, w=1.5, label = "Bolus")
  plot!(pl1, datatsteps, ode_data[4,:], color = :purple, w=1.5, label = "Meals")

  pl2 = plot(datatsteps, ode_data[1,:], seriestype = :scatter, color = :red, w=1.5, label = "CGM", xlabel = "t", title = "CGM, Steps, Bolus, Meals", legend = :top)
  plot!(pl2, datatsteps, ode_data[2,:], seriestype = :scatter, color = :blue, w=1.5, label = "Steps")
  plot!(pl2, datatsteps, ode_data[3,:], seriestype = :scatter, color = :green, w=1.5, label = "Bolus")
  plot!(pl2, datatsteps, vcat(new_exo_data[2,1:7]..., 3.5, 3.5, 3.5, new_exo_data[2, 11:end]...), color = :purple, w=1.5, label = "Data: Meals")

  for i = 1:nsamples
    resol1 = resol1s[i]
    resol2 = resol2s[i]
    @show size(resol1)
    @show length(tsteps)
    @show typeof(resol1)
    plot!(pl1, tsteps,resol1[1,:], alpha=0.1, color = :red, label = "")
    plot!(pl1, tsteps,resol1[2,:], alpha=0.1, color = :blue, label = "")

    plot!(pl2, tsteps,resol2[1,:], alpha=0.1, color = :red, label = "")
    plot!(pl2, tsteps,resol2[2,:], alpha=0.1, color = :blue, label = "")
  end
  pl1, pl2
end

# ╔═╡ 850713d5-9103-4066-8168-dcb5851c259c
function plotexample()
  resol1s, resol2s, new_exo_data, tsteps, exo_data, datatsteps = examplep()
  nsamples = length(resol1s)
  # @show tsteps
  genplot(resol1s, resol2s, tsteps, ode_data, nsamples, exo_data, new_exo_data, datatsteps)
end

# ╔═╡ 50203b94-754c-430b-b999-a2606c135a48
p1, p2 = plotexample();

# ╔═╡ 0639a05b-222c-4926-9db2-db32fcd85e67
p1

# ╔═╡ 313b31cd-76a4-445f-8bb4-dd6d040546cc
p2

# ╔═╡ 9ff58d29-9eb3-4f57-bd97-0a6be1816f4c
network = gennet()

# ╔═╡ e8291111-4d8f-4484-a468-09d91b4fac79
function examplex(ω)
  exo_data = Data(ω)
  tspan = (0.0, 1.0)
  tsteps =  range(tspan[1], tspan[2], length = 500)
  datatsteps = range(tspan[1], tspan[2], length=size(ode_data)[2])
  datasize = size(ode_data)[2]
  u0 = ode_data[1:2, 1]

  function f_int(f, u, p, t)
    scaled_t = scaleidx(t, tspan[1], tspan[2], datasize)
    if scaled_t in [8, 9, 10]
        f(vcat(u..., exo_data[1, scaled_t], 3.5), p) # intervening on meals!
    else
        f(vcat(u, exo_data[:, scaled_t]...), p)
    end
  end

  intervention = I(ω)

  nsamples = 100
  noise = (@~ Normal(0, 0.05))(ω)
  noisy_theta = theta .+ noise
  generate_series(network, u0, tspan, tsteps, noisy_theta;
	  			  exo_data = exo_data, intervention = intervention)
end

# ╔═╡ 4ca22489-90b4-47cc-bb1b-4c66a2a221a4
randsample(examplex)

# ╔═╡ b453e716-c5c8-4bca-9131-26d8df68eb8f
exo_data = ode_data[3:4, :]

# ╔═╡ ba89f0ad-1fc3-48c3-b578-44472e2ae78b
function DataAddMeals(ω)
	new_exo_data = deepcopy(exo_data)
	new_exo_data[2, 8] = 3.5 # 13, 3.5
  	new_exo_data[2, 9] = 3.5 # 14, 3.5
  	new_exo_data[2, 10] = 3.5 # 15, 3.5
	new_exo_data
end

# ╔═╡ 4ecdc9aa-7c0b-4312-818d-134b2ed52a41
examples |ᵈ (Data => DataAddMeals)

# ╔═╡ 281efedb-aaea-4e48-8840-3b5b588dcae3
tspan = (0.0, 1.0)

# ╔═╡ 35560947-58f2-4a51-b926-22a74fdab58a
tsteps =  range(tspan[1], tspan[2], length = 500)

# ╔═╡ 9b128cdd-5e94-4639-a097-91d0f98e131d
datatsteps = range(tspan[1], tspan[2], length=size(ode_data)[2])

# ╔═╡ 6be8b340-781b-42ee-8e80-c78da347c5dc
datasize = size(ode_data)[2]

# ╔═╡ 39eb16c8-a93d-442c-8dc9-1982cd0ae79f
u0 = ode_data[1:2, 1]

# ╔═╡ 8a47a4d9-dac3-472c-9f7f-776e47bdb1a9
T = @~ DiscreteUniform(1, length(datatsteps))

# ╔═╡ 89c0ca23-5446-4093-89a7-a16e02c8cf55
randsample(T)

# ╔═╡ 8171ca8d-16f7-423c-be29-07b82d64dca9
glevel(data) = minimum(data[1,:] .* glucose_std .+ glucose_mean)

# ╔═╡ 141ad1e5-3abb-4db8-9d5d-9341839a28d4
function genmodel(ω)
  t_ = T(ω)

  # Generate new exogenous data
  new_exo_data = copy(exo_data)
  new_exo_data[2, 8] = 3.5
  new_exo_data[2, 9] = 3.5
  new_exo_data[2, 10] = 3.5

  function f_int(f, u, p, t)
    scaled_t = scaleidx(t, tspan[1], tspan[2], datasize)
      # scaled_t = Int(round(t*(datasize - 1)/(tspan[2] - tspan[1]) + (1 - (datasize - 1)*tspan[1]/(tspan[2] - tspan[1]))))
    if scaled_t in [t_-1, t_, t_+1]
        f(vcat(u..., exo_data[1, scaled_t], 3.5), p) # intervening on bolus!
    else
        f(vcat(u, exo_data[:, scaled_t]...), p)
    end
  end

  # dist = Normal(0, 0.05)
  # noise = rand(dist, length(theta))
  # noise = rand(dist, 1)
  noise = (@~ Normal(0.0, 0.05))(ω)
  noisy_theta = theta .+ noise
  resol1 = generate_series(network, u0, tspan, tsteps, noisy_theta; intervene = f_int, exo_data = exo_data)
  # @show glevel(resol1) > 100.0
  Omega.cond!(ω, glevel(resol1) > 100.0)
  t_
end

# ╔═╡ 6b01d502-73e2-4df7-9065-c35b1f5d9f2d
samplse = randsample(genmodel, 1000)

# ╔═╡ 2a20b0da-65d1-4d03-868e-e410e74effd4
histogram(samplse)

# ╔═╡ 2467b923-fbcf-4495-a002-3e2db55299c9
glevelrv(ω) = glevel(genmodel(ω))

# ╔═╡ 2b3dc1a7-5ea2-43df-8a64-c36fa2b52bed
gm2 = ~genmodel

# ╔═╡ 22fb08da-474a-4614-9e89-f45660ba60d2
function fscale(f, alpha = 4)
  Plots.scalefontsizes(alpha)
  try
    x = f()
    Plots.scalefontsizes(1/alpha)
    x
  catch e
    Plots.scalefontsizes(1/alpha)
    throw(e)
  end
end

# ╔═╡ Cell order:
# ╠═9b06e82a-83d9-11ec-2bec-45900ccb18c2
# ╠═dc3cc130-aa7d-417e-ac66-b3b38d4a0036
# ╠═0dd8222c-a031-405e-b1fb-523b4fae908f
# ╠═28566c39-a010-4f61-864d-4f052b1cb784
# ╠═ab98112d-6ddc-4b51-aaed-c46711b745e1
# ╠═746fdb88-587f-417b-bbf9-6dca06110730
# ╠═2cb2b71b-e327-49a4-852c-84768a379326
# ╠═95ff457e-a43e-4429-8289-7a86fbe44b72
# ╠═87467485-7e2f-454b-95f6-89996703890d
# ╠═6582679c-f481-4359-bb31-89eb23cd361b
# ╠═4ebd0515-b8c5-4069-83bd-a8bd17a4b90f
# ╠═99397fa7-0cf3-447a-b03e-711514b5a505
# ╠═5bfbf974-b7b5-4aac-ae8c-3c2b5d61eb81
# ╠═6bc373f4-1964-4191-a283-735277652912
# ╠═5ef4dd73-79a3-49c3-95b0-936865aaff58
# ╠═615ca691-a2c1-423b-8199-f66b1d07cd19
# ╠═ec7bdeff-b05d-4092-ac3e-bc334c1b5d4d
# ╠═8fbaec56-660d-4148-bd89-1b81f0f382a3
# ╠═ba89f0ad-1fc3-48c3-b578-44472e2ae78b
# ╠═e8291111-4d8f-4484-a468-09d91b4fac79
# ╠═4ecdc9aa-7c0b-4312-818d-134b2ed52a41
# ╠═4ca22489-90b4-47cc-bb1b-4c66a2a221a4
# ╠═79d0e2cc-fea3-42b6-b8bd-e4b36fee2bd9
# ╠═850713d5-9103-4066-8168-dcb5851c259c
# ╠═50203b94-754c-430b-b999-a2606c135a48
# ╠═0639a05b-222c-4926-9db2-db32fcd85e67
# ╠═313b31cd-76a4-445f-8bb4-dd6d040546cc
# ╠═43e1fed4-ad83-4d6a-ba37-ed9ff68eceeb
# ╠═9ff58d29-9eb3-4f57-bd97-0a6be1816f4c
# ╠═b453e716-c5c8-4bca-9131-26d8df68eb8f
# ╠═281efedb-aaea-4e48-8840-3b5b588dcae3
# ╠═35560947-58f2-4a51-b926-22a74fdab58a
# ╠═9b128cdd-5e94-4639-a097-91d0f98e131d
# ╠═6be8b340-781b-42ee-8e80-c78da347c5dc
# ╠═39eb16c8-a93d-442c-8dc9-1982cd0ae79f
# ╠═8a47a4d9-dac3-472c-9f7f-776e47bdb1a9
# ╠═141ad1e5-3abb-4db8-9d5d-9341839a28d4
# ╠═89c0ca23-5446-4093-89a7-a16e02c8cf55
# ╠═6b01d502-73e2-4df7-9065-c35b1f5d9f2d
# ╠═2a20b0da-65d1-4d03-868e-e410e74effd4
# ╠═8171ca8d-16f7-423c-be29-07b82d64dca9
# ╠═2467b923-fbcf-4495-a002-3e2db55299c9
# ╠═2b3dc1a7-5ea2-43df-8a64-c36fa2b52bed
# ╠═22fb08da-474a-4614-9e89-f45660ba60d2

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using NCDatasets
using Test

sz = (4,5)
filename = tempname()

# The mode "c" stands for creating a new file (clobber)
NCDatasets.NCDataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,Char,String]
        fv = fillvalue(T)

        data, scalar_data =
            if T == String
                [Char(i+60) * Char(j+60) for i = 1:sz[1], j = 1:sz[2]], "abcde"
            else
                [T(i+2*j+62) for i = 1:sz[1], j = 1:sz[2]], T(100)
            end

        # mask the frist element
        datam = Array{Union{T,Missing}}(data)
        datam[1] = missing

        #for define_type in [:with_kw,:with_attrib]
        for define_type in [:with_kw]

            kwargs =
                if define_type == :with_kw
                    Dict(:fillvalue => fv)
                else
                    Dict(:attrib => Dict("_FillValue" => fv))
                end

            v = NCDatasets.defVar(ds,"var-$T-$define_type",T,("lon","lat"); kwargs...)

            @test fillvalue(v) == fv
            @test v.attrib["_FillValue"] == fv
            @test fillmode(v) == (false,fv)

            v[:,:] = datam
            @test ismissing(v[1,1])
            @test isequal(v[:,:],datam)

            # load without transformation
            @test v.var[1,1] == fv

            # write/read without transformation
            v.var[:,:] = data
            @test v.var[:,:] == data
        end
    end
end

# NaN as fillvalue

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDatasets.NCDataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
ds.dim["lon"] = 10
ds.dim["lat"] = 11

fv = NaN32
v = NCDatasets.defVar(ds,"var_with_missing_data",Float32,("lon","lat"), fillvalue = fv)

data = [Float32(i+j) for i = 1:10, j = 1:11]
# mask the frist element
datam = Array{Union{Float32,Missing}}(data)
datam[1] = missing

v[:,:] = datam
@test ismissing(v[1,1])
@test isequal(v[:,:],datam)

# load without transformation
@test isnan(v.var[1,1])

NCDatasets.close(ds)

# all fill-values
filename = tempname()
ds = NCDatasets.NCDataset(filename,"c")
ds.dim["lon"] = 3

v = NCDatasets.defVar(ds,"var_with_all_missing_data",Float32,("lon",), fillvalue = fv)
data = [missing, missing, missing]

v[:] = data
@test all(ismissing.(v[:]))

v[1] = 1234.
@test !ismissing(v[1])

v[1] = missing
@test ismissing(v[1])

v[1] = 1234.
v[1:1] .= missing
@test ismissing(v[1])

NCDatasets.close(ds)

# test nomissing

data = [missing, Float64(1.), Float64(2.)]
@test_throws ErrorException NCDatasets.nomissing(data)

dataf = NCDatasets.nomissing(data,-9999.)
@test eltype(dataf) == Float64
@test dataf == [-9999., 1., 2.]

data = Union{Float64,Missing}[1., 2.]
dataf = NCDatasets.nomissing(data)
@test eltype(dataf) == Float64
@test dataf == [1., 2.]

@test nomissing(Union{Int64,Missing}[]) == []

@test nomissing([1,2,3,4],-9999) == [1,2,3,4]
@test nomissing([missing,2,3,4],-9999) == [-9999,2,3,4]

# issue 39
using NCDatasets

filename = tempname()

NCDataset(filename, "c") do ds
    defDim(ds, "lon", 2)
    defDim(ds, "lat", 2)
    v = defVar(ds, "Char variable", Char, ("lon","lat"), fillvalue = ' ')
    v[:,:] = ['a' 'b'; 'c' 'd']
end

NCDataset(filename, "r") do ds
    @test ds["Char variable"][:,:] == ['a' 'b'; 'c' 'd']
end

rm(filename)

# String fillvalues

filename = tempname()
ds = Dataset(filename,"c")
defDim(ds, "mode_items", 10)
v = defVar(ds,"instrument_mode",String,("mode_items",),
       attrib = Dict("_FillValue" => "UNDEFINED MODE"))

@test fillvalue(v) == "UNDEFINED MODE"
close(ds)
rm(filename)

filename = tempname()
ds = Dataset(filename,"c")
defDim(ds, "mode_items", 10)
v = defVar(ds,"instrument_mode",String,("mode_items",),fillvalue = "UNDEFINED MODE")
@test fillvalue(v) == "UNDEFINED MODE"
close(ds)
rm(filename)


text:
using Plots
using LaTeXStrings
using LinearAlgebra
using HDF5
using DelimitedFiles
using LsqFit
using Printf

pgfplotsx()
theme(
    :solarized_light;
    fg = RGB(88 / 255, 110 / 255, 117 / 255),
    fgtext = RGB(7 / 255, 54 / 255, 66 / 255),
    fgguide = RGB(7 / 255, 54 / 255, 66 / 255),
    fglegend = RGB(7 / 255, 54 / 255, 66 / 255),
)

function extract_time(file::HDF5.File, i::Int; group = "/χ")
    i = string(i, pad = 4)
    read(attributes(file["$group/$i"]), "time")
end
extract_time(file::AbstractString, i::Int; kwargs...) =
    h5open(io -> extract_time(io, i; kwargs...), file, "r")

function extract_time(sizes::AbstractVector{Int}; suffix)
    map(n -> extract_time("output_square_$(n)x$(n)_$suffix.h5", 2), sizes)
end

function save_timings(sizes)
    sizes = collect(sizes)
    cpu = extract_time(sizes; suffix = "cpu")
    cuda = extract_time(sizes; suffix = "cuda")
    writedlm("benchmark_result_square.dat", [sizes cpu cuda])
end

function determine_scaling(x, y)
    @. model(x, p) = p[1] + p[2] * x
    p0 = [-10.0, 4.0]
    fit = curve_fit(model, log.(x), log.(y), p0)

    f = n -> exp(fit.param[1]) * n .^ fit.param[2]
    return fit.param[2], f
end

function plot_times()
    table = readdlm("benchmark_result_square.dat")
    sizes, cpu, cuda = table[:, 1], table[:, 2], table[:, 3]
    sizes = sizes .^ 2 # benchmark.sh stores width of the sample rather
    cpu_exp, cpu_line = determine_scaling(sizes[5:end], cpu[5:end])
    cuda_exp, cuda_line = determine_scaling(sizes[5:end], cuda[5:end])
    continuous_size = collect(sizes[begin]:0.1:sizes[end])

    g = plot(xlabel = "System size", ylabel = "Time per frequency, seconds", legend = :topleft)
    plot!(
        g,
        continuous_size,
        cpu_line(continuous_size),
        linewidth = 2,
        linecolor = :black,
        linealpha = 0.6,
        label = nothing,
    )
    plot!(
        g,
        continuous_size,
        cuda_line(continuous_size),
        linewidth = 2,
        linecolor = :black,
        linealpha = 0.6,
        label = nothing,
    )
    scatter!(g, sizes, cpu, label = @sprintf("CPU (Xeon E5-2450 v2), ν = %.2f", cpu_exp))
    scatter!(g, sizes, cuda, label = @sprintf("GPU (Tesla K40), ν = %.2f", cuda_exp))
    savefig(g, "benchmark_result_square.png")
    g
end


text:
for i = 1:10
  print(i)
  i == 10 && break
  print(", ")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Dojo

# Parameters
origin = Origin{Float64}()
pbody = Box(1.0, 1.0, 1.0, 1.0)
cbody = Sphere(0.5, 1.0)
joint1 = JointConstraint(Floating(origin, pbody))
joint2 = JointConstraint(Floating(origin, cbody))

bodies = [pbody, cbody]
joints = [joint1, joint2]

side = 1.0
corners = [
            [[ side / 2.0;  side / 2.0; -side / 2.0]]
            [[ side / 2.0; -side / 2.0; -side / 2.0]]
            [[-side / 2.0;  side / 2.0; -side / 2.0]]
            [[-side / 2.0; -side / 2.0; -side / 2.0]]
            [[ side / 2.0;  side / 2.0;  side / 2.0]]
            [[ side / 2.0; -side / 2.0;  side / 2.0]]
            [[-side / 2.0;  side / 2.0;  side / 2.0]]
            [[-side / 2.0; -side / 2.0;  side / 2.0]]
        ]

normal = [[0.0, 0.0, 1.0] for i = 1:8]
contact_radius = [0.0 for i = 1:8]
friction_coefficient = 0.5 * ones(8)

box_contacts = contact_constraint(pbody, normal,
    friction_coefficient=friction_coefficient,
    contact_origins=corners,
    contact_radius=contact_radius,
    contact_type=:nonlinear)

sphere_contact = contact_constraint(cbody, [0.0; 0.0; 1.0],
    friction_coefficient=0.5,
    contact_origin=SA[0.0; 0.0; 0.0],
    contact_radius=0.5)

collision = SphereBoxCollision{Float64,2,3,6}(
    szeros(3),
    SA[0.5; 0.0; 0.0],
    SA[0.0; 0.5; 0.0],
    SA[0.0; 0.0; 0.5],
    0.5,
)

friction_parameterization = SA{Float64}[
    1.0  0.0
    0.0  1.0
]
body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

contacts = [
        # corner_contact1,
        # corner_contact2,
        # corner_contact3,
        # corner_contact4,
        # corner_contact5,
        # corner_contact6,
        # corner_contact7,
        # corner_contact8,
        box_contacts...,
        sphere_contact,
        ContactConstraint((body_body_contact, cbody.id, pbody.id), name=:body_body)
]

mech = Mechanism(origin, bodies, joints, contacts,
            gravity=1.0 * -9.81,
            timestep=0.01)

mech.bodies[1].state.x2 = [0.0, 0.0, 1.5]
q = rand(4)
q ./= norm(q)
mech.bodies[1].state.q2 = Quaternion(q..., true)
mech.bodies[2].state.x2 = [0.3, 0.25, 3.0]

storage = simulate!(mech, 5.0,
    verbose=true,
    record=true,
    opts=SolverOptions(verbose=true, rtol=1.0e-5, btol=1.0e-5))

vis = Visualizer()
visualize(mech, storage,
    vis=vis)
render(vis)


text:
using AutoOptimize
using Test

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


text:
abstract type AbstractCFRPolicy <: AbstractPolicy end

function Base.run(
    p::AbstractCFRPolicy,
    env::AbstractEnv,
    stop_condition = StopAfterStep(1),
    hook = EmptyHook(),
)
    @assert NumAgentStyle(env) isa MultiAgent
    @assert DynamicStyle(env) === SEQUENTIAL
    @assert RewardStyle(env) === TERMINAL_REWARD
    @assert ChanceStyle(env) === EXPLICIT_STOCHASTIC
    @assert DefaultStateStyle(env) isa InformationSet

    RLBase.reset!(env)

    while true
        update!(p, env)
        hook(POST_ACT_STAGE, p, env)
        stop_condition(p, env) && break
    end
    update!(p)
end


text:
Tables.istable(::Type{<:Result}) = true
Tables.rowaccess(::Type{<:Result}) = true
Tables.rows(jl_result::Result) = jl_result

Base.eltype(jl_result::Result) = Row
Base.length(jl_result::Result) = num_rows(jl_result)

function Tables.schema(jl_result::Result)
    types = map(jl_result.not_null, column_types(jl_result)) do not_null, col_type
        not_null ? col_type : Union{col_type, Missing}
    end
    return Tables.Schema(map(Symbol, column_names(jl_result)), types)
end

function Base.iterate(jl_result::Result, (len, row)=(length(jl_result), 1))
    row > len && return nothing
    return Row(jl_result, row), (len, row + 1)
end

struct Row
    result::Result
    row::Int
end

Base.propertynames(r::Row) = column_names(getfield(r, :result))

function Base.getproperty(pqrow::Row, name::Symbol)
    jl_result = getfield(pqrow, :result)
    row = getfield(pqrow, :row)
    col = column_number(jl_result, name)
    if libpq_c.PQgetisnull(jl_result.result, row - 1, col - 1) == 1
        return missing
    else
        oid = jl_result.column_oids[col]
        T = jl_result.column_types[col]
        return jl_result.column_funcs[col](PQValue{oid}(jl_result, row, col))::T
    end
end

# sink
function load!(table::T, connection::Connection, query::AbstractString) where {T}
    Tables.istable(T) || throw(ArgumentError("$T doesn't support the required Tables.jl interface"))
    stmt = prepare(connection, query)
    rows = Tables.rows(table)
    state = iterate(rows)
    state === nothing && return
    st, row = state
    names = propertynames(row)
    sch = Tables.Schema(names, nothing)
    parameters = Vector{Parameter}(undef, length(names))
    while true
        Tables.eachcolumn(sch, row) do val, col, nm
            parameters[col] = if ismissing(val)
                missing
            elseif val isa AbstractString
                convert(String, val)
            else
                string(val)
            end
        end
        close(execute(stmt, parameters; throw_error=true))
        state = iterate(rows, st)
        state === nothing && break
        row, st = state
    end
    return
end


text:
Base.promote_rule(::Type{Decimal}, ::Type{<:Real}) = Decimal
Base.promote_rule(::Type{BigFloat}, ::Type{Decimal}) = Decimal
Base.promote_rule(::Type{BigInt}, ::Type{Decimal}) = Decimal

const RoundUnnecessary = RoundingMode{:RoundUnnecessary}()
const RoundFloor       = RoundingMode{:RoundFloor}()
const RoundCeiling     = RoundingMode{:Ceiling}()

function divide(x::BigInt, y::Union{Integer, BigInt}, round::RoundingMode)
    c, m = divrem(x, y)
    if round === RoundUnnecessary
        m == 0 || throw(DivideError())
        return c
    elseif round === RoundFloor
        c > 0 && return c
        return m == 0 ? c : c - 1
    elseif round === RoundCeiling
        c < 0 && return c
        return m == 0 ? c : c + 1
    elseif round === RoundDown
        return c
    elseif round === RoundUp
        m == 0 && return c
        return c >= 0 ? c + 1 : c - 1
    end
end

function setexponent(x::Decimal, q::Integer, round::RoundingMode=RoundUnnecessary)
    # c*10^q = c*(10^p*10^(q-p)) = (c*10^(q-p))*10^p
    p = x.q - q

    p == 0 && return x

    if p > 0
        c = x.c * big"10"^p
    else
        c = divide(x.c, big"10"^-p, round)
    end
    Decimal(c, q)
end

function strip_trailing_zeros(x::Decimal)

    iszero(x.c) && return Zero

    c, q = x.c, x.q
    while true
        d, m = divrem(c, big"10")
        if !iszero(m)
            break
        end
        c = d
        q += 1
    end
    Decimal(c, q)
end

function +(x::Decimal, y::Decimal)
    q = min(x.q, y.q)
    x = setexponent(x, q)
    y = setexponent(y, q)
    Decimal(x.c + y.c, q)
end

-(x::Decimal) = Decimal(-x.c, x.q)
-(x::Decimal, y::Decimal) = x + (-y)

function mul(x::Decimal, y::Decimal, precision::Integer=PRECISION, round::RoundingMode=ROUNDING)
    p = Decimal(x.c * y.c, x.q + y.q)
    p = p.q >= -precision ? p : strip_trailing_zeros(p)
    p.q >= -precision ? p : setexponent(p, -precision, round)
end

*(x::Decimal, y::Decimal) = mul(x, y)

function div(x::Decimal, y::Decimal, precision::Integer=PRECISION, round::RoundingMode=ROUNDING)
    if x.q > (p = y.q - precision)
        x = setexponent(x, p)
    else
        y = setexponent(y, x.q + precision)
    end
    c = divide(x.c, y.c, round)
    Decimal(c, -precision)
end

/(x::Decimal, y::Decimal) = div(x, y)

const Zero = Decimal(0, 0)
const One  = Decimal(1, 0)
zero(::Type{Decimal}) = Zero
one(::Type{Decimal}) = One
eps(x::Decimal) = Decimal(1, x.q)

function inv(x::Decimal, precision::Integer=PRECISION, round::RoundingMode=ROUNDING)
    div(One, x, precision, round)
end

for comp in (:(==), :<, :>, :<=, :>=)
    @eval function $comp(x::Decimal, y::Decimal)
        q = min(x.q, y.q)
        x = setexponent(x, q)
        y = setexponent(y, q)
        $comp(x.c, y.c)
    end
end

number(x::Decimal) = x.c * 10.0^x.q

signbit(x::Decimal) = signbit(x.c)
sign(x::Decimal) = sign(x.c)


text:
@testset "1012.numbers-with-repeated-digits.jl" begin
    @test num_dup_digits_at_most_n(3562) == 1561
    @test num_dup_digits_at_most_n(20) == 1
    @test num_dup_digits_at_most_n(8315) == 3870
    @test num_dup_digits_at_most_n(1000) == 262
end


text:
module ChargedParticleDynamics

    include("initialization/initial_conditions.jl")

    include("ChargedParticle3d.jl")
    include("GuidingCenter4d.jl")
    include("GyroKinetics4d.jl")
    include("PauliParticle3d.jl")
    include("ChargedParticlePlots.jl")

end