text:
module HEP

import Base.convert
import Base.eta

import Base.+
import Base.-
import Base.*

#A 4-tuple representing a 4 vector in the natural system of units
abstract AbstractFourVector

#Cartesian
immutable FourVector <: AbstractFourVector
    t::Float64
    x::Float64
    y::Float64
    z::Float64
end

#Spherical
immutable FourVectorSph <: AbstractFourVector
    perp::Float64 #perpendicular component (e.g. in case of 4-momentum -> pt or transverse momentum)
    eta::Float64
    phi::Float64
    l::Float64 #length
end

+(v1::FourVector, v2::FourVector) = FourVector(v1.t+v2.t, v1.x+v2.x, v1.y+v2.y, v1.z+v2.z)
*(v1::FourVector, a::FourVector) = FourVector(a * v1.t, a * v1.x, a * v1.y, a*v1.z)
-(v1::FourVector, v2::FourVector) = v1 + -1.0 * v2

#magnitude
l2(v::FourVector) = v.t^2 - v.x^2 - v.y^2 - v.z^2
function l(v::FourVector)
    sq = l2(v)
    return sq > 0 ? sqrt(sq) : -sqrt(-sq)
end

#polar
phi(v::FourVector) = atan2(v.y, v.x)
rho(v::FourVector) = sqrt(v.x^2 + v.y^2 + v.z^2)
theta(v::FourVector) = atan2(perp(v), v.z)

#pseudorapidity
function eta(v::FourVector)
    ct = cos(theta(v))
    return ct^2 < 1 ? -0.5 * log((1.0-ct)/(1.0+ct)) : sign(v.z)*Inf
end

#transverse component (perpendicular)
perp(v::FourVector) = sqrt(v.x^2 + v.y^2)

phi(v::FourVectorSph) = v.phi
eta(v::FourVectorSph) = v.eta
perp(v::FourVectorSph) = v.perp
l(v::FourVectorSph) = v.l
rho(v::FourVectorSph) = rho(convert(v, FourVector))
theta(v::FourVectorSph) = theta(convert(v, FourVector))

x(v::FourVector) = v.x
y(v::FourVector) = v.y
z(v::FourVector) = v.z
t(v::FourVector) = v.t

x(v::FourVectorSph) = perp(v) * cos(phi(v))
y(v::FourVectorSph) = perp(v) * sin(phi(v))
z(v::FourVectorSph) = perp(v) * sinh(eta(v))
t(v::FourVectorSph) = sqrt(sign(l(v)) * l(v)^2 + perp(v)^2 * (1 + sinh(eta(v))^2))

# #spherical components
# function FourVectorSph(perp, eta, phi, l)
#     perp = abs(perp)
#     x = perp * cos(phi)
#     y = perp * sin(phi)
#     z = perp * sinh(eta)
#     t = l>0 ? sqrt(l^2+x^2+y^2+z^2) : sqrt(-l^2+x^2+y^2+z^2)
#     return FourVector(t, x, y, z)
# end

function convert(a::FourVectorSph, b::Type{FourVector})
    FourVector(t(a), x(a), y(a), z(a))
end

function convert(a::FourVector, b::Type{FourVectorSph})
    FourVector(perp(a), eta(a), phi(a), l(a))
end

#calculates the delta R between 2 four momenta
function deltar{T1 <: AbstractFourVector, T2 <: AbstractFourVector}(v1::T1, v2::T2)
    deta = eta(v1) - eta(v2)
    dphi = phi(v1) - phi(v2)

    #bring dphi to -pi...+pi
    while dphi >= pi
        dphi -= 2 * pi
    end

    while dphi < -pi
        dphi += 2 * pi
    end

    return sqrt(deta^2 + dphi^2)
end

export FourVector, FourVectorSph
export l2, l, phi, rho, theta, eta, perp
export x, y, z, t
export deltar
export +, -, *

#
# function to_jl(o::TH1D; error_type=:errors)
#
#     nb = Int64(GetNbinsX(o))
#     #+3 = underflow, overflow, superfluous extra bin
#     conts = zeros(Float64, nb+2)
#     errs = zeros(Float64, nb+2)
#     ents = zeros(Float64, nb+2)
#
#     #underflow low, overflow low, overflow high
#     edges = zeros(Float64, nb+3)
#     for n=0:nb+1
#         conts[n+1] = GetBinContent(o, Int32(n)) |> Float64
#         errs[n+1] = GetBinError(o, Int32(n)) |> Float64
#         #entries[n+1] = GetEntries(h) * conts[n+1]
#         edges[n+1] = GetBinLowEdge(o, Int32(n))
#     end
#
#     if error_type == :errors
#         ents = (conts ./ errs).^2
#         ents[isnan(ents)] = 0.0
#         ents[ents .== Inf] = 1.0
#         #println(hcat(conts, errs, ents, conts./sqrt(ents)))
#     end
#
#     #this works for Poisson bin errors
#     if error_type == :contents
#         ents = conts ./ sum(conts) .* Float64(GetEntries(o))
#         ents[isnan(ents)] = 0.0
#     end
#
#     edges[1] = -Inf
#     edges[nb+2] = edges[nb+1] + GetBinWidth(o, Int32(nb))
#     edges[nb+3] = Inf
#
#     return (edges,), conts, errs
#
# end
#
# function to_jl(o::TH2D)
#
#     nx = Int64(GetNbinsX(o))
#     ny = Int64(GetNbinsY(o))
#
#     arr = zeros(Float64, nx+2, ny+2)
#     errs = zeros(Float64, nx+2, ny+2)
#
#     edges_x = zeros(Float64, nx + 3)
#     edges_y = zeros(Float64, ny + 3)
#
#     for n=0:nx+1
#         #edges_x[n+1] = GetBinLowEdge(o, int32(n))
#         edges_x[n+1] = ROOT.GetBinLowEdgeX(o, Int32(n))
#     end
#     for n=0:ny+1
#         #edges_x[n+1] = GetBinLowEdge(o, int32(n))
#         edges_y[n+1] = ROOT.GetBinLowEdgeY(o, Int32(n))
#     end
#     for x=0:nx+1
#         for y=0:ny+1
#             arr[x+1, y+1] = GetBinContent(root_cast(TH1, o), Int32(x), Int32(y))
#             errs[x+1, y+1] = GetBinError(root_cast(TH1, o), Int32(x), Int32(y))
#         end
#     end
#     edges_x[1] = -Inf
#     edges_x[end] = +Inf
#
#     edges_y[1] = -Inf
#     edges_y[end] = +Inf
#
#     ents = arr .* GetEntries(o) ./ sum(arr)
#     ents[isnan(ents)] = 0
#     ents = round(ents)
#     return (edges_x, edges_y), arr, errs
# end
#
# function to_jl(o::TH3D)
#
#     nx = Int64(GetNbinsX(o))
#     ny = Int64(GetNbinsY(o))
#     nz = Int64(GetNbinsZ(o))
#
#     arr = zeros(Float64, nx+2, ny+2, nz+2)
#     errs = zeros(Float64, nx+2, ny+2, nz+2)
#
#     edges_x = zeros(Float64, nx + 3)
#     edges_y = zeros(Float64, ny + 3)
#     edges_z = zeros(Float64, nz + 3)
#
#     for n=0:nx+1
#         #edges_x[n+1] = GetBinLowEdge(o, int32(n))
#         edges_x[n+1] = ROOT.GetBinLowEdgeX(o, Int32(n))
#     end
#     for n=0:ny+1
#         #edges_x[n+1] = GetBinLowEdge(o, int32(n))
#         edges_y[n+1] = ROOT.GetBinLowEdgeY(o, Int32(n))
#     end
#     for n=0:nz+1
#         #edges_x[n+1] = GetBinLowEdge(o, int32(n))
#         edges_z[n+1] = ROOT.GetBinLowEdgeZ(o, Int32(n))
#     end
#     for x=0:nx+1
#         for y=0:ny+1
#             for z=0:nz+1
#                 arr[x+1, y+1, z+1] = GetBinContent(root_cast(TH1, o), Int32(x), Int32(y), Int32(z))
#                 errs[x+1, y+1, z+1] = GetBinError(root_cast(TH1, o), Int32(x), Int32(y), Int32(z))
#             end
#         end
#     end
#     edges_x[1] = -Inf
#     edges_x[end] = +Inf
#
#     edges_y[1] = -Inf
#     edges_y[end] = +Inf
#
#     edges_z[1] = -Inf
#     edges_z[end] = +Inf
#
#     ents = arr .* GetEntries(o) ./ sum(arr)
#     ents[isnan(ents)] = 0
#     ents = round(ents)
#     return (edges_x, edges_y, edges_z), arr, errs
# end
#
# function read_file(fn)
#     tf = TFile(fn)
#     kl = GetListOfKeys(tf)
#
#     hd = Dict()
#     for i=1:length(kl)
#         k = root_cast(TKey, kl[i])
#         kn = k|>GetName|>bytestring
#         o = ReadObj(k);
#         clname = k|>GetClassName|>bytestring|>symbol
#
#         o = root_cast(eval(clname), o)
#         jo = to_jl(o)
#
#         #println(jo)
#
#         #h1 = Histogram(jo[1], jo[2])
#         #h2 = Histogram(jo[1], jo[3])
#
#     #    println(h2)
#         hd[symbol(kn)] = ErrorHistogram{Float64, length(jo[1]), typeof(jo[1])}(
#             jo[1], jo[2], jo[3].^2, :right
#         )
#     #    println(jo)
#         #println(kn, " ", clname, " ", o,)
#     end
#     Close(tf)
#     return hd
#end

quadsum(arr::AbstractArray) = arr.^2 |> sum |> sqrt
quadsum(arr::AbstractArray, n::Int64) = sum(arr.^2, n) |> sqrt
#
# import Base.sum
# function sum{T<:Real, N, E}(h::ErrorHistogram{T,N,E}, n::Int64)
#     inds = [1:N]
#     splice!(inds, n)
#     newedges = deepcopy(h.edges[inds])
#     dims = map(x->size(x)[1]-1, newedges)
#     #println(newedges, " ", dims)
#     return ErrorHistogram{T, N-1, typeof(newedges)}(
#         newedges, reshape(sum(h.weights, n), dims...), reshape(sum(sqrt(h.weights_sq), n).^2, dims...), :right
#     )
# end
#
# import Base.getindex
# function getindex{T<:Real, N, E}(h::ErrorHistogram{T, N, E}, inds...)
#     newinds = Any[]
#     for i in inds
#         push!(newinds, i.start:i.stop+1)
#     end
#     edgs = tuple([h.edges[i][newinds[i]] for i=1:length(newinds)]...)
#     return ErrorHistogram{T, N, E}(edgs, h.weights[inds...], h.weights_sq[inds...], :right)
# end

#export to_jl, read_file

end # module


text:
"""
This code is based on the original work on the Scilab DoE package by the
following individuals
    Copyright (C) 2012 - 2013 - Michael Baudin
    Copyright (C) 2012 - Maria Christopoulou
    Copyright (C) 2010 - 2011 - INRIA - Michael Baudin
    Copyright (C) 2009 - Yann Collette
    Copyright (C) 2009 - CEA - Jean-Marc Martinez
and
    Copyright (C) 2014, Abraham D. Lee
who converted the package to Python.
"""

"""
$(TYPEDSIGNATURES)

Constructs a central composite design.

```jldoctest
julia> ccdesign(3)
22×3 Matrix{Float64}:
 -1.0      -1.0      -1.0
  1.0      -1.0      -1.0
 -1.0       1.0      -1.0
  1.0       1.0      -1.0
 -1.0      -1.0       1.0
  1.0      -1.0       1.0
 -1.0       1.0       1.0
  1.0       1.0       1.0
  0.0       0.0       0.0
  0.0       0.0       0.0
  ⋮
  1.82574   0.0       0.0
  0.0      -1.82574   0.0
  0.0       1.82574   0.0
  0.0       0.0      -1.82574
  0.0       0.0       1.82574
  0.0       0.0       0.0
  0.0       0.0       0.0
  0.0       0.0       0.0
  0.0       0.0       0.0

```
"""
function ccdesign(n::Int, center::Array{Int}=[4, 4], alpha::Symbol=:orthogonal, face::Symbol=:circumscribed)
    @assert n>1
    @assert alpha in [:orthogonal, :rotatable]
    @assert face in [:circumscribed, :inscribed, :faced]
    @assert length(center) == 2

    H2, a = star(n, alpha, center)
    if face == :inscribed
        H1 = explicit_fullfactorial(fullfactorial(Tuple(fill([-1, 1], n))))
        H1 = H1/a
        H2, a = star(n, :faced, [1,1])
    elseif face == :faced
        H2, a = star(n, :faced, [1,1])
        H1 = explicit_fullfactorial(fullfactorial(Tuple(fill([-1, 1], n))))
    elseif face == :circumscribed
        H1 = explicit_fullfactorial(fullfactorial(Tuple(fill([-1, 1], n))))
    else
        print("Invalid value of face")
    end

    C1 = zeros(center[1], n)
    C2 = zeros(center[2], n)

    H1 = vcat(H1, C1)
    H2 = vcat(H2, C2)
    H = vcat(H1, H2)
end

function star(n, alpha, center)
    if alpha  == :faced
        a = 1
    elseif alpha == :orthogonal
        nc = 2.0^n
        nco = center[1]
        na = 2.0*n
        nao = center[2]
        a = (n*(1 + nao/na)/(1+nco/nc))^0.5
    elseif alpha == :rotatable
        nc = 2.0^n
        a = nc^0.25
    else
        print("Invalid value of alpha")
    end

    H = zeros(2*n, n)
    for i in 1:n
         H[2*i-1:2*i, i] = [-1, 1]
    end
    H *= a
    return H, a
end

ext: jl


text:
module AISIS10016

using CSV, DataFrames

export table_g53, f322, app2C2262

function e2(Fcre, Fy, Ag, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.80
    elseif ASDorLRFD==1
        StrengthFactor=0.85
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    λc = sqrt(Fy/Fcre)

    if λc <= 1.5
        Fn = 0.658^((λc)^2) * Fy
    else
        Fn = (0.877 /(λc)^2) * Fy
    end

    Pne = Ag * Fn

    ePne = Pne * StrengthFactor

    return Pn, ePne

end

function f321(Mne, Mcrℓ, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.67
    elseif ASDorLRFD==1
        StrengthFactor=0.90
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    λℓ=sqrt(Mne/Mcrℓ)

    if λℓ <= 0.776
        Mnℓ=Mne
    else
        Mnℓ=(1-0.15*(Mcrℓ/Mne)^0.4)*(Mcrℓ/Mne)^0.4*Mne
    end

    eMnℓ = Mnℓ * StrengthFactor

    return Mnℓ, eMnℓ

end

function f322(Mne, Mcrℓ, My_net, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.67
    elseif ASDorLRFD==1
        StrengthFactor=0.90
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    #Mcrℓ includes the influence of holes here.
    Mnℓ, eMnℓ = f321(Mne, Mcrℓ, ASDorLRFD)

    if Mnℓ > My_net

        Mnℓ = My_net
        eMnℓ = Mnℓ * StrengthFactor

    end

    return Mnℓ, eMnℓ

end

function f411(My, Mcrd, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.67
    elseif ASDorLRFD==1
        StrengthFactor=0.90
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    λd=sqrt(My/Mcrd)

    if λd <= 0.673
        Mnd=My
    else
        Mnd=(1-0.22*(Mcrd/My)^0.4)*(Mcrd/My)^0.4*My
    end

    eMnd=Mnd*StrengthFactor

    return Mnd, eMnd

end

function g21(E, h, t, Fy, Vcr, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.60
    elseif ASDorLRFD==1
        StrengthFactor=0.95
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    Aw = h.*t

    Vy=0.6 .*Aw .*Fy
    λv=sqrt.(Vy./Vcr)

    if λv <= 0.815
        Vn = Vy
    elseif (λv>0.815) & (λv<=1.227)
        Vn = 0.815 .*sqrt.(Vcr.*Vy)
    elseif λv > 1.227
        Vn = Vcr
    end

    eVn=Vn*StrengthFactor

    return Vn, eVn

end

function g231(h, t, Fcr)

    Aw = h.*t
    Vcr = Aw.*Fcr

    return Vcr

end

function g232(E, μ, kv, h, t)

    Fcr = (π^2 .*E.*kv)./(12 .*(1-μ.^2).*(h./t).^2)

end

function g233(a, h)

    if a./h <= 1.0
        kv = 4.00 .+ 5.34./(a./h).^2
    elseif a./h > 1.0
        kv = 5.34 .+ 4.00./(a./h).^2
    end

    return kv

end

function g51(t, h, Fy, θ, C, C_R, R, C_N, N, C_h, ϕ_w, Ω_w, ASDorLRFD)

    Pn = C * t^2 * Fy * sin(deg2rad(θ)) * (1-C_R*sqrt(R/t)) * (1+C_N*sqrt(N/t)) * (1-C_h*sqrt(h/t))

    if ASDorLRFD == 0
        ePn = Pn / Ω_w
    elseif ASDorLRFD == 1
        ePn = ϕ_w * Pn
    else   #nominal
        ePn = Pn
    end

    return Pn, ePn

end

function table_g53()

    filename = string(@__DIR__, "/assets/AISI_S100_16_Table_G5_3.csv")
    data = CSV.File(filename)

    table = DataFrame(support_condition = Vector(data.support_condition), flange_condition = Vector(data.flange_condition), load_case = Vector(data.load_case), load_location=Vector(data.load_location), C=Vector(data.C), C_R = Vector(data.C_R), C_N = Vector(data.C_N), C_h = Vector(data.C_h), ASD = Vector(data.ASD), LRFD = Vector(data.LRFD), LSD = Vector(data.LSD), limits = Vector(data.limits))

    return table

end

function h21(Mbar, Vbar, Maℓo, Va)

    Interaction=sqrt.((Mbar./Maℓo).^2 .+ (Vbar./Va).^2)

end

function h121(Pbar, Mxbar, Mybar, Pa, Max, May)

    ActionP = abs.(Pbar./Pa)
    ActionMx = abs.(Mxbar./Max)
    ActionMy = abs.(Mybar./May)

    Interaction = ActionP .+ActionMx .+ ActionMy

    return ActionP, ActionMx, ActionMy, Interaction

end

function h31(Pbar, Mbar, Pn, Mnℓo, design_code)

    actionP = 0.91 * (Pbar/Pn)
    actionM = Mbar / Mnℓo
    interaction = actionP + actionM

    Ω = 1.70
    ϕ_LRFD = 0.95
    ϕ_LSD = 0.75

    if design_code == "AISI S100-16 ASD"
        DC = interaction / (1.33/Ω)
    elseif design_code == "AISI S100-16 LRFD"
        DC = interaction / (1.33*ϕ_LRFD)
    elseif design_code == "AISI S100-16 LSD"
        DC = interaction / (1.33*ϕ_LSD)
    end

    return actionP, actionM, interaction, DC

end

function h33(Pbar, Mbar, Pn, Mnℓo, design_code)

    actionP = 0.86 * (Pbar/Pn)
    actionM = Mbar / Mnℓo
    interaction = actionP + actionM

    Ω = 1.70
    ϕ_LRFD = 0.90
    ϕ_LSD = 0.80

    if design_code == "AISI S100-16 ASD"
        DC = interaction / (1.65/Ω)
    elseif design_code == "AISI S100-16 LRFD"
        DC = interaction / (1.65*ϕ_LRFD)
    elseif design_code == "AISI S100-16 LSD"
        DC = interaction / (1.65*ϕ_LSD)
    end

    return actionP, actionM, interaction, DC

end

function table23131(CorZ,t,b,d,θ)

    CorZ = convert(Int8, CorZ)

    Af=(b+d)*t
    Jf=1/3*b*t^3+1/3*d*t^3
    Cwf=0.0

    θ = deg2rad(θ)  #convert degrees to radians

    if CorZ==0

        Ixf=t*(t^2*b^2+4*b*d^3+t^2*b*d+t^2*b*d+d^4)/(12*(b+d))
        Iyf=t*(b^4+4*d*b^3)/(12*(b+d))
        Ixyf=(t*b^2*d^2)/(4*(b+d))
        xof=b^2/(2*(b+d))
        hxf=-(b^2+2*d*b)/(2*(b+d))
        hyf=-d^2/(2*(b+d))
        yof=hyf

    else

        Ixf=t*(t^2*b^2+4*b*d^3-4*b*d^3*cos(θ)^2+t^2*b*d+d^4-d^4*cos(θ)^2)/(12*(b+d))
        Iyf=t*(b^4+4*d*b^3+6*d^2*b^2*cos(θ)+4*d^3*b*cos(θ)^2+d^4*cos(θ)^2)/(12*(b+d))
        Ixyf=t*b*d^2*sin(θ)*(b+d*cos(θ))/(4*(b+d))
        xof=(b^2-d^2*cos(θ))/(2*(b+d))
        hxf=-(b^2+2*d*b+d^2*cos(θ))/(2*(b+d))
        hyf=-d^2*sin(θ)/(2*(b+d))
        yof=hyf
    end

    return Af,Jf,Ixf,Iyf,Ixyf,Cwf,xof,hxf,hyf,yof

end

function app23331(CorZ, t, ho, b, d, θc, E, μ, G, f1, f2, M1, M2, CurvatureSign, Lm, kϕ, Sf)

    θc=deg2rad(θc)
    bc=b-t/2-t/2*tan(θc/2)
    dc=d-t/2*tan(θc/2)

    Af,Jf,Ixf,Iyf,Ixyf,Cwf,xof,hxf,hyf,yof=table23131(CorZ,t,bc,dc,θc)

    L = app23334(ho, μ, t, Ixf, xof, hxf, Cwf, Ixyf, Iyf, Lm)

    β=app23333(L, Lm, M1, M2, CurvatureSign)

    kϕfe=app23133(E,Ixf,xof,hxf,Cwf,Ixyf,Iyf,L,G,Jf)

    kϕwe=app23335(E,t,μ,ho,L)

    kϕfg=app23135(L,Af,xof,hxf,Ixyf,Iyf,yof,Ixf)

    kϕwg=app23336(f1,f2,ho,t,L)

    Fcrd=app23332(β, kϕfe, kϕwe, kϕ, kϕfg, kϕwg)

    Mcrd=Sf*Fcrd

    return Mcrd

end

function app23333(L, Lm, M1, M2)

    β=1+0.4(L/Lm)^0.7*(1+M1/M2)^0.7

    if β>=1.3
        β=1.3
    end

    if β<=1.0
        β=1.0
    end

    return β

end

function app23334(ho, μ, t, Ixf, xof, hxf, Cwf, Ixyf, Iyf, Lm)

    Lcrd=(4*π^4*ho*(1-μ^2)/t^3*(Ixf*(xof-hxf)^2+Cwf-Ixyf^2/Iyf*(xof-hxf)^2)+π^4*ho^4/720)^(1/4)

    L=minimum([Lcrd, Lm])

    return Lcrd, L

end

function app23133(E,Ixf,xof,hxf,Cwf,Ixyf,Iyf,L,G,Jf)

    kϕfe=(π/L)^4*(E*Ixf*(xof-hxf)^2+E*Cwf-E*Ixyf^2/Iyf*(xof-hxf)^2)+(π/L)^2*G*Jf

end

function app23335(E,t,μ,ho,L)

    kϕwe=E*t^3/(12*(1-μ^2))*(3/ho+(π/L)^2*19/60*ho+(π/L)^4*ho^3/240)

end

function app23135(L,Af,xof,hxf,Ixyf,Iyf,yof,Ixf)

    kϕfg=(π/L)^2*(Af*((xof-hxf)^2*(Ixyf/Iyf)^2-2*yof*(xof-hxf)*(Ixyf/Iyf)+hxf^2+yof^2)+Ixf+Iyf)

end

function app23336(f1,f2,ho,t,L)

    ξweb=(f1-f2)/f1

    kϕwg=ho*t*π^2/13440*((((45360*(1-ξweb)+62160)*(L/ho)^2)+448*π^2+(ho/L)^2*(53+3*(1-ξweb))*π^4)/(π^4+28*π^2*(L/ho)^2+420*(L/ho)^4))

    return kϕwg

end

function app23332(β, kϕfe, kϕwe, kϕ, kϕfg, kϕwg)

    Fcrd=β*(kϕfe+kϕwe+kϕ)/(kϕfg+kϕwg)

end

function app2C2262(t, Lh, Lcrd)

    tr = t * (1 - Lh/Lcrd)^(1/3)

end

end #module


text:
using DataStructures;

#Heuristic misplaced

function misplaced(goal, state)
    x, y = size(state)
    z = 0
    for i = 1:x
        for j = 1:y
            if goal[i, j] == 0
                continue
            elseif goal[i, j] == state[i, j]
                continue
            else
                z += 1
            end
        end
    end
    return z
end

#Heuristic manhattan

function manhattan(goal, state)
    x, y = size(state)
    z = 0
    for i = 1:x
        for j = 1:y
            if goal[i, j] == 0
                continue
            end
            cord = findfirst(isequal(goal[i, j]), state)
            temp = abs(cord[1] - i)
            temp1 = abs(cord[2] - j)
            z += temp + temp1
        end
    end
    return z
end

#  Perform A* search with heuristic.
function astar(goal, heuristic, start)
    reset_timer!(to::TimerOutput)

    node = createtree(start)
    frontier = PriorityQueue{TreeNode,Integer}()
    frontier[node] = 0
    explored = []

    while !isempty(frontier)
        @timeit to "Moves" acts = moves(node.state)
        @timeit to "Dequeue" node = dequeue!(frontier)
        @timeit to "Push" push!(explored, node.state)
        @timeit to "Expand" states = expansion(node.state)

        for (a, s) in zip(acts, states)
            if in(explored, s)
                continue
            end
            @timeit to "Node Creation" cnode = addnode(a, node, s)
            @timeit to "Heuristic Score" score =
                node.pathcost + heuristic(goal, cnode.state)
            if s == goal
                displaysolution(cnode)
                return cnode
            end
            if !(cnode in keys(frontier))
                frontier[cnode] = score
            elseif frontier[cnode] > score
                frontier[cnode] = score
            end
        end
    end
    println("No solution found!"^10)
    return nothing
end

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
function counting_sort(A, k)
    B = Array(zeros(Int, (length(A),))')
    C = zeros(Int, (k+1,))

    for j in 1:length(A)
        C[A[j]+1] += 1
    end # C[i+1] now contains number of elems equal to i

    for i in 1:k
        C[i+1] += C[i]
    end # C[i+1] now contains number of elems ≤ i

    for j in range(length(A), 1, step = -1)
        B[C[A[j]+1]] = A[j]
        C[A[j]+1] -= 1
    end # B is sorted version of A

    return B
end

function radix_sort!(A)
    max = maximum(A)

    place = 1
    while max // place > 0
        counting_sort!(A, place)
        place *= 10
    end
end

function counting_sort!(A, place)
    B = Array(zeros(Int, (length(A),))')
    C = zeros(Int, (10,))

    for j in 1:length(A)
        index = Int(floor(A[j] / place))
        digit = (index % 10) + 1
        C[digit] += 1
    end # C[i+1] now contains number of elems equal to i

    for i in 1:9
        C[i+1] += C[i]
    end # C[i+1] now contains number of elems ≤ i

    for j in range(length(A), 1, step = -1)
        index = Int(floor(A[j] / place))
        digit = (index % 10) + 1
        B[C[digit]] = A[j]
        C[digit] -= 1
    end # B is sorted version of A

    # copy
    A .= B
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
##########################################################
#FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
#Test problem is chosen to give an exact solution at all nodes of the mesh.
#  -Laplace(u) = f    in the unit square
#            u = u_D  on the boundary
#  u_D = 1 + x^2 + 2y^2
#    f = -6
###########################################################

module ft01

using FenicsPy

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "P", 1)

# Define boundary condition
u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)

bc = DirichletBC(V, u_D, "on_boundary")

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = FeFunction(V)
solve(a == L, u, bc)

# Save solution to file in VTK format
vtkfile = File("poisson/solution.pvd")
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, "L2")

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = max(abs.(vertex_values_u_D - vertex_values_u)...)

# Print errors
println("error_L2  =", error_L2)
println("error_max =", error_max)

# Plot solution and mesh
#plot(u)
#plot(mesh)

end # module ft01


text:
#=
MIT License

Copyright (c) 2017 Joseph R. Peterson

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

# Define module name
module PlogP

# External modules
using DataTables

# Export module functions
export plogp

"""
	plogp(dt, cv_threshold)

	dt - DataFrame containing the expression value with columns as experiments and rows as genes. The first column should contain the gene identifiers.
	cv_threshold - Throw away genes with coefficients of vartion below this value (Default: 0).
	eig_threshold - Threshold for minimum eigenvalue to consider significant (Default: 1e-12; works well for E. coli example data)

Computes the gene interaction matrix (inverse of covariance matrix) using the principle of entropy maximization[^1]:
``\rho(\vec{x})=Ae^{\frac{\vec{x}M\vec{x}}{2}}``

[^1]T.R. Lezon, J.R. Banavar, M. Cieplak, A. Maritan, N.V. Fedoroff (2006) Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns, PNAS, 103(50):19033-19038.
"""
function plogp(dt::DataTable, cv_threshold::Float64=0.0, eig_threshold::Float64=1.0e-12)
	# Extract data
	geneNames = dt[:,1]
	rawData   = convert(Array,dt[:,2:end])

	geneCount   = size(rawData)[1]
	sampleCount = size(rawData)[2]

	# Normalize rows to unit variance
	means = mean(rawData,2)
	stds  = std(rawData,2)
	transformedData = (rawData.-means)./stds

	# Prefilter data if necessary
	cvs = stds./means
	filteredData = transformedData[(cvs.>cv_threshold)[:,1],1]
	sigGeneCount = size(filteredData)[1]
	sigGenes = geneNames[(cvs.>cv_threshold)[:,1]]

	# Compute covariance matrix
	covariance = cov(filteredData',filteredData')

	# Compute spectral decomposition
	eigenValues, eigenVectors = eig(covariance)

	# Get positive guaranteed eigenvalues/eigenvectors
	countPos = length(eigenValues[real(eigenValues).>eig_threshold])
	posEV   = real(eigenValues[end-(countPos-1):end])
	posEVec = real(eigenVectors[end-(countPos-1):end])

	# Compute partial inverse
	interaction = zeros(sigGeneCount, sigGeneCount)
	for i in 1:sigGeneCount
		for j in 1:sigGeneCount
			for k in 1:sampleCount
				interaction[i,j] += 1.0/posEV[k] * posEVec[k,i]*posEVec[k,j]
			end
		end
	end

	return sigGenes, interaction
end

"""
	createNetwork
"""
function createNetwork()

end

# Module end
end


text:
# -------------------- Dreaming Utilities --------------------

"""
    make_step(img, iterations, η, solo_call = false; path = "")

The core function of the deepdream. This is responsible for the forward
pass of the image through the model and generate the dream. At its core it
maximizes the `L2 Norm` of the image.

Arguments:
1. `img`: An array of the image on which to dream.
2. `iteratons`: Number of iterations. Ideally should be close to 10
3. `η`: Learning Rate. Taking a very high η will lead to divergence.
        Something close to 0.03 is recommended.
4. `solo_call`: Set it to `true` if this function is called indepently.
5. `path`: Keyword argument which **MUST** be specified if `solo_call`
           is set to `true`. This is the path where the final image is
           saved.
"""

function make_step(img, iterations, η, solo_call = false; path = "")
    if(solo_call && path == "")
        error("Image Save Path must be specified for solo calls")
    end
    input = param(img)
    for i in 1:iterations
        out = model(input)
        Flux.back!(out, out.data)
        input.data .= input.data + η * input.grad / mean(abs.(input.grad))
        info("$i iterations complete")
    end
    if(solo_call)
        save_image(path, input.data)
    end
    return input.data
end

"""
    deepdream(base_img, iterations, η, octave_scale, num_octaves, path, guide = 1.0)

Iteratively calls the `make_step` or the `guided_step` function depending on the
arguments. First it generates the `octaves` which are essentially `centre-zoomed`
parts of the `base_img`.

Arguments:
1. `base_img`: An array of the image on which to dream.
2. `iteratons`: Number of iterations. Ideally should be close to 10.
3. `η`: Learning Rate. Taking a very high η will lead to divergence.
        Something close to 0.03 is recommended.
4. `octave_scale`: Amount of zoom to be applied to the base image.
        Ideally taken to abe around 2.0.
5. `num_octaves` : Total number of octaves to be generated.
6. `path`: This is the path where the final image is saved.
7. `guide`: Pass the guiding image array.
"""

function deepdream(base_img, iterations, η, octave_scale, num_octaves, path, guide = 1.0)
    octaves = [copy(base_img)]
    for i in 1:(num_octaves-1)
        push!(octaves, zoom_image(octaves[end], octave_scale, octave_scale))
    end
    detail = zeros(octaves[end])
    out = base_img
    for (i, oct) in enumerate(octaves[length(octaves):-1:1])
        info("OCTAVE NUMBER = $i")
        w, h = size(oct)[1:2]
        if(i > 1)
            w1, h1 = size(detail)[1:2]
            detail = zoom_image(detail, w1 / w, h1 / h)
        end
        input_oct = (oct + detail) |> gpu
        if(guide != 1.0)
            out = guided_step(input_oct, guide, iterations, η)
        else
            out = make_step(input_oct, iterations, η)
        end
        detail = out - oct
    end
    save_image(path, out)
end

"""
    guided_step(img, guide, iterations, η, solo_call; path = "")

Core function to generate guided dreams. The `deepdream` function falls back
to this function to generate the dream using the guide. Internally it calculates
the similarity between the img feature vector and the guide feature vector by
computing their dot product and then maximizes the activation corresponding to
those features.

Arguments:
1. `img`: An array of the image on which to dream.
2. `guide`: The guiding image array.
3. `iteratons`: Number of iterations. Ideally should be close to 10
4. `η`: Learning Rate. Taking a very high η will lead to divergence.
        Something close to 0.03 is recommended.
5. `solo_call`: Set it to `true` if this function is called indepently.
6. `path`: Keyword argument which **MUST** be specified if `solo_call`
           is set to `true`. This is the path where the final image is
           saved.
"""

function guided_step(img, guide, iterations, η, solo_call = false; path = "")
    if(solo_call && path == "")
        error("Image Save Path must be specified for solo calls")
    end
    input = param(img)
    out = model(input)
    guide_features = model(guide)
    ch = size(out.data, 3)
    y = reshape(guide_features, :, ch)
    for i in 1:iterations
        out = model(input)
        x = reshape(out.data , :, ch)
        dot_prod = (x * y') |> cpu
        result = y[argmax(dot_prod, 2), :]
        result = reshape(result, size(out.data))
        Flux.back!(out, result)
        input.data .= input.data + η * input.grad / mean(abs.(input.grad))
        info("$i iterations complete")
    end
    if(solo_call)
        save_image(path, input.data)
    end
    return input.data
end

"""
    recurdream(base_img, iterations, η, octave_scale, num_octaves, num_frames, path, guide = 1.0; s = 0.05)

Iteratively calls the `deepdream` function. Also creates a global array
`frames` which stores the generated frames.

Arguments:
1. `base_img`: An array of the image on which to dream.
2. `iteratons`: Number of iterations. Ideally should be close to 10.
3. `η`: Learning Rate. Taking a very high η will lead to divergence.
        Something close to 0.03 is recommended.
4. `octave_scale`: Amount of zoom to be applied to the base image.
        Ideally taken to abe around 2.0.
5. `num_octaves` : Total number of octaves to be generated.
6. `num_frames`: Total number of frames to be generated
7. `path`: This is the path where the final image is saved.
8. `guide`: Pass the guiding image array.
9. `s`: Zoom scale
"""

function recurdream(base_img, iterations, η, octave_scale, num_octaves, num_frames, path, guide = 1.0; s = 0.05)
    global frames = []
    for i in 1:num_frames
        path_save = "." * *(split(path, ".")[1:end-1]...) * "_$i." * split(path, ".")[end]
        base_img = deepdream(base_img, iterations, η, octave_scale, num_octaves, path_save, 1.0)
        push!(frames, imresize(load(path_save), 256, 256))
        η /= 1.005
        info("Image saved at $path_save")
        base_img = zoom_image(base_img, s, s)
    end
end


text:
#include("ER_graphs.jl")
#include("utils.jl")
#using LoopVectorization
using Combinatorics
using Serialization

function init_mp(PG::Pair_ER)
    ones(2*ne(PG.N1.G),2*ne(PG.N2.G))
end

function sum_combi(neig_i1::Array{Int64,1},neig_i2::Array{Int64,1},k::Int64,mp_previous::Array{Float64,2})
    if  k == 0
        return Float64(1)
    end
    σ1 = combinations(neig_i1,k)
    σ2 = combinations(neig_i2,k)
    sum_inj = Float64(0)
    for s1 in σ1
        for s2 in σ2
            for ps2 in permutations(s2)
                sum_inj += prod(mp_previous[s1[j],ps2[j]] for j in 1:k)
            end
        end
    end
    sum_inj
end

function m_passing(PG::Pair_ER,e1::Int64,e2::Int64,mp_previous::Array{Float64,2})
    # trees T_1(e1) and T_2(e2)
    d1 = PG.N1.deg[e1]
    d2 = PG.N2.deg[e2]
    deg_min = min(d1,d2)
    if deg_min == 1
        output = psi(PG,0,d1-1,d2-1)
    else
        neig_i1 = PG.N1.neig[e1]
        neig_i2 = PG.N2.neig[e2]
        output = sum((psi(PG,k,d1-1,d2-1) * sum_combi(neig_i1,neig_i2,k,mp_previous)) for k in 0:(deg_min-1))
    end
    output
end

@inline function compute_next!(PG::Pair_ER,mp_previous::Array{Float64,2},mp_next::Array{Float64,2},ne1::Int64,ne2::Int64)
    @inbounds for e1 in 1:ne1
        @inbounds for e2 in 1:ne2
            mp_next[e1,e2] = m_passing(PG,e1,e2,mp_previous)
        end
    end
    nothing
end

function run_bp(PG::Pair_ER,n_iter::Int64;mp_previous=nothing,verbose=false)
    if mp_previous === nothing
        mp_previous = init_mp(PG)
    end
    mp_next = similar(mp_previous)
    mp_opt = similar(mp_previous)
    (ne1,ne2) = size(mp_previous)
    res = zeros(6,n_iter)
    perf = 0
    last_iter = 0
    ov_opt = 0
    for i in 1:n_iter
        compute_next!(PG,mp_previous,mp_next,ne1,ne2)
        (mp_previous, mp_next) = (mp_next, mp_previous)
        M_lr = create_matrix_lr(PG,mp_previous)
        M_loglr = log.(M_lr)
        ov1, ov2, v1,v2 = eval_M(PG, M_loglr)
        m1,m2 = eval_edges(PG, M_loglr)
        new_perf = (m1+m2)/2
        if new_perf > perf
            perf = new_perf
            mp_opt = mp_previous
            last_iter = i
            ov_opt = (ov1+ov2)/2
        end
        res[1,i] = ov1
        res[2,i] = ov2
        res[3,i] = v1
        res[4,i] = v2
        res[5,i] = m1
        res[6,i] = m2
        if verbose
            println(i, " | ", m1, " | ", m2, " | " ,ov1," | ",ov2, " | ", v1, " | ", v2)
        end
        if max(v1,v2) == Inf
            break
        end
    end
    res, mp_opt, last_iter, perf, ov_opt/PG.n
end

function lr_root(PG::Pair_ER,mp,d1,d2,neig_i1,neig_i2)
    sum = Float64(0)
    for k in 0:min(d1,d2)
        sum += psi(PG,k,d1,d2) * sum_combi(neig_i1,neig_i2,k,mp)
    end
    sum
end

function create_matrix_lr(PG::Pair_ER,mp)
    v1 = vertices(PG.N1.G)
    v2 = vertices(PG.N2.G)
    matrix_lr = zeros(length(v1),length(v2))
    deg1, neig1 = neighbors_edges(PG.N1.G)
    deg2, neig2 = neighbors_edges(PG.N2.G)
    for i in v1
        for j in v2
            matrix_lr[i,j] = lr_root(PG,mp,deg1[i],deg2[j],neig1[i],neig2[j])
        end
    end
    matrix_lr
end

function init_mp_biased(PG::Pair_ER,bias=10)
    mp_init = ones(2*ne(PG.N1.G),2*ne(PG.N2.G))
    for (i,j) in eachrow(PG.edge_G1G2)
        mp_init[i,j] = bias
    end
    mp_init
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 438caa30-66e8-11eb-31e9-917e458e4d33
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="StatsBase", version="0.33"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="HypertextLiteral", version="0.6"),
    ])
    using Plots, Random, StatsBase, PlutoUI, LinearAlgebra, SparseArrays, Markdown, InteractiveUtils, HypertextLiteral
    md"Packages"
end

# ╔═╡ f36826be-93cd-11eb-3cd4-278a16171c91
	TableOfContents(aside = true)

# ╔═╡ 24d68cf0-7ed0-11eb-004b-7b2702a5de1a
md"
# Odyssey across the ocean ⛵
"
# CC-BY-4.0 (IMOOX, TU Graz, Institute of Theoretical and Computational Physics)
# Authors: Johanna Moser, Gerhard Dorn, Wolfgang von der Linden

# ╔═╡ b4dc3a2f-f407-4b8c-a448-4849a914be99
md"""
# About the creators

We, **Johanna Moser, Prof. Wolfgang von der Linden** and **Gerhard Dorn** created this notebook in the context of teaching **Bayesian probability theory** using story telling.

In order to make the content more approachable we invented the story of **Captain Bayes and her crew** who navigate on the ocean of uncertainty and experience a lot of adventures in probability theory.

The course is a free massive open online course (MOOC) available on the platform [`IMOOX`](https://imoox.at/mooc/local/landingpage/course.php?shortname=bayes&lang=en)

$(Resource("https://raw.githubusercontent.com/Captain-Bayes/images/main/adventure_map.gif"))

We even have a **trailer**:
"""

# ╔═╡ 97565f62-8129-4515-852e-2a47748dacce
html"""

"""

# ╔═╡ c6f9c954-54c4-48ce-8465-fc85202b79f4
begin

	md"""Switch to the next 8 days 📅 by clicking the up 🔼 button: 👉 $(@bind dial NumberField(0:8; default=0))"""
	#md"""$(@bind dial ClickCounterWithReset("Dial!", "Start over!"))"""
end

# ╔═╡ 89c46920-7ed0-11eb-3c5f-574d525a9f1f
md" Here again you can see our compass 🧭! Now you can change the probabilities by increasing the number next to the directions ⬅⬆➡⬇ and thus simulate wind or currents. Let's see how this might affect our journey!
"

# ╔═╡ 66966ea0-96d4-11eb-1c5c-2b06dfd9313b
begin
	hint(text, headline=md"Law of large numbers") = Markdown.MD(Markdown.Admonition("hint", string(headline), [text]));

	hint(md"You might have noticed that the more days our crew is sailing, the more similar the simulation is to the theoretical values. This is called the **law of large numbers**, and it tells us, that for an infinite number of days, the results of the experiment will perfectly coincide with the theoretical distribution.")
end

# ╔═╡ 2ea667f0-6f26-11eb-02fb-1335862dc98e
md" Where do you spot turtle island?

x coordinate = $(@bind x_island NumberField(0:6, default=5))

y coordinate = $(@bind y_island NumberField(0:6, default=0))"

# ╔═╡ 4ae39013-0e68-44a2-a124-402dc1f76557
md"**Question:** But before you manipulate the compass 🧭, what do you think, how do we have to change the compass probabilities to randomly reach turtle islans 🐢 (which lies 5 days in east direction) probably sooner?

👉 Increase the weight in East direction? $(@bind answer_1 CheckBox())

👉 Decrease the weights in North and South direction? $(@bind answer_2 CheckBox())

👉 Set West, North and South to zero? $(@bind answer_3 CheckBox())

👉 Increase West to the maximum? $(@bind answer_4 CheckBox())

"

# ╔═╡ 4e52d0f9-c2fd-4de4-b0c8-05a6b332c9bd
md"""Click here if you want to see the analytic solution using a Markov process 👉 $(@bind plot_exact CheckBox())"""

# ╔═╡ e4bbb3f0-e6f2-44a4-aa63-3bad12d14dc0
begin
	if plot_exact  == false
md" ## An exact solution for Turtle island!🐢 (hidden)"
	else md" ## An exact solution for Turtle island!🐢"
	end
end

# ╔═╡ ad8dfe0a-1184-4831-b33b-fcc42eb6e0c4
md"""Just click in the box to show Bernoulli

the cummulative probability distribution of first return to turtle island🐢 👉$(@bind cummulative_probability CheckBox())

and the probability for the position after one year 👉$(@bind prob_dist_one_year CheckBox())"""

# ╔═╡ 59a47ad8-ed45-4765-a8e2-28cd49ae0ab2
md"## Program code"

# ╔═╡ b34109a1-5762-46ac-b95c-5c103b1552bc
function sub2ind(siz, ix,iy)

	return siz[1]*(iy-1) + ix
end

# ╔═╡ bdc92620-ab48-4587-a333-1578ce8a4e17
begin

	q_left = 0.25
	q_right = 0.25
	q_up = 0.25
	q_down = 0.25
local N = 100
diag_up = [repeat([ones(N-1)*q_up;0],(N-1));ones(N-1)*q_up]
diag_down = [repeat([ones(N-1)*q_down;0],(N-1));ones(N-1)*q_down]

local C = spdiagm( 1=>diag_up, -1=>diag_down, N=> ones(N^2-N)*q_right, -N => ones(N^2-N)*q_left)

	#index of turtle island:
d = CartesianIndex(N÷2 + 5, N÷2)

local turtle_island = sub2ind([N,N], (N+1)÷2 + 5, (N+1)÷2)
origin = sub2ind([N,N], (N+1)÷2 , (N+1)÷2)
	# make turtle island a stop position:
	C[turtle_island, :] .= 0
	C[turtle_island, turtle_island] = 1

	# lost probability weights could be sent to a special state but this is not necessary, the border / edge is given by 1-sum(pi)
	steps = 365
	time_prob = zeros(steps,1)
	out_of_simulation = zeros(steps,1)

	let
		pi_vec= zeros(1,N^2)
		pi_vec[origin] = 1

		local C10 = C^10
		for t = 1:steps
			pi_vec = pi_vec * C

			time_prob[t] = pi_vec[turtle_island]
			out_of_simulation[t] = 1-sum(pi_vec)
		end
		global distri_markov = reshape(pi_vec, N,N)
	end

	#=end

	pdf = diff(time_prob,dims=1)
	md"theoretical values turtle island"

	=#

	md"Theoretical values turtle island Markov Chain"

end

# ╔═╡ 0838840f-c652-48ec-8a7c-5f01c6aaded5
if plot_exact == true
	pdf_return = diff(time_prob; dims=1)
	days_of_pos_return = [2:steps;]
	day_max_prob_return = days_of_pos_return[pdf_return[:] .== maximum(pdf_return)][1]
	max_prob_return = maximum(pdf_return)

	plot((2:1:steps) , pdf_return , line = (1.2, 0, :bar), xlim = [0,30], title="Reaching turtle island", xlabel="days", ylabel="probability of first arrival", label=:none)

		#plot!((1:1:steps) , out_of_simulation)
	end

# ╔═╡ 979bab1b-8463-4ca0-a506-d97ed2c29871
if cummulative_probability
	plot(1:steps, time_prob, line = (1, 1.0, :path), label=:none, xlabel="t | days", ylabel="probability of first arrival after t days", title="Cummulative probability")

end

# ╔═╡ d69ab95f-bd5c-4507-b24d-f12fbde4cae8
if prob_dist_one_year
heatmap(-49:50, -49:50, permutedims(distri_markov,[2,1]), clim= (0,0.0019), c = :dense, title="probability distribution after one year")
end

# ╔═╡ 264b06d9-6fea-4260-9b42-66feca66a654
begin
	# define images
bayes = "https://raw.githubusercontent.com/Captain-Bayes/images/main/bayes_50px.gif"
	ernesto_short = "https://raw.githubusercontent.com/Captain-Bayes/images/main/Ernesto_animated.gif"
	ernesto_completed = "https://raw.githubusercontent.com/Captain-Bayes/images/main/Ernesto_completed.gif"
	bernoulli = "https://raw.githubusercontent.com/Captain-Bayes/images/main/bernoulli_100px.gif"
	turtle_island = "https://github.com/Captain-Bayes/images/blob/main/island_in_sight.gif?raw=true"
	bottle = "https://raw.githubusercontent.com/Captain-Bayes/images/main/Flaschenpost_schwimmend.gif"
	desparate_bernoulli = "https://raw.githubusercontent.com/Captain-Bayes/images/main/Bernoulli_desperate.png"

	bayes_large = "https://raw.githubusercontent.com/Captain-Bayes/images/main/bayes_100px.gif"

	md"""Images"""

	#compass = Resource("https://raw.githubusercontent.com/Captain-Bayes/images/main/Kompass_empty.png")
end

# ╔═╡ 1aa40340-9382-11eb-1031-5fa5b0f133ee
begin
local	rng1 = MersenneTwister(32) #seperate rng so that it won't actualize again when you're further down
	dir_index = 5 .-sum(rand(rng1, 1, 8).<=cumsum([0.25, 0.25, 0.25, 0.25]), dims = 1)
local	K_start = [sum(dir_index[1:j] .== i) for i in 1:4, j in 1:8]
	pos_start = [0 0; K_start[1,:] - K_start[3,:] K_start[2,:] - K_start[4,:]]

	dial_index = Integer(round(dial,digits=0))

	md" ## Welcome
to this random journey! Each day we choose the sailing direction randomly. Hmm, let us try to use my famous random compass for this navigation!

Can you help me with the dial process? Just switch to the next days by clicking the up button 🔼, and see what happens!

$(Resource(bayes, :width => 200))"
	#Can you help me with to dial? Click the button below, and see what happens!

end

# ╔═╡ 5d25f3a0-93a4-11eb-3da6-c96ae54a0d70
begin

	if dial_index > 0 && dial_index < 9

	lim_start = maximum(abs.(pos_start[:]))
	plot(
			pos_start[1:dial_index+1,1], pos_start[1:dial_index+1,2],
			linecolor   = :green,
			linealpha = 0.2,
			linewidth = 2, aspect_ratio =:equal,
			marker = (:dot , 5, 0.2, :green),
			label=false,
			xlim =[-lim_start, lim_start],
		    ylim =[-lim_start, lim_start],
			legend = :bottom
			)
	plot!(
			[0],[0],
			marker = (:dot, 10, 1.0, :red),
			label = "initial position"
			)

	plot!(
			[pos_start[dial_index+1,1]], [pos_start[dial_index+1,2]],
			marker = (:circle, 10, 1.0, :green),
			label = "current position"
			)

	end
end

# ╔═╡ 38fd06d0-93cc-11eb-030e-a7888d7d7eee
begin
#dial first steps

local dir = ["E", "N", "W",  "S"]
local word = [ "east", "north", "west" ,"south"]
local url = ["https://raw.githubusercontent.com/Captain-Bayes/images/main/Kompass_east-export.gif", "https://raw.githubusercontent.com/Captain-Bayes/images/main/Kompass_north-export.gif",  "https://raw.githubusercontent.com/Captain-Bayes/images/main/Kompass_west-export.gif",  "https://raw.githubusercontent.com/Captain-Bayes/images/main/Kompass_south-export.gif"]

if dial_index < 8	&& dial_index > 0
	md"Well done! The compass needle landed on **$(dir[dir_index[dial_index]])**. Seems like we'll be heading **$(word[dir_index[dial_index]])wards** today! $(Resource(url[dir_index[dial_index]], :width => 200))"
elseif dial_index >= 8
		md"""
		$(Resource(ernesto_short, :width => 30))

		**Thank you for helping Captain Bayes dial the compass! From now on, she can handle it on her own. Scroll down further to see the whole journey of our crew. You can also change the seed to see different possible journeys!** """
	end

end

# ╔═╡ 9d9726dd-3456-4988-ae96-c25129092c39
md"""
That's cool! But how will we know where we'll (most likely) end up?

$(Resource(bernoulli, :width => 200))

Click here to find out! 👉 $(@bind see_distribution CheckBox())
"""

# ╔═╡ 702cce30-8047-11eb-015b-c5fb0509a38f
begin
	if see_distribution == false
	md" ## Probability distribution (hidden)"
	else
		md" ## Probability distribution"
	end
end

# ╔═╡ 299d527e-96e7-11eb-32c5-05fdf76bb79f
md"$(Resource(ernesto_completed, ))

Congrats on completing this section! You really deserve a pause! 🍍🍎🥕🥛🍵

When you're ready again, click here 👉 $(@bind see_turtle_island CheckBox())!
"

# ╔═╡ 9579e880-8047-11eb-25cc-1710d87cbd23
begin
	if see_turtle_island  == false
md" ## Let's go to Turtle island!🐢 (hidden)"
	else md" ## Let's go to Turtle island!🐢"
	end
end

# ╔═╡ 11336f00-96e6-11eb-2434-8bac971b9849
md" $(Resource(turtle_island, :width => 700))

 **Bernoulli:** Look, there is a giant turtle over there!

**Bayes:** This must be turtle island! I  hear they have the most delicious fresh lemonade there!

**Bernoulli:** Oh how wonderful, I really long for some fresh lemonade! Let's just go there!

**Bayes:** Wait ... we are currently  analyzing our journey with the random compass 🧭, and I don't want to stop this experiment early. And sooner or later we will reach the island for sure. Bernoulli, why don't you try to figure out how long it will take to reach that island 🐢🏝 on average?
Just perform **random walks** 🤪🌀 on your ship map and count how often we hit the island.
You could then produce a histogram showing us the statistical likelihood of reaching the island within a certain amount of time. The median of the arrival times of your simulated journeys should be an appropriate benchmark." #formulierung

# ╔═╡ f16d2de5-2473-432d-9b8e-c0482fa74031
md" $(Resource(desparate_bernoulli, :width=>100)) Oh that's sooo long 😰. But if the wind and the currents were different, we could get there a lot faster! Oh dear reader, will you tinker once again with the probabilities of our compass, so I can get my lemonade faster? Don't tell the Captain though!🤫
"

# ╔═╡ 5e7487fb-cab0-4bae-ae33-7537469d530e
md" $(Resource(bayes_large, :width=>180)) **Bernoulli, there must be a more accurate solution!**

A friend of mine - Markov - has some brilliant ideas - let's see if we can apply them.
The probabilities for the next day should depend on those of the previous. If we just multiply this vector with this transfer matrix...

 >$\vec{P}^{(t+1)} = M \cdot \vec{P}^{(t)}$

Bernoulli, just have a look 😀

"

# ╔═╡ 454efd23-d99a-4d99-b348-1336dfe294a5
md"$(Resource(bernoulli, :width=>180)) Brilliant, so it seams the highest probability to reach turtle is on **the $(day_max_prob_return) day** with a **probability of $(round(max_prob_return; digits=4))**...

Puh, that's not much! Oh, I what may be the chances that I get there within a year?

Can you show me the cummulative probabilities and where we might be after one year?
"

# ╔═╡ e130ac04-e3eb-4be5-ae6c-c87eaf7064a5
begin
	days_max_first_journey = 200
	days_max = 400
	n_reps_max = 2000
	days_slider = @bind days Slider(1:1:days_max_first_journey, show_value = true, default = 100)
	days_slider_2 = @bind days_2 Slider(1:25:days_max, show_value = true, default = 100)
	seed_slider = @bind seed NumberField(1:100, default = 20)
	seed_slider_2 = @bind seed_2 NumberField(1:100, default = 20)

	md"define sliders"
end

# ╔═╡ 4d6c28d0-70a9-11eb-0a98-5583ef673517
md"Now you found an old log book 📘 of a previous random walk noted by Bernoulli. Just choose the number of 👉 **$(days_slider) days** to see the path of the odyssey of our crew on the ocean.

The darker a green circle 🟢 the more often it has been visited.
"

#@bind days_clock Clock(0.3, true)

# ╔═╡ ce4ef0c2-81d1-11eb-05fe-e590b8cf2191
md" You can also examine another random walk/odyssey by taking another logbook 📕📗📙 out of the shelf in the captain's cabin. The seed  👉 $(seed_slider) is like the ID of the book"

# ╔═╡ a3a332d0-6d8d-11eb-2157-61140f731b59
md"""**Where do we end?** By repeating our odyssee again and again we can find a probability distribution about where the ship will end up after a certain number of days. The starting point will always be [0,0] and we'll mark the endpoint after each journey. See how the distribution changes depending on the number of **repetitions** 👉 $(@bind nr_reps Slider(100:100:n_reps_max, show_value = true, default = 100)) times,

and the number of **days** 👉 $days_slider_2 days spent travelling!

Does it look similar to the spilled **ink stain** or rather a **chessboard pattern** ♟?
Since on even and odd days we can just reach half of the possitions you can chosse the option to average two consecutive days 👉 $(@bind averaged_final CheckBox()).

You can also simulate wind or currents by changing the probabilities on the compass 🧭."""

# ╔═╡ 9d704724-4ddc-4ef2-a8f2-ce58af8f2339
md"""So if you wish you can change the seed 👉 $(seed_slider_2) and check what median arrival time ⏱ another random sample yields"""

# ╔═╡ ae845910-8109-11eb-39a8-0182f17e791e
begin
hide_everything_below =
	html"""

	""";

md"definition hide everything below"
end

# ╔═╡ 4f304620-93cb-11eb-1da6-739664f2a105
begin
	if dial < 8
			hide_everything_below
	end
end

# ╔═╡ 81cfef50-93d4-11eb-3448-975c908bd1a2
begin
if see_distribution == false
		hide_everything_below
	end
end

# ╔═╡ e527cc7e-96e5-11eb-029c-277bfb1730dc
begin
if see_turtle_island == false
	hide_everything_below
	end
end

# ╔═╡ 0446b6c1-1a7b-4ac0-b50e-6bc7e3b60c72
begin
if plot_exact == false
		hide_everything_below
	end
end

# ╔═╡ 49b1e2ad-2129-4562-911f-a81976a6bd55
html"""

"""

# ╔═╡ 69d1a4d0-96c6-11eb-002f-9138e617a1c2
begin
	see_distribution
	# used to reset the compass to make it fair again, when entering the next section

	W1 = @bind W Scrubbable(0:1:3, default=1)
	N1 = @bind N Scrubbable(0:1:3, default=1)
	E1 = @bind E Scrubbable(0:1:3, default=1)
	S1 = @bind S Scrubbable(0:1:3, default=1)

	md"define tablestyle"
end

# ╔═╡ 8510bdd0-96c6-11eb-3a9a-bd311edac8f4
@htl("""

            	$(N1)

            $(W1)

            $(E1)

            	$(S1)

""")

# ╔═╡ 815094e0-93ce-11eb-3878-3be919148949
begin
	angles = [0.0, pi/2, pi, 3*pi/2, 0.0]
	numbers = [E, S, N, W]
	su = sum(numbers)
	weighted = [E/su, N/su, W/su, S/su, E/su]
	plot(angles, weighted, proj=:polar, m=2, label = "weights")
	#warum so wahnsinnig langsam?
end

# ╔═╡ c3a505f0-66e8-11eb-3540-69ce34959d64
begin

#probabilities for the four directions
prob = weighted[1:4]
# create all random variables using the seed defined above
rng = MersenneTwister(seed_2)

local sample_temp = 5 .- sum(rand(rng, 1,days_max*n_reps_max) .<= cumsum(prob), dims=1)
local sample_array = reshape(sample_temp, (days_max,:))

# first calculation: positions of first run -> moved to first section!
# K_first_run = [sum(sample_array[1:i,1] .== j) for j ∈ 1:4, i ∈ 1:days_max]
#pos_first_run = [0 0; K_first_run[1,:] - K_first_run[3,:]  K_first_run[2,:] - K_first_run[4,:]]

# second calculation: final positions
local K = [sum(sample_array[1:days_2,i] .== j ) for j ∈ 1:4, i ∈ 1:n_reps_max]
final_pos = [K[1,:] - K[3,:]  K[2,:] - K[4,:]]

local K2 = [sum(sample_array[1:days_2-1,i] .== j ) for j ∈ 1:4, i ∈ 1:n_reps_max]
final_pos_2 = [K2[1,:] - K2[3,:]  K2[2,:] - K2[4,:]]

# third calculation: return to turtle island

#flip days_max and repetitions, we only take 1000 repetitions but let the walker go 10000 days (overflow)
overfl = n_reps_max

# since approx 50% of the walker return before 10000 days we will get 2*days_max repetitions
max_runs = Integer(days_max*1.5)

local x_logic = (sample_temp.==1) - (sample_temp.==3)
local y_logic = (sample_temp.==2) - (sample_temp.==4)

# those vectors	will be used to display histogram
reached = zeros(max_runs)
not_reached = zeros(max_runs)
let
	ind = 0
	step = 1
	for k= 1:max_runs

L_reach_island = findall((cumsum(x_logic[ind .+ step .* (1:overfl)]) .== x_island) .& (cumsum(y_logic[ind .+ step .* (1:overfl)]) .== y_island))
		if isempty(L_reach_island) # not reached (overflow)
			not_reached[k] = 1
			ind = ind + overfl
		else 			# reached
			reached[k] =  minimum(L_reach_island)
			ind = Integer(ind + reached[k])
		end
		# if more than 50% of walkers for some random seed do not come back we set index to zero and take only every second value to generate a new sample set (backup to prevent error)
		if ind + overfl > length(sample_temp)
			ind = 0
			step = 2
		end

	end
end

	days_histogram = reached
	days_histogram[reached .== 0] .= overfl

	md" probability distribution and turtle island walker"
end

# ╔═╡ 8aeb5aa0-93d1-11eb-1935-9b20967fb2e3
begin
	rng2 = MersenneTwister(seed)

local sample_array = 5 .- sum(rand(rng2, 1,days_max) .<= cumsum(prob), dims=1)

# first calculation: positions of first run:
 K_first_run = [sum(sample_array[1:i] .== j) for j ∈ 1:4, i ∈ 1:days_max]
pos_first_run = [0 0; K_first_run[1,:] - K_first_run[3,:]  K_first_run[2,:] - K_first_run[4,:]]

	md""" Now let's compare the directions chosen by our compass with the theoretical distribution we chose. What do you notice?"""

#how to blend in law of large numbers
end

# ╔═╡ d4a2ab30-7f4a-11eb-0a76-d50b21c3217b
begin
	median_island = StatsBase.median(days_histogram)

	plot1 = histogram(days_histogram, bins = 100, label = :none, xlabel = "days", ylabel = "occurrences")

	plot!([StatsBase.median(days_histogram)],[0], marker = "red", label = string("Median =",StatsBase.median(days_histogram)), legend=:top)
# calculate median manually:
	hist_return = [sum(days_histogram .== i) for i in 1:overfl]

	#med = minimum(findall(cumsum(hist_return./max_runs) .>= 0.5))

	#einfügen: theoretische kurve?

	plot1
end

# ╔═╡ 6cbe3250-805d-11eb-0f34-43f4c1669537
md"$(Resource(bernoulli, :width=>180)) The simulation is finished, Captain! Here's the histogram. Taking the median as a measure of average, it looks like we'll arrive in about $(median_island) days...

"

# ╔═╡ 10b945c9-4946-49bf-9e16-62b27b3766d7
md"""If you are still not happy with the expected time to reach turtle island, you could say our simulation was just bad luck, we only use **$(max_runs) runs** and **$(overfl) days** for each run to get the average return statistic, so maybe the true value is different 🤔."""

# ╔═╡ 76380dec-000c-43d6-957f-4fb156846ff9
@htl("""

            	$(N1)

            $(W1)

            $(E1)

            	$(S1)

""")

# ╔═╡ 55e65cf6-7f5c-4faa-b271-cb78761300aa
begin
#define variables

first_steps_x = [0]
first_steps_y = [0]

	x0 = [0,0]

compass_dict = Dict("N"=> 1, "E" => 2, "S" => 3, "W" => 4)

compass = ["E", "N", "W", "S"] #possible directions
compass_numbers = [1, 2, 3, 4]
times_compass = [0, 0, 0, 0] #counts times every direction NESW is chosen
actual_directions = [[1, 0], [0, 1], [-1, 0],  [0, -1] ]
#calculate changes to probability:
norm = N + S + E + W
n = N/norm
e = n + E/norm
s = e + S/norm
w = s + W/norm
weights = [n, E/norm, S/norm, W/norm]
	md"variables"
end

# ╔═╡ 5cea77d0-93d1-11eb-1508-aff9495e46d8
begin
	dial
if length(first_steps_x) < 7
md"""
## One journey (hidden)"""
	else
		md"""
## One journey"""
end
end

# ╔═╡ 5b711a00-6d8c-11eb-00ac-4dd20bc3dcc6
begin
	max_x_1_run = maximum(abs.(pos_first_run[1:days,1]))
	max_y_1_run = maximum(abs.(pos_first_run[1:days,2]))
	plot(
			pos_first_run[1:days,1], pos_first_run[1:days,2], linecolor   = :green,
			linealpha = 0.2,
			linewidth = 2, aspect_ratio =:equal,
			marker = (:dot , 5, 0.2, :green),
			label=false,
			xlim = [-max_x_1_run, max_x_1_run],
			ylim = [-maximum([max_x_1_run*0.6, max_y_1_run]), maximum([max_x_1_run*0.6, max_y_1_run])]
			)
	plot!(
			[x0[1]],[x0[2]],
			marker = (:dot, 10, 1.0, :red),
			label = "initial position"
			)

	plot!(
			[pos_first_run[days,1]], [pos_first_run[days,2]],
			marker = (:circle, 10, 1.0, :green),
			label = "current position"
			)
end

# ╔═╡ 29523c50-7554-11eb-25c1-1b56caf928c5
begin
plot(compass, K_first_run[:,days]/days,
		line = (1., 1., :bar), label = "simulation", title = "cardinal directions chosen")
plot!(compass, weights, line = (1.0, 0.0, :bar), bar_width = 0.02,
    marker = (:circle, 50, 1), color = [:red], label = :none, legend = :right)

plot!(compass, weights,
		line = (0., 0, :path),
    normalize = false,
    bins = 10,
	bar_width = 0.2,
    marker = (7, 1., :o),
    markerstrokewidth = 1,
    color = [:red],
    fill = 1.,
    orientation = :v,
	ylabel = "Relative frequency",
	xlabel = "Directions",
	label = "theory")
end

# ╔═╡ 9600df90-7f46-11eb-2d6f-953d8166854e
begin

	#shifts the position [0,0] to the middle of a matrix hist_data
	max_distance_1 = maximum(abs.([final_pos[1:nr_reps,1]; final_pos[1:nr_reps,2]]))
	max_distance_2 = maximum(abs.([final_pos_2[1:nr_reps,1]; final_pos_2[1:nr_reps,2]]))
	max_distance = maximum([max_distance_1, max_distance_2])
	# initialize with a sparse matrix, plot with a full matrix - Array command
	hist_data = sparse(final_pos[1:nr_reps,2].+max_distance.+1, final_pos[1:nr_reps,1].+max_distance.+1, 	 	 ones(size(final_pos[1:nr_reps],1)), 2*max_distance + 1, 2*max_distance + 1)

	hist_data_2 = sparse(final_pos_2[1:nr_reps,2].+max_distance.+1, final_pos_2[1:nr_reps,1].+max_distance.+1, 	 	 ones(size(final_pos_2[1:nr_reps],1)), 2*max_distance + 1, 2*max_distance + 1)

	if averaged_final
	heatmap(-max_distance:1:max_distance, -max_distance:1:max_distance, Array(hist_data + hist_data_2)./2, seriestype = :bar, c=:dense, axis = :equal)
	else
	heatmap(-max_distance:1:max_distance, -max_distance:1:max_distance, Array(hist_data), seriestype = :bar, c=:dense)
	end

	L_phi  = [0:.01:2*pi;]
	radius = sqrt(days_2) * sqrt(pi)/2
	Lc_x    = radius * cos.(L_phi) .+ days_2*(weights[2] - weights[4])
	Lc_y    = radius * sin.(L_phi) .+ days_2*(weights[1] - weights[3])
	plot!(Lc_x,Lc_y, linewidth = 2, label = "measure for mean distance")
#here are different color schemes: https://docs.juliaplots.org/latest/generated/colorschemes/

	#heatmap(x_array, y_array, final_position, seriestype = :bar)
	#plot!([0],[0], marker = (:dot, 10, 1.0, :red))
end

# ╔═╡ 073ac878-52dd-4112-9b42-ad08649fe927
ClickCounterWithReset(text="Click", reset_text="Reset") = HTML("""

<button>$(text)</button>
$(reset_text)

""")

# ╔═╡ 8e804243-9123-497d-a4b2-552f04c1d9d5
begin
almost(text, headline=md"Almost there!") = Markdown.MD(Markdown.Admonition("warning", string(headline), [text]));
#brown

correct(text=md"Great! You got the right answer!", headline=md"Got it!") = Markdown.MD(Markdown.Admonition("correct", string(headline), [text]));
#green

keep_working(text=md"The answer is not quite right.", headline=md"Keep working on it!") = Markdown.MD(Markdown.Admonition("danger", string(headline), [text]));
#red
md"admonitions"
end

# ╔═╡ 473d0dab-ca56-4ce7-8f0e-7a8436ea5833
begin
	if answer_1
		keep_working(md"""Rethink your answer or try it out, you will see that the "ink stain" moves far to the right ➡ missing turtle island!""", md"Not really!")
	elseif answer_2
		correct( md"You surly realized that setting North and South to zero will lead to a one dimensional random walk which will highly increase the probability to reach turtle island soon 🐢", md"""Clever!""")
	elseif answer_3
		almost(md"Well, you can try it but Captain Bayes will find out, you killed all randomness. Your manipulation will directly steer the ship onto turtle island, so you definitly will reach it in 5 days, but that's no random walk anymore!.", "This is kind of cheating")
	elseif answer_4
		keep_working(md"Maybe you are like Magellan and hope to circumnavigate the globe to tackle turtle island from a direction it would not expect^^. Well played, but this definitly takes tooo long. And be careful in a pure Eucledian space setting East to zero would destroy your dream from ever drinking lemonade on turtle island beach.", md"Wrong direction!")
	end
	#almost, correct, keep_working

end

# ╔═╡ Cell order:
# ╟─f36826be-93cd-11eb-3cd4-278a16171c91
# ╟─24d68cf0-7ed0-11eb-004b-7b2702a5de1a
# ╟─b4dc3a2f-f407-4b8c-a448-4849a914be99
# ╟─97565f62-8129-4515-852e-2a47748dacce
# ╟─1aa40340-9382-11eb-1031-5fa5b0f133ee
# ╟─c6f9c954-54c4-48ce-8465-fc85202b79f4
# ╟─38fd06d0-93cc-11eb-030e-a7888d7d7eee
# ╟─5d25f3a0-93a4-11eb-3da6-c96ae54a0d70
# ╟─4f304620-93cb-11eb-1da6-739664f2a105
# ╟─5cea77d0-93d1-11eb-1508-aff9495e46d8
# ╟─4d6c28d0-70a9-11eb-0a98-5583ef673517
# ╟─ce4ef0c2-81d1-11eb-05fe-e590b8cf2191
# ╟─5b711a00-6d8c-11eb-00ac-4dd20bc3dcc6
# ╟─89c46920-7ed0-11eb-3c5f-574d525a9f1f
# ╟─8510bdd0-96c6-11eb-3a9a-bd311edac8f4
# ╟─815094e0-93ce-11eb-3878-3be919148949
# ╟─8aeb5aa0-93d1-11eb-1935-9b20967fb2e3
# ╟─29523c50-7554-11eb-25c1-1b56caf928c5
# ╟─66966ea0-96d4-11eb-1c5c-2b06dfd9313b
# ╟─9d9726dd-3456-4988-ae96-c25129092c39
# ╟─81cfef50-93d4-11eb-3448-975c908bd1a2
# ╟─702cce30-8047-11eb-015b-c5fb0509a38f
# ╟─a3a332d0-6d8d-11eb-2157-61140f731b59
# ╟─9600df90-7f46-11eb-2d6f-953d8166854e
# ╟─299d527e-96e7-11eb-32c5-05fdf76bb79f
# ╟─e527cc7e-96e5-11eb-029c-277bfb1730dc
# ╟─9579e880-8047-11eb-25cc-1710d87cbd23
# ╟─11336f00-96e6-11eb-2434-8bac971b9849
# ╟─2ea667f0-6f26-11eb-02fb-1335862dc98e
# ╟─6cbe3250-805d-11eb-0f34-43f4c1669537
# ╟─d4a2ab30-7f4a-11eb-0a76-d50b21c3217b
# ╟─f16d2de5-2473-432d-9b8e-c0482fa74031
# ╟─4ae39013-0e68-44a2-a124-402dc1f76557
# ╟─473d0dab-ca56-4ce7-8f0e-7a8436ea5833
# ╟─76380dec-000c-43d6-957f-4fb156846ff9
# ╟─10b945c9-4946-49bf-9e16-62b27b3766d7
# ╟─9d704724-4ddc-4ef2-a8f2-ce58af8f2339
# ╟─5e7487fb-cab0-4bae-ae33-7537469d530e
# ╟─4e52d0f9-c2fd-4de4-b0c8-05a6b332c9bd
# ╟─0446b6c1-1a7b-4ac0-b50e-6bc7e3b60c72
# ╟─e4bbb3f0-e6f2-44a4-aa63-3bad12d14dc0
# ╟─0838840f-c652-48ec-8a7c-5f01c6aaded5
# ╟─454efd23-d99a-4d99-b348-1336dfe294a5
# ╟─ad8dfe0a-1184-4831-b33b-fcc42eb6e0c4
# ╟─979bab1b-8463-4ca0-a506-d97ed2c29871
# ╟─d69ab95f-bd5c-4507-b24d-f12fbde4cae8
# ╟─59a47ad8-ed45-4765-a8e2-28cd49ae0ab2
# ╟─b34109a1-5762-46ac-b95c-5c103b1552bc
# ╟─bdc92620-ab48-4587-a333-1578ce8a4e17
# ╟─264b06d9-6fea-4260-9b42-66feca66a654
# ╟─55e65cf6-7f5c-4faa-b271-cb78761300aa
# ╟─e130ac04-e3eb-4be5-ae6c-c87eaf7064a5
# ╟─c3a505f0-66e8-11eb-3540-69ce34959d64
# ╟─438caa30-66e8-11eb-31e9-917e458e4d33
# ╟─ae845910-8109-11eb-39a8-0182f17e791e
# ╟─49b1e2ad-2129-4562-911f-a81976a6bd55
# ╟─69d1a4d0-96c6-11eb-002f-9138e617a1c2
# ╟─073ac878-52dd-4112-9b42-ad08649fe927
# ╟─8e804243-9123-497d-a4b2-552f04c1d9d5

max_stars_repo_path: Interactive Story Telling.jl

max_issues_repo_path: Interactive Story Telling.jl

max_forks_repo_path: Interactive Story Telling.jl


text:
module TestGraphGenerators
using Test

using Catlab.Graphs.BasicGraphs, Catlab.Graphs.GraphGenerators
using Catlab.CategoricalAlgebra: homomorphisms

# Path graphs
#------------

n = 5
g = path_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n-1)
g = path_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))

# Cycle graphs
#-------------

g = cycle_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n)
g = cycle_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2n)

# Complete graphs
#----------------

for T in (Graph, SymmetricGraph)
  g = complete_graph(T, n)
  @test (nv(g), ne(g)) == (n, n*(n-1))
end
for T in (ReflexiveGraph, SymmetricReflexiveGraph)
  g = complete_graph(T, n)
  @test (nv(g), ne(g)) == (n, n*n)
end

# Star graphs
#------------

g = star_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n-1)
@test length(unique(src(g,e) for e in edges(g))) == 1
g = star_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))

# Wheel graphs
#-------------

g = wheel_graph(Graph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))
triangle = Graph(3)
add_edges!(triangle, [1,2,1], [2,3,3])
@test length(homomorphisms(triangle, g)) == n-1

g = wheel_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 4(n-1))
triangle = cycle_graph(SymmetricGraph, 3)
@test length(homomorphisms(triangle, g)) == 6(n-1) # == 3! * (n-1)

end


text:
using GenderInference
import GenderInference: RawDataSet, NameDataset, parsedataset
using Test

@testset "Data Handling" begin
    @test typeof(GenderInference.NAMES) <: NameDataset

    unknown_raw = RawDataSet(:test)
    @test unknown_raw == RawDataSet{:test}()
    @test unknown_raw == RawDataSet("test")

    @test_throws ArgumentError parsedataset("a/path", unknown_raw)

end

@testset "Gender Inference" begin

    @test gender("Kevin") == :male
    @test gender("rachel") == :female
    @test gender("foo") === missing
    @test gender("casey") == :male
    @test gender("casey", threshold=0.8) === missing

    @test typeof(gendercount("stefan")) <: NamedTuple{(:female,:male), Tuple{T,T}} where T <: Integer
    @test typeof(gendercount("bar")) <: NamedTuple{(:female,:male), Tuple{T,T}} where T <: Integer

    @test typeof(proportionfemale("Viral")) <: AbstractFloat
    @test typeof(proportionmale("Jane")) <: AbstractFloat
    @test proportionfemale("Foo") === missing
    @test proportionmale("Gadz") === missing

    @test gendercount("Julia", 1850:1900) === (female=missing, male=missing)
end


text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport Compute_Jacobian_MA_journal
@pyimport GLS_Apr_weekday_PM_journal

numNodes = Compute_Jacobian_MA_journal.numNodes;
numLinks = Compute_Jacobian_MA_journal.numLinks;
numODpairs = Compute_Jacobian_MA_journal.numODpairs;
numZones = Compute_Jacobian_MA_journal.numZones
od_pairs = Compute_Jacobian_MA_journal.od_pairs;
link_list_js = Compute_Jacobian_MA_journal.link_list_js;
link_length_list = Compute_Jacobian_MA_journal.link_length_list;

flow_observ = GLS_Apr_weekday_PM_journal.x_

include("../Julia_files/initia_data.jl");
include("prepare_data_MA_journal.jl");
include("extract_data_MA_journal.jl");
include("../Julia_files/tap_MSA.jl");
include("../Julia_files/inverseVI.jl");
include("../Julia_files/demands_adjustment_gradi.jl");

# demandsDict_ = readall("../08_develop_new_OD_demand_estimator_MA_Dijkstra_uni_class_Apr_weekend_journal/results/demandsDict.json");
# demandsDict_ = JSON.parse(demandsDict_);

# demandsDict_

# demandsDict__ = demandsDict_["24"]

# _demandsDict__ = Dict()

# for key_ in keys(demandsDict__)
#     key = (int(split(split(key_, ",")[1], "(")[2]),int(split(split(key_, ",")[2], ")")[1]))
#     _demandsDict__[key] = demandsDict__[key_]
# end

# _demandsDict__

# demandsDict[0] = _demandsDict__;
# demandsDict[1] = _demandsDict__;

demandsDict[0] == demandsDict[1]

demandsDiffDict[1] = norm(demandsDicToVec(demandsDict[1]) - demandsDicToVec(demandsDict[0]))/
                     norm(demandsDicToVec(demandsDict[0]));

demandsDiffDict[1]

function demandsDictFixed(day)
#     day = 4  # day of April

    # observed flow vector
    xl = flow_observ[:, day]

    tapFlows = Dict{}()

    for i = 1:length(ta_data.start_node)
        key = (ta_data.start_node[i], ta_data.end_node[i])
        tapFlows[key] = xl[i]
    end

    tapFlowVect = xl;

    # get observed flow vector (corresponding to ground truth demands and ground truth costs)
    tapFlowDicDict[0], tapFlowVecDict[0] = tapFlows, tapFlowVect;

    # get arcs data corresponding to ground truth demands and flows
    arcsDict[0] = observFlow("./benchmark_data/journal_net_Apr_weekend.txt", tapFlowDicDict[0]);

    arcsDict[0]

    deg = 6
    c = 0.5
    lam = 10000.0

    coeffs_dict_Apr_weekend_ = readall("../temp_files/coeffs_dict_Apr_weekend.json")
    coeffs_dict_Apr_weekend_ = JSON.parse(coeffs_dict_Apr_weekend_)
    fcoeffs = coeffs_dict_Apr_weekend_["($(deg),$(c),$(lam),1)"]

    # fcoeffs = [1, 0, 0, .15]

    demandsVecDict[1] = demandsDicToVec(demandsDict[1]);
    objFunDict[1] = objF(demandsVecDict[1], fcoeffs);

    # get initial flow vector (corresponding to initial demands)
    tapFlowDicDict[1], tapFlowVecDict[1] = tapMSA(demandsDict[1], fcoeffs);

    demandsVecDict[0] = demandsDicToVec(demandsDict[0]);

    # get arcs data corresponding to initial demands and flows
    arcsDict[1] = observFlow("./benchmark_data/journal_net_Apr_weekend.txt", tapFlowDicDict[1]);

    linkCostDicDict[1] = tapFlowVecToLinkCostDict(tapFlowVecDict[1], fcoeffs);

    linkCostDicDict[1]["0"], link_length_list[1]

    jacobiSpiessDict[1] = Compute_Jacobian_MA_journal.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                                  link_list_js, [linkCostDicDict[1]["$(i)"] for i=0:numLinks-1]);

    # maximum number of iterations
    N = 100;

    # Armijo rule parameters
    rho = 2;
    M = 10;

    # search direction parameter
    epsilon_1 = 0;

    # stop criterion parameter
    epsilon_2 = 1e-20;

    for l = 1:N

        jacobDict[l] = jacobiSpiessDict[l]

        descDirecDict[l] = descDirec(tapFlowVecDict[l], tapFlowVecDict[0], jacobDict[l], numODpairs, numLinks);

        demandsVecDict[l] = demandsDicToVec(demandsDict[l]);

        searchDirecDict[l] = searchDirec(demandsVecDict[l], descDirecDict[l], epsilon_1);

        thetaMaxDict[l] = thetaMax(demandsVecDict[l], searchDirecDict[l]);

        demandsVecDict[l+1] = similar(demandsVecDict[0]);

        demandsVecDict[l+1], objFunDict[l+1] = armijo(objFunDict[l], demandsVecDict[l], fcoeffs, searchDirecDict[l],
        thetaMaxDict[l], rho, M);

        demandsDict[l+1] = demandsVecToDic(demandsVecDict[l+1]);

        tapFlowDicDict[l+1], tapFlowVecDict[l+1] = tapMSA(demandsDict[l+1], fcoeffs);

        arcsDict[l+1] = observFlow("./benchmark_data/journal_net_Apr_weekend.txt", tapFlowDicDict[l+1]);

        linkCostDicDict[l+1] = tapFlowVecToLinkCostDict(tapFlowVecDict[l+1], fcoeffs);

        jacobiSpiessDict[l+1] = Compute_Jacobian_MA_journal.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                                  link_list_js, [linkCostDicDict[l+1]["$(i)"] for i=0:numLinks-1]);

        demandsDiffDict[l+1] = norm(demandsVecDict[l+1] - demandsVecDict[0]) / norm(demandsVecDict[0]);

        # stopping criterion
        if (objFunDict[l] - objFunDict[l+1]) / objFunDict[1] < epsilon_2
            break
        end

        println("iteration $(l) finished...")

    end

    # normalize objective function value
    for l = 1:(length(objFunDict))
        norObjFunDict[l] = objFunDict[l] / objFunDict[1];
    end

    outfile = open("./results/demandsVecDict$(day)_journal.json", "w")

    JSON.print(outfile, demandsVecDict)

    close(outfile)

    outfile = open("./results/demandsDict$(day)_journal.json", "w")

    JSON.print(outfile, demandsDict)

    close(outfile)

    outfile = open("./results/tapFlowDicDict$(day)_journal.json", "w")

    JSON.print(outfile, tapFlowDicDict)

    close(outfile)

    outfile = open("./results/tapFlowVecDict$(day)_journal.json", "w")

    JSON.print(outfile, tapFlowVecDict)

    close(outfile)

    demandsDict[length(demandsDict)-1]

    demandsDict_ = Dict{}()

    for key in keys(demandsDict[length(demandsDict)-1])
        demandsDict_[key] = demandsDict[length(demandsDict)-1][key]
    end

    outfile = open("./results/demandsDictFixed$(day)_journal.json", "w")

    JSON.print(outfile, demandsDict_)

    close(outfile)
end

weekend_Apr_list = [1, 7, 8, 14, 15, 21, 22, 28, 29]

# weekend_Apr_list = [6]

for day in weekend_Apr_list
    demandsDictFixed(day)
    println("day $(day) finished...")
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
using SafeTestsets

@safetestset "My_f tests " begin include("my_f_tests.jl") end


text:
# Three-state Bioreactor
id = :PerezGalvan2017e

pL = [6.45; 0.46; 1.05]
pU = [6.55; 0.47; 1.1]
PerezGalvan2017e_x0(p::Vector{T}) where T = T[p[1]; 5.0*one(T); 15.0*one(T)]
function PerezGalvan2017e_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    x3m = 50.0; x2f = 20.0; D = 0.202; Y = 0.4; β = 0.2; α = 0.5
    μ = p[2]*(1.0 - u[3]/x3m)*u[2]/(p[3] + u[2])
    du[1] = (μ - D)*u[1]
    du[2] = D*(x2f - u[2]) - μ*u[1]/Y
    du[3] = D*u[3] + (α*μ + β)*u[1]
    return
end
tspan = (0.0, 7.7)

prob = ODERelaxProb(PerezGalvan2017e_f!, tspan, PerezGalvan2017e_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0098135417300923"
source = "PerezGalvan2017"
desc = "Three-state Bioreactor"


text:
# https://github.com/JuliaLang/julia/blob/b166d0d15f616f4a025ac5905a115399cc932ddc/base/compiler/abstractinterpretation.jl#L226-L314

function abstract_call_method_with_const_args(interp::AbstractInterpreter, @nospecialize(rettype), @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch, sv::InferenceState, edgecycle::Bool)
    method = match.method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    length(argtypes) >= nargs || return Any
    haveconst = false
    allconst = true
    # see if any or all of the arguments are constant and propagating constants may be worthwhile
    for a in argtypes
        a = widenconditional(a)
        if allconst && !isa(a, Const) && !isconstType(a) && !isa(a, PartialStruct)
            allconst = false
        end
        if !haveconst && has_nontrivial_const_info(a) && const_prop_profitable(a)
            haveconst = true
        end
        if haveconst && !allconst
            break
        end
    end
    haveconst || improvable_via_constant_propagation(rettype) || return Any
    if nargs > 1
        if istopfunction(f, :getindex) || istopfunction(f, :setindex!)
            arrty = argtypes[2]
            # don't propagate constant index into indexing of non-constant array
            if arrty isa Type && arrty <: AbstractArray && !issingletontype(arrty)
                return Any
            elseif arrty ⊑ Array
                return Any
            end
        elseif istopfunction(f, :iterate)
            itrty = argtypes[2]
            if itrty ⊑ Array
                return Any
            end
        end
    end
    if !allconst && (istopfunction(f, :+) || istopfunction(f, :-) || istopfunction(f, :*) ||
                     istopfunction(f, :(==)) || istopfunction(f, :!=) ||
                     istopfunction(f, :<=) || istopfunction(f, :>=) || istopfunction(f, :<) || istopfunction(f, :>) ||
                     istopfunction(f, :<<) || istopfunction(f, :>>))
        return Any
    end
    force_inference = allconst || InferenceParams(interp).aggressive_constant_propagation
    if istopfunction(f, :getproperty) || istopfunction(f, :setproperty!)
        force_inference = true
    end
    mi = specialize_method(match, !force_inference)
    mi === nothing && return Any
    mi = mi::MethodInstance
    # decide if it's likely to be worthwhile
    if !force_inference && !const_prop_heuristic(interp, method, mi)
        return Any
    end
    inf_cache = get_inference_cache(interp)
    inf_result = cache_lookup(mi, argtypes, inf_cache)
    if inf_result === nothing
        if edgecycle
            # if there might be a cycle, check to make sure we don't end up
            # calling ourselves here.
            infstate = sv
            cyclei = 0
            while !(infstate === nothing)
                if method === infstate.linfo.def && any(infstate.result.overridden_by_const)
                    return Any
                end
                if cyclei < length(infstate.callers_in_cycle)
                    cyclei += 1
                    infstate = infstate.callers_in_cycle[cyclei]
                else
                    cyclei = 0
                    infstate = infstate.parent
                end
            end
        end
        inf_result = InferenceResult(mi, argtypes)
        frame = InferenceState(inf_result, #=cache=#false, interp)
        frame === nothing && return Any # this is probably a bad generated function (unsound), but just ignore it
        frame.limited = true
        frame.parent = sv
        push!(inf_cache, inf_result)
        typeinf(interp, frame) || return Any
    end
    result = inf_result.result
    # if constant inference hits a cycle, just bail out
    isa(result, InferenceState) && return Any
    add_backedge!(inf_result.linfo, sv)
    return result
end


text:
using InstantiateFromURL
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.7.0")
# uncomment to force package installation and precompilation
# github_project("QuantEcon/quantecon-notebooks-julia", version="0.7.0", instantiate=true, precompile = true)

max_stars_repo_name: artemsolod/lecture-source-jl

max_issues_repo_name: artemsolod/lecture-source-jl

max_forks_repo_name: artemsolod/lecture-source-jl


text:
function sim_logP!(gpT::GPE, gpC::GPE, gpNull::GPE, treat::BitVector)
    Ysim = prior_rand(gpNull)

    gpT.y = Ysim[treat]
    gpC.y = Ysim[.!treat]
    gpNull.y = Ysim

    update_alpha!(gpT)
    update_alpha!(gpC)
    update_alpha!(gpNull)

    gpT.mll = mll(gpT)
    gpC.mll = mll(gpC)
    gpNull.mll = mll(gpNull)

    mll_altv = gpT.mll + gpC.mll
    mll_null = gpNull.mll
    return mll_null, mll_altv
end

function nsim_logP(gpT::GPE, gpC::GPE, gpNull::GPE, nsim::Int)
    gpT_mod = modifiable(gpT)
    gpC_mod = modifiable(gpC)
    treat = BitVector(gpNull.nobs)
    treat[:] = false
    treat[1:gpT.nobs] = true
    mll_sims = [sim_logP!(gpT_mod, gpC_mod, gpNull, treat)
        for _ in 1:nsim];
    return mll_sims
end

function boot_mlltest(gpT::GPE, gpC::GPE, nsim::Int)
    mll_alt = gpT.mll + gpC.mll

    gpNull = make_null(gpT, gpC)
    mll_null = gpNull.mll
    mll_sims = GeoRDD.nsim_logP(gpT, gpC, gpNull, nsim)

    mll_sim_null = [sim[1] for sim in mll_sims]
    mll_sim_altv = [sim[2] for sim in mll_sims]

    Δmll_obs = mll_alt - mll_null
    Δmll_sim = mll_sim_altv .- mll_sim_null
    return mean(Δmll_sim .> Δmll_obs)
end

function placebo_mll(angle::Float64, X::Matrix, Y::Vector,
                 kern::Kernel, m::Mean, logNoise::Float64,
                 nsim::Int)
    shift = shift_for_even_split(angle, X)
    left = left_points(angle, shift, X)
    gp_left  = GPE(X[:,left],  Y[left],  m, kern, logNoise)
    gp_right = GPE(X[:,.!left], Y[.!left], m, kern, logNoise)
    pval = boot_mlltest(gp_left, gp_right, nsim)
    return pval
end
function placebo_mll(angle::Float64, gp::GPE, nsim::Int; kwargs...)
    return placebo_mll(angle, gp.x, gp.y, gp.kernel, gp.m, gp.logNoise, nsim; kwargs...)
end


text:
using Flux, Flux.Data.MNIST, Statistics, BSON, Random, StatsBase, Plots; pyplot()
using Flux: params, onehotbatch, crossentropy, update!
Random.seed!(0)

nTrain = 20000
miniBatchSize = 1000
imgs   = Flux.Data.MNIST.images()[1:nTrain]
labels = Flux.Data.MNIST.labels()[1:nTrain]

trainData = hcat([vcat(float.(imgs[i])...) for i in 1:nTrain]...)
trainLabels = labels

testImgs = Flux.Data.MNIST.images(:test)
testLabels = Flux.Data.MNIST.labels(:test)
nTest = length(testImgs)
testData = hcat([vcat(float.(testImgs[i])...) for i in 1:nTest]...)

W = randn(10,28*28)
b = randn(10)

logisticM(imgVec) = softmax(W*imgVec .+ b)
logisticMclassifier(imgVec) = argmax(logisticM(imgVec))-1
loss(x,y) = crossentropy(logisticM(x),onehotbatch(y,0:9))
opt = ADAM(0.01)

lossValue = 0.0
lossArray = []
epochNum = 0
while true
    global lossValue
    prevLossValue = lossValue
    for batch in Iterators.partition(1:nTrain,miniBatchSize)
      gs = gradient(()->loss(trainData[:,batch],trainLabels[batch]),params(W,b))
      for p in (W,b)
          update!(opt,p,gs[p])
      end
    end
    global epochNum += 1
    lossValue = loss(trainData,trainLabels)
    push!(lossArray,lossValue)
    print(".")
    abs(prevLossValue-lossValue) < 5e-4 && break
end

println("\nNumber of epochs: ", epochNum)
acccuracy = mean([logisticMclassifier(testData[:,k]) for k in 1:nTest]
                .== testLabels)
println("Accuracy: ", acccuracy)

plot(lossArray, xlabel = "Epoch", ylabel = "Cross Entropy Loss", legend = false)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using FeatureDescriptors
using FeatureDescriptors.TestUtils
using Test

# used for testing
abstract type NaiveFakeDescriptor <: FakeDescriptor end
abstract type SpecialFakeDescriptor <: FakeDescriptor end

# used to test isless
abstract type FakeDescriptorA <: FakeDescriptor end
abstract type FakeDescriptorB <: FakeDescriptor end
abstract type FakeDescriptorC{T} <: FakeDescriptor end

@testset "FeatureDescriptors.jl" begin

    @testset "FakeDescriptor" begin

        test_interface(FakeDescriptor)

        @testset "naive subtyping" begin
            # Inherits the supertype behaviour
            test_interface(NaiveFakeDescriptor)
            @test sources(NaiveFakeDescriptor) == ["fake_table"]
            @test quantity_key(NaiveFakeDescriptor) == :quantity
            @test categorical_keys(NaiveFakeDescriptor) == [:category1, :category2]
            @test label(NaiveFakeDescriptor) == :naivefakedescriptor
        end

        @testset "special subtyping" begin
            # Define different quantitative variable and no categories
            FeatureDescriptors.quantity_key(::Type{<:SpecialFakeDescriptor}) = :foo
            FeatureDescriptors.categorical_keys(::Type{<:SpecialFakeDescriptor}) = []

            test_interface(SpecialFakeDescriptor)
            @test sources(SpecialFakeDescriptor) == ["fake_table"]
            @test quantity_key(SpecialFakeDescriptor) == :foo
            @test isempty(categorical_keys(SpecialFakeDescriptor))
            @test label(SpecialFakeDescriptor) == :specialfakedescriptor
        end

    end

    @testset "isless" begin
        @test isless(FakeDescriptorA, FakeDescriptorB)
        @test isless(FakeDescriptorA, FakeDescriptorC{Float64})
        @test !isless(FakeDescriptorC, FakeDescriptorB)
        unordered = [FakeDescriptorC{Int64}, FakeDescriptorB, FakeDescriptorC{Float64}, FakeDescriptorA]
        ordered = [FakeDescriptorA, FakeDescriptorB, FakeDescriptorC{Float64}, FakeDescriptorC{Int64}]
        @test sort(unordered) == ordered
    end

end

max_stars_repo_name: invenia/FeatureDescriptors

max_issues_repo_name: invenia/FeatureDescriptors

max_forks_repo_name: invenia/FeatureDescriptors


text:
"""

```
Residue(atoms::AbstractVector{Atom}, range::UnitRange{Int})
```

Residue data structure. It contains two fields: `atoms` which is a vector of
`Atom` elements, and `range`, which indicates which atoms of the `atoms` vector
compose the residue.

The Residue structure carries the properties of the residue or molecule of the atoms
it contains, but it does not copy the original vector of atoms, only the residue
meta data for each residue.

### Example

```julia-repl
julia> pdb = wget("1LBD");

julia> residues = collect(eachresidue(pdb))
   Array{Residue,1} with 238 residues.

julia> resnum.(residues[1:3])
3-element Vector{Int64}:
 225
 226
 227

julia> residues[5].chain
"A"

julia> residues[8].range
52:58

```

"""
struct Residue{T<:AbstractVector{Atom}}
  atoms::T
  range::UnitRange{Int}
  name::String
  resname::String
  residue::Int
  resnum::Int
  chain::String
  model::Int
  segname::String
end
name(residue::Residue) = residue.name
resname(residue::Residue) = residue.resname
residue(residue::Residue) = residue.residue
resnum(residue::Residue) = residue.resnum
chain(residue::Residue) = residue.chain
model(residue::Residue) = residue.model
segname(residue::Residue) = residue.segname
mass(residue::Residue) = mass(@view residue.atoms[residue.range])

function Residue(atoms::AbstractVector{Atom}, range::UnitRange{Int})
  i = range[begin]
  # Check if the range effectivelly corresponds to a single residue (unsafe check)
  for j in range[begin]+1:range[end]
    if atoms[j].residue != atoms[i].residue
      error("Range $range does not correspond to a single residue or molecule.")
    end
  end
  Residue(atoms,range,
          atoms[i].resname, atoms[i].resname, atoms[i].residue,
          atoms[i].resnum, atoms[i].chain, atoms[i].model, atoms[i].segname)
end
Residue(atoms::AbstractVector{Atom}) = Residue(atoms,1:length(atoms))

function Base.getindex(residue::Residue,i)
  @assert i > 0 "Index must be in 1:$(residue.range[end]-residue.range[begin])"
  @assert (i <= length(residue.range)) "Residue has $(residue.range[end]-residue.range[begin]+1) atoms."
  i = residue.range[begin] + i - 1
  residue.atoms[i]
end

#
# Structure and function to define the eachresidue iterator
#
struct EachResidue{T<:AbstractVector{Atom}}
  atoms::T
end

"""

```
eachresidue(atoms::AbstractVector{Atom})
```

Iterator for the residues (or molecules) of a selection.

### Example

```julia-repl
julia> atoms = wget("1LBD");

julia> length(eachresidue(atoms))
238

julia> for res in eachresidue(atoms)
         println(res)
       end
 Residue of name SER with 6 atoms.
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       1    N     SER     A      225        1   45.228   84.358   70.638 67.05  1.00     1       -         1
       2   CA     SER     A      225        1   46.080   83.165   70.327 68.73  1.00     1       -         2
       3    C     SER     A      225        1   45.257   81.872   70.236 67.90  1.00     1       -         3
       4    O     SER     A      225        1   45.823   80.796   69.974 64.85  1.00     1       -         4
       5   CB     SER     A      225        1   47.147   82.980   71.413 70.79  1.00     1       -         5
       6   OG     SER     A      225        1   46.541   82.639   72.662 73.55  1.00     1       -         6

 Residue of name ALA with 5 atoms.
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       7    N     ALA     A      226        2   43.940   81.982   70.474 67.09  1.00     1       -         7
       8   CA     ALA     A      226        2   43.020   80.825   70.455 63.69  1.00     1       -         8
       9    C     ALA     A      226        2   41.996   80.878   69.340 59.69  1.00     1       -         9
                                                      ...

```

"""
eachresidue(atoms::AbstractVector{Atom}) = EachResidue(atoms)

# Collect residues default constructor
Base.collect(r::EachResidue) = collect(Residue,r)

#
# Iterate over the resiudes
#
function Base.iterate(residues::EachResidue, state=1)
  r0 = state
  r0 > length(residues.atoms) && return nothing
  residue0 = residues.atoms[r0].residue
  r1 = r0
  while r1 <= length(residues.atoms)
    if residues.atoms[r1].residue != residue0
      return (Residue(residues.atoms,r0:r1-1),r1)
    end
    r1 += 1
  end
  return (Residue(residues.atoms,r0:r1-1),r1)
end

#
# Iterate over atoms of one residue
#
function Base.iterate(residue::Residue, state=1)
  i1 = residue.range[begin] + state - 1
  if i1 <= residue.range[end]
    return (residue.atoms[i1],state+1)
  else
    return nothing
  end
end

#
# Length of the Residue struct and eachresidue iterator (number of residues)
#
Base.length(residues::EachResidue) = sum( 1 for residue in residues )
Base.length(residue::Residue) = length(residue.range)

#
# io show functions
#
function Base.show(io::IO, residue::Residue)
  natoms = residue.range[end]-residue.range[begin]+1
  println(" Residue of name $(name(residue)) with $natoms atoms.")
  print_short_atom_list(@view residue.atoms[residue.range])
end

function Base.show(io::IO, residues::EachResidue)
  println(" Iterator with $(length(residues)) residues.")
end

function Base.show(io::IO,::MIME"text/plain", residues::AbstractVector{Residue} )
  println("   Array{Residue,1} with $(length(residues)) residues.")
end


text:
using FunctionalStateMachine
using Graphs
using Dates
using Test

## User state functions

function bar!(usrdata)
  println("do bar!")
  return FunctionalStateMachine.exitStateMachine
end

function foo!(usrdata)
  println("do foo!")
  return bar!
end

@testset "Test generic state machine..." begin

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, verbose=true); end

@test statemachine.next == emptyState

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, verbose=false); end

@test statemachine.next == emptyState

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, breakafter=foo!); end

@test statemachine.next == bar!

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, iterlimit=1); end

@test statemachine.next == bar!

statemachine = StateMachine{Nothing}(next=bar!)
while statemachine(nothing, verbose=true); end

@test statemachine.next == emptyState

end

@testset "test watchdog timeout" begin

function longwait(x)
  @info "starting stagnant function call, but should not see it's end (but watchdog timeout)"
  while true
    print(".")
    sleep(0.5)
  end
  @info "done with longwait"
  return exitStateMachine
end

statemachine = StateMachine{Nothing}(next=longwait)
try
  while statemachine(nothing, 2.0, verbose=true); end
catch e
  @info " watchdog test, successfully caught exception for stagnant FSM step"
  @test_throws InterruptException throw(e)
end

end

@testset "test recording and rendering of an FSM run" begin

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, recordhistory=true); end

hists = Dict{Symbol,Vector{Tuple{DateTime,Int,Function,Nothing}}}(:first => statemachine.history)

animateStateMachineHistoryIntervalCompound(hists, interval=1)

end

#