text:
using BQCESubroutine.Utilities
using BenchmarkTools
using LoopVectorization
using PaddedMatrices
using YaoLocations

function subspace_mul!(st::AbstractMatrix, comspace, U::AbstractMatrix{T}, subspace) where T
    y = similar(st, (size(U, 1), ))
    indices = [b + 1 for b in comspace]::Vector{Int}
    idx = similar(indices)
    @inbounds for k in subspace, b in 1:size(st, 2)
        for i in 1:size(U, 1)
            idx[i] = k + indices[i]
        end

        for i in 1:size(U, 1)
            y[i] = zero(T)
            for j in 1:size(U, 2)
                y[i] += U[i, j] * st[idx[j], b]
            end
        end

        for i in 1:size(U, 1)
            st[idx[i], b] = y[i]
        end
    end
    return st
end

function subspace_mul_avx!(st::AbstractMatrix, comspace, U::AbstractMatrix{Complex{T}}, subspace) where T
    indices = StrideArray{Int}(undef, (PaddedMatrices.static_length(comspace), ))
    @simd ivdep for i in eachindex(indices)
        indices[i] = comspace[i] + 1
    end

    C_re = StrideArray{T}(undef, (PaddedMatrices.static_length(indices), StaticInt{32}()))
    C_im = StrideArray{T}(undef, (PaddedMatrices.static_length(indices), StaticInt{32}()))
    U_re = StrideArray{T}(undef, (PaddedMatrices.size(U, StaticInt{1}()), PaddedMatrices.size(U, StaticInt{2}())))
    U_im = StrideArray{T}(undef, (PaddedMatrices.size(U, StaticInt{1}()), PaddedMatrices.size(U, StaticInt{2}())))
    Ur = reinterpret(reshape, T, U)

    @inbounds @simd ivdep for i ∈ eachindex(U_re)
        U_re[i] = Ur[1,i]
        U_im[i] = Ur[2,i]
    end
    str = reinterpret(reshape, T, st)
    Bmax = size(st,2)
    for k in subspace
        # idx = k .+ indices
        # _k = k - 1
        for _b ∈ 0:(Bmax-1) >>> 5
            b =    _b << 5;
            bmax = b + 32
            if bmax ≤ Bmax # full static block
                @avx for n ∈ 1:32, m ∈ axes(U, 1)
                    C_re_m_n = zero(T)
                    C_im_m_n = zero(T)
                    for i ∈ axes(U, 2)
                        str_i = k + indices[i]
                        C_re_m_n += U_re[m,i] * str[1,str_i,n+b] - U_im[m,i] * str[2,str_i,n+b]
                        C_im_m_n += U_re[m,i] * str[2,str_i,n+b] + U_im[m,i] * str[1,str_i,n+b]
                    end
                    C_re[m,n] = C_re_m_n
                    C_im[m,n] = C_im_m_n
                end
                @avx for n ∈ 1:32, m ∈ axes(U, 1)
                    str_m = k + indices[m]
                    str[1,str_m,n+b] = C_re[m,n]
                    str[2,str_m,n+b] = C_im[m,n]
                end
                # AmulB!(C_re, C_im, U_re, U_im,
            else # dynamic block
                Nmax = 32 + Bmax - bmax
                @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
                    C_re_m_n = zero(T)
                    C_im_m_n = zero(T)
                    for i ∈ axes(U, 2)
                        str_i = k + indices[i]
                        C_re_m_n += U_re[m,i] * str[1,str_i,n+b] - U_im[m,i] * str[2,str_i,n+b]
                        C_im_m_n += U_re[m,i] * str[2,str_i,n+b] + U_im[m,i] * str[1,str_i,n+b]
                    end
                    C_re[m,n] = C_re_m_n
                    C_im[m,n] = C_im_m_n
                end
                @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
                    str_m = k + indices[m]
                    str[1,str_m,n+b] = C_re[m,n]
                    str[2,str_m,n+b] = C_im[m,n]
                end
            end
        end
    end
    return st
end

function subspace_mul4x4!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        for b in 1:size(st, 2)
            Base.Cartesian.@nexprs 4 i -> begin
                y_i = zero(T)
                Base.Cartesian.@nexprs 4 j -> begin
                    y_i += U_i_j * st[idx_j, b]
                end
            end

            Base.Cartesian.@nexprs 4 i -> begin
                st[idx_i, b] = y_i
            end
        end
    end
    return st
end

function subspace_mul8x8!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    Base.Cartesian.@nextract 8 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 8 U i->begin
        Base.Cartesian.@nextract 8 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 8 idx i-> k + indices_i + offset

        for b in 1:size(st, 2)
            Base.Cartesian.@nexprs 8 i -> begin
                y_i = zero(T)
                Base.Cartesian.@nexprs 8 j -> begin
                    y_i += U_i_j * st[idx_j, b]
                end
            end

            Base.Cartesian.@nexprs 8 i -> begin
                st[idx_i, b] = y_i
            end
        end
    end
    return st
end

using StaticArrays
n = 8
locs = Locations((1, 3))
st = rand(ComplexF64, 1 << n, 100)
U = rand(ComplexF64, 1 << 2, 1 << 2)
SU = @SMatrix rand(ComplexF64, 1 << 2, 1 << 2)
subspace = bsubspace(n, locs)
comspace = bcomspace(n, locs)

subspace_mul4x4!(deepcopy(st), comspace, U, subspace) ≈ subspace_mul!(deepcopy(st), comspace, U, subspace)
subspace_mul8x8!(deepcopy(st), comspace, U, subspace) ≈ subspace_mul!(deepcopy(st), comspace, U, subspace)
subspace_mul_avx!(deepcopy(st), comspace, U, subspace) ≈ subspace_mul!(deepcopy(st), comspace, U, subspace)

@benchmark subspace_mul4x4!($(deepcopy(st)), $comspace, $SU, $subspace)
@benchmark subspace_mul8x8!($(deepcopy(st)), $comspace, $SU, $subspace)
@benchmark subspace_mul!($(deepcopy(st)), $comspace, $SU, $subspace)
@benchmark subspace_mul_avx!($(deepcopy(st)), $comspace, $SU, $subspace)


text:
module TestLib
  include("lens.jl")

  module Omega
    include("omega/space.jl")
  end
end

ext: jl


text:
# Copyright (c) 2019 Harshavardhan P K
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

function getcoords(arr::Array{Tuple{Char,Int64}})::Set{Tuple{Int64,Int64}}
    coords::Set{Tuple{Int64,Int64}} = Set([])
    x, y = 0, 0
    for (dir, steps) in arr
        for _ ∈ 1:steps
            if dir == 'R'
                x += 1
            elseif dir == 'L'
                x -= 1
            elseif dir == 'U'
                y += 1
            elseif dir == 'D'
                y -= 1
            end
            push!(coords, (x, y))
        end
    end
    return coords
end

function closestcoods(arr1::Set{Tuple{Int64,Int64}}, arr2::Set{Tuple{Int64,Int64}})::Int64
    return minimum([abs(x) + abs(y) for (x, y) in arr1 ∩ arr2])
end

function getData(fname::String)
    f = open(fname)
    line = readline(f)
    line = split(strip(line), ',')
    w1 = map(x -> (x[1], parse(Int64, x[2:end])), line)
    line = readline(f)
    line = split(strip(line), ',')
    w2 = map(x -> (x[1], parse(Int64, x[2:end])), line)
    return w1, w2
end

w1, w2 = getData("input")
c1, c2 = getcoords(w1), getcoords(w2)
println(closestcoods(c1, c2))


text:
module AlarmMerge

    include( "AlrmReport.jl")
    using Printf
    mutable struct AlarmInfo
        AlarmID::String
        S1::UInt
        S2::UInt
        SC::UInt
        t0::UInt
        tfresh::UInt
        Tc::UInt
        function AlarmInfo()
            new()
        end
    end

    AlarmInfoHistory = Dict()
    MergeWidth = 0 #In Grain or Seg
    MergeTime  = 0

    ####################################################################
    function FreshAlarmInfoDict( )
        for key1   in keys(AlarmInfoHistory)
            data  = AlarmInfoHistory[ key1 ]
            if data.tfresh - curtime >= 30*60 # 至少保存30分钟的历史数据
                delete!( AlarmInfoHistory ,key1 )
            end
        end
    end

    function getAlarmIDPre()
        a=Base.Libc.time();
        b=floor(UInt,100a);
        return b
    end

    function IsinSpaceRange(a::AlarmInfo,es1::UInt,es2::UInt,timetick::UInt )::Bool
        chn2Beg = 4096
        che2End = 8192
        try
            chn2Beg = Main.SegParaAlarmProc.getchn2BegID() #4096
            che2End = Main.SegParaAlarmProc.getSegNUmbers() #8192
        catch(e)
                println("-----:",e)
        end
        as1 = (a.S1   >= chn2Beg ) ? max( a.S1 - MergeWidth ,chn2Beg) : max( a.S1 - MergeWidth ,1)
        as2 = (a.S2 <  chn2Beg ) ? min( a.S2 + MergeWidth ,chn2Beg ) : min(  a.S2 + MergeWidth ,che2End )
        inrng = ( es2>= as1 && es2 <= as2 ) ||
                ( es1>= as1 && es1 <= as2 ) ||
                (es1 < as1 && es2 > as2 )

        return inrng
    end

    function IsinTimeRange(a::AlarmInfo,es1::UInt,es2::UInt,timetick::UInt )::Bool
        inrng =  ( timetick - a.tfresh ) <= MergeTime
        return inrng
    end

    function freshEvent2AlarmHIstory( a::AlarmInfo ,es1::UInt,es2::UInt,tick::UInt  )
        a.Tc  = getAlarmIDPre()
        a.S1 = es1
        a.S2 = es2
        a.SC = floor( UInt,(a.S2+a.S1)/2*10 )
        a.tfresh = tick
        AlarmReport.LogAlarm(a,false)
    end

    function addEvent2AlarmHistory( es1::UInt,es2::UInt,tick::UInt  )
        alrmobj = AlarmInfo()
        alrmobj.Tc  = getAlarmIDPre()
        # alrmobj.AlarmID = join( (Tc,string(es1),string(es2) ),"-" ) )
        alrmobj.S1 = es1
        alrmobj.S2 = es2
        alrmobj.SC = floor( UInt,(es1+es2)/2*10 )
        alrmobj.t0 = tick
        alrmobj.tfresh =tick
        mainID=@sprintf(":%013u-%04u-%04u",alrmobj.Tc ,es1, es2)
        # mainID = join( (string( alrmobj.Tc ),string(es1),string(es2) ) ,"-" )
        for i =1 : 10000
            alrmobj.AlarmID  = @sprintf("%s-%05u",mainID,i)
            # alrmobj.AlarmID = join( ( mainID, string(i) ),"-"  )
            if  get(AlarmInfoHistory,alrmobj.AlarmID,nothing ) == nothing
                AlarmInfoHistory[ alrmobj.AlarmID  ] = alrmobj
                break
            end
        end
        # writeLog
        AlarmReport.LogAlarm(alrmobj,true)
    end

    function mergeProc( es1::UInt,es2::UInt,tick::UInt )
        for data in values( AlarmInfoHistory )
            intimeRng = false
            inSpaceWidth = IsinSpaceRange(data,es1,es2,tick )
            if inSpaceWidth == true
                intimeRng = IsinTimeRange(data,es1,es2,tick )
            end
            if intimeRng == true
                 return freshEvent2AlarmHIstory(data,es1,es2,tick  )
            end
        end
        addEvent2AlarmHistory( es1,es2,tick  )
    end
    function initMergeCfg( ;MergeWidth1::UInt , MergeTime1::UInt )
        global MergeWidth , MergeTime
        MergeWidth =MergeWidth1;
        MergeTime = MergeTime1;
    end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# Dati di input

V = [0.6540618 0.2054992 0.2972308; 0.7142365 0.1455625 0.969203; 0.5941251 0.8769965 0.3624924; 0.6542998 0.8170598 1.0344647; 1.3260341 0.2707609 0.2428771; 1.3862088 0.2108241 0.9148494; 1.2660973 0.9422582 0.3081388; 1.326272 0.8823214 0.980111; -0.3874063 0.4902226 0.4536339; 0.3249123 0.707347 0.5231232; -0.1702819 -0.0864242 0.0297177; 0.5420367 0.1307001 0.099207; -0.317917 0.0663064 1.0658723; 0.3944016 0.2834308 1.1353616; -0.1007926 -0.5103404 0.6419561; 0.611526 -0.2932161 0.7114454; 0.7899026 0.0605793 0.6679889; 0.46601 0.0749997 0.6686316; 0.804323 0.3844725 0.6679746; 0.4804304 0.3988929 0.6686173; 0.7905452 0.060565 0.9922023; 0.4666527 0.0749854 0.992845; 0.8049656 0.3844582 0.992188; 0.4810731 0.3988786 0.9928307; -0.2261907 -0.0720455 0.4715635; -0.0499888 0.0863489 0.7965885; -0.0677963 0.219164 0.24378; 0.1084056 0.3775584 0.5688049; 0.0988343 -0.2998291 0.4063673; 0.2750362 -0.1414347 0.7313923; 0.2572286 -0.0086196 0.1785838; 0.4334306 0.1497748 0.5036087];
EV = SparseArrays.sparse([1, 5, 9, 1, 6, 10, 2, 5, 11, 2, 6, 12, 3, 7, 9, 3, 8, 10, 4, 7, 11, 4, 8, 12, 13, 17, 21, 13, 18, 22, 14, 17, 23, 14, 18, 24, 15, 19, 21, 15, 20, 22, 16, 19, 23, 16, 20, 24, 25, 29, 33, 25, 30, 34, 26, 29, 35, 26, 30, 36, 27, 31, 33, 27, 32, 34, 28, 31, 35, 28, 32, 36], [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24], Int8[-1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1]);
FE = SparseArrays.sparse([1, 3, 1, 4, 2, 3, 2, 4, 1, 5, 1, 6, 2, 5, 2, 6, 3, 5, 3, 6, 4, 5, 4, 6, 7, 9, 7, 10, 8, 9, 8, 10, 7, 11, 7, 12, 8, 11, 8, 12, 9, 11, 9, 12, 10, 11, 10, 12, 13, 15, 13, 16, 14, 15, 14, 16, 13, 17, 13, 18, 14, 17, 14, 18, 15, 17, 15, 18, 16, 17, 16, 18], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 35, 35, 36, 36], Int8[1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1]);

function arrange3Dfaces(V, copEV, copFE)
	EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments
	triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
	FVs = convert(Array{Lar.Cells}, triangulated_faces)
	V = convert(Lar.Points,V')
	return V,FVs,EVs
end

# Dati di input visualizzati prima della frammentazione

V,FVs,EVs = arrange3Dfaces(V, EV, FE);
GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));

# Dati visualizzati dopo la frammentazione

V = [0.6540617999999998 0.2054992 0.2972308000000001; 0.7142364999999998 0.1455625 0.969203; 0.6872816200343089 0.1724107693563844 0.6681972618528844; 0.5941250999999997 0.8769965 0.36249240000000005; 0.6542998079614951 0.8170598007732168 1.0344646993560203; 0.6929229746417847 0.3843473621810608 0.9924101013070119; 0.6660927570227951 0.3906268097799347 0.66824889023436; 0.6924256885846046 0.3894687221287142 0.9924113129852241; 1.3260341 0.2707609000000001 0.2428771; 1.3862088 0.21082410000000004 0.9148494000000003; 1.2660973 0.9422582000000003 0.3081388; 1.3262719920385087 0.8823213992267827 0.9801110006439793; 0.7946185131077457 0.1533691699791747 0.9627011677876783; 0.795311803386693 0.18207414025178234 0.6679835359595829; 0.8049461769791598 0.3844586322271989 0.9823884300888783; 0.6978893364451173 0.3892254686183435 0.9924004714668809; -0.38740630000000004 0.49022260000000006 0.45363390000000015; 0.3249123 0.707347 0.5231232000000001; -0.17028190000000007 -0.08642419999999995 0.029717700000000125; 0.542036694733776 0.13070013212634024 0.09920695360173326; -0.317917 0.06630639999999999 1.0658723; 0.3944015999999999 0.28343080000000004 1.1353616; -0.10079260000000007 -0.5103404 0.6419561000000001; 0.6115259947337761 -0.2932160678736597 0.7114453536017331; 0.4729938678063347 0.07470259548304928 0.9819171231357612; 0.4673971777131105 0.08956648867512645 0.992844167328483; 0.4666328572520818 0.07498584149882567 0.9828352213039846; 0.5039882186318663 0.07330882745780876 0.6685562398308738; 0.4707836778272794 0.18222015069807973 0.6686268661789595; 0.4673140236772549 0.08983923498806379 0.9928443441979011; 0.7899025999999999 0.060579300000000016 0.6679888999999999; 0.46601 0.0749997 0.6686315999999999; 0.8043229999999999 0.38447249999999994 0.6679745999999999; 0.4804303999999999 0.3988929 0.6686172999999999; 0.7905452000000002 0.060564999999999994 0.9922023000000001; 0.46665270000000003 0.0749854 0.9928450000000001; 0.8049656 0.3844582 0.9921880000000001; 0.4810731000000001 0.3988786 0.9928307000000001];
copEV = SparseArrays.sparse([1, 4, 14, 2, 5, 15, 1, 2, 7, 17, 3, 4, 19, 3, 6, 20, 5, 6, 9, 23, 7, 8, 46, 47, 8, 9, 54, 55, 10, 12, 14, 10, 13, 16, 11, 12, 19, 11, 13, 20, 15, 16, 18, 21, 17, 18, 48, 49, 21, 22, 62, 63, 22, 23, 53, 54, 24, 26, 36, 24, 27, 37, 25, 26, 38, 25, 27, 39, 28, 30, 36, 28, 31, 37, 29, 30, 38, 29, 33, 39, 32, 33, 34, 41, 31, 32, 35, 43, 65, 66, 67, 34, 35, 60, 61, 40, 41, 44, 45, 40, 42, 50, 51, 67, 42, 43, 57, 58, 44, 48, 59, 45, 50, 60, 46, 49, 62, 47, 51, 64, 52, 56, 59, 52, 57, 61, 65, 53, 56, 63, 55, 58, 64, 66], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 38], Int8[-1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 18, 20, 1, 2, 25, 26, 1, 2, 21, 22, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 18, 20, 4, 5, 28, 29, 6, 7, 6, 9, 8, 9, 28, 29, 8, 9, 25, 27, 8, 9, 21, 22, 10, 13, 10, 14, 10, 15, 10, 16, 11, 13, 11, 14, 11, 15, 11, 16, 12, 17, 11, 16, 11, 12, 23, 24, 11, 12, 30, 32, 13, 15, 13, 16, 14, 15, 14, 16, 16, 17, 18, 19, 16, 17, 23, 24, 16, 17, 16, 17, 18, 23, 19, 24, 20, 25, 18, 26, 18, 28, 20, 29, 19, 30, 18, 31, 21, 23, 21, 27, 22, 25, 21, 26, 21, 28, 21, 21, 23, 28, 24, 30, 23, 32, 25, 29, 27, 28, 26, 31, 32, 31, 30, 31], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 66, 67, 67], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1]);

V,FVs,EVs = arrange3Dfaces(V, copEV, copFE);
GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.,1.,1.,99,1));

copEF = convert(Lar.ChainOp, LinearAlgebra.transpose(copFE))
println(">>>>>>>> going to loop?")

# dopo stoppato il loop

ld_bounds=copFE   # cobordo_1

println(">>>>>>>>>>>>> eccomi !!!!")

lld_cellsnum, ld_cellsnum = size(ld_bounds) # !! i nomi lowerlower e lower sembrano invertirti ...
(32, 67)

count_marks = zeros(Int64, ld_cellsnum)
67-element Vector{Int64}: # il conto torna i marker 'uso sono associati agli spigoli !!

angles = Array{Array{Int64, 1}, 1}(undef, lld_cellsnum)
32-element Vector{Vector{Int64}}: #undef, #undef, #undef, ....
# array di array associati alle faccie, ma non mi convincono (dovrebbero essere associati agli spigoli)
# !! potrebbe diventare causa di errore quando gli spigoli diventano più numerosi delle facce .... (in esempi complessi)

julia> @show collect(enumerate(angles));
collect(enumerate(angles)) = [(1, [4, 1]), (2, [2, 5]), (3, [6, 1]), (4, [7, 1]), (5, [2, 8]), (6, [9, 1]), (7, [1, 20, 2, 18]), (8, [2, 22, 1, 21]), (9, [25, 2, 26, 1]), (10, [4, 3]), (11, [3, 6]), (12, [3, 7]), (13, [3, 9]), (14, [4, 7]), (15, [5, 8]), (16, [9, 4]), (17, [4, 18, 5, 20]), (18, [5, 28, 4, 29]), (19, [7, 6]), (20, [6, 9]), (21, [8, 21, 9, 22]), (22, [25, 8, 27, 9]), (23, [29, 9, 28, 8]), (24, [10, 13]), (25, [14, 10]), (26, [15, 10]), (27, [10, 16]), (28, [13, 11]), (29, [14, 11]), (30, [11, 15]), (31, [11, 16]), (32, [17, 12]), (33, [16, 11]), (34, [11, 32, 12, 30]), (35, [12, 23, 11, 24]), (36, [13, 15]), (37, [16, 13]), (38, [15, 14]), (39, [16, 14]), (40, [16, 17]), (41, [17, 24, 16, 23]), (42, [19, 17, 18, 16]), (43, [17, 16]), (44, [23, 18]), (45, [19, 24]), (46, [20, 25]), (47, [26, 18]), (48, [18, 28]), (49, [29, 20]), (50, [30, 19]), (51, [18, 31]), (52, [21, 23]), (53, [27, 21]), (54, [25, 22]), (55, [21, 26]), (56, [28, 21]), (57, [21]), (58, [21]), (59, [28, 23]), (60, [24, 30]), (61, [32, 23]), (62, [25, 29]), (63, [28, 27]), (64, [26, 31]), (65, [32]), (66, [31]), (67, [30, 31])]

@show collect(enumerate(Lar.cop2lar(sparse(ld_bounds'))))
collect(enumerate(Lar.cop2lar(ld_bounds))) = [(1, [1, 4]), (2, [2, 5]), (3, [1, 6]), (4, [1, 7]), (5, [2, 8]), (6, [1, 9]), (7, [1, 2, 18, 20]), (8, [1, 2, 21, 22]), (9, [1, 2, 25, 26]), (10, [3, 4]), (11, [3, 6]), (12, [3, 7]), (13, [3, 9]), (14, [4, 7]), (15, [5, 8]), (16, [4, 9]), (17, [4, 5, 18, 20]), (18, [4, 5, 28, 29]), (19, [6, 7]), (20, [6, 9]), (21, [8, 9, 21, 22]), (22, [8, 9, 25, 27]), (23, [8, 9, 28, 29]), (24, [10, 13]), (25, [10, 14]), (26, [10, 15]), (27, [10, 16]), (28, [11, 13]), (29, [11, 14]), (30, [11, 15]), (31, [11, 16]), (32, [12, 17]), (33, [11, 16]), (34, [11, 12, 30, 32]), (35, [11, 12, 23, 24]), (36, [13, 15]), (37, [13, 16]), (38, [14, 15]), (39, [14, 16]), (40, [16, 17]), (41, [16, 17, 23, 24]), (42, [16, 17, 18, 19]), (43, [16, 17]), (44, [18, 23]), (45, [19, 24]), (46, [20, 25]), (47, [18, 26]), (48, [18, 28]), (49, [20, 29]), (50, [19, 30]), (51, [18, 31]), (52, [21, 23]), (53, [21, 27]), (54, [22, 25]), (55, [21, 26]), (56, [21, 28]), (57, [21]), (58, [21]), (59, [23, 28]), (60, [24, 30]), (61, [23, 32]), (62, [25, 29]), (63, [27, 28]), (64, [26, 31]), (65, [32]), (66, [31]), (67, [30, 31])]

# ERRORI TOPOLOGICI NEI DATI D'INGRESSO: (edge, [incident-faces])
(57, [21])
(58, [21])
(65, [32])
(66, [31])

# visualizziamo il sotto-complesso delle facce [21, 31, 22]

julia> V[:,30], V[:,26]
([0.4673140236772549, 0.08983923498806379, 0.9928443441979011], [0.4673971777131105, 0.08956648867512645, 0.992844167328483])

julia> (V[:,30] - V[:,26]).^2
3-element Vector{Float64}:
 6.91459367907873e-9
 7.439055122090963e-8
 3.1282791063085847e-14

julia> sum((V[:,30] - V[:,26]).^2)
8.130517618277942e-8

julia> sqrt(sum((V[:,30] - V[:,26]).^2))
0.0002851406252759845


text:
@testset "file3 gets run" begin
    global file3_run
    @test !file3_run
    file3_run = true
    @test file3_run
end


text:
#!/usr/bin/env julia
# ARGS[1] points to the problem-specifications repo
# call with --project=generators

using JSON

p = joinpath(ARGS[1], "exercises", "clock", "canonical-data.json")
@debug "loading canonical data from $p"
const data = JSON.parsefile(p)

s  = "# canonical data version: $(data["version"]) (auto-generated)\n\n"
s *= "using Test\n\n"
s *= "include(\"$(data["exercise"]).jl\")\n\n"

parsetime(timestr) = Tuple(parse.(Int, split(timestr, ":")))

function parse_create(desc, input, expected)
    h, min = parsetime(expected)
    "    # $desc\n    @test Clock($(input["hour"]), $(input["minute"])) == Clock($h, $min)"
end

function parse_add(desc, input, expected; op = "+")
    h, min = parsetime(expected)
    "    # $desc\n    @test Clock($(input["hour"]), $(input["minute"])) $op Dates.Minute($(input["value"])) == Clock($h, $min)"
end

parse_sub(desc, input, expected) = parse_add(desc, input, expected; op = "-")

function parse_equal(desc, input, expected)
    c1 = (h=input["clock1"]["hour"], min=input["clock1"]["minute"])
    c2 = (h=input["clock2"]["hour"], min=input["clock2"]["minute"])

    neg = expected ? "=" : "!"

    "    # $desc\n    @test Clock($(c1.h), $(c1.min)) $neg= Clock($(c2.h), $(c2.min))"
end

function parse_show(desc, input, expected)
    "    @test sprint(show, Clock($(input["hour"]), $(input["minute"]))) == \"\\\"$expected\\\"\""
end

const mappings = Dict(
    "create" => parse_create,
    "add" => parse_add,
    "subtract" => parse_sub,
    "equal" => parse_equal,
    "show" => parse_show,
)

for case in data["cases"]
    global s
    s *= "@testset \"$(case["description"])\" begin\n"
    for c in case["cases"]
        s *= mappings[c["property"]](c["description"], c["input"], c["expected"]) * "\n\n"
    end
    s = s[1:end - 1] # Remove trailing newline
    s *= "end\n\n"
end

# showing clocks
s *= "@testset \"displaying clocks\" begin\n"
for c in data["cases"][1]["cases"]
    global s
    s *= mappings["show"](c["description"], c["input"], c["expected"]) * "\n"
end
s *= "end\n"

print(s)


text:
using PyPlot
using HDF5, JLD
using Statistics
using Base.Iterators

using SubmodularMaximization
using RosDataProcess

close("all")

experiment_name = "connectivity_test"
data_folder = "./data"

reprocess = false

num_sensors = 10
nominal_area = 2.0

num_hops = 1

num_agents = 2 .^ (3:8)
trials = 1:25

num_steps = 2 .^ (2:4)

sensor_radius(num_agents) = sqrt(nominal_area / (num_agents * pi))
station_radius(num_agents) = 2 * sensor_radius(num_agents)

# works out to the maximum radius for all interactions
communication_radius(num_agents) = 3*station_radius(num_agents)

#
# Define the solver factories
#

function make_hop_solver(p, num_partitions)
  num_agents = length(p.partition_matroid)
  radius = communication_radius(num_agents)

  partition_solver = generate_by_global_partition_size(num_agents,
                                                       num_partitions)

  MultiHopSolver(p, solver=partition_solver,
                 communication_range=radius,
                 num_hops=num_hops)
end

function make_sequential_solver(p)
  num_agents = length(p.partition_matroid)
  radius = communication_radius(num_agents)

  SequentialCommunicationSolver(p, radius)
end

# Converging auction solver
# (should converge in the same time as the sequential solver or faster)
function make_auction_solver(p, x...)
  num_agents = length(p.partition_matroid)
  radius = communication_radius(num_agents)

  AuctionSolver(p, radius, x...)
end

# Converging auction solver
# (should converge in the same time as the sequential solver or faster)
function make_local_auction_solver(p, x...)
  num_agents = length(p.partition_matroid)
  radius = communication_radius(num_agents)

  LocalAuctionSolver(p, radius, x...)
end

#
# Assemble solvers and test specifications
#

solvers_strings =
  [map(n -> (p->make_hop_solver(p, n), string("RSP-",n)), num_steps)...,
   (make_sequential_solver, "Sequential"),
   map(n -> (p->make_auction_solver(p, n), string("Auction-",n)), num_steps)...,
   (p->make_auction_solver(p), string("Auction")),
   map(n -> (p->make_local_auction_solver(p, n), string("Local-Auction-",n)),
       num_steps)...,
   (p->make_local_auction_solver(p), string("Local-Auction")),
  ]

solver_factories = map(first, solvers_strings)
solver_strings = map(last, solvers_strings)

solver_inds = 1:length(solver_factories)
all_tests = product(solver_inds, num_agents, trials)

#
# Construct and run trials
#

function trial_fun(x)
  (solver_ind, num_agents, trial) = x
  make_solver = solver_factories[solver_ind]

  agent_specification = CircleAgentSpecification(sensor_radius(num_agents),
                                                 station_radius(num_agents),
                                                 num_sensors)

  problem = generate_connected_problem(agent_specification, num_agents,
                                       objective =
                                        x -> mean_area_coverage(x, 100),
                                       communication_radius =
                                        communication_radius(num_agents))

  solver = make_solver(problem)

  solution = solve_problem(solver, problem, threaded=false)

  (solver=solver, solution=solution)
end
function print_summary(x)
  (solver_ind, num_agents, trial) = x

  println("Num. Agents: ", num_agents,
          " Solver: ", solver_strings[solver_ind],
          " Trial: ", trial)
end

println("Loading data / Running experiments")

@time results = run_experiments(all_tests,
                                trial_fun=trial_fun,
                                print_summary=print_summary,
                                experiment_name=experiment_name,
                                data_folder=data_folder,
                                reprocess=reprocess
                               )

println("Plotting")

solver_groups = 3
solver_configurations = length(num_steps) + 1
# One color for each solver group
colors = repeat(RosDataProcess.generate_colors(solver_groups), inner=num_steps)
# One style for each solver configuration
style_atoms = [(0, (1, 1)), # densely dotted
               (0, (2, 1, 1, 0)), # densely dashed
               (0, (2, 1, 1, 1, 1, 0)), # dashdotted
               "-"
              ]
styles = repeat(style_atoms, outer=solver_groups)

series_by_solver = Any[]
for solver_ind in solver_inds
  solver_data = map(product(num_agents, trials)) do (num_agents, trial)
    results[solver_ind, num_agents, trial]
  end
  push!(series_by_solver, TimeSeries(num_agents, solver_data))
end

plot_specs = [
              (get = x -> x.solution.value,
               ylabel = "Objective (fraction of unit area)",
               name = "objective"),
              (get = x -> communication_span(x.solver),
               ylabel = "Messaging span",
               name = "span",
               after = [loglog]
              ),
              (get = x -> communication_messages(x.solver),
               ylabel = "Messages sent",
               name = "messages",
               after = [loglog]
              ),
              (get = x -> communication_volume(x.solver),
               ylabel = "Messaging volume",
               name = "volume",
               after = [loglog]
              ),
             ]

fig_dir = "fig/connectivity_test"
if !isdir(fig_dir)
  mkpath(fig_dir)
end

for spec in plot_specs
  println("  Plot: " * spec.name)

  figure()
  for solver_ind in reverse(solver_inds)
    data = map(spec.get, series_by_solver[solver_ind])
    plot_trials(data,
                label=solver_strings[solver_ind],
                linestyle=styles[solver_ind],
                color=colors[solver_ind],
                mean=true,
                standard_error=true)
  end
  legend(ncol=3)
  xlabel("Num. Agents")
  ylabel(spec.ylabel)

  if in(:after, keys(spec))
    foreach(x->x(), spec.after)
  end

  println("  Saving")
  save_latex(fig_dir, spec.name)
  println("  Saved")
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
## Vector of fun routines

function coefficientmatrix{N,F}(::Type{N},f::Vector{F},o...)
    if isempty(f)
        return Array(N,0,0)
    end

    n=mapreduce(ncoefficients,max,f)
    m=length(f)
    R=zeros(N,n,m)
    for k=1:m
        R[1:ncoefficients(f[k]),k]=coefficients(f[k],o...)
    end
    R
end

scalarorfuntype{S,T<:Number}(::Fun{S,T})=T
scalarorfuntype{T<:Number}(::T)=T
scalarorfuntype{T<:Number}(b::Vector{T})=T
scalarorfuntype(b::Vector{Any})=promote_type(map(scalarorfuntype,b)...)
scalarorfuntype{F<:Fun}(b::Vector{F})=promote_type(map(scalarorfuntype,b)...)

coefficientmatrix{F<:Fun}(Q::Vector{F},o...)=coefficientmatrix(scalarorfuntype(Q),Q,o...)
coefficientmatrix(Q::Vector{Any})=(@assert isempty(Q); zeros(0,0))

function values{D,N}(f::Vector{Fun{D,N}})
    n=mapreduce(ncoefficients,max,f)
    m=length(f)
    R=zeros(N,n,m)
    for k=1:m
        R[:,k] = values(pad(f[k],n))
    end
    R
end

function values{D,T}(p::Array{Fun{D,T},2})
    @assert size(p)[1] == 1

   values(vec(p))
end

## evaluation

#TODO: fix for complex
evaluate{T<:Fun}(A::AbstractArray{T},x::Number)=typeof(first(A)(x))[Akj(x) for Akj in A]

function evaluate{T<:Fun}(A::AbstractVector{T},x::AbstractVector)
    n=length(x)
    ret=Array(promote_type(eltype(x),mapreduce(eltype,promote_type,A)),length(A),n)

    for k=1:length(A)
        bkr=evaluate(A[k],x)

        for j=1:n
            ret[k,j]=bkr[j]
        end
    end

    ret
end

## Algebra

## scalar fun times vector

 for op in (:*,:(Base.Ac_mul_B),:(Base.At_mul_B))
     @eval begin
         function ($op){T<:Number,V<:Number,D}(A::Array{T,2}, p::Vector{Fun{D,V}})
             cfs=$op(A,coefficientmatrix(p).')
             ret = Array(Fun{D,promote_type(T,V)},size(cfs,1))
             for i = 1:size(cfs,1)
                 ret[i] = chop!(Fun(vec(cfs[i,:]),first(p).space),eps())
             end
             ret
         end

         function ($op){T<:Number,D}(p::Vector{Fun{D,T}},A::Array{T,2})
             cfs=$op(A,coefficientmatrix(p).')
             ret = Array(Fun{D,T},size(cfs,1))
             for i = 1:size(cfs,1)
                 ret[i] = chop!(Fun(vec(cfs[i,:]),first(p).space),eps())
             end
             ret
         end
     end
 end

#Allow vecfun + constvec, etc
#can't just promote constant vector to a vector-valued fun because don't know the domain.
for op = (:+,:-,:.*,:./)
    @eval begin
        ($op){T<:Number,S,V}(f::Fun{S,V},c::AbstractArray{T})=devec($op(vec(f),c))
        ($op){T<:Number,S,V}(c::AbstractArray{T},f::Fun{S,V})=devec($op(c,vec(f)))
    end
end


text:
observe(a::Agent, reward::Float64) = 0
function choose(a::Agent)
  action_idx = choose(a.policy, a)
  a.last_action = action_idx
  return action_idx
end

mutable struct BasicAgent <: Agent
  empirical_mean::AbstractVector{Float64}
  policy::Policy
  bandit::Bandit
  period::Integer
  action_attempts::AbstractVector{Float64}
  total_reward::AbstractVector{Float64}
  last_action::Integer
end
function BasicAgent(policy::Policy, bandit::Bandit, initial_mean::Integer, belief::AbstractVector{Float64})
  BasicAgent(belief, policy, bandit, 0, zeros(bandit.num_arms), zeros(bandit.num_arms), -1)
end
BasicAgent(policy::Policy, bandit::Bandit, belief::AbstractVector{Float64}) = BasicAgent(policy, bandit, 0, belief)

function observe(agent::BasicAgent, reward::Float64)
  agent.period += 1
  agent.action_attempts[agent.last_action] += 1
  agent.total_reward[agent.last_action] += reward
  agent.empirical_mean[agent.last_action] = agent.total_reward[agent.last_action] / agent.action_attempts[agent.last_action]
end

mutable struct BetaBernoulliAgent <: Agent
  α::AbstractVector{Float64}
  β::AbstractVector{Float64}
  empirical_mean::AbstractVector{Float64}
  policy::Policy
  bandit::Bandit
  period::Integer
  action_attempts::AbstractVector{Float64}
  total_reward::AbstractVector{Float64}
  last_action::Integer
end
function BetaBernoulliAgent(α::AbstractVector{Float64}, β::AbstractVector{Float64}, policy::Policy, bandit::Bandit)
  BetaBernoulliAgent(α, β, α ./ (α + β), policy, bandit, 0, zeros(bandit.num_arms), zeros(bandit.num_arms), -1)
end
BetaBernoulliAgent(means::AbstractVector{Float64}, policy::Policy, bandit::Bandit) = BetaBernoulliAgent(means, 1 .- means, policy, bandit)
BetaBernoulliAgent(policy::Policy, bandit::Bandit) = BetaBernoulliAgent(ones(bandit.num_arms), ones(bandit.num_arms), policy, bandit)

observe(agent::BetaBernoulliAgent, reward::Integer) = observe(agent, convert(AbstractFloat, reward))
function observe(agent::BetaBernoulliAgent, reward::Float64)
  agent.period += 1
  agent.action_attempts[agent.last_action] += 1
  agent.total_reward[agent.last_action] += reward
  if reward == 1
    agent.α[agent.last_action] += 1
  else
    agent.β[agent.last_action] += 1
  end
  α = agent.α[agent.last_action]
  β = agent.β[agent.last_action]
  agent.empirical_mean[agent.last_action] = α / (α + β)
end

randomSample(agent::BetaBernoulliAgent) = [rand(Beta(agent.α[i], agent.β[i])) for i in 1:length(agent.empirical_mean)]

mutable struct NormalAgent <: Agent
  prior_mean::AbstractVector{Float64}
  prior_variance::AbstractVector{Float64}
  empirical_mean::AbstractVector{Float64}
  policy::Policy
  bandit::Bandit
  period::Integer
  action_attempts::AbstractVector{Float64}
  total_reward::AbstractVector{Float64}
  last_action::Integer
end
function NormalAgent(prior_mean::AbstractVector{Float64}, prior_variance::AbstractVector{Float64}, policy::Policy, bandit::Bandit)
  NormalAgent(prior_mean, prior_variance, zeros(length(prior_mean)), policy, bandit, -1, zeros(length(prior_mean)), zeros(length(prior_mean)), -1)
end
NormalAgent(policy::Policy, bandit::Bandit) = NormalAgent(zeros(bandit.num_arms), ones(bandit.num_arms), policy, bandit)

get_variances(a::NormalAgent) = a.prior_variance ./ (1 .+ a.action_attempts)

function observe(agent::NormalAgent, reward::Float64)
  agent.period += 1
  agent.action_attempts[agent.last_action] += 1
  agent.total_reward[agent.last_action] += reward
  agent.empirical_mean[agent.last_action] = agent.total_reward[agent.last_action] / agent.action_attempts[agent.last_action]
end

function randomSample(agent::NormalAgent)
  variances = get_variances(agent)
  rands = [rand(Normal(agent.empirical_mean[i], variances[i])) for i in 1:length(agent.empirical_mean)]
  return rands
end


text:
module AC274

export neighbor, has_neighbor, DG1D, DG2D, CG2D, solve, generateMesh, Mesh1D,
    plotSolution, invchi, chi, 𝜒⁻¹, 𝜒, plotMesh, evaluate, evaluate_ref, interpolate,
    chi, invchi, vertices, mesh, nbf, 𝖓, cid, do_quad_ref, ∂, rewriteMesh, calciᴰ, ℳ,
    𝒩, createDualMesh, apply∂D!, applyA∂D!, drawSupport, generateMatrices, stiffness,
    loadV

include("simpletypes.jl")
include("mesh.jl")
include("simulation.jl")
include("quadrature.jl")
include("neighbors.jl")
include("nodalbasis.jl")
include("RHS.jl")
include("solve.jl")
include("plot1d.jl")
include("plot2d.jl")
include("cg.jl")

end


text:
abstract type FiniteMDP end

struct DeterministicFiniteMDP <: FiniteMDP
    next_state::Function
    reward::Function
    actions_no::Int
    states_no::Int
    terminal_state::Int
end # struct

struct StochasticFiniteMDP{T<:Real} <: FiniteMDP
    probabilities::Function
    states_no::Int
    actions_no::Int
    rewards::Array{T}
    terminal_state::Int
end # struct

states_no(fmdp::DeterministicFiniteMDP) = fmdp.states_no
states_no(fmdp::StochasticFiniteMDP) = fmdp.states_no

actions_no(fmdp::DeterministicFiniteMDP) = fmdp.actions_no
actions_no(fmdp::StochasticFiniteMDP) = fmdp.actions_no

terminal_state(fmdp::DeterministicFiniteMDP) = fmdp.terminal_state
terminal_state(fmdp::StochasticFiniteMDP) = fmdp.terminal_state

act_once(fmdp::DeterministicFiniteMDP, s::Integer, a::Integer) = (fmdp.next_state(s, a), fmdp.reward(s, a))

random_action(fmdp::FiniteMDP) = rand(1:actions_no(fmdp))


text:
function fib(n)
    n < 3 ? 1 : fib(n-2) + fib(n-1)
end

"""
   This is code text.
"""
function memoize(f)
    memo = Dict()
    (args...; kwargs...) -> begin
        x = (args, kwargs)
        haskey(memo, x) || (memo[x] = f(args..., kwargs...))
        memo[x]
    end
end

fib! = memoize(fib)

fib!(20)

max_stars_repo_path: Chapter 12/fibonacci.jl

max_issues_repo_path: Chapter 12/fibonacci.jl

max_forks_repo_path: Chapter 12/fibonacci.jl


text:
using Memento
logger = get_logger("bot-logger")

include("utils.jl")

struct GameMap
    id::String
    width::Int
    height::Int
    players::Dict{String, Player}
    planets::Dict{String, Planet}

    function GameMap(id, width, height, map_string::String)
        game_map = new(id, width, height, Dict{String, Player}(), Dict{String, Planet}())
        map_parse!(game_map, map_string)

        game_map
    end
end

"""
    map_parse!(game_map, map_string)

Parse the map description from the game. Changes definition of game_map.players and game_map.planets, so it should be used only on newly created GameMap. `map_string` contains the Halite engine outputs.
"""
function map_parse!(game_map::GameMap, map_string::String)
    tokens = Vector{String}(split(map_string))
    num_players = parse(Int, shift!(tokens))
    debug(logger, @sprintf("Number of players: %s", num_players))
    for i in 1:num_players
        player_id = shift!(tokens)
        player = Player(player_id)
        num_ships = parse(Int, shift!(tokens))
        for j in 1:num_ships
            ship = Ship(player.id, tokens)
            player.ships[ship.id] = ship
            debug(logger, @sprintf("Owner: %s, Ship: %s, x: %d, y: %d", ship.owner_id, ship.id, ship.x, ship.y))
        end
        game_map.players[player_id] = player
    end
    num_planets = parse(Int, shift!(tokens))
    for i in 1:num_planets
        planet = Planet(tokens)
        game_map.planets[planet.id] = planet
    end
    @assert isempty(tokens)
end

# Utility functions
"""
    get_player(game_map, player_id)
    get_player(game_map, entity)

Returns player defined by `player_id` or `entity.owner_id`.
"""
get_player(game_map::GameMap, player_id::String) = game_map.players[player_id]
get_player(game_map::GameMap, entity::Entity) = get_player(game_map, entity.owner_id)

"""
    get_me(game_map)

Returns the user's player.
"""
get_me(game_map::GameMap) = get_player(game_map, game_map.id)

"""
    all_players(game_map)

List of all players.
"""
all_players(game_map::GameMap) = values(game_map.players)

"""
    all_planets(game_map)

List of all planets.
"""
all_planets(game_map::GameMap) = values(game_map.planets)

"""
    get_planet(game_map, planet_id)
    get_planet(game_map, ship)

Returns planet defined by `planet_id` or `ship.docked_planet_id`.
"""
get_planet(game_map::GameMap, planet_id::String) = game_map.planets[planet_id]
get_planet(game_map::GameMap, ship::Ship) = get_planet(game_map, ship.docked_planet_id)

"""
    all_docked_ships(game_map, planet)

List of all docked ships, empty if planet is not owned.
"""
function all_docked_ships(game_map::GameMap, planet::Planet)
    if planet.owned
        return map(x -> get_ship(get_player(game_map, planet), x), planet.docked_ships_ids)
    else
        return Vector{Ship}()
    end
end

"""
    get_docked_ship(game_map, planet, ship_id)

Return docked ship. Throw error if ship is not docked on this planet.
"""
get_docked_ship(game_map::GameMap, planet::Planet, ship_id::String) = get_ship(get_player(game_map, planet), ship_id)

"""
    all_ships(game_map)

List of all ships in current game map.
"""
function all_ships(game_map::GameMap)
    all_ships = Vector{Ship}()
    for player in all_players(game_map)
        all_ships = vcat(all_ships, all_ships(player, true))
    end

    return all_ships
end

"""
    obstacles_between(game_map, entity1, entity2, ignore_ships, ignore_planets)

Check whether there is a straight-line path between two entities, without planetary or ships obstacles in between. Returns list of obstacles, empty if there is no obstacles were found.
"""
function obstacles_between(game_map::GameMap, e1::Entity, e2::Entity, ignore_ships::Bool, ignore_planets::Bool)
    obstacles = Vector{Entity}()
    if !ignore_ships
        for player in all_players(game_map)
            for entity in all_ships(player)
                if e1 == entity || e2 == entity
                    continue
                end
                if intersect_segment_circle(e1, e2, entity, radius(e1) + 0.1)
                    push!(obstacles, entity)
                end
            end
        end
    end

    if !ignore_planets
        for entity in all_planets(game_map)
            if e1 == entity || e2 == entity
                continue
            end
            if intersect_segment_circle(e1, e2, entity, radius(e1) + 0.1)
                push!(obstacles, entity)
            end
        end
    end

    return obstacles
end

"""
    navigate(game_map, ship, target)

Move a `ship` to a specific `target` position (Entity). It is recommended to place the position itself here, else navigate will crash into the target. If `avoid_obstacles` is set to true (default) will avoid obstacles on the way, with up to `max_corrections` corrections. Note that each correction accounts for `angular_step` degrees difference, meaning that the algorithm will naively try `max_correction` degrees before giving up (and returning empty string). The navigation will only consist of up to one command; call this method again in the next turn to continue navigating to the position.

# Arguments
- `game_map::GameMap`: The map of the game, from which obstacles will be extracted
- `ship::Ship`: ship to navigate
- `target::Entity`: The entity to which you will navigate
- `speed::Int`: The (max) speed to navigate. If the obstacle is nearer, will adjust accordingly.
- `avoid_obstacles::Bool`: Whether to avoid the obstacles in the way (simple pathfinding).
- `max_corrections::Int`: The maximum number of degrees to deviate per turn while trying to pathfind. If exceeded returns empty string.
- `angular_step::Float64`: The degree difference to deviate if the original destination has obstacles
- `ignore_ships::Bool`: Whether to ignore ships in calculations (this will make your movement faster, but more precarious)
- `ignore_planets::Bool`: Whether to ignore planets in calculations (useful if you want to crash onto planets)
"""
function navigate(game_map::GameMap, ship::Ship, target::Entity;
                 speed::Int = 0,
                 avoid_obstacles::Bool = true,
                 max_corrections::Int = 90,
                 angular_step::Float64 = 1.0,
                 ignore_ships::Bool = false,
                 ignore_planets::Bool = false)
    if max_corrections <= 0
        return ""
    end
    distance = calculate_distance_between(ship, target)
    angle = calculate_angle_between(ship, target)
    if avoid_obstacles && !isempty(obstacles_between(game_map, ship, target, ignore_ships, ignore_planets))
        new_target_dx = cos(deg2rad(angle + angular_step)) * distance
        new_target_dy = sin(deg2rad(angle + angular_step)) * distance
        new_target = Position(ship.x + new_target_dx, ship.y + new_target_dy)

        return navigate(game_map, ship, new_target,
                        speed = speed,
                        avoid_obstacles = avoid_obstacles,
                        max_corrections = max_corrections - 1,
                        angular_step = angular_step,
                        ignore_ships = ignore_ships,
                        ignore_planets = ignore_planets)
    end

    speed = distance >= speed ? speed : floor(Int, distance)
    return thrust(ship, speed, angle)
end

max_stars_repo_name: xinnosuke/Halite-II

max_issues_repo_name: xinnosuke/Halite-II

max_forks_repo_name: xinnosuke/Halite-II


text:
#################### Heidelberger and Welch Diagnostic ####################

function heideldiag(x::Vector{<:Real}; alpha::Real=0.05, eps::Real=0.1,
                             etype=:imse, start::Integer=1, args...)
  n = length(x)
  delta = trunc(Int, 0.10 * n)
  y = x[trunc(Int, n / 2):end]
  S0 = length(y) * mcse(y, etype; args...)^2
  i, pvalue, converged, ybar = 1, 1.0, false, NaN
  while i < n / 2
    y = x[i:end]
    m = length(y)
    ybar = mean(y)
    B = cumsum(y) - ybar * collect(1:m)
    Bsq = (B .* B) ./ (m * S0)
    I = sum(Bsq) / m
    pvalue = 1.0 - pcramer(I)
    converged = pvalue > alpha
    if converged
      break
    end
    i += delta
  end
  halfwidth = sqrt(2.0) * erfinv(1.0 - alpha) * mcse(y, etype; args...)
  passed = halfwidth / abs(ybar) <= eps
  [i + start - 2, converged, pvalue, ybar, halfwidth, passed]
end

function heideldiag(
    chains::Chains;
    sections = _default_sections(chains),
    alpha = 0.05,
    eps = 0.1,
    etype = :imse,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Preallocate.
    _, p, m = size(_chains.value)
    vals = [Array{Float64}(undef, p, 6) for i in 1:m]

    # Perform tests.
    for j in 1:p, k in 1:m
        vals[k][j, :] = heideldiag(
                            collect(skipmissing(_chains.value[:, j, k])),
                            alpha=alpha,
                            eps=eps,
                            etype=etype,
                            start=first(_chains);
                            kwargs...
                           )
    end

    # Retrieve columns.
    data = [[vals[k][:, i] for i in 1:6] for k in 1:m]

    # Obtain names of parameters.
    names_of_params = names(_chains)

    # Compute data frames.
    vector_of_df = [
        ChainDataFrame(
            "Heidelberger and Welch Diagnostic - Chain $i",
            (parameters = names_of_params, burnin = columns[1], stationarity = columns[2],
             pvalue = columns[3], mean = columns[4], halfwidth = columns[5],
             test = columns[6])
        )
        for (i, columns) in enumerate(data)
    ]

    return vector_of_df
end


text:
module constants

const g = 9.81             #: Gravitational constant, m/s^2
const Cp = 1004.0          #: Specific heat of dry air at constant pressure, J/kg
const L = 2.25e6           #: Latent heat of condensation, J/kg
const rho_w = 1e3          #: Density of water, kg/m^3
const R = 8.314            #: Universal gas constant, J/(mol K)
const Mw = 18.0/1e3        #: Molecular weight of water, kg/mol
const Ma = 28.9/1e3        #: Molecular weight of dry air, kg/mol
const Rd = R/Ma            #: Gas constant for dry air, J/(kg K)
const Rv = R/Mw            #: Gas constant for water vapor, J/(kg K)
const Dv = 3.e-5           #: Diffusivity of water vapor in air, m^2/s
const ac = 1.0             #: condensation constant
const Ka = 2.e-2           #: Thermal conductivity of air, J/m/s/K
const at = 0.96            #: thermal accomodation coefficient
const epsilon = 0.622      #: molecular weight of water / molecular weight of dry air

# Additional fixed model parameters
const N_STATE_VARS = 7
const STATE_VARS = ["z", "P", "T", "wv", "wc", "wi", "S"]
const STATE_VAR_MAP = {var: i for i, var in enumerate(STATE_VARS)}

end # module constants


text:
# ---
# title: 125. Valid Palindrome
# id: problem125
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Two Pointers, String
# link: <https://leetcode.com/problems/valid-palindrome/description/>
# hidden: true
# ---
#
# Given a string, determine if it is a palindrome, considering only alphanumeric
# characters and ignoring cases.
#
# **Note:**  For the purpose of this problem, we define empty string as valid
# palindrome.
#
# **Example 1:**
#
#
#
#     Input: "A man, a plan, a canal: Panama"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: "race a car"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `s` consists only of printable ASCII characters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Maximum flow" begin
    #### Graphs for testing
    graphs = [
    # Graph with 8 vertices
    (8,
    [
    (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
    (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
    (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
    ],
    1, 8, # source/target
    3,    # answer for default capacity
    28,   # answer for custom capacity
    15, 5 # answer for restricted capacity/restriction
    ),

    # Graph with 6 vertices
    (6,
    [
    (1, 2, 9), (1, 3, 9), (2, 3, 10), (2, 4, 8), (3, 4, 1),
    (3, 5, 3), (5, 4, 8), (4, 6, 10), (5, 6, 7)
    ],
    1, 6, # source/target
    2,   # answer for default capacity
    12,  # answer for custom capacity
    8, 5  # answer for restricted capacity/restriction
    )
    ]

    for (nvertices, flow_edges, s, t, fdefault, fcustom, frestrict, caprestrict) in graphs
        flow_graph = Graphs.DiGraph(nvertices)
        for g in testdigraphs(flow_graph)
            capacity_matrix = zeros(Int, nvertices, nvertices)
            for e in flow_edges
                u, v, f = e
                Graphs.add_edge!(g, u, v)
                capacity_matrix[u, v] = f
            end

            # Test DefaultCapacity
            d = @inferred(GraphsFlows.DefaultCapacity(g))
            T = eltype(d)
            @test typeof(d) <: AbstractMatrix{T}
            @test d[s, t] == 0
            @test size(d) == (nvertices, nvertices)
            @test typeof(transpose(d)) <: GraphsFlows.DefaultCapacity
            @test typeof(adjoint(d)) <: GraphsFlows.DefaultCapacity

            # Test all algorithms - type instability in PushRelabel #553
            for ALGO in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
                @test maximum_flow(g, s, t, algorithm=ALGO())[1] == fdefault
                @test maximum_flow(g, s, t, capacity_matrix, algorithm=ALGO())[1] == fcustom
                @test maximum_flow(g, s, t, capacity_matrix, algorithm=ALGO(), restriction=caprestrict)[1] == frestrict
            end
        end
    end
end


text:
function get_dft_interval(max_val::Float64,N::Int64)
    #Returns the interval over which we obtain the Fourier transform
	(-N÷2:N-1-N÷2)*(2*max_val/N)
end

function get_reciprocal_dft_interval(max_val::Float64,N::Int64)
    #Returns the interval over which we have to sample our function to perform the DFT
	(-N÷2:N-1-N÷2)*(π/max_val)
end

function get_fourier_transform(k_max::Float64,N::Int64,ψ0)
	#Calculates the Fourier transform of the given initial profile ψ0

    ks = get_dft_interval(k_max,N)
    rs = get_reciprocal_dft_interval(k_max,N)

	#This is a matrix containing the ψ0s
	ψ0s = convert(Matrix{ComplexF64},[ ψ0(x,y) for x in rs, y in rs])

	#Then, we calculate the Fourier transform of ψ, which returns a matrix
	ϕs =  (0.5*π/k_max^2)*ifftshift( fft( fftshift(  ifftshift(fft(fftshift(ψ0s,1),1),1)  ,2),2),2)

	#Now, we interpolate this matrix, to be able to calculate the function at an arbitrary point.
	CubicSplineInterpolation( (ks , ks),ϕs,extrapolation_bc=0.0)
end

function display_fourier_transform(plot_range::Float64,plot_N::Int64,clim::Float64,ϕ)
	#Displays the fourier transform ϕ. It is useful to see if the chosen parameters are adequate
	ks = LinRange(-plot_range,plot_range,plot_N)

	I = [ real(ϕ(kx,ky)*conj(ϕ(kx,ky))) for ky in ks, kx in ks ]

	heatmap(ks,ks,I,clims=(0,clim),aspect_ratio=:equal,
	xlabel=L"k_x",ylabel=L"k_y",xlims=(-plot_range,plot_range),ylims=(-plot_range,plot_range),
	title="Fourier transform of initial profile")
end

function propagate_beam(r_max::Float64,N::Int64,z::Float64,ϕ)
	#Propagates a beam to a distance z given the fourier transform of its initial profile

	rs = get_dft_interval(r_max,N)
    ks = get_reciprocal_dft_interval(r_max,N)

	s = zeros(ComplexF64,N,N)

	Threads.@threads for n in 1:N
		Threads.@threads for m in 1:N
			s[n,m] = ϕ(ks[m],ks[n])*exp(-im*z*(ks[n]^2+ks[m]^2)/2)
		end
	end

	( 0.5*π/r_max^2 )*ifftshift( bfft( fftshift(  ifftshift(bfft(fftshift(s,1),1),1)  ,2),2),2)
end

function propagate_beam(r_max::Float64,N::Int64,z::Float64,ϕ,bplan1,bplan2)
	#Same as before, but accepts plans for the FFT.

	rs = get_dft_interval(r_max,N)
    ks = get_reciprocal_dft_interval(r_max,N)

	s = zeros(ComplexF64,N,N)

	Threads.@threads for n in 1:N
		Threads.@threads for m in 1:N
			s[n,m] = ϕ(ks[m],ks[n])*exp(-im*z*(ks[n]^2+ks[m]^2)/2)
		end
	end

	( 0.5*π/r_max^2 )*ifftshift(bplan2*fftshift(ifftshift(bplan1*fftshift(s,1),1),2),2)
end

function decrease_size(a,rs,I)
	#This function decreses the size of a range and corresponding intensities.
	#It is analogous to setting xlims=(-a,a), ylims=(-a,a) but somehow is faster
	K=1
	for k in eachindex(rs)
		if rs[k]>-a
			K = k
			break
		end
	end
	rs[ K:(length(rs)+1-K) ],I[K:(length(rs)+1-K),K:(length(rs)+1-K)]
end;

function display_beam(r_max::Float64,N::Int64,z::Float64,plot_range::Float64,clim::Float64,ϕ)
	#Displays a propageted beam
	rs = get_dft_interval(r_max,N)
	ψs = propagate_beam(r_max,N,z,ϕ)

	I = real(ψs.*conj(ψs))

	drs,dI = decrease_size(plot_range,rs,I)

	heatmap(drs,drs,dI,clims=(0,clim),aspect_ratio=:equal,title="z="*(@sprintf "%.2f" z),xlabel=L"x",ylabel=L"y",xlims=(-plot_range,plot_range),ylims=(-plot_range,plot_range))
end

function animate_beam(r_max::Float64,N::Int64,zmax::Float64,
	plot_range::Float64,clim::Float64,ϕ,fps::Int64,nframes::Int64,file_name::String)
	#Animates a beam propagation
	#The file_name should end in .gif

	rs = get_dft_interval(r_max,N)
	ks = get_reciprocal_dft_interval(r_max,N)

	plan1 = plan_bfft( zeros(ComplexF64, N,N ) , 1)
	plan2 = plan_bfft( zeros(ComplexF64, N,N ) , 2)

	anim = @animate for z in LinRange(0,zmax,nframes)
		ψs = propagate_beam(r_max,N,z,ϕ,plan1,plan2)

		I = real(ψs.*conj(ψs))

		drs,dI = decrease_size(plot_range,rs,I)

		heatmap(drs,drs,dI,clims=(0,clim),aspect_ratio=:equal,
		title=L"z="*(@sprintf "%.2f" z),xlabel=L"x",ylabel=L"y",
		xlims=(-plot_range,plot_range),size=(1200,900),labelfontsize=20,titlefontsize=20,tickfontsize=14)
	end

	gif(anim,file_name,fps=fps)
end;


text:
using HSA
include("helpers.jl")
using FactCheck

facts("Code Objects and Executables") do
    @use rt = HSA.Runtime() begin

        @with_agents context("The ISAs") do
            context("can be queried") do
                a = get_testagent()

                i_s_a = HSA.agent_info_isa(a)

                @fact isa(i_s_a, HSA.hsa_isa_t) --> true
            end
        end

        context("Executables") do
            context("can be created") do
                @use ex = Executable() begin

                    @fact isa(ex, Executable) --> true
                    @fact HSA.executable_info_state(ex) --> HSA.ExecutableStateUnfrozen
                    @fact HSA.executable_info_profile(ex) --> HSA.ProfileFull
                end
            end

            context("can be created with a profile") do
                @use ex = Executable(profile = HSA.ProfileBase) begin

                    @fact HSA.executable_info_profile(ex) --> HSA.ProfileBase
                end
            end

            context("can be frozen") do
                @use ex = Executable() begin

                    HSA.executable_freeze(ex)

                    @fact HSA.executable_info_state(ex) --> HSA.ExecutableStateFrozen
                end
            end
        end

    end
end


text:
module Record
using MLStyle
using MLStyle.Toolz: isCapitalized, ($), cons, nil
using MLStyle.MatchCore
using MLStyle.Qualification
using MLStyle.Infras
using MLStyle.Pervasives
using MLStyle.Render: render

export @as_record
export def_record

function def_record(ctor, record_fields, qualifier :: Qualifier, mod)
    ctor_name = string(ctor)
    n_destructor_args = length(record_fields)
    mk_match(tag, hd_obj, destruct_fields, mod) = begin
        check_if_given_field_names = map(destruct_fields) do field
            @match field begin
            Expr(:kw, _...) => true
            _               => false
            end
        end
        TARGET = mangle(mod)

        if all(check_if_given_field_names) # begin if
            map(destruct_fields) do field_
            @match field_ begin
                Expr(:kw, field::Symbol, pat) => begin
                    let ident = mangle(mod), field = field
                        function(body)
                            @format [TARGET, body, ident] quote
                                ident = TARGET.$field
                                body
                            end
                        end ∘ mk_pattern(ident, pat, mod)
                    end
                end
                _ => @syntax_err "The field name of destructor must be a Symbol!"
            end
            end
        elseif all(map(!, check_if_given_field_names))
            n_d = length(destruct_fields)
            if n_d == 1 && destruct_fields[1] == :(_)
                []
                # ignore fields
            else
                @assert n_d == n_destructor_args "Malformed destructing for case class $ctor_name(from module $(nameof(mod)))."
                map(zip(destruct_fields, record_fields)) do (pat, field)
                        let ident = mangle(mod)
                            function (body)
                                @format [TARGET, body, ident] quote
                                    ident = TARGET.$field
                                    body
                                end
                            end ∘ mk_pattern(ident, pat, mod)
                        end
                end
            end
        else
            @syntax_err "Destructor should be used in the form of `C(a, b, c)` or `C(a=a, b=b, c=c)` or `C(_)`"
        end |> x -> (TARGET, reduce(∘, x, init=identity))

    end

    def_app_pattern(mod,
                predicate = (hd_obj, args) -> hd_obj === ctor,
                rewrite   = (tag, hd_obj, destruct_fields, mod) -> begin
                    TARGET, match_fields = mk_match(tag, hd_obj, destruct_fields, mod)
                    (@typed_as hd_obj) ∘ match_fields
                end,
                qualifiers = Set([qualifier]))

    # GADT syntax support!!!
    def_gapp_pattern(mod,
                predicate = (spec_vars, hd_obj, args) -> hd_obj === ctor,
                rewrite   = (tag, forall, spec_vars, hd_obj, destruct_fields, mod) -> begin
                    hd = :($hd_obj{$(spec_vars...)})
                    TARGET, match_fields = mk_match(tag, hd, destruct_fields, mod)
                    if isempty(forall)
                        @typed_as hd
                    else
                        function (body)
                            NAME = mangle(mod)
                            @format [TARGET, tag, body, hd] quote
                                @inline __L__ function NAME(TARGET :: hd) where {$(forall...)}
                                    body
                                end
                                @inline L function NAME(_)
                                    failed
                                end
                                NAME(tag)
                            end
                        end
                    end ∘ match_fields
                end,
                qualifiers = Set([qualifier]))
    ctor
end

macro as_record(qualifier, ctor)
    let mod = __module__, ctor = mod.eval(ctor)
        def_record(ctor, fieldnames(ctor), get_qualifier(qualifier, mod), mod)
    end
end

macro as_record(ctor)
    let mod = __module__, ctor = mod.eval(ctor)
        def_record(ctor, fieldnames(ctor), get_qualifier(:public, mod), mod)
    end
end

end


text:
# This file was generated by JuDoc, do not modify it. # hide
@show round.(randn(3), digits=2)


text:
import GitHub

function get_all_pull_requests(repo::GitHub.Repo,
                               state::String;
                               auth::GitHub.Authorization,
                               per_page::Integer = 100,
                               page_limit::Integer = 100)
    all_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    myparams = Dict("state" => state,
                    "per_page" => per_page,
                    "page" => 1)
    prs, page_data = GitHub.pull_requests(repo;
                                          auth=auth,
                                          params = myparams,
                                          page_limit = page_limit)
    append!(all_pull_requests, prs)
    while haskey(page_data, "next")
        prs, page_data = GitHub.pull_requests(repo;
                                              auth=auth,
                                              page_limit = page_limit,
                                              start_page = page_data["next"])
        append!(all_pull_requests, prs)
    end
    unique!(all_pull_requests)
    return all_pull_requests
end

_repos_are_the_same(::GitHub.Repo, ::Nothing) = false
_repos_are_the_same(::Nothing, ::GitHub.Repo) = false
_repos_are_the_same(::Nothing, ::Nothing) = false
function _repos_are_the_same(x::GitHub.Repo, y::GitHub.Repo)
    if x.name == y.name && x.full_name == y.full_name &&
                           x.owner == y.owner &&
                           x.id == y.id &&
                           x.url == y.url &&
                           x.html_url == y.html_url &&
                           x.fork == y.fork
       return true
    else
        return false
    end
end

function exclude_pull_requests_from_forks(repo::GitHub.Repo, pr_list::Vector{GitHub.PullRequest})
    non_forked_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    for pr in pr_list
        always_assert(_repos_are_the_same(repo, pr.base.repo))
        if _repos_are_the_same(repo, pr.head.repo)
            push!(non_forked_pull_requests, pr)
        end
    end
    return non_forked_pull_requests
end

function only_my_pull_requests(pr_list::Vector{GitHub.PullRequest}; my_username::String)
    _my_username_lowercase = lowercase(strip(my_username))
    n = length(pr_list)
    pr_is_mine = BitVector(undef, n)
    for i = 1:n
        pr_user_login = pr_list[i].user.login
        if lowercase(strip(pr_user_login)) == _my_username_lowercase
            pr_is_mine[i] = true
        else
            pr_is_mine[i] = false
        end
    end
    my_pr_list = pr_list[pr_is_mine]
    return my_pr_list
end

function create_new_pull_request(repo::GitHub.Repo;
                                 base_branch::String,
                                 head_branch::String,
                                 title::String,
                                 body::String,
                                 auth::GitHub.Authorization)
    params = Dict{String, String}()
    params["title"] = title
    params["head"] = head_branch
    params["base"] = base_branch
    params["body"] = body
    @info params
    @info repo
    @info auth
    result = GitHub.create_pull_request(repo; params = params, auth = auth)
    return result
end


text:
module Experiment

using JuLIP: Atoms, AbstractCalculator, neighbourlist,
             get_data, set_data!, has_data, JVecF, mat
import JuLIP: energy, forces

struct FastLJ <: AbstractCalculator
    rcut::Float64        # cut-off radius
    rcut2::Float64       # rcut^2
    r0::Float64          # spatial scale
    fcut::Float64        # ϕlj(rcut)
    dfcut::Float64       # ∂ϕlj(rcut) / ∂(rcut^2)
    rbuf::Float64   # buffer for the neighbourlist
end

FastLJ(; r0=1.0, rcut=1.01*r0, rbuf=0.5*rcut) =
   FastLJ( rcut,
           rcut^2,
           r0,
           (r0/rcut)^12 - 2 * (r0/rcut)^6,
           6 / r0^2 * (-(r0/rcut)^14 + (r0/rcut)^8),
           rbuf )

function update!(V::FastLJ, at::Atoms{Float64, Int})
    # first check whether we already have a neighbourlist that is still usable
    if has_data(at, :nlist_fastlj_X)
        Xold = get_data(at, :nlist_fastlj_X)::Vector{JVecF}
        d2 = 0.0
        for n = 1:length(Xold)
            d2 = max(d2, sum(abs2, at.X[n] - Xold[n]))
        end
        if d2 < V.rbuf^2
            i = get_data(at, :nlist_fastlj_i)::Vector{Int}
            j = get_data(at, :nlist_fastlj_j)::Vector{Int}
            return i, j
        end
    end
    # if not, then assemble a new neighbourlist
    nlist = neighbourlist(at, V.rcut)
    set_data!(at, :nlist_fastlj_X, copy(at.X))
    set_data!(at, :nlist_fastlj_i, nlist.i)
    set_data!(at, :nlist_fastlj_j, nlist.j)
    return nlist.i, nlist.j
end

function energy(V::FastLJ, at::Atoms{Float64, Int})
    i, j = update!(V, at)
    return energy_inner(V, at.X, i, j)
end

function energy_inner(V::FastLJ, X, i, j)
    E = 0.0
    for n = 1:length(i)
        @inbounds r2 = sum(abs2, X[i[n]]-X[j[n]])
        # evaluate LJ
        r2inv = (V.r0*V.r0)/r2
        r4inv = r2inv * r2inv
        r6inv = r2inv * r4inv
        r12inv = r6inv * r6inv
        lj1 = r12inv - 2 * r6inv
        # evaluate a cutoff
        lj2 = lj1 - V.fcut - V.dfcut * (r2 - V.rcut2)
        lj3 = lj2 * (r2 < V.rcut2)
        E += 0.5 * lj3
    end
    return E
end

forces(V::FastLJ, at::Atoms{Float64, Int}) =
    forces!(zeros(JVecF, length(at)), V, at)

function forces!(F, V::FastLJ, at::Atoms{Float64, Int})
    i, j = update!(V, at)
    return forces_inner!(F, V, at.X, i, j)
end

function forces_inner!(F, V::FastLJ, X, i, j)
    @assert length(F) == length(X)
    r02inv = 1.0 / V.r0^2
    for n = 1:length(i)
        R = X[i[n]]-X[j[n]]
        r2 = sum(abs2, R)
        r = sqrt(r2)
        # evaluate LJ
        r2inv = (V.r0*V.r0)/r2
        r4inv = r2inv * r2inv
        r6inv = r2inv * r4inv
        r8inv = r4inv * r4inv
        r14inv = r8inv * r6inv
        dlj1 = (6 * r02inv) * ( - r14inv + r8inv)
        # evaluate a cutoff
        dlj2 = (dlj1 - V.dfcut) * (r2 < V.rcut2)
        F[i[n]] -= dlj2 * R
        F[j[n]] += dlj2 * R
    end
    return F
end

end


text:
# Gaver-Stehfest method "Post-Widder"

# 1.
# Stehfest, Harald.
# Algorithm 368: Numerical inversion of Laplace transforms [D5].
# Communications of the ACM 13.1 (1970): 47-49.

# 2.
# Kuhlman, Kristopher L.
# Review of inverse Laplace transform algorithms for Laplace-space numerical approaches.
# Numerical Algorithms 63.2 (2013): 339-355.

# 3.
# Al-Shuaibi, Abdulaziz.
# Inversion of the laplace transform via post—widder formula.
# Integral Transforms and Special Functions 11.3 (2001): 225-232.

#= Compute Post-Widder coefficients Vk for N terms. Also referred to as Salzer summation weights or Stehfest coefficients.
Only depend on the approximation order (N) and the precision. M must be even but can be precomputed and used only once.
Coefficients get very large in magnitude and oscillate rapidly. Must be careful to avoid catastrophic cancellation for lower precision.
For  N < 18 double precision is usually ok but less accurate, utilizng defualt BigFloat precision and N = 36 usually provides sufficient accuracy.
=#
function _PWcoeffs(N::Integer)
    if isodd(N)
        N += 1
        @warn "N must be even... increment N += 1"
    end
    v = zeros(BigFloat, N)
    aux = zero(eltype(v))
    for k in 1:N
        for j in fld(k + 1, 2):minimum(Int, [k, div(N, 2)])
            aux = big(j)^(div(N, 2)) * factorial(big(2 * j))
            aux /= factorial(big(div(N, 2) - j)) * factorial(big(j)) * factorial(big(j - 1))
            aux /= factorial(big(k - j)) * factorial(big(2 * j - k))
            v[k] += aux
        end
        v[k] *= (-1)^(k + div(N, 2))
    end
    return v
end
"""
    postwid(func::Function, t::AbstractFloat, v::Array{AbstractFloat,1}=_PWcoeffs(N))

Evaluate the inverse Laplace transform of `func` at the point `t` using the Gaver-Stehfest "Post-Widder" algorithm.
Vk coefficients only depend on the number of expansion terms so can be computed once. N (which must be even) defaults to 18.
In general the accuracy in double precision is poor, so arbitrary precision is used throughout.

Increasing precision should be accompanied by an increase in the number of coefficients used.
Increasing precision without increasing number of coefficients will not yield better accuracy. The inverse is generally true as well.

This method is not robust to oscillating F(t) and must be smooth.
# Example

```jldoctest
julia> setprecision(53); InverseLaplace.postwid(s -> 1/(s + 1), 2.0)
0.1353356835639731
```
"""
function postwid(func::Function, t::AbstractFloat; v = _PWcoeffs(18))
    N = length(v)
    a = zero(eltype(v))
    for k in 1:N
        bk = convert(eltype(v), k)
        a += v[k] * func(bk * log(convert(eltype(a), 2)) / t)
    end
    return a * log(convert(eltype(a), 2)) / t
end
"""
    postwid(func::Function, t::AbstractArray, v::Array{AbstractFloat,1}=_PWcoeffs(N))

Evaluate the inverse Laplace transform of `func` over an array of `t` using the Gaver-Stehfest "Post-Widder" algorithm.
Computes coefficients once and calculates f(t) across available threads.
# Example

```jldoctest
julia> setprecision(53); postwid(s -> 1/(s + 1), 2.0:4.0)
 0.1353356835639731
 0.049786688998390505
 0.01831327193414956
```
"""
function postwid(f::Function, t::AbstractArray; v = _PWcoeffs(18))
    N = length(v)
    a = zeros(eltype(v), length(t))
    Threads.@threads for ind in eachindex(t)
        a[ind] = postwid(f, t[ind], v = v)
    end
    return a
end


text:
"""
    fuzzify(df::DataFrames.DataFrame; kwargs...)
This function adds "fuzziness" to an input `df`

# Keyword Arguments
- value::Symbol=VALUE_COL
- nsample::Integer=50
- distribution::Symbol=:normal
- factor
"""
function fuzzify(df;
    value = VALUE_COL,
    sample = SAMPLE_COL,
    units = :Units,
    label,
    sublabel,
    kwargs...,
)
    # If samples are not already included in the DataFrame, add them.
    if !(:Sample in propertynames(df))
        val = random_samples(df[:,value]; kwargs...)
        update_stop!(val)

        df = DataFrames.crossjoin(df, DataFrames.DataFrame(sample=>1:size(val,2)))
        df[!,value] .= vcat(val'...)
    end

    idx = DataFrames.Not([value, sample])
    return DataFrames.groupby(df, idx)
end

"""
"""
_distribution_error(args...) =  _option_error("distribution", [:normal,:uniform], args...)

""
function _option_error(option, allowed::String)
    throw(ArgumentError("Allowed values for $option: $allowed"))
end

function _option_error(option, allowed::AbstractArray)
    return _option_error(option, string((string.(allowed).*", ")...,)[1:end-2])
end

function _option_error(option, allowed, value)
    return if isempty(value)
        _option_error(option, allowed)
    else
        throw(ArgumentError("$option = $value. Allowed values: $allowed"))
    end
end


text:
using BinaryBuilder

json = String(read(ARGS[1]))
buff = IOBuffer(strip(json))
objs = []
while !eof(buff)
    push!(objs, BinaryBuilder.JSON.parse(buff))
end
json_obj = objs[1]

BinaryBuilder.cleanup_merged_object!(json_obj)
json_obj["dependencies"] = Dependency[dep for dep in json_obj["dependencies"] if !isa(dep, BuildDependency)]
BinaryBuilder.rebuild_jll_package(json_obj, download_dir="products", upload_prefix="terasakisatoshi/hellowrap_jll")


text:
using CairoMakie

"""
plot_mcmc_densities!(chain_X; n_columns=3)

Plots marginal density distributions for each parameter in the
constituent parameter vectors of `chain_X`, each representing a sample
generated by MCMC.
"""
function plot_mcmc_densities2!(fig, axes, chain_X, parameter_names;
                                n_columns = 3,
                                directory = pwd(),
                                filename = "mcmc_densities.png",
                                label = false,
                                color = (:blue, 0.5),
                                type = "hist",
                                show_means = false,
                                kwargs...)
    @assert type in ["hist", "density"]

    sample_means = []
    for (i, param_name) in enumerate(parameter_names)
        samples = getindex.(chain_X, i)

        ax = axes[i+1]
        if type == "hist"
            hist!(ax, samples; label, bins = 50, color, normalization = :pdf, kwargs...)
        elseif type == "density"
            density!(ax, samples; label, color, bandwidth=sqrt(var(samples))/15, kwargs...)
        end

        push!(sample_means, mean(samples))
        show_means && vlines!(ax,  [mean(samples)]; color=(color[1], 1.0), label="Sample mean", linestyle=nothing, linewidth=5)
    end

    fig[1,1] = Legend(fig, axes[end], nothing; framevisible=true,
                                                            tellheight = false,
                                                            tellwidth = false,
                                                            nbanks = 2)
    save(joinpath(directory, filename), fig)

    return sample_means
end

function plot_mcmc_densities2(chain_X, parameter_names;
                                n_columns = 3, kwargs...)

    N_axes = length(first(chain_X)) + 1
    n_rows = Int(ceil(N_axes / n_columns))

    fig = Figure(resolution = (500n_columns, 200n_rows))

    ax_coords = [(i, j) for i = 1:n_rows, j = 1:n_columns]

    # Reserve position [1,1] for the legend
    ax1 = Axis(fig[1,1])
    hidedecorations!(ax1)
    hidespines!(ax1)
    axes = [ax1]
    for (i, param_name) in enumerate(parameter_names)
        coords = ax_coords[i+1]
        ax = Axis(fig[coords...],
            xlabel = string(param_name),
            ylabel = "Density")
        ax.xlabelsize = 20
        push!(axes, ax)
    end

    sample_means = plot_mcmc_densities2!(fig, axes, chain_X, parameter_names; n_columns, kwargs...)

    return fig, axes
end


text:
# Unit tests for functions in src/algebra.jl file
using COSMO, Test, Random, LinearAlgebra, SparseArrays

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Algebra" begin

for T in UnitTestFloats

  rng = Random.MersenneTwister(9191)

  @testset "Algebra (T = $(T))" begin

    # clip()
    n = 50
    x = 10 * rand(rng, T, n) .- 5
    x[10] = -8
    x[20] = 8
    x_c = zeros(T, n)
    x_c2 = zeros(T, n)
    @. x_c = COSMO.clip(x, -2, 2)
    @test (maximum(x_c) <= 2. && minimum(x_c) >= -2.)
    @. x_c2 = COSMO.clip(x, -2, 2, -100, 100)
    @test (x_c2[10] == -100 && x_c2[20] == 100)

    # scaled_norm
    E = Diagonal(T[1, 2, 3.])
    v = T[-1, 5, 4]
    @test COSMO.scaled_norm(E, v, T(1)) == 23
    @test COSMO.scaled_norm(E, v, T(2)) == norm(T[-1; 10; 12], 2)
    @test COSMO.scaled_norm(E, v, T(Inf)) == 12
    @test_throws ArgumentError COSMO.scaled_norm(E, v, T(3));

    # col_norms!()
    v = zeros(T, 4)
    A = T[1. 2 3 4; 2 -1 30 4.1]
    COSMO.col_norms!(v, A)
    @test v == T[2; 2; 30; 4.1]

    A = sparse(T[1. 2 3 4; 2 -1 30 4.1])
    B = sparse(T[1 3 -2 4.2; -100 -1 -2 -100])
    COSMO.col_norms!(v, A, reset = true)
    COSMO.col_norms!(v, B, reset = false)
    @test v == T[100; 3; 30; 100]

    # row_norms!()
    A = Matrix(A')
    B = Matrix(B')
    COSMO.row_norms!(v, A, reset = true)
    COSMO.row_norms!(v, B, reset = false)
    @test v == T[100; 3; 30; 100]
    A = sparse(A)
    B = sparse(B)
    COSMO.row_norms!(v, A, reset = true)
    COSMO.row_norms!(v, B, reset = false)
    @test v == T[100; 3; 30; 100]

    # lmul!()
    L = Diagonal(rand(rng, T, 5))
    Mo = sprand(rng, T, 5, 5, 0.8)
    M = deepcopy(Mo)
    COSMO.lmul!(L, M)
    @test norm(M -  L * Mo, Inf) <= 1e-10

    M = deepcopy(Mo)
    COSMO.lmul!(I, M)
    @test M == Mo

    x = rand(T, 5)
    y = copy(x)
    COSMO.lmul!(L, x)
    @test norm(x - L * y, Inf ) <= 1e-10

    # rmul!()
    R = Diagonal(rand(rng, T, 5))
    Mo = sprand(rng, T, 5, 5, 0.8)
    M = deepcopy(Mo)
    COSMO.rmul!(M, R)
    @test norm(M -  Mo * R, Inf) <= 1e-10

    M = deepcopy(Mo)
    COSMO.rmul!(M, I)
    @test M == Mo

    # lrmul!()
    L = Diagonal(rand(rng, T, 5))
    Mo = rand(rng, T, 5, 5)
    M = deepcopy(Mo)
    R = Diagonal(rand(rng, T, 5))
    COSMO.lrmul!(L, M, R)
    @test norm(M -  L * Mo * R, Inf) <= 1e-7

    M = deepcopy(Mo)
    COSMO.lrmul!(L, M, I)
    @test norm(M -  L * Mo, Inf) <= 1e-10

    M = deepcopy(Mo)
    COSMO.lrmul!(I, M, I)
    @test M == Mo

    # symmetrize!()
    A = rand(rng, T, 3, 3)
    B = deepcopy(A)
    COSMO.symmetrize_upper!(A)
    @test Symmetric(A) == (0.5 * (B + B'))
    @test_throws AssertionError COSMO.symmetrize_upper!(rand(T, 2, 3));

    # is_pos_def()
    A = Symmetric(rand(rng, T, 4, 4))
    Apos = A + 4 * diagm(0 => ones(T, 4))
    Aneg = A - 4 * diagm(0 => ones(T, 4))
    @test COSMO.is_pos_def(Apos, T(1e-6))
    @test !COSMO.is_pos_def(Aneg, T(1e-6))
    @test COSMO.is_neg_def(Aneg, T(1e-6))
    @test !COSMO.is_neg_def(Apos, T(1e-6))

  end
end
end
nothing


text:
pts = rand(3, 3000)
μF = [5, 5, 5]

# Measure induced at resolution μF by the same resolution

# Accepted discrepancy between induced and original measures is sometimes
# nonzero, because the random starting distribution is *random* every
# time, and will yield slightly different measures every time.

μ_induced1a = μϵF_induced_by_ϵj(pts, μF , μF)
μ_induced1b = induced_measure(pts, μF , μF)

@test all(μ_induced1a .≈ μ_induced1b)

μ_induced2a = μϵF_induced_by_ϵj(pts, μF , [2, 6, 4])
μ_induced2b = induced_measure(pts, μF , [2, 6, 4])

@test all(μ_induced2a .≈ μ_induced2b)

@testset "Induced measure" begin
	@testset "Induced measures sum to 1" begin
		# Test both the precise-name function and the alias.
		@test abs(sum(μ_induced1a) - 1) < 1e-7
		@test abs(sum(μ_induced1b) - 1) < 1e-7
		@test abs(sum(μ_induced2a) - 1) < 1e-7
		@test abs(sum(μ_induced2b) - 1) < 1e-7
	end

	@testset "Induced measure at resultions differ from regular measure" begin
		# If we obtain the same measure by inducing the measure to a final
		# parition from a different partition, then something is very wrong.
		# Check that this is not the case.
		@test !(all(μ_induced1a .≈ μ_induced2a))
		@test !(all(μ_induced1b .≈ μ_induced2b))
	end

	@testset "Comparing regular to induced measure" begin
		# Measure computed regularly, either manually or using the
		# constructor that also stored info.
		μ_orig1 = invariantmeasure(transferoperator(pts, RectangularBinning(μF)))
		μ_orig2 = rectangularinvariantmeasure(pts, μF).measure;

		@testset "The same" begin
			@test sum(μ_orig1.dist) ≈ 1
			@test sum(μ_orig2.dist) ≈ 1
			@test length(μ_orig1.dist) == length(μ_orig2.dist)

			# We need to subset the induced measure at nonzero entries, because
			# the regularly computed measures are computed only on the bins that
			# are visited by the orbit, while the induced measure considers all
			# potential bins within the hypercube enclosing the points.
			# Next, we need to sort the points, because (for the same reason) the
			# order in which the bins are considered is different.
			#@show sort(μ_orig1.dist) .- sort(μ_induced1a[μ_induced1a .> 1e-7]) .< 1e-4
			#@test all(sort(μ_orig1.dist) .- sort(μ_induced1a[μ_induced1a .> 1e-7]) .< 1e-4)
			#@test all(sort(μ_orig1.dist) .- sort(μ_induced1b[μ_induced1b .> 1e-7]) .< 1e-4)
		end

		@testset "Different" begin
			# If induced from different resolution, the distributions will not be
			# identical.
			#sort(μ_orig1.dist) .≈ sort(μ_induced2a[μ_induced2a .> 0])
			@test !all(sort(μ_orig1.dist) .≈ sort(μ_induced2a[μ_induced2a .> 1e-7]))
			@test !all(sort(μ_orig1.dist) .≈ sort(μ_induced2b[μ_induced2b .> 1e-7]))
		end

	end
end


text:
using DataFrames, XLSX

# Files:
studentFile = "/Users/cfranken/GDrive/work/Caltech/OptionRepWork/TA2021/TA Preferences 2021-22 (Students).xlsx"
classFile   = "/Users/cfranken/GDrive/work/Caltech/OptionRepWork/TA2021/ClassList.xlsx"
facultyFile = "/Users/cfranken/GDrive/work/Caltech/OptionRepWork/TA2021/TA Preferences 2021-22 (Faculty)_May16.xlsx"

# Read into DataFrames:
df_students = DataFrame(XLSX.readtable(studentFile , "Form Responses 1")...)
df_faculty  = DataFrame(XLSX.readtable(facultyFile , "Form Responses 1")...)
df_classes  = DataFrame(XLSX.readtable(classFile, "Sheet1")...)

# Dicts that are helpful (and can be changed later)
const global StuClass       = Dict([("C1", "First Choice",), ("C2", "Second Choice"), ("C3", "Third Choice"), ("C4", "Fourth Choice"), ("C5", "Fifth Choice")])
const global StuInterest    = Dict([("C1", :"Relative interest in first choice"), ("C2", :"Relative interest in second choice"),("C3", :"Relative interest in third choice"),("C4", :"Relative interest in fourth choice"),("C5", :"Relative interest in fifth choice")])
const global FacTA_student  = Dict([("Class1", :"Course 1 requested TA's:  please list 3, comma separated, top choice first"), ("Class2", :"Course 2 requested TA's:  please list up to 3, comma separated, top choice first"), ("Class3", :"Course 3 requested TA's:  please list up to 3, comma separated, top choice first"), ("Class4", :"Course 4 requested TA's:  please list up to 3, comma separated, top choice first")])
const global FacTA_class    = Dict([("Class1", :"First course"), ("Class2", :"Second course"), ("Class3", :"Third course"), ("Class4", :"Fourth course")])

df_classes."Matched"   .= ""
df_classes."Unmatched" .= ""

max_stars_repo_name: cfranken/TA-matching

max_issues_repo_name: cfranken/TA-matching

max_forks_repo_name: cfranken/TA-matching


text:
using FiniteDifferences
import LinearAlgebra
import ChainRulesCore: frule, rrule
using ChainRulesCore
const RealOrComplex = Union{Real, Complex}

#Required for ChainRulesTestUtils
function FiniteDifferences.to_vec(M::GBMatrix)
    x, back = FiniteDifferences.to_vec(Matrix(M))
    function backtomat(xvec)
        M2 = GBMatrix(back(xvec))
        return mask(M2, M; structural=true)
    end
    return x, backtomat
end

function FiniteDifferences.to_vec(v::GBVector)
    x, back = FiniteDifferences.to_vec(Vector(v))
    function backtovec(xvec)
        v2 = GBVector(back(xvec))
        return mask(v2, v; structural=true)
    end
    return x, backtovec
end

# LinearAlgebra.norm doesn't like the nothings.
LinearAlgebra.norm(A::GBArray, p::Real=2) = norm(nonzeros(A), p)

# Broadcast b into the rows of A. WARNING: THIS DOES NOT MATCH JULIA.
function broadcast_emul!(C, A, b, op; mask = nothing, accum = nothing, desc = nothing)
    B = diagm(b)
    mul!(C, A, B, (any, op); mask, accum, desc)
    return C
end

function broadcast_emul(A, b, op; mask = nothing, accum = nothing, desc = nothing)
    B = diagm(b)
    mul(A, B, (any, op); mask, accum, desc)
end

function find_k(A, B::GBArray, op, minmax; mask = nothing, accum = nothing, desc = nothing)
    K = []
    for col ∈ axes(B, 2)
        intermediate = broadcast_emul(A, B[:, col], op; mask, accum, desc)
        push!(K, argminmax(intermediate, minmax, 2))
    end
    return hcat(K...)
end


text:
doc"""
    Weibull(α,θ)

The *Weibull distribution* with shape `α` and scale `θ` has probability density function

$f(x; \alpha, \theta) = \frac{\alpha}{\theta} \left( \frac{x}{\theta} \right)^{\alpha-1} e^{-(x/\theta)^\alpha},
    \quad x \ge 0$

```julia
Weibull()        # Weibull distribution with unit shape and unit scale, i.e. Weibull(1.0, 1.0)
Weibull(a)       # Weibull distribution with shape a and unit scale, i.e. Weibull(a, 1.0)
Weibull(a, b)    # Weibull distribution with shape a and scale b

params(d)        # Get the parameters, i.e. (a, b)
shape(d)         # Get the shape parameter, i.e. a
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Weibull distribution on Wikipedia](http://en.wikipedia.org/wiki/Weibull_distribution)

"""
immutable Weibull <: ContinuousUnivariateDistribution
    α::Float64   # shape
    θ::Float64   # scale

    function Weibull(α::Real, θ::Real)
    	@check_args(Weibull, α > zero(α) && θ > zero(θ))
    	new(α, θ)
    end
    Weibull(α::Real) = Weibull(α, 1.0)
    Weibull() = new(1.0, 1.0)
end

@distr_support Weibull 0.0 Inf

#### Parameters

shape(d::Weibull) = d.α
scale(d::Weibull) = d.θ

params(d::Weibull) = (d.α, d.θ)

#### Statistics

mean(d::Weibull) = d.θ * gamma(1.0 + 1.0 / d.α)
median(d::Weibull) = d.θ * logtwo ^ (1.0 / d.α)
mode(d::Weibull) = d.α > 1.0 ? (iα = 1.0 / d.α; d.θ * (1.0 - iα) ^ iα) : 0.0

var(d::Weibull) = d.θ^2 * gamma(1.0 + 2.0 / d.α) - mean(d)^2

function skewness(d::Weibull)
    μ = mean(d)
    σ2 = var(d)
    σ = sqrt(σ2)
    r = μ / σ
    gamma(1.0 + 3.0 / d.α) * (d.θ / σ)^3 - 3.0 * r - r^3
end

function kurtosis(d::Weibull)
    α, θ = params(d)
    μ = mean(d)
    σ = std(d)
    γ = skewness(d)
    r = μ / σ
    r2 = r^2
    r4 = r2^2
    (θ / σ)^4 * gamma(1.0 + 4.0 / α) - 4.0 * γ * r - 6.0 * r2 - r4 - 3.0
end

function entropy(d::Weibull)
    α, θ = params(d)
    0.5772156649015328606 * (1.0 - 1.0 / α) + log(θ / α) + 1.0
end

#### Evaluation

function pdf(d::Weibull, x::Float64)
    if x >= 0.0
        α, θ = params(d)
        z = x / θ
        (α / θ) * z^(α - 1.0) * exp(-z^α)
    else
        0.0
    end
end

function logpdf(d::Weibull, x::Float64)
    if x >= 0.0
        α, θ = params(d)
        z = x / θ
        log(α / θ) + (α - 1.0) * log(z) - z^α
    else
        -Inf
    end
end

zv(d::Weibull, x::Float64) = (x / d.θ) ^ d.α
xv(d::Weibull, z::Float64) = d.θ * z ^ (1.0 / d.α)

cdf(d::Weibull, x::Float64) = x > 0.0 ? -expm1(-zv(d, x)) : 0.0
ccdf(d::Weibull, x::Float64) = x > 0.0 ? exp(-zv(d, x)) : 1.0
logcdf(d::Weibull, x::Float64) = x > 0.0 ? log1mexp(-zv(d, x)) : -Inf
logccdf(d::Weibull, x::Float64) = x > 0.0 ? -zv(d, x) : 0.0

quantile(d::Weibull, p::Float64) = xv(d, -log1p(-p))
cquantile(d::Weibull, p::Float64) = xv(d, -log(p))
invlogcdf(d::Weibull, lp::Float64) = xv(d, -log1mexp(lp))
invlogccdf(d::Weibull, lp::Float64) = xv(d, -lp)

function gradlogpdf(d::Weibull, x::Float64)
    if insupport(Weibull, x)
        α, θ = params(d)
        (α - 1.0) / x - α * x^(α - 1.0) / (θ^α)
    else
        0.0
    end
end

#### Sampling

rand(d::Weibull) = xv(d, randexp())


text:
module MemoryBasedCF
    include("core.jl")
    export Config, Memory
    export memorize
    export userbased_scores, itembased_scores
    export userbased_rankings, itembased_rankings
end


text:
const _libsb = joinpath(dirname(@__FILE__),"..","deps","usr","lib", "libstemmer." * Libdl.dlext)
#@BinDeps.load_dependencies [:libstemmer=>:_libsb]
##
# character encodings supported by libstemmer
const UTF_8         = "UTF_8"

##
# lists the stemmer algorithms loaded
function stemmer_types()
    cptr = ccall((:sb_stemmer_list, _libsb), Ptr{Ptr{UInt8}}, ())
    (C_NULL == cptr) && error("error getting stemmer types")

    stypes = AbstractString[]
    i = 1
    while true
        name_ptr = unsafe_load(cptr, i)
        (C_NULL == name_ptr) && break
        push!(stypes, unsafe_string(name_ptr))
        i += 1
    end
    stypes
end

mutable struct Stemmer
    cptr::Ptr{Cvoid}
    alg::String
    enc::String

    function Stemmer(stemmer_type::String, charenc::String=UTF_8)
        cptr = ccall((:sb_stemmer_new, _libsb),
                    Ptr{Cvoid},
                    (Ptr{UInt8}, Ptr{UInt8}),
                    stemmer_type, charenc)

        if cptr == C_NULL
            if charenc == UTF_8
                error("stemmer '$(stemmer_type)' is not available")
            else
                error("stemmer '$(stemmer_type)' is not available for encoding '$(charenc)'")
            end
        end

        stm = new(cptr, stemmer_type, charenc)
        finalizer(release, stm)
        stm
    end
end

show(io::IO, stm::Stemmer) = println(io, "Stemmer algorithm:$(stm.alg) encoding:$(stm.enc)")

function release(stm::Stemmer)
    (C_NULL == stm.cptr) && return
    ccall((:sb_stemmer_delete, _libsb), Nothing, (Ptr{Cvoid},), stm.cptr)
    stm.cptr = C_NULL
    nothing
end

function stem(stemmer::Stemmer, bstr::String)::String
    sres = ccall((:sb_stemmer_stem, _libsb),
                Ptr{UInt8},
                (Ptr{UInt8}, Ptr{UInt8}, Cint),
                stemmer.cptr, bstr, sizeof(bstr))
    (C_NULL == sres) && error("error in stemming")
    slen = ccall((:sb_stemmer_length, _libsb), Cint, (Ptr{Cvoid},), stemmer.cptr)
    bytes = unsafe_wrap(Array, sres, Int(slen), own=false)
    String(copy(bytes))
end

function stem(stemmer::Stemmer, word::SubString{String})::String
    sres = ccall((:sb_stemmer_stem, _libsb),
                Ptr{UInt8},
                (Ptr{UInt8}, Ptr{UInt8}, Cint),
                stemmer.cptr, pointer(word.string)+word.offset, lastindex(word))
    (C_NULL == sres) && error("error in stemming")
    unsafe_string(sres)
    #slen = ccall((:sb_stemmer_length, _libsb), Cint, (Ptr{Cvoid},), stemmer.cptr)
    #bytes = pointer_to_array(sres, @compat(Int(slen)), false)
    #bytestring(bytes)
end


text:
module Queues

export Queue, enqueue!, dequeue!, peek

include("LinkedLists.jl")
using .LinkedLists

struct Queue{T}
  list::LinkedList{T}
end

Queue(t::Type{T}) where {T} = Queue{T}(LinkedList{T}())
Queue{T}() where {T} = Queue{T}(LinkedList{T}())
Queue() = Queue{Any}(LinkedList{Any}())

function enqueue!(q::Queue{T}, elem::T) where {T}
  pushfirst!(q.list, Node(elem))
  q
end

function dequeue!(q::Queue{T}) where {T}
  isempty(q) && throw(ArgumentError("Queue must not be empty"))
  pop!(q.list).data
end

function peek(q::Queue{T}) where {T}
  isempty(q) && throw(ArgumentError("Queue must not be empty"))
  (q.list[end]).data
end

function Base.length(q::Queue{T}) where {T}
  length(q.list)
end

function Base.isempty(q::Queue{T}) where {T}
  length(q) == 0
end

function Base.empty!(q::Queue{T}) where {T}
  empty!(q.list)

  q
end

function Base.eltype(q::Queue{T}) where {T}
  eltype(q.list)
end

function Base.iterate(q::Queue{T}, state::Union{Queue,Nothing} = q) where {T}
  (state == nothing || isempty(q)) && return nothing
  dequeue!(q), q
end

end

max_stars_repo_path: Chapter 6/Queues.jl

max_issues_repo_path: Chapter 6/Queues.jl

max_forks_repo_path: Chapter 6/Queues.jl


text:
# Gurobi environment and other supporting facilities

mutable struct Env
    ptr_env::Ptr{Cvoid}
end

# Note: `Env()` is an outer constructor because sometimes we want to
# make an `Env` object from a pointer returned from `GRBgetenv`, and
# sometimes we want to create a new `Env`.
function Env()
    a = Ref{Ptr{Cvoid}}()
    ret = @grb_ccall(
        loadenv,
        Cint,
        (Ptr{Ptr{Cvoid}}, Ptr{UInt8}),
        a, C_NULL
    )
    if ret == 10009
        error("Invalid Gurobi license")
    elseif ret != 0
        error("Failed to create environment (error $ret).")
    end
    # TODO(odow): no finalizer is set for Env, because when a model + Env
    # falls out of scope, Julia's GC will sometimes GC the env first, and
    # then the model. This causes an error. Thus, `Env`'s are finalized in
    # the `Model` constructor.
    # We should probably add a flag to differentiate between envs created
    # manually, and envs created by `Model()`. CPLEX.jl does something
    # similar.
    # finalizer(env, free_env)
    return Env(a[])
end

Base.unsafe_convert(ty::Type{Ptr{Cvoid}}, env::Env) = env.ptr_env::Ptr{Cvoid}

function is_valid(env::Env)
    env.ptr_env != C_NULL
end

function free_env(env::Env)
    if env.ptr_env != C_NULL
        @grb_ccall(freeenv, Nothing, (Ptr{Cvoid},), env.ptr_env)
        env.ptr_env = C_NULL
    end
end

function get_error_msg(env::Env)
    @assert env.ptr_env != C_NULL
    sz = @grb_ccall(geterrormsg, Ptr{UInt8}, (Ptr{Cvoid},), env.ptr_env)
    unsafe_string(sz)
end

# error

mutable struct GurobiError <: Exception
    code::Int
    msg::String

    function GurobiError(env::Env, code::Integer)
        new(convert(Int, code), get_error_msg(env))
    end
end


text:
struct FastWriteBuffer{V<:AbstractVector{UInt8}} <: IO
    data::V
    position::Base.RefValue{Int}
end

FastWriteBuffer(data::AbstractVector{UInt8}) = FastWriteBuffer(data, Ref(0))
FastWriteBuffer() = FastWriteBuffer(Vector{UInt8}())

@inline function ensureroom!(buf::FastWriteBuffer, n::Integer)
    resize!(buf.data, buf.position[] + n)
end

function Base.take!(buf::FastWriteBuffer)
    resize!(buf.data, buf.position[])
    buf.position[] = 0
    buf.data
end

@inline function Base.unsafe_write(buf::FastWriteBuffer, p::Ptr{UInt8}, n::UInt)
    ensureroom!(buf, n)
    unsafe_copyto!(pointer(buf.data, buf.position[] + 1), p, n)
    buf.position[] += n
    Int(n)
end

@inline function Base.write(buf::FastWriteBuffer, x::UInt8)
    ensureroom!(buf, 1)
    position = buf.position[] += 1
    buf.data[position] = x
    1
end

@inline Base.write(buf::FastWriteBuffer, x::Int8) = write(buf, reinterpret(UInt8, x))

@inline function Base.write(
            buf::FastWriteBuffer,
            x::Union{Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}) # TODO: more?
    n = Core.sizeof(x)
    ensureroom!(buf, n)
    position = buf.position[]
    @inbounds for i = Base.OneTo(n) # LLVM unrolls this loop on Julia 0.6.4
        position += 1
        buf.data[position] = x % UInt8
        x = x >>> 8
    end
    buf.position[] = position
    n
end

@inline Base.write(buf::FastWriteBuffer, x::Union{Float16, Float32, Float64}) = write(buf, reinterpret(Unsigned, x))

Base.position(buf::FastWriteBuffer) = buf.position[]
Base.isopen(buf::FastWriteBuffer) = true
Base.close(buf::FastWriteBuffer) = throw(MethodError("Unsupported operation"))
Base.eof(buf::FastWriteBuffer) = true
Base.iswritable(buf::FastWriteBuffer) = true
Base.isreadable(buf::FastWriteBuffer) = false


text:
using Traceur
using Tmatrix

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 3
k1_r = 1e7; k1_i = 0.0
k2_r = 1.5e7; k2_i = 1e3

# Using Complex numbers
@trace T_complex = calculate_Tmatrix_for_spheroid(rx, rz, n_max, Complex(k1_r, k1_i), Complex(k2_r, k2_i))


text:
function treepositions(hc::Hclust, useheight::Bool, orientation=:vertical)

    order = StatsBase.indexmap(hc.order)
    nodepos = Dict(-i => (float(order[i]), 0.0) for i in hc.order)

    xs = Array{Float64}(undef, 4, size(hc.merges, 1))
    ys = Array{Float64}(undef, 4, size(hc.merges, 1))

    for i in 1:size(hc.merges, 1)
        x1, y1 = nodepos[hc.merges[i, 1]]
        x2, y2 = nodepos[hc.merges[i, 2]]

        xpos = (x1 + x2) / 2
        ypos = useheight ?  hc.heights[i] : (max(y1, y2) + 1)

        nodepos[i] = (xpos, ypos)
        xs[:, i] .= [x1, x1, x2, x2]
        ys[:, i] .= [y1, ypos, ypos, y2]
    end
    if orientation == :horizontal
        return ys, xs
    else
        return xs, ys
    end
end

@recipe function f(hc::Hclust; useheight=true, orientation=:vertical )
    typeof(useheight) <: Bool || error("'useheight' argument must be true or false")

    legend --> false
    linecolor --> :black

    if orientation == :horizontal
        yforeground_color_axis --> :white
        ygrid --> false
        ylims --> (0.5, length(hc.order) + 0.5)
        yticks --> (1:nnodes(hc), string.(1:nnodes(hc))[hc.order])
        if useheight
            hs=sum(hc.heights)
            xlims --> (0, hs+hs*0.01)
        else
            xlims --> (0, Inf)
        end
        xshowaxis --> useheight
    else
        xforeground_color_axis --> :white
        xgrid --> false
        xlims --> (0.5, length(hc.order) + 0.5)
        xticks --> (1:nnodes(hc), string.(1:nnodes(hc))[hc.order])
        ylims --> (0, Inf)
        yshowaxis --> useheight
    end

    treepositions(hc, useheight, orientation)
end


text:
# PolyDerivEvaluator: used to compute the derivatives once instead of every call
struct PolyDerivEvaluator{F, D<:Tuple{Vararg{Polynomial}}}
    f::F
    derivs::D
end

@inline function PolyDerivEvaluator(p::Polynomial, ::Val{num_derivs}) where num_derivs
    PolyDerivEvaluator(p, make_poly_derivs(p, Val(num_derivs)))
end

make_poly_derivs(p::Polynomial, ::Val{0}) = ()
@inline function make_poly_derivs(p::Polynomial, ::Val{num_derivs}) where num_derivs
    p′ = SUP.derivative(p)
    (p′, make_poly_derivs(p′, Val(num_derivs - 1))...)
end

@generated function (pde::PolyDerivEvaluator)(θ, ::Val{num_derivs}) where num_derivs
    num_derivs >= 0 || return :(throw(ArgumentError("num_derivs must be nonnegative")))
    exprs = [:(pde.f(θ))]
    for i = 1 : num_derivs
        push!(exprs, :(pde.derivs[$i](θ)))
    end
    return :(tuple($(exprs...)))
end