text:
function elementary_matrix!(A,el::AbstractElement, u::PolynomialBasis,
                           v::PolynomialBasis, q::AbstractQuadratureRule;
                           f=(u,v)->(i,j,x̂,_,_)->u(x̂)[j]*v(x̂)[i])
    x̂,ŵ = q()
    x,w = q(el)
    m, n = length(v), length(u)
    for k in 1:length(w)
        for i in 1:m
            for j in 1:n
                A[i,j] += f(u, v)(i, j, x̂[k], el, x[k]) * w[k]
            end
        end
    end
    return A
end

function elementary_matrix_new!(A,el::AbstractElement, u::PolynomialBasis,
                           v::PolynomialBasis, q::AbstractQuadratureRule;
                           f=(u,v)->(i,j,x̂,_,_)->u(x̂)[j]*v(x̂)[i])
    x̂,ŵ = q()
    x,w = q(el)
    m, n = length(v), length(u)
    for k in 1:length(w)
        for i in 1:m
            for j in 1:n
                A[i,j] += f(u, v)(i, j, x̂[k], el, x[k]) * w[k]
            end
        end
    end
    return A
end

function elementary_matrix(el::AbstractElement, u::PolynomialBasis,
    v::PolynomialBasis, q::AbstractQuadratureRule;
    f=(u,v)->(i,j,x̂,_,_)->u(x̂)[j]*v(x̂)[i])
    x̂,ŵ = q()
    x,w = q(el)
    m, n = length(v), length(u)
    A = zeros(Float64,m,n)
    elementary_matrix!(A,el,u,v,q,f=f)
end

function elementary_matrix(el::AbstractElement, v::PolynomialBasis,
                           q::AbstractQuadratureRule;
                           f=(v)->(i,x̂,_,_)->v(x̂)[i])
    x̂,ŵ = q()
    x,w = q(el)
    m = length(v)
    b = zeros(Float64,m)
    for k in 1:length(w)
        for i in 1:m
            b[i] += f(v)(i, x̂[k], el, x[k]) * w[k]
        end
    end
    return b
end

# NOTE: for testing purposes regarding speed differences between generic and a
# specific elementary matrix kernel
function (u::LagrangeBasis)(q::AbstractQuadratureRule)
    x̂,_ = q()
    mapreduce(x->u(x),hcat,x̂)
end

function mass_matrix(el::AbstractElement,u::PolynomialBasis,v::PolynomialBasis,q::AbstractQuadratureRule)
    x̂,ŵ = q()
    U = u(q)
    V = v(q)
    μ = map(u->measure(el,u),x̂)
    Σ = diagm(μ)*diagm(ŵ)
    return V*Σ*transpose(U)
end

function mass_matrix_unrolled!(A,el::AbstractElement,u::PolynomialBasis,v::PolynomialBasis,q::AbstractQuadratureRule)
    x̂,ŵ = q()
    U = u(q)
    V = v(q)
    U *= diagm(ŵ) # put precomputed weights here. Should be computed only once.
    for k in 1:length(ŵ)
        J = @inbounds jacobian(el,x̂[k])
        μ = abs(det(J))
        for j in 1:size(U,2)
            @inbounds uj = U[j,k]*μ
            for i in 1:size(V,2)
                @inbounds A[i,j] += uj*V[i,k]
            end
        end
    end
    return A
end

function (gu::gradLagrangeBasis)(q::AbstractQuadratureRule)
    x̂,_ = q()
    mapreduce(x->gu(x),hcat,x̂)
end

function stiffness_matrix_unrolled!(A,el::AbstractElement,u::PolynomialBasis,v::PolynomialBasis,q::AbstractQuadratureRule)
    x̂,ŵ = q()
    ∇U = grad(u)(q)
    ∇V = grad(v)(q)
    ∇U *= diagm(ŵ) # put the reference weights in one of the precomputed matrices
    for k in 1:length(ŵ)
        J  = jacobian(el,x̂[k])
        iJ = inv(J)
        μ  = abs(det(J))
        Σ  = iJ*transpose(iJ)*μ
        for i in 1:size(∇V,2)
            for j in 1:size(∇U,2)
                A[i,j] += transpose(∇V[i,k]) * Σ * ∇U[j,k]
            end
        end
    end
    return A
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Autosampler, ImagineFormat, CSV, DataFrames
using Test, Random

if !isfile("recording.ai")
    include("create_ai.jl")
end

@testset "Autosampler.jl" begin
    rows = ["chemical1_10uM,RA1", "chemical1_01uM,RA2"]
    rowid = [split(row, ',')[1] for row in rows]
    # The sequence of random numbers might change across Julia versions,
    # so determine what the sequence is
    Random.seed!(1234321)
    p = vcat(randperm(2), randperm(2), randperm(2), randperm(2))
    for stimfile in ("stimuli_tabs.csv", "stimuli_comma_quotes.csv")
        Random.seed!(1234321)
        mktemp() do pathcsv, iocsv
            # Generate sequence file
            stim_randomize(iocsv, stimfile, 4)
            str = readlines(pathcsv)
            @test str == rows[p]
            # Test insertion into .imagine file
            pathi, io = mktemp()
            close(io)
            update_imagine("recording.imagine", pathcsv; updatedfile=pathi, um_per_pixel=0.33)
            header = ImagineFormat.parse_header(pathi)
            @test split(header["stimulus sequence"], '\$') == rowid[p]
            @test header["stimulus scan hi"] == 1000:1000:8000
            @test length(header["stimulus frame hi"]) == 8
            rm(pathi)
        end
    end
end

ext: jl


text:
module STBParameters

   import TBComponents: RateBands, WEEKS_PER_YEAR
   import Parameters: @with_kw # FIXME maybe use Base.kwdef; see: https://github.com/JuliaLang/julia/pull/29316

   import Utils
   import JSON

   using Definitions

   export IncomeTaxSys, NationalInsuranceSys, TaxBenefitSystem
   export weeklyise!, annualise!, fromJSON, get_default_it_system

   ## TODO Use Unitful to have currency weekly monthly annual counts as annotations
   # using Unitful

   Fuel_Dict = Dict{Fuel_Type,Real}
   Default_Fuel_Dict_2020_21 = Fuel_Dict(
         Missing_Fuel_Type=>0.1,
         No_Fuel=>0.1,
         Other=>0.1,
         Dont_know=>0.1,
         Petrol=>0.25, # dunno
         Diesel=>0.37,
         Hybrid_use_a_combination_of_petrol_and_electricity=>0.16,
         Electric=>0.02,
         LPG=>0.02,
         Biofuel_eg_E85_fuel=>0.02 )

   @with_kw mutable struct IncomeTaxSys
      non_savings_rates :: RateBands = [19.0,20.0,21.0,41.0,46.0]
      non_savings_thresholds :: RateBands = [2_049.0, 12_444.0, 30_930.0, 150_000.0]
      non_savings_basic_rate :: Integer = 2 # above this counts as higher rate

      savings_rates  :: RateBands = [0.0, 20.0, 40.0, 45.0]
      savings_thresholds  :: RateBands = [5_000.0, 37_500.0, 150_000.0]
      savings_basic_rate :: Integer = 2 # above this counts as higher rate

      dividend_rates :: RateBands = [0.0, 7.5,32.5,38.1]
      dividend_thresholds :: RateBands = [2_000.0, 37_500.0, 150_000.0]
      dividend_basic_rate :: Integer = 2 # above this counts as higher rate

      personal_allowance          = 12_500.00
      personal_allowance_income_limit = 100_000.00
      personal_allowance_withdrawal_rate = 50.0
      blind_persons_allowance     = 2_450.00

      married_couples_allowance   = 8_915.00
      mca_minimum                 = 3_450.00
      mca_income_maximum          = 29_600.00
      mca_credit_rate             = 10.0
      mca_withdrawal_rate         = 50.0

      marriage_allowance          = 1_250.00
      personal_savings_allowance  = 1_000.00

      # FIXME better to have it straight from
      # the book with charges per CO2 range
      # and the data being an estimate of CO2 per type
      company_car_charge_by_CO2_emissions = Default_Fuel_Dict_2020_21
      fuel_imputation = 24_100.00

      #
      # pensions
      #
      pension_contrib_basic_amount = 3_600.00
      pension_contrib_annual_allowance = 40_000.00
      pension_contrib_annual_minimum = 10_000.00
      pension_contrib_threshold_income = 110_000.00
      it.pension_contrib_withdrawal_rate = 50.0
   end

   function annualise!( it :: IncomeTaxSys )
      it.non_savings_rates .*= 100.0
      it.savings_rates .*= 100.0
      it.dividend_rates .*= 100.0
      it.personal_allowance_withdrawal_rate *= 100.0
      it.non_savings_thresholds .*= WEEKS_PER_YEAR
      it.savings_thresholds .*= WEEKS_PER_YEAR
      it.dividend_thresholds .*= WEEKS_PER_YEAR
      it.personal_allowance *= WEEKS_PER_YEAR
      it.blind_persons_allowance *= WEEKS_PER_YEAR
      it.married_couples_allowance *= WEEKS_PER_YEAR
      it.mca_minimum *= WEEKS_PER_YEAR
      it.marriage_allowance *= WEEKS_PER_YEAR
      it.personal_savings_allowance *= WEEKS_PER_YEAR
      it.pension_contrib_basic_amount *= WEEKS_PER_YEAR

      it.mca_income_maximum       *= WEEKS_PER_YEAR
      it.mca_credit_rate             *= 100.0
      it.mca_withdrawal_rate         *= 100.0
      for k in it.company_car_charge_by_CO2_emissions
         it.company_car_charge_by_CO2_emissions[k.first] *= WEEKS_PER_YEAR
      end
   end

   function weeklyise!( it :: IncomeTaxSys )

      it.non_savings_rates ./= 100.0
      it.savings_rates ./= 100.0
      it.dividend_rates ./= 100.0
      it.personal_allowance_withdrawal_rate /= 100.0
      it.non_savings_thresholds ./= WEEKS_PER_YEAR
      it.savings_thresholds ./= WEEKS_PER_YEAR
      it.dividend_thresholds ./= WEEKS_PER_YEAR
      it.personal_allowance /= WEEKS_PER_YEAR
      it.blind_persons_allowance /= WEEKS_PER_YEAR
      it.married_couples_allowance /= WEEKS_PER_YEAR
      it.mca_minimum /= WEEKS_PER_YEAR
      it.marriage_allowance /= WEEKS_PER_YEAR
      it.personal_savings_allowance /= WEEKS_PER_YEAR
      it.mca_income_maximum       /= WEEKS_PER_YEAR
      it.mca_credit_rate             /= 100.0
      it.mca_withdrawal_rate         /= 100.0
      for k in it.company_car_charge_by_CO2_emissions
         it.company_car_charge_by_CO2_emissions[k.first] /= WEEKS_PER_YEAR
      end
      it.pension_contrib_basic_amount /= WEEKS_PER_YEAR

   end

   function get_default_it_system(
      ;
      year     :: Integer=2019,
      scotland :: Bool = true,
      weekly   :: Bool = true )::Union{Nothing,IncomeTaxSys}
      it = nothing
      if year == 2019
         it = IncomeTaxSys()
         if ! scotland
            it.non_savings_rates = [20.0,40.0,45.0]
            it.non_savings_thresholds = [37_500, 150_000.0]
            it.non_savings_basic_rate = 1
         end
         if weekly
            weeklyise!( it )
         end
      end
      it
   end

   function to_rate_bands( a :: Vector ) :: RateBands
      n = size( a )[1]
      rb :: RateBands = zeros(n)
      for i in 1:n
         rb[i] =  Real(a[i])
      end
      rb
   end

   function to_fuel_charges( d :: Dict ) :: Fuel_Dict
      fd = Fuel_Dict()
      for i in instances(Fuel_Type)
         k = String(Symbol(i))
         fd[i] = d[k]
      end
      fd
   end

   """
   Map from
   """
   function fromJSON( json :: Dict ) :: IncomeTaxSys
      it = IncomeTaxSys()
      println( typeof(json["non_savings_thresholds"]))
      it.non_savings_rates = to_rate_bands( json["non_savings_rates"] )
      it.non_savings_thresholds  = to_rate_bands( json["non_savings_thresholds"] )
      it.non_savings_basic_rate = json["non_savings_basic_rate"]

      it.savings_rates = to_rate_bands( json["savings_rates"] )
      it.savings_thresholds = to_rate_bands( json["savings_thresholds"] )
      it.savings_basic_rate = json["savings_basic_rate"]

      it.dividend_rates = to_rate_bands( json["dividend_rates"] )
      it.non_savings_thresholds = to_rate_bands( json["non_savings_thresholds"] )
      it.dividend_basic_rate = json["dividend_basic_rate"]

      it.savings_thresholds = to_rate_bands( json["savings_thresholds"] )
      it.dividend_thresholds = to_rate_bands( json["dividend_thresholds"] )
      it.personal_allowance = json["personal_allowance"]
      it.personal_allowance_income_limit = json["personal_allowance_income_limit"]
      it.personal_allowance_withdrawal_rate = json["personal_allowance_withdrawal_rate"]
      it.blind_persons_allowance = json["blind_persons_allowance"]
      it.married_couples_allowance = json["married_couples_allowance"]
      it.mca_minimum = json["mca_minimum"]
      it.marriage_allowance = json["marriage_allowance"]
      it.personal_savings_allowance = json["personal_savings_allowance"]

      it.mca_income_maximum = json["mca_income_maximum"]
      it.mca_credit_rate = json["mca_credit_rate"]
      it.mca_withdrawal_rate = json["mca_withdrawal_rate"]
      ## CAREFUL!
      it.fuel_imputation = json["fuel_imputation"]
      it.company_car_charge_by_CO2_emissions =
         to_fuel_charges(json["company_car_charge_by_CO2_emissions"])
      it.pension_contrib_basic_amount = json["pension_contrib_basic_amount"]
      it.pension_contrib_annual_allowance = json["pension_contrib_annual_allowance"]
      it.pension_contrib_annual_minimum = json["pension_contrib_annual_allowance"]
      it.pension_contrib_threshold_income = json["pension_contrib_threshold_income"]
      it.pension_contrib_withdrawal_rate = json["pension_contrib_withdrawal_rate"]
      it
   end

   @with_kw mutable struct NationalInsuranceSys
      rates :: RateBands = [1.0]
   end

   @with_kw mutable struct TaxBenefitSystem
      name :: AbstractString = "Scotland 2919/20"
      it   :: IncomeTaxSys = IncomeTaxSys()
      ni   :: NationalInsuranceSys = NationalInsuranceSys()
   end

   function save( filename :: AbstractString, sys :: TaxBenefitSystem )
       JSON.print( filename, sys )
   end

   # include( "../default_params/default2019_20.jl")
   # defsys = load()

end


text:
# -----------------------------------------------------------
#
# Constructing the Hierarchical Discontinuous Galerkin Basis
#
# -----------------------------------------------------------

# Efficiency criticality: LOW
# Computations only performed once

# -----------------------------------------------------
# Defining the Inner Product
# -----------------------------------------------------

# The coordinate convention here is to have vectors of length n even
# The first half has component i is the respective coefficient of x^i
# THe second half has component n/2 + i is the coefficient of sgn(x)*x^i

function product_matrix(i::Int, j::Int, n::Int)
    # This is the inner product of <x^i, x^j> when i,j < n/2
    # or the appropriate reflection when they are >=
    # Exact expressions only. No NIntegrations
    k= Int(round(n/2))
    if i < k && j < k
        return (1 + (-1)^(i + j))/(1 + i + j)
    elseif i >= k && j< k
        return (1 - (-1)^((i-k) + (j)))/(1 + (i-k) + (j))
    elseif i < k && j >= k
        return product_matrix(j, i, n)
    else
        return (1 + (-1)^((i-k) + (j-k)))/(1 + (i-k) + (j-k))
    end
end

function inner_product(v1::AbstractArray{T}, v2::AbstractArray{T}) where T <: Real
    value=zero(eltype(v1)) #Get 0 of the same type as v1
    n=length(v1)
    for i in 1:n
        for j in 1:n
            if v1[i]==0 || v2[j]==0
                continue
            else
                value += product_matrix(i-1,j-1,n)*v1[i]*v2[j]
            end
        end
    end
    return value
end

function inner_product(v1::AbstractArray{T}, j::Int) where T <: Real #we consider x^(j-1)
    value=zero(eltype(v1)) #Get 0 of the same type as v1
    n=length(v1)
    for i in 1:n
        value += product_matrix(i-1,j-1,n)*v1[i]
        #note we need to shift appropriately
    end
    return value
end

# ------------------------------------------------------
# Defining Gram-Schmidt and forming Legendre Polynomials
# ------------------------------------------------------

# The gram schmidt process on a set of vectors
# I think using an array of arrays is easiest here
function gram_schmidt(Q_initial::Array{Array{T, 1}, 1}) where T <: Real
    n = length(Q_initial[1])
    k = Int(round(n/2))
    # n = 2k, k polys of degree up to k-1 and then k polys according to f
    Q_final = deepcopy(Q_initial)
    # we'll return a modified version of the initial without modifying Q_initial
    for i = 1:k
        for j = 1:i-1
            proj = inner_product(Q_initial[i],Q_final[j])/inner_product(Q_final[j], Q_final[j])
            #This gives us the number corresponding to the projection along Q_final[j]
            Q_final[i] -= proj * Q_final[j]
            #This subtracts it out from Q_final[i] to orthogonalize
        end
        Q_final[i] /= sqrt(inner_product(Q_final[i], Q_final[i]))
        #Normalize every time
    end
    return Q_final
end

#We can now form the Legendre polynomials
function legendre(k::Int)
    Q = [[i==j ? 1.0 : 0.0 for i = 1:2*(k+1)] for j = 1:(k+1)]
    #start with just the x^i basis for i = 0 to k
    Q = gram_schmidt(Q)
    #Orthgonalize, and now we have the Legendre polynomials
    return Q
end

# I want to be able to do this completely analytically, using fractions
# so that I can avoid any numerical error in this process

# -----------------------------------------------------------
# Making a basis of functions have the first k moments vanish
# -----------------------------------------------------------

function orthogonalize_1(Q_initial::Array{Array{T, 1}, 1}) where T <: Real
    n = length(Q_initial[1])
    k = Int(round(n/2))
    # n = 2k, because basis includes the f(j,x) functions in addition to x^j
    Q_final = deepcopy(Q_initial)
    # As before
    legendre_polys=legendre(k-1)
    # We need an orthogonal basis for 1..x^(k-1) for the projection to work
    for i = 1:k
        for j in 1:k
            proj =  inner_product(Q_initial[i],legendre_polys[j])/
                        inner_product(legendre_polys[j],legendre_polys[j])
            Q_final[i] -= proj*legendre_polys[j] #subtract projection
            #project out v[i] by each element of this orthogonal basis for x^j
        end
    end
    return Q_final
end

# -----------------------------------------------------
# Make some functions have higher vanishing moments
# -----------------------------------------------------

# This function will make k-1 of the basis vectors orth to x^k,
# k-2 to x^k+1  all the way to 1 vector orth to x^2k-2

function orthogonalize_2(Q_initial::Array{Array{T, 1}, 1}) where T <: Real
    n = length(Q_initial[1])
    k = Int(round(n/2))
    Q_final = deepcopy(Q_initial)
    for i = 1:k-1
        fi=copy(Q_initial[i])
        # We assume this isn't perp to x^(k+i-1) and subtract it from the next ones.
        # We know that it isn't perp by parity considerations on the degree! :)
        # (Thank goodness, cuz otherwise each time, we'd have to rearrange the basis
        # until we found a non-perp function and subtract by that one)
        for j = i+1:k
            a= inner_product(Q_final[j],k+i)/inner_product(fi, k+i)
            Q_final[j] -= a * fi
        end
    end
    return Q_final
end

# Standard gram-schmidt process, starting at the end
# (this is important, because the last function is orthogonal to a lot of higher
# polynomials, and we don't want to do anything other than normalize it,
# with similar reasoning for the penultimate, etc. functions)

function gram_schmidt_rev(Q_initial::Array{Array{T, 1}, 1}) where T <: Real
    n = length(Q_initial[1])
    k = Int(round(n/2))
    Q_final = [[0.0 for i in 1:n] for j in 1:k]
    for i = k:-1:1    #note we're going in reverse
        fi = copy(Q_initial[i])  #we won't modify the original, so we copy
        Q_final[i] = fi          #start with Q_final = f_i, and we'll subtract projections
        for j = i+1:k #because we're going in reverse
            proj = inner_product(fi,Q_final[j])/inner_product(Q_final[j], Q_final[j])
            #find the projection
            Q_final[i] -= proj * Q_final[j]    #project out direction j
        end
        Q_final[i] /= sqrt(inner_product(Q_final[i], Q_final[i])) #now normalize
    end
    return Q_final
end

# -----------------------------------------------------
# TODO: Lastly, perform a rotation so that there
# is only one 'dicontinuous' basis element,
# namely the last one
# -----------------------------------------------------

function rotate_discontinuity(Q_initial::Array{Array{T, 1}, 1}) where T <: Real
    # To be implemented
end

# -----------------------------------------------------
# All together, for the final result:
# -----------------------------------------------------

function dg_basis(k::Int)
    Q = [[j==(i-k) ? 1.0 : 0.0 for i in 1:2*k] for j in 1:k]
    Q = orthogonalize_1(Q)
    Q = orthogonalize_2(Q)
    Q = gram_schmidt_rev(Q)
    return Q
end


text:
function sdplib_data(path)
    # Read data from file
    data = DelimitedFiles.readdlm(path, use_mmap=true)

    # Parse SDPLIB data
    m = data[1, 1]
    if isa(data[3, 1], Float64) || isa(data[3, 1], Int64)
        blks = data[3, :]
        for elem = 1:length(blks)
            if blks[elem] == ""
                blks = blks[1:elem-1]
                break
            end
        end
    else
        blks = parse.(Float64, split(data[3, 1][2:end - 1], ","))
    end
    cum_blks = pushfirst!(cumsum(blks), 0)
    if isa(data[4, 1], Float64) || isa(data[4, 1], Int64)
        c = data[4, :]
    else
        c = [parse(Float64,string) for string in split(data[4, 1][2:end - 1], ",")]
    end
    n = abs(cum_blks[end])
    n = length(c)
    F = Dict(i => SparseArrays.spzeros(n, n) for i = 0:m)
    for k=5:size(data)[1]
        idx = cum_blks[data[k, 2]]
        # if data[k, 2] == 1
        #     idx = 0
        # else
        #     idx = 161
        # end

        i, j = data[k, 3] + idx, data[k, 4] + idx
        if data[k, 1] == 0
            F[0][i, j] = - data[k, 5]
            F[0][j, i] = - data[k, 5]
        else
            F[data[k, 1]][i, j] = data[k, 5]
            F[data[k, 1]][j, i] = data[k, 5]
        end
    end
    return n, m, F, c
end
function sdplib_eval(F,c,n,m,XX)
    rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-10])
    @show rank
    @show tr(F[0] * XX)
    nothing
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libnl_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libnl")
JLLWrappers.@generate_main_file("libnl", UUID("7c700954-19d3-5208-81e2-8fa5fe7c0bd8"))
end  # module libnl_jll


text:
struct MyComplex
    real::Float64
    imag::Float64
end

a = MyComplex(1.0, 2.0)

println(a.real, a.imag)

max_stars_repo_name: liushooter/OneDayOneCommit

max_issues_repo_name: liushooter/OneDayOneCommit

max_forks_repo_name: liushooter/OneDayOneCommit


text:
# test for variations in DOF solving,  #227 #316 #430

# using Revise

using IncrementalInference
using Test

##

@testset "Test for variations in N while solving" begin

##

fg = generateCanonicalFG_CaesarRing1D(graphinit=true)

pts_ = approxConv(fg, :x0x1f1, :x1, N=101)
@test length(pts_) == 101
# if length(pts_) != 101
#   @warn "approxConv not adhering to N=101 != $(length(pts_)), see issue #105"
# end

##

# Change to N=150 AFTER constructing the graph, so solver must update the belief sample values during inference
getSolverParams(fg).N = 150
# getSolverParams(fg).multiproc = false
# getSolverParams(fg).async = false
smtasks = Task[]
tree = solveTree!(fg; smtasks, recordcliqs=ls(fg));

##

pts_ = getBelief(fg, :x1) |> getPoints
@test length(pts_) == 150
@test length(pts_[1]) == 1

##

# test with change for second solve
getSolverParams(fg).N = 200
tree = solveTree!(fg)

pts_ = getBelief(fg, :x1) |> getPoints
@test length(pts_) == 200
@test length(pts_[1]) == 1

println("test making N smaller than current")
getSolverParams(fg).N = 99
tree = solveTree!(fg)

pts_ = getBelief(fg, :x1) |> getPoints
@warn "removing older solve N size test, likely to be reviewed and updated to new workflow in the future"
@test length(pts_) == 99
@test length(pts_[1]) == 1

##

end

#


text:
using Slomo.CosmologyTools: Ωm, ρcrit, ρm, default_cosmo
using Slomo.Constants: rho_crit_h2

using Test

cosmo = default_cosmo

z = [0.0, 0.1, 0.5]
@test size(Ωm.(z; cosmo = cosmo)) == size(z)
@test size(ρcrit.(z; cosmo = cosmo)) == size(z)
@test size(ρm.(z; cosmo = cosmo)) == size(z)

z = 0
@test Ωm(z; cosmo = cosmo) == cosmo.Ω_m
@test ρcrit(z; cosmo = cosmo) == rho_crit_h2 * cosmo.h ^ 2
@test ρm(z; cosmo = cosmo) == cosmo.Ω_m * rho_crit_h2 * cosmo.h ^ 2

z = 0.0
@test Ωm(z; cosmo = cosmo) == cosmo.Ω_m
@test ρcrit(z; cosmo = cosmo) == rho_crit_h2 * cosmo.h ^ 2
@test ρm(z; cosmo = cosmo) == cosmo.Ω_m * rho_crit_h2 * cosmo.h ^ 2


text:
using AzureBlob
using Dates
using TimeZones

using Test

@testset "Black box signature" begin
	timestamp = DateTime(2019, 1, 1, 12, 34, 56) |> AzureBlob.http_date

	getsig = storage_signature("blob", "directory", "container", "storageaccount", "storagekeyxx", timestamp)
	@test getsig == "QDFeyZXEbQXqHJ3Pu4WByyWo1zU52SRC/GBktr+8SdU="

	putsig = storage_signature("content", "blob", "directory", "container", "storageaccount", "storagekeyxx", timestamp, "application/json")
	@test putsig == "Zxe1fAisDpvx84CpCfV4u4ptem78d0ikwOKsmvNe4lg="
end

@testset "Time zones" begin
	timestamp = DateTime(2019, 1, 1, 12, 34, 56)
	@test http_date(timestamp) == "Tue, 01 Jan 2019 12:34:56 GMT"

	timestamp = ZonedDateTime(2019, 1, 1, 12, 34, 56, tz"UTC")
	@test http_date(timestamp) == "Tue, 01 Jan 2019 12:34:56 GMT"

	timestamp = ZonedDateTime(2019, 1, 1, 12, 34, 56, tz"America/New_York")
	@test http_date(timestamp) == "Tue, 01 Jan 2019 17:34:56 GMT"

	timestamp = ZonedDateTime(2019, 1, 1, 12, 34, 56, tz"Europe/Copenhagen")
	@test http_date(timestamp) == "Tue, 01 Jan 2019 11:34:56 GMT"
end

max_stars_repo_name: robertdj/AzureBlob

max_issues_repo_name: robertdj/AzureBlob

max_forks_repo_name: robertdj/AzureBlob


text:
using ReinforcementLearningZoo
using Test
using ReinforcementLearningBase
using ReinforcementLearningCore
using ReinforcementLearningEnvironments
using Flux
using Statistics
using Random
using StableRNGs
using OpenSpiel

# used for OpenSpielEnv("kuhn_poker")
function get_optimal_kuhn_policy(α = 0.2)
    TabularRandomPolicy(
        table = Dict(
            "0" => [1 - α, α],
            "0pb" => [1.0, 0.0],
            "1" => [1.0, 0.0],
            "1pb" => [2.0 / 3.0 - α, 1.0 / 3.0 + α],
            "2" => [1 - 3 * α, 3 * α],
            "2pb" => [0.0, 1.0],
            "0p" => [2.0 / 3.0, 1.0 / 3.0],
            "0b" => [1.0, 0.0],
            "1p" => [1.0, 0.0],
            "1b" => [2.0 / 3.0, 1.0 / 3.0],
            "2p" => [0.0, 1.0],
            "2b" => [0.0, 1.0],
        ),
    )
end

# used for julia version KuhnPokerGame
function get_optimal_kuhn_policy(env::KuhnPokerEnv; α = 0.2)
    TabularRandomPolicy(
        table = Dict(
            (:J,) => [1 - α, α],
            (:J, :pass, :bet) => [1.0, 0.0],
            (:Q,) => [1.0, 0.0],
            (:Q, :pass, :bet) => [2.0 / 3.0 - α, 1.0 / 3.0 + α],
            (:K,) => [1 - 3 * α, 3 * α],
            (:K, :pass, :bet) => [0.0, 1.0],
            (:J, :pass) => [2.0 / 3.0, 1.0 / 3.0],
            (:J, :bet) => [1.0, 0.0],
            (:Q, :pass) => [1.0, 0.0],
            (:Q, :bet) => [2.0 / 3.0, 1.0 / 3.0],
            (:K, :pass) => [0.0, 1.0],
            (:K, :bet) => [0.0, 1.0],
        ),
    )
end

@testset "ReinforcementLearningZoo.jl" begin
    include("cfr/cfr.jl")
end


text:
export @class

macro class(Type, MethodList = [])
    esc(quote
        local methodMap = Set{Symbol}($MethodList)

        function Base.propertynames(obj::$Type)
            return vcat(collect(keys(methodMap)), collect(fieldnames($Type)))
        end

        function Base.hasproperty(obj::$Type, sym::Symbol)
            if (hasfield($Type, sym))
                return true
            else
                return in(sym, methodMap)
            end
        end

        function Base.getproperty(obj::$Type, sym::Symbol)
            #if in(sym, methodMap)
            if hasfield($Type, sym)
                return getfield(obj, sym)
            else
                # assumes unknown properties are methods
                return (args...;kwargs...) -> begin
                    getfield(@__MODULE__, sym)(obj, args...;kwargs...)
                end
            end
            # else
            #     return getfield(obj, sym)
            # end
        end
    end)
end


text:
using BHAtp, Statistics

ProjDir = @__DIR__
!isdir(joinpath(ProjDir, "plots")) && mkdir(joinpath(ProjDir, "plots"))

# Actual data

incs = [5,10, 20, 30, 40, 50]
tpf = DataFrame(incl=incs)
tpf[!, :tp] = [1.95, 2.70, 3.60, 4.30, 4.80, 5.05]

actual = DataFrame(
  incls = [4, 6.5, 7.5, 12, 13, 15,17.5, 20.1, 23, 27, 30.5, 34.4, 40.1, 45],
  atp = [1.75, 1.8, 1.95, 2.05, 2.45, 2.5, 2.6, 2.7, 3, 3.1, 3.4, 3.8, 3.8, 4]
)

fcf = DataFrame(
  incl = [5, 10, 20, 30, 40, 50],
  ft = [750, 1200, 2000, 2500, 3000, 3500],
  fcf = [0.18, 0.7, 0.9, 0.95, 1, 0.95]
)

name = "figure_08"

figure_08 = @pgf GroupPlot({
    group_style = {
      group_size = "1 by 2",
      vertical_sep = "0pt",
      xticklabels_at = "edge bottom"
    }
  },
  Axis({
      height = "6cm", width = " 8cm",
      xmin = 0, xmax = 55,
      ymin = 0, ymax = 5,
      grid = "major",
      xlabel = "Inclination [°]",
      ylabel = "TP [°/100ft]",
      title = "Figure 8: Actual vs. TP",
      legend_pos = "south east",
    },
    Plot({ color="blue", mark = "." },
      Table(tpf.incl, tpf.tp)),
    LegendEntry("Isotropic TP"),
    Plot({ color="green", mark = "." },
      Table(tpf.incl, tpf.tp - fcf.fcf)),
    LegendEntry("Predicted TP"),
    Plot({ thick, color="black", mark = "+", },
      Table(actual.incls, actual.atp)),
    LegendEntry("TP, actual")),
    LegendEntry("RPM"),
  Axis({
        height = "6cm", width = "8cm",
        ylabel = "Sideforce NBS",
        xmin = 0, xmax = 55,
        ymin = 0, ymax = 4000,
        grid = "major",
        xlabel = "Inclination [°]",
        ylabel = "Sideforce",
        legend_pos = "south east",
      },
      Plot({ color="red", mark = "*", title ="TP" },
        Table(fcf.incl, fcf.ft)),
      LegendEntry("sideforce"),
  )
)

savefig(joinpath(ProjDir, "plots", name), figure_08, ProjDir)


text:
function show(io::IO, geo::RectRegion)
    print(
		io,
		"The Rectilinear Region $(geo.regID) has the following properties:\n",
		"    Region ID (regID) : ", geo.regID, '\n',
		"    Parent ID (parID) : ", geo.parID, '\n',
		"    Name      (name)  : ", geo.name, '\n',
		"    Bounds  (N,S,E,W) : ",[geo.N,geo.S,geo.E,geo.W], '\n',
        "        (is180,is360) : ",(geo.is180,geo.is360),"\n",
	)
end

function show(io::IO, geo::PolyRegion)
    print(
		io,
		"The Polygonal Region $(geo.regID) has the following properties:\n",
		"    Region ID (regID) : ", geo.regID, '\n',
		"    Parent ID (parID) : ", geo.parID, '\n',
		"    Name       (name) : ", geo.name,  '\n',
		"    Bounds  (N,S,E,W) : ",[geo.N,geo.S,geo.E,geo.W], '\n',
		"    Shape     (shape) : ", geo.shape, '\n',
        "        (is180,is360) : ",(geo.is180,geo.is360),"\n",
	)
end

function show(io::IO, grd::PolyGrid)
	nlon = length(grd.ilon)
	nlat = length(grd.ilat)
	mask = grd.mask
    print(
		io,
		"The Polygonal Grid has the following properties:\n",
		"    Grid Bounds       (igrid) : ", grd.igrid, '\n',
		"    Longitude Indices  (ilon) : ", grd.ilon,  '\n',
		"    Latitude Indices   (ilat) : ", grd.ilat,  '\n',
		"    Longitude Points   (glon) : ", grd.glon,  '\n',
		"    Latitude Points    (glat) : ", grd.glat,  '\n',
		"    Region Size (nlon * nlat) : $(nlon) lon points x $(nlat) lat points\n",
		"    Region Mask (sum(mask) / (nlon * nlat)) : $(sum(isone.(mask))) / $(nlon*nlat)\n"
	)
end

function show(io::IO, grd::RectGrid)
	nlon = length(grd.ilon)
	nlat = length(grd.ilat)
    print(
		io,
		"The Rectilinear Grid has the following properties:\n",
		"    Grid Bounds       (igrid) : ", grd.igrid, '\n',
		"    Longitude Indices  (ilon) : ", grd.ilon,  '\n',
		"    Latitude Indices   (ilat) : ", grd.ilat,  '\n',
		"    Longitude Points   (glon) : ", grd.glon,  '\n',
		"    Latitude Points    (glat) : ", grd.glat,  '\n',
		"    Region Size (nlon * nlat) : $(nlon) lon points x $(nlat) lat points\n",
	)
end


text:
using Mapper2
using Test

using Logging
disable_logging(Logging.BelowMinLevel)

# Add path to example architecture to the LOAD_PATH variable
for i in 1:2
    include(joinpath(Mapper2.PKGDIR, "example", "ex$i", "Example$i.jl"))
end
include(joinpath(Mapper2.PKGDIR, "example", "chessboard", "Chessboard.jl"))
using .Chessboard
using .Example1
using .Example2

# Taskgraph
include("Taskgraph.jl")

# Architecture
include("Architecture/Architecture.jl")

# Placement
include("Place/MoveGenerators.jl")
include("Place/InitialPlacement.jl")

# Maps
include("Map/IntegrationTest.jl")
include("Map/SaveLoad.jl")

# Long running placement
include("Place/Placement.jl")


text:
"""
    add_stream(conn::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}

Add and return a listener for a given set of calls to KRPC. The returned listener will produce a combined
value for *any* update of the streamed calls.

# Examples

```
using KRPC.Interface.SpaceCenter
stream = add_stream(conn, (KRPC.Interface.SpaceCenter.get_ActiveVessel(), KRPC.Interface.SpaceCenter.get_GameMode()))
for (vessel, game_mode) in stream
    println("The current vessel is \$vessel in game mode \$game_mode")
end
```
"""
function add_stream(conn::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}
    request_map = Dict{UInt64, Pair{Request, Int}}()
    handles = kerbal(conn, ((AddStream_Phantom.(calls, true))..., ))
    for ((index, call), handle) in zip(enumerate(calls), handles)
        request_map[handle.id] = call=>index
    end
    default_value = Any[Nothing() for i=1:K]
    out = Channel{Vector}(0)
    bind(out, conn.str_listener)
    listener = Listener(request_map, conn, default_value, out)
    setindex!.((conn.listeners, ), (listener, ), getproperty.(handles, :id))
    return listener
end

"""
    add_stream(f::Function, conn::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}

# Examples

Do-notation version of add_stream.
```
add_stream(conn, (KRPC.Interface.SpaceCenter.get_ActiveVessel(), KRPC.Interface.SpaceCenter.get_GameMode())) do stream
    nmessages = 0
    for (vessel, game_mode) in stream
        println("The current vessel is \$vessel in game mode \$game_mode")
    end
end
```
"""
function add_stream(f::Function, conn::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}
    listener = add_stream(conn, calls)
    try
        f(listener)
    finally
        close(listener)
    end
end

"""
    next_value(channel::KRPC.Listener)

Returns the next value is sent from the server for any of the RPC calls. `next_value` will block
until a value is produced.
"""
function next_value(channel::Listener)
    return take!(channel.channel)
end

Base.iterate(channel::Listener) = iterate(channel.channel)
Base.iterate(channel::Listener, state) = iterate(channel.channel, state)

function Base.close(channel::Listener)
    req = Request[]
    for id in keys(channel.streams)
        if !(id in keys(channel.streams))
            error("Attempted to remove unbound stream. Check if the stream has already been removed.")
        end
        delete!(channel.connection.listeners, id)
        push!(req, RemoveStream_Phantom(id))
    end
    kerbal(channel.connection, (req..., ))
end

function update_value(listener::Listener, stream::UInt64, value)
    req, arg_index = listener.streams[stream]
    processed_value = process_return(listener.connection, value, req)
    listener.current_value[arg_index] = processed_value
end

function send_current_values(listener::Listener)
    try_put!(listener.channel, listener.current_value)
end

# writes to c only if there's a task waiting on it
# we need this because KRPC will just yammer away on the
# socket and then crash if there's nothing listening on the
# other end, so can't wait for there to be a listener on ours
# assumes c is unbuffered
function try_put!(c::Channel, value)
    lock(c)
    taker = try
        Base.check_channel_state(c)
        notify(c.cond_wait)
        if isempty(c.cond_take.waitq)
            nothing
        else
            popfirst!(c.cond_take.waitq)
        end
    finally
        unlock(c)
    end
    if !isnothing(taker)
        schedule(taker, value)
        yield()
    end
    return value
end

function stream_listener(conn::KRPCConnection)
    try
        while (!isready(conn.active) || take!(conn.active))
            updated = Set{Listener}()
            data = readproto(RecvRawProto(conn.stream_conn), krpc.schema.StreamUpdate())
            for result in data.results
                if !(result.id in keys(conn.listeners))
                    continue
                end
                listener = conn.listeners[result.id]
                update_value(listener, result.id, result.result.value)
                push!(updated, listener)
            end
            send_current_values.(updated)
        end
    finally
        Base.close(conn.stream_conn) # we're about to die, so close the channel
    end
end


text:
using Turing, TuringBenchmarks.TuringTools
using LinearAlgebra

include("lr_helper.jl")

x, y = readlrdata()

@model hlr_nuts(x, y, θ) = begin

    N,D = size(x)

    σ² ~ Exponential(θ)
    α ~ Normal(0, sqrt(σ²))
    β ~ MvNormal(zeros(D), ones(D)*sqrt(σ²))

    for n = 1:N
        y[n] ~ BinomialLogit(1, dot(x[n,:], β) + α)
    end
end

# Sampling parameter settings
n_samples = 10_000

# Sampling
@tbenchmark chain = sample(hlr_nuts(x, y, 1/0.1), NUTS(0.65), n_samples);


text:
# Good starting points for OSM overpass are:
# https://wiki.openstreetmap.org/wiki/Overpass_API/Language_Guide
# https://wiki.openstreetmap.org/wiki/Map_Features

"""
Generate a Overpass QL query. See:
https://wiki.openstreetmap.org/wiki/Overpass_API/Overpass_QL

# Arguments
- `latitudenorth::Float64`: The northern boundary of the region to query.
- `longitudeeast::Float64`: The eastern boundary of the region to query.
- `latitudesouth::Float64`: The southern boundary of the region to query.
- `longitudewest::Float64`: The western boundary of the region to query.
- `queryname::String = "cardrive"`: Specifeis the query content, so far only "cardrive" implemented.
- `timeout::Integer = 600`: Overpass timeout.
- `output::String = "json"`: Overpass output format, so far only "json" implemented.
"""
function overpassquery(latitudesouth::Float64, longitudewest::Float64, latitudenorth::Float64, longitudeeast::Float64;
                       queryname::String="cardrive", timeout::Integer=600, output::String = "json")::String

    if queryname === "cardrive" # a overpass query optimized for car routing
        query = """
        [out:$output]
        [timeout:$timeout]
        [bbox:$latitudesouth,$longitudewest,$latitudenorth,$longitudeeast];
        (
            way["highway"]
                ["area"!~"yes"]
                ["highway"!~"cycleway|footway|path|pedestrian|steps|track|corridor|
                    proposed|construction|bridleway|abandoned|platform|raceway|service"]
                ["motor_vehicle"!~"no|privat|delivery|destination|customers|agricultural|forestry"]
                ["motorcar"!~"no|privat|designated|delivery|destination|customers|agricultural|forestry"]
                ["vehicle"!~"no|privat|designated|delivery|destination|customers|agricultural|forestry"]
                ["access"!~"no|privat|designated|delivery|destination|customers|agricultural|forestry"];
            way["highway"]
                ["service"~"parking|parking_aisle"];
        )->.streets_network;
        (
            rel(bw)["restriction"];
            rel(bw)["restriction:motor_vehicle"];
            rel(bw)["restriction:motorcar"];
            rel(bw)["restriction:vehicle"];
            rel(bw)["restriction:conditional"];
            rel(bw)["restriction:motor_vehicle:conditional"];
            rel(bw)["restriction:motorcar:conditional"];
            rel(bw)["restriction:vehicle:conditional"];
        )->.turn_restrictions;
        ((.streets_network; .turn_restrictions;);>;);
        out body;
        """
    else
        query = ""
        @warn "Not yet implemented"
    end

    return query
end

"""
Query OSM data using Overpass QL.

# Arguments
- `query::String`: A Overpass QL query. See also: [`overpassquery`](@ref).
"""
function overpass(query::String)::Dict
    if !occursin("[out:json]", query)
        error("Use '[out:json]' as output flag, other formats are not yet implemented")
    end
    str =  String(HTTP.get("https://overpass-api.de/api/interpreter", body=query).body)
    data = JSON.parse(str)

    return data
end


text:
module Simple

using ContextualDispatch
import ContextualDispatch: Context, overdub

foo(x, y) = begin
    q = 10 + 15
    x + y + 10 + q
end

mutable struct SimpleCtx <: Context
    call_num::Int
    SimpleCtx() = new(0)
end
function overdub(ctx::SimpleCtx, ::typeof(+), args...)
    ctx.call_num += 1
    ret = *(args...)
    ret
end

ContextualDispatch.@load_call()

r, c = call(SimpleCtx(), foo, 5, 10)
display((r, c))

end # module


text:
function detect_data_sources(dir::String = joinpath(Pkg.dir(), "OpenSecrets/data/"))
    data_sources = Dict()

    for thing in readdir(dir)
        thing_path = abspath(joinpath(dir, thing))

        !(isdir(thing_path) && ismatch(r"^CampaignFin[0-9]{2}$", thing)) && continue

        year = int(match(r"[0-9]{2}$", thing).match)
        data_sources[year] = abspath(joinpath(dir, thing))
    end

    isdefined(:_campaign_finance_data_sources) && merge!(_campaign_finance_data_sources, data_sources)
    data_sources
end


text:
# Define types first
struct substance
    name::String
    density::Function
    enthalpy::Function
    heat_capacity::Function
    viscosity::Function
    thermal_conductivity::Function
    alt_heat_capacity::Function
    visc_kin::Function
    exp_coeff::Function
    therm_diff::Function
    Pr::Function
end

function substance(name, x::Vector{Function})
    dens, enth, Cp, visc, thermCond = x
    function alt_heat_capacity(p, T)
        return derivative(T_ -> enth(p, T_), T)
    end
    function visc_kin(p, T)
        return visc(p, T)/dens(p, T)
    end
    function exp_coeff(p, T)
        dRdT = derivative(T_ -> dens(p, T_), T)
        return (-1/dens(p, T))*dRdT
    end
    function thermDiff(p, T)
        if Cp == nothing
            Cp_ = alt_heat_capacity(p, T)
        else
            Cp_ = Cp(p, T)
        end
        return thermCond(p, T) / (dens(p,T)*Cp_)
    end
    function Pr(p, T)
        if Cp == nothing
            Cp_ = alt_heat_capacity(p, T)
        else
            Cp_ = Cp(p, T)
        end
        return Cp_ * visc(p, T) / thermCond(p, T)
    end
    return substance(name, dens, enth, Cp, visc, thermCond, alt_heat_capacity, visc_kin, exp_coeff, thermDiff, Pr)
end

@enum hx_type conc_parallel=0 conc_counter=1 st_oneshell=2 st_twoshell=3 Cr_0=4

function ϵ_Cp(NTU, Cr)
    top = 1 - exp(-NTU *(1+Cr))
    return top/(1+Cr)
end

function ϵ_Cc(NTU, Cr)

    if Cr==1
        return NTU/(1+NTU)
    end

    top = 1 - exp(-NTU *(1+Cr))
    bottom = 1 - Cr*exp(-NTU *(1+Cr))
    return top/bottom

end

function ϵ_St₁(NTU, Cr)
    term = exp(-NTU * √(1+Cr^2))
    return 2*(1 + Cr + √(1+Cr^2) * (1+term)/(1-term))
end

function ϵ_Stₙ(NTU, Cr, n)
    term = ((1 - ϵ_St₁(NTU, Cr)*Cr) / (1 - ϵ_St₁(NTU, Cr)))^n
    return (term - 1)/(term - Cr)
end

function ϵ_Cr0(NTU, Cr=0)
    return 1 - exp(-NTU)
end

ϵ_functions = Dict(conc_parallel => ϵ_Cp, conc_counter => ϵ_Cc, st_oneshell => ϵ_St₁, st_twoshell => ϵ_Stₙ, Cr_0 => ϵ_Cr0)

function ϵ(NTU, Cr, hx::hx_type)
    return ϵ_functions[hx](NTU, Cr)
end


text:
include("sand_01.jl")
using .Sand_01
using Plots
fig = Plots.plot()

NN = 200 # no of cells
tmax = 10.0 # max time
dt = 1.0 # how often to plot
tp = 0.0#5000 # when to start plotting

r=Sand_01.ring(NN; U=0.5)
for t in 0.0:dt:tmax
    if t > tp
        t = r.t
        x0_ = 0.5 * (r.X_[1:end-1] + r.X_[2:end])
        Plots.plot!(x0_[1:end],r.Σ_[1:end],label="t: $t",legend=:bottomright)#,title="Σ")#;legend=false)
        display(fig)
    end
    Sand_01.tstep!(r,dt)
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
function add_project_toml() :: Int64
    touch("Project.toml")

    open("Project.toml", "w") do io_project
        open("../sample/Project.toml") do io_sample
            write(io_project, io_sample)
        end
    end
end

function compare_http_header(headers::Array, key::String, value::String) :: Nothing
    for item in headers
        if item[1]==key
            @test item[2]==value
        end
    end
end

function delete_project() :: Nothing
    cd("../..")
    rm("demo", recursive=true)
    Dance.Router.delete_routes!()
    nothing
end

function extract_html_body_content(html_body::Array{UInt8,1}) :: String
    return split(
        split(String(html_body), "
")[2],
        "
"
    )[1]
end

function get_webpack_bundle_id() :: String
    bundle_id::String = ""

    for line in eachline(open("../html/base.html", "r"))
        if occursin("")[1]
        end
    end

    cd("..")
    return bundle_id
end

function project_settings() :: Nothing
    cd("settings")

    touch("dev.jl")
    open("dev.jl", "w") do io
        write(io, ":dev = true \n")
    end

    open("Global.jl", "a+") do io
        write(io, "include(\"dev.jl\")")
    end

    cd("..")
end

function static_files_append_content() :: Int64
    open("src/js/main.js", "a") do io_js
        write(io_js, "\$('#js-button').click(function() {\n  console.log(\"Button clicked\");\n});")
    end

    open("src/css/main.css", "a") do io_css
        write(io_css, "#js-button {\n  color: red;\n}")
    end
end


text:
# %% imports

include("../../Modules/CMAES/cmaes.jl")
include("../../Modules/save_param.jl")
include("critic_sim.jl")

# %%

function optim()
    cfg = YAML.load(open("Julia/CartPole/Critic/cfg_cmaes.yml"))
    N = length(cfg["tunekeys"])
    hyp = String[]
    for k in keys(cfg["tunekeys"])
        push!(hyp, k)
    end
    best = Array{Float64}[]
    for i in 1:20
        push!(best,[])
    end
    best_fit = [1000.0 for i in 1:20]
    error = false

    # Create initial point (optional)
    # results = load_param("Julia/CartPole/LIF_base/results_LIF88888EI.jld")
    # genes_init = results["weights"][20]

    lambda = 4+floor(Int64,3*log(N))
    mu = floor(Int64,lambda/2)
    c = CMAES( N=N, μ=mu, λ=lambda , τ=sqrt(N), τ_c=N^2, τ_σ=sqrt(N))
    for i in 1:1000
        @time step!(c)
        bestind = argmin(c.F_λ)
        minfit = c.F_λ[bestind]
        println(i," ", minfit, " (current best = $(minimum(best_fit)))")
        if minfit < best_fit[1]
            best_fit[1] = minfit
            best[1] = copy(c.offspring[bestind])
            idx = sortperm(best_fit, rev = true)
            best_fit = best_fit[idx]
            best = best[idx]
        end

        if length(best[1]) > 0 && i%20 == 0
            save_param(best_fit,best,hyp)
        end
    end

    return best_fit, best
end

# %%

optim()

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using PiezoSerial
using Test


text:
# helper test functions to check for SDD properties holding

function validate(sdd::Sdd)
   for node in linearize(sdd) # linearize first so that callee can use bit field
      validate_node(node)
   end
    #TODO make one of these for structured decomposability
    @test isdecomposable(sdd)
    @test iscanonical(sdd, 5; verbose = true)
 end

function validate_node(n::Sdd)
   validate_node(GateType(n), n)
end

function validate_node(::⋁Gate, n::Sdd)
   size = num_children(n)
   primes = compile(mgr(n),false)
   for i = 1:size
      element = children(n)[i]
      # has alternating layers
      @test GateType(element) isa ⋀Gate
      for j = i+1:size
         other_element = children(n)[j]
         # is deterministic
         @test isfalse(prime(element) & prime(other_element))
         # is compressed
         @test sub(element) !== sub(other_element)
      end
      primes = primes | prime(element)
   end
   # is exhaustive
   @test istrue(primes)
   # cannot be trimmed to the sub
   @test size >= 1
   # cannot be trimmed to the prime
   if size == 2
      e1 = children(n)[1]
      e2 = children(n)[2]
      has_false_sub = (isfalse(sub(e1)) || isfalse(sub(e2)))
      has_true_sub = (istrue(sub(e1)) || istrue(sub(e2)))
      @test !(has_false_sub && has_true_sub)
   end
   @test NodeType(vtree(n)) isa Inner
end

function validate_node(::⋀Gate, n::Sdd)
   @test num_children(n) == 2
   @test !(GateType(prime(n)) isa ⋀Gate)
   @test !(GateType(sub(n)) isa ⋀Gate)
   # has no false prime
   @test !isfalse(prime(n))
   @test NodeType(vtree(n)) isa Inner
   @test GateType(prime(n)) isa ConstantGate || varsubset_left(prime(n), n)
   @test GateType(sub(n)) isa ConstantGate || varsubset_right(sub(n), n)
end

function validate_node(::LiteralGate, l::Sdd)
   @test variable(l) == first(variables(vtree(l)))
end

function validate_node(::ConstantGate, ::Sdd)
   # nothing to check?
end


text:
mutable struct SimpleBush{T<:Integer} <: AbstractBush{T}
    net::AbstractNetwork{T}
    src::T
    linkactivity::Vector{Bool}
    orderednodes::Union{Nothing, Vector{T}} # `nothing` when not updated
end

function SimpleBush{T}(net::AbstractNetwork{T}, src::U, costfn::CostFunction) where {T<:Integer, U<:Integer}
    nlinks = numlinks(net)
    linkactivity = fill(false, nlinks)
    _, parentvec = dijkstra(net, src, costfn(zeros(nlinks)).costs)

    for l in links(net)
        i,j,lidx = upn(l), dwn(l), idx(net,l)
        if parentvec[j] == i
            linkactivity[lidx] = true
        end
    end

    SimpleBush{T}(net, src, linkactivity, nothing)
end

src(b::SimpleBush) = b.src
net(b::SimpleBush) = b.net

for fn in (:numnodes, :numzones, :numsources, :numsinks, :sources, :sinks, :zones)
    @eval $fn(b::SimpleBush) = $fn(b.net)
end

for fn in (:issource, :issink, :throughflowallowed, :props)
    @eval $fn(b::SimpleBush{T}, i::U) where {T<:Integer, U<:Integer} = $fn(b.net, i)
end

function idx(b::SimpleBush, args...; warn=true, kwargs...)
    i = idx(b.net, args...; kwargs...)
    if !b.linkactivity[i]
        @warn "Link $i inactive in this bush!"
    end
    return i
end

numlinks(b::SimpleBush) = sum(b.linkactivity)
links(b::SimpleBush) = links(b.net)[b.linkactivity]

function outneighbors(b::SimpleBush, i::U) where {U<:Integer}
    adjlist, points = fstar(b.net)
    rv = adjlist[i]
    activeones = b.linkactivity[points[i]:(points[i+1]-1)]
    return rv[activeones]
end

function inneighbors(b::SimpleBush, j::U) where {U<:Integer}
    adjlist, points = bstar(b.net)
    rv = adjlist[j]
    activeones = b.linkactivity[idx.(Ref(b.net), adjlist[j], Ref(j))]
    return rv[activeones]
end

function outdegree(b::SimpleBush, i::U) where {U<:Integer}
    _,points = fstar(b.net)
    sum(b.linkactivity[points[i]:(points[i+1]-1)])
end

function indegree(b::SimpleBush, j::U) where {U<:Integer}
    adjlist,_ = bstar(b.net)
    return sum(b.linkactivity[idx.(Ref(b.net), adjlist[j], Ref(j))])
end

add_link!(b::SimpleBush, i::U) where {U<:Integer} = (b.orderednodes = nothing; b.linkactivity[i] = true;)
rem_link!(b::SimpleBush, j::U) where {U<:Integer} = (b.orderednodes = nothing; b.linkactivity[j] = false;)

add_link!(b::SimpleBush, l::AbstractLink) = add_link!(b, idx(b.net, l))
rem_link!(b::SimpleBush, l::AbstractLink) = rem_link!(b, idx(b.net, l))

add_link!(b::SimpleBush, i::U, j::U) where {U<:Integer} = add_link!(b, idx(b.net, i, j))
rem_link!(b::SimpleBush, i::U, j::U) where {U<:Integer} = rem_link!(b, idx(b.net, i, j))

add_link!(b::SimpleBush, args...; kwargs...) = error("Can only activate an existing link for a bush!")
rem_link!(b::SimpleBush, args...; kwargs...) = error("Can only deactivate an existing link for a bush!")

add_zone!(b::SimpleBush, args...; kwargs...) = error("Cannot add a zone to a bush!")
rem_zone!(b::SimpleBush, args...; kwargs...) = error("Cannot remove a zone from a bush!")

has_link(b::SimpleBush, i, j) = has_link(b.net, i, j) && b.linkactivity[idx(b.net, i, j)]

function orderednodes(b::SimpleBush) # NOTE: it updates bush, when b.orderednodes is not updated
    if b.orderednodes == nothing
        b.orderednodes = invperm(topologicalorder(b))
    end
    b.orderednodes
end


text:
# Usage:
#     julia benchmark/run_benchmarks.jl
using BenchmarkTools
using ImageContrastAdjustment
using ImageTransformations
using TestImages
using ImageCore

on_CI = haskey(ENV, "GITHUB_ACTIONS")

img = testimage("mandril_gray")
tst_sizes = on_CI ? (64, ) : (64, 128)
tst_types = (N0f8, Float32, Gray{N0f8}, Gray{Float32}, RGB{N0f8}, RGB{Float32})

const SUITE = BenchmarkGroup()

# TODO: add MidwayEqualization
alg_list = (( "LinearStretching", LinearStretching()),
            ( "AdaptiveEqualization", AdaptiveEqualization()),
            ( "ContrastStretching", ContrastStretching()),
            ( "Equalization", Equalization()),
            ( "GammaCorrection", GammaCorrection()),
            ( "Matching", Matching))

function add_algorithm_benchmark!(suite, img, alg_name, alg;
                                  tst_sizes, tst_types)
    haskey(suite, alg_name) || (suite[alg_name] = BenchmarkGroup())

    for T in tst_types
        haskey(suite[alg_name], T) || (suite[alg_name][T] = BenchmarkGroup())
        for sz in tst_sizes
            tst_img = imresize(T.(img), (sz, sz))
            if alg === Matching
                tst_alg = alg(targetimg = tst_img)
            else
                tst_alg = alg
            end

            suite[alg_name][T]["$sz×$sz"] = @benchmarkable adjust_histogram($tst_img, $tst_alg)
        end
    end
end

for (alg_name, alg) in alg_list
    add_algorithm_benchmark!(SUITE, img, alg_name, alg;
                             tst_sizes=tst_sizes,
                             tst_types=tst_types)
end


text:
module Testsanitise
    using Pkg
    Pkg.activate(".")
    using Test
    using st16SseqJuliaTools

    #include("AlignmentMatrix.jl")
    include("UCparsing.jl")
    include("HMMutils.jl")
    include("BLASTparsing.jl")
    include("TaxonKitWrap.jl")
    include("Utils.jl")

end


text:
using LightGraphs

include("tree_helpers.jl")

@enum Int64 S I R

mutable struct Person
    SIR
    is_infected::Bool
    Person(SIR) = new(SIR, false)
end

function infect(people, graph, infection_pt)
    if people[infection_pt].SIR == S
        people[infection_pt].is_infected = true
        people[infection_pt].SIR = I
        for neighbor in graph.fadjlist[infection_pt]
            infect(people, graph, neighbor)
        end
    else
        return
    end

end

function sir_test()
    people = [Person(S) for i = 1:6]

    # setting the 4th person to be recovered
    people[4].SIR = R

    # creating a graph to work with
    g = SimpleGraph(6)

    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 4)
    add_edge!(g, 3, 5)
    add_edge!(g, 5, 6)

    infect(people, g, 1)

    for person in people
        println("infection status: ", person.is_infected)
    end

    return g
end


text:
# ---
# title: 1418. Display Table of Food Orders in a Restaurant
# id: problem1418
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/display-table-of-food-orders-in-a-restaurant/description/>
# hidden: true
# ---
#
# Given the array `orders`, which represents the orders that customers have done
# in a restaurant. More specifically
# `orders[i]=[customerNamei,tableNumberi,foodItemi]` where `customerNamei` is
# the name of the customer, `tableNumberi` is the table customer sit at, and
# `foodItemi` is the item customer orders.
#
# _Return the restaurant 's " **display table** "_. The " **display table** " is
# a table whose row entries denote how many of each food item each table
# ordered. The first column is the table number and the remaining columns
# correspond to each food item in alphabetical order. The first row should be a
# header whose first column is "Table", followed by the names of the food items.
# Note that the customer names are not part of the table. Additionally, the rows
# should be sorted in numerically increasing order.
#
#
#
# **Example 1:**
#
#
#
#     Input: orders = [["David","3","Ceviche"],["Corina","10","Beef Burrito"],["David","3","Fried Chicken"],["Carla","5","Water"],["Carla","5","Ceviche"],["Rous","3","Ceviche"]]
#     Output: [["Table","Beef Burrito","Ceviche","Fried Chicken","Water"],["3","0","2","1","0"],["5","0","1","0","1"],["10","1","0","0","0"]]
#     Explanation: The displaying table looks like:
#     **Table,Beef Burrito,Ceviche,Fried Chicken,Water**
#     3    ,0           ,2      ,1            ,0
#     5    ,0           ,1      ,0            ,1
#     10   ,1           ,0      ,0            ,0
#     For the table 3: David orders "Ceviche" and "Fried Chicken", and Rous orders "Ceviche".
#     For the table 5: Carla orders "Water" and "Ceviche".
#     For the table 10: Corina orders "Beef Burrito".
#
#
# **Example 2:**
#
#
#
#     Input: orders = [["James","12","Fried Chicken"],["Ratesh","12","Fried Chicken"],["Amadeus","12","Fried Chicken"],["Adam","1","Canadian Waffles"],["Brianna","1","Canadian Waffles"]]
#     Output: [["Table","Canadian Waffles","Fried Chicken"],["1","2","0"],["12","0","3"]]
#     Explanation:
#     For the table 1: Adam and Brianna order "Canadian Waffles".
#     For the table 12: James, Ratesh and Amadeus order "Fried Chicken".
#
#
# **Example 3:**
#
#
#
#     Input: orders = [["Laura","2","Bean Burrito"],["Jhon","2","Beef Burrito"],["Melissa","2","Soda"]]
#     Output: [["Table","Bean Burrito","Beef Burrito","Soda"],["2","1","1","1"]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= orders.length <= 5 * 10^4`
#   * `orders[i].length == 3`
#   * `1 <= customerNamei.length, foodItemi.length <= 20`
#   * `customerNamei` and `foodItemi` consist of lowercase and uppercase English letters and the space character.
#   * `tableNumberi `is a valid integer between `1` and `500`.
#
#
## @lc code=start
using LeetCode

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


text:
# Create example data
x = randn(100_000)
y = randn(100_000)
# Draw an hexagonal plot in the bigger bottom-right region
subplot(3, 3, (5, 9))
hexbin(x, y, colorbar = false)
# Draw marginal histograms
subplot(3, 3, (2, 3))
histogram(x)
subplot(3, 3, (4, 7))
histogram(y, horizontal = true, xflip = true)
# Draw a shade plot in the smaller top-left region
subplot(3, 3, 1)
shade(x, y)


text:
# the SIR with age-structure example but where the rules are programatically generated.

using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq
using Random
using StatsBase: sample

using Distributions: Exponential, Geometric, cdf

using Plots
using LinearAlgebra
using ProgressBars

include("utils.jl")

# age structure ------------------------------------------------------------
@present ThAgeSIR(FreeSchema) begin
    (S,I,R,Agent,Age)::Ob
    s::Hom(S, Agent)
    i::Hom(I, Agent)
    r::Hom(R, Agent)
    age::Hom(Age, Agent)

    AgeValue::AttrType
    agevalue::Attr(Age, AgeValue)
end

@acset_type AgeSIR(ThAgeSIR, index = [:s, :i, :r, :age])

getage(st::AgeSIR, sir::Symbol) = st[vcat(incident(st, st[sir], :age)...), :agevalue]

N = 1000

pop_TW_actual = [1011137, 1041749, 975803, 1213008, 1502279, 1618075, 1621625, 1926961, 2032452, 1755391, 1806638, 1852580, 1684376, 1416638, 903130, 1454933]
pop_TW = N * (pop_TW_actual / sum(pop_TW_actual))
pop_TW = Int.(round.(pop_TW))
N_ages = Float64.(pop_TW)

contact_unnorm = [1.13133995414878	0.557101193757009	0.277403588569526	0.166283484217312	0.243637116538954	0.501429856885401	0.773456505042611	0.664734678057495	0.373030323264845	0.207403718009079	0.219768245286902	0.183713318566528	0.100598824431586	0.0723635363224441	0.0451987854057666	0.0264190015627899
0.522116359073187	5.13993582380529	0.999881750990576	0.256872404884024	0.152965888814777	0.440037571740546	0.81718457664922	0.901816339582795	0.776770796127256	0.324978689513501	0.203946781656598	0.170308924697168	0.114410112258284	0.0766043987718478	0.0377257124824009	0.0309979523804595
0.233720556183555	1.92985936809443	11.0692392115508	0.983950046684279	0.331933259964229	0.36113787699594	0.563341231085396	0.903032808933502	1.17790738474796	0.651636521136264	0.342842258726613	0.173034004951767	0.085317189948624	0.0792891389004607	0.0627981424580935	0.0524757101169723
0.117682987086876	0.345910173819922	3.30103123371644	10.0791939893023	1.50916450884044	0.791285884658043	0.640497694236322	0.824953404104752	1.08583985409548	1.09515652926247	0.535792916941561	0.197364288066501	0.0713874830715285	0.0516466911367817	0.0315463617042823	0.0232625811949096
0.184283846973784	0.180619839511097	0.291480297355471	2.42727014100382	4.36737726042213	2.04934997728767	1.39723781506946	1.20032801881538	1.01070904475343	1.23524796447696	0.800832127889001	0.39902019518293	0.0944429030797382	0.0416382822716035	0.0470826453808244	0.0394951197797532
0.56156192817194	0.342992537769556	0.200122810772951	0.916989232128664	2.3980129020946	4.30977975618343	2.22681401915761	1.59261033191491	1.30238635661337	1.05207634128304	1.00949814315475	0.536928086612407	0.168547286689823	0.0624501078266029	0.0317257423796732	0.0269431049949397
0.817794867047499	1.1441520682733	0.841196707142586	0.551555851206977	1.22865834990464	2.17367993548303	3.72979423098863	2.14929466728635	1.55306730764128	1.11703404262428	0.859302034849221	0.601599257124208	0.256942068421412	0.106203837431624	0.0555526815130332	0.0511357782774989
0.665898107120277	1.15314199919716	0.9740043186808	0.697233734859189	0.825555277703876	1.54165452016297	1.99243882433986	3.25956559702894	2.06558139725526	1.22850364764608	0.791670106907652	0.438605150114473	0.244794735730732	0.148777023280774	0.0869686510606626	0.0349673148163165
0.344192143343282	0.807995410469339	1.14415043227526	1.22438695056131	0.988030619038175	1.30424464223099	1.72028688625099	1.87814450829015	2.90470388452872	1.49758830410338	0.934324050254837	0.330600903622231	0.182537686520201	0.114580140863156	0.0813846406012246	0.0400784886603155
0.235791530653026	0.585194161111923	0.776189668181761	1.67596626847397	0.98479644076933	1.01459833535067	1.25482711892082	1.36074838781177	1.44439037412333	2.06915205205562	0.978401494457824	0.39999986541145	0.147110687316517	0.0835852658511084	0.0775170890988546	0.0761516911425131
0.247247471139332	0.643310594857866	1.05811260803268	1.31212557944165	1.09555967723574	1.35073741525286	1.16660463514035	1.06882688561953	1.37118566862	1.49192984016939	1.73933382380217	0.740193371293642	0.243542236335681	0.109266882350377	0.0798121160310869	0.0833033422190585
0.442684358021657	0.702087727411373	0.662177052932204	0.816859526244507	0.738982978496004	1.2358668685183	1.27567127460055	0.896236373088294	0.950326986928672	0.797183567958301	1.02749599715298	1.51681655407295	0.451707898396268	0.195099577202936	0.0961333966552319	0.0769572529007728
0.327638342195265	0.328285016941782	0.238165901349415	0.354821937305308	0.320966945158153	0.529029787102833	0.677350792373831	0.650949037046765	0.501900491233428	0.407240319726415	0.395920768113556	0.541100083728198	0.836769891969799	0.267847588746602	0.145638827229407	0.059691205799689
0.193568331575816	0.313138348784034	0.2640962452464	0.160094097408356	0.204553272242237	0.325505767060564	0.516407578574901	0.473341415992634	0.438186035505819	0.281001067157269	0.296302451590512	0.375946282497057	0.331749859023468	0.679792406234649	0.173000119847818	0.0702778253123571
0.0963863876721027	0.282393513167451	0.295599962015119	0.30895124164141	0.1384465024105	0.257953175229946	0.27963826844241	0.436956466193058	0.509042820432503	0.403187215044375	0.316140357668137	0.277393058839257	0.385660305905714	0.370945733661385	0.606083091791185	0.185944796718132
0.188903590501605	0.259894019899239	0.37803587788717	0.31873977192663	0.125868047138266	0.160125548829692	0.269195296463481	0.305611946209338	0.348244256184042	0.390220623327373	0.393237304756271	0.240097073151863	0.144076910885202	0.19701354194507	0.184345002762545	0.325239463933643];

contact = (contact_unnorm + transpose(contact_unnorm) .* ( pop_TW * transpose(1 ./ pop_TW))) ./ 2

C = contact
for i = 1:16
    for j = 1:16
        C[i,j] = contact[i,j]*N_ages[i]/N_ages[j]
    end
end

# parameters
N = 1000
I0 = 10
S0 = N - I0
Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5
β = R0 * γ / max(eigvals(C)...);

# rules are the dynamics of the model
struct rule_span
    l::ACSetTransformation # I->L
    r::ACSetTransformation # I->R
end

mutable struct rule_queue
    r::rule_span
    m::Vector{ACSetTransformation}
    fire::BitVector
    function rule_queue(l::ACSetTransformation, r::ACSetTransformation)
        return new(rule_span(l, r), Vector{ACSetTransformation}(), BitVector())
    end
end

function make_age_infection_prob(i,j, β, C, N_ages, Δt)
    Cij = C[i,j]
    Nj = N_ages[j]
    function(sr::rule_queue)
        p = cdf(Exponential(), β * Cij * (1 / Nj) * Δt)
        willfire = randsubseq(1:length(sr.fire), p)
        sr.fire[willfire] .= true
    end
end

function find_rule_matches(rules::Vector{rule_queue}, state::ACSet)
    for rule in rules
        rule.m = homomorphisms(codom(rule.r.l), state; monic = true)
        rule.fire = falses(length(rule.m))
    end
end

# apply rule matches
function fire_events(state::ACSet, events::Vector{rule_queue})
    newstate = state
    # for each event type
    for i in 1:length(events)
      if sum(events[i].fire) == 0
        continue
      end
      # for each match
      for j in 1:length(events[i].m)
          if events[i].fire[j] === false
              continue
          else
              # apply rewrite
              _, kg, _, kh = rewrite_match_maps(events[i].r.l, events[i].r.r, events[i].m[j])
              newstate = codom(kh)
              # update remaining matches post rewrite
              if j < length(events[i].m)
                  for k in j+1:length(events[i].m)
                  # if the event wasn't going to happen anyway, continue
                  if events[i].fire[k] === false
                      continue
                  else
                      m = postcompose_partial(kg, kh, events[i].m[k])
                      if isnothing(m)
                          events[i].fire[k] = false
                      else
                          events[i].m[k] = m
                      end
                  end
                  end
              end
              # update remaining matches for other event types post rewrite
              if i < length(events)
                  for k in i+1:length(events)
                      for l in 1:length(events[k].m)
                      # if the event wasn't going to happen anyway, continue
                      if events[k].fire[l] === false
                          continue
                      else
                          m = postcompose_partial(kg, kh, events[k].m[l])
                          if isnothing(m)
                              events[k].fire[l] = false
                          else
                              events[k].m[l] = m
                          end
                      end
                      end
                  end
              end
          end
      end
    end
    return newstate
end

# i=1
# j=2
# L = @acset AgeSIR{Int64} begin Agent=2; I=1; S=1; Age=2; agevalue=[i,j]; age=[1,2]; s=[1]; i=[2] end
# I = @acset AgeSIR{Int64} begin Agent=2; I=1; Age=2; agevalue=[i,j]; age=[1,2]; i=[2] end
# R = @acset AgeSIR{Int64} begin Agent=2; I=2; Age=2; agevalue=[i,j]; age=[1,2]; i=[1,2] end
# l = ACSetTransformation(I, L; Agent = [1,2], I = [1], Age = [1,2])
# r = ACSetTransformation(I, R; Agent = [1,2], I = [2], Age = [1,2])

# x = rule_span(l, r)
# x = rule_queue(l, r)

# infection rules
N_age = size(C)[1]

age_sir_rules = Vector{rule_queue}(undef, N_age^2 + 1)
age_sir_probability = Vector{Function}(undef, N_age^2 + 1)

for i in 1:16 # infectee
    for j in 1:16 # infector
        L = @acset AgeSIR{Int64} begin Agent=2; I=1; S=1; Age=2; agevalue=[i,j]; age=[1,2]; s=[1]; i=[2] end
        I = @acset AgeSIR{Int64} begin Agent=2; I=1; Age=2; agevalue=[i,j]; age=[1,2]; i=[2] end
        R = @acset AgeSIR{Int64} begin Agent=2; I=2; Age=2; agevalue=[i,j]; age=[1,2]; i=[1,2] end
        l = ACSetTransformation(I, L; Agent = [1,2], I = [1], Age = [1,2])
        r = ACSetTransformation(I, R; Agent = [1,2], I = [2], Age = [1,2])
        age_sir_rules[(i-1)*16+j] = rule_queue(l, r)
        age_sir_probability[(i-1)*16+j] = make_age_infection_prob(i, j, β, C, N_ages, Δt)
    end
end

# recovery rule
L = @acset AgeSIR{Int64} begin Agent=1; I=1; i=[1] end
I = @acset AgeSIR{Int64} begin Agent=1 end
R = @acset AgeSIR{Int64} begin Agent=1; R=1; r=[1] end

l = ACSetTransformation(I, L; Agent = [1])
r = ACSetTransformation(I, R; Agent = [1])

age_sir_rules[257] = rule_queue(l, r)
age_sir_probability[257] = function(sr::rule_queue)
    p = cdf(Exponential(), γ * Δt)
    willfire = randsubseq(1:length(sr.fire), p)
    sr.fire[willfire] .= true
end

# find homomorphisms into the current world state

# age_sir_rules[3].m = homomorphisms(codom(age_sir_rules[3].r.l), state; monic = true)
# age_sir_rules[3].fire = trues(length(age_sir_rules[3].m))

# new_state = fire_events(age_sir_rules)
# Debugger.@enter fire_events(age_sir_rules)

# find_rule_matches(age_sir_rules, state)

# for i in 1:length(age_sir_rules)
#     age_sir_probability[i](age_sir_rules[i])
# end

# make world state
ages = vcat([fill(i, pop_TW[i]) for i in 1:16]...)
shuffle!(ages)

sir_index = vcat(fill.(["S", "I", "R"],[S0, I0, N-S0-I0])...)
shuffle!(sir_index)

state = AgeSIR{Int64}()
add_parts!(state, :S, sum(sir_index .== "S"))
add_parts!(state, :I, sum(sir_index .== "I"))
add_parts!(state, :R, sum(sir_index .== "R"))

add_parts!(state, :Agent, N)
add_parts!(state, :Age, N)

set_subpart!(state, 1:nparts(state, :S), :s, findall(sir_index .== "S"))
set_subpart!(state, 1:nparts(state, :I), :i, findall(sir_index .== "I"))
set_subpart!(state, 1:nparts(state, :R), :r, findall(sir_index .== "R"))

set_subpart!(state, 1:N, :age, 1:N)
set_subpart!(state, 1:N, :agevalue, ages)

#  --------------------------------------------------------------------------------

# Simulation loop
out = fill(-1, (steps, 3))

for t = ProgressBar(1:steps)

    # enabled events
    find_rule_matches(age_sir_rules, state)

    # probability sample
    for i in 1:length(age_sir_rules)
        age_sir_probability[i](age_sir_rules[i])
    end

    global state = fire_events(state, age_sir_rules)

    # write output
    out[t, :] = [nparts(state, x) for x in [:S, :I, :R]]
end

# plots
plot(
    (1:steps) * Δt,
    out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

ever_infected_ages = [getage(state, :i); getage(state, :r)]
ever_infected_ages = [sum(ever_infected_ages .== i) for i = 1:16]

bar(pop_TW, label = false, xlabel = "Age bins", ylabel = "Number")
bar!(ever_infected_ages, label = false)


text:
"""
    NexPage(nexSerial, name; pid=pid, cid=cid)

A Nextion Page UI control.

A `NexPage` is the main element of a Nextion UI. It's where controls are put on.
"""
struct NexPage <: AbstractNexObject
    _nid::NexID

    function NexPage(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        new(NexID(nexSerial, name, pid, cid))
    end
end

"""
    show(page, Name)

Show Nextion page named `page` using its name.
"""
function show(page::NexPage, ::Type{Name})
    _nid = NexID(page)
    _name = String(Name(_nid))
    cmd = "page $_name"
    nexSerial = NexSerial(_nid)
    send(nexSerial, cmd)
end

"""
    show(page, PageID)

Show Nextion page named `page` using its page id.
"""
function show(page::NexPage, ::Type{PageID})
    _nid = NexID(page)
    _pid = String(PageID(_nid))
    cmd = "page $_pid"
    nexSerial = NexSerial(_nid)
    send(nexSerial, cmd)
end

"""
    show(page)

Show Nextion page named `page`.

It calls `show(page, Name)` so it's showing
page using its name.
"""
function show(page::NexPage)
    show(page, Name)
end

"""
    isshown(page) -> Bool

Return if a page is currently shown or not.
"""
function isshown(page::NexPage)::Bool
    _nid = NexID(page)
    pid1 = PageID(_nid._nexSerial)  #ToDo: should be implemented
    pid2 = PageID(NexID(page))
    pid1 == pid2
end


text:
@testset "SS-HOPM" begin

Random.seed!(1)
eigvec = randn(20)
eigval = norm(eigvec)
eigvec ./= eigval
T = eigval * eigvec .* eigvec' .* reshape(eigvec, 1, 1, 20)

@testset "Dense representation" begin
    @time (lbd, x) = sshopm(T, 1, verbose=true)
    @test norm(eigvec .- x) < 1e-5
    @test norm(eigval .- lbd) < 1e-5
end

@testset "Sparse representation" begin
    spT = SparseArray(T)

    @time (lbd, x) = sshopm(spT, 1, verbose=true)
    @test norm(eigvec .- x) < 1e-5
    @test norm(eigval .- lbd) < 1e-5
end

end


text:
using Blue
using Test


text: println(PROGRAM_FILE); for x in ARGS; println(x); end


text:
struct UniformScalingMap{T} <: LinearMap{T}
    λ::T
    M::Int
    function UniformScalingMap(λ::Number, M::Int)
        M < 0 && throw(ArgumentError("size of UniformScalingMap must be non-negative, got $M"))
        return new{typeof(λ)}(λ, M)
    end
end
UniformScalingMap(λ::Number, M::Int, N::Int) =
    (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
UniformScalingMap(λ::T, sz::Dims{2}) where {T} =
    (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))

MulStyle(::UniformScalingMap) = FiveArg()

# properties
Base.size(A::UniformScalingMap) = (A.M, A.M)
Base.isreal(A::UniformScalingMap) = isreal(A.λ)
LinearAlgebra.issymmetric(::UniformScalingMap) = true
LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A)
LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ)

# comparison of UniformScalingMap objects
Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (A.λ == B.λ && A.M == B.M)

# special transposition behavior
LinearAlgebra.transpose(A::UniformScalingMap) = A
LinearAlgebra.adjoint(A::UniformScalingMap)   = UniformScalingMap(conj(A.λ), size(A))

# multiplication with scalar
Base.:(*)(α::Number, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J))
Base.:(*)(J::UniformScalingMap, α::Number) = UniformScalingMap(J.λ * α, size(J))

# multiplication with vector
Base.:(*)(A::UniformScalingMap, x::AbstractVector) =
    length(x) == A.M ? A.λ * x : throw(DimensionMismatch("A_mul_B!"))

if VERSION < v"1.3.0-alpha.115"
function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector)
    (length(x) == length(y) == A.M || throw(DimensionMismatch("A_mul_B!")))
    if iszero(A.λ)
        return fill!(y, zero(eltype(y)))
    elseif isone(A.λ)
        return copyto!(y, x)
    else
        # call of LinearAlgebra.generic_mul! since order of arguments in mul! in
        # stdlib/LinearAlgebra/src/generic.jl reversed
        return LinearAlgebra.generic_mul!(y, A.λ, x)
    end
end
else # 5-arg mul! exists and order of arguments is corrected

A_mul_B!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector) = mul!(y, J, x)

end # VERSION

@inline function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false)
    @boundscheck (length(x) == length(y) == J.M || throw(DimensionMismatch("mul!")))
    _scaling!(y, J, x, α, β)
    return y
end

@inline function LinearAlgebra.mul!(Y::AbstractMatrix, J::UniformScalingMap, X::AbstractMatrix, α::Number=true, β::Number=false)
    @boundscheck size(X) == size(Y) || throw(DimensionMismatch("mul!"))
    @boundscheck size(X,1) == J.M || throw(DimensionMismatch("mul!"))
    _scaling!(Y, J, X, α, β)
    return Y
end

function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false)
    λ = J.λ
    if isone(α)
        if iszero(β)
            iszero(λ) && return fill!(y, zero(eltype(y)))
            isone(λ) && return copyto!(y, x)
            y .= λ .* x
            return y
        elseif isone(β)
            iszero(λ) && return y
            isone(λ) && return y .+= x
            y .+= λ .* x
            return y
        else # β != 0, 1
            iszero(λ) && (rmul!(y, β); return y)
            isone(λ) && (y .= y .* β .+ x; return y)
            y .= y .* β .+ λ .* x
            return y
        end
    elseif iszero(α)
        iszero(β) && (fill!(y, zero(eltype(y))); return y)
        isone(β) && return y
        # β != 0, 1
        rmul!(y, β)
        return y
    else # α != 0, 1
        if iszero(β)
            iszero(λ) && return fill!(y, zero(eltype(y)))
            isone(λ) && return y .= x .* α
            y .= λ .* x .* α
            return y
        elseif isone(β)
            iszero(λ) && return y
            isone(λ) && return y .+= x .* α
            y .+= λ .* x .* α
            return y
        else # β != 0, 1
            iszero(λ) && (rmul!(y, β); return y)
            isone(λ) && (y .= y .* β .+ x .* α; return y)
            y .= y .* β .+ λ .* x .* α
            return y
        end
    end
end

At_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
Ac_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

# combine LinearMap and UniformScaling objects in linear combinations
Base.:(+)(A₁::LinearMap, A₂::UniformScaling) = A₁ + UniformScalingMap(A₂.λ, size(A₁, 1))
Base.:(+)(A₁::UniformScaling, A₂::LinearMap) = UniformScalingMap(A₁.λ, size(A₂, 1)) + A₂
Base.:(-)(A₁::LinearMap, A₂::UniformScaling) = A₁ - UniformScalingMap(A₂.λ, size(A₁, 1))
Base.:(-)(A₁::UniformScaling, A₂::LinearMap) = UniformScalingMap(A₁.λ, size(A₂, 1)) - A₂


text:
using LinearAlgebra
using MLDatasets.MNIST
using ImageCore
using Plots
using Random
using PyCall
@pyimport pickle

function h(v)
    map(x -> 1.0 / (1.0 + exp(-x)), v)
end

function σ(a)
    c = maximum(a)
    exp_a = map(x -> exp(x .- c), a)
    exp_a / sum(exp_a)
end

function init_network()::Matrix{Float64}
    rand(Float64, (2,3))
end

function loss(
    x::Matrix{Float64},
    t::Matrix{Float64},
    W::Matrix{Float64},
    )::Float64
    z = x * W
    y = σ(z)
    cross_entropy_error(y, t)
end

function forward(network, x)
    W1, W2, W3 = network["W1"], network["W2"], network["W3"]
    b1, b2, b3 = network["b1"], network["b2"], network["b3"]

    a1 = x * W1 + b1
    z1 = h(a1)
    a2 = z1 * W2 + b2
    z2 = h(a2)
    a3 = z2 * W3 + b3
end

function prepare_testdata(;batch=1)
    data, labels = MNIST.testdata()
    batch_num = div(length(labels), batch)
    data0 = eachslice(reshape(data, (28 * 28 * batch, batch_num)), dims=2)
    labels0 = eachslice(transpose(reshape(labels, (batch, batch_num))), dims=1)
    collect(zip(data0, labels0))
end

function cross_entropy_error(
    y::AbstractArray{T, N},
    t::AbstractArray{T, N},
    )::AbstractFloat where{T <: Number, N}
    - sum(float.(t) .* log.(float.(y)))
end

function numerical_gradient(
    f::Function,
    xs::AbstractArray{T, N};
    h=1e-4::T,
    )::AbstractArray{T, N} where{T <: Number, N}
    grad = zeros(size(xs))
    for (i, x) in enumerate(xs)
        xs[i] = x + h
        fxh1 = f(xs)
        xs[i] = x - h
        fxh2 = f(xs)
        grad[i] = (fxh1 - fxh2) / (2h)
        xs[i] = x
    end
    grad
end

function gradient_descent(
    f::Function,
    init_x::AbstractArray{T, N};
    learning_rate=0.01::T,
    steps=100::UInt128
    )::AbstractArray{T, N} where{T <: Number, N}
    x = init_x
    for _ = 1:steps
        x -= learning_rate * numerical_gradient(f, x)
    end
    x
end

function make_predict(fname="./sample_weight.pkl")
    f = pybuiltin("open")(fname, "rb")
    network = pickle.load(f)
    f.close()
    W₁ ,W₂, W₃ = network["W1"], network["W2"], network["W3"]
    b₁ ,b₂, b₃ = network["b1"], network["b2"], network["b3"]
    function (x::AbstractArray)
        a₁ = x * W₁ .+ transpose(b₁)
        z₁ = h(a₁)
        a₂ = z₁ * W₂ .+ transpose(b₂)
        z₂ = h(a₂)
        a₃ = z₂ * W₃ .+ transpose(b₃)
        σ(a₃)
    end
end

function bench01()
    predict = make_predict()
    data = prepare_testdata()
    count = 0
    for (x, l) in data
        if argmax(predict(x))[2] - 1 == l[1]
            count += 1
        end
    end
    count / length(data)
end
bench01()

function bench01_batch(;batch=100)
    predict = make_predict()
    data = prepare_testdata(batch=batch)
    count = 0
    for (xs, ls) in data
        for (p, l) in zip(eachslice(predict(transpose(reshape(xs, (28 * 28, batch)))), dims=1), ls)
            if argmax(p) - 1 == l
                count += 1
            end
        end
    end
    count / (length(data) * batch)
end

max_stars_repo_name: kirisaki/deep-julia

max_issues_repo_name: kirisaki/deep-julia

max_forks_repo_name: kirisaki/deep-julia


text:
__precompile__()

module MendelImpute

    using LinearAlgebra
    using GeneticVariation
    using VCFTools
    using SnpArrays
    using GroupSlices
    using Random
    using ProgressMeter
    using JLD2, FileIO, JLSO
    using Distances
    using LazyArrays

    # main functions that users are exposed to
    export phase               # main function for imputation and phasing
    export compress_haplotypes # for compressing haplotype panels
    export paint, composition, unique_populations # chromosome painting and admix proportions
    export convert_compressed  # for ultra-compression outputs

    include("compress.jl")
    include("data_structures.jl")
    include("phasing.jl")
    include("haplotype_pair.jl")
    include("haplotype_rescreen.jl")
    include("haplotype_thinning.jl")
    include("haplotype_stepscreen.jl")
    include("dynamic_programming.jl")
    include("impute.jl")
    include("breakpoints.jl")
    include("intersect.jl")
    include("painting.jl")
    include("ultra_compress.jl")

end # module


text:
# strings can be converted using [float](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.float) and [int](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.int):
e_str1 = "2.718"
e = float(e_str1)
println(5e)
#> 13.5914
num_15 = int("15")
println(3num_15)
#> 45

# numbers can be converted to strings and formatted using [printf](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.@printf)
@printf "e = %0.2f\n" e
#> 2.718
# or to create another string [sprintf](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.@sprintf)
e_str2 = @sprintf("%0.3f", e)

# to show that the 2 strings are the same
println("e_str1 == e_str2: $(e_str1 == e_str2)")
#> e_str1 == e_str2: true

# available number format characters are [f, e, g, c, s, p, d](https://github.com/JuliaLang/julia/blob/master/base/printf.jl#L15):
# (pi is a predefined constant; however, since its type is
# "MathConst" it has to be converted to a float to be formatted)
@printf "fix trailing precision: %0.3f\n" float(pi)
#> fix trailing precision: 3.142
@printf "scientific form: %0.6e\n" 1000pi
#> scientific form: 3.141593e+03
# g is not implemented yet
@printf "a character: %c\n" 'α'
#> a character: α
@printf "a string: %s\n" "look I'm a string!"
#> a string: look I'm a string!
@printf "right justify a string: %50s\n" "width 50, text right justified!"
#> right justify a string:                    width 50, text right justified!
@printf "a pointer: %p\n" 1e10
#> a pointer: 0x00000002540be400
@printf "print a integer: %d\n" 1e10
#> print an integer: 10000000000

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
function loadYears(fileName::String)::Vector{Int}
    file = open(fileName)
    years = Vector{Int}()
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line, ". ")
        try
            year = parse(Int, fields[1])
            push!(years, year)
        catch
            println(line)
        end
    end
    close(file)
    years
end

function loadNonStopWords(fileName::String)::Tuple{Vector{String},Dict{String,Int}}
    file = open(fileName)
    nonStopWords = Vector{String}()
    nonStopWordIDs = Dict{String,Int}()
    for phrase in readlines(file)
        if phrase != ""
            push!(nonStopWords, phrase)
            nonStopWordIDs[phrase] = length(nonStopWords)
        end
    end
    close(file)
    (nonStopWords, nonStopWordIDs)
end

function loadDualAspectPairFreqs(fileName::String)::Vector{Tuple{Tuple{Int,Int},Int}}
    file = open(fileName)
    dualAspectPairFreqs = Vector{Tuple{Tuple{Int,Int},Int}}()
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line)
        @assert length(fields) == 3
        i = parse(Int, fields[1])
        j = parse(Int, fields[2])
        freq = parse(Int, fields[3])
        push!(dualAspectPairFreqs, ((i, j), freq))
    end
    close(file)
    dualAspectPairFreqs
end

function loadPhraseIDGroups(fileName::String, nonStopWordIDs::Dict{String,Int})::Vector{Vector{Int}}
    file = open(fileName)
    phraseIDGroups = Vector{Vector{Int}}()
    for line in readlines(file)
        if length(line) == 0
            continue
        end
        phrases = split(line)
        idGroup = Vector{Int}()
        for phrase in phrases
            if haskey(nonStopWordIDs, phrase)
                push!(idGroup, nonStopWordIDs[phrase])
            end
        end
        push!(phraseIDGroups, idGroup)
    end
    close(file)
    phraseIDGroups
end

function loadTopics(fileName::String, nonStopWordIDs::Dict{String,Int}, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}})::Tuple{Vector{Int},Vector{Set{Int}}}
    pairIDs = Dict{Tuple{Int,Int},Int}()
    for (pairID, pairFreq) in enumerate(dualAspectPairFreqs)
        pair = pairFreq[1]
        pairIDs[pair] = pairID
    end
    numPairs = length(dualAspectPairFreqs)
    topicIDs = Vector{Int}(undef, numPairs)
    topics = Vector{Set{Int}}()
    file = open(fileName)
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line)
        @assert length(fields) == 4
        topicID = parse(Int, fields[1])
        freq = parse(Int, fields[2])
        phrase1 = fields[3]
        phrase2 = fields[4]
        phraseID1 = nonStopWordIDs[phrase1]
        phraseID2 = nonStopWordIDs[phrase2]
        pairID = pairIDs[(phraseID1, phraseID2)]
        @assert freq == dualAspectPairFreqs[pairID][2]
        topicIDs[pairID] = topicID
        if topicID > length(topics)
            @assert topicID == length(topics) + 1
            push!(topics, Set{Int}())
        end
        push!(topics[topicID], pairID)
    end
    close(file)

    (topicIDs, topics)
end

function getPairYearDistributions(dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, phraseIDGroups::Vector{Vector{Int}}, years::Vector{Int})::Vector{Dict{Int,Int}}
    numPairs = length(dualAspectPairFreqs)
    pairIDs = Dict{Tuple{Int,Int},Int}()
    for (pairID, pairFreq) in enumerate(dualAspectPairFreqs)
        pair = pairFreq[1]
        pairIDs[pair] = pairID
    end

    pairYearDistributions = [Dict{Int,Int}() for i in 1:numPairs]
    for (groupID, phraseIDGroup) in enumerate(phraseIDGroups)
        year = years[groupID]
        groupSize = length(phraseIDGroup)
        groupPairIDs = Set{Int}()
        for i in 1:groupSize
            for j in 1:groupSize
                if i == j
                    continue
                end
                pair = (phraseIDGroup[i], phraseIDGroup[j])
                if haskey(pairIDs, pair)
                    pairID = pairIDs[pair]
                    push!(groupPairIDs, pairID)
                end
            end
        end
        for pairID in groupPairIDs
            pairYearDistributions[pairID][year] = get(pairYearDistributions[pairID], year, 0) + 1
        end
    end

    pairYearDistributions
end

function getTopicYearDistributions(topicIDs::Vector{Int}, topics::Vector{Set{Int}}, pairYearDistributions::Vector{Dict{Int,Int}})::Vector{Dict{Int,Int}}
    numTopics = length(topics)
    topicYearDistributions = [Dict{Int,Int}() for i in 1:numTopics]
    for (pairID, yearDistributions) in enumerate(pairYearDistributions)
        topicID = topicIDs[pairID]
        for (year, freq) in yearDistributions
            topicYearDistributions[topicID][year] = get(topicYearDistributions[topicID], year, 0) + freq
        end
    end
    topicYearDistributions
end

function getFiveYearSlopeSums(topicYearDistributions::Vector{Dict{Int,Int}}, toYear::Int)::Vector{Float64}
    numTopics = length(topicYearDistributions)
    slopeSums = Vector{Float64}(undef, numTopics)
    for (topicID, yearDistributions) in enumerate(topicYearDistributions)
        topicSlopeSum = 0.0
        for year in toYear - 4:toYear
            slope = 0.3 * (get(yearDistributions, year, 0) - get(yearDistributions, year - 3, 0)) + 0.1 *  (get(yearDistributions, year - 1, 0) - get(yearDistributions, year - 2, 0))
            topicSlopeSum += slope
        end
        slopeSums[topicID] = topicSlopeSum
    end
    slopeSums
end

function loadConns(fileName::String, n::Int)::Vector{Dict{Int,Float64}}
    sparseConns = [Dict{Int,Float64}() for i in 1:n]
    file = open(fileName)
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line)
        @assert length(fields) == 3
        i = parse(Int, fields[1])
        j = parse(Int, fields[2])
        conn = parse(Float64, fields[3])
        sparseConns[i][j] = conn
    end
    close(file)
    sparseConns
end

function findTop3Representations(topics::Vector{Set{Int}}, topicID::Int, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, sparseConns::Vector{Dict{Int,Float64}})::NTuple{3, Int}
    topic = topics[topicID]
    numPairsInTopic = length(topic)
    topicPairIDs = [pairID for pairID in topic]
    pairScores = Vector{Float64}(undef, numPairsInTopic)
    for (i, pairID) in enumerate(topicPairIDs)
        freq = dualAspectPairFreqs[pairID][2]
        pairConns = sparseConns[pairID]
        sumPairConns = 0.0
        for (neighbor, conn) in pairConns
            sumPairConns += conn
        end
        pairScores[i] = freq^2 + 0.5 * sumPairConns
    end
    pairRanks = sortperm(pairScores, rev=true)
    (topicPairIDs[pairRanks[1]], topicPairIDs[pairRanks[2]], topicPairIDs[pairRanks[3]])
end

function saveTopKTopics(fileName::String, K::Int, topics::Vector{Set{Int}}, topicSlopeSums::Vector{Float64}, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, sparseConns::Vector{Dict{Int,Float64}}, nonStopWords::Vector{String})
    file = open(fileName, "w")
    topicRanks = sortperm(topicSlopeSums, rev=true)
    for k = 1:K
        topicID = topicRanks[k]
        score = topicSlopeSums[topicID]
        pairID1, pairID2, pairID3 = findTop3Representations(topics, topicID, dualAspectPairFreqs, sparseConns)
        pair1 = dualAspectPairFreqs[pairID1][1]
        phrase11 = nonStopWords[pair1[1]]
        phrase12 = nonStopWords[pair1[2]]
        pair2 = dualAspectPairFreqs[pairID2][1]
        phrase21 = nonStopWords[pair2[1]]
        phrase22 = nonStopWords[pair2[2]]
        pair3 = dualAspectPairFreqs[pairID3][1]
        phrase31 = nonStopWords[pair3[1]]
        phrase32 = nonStopWords[pair3[2]]
        println(file, "$topicID $score $phrase11 $phrase12 $phrase21 $phrase22 $phrase31 $phrase32")
    end
    close(file)
end

years = loadYears("DualAspectTopicRepresentation/data/NIPS.txt")
nonStopWords, nonStopWordIDs = loadNonStopWords("DualAspectTopicRepresentation/workspace/NIPS.nonStopWords.txt")
dualAspectPairFreqs = loadDualAspectPairFreqs("DualAspectTopicRepresentation/workspace/NIPS.DAPairFreq.txt")
phraseIDGroups = loadPhraseIDGroups("DualAspectTopicRepresentation/workspace/NIPS.w2v.training.txt", nonStopWordIDs)
topicIDs, topics = loadTopics("DualAspectTopicRepresentation/workspace/NIPS.topics.txt", nonStopWordIDs, dualAspectPairFreqs)
pairYearDistributions = getPairYearDistributions(dualAspectPairFreqs, phraseIDGroups, years)
topicYearDistributions = getTopicYearDistributions(topicIDs, topics, pairYearDistributions)
topicSlopeSums = getFiveYearSlopeSums(topicYearDistributions, 2020)
sparseConns = loadConns("DualAspectTopicRepresentation/workspace/NIPS.DAPairConns.txt", length(dualAspectPairFreqs))
saveTopKTopics("DualAspectTopicRepresentation/workspace/NIPS.Top100Topics.txt", 100, topics, topicSlopeSums, dualAspectPairFreqs, sparseConns, nonStopWords)

max_stars_repo_path: src/06 rank_topics.jl

max_issues_repo_path: src/06 rank_topics.jl

max_forks_repo_path: src/06 rank_topics.jl


text:
# Julia wrapper for header: /usr/include/scip/pub_misc_linear.h
# Automatically generated using Clang.jl wrap_c

function SCIPconsGetRhs(scip, cons, success)
    ccall((:SCIPconsGetRhs, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{UInt32}), scip, cons, success)
end

function SCIPconsGetLhs(scip, cons, success)
    ccall((:SCIPconsGetLhs, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{UInt32}), scip, cons, success)
end

function SCIPgetConsVals(scip, cons, vals, varssize, success)
    ccall((:SCIPgetConsVals, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Cdouble}, Cint, Ptr{UInt32}), scip, cons, vals, varssize, success)
end

function SCIPconsGetDualfarkas(scip, cons, dualfarkas, success)
    ccall((:SCIPconsGetDualfarkas, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Cdouble}, Ptr{UInt32}), scip, cons, dualfarkas, success)
end

function SCIPconsGetDualsol(scip, cons, dualsol, success)
    ccall((:SCIPconsGetDualsol, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Cdouble}, Ptr{UInt32}), scip, cons, dualsol, success)
end

function SCIPconsGetRow(scip, cons)
    ccall((:SCIPconsGetRow, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end