text:
using Interpolations

@inline function calcpoint(blockIdx, blockDim, threadIdx, size)
    i = (blockIdx - 1) * blockDim + threadIdx
    return i, Float32(i)
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,1}, texture::CuDeviceTexture{<:Any,1})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    @inbounds dst[i] = texture[u]
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,2}, texture::CuDeviceTexture{<:Any,2})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    @inbounds dst[i,j] = texture[u,v]
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,3}, texture::CuDeviceTexture{<:Any,3})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    k, w = calcpoint(blockIdx().z, blockDim().z, threadIdx().z, size(dst)[3])
    @inbounds dst[i,j,k] = texture[u,v,w]
    return nothing
end

function fetch_all(texture)
    dims = size(texture)
    d_out = CuArray{eltype(texture)}(undef, dims...)

    kernel = @cuda launch=false kernel_texture_warp_native(d_out, texture)
    config = launch_configuration(kernel.fun)

    dim_x, dim_y, dim_z = size(texture, 1), size(texture, 2), size(texture, 3)
    threads_x = min(dim_x, config.threads)
    blocks_x = cld(dim_x, threads_x)

    kernel(d_out, texture; threads=threads_x, blocks=(blocks_x, dim_y, dim_z))
    d_out
end

@testset "CuTextureArray(::CuArray)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)
    d_a3D = CuArray(a3D)

    texarr1D = CuTextureArray(d_a1D)
    tex1D = CuTexture(texarr1D)
    @test fetch_all(tex1D) == d_a1D

    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D

    texarr3D = CuTextureArray(d_a3D)
    tex3D = CuTexture(texarr3D)
    @test fetch_all(tex3D) == d_a3D
end

@testset "CuTextureArray(::Array)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end

    texarr1D = CuTextureArray(a1D)
    copyto!(texarr1D, a1D)
    tex1D = CuTexture(texarr1D)
    @test Array(fetch_all(tex1D)) == a1D

    texarr2D = CuTextureArray(a2D)
    tex2D = CuTexture(texarr2D)
    @test Array(fetch_all(tex2D)) == a2D

    tex2D_dir = CuTexture(CuTextureArray(a2D))
    @test Array(fetch_all(tex2D_dir)) == a2D

    texarr3D = CuTextureArray(a3D)
    tex3D = CuTexture(texarr3D)
    @test Array(fetch_all(tex3D)) == a3D
end

@testset "CuTexture(::CuArray)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)

    # NOTE: tex1D is not supported for linear memory

    # This works as long as d_a2D is well pitched
    texwrap2D = CuTexture(d_a2D)
    @test fetch_all(texwrap2D) == d_a2D
end

@testset "type support" begin
    for T in (Int32, UInt32, Int16, UInt16, Int8, UInt8, Float32, Float16)
        testheight, testwidth, testdepth = 32, 32, 4
        a2D = rand(T, testheight, testwidth)
        d_a2D = CuArray(a2D)

        # Using CuTextureArray
        tex_2D = CuTexture(d_a2D)
        @test fetch_all(tex_2D) == d_a2D

        # Wrapping CuArray
        # This works as long as d_a2D is well pitched
        texwrap_2D = CuTexture(d_a2D)
        @test fetch_all(texwrap_2D) == d_a2D
    end
end

@testset "multiple channels" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int32(i), Int32(j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D

    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int16(i), Int16(j), Int16(i + j), Int16(i - j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D
end

@testset "interpolations" begin
    @testset "$interpolate $T" for T in (Float16, Float32,)
        @testset "$(N)D" for N in 1:3
            cpu_src = rand(T, fill(10, N)...)
            cpu_idx = [tuple(rand(1:0.1:10, N)...) for _ in 1:10]

            gpu_src = CuTextureArray(CuArray(cpu_src))
            gpu_idx = CuArray(cpu_idx)

            @testset "nearest neighbour" begin
                cpu_dst = similar(cpu_src, size(cpu_idx))
                cpu_int = interpolate(cpu_src, BSpline(Constant()))
                broadcast!(cpu_dst, cpu_idx, Ref(cpu_int)) do idx, int
                    int(idx...)
                end

                gpu_dst = CuArray{T}(undef, size(cpu_idx))
                gpu_tex = CuTexture(gpu_src; interpolation=CUDA.NearestNeighbour())
                broadcast!(gpu_dst, gpu_idx, Ref(gpu_tex)) do idx, tex
                    tex[idx...]
                end

                @test cpu_dst ≈ Array(gpu_dst)
            end

            @testset "linear interpolation" begin
                cpu_dst = similar(cpu_src, size(cpu_idx))
                cpu_int = interpolate(cpu_src, BSpline(Linear()))
                broadcast!(cpu_dst, cpu_idx, Ref(cpu_int)) do idx, int
                    int(idx...)
                end

                gpu_dst = CuArray{T}(undef, size(cpu_idx))
                gpu_tex = CuTexture(gpu_src; interpolation=CUDA.LinearInterpolation())
                broadcast!(gpu_dst, gpu_idx, Ref(gpu_tex)) do idx, tex
                    tex[idx...]
                end

                @test cpu_dst ≈ Array(gpu_dst) rtol=0.01
            end

            N<3 && @testset "cubic interpolation" begin
                cpu_dst = similar(cpu_src, size(cpu_idx))
                cpu_int = interpolate(cpu_src, BSpline(Cubic(Line(OnGrid()))))
                broadcast!(cpu_dst, cpu_idx, Ref(cpu_int)) do idx, int
                    int(idx...)
                end

                gpu_dst = CuArray{T}(undef, size(cpu_idx))
                gpu_tex = CuTexture(gpu_src; interpolation=CUDA.CubicInterpolation())
                broadcast!(gpu_dst, gpu_idx, Ref(gpu_tex)) do idx, tex
                    tex[idx...]
                end

                # FIXME: these results, although they look OK in an image,
                #        do not match the output from Interpolations.jl
                @test_skip cpu_dst ≈ Array(gpu_dst)
            end
        end
    end
end


text:
module ModalSSM

	include("modalparam.jl")
	include("modalparami.jl")
	include("modalparam_test.jl")
	include("mac.jl")

	export modalparam, modalparami, modalparam_test, mac

end # module

ext: jl


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

mutable struct RouteFilterPropertiesFormat <: SwaggerModel
    rules::Any # spec type: Union{ Nothing, Vector{RouteFilterRule} } # spec name: rules
    peerings::Any # spec type: Union{ Nothing, Vector{ExpressRouteCircuitPeering} } # spec name: peerings
    ipv6Peerings::Any # spec type: Union{ Nothing, Vector{ExpressRouteCircuitPeering} } # spec name: ipv6Peerings
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function RouteFilterPropertiesFormat(;rules=nothing, peerings=nothing, ipv6Peerings=nothing, provisioningState=nothing)
        o = new()
        validate_property(RouteFilterPropertiesFormat, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(RouteFilterPropertiesFormat, Symbol("peerings"), peerings)
        setfield!(o, Symbol("peerings"), peerings)
        validate_property(RouteFilterPropertiesFormat, Symbol("ipv6Peerings"), ipv6Peerings)
        setfield!(o, Symbol("ipv6Peerings"), ipv6Peerings)
        validate_property(RouteFilterPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type RouteFilterPropertiesFormat

const _property_map_RouteFilterPropertiesFormat = Dict{Symbol,Symbol}(Symbol("rules")=>Symbol("rules"), Symbol("peerings")=>Symbol("peerings"), Symbol("ipv6Peerings")=>Symbol("ipv6Peerings"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_RouteFilterPropertiesFormat = Dict{Symbol,String}(Symbol("rules")=>"Vector{RouteFilterRule}", Symbol("peerings")=>"Vector{ExpressRouteCircuitPeering}", Symbol("ipv6Peerings")=>"Vector{ExpressRouteCircuitPeering}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ RouteFilterPropertiesFormat }) = collect(keys(_property_map_RouteFilterPropertiesFormat))
Swagger.property_type(::Type{ RouteFilterPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RouteFilterPropertiesFormat[name]))}
Swagger.field_name(::Type{ RouteFilterPropertiesFormat }, property_name::Symbol) =  _property_map_RouteFilterPropertiesFormat[property_name]

function check_required(o::RouteFilterPropertiesFormat)
    true
end

function validate_property(::Type{ RouteFilterPropertiesFormat }, name::Symbol, val)
end


text:
struct Operation <: Expression
    op::Function
    args::Vector{Expression}
end

# Recursive ==
function Base.isequal(x::Operation,y::Operation)
    x.op == y.op && length(x.args) == length(y.args) && all(isequal.(x.args,y.args))
end
Base.isequal(::Operation, ::Number   ) = false
Base.isequal(::Number   , ::Operation) = false
Base.isequal(::Operation, ::Variable ) = false
Base.isequal(::Variable , ::Operation) = false
Base.isequal(::Operation, ::Constant ) = false
Base.isequal(::Constant , ::Operation) = false

Base.convert(::Type{Expr}, O::Operation) =
    build_expr(:call, Any[Symbol(O.op); convert.(Expr, O.args)])
Base.show(io::IO, O::Operation) = print(io, convert(Expr, O))

# For inv
Base.convert(::Type{Operation}, x::Number) = Operation(identity, Expression[Constant(x)])
Base.convert(::Type{Operation}, x::Operation) = x
Base.convert(::Type{Operation}, x::Expression) = Operation(identity, Expression[x])
Operation(x) = convert(Operation, x)


text:
using ProgressMeter
using Printf

function main()
    PARAMDAT = "PARAMDAT.json"

    # 膨張率 [%]
    expantion = 10

    # 構造体メモリ
    memori = Int(5.0e6)

    # ----------------------------------------------------------
    # 変数はないはず
    # ----------------------------------------------------------
    out_file,in_file,n_div,s_x,e_x,n_x,s_y,e_y,n_y,s_z,e_z,n_z = input_para(PARAMDAT)
    main_pre(n_x,n_y,s_x,e_x,s_y,e_y,s_z,e_z,n_z)

    nodes, nodes_num, element, element_shape = read_allgrid(in_file)

    println("nodes    : "*string(size(nodes)[1]))
    println("elements : "*string(size(element)[1]))

    # 空間分割
    divnum = 2^n_div

    # 辺の長さ
    lx = abs(e_x - s_x)
    ly = abs(e_y - s_y)
    lz = abs(e_z - s_z)
    lx_exp = lx *expantion/200
    ly_exp = ly *expantion/200
    lz_exp = lz *expantion/200

    # モートン単位長
    unitx = abs(e_x - s_x) / divnum *(1 + expantion/100)
    unity = abs(e_y - s_y) / divnum *(1 + expantion/100)
    unitz = abs(e_z - s_z) / divnum *(1 + expantion/100)

    # 空間分割中心算出
    dx = (e_x - s_x) / n_x
    dy = (e_y - s_y) / n_y
    dz = (e_z - s_z) / n_z

    new_morton_num = zeros(n_x,n_y,n_z,2)
    new_cell_center = zeros(n_x,n_y,n_z,3)
    for i in 1:n_x
        for j in 1:n_y
            for k in 1:n_z
                # 左上の点のモートン空間番号
                ulx = floor((s_x + dx * (i-1) -s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
                uly = floor((s_y + dy *  j    -s_y + ly_exp) / unity)
                ulz = floor((s_z + dz *  k    -s_z + lz_exp) / unitz)

                # 右下の点のモートン空間番号
                lrx = floor((s_x + dx *  i    -s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
                lry = floor((s_y + dy * (j-1) -s_y + ly_exp) / unity)
                lrz = floor((s_z + dz * (k-1) -s_z + lz_exp) / unitz)

                # 0<=s, 0<=t
                #new_morton_num[i,j,1],new_morton_num[i,j,2] = cal_morton(ulx,uly,lrx,lry)
                s,t = cal_morton(ulx,uly,ulz,lrx,lry,lrz,n_div)
                new_morton_num[i,j,k,1] = s
                new_morton_num[i,j,k,2] = t

                new_cell_center[i,j,k,1] = s_x +dx/2+dx*(i-1)
                new_cell_center[i,j,k,2] = s_y +dy/2+dy*(j-1)
                new_cell_center[i,j,k,3] = s_z +dz/2+dz*(k-1)
            end
        end
    end
    println(" fin kuukan bunkatu\n")

    # 登録管理構造体
    s = Int64(0)
    for i in 1:(n_div+1)
        s = s + 8^(i-1)
    end
    morton_tree = zeros(Int64,s,memori)
    morton_tree_num = ones(Int64,s)

    # セルを囲む直方体
    # これでモートン空間番号を算出
    cellnum = size(element)[1]
    cell_center = zeros(cellnum,3)

    for i in 1:cellnum
        loop = element_shape[i]
        x = zeros(loop)
        y = zeros(loop)
        z = zeros(loop)

        for j in 1:loop
            x[j] = nodes[element[i,j],1]
            y[j] = nodes[element[i,j],2]
            z[j] = nodes[element[i,j],3]

            cell_center[i,1] += x[j] / loop
            cell_center[i,2] += y[j] / loop
            cell_center[i,3] += z[j] / loop
        end

        maxx = maximum(x)
        minx = minimum(x)
        maxy = maximum(y)
        miny = minimum(y)
        maxz = maximum(z)
        minz = minimum(z)

        # 左上の点のモートン空間番号
        ulx = floor((minx-s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
        uly = floor((maxy-s_y + ly_exp) / unity)
        ulz = floor((maxz-s_z + lz_exp) / unitz)

        # 右下の点のモートン空間番号
        lrx = floor((maxx-s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
        lry = floor((miny-s_y + ly_exp) / unity)
        lrz = floor((minz-s_z + lz_exp) / unitz)

        # =-------------異常値をはじく事-----------
        if (ulx < 0 || lrx < 0) || (uly < 0 || lry < 0) || (ulz < 0 || lrz < 0)
            continue
        elseif (ulx > divnum || lrx > divnum) || (uly > divnum || lry > divnum) || (ulz > divnum || lrz > divnum)
            continue
        else
            #morton_num[i,1],morton_num[i,2] = cal_morton(ulx,uly,lrx,lry)
            # 0<=s, 0<=t
            s,t = cal_morton(ulx,uly,ulz,lrx,lry,lrz,n_div)

            num = Int64(t + (8^s-1)/7   +1)     # 線形空間番号

            s = morton_tree_num[num]
            morton_tree[num,s] = i
            morton_tree_num[num] += 1
        end
    end
    println(" fin kuukan bunkatu2 \n")

    # 線形探索用
    s,search0,search1,search2,search3,search4 = liner_morton(n_div)

    # 分割したセル中心と近い点を四点計算
    point = 4
    close_point_num = zeros(n_x,n_y,n_z,point)
    distance = zeros(n_x,n_y,n_z,point)
    distance[:,:,:,:] .= 100

    prog = Progress(n_x,1)
    for i in 1:n_x
        next!(prog)
        for j in 1:n_y
            for k in 1:n_z
                if new_morton_num[i,j,k,1] ==0      # 所属空間
                    temp = copy(search0)
                    loop = Int(s[1] * 8^0)
                elseif new_morton_num[i,j,k,1] ==1      # 所属空間
                    temp = copy(search1)
                    loop = Int(s[2] * 8^1)
                elseif new_morton_num[i,j,k,1] ==2      # 所属空間
                    temp = copy(search2)
                    loop = Int(s[3] * 8^2)
                elseif new_morton_num[i,j,k,1] ==3      # 所属空間
                    temp = copy(search3)
                    loop = Int(s[4] * 8^3)
                elseif new_morton_num[i,j,k,1] ==4      # 所属空間
                    temp = copy(search4)
                    loop = Int(s[5] * 8^4)
                else
                    "\n check new_morton programm \n"
                    throw(UndefVarError(:x))
                end

                # 空間事に探査
                for pointi in 1:loop
                    for l in 1:memori

                        if morton_tree[temp[pointi],l] == 0
                            break
                        end
                        bangou = morton_tree[temp[pointi],l]

                        x1 = new_cell_center[i,j,k,1]
                        y1 = new_cell_center[i,j,k,2]
                        z1 = new_cell_center[i,j,k,3]
                        x2 = cell_center[bangou,1]
                        y2 = cell_center[bangou,2]
                        z2 = cell_center[bangou,3]
                        temp_distance = ((x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 )^0.5

                        for ii in 1:point
                            if temp_distance < distance[i,j,k,ii]
                                for m in 1:(point+1-ii)-1
                                    distance[i,j,k,point+1-m] = distance[i,j,k,point-m]
                                    close_point_num[i,j,k,point+1-m] = close_point_num[i,j,k,point-m]
                                end
                                distance[i,j,k,ii] = temp_distance
                                close_point_num[i,j,k,ii] = morton_tree[temp[pointi],l]
                                break
                            end
                        end
                    end
                end
            end
        end
    end

    println(" fin kuukan bunkatu3 \n")

    #=
    println("\n ---------------------------------- \n")
    println(" start interpolation ")
    println("\n ---------------------------------- \n")

    # 補間 loop
    Qcell = zeros(n_x,n_y,5)
    for i in 1:n_x
        for j in 1:n_y
            # cell中心ベクトル p
            p  = [new_cell_center[i,j,1], new_cell_center[i,j,2]]

            # 四点ベクトル p
            x0 = cell_center[Int64(close_point_num[i,j,1]),1]
            y0 = cell_center[Int64(close_point_num[i,j,1]),2]
            x1 = cell_center[Int64(close_point_num[i,j,2]),1]
            y1 = cell_center[Int64(close_point_num[i,j,2]),2]
            x2 = cell_center[Int64(close_point_num[i,j,3]),1]
            y2 = cell_center[Int64(close_point_num[i,j,3]),2]

            p0 = [x0, y0]
            p1 = [x1, y1]
            p2 = [x2, y2]

            # 局所座標ベクトル e
            eu = [p1[1] - p0[1], p1[2] - p0[2]]
            ev = [p2[1] - p0[1], p2[2] - p0[2]]
            inv_A = inverse_matrix(eu,ev)

            # 補間距離ベクトル u
            up = zeros(2)
            for l in 1:2
                up[l] = inv_A[l,1]*(p[1]-p0[1]) + inv_A[l,2]*(p[2]-p0[2])
            end
            for k in 1:5
                Qcell[i,j,k] = readQ[Int64(close_point_num[i,j,1]),k]
                            + up[1]*readQ[Int64(close_point_num[i,j,2]),k]
                            + up[2]*readQ[Int64(close_point_num[i,j,3]),k]
                #Qcell[i,j] = Qcell[i,j]*avoga
            end
        end
    end
    =#

    println("\n ---------------------------------- \n")
    println(" start atari ")
    println("\n ---------------------------------- \n")

    #println(nodes)
    #println(close_point_num)

    # 分割した空間がセルと当たるか否かの判定
    atari = zeros(Int64,n_x,n_y,n_z)

    for i in 1:n_x
        for j in 1:n_y
            for k in 1:n_z
                for l in 1:point
                    # ---------- 過大評価 ---------------
                    cnum = Int64(close_point_num[i,j,k,l])

                    loop = element_shape[cnum]
                    x = zeros(loop)
                    y = zeros(loop)
                    z = zeros(loop)

                    for m in 1:loop
                        x[m] = nodes[element[cnum,m],1]
                        y[m] = nodes[element[cnum,m],2]
                        z[m] = nodes[element[cnum,m],3]
                    end

                    maxx = maximum(x)
                    minx = minimum(x)
                    maxy = maximum(y)
                    miny = minimum(y)
                    maxz = maximum(z)
                    minz = minimum(z)

                    x = new_cell_center[i,j,k,1]
                    y = new_cell_center[i,j,k,2]
                    z = new_cell_center[i,j,k,3]

                    if minx < x && x < maxx
                        if miny < y && y < maxy
                            if minz < z && z < maxz
                                atari[i,j,k] += 1
                            end
                        end
                    end
                end
            end
        end
    end

    # 修正当たり判定
    #=
    atari_temp = zeros(n_x,n_y,n_z)
    for i in 2:n_x-1
        for j in 2:n_y-1
            for k in 2:n_z-1
                if atari[i,j,k] == 0 && atari[i+1,j,k] ==1 && atari[i,j+1,k] ==1 && atari[i-1,j,k] ==1 && atari[i,j-1,k] ==1 && atari[i,j,k+1] ==1 && atari[i,j,k-1] ==1
                    atari_temp[i,j,k] = 1
                else atari[i,j,k] == 1 && atari[i+1,j,k] ==0 && atari[i,j+1,k] ==0 && atari[i-1,j,k] ==0 && atari[i,j-1,k] ==0 && atari[i,j,k+1] ==0 && atari[i,j,k-1] ==0
                    atari_temp[i,j,k] = 0
                end
            end
        end
    end

    for i in 2:n_x-1
        for j in 2:n_y-1
            for k in 2:n_z-1
                atari[i,j,k] = atari[i,j,k] * atari_temp[i,j,k]
            end
        end
    end
    =#

    # output
    fff = "grid_morton/atari"
    open(fff,"w") do f
        write(f,"result:atari\n")
        for i in 1:n_x
            for j in 1:n_y
                for k in 1:n_z
                    a1 = string(atari[i,j,k])
                    write(f, a1*"\n")
                end
            end
        end
    end
    println("\nwrite "*fff)

    # ダミー
    Qcell = ones(n_x,n_y,n_z,5)

    fff = "result_morton/"* out_file
    open(fff,"w") do f
        write(f,"result:rho[kg/m^3], u[m/s], v[m/s], p[Pa], T[K]\n")
        for i in 1:n_x
            for j in 1:n_y
                for k in 1:n_z

                    a1 = @sprintf("%8.8f", Qcell[i,j,k,1])
                    a2 = @sprintf("%8.8f", Qcell[i,j,k,2])
                    a3 = @sprintf("%8.8f", Qcell[i,j,k,3])
                    a4 = @sprintf("%8.8f", Qcell[i,j,k,4])
                    a5 = @sprintf("%8.8f", Qcell[i,j,k,5])

                    #=
                    a1 = string(Int64(close_point_num[i,j,1]))
                    a2 = string(Int64(close_point_num[i,j,2]))
                    a3 = string(Int64(close_point_num[i,j,3]))
                    =#

                    write(f, a1*" "*a2*" "*a3*" "*a4* " "*a5*"\n")
                end
            end
        end
    end
    println("\nwrite "*fff)

end

# ----------------
main()
# ----------------


text:
module StackOverflow

using HTTP
using CodecZlib
using JSON

mutable struct Answer
    owner::Dict
    is_accepted::Bool
    score::Int
    last_activity_date::Int
    last_edit_date::Int
    community_owned_date::Int
    locked_date::Int
    creation_date::Int
    answer_id::Int
    question_id::Int
end

mutable struct Question
    link::String
    view_count::Int
    creation_date::Int
    is_answered::Bool
    owner::Dict
    last_activity_date::Int
    score::Int
    accepted_answer_id::Int
    question_id::Int
    tags::Vector
    title::String
    answer_count::Int
end

include("questions.jl")
export
    getrecentquestionsfortag
    getquestionfromanswer

include("privileges.jl")
export
    checkmyprivileges,
    getprivileges

include("utils.jl")
export
    convert_HTTP_Response_To_JSON

include("answers.jl")
export
    getanswers
    getanswerfromquestion

include("badges.jl")
export
    getrecentbadgerecipients

include("siteinfo.jl")
export
    getsiteinfo

include("search.jl")
export
    searcherror

include("tags.jl")
export
    searchtag

include("users.jl")
export
    getusers

include("comments.jl")
export
    getcomments

include("suggestededits.jl")
export
    getsuggested_edits

include("errors.jl")
export
    geterrors

include("unanswered.jl")
export
    getopenquestions

include("noanswers.jl")
export
    getnoanswerquestions

include("featuredquestions.jl")
export
    getfeaturedquestions

include("moderators.jl")
export
    getmoderators

end # module


text:
using MeshPlotter
using Test

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


text:
using StanModels

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "WaffleDivorce.csv"), delim=';')
mean_ma = mean(df[!, :MedianAgeMarriage])
df[!, :MedianAgeMarriage_s] =
  convert(Vector{Float64},  (df[!, :MedianAgeMarriage]) .-
    mean_ma)/std(df[!, :MedianAgeMarriage]);

# Define the Stan language model

m5_1s = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] divorce; // Predictor
 vector[N] median_age; // Outcome
}

parameters {
 real a; // Intercept
 real bA; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
  # priors
  a ~ normal(10, 10);
  bA ~ normal(0, 1);
  sigma ~ uniform(0, 10);

  # model
  divorce ~ normal(a + bA*median_age , sigma);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m5.1s", m5_1s);

# Input data for cmdstan

m5_1_data = Dict("N" => length(df[!, :Divorce]), "divorce" => df[!, :Divorce],
    "median_age" => df[!, :MedianAgeMarriage_s]);

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m5_1_data);

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.69 0.22  9.34 10.03  2023    1
bA    -1.04 0.21 -1.37 -0.71  1882    1
sigma  1.51 0.16  1.29  1.79  1695    1
"

# Describe the draws
if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
# structs to represent state transition/occupation priors.

abstract type StatePrior end

Base.broadcastable(s::StatePrior) = Ref(s)

"""
    add(::StatePrior, state, n)

Update the sufficient statistics of the state prior, returning a new version of
the state prior and the index corresponding to the state (possibly changed in
the update).
"""
function add(::StatePrior, state, n) end

"""
    state_to_index(::StatePrior, state)

Convert a state to the index of the corresponding mixture component.  Defaults
to the identity mapping.
"""
state_to_index(::StatePrior, state) = state

function log_prior(p::StatePrior)
    cands = candidates(p)
    log_weights = log_prior.(p, cands)
    log_weights .-= StatsFuns.logsumexp(log_weights)
    return log_weights
end

function Base.rand(p::StatePrior)
    weights = Weights(exp.(log_prior(p)))
    sample(candidates(p), weights)
end

StatsBase.entropy(p::StatePrior) = entropy(exp.(log_prior(p)))

function simulate(p::StatePrior, n::Int)
    states = Int[]
    for _ in 1:n
        p, x = add(p, rand(p))
        push!(states, x)
    end
    p, states
end

struct ChineseRestaurantProcess <: StatePrior
    α::Float64
    N::Vector{Float64}
end

ChineseRestaurantProcess(α::Float64) = ChineseRestaurantProcess(α, Vector{Float64}())

candidates(crp::ChineseRestaurantProcess) = 1:length(crp.N)+1
function add(crp::ChineseRestaurantProcess, x::Int, n::Float64=1.)
    N = copy(crp.N)
    if x == length(N)+1
        push!(N, n)
    else
        N[x] += n
    end
    return ChineseRestaurantProcess(crp.α, N), x
end

log_prior(crp::ChineseRestaurantProcess, x::Int) =
    x == length(crp.N)+1 ? log(crp.α) : log(crp.N[x])

marginal_log_prior(crp::ChineseRestaurantProcess) =
    log_prior = sum(lgamma(n) for n in crp.N) + length(crp.N) * log(crp.α)

"""
    StickyCRP <: StatePrior

Represents a "sticky" Chinese Restaurant Process, where there's a constant,
additional probability that the state at time t is the same as time t-1.  That
is, with probability ρ, ``x_t = x_{t-1}``, and with probability ``(1-\rho)``, a
new ``x`` is sampled according to a CRP.

There's some additional book-keeping that's necessary, as "sticking" transitions
are not informative for the CRP prior, but the proportion of transitions that
are same-state are informative when re-sampling the stickiness parameter.  So we
need to track
1. The number of times a state ``j`` was visited when ``s=0`` (non-sticky).
2. The number of times a state ``j`` was visited when the previous state was
   also ``j``, _regardless_ of whether sticky or not.

The first is stored on `N`, and the second on `Nsame` (and the value of the last
`x` on `last`).

When generating candidates the value `x=0` is used as a sentinal for sticking.
"""
struct StickyCRP <: StatePrior
    α::Float64
    ρ::Float64                  # probability of "sticking"
    last::Int64
    N::Vector{Float64}          # number of non-sticky occuptations of state
    Nsame::Vector{Float64}      # number of transitions to same state (sticky or not)
end

StickyCRP(α::Float64, κ::Float64) = StickyCRP(α, κ, 1, Vector{Float64}(), Vector{Float64}())
Base.isempty(crp::StickyCRP) = isempty(crp.N)

state_to_index(crp::StickyCRP, state::Int) = state == 0 ? crp.last : state

candidates(crp::StickyCRP) = isempty(crp) ? (1:1) : 0:length(crp.N)+1
add(crp::StickyCRP, x::Int, n::Float64=1.0) = add(crp::StickyCRP, state_to_index(crp, x), x==0, n)
function add(crp::StickyCRP, x::Int, sticky::Bool, n::Float64)
    N = copy(crp.N)
    Nsame = copy(crp.Nsame)

    # new component
    if x == length(N)+1
        N = push!(N, 0.)
        Nsame = push!(Nsame, 0.)
    end

    # same transition -> add to Nsame
    if x == crp.last
        Nsame[x] += n
    end

    # sticky -> don't update N
    if !sticky
        N[x] += n
    end

    return StickyCRP(crp.α, crp.ρ, x, N, Nsame), x
end

function log_prior(crp::StickyCRP, x::Int)
    if x == 0
        # sticky: ∝ ρ / (1-ρ) * (∑N+α)
        return logit(crp.ρ) + log(sum(crp.N)+crp.α)
    elseif 0 < x ≤ length(crp.N)
        return log(crp.N[x])
    elseif x == length(crp.N)+1
        return log(crp.α)
    else
        throw(ArgumentError("state $x is invalid for sticky CRP with " *
                            "$(length(crp.N)) components (valid values are " *
                            "0..$(length(crp.N)+1)"))
    end
end

################################################################################
# a changepoint state prior

"""
    struct ChangePoint <: StatePrior

A changepoint prior.

# Fields

* logp::Float64 - log-probability of a change
* k::Int - number of states.
* n::Int - number of data points seen so far.
"""
struct ChangePoint <: Particles.StatePrior
    logp::Float64
    k::Int
    n::Int
end

function ChangePoint(p::Float64)
    0 ≤ p ≤ 1 || throw(ArgumentError("ChangePoint probability p must be between 0 and 1 (got p=$p)"))
    ChangePoint(log(p), 0, 0)
end

function add(cp::ChangePoint, state, n=1.)
    state == cp.k || state == cp.k+1 ||
        throw(ArgumentError("ChangePoint state must be same as last ($(cp.k)) or last+1, got $state"))
    ChangePoint(cp.logp, cp.k + (state>cp.k), cp.n+1), state
end

candidates(cp::ChangePoint) =
    max(cp.k,1):cp.k+1

# prior is 1 for starting transition (k=0).  then p for a change, 1-p for stay
log_prior(cp::ChangePoint, state::Int) =
    cp.k == 0 ? 0. :
    state == cp.k ? StatsFuns.log1mexp(cp.logp) : cp.logp

# marginal prior is 1 if there's just one data point.  if > 1, there are k-1
# changepoints, and n-1 opportunities for a changepoint.
marginal_log_prior(cp::ChangePoint) =
    cp.n == 1 ? 0. : cp.logp*(cp.k-1) + StatsFuns.log1mexp(cp.logp)*(cp.n-cp.k)

"""
    struct NStatePrior <: StatePrior

A state prior with a fixed, finite number of states.

# Fields
* `counts::Vector{Float64}`
* `total_count::Float64`
* `n::Int` number of states
"""
struct NStatePrior <: StatePrior
    counts::Vector{Float64}
    total_count::Float64
    n::Int
end

NStatePrior(counts::Vector{<:Real}) = NStatePrior(counts, sum(counts), length(counts))

function add(ns::NStatePrior, state, n=1.)
    state ∈ 1:ns.n || throw(ArgumentError("Invalid state for NStatePrior: $state " *
                                          "(valid states are 1:$(ns.n))"))
    counts = copy(ns.counts)
    counts[state] += n
    NStatePrior(counts, ns.total_count + n, ns.n), state
end

candidates(ns::NStatePrior) = 1:ns.n
log_prior(ns::NStatePrior, state) = log(ns.counts[state]) - log(ns.total_count)
marginal_log_prior(ns::NStatePrior) = 0.


text:
# theobins
function theobins_test()
	d = Distributions.Normal()
	b = collect(-3.:1:3);
	n = 100;
	bc = HypoTest.theobins(d,b,n);
	@test trunc.(Int,round.(bc)) == [2,14,34,34,14,2];
	q = [Distributions.quantile(d,0.99),Distributions.quantile(d,0.99999)]
	@test HypoTest.theobins(d,q,n) ≈ [0.999] atol = 0.001
end

# countbins
function countbins_test()
	x = collect(1.:1.:10);
	bc = HypoTest.countbins(x,[0.,1.,9.,10.]);
	@test bc == [0,8,2]; # see PyPlot.plt[:hist](x,[0.,1.,9.,10.])
	@test HypoTest.countbins(x,[1.,1.,9.]) == [0,9];
	@test HypoTest.countbins([-1.001,0.,0.999],[-2.,-1.,0.,1]) == [1,0,2]
end

# histbins
function histbins_test()
	hb = HypoTest.histbins([0.,1.],method=10)
	@test hb == collect(0:0.1:1.)
	hb = HypoTest.histbins([-1.,0.,1.,2.],method="square")
	@test hb == collect(-1:1.5:2)
end

# prepdata
function prepdata_test()
	@test [1,2,3.] == HypoTest.prepdata([1.,2.,3.])
	@test [-1.,0.99] == HypoTest.prepdata([NaN,-1.,NaN,0.99,NaN])
	# two vectors
	x,y = HypoTest.prepdata([1.,2,NaN],[2.,NaN,9]);
	@test x == [1.]
	@test y == [2.];
	x,y = HypoTest.prepdata([1.,2.,3.],[2.,3.,4.]);
	@test x == [1.,2.,3.];
	@test y == [2.,3.,4.];
end

# Run functions
theobins_test();
countbins_test();
histbins_test();
prepdata_test();

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
## set up McMC
nsamples, nchains, nchainsatone = 20001, 8, 1
Tmax = 2.50
addprocs(nchains)
@info "workers are $(workers())"
@everywhere using Distributed
@everywhere using transD_GP
## run McMC
@time transD_GP.main(opt, csem, Tmax=Tmax, nsamples=nsamples, nchains=nchains, nchainsatone=nchainsatone)
rmprocs(workers())


text:
function u_pcm(x, c, η; m)
    d_squared = FuzzySets.euclidean_distance(x, c; squared=true)
    u = 1 / (1 + (d_squared / η)^(1 / (m - 1)))
    u
end


text:
# Taken from julia/base/util.jl
const _mem_units = ["byte", "KiB", "MiB", "GiB", "TiB", "PiB"]
const _cnt_units = ["", " k", " M", " G", " T", " P"]
function prettyprint_getunits(value, numunits, factor)
    if value == 0 || value == 1
        return (value, 1)
    end
    unit = ceil(Int, log(value) / log(factor))
    unit = min(numunits, unit)
    number = value/factor^(unit-1)
    return number, unit
end

function format_bytes(bytes)
    bytes, mb = prettyprint_getunits(bytes, length(_mem_units), Int64(1024))
    if mb == 1
        @sprintf("%d %s%s", bytes, _mem_units[mb], bytes==1 ? "" : "s")
    else
        @sprintf("%.3f %s", bytes, _mem_units[mb])
    end
end

format_number(n::Int) = Formatting.format(n, commas=true)

format_observations(x::Array{T, 2}) where T <: Real = format_number(size(x, 2))

function log_experiment_info(experiment)
    split = experiment[:split_strategy]
    debug(LOGGER, "Experiment info: \n\tdata file: $(experiment[:data_file]),\n" *
        "\tparams: adjust_K = $(experiment[:param][:adjust_K]), num_al_iterations = $(experiment[:param][:num_al_iterations])\n" *
        "\tquery_strategy: $(experiment[:query_strategy][:type])\n" *
        "\tmodel: $(experiment[:model][:type])\n" *
        "\tinit_strategy: $(experiment[:model][:init_strategy])\n" *
        "\tinitial_pools: $(countmap(experiment[:param][:initial_pools]))\n" *
        "\tsplit_strategy: \n\t\ttrain: $(split.train_strat); \n\t\ttest:  $(split.test_strat); \n\t\ttrain: $(split.query_strat)")
end

function Base.push!(logger::Logger, handler::Handler, name::String)
    logger.handlers[name] = handler
end


text:
using Test

@testset "Atmos Model" begin
    include("ref_state.jl")
end


text:
using VegaLite
using RDatasets

mpg = dataset("ggplot2", "mpg")

data_values(mpg) +
  mark_point() +
  encoding_x_quant(:Cty) +
  encoding_y_quant(:Hwy)

data_values(mpg) +
  mark_point() +
  encoding_x_quant(:Cty, axis=false) +
  encoding_y_quant(:Hwy, scale=scale(zero=false)) +
  encoding_color_nominal(:Manufacturer) +
  config_cell(width=350, height=400)

data_values(mpg) +
  mark_line() +
  encoding_x_ord(:Year,
                 axis  = axis(labelAngle=-45, labelAlign="right"),
                 scale = scale(bandSize=50)) +
  encoding_y_quant(:Hwy, aggregate="mean") +
  encoding_color_nominal(:Manufacturer)

data_values(mpg) +
  mark_point() +
  encoding_column_ord(:Cyl) +
  encoding_row_ord(:Year) +
  encoding_x_quant(:Displ) +
  encoding_y_quant(:Hwy) +
  encoding_size_quant(:Cty) +
  encoding_color_nominal(:Manufacturer)

data_values(mpg) +
  mark_text() +
  encoding_column_ord(:Cyl) +
  encoding_row_ord(:Year) +
  encoding_text_quant(:Displ, aggregate="mean") +
  config_mark(fontStyle="italic", fontSize=12, font="courier")


text:
# A line segment defined by its two endpoints.
struct LineSegment{N,T} <: Edge{N,T}
    points::SVector{2, Point{N,T}}
end

const LineSegment_2D = LineSegment{2}
const LineSegment_3D = LineSegment{3}

# Constructors & Conversions
# -------------------------------------------------------------------------------------------------
LineSegment(p₁::Point{N,T}, p₂::Point{N,T}) where {N,T} = LineSegment{N,T}(SVector(p₁, p₂))
LineSegment{N}(p₁::Point{N,T}, p₂::Point{N,T}) where {N,T} = LineSegment{N,T}(SVector(p₁, p₂))
LineSegment{N,T}(p₁::Point{N,T}, p₂::Point{N,T}) where {N,T} = LineSegment{N,T}(SVector(p₁, p₂))
function LineSegment{N,T₁}(p₁::Point{N,T₂}, p₂::Point{N,T₂}) where {N,T₁,T₂}
    return LineSegment{N,T₁}(SVector(Point{N, T₁}(p₁), Point{N, T₁}(p₂)))
end

# Methods
# -------------------------------------------------------------------------------------------------
# Interpolation
# l(0) yields points[1], and l(1) yields points[2]
@inline (l::LineSegment{N,T})(r) where {N,T} = l[1] + (l[2] - l[1])T(r)
@inline arclength(l::LineSegment_2D) = distance(l[1], l[2])

function Base.intersect(𝐥₁::LineSegment_3D{T}, 𝐥₂::LineSegment_3D{T}) where {T}
    # NOTE: Doesn't work for colinear/parallel lines. (𝐯 × 𝐮 = 0⃗).
    # Using the equation of a line in parametric form
    # For 𝐥₁ = 𝐱₁ + r𝐯 and 𝐥₂ = 𝐱₂ + s𝐮
    # 1) 𝐱₁ + r𝐯 = 𝐱₂ + s𝐮                  subtracting 𝐱₁ from both sides
    # 2) r𝐯 = (𝐱₂-𝐱₁) + s𝐮                  𝐰 = 𝐱₂-𝐱₁
    # 3) r𝐯 = 𝐰 + s𝐮                        cross product with 𝐮 (distributive)
    # 4) r(𝐯 × 𝐮) = 𝐰 × 𝐮 + s(𝐮 × 𝐮)        𝐮 × 𝐮 = 0
    # 5) r(𝐯 × 𝐮) = 𝐰 × 𝐮                   let 𝐰 × 𝐮 = 𝐚 and 𝐯 × 𝐮 = 𝐜
    # 6) r𝐜 = 𝐚                             dot product 𝐜 to each side
    # 7) r𝐜 ⋅ 𝐜 = 𝐚 ⋅ 𝐜                     divide by 𝐜 ⋅ 𝐜
    # 8) r = 𝐚 ⋅ 𝐜/𝐜 ⋅ 𝐜                    definition of 2-norm
    # 9) r = 𝐚 ⋅ 𝐜/‖𝐜‖
    # Note that if the lines are parallel or collinear, 𝐜 = 𝐯 × 𝐮 = 0⃗
    # We need to ensure r, s ∈ [0, 1].
    # 𝐱₂ + s𝐮 = 𝐱₁ + r𝐯                     subtracting 𝐱₂ from both sides
    # s𝐮 = -𝐰 + r𝐯                          cross product with 𝐰
    # s(𝐮 × 𝐰) = -𝐰 × 𝐰 + r(𝐯 × 𝐰)          𝐰 × 𝐰 = 0 and substituting for r
    # s(𝐮 × 𝐰) = (𝐯 × 𝐰)[𝐚 ⋅ 𝐜/‖𝐜‖]         using 𝐮 × 𝐰 = -(𝐰 × 𝐮), likewise for 𝐯 × 𝐰
    # s(𝐰 × 𝐮) = (𝐰 × 𝐯)[𝐚 ⋅ 𝐜/‖𝐜‖]         let 𝐰 × 𝐯 = 𝐛. use 𝐰 × 𝐮 = 𝐚
    # s𝐚 = 𝐛[𝐚 ⋅ 𝐜/‖𝐜]                      dot product 𝐚 to each side
    # s(𝐚 ⋅ 𝐚) = (𝐛 ⋅ 𝐚)[𝐚 ⋅ 𝐜/‖𝐜‖]         definition of 2-norm and divide
    # s = (𝐚 ⋅ 𝐛)(𝐚 ⋅ 𝐜)/(‖𝐚‖‖𝐜‖)           substitute for r
    # s = r𝐚 ⋅ 𝐛/‖𝐚‖
    ϵ = T(5e-6)
    𝐰 = 𝐥₂[1] - 𝐥₁[1]
    𝐯 = 𝐥₁[2] - 𝐥₁[1]
    𝐮 = 𝐥₂[2] - 𝐥₂[1]
    𝐜 = 𝐯 × 𝐮
    # Note: 0 ≤ 𝐜 ⋅ 𝐰, and the minimum distance between two lines is d = (𝐜 ⋅ 𝐰 )/‖𝐜‖.
    # Hence 𝐜 ⋅𝐰 ≈ 0 for the lines to intersect
    # (https://math.stackexchange.com/questions/2213165/find-shortest-distance-between-lines-in-3d)
    𝐜 ⋅𝐰  ≤ T(1e-8) || return (false, Point_3D{T}(0,0,0))
    𝐚 = 𝐰 × 𝐮
    𝐛 = 𝐰 × 𝐯
    r = (𝐚 ⋅ 𝐜)/(𝐜 ⋅ 𝐜)
    s = r*(𝐚 ⋅ 𝐛)/(𝐚 ⋅ 𝐚)
    return (T(1e-8)^2 < abs(𝐜 ⋅ 𝐜) && -ϵ ≤ r && r ≤ 1 + ϵ && -ϵ ≤ s && s ≤ 1 + ϵ  , 𝐥₂(s)) # (hit, point)
end

function Base.intersect(𝐥₁::LineSegment_2D{T}, 𝐥₂::LineSegment_2D{T}) where {T}
    # NOTE: Doesn't work for colinear/parallel lines. (𝐯 × 𝐮 = 0⃗). Also, the cross product
    # operator for 2D points returns a scalar (the 2-norm of the cross product).
    #
    # From the 3D intersection routine we know:
    # r = 𝐚 ⋅ 𝐜/𝐜 ⋅ 𝐜
    # s = (𝐚 ⋅ 𝐛)(𝐚 ⋅ 𝐜)/(‖𝐚‖‖𝐜‖)
    # Since the 2D cross product returns a scalar
    # r = 𝐚 ⋅ 𝐜/𝐜 ⋅ 𝐜 = 𝐚/𝐜 = a/c
    # s = (𝐚 ⋅ 𝐛)(𝐚 ⋅ 𝐜)/(‖𝐚‖‖𝐜‖) = 𝐛/𝐜 = b/c
    #
    # Simply evaluating everything removes branches and is faster than failing early with
    # 1e-8 < abs(c) or delaying division by vxu and testing against r and s's numerators.
    # This has been tested.
    ϵ = T(5e-6)
    𝐰 = 𝐥₂[1] - 𝐥₁[1]
    𝐯 = 𝐥₁[2] - 𝐥₁[1]
    𝐮 = 𝐥₂[2] - 𝐥₂[1]
    c = 𝐯 × 𝐮
    r = 𝐰 × 𝐮/c
    s = 𝐰 × 𝐯/c
    # -ϵ ≤ r ≤ 1 + ϵ introduces a branch, but -ϵ ≤ r && r ≤ 1 + ϵ doesn't for some reason.
    return (T(1e-8) < abs(c) && -ϵ ≤ r && r ≤ 1 + ϵ && -ϵ ≤ s && s ≤ 1 + ϵ  , 𝐥₂(s)) # (hit, point)
end

# Return if the point is left of the line segment
#   p    ^
#   ^   /
# v⃗ |  / u⃗
#   | /
#   o
#   We rely on v⃗ × u⃗ = |v⃗||u⃗|sin(θ). We may determine if θ ∈ (0, π] based on the sign of v⃗ × u⃗
@inline function isleft(p::Point_2D, l::LineSegment_2D)
    u⃗ = l[2] - l[1]
    v⃗ = p - l[1]
    return u⃗ × v⃗ >= 0
end

# A random line within [0, 1] × [0, 1]
function Base.rand(::Type{LineSegment{N,F}}) where {N,F}
    return LineSegment{N,F}(rand(Point{N,F}, 2))
end

# N random lines within [0, 1] × [0, 1]
function Base.rand(::Type{LineSegment{N,F}}, NP::Int64) where {N,F}
    return [ rand(LineSegment{N,F}) for i ∈ 1:NP ]
end

# # Plot
# # -------------------------------------------------------------------------------------------------
# if enable_visualization
#     function convert_arguments(LS::Type{<:LineSegments}, l::LineSegment_2D)
#         return convert_arguments(LS, [l[1], l[2]])
#     end
#
#     function convert_arguments(LS::Type{<:LineSegments}, L::Vector{<:LineSegment_2D})
#         return convert_arguments(LS, reduce(vcat, [[l[1], l[2]] for l in L]))
#     end
# end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
### Sparse-array iterators

const SubSparseArray{I,T,N,P<:AbstractSparseArray} = SubArray{T,N,P,I,false}

## SparseMatrixCSC

const SubSparseMatrixCSC{I,T,N,P<:SparseMatrixCSC} = SubArray{T,N,P,I,false}
const ContiguousCSC{I<:Union{Tuple{CSlice,Union{Int,UnitRange{Int}}},Tuple{UnitRange{Int},Int}},T,N,P<:SparseMatrixCSC} = Union{<:SparseMatrixCSC{T},SubSparseMatrixCSC{I,T,N,P}}

indextype(::Type{SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = Ti
indextype(A::SparseMatrixCSC) = indextype(typeof(A))

# Indexing along a particular column
struct IndexCSC
    row::Int; col::Int  # where you are currently (might not be a stored value)
    stored::Bool        # true if this location corresponds to a stored value
    cscindex::Int       # for stored value, the index into rowval & nzval
end

function getindex(I::IndexCSC, d)
    @boundscheck d==1 || d==2 || Base.throw_boundserror(I, d)
    ifelse(d == 1, I.row, I.col)
end

@inline getindex(A::SparseMatrixCSC, i::IndexCSC) = (@inbounds ret = i.stored ? A.nzval[i.cscindex] : zero(eltype(A)); ret)
@inline getindex(A::SubSparseMatrixCSC, i::IndexCSC) = A.parent[i]
# @inline function getindex(a::AbstractVector, i::IndexCSC)
#     @boundscheck 1 <= i.rowval <= length(a)
#     @inbounds ret = a[i.rowval]
#     ret
# end

@inline setindex!(A::SparseMatrixCSC, val, i::IndexCSC) = (@inbounds A.nzval[i.cscindex] = val; val)
@inline setindex!(A::SubSparseMatrixCSC, val, i::IndexCSC) = A.parent[i] = val
# @inline function setindex!(a::AbstractVector, val, i::IndexCSC)
#     @boundscheck 1 <= i.rowval <= length(a) || throw(BoundsError(a, i.rowval))
#     @inbounds a[i.rowval] = val
#     val
# end

struct IteratorCSC{isstored,S<:ContiguousCSC}
    A::S
    cscrange::UnitRange{Int}

    IteratorCSC{isstored,S}(A::S, ::CSlice, ::CSlice) where {isstored,S<:SparseMatrixCSC} = new{isstored,S}(A, 1:length(A.nzval))
end
function IteratorCSC{isstored,S}(A::S, c::CSlice, j) where {isstored,S<:ContiguousCSC}
    IteratorCSC{isstored,S}(A, c, j)
end
function IteratorCSC{E}(A::S, c::CSlice, j) where {E,S<:ContiguousCSC}
    IteratorCSC{E,S}(A, c, j)
end

struct ColIteratorCSC{isstored,S<:ContiguousCSC}
    A::S
    col::Int
    cscrange::UnitRange{Int}

    function ColIteratorCSC{E}(A::S, ::CSlice, col::Integer) where {E,S<:SparseMatrixCSC}
        @boundscheck 1 <= col <= size(A, 2) || throw(BoundsError(A, (:,col)))
        @inbounds r = A.colptr[col]:A.colptr[col+1]-1
        new{E,S}(A, col, r)
    end
    function ColIteratorCSC{E}(A::S, ::CSlice, col::Integer) where {E,I<:Tuple{CSlice,Any},S<:SubSparseMatrixCSC{I}}
        @boundscheck 1 <= col <= size(A, 2) || throw(BoundsError(A, (:,col)))
        @inbounds j = A.indexes[2][col]
        @inbounds r = A.parent.colptr[j]:A.parent.colptr[j+1]-1
        new{E,I}(A, col, r)
    end
    function ColIteratorCSC{E}(A::S, ::CSlice, col::Integer) where {E,I<:Tuple{UnitRange{Int},Any},S<:SubSparseMatrixCSC{I}}
        @boundscheck 1 <= col <= size(A, 2) || throw(BoundsError(A, (:,col)))
        @inbounds j = A.indexes[2][col]
        @inbounds r1, r2 = Int(A.parent.colptr[j]), Int(A.parent.colptr[j+1]-1)
        rowval = A.parent.rowval
        i = A.indexes[1]
        r1 = searchsortedfirst(rowval, first(i), r1, r2, Forward)
        r1 <= r2 && (r2 = searchsortedlast(rowval, last(i), r1, r2, Forward))
        new{E,I}(A, col, r1:r2)
    end
    function ColIteratorCSC{E}(A::S, i::UnitRange, col::Integer) where {E,S<:SparseMatrixCSC}
        @boundscheck 1 <= col <= size(A, 2) || throw(BoundsError(A, (i,col)))
        @boundscheck (1 <= first(i) && last(i) <= size(A, 1)) || throw(BoundsError(A, (i,col)))
        @inbounds r1, r2 = Int(A.parent.colptr[j]), Int(A.parent.colptr[j+1]-1)
        rowval = A.parent.rowval
        r1 = searchsortedfirst(rowval, first(i), r1, r2, Forward)
        r1 <= r2 && (r2 = searchsortedlast(rowval, last(i), r1, r2, Forward))
        new{E,S}(A, col, r1:r2)
    end
end
# Default is to visit each site, not just the stored sites
ColIteratorCSC(A::ContiguousCSC, i, col::Integer) = ColIteratorCSC{false}(A, i, col)
# ...but you can choose with ColIteratorCSC{true/false}(A, col)
#(::Type{ColIteratorCSC{E}})(A::ContiguousCSC, i, col::Integer) where E = ColIteratorCSC{E,typeof(A)}(A, i, col)

# Iteration when you're visiting every entry
# The column iterator state has the following structure:
#    (row::Int, nextrowval::Ti<:Integer, cscindex::Int)
# nextrowval = A.rowval[cscindex], but we cache it in the state to
# avoid looking it up each time. We use it to decide when the cscindex
# needs to be incremented.
# The full-matrix iterator is similar, except it adds the column:
#    (row::Int, col::Int, nextrowval::Ti<:Integer, nextcolval::Ti, cscindex::Int)
length(iter::IteratorCSC{false})  = length(iter.A)
length(iter::ColIteratorCSC{false})  = size(iter.A, 1)
function start(iter::IteratorCSC{false})
    cscindex = first(iter.cscrange)
    nextrow, nextcol = _nextrowcolval(iter, 0, cscindex)
    (1, 1, nextrow, nextcol, cscindex)
end
function start(iter::ColIteratorCSC{false})
    cscindex = first(iter.cscrange)
    nextrowval = _nextrowval(iter, cscindex)
    (1, nextrowval, cscindex)
end
done(iter::IteratorCSC{false}, s) = s[2] > size(iter.A, 2)
done(iter::ColIteratorCSC{false}, s) = s[1] > size(iter.A, 1)
function iterate(iter::IteratorCSC{false,S}, s=start(iter)) where S<:SparseMatrixCSC
    done(iter, s) && return nothing
    row, col, nextrowval, nextcolval, cscindex = s
    item = IndexCSC(row, col, row==nextrowval && col==nextcolval, cscindex)
    newrow = row+1
    newcol = col
    if newrow > size(iter.A, 1)
        newrow = 1
        newcol += 1
    end
    if item.stored
        nrv, ncv = _nextrowcolval(iter, col, cscindex+1)
        return (item, (newrow, newcol, nrv, ncv, cscindex+1))
    end
    return (item, (newrow, newcol, nextrowval, nextcolval, cscindex))
end
function iterate(iter::ColIteratorCSC{false,S}, s=start(iter)) where S<:SparseMatrixCSC
    done(iter, s) && return nothing
    row, nextrowval, cscindex = s
    item = IndexCSC(row, iter.col, row==nextrowval, cscindex)
    item.stored ? (item, (row+1, _nextrowval(iter, cscindex+1), cscindex+1)) :
                  (item, (row+1, nextrowval, cscindex))
end
function _nextrowval(iter::ColIteratorCSC, cscindex)
    if cscindex <= last(iter.cscrange)
        return iter.A.rowval[cscindex]
    end
    convert(indextype(iter.A), size(iter.A, 1)+1)  # out-of-bounds fallback
end
function _nextrowcolval(iter::IteratorCSC, col, cscindex)
    if cscindex <= last(iter.cscrange)
        nextcol = col
        nextcscindex = iter.A.colptr[col+1]
        if cscindex >= nextcscindex
            nextcol = findnext(j->j!=nextcscindex, iter.A.colptr, col+1)-1
        end
        return (iter.A.rowval[cscindex], nextcol)
    end
    # out-of-bounds fallback
    convert(indextype(iter.A), size(iter.A, 1)+1), size(iter.A, 2)+1
end

# Iteration when you're visting just the stored entries
# We use similar caching tricks with nextcol and nextcolptrindex for IteratorCSC
length(iter::IteratorCSC{true}) = length(iter.cscrange)
length(iter::ColIteratorCSC{true}) = length(iter.cscrange)
function start(iter::IteratorCSC{true})
    nextcol = findfirst(j->j!=1, iter.A.colptr)-1
    nextcolptrindex = iter.A.colptr[nextcol+1]
    (nextcol, nextcolptrindex, first(iter.cscrange))
end
start(iter::ColIteratorCSC{true}) = first(iter.cscrange)
done(iter::IteratorCSC{true}, s) = last(iter.cscrange) < s[3]
done(iter::ColIteratorCSC{true}, s) = last(iter.cscrange) < s
function iterate(iter::IteratorCSC{true,S}, s=start(iter)) where S<:SparseMatrixCSC
    done(iter, s) && return nothing
    @inbounds begin
        col, nextcolptrindex, cscindex = s
        row = iter.A.rowval[cscindex]
    end
    nextcol = col
    if s == nextcolptrindex
        tmp = nextcolptrindex  # work around julia #15276
        nextcol = findnext(j->j!=tmp, iter.A.colptr, col)-1
        nextcolptrindex = iter.A.colptr[nextcol+1]
    end
    idx = IndexCSC(row, col, true, cscindex)
    (idx, (nextcol, nextcolptrindex, cscindex+1))
end
function iterate(iter::ColIteratorCSC{true,S}, s=start(iter)) where S<:SparseMatrixCSC
    done(iter, s) && return nothing
    @inbounds row = iter.A.rowval[s]
    idx = IndexCSC(row, iter.col, true, s)
    (idx, s+1)
end
function iterate(iter::ColIteratorCSC{true,S}, s=start(iter)) where S<:SubSparseMatrixCSC
    done(iter, s) && return nothing
    @inbounds row = iter.A.parent.rowval[s]
    idx = IndexCSC(row, iter.col, true, s)
    (idx, s+1)
end

# nextstored{S<:SparseMatrixCSC}(iter::ColIteratorCSC{S}, s, index::Integer) =

each(w::ArrayIndexingWrapper{A,<:NTuple{N,<:CSlice},true,isstored}) where {A<:SparseMatrixCSC,N,isstored} = IteratorCSC{isstored}(w.data, w.indexes...)
each(w::ArrayIndexingWrapper{A,I,true,isstored}) where {A<:SparseMatrixCSC,I,isstored} = ColIteratorCSC{isstored}(w.data, w.indexes...)
each(w::ArrayIndexingWrapper{A,<:NTuple{N,<:CSlice},false,isstored}) where {A<:SparseMatrixCSC,N,isstored} = ValueIterator(w.data, IteratorCSC{isstored}(w.data, w.indexes...))
each(w::ArrayIndexingWrapper{A,I,false,isstored}) where {A<:SparseMatrixCSC,I,isstored} = ValueIterator(w.data, ColIteratorCSC{isstored}(w.data, w.indexes...))


text:
using CLIMA.PlanetParameters
export PeriodicBC, NoFluxBC, InitStateBC, DYCOMS_BC, RayleighBenardBC

#TODO: figure out a better interface for this.
# at the moment we can just pass a function, but we should do something better
# need to figure out how subcomponents will interact.
function atmos_boundary_state!(::Rusanov, f::Function, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t, _...)
  f(stateP, auxP, nM, stateM, auxM, bctype, t)
end

function atmos_boundary_state!(::CentralNumericalFluxDiffusive, f::Function,
                               m::AtmosModel, stateP::Vars, diffP::Vars,
                               auxP::Vars, nM, stateM::Vars, diffM::Vars,
                               auxM::Vars, bctype, t, _...)
  f(stateP, diffP, auxP, nM, stateM, diffM, auxM, bctype, t)
end

# lookup boundary condition by face
function atmos_boundary_state!(nf::Rusanov, bctup::Tuple, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t, _...)
  atmos_boundary_state!(nf, bctup[bctype], m, stateP, auxP, nM, stateM, auxM,
                        bctype, t)
end

function atmos_boundary_state!(nf::CentralNumericalFluxDiffusive,
                               bctup::Tuple, m::AtmosModel, stateP::Vars,
                               diffP::Vars, auxP::Vars, nM, stateM::Vars,
                               diffM::Vars, auxM::Vars, bctype, t, _...)
  atmos_boundary_state!(nf, bctup[bctype], m, stateP, diffP, auxP, nM, stateM,
                        diffM, auxM, bctype, t)
end

abstract type BoundaryCondition
end

"""
    PeriodicBC <: BoundaryCondition

Assume that the topology is periodic and hence nothing special needs to be done at the boundaries.
"""
struct PeriodicBC <: BoundaryCondition end

# TODO: assert somewhere that the topology is actually periodic when using those
atmos_boundary_state!(_, ::PeriodicBC, _...) = nothing

"""
    NoFluxBC <: BoundaryCondition

Set the momentum at the boundary to be zero.
"""
# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)
struct NoFluxBC <: BoundaryCondition
end

function atmos_boundary_state!(::Rusanov, bc::NoFluxBC, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t, _...)
  FT = eltype(stateM)
  stateP.ρ = stateM.ρ
  stateP.ρu -= 2 * dot(stateM.ρu, nM) * SVector(nM)
end

function atmos_boundary_state!(::CentralNumericalFluxDiffusive, bc::NoFluxBC,
                               m::AtmosModel, stateP::Vars, diffP::Vars,
                               auxP::Vars, nM, stateM::Vars, diffM::Vars,
                               auxM::Vars, bctype, t, _...)
  FT = eltype(stateM)
  stateP.ρ = stateM.ρ
  stateP.ρu -= 2 * dot(stateM.ρu, nM) * SVector(nM)
  diffP.ρτ = SVector(FT(0), FT(0), FT(0), FT(0), FT(0), FT(0))
  diffP.ρd_h_tot =  SVector(FT(0), FT(0), FT(0))
end

"""
    InitStateBC <: BoundaryCondition

Set the value at the boundary to match the `init_state!` function. This is
mainly useful for cases where the problem has an explicit solution.
"""
# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)
struct InitStateBC <: BoundaryCondition
end
function atmos_boundary_state!(::Rusanov, bc::InitStateBC, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t, _...)
  init_state!(m, stateP, auxP, auxP.coord, t)
end
function atmos_boundary_state!(::CentralNumericalFluxDiffusive, bc::InitStateBC,
                               m::AtmosModel, stateP::Vars, diffP::Vars,
                               auxP::Vars, nM, stateM::Vars, diffM::Vars,
                               auxM::Vars, bctype, t, _...)
  init_state!(m, stateP, auxP, auxP.coord, t)
end

"""
  DYCOMS_BC <: BoundaryCondition
  Prescribes boundary conditions for Dynamics of Marine Stratocumulus Case
"""
struct DYCOMS_BC{FT} <: BoundaryCondition
  C_drag::FT
  LHF::FT
  SHF::FT
end
function atmos_boundary_state!(::Rusanov, bc::DYCOMS_BC, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t, state1::Vars, aux1::Vars)
  # stateM is the 𝐘⁻ state while stateP is the 𝐘⁺ state at an interface.
  # at the boundaries the ⁻, minus side states are the interior values
  # state1 is 𝐘 at the first interior nodes relative to the bottom wall
  FT = eltype(stateP)
  # Get values from minus-side state
  ρM = stateM.ρ
  UM, VM, WM = stateM.ρu
  EM = stateM.ρe
  QTM = stateM.moisture.ρq_tot
  uM, vM, wM  = UM/ρM, VM/ρM, WM/ρM
  q_totM = QTM/ρM
  UnM = nM[1] * UM + nM[2] * VM + nM[3] * WM

  # Assign reflection wall boundaries (top wall)
  stateP.ρu = SVector(UM - 2 * nM[1] * UnM,
                      VM - 2 * nM[2] * UnM,
                      WM - 2 * nM[3] * UnM)

  # Assign scalar values at the boundaries
  stateP.ρ = ρM
  stateP.moisture.ρq_tot = QTM

  if bctype == 1 # bctype identifies bottom wall
    stateP.ρu = SVector(0,0,0)
  end
end
function atmos_boundary_state!(::CentralNumericalFluxDiffusive, bc::DYCOMS_BC,
                               m::AtmosModel, stateP::Vars, diffP::Vars,
                               auxP::Vars, nM, stateM::Vars, diffM::Vars,
                               auxM::Vars, bctype, t, state1::Vars, diff1::Vars,
                               aux1::Vars)
  # stateM is the 𝐘⁻ state while stateP is the 𝐘⁺ state at an interface.
  # at the boundaries the ⁻, minus side states are the interior values
  # state1 is 𝐘 at the first interior nodes relative to the bottom wall
  FT = eltype(stateP)
  # Get values from minus-side state
  ρM = stateM.ρ
  UM, VM, WM = stateM.ρu
  EM = stateM.ρe
  QTM = stateM.moisture.ρq_tot
  uM, vM, wM  = UM/ρM, VM/ρM, WM/ρM
  q_totM = QTM/ρM
  UnM = nM[1] * UM + nM[2] * VM + nM[3] * WM

  # Assign reflection wall boundaries (top wall)
  stateP.ρu = SVector(UM - 2 * nM[1] * UnM,
                      VM - 2 * nM[2] * UnM,
                      WM - 2 * nM[3] * UnM)

  # Assign scalar values at the boundaries
  stateP.ρ = ρM
  stateP.moisture.ρq_tot = QTM
  # Assign diffusive fluxes at boundaries
  diffP = diffM
  xvert = auxM.coord[3]

  if bctype == 1 # bctype identifies bottom wall
    # ------------------------------------------------------------------------
    # (<var>_FN) First node values (First interior node from bottom wall)
    # ------------------------------------------------------------------------
    z_FN             = aux1.coord[3]
    ρ_FN             = state1.ρ
    U_FN, V_FN, W_FN = state1.ρu
    E_FN             = state1.ρe
    u_FN, v_FN, w_FN = U_FN/ρ_FN, V_FN/ρ_FN, W_FN/ρ_FN
    windspeed_FN     = sqrt(u_FN^2 + v_FN^2 + w_FN^2)
    q_tot_FN         = state1.moisture.ρq_tot / ρ_FN
    e_int_FN         = E_FN/ρ_FN - windspeed_FN^2/2 - grav*z_FN
    TS_FN            = PhaseEquil(e_int_FN, q_tot_FN, ρ_FN)
    T_FN             = air_temperature(TS_FN)
    q_vap_FN         = q_tot_FN - PhasePartition(TS_FN).liq
    # --------------------------
    # Bottom boundary quantities
    # --------------------------
    zM          = auxM.coord[3]
    q_totM      = QTM/ρM
    windspeed   = sqrt(uM^2 + vM^2 + wM^2)
    e_intM      = EM/ρM - windspeed^2/2 - grav*zM
    TSM         = PhaseEquil(e_intM, q_totM, ρM)
    q_vapM      = q_totM - PhasePartition(TSM).liq
    TM          = air_temperature(TSM)
    # ----------------------------------------------------------
    # Extract components of diffusive momentum flux (minus-side)
    # ----------------------------------------------------------
    ρτM = diffM.ρτ

    # ----------------------------------------------------------
    # Boundary momentum fluxes
    # ----------------------------------------------------------
    # Case specific for flat bottom topography, normal vector is n⃗ = k⃗ = [0, 0, 1]ᵀ
    # A more general implementation requires (n⃗ ⋅ ∇A) to be defined where A is replaced by the appropriate flux terms
    C_drag = bc.C_drag
    ρτ13P  = -ρM * C_drag * windspeed_FN * u_FN
    ρτ23P  = -ρM * C_drag * windspeed_FN * v_FN
    # Assign diffusive momentum and moisture fluxes
    # (i.e. ρ𝛕 terms)
    stateP.ρu = SVector(0,0,0)
    diffP.ρτ = SHermitianCompact{3,FT,6}(SVector(FT(0),ρτM[2,1],ρτ13P, FT(0), ρτ23P,FT(0)))

    # ----------------------------------------------------------
    # Boundary moisture fluxes
    # ----------------------------------------------------------
    diffP.moisture.ρd_q_tot  = SVector(FT(0),
                                       FT(0),
                                       bc.LHF/(LH_v0))
    # ----------------------------------------------------------
    # Boundary energy fluxes
    # ----------------------------------------------------------
    # Assign diffusive enthalpy flux (i.e. ρ(J+D) terms)
    diffP.ρd_h_tot  = SVector(FT(0),
                              FT(0),
                              bc.LHF + bc.SHF)
  end
end

"""
  RayleighBenardBC <: BoundaryCondition

# Fields
$(DocStringExtensions.FIELDS)
"""
struct RayleighBenardBC{FT} <: BoundaryCondition
  "Prescribed bottom wall temperature [K]"
  T_bot::FT
  "Prescribed top wall temperature [K]"
  T_top::FT
end
# Rayleigh-Benard problem with two fixed walls (prescribed temperatures)
function atmos_boundary_state!(::Rusanov, bc::RayleighBenardBC, m::AtmosModel,
                               stateP::Vars, auxP::Vars, nM, stateM::Vars,
                               auxM::Vars, bctype, t,_...)
  # Dry Rayleigh Benard Convection
  @inbounds begin
    FT = eltype(stateP)
    ρP = stateM.ρ
    stateP.ρ = ρP
    stateP.ρu = SVector{3,FT}(0,0,0)
    if bctype == 1
      E_intP = ρP * cv_d * (bc.T_bot - T_0)
    else
      E_intP = ρP * cv_d * (bc.T_top - T_0)
    end
    stateP.ρe = (E_intP + ρP * auxP.coord[3] * grav)
    nothing
  end
end
function atmos_boundary_state!(::CentralNumericalFluxDiffusive, bc::RayleighBenardBC,
                               m::AtmosModel, stateP::Vars, diffP::Vars,
                               auxP::Vars, nM, stateM::Vars, diffM::Vars,
                               auxM::Vars, bctype, t, _...)
  # Dry Rayleigh Benard Convection
  @inbounds begin
    FT = eltype(stateM)
    ρP = stateM.ρ
    stateP.ρ = ρP
    stateP.ρu = SVector{3,FT}(0,0,0)
    if bctype == 1
      E_intP = ρP * cv_d * (bc.T_bot - T_0)
    else
      E_intP = ρP * cv_d * (bc.T_top - T_0)
    end
    stateP.ρe = (E_intP + ρP * auxP.coord[3] * grav)
    diffP.ρd_h_tot = SVector(diffP.ρd_h_tot[1], diffP.ρd_h_tot[2], FT(0))
    nothing
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using DataStructures
using Statistics

"""
PS with PF
----------

Prioritized Sweeping with: ParticleFilterJump for the Tsas' estimation
and Perfect Integration for Rsas' estimation.

See learner/prioritizedsweepinglearner.jl for the core implementation and more info.

Note 1
------
With particle filtering we need to run multiple times with different seeds.
As final LL we report the MEAN across seeds.
For signals/results we keep track of the ones of all seeds
and in the end report the one with the MAXIMUM LL.
Note 2
------
get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl
"""
function getLL_Prio_Sweep_TParticleFilterJump(data::Array{Int,2},
                                                gamma::Float64, temp::Float64,
                                                surpriseprobability::Float64,
                                                stochasticity::Float64,
                                                nr_upd_cycles::Union{Int, Float64},
                                                nparticles::Union{Int, Float64};
                                                particlefilterseeds = [1,2,3],
                                                init_R = 0.,
                                                do_debug_print = false
                                                )
    nr_upd_cycles = converttoInt(nr_upd_cycles)
    nparticles = converttoInt(nparticles)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_Prio_Sweep_TParticleFilter supports only one subject")
    end
    all_puzzles = get_all_puzzles(data)

    LL_seeds = zeros(Float64, length(particlefilterseeds))
    action_selection_probs_seeds = Array{Array{Float64, 1}, 1}(undef, length(particlefilterseeds))
    [action_selection_probs_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    # signals:
    Surprise_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [Surprise_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    γ_surprise_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [γ_surprise_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]

    QSA_seeds = Array{Any, 1}(undef, length(particlefilterseeds))
    Tsas_seeds = Array{Any, 1}(undef, length(particlefilterseeds))
    learner = []
    for i_seed in 1:length(particlefilterseeds)
        action_counter = 0
        for i_puzzle in all_puzzles
            # Puzzle = independent block of the experiment.
            # The code was made for general purposes, in this experiment we
            # have only ONE puzzle.
            learner = PrioritizedSweeping(nr_states = nr_states,
                                        nr_actions = nr_actions,
                                        γ = gamma, init_R = init_R,
                                        nr_upd_cycles = nr_upd_cycles,
                                        Testimatetype = TParticleFilterJump,
                                        nparticles = nparticles,
                                        surpriseprobability = surpriseprobability,
                                        stochasticity = stochasticity,
                                        seedlearner = particlefilterseeds[i_seed],
                                        do_debug_print = do_debug_print
                                        )
            all_epi_in_puzzle = get_all_episodes_in_puzzle(data, i_puzzle)
            for i_episode in all_epi_in_puzzle
                episode_data = get_data_by_puzzle_episode(data, i_puzzle, i_episode)
                path_length = size(episode_data)[1]  # = T
                if do_debug_print
                    println("i_puzzle:$i_puzzle | i_episode:$i_episode | path_length:$path_length")
                end
                for t in 1:path_length
                    action_counter += 1
                    state = episode_data[t, Int(c_state)]
                    action = episode_data[t, Int(c_action)]
                    reward = episode_data[t, Int(c_reward)]
                    next_state = episode_data[t, Int(c_nextState)]
                    done = reward>0 ? true : false
                    # println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | done: $done")
                    if do_debug_print
                        println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | done: $done")
                    end

                    p_actions = getSoftmaxProbs(learner.QSA[state,:], temp, [0., 0.])
                    action_selection_probs_seeds[i_seed][action_counter] = copy(p_actions[action])
                    LL_seeds[i_seed] += log(p_actions[action])
                    if do_debug_print
                        println("p(s,a)=$p_actions, action=$action")
                    end

                    # Prioritized Sweeping with Particle Filtering
                    update!(learner, state, action, next_state, reward, done)

                    Surprise_seeds[i_seed][action_counter] = copy(learner.Testimate.Sgm[1])
                    γ_surprise_seeds[i_seed][action_counter] = copy(learner.Testimate.γ_surprise[1])
                    QSA_seeds[i_seed] = copy(learner.QSA)
                    Tsas_seeds[i_seed]= deepcopy(learner.Testimate.Ps1a0s0)

                end # end t
            end # end i_episode
        end # end i_puzzle
    end # end i_seed

    LL = get_meanLL_acrossseeds(LL_seeds)

    action_selection_probs_cat = hcat(action_selection_probs_seeds...)
    action_selection_probs = mean(action_selection_probs_cat, dims=2)[:]
    trial_by_trial_p_sa = (action_selection_probs)

    # For signals: return the ones corresponding to the Max LL across seeds
    LL_max, seed_idx = findmax(LL_seeds)
    trial_by_trial_p_sa_max = (action_selection_probs_seeds[seed_idx])
    if do_debug_print
        println(finalQSA_record)
        println("action_selection_probs: $action_selection_probs")
        println("LL_mean: $LL_mean")
        println("LL_max: $LL_max")
        println("trial_by_trial_p_sa: $trial_by_trial_p_sa")
    end
    algo = "prio_sweep_particlefilter"
    param_v = [nr_upd_cycles, nparticles, gamma, temp, surpriseprobability, stochasticity]
    param_str = "nr_upd_cycles=$nr_upd_cycles, nparticles=$nparticles,
                gamma=$(round(gamma,digits=4)), temp=$(round(temp,digits=4)),
                surpriseprobability=$(round(surpriseprobability,digits=4)),
                stochasticity=$(round(stochasticity,digits=4))"
    results_dict = Dict{String,Any}("final_QSA" => QSA_seeds[seed_idx],
                                    "final_Tsas" => Tsas_seeds[seed_idx],
                                    "Surprise" => Surprise_seeds[seed_idx],
                                    "gamma_surprise" => γ_surprise_seeds[seed_idx],
                                    "nr_upd_cycles" => nr_upd_cycles,
                                    "nparticles" => nparticles,
                                    "particlefilterseeds" => particlefilterseeds,
                                    "trial_by_trial_p_sa_max" => trial_by_trial_p_sa_max)

    signals = Step_by_step_signals(data_stats.all_subject_ids[1],
    [], #RPE
    [], #SPE
    trial_by_trial_p_sa,
    LL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
    return LL, signals
end


text:
using ColorPrinting
using Documenter

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

makedocs(;
    modules=[ColorPrinting],
    authors="Andrey Oskin",
    repo="https://github.com/Arkoniak/ColorPrinting.jl/blob/{commit}{path}#{line}",
    sitename="ColorPrinting.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Arkoniak.github.io/ColorPrinting.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Arkoniak/ColorPrinting.jl",
    devbranch="main",
)


text:
using Test
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using jInv.LinearSolvers

println("\t\t-- test ConvDiffFEMParam --")
domain = [-1 1. -1 1]
n      = [32 32]
M      = getRegularMesh(domain,n)
Ainv   = getJuliaSolver()
pFor   = getConvDiffFEMParam(M,sig=.0001,Ainv=Ainv)
m0     = zeros(tuple(n...))
m0[10:20,15:25] .=1.;
isOK, his = checkDerivative(vec(m0),pFor,out=true)
@test isOK

v = randn(M.nc)
w = randn(prod(M.n.+1))
t1 = dot(w, getSensMatVec(v,vec(m0),pFor))
t2 = dot(v, getSensTMatVec(w,vec(m0),pFor))
@test abs(t1-t2)/abs(t1) < 1e-8

domain = [-1 1. 0 1 -1 1]
n      = [6 8 10]
M      = getRegularMesh(domain,n)
clear!(Ainv)
pFor   = getConvDiffFEMParam(M,sig=.5,Ainv=Ainv)

m0              = zeros(tuple(n...))
m0[3:4,4:6,4:7] .=1.;
isOK, his = checkDerivative(vec(m0),pFor,out=true)
@test isOK

v = randn(M.nc)
w = randn(prod(M.n.+1))
t1 = dot(w, getSensMatVec(v,vec(m0),pFor))
t2 = dot(v, getSensTMatVec(w,vec(m0),pFor))
@test abs(t1-t2)/abs(t1) < 1e-8

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
# This file was generated, do not modify it. # hide
coerce!(df, autotype(df, :few_to_finite))
schema(df)


text:
using Documenter, LinearAlgebra, SuiteSparse

makedocs(
    format = :html,
    sitename = "SuiteSparse.jl",
    modules = [SuiteSparse.CHOLMOD, SuiteSparse.UMFPACK, SuiteSparse.SPQR],
    pages = [
        "index.md",
        "cholmod.md",
        "umfpack.md",
        "spqr.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaLinearAlgebra/SuiteSparse.jl.git",
    target = "build",
    julia  = "nightly",
    deps = nothing,
    make = nothing
)


text:
# macro for argument checking (taken from Distributions.jl)

macro check_args(D, cond)
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)), " is not satisfied.")))
        end
    end
end

"""
    BIC(fm::ContinuousUnivariateDistribution, data::Array{<:Real,1}; censoring::Real=0)

Calculates the Bayesian information criterion (BIC) of the specified model.

```math
\\textrm{BIC} = \\ln f_{(Y|\\hat{\\theta}_{MV})}(y) - \\frac{k}{2}\\ln{n}
```

"""
function BIC(fm::ContinuousUnivariateDistribution, data::Array{<:Real,1}; censoring::Real=0)
    n = size(data,1)
    k = size(params(fm),1)

    r_data = data[data .>= censoring]
    l_data = fill(censoring, count(data .< censoring))

    ll = sum(logcdf.(fm, l_data)) + sum(logpdf.(fm, r_data))

    return ll - k/2*log(n)
end

function BIC(fm::ContinuousUnivariateDistribution, data::TimeArray; censoring::Real=0)
    return BIC(fm, values(data), censoring=censoring)
end

function window(ts::TimeArray, year_from::Int64, year_to::Int64)
    return to(from(ts, DateTime(year_from, 01, 01)), DateTime(year_to-1, 12, 31))
end


text:
@testset "GetIndex(Constant)" begin
    func = Constant([1.0 2.0 3.0; 4.0 5.0 6.0])[:, 2]
    @test typetree(func) == [:Constant]
end

@testset "GetIndex(Inflate)" begin
    func = Inflate(Constant([1.0, 2.0, 3.0]), (5,), 1 => [1, 2, 3])[:]
    @test typetree(func) == [:Inflate, [:Constant], [:Constant]]

    func = Inflate(Constant([1.0, 2.0, 3.0]), (3,))[:]
    @test typetree(func) == [:Inflate, [:Constant]]

    func = Inflate(Constant([1.0, 2.0, 3.0]), (3,))[1]
    @test typetree(func) == [:Inflate, [:Constant]]
end

@testset "InsertAxis(Constant)" begin
    func = InsertAxis(Constant([1.0, 2.0, 3.0]), 1)
    @test typetree(func) == [:Constant]
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
function demoscaletoChebyinterval()

    # ground truth: we know u for x[i] is i.
    x = LinRange(-1,1,5)
    println(x)

    c = collect(1:5)
    u = coordinatetolattice1D(collect(x),x[1],x[end], length(x))

    x_rec = lattice1Dtocoordinate(u,x[1],x[end], length(x))
    println("u is ", u)
    println("x_rec is ", xx)
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
# Michelle Medina
# July 12, 2020
# Utilities for Nodal / Hierarchical Bases

export calcNodal, getMaximum, evaluate_upto

function calcNodal(u::Function, x::Array{T,1})::NodalBasis{T} where {T}
	coefficients = zeros(length(x))
 	for index in CartesianIndices(x)
 		coefficients[index] = u(x[index])
 	end
 	return NodalBasis(coefficients)
end

function getMaximum(values::HierarchicalBasis{T})::Array{T,1} where {T}
	M = []
	for l in 1:maxlevel(values)
		absolute = abs.(values.levels[l].coefficients)
		max = maximum(absolute)
		append!(M,max)
	end
	return M
end

function evaluate_upto(H::HierarchicalBasis{T}, level::Int, x::Float64)::T where {T}
	value = T(0)
	for l in 0:level
		for j in 0:2^l
			value = value + H[l, j]*basis(l, j, x)
		end
	end
	return value
end

function evaluate_upto(H::HierarchicalBasis{T}, level::Int, x::Array{T,1})::Array{T,1} where {T}
	A = zeros(size(x))
	for index in CartesianIndices(x)
		A[index] = evaluate_upto(H, level, x[index])
	end
	return A
end


text:
function FESpaces._change_argument(
  op::typeof(jacobian),f,trian,uh::MultiFieldFEFunction)

  U = get_fe_space(uh)
  function g(cell_u)
    single_fields = GenericCellField[]
    nfields = length(U.spaces)
    for i in 1:nfields
      cell_values_field = lazy_map(a->a.array[i],cell_u)
      cf = CellField(U.spaces[i],cell_values_field)
      cell_data = lazy_map(BlockMap((1,nfields),i),get_data(cf))
      uhi = GenericCellField(cell_data,get_triangulation(cf),DomainStyle(cf))
      push!(single_fields,uhi)
    end
    xh = MultiFieldCellField(single_fields)
    cell_grad = f(xh)
    get_contribution(cell_grad,trian)
  end
  g
end

function FESpaces._change_argument(
  op::typeof(gradient),f,trian,uh::MultiFieldFEFunction)

  U = get_fe_space(uh)
  function g(cell_u)
    single_fields = GenericCellField[]
    nfields = length(U.spaces)
    for i in 1:nfields
      cell_values_field = lazy_map(a->a.array[i],cell_u)
      cf = CellField(U.spaces[i],cell_values_field)
      cell_data = lazy_map(BlockMap(nfields,i),get_data(cf))
      uhi = GenericCellField(cell_data,get_triangulation(cf),DomainStyle(cf))
      push!(single_fields,uhi)
    end
    xh = MultiFieldCellField(single_fields)
    cell_grad = f(xh)
    get_contribution(cell_grad,trian)
  end
  g
end

function Algebra.hessian(f::Function,uh::MultiFieldFEFunction)
  @notimplemented
end

#function FESpaces._change_argument(
#  op::typeof(hessian),f,trian,uh::MultiFieldFEFunction)
#
#  U = get_fe_space(uh)
#  function g(cell_u)
#    single_fields = GenericCellField[]
#    nfields = length(U.spaces)
#    for i in 1:nfields
#      cell_values_field = lazy_map(a->a.array[i],cell_u)
#      cf = CellField(U.spaces[i],cell_values_field)
#      cell_data = lazy_map(BlockMap((nfields,nfields),i),get_data(cf))
#      uhi = GenericCellField(cell_data,get_triangulation(cf),DomainStyle(cf))
#      push!(single_fields,uhi)
#    end
#    xh = MultiFieldCellField(single_fields)
#    cell_grad = f(xh)
#    get_contribution(cell_grad,trian)
#  end
#  g
#end


text:
using Test
using NumberIntervals
using Turmeric
using Base.Threads
using StaticArrays

import Turmeric: strict_isinterior, contains_root

@testset "guarded_reduce()" begin
    nothing
end

@testset "explore_region!()" begin
    nothing
end

@testset "ThreadedRootSearch iterator interface" begin
    nothing
end

@testset "pop_unfinished!()" begin
    nothing
end

@testset "collect_regions()" begin
    nothing
end

@testset "estimate_generations()" begin
    nothing
end

@testset "chunks(regions, chunksize)" begin
    nothing
end

@testset "bisect()" begin
    x1 = NumberInterval(-1., 1.)
    x2 = NumberInterval(1., 10.)
    c = SVector(x1, x2)
    a, b = bisect(c)
    @test all(a .∪ b .=== c)
    @test !all(a .=== c)
end

@testset "contains_root()" begin
    x1 = NumberInterval(-1., 1.)
    x2 = NumberInterval(1., 10.)
    @test contains_root(x1)
    @test !contains_root(x2)
    @test contains_root(SVector(x1, x1))
    @test !contains_root(SVector(x1, x2))
end

@testset "isinterior()" begin
    x1 = NumberInterval(1. ,2.)
    x2 = NumberInterval(0., 5.)
    @test isinterior(SVector(x1, x1), SVector(x2, x2))
    @test !isinterior(SVector(x1, x1), SVector(x2, x1))
end

@testset "strict_isinterior()" begin
    x = entireinterval()
    @test !strict_isinterior(x, x)
end

@testset "diam()" begin
    x1 = NumberInterval(1. ,2.)
    x2 = NumberInterval(-3., 5.)
    @test diam(SVector(x1, x2)) >= max(diam(x1), diam(x2))
    x3 = NumberInterval(1)
    @test diam(SVector(x3, x3)) == 0
end

@testset "Krawczyk contractor" begin
    nothing
end

@testset "GradientContractor" begin
    x = NumberInterval(-10.0, -1.0)
    contractor = GradientContractor(log, Newton(), x)
    @test_throws DomainError contractor(x)
end

@testset "trivial contractor" begin
    nothing
end

@testset "fullcontract()" begin
    nothing
end


text:
using Agents, Test

mutable struct Agent1 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
end
model1 = ABM(Agent1, Space((3,3)))

agent = add_agent!((1,1), model1)
@test agent.pos == (1, 1)
@test agent.id == 1
pos1 = model1.space.agent_positions[coord2vertex((1,1), model1)]
@test length(pos1) == 1
@test pos1[1] == 1

move_agent!(agent, (2,2), model1)

@test agent.pos == (2,2)
pos1 = model1.space.agent_positions[coord2vertex((1,1), model1)]
@test length(pos1) == 0
pos2 = model1.space.agent_positions[coord2vertex((2,2), model1)]
@test pos2[1] == 1

# %% Scheduler tests
N = 1000
mutable struct Agent0 <: AbstractAgent
  id::Int
end
# fastest
model = ABM(Agent0)
for i in 1:N; add_agent!(model); end
@test sort!(collect(keys(model.agents))) == 1:N

# by_id
model = ABM(Agent0; scheduler = by_id)
for i in 1:N; add_agent!(model); end

@test sort!(collect(keys(model.agents))) == 1:N

@test model.scheduler(model) == 1:N

using Random
Random.seed!(12)

# random
model = ABM(Agent0; scheduler = random_activation)
for i in 1:N; add_agent!(model); end

@test model.scheduler(model)[1:3] == [913, 522, 637] # reproducibility test

# partial
Random.seed!(12)
model = ABM(Agent0; scheduler = partial_activation(0.1))
for i in 1:N; add_agent!(model); end

a = model.scheduler(model)
@test length(a) < N
@test a[1] == 74 # reproducibility test

# by property
mutable struct Agentw <: AbstractAgent
  id::Int
  w::Float64
end
model = ABM(Agentw; scheduler = property_activation(:w))
for i in 1:N; add_agent!(model, rand()/rand()); end

Random.seed!(12)
a = model.scheduler(model)

ids = collect(keys(model.agents))
properties = [model.agents[id].w for id in ids]

@test ids[sortperm(properties)] == a


text:
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4

##

dense_coefs = zeros(9,9,9)
dense_coefs[1:9,1:9,5] .= 1
for i = parse(Int,ARGS[1])
    global dense_coefs
    for x = (i-1):-1:0
        for y = (i-1):-1:0
                dense_coefs[1+x,1+y,5] = 0
                dense_coefs[1+x,9-y,5] = 0
                dense_coefs[9-x,1+y,5] = 0
                dense_coefs[9-x,9-y,5] = 0
        end
    end
    @show sum(dense_coefs .== 1)/81
    st_def = CreateStencilDefinition(dense_coefs)
    bdim = 32
    if length(ARGS) == 2
        bdim = 16
    end
    st_inst1 = NewStencilInstance(st_def, bdim=bdim, m_step=false)
    println(st_inst1.kernel_expr)
    radius = 4
    nx = 12
    ny = 12
    nz = 0
    dx = 1<<(nx)
    dy = 1<<(ny)
    dz = 1<<(nz)

    data = CreateData(dx,dy,dz)

    t_steps = 16
    gpu_out = ApplyStencil(st_inst1, data, t_steps)
    CUDA.cuProfilerStart()
    NVTX.@range "standard" begin
        gpu_out = ApplyStencil(st_inst1, data, t_steps) end
    #exit()

end


text:
#!/usr/bin/env julia

const DATA_DIR = "../09_data"
const IS_TEST = false

if IS_TEST
    const INPUT_FILE    = "input_test.txt"
    const PREAMBLE_SIZE = 5
else
    const INPUT_FILE    = "input.txt"
    const PREAMBLE_SIZE = 25
end

macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end

function does_contain_distinct_sum(window, val)
    for i in 1:length(window)-1
        for j in i+1:length(window)
            if val == window[i] + window[j]
                return (true, i, j)
            end
        end
    end
    return (false, 0, 0)
end

function get_encryption_weakness(vals, val)
    for width in 2 : length(vals) - 1
        for start in 1 : length(vals) - width + 1
            # println("\t(width, start) = ($(width), $(start))")
            src_range = vals[start : start + width - 1]
            if sum(src_range) == val
                return (true, start, width)
            end
        end
    end
    return (false, 0, 0)
end

function main(filename::String, preamble_size::Int64)
    lines = readlines(filename)
    vals = map(s->parse(Int64, s), lines)
    for k in PREAMBLE_SIZE + 1 : 1 : length(vals)
        # println("size of vals=$(length(vals)), k=$(k)")
        (is_valid, i, j) = does_contain_distinct_sum(
            view(vals, k-PREAMBLE_SIZE : k), vals[k]
            )
        if is_valid
            val_i = vals[k - PREAMBLE_SIZE + i]
            val_j = vals[k - PREAMBLE_SIZE + j]
            # println("At position $(k): $(vals[k]) = $(val_i) + $(val_j)")
        else
            println("========== PART 1 ==========")
            println("At position $(k): $(vals[k]) is not a distinct sum from the lagging window")
            (did_find_sum, start, width) = get_encryption_weakness(vals, vals[k])
            println()
            println("========== PART 2 ==========")
            if did_find_sum
                items = @view vals[start : start + width - 1]
                mn = minimum(items)
                mx = maximum(items)
                println("Target is the sum of numbers ranging from $(mn) to $(mx)")
                println("items=$(items)")
                println("$(mn) + $(mx) = $(mn + mx)")
            else
                println("Did not find range that sums to target")
            end
            return 0
        end
    end
    println("All values after preamble are distinct sums from lagging window")
end

main("$DATA_DIR/$INPUT_FILE", PREAMBLE_SIZE)

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
"""
    multinomial_hyper(α::AbstractArray{Float32,1})

[Dirichlet Distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution)
"""
struct multinomial_hyper <: distribution_hyper_params
    α::AbstractArray{Float32,1}
end

mutable struct multinomial_sufficient_statistics <: sufficient_statistics
    N::Float32
    points_sum::AbstractArray{Float32,1}
end

function calc_posterior(prior:: multinomial_hyper, suff_statistics::Vector{Tuple{sufficient_statistics,Number}})
    if length(suff_statistics) == 0
        return prior
    end
    points_num = sum([post_kernel(x[2],global_time)*x[1].N for x in suff_statistics])
    if points_num == 0
        return prior
    end

    suff_sum = sum([post_kernel(x[2],global_time)*x[1].points_sum for x in suff_statistics],dims=1)[1]
    return multinomial_hyper(prior.α + suff_sum)
end

function sample_distribution(hyperparams::multinomial_hyper)
    return multinomial_dist(log.(rand(Dirichlet(Float64.(hyperparams.α)))))
end

function create_sufficient_statistics(hyper::multinomial_hyper,posterior::multinomial_hyper,points::AbstractArray{Float32,2}, pts_to_group = 0)
    pts = copy(points)
    points_sum = sum(pts, dims = 2)[:]
    S = pts * pts'
    return multinomial_sufficient_statistics(size(points,2),points_sum)
end

function log_marginal_likelihood(hyper::multinomial_hyper, posterior_hyper::multinomial_hyper, suff_stats::Vector{Tuple{sufficient_statistics,Number}})
    D = length(suff_stats[1][1].points_sum)
    logpi = log(pi)
    val = logabsgamma(sum(hyper.α))[1] -logabsgamma(sum(posterior_hyper.α))[1] + sum((x-> logabsgamma(x)[1]).(posterior_hyper.α) - (x-> logabsgamma(x)[1]).(hyper.α))
    return val
end

function aggregate_suff_stats(suff_l::multinomial_sufficient_statistics, suff_r::multinomial_sufficient_statistics)
    return multinomial_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum)
end

function posterior_predictive!(r::AbstractArray,x::AbstractMatrix,posterior_hyper::multinomial_hyper)
    v_sum = log.(posterior_hyper.α/sum(posterior_hyper.α))
    r .= (v_sum' * x)[1,:]
end


text:
using Test

import TLOU

const pr = TLOU.Pricing(
    0.6, # K
    (0.1,0.4,0.6,0.8), #CL
    (0.2,0.5,0.7,0.9, 1.0), #CH
    (0.9,0.7,0.5,0.3), #piL
    (1.1,1.2,1.4,1.6, 2.0)  #piH
)

const test_table = [
    (0.05, 0, 0),
    (0.15, 1, 0),
    (0.20, 1, 1),
    (0.40, 2, 1),
    (0.50, 2, 2),
    (0.60, 3, 2),
    (0.90, 4, 4),
    (1.00, 4, 5),
]

@testset "TLOU.Pricing indexing" begin
    @test pr isa TLOU.Pricing{4,5}
    for (c, idxL, idxH) in test_table
        (iL, iH) = TLOU.price_segment(pr, c, :both)
        @test iL == idxL
        @test iH == idxH
        (pL, pH) = TLOU.price_value(pr, c, :both)
        @test (iL == 0 && pL ≈ 1.0) || pr.piL[iL] ≈ pL
        @test (iH == 0 && pH ≈ 1.0) || pr.piH[iH] ≈ pH
    end
end

@testset "Pricing computation" begin
    for (c, idxL, idxH) in test_table
        xlow =  [0.0, 0.5 * c, 0.9 * c, 1.0 * c]
        xhigh = [1.05c, 1.5 * c, 1.9 * c, 3.0 * c]
        for x in xlow
            cost = pr(c, x)

            @test cost ≈ pr.K*c + x * TLOU.price_value(pr, c, Val{:low}())
        end
        for x in xhigh
            cost = pr(c, x)
            @test cost ≈ pr.K*c + x * TLOU.price_value(pr, c, Val{:high}())
        end
    end
end


text:
function CHISQUAREHEUR(trans::SparseMatrixCSC{Float64, Int64},
                       capacity::Array{Int64,1},
                       Q::Array{Int64,2},
                       sig::Float64)
    if CHECKCAPACITY(Q,capacity) == 1
        # Number of transactions
        J = size(trans,1)

        # Number of products
        I = size(trans,2)

        # Chi Square Test
        ## We start the chi-square heuristic by performing chi-square tests
        ## for all unique SKU-combinations. For each undirected SKU-pair we
        ## evaluate a chi-square test and either accept or reject the null
        ## hypothesis of independence. Our goal is a matrix dep with the
        ## “net-benefit” of all unique positive dependent SKU-combinations
        ## and a matrix nor with all "independent contributions" of the unique
        ## SKU-combinations. More details in our article.
        sum_cond_sku = Int.(vec(sum(trans,dims=1)))
        dep, nor = HYOPTHESISCHI(Q::Array{Int64,2},
                                 I::Int64,
                                 J::Int64,
                                 sig::Float64,
                                 sum_cond_sku::Array{Int64,1})

        ## Create the sku-warehouse allocation matrix
        X = Array{Bool,2}(undef,I,size(capacity,1)) .= false

        ## Determine the sum of the coappearances for each SKU on the base of
        ## the matrices dep and nor
        sum_dep = vec(sum(dep,dims = 1))
        sum_nor = vec(sum(nor,dims = 1))

        ## Calculate the weight of each warehouse. It shows us the density of
        ## the independent coappearances in each warehouse if we were to allocate
        ## all SKUs simply according to the highest independent coappearances.
        weight = WHWEIGHT(capacity::Array{Int64,1},sum_nor::Array{Float64,1})

        ## Copy the capacity to have an array that keeps the left-over capacity
        cap_left = copy(capacity)

        ## In the heuristic we compare the split-delivery minimising potential
        ## of SKU allocations to maximise dependent coappearances (SKUs with
        ## positive dependencies are stored in the same warehouse) with the
        ## maximisation of independent coappearances (SKUs with high independent
        ## coappearances are stored in warehouses with high capacities).
        ## Start the allocation:
        ## If all SKUs are allocated once, continue after the while loop, else
        while sum(X) < I
        ## select the SKU i with the highest coappearance not being allocated
        ## to the warehouses yet. In addition, select the first sorted
        ## warehouse k which has unused capacity left.
        ## Then, check for each warehouse with free space whether SKU i has
        ## significant dependencies to other already allocated SKUs.
        ## If dependencies exist check whether they are expected to dominate.
        ## If so, allocate SKU i to to the warehouse with the dependencies
        ## to complement the positive dependent SKUs in the warehouse.
        ## Otherwise allocate SKU i to warehouse k as the independent
        ## coappearances are expected to dominate.
        ## If no dependencies to yet allocated SKUs exist, select the warehouse
        ## with the highest unused capacity. If the dependet coappearances are
        ## expected to to dominate later, we allocate SKU i to the warehouse
        ## with the highest unused capacity. This maximises the likelihood to
        ## allocate all significant SKUs from a SKU-cluster into one warehouse.
        ## Otherwise we allocate it to warehouse k to maximise the independent
        ## coappearances in the allocation.
        i,k  = SELECTIK(sum_dep::Array{Float64,1},
                            sum_nor::Array{Float64,1},
                            weight,
                            cap_left::Array{Int64,1},
                            X::Array{Bool,2},
                            dep::Array{Float64,2})
        ALLOCATEONE!(X::Array{Bool,2},
                        sum_dep::Array{Float64,1},
                        sum_nor::Array{Float64,1},
                        cap_left::Array{Int64,1},
                        i::Int64,
                        k::Int64)

        ## Check for all unallocated SKUs whether they have positive
        ## dependencies to the SKUs in the warehouse k the last SKU was allocated
        ## to. If so, check whether the dependencies are expected to dominate the
        ## independent coapperances. If yes, allocate the corresponding SKUs to
        ## the warehouse k.
        ADDDEPENDENT!(X::Array{Bool,2},
                        cap_left::Array{Int64,1},
                        k::Int64,
                        dep::Array{Float64,2},
                        nor::Array{Float64,2},
                        sum_dep::Array{Float64,1},
                        sum_nor::Array{Float64,1})
        FILLLAST!(X::Array{Bool,2},cap_left::Array{Int64,1})
        end
        if sum(cap_left) > 0
            ## First, remove every so far assigned dependent SKU-pair from the coappearance
            ## matrix dep to prevent the allocation bias described in our article.
            Qs = REMOVEALLOC(X::Array{Bool,2},
                            cap_left::Array{Int64,1},
                            nor::Array{Float64,2},
                            dep::Array{Float64,2})

            ## Afterwards allocate the SKUs with the highest potential allocation value to
            ## each warehouse with leftover storage space until it is full. If no SKU
            ## with coappearances is found and there is still storage space left, terminate the
            ## algorithm as further allocations pose no benefit.
            FILLUP!(X::Array{Bool,2},
                        Qs::Array{Float64,2},
                        cap_left::Array{Int64,1})
        end
    end
    ## return the resulting allocation matrix
    X = convert(Matrix{Int64},X)
    return X::Array{Int64,2}
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
using Krylov, Stopping

include("random_coordinate_descent_method.jl")
include("stop_random_coordinate_descent_method.jl")
include("stop_random_coordinate_descent_method_LS.jl")

n = 5000;
A = rand(n, n)
sol = rand(n)
b = A * sol

#@time stp = StopRandomizedCD(A,b, max_iter = 1000)

x0 = zeros(size(A, 2))
pb = LinearSystem(A, b)
s0 = GenericState(x0, similar(b))
mcnt = Stopping._init_max_counters_linear_operators()
#
@time meta = StoppingMeta(
  max_cntrs = mcnt,
  atol = 1e-7,
  rtol = 1e-15,
  max_iter = 99,
  retol = false,
  tol_check = (atol, rtol, opt0) -> (atol + rtol * opt0),
  optimality_check = (pb, state) -> state.res,
)
s0 = GenericState(x0, similar(b))
@time stp = LAStopping(pb, meta, cheap_stop_remote_control(), s0)
@time StopRandomizedCD(stp)

@time x, OK = RandomizedCD(A, b, max_iter = 100)

###############################################################################
#
# We compare the stopping_randomized_CD with the same version with an
# artificial substopping created and called at each iteration.
# It appears that the expansive part is to create the SubStopping, and in
# particular create the StoppingMeta.
#
@time meta = StoppingMeta(
  max_cntrs = mcnt,
  atol = 1e-7,
  rtol = 1e-15,
  max_iter = 99,
  retol = false,
  tol_check = (atol, rtol, opt0) -> (atol + rtol * opt0),
  optimality_check = (pb, state) -> state.res,
)
s0 = GenericState(x0, similar(b))
@time stp2 = LAStopping(pb, meta, cheap_stop_remote_control(), s0)
@time StopRandomizedCD_LS(stp2)

;


text:
using LocalFunctionApproximation
using Test
using NearestNeighbors
using StaticArrays

# Only testing NNFA as the GIFA is just a wrapper
# for the already tested GridInterpolations package
points = [SVector(0.,0.), SVector(0.,1.), SVector(1.,1.), SVector(1.,0.)]
nntree = KDTree(points)
vals = [1., 1., -1., -1]
k = 2
r = 0.5*sqrt(2)

knnfa = LocalNNFunctionApproximator(nntree, points, k)
set_all_interpolating_values(knnfa, vals)

@test n_interpolating_points(knnfa) == 4

for i = 1:10
    point = [rand()/2, 0.5]
    @test compute_value(knnfa, point) == 1.0
end

for i = 1:10
    point = [1.0 - rand()/2., 0.5]
    @test compute_value(knnfa, point) == -1.0
end

rnnfa = LocalNNFunctionApproximator(nntree, points, r)
set_all_interpolating_values(rnnfa,vals)
for i = 1:10
    point = [rand()/2, 0.5]
    @test compute_value(rnnfa, point) == 1.0
end

for i = 1:10
    point = [1.0 - rand()/2., 0.5]
    @test compute_value(rnnfa, point) == -1.0
end