text:
__precompile__()
module PackedASCIITypes

export PackedASCII, @pasc_str

using Base: StringVector
import Base:
    length, ==, hash, <, typemin, show, String, convert, promote_rule, write

abstract type PackedASCII end

primitive type PackedASCII1 <: PackedASCII 8 end
primitive type PackedASCII2 <: PackedASCII 16 end
primitive type PackedASCII4 <: PackedASCII 32 end
primitive type PackedASCII8 <: PackedASCII 64 end
primitive type PackedASCII17 <: PackedASCII 128 end

"""
    maxlength(::Type{<:PackedASCII})

Maximum number of characters allowed by a type.
"""
function maxlength end

for (T, S, maxlen, lenbits) in [(PackedASCII1, UInt8, 1, 1),
                                (PackedASCII2, UInt16, 2, 2),
                                (PackedASCII4, UInt32, 4, 3),
                                (PackedASCII8, UInt64, 8, 4),
                                (PackedASCII17, UInt128, 17, 5)]
    # sanity checks and useful quantities
    @assert maxlen+1 ≤ 2^lenbits
    nbits = sizeof(S)*8
    @assert nbits ≥ maxlen*7 + lenbits
    firstoffset = lenbits + (maxlen - 1) * 7 # first char starts here
    lenmask = S(2^lenbits - 1)              # extract length bits

    @eval length(p::$T) = Int(reinterpret($S, p) & $lenmask)

    # reflection
    @eval maxlength(::Type{$T}) = $maxlen

    # compose and decompose
    #
    # The decomposed representation is the character bits shifted to the right,
    # and the length. These functions are used internally.
    @eval function _compose(::Type{$T}, s::$S, len::Integer)
        u = s << ($lenbits + 7 * ($maxlen - len))
        u |= $S(len)
        reinterpret($T, u)
    end
    @eval _compose(::Type{$T}, s, len) = _compose($T, convert($S, s), len)

    @eval function _decompose(p::$T)
        u = reinterpret($S, p)
        len = u & $lenmask
        s = u >> ($lenbits + 7*($maxlen - len))
        s, len
    end

    # conversion to and from String
    Tname = nameof(T)
    @eval export $(Tname)

    @eval function ($Tname)(str::AbstractString)
        len = length(str)
        (len ≤ $maxlen && isascii(str)) || throw(InexactError(Symbol($(Tname)), $(T), str))
        s = zero($S)
        for c in str
            s <<= 7
            s |= $S(c)
        end
        _compose($T, s, len)
    end
    @eval ($Tname)(p::$T) = p
    @eval ($Tname)(p::PackedASCII) = convert($T, p)
    @eval function String(p::$T)
        len = length(p)
        sv = StringVector(len)
        u = reinterpret($S, p)
        for i in 1:len
            sv[i] = UInt8((u >> $firstoffset) & 0x7f)
            u <<= 7
        end
        String(sv)
    end

    # read macro
    readT = Symbol("pasc$(maxlen)_str")
    readTmacro = Symbol("@" * string(readT))
    @eval export $(readTmacro)
    @eval macro ($readT)(str)
        $T(str)
    end

    # methods
    @eval (==)(p::$T, q::$T) = reinterpret($S, p) == reinterpret($S, q)
    @eval function hash(p::$T, h::UInt)
        s, len = _decompose(p)
        hash(s, hash(len, h))
    end
    @eval (<)(p::$T, q::$T) = reinterpret($S, p) < reinterpret($S, q)
    @eval typemin(::Type{$T}) = reinterpret($T, zero($S))
    @eval show(io::IO, p::$T) = print(io, "pasc$($(maxlen))\"", String(p), "\"")
    @eval write(io::IO, p::$T) = write(io, reinterpret($S, p))
end

function convert(::Type{T}, p::S) where {T <: PackedASCII, S <: PackedASCII}
    shiftedbits, len = _decompose(p)
    _compose(T, shiftedbits, len)
end

function (==)(p::PackedASCII, q::PackedASCII)
    sa, lena = _decompose(p)
    sb, lenb = _decompose(q)
    lena == lenb && sa == sb
end

"""
    PackedASCII(str)

Construct the shortest `PackedASCII` representation that can fit the string `str`.

The maximum length is 17 characters, and the string has to be ASCII.

See the `README.md` of the package for details.
"""
function PackedASCII(str::AbstractString)
    len = length(str)
    if len ≤ 1               # FIXME ugly implementation, worth code generation?
        PackedASCII1(str)
    elseif len ≤ 2
        PackedASCII2(str)
    elseif len ≤ 4
        PackedASCII4(str)
    elseif len ≤ 8
        PackedASCII8(str)
    elseif len ≤ 17
        PackedASCII17(str)
    else
        throw(InexactError(:PackedASCII, PackedASCII, str))
    end
end

macro pasc_str(str)
    PackedASCII(str)
end

promote_rule(::Type{T}, ::Type{S}) where {T <: PackedASCII, S <: PackedASCII} =
    maxlength(T) ≥ maxlength(S) ? T : S

(<)(a::PackedASCII, b::PackedASCII) = (<)(promote(a, b)...)

end # module


text:
using LinearAlgebra
using Oceananigans.Architectures: array_type
using Oceananigans.TimeSteppers: _compute_w_from_continuity!

function can_solve_single_tridiagonal_system(arch, N)
    ArrayType = array_type(arch)

    a = rand(N-1)
    b = 3 .+ rand(N) # +3 to ensure diagonal dominance.
    c = rand(N-1)
    f = rand(N)

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    M = Tridiagonal(a, b, c)
    ϕ_correct = M \ f

    # Convert to CuArray if needed.
    a, b, c, f = ArrayType.([a, b, c, f])

    ϕ = reshape(zeros(N), (1, 1, N)) |> ArrayType

    grid = RegularCartesianGrid(size=(1, 1, N), length=(1, 1, 1))
    btsolver = BatchedTridiagonalSolver(arch; dl=a, d=b, du=c, f=f, grid=grid)

    solve_batched_tridiagonal_system!(ϕ, arch, btsolver)

    return Array(ϕ[:]) ≈ ϕ_correct
end

function can_solve_single_tridiagonal_system_with_functions(arch, N)
    ArrayType = array_type(arch)

    grid = RegularCartesianGrid(size=(1, 1, N), length=(1, 1, 1))

    a = rand(N-1)
    c = rand(N-1)

    @inline b(i, j, k, grid, p) = 3 .+ cos(2π*grid.zC[k])  # +3 to ensure diagonal dominance.
    @inline f(i, j, k, grid, p) = sin(2π*grid.zC[k])

    bₐ = [b(1, 1, k, grid, nothing) for k in 1:N]
    fₐ = [f(1, 1, k, grid, nothing) for k in 1:N]

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    M = Tridiagonal(a, bₐ, c)
    ϕ_correct = M \ fₐ

    # Convert to CuArray if needed.
    a, c = ArrayType.([a, c])

    ϕ = reshape(zeros(N), (1, 1, N)) |> ArrayType

    btsolver = BatchedTridiagonalSolver(arch; dl=a, d=b, du=c, f=f, grid=grid)

    solve_batched_tridiagonal_system!(ϕ, arch, btsolver)

    return Array(ϕ[:]) ≈ ϕ_correct
end

function can_solve_batched_tridiagonal_system_with_3D_RHS(arch, Nx, Ny, Nz)
    ArrayType = array_type(arch)

    a = rand(Nz-1)
    b = 3 .+ rand(Nz) # +3 to ensure diagonal dominance.
    c = rand(Nz-1)
    f = rand(Nx, Ny, Nz)

    M = Tridiagonal(a, b, c)
    ϕ_correct = zeros(Nx, Ny, Nz)

    # Solve the systems with backslash on the CPU to avoid scalar operations on the GPU.
    for i = 1:Nx, j = 1:Ny
        ϕ_correct[i, j, :] .= M \ f[i, j, :]
    end

    # Convert to CuArray if needed.
    a, b, c, f = ArrayType.([a, b, c, f])

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(1, 1, 1))
    btsolver = BatchedTridiagonalSolver(arch; dl=a, d=b, du=c, f=f, grid=grid)

    ϕ = zeros(Nx, Ny, Nz) |> ArrayType

    solve_batched_tridiagonal_system!(ϕ, arch, btsolver)

    return Array(ϕ) ≈ ϕ_correct
end

function can_solve_batched_tridiagonal_system_with_3D_functions(arch, Nx, Ny, Nz)
    ArrayType = array_type(arch)

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(1, 1, 1))

    a = rand(Nz-1)
    c = rand(Nz-1)

    @inline b(i, j, k, grid, p) = 3 + grid.xC[i]*grid.yC[j] * cos(2π*grid.zC[k])
    @inline f(i, j, k, grid, p) = (grid.xC[i] + grid.yC[j]) * sin(2π*grid.zC[k])

    ϕ_correct = zeros(Nx, Ny, Nz)

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    for i = 1:Nx, j = 1:Ny
        bₐ = [b(i, j, k, grid, nothing) for k in 1:Nz]
        M = Tridiagonal(a, bₐ, c)

        fₐ = [f(i, j, k, grid, nothing) for k in 1:Nz]
        ϕ_correct[i, j, :] .= M \ fₐ
    end

    # Convert to CuArray if needed.
    a, c = ArrayType.([a, c])

    btsolver = BatchedTridiagonalSolver(arch; dl=a, d=b, du=c, f=f, grid=grid)

    ϕ = zeros(Nx, Ny, Nz) |> ArrayType
    solve_batched_tridiagonal_system!(ϕ, arch, btsolver)

    return Array(ϕ) ≈ ϕ_correct
end

function vertically_stretched_poisson_solver_correct_answer(arch, Nx, Ny, zF)
    Lx, Ly, Lz = 1, 1, zF[end]
    Δx, Δy = Lx/Nx, Ly/Ny

    #####
    ##### Vertically stretched operators
    #####

    @inline δx_caa(i, j, k, f) = @inbounds f[i+1, j, k] - f[i, j, k]
    @inline δy_aca(i, j, k, f) = @inbounds f[i, j+1, k] - f[i, j, k]
    @inline δz_aac(i, j, k, f) = @inbounds f[i, j, k+1] - f[i, j, k]

    @inline ∂x_caa(i, j, k, Δx,  f) = δx_caa(i, j, k, f) / Δx
    @inline ∂y_aca(i, j, k, Δy,  f) = δy_aca(i, j, k, f) / Δy
    @inline ∂z_aac(i, j, k, ΔzF, f) = δz_aac(i, j, k, f) / ΔzF[k]

    @inline ∂x²(i, j, k, Δx, f)       = (∂x_caa(i, j, k, Δx, f)  - ∂x_caa(i-1, j, k, Δx, f))  / Δx
    @inline ∂y²(i, j, k, Δy, f)       = (∂y_aca(i, j, k, Δy, f)  - ∂y_aca(i, j-1, k, Δy, f))  / Δy
    @inline ∂z²(i, j, k, ΔzF, ΔzC, f) = (∂z_aac(i, j, k, ΔzF, f) - ∂z_aac(i, j, k-1, ΔzF, f)) / ΔzC[k]

    @inline div_ccc(i, j, k, Δx, Δy, ΔzF, u, v, w) = ∂x_caa(i, j, k, Δx, u) + ∂y_aca(i, j, k, Δy, v) + ∂z_aac(i, j, k, ΔzF, w)

    @inline ∇²(i, j, k, Δx, Δy, ΔzF, ΔzC, f) = ∂x²(i, j, k, Δx, f) + ∂y²(i, j, k, Δy, f) + ∂z²(i, j, k, ΔzF, ΔzC, f)

    #####
    ##### Generate "fake" vertically stretched grid
    #####

    function grid(zF)
        Nz = length(zF) - 1
        ΔzF = [zF[k+1] - zF[k] for k in 1:Nz]
        zC = [(zF[k] + zF[k+1]) / 2 for k in 1:Nz]
        ΔzC = [zC[k+1] - zC[k] for k in 1:Nz-1]
        return zF, zC, ΔzF, ΔzC
    end

    Nz = length(zF) - 1
    zF, zC, ΔzF, ΔzC = grid(zF)

    # Need some halo regions.
    ΔzF = OffsetArray([ΔzF[1], ΔzF...], 0:Nz)
    ΔzC = [ΔzC..., ΔzC[end]]

    # Useful for broadcasting z operations
    ΔzC = reshape(ΔzC, (1, 1, Nz))

    # Temporary hack: Useful for reusing fill_halo_regions! and
    # BatchedTridiagonalSolver which only need Nx, Ny, Nz.
    fake_grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))

    #####
    ##### Generate batched tridiagonal system coefficients and solver
    #####

    function λi(Nx, Δx)
        is = reshape(1:Nx, Nx, 1, 1)
        @. (2sin((is-1)*π/Nx) / Δx)^2
    end

    function λj(Ny, Δy)
        js = reshape(1:Ny, 1, Ny, 1)
        @. (2sin((js-1)*π/Ny) / Δy)^2
    end

    kx² = λi(Nx, Δx)
    ky² = λj(Ny, Δy)

    # Lower and upper diagonals are the same
    ld = [1/ΔzF[k] for k in 1:Nz-1]
    ud = copy(ld)

    # Diagonal (different for each i,j)
    @inline δ(k, ΔzF, ΔzC, kx², ky²) = - (1/ΔzF[k-1] + 1/ΔzF[k]) - ΔzC[k] * (kx² + ky²)

    d = zeros(Nx, Ny, Nz)
    for i in 1:Nx, j in 1:Ny
        d[i, j, 1] = -1/ΔzF[1] - ΔzC[1] * (kx²[i] + ky²[j])
        d[i, j, 2:Nz-1] .= [δ(k, ΔzF, ΔzC, kx²[i], ky²[j]) for k in 2:Nz-1]
        d[i, j, Nz] = -1/ΔzF[Nz-1] - ΔzC[Nz] * (kx²[i] + ky²[j])
    end

    #####
    ##### Random right hand side
    #####

    # Random right hand side
    Ru = CellField(Float64, arch, fake_grid, UVelocityBoundaryConditions(fake_grid))
    Rv = CellField(Float64, arch, fake_grid, VVelocityBoundaryConditions(fake_grid))
    Rw = CellField(Float64, arch, fake_grid, WVelocityBoundaryConditions(fake_grid))

    interior(Ru) .= rand(Nx, Ny, Nz)
    interior(Rv) .= rand(Nx, Ny, Nz)
    interior(Rw) .= zeros(Nx, Ny, Nz)

    U = (u=Ru, v=Rv, w=Rw)
    fill_halo_regions!(U, arch)

    _compute_w_from_continuity!(U, fake_grid)

    fill_halo_regions!(Rw, arch)

    R = zeros(Nx, Ny, Nz)
    for i in 1:Nx, j in 1:Ny, k in 1:Nz
        R[i, j, k] = div_ccc(i, j, k, Δx, Δy, ΔzF, Ru.data, Rv.data, Rw.data)
    end

    # @show sum(R)  # should be zero by construction.

    F = zeros(Nx, Ny, Nz)
    F = ΔzC .* R  # RHS needs to be multiplied by ΔzC

    #####
    ##### Solve system
    #####

    F̃ = fft(F, [1, 2])

    btsolver = BatchedTridiagonalSolver(arch, dl=ld, d=d, du=ud, f=F̃, grid=fake_grid)

    ϕ̃ = zeros(Complex{Float64}, Nx, Ny, Nz)
    solve_batched_tridiagonal_system!(ϕ̃, arch, btsolver)

    ϕ = CellField(Float64, arch, fake_grid, PressureBoundaryConditions(fake_grid))
    interior(ϕ) .= real.(ifft(ϕ̃, [1, 2]))
    ϕ.data .= ϕ.data .- mean(interior(ϕ))

    #####
    ##### Compute Laplacian of solution ϕ to test that it's correct
    #####

    fill_halo_regions!(ϕ, arch)

    ∇²ϕ = CellField(Float64, arch, fake_grid, PressureBoundaryConditions(fake_grid))
    for i in 1:Nx, j in 1:Ny, k in 1:Nz
        ∇²ϕ.data[i, j, k] = ∇²(i, j, k, Δx, Δy, ΔzF, ΔzC, ϕ.data)
    end

    return interior(∇²ϕ) ≈ R
end

@testset "Solvers" begin
    @info "Testing Solvers..."

    for arch in archs
        @testset "Batched tridiagonal solver [$arch]" begin
            for Nz in [8, 11, 18]
                @test can_solve_single_tridiagonal_system(arch, Nz)
                @test can_solve_single_tridiagonal_system_with_functions(arch, Nz)
            end

            for Nx in [3, 8], Ny in [5, 16], Nz in [8, 11]
                @test can_solve_batched_tridiagonal_system_with_3D_RHS(arch, Nx, Ny, Nz)
                @test can_solve_batched_tridiagonal_system_with_3D_functions(arch, Nx, Ny, Nz)
            end
        end
    end

    for arch in [CPU()]
        @testset "Vertically stretched Poisson solver [FACR, $arch]" begin
            @info "  Testing vertically stretched Poisson solver [FACR, $arch]..."

            Nx = Ny = 8
            zF = [1, 2, 4, 7, 11, 16, 22, 29, 37]
            @test vertically_stretched_poisson_solver_correct_answer(arch, Nx, Ny, zF)
        end
    end
end

ext: jl


text:
using Plots
using ForwardDiff
using LinearAlgebra
using Random
using Statistics

# https://openai.com/blog/evolution-strategies/

sol = [0.5; 0.1]
f(w) = -sum((w - sol).^2.0)
f(a, b) = f([a; b])

n = 1000
x_limits = range(-2.0, stop = 2.0, length = n)
y_limits = range(-2.0, stop = 2.0, length = n)

npop = 50
sigma = 0.1
alpha = 0.001

w = Random.randn(2)

plt = Plots.contour(x_limits, y_limits, f)
plt = plot!((w[1], w[2]), markershape = :circle, color = :cyan, outline = :cyan)

for i = 1:300
    N = Random.randn(npop, 2)
    R = zeros(npop)
    i % 10 == 0 && (plt = Plots.contour(x_limits, y_limits, f))

    for j = 1:npop
        w_try = w + sigma * N[j, :]
        R[j] = f(w_try)

        if i % 10 == 0
            plt = plot!((w_try[1], w_try[2]),
                markershape = :circle, color = :cyan, outline = :cyan,
                label = "")
        end
    end
    i % 10 == 0 && display(plt)

    A = (R .- mean(R)) ./ std(R)
    w = w + alpha / (npop * sigma) * N' * A

    # plt = plot!((w[1], w[2]),
    #     markershape = :circle, color = :cyan, outline = :cyan,
    #     label = "")
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# Purpose: This provides a type of array that is representable as a
# multiplication by two arrays. For example, the array `S`
#
#   S[1,2] = d1[1] * d2[2]
#
# is representable as the multiplication of elements from arrays `d1` and `d2`.
# Then, instead of saving the entire 2D matrix `S`, we can simply save the
# vectors `d1` and `d2`, saving space, and saving computation time as we can
# specialize some algorithms.
#
# It is not possible to write to a SeparableArray. This is because the meaning
# of an operation such as `S[3,4] = 5.0` is undefined: do you want to change
# the first array or the second array?

module SeparableArrays

export SeparableArray, @SeparableArray
export exponentiate

struct SeparableArray{T,N,A1,A2} <: AbstractArray{T,N}
    arr1::A1
    arr2::A2
    name1::Symbol
    name2::Symbol
end

function SeparableArray(arr1, arr2; name1=:arr1, name2=:arr2)
    T = promote_type(eltype(arr1), eltype(arr2))
    N = ndims(arr1) + ndims(arr2)
    A1 = typeof(arr1)
    A2 = typeof(arr2)
    return SeparableArray{T,N,A1,A2}(arr1, arr2, name1, name2)
end

macro SeparableArray(arr1, arr2)
    sym1 = "$arr1"
    sym2 = "$arr2"
    return :( SeparableArray($(esc(arr1)), $(esc(arr2)); name1=Symbol($sym1), name2=Symbol($sym2)) )
end

Base.size(A::SeparableArray) = size(A.arr1)..., size(A.arr2)...

function Base.getindex(A::SeparableArray, i::Int)
    N1 = length(A.arr1)
    m = (i-1) % N1 + 1
    n = (i-1) ÷ N1 + 1
    return A.arr1[m] * A.arr2[n]
end

function calc_index(I, arr, off, A)
    idx = I[1+off]
    sz = size(arr, 1)
    L = 1
    1 <= idx <= sz || throw(BoundsError(A, I))
    for i=2:ndims(arr)
        L *= sz
        sz = size(arr, i)
        1 <= I[i+off] <= sz || throw(BoundsError(A, I))
        idx += (I[i+off] - 1) * L
    end
    return arr[idx]
end

function Base.getindex(A::SeparableArray{T,N,A1,A2}, I::Vararg{Int,N}) where {T,N,A1,A2}
    arr1 = getfield(A, :arr1)
    D1 = ndims(arr1)
    V1 = calc_index(I, arr1, 0, A)
    arr2 = getfield(A, :arr2)
    V2 = calc_index(I, arr2, D1, A)
    return V1 * V2
end

function Base.getproperty(s::SeparableArray, f::Symbol)
    #@show "getproperty",f
    if f == getfield(s, :name1)
        return getfield(s, :arr1)
    elseif f == getfield(s, :name2)
        return getfield(s, :arr2)
    else
        return getfield(s, f)
    end
end

function Base.propertynames(s::SeparableArray, private::Bool=false)
    tuple = (s.name1, s.name2)
    if private
        tuple = tuple..., fieldnames(typeof(s))...
    end
    return Tuple(unique(tuple))
end

# Optimizations
@doc raw"""
    exponentiate(s::SeparableArray, exp::Number)

This function is used for elementwise exponentiation of the array 's'. It could
be made more elegant by extending the broadcast syntax. PRs welcome.
"""
function exponentiate(s::SeparableArray, exponent::Number)
    a1 = s.arr1 .^ exponent
    a2 = s.arr2 .^ exponent
    return SeparableArray(a1, a2, name1=s.name1, name2=s.name2)
end
#Base.broadcasted(::typeof(^), s::SeparableArray, e::Number) = exponentiate(s, e)

end

# vim: set sw=4 et sts=4 :


text:
########################################################################
# Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Submitted to Science, 2016
# Code in this file by Michael K. Borregaard
#########################################################################

#Load the necessary library
using DataFrames
include("GD_summary_functions.jl")

# read in the mammals and subset to valid sequence comparisons
mamcyt = readtable("pairwise_equalarea_mammals_cytb.csv")
mamcyt = mamcyt[!isna(mamcyt[:seq1]),:]
mamcyt = mamcyt[mamcyt[:overlap] .>= 0.5,:]

# read in the amphibians and subset to valid sequence comparisons
amphs = readtable("pairwise_equalarea_amphibians_cytb.csv")
amphs = amphs[!isna(amphs[:seq1]),:]
amphs = amphs[amphs[:overlap] .>= 0.5,:]

# Calculate the values and write them to disk
amp = totalbasepairs(amphs)
amp[:genespecies] = totalspecies(amphs)[:x1]
writetable("amphfig2.csv", amp)

mam = totalbasepairs(mamcyt)
mam[:genespecies] = totalspecies(mamcyt)[:x1]
writetable("mamfig2.csv", mam)

total = vcat(mamcyt, amphs)
fig2 = totalbasepairs(total)
fig2[:genespecies] = totalspecies(total)[:x1]
writetable("totalfig2.csv", fig2)

max_stars_repo_path: Miraldo et al. 2016/Calculate_Values_for_Figure_S10_Fig_2_cytb.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
include("header.jl")

#=
# Example calculation in a pipe system
In this notebook we will see an example of how to solve a problem in a pipe
system, using the 1-d steady-flow energy equation.

In our calculations, we will use the units that we discussed in notebook 1.0
in order to take advantage of the useful tools that come with them.
=#

# ### Set up the module
using MAE103

#=
## Example: Find the flow rate
In this example (Example 8.10 in the textbook), we are to determine the change of
flow rate that occurs after we modify a design for a fume hood. The fume hood
involves a fan that pumps air and noxious gas out of an enclosed region,
through a duct, and out to the atmosphere where it can safely mix with the
ambient air. The existing design involves a fan with a short duct of diameter
8 inches, and the losses correpond to a loss coefficient of $K_L = 5$. This
leads to a flow rate of $\dot{Q}_0 = 9$ ft$^3$/s, which is safely within
the regulated range of 6 to 12 ft$^3$/s.

The new design involves a long 100 ft pipe of galvanized iron (leading to
major loss from viscosity) and a new total loss coefficient of $K_L = 10$.
We wish to see if the new design's volume flow rate falls within the
regulated range. We will assume that the fan's head does not change in
the redesign.
=#
#-
#=
First, let's put in the parameters. We treat the gas as air, since any noxious
gas species are assumed to be at low concentrations.

Galvanized iron has a typical roughness of 0.0005 ft.
=#
Q0 = VolumeFlowRate(9u"ft^3/s")
KL0 = 5
D = Diameter(8u"inch")
μ = Viscosity(Air)
ρ = Density(Air)
ϵ = Height(0.0005u"ft")
g = Gravity()
#=
Now let's calculate the flow velocity in the original design
=#
A = Area(π/4*D^2)
V0 = Velocity(Q0/A)
value(V0,u"ft/s")

#=
And now, we can calculate the fan head $h_p$, since this is equal
to the kinetic energy of the exiting flow plus the head loss. (Pressures
are equal and ambient in the fume hood and at the exit, and elevation
changes are negligible.

$$h_p = \dfrac{V^2}{2g} + K_L\dfrac{V^2}{2g}$$
=#
hp = Head((KL0+1)*V0^2/(2*g))
value(hp,u"ft")

#=
Now we will analyze the modified system. We add 100 ft of duct after the
fan and increase the loss coefficient to 10:
=#
L = Length(100u"ft")
KL = 10.0

#=
The energy equation now becomes

$$h_p = \dfrac{V^2}{2g} + K_L\dfrac{V^2}{2g} + \dfrac{fL}{D}\dfrac{V^2}{2g}$$

The unknown in this equation is $V$, but $f$ (the friction factor) is also
unknown. Let's solve for $V$, pretending that we know $f$:

$$V = \left(\dfrac{2gh_p}{1 + K_L + fL/D}\right)^{1/2}$$

To get $f$, we need Reynolds number and roughness coefficient.
And to get Reynolds number, we need $V$. So, we clearly must iterate.

Our approach will be
0. Guess a value for $f$, based on $\epsilon/D$ and a Reynolds number close to $\infty$.
1. Calculate $V$ from the energy equation with the current guess of $f$.
2. Calculate $Re_D = \rho V D/\mu$ from $V$.
3. Calculate a new value of $f$. Check if it is equal to the guess used in 1. If
yes, we stop. If not, then return to step 1 with this $f$ as our new guess.

Let's prepare ourselves for this iteration by defining some equations. The
roughness coefficient:
=#
eD = RoughnessCoefficient(ϵ/D)

# Here is $V$ as a function of $f$ (and the other quantities, which don't change).
Vfromf(f) = Velocity(sqrt(2*g*hp/(1+KL+f*L/D)))

# Here is $Re_D$ as a function of $V$ (and other quantities that don't change)
ReD(V) = ReynoldsNumber(ρ*V*D/μ)

#=
To do the iteration, we will define a simple function that takes in
an initial guess for $f$, iterates using the procedure above, and then returns
the correct value.

In the function, we create a quantity `f_old`, which simply hold our old guess
for $f$, and `f` will hold our next guess. We will compare these in each iteration,
checking whether they are nearlyequal. We'll say they're equal if they are closer
than $10^{-8}$. To ensure that the test fails at first, we initialize `f_old` to
infinity. This will force it to proceed into the iteration loop:
=#
function iterate_f(f0::FrictionFactor)
  f = f0
  f_old = Inf
  while abs(f - f_old) > 1e-8
      f_old = f
      f = FrictionFactor(ReD(Vfromf(f_old)),eD)
  end
  return f
end

#=
Now, to get $f$ from $Re_D$ and $\epsilon/D$, we use the function
`FrictionFactor(Re,eD)`, which solves the Colebrook equation for $f$.
Here, we use it to find our initial guess for $f$, treating $Re_D$ as very large.
=#
f0 = FrictionFactor(ReynoldsNumber(1e10),eD)

f = iterate_f(f0)

# We converged on the final $f$ value! The actual velocity is thus
V = Vfromf(f)
value(Vfromf(f),u"ft/s")

# And the final Reynolds number is
ReD(V)

# And finally, the volume flow rate after making the change to the setup is
Q = VolumeFlowRate(V*A)
value(Q,u"ft^3/s")

#=
So this shows that the flow rate is too small to meet the specifications.
We would need to decrease the length, $L$, or reduce the losses in $K_L$
to make it meet the specifications.
=#


text:
printstyled("  time\n", color=:light_green)

# Timestamp
t0 = time()
ts0 = String(split(string(u2d(t0)), '.')[1])
ts1 = String(split(timestamp(t0), '.')[1])
ts2 = String(split(timestamp(u2d(t0)), '.')[1])
ts3 = String(split(timestamp(string(u2d(t0))), '.')[1])
ts4 = String(split(timestamp(DateTime(string(u2d(t0)))), '.')[1])
@test ts1 == ts2 == ts3 == ts4 == ts0

# j2md
@test j2md(1, 1) == (1,1)
@test j2md(2000, 1) == (1,1)
@test j2md(2000, 60) == (2,29)
@test j2md(2000, 366) == (12,31)
@test_throws BoundsError j2md(2000, 367)
@test j2md(2012, 366) == (12,31)
@test j2md(2013, 60) == (3,1)
@test j2md(2015, 60) == (3,1)
@test j2md(2015, 365) == (12, 31)
@test j2md(2016, 61) == (3,1)
@test_throws BoundsError j2md(2013, 366)
@test_throws BoundsError j2md(2015, 366)
@test j2md(2100, 60) == (3,1)

# md2j
@test md2j(2000,1,1) == 1
@test md2j(2000,2,29) == 60
@test md2j(2013,3,1) == 60
@test md2j(2100,3,1) == 60
@test md2j(2001,1,1) == 1
@test 60 == md2j(2015,3,1)
@test 61 == md2j(2016,3,1)
@test 365 == md2j(2015,12,31)
@test 365 == md2j(1900,12,31)

# parsetimewin
t1 = t0 + 86400.0
dt0 = u2d(t0)
dt1 = u2d(t1)
st0 = string(dt0)
st1 = string(dt1)

d0, d1 = parsetimewin(dt0, dt1);  @test d1 > d0                      # dt, dt
d0, d1 = parsetimewin(dt0, t1);   @test d1 > d0                      # dt, r
d0, d1 = parsetimewin(dt0, st1);  @test d1 > d0                      # dt, s
d0, d1 = parsetimewin(t0, dt1);   @test d1 > d0                      # r, dt
d0, d1 = parsetimewin(t0, t1);    @test d1 > d0                      # r, r
d0, d1 = parsetimewin(t0, st1);   @test d1 > d0                      # r, s
d0, d1 = parsetimewin(st0, dt1);  @test d1 > d0                      # s, dt
d0, d1 = parsetimewin(st0, t1);   @test d1 > d0                      # s, r
d0, d1 = parsetimewin(st0, st1);  @test d1 > d0                      # s, s

# Checking that parsetimewin sorts correctly; type mismatch intentional
d0, d1 = parsetimewin(0, -600);  @test d1 > d0
d0, d1 = parsetimewin(-600.0, 0);  @test d1 > d0
d0, d1 = parsetimewin(600, 0.0);  @test d1 > d0
d0, d1 = parsetimewin(0.0, 600.0);  @test d1 > d0

@test ≈(600000, (DateTime(d1)-DateTime(d0)).value)
d0, d1 = parsetimewin("2016-02-29T23:30:00", "2016-03-01T00:30:00")
@test ≈(3600000, (DateTime(d1)-DateTime(d0)).value)
t = DateTime(now())
s = t-Hour(2)
d0, d1 = parsetimewin(s, t)
@test ≈(7200000, (DateTime(d1)-DateTime(d0)).value)

# t_collapse, t_expand
ts = 1451606400000000
T = Int64[1 ts; 100001 30000000; 250001 12330000; 352303 99000000; 360001 0]
fs = 100.0
Δ = round(Int64, sμ/fs)
t_long = t_expand(T, fs)
@test ≈(T, t_collapse(t_long, fs))

T1 = hcat(cumsum(ones(Int64,size(T,1))), cumsum(T[:,2]))
fs1 = 0.0
@test ≈(T1, t_collapse(t_expand(T1, fs1), fs1))

# starttime
printstyled(stdout, "    starttime\n", color=:light_green)
@test starttime(T, fs) == first(t_long)
@test starttime(T, Δ) == starttime(T, fs)
@test starttime(Array{Int64,2}(undef, 0, 0), Δ) == 0
@test starttime(Array{Int64,2}(undef, 0, 0), 100Δ) == 0
@test starttime(Array{Int64,2}(undef, 0, 0), fs) == 0

# why starttime exists: it handles segments that aren't in chronological order
printstyled(stdout, "      non-chronological\n", color=:light_green)
Tnco = Int64[1 ts; 1000 -2000Δ; 120000 0]
t_long = t_expand(Tnco, fs);
@test minimum(t_long) != first(t_long)
@test starttime(Tnco, fs) == starttime(Tnco, Δ) == ts - 1001Δ == minimum(t_expand(Tnco, fs))

# endtime
printstyled(stdout, "    endtime\n", color=:light_green)
T = Int64[1 ts; 100001 30000000; 250001 12330000; 352303 99000000; 360001 0]
t_long = t_expand(T, fs)
@test endtime(T, fs) == last(t_long)
@test endtime(T, Δ) == endtime(T, fs)
@test endtime(Array{Int64,2}(undef, 0, 0), Δ) == 0
@test endtime(Array{Int64,2}(undef, 0, 0), 100Δ) == 0
@test endtime(Array{Int64,2}(undef, 0, 0), fs) == 0

# endtime must also handle segments that aren't in chronological order
printstyled(stdout, "      non-chronological\n", color=:light_green)
nx = 120000
igap = 1000
tgap = 10000
Tnco = Int64[1 ts; nx-2000 tgap*Δ; nx-(igap-1) -8000Δ; nx 0]
t_long = t_expand(Tnco, fs);
@test maximum(t_long) != last(t_long)
@test endtime(Tnco, fs) == endtime(Tnco, Δ) == ts + (nx-igap)*Δ + (tgap-1)*Δ == maximum(t_long)

printstyled(stdout, "    t_win, w_time\n", color=:light_green)
printstyled(stdout, "      faithful representation of gaps\n", color=:light_green)
fs = 100.0
Δ = round(Int64, sμ/fs)
t = [1 0; 6 980000; 8 100000; 10 0]
# t1 = t_expand(t, fs)
 #       0
 #   10000
 #   20000
 #   30000
 #   40000
 # 1030000
 # 1040000
 # 1150000
 # 1160000
 # 1170000

t2 = t_win(t, Δ)
 # first      last
 # sample     sample
 # time       time
 #       0    40000 ===> so that collect(w[i,1]:Δ:w[i,2]) for each window == t_expand
 # 1030000  1040000
 # 1150000  1170000
 #
@test t2[1,2] == 40000
@test t2[2,:] == [1030000, 1040000]
@test t2[3,:] == [1150000, 1170000]

printstyled(stdout, "      arbitrary windows containing gaps\n", color=:light_green)
t = [1 999999998990000; 101 10000; 297 0]
@test w_time(t_win(t, Δ), Δ) == t

t = [1 999999998990000; 101 10000; 297 1000000; 303 40000; 500 1000000; 10000 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      length-0 gap\n", color=:light_green)
t = [1 0; 6 0; 8 0; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      negative gap\n", color=:light_green)
t = [1 0; 6 -2Δ; 8 -10Δ; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      single-point gap\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 0]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      non-null gap at end\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 Δ]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "      negative gap at end\n", color=:light_green)
t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 -5Δ]
@test w_time(t_win(t, Δ), Δ) == t

printstyled(stdout, "    t_bounds\n", color=:light_green)
function test_t_bounds(t::Array{Int64, 2}, Δ::Int64)
  (t0, t1) = t_bounds(t, Δ)
  W = t_win(t, Δ)
  @assert t0 == minimum(W)
  @assert t1 == maximum(W)
  return nothing
end

fs = 100.0
Δ = round(Int64, sμ/fs)
t = [1 0; 6 980000; 8 100000; 10 0]
test_t_bounds(t, Δ)

t = [1 999999998990000; 101 10000; 297 0]
test_t_bounds(t, Δ)

t = [1 999999998990000; 101 10000; 297 1000000; 303 40000; 500 1000000; 10000 0]
test_t_bounds(t, Δ)

t = [1 0; 6 0; 8 0; 10 0]
test_t_bounds(t, Δ)

t = [1 0; 6 -2Δ; 8 -10Δ; 10 0]
test_t_bounds(t, Δ)

t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 0]
test_t_bounds(t, Δ)

t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 Δ]
test_t_bounds(t, Δ)

t = [1 0; 6 2Δ; 7 4Δ; 8 Δ; 10 -5Δ]
test_t_bounds(t, Δ)

printstyled(stdout, "    mktime\n", color=:light_green)
fv = 0.0005
iv = Array{Int32,1}([1980, 082, 10, 35, 39, 890])
(m,d) = j2md(iv[1], iv[2])

ts_0 = round(Int64, d2u(DateTime(iv[1], m, d, iv[3], iv[4], iv[5]))*sμ) +
       iv[6]*1000 +
       round(Int64,fv*1000.0)
ts_1 = Date(iv[1], m, d).instant.periods.value * 86400000000 +
       div(Time(iv[3], iv[4], iv[5]).instant.value, 1000) +
       iv[6]*1000 +
       round(Int64,fv*1000.0) -
       SeisIO.dtconst
ts_2 = mktime(iv[1], iv[2], iv[3], iv[4], iv[5], iv[6]*Int32(1000)) +
       round(Int64, fv*1000.0)
iv[6] *= Int32(1000)
ts_3 = mktime(iv) + round(Int64, fv*1000.0)
@test ts_0 == ts_1 == ts_2 == ts_3
# timespec()

printstyled(stdout, "    int2tstr, tstr2int\n", color=:light_green)
s = "2018-01-01T00:00:00.000001"
t = "2018-01-04T00:00:00.003900"
si = tstr2int(s)
ti = tstr2int(t)
j = loop_time(si, ti)

t = "2018-01-04T00:00:00.39"
for (n,s) in enumerate(["2018-01-01T00:00:00.000001",
                        "2018-01-01T00:00:00",
                        "2018-01-01T00:00:00.035",
                        "2016-02-29T00:00:00.02",
                        "2018-02-28T00:00:00.33"])

  s_str = identity(s)
  if length(s) == 19
    s_str *= "."
  end
  s_str = rpad(s_str, 26, '0')
  @test (int2tstr(tstr2int(s))) == s_str
  if n < 4
    @test loop_time(tstr2int(s), tstr2int(t)) == 4
    @test loop_time(tstr2int(s), tstr2int(t), ti=43200000000) == 7
    @test loop_time(tstr2int(s), tstr2int(t), ti=3600000000) == 73
  elseif n == 4
    @test loop_time(tstr2int(s), tstr2int("2016-03-02T00:00:02")) == 3
  elseif n == 5
    @test loop_time(tstr2int(s), tstr2int("2018-03-02T00:00:02")) == 3
  end
end

printstyled(stdout, "    sort_segs!\n", color=:light_green)
Δ = 20000
ts = 1583455810004000
nx = 40000
gi = 10
gl = 6
W = ts .+ Δ.*[    0     gi-1
                -gl       -1
                 gi  nx-gl-1]
W0 = deepcopy(W)
sort_segs!(W)
@test W == W0[[2,1,3], :]

t = sort_segs(w_time(W0, Δ), Δ)
@test t == w_time(W, Δ)

printstyled(stdout, "    t_extend\n", color=:light_green)
printstyled(stdout, "      time-series\n", color=:light_green)
nx = 6000
Δ = 10000
g = 300000
fs = 100.0
t0 = 1411776000000000
t = [1 t0; nx 0]
ts = t[1,2] + nx*Δ
t2 = [1 ts; nx 0]

# should have only 2 rows
t1 = deepcopy(t)
@test t_extend(t1, ts, nx, fs) == nothing
@test t1 == [1 t0; 2nx 0]

# should have a gap of Δ at point nx+1
@test t_extend(t, ts + Δ, nx, fs) == [1 t0; nx+1 Δ; 2nx 0]

# should correctly log the gaps at nx and nx+1
t1 = [1 t0; nx g]
@test t_extend(t1, ts+g+Δ, nx, fs) == [1 t0; nx g; nx+1 Δ; 2nx 0]

# should only extend the expected length of nx to incorporate new start time
for ts1 in ts : 100Δ : ts + nx
  t1 = deepcopy(t)
  # @test t_extend(t, ts1, 0, fs) == [1 t[1,2]; nx + div(ts1-endtime(t, Δ), Δ) 0]
  @test t_extend(t1, ts1, 0, fs) == nothing
  @test t1 == [1 t[1,2]; nx-1 + div(ts1-endtime(t, Δ), Δ) 0]
end

# should initialize a new starter time array at ts
@test t_extend(Array{Int64,2}(undef,0,2), ts, 0, fs) == [1 ts]

t = [1 12356; 1231 333; 14134 0]
ts_new = 8348134123
nx_new = 65536
dt = 20000
fs = 50.0
t2 = deepcopy(t)
t_extend(t2, ts_new, 0, dt)
t1 = t_extend(t, ts_new, nx_new, dt)
 #     1       12356
 #  1231         333
 # 14135  8065441434
 # 79670           0
@test size(t1) == (4,2)
@test t1[end,1] == 79670
@test endtime(t1, dt) == t_expand(t1, fs)[end]
@test t1[1:end-2,:] == t2[1:end-1,:]

t = [1 3301; 505 1200; 1024 3]
ts_new = 1181381433
nx_new = 3000
t2 = t_extend(t, ts_new, 1, dt)
t1 = t_extend(t, ts_new, nx_new, dt)
# 1           3301
# 505         1200
# 1024           3
# 1025  1160896929
# 4024           0
@test endtime(t1, dt) == t_expand(t1, fs)[end]
@test t1[1:end-1,:] == t2

t1 = Array{Int64,2}(undef, 0, 0)
t2 = Array{Int64,2}(undef, 0, 2)
@test t_extend(t1, ts_new, nx_new, dt) == t_extend(t2, ts_new, nx_new, dt) == [1 ts_new; nx_new 0]
@test t_extend(t1, ts_new, 0, dt) == t_extend(t2, ts_new, 0, dt) == [1 ts_new]

printstyled(stdout, "      irregular\n", color=:light_green)
nx = 200
nn = 120
t1 = zeros(Int64, nx, 2)
t1[1:nx, 1] .= 1:nx
t1[1:nx, 2] .= sort(abs.(rand(1262304000000000:t0, nx)))
t2 = zeros(Int64, nn, 2)
t2[1:nn, 1] .= 1:nn
t2[1:nn, 2] .= sort(abs.(rand(t0+1:1583020800000000, nn)))
t1o = deepcopy(t1)
t2o = deepcopy(t2)
t3 = t_extend(t1, t2[1,2], nn, 0.0)
@test t1 == t1o
@test t2 == t2o
@test size(t3, 1) == nx+nn
@test t3[:, 1] == collect(1:1:nx+nn)
@test t3[nx+1, 2] == t2[1,2]

# These were in test_time_utils.jl
buf = BUF.date_buf
dstr = "2019-06-01T03:50:04.02"
dt = DateTime(dstr)
t = round(Int64, d2u(dt)*sμ)
nx = 12345

# Tests for mk_t
printstyled("    mk_t\n", color=:light_green)
C = randSeisChannel(s=true)
C.x = randn(nx)
mk_t!(C, nx, t)
@test C.t == [1 t; nx 0]
@test C.t == mk_t(nx, t)

# Tests for t_arr
printstyled("    t_arr!\n", color=:light_green)
t_arr!(buf, t)
@test buf[1:6] == Int32[2019, 152, 3, 50, 4, 20]
t = round(Int64, d2u(DateTime("2020-03-01T13:49:00.3"))*sμ)
t_arr!(buf, t)
@test buf[1:6] == Int32[2020, 61, 13, 49, 0, 300]
t = round(Int64, d2u(DateTime("2020-03-01T13:49:00.030"))*sμ)
t_arr!(buf, t)
@test buf[1:6] == Int32[2020, 61, 13, 49, 0, 30]
t = round(Int64, d2u(DateTime("2020-03-01T13:49:00.003"))*sμ)
t_arr!(buf, t)
@test buf[1:6] == Int32[2020, 61, 13, 49, 0, 3]

# Tests for x_inds
printstyled("    x_inds\n", color=:light_green)
function test_xinds(t::Array{Int64, 2})
  # Check that the sets of indices match what's in :t
  xi = x_inds(t)
  nt = size(t, 1)
  @test xi[1:nt-1,1] == t[1:nt-1,1]

  # Check that the length is what we expect
  @test size(xi, 1) == nt - (t[nt,2] == 0 ? 1 : 0)

  # Check xi and t_win(t, dt) have the same content
  for i in 1:size(xi,1)-1
    @test xi[i,2] == xi[i+1,1]-1
    for dt in ([250, 500, 1000, 10000, 20000, 100000])
      w = t_win(t, dt)
      @test xi[i,2]-xi[i,1] == div(w[i,2]-w[i,1], dt)
    end
  end
end

ts = 1583455810004000
for t in [[1  ts; 639000 0],
          [1  ts; 639000 57425593],
          [1  ts; 638999 300; 639000 57425593],
          [1  ts; 638998 1234; 638999 300; 639000 57425593],
          [1  ts; 638998 -1234525827; 639000 0],
          [1  ts; 638998 1234; 639000 0],
          [1 ts; 8 1234; 9 100000; 10 134123; 11 12; 2400 -3; 3891 3030303; 3892 -30000000; 3893 1234; 3894 57425593],
          [1 ts; 3894 0],
          [1 ts; 10 13412300; 11 123123123; 184447 3030303; 184448 -30000000; 184449 12300045; 184450 57425593]]
  test_xinds(t)
end
for i in 1:100
  t = breaking_tstruct(ts, 39000, 100.0)
  test_xinds(t)
end


text:
function juliacon2021(::Val{:terminal})
    if myid() == 1
        return println(
            "Welcome to JuliaCon 2021! Find more information on https://juliacon.org/2021/."
        )
    else
        return println("Greetings from ", rand(countries), "!")
    end
    return nothing
end

# TODO: needs love for a distributed version based on the :terminal method (no hurry though)
function juliacon2021(::Val{:text})
    return "Welcome to JuliaCon 2021! Find more information on https://juliacon.org/2021/."
end

juliacon2021(; output=:terminal) = juliacon2021(Val(output))


text:
function straight_line()

    g = 0:0.02:1

    M_bf = 2

    ζ_vec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]

    # For each m there are many non-unique solution
    I_sols =
    [[[1, length(g)]];
     [Vector{Int}() for m=1:9]]

    f_sols = zeros(Float64, 10)

    return g, ζ_vec, I_sols, f_sols
end


text:
0 < incx || throw(ArgumentError("non-positive increment"))
0 < ix || throw(BoundsError())
ix+(n-1)*incx <= length(x) || throw(BoundsError())


text:
function save_loc()
	loc=""
	if haskey(ENV,"FW_SAVE_LOC")
		loc=ENV["FW_SAVE_LOC"]
	else
		loc=ENV["APPDATA"]*"/FoamWorld/sav/"
		ENV["FW_SAVE_LOC"]=loc
	end
end
function save_game(g::Game,path::String=save_loc()*g.set["name"])
	if !isdir(path)
		mkdir(path)
	end
	cd(path)
	write_code("setting.fw",g.set)
end
function load_game(path::String)
	g=Game()
	cd(path)
	g.set=read_code("setting.fw")
end
function get_loadedimg(s::String)
	if haskey(loadedimgs,s)
		return @inbounds loadedimgs[s]
	else
		return loadedimgs["notexture"]
	end
end


text:
import Pkg
Pkg.activate("./")

using PDMesh

obj = Cuboid([-10 10; 0 3; -2 2])

f = out -> begin
x=out[:x];
mask = (x[1, :] .< 0.1) .& (x[1, :] .> -0.1) .& (x[2, :] .> 2)
mask
end

obj = delete(obj,  f)

x, v, y, vol, type = create(obj, resolution=0.1, rand_=0.0, type=1)

write_data("./examples/ovito_ref.data", x, type, vol)

write_data("./examples/ovito.data", y, v, type, vol)

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
prices = Dict(
    "milk" => 1,
    "eggs" => 2,
    "cheese" => 3.5,
    "bread" => 0.9,
    "chocolate" => 2.2,
    "flour" => 0.5,
    "tomatoes" => 1.5
)

stock = Dict(
    "milk" => 10,
    "eggs" => 20,
    "cheese" => 35,
    "bread" => 20,
    "chocolate" => 22,
    "flour" => 15,
    "tomatoes" => 15
)

mutable struct Store
    stock::Dict
    prices::Dict
    funds::Float64
end

store = Store(prices, stock, 100)
println(store.stock["milk"])
println(store.prices["tomatoes"])
println(sum(values(store.prices)))

map(sum, collect(values(prices)))
map(sum, collect(values(prices)))
println(keys(store.prices))


text:
using ArtGallery
using Test

@test album("computer") == nothing


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays
using LinearMaps: LinearMapVector, LinearMapTuple

@testset "composition" begin
    F = @inferred LinearMap(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
    FC = @inferred LinearMap{ComplexF64}(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
    FCM = LinearMaps.CompositeMap{ComplexF64}((FC,))
    L = LowerTriangular(ones(10,10))
    @test_throws DimensionMismatch F * LinearMap(rand(2,2))
    @test_throws ErrorException LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10))))
    A = 2 * rand(ComplexF64, (10, 10)) .- 1
    B = rand(size(A)...)
    H = LinearMap(Hermitian(A'A))
    S = LinearMap(Symmetric(real(A)'real(A)))
    M = @inferred 1 * LinearMap(A)
    N = @inferred LinearMap(B)
    v = rand(ComplexF64, 10)
    α = rand(ComplexF64)
    @test FCM * v == F * v
    @test @inferred (F * F) * v == @inferred F * (F * v)
    @test @inferred (F * A) * v == @inferred F * (A * v)
    @test @inferred (A * F) * v == @inferred A * (F * v)
    @test @inferred A * (F * F) * v == @inferred A * (F * (F * v))
    F2 = F*F
    FC2 = FC*FC
    F4 = FC2 * F2
    F4v = @inferred LinearMaps.CompositeMap{ComplexF64}([FC2, F2])
    @test occursin("10×10 LinearMaps.CompositeMap{$(eltype(F4))}", sprint((t, s) -> show(t, "text/plain", s), F4))
    @test length(F4.maps) == 4
    @test (@inferred F4 * v) == F4v * v == (@inferred F * (F * (F * (F * v))))
    @test @inferred Matrix(M * transpose(M)) ≈ A * transpose(A)
    @test @inferred !isposdef(M * transpose(M))
    @test @inferred isposdef(LinearMap(M * M', isposdef=true))
    @test @inferred issymmetric(N * N')
    @test @inferred ishermitian(N * N')
    @test @inferred !issymmetric(M' * M)
    @test @inferred ishermitian(M' * M)
    @test @inferred issymmetric(F'F)
    @test @inferred issymmetric(F'*S*F)
    @test @inferred ishermitian(F'F)
    @test @inferred ishermitian(LinearMaps.CompositeMap{ComplexF64}([F, F']))
    @test @inferred ishermitian(F'*H*F)
    @test @inferred ishermitian(LinearMaps.CompositeMap{ComplexF64}([F, H, F']))
    @test @inferred !issymmetric(FC'FC)
    @test @inferred ishermitian(FC'FC)
    @test @inferred ishermitian(FC'*H*FC)
    @test @inferred issymmetric(transpose(F) * F * 3)
    @test @inferred issymmetric(transpose(F) * 3 * F)
    @test @inferred !isposdef(-5*transpose(F) * F)
    @test @inferred ishermitian((M * F)' * M * 4 * F)
    @test @inferred transpose(M * F) == @inferred transpose(F) * transpose(M)
    @test @inferred (4*((-3*M)*2)) == @inferred -12M*2
    @test @inferred (4*((3*(-M))*2)*(-5)) == @inferred -12M*(-10)
    L = @inferred 3 * F + 1im * A + F * M' * F
    LF = 3 * Matrix(F) + 1im * A + Matrix(F) * Matrix(M)' * Matrix(F)
    @test Array(L) ≈ LF
    R1 = rand(ComplexF64, 10, 10); L1 = LinearMap(R1)
    R2 = rand(ComplexF64, 10, 10); L2 = LinearMap(R2)
    R3 = rand(ComplexF64, 10, 10); L3 = LinearMap(R3)
    CompositeR = *(R1, R2, R3)
    CompositeL = prod(LinearMap, [R1, R2, R3])
    @test prod([L1, L2, LinearMap(randn(10, 10))]) isa LinearMaps.CompositeMap{ComplexF64,<:LinearMapVector}
    @test (@inferred prod([L1, L2, L3])) isa LinearMaps.CompositeMap{<:Any,<:LinearMapVector}
    @test prod([L1, L2, L3]) == (@inferred L1 * prod([L2, L3])) == (@inferred prod([L1, L2]) * L3)
    @test (@inferred prod((L1, L2, L3))) isa LinearMaps.CompositeMap{<:Any,<:LinearMapTuple}
    @test (@inferred L1 * L2 * L3) == CompositeL == prod([L1, L2, L3])
    @test Matrix(L1 * L2) ≈ sparse(L1 * L2) ≈ R1 * R2
    @test Matrix(@inferred((α * L1) * (L2 * L3))::LinearMaps.ScaledMap) ≈ α * CompositeR
    @test Matrix(@inferred((L1 * L2) * (L3 * α))::LinearMaps.ScaledMap) ≈ α * CompositeR
    @test @inferred transpose(CompositeL) == transpose(L3) * transpose(L2) * transpose(L1)
    @test @inferred adjoint(CompositeL) == L3' * L2' * L1'
    @test @inferred adjoint(adjoint((CompositeL))) == CompositeL
    @test transpose(transpose((CompositeL))) == CompositeL
    Lt = @inferred transpose(LinearMap(CompositeL))
    @test Lt * v ≈ transpose(CompositeR) * v
    Lc = @inferred adjoint(LinearMap(CompositeL))
    @test Lc * v ≈ adjoint(CompositeR) * v

    # convert to AbstractMatrix
    for A in (LinearMap(sprandn(10, 10, 0.3)), LinearMap(rand()*I, 10))
        for B in (LinearMap(sprandn(10, 10, 0.3)), LinearMap(rand()*I, 10))
            AA = convert(AbstractMatrix, A*B)
            if A isa LinearMaps.UniformScalingMap && B isa LinearMaps.UniformScalingMap
                @test isdiag(AA)
            else
                @test issparse(AA)
            end
        end
    end

    # test inplace operations
    w = similar(v)
    mul!(w, L, v)
    @test w ≈ LF * v

    # test composition of several maps with shared data #31
    global sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
    N = length(sizes) - 1
    global Lf = []
    global Lt = []
    global Lc = []

    # build list of operators [LN, ..., L2, L1] for each mode
    for (fi, i) in [ (Symbol("f$i"), i) for i in 1:N]
        @eval begin
            function ($fi)(source)
                dest = ones(prod(sizes[$i + 1]))
                tmp = reshape(source, sizes[$i])
                return conj.($i * dest)
            end
            insert!(Lf, 1, LinearMap($fi, prod(sizes[$i + 1]), prod(sizes[$i])))
            insert!(Lt, 1, LinearMap(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
            insert!(Lc, 1, LinearMap{ComplexF64}(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
        end
    end
    @test size(prod(Lf[1:N])) == (prod(sizes[end]), prod(sizes[1]))
    @test isreal(prod(Lf[1:N]))
    # multiply as composition and as recursion
    v1 = ones(prod(sizes[1]))
    u1 = ones(prod(sizes[1]))
    w1 = im.*ones(ComplexF64, prod(sizes[1]))
    for i = N:-1:1
        v2 = prod(Lf[i:N]) * ones(prod(sizes[1]))
        u2 = transpose(LinearMap(prod(Lt[N:-1:i]))) * ones(prod(sizes[1]))
        w2 = adjoint(LinearMap(prod(Lc[N:-1:i]))) * ones(prod(sizes[1]))

        v1 = Lf[i] * v1
        u1 = transpose(Lt[i]) * u1
        w1 = adjoint(Lc[i]) * w1

        @test v1 == v2
        @test u1 == u2
        @test w1 == w2
    end
    L1 = LinearMap(rand(2,3))
    L2 = LinearMap(rand(4,2))
    L3 = LinearMap(rand(3, 4))
    L4 = LinearMap(rand(5, 3))
    Ls = L4*L3*L2*L1
    X = rand(size(Ls, 2), 10)
    Y = similar(X, (size(Ls, 1), size(X, 2)))
    @test mul!(Y, Ls, X) ≈ L4.lmap * L3.lmap * L2.lmap * L1.lmap * X

    # test isposdef on a case where sufficient conditions work
    B = LinearMap([1 0; 0 1], isposdef=true) # isposdef!
    C = B' * B * B * B * B # no B' at end on purpose
    @test @inferred isposdef(C)
    @test @inferred isposdef(LinearMaps.CompositeMap{Float64}([B, B, B, B, B']))
    @test @inferred isposdef(B * B) # even case for empty tuple test
    @test @inferred isposdef(LinearMaps.CompositeMap{Float64}([B, B]))

    M = LinearMap(cumsum, 3)
    for i in 1:4
        P = prod(fill(M, i))
        @test P isa LinearMaps.CompositeMap{<:Any,<:LinearMapVector}
        @test P * ones(3) == (LowerTriangular(ones(3,3))^i) * ones(3)
    end
end


text:
module MyExample

# Write your package code here.

include("extra_file.jl")

export my_f

end


text:
using Flux: @treelike

"""
    TransformerModel(embed::AbstractEmbed, transformers::AbstractTransformer)
    TransformerModel(
                      embed::AbstractEmbed,
                      transformers::AbstractTransformer,
                      classifier
                     )

a structure for putting everything together
"""
struct TransformerModel{E <: AbstractEmbed, T <: AbstractTransformer, C}
    embed::E
    transformers::T
    classifier::C
end

TransformerModel(embed, transformers) = TransformerModel(embed, transformers, identity)

@treelike TransformerModel

"""
    set_classifier(model::TransformerModel, classifier)

return a new TransformerModel whose classifier is set to `classifier`.
"""
set_classifier(model::TransformerModel, classifier) = TransformerModel(model.embed, model.transformers, classifier)

"""
    clear_classifier(model::TransformerModel)

return a new TransformerModel without classifier.
"""
clear_classifier(model::TransformerModel) = TransformerModel(model.embed, model.transformers, identity)

Base.print(io::IO, model::TransformerModel) = Base.summary(io, model)

function recursive_print(io::IO, x, i=1; tabchar = '\t')
    print(io, x)
    print(io, '\n')
end

function recursive_print(io::IO, x::AbstractArray{T}, i=1; tabchar = '\t') where T
    if T == Any || T <: Number
        recursive_print(io, typeof(x), i; tabchar = tabchar)
    else
        recursive_print(io, Tuple(x), i; tabchar = tabchar)
    end
end

function recursive_print(io::IO, x::Union{Tuple, NamedTuple}, i=1; tabchar = '\t')
    print(io, "(\n$(tabchar)")
    for k ∈ keys(x)
        print(io, tabchar^i)
        print(io, "$k => ")
        recursive_print(io, x[k], i+1; tabchar = tabchar)
        print(io, "$(tabchar)")
    end
    print(io, tabchar^(i-1))
    print(io, ")\n")
end

function Base.show(io::IO, model::TransformerModel)
    tabchar = "  "
    print(io, "TransformerModel{")
    print(io, typeof(model.transformers))
    print(io, "}(\n$(tabchar)")
    print(io, "embed = ")
    print(io, model.embed)
    print(io, ",\n$(tabchar)")
    print(io, "transformers = ")
    print(io, model.transformers)
    if model.classifier !== identity
        print(io, ",\n$(tabchar)")
        print(io, "classifier = \n$(tabchar^2)")
        recursive_print(io, model.classifier, 2; tabchar = tabchar)
    else
        print(io, '\n')
    end
    print(io, ")")
end


text:
struct Events{T <: EpidemicModel}
  exposure::Union{Nothing, Vector{Float64}}
  infection::Vector{Float64}
  removal::Union{Nothing, Vector{Float64}}
  individuals::Int64

  function Events{T}(e::V, i::V, r::V) where {T <: SEIR, V <: Vector{Float64}}
    if length(unique((length.([e; i; r])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(e, i, r, length(i))
  end

  function Events{T}(n::Integer) where T <: SEIR
    return new{T}(fill(NaN, n), fill(NaN, n), fill(NaN, n), n)
  end

  function Events{T}(e::V, i::V) where {T <: SEI, V <: Vector{Float64}}
    if length(unique((length.([e; i])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(e, i, nothing, length(i))
  end

  function Events{T}(n::Integer) where T <: SEI
    return new{T}(fill(NaN, n), fill(NaN, n), nothing, n)
  end

  function Events{T}(i::V, r::V) where {T <: SIR, V <: Vector{Float64}}
    if length(unique((length.([i; r])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(nothing, i, r, length(i))
  end

  function Events{T}(n::Integer) where T <: SIR
    return new{T}(nothing, fill(NaN, n), fill(NaN, n), n)
  end

  function Events{T}(i::V) where {T <: SI, V <: Vector{Float64}}
    return new{T}(nothing, i, nothing, length(i))
  end

  function Events{T}(n::Integer) where T <: SI
    return new{T}(nothing, fill(NaN, n), nothing, n)
  end
end

function Events{T}(a::Array{Float64,2}) where T <: EpidemicModel
  if size(a, 2) == 3
    return Events{T}(a[:,1], a[:,2], a[:,3])
  elseif size(a, 2) == 2
    return Events{T}(a[:,1], a[:,2])
  elseif size(a, 2) == 1
    return Events{T}(a[:,1])
  else
    @error "Invalid array size for construction of an $(Events{T}) object"
  end
end

function Events{T}(x::Vector{DiseaseState}) where T <: EpidemicModel
  events = Events{T}(length(x))
  for i = 1:length(x)
    if x[i] in _state_progressions[T]
      for j = _state_progressions[T][2:findfirst(Ref(x[i]) .== _state_progressions[T])]
        events[j][i] = -Inf
      end
    else
      @error "Invalid initial state for individual $i"
    end
  end
  return events
end

function Base.show(io::IO, x::Events{T}) where T <: EpidemicModel
  return print(io, "$T model event times (n=$(x.individuals))")
end

function Base.getindex(x::Events{T}, new_state::DiseaseState) where T <: EpidemicModel
  if new_state == State_E
    return x.exposure
  elseif new_state == State_I
    return x.infection
  elseif new_state == State_R
    return x.removal
  else
    @error "Unrecognized indexing disease state"
  end
end

function Base.getindex(x::Events{T}, states::Vector{DiseaseState}) where T <: EpidemicModel
  y = x[states[1]]
  for i = 2:length(states)
    y = hcat(y, x[states[i]])
  end
  return y
end

function Base.convert(::Type{Array{Float64, 2}}, x::Events{T}) where T <: EpidemicModel
  return x[_state_progressions[T][2:end]]
end

function Base.convert(::Type{Vector{Float64}}, x::Events{T}) where T <: EpidemicModel
  return x[_state_progressions[T][2:end]][:]
end

function Base.convert(::Type{Array{Float64, 2}}, x::Array{Events{T}, 1}) where T <: EpidemicModel
  y = convert(Vector{Float64}, x[1])'
  for i = 2:length(x)
    y = vcat(y, convert(Vector{Float64}, x[i])')
  end
  return y
end

function Base.minimum(x::Events{T}) where T <: EpidemicModel
  y = convert(Array{Float64, 2}, x)
  return minimum(y[.!isnan.(y)])
end

function Base.maximum(x::Events{T}) where T <: EpidemicModel
  y = convert(Array{Float64, 2}, x)
  return maximum(y[.!isnan.(y)])
end

function Statistics.mean(x::Vector{Events{T}}) where T <: EpidemicModel
  return Events{T}(mean([convert(Array{Float64, 2}, i) for i in x]))
end


text:
export MonteCarloTreeSearch, MCTS, MCTSTreeSolution
export MCTSNodeSelector, MaxUCBSelector, BoltzmannUCBSelector

## MCTS solution type ##

"MCTS policy solution."
@kwdef mutable struct MCTSTreeSolution <: PolicySolution
	Q::Dict{UInt64,Dict{Term,Float64}} = Dict()
	s_visits::Dict{UInt64,Int} = Dict()
	a_visits::Dict{UInt64,Dict{Term,Float64}} = Dict()
end

get_action(sol::MCTSTreeSolution, state::State) =
    best_action(sol, state)
best_action(sol::MCTSTreeSolution, state::State) =
    argmax(sol.Q[hash(state)])
rand_action(sol::MCTSTreeSolution, state::State) =
    best_action(sol, state)
get_value(sol::MCTSTreeSolution, state::State) =
    maximum(sol.Q[hash(state)])
get_value(sol::MCTSTreeSolution, state::State, action::Term) =
    sol.Q[hash(state)][action]
get_action_values(sol::MCTSTreeSolution, state::State) =
    pairs(sol.Q[hash(state)])

has_state_node(sol::MCTSTreeSolution, state::State) =
	hash(state) in keys(sol.s_visits)

## Selection strategies for leaf nodes ##

"Max upper-confidence bound (UCB) action policy."
struct MaxUCBPolicy <: PolicySolution
	tree::MCTSTreeSolution
	confidence::Float64
end

function get_action_values(sol::MaxUCBPolicy, state::State)
	state_id = hash(state)
	s_visits = sol.tree.s_visits[state_id]
	vals = map(collect(sol.tree.Q[state_id])) do (act, val)
		a_visits = sol.tree.a_visits[state_id][act]
		if s_visits != 0 && !(s_visits == 1 && a_visits == 0)
			val += sol.confidence * sqrt(log(s_visits) / a_visits)
		end
		return act => val
	end
	return vals
end

function best_action(sol::MaxUCBPolicy, state::State)
    a_vals = get_action_values(sol, state)
	best_idx = argmax(last.(a_vals))
	return first(a_vals[best_idx])
end

get_action(sol::MaxUCBPolicy, state::State) =
    best_action(sol, state)
rand_action(sol::MaxUCBPolicy, state::State) =
    best_action(sol, state)
get_value(sol::MaxUCBPolicy, state::State) =
    maximum(last, get_action_values(sol, state))
get_value(sol::MaxUCBPolicy, state::State, action::Term) =
    findfirst(x -> first(x) == action, get_action_values(sol, state))

BoltzmannUCBPolicy(tree::MCTSTreeSolution, confidence, temperature) =
	BoltzmannPolicy(MaxUCBPolicy(tree, confidence), temperature)

"Abstract type for MCTS node selection strategies."
abstract type MCTSNodeSelector end

(sel::MCTSNodeSelector)(tree::MCTSTreeSolution, domain::Domain, state::State) =
	error("Not implemented.")

"Max UCB selection strategy."
@kwdef struct MaxUCBSelector <: MCTSNodeSelector
	confidence::Float64 = 2.0
end

(sel::MaxUCBSelector)(tree::MCTSTreeSolution, ::Domain, state::State) =
 	get_action(MaxUCBPolicy(tree, sel.confidence), state)

"Boltzmann UCB selection strategy."
@kwdef struct BoltzmannUCBSelector <: MCTSNodeSelector
	confidence::Float64 = 2.0
	temperature::Float64 = 1.0
end

(sel::BoltzmannUCBSelector)(tree::MCTSTreeSolution, ::Domain, state::State) =
 	get_action(BoltzmannUCBPolicy(tree, sel.confidence, sel.temperature), state)

## Leaf node estimators ##

"Abstract type for MCTS leaf value estimators."
abstract type MCTSLeafEstimator end

(e::MCTSLeafEstimator)(::Domain, ::State, ::Specification, depth::Int) =
	error("Not implemented.")

"Estimates value as a constant."
struct ConstantEstimator{C <: Real} <: MCTSLeafEstimator
	value::C
end

(e::ConstantEstimator)(::Domain, ::State, ::Specification, ::Int) = e.value

"Estimates value via uniform random rollouts."
struct RandomRolloutEstimator <: MCTSLeafEstimator end

function (e::RandomRolloutEstimator)(domain::Domain, state::State,
								     spec::Specification, depth::Int)
	sim = RewardAccumulator(depth)
	policy = RandomPolicy(domain)
	return sim(policy, domain, state, spec)
end

"Estimates value via policy rollouts."
struct PolicyRolloutEstimator{P <: PolicySolution} <: MCTSLeafEstimator
	policy::P
end

function (e::PolicyRolloutEstimator)(domain::Domain, state::State,
									 spec::Specification, depth::Int)
	sim = RewardAccumulator(depth)
	return sim(e.policy, domain, state, spec)
end

## Main algorithm ##

"Planner that uses Monte Carlo Tree Search (MCTS)."
@kwdef mutable struct MonteCarloTreeSearch{S,E} <: Planner
	n_rollouts::Int64 = 50
	max_depth::Int64 = 50
	heuristic::Heuristic = NullHeuristic() # Initial value heuristic
	selector::S = BoltzmannUCBSelector() # Node selection strategy
	estimator::E = RandomRolloutEstimator() # Leaf node value estimator
end

const MCTS = MonteCarloTreeSearch

function solve(planner::MonteCarloTreeSearch,
			   domain::Domain, state::State, spec::Specification)
	@unpack n_rollouts, max_depth = planner
	@unpack heuristic, selector, estimator = planner
	discount = get_discount(spec)
	# Precompute heuristic information
    precompute!(heuristic, domain, state, spec)
    # Initialize solution
	sol = MCTSTreeSolution()
	sol = insert_node!(planner, sol, domain, state, spec)
	# Perform rollouts from initial state
	a_visited, s_visited = Term[], State[]
	initial_state = state
	for n in 1:n_rollouts
		state = initial_state
		value = 0.0
        # Rollout until maximum depth
        for t in 1:max_depth
			# Terminate if rollout reaches goal
            if is_goal(spec, domain, state) break end
			# Select action
			act = selector(sol, domain, state)
			push!(s_visited, state)
			push!(a_visited, act)
			# Transition to next state
            state = transition(domain, state, act)
			# Insert leaf node and evaluate
			if !has_state_node(sol, state)
				insert_node!(planner, sol, domain, state, spec)
				value = get_discount(spec)^t *
						estimator(domain, state, spec, max_depth-t)
				break
			end
        end
        # Backpropagate value estimates
		next_state = state
        while length(s_visited) > 0
            state, act = pop!(s_visited), pop!(a_visited)
			state_id = hash(state)
			# Update visitation counts
			sol.s_visits[state_id] += 1
			sol.a_visits[state_id][act] += 1
			# Compute value (accumulated discounted reward)
			reward = get_reward(spec, domain, state, act, next_state)
			value = get_discount(spec) * value + reward
			# Take weighted average with existing Q value
			sol.Q[state_id][act] +=
				(value-sol.Q[state_id][act]) / sol.a_visits[state_id][act]
			state = next_state
        end
	end
	return sol
end

function insert_node!(planner::MonteCarloTreeSearch, sol::MCTSTreeSolution,
			   		  domain::Domain, state::State, spec::Specification)
	actions = available(state, domain)
	qs = map(actions) do act
        next_state = transition(domain, state, act)
        r = get_reward(spec, domain, state, act, next_state)
        h_val = planner.heuristic(domain, next_state, spec)
        return get_discount(spec) * (-h_val) + r
    end
	state_id = hash(state)
    sol.Q[state_id] = Dict{Term,Float64}(zip(actions, qs))
	sol.a_visits[state_id] = Dict{Term,Int}(a => 0 for a in actions)
	sol.s_visits[state_id] = 0
	return sol
end


text:
using HMMs
using Brobdingnag
using PyPlot
using LinearAlgebra

m = 3
graph1 = HMMs.fullyconnected( m )

hmm1 = HMMs.randomhmm( graph1, seed=1 )
y1 = HMMs.rand( hmm1, 100_000 );

hmm2 = HMMs.randomhmm( graph1, calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y1 );
HMMs.em( hmm2, debug=2, keepintermediates=true )
# 631 iterations, 232 seconds, exp(117683.72053268638) likelihood
hmm2a = HMMs.randomhmm( graph1, calc=Brob, seed=2 )
HMMs.setobservations( hmm2a, y1 );
HMMs.em( hmm2a, debug=2, keepintermediates=true, acceleration=10 )
# 296 iterations, 108 seconds, exp(117681.20996505348) likelihood

C = convert( Matrix{Float64}, HMMs.sandwich( hmm2a ) )
convert(Matrix{Float64}, [HMMs.getparameters( hmm2a ) sqrt.(diag(C))] )
convert(Vector{Float64}, HMMs.getparameters( hmm2a )./sqrt.(diag(C)) )

graph3 = HMMs.Digraph( [1,1,2,2,2,3,3], [1,2,1,2,3,2,3] )
hmm3 = HMMs.randomhmm( graph3, seed=1 )
y3 = rand( hmm3, 100_000 );

hmm4 = HMMs.randomhmm( graph1, calc=Brob, seed=2 )
HMMs.setobservations( hmm4, y3 );
HMMs.em( hmm4, debug=2, keepintermediates=true )
# 7406 iterations, 2792 seconds

hmm4a = HMMs.randomhmm( graph1, calc=Brob, seed=2 )
HMMs.setobservations( hmm4a, y3 );
HMMs.em( hmm4a, debug=2, keepintermediates=true, acceleration=10 )
# 960 iterations, 394 seconds, exp(-17394.489088597584) likelihood
# old version: 4512 iterations, 1782 seconds, exp(-17394.488432664482) likelihood

C = convert( Matrix{Float64}, HMMs.sandwich( hmm4a ) )
convert(Matrix{Float64}, [HMMs.getparameters( hmm4a ) sqrt.(diag(C))] )

hmm4a3 = HMMs.randomhmm( graph1, calc=Brob, seed=3 )
HMMs.setobservations( hmm4a3, y3 );
HMMs.em( hmm4a3, debug=2, keepintermediates=true, acceleration=10 )
# 551 iterations, 223 seconds, exp(-17394.490665902726) likelihood

C = convert( Matrix{Float64}, HMMs.sandwich( hmm4a3 ) )
convert(Vector{Float64}, [HMMs.getparameters( hmm4a3 ); sqrt.(diag(C))] )


text:
export AbstractLearner

"""
A learner is used to define:

- How to generate necessary training data?
- How to update the inner [Approximators](@ref)?
"""
abstract type AbstractLearner end

approximator(x::AbstractLearner) = x.approximator

(learner::AbstractLearner)(obs::Observation) = approximator(learner)(get_state(obs))
(learner::AbstractLearner)(obs::Observation, a) = approximator(learner)(get_state(obs), a)
(learner::AbstractLearner)(s) = approximator(learner)(s)
(learner::AbstractLearner)(s, a) = approximator(learner)(s, a)


text:
function add_view_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(view_pole_cb,b["select_pole_button"],"clicked",Void,(),false,(handles,))
    signal_connect(view_discrete_cb,b["view_discrete_points"],"clicked",Void,(),false,(handles,))
    signal_connect(view_whiskers_cb,b["view_tracked"],"clicked",Void,(),false,(handles,))

    nothing
end

view_pole(b::Gtk.GtkBuilder)=getproperty(b["select_pole_button"],:active,Bool)
view_discrete(b::Gtk.GtkBuilder)=getproperty(b["view_discrete_points"],:active,Bool)
show_tracked(b::Gtk.GtkBuilder)=getproperty(b["view_tracked"],:active,Bool)

function view_pole_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data
    redraw_all(han)
    nothing
end

function view_discrete_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data
    redraw_all(han)
    nothing
end

function view_whiskers_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data
    redraw_all(han)
    nothing
end


text:
using DataFrames
using CSV

function unrate_full(path)
    df = DataFrame(CSV.File(path*"unrate_full.csv", header=0))
    #data = convert(Matrix, df[:, 1])
    g = x -> parse(Float64, x)
    target = convert(Array{}, g.(df[2:end, 2]))
    data = collect(1:1:length(target))
    return data, target
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
A unified interface for log density problems, for

1. defining mappings to a log density (eg Bayesian for inference),

2. optionally obtaining a gradient using automatic differentiation,

3. defining a common interface for working with such log densities and gradients (eg MCMC).

These use cases utilize different parts of this package, make sure you read the
documentation.
"""
module LogDensityProblems

export TransformedLogDensity, ADgradient

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES, TYPEDEF
using UnPack: @unpack
import Random
using Requires: @require
import TransformVariables

####
#### interface for problems
####

"""
$(TYPEDEF)

A trait that means that a log density supports evaluating derivatives up to order `K`.

Typical values for `K` are `0` (just the log density) and `1` (log density and gradient).
"""
struct LogDensityOrder{K}
    function LogDensityOrder{K}() where K
        _K = Int(K)
        @argcheck _K ≥ 0
        new{_K}()
    end
end

LogDensityOrder(K::Integer) = LogDensityOrder{K}()

Base.isless(::LogDensityOrder{A}, ::LogDensityOrder{B}) where {A, B} = A ≤ B

"""
$(SIGNATURES)

Test if the type (or a value, for convenience) supports the log density interface.

When `nothing` is returned, it doesn't support this interface.  When `LogDensityOrder{K}()`
is returned (typically with `K == 0` or `K = 1`), derivatives up to order `K` are supported.
*All other return values are invalid*.

# Interface description

The following methods need to be implemented for the interface:

1. [`dimension`](@ref) returns the *dimension* of the domain,

2. [`logdensity`](@ref) evaluates the log density at a given point.

3. [`logdensity_and_gradient`](@ref) when `K ≥ 1`.

See also [`LogDensityProblems.stresstest`](@ref) for stress testing.
"""
capabilities(T::Type) = nothing

capabilities(x) = capabilities(typeof(x)) # convenience function

"""
    dimension(ℓ)

Dimension of the input vectors `x` for log density `ℓ`. See [`logdensity`](@ref),
[`logdensity_and_gradient`](@ref).

!!! note
    This function is *distinct* from `TransformedVariables.dimension`.
"""
function dimension end

"""
    logdensity(ℓ, x)

Evaluate the log density `ℓ` at `x`, which has length compatible with its
[`dimension`](@ref).

Return a real number, which may or may not be finite (can also be `NaN`). Non-finite values
other than `-Inf` are invalid but do not error, caller should deal with these appropriately.

# Note about constants

Log densities can be shifted by *the same constant*, as long as it is consistent between
calls. For example,

```julia
logdensity(::StandardMultivariateNormal) = -0.5 * sum(abs2, x)
```

is a valid implementation for some callable `StandardMultivariateNormal` that would
implement the standard multivariate normal distribution (dimension ``k``) with pdf
```math
(2\\pi)^{-k/2} e^{-x'x/2}
```
"""
function logdensity end

"""
    logdensity_and_gradient(ℓ, x)

Evaluate the log density `ℓ` and its gradient at `x`, which has length compatible with its
[`dimension`](@ref).

Return two values:

- the log density as real number, which equivalent to `logdensity(ℓ, x)`

- *if* the log density is finite, the gradient, a vector of real numbers, otherwise this
   value is arbitrary and should be ignored.

The first argument (the log density) can be shifted by a constant, see the note for
[`logdensity`](@ref).
"""
function logdensity_and_gradient end

#####
##### Transformed log density (typically Bayesian inference)
#####

"""
    TransformedLogDensity(transformation, log_density_function)

A problem in Bayesian inference. Vectors of length compatible with the dimension (obtained
from `transformation`) are transformed into a general object `θ` (unrestricted type, but a
named tuple is recommended for clean code), correcting for the log Jacobian determinant of
the transformation.

`log_density_function(θ)` is expected to return *real numbers*. For zero densities or
infeasible `θ`s, `-Inf` or similar should be returned, but for efficiency of inference most
methods recommend using `transformation` to avoid this. It is recommended that
`log_density_function` is a callable object that also encapsulates the data for the problem.

Use the property accessors `ℓ.transformation` and `ℓ.log_density_function` to access the
arguments of `ℓ::TransformedLogDensity`, these are part of the public API.

# Usage note

This is the most convenient way to define log densities, as `capabilities`, `logdensity`,
and `dimension` are automatically defined. To obtain a type that supports derivatives, use
[`ADgradient`](@ref).
"""
struct TransformedLogDensity{T <: TransformVariables.AbstractTransform, L}
    transformation::T
    log_density_function::L
end

function Base.show(io::IO, ℓ::TransformedLogDensity)
    print(io, "TransformedLogDensity of dimension $(dimension(ℓ))")
end

capabilities(::Type{<:TransformedLogDensity}) = LogDensityOrder{0}()

dimension(p::TransformedLogDensity) = TransformVariables.dimension(p.transformation)

function logdensity(p::TransformedLogDensity, x::AbstractVector)
    @unpack transformation, log_density_function = p
    TransformVariables.transform_logdensity(transformation, log_density_function, x)
end

#####
##### AD wrappers --- interface and generic code
#####

"""
An abstract type that wraps another log density for calculating the gradient via AD.

Automatically defines the methods `capabilities`, `dimension`, and `logdensity` forwarding
to the field `ℓ`, subtypes should define a [`logdensity_and_gradientent`](@ref).

This is an implementation helper, not part of the API.
"""
abstract type ADGradientWrapper end

logdensity(ℓ::ADGradientWrapper, x::AbstractVector) = logdensity(ℓ.ℓ, x)

capabilities(::Type{<:ADGradientWrapper}) = LogDensityOrder{1}()

dimension(ℓ::ADGradientWrapper) = dimension(ℓ.ℓ)

Base.parent(ℓ::ADGradientWrapper) = ℓ.ℓ

"""
$(SIGNATURES)

Wrap `P` using automatic differentiation to obtain a gradient.

`kind` is usually a `Val` type with a symbol that refers to a package, for example
```julia
ADgradient(Val(:ForwardDiff), P)
```

The symbol can also be used directly as eg

```julia
ADgradient(:ForwardDiff, P)
```

and should mostly be equivalent if the compiler manages to fold the constant.

`parent` can be used to retrieve the original argument.

See `methods(ADgradient)`. Note that **some methods are defined conditionally on the
relevant package being loaded.**
"""
ADgradient(kind::Symbol, P; kwargs...) = ADgradient(Val{kind}(), P; kwargs...)

function ADgradient(v::Val{kind}, P; kwargs...) where kind
    @info "Don't know how to AD with $(kind), consider `import $(kind)` if there is such a package."
    throw(MethodError(ADgradient, (v, P)))
end

####
#### AD wrappers - specific
####

function __init__()
    @require DiffResults="163ba53b-c6d8-5494-b064-1a9d43ac40c5" include("DiffResults_helpers.jl")
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("AD_ForwardDiff.jl")
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("AD_Tracker.jl")
    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include("AD_Zygote.jl")
end

####
#### stress testing
####

"""
$(SIGNATURES)

Test `ℓ` with random values.

`N` random vectors are drawn from a standard multivariate Cauchy distribution, scaled with
`scale` (which can be a scalar or a conformable vector).

Each random vector is then used as an argument in `f(ℓ, ...)`. `logdensity` and
`logdensity_and_gradient` are  recommended for `f`.

In case the call produces an error, the value is recorded as a failure, which are returned
by the function.

Not exported, but part of the API.
"""
function stresstest(f, ℓ; N = 1000, rng::Random.AbstractRNG = Random.GLOBAL_RNG, scale = 1)
    failures = Vector{Float64}[]
    d = dimension(ℓ)
    for _ in 1:N
        x = TransformVariables.random_reals(d; scale = scale, cauchy = true, rng = rng)
        try
            f(ℓ, x)
        catch e
            push!(failures, x)
        end
    end
    failures
end

end # module


text:
using Instruments, Currencies

using Test

@testset "All currencies" begin
    for sym in Currencies.allsymbols()
        ccy = Currency{sym}()
        # This needs some new tests, since Cash is not yet defined
    end
end


text:
using ControlSystems;

A=[0 1;1 0];
B=[0;1];
Q=[2 2;2 2];
R=8;

l=lqr(A,B,Q,R)

care(A,B,Q,R)

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design


text:
push!(LOAD_PATH, pwd())
include("startup.jl")

using HDF5
#using PyCall
#pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
using PyPlot, Colors
import .InterpolationFunctions
import  .MasslistFunctions
import .ResultFileFunctions
using Statistics
using DataFrames
using CSV

#fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/DryAtmo/results/"
#fn = "_result_withDMS.hdf5"

stagesfile = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/DMS_TME_BAG_List.csv"		# add filename, when plotstages = true

ramping = "O3"
plotSimulatedData = true

if ramping == "DMS"
	fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/DMSramp/results/"
	simulationfile = string(fp, "combiningAllOutputs_DMSramp_O314e12_TME3e10.csv")
	plotHighTimeRes = true # Plot every datapoint or only file averages
elseif ramping == "RO2"
	fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/RO2ramp/results/"
	simulationfile = string(fp, "combiningAllOutputs_TMEramp_O3_14e12_DMS_2e12.csv")
	plotHighTimeRes = true # Plot every datapoint or only file averages
elseif ramping == "O3"
	fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/O3ramp/results/"
	simulationfile = string(fp, "combiningAllOutputs_O3ramp_DMS1e12_TME3e10.csv")
	plotHighTimeRes = false # only file averages
end

fn = "_result.hdf5"
file = string(fp,fn)

plotFittedInsteadOfSummed = true # Use multi peak fitted data instead of raw
smoothing = 100 # Average n samples, 1 for raw
plotsymbol = "-"
isobarToPlot = 0

backgroundSubstractionMode = 2 # 0=no correction, 1=minimum found in all data,  2=mean of range given below
backgroundStart = DateTime(2020,5,28,8,02)
backgroundEnd = DateTime(2020,5,28,8,20)

calibrate = false
calibrationfactor = 4.5 * 2.47e7	# cps / molec/cm³

normCPS = true

massrange2plot = false
plotWronglyAssigned = false

include("manualMassLibrary.jl")

massesToPlot = [
# Examples for selecting what to plot:

#H3O[1]
#H3OH2O[1]
#H3OH2OH2O[1]
#H3OH2OH2OH2O[1]
#NH3H2O[1]
#NH3[1]
#NH3NH3[1]
#O2[1]
#NO[1]

#DMS[1]
#DMSO[1]
#MSIA[1]
#MasslistFunctions.massFromComposition(C=2,H=4,S=1,O=3)
#MasslistFunctions.massFromComposition(C=2,H=8,S=1,O=2)
#MasslistFunctions.massFromComposition(C=2,H=10,S=1,O=3)
#MasslistFunctions.massFromComposition(C=2,H=4,S=1,O=1)
#MasslistFunctions.massFromComposition(C=1, S=1,O=1, Hplus = 0)
#MasslistFunctions.massFromComposition(C=1,H=2,O=1)
# MasslistFunctions.massFromComposition(C=3, O=3, H=5, Hplus=1)		#non existent
#MasslistFunctions.massFromComposition(C=6, O=3, H=13, Hplus=1)
#MasslistFunctions.massFromComposition(C=2, O=2, H=5, S=1, Hplus=1)
#MasslistFunctions.massFromComposition(C=2, O=4, H=5, S=1, Hplus=1)
#MasslistFunctions.massFromComposition(C=2, O=3, H=7, S=1, Hplus=1)
#MasslistFunctions.createMassList(; C=4:5, O=2:6, N=0:0, S=1:2, nHplus=1, allowRadicals=false)[1]
# MasslistFunctions.createMassList(; C=3, S=1:2, nHplus=1, allowRadicals=false)[1]		#non existent
# MasslistFunctions.createMassList(; C=8, S=1:2,O=3:5, nHplus=1, allowRadicals=false)[1]		#non existent
#MasslistFunctions.massFromComposition(C=9, H=18, O=4, Hplus=1)
#MasslistFunctions.massFromComposition(C=6, H=10, O=4, Hplus=1)
# MasslistFunctions.massFromComposition(C=12, H=26, O=4, Hplus=1)		#non existent

# DMS + TME radicals
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=2)
#MasslistFunctions.massFromComposition(C=1,H=3,O=2)
#MasslistFunctions.massFromComposition(C=3,H=5,O=3)
  #MasslistFunctions.massFromComposition(C=6,H=13,O=3)
#MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=2)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=3)
MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=4)
MasslistFunctions.massFromComposition(C=2,H=7,S=1,O=3)
MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=4)

#MasslistFunctions.massFromComposition(C=4,H=10,S=2,O=2)

MasslistFunctions.massFromComposition(C=3,H=8,S=1,O=4)
MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=4)
#MasslistFunctions.massFromComposition(C=2,H=8,O=3)

# DMS + TME dimers
#=
  #MasslistFunctions.massFromComposition(C=2,H=6,O=2)
MasslistFunctions.massFromComposition(C=4,H=8,O=3)
MasslistFunctions.massFromComposition(C=6,H=10,O=4)
  #MasslistFunctions.massFromComposition(C=7,H=16,O=3)
  #MasslistFunctions.massFromComposition(C=9,H=18,O=4)
  #MasslistFunctions.massFromComposition(C=12,H=26,O=4)
MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=2)
  #MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=3)
MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=4)

  #MasslistFunctions.massFromComposition(C=3,H=8,S=1,O=2)
MasslistFunctions.massFromComposition(C=3,H=8,S=1,O=4)
  #MasslistFunctions.massFromComposition(C=3,H=10,S=1,O=3)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=3)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=4)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=5)
MasslistFunctions.massFromComposition(C=5,H=10,S=1,O=3)
MasslistFunctions.massFromComposition(C=5,H=10,S=1,O=5)
MasslistFunctions.massFromComposition(C=5,H=12,S=1,O=4)
MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=3)
MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=4)
  #MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=5)
MasslistFunctions.massFromComposition(C=8,H=18,S=1,O=3)
  #MasslistFunctions.massFromComposition(C=8,H=18,S=1,O=5)
MasslistFunctions.massFromComposition(C=8,H=20,S=1,O=4)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=2)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=3)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=4)
MasslistFunctions.massFromComposition(C=4,H=10,S=2,O=2)
MasslistFunctions.massFromComposition(C=4,H=10,S=2,O=4)
MasslistFunctions.massFromComposition(C=4,H=12,S=2,O=3)
=#
]

if massrange2plot
	result4masses = ResultFileFunctions.loadResults(file; useAveragesOnly= true, masslistOnly = true)
	massesToPlot = result4masses.MasslistMasses[59.8 .< result4masses.MasslistMasses .< 60.0]
end
measResult = ResultFileFunctions.loadResults(file, massesToLoad=massesToPlot, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed, massMatchTolerance=0.001)

if (isobarToPlot != 0)
  	isobarResult = ResultFileFunctions.loadResults(file, massesToLoad=[isobarToPlot+0.3], massMatchTolerance=0.5, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed)
  	measResult=ResultFileFunctions.joinResultsMasses(measResult, isobarResult)
end

if (backgroundSubstractionMode == 0)
  background=0
elseif (backgroundSubstractionMode == 1)
  background = minimum(InterpolationFunctions.averageSamples(measResult.Traces,smoothing),dims=1)
elseif backgroundSubstractionMode == 2
  import Statistics
  background = Statistics.mean(measResult.Traces[(measResult.Times.>backgroundStart) .& (measResult.Times.<backgroundEnd),:],dims=1)
end
bgCorrectedTraces = measResult.Traces .- background

stagestimes, stagesdata = ImportFunctions.importFromCSV(stagesfile; headerrow = 9)
meanTraces = DatasetFunctions.calculateMeanTraces(stagestimes, measResult)

showThis = meanTraces[:,1] .!= 0.0

if ramping == "DMS"
	if plotSimulatedData
		simResults = DataFrame(CSV.File(simulationfile, delim=' '))
		DMSinit = simResults.DMS[1:2:end]
		simResults = simResults[2:2:end,:]
	end
	showThis[end-1:end] .= false
	tme = median(stagesdata.TME_percm3[showThis])
	o3 = median(stagesdata.O3_percm3[showThis])

	for name in names(meanTraces)
		xdata = stagesdata.DMS_percm3[showThis]
		colordata = stagesdata.Uextr_V[showThis]
		if name == "COS"
			ydata = meanTraces[showThis, name]*2.47e10/170
		elseif name == "CH2O"
			ydata = meanTraces[showThis, name]*2.47e10/700
		else
			ydata = meanTraces[showThis, name]*2.47e10/4500
			ydata[colordata .> 70] = ydata[colordata .> 70]*2
		end
		fig=figure()
		ax = subplot(111)
		ax.scatter(xdata, ydata, c = colordata, cmap = "rainbow")
		ax.set_yscale("log")
		ax.set_ylabel("estimated concentration [molec / cm³]")
		ax.set_xscale("log")
		ax.set_xlim(minimum(xdata[xdata .> 0])/2, maximum(xdata)*2)
		ax.set_ylim(minimum(ydata[ydata .> 0])/2, maximum(ydata)*2)
		title("$(name) - col. by Extr. Voltage, TME = $(tme), O3 = $(o3)")
		savefig("$(fp)_$(name)_DMSramp.png")
	end

	fig = figure()
	ax = subplot(111)
	showThis2 = copy(showThis)
	a = findall(showThis .== true)
	b = findall(stagesdata.Uextr_V .> 70)
	showThis2[a[findall(x -> x in b, a)]] .= false
	xdata = stagesdata.DMS_percm3[showThis2]
	colordata = stagesdata.Uextr_V[showThis2]
	for name in names(meanTraces)
		if occursin("C", name) && (name != "C2H8O2S") && (name != "C2H6S")
			ydata = meanTraces[showThis2, name]*2.47e10/4500
			ydata[7] = NaN
			ydata[stagesdata.DMS_percm3[showThis2] .== 0] .= NaN
			ax.plot(xdata, ydata, "o", markersize = 10, label = name)
		end
	end
	ax.set_yscale("log")
	ax.set_xlabel("DMS [molec / cm³]")
	ax.set_ylabel("estimated concentration [molec / cm³]")
	ax.set_xscale("log")
	ax.set_xlim(minimum(xdata[xdata .> 0])/2, maximum(xdata)*2)
	ax.set_ylim(1e3, 3e10)
	title("DMS ramp: soft mode, TME = $(tme), O3 = $(o3)")
	ax.legend(loc = 3)
	if plotSimulatedData
		ax2 = ax.twinx()
		ax2.semilogy(DMSinit, simResults[:,"HCHO"], "P:", color = "blue", label = "HCHO")
		ax2.semilogy(DMSinit, simResults[:,"OCS"], "P:", color = "orange", label = "OCS")
		ax2.semilogy(DMSinit, simResults[:,"CH3SCHO"], "P:", color = "green", label = "CH3SCHO")
		ax2.semilogy(DMSinit, simResults[:,"DMSO"], "P:", color = "red", label = "DMSO")
		ax2.semilogy(DMSinit, simResults[:,"MSIA"], "P:", color = "purple", label = "MSIA")
		ax2.semilogy(DMSinit, simResults[:,"HPMTF"], "P:", color = "brown", label = "HPMTF")
		ax2.legend(loc = 4)
		ax2.set_ylim(1e3, 3e10)
		ax2.set_ylabel("simulated concentrations")
	end
	savefig("$(fp)ramp_soft.png")

	fig = figure()
	ax = subplot(111)
	showThis3 = copy(showThis)
	a = findall(showThis .== true)
	b = findall(stagesdata.Uextr_V .< 70)
	showThis3[a[findall(x -> x in b, a)]] .= false
	xdata = stagesdata.DMS_percm3[showThis3]
	colordata = stagesdata.Uextr_V[showThis3]
	for name in names(meanTraces)
		if occursin("C", name) && (name != "C2H6S") && (name != "C2H8O2S")
			if name == "COS"
				ydata = meanTraces[showThis3, name]*2.47e10/170
			elseif name == "CH2O"
				ydata = meanTraces[showThis3, name]*2.47e10/700
			else
				ydata = meanTraces[showThis3, name]*2.47e10*2/4500
			end
			ydata[7] = NaN			# problem with 7th!
			ydata[stagesdata.DMS_percm3[showThis3] .== 0] .= NaN
			ax.plot(xdata, ydata,  "o", markersize = 10, label = name)
		end
	end
	ax.set_yscale("log")
	ax.set_xlabel("DMS [molec / cm³]")
	ax.set_ylabel("estimated concentration [molec / cm³]")
	ax.set_xscale("log")
	ax.set_xlim(minimum(xdata[xdata .> 0])/2, maximum(xdata)*2)
	ax.set_ylim(1e3, 3e10)
	title("DMS ramp: hard mode, TME = $(tme), O3 = $(o3)")
	ax.legend(loc = 3)
	if plotSimulatedData
		ax2 = ax.twinx()
		ax2.semilogy(DMSinit, simResults[:,"HCHO"], "P:", color = "blue", label = "HCHO")
		ax2.semilogy(DMSinit, simResults[:,"OCS"], "P:", color = "orange", label = "OCS")
		ax2.semilogy(DMSinit, simResults[:,"CH3SCHO"], "P:", color = "green", label = "CH3SCHO")
		ax2.semilogy(DMSinit, simResults[:,"DMSO"], "P:", color = "red", label = "DMSO")
		ax2.semilogy(DMSinit, simResults[:,"MSIA"], "P:", color = "purple", label = "MSIA")
		ax2.semilogy(DMSinit, simResults[:,"HPMTF"], "P:", color = "brown", label = "HPMTF")
		ax2.legend(loc = 4)
		ax2.set_ylim(1e3, 3e10)
		ax2.set_ylabel("simulated concentrations")
	end
	savefig("$(fp)ramp_hard.png")

elseif ramping == "RO2"

	if plotSimulatedData
		simResults = DataFrame(CSV.File(simulationfile, delim=','))
		DMSinit = simResults.DMS[1:2:end]
		TMEinit = simResults.DM23BU2ENE[1:2:end]
		O3init = simResults.O3[1:2:end]
		simResults = simResults[2:2:end,:]
	end

	# showThis[end-1:end] .= false
	dms = median(stagesdata.DMS_percm3[showThis])

	BGstep = 3

	fig = figure()
	ax = subplot(111)
	xdata = stagesdata.TME_percm3[showThis].*stagesdata.O3_percm3[showThis].*1e-15
	colordata = stagesdata.O3_percm3[showThis]
	for name in names(meanTraces)
		if occursin("C", name) && (name != "C2H8O2S") && (name != "C2H6S")
			ydata = (meanTraces[showThis, name] .- meanTraces[showThis, name][BGstep])*2.47e10/4500
			# ydata[7] = NaN
			ydata[stagesdata.TME_percm3[showThis] .== 0] .= NaN
			ax.plot(xdata, ydata, "o", markersize = 10, label = name)
		end
	end
	ax.set_yscale("log")
	ax.set_xlabel("TME * O3 * 1e-15")
	ax.set_ylabel("estimated concentration [molec / cm³]")
	ax.set_xscale("log")
	ax.set_xlim(minimum(xdata[xdata .> 0])/2, 4e13)
	ax.set_ylim(1e3, 3e10)
	title("RO2 ramp: soft mode, DMS = $(dms)")
	ax.legend(loc = 3)
	if plotSimulatedData
		ax2 = ax.twinx()
		ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"HCHO"]./1000, ":", color = "blue", label = "HCHO/1000")
		#ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"OCS"], ":", color = "orange", label = "OCS")
		ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"CH3SCHO"], ":", color = "orange", label = "CH3SCHO")
		ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"DMSO"]./10, ":", color = "green", label = "DMSO/10")
		ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"MSIA"]./10, ":", color = "red", label = "MSIA/10")
		ax2.semilogy(TMEinit.*O3init.*1e-15, simResults[:,"HPMTF"]./10, ":", color = "purple", label = "HPMTF/10")
		ax2.legend(loc = 4)
		ax2.set_ylim(1e3, 3e10)
		ax2.set_ylabel("simulated concentration")
	end
	savefig("$(fp)RO2ramp_soft_a.png")
elseif ramping == "O3"
	if plotSimulatedData
		simResults = DataFrame(CSV.File(simulationfile, delim=','))
		DMSinit = simResults.DMS[1:2:end]
		TMEinit = simResults.DM23BU2ENE[1:2:end]
		O3init = simResults.O3[1:2:end]
		simResults = simResults[2:2:end,:]
	end
	dms = median(stagesdata.DMS_percm3[showThis])
	tme = median(stagesdata.TME_percm3[showThis])

	BGstep = 1
	fig = figure()
	ax = subplot(111)
	xdata = stagesdata.O3_percm3[showThis]
	# colordata = stagesdata.O3_percm3[showThis]
	for name in names(meanTraces)
		if occursin("C", name) && (name != "C2H8O2S") && (name != "C2H6S")
			ydata = (meanTraces[showThis, name] .- meanTraces[showThis, name][BGstep])*2.47e10/4500
			# ydata[7] = NaN
			ydata[stagesdata.O3_percm3[showThis] .== 0] .= NaN
			ydata[stagesdata.Uextr[showThis] .!= 53] .= NaN
			ydata[stagesdata.Uextr[showThis] .!= 53] .= NaN
			ydata[stagesdata.EN_Vpp[showThis] .!= 0] .= NaN
			ax.plot(xdata, ydata, "o", markersize = 10, label = name)
		end
	end
	ax.set_yscale("log")
	ax.set_xlabel("O3")
	ax.set_ylabel("estimated concentration [molec / cm³]")
	ax.set_xscale("log")
	ax.set_xlim(minimum(xdata[xdata .> 0])/2, 1e13)
	ax.set_ylim(1e5, 3e10)
	title("O3 ramp: soft mode, DMS = $(dms), TME = $(tme)")
	ax.legend(loc = 3)
	if plotSimulatedData
		ax2 = ax.twinx()
		ax2.semilogy(O3init, simResults[:,"HCHO"]./100, ":", color = "blue", label = "HCHO/100")
		#ax2.semilogy(TMEinit.*O3init.*1e-12, simResults[:,"OCS"], ":", color = "orange", label = "OCS")
		ax2.semilogy(O3init, simResults[:,"CH3SCHO"], ":", color = "orange", label = "CH3SCHO")
		ax2.semilogy(O3init, simResults[:,"DMSO"], ":", color = "green", label = "DMSO")
		ax2.semilogy(O3init, simResults[:,"MSIA"], ":", color = "red", label = "MSIA")
		ax2.semilogy(O3init, simResults[:,"HPMTF"], ":", color = "purple", label = "HPMTF")
		ax2.legend(loc = 4)
		ax2.set_ylim(1e5, 3e10)
		ax2.set_ylabel("simulated concentration")
	end
	savefig("$(fp)O3ramp_soft.png")
end


text:
#=
if_while:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-08
=#

module ForWhile
    # methods
    function while_under_val(val)
        i = 1
        while i < val
            i += 1
        end
    end

    function for_val_times(val)
        j = 0
        for i = 1:val
            j += 1
        end
    end

    function main()
        println("While 100000000 times")
        @time while_under_val(100000001)
        println("")
        println("For 100000000 times")
        @time for_val_times(100000000)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .ForWhile
    ForWhile.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
Pkg.add(“DataFrames”)
Pkg.add("CSV")
Pkg.add("ScikitLearn")
using DataFrames, CSV, ScikitLearn: fit!, predict, @sk_import, fit_transform!
@sk_import preprocessing: LabelEncoder
 @sk_import model_selection: cross_val_score
 @sk_import metrics: accuracy_score
 @sk_import linear_model: LogisticRegression
 @sk_import ensemble: RandomForestClassifier
 @sk_import tree: DecisionTreeClassifier

##change directory to the folder containing test and train files
;cd /path to the train file
train = CSV.read("........train file..........") // readtable has been deprecated.
test  = CSV.read("........test file ..........")
describe(train) // gives us the "missing" values info ; showcols() has been deprecated
##scikitlearn uses PyCall.jl and pyhton doesn't support the data type "missing" , so we fill out the missing data
train[ismissing.train[://any column]), ://column] = mode(dropmissing(train[://any column]))    //isna and dropna didn't work for me

##Another problem I faced was, version mismatch of python.
##ScikitLearn.jl installed  and minimal Python 3.7 while I had my default python enivronment set to Pyhton 3 which was recognised by PyCall.
##You can fix this:
   In julia terminal type
                ENV["PYTHON"] = ""   //This will set the python environment to Julia's i.e. the one installed by ScikitLearn.jl
                restart julia and in the terminal type
                using PyCall
                PyCall.conda
                If it returns true you're good to go.

function classification_model(model, predictors)
     y = convert(Array, train[:13])
     X = convert(Array, train[predictors])
     X2 = convert(Array, test[predictors])

    #Fit the model:
     fit!(model, X, y)

     #Make predictions on training set:
     predictions = predict(model, X)

     #Print accuracy
     accuracy = accuracy_score(predictions, y)
     println("\naccuracy: ",accuracy)

     #5 fold cross validation
     cross_score = cross_val_score(model, X, y, cv=5)

     #print cross_val_score
     println("cross_validation_score: ", mean(cross_score))

     #return predictions
     fit!(model, X, y)
     global pred = predict(model, X2)  //declaring scope of pred as global so that it can be written into a csv file
     return pred
 end

model = //any
predictors =[:Gender, :Married, :Dependents, :Education,
 :Self_Employed, :Loan_Amount_Term, :Credit_History, :Property_Area,
 :LoanAmount]  //you may add any nuber of features
classification_model(model, predictors)

sample_submit=CSV.read("sample_submission.csv", copycols = true)

submission = [ test[: , 1] , pred]
CSV.write("submission.csv" , submission )

max_stars_repo_path: Bank Loan challange/Bank load code.jl

max_stars_repo_name: Sov-trotter/Bank-Loan-challange

max_issues_repo_path: Bank Loan challange/Bank load code.jl

max_issues_repo_name: Sov-trotter/Bank-Loan-challange

max_forks_repo_path: Bank Loan challange/Bank load code.jl

max_forks_repo_name: Sov-trotter/Bank-Loan-challange