text:
using Statistics

mutable struct Acao
    DATAPREGAO::Vector{String}
    CODNEG::String
    QUATOT::Vector{Int}
    PREABE::Vector{Float64}
    PREMAX::Vector{Float64}
    PREMIN::Vector{Float64}
    PREMED::Vector{Float64}
    PREULT::Vector{Float64}
    variancia::Float64
    media::Float64
    desvioPadrao::Float64
    mediaMovelSimplesCurtoPrazo::Vector{Float64}
    mediaMovelSimplesLongoPrazo::Vector{Float64}
    mediaMovelSimplesRetorno5::Vector{Float64}
end

function buildAcao(DATAPREGAO::String, CODNEG::String, QUATOT::Int, PREABE::Float64, PREMAX::Float64, PREMIN::Float64, PREMED::Float64, PREULT::Float64)
    return Acao([DATAPREGAO], CODNEG, [QUATOT], [PREABE], [PREMAX], [PREMIN], [PREMED], [PREULT], 0, 0, 0, Float64[], Float64[], Float64[])
end

function calcularRetorno(serie::Vector{Float64})::Vector{Float64}
    retorno::Vector{Float64} = Float64[]
    anterior = serie[1]

    for valor in serie[2:lastindex(serie)]
        append!(retorno, (valor - anterior) / anterior)
    end

    return retorno
end

function mediaMovelSimples(serieAnterior::Union{Vector{Float64},Nothing}, serie::Vector{Float64}, periodos::Int64)::Vector{Float64}
    _serie = serieAnterior !== nothing ? append!(serieAnterior, serie) : serie
    _mediaMovel::Vector{Float64} = Float64[]
    _ultimoIndice = lastindex(_serie)
    periodoInicial = 1
    periodoAtual = periodos

    while periodoAtual <= _ultimoIndice
        media = 0
        for valor in _serie[periodoInicial:periodoAtual]
            media += valor
        end
        append!(_mediaMovel, media / periodos)
        periodoAtual += 1
        periodoInicial += 1
    end

    return _mediaMovel
end

function mediaMovelExponencial(serieAnterior::Vector{Float64}, serie::Vector{Float64}, periodos::Int16)::Vector{Float64}
    _serie = length(serieAnterior) > 0 ? append!(serieAnterior, serie) : serie
    _mediaMovel::Vector{Float64} = Float64[]
    _ultimoIndice = lastindex(_serie)
    periodoInicial = 1
    periodoAtual = periodos

    while periodoAtual <= _ultimoIndice
        media = 0
        for valor in _serie[periodoInicial:periodoAtual]
            media += valor
        end
        append!(_mediaMovel, media / periodos)
        periodoAtual += 1
        periodoInicial += 1
    end

    return _mediaMovel
end
@time begin
    println("Portifolio de Ações")
    _dados = Float64[]

    arquivo_cotacoes = "/Users/rafaeldias/Downloads/COTAHIST_A2021 3.TXT"
    println("Carregar cotações do arquivo $arquivo_cotacoes")
    _arquivo = open(arquivo_cotacoes, "r")
    _dados = readlines(_arquivo)
    close(_arquivo)

    arquivo_cotacoes = "/Users/rafaeldias/Downloads/COTAHIST_A2022 5.TXT"
    println("Carregar cotações do arquivo $arquivo_cotacoes")
    _arquivo = open(arquivo_cotacoes, "r")
    append!(_dados, readlines(_arquivo))
    close(_arquivo)

    global acoes = Dict{String,Acao}()
    global carteiraAtual = Dict{String,Acao}()

    for _cotacao in _dados
        if length(_cotacao) > 0
            TIPREG = _cotacao[1:2]
            DATAPREGAO = _cotacao[3:10]
            CODNEG = replace(_cotacao[13:24], " " => "")
            TPMERC = _cotacao[25:27]
            QUATOT = _cotacao[153:170]
            PREABE = _cotacao[57:69]
            PREMAX = _cotacao[70:82]
            PREMIN = _cotacao[83:95]
            PREMED = _cotacao[96:108]
            PREULT = _cotacao[109:121]

            # Mercado Fracionario
            if TIPREG == "01" && TPMERC == "020"
                QUATOT = parse(Int, QUATOT)
                PREABE = parse(Int, PREABE) / 100.00 # DIVIDIR POR 100
                PREMAX = parse(Int, PREMAX) / 100.00 # DIVIDIR POR 100
                PREMIN = parse(Int, PREMIN) / 100.00 # DIVIDIR POR 100
                PREMED = parse(Int, PREMED) / 100.00 # DIVIDIR POR 100
                PREULT = parse(Int, PREULT) / 100.00 # DIVIDIR POR 100
                if haskey(acoes, CODNEG)
                    append!(acoes[CODNEG].QUATOT, QUATOT)
                    append!(acoes[CODNEG].PREABE, PREABE)
                    append!(acoes[CODNEG].PREMAX, PREMAX)
                    append!(acoes[CODNEG].PREMIN, PREMIN)
                    append!(acoes[CODNEG].PREMED, PREMED)
                    append!(acoes[CODNEG].PREULT, PREULT)
                else
                    acoes[CODNEG] = buildAcao(DATAPREGAO, CODNEG, QUATOT, PREABE, PREMAX, PREMIN, PREMED, PREULT)
                end
            end

            # Mercado a Vista
            if TIPREG == "01" && TPMERC == "010" && (startswith(CODNEG, "BOVA11"))
                QUATOT = parse(Int, QUATOT)
                PREABE = parse(Int, PREABE) / 100.00 # DIVIDIR POR 100
                PREMAX = parse(Int, PREMAX) / 100.00 # DIVIDIR POR 100
                PREMIN = parse(Int, PREMIN) / 100.00 # DIVIDIR POR 100
                PREMED = parse(Int, PREMED) / 100.00 # DIVIDIR POR 100
                PREULT = parse(Int, PREULT) / 100.00 # DIVIDIR POR 100

                if haskey(acoes, CODNEG)
                    append!(acoes[CODNEG].QUATOT, QUATOT)
                    append!(acoes[CODNEG].PREABE, PREABE)
                    append!(acoes[CODNEG].PREMAX, PREMAX)
                    append!(acoes[CODNEG].PREMIN, PREMIN)
                    append!(acoes[CODNEG].PREMED, PREMED)
                    append!(acoes[CODNEG].PREULT, PREULT)
                else
                    acoes[CODNEG] = buildAcao(DATAPREGAO, CODNEG, QUATOT, PREABE, PREMAX, PREMIN, PREMED, PREULT)
                end
            end
        end
    end
    _dados = nothing
    # Carteira atual
    # ITUB4 6 23.2
    # SBSP3 3 36.17
    # ABEV3 7 15.72
    # BOVA11 1 103.49
    # LUPA3F 1

    somaVariancia = 0.0
    alocacoes = Dict{String,Vector{Float64}}()
    alocacoesAux = Vector{Float64}()
    codigos = Vector{String}()
    for acao::Acao in sort!(collect(values(acoes)), by = x -> var(x.PREULT))
        if acao.CODNEG == "CBAV3F" || acao.CODNEG == "ITUB4F" || acao.CODNEG == "SBSP3F" || acao.CODNEG == "_ABEV3F" || acao.CODNEG == "_BOVA11" || acao.CODNEG == "_LUPA3F" || acao.CODNEG == "ENAT3F"
            global somaVariancia += var(acao.PREULT)
        end
    end
    f = open("cotacoes.json", "w")
    write(f, "[")
    for acao::Acao in sort!(collect(values(acoes)), by = x -> var(x.PREULT))
        acao.variancia = var(acao.PREULT)
        acao.desvioPadrao = std(acao.PREULT)
        acao.media = mean(acao.PREULT)
        acao.mediaMovelSimplesCurtoPrazo = mediaMovelSimples(nothing, acao.PREULT, 5)
        acao.mediaMovelSimplesLongoPrazo = mediaMovelSimples(nothing, acao.PREULT, 30)
        acao.mediaMovelSimplesRetorno5 = mediaMovelSimples(nothing, calcularRetorno(acao.PREULT), 5)

        mediaRetorno = mean(last(acao.mediaMovelSimplesRetorno5, 5))
        mediaVolume = mean(last(acao.QUATOT, 5))

        acimaMedia = 0
        abaixoMedia = 0

        for retorno in acao.mediaMovelSimplesRetorno5
            if retorno > mediaRetorno
                acimaMedia += 1
            else
                abaixoMedia += 1
            end
        end

        if acao.CODNEG == "CBAV3F" || acao.CODNEG == "ITUB4F" || acao.CODNEG == "SBSP3F" || acao.CODNEG == "_BOVA11" || acao.CODNEG == "ENAT3F" || acao.CODNEG == "_LUPA3F" || acao.CODNEG == "_ABEV3F"
            global alocacoes[acao.CODNEG] = [1. / (acao.variancia / somaVariancia)]
            push!(alocacoesAux, 1. / (acao.variancia / somaVariancia))
            push!(codigos, acao.CODNEG)
            write(f, replace("{\"fatorAlocacao\":$(1. / (acao.variancia / somaVariancia)),\"portifolio\":\"SIM\",\"codigo\":\"$(acao.CODNEG)\",\"variancia\":$(acao.variancia),\"desviopadrao\":$(acao.desvioPadrao),\"media\":$(acao.media),\"volume\":$(mean(acao.QUATOT)),\"mmscurto\":$(last(acao.mediaMovelSimplesCurtoPrazo, 5)),\"mmslongo\":$(last(acao.mediaMovelSimplesLongoPrazo, 5)),\"mmsretorno\":$(last(acao.mediaMovelSimplesRetorno5, 5)), \"mediaretorno\":$(mediaRetorno)},\n", "Float64[]" => "[]", "NaN" => "0"))
        end

        # if mediaRetorno > 0 && acimaMedia > abaixoMedia && last(acao.mediaMovelSimplesRetorno5) > mediaRetorno && mediaVolume > 1000
        #     global alocacoes[acao.CODNEG] = [1. / (acao.variancia / somaVariancia)]
        #     push!(alocacoesAux, 1. / (acao.variancia / somaVariancia))
        #     push!(codigos, acao.CODNEG)
        #     write(f, replace("{\"portifolio\":\"NAO\",\"codigo\":\"$(acao.CODNEG)\",\"variancia\":$(acao.variancia),\"desviopadrao\":$(acao.desvioPadrao),\"media\":$(acao.media),\"volume\":$(mean(acao.QUATOT)),\"mmscurto\":$(last(acao.mediaMovelSimplesCurtoPrazo, 5)),\"mmslongo\":$(last(acao.mediaMovelSimplesLongoPrazo, 5)),\"mmsretorno\":$(last(acao.mediaMovelSimplesRetorno5, 5)), \"mediaretorno\":$(mediaRetorno)},\n", "Float64[]" => "[]", "NaN" => "0"))
        # end
    end
    write(f, "{\"codigo\":\"AAAAAA\"}]")
    flush(f)
    close(f)
    acoes = nothing
end

b = 1600.
a = transpose(alocacoesAux) \ [b]
for i in 1:length(codigos)
    print("$(codigos[i]) $(a[i] * alocacoesAux[i]) -> $(100. * ((a[i] * alocacoesAux[i]) / b))% => $(b * (((a[i] * alocacoesAux[i]) / b)))\n")
end

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
import Chemfiles
using CellListMap
using FastPow
using StaticArrays
using Printf
using Base.Threads
using Parameters
using Statistics: mean
using LinearAlgebra: norm_sqr

#
# Simulation setup
#
@with_kw struct Params{V,N,T,M,UnitCellType}
    x0::V = getcoor("./ne10k_initial.pdb")
    temperature::T = 300.
    nsteps::Int = 2_000
    dt::T = 1.0 # fs
    ibath::Int = 10
    print_energy::Int = 50
    print_traj::Int = 100
    trajfile::String = "ne10k_traj.xyz"
    cutoff::T = 12.
    box::Box{UnitCellType,N,T,M} = Box([ 46.37, 46.37, 46.37 ], cutoff, lcell=2)
    # Parameters for Neon
    mass::T = 20.17900 # g/mol
    ε::T = 0.0441795 # kcal/mol
    σ::T = 2*1.64009 # Å
    kB::T = 0.001985875 # Boltzmann constant kcal / mol K
end

function potential_energy(d2,ε,σ,u)
    @fastpow u += ε*( σ^12/d2^6 - 2*σ^6/d2^3 )
    return u
end

function forces(x,y,i,j,d2,ε,σ,f)
    r = y - x
    @fastpow dudr = -12*ε*(σ^12/d2^7 - σ^6/d2^4)*r
    @inbounds f[i] = f[i] + dudr
    @inbounds f[j] = f[j] - dudr
    return f
end

# Kinetic energy and temperature
compute_kinetic(v::AbstractVector,m) = (m/2)*sum(x -> norm_sqr(x), v)
compute_temp(kinetic,kB,n) = 2*kinetic/(3*kB*n)
compute_temp(v::AbstractVector,m,kB) = 2*compute_kinetic(v,m)/(3*kB*length(v))

# Remove drift from velocities
function remove_drift!(v)
    vmean = mean(v)
    v .= v .- Ref(vmean)
end

# Function to print output data
function print_data(istep,x,params,cl,kinetic,trajfile;parallel=parallel)
    @unpack print_energy, print_traj, kB, box, ε, σ = params
    if istep%print_energy == 0
        u = map_pairwise!(
            (x,y,i,j,d2,output) -> potential_energy(d2,ε,σ,output),
            0., box, cl, parallel=parallel,
        )
        temp = compute_temp(kinetic,kB,length(x))
        @printf(
            "STEP = %8i U = %12.5f K = %12.5f TOT = %12.5f TEMP = %12.5f\n",
            istep, u, kinetic, u+kinetic, temp
        )
    end
    if istep%print_traj == 0 && istep > 0
        println(trajfile,length(x))
        println(trajfile," step = ", istep)
        for i in 1:length(x)
           @printf(trajfile,"Ne %12.5f %12.5f %12.5f\n", ntuple(j -> x[i][j], 3)...)
        end
    end
    return nothing
end

# Read coordinates from NAMD-DCD file
function getcoor(file)
    traj = redirect_stdout(() -> Chemfiles.Trajectory(file), devnull)
    frame = Chemfiles.read_step(traj,0)
    Chemfiles.close(traj)
    return copy(reinterpret(reshape,SVector{3,Float64},Chemfiles.positions(frame)))
end

#
# Simulation
#
function simulate(params::Params{V,N,T,UnitCellType}; parallel=true) where {V,N,T,UnitCellType}
    @unpack x0, temperature, nsteps, box, dt, ε, σ, mass, kB = params
    trajfile = open(params.trajfile,"w")

    # To use coordinates in Angstroms, dt must be in 10ps. Usually packages
    # use ps and nm internally (thus multiply coordinates by 10 and divide
    # the timestep given in fs by 1000)
    dt = dt/100

    # Initial arrays
    x = copy(x0)
    f = similar(x)
    flast = similar(x)

    # Initial velocities
    v = randn(eltype(x),size(x))
    remove_drift!(v)
    # Adjust average to desidred temperature
    t0 = compute_temp(v,mass,kB)
    @. v = v * sqrt(temperature/t0)

    # Build cell lists for the first time
    cl = CellList(x,box,parallel=parallel)

    # preallocate threaded output, since it contains the forces vector
    f .= Ref(zeros(eltype(f)))
    f_threaded = [ deepcopy(f) for _ in 1:nthreads() ]
    aux = CellListMap.AuxThreaded(cl)

    # Print data at initial point
    kinetic = compute_kinetic(v,mass)
    print_data(0,x,params,cl,kinetic,trajfile,parallel=parallel)

    # Simulate
    for istep in 1:nsteps

        # Update positions (velocity-verlet)
        @. x = x + v*dt + 0.5*(f/mass)*dt^2

        # Update cell lists
        cl = UpdateCellList!(x,box,cl,aux,parallel=parallel)

        # Reset forces
        flast .= f
        f .= Ref(zeros(eltype(f)))
        @threads for it in 1:nthreads()
            f_threaded[it] .= Ref(zeros(eltype(f)))
        end

        # Update forces
        map_pairwise!(
            (x,y,i,j,d2,output) -> forces(x,y,i,j,d2,ε,σ,output),
            f, box, cl, parallel=parallel,
            output_threaded=f_threaded
        )

        # Update velocities
        @. v = v + 0.5*((flast + f)/mass)*dt

        # Print data and output file
        kinetic = compute_kinetic(v,mass)
        print_data(istep,x,params,cl,kinetic,trajfile,parallel=parallel)

        # Isokinetic bath
        if istep%params.ibath == 0
            remove_drift!(v)
            temp = compute_temp(kinetic,kB,length(v))
            @. v = v * sqrt(temperature/temp)
        end

   end
   close(trajfile)

end

params = Params()
@time simulate(params)

ext: jl


text:
export Passive,HH

const power_exp=Base.exp(series(0.0,1.0,zeros(Float64,6)...))

exp(x::Float64)=polyval(power_exp,x)::Float64

#implicit_euler(x,dt,xtau,xinf)=x +(1.0-exp(-dt/xtau))*(-(xinf/xtau)/(-1.0/xtau)-x) slower?

function implicit_euler(x::Float64, dt::Float64,xtau::Float64,xinf::Float64)
    x = x + (1.0 - exp(dt*(( ( ( - 1.0 ) ) ) / xtau)))*(- ( ( ( xinf ) ) / xtau ) / ( ( ( ( - 1.0) ) ) / xtau ) - x)
end

speed_init(prop::Channel,l::Int64)=(:(begin; end))

#=
Passive
still needs range variables
adapted for Julia by PMT
=#

type Passive <: Channel
    nodevar::Array{ASCIIString,1}
    g::Float64
    e::Float64
end

Passive()=(myvars=Array(ASCIIString,0); Passive(myvars,.001,-70.0))

function speed_init(prop::Passive,l::Int64)
    :(begin
      p.p[$(1+l),i]=.001
      p.p[$(2+l),i]=-70.0
      end)
end

speed_con(myprop::Passive,l::Int64)=(:(begin; end))

function speed_cur(myprop::Passive,l::Int64)
    :(begin
      im[k]+= p.p[$(1+l),i]*(v[k]- p.p[$(2+l),i])
      end)
end

#=
HH
 This is the original Hodgkin-Huxley treatment for the set of sodium,
  potassium, and leakage channels found in the squid giant axon membrane.
  ("A quantitative description of membrane current and its application
  conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
 Membrane voltage is in absolute mV and has been reversed in polarity
  from the original HH convention and shifted to reflect a resting potential
  of -65 mV.
SW Jaslove  6 March, 1992
adapted for Julia by PMT
=#

const q10 = 3^((6.3 - 6.3)/10)

type HH <: Channel
    nodevar::Array{ASCIIString,1}
    gnabar::Float64
    gkbar::Float64
    gl::Float64
    el::Float64
    gna::Float64
    gk::Float64
    m::Float64
    h::Float64
    n::Float64 #9
    minf::Float64
    hinf::Float64
    ninf::Float64
    mtau::Float64
    htau::Float64
    ntau::Float64 #15
end

HH()=(myvars=["ena", "ek"];HH(myvars,.12,.036,.0003,-54.3,zeros(Float64,11)...))

function speed_cur(myprop::HH,l::Int64)

    :(begin
      @inbounds ina = p.p[$(5+l),i] * (v[k] - ena)
      @inbounds ik  = p.p[$(6+l),i] * (v[k] - ek)
      @inbounds il  = p.p[$(3+l),i] * (v[k] - p.p[$(4+l),i])

      @inbounds im[k]+=ina+ik+il
      end)
end

function speed_con(myprop::HH,l::Int64)

    :(begin

        #"m" sodium activation system
        @inbounds alpha = .1 * vtrap(-(v[k]+40.0),10.0)
        @inbounds beta =  4.0 * exp(-(v[k]+65.0)/18.0)
        mysum = alpha + beta
        @inbounds p.p[$(13+l),i] = 1.0/(q10*mysum)
        @inbounds p.p[$(10+l),i] = alpha/mysum

        #"h" sodium inactivation system
        @inbounds alpha = .07 * exp(-(v[k]+65.0)/20.0)
        @inbounds beta = 1.0 / (exp(-(v[k]+35.0)/10.0)+1)
        mysum = alpha + beta
        @inbounds p.p[$(14+l),i] = 1.0/(q10*mysum)
        @inbounds p.p[$(11+l),i] = alpha/mysum

        #"n" potassium activation system
        @inbounds alpha = .01*vtrap(-(v[k]+55.0),10.0)
        @inbounds beta = .125*exp(-(v[k]+65.0)/80.0)
        mysum = alpha + beta
        @inbounds p.p[$(15+l),i] = 1/(q10*mysum)
        @inbounds p.p[$(12+l),i] = alpha/mysum

        @inbounds p.p[$(7+l),i]=implicit_euler(p.p[$(7+l),i],dt,p.p[$(13+l),i],p.p[$(10+l),i])
        @inbounds p.p[$(8+l),i]=implicit_euler(p.p[$(8+l),i],dt,p.p[$(14+l),i],p.p[$(11+l),i])
        @inbounds p.p[$(9+l),i]=implicit_euler(p.p[$(9+l),i],dt,p.p[$(15+l),i],p.p[$(12+l),i])

        @inbounds p.p[$(5+l),i] = p.p[$(1+l),i] * p.p[$(7+l),i] * p.p[$(7+l),i] * p.p[$(7+l),i] * p.p[$(8+l),i]
        @inbounds p.p[$(6+l),i] = p.p[$(2+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i]

    end)
end

function speed_init(myprop::HH,l::Int64)

    :(begin

      @inbounds p.p[$(1+l),i]=.12
      @inbounds p.p[$(2+l),i]=.036
      @inbounds p.p[$(3+l),i]=.0003
      @inbounds p.p[$(4+l),i]=-54.3

        #"m" sodium activation system
        @inbounds alpha = .1 * vtrap(-(v[k]+40.0),10.0)
        @inbounds beta =  4.0 * exp(-(v[k]+65.0)/18.0)
        mysum = alpha + beta
        @inbounds p.p[$(13+l),i] = 1.0/(q10*mysum)
        @inbounds p.p[$(10+l),i] = alpha/mysum

        #"h" sodium inactivation system
        @inbounds alpha = .07 * exp(-(v[k]+65.0)/20.0)
        @inbounds beta = 1.0 / (exp(-(v[k]+35.0)/10.0)+1)
        mysum = alpha + beta
        @inbounds p.p[$(14+l),i] = 1.0/(q10*mysum)
        @inbounds p.p[$(11+l),i] = alpha/mysum

        #"n" potassium activation system
        @inbounds alpha = .01*vtrap(-(v[k]+55.0),10.0)
        @inbounds beta = .125*exp(-(v[k]+65.0)/80.0)
        mysum = alpha + beta
        @inbounds p.p[$(15+l),i] = 1/(q10*mysum)
        @inbounds p.p[$(12+l),i] = alpha/mysum

        @inbounds p.p[$(7+l),i]=p.p[$(10+l),i]
        @inbounds p.p[$(8+l),i]=p.p[$(11+l),i]
        @inbounds p.p[$(9+l),i]=p.p[$(12+l),i]

        @inbounds p.p[$(5+l),i] = p.p[$(1+l),i] * p.p[$(7+l),i] * p.p[$(7+l),i] * p.p[$(7+l),i] * p.p[$(8+l),i]
        @inbounds p.p[$(6+l),i] = p.p[$(2+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i] * p.p[$(9+l),i]
    end)
end

function vtrap(x::Float64,y::Float64)
    if abs(x/y) < 1e-6
        trap = y*(1 - x/y/2)
    else
        trap = x/(exp(x/y) - 1)
    end
    trap
end

#=
Ca_HVA

Reuveni, Friedman, Amitai, and Gutnick, J.Neurosci. 1993
=#

type Ca_HVA <: Channel
    nodevar::Array{ASCIIString,1}
    gCa_HVAbar::Float64
    gCA_HVA::Float64
    m::Float64
    h::Float64
    minf::Float64
    hinf::Float64
    ica::Float64
    mtau::Float64
    htau::Float64
end

Ca_HVA()=(myvars=["eca"];Ca_HVA(myvars,.00001,zeros(Float64,7)...))

function con_calc(prop::Ca_HVA,v::Float64,dt::Float64)

    rates(prop,v)

    prop.m=implicit_euler(prop.m,dt,prop.mtau,prop.minf)
    prop.h=implicit_euler(prop.h,dt,prop.htau,prop.hinf)

    prop.gCa_HVA = prop.gCa_HVAbar * prop.m * prop.m * prop.h

    nothing
end

function cur_calc(prop::Ca_HVA,vars::Dict{ASCIIString,Float64},v::Float64)

    prop.ica = prop.gCA_HVA*(v-vars["eca"])

end

function prop_init(prop::Ca_HVA,v::Float64)
    rates(prop,v)
    prop.m=prop.mInf
    prop.h=prop.hInf
    nothing
end

function rates(prop::Ca_HVA,v::Float64)
    if v == -27.0
        v+=.0001
    end

    mAlpha =  (0.055*(-27-v))/(exp((-27-v)/3.8) - 1)
    mBeta  =  (0.94*exp((-75-v)/17))

    prop.mInf = mAlpha/(mAlpha + mBeta)
    prop.mTau = 1/(mAlpha + mBeta)

    hAlpha =  (0.000457*exp((-13-v)/50))
    hBeta  =  (0.0065/(exp((-v-15)/28)+1))

    prop.hInf = hAlpha/(hAlpha + hBeta)
    prop.hTau = 1/(hAlpha + hBeta)
    nothing
end

#=
LVA ca channel

Note: mtau is an approximation from the plots
Avery and Johnston 1996, tau from Randall 1997
shifted by -10 mv to correct for junction potential
corrected rates using q10 = 2.3, target temperature 34, orginal 21
=#

type Ca_LVAst <: Channel
    nodevar::Array{ASCIIString,1}
    gCa_LVAstbar::Float64
    gCA_LVAst::Float64
    m::Float64
    h::Float64
    minf::Float64
    hinf::Float64
    ica::Float64
    mtau::Float64
    htau::Float64
    qt::Float64
end

Ca_LVAst()=(myvars=["eca"];Ca_LVAst(myvars,.00001,zeros(Float64,7)...,2.3^((34-21)/10)))

function con_calc(prop::Ca_LVAst,v::Float64,dt::Float64)

    rates(prop,v)

    prop.m=implicit_euler(prop.m,dt,prop.mtau,prop.minf)
    prop.h=implicit_euler(prop.h,dt,prop.htau,prop.hinf)

    prop.gCa_LVAst = prop.gCa_LVAstbar * prop.m * prop.m * prop.h

    nothing
end

function cur_calc(prop::Ca_LVAst,vars::Dict{ASCIIString,Float64},v::Float64)

    prop.ica = prop.gCA_LVAst*(v-vars["eca"])

end

function prop_init(prop::Ca_LVAst,v::Float64)
    rates(prop,v)
    prop.m=prop.mInf
    prop.h=prop.hInf
    nothing
end

function rates(prop::Ca_LVAst,v::Float64)

    prop.mInf = 1.0000/(1+ exp((v - -30.000)/-6))
    prop.mTau = (5.0000 + 20.0000/(1+exp((v - -25.000)/5)))/prop.qt
    prop.hInf = 1.0000/(1+ exp((v - -80.000)/6.4))
    prop.hTau = (20.0000 + 50.0000/(1+exp((v - -40.000)/7)))/prop.qt

    nothing
end

type Ih <: Channel
    nodevar::Array{ASCIIString,1}
    gIhbar::Float64
    ehcn::Float64
    gIh::Float64
    m::Float64
    minf::Float64
    ihcn::Float64
    mtau::Float64
end

Ih()=(myvars=Array(ASCIIString,0);Ih(myvars,.00001, -45, zeros(Float64, 5)...))

function con_calc(prop::Ih,v::Float64,dt::Float64)

    rates(prop,v)

    prop.m=implicit_euler(prop.m,dt,prop.mtau,prop.minf)

    prop.gIh = prop.gIhbar * prop.m

    nothing
end

function cur_calc(prop::Ih,vars::Dict{ASCIIString,Float64},v::Float64)

    prop.ica = prop.gIh*(v-prop.ehcn)

end

function prop_init(prop::Ih,v::Float64)
    rates(prop,v)
    prop.m=prop.mInf
    nothing
end

function rates(prop::Ih,v::Float64)

    if v == -154.9
        v=+ .0001
    else
    end

    mAlpha =  0.001*6.43*(v+154.9)/(exp((v+154.9)/11.9)-1)
    mBeta  =  0.001*193*exp(v/33.1)
    prop.mInf = mAlpha/(mAlpha + mBeta)
    prop.mTau = 1/(mAlpha + mBeta)

    nothing
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
import MultivariatePolynomials: AbstractVariable, similarvariable, @similarvariable

@testset "Variable" begin
    @testset "polyvar macro index set" begin
        Mod.@polyvar x y z
        Mod.@polyvar x[1:3] y z[1:2]
        @test length(x) == 3
        @test length(z) == 2
        @test x[1] > x[2] > x[3] > y > z[1] > z[2]
    end
    @testset "PolyVar" begin
        Mod.@polyvar x
        @test 1 != x
        @test x != 0
        @test copy(x) == x
        @test nvariables(x) == 1
        @test !isapproxzero(x)
        @test !iszero(x)
        @test zero(x) == 0
        @test iszero(zero(x))
        @test zero(x) isa AbstractPolynomial{Int}
        @inferred zero(x)

        typetests(x)
        @test (@inferred polynomial(x)) isa AbstractPolynomial{Int}
        @test (@inferred polynomial(x, Float64)) isa AbstractPolynomial{Float64}

        @test nterms(x) == 1
        @test @inferred(terms(x)) == [x]

        Mod.@polyvar y
        @test degree(x, x) == 1
        @test degree(x, y) == 0
        @test length(exponents(x)) == 1
        @test first(exponents(x)) == 1
        @test isconstant(x) == false

        @test divides(x, x) == true
        @test divides(x, y) == false

        @testset "Issue #82" begin
            for v in (x, typeof(x))
                for fun in (oneunit, one)
                    @test fun(v) == 1
                    @test 1 == fun(v)
                    #@test isone(fun(v)) # Enable when Julia v0.6 is dropped
                    @test fun(v) isa AbstractMonomial
                    @inferred fun(v)
                end
            end
        end
    end
    @testset "Create similar variable" begin
        Mod.@polyvar x y
        f = x^2 + y

        z = similarvariable(f, Val{:z})
        @test z isa AbstractVariable

        z = similarvariable(typeof(f), Val{:z})
        @test z isa AbstractVariable

        @inferred similarvariable(f, Val{:z})
        @inferred similarvariable(typeof(f), Val{:z})

        w = similarvariable(f, :w)
        @test w isa AbstractVariable

        @similarvariable f o
        @test o isa AbstractVariable

        m = @similarvariable f u
        @test m isa AbstractVariable
    end
end


text:
using .CholeskyExperiment: get_indexes,
                           process_decomposition,
                           process_instances,
                           mpirun

@testset "mpi.jl" begin
    manager = ExperimentManager(db="choleskytest")
    init_manager(manager)
    @testset "get_indexes" begin
        @test get_indexes(manager, :gen, 2) == [[1], [2]]
        @test get_indexes(manager, :gen, 4) == [[1], [2]]
        @test get_indexes(manager, :gen, 1) == [[1, 2]]
        # TODO: Add more test for `get_indexes`
    end
    @testset "mpirun" begin
        mpirun(:genoriginal; dec_type=:cholesky, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 2
        mpirun(:genoriginal; dec_type=:MD, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 4
        mpirun(:genoriginal; dec_type=:EG, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 6
        mpirun(:genoriginal; dec_type=:LEX_M, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 8
        mpirun(:genoriginal; dec_type=:MCS_M, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 10
        mpirun(:genoriginal; dec_type=:minLC, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 12
        mpirun(:genoriginal; dec_type=:maxLC, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 14
        mpirun(:genoriginal; dec_type=:minTrig, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 16
        mpirun(:genoriginal; dec_type=:maxTrig, db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 18
        mpirun(:gen; db_name="choleskytest", finalize=false)
        @test length(manager.decompositions) == 20
        mpirun(:merge; db_name="choleskytest", merge_type=:molzahn, finalize=false, unsolved=[true, false], original=[true, false], merged=[false])
        @test length(manager.decompositions) == 40
        mpirun(:merge; db_name="choleskytest", merge_type=:sliwak, finalize=false, unsolved=[true, false], original=[true, false], merged=[false])
        @test length(manager.decompositions) == 60
        mpirun(:solve, db_name="choleskytest", finalize=false, original=[true])
        for d in manager.decompositions
            if d["nb_added_edges"] == 0
                @test haskey(d, "features") && haskey(d["features"], "solver")
            else
                @test haskey(d, "features") && !haskey(d["features"], "solver")
            end
        end
        mpirun(:solve, db_name="choleskytest", finalize=false, original=[false])
        for d in manager.decompositions
            @test haskey(d, "features") && haskey(d["features"], "solver")
        end
        CholeskyExperiment.db.compute_cholesky_diff(manager.decompositions)
        mpirun(:combine; db_name="choleskytest", finalize=false, unsolved=[true, false], original=[true, false], merged=[true, false])
        @test length(manager.decompositions) > 60
        mpirun(:solve; db_name="choleskytest", unsolved=[true])
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
function _conv(h, x)
    m = length(h)
    n = length(x)
    zero_pad_x(i) = 1 <= i <= n ? x[i] : 0
    return [sum(h[j] * zero_pad_x(i - j + 1) for j in 1:m) for i in 1:m+n-1]
end

@testset "Conv (issue #364)" begin
    n = 3
    m = 11
    h = rand(m)
    x = rand(n)
    hvar = Variable(m)
    hvar.value = h
    @test evaluate(conv(hvar, x)) ≈ _conv(h, x)
end

@testset "`conj` (issue #416)" begin
    A = [1 1im; -1im 1]
    X = ComplexVariable(2, 2)
    p = minimize(real(tr(conj(X))), [X == A])
    solve!(
        p,
        MOI.OptimizerWithAttributes(
            SCS.Optimizer,
            "verbose" => 1,
            "eps_abs" => 1e-6,
        ),
    )
    @test evaluate(X) ≈ A atol = 1e-3
end

@testset "`logisticloss` (issue #458)" begin
    x = Variable()
    expr = logisticloss(x)
    @test expr isa Convex.AbstractExpr

    set_value!(x, 1.5)
    @test evaluate(expr) ≈ log(1 + exp(1.5))
end


text:
using DashBootstrapComponents, DashHtmlComponents

collapse = html_div([
    dbc_button(
        "Open collapse",
        id="collapse-button",
        className="mb-3",
        color="primary",
        n_clicks=0,
    ),
    dbc_collapse(
        dbc_card(dbc_cardbody("This content is hidden in the collapse")),
        id="collapse",
        is_open=false,
    ),
]);

callback!(
    app,
    Output("collapse", "is_open"),
    Input("collapse-button", "n_clicks"),
    State("collapse", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
isdefined(:AbstractTrees) && workspace()
using AbstractTrees
using Base.Test
import Base: ==

AbstractTrees.children(x::Array) = x
tree = Any[1,Any[2,3]]

AbstractTrees.print_tree(STDOUT, tree)
@test collect(Leaves(tree)) == [1,2,3]
@test collect(PostOrderDFS(tree)) == Any[1,2,3,Any[2,3],Any[1,Any[2,3]]]
@test collect(StatelessBFS(tree)) == Any[Any[1,Any[2,3]],1,Any[2,3],2,3]

tree2 = Any[Any[1,2],Any[3,4]]
@test collect(PreOrderDFS(tree2)) == Any[tree2,Any[1,2],1,2,Any[3,4],3,4]

"""
    A tree in which every node has 0 or 1 children
"""
immutable OneTree
    nodes::Vector{Int}
end
AbstractTrees.treekind(::Type{OneTree}) = AbstractTrees.IndexedTree()
AbstractTrees.siblinglinks(::Type{OneTree}) = AbstractTrees.StoredSiblings()
AbstractTrees.relative_state(t::OneTree, _, __::Int) = 1
Base.getindex(t::OneTree, idx) = t.nodes[idx]
AbstractTrees.childindices(tree::OneTree, node::Int) =
    (ret = (node == 0 || tree[node] == 0) ? () : (tree[node],))
AbstractTrees.children(tree::OneTree) = AbstractTrees.children(tree, tree)
AbstractTrees.rootstate(tree::OneTree) = 1
AbstractTrees.printnode(io::IO, t::OneTree) =
    AbstractTrees.printnode(io::IO, t[AbstractTrees.rootstate(t)])

ot = OneTree([2,3,4,0])
AbstractTrees.print_tree(STDOUT, ot)
@test collect(AbstractTrees.Leaves(ot)) == [0]
@test collect(AbstractTrees.PreOrderDFS(ot)) == [2,3,4,0]
@test collect(AbstractTrees.PostOrderDFS(ot)) == [0,4,3,2]

"""
    Stores an explicit parent for some other kind of tree
"""
immutable ParentTree{T}
    tree::T
    parents::Vector{Int}
end
AbstractTrees.treekind{T}(::Type{ParentTree{T}}) = AbstractTrees.treekind(T)
AbstractTrees.parentlinks{T}(::Type{ParentTree{T}}) = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks{T}(::Type{ParentTree{T}}) = AbstractTrees.siblinglinks(T)
AbstractTrees.relative_state(t::ParentTree, _, __::Int) =
    AbstractTrees.relative_state(t.tree, _, __)
Base.getindex(t::ParentTree, idx) = t.tree[idx]
AbstractTrees.childindices(tree::ParentTree, node::Int) = AbstractTrees.childindices(tree.tree, node)
AbstractTrees.children(tree::ParentTree) = AbstractTrees.children(tree, tree)
AbstractTrees.rootstate(tree::ParentTree) = AbstractTrees.rootstate(tree.tree)
AbstractTrees.parentind(tree::ParentTree, node::Int) = tree.parents[node]
AbstractTrees.printnode(io::IO, t::ParentTree) =
    AbstractTrees.printnode(io::IO, t[AbstractTrees.rootstate(t)])

pt = ParentTree(ot,[0,1,2,3])
AbstractTrees.print_tree(STDOUT, pt)
@test collect(AbstractTrees.Leaves(pt)) == [0]
@test collect(AbstractTrees.PreOrderDFS(pt)) == [2,3,4,0]
@test collect(AbstractTrees.PostOrderDFS(pt)) == [0,4,3,2]

#=
immutable IntTree
    num::Int
    children::Vector{IntTree}
end
==(x::IntTree,y::IntTree) = x.num == y.num && x.children == y.children

@test treemap(PostOrderDFS(tree)) do ind, x, children
    IntTree(isa(x,Int) ? x : mapreduce(x->x.num,+,0,children),
        isempty(children) ? IntTree[] : children)
end == IntTree(6,[IntTree(1,IntTree[]),IntTree(5,[IntTree(2,IntTree[]),IntTree(3,IntTree[])])])

@test collect(PostOrderDFS([])) == Any[[]]
=#


text: abstract type AbstractMolecule end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
function mpr_delays(nt, dt, r, V, weights,idelays,g,Delta,tau,eta,J,I)
    nl, nt_, nn, nk = size(r)
    rtau = 1 / tau
    Delta_rpitau = Delta / (pi * tau)
    Threads.@threads for k=1:nk
        for t=1:nt
            for i=1:nn
                @inbounds for l=1:nl
                    acc = 0.0
                    for j=1:nn
                        acc += weights[j, i]*r[l, 256 + t - idelays[j, i], j, k]
                    end
                    r_cl = g[l, k] * acc
                    rl = r[l, 256 + t, i, k]
                    Vl = V[l, 256 + t, i, k]
                    r_noisel = r[l, 256 + t + 1, i, k]
                    V_noisel = V[l, 256 + t + 1, i, k]
                    drl = rtau * (Delta_rpitau + 2 * Vl * rl)
                    dVl = 1/tau * ( Vl^2 - pi^2 * tau^2 * rl^2 + eta + J * tau * rl + I + r_cl )
                    r[l, 256 + t + 1, i, k] = rl + dt*drl + r_noisel
                    V[l, 256 + t + 1, i, k] = Vl + dt*dVl + V_noisel
                end
            end
        end
    end
end

function mpr_benchmark()
    using BenchmarkTools
    dt = 1e-4
    nt = 1000
    k = 16
    nl = 8
    g = reshape(range(0,stop=1,length=k*nl),(nl,k))

    r = randn(Float32, nl, 2000, 76, k) * 1e-5;
    V = randn(Float32, nl, 2000, 76, k) * 1e-5;
    @. V -= 2;

    pars = [1.0,1.0,-5.0,15.0,0.0]

    weights, lengths = VirtualBrain.conn76()

    weights = Float32.(weights);
    idelays = trunc.(Int32, lengths / 10.0 / 0.1);
    maximum(idelays)

    delays(1000, dt, r, V, weights, idelays, g, pars...)

    @time delays(1000, dt, r, V, weights, idelays, g, pars...)
    @time delays(1000, dt, r, V, weights, idelays, g, pars...)
    @time delays(1000, dt, r, V, weights, idelays, g, pars...)

    1000 * k * nl / 0.138 * 1e-6
    # matches numbers for numba until boundscheck=False
end


text:
export directConv, directConv!, directConv2D!, directCrossCorrelation,directCrossCorrelation2D
export BoundaryExtension, ZeroPaddingBE, ConstantBE, PeriodicBE, MirrorBE
export boundaryExtension

# first index
const tilde_i0 = Int(1)

# ================================================================
# Used for tag dispatching, parent of available boundary extensions
# ================================================================

"""
    abstract type BoundaryExtension end

Abstract type associated to boundary extension.
"""
abstract type BoundaryExtension end

"""
    struct ZeroPaddingBE <: BoundaryExtension end

A type used to tag zero padding extension
"""
struct ZeroPaddingBE <: BoundaryExtension end

"""
    struct ConstantBE <: BoundaryExtension end

A type used to tag constant constant extension
"""
struct ConstantBE <: BoundaryExtension end

"""
    struct PeriodicBE <: BoundaryExtension end

A type used to tag periodic extension
"""
struct PeriodicBE <: BoundaryExtension end

"""
    struct MirrorBE <: BoundaryExtension end

A type used to tag mirror extension
"""
struct MirrorBE <: BoundaryExtension end

"""
    scale(λ::Int,Ω::UnitRange{Int})

Range scaling.

**Caveat:**

We do not use Julia `*` operator as it returns a step range:
```jldoctest
julia> r=6:8
6:8

julia> -2*r
-12:-2:-16
```

What we need is:

```jldoctest
julia> r=6:8
6:8

julia> scale(-2,r)
-16:-12
```
"""
function scale(λ::Int,Ω::UnitRange{Int})
    ifelse(λ>0,
           UnitRange(λ*first(Ω),λ*last(Ω)),
           UnitRange(λ*last(Ω),λ*first(Ω)))
end

#+BoundaryExtension,Internal
#
# In
# $$
# \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^*
# $$
# the computation of $\gamma[k],\ k\in\Omega^\gamma$ is splitted into two parts:
#  - one part $\Omega^\gamma \cap \Omega^\gamma_1$ *free of boundary effect*,
#  - one part $\Omega^\gamma \setminus \Omega^\gamma_1$ *that requires boundary extension* $\tilde{\beta}=\Phi(\beta,k)$
#
# *Example:*
#!DirectConvolution.compute_Ωγ1(-1:2,-2,1:20)
function compute_Ωγ1(Ωα::UnitRange{Int},
                     λ::Int,
                     Ωβ::UnitRange{Int})

    λΩα = scale(λ,Ωα)

    UnitRange(first(Ωβ)-first(λΩα),
              last(Ωβ)-last(λΩα))
end

#+BoundaryExtension,Internal
#
# Left relative complement
#
# $$
# (A\setminus B)_{\text{Left}}=[  l(A), \min{(u(A),l(B)-1)} ]
# $$
#
# *Example:*
#!DirectConvolution.relativeComplement_left(1:10,-5:5)
#
# $(A\setminus B)=\{6,7,8,9,10\}$ and the left part (elements that are
# $\in A$ but on the left side of $B$) is *empty*.
function relativeComplement_left(A::UnitRange{Int},
                                 B::UnitRange{Int})
    UnitRange(first(A),
              min(last(A),first(B)-1))
end

#+BoundaryExtension,Internal
#
# Left relative complement
#
# $$
# (A\setminus B)_{\text{Right}}=[ \max{(l(A),u(B)+1)}, u(A) ]
# $$
#
# *Example:*
#!DirectConvolution.relativeComplement_right(1:10,-5:5)
#
# $(A\setminus B)=\{6,7,8,9,10\}$ and the right part (elements that are
# $\in A$ but on the right side of $B$) is $\{6,7,8,9,10\}$
function relativeComplement_right(A::UnitRange{Int},
                                  B::UnitRange{Int})
    UnitRange(max(first(A),last(B)+1),
              last(A))
end

# ================================================================

"""
    boundaryExtension(β::AbstractArray{T,1},
                      k::Int,
                      ::Type{BOUNDARY_EXT_TYPE})

Computes extended boundary value `β[k]` given boundary extension type
`BOUNDARY_EXT_TYPE`

Available `BOUNDARY_EXT_TYPE` are:
* `ZeroPaddingBE`: zero padding
* `ConstantBE`: constant boundary extension padding
* `PeriodicBE`: periodic boundary extension padding
* `MirrorBE`: mirror symmetry boundary extension padding

The routine is robust in the sens that there is no restriction on `k`
value.

```jldoctest
julia> dom = [-5:10;];

julia> hcat(dom,map(x->boundaryExtension([1; 2; 3],x,ZeroPaddingBE),dom))'
2×16 adjoint(::Matrix{Int64}) with eltype Int64:
 -5  -4  -3  -2  -1  0  1  2  3  4  5  6  7  8  9  10
  0   0   0   0   0  0  1  2  3  0  0  0  0  0  0   0

julia> hcat(dom,map(x->boundaryExtension([1; 2; 3],x,ConstantBE),dom))'
2×16 adjoint(::Matrix{Int64}) with eltype Int64:
 -5  -4  -3  -2  -1  0  1  2  3  4  5  6  7  8  9  10
  1   1   1   1   1  1  1  2  3  3  3  3  3  3  3   3

julia> hcat(dom,map(x->boundaryExtension([1; 2; 3],x,PeriodicBE),dom))'
2×16 adjoint(::Matrix{Int64}) with eltype Int64:
 -5  -4  -3  -2  -1  0  1  2  3  4  5  6  7  8  9  10
  1   2   3   1   2  3  1  2  3  1  2  3  1  2  3   1

julia> hcat(dom,map(x->boundaryExtension([1; 2; 3],x,MirrorBE),dom))'
2×16 adjoint(::Matrix{Int64}) with eltype Int64:
 -5  -4  -3  -2  -1  0  1  2  3  4  5  6  7  8  9  10
  3   2   1   2   3  2  1  2  3  2  1  2  3  2  1   2

```

"""
function boundaryExtension end

function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{ZeroPaddingBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    if (k>=kmin)&&(k<=kmax)
        β[k]
    else
        T(0)
    end
end

function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{ConstantBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    if k<kmin
        β[kmin]
    elseif k<=kmax
        β[k]
    else
        β[kmax]
    end
end

function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{PeriodicBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmin+mod(k-kmin,1+kmax-kmin)]
end

function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{MirrorBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmax-abs(kmax-kmin-mod(k-kmin,2*(kmax-kmin)))]
end

# ================================================================

#+Convolution,Internal
function directConv!(tilde_α::AbstractArray{T,1},
                     α_offset::Int,
                     λ::Int,
                     β::AbstractArray{T,1},
                     γ::AbstractArray{T,1},
                     Ωγ::UnitRange{Int},
                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE;
                     accumulate::Bool=false)::Nothing where {T <: Number,
                                                             LeftBE <: BoundaryExtension,
                                                             RightBE <: BoundaryExtension}
    # Sanity check
    @assert λ!=0
    @assert (first(Ωγ)>=1)&&(last(Ωγ)<=length(γ))

    # Initialization
    Ωα = filter_range(length(tilde_α),α_offset)
    Ωβ = UnitRange(1,length(β))
    tilde_Ωα = 1:length(Ωα)

    if !accumulate
        for k in Ωγ
            γ[k]=0
        end
    end

    rΩγ1=intersect(Ωγ,compute_Ωγ1(Ωα,λ,Ωβ))

    # rΩγ1 part: no boundary effect
    #
    β_offset = λ*(first(Ωα)-tilde_i0)
    for k in rΩγ1
        @simd for i in tilde_Ωα
            @inbounds γ[k]+=tilde_α[i]*β[k+λ*i+β_offset]
        end
    end

    # Left part
    #
    rΩγ1_left = relativeComplement_left(Ωγ,rΩγ1)

    for k in rΩγ1_left
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*boundaryExtension(β,k+λ*i+β_offset,LeftBE)
        end
    end

    # Right part
    #
    rΩγ1_right = relativeComplement_right(Ωγ,rΩγ1)

    for k in rΩγ1_right
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*boundaryExtension(β,k+λ*i+β_offset,RightBE)
        end
    end

    nothing
end

# """
#     directConv!

# Computes a convolution.

# Inplace modification of ``\\gamma[k], k\\in\\Omega_\\gamma``

# ```math
# \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^*
# ```

# If ``k\\notin \\Omega_\\gamma``, ``\\gamma[k]`` is unmodified.

# If *accumulate=false* then an erasing step ``\\gamma[k]=0, k\\in\\Omega_\\gamma`` is performed before computation.

# If ``\\lambda=-1`` you compute a convolution, if ``\\lambda=+1`` you
# compute a cross-correlation.

# **Example:**

# ```@jldoctest

# julia> β=[1:15;]

# julia> γ=ones(Int,15)

# julia> α=LinearFilter([0,0,1],0)

# julia> directConv!(α,1,β,γ,5:10)

# julia> hcat([1:length(γ);],γ)

# ```

# """

"""
    directConv!

Computes a convolution.
```math
\\gamma[k]=\\sum\\limits_{i\\in\\Omega^\\alpha}\\alpha[i]\\beta[k+\\lambda i]
```

The following example shows how to apply inplace the `[0 0 1]` filter on `γ[5:10]`
```jldoctest
julia> β=[1:15;];

julia> γ=ones(Int,15);

julia> α=LinearFilter([0,0,1],0);

julia> directConv!(α,1,β,γ,5:10)

julia> hcat([1:length(γ);],γ)
15×2 Matrix{Int64}:
  1   1
  2   1
  3   1
  4   1
  5   7
  6   8
  7   9
  8  10
  9  11
 10  12
 11   1
 12   1
 13   1
 14   1
 15   1

```

"""
function directConv!(α::LinearFilter{T},
                     λ::Int,

                     β::AbstractArray{T,1},

                     γ::AbstractArray{T,1},
                     Ωγ::UnitRange{Int},

                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE;

                     accumulate::Bool=false)::Nothing where {T <: Number,
                                                             LeftBE <: BoundaryExtension,
                                                             RightBE <: BoundaryExtension}

    directConv!(fcoef(α),
                offset(α),
                λ,

                β,

                γ,
                Ωγ,

                LeftBE,
                RightBE,

                accumulate=accumulate)

    nothing
end

#+Convolution
#
# Computes a convolution.
#
# Convenience function that allocate $\gamma$ and compute all its
# component using [[directConv_details][]]
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
# *Example:*
#!β=[1:15;];
#!γ=ones(Int,15);
#!α=LinearFilter([0,0,1],0);
#!γ=directConv(α,1,β);
#!hcat([1:length(γ);],γ)'
#
function directConv(α::LinearFilter{T},

                    λ::Int64,

                    β::AbstractArray{T,1},

                    ::Type{LeftBE}=ZeroPaddingBE,
                    ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}
    γ = Array{T,1}(undef,length(β))

    directConv!(α,
                λ,

                β,

                γ,
                UnitRange(1,length(γ)),

                LeftBE,
                RightBE,

                accumulate=false)

    return γ
end

#+Convolution
#
# Computes a convolution.
#
# This is a convenience function where $\lambda=-1$
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
function directConv(α::LinearFilter{T},
                    β::AbstractArray{T,1},

                    ::Type{LeftBE}=ZeroPaddingBE,
                    ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}

    return directConv(α,-1,β,LeftBE,RightBE)
end

#+Convolution
#
# Computes a cross-correlation
#
# This is a convenience function where $\lambda=+1$
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
function directCrossCorrelation(α::LinearFilter{T},
                                β::AbstractArray{T,1},

                                ::Type{LeftBE}=ZeroPaddingBE,
                                ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                                      LeftBE <: BoundaryExtension,
                                                                      RightBE <: BoundaryExtension}

    return directConv(α,+1,β,LeftBE,RightBE)
end

# 2D

# +Convolution L:directConv2D_inplace
# Computes a 2D (separable) convolution.
#
# For general information about parameters, see [[directConv_details][]]
#
# α_I must be interpreted as filter for *running index I*
#
# CAVEAT: the result overwrites β
#
# TODO: @parallel
function directConv2D!(α_I::LinearFilter{T},
                       λ_I::Int,
                       α_J::LinearFilter{T},
                       λ_J::Int,

                       β::AbstractArray{T,2},

                       min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Nothing where {T<:Number}

    γ=similar(β)

    (n,m)=size(β)
    α_I_coef=fcoef(α_I)
    α_I_offset=offset(α_I)
    α_J_coef=fcoef(α_J)
    α_J_offset=offset(α_J)
    Ωγ_I = 1:n
    Ωγ_J = 1:m

    # i running (for filter)
    for j in 1:m

        directConv!(α_I_coef,
                    α_I_offset,
                    λ_I,

                    view(β,:,j),

                    view(γ,:,j),
                    Ωγ_I,

                    min_I_BE,
                    max_I_BE,

                    accumulate=false)
    end

    # j running (for filter)
    for i in 1:n

        directConv!(α_J_coef,
                    α_J_offset,
                    λ_J,

                    view(γ,i,:),

                    view(β,i,:),
                    Ωγ_J,

                    min_J_BE,
                    max_J_BE,

                    accumulate=false)
    end

    nothing
end

# +Convolution
# Computes a 2D cross-correlation
#
# This is a wrapper that calls [[directConv2D_inplace][]]
#
# *Note:* β is not modified, instead the function returns the result.
#
function directCrossCorrelation2D(α_I::LinearFilter{T},
                                  α_J::LinearFilter{T},

                                  β::AbstractArray{T,2},

                                  min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Array{T,2} where {T<:Number}

    γ=similar(β)
    γ.=β

    directConv2D!(α_I,+1,α_J,+1,γ,
                  min_I_BE,
                  max_I_BE,
                  min_J_BE,
                  max_J_BE)

    return γ
end

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
mutable struct MultiCarRacingEnv{A,T,R<:AbstractRNG} <: AbstractEnv
    N::Int
    envs::Vector{CarRacingEnv}
    action_space::A
    observation_space::Space{Vector{A}}
    state::Vector{T}
    done::Bool
    t::Int
    dt::T
    δt::T
    track::Track
    rng::R
end

"""
    MultiCarRacingEnv(;kwargs...)

# Keyword arguments
- `T = Float64`
- `N = 2                                                     # Number of cars`
- `dt = 0.1`,                                                # Time step between actions
- `δt = 0.01,                                                # Time step used for integration`
- `track = string(@__DIR__, "/car_racing_tracks/curve.csv"), # Track to load`
- 'car_params = []                                           # Vector of car parameters to use
- `rng = Random.GLOBAL_RNG`
"""
function MultiCarRacingEnv(N=2;
    T = Float64,
    dt = 0.1,
    δt = 0.01,
    track = string(@__DIR__, "/car_racing_tracks/curve.csv"),
    car_params = [],
    rng = Random.GLOBAL_RNG,
)

    length(car_params) <= N || error("# Car parameters must be ≤ # cars")

    envs = Vector{CarRacingEnv}(undef, N)
    for ii in 1:N
        if length(car_params) >= ii
            cre = CarRacingEnv(car_params, T=T, dt=dt, δt=δt, track=track, rng=rng)
        else
            cre = CarRacingEnv(T=T, dt=dt, δt=δt, track=track, rng=rng)
        end
        envs[ii] = cre
    end

    endpts_l = []
    endpts_r = []
    single_state_size = length(RLBase.state_space(envs[1]))
    obs_space_vec = Vector{ClosedInterval}(undef, N*single_state_size)
    state = zeros(T, N*single_state_size)
    for (idx, en) in enumerate(envs)
        endpts_l = [endpts_l; leftendpoint(RLBase.action_space(en))]
        endpts_r = [endpts_r; rightendpoint(RLBase.action_space(en))]
        start_idx = single_state_size*(idx-1) + 1
        end_idx = single_state_size*idx
        obs_space_vec[start_idx:end_idx] =  RLBase.state_space(en)[:]
    end
    action_space = ClosedInterval{Vector{T}}(endpts_l, endpts_r)
    observation_space = Space(obs_space_vec)

    env = MultiCarRacingEnv(
        N,
        envs,
        action_space,
        observation_space,
        state,
        false,
        0,
        dt,
        δt,
        Track(track),
        rng,
    )

    reset!(env)
    env
end

MultiCarRacingEnv{T}(; kwargs...) where {T} = MultiCarRacingEnv(; T = T, kwargs...)

# Might not want to have the same seed for every envrionment here
function Random.seed!(env::MultiCarRacingEnv, seed)
    for en in env.envs
        Random.seed!(en.rng, seed)
    end
end

RLBase.action_space(env::MultiCarRacingEnv) = env.action_space
RLBase.state_space(env::MultiCarRacingEnv{T}) where {T} = env.observation_space
RLBase.is_terminated(env::MultiCarRacingEnv) = env.done
RLBase.state(env::MultiCarRacingEnv) = env.state

function _update_states_env2envs(env::MultiCarRacingEnv)
    for (idx, en) in enumerate(env.envs)
        ss_size = length(RLBase.state_space(en))
        start_idx = ss_size*(idx-1) + 1
        end_idx = ss_size*idx
        en.state = env.state[start_idx:end_idx]
    end
end

function _update_states_envs2env(env::MultiCarRacingEnv)
    for (idx, en) in enumerate(env.envs)
        ss_size = length(RLBase.state_space(en))
        start_idx = ss_size*(idx-1) + 1
        end_idx = ss_size*idx
        env.state[start_idx:end_idx] = en.state
    end
end

function within_track(env::MultiCarRacingEnv)
    within = true
    for en in env.envs
        within = within && within_track(en).within
    end
    return within
end

function exceed_β(env::MultiCarRacingEnv)
    exceed = false
    for en in env.envs
        exceed = exceed || exceed_β(en)
    end
    return exceed
end

"""
    reward(env::MultiCarRacingEnv)
    - Sum over rewards of the ohter envrionments
    - Collision (≤ 4m): -7000
    - Penalize distance away by the distance amount
"""
function RLBase.reward(env::MultiCarRacingEnv{T}) where {T}
    rew = 0.0
    for (ii, en) in enumerate(env.envs)
        rew += reward(en)
        for jj ∈ (ii+1):env.N
            Δd = norm(env.envs[jj].state[1:2] - en.state[1:2])
            rew += -Δd
            if Δd ≤ 4.0
                rew += -11000.0
            end
        end
    end
    return rew
end

function RLBase.reset!(env::MultiCarRacingEnv{A,T}) where {A,T}
    ss_size = length(env.state)
    ind_ss_size = round(Int, length(env.state)/env.N)
    env.envs[1].state = zeros(T, ind_ss_size)
    env.envs[1].state[3] = deg2rad(90)
    env.envs[1].state[4] = 10.0
    for ii ∈ 2:env.N
        if mod(ii,2) == 0
            env.envs[ii].state[1] = ii/2*5.0
        else
            env.envs[ii].state[1] = (1-ii)/2*5.0
        end
        env.envs[ii].state[3] = deg2rad(90)
        env.envs[ii].state[4] = 10.0
    end

    _update_states_envs2env(env)
    env.t = 0
    env.done = false
    nothing
end

function RLBase.reset!(env::MultiCarRacingEnv{A,T}, state::Vector{T}) where {A,T}
    env.state = state
    _update_states_env2envs(env)
    env.t = 0
    env.done = false
    nothing
end

"""
    a::Vector{Float64}
    size = (2N,)
    a[1] = Turn angle [-max turn angle, max turn angle] (-1 right turn, +1 left turn)
    a[2] = Pedal amount (-1 = full brake, 1 = full throttle)
    a[3] = Turn angle Car 2
    a[4] = Pedal amount Car 2
    a[2N-1] = turn angle Car N
    a[2N] = Pedal amount Car N
"""
function (env::MultiCarRacingEnv)(a::Vector{Float64})
    length(a) == env.N * 2 || error("Action space of each car is of size 2")
    for (ii, en) in enumerate(env.envs)
        aᵢ = a[(2*ii-1):(2*ii)]
        _step!(en, aᵢ)
    end
    _update_states_envs2env(env)
end

function (env::MultiCarRacingEnv)(a::Vector{Int})
    env(Float64.(a))
end

function (env::MultiCarRacingEnv)(a::Matrix{Float64})
    size(a)[2] == 1 || error("Only implented for one step")
    env(vec(a))
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
# This file was generated, do not modify it. # hide
mytrace = mysim.trace
@with mytrace begin
    Ey = α .+ β .* x
    r = (y - Ey) / σ
    (;Ey, r)
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# Tests for Enhanced Russell Graph DEA Models
@testset "EnhancedRussellGraphDEAModel" begin

    # Data for ERG test
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    # Enhanced Russell Graph CRS
    deaergcrs = deaerg(X, Y, rts = :CRS)

    @test typeof(deaergcrs) == EnhancedRussellGraphDEAModel

    @test nobs(deaergcrs) == 8
    @test ninputs(deaergcrs) == 1
    @test noutputs(deaergcrs) == 1
    @test efficiency(deaergcrs) ≈ [0.4; 1.0; 0.8; 0.6; 0.4; 0.4; 0.514285714285650175; 0.2]
    @test efficiency(deaergcrs, :beta) ≈ [0.4; 1.0; 1.0; 0.6; 0.4; 0.4; 0.514285714285650175; 0.2]
    @test convert(Matrix, peers(deaergcrs)) ≈ [
        0 0.5     0    0    0    0    0    0
        0 1.0     0    0    0    0    0    0
        0 1.6     0    0    0    0    0    0
        0 3.0     0    0    0    0    0    0
        0 1.5     0    0    0    0    0    0
        0 3.5     0    0    0    0    0    0
        0 3.5     0    0    0    0    0    0
        0 2.353   0    0    0    0    0    0]

    @test slacks(deaergcrs, :X) ≈ [0; 0; 1.6; 0; 0; 0; 0; 0]
    @test slacks(deaergcrs, :Y) ≈ [1.5; 0; 0; 6.0; 4.5; 10.5; 8.5; 9.412]

    @test efficiency(deaerg(targets(deaergcrs, :X), targets(deaergcrs, :Y), rts = :CRS)) ≈ ones(8,1) atol = 1e-6

    # Test default is CRS
    @test efficiency(deaerg(X, Y)) ≈ efficiency(deaergcrs)

    # Enhanced Russell Graph VRS
    deaergvrs = deaerg(X, Y, rts = :VRS)

    @test efficiency(deaergvrs) ≈ [1.0; 1.0; 1.0; 1.0; 0.4; 0.4761904761904763; 0.8571428571428572; 0.2] atol = 1e-10
    @test efficiency(deaergvrs, :beta) ≈ [1.0; 1.0; 1.0; 1.0; 0.6; 1.0; 1.0; 0.4706]
    @test convert(Matrix, peers(deaergvrs)) ≈ [
        1.0  0     0    0    0    0    0    0
        0    1.0   0    0    0    0    0    0
        0    0     1.0  0    0    0    0    0
        0    0.0   0    1    0    0    0    0
        0    1.0   0    0    0    0    0    0
        0    1/3   2/3  0    0    0    0    0
        0    0     0    1    0    0    0    0
        0    1     0    0    0    0    0    0]

    @test slacks(deaergvrs, :X) ≈ [0; 0; 0; 0; 2.0; 22/3; 2; 5.412]
    @test slacks(deaergvrs, :Y) ≈ [0; 0; 0; 0; 2.0; 0; 0.0; 2.647]

    @test efficiency(deaerg(targets(deaergvrs, :X), targets(deaergvrs, :Y), rts = :VRS)) ≈ ones(8,1) atol = 1e-6

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaerg_crs_ref_eff = zeros(size(X, 1))
    deaerg_vrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deaerg_crs_ref_eff[i] = efficiency(deaerg(Xeval, Yeval, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaerg_vrs_ref_eff[i] = efficiency(deaerg(Xeval, Yeval, rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test deaerg_crs_ref_eff ≈ efficiency(deaergcrs)
    @test deaerg_vrs_ref_eff ≈ efficiency(deaergvrs)

    # Print
    show(IOBuffer(), deaergcrs)
    show(IOBuffer(), deaergvrs)

    # ------------------
    # Test errors
    # ------------------
    @test_throws DimensionMismatch deaerg([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaerg([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaerg([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deaerg([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deaerg([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale

    @test_throws ArgumentError efficiency(deaergcrs, :Error) # Invalid efficiency type

end


text:
using Lint.ExpressionUtils

@testset "Expressions" begin
    @test expand_trivial_calls(:(A')) == :(ctranspose(A))
end


text:
using Test
using MPI
using CLIMA.Mesh.Topologies

function main()
  MPI.Init()
  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  @assert csize == 3

  topology = StackedBrickTopology(comm, (2:5,4:6), periodicity=(false,true),
                                  boundary=((1,2),(3,4)))

  elems       = topology.elems
  realelems   = topology.realelems
  ghostelems  = topology.ghostelems
  sendelems   = topology.sendelems
  elemtocoord = topology.elemtocoord
  elemtoelem  = topology.elemtoelem
  elemtoface  = topology.elemtoface
  elemtoordr  = topology.elemtoordr
  elemtobndy  = topology.elemtobndy
  nabrtorank  = topology.nabrtorank
  nabrtorecv  = topology.nabrtorecv
  nabrtosend  = topology.nabrtosend

  globalelemtoface = [1 1 2 2 2 2
                      1 1 1 1 2 2
                      4 4 4 4 4 4
                      3 3 3 3 3 3]

  globalelemtoordr = ones(Int, size(globalelemtoface))

  globalelemtocoord = Array{Int}(undef, 2, 4, 6)
  globalelemtocoord[:,:, 1] = [2 3 2 3; 4 4 5 5]
  globalelemtocoord[:,:, 2] = [2 3 2 3; 5 5 6 6]
  globalelemtocoord[:,:, 3] = [3 4 3 4; 4 4 5 5]
  globalelemtocoord[:,:, 4] = [3 4 3 4; 5 5 6 6]
  globalelemtocoord[:,:, 5] = [4 5 4 5; 4 4 5 5]
  globalelemtocoord[:,:, 6] = [4 5 4 5; 5 5 6 6]

  globalelemtobndy = [1 1 0 0 0 0
                      0 0 0 0 2 2
                      0 0 0 0 0 0
                      0 0 0 0 0 0]

  if crank == 0
    nrealelem = 2
    globalelems = [1,2,3,4]
    elemtoelem_expect = [1 2 3 4
                         3 4 3 4
                         2 1 3 4
                         2 1 3 4]
    nabrtorank_expect = [1]
    nabrtorecv_expect = UnitRange{Int}[1:2]
    nabrtosend_expect = UnitRange{Int}[1:2]
  elseif crank == 1
    nrealelem = 2
    globalelems = [3,4,1,2,5,6]
    elemtoelem_expect = [3 4 3 4 5 6
                         5 6 3 4 5 6
                         2 1 3 4 5 6
                         2 1 3 4 5 6]
    nabrtorank_expect = [0, 2]
    nabrtorecv_expect = UnitRange{Int}[1:2, 3:4]
    nabrtosend_expect = UnitRange{Int}[1:2, 3:4]
  elseif crank == 2
    nrealelem = 2
    globalelems = [5,6,3,4]
    elemtoelem_expect = [3 4 3 4
                         1 2 3 4
                         2 1 3 4
                         2 1 3 4]
    nabrtorank_expect = [1]
    nabrtorecv_expect = UnitRange{Int}[1:2]
    nabrtosend_expect = UnitRange{Int}[1:2]
  end

  @test elems == 1:length(globalelems)
  @test realelems == 1:nrealelem
  @test ghostelems == nrealelem+1:length(globalelems)
  @test elemtocoord == globalelemtocoord[:,:,globalelems]
  @test elemtoface[:,realelems] == globalelemtoface[:,globalelems[realelems]]
  @test elemtoelem == elemtoelem_expect
  @test elemtobndy == globalelemtobndy[:,globalelems]
  @test elemtoordr == ones(eltype(elemtoordr), size(elemtoordr))
  @test nabrtorank == nabrtorank_expect
  @test nabrtorecv == nabrtorecv_expect
  @test nabrtosend == nabrtosend_expect
end

main()

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
module Kyulacs

using PyCall
using Reexport

include("pywraputils.jl")

# qulacs
include("qulacs.jl")
@reexport using .Qulacs

include("qulacs_gpu.jl")

# qulacs.circuit
include("circuit.jl")

# qulacs.gate
include("gate.jl")

# qulacs.observable
include("observable.jl")

# qulacs.quantum_operator
include("quantum_operator.jl")

# qulacs.state
include("state.jl")

end # module Kyulacs


text:
@testset "YAS style" begin
    @testset "basic" begin
        str_ = "foo(; k =v)"
        str = "foo(; k = v)"
        @test yasfmt(str_, 4, 80) == str

        str_ = "a = (arg1, arg2, arg3)"
        str = """
        a = (arg1, arg2,
             arg3)"""
        @test yasfmt(str_, 4, length(str_) - 1) == str
        @test yasfmt(str_, 4, 16) == str

        str = """
        a = (arg1,
             arg2,
             arg3)"""
        @test yasfmt(str_, 4, 15) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = [arg1, arg2, arg3]"
        str = """
        a = [arg1, arg2,
             arg3]"""
        @test yasfmt(str_, 4, length(str_) - 1) == str
        @test yasfmt(str_, 4, 16) == str

        str = """
        a = [arg1,
             arg2,
             arg3]"""
        @test yasfmt(str_, 4, 15) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = {arg1, arg2, arg3}"
        str = """
        a = {arg1,arg2,arg3}"""
        @test yasfmt(str_, 4, 20) == str

        str = """
        a = {arg1,arg2,
             arg3}"""
        @test yasfmt(str_, 4, 19) == str
        @test yasfmt(str_, 4, 15) == str

        str = """
        a = {arg1,
             arg2,
             arg3}"""
        @test yasfmt(str_, 4, 14) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = Union{arg1, arg2, arg3}"
        str = """
        a = Union{arg1,arg2,arg3}"""
        @test yasfmt(str_, 4, 25) == str

        str = """
        a = Union{arg1,arg2,
                  arg3}"""
        @test yasfmt(str_, 4, 24) == str
        @test yasfmt(str_, 4, 20) == str

        str = """
        a = Union{arg1,
                  arg2,
                  arg3}"""
        @test yasfmt(str_, 4, 19) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = fcall(arg1,arg2,arg3)"
        str = """
        a = fcall(arg1, arg2, arg3)"""
        @test yasfmt(str_, 4, length(str)) == str

        str = """
        a = fcall(arg1, arg2,
                  arg3)"""
        @test yasfmt(str_, 4, 26) == str
        @test yasfmt(str_, 4, 21) == str

        str = """
        a = fcall(arg1,
                  arg2,
                  arg3)"""
        @test yasfmt(str_, 4, 20) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = @call(arg1,arg2,arg3)"
        str = """
        a = @call(arg1, arg2, arg3)"""
        @test yasfmt(str_, 4, length(str)) == str

        str = """
        a = @call(arg1, arg2,
                  arg3)"""
        @test yasfmt(str_, 4, 26) == str
        @test yasfmt(str_, 4, 21) == str

        str = """
        a = @call(arg1,
                  arg2,
                  arg3)"""
        @test yasfmt(str_, 4, 20) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = "a = array[arg1,arg2,arg3]"
        str = """
        a = array[arg1, arg2, arg3]"""
        @test yasfmt(str_, 4, length(str)) == str

        str = """
        a = array[arg1, arg2,
                  arg3]"""
        @test yasfmt(str_, 4, 26) == str
        @test yasfmt(str_, 4, 21) == str

        str = """
        a = array[arg1,
                  arg2,
                  arg3]"""
        @test yasfmt(str_, 4, 20) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = """
        using Cassette: A, B, C"""
        str = """
        using Cassette: A,
                        B,
                        C"""
        @test yasfmt(str_, 4, 1) == str
    end

    # more complicated samples
    @testset "pretty" begin
        str_ = "comp = [a * b for a in 1:10, b in 11:20]"
        str = """
        comp = [a * b
                for a in 1:10, b in 11:20]"""
        @test yasfmt(str_, 2, length(str_) - 1, always_for_in = true) == str
        @test yasfmt(str_, 2, 34, always_for_in = true) == str

        str = """
        comp = [a * b
                for a in 1:10,
                    b in 11:20]"""
        @test yasfmt(str_, 2, 33, always_for_in = true) == str

        str = """
        comp = [a *
                b
                for a in
                    1:10,
                    b in
                    11:20]"""
        @test yasfmt(str_, 2, 1, always_for_in = true) == str

        str_ = "comp = Typed[a * b for a in 1:10, b in 11:20]"
        str = """
        comp = Typed[a * b
                     for a in 1:10, b in 11:20]"""
        @test yasfmt(str_, 2, length(str_) - 1, always_for_in = true) == str
        @test yasfmt(str_, 2, 39, always_for_in = true) == str

        str = """
        comp = Typed[a * b
                     for a in 1:10,
                         b in 11:20]"""
        @test yasfmt(str_, 2, 38, always_for_in = true) == str

        str = """
        comp = Typed[a *
                     b
                     for a in
                         1:10,
                         b in
                         11:20]"""
        @test yasfmt(str_, 2, 1, always_for_in = true) == str

        str_ = "foo(arg1, arg2, arg3) == bar(arg1, arg2, arg3)"
        str = """
        foo(arg1, arg2, arg3) ==
        bar(arg1, arg2, arg3)"""
        @test yasfmt(str, 2, length(str_)) == str_
        @test yasfmt(str_, 2, length(str_) - 1) == str
        @test yasfmt(str_, 2, 24) == str

        str = """
        foo(arg1, arg2,
            arg3) ==
        bar(arg1, arg2, arg3)"""
        @test yasfmt(str_, 2, 23) == str
        @test yasfmt(str_, 2, 21) == str

        str = """
        foo(arg1, arg2,
            arg3) ==
        bar(arg1, arg2,
            arg3)"""
        @test yasfmt(str_, 2, 20) == str
        @test yasfmt(str_, 2, 15) == str

        str = """
        foo(arg1,
            arg2,
            arg3) ==
        bar(arg1,
            arg2,
            arg3)"""
        @test yasfmt(str_, 2, 14) == str
        @test yasfmt(str_, 2, 1) == str

        str_ = """
        function func(arg1::Type1, arg2::Type2, arg3) where {Type1,Type2}
          body
        end"""
        str = """
        function func(arg1::Type1, arg2::Type2,
                      arg3) where {Type1,Type2}
          body
        end"""
        @test yasfmt(str_, 2, 64) == str
        @test yasfmt(str_, 2, 39) == str

        str = """
        function func(arg1::Type1,
                      arg2::Type2,
                      arg3) where {Type1,
                                   Type2}
          body
        end"""
        @test yasfmt(str_, 2, 31) == str
        @test yasfmt(str_, 2, 1) == str

        str_ = """
        @test TimeSpan(spike_annotation) == TimeSpan(first(spike_annotation), last(spike_annotation))"""
        str = """
        @test TimeSpan(spike_annotation) ==
              TimeSpan(first(spike_annotation), last(spike_annotation))"""
        @test yasfmt(str_, 4, length(str_) - 1) == str
        @test yasfmt(str_, 4, 63) == str
        str_ = """
        @test TimeSpan(spike_annotation) == TimeSpan(first(spike_annotation), last(spike_annotation))"""
        str = """
        @test TimeSpan(spike_annotation) ==
              TimeSpan(first(spike_annotation),
                       last(spike_annotation))"""
        @test yasfmt(str_, 4, 62) == str

        str_ = raw"""ecg_signal = signal_from_template(eeg_signal; channel_names=[:avl, :avr], file_extension=Symbol("lpcm.zst"))"""
        str = raw"""
        ecg_signal = signal_from_template(eeg_signal; channel_names = [:avl, :avr],
                                          file_extension = Symbol("lpcm.zst"))"""
        @test yasfmt(str_, 4, length(str_) - 1) == str
        @test yasfmt(str_, 4, 75) == str

        str = raw"""
        ecg_signal = signal_from_template(eeg_signal;
                                          channel_names = [:avl, :avr],
                                          file_extension = Symbol("lpcm.zst"))"""
        @test yasfmt(str_, 4, 73) == str
        str = raw"""
        ecg_signal = signal_from_template(eeg_signal;
                                          channel_names = [:avl,
                                                           :avr],
                                          file_extension = Symbol("lpcm.zst"))"""
        @test yasfmt(str_, 4, 1) == str

    end

    @testset "inline comments with arguments" begin
        str_ = """
        var = fcall(arg1,
            arg2, arg3, # comment
                            arg4, arg5)"""
        str = """
        var = fcall(arg1, arg2, arg3, # comment
                    arg4, arg5)"""
        @test yasfmt(str_, 4, 80) == str
        @test yasfmt(str_, 4, 29) == str

        str = """
        var = fcall(arg1, arg2,
                    arg3, # comment
                    arg4, arg5)"""
        @test yasfmt(str_, 4, 28) == str
        @test yasfmt(str_, 4, 23) == str

        str = """
        var = fcall(arg1,
                    arg2,
                    arg3, # comment
                    arg4,
                    arg5)"""
        @test yasfmt(str_, 4, 22) == str
        @test yasfmt(str_, 4, 1) == str

        str_ = """
        comp = [
        begin
                    x = a * b + c
                    y = x^2 + 3x # comment 1
            end
                       for a in 1:10,  # comment 2
                    b in 11:20,
           c in 300:400]"""

        str = """
        comp = [begin
                  x = a * b + c
                  y = x^2 + 3x # comment 1
                end
                for a = 1:10,  # comment 2
                    b = 11:20, c = 300:400]"""
        @test yasfmt(str_, 2, 80) == str
        @test yasfmt(str_, 2, 35) == str

        str = """
        comp = [begin
                  x = a * b + c
                  y = x^2 + 3x # comment 1
                end
                for a = 1:10,  # comment 2
                    b = 11:20,
                    c = 300:400]"""
        @test yasfmt(str_, 2, 34) == str

        str_ = """
        ys = ( if p1(x)
                 f1(x)
        elseif p2(x)
            f2(x)
        else
            f3(x)
        end for    x in xs)
        """
        str = """
        ys = (if p1(x)
                f1(x)
              elseif p2(x)
                f2(x)
              else
                f3(x)
              end
              for x in xs)
        """
        @test yasfmt(str_, 2, 80) == str

        str_ = """spike_annotation = first(ann for ann in recording.annotations if ann.value == "epileptiform_spike")"""
        str = """
        spike_annotation = first(ann
                                 for ann in recording.annotations
                                 if ann.value == "epileptiform_spike")"""
        @test yasfmt(str_, 2, 80) == str

        # only that
        str_ = "foo(a, b) = (arg1, arg2, arg3)"
        str = """
        foo(a, b) = (arg1, arg2,
                     arg3)"""
        @test yasfmt(str_, 2, length(str_) - 1) == str

        str = """
        foo(a, b) = (arg1,
                     arg2,
                     arg3)"""
        @test yasfmt(str_, 2, 1) == str

        str_ = """
        fooooooooooooooooooo(arg1, arg2,
        x -> begin
        body
        end
        )"""
        str = """
        fooooooooooooooooooo(arg1, arg2,
                             x -> begin
                                 body
                             end)"""
        @test yasfmt(str_, 4, 32) == str

    end

    @testset "inline comments with arguments" begin
        # parsing error is newline is placed front of `for` here
        str_ = "var = (x, y) for x = 1:10, y = 1:10"
        str = """
        var = (x, y) for x = 1:10,
                         y = 1:10"""
        @test yasfmt(str_, 4, length(str_) - 1) == str
    end

    @testset "invisbrackets" begin
        str_ = """
        if ((
          aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa ||
          aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa ||
          aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
        ))
          nothing
        end"""
        str = """
        if ((aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa ||
             aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa ||
             aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa))
          nothing
        end"""
        @test yasfmt(str_, 2, 92) == str
    end

    @testset "issue 189" begin
        str_ = """
    D2 = [
            (b_hat * y - delta_hat[i] * y) * gamma[i] + (b * y_hat - delta[i] * y_hat) *
                                                            gamma_hat[i] + (b_hat - y_hat) *
                                                                           delta[i] + (b - y) *
                                                                                delta_hat[i] - delta[i] * delta_hat[i]
            for i = 1:8
        ]"""
        str = """
        D2 = [(b_hat * y - delta_hat[i] * y) * gamma[i] +
              (b * y_hat - delta[i] * y_hat) * gamma_hat[i] +
              (b_hat - y_hat) * delta[i] +
              (b - y) * delta_hat[i] - delta[i] * delta_hat[i]
              for i = 1:8]"""
        @test yasfmt(str_, 2, 60) == str

    end

    @testset "issue 237" begin
        str_ = """
        for x in (arg1, arg2,)
            @info "Test"
        end"""
        str = """
        for x in (arg1, arg2)
            @info "Test"
        end"""
        @test yasfmt(str_, 4, 92) == str
    end

end


text:
#=
JoiceTypeDefinitions:
- Julia version: 1.6
- Author: connorforsythe
- Date: 2021-11-05
=#

abstract type JoiceModel end
abstract type JoiceModelType end
abstract type JoiceModelSpace end
abstract type JoiceDataStructure end
export EstimateModel, f

struct WTP <: JoiceModelSpace end
struct Preference <: JoiceModelSpace end

struct Panel <: JoiceDataStructure end
struct Cross <: JoiceDataStructure end

ColType = Union{String, Symbol, Nothing}
ColsType = Union{Vector{String}, Vector{Symbol}, Nothing, ColType}
RandColsType = Union{Dict{String, String}, Dict{Symbol, String}, Nothing}
NameType = Union{String, Nothing}

struct ModelInfo

    Data::DataFrame
    ChoiceColumn::Symbol
    Parameters::Vector{Symbol}
    ChoiceSetIDColumn::Symbol
    RandomParameters::Union{Dict{Symbol, String}, Nothing}
    RandomDraws::Union{Matrix{T}, Nothing} where {T<:Real}
    PanelIDColumn::Union{Symbol, Nothing}
    ObservationIDColumn::Symbol
    DataStructure::JoiceDataStructure
    ModelSpace::JoiceModelSpace
    Name::Union{String, Nothing}
    InterDrawCount::Int
    WeightColumn::Union{Symbol, Nothing}
    TimeCreated::DateTime
    PreferenceSearchBounds::Vector{T} where {T<:Real}

end

struct LogitPrefModel <: JoiceModel

    Data::DataFrame
    ChoiceColumn::Symbol
    Parameters::Vector{Symbol}
    ChoiceSetIDColumn::Symbol
    RandomParameters::Union{Dict{Symbol, String}, Nothing}
    PanelIDColumn::Union{Symbol, Nothing}
    Space::JoiceModelSpace
    Name::Union{String, Nothing}
    InterDrawCount::Int
    WeightColumn::Union{Symbol, Nothing}
    TimeCreated::DateTime

end

max_stars_repo_name: crforsythe/Joice-Releases

max_issues_repo_name: crforsythe/Joice-Releases

max_forks_repo_name: crforsythe/Joice-Releases


text:
using JSON

function Model_DF()::DataFrame
	DataFrame(
		bitarr = String[],
		rtol = Number[],
		objective = Number[],
		dettime = Number[],
		solve_time = Number[],
		elapsed_time = Number[],
		time_limit = Number[],
		solution_status = String[],
		term_stat = String[],
		gap = Number[],

		#this snapshot's index
		index = Int[]
	)
end

"""Record attributes of interest from a model"""
function get_DF_row(m::Model;
		index::Int=0,
		elapsed_time::Number=-1,
		getdettime=nothing)

	sol = json(value.(all_variables(m)))
	objective = objective_value(m)
	gap = MOI.get(m, MOI.RelativeGap())

	att = is_gurobi(m) ? "MIPGap" : "CPXPARAM_MIP_Tolerances_MIPGap"
	rtol = get_optimizer_attribute(m, att)

	solve_time = MOI.get(m, MOI.SolveTime())

	if !isnothing(getdettime)
		env = m.moi_backend.optimizer.model.env
		dettime_P = Ref{Float64}()
		getdettime(env, dettime_P)
		dettime = dettime_P[]
		dettime_P = nothing
		env = nothing
	else
		dettime = -1.0
	end

	att = is_gurobi(m) ? "TimeLimit" : "CPXPARAM_TimeLimit"
	time_limit = get_optimizer_attribute(m, att)
	solution_status = "$(primal_status(m))"
	term_stat = "$(termination_status(m))"
	objective = objective_value(m)

	[sol, rtol, objective, dettime, solve_time, elapsed_time, time_limit,
		solution_status, term_stat, gap, index]
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
"""
Provides storage of simulation data
"""
abstract type SimulationStore end

# Required methods:
# - open_store
# - Base.isopen(store::SimulationStore)
# - Base.close(store::SimulationStore)
# - Base.flush(store::SimulationStore)
# - get_params(store::SimulationStore)
# - initialize_problem_storage!
# - list_decision_model_keys(store::SimulationStore, problem::Symbol, container_type::Symbol)
# - list_emulation_model_keys(store::SimulationStore, container_type::Symbol)
# - list_decision_models(store::SimulationStore)
# - log_cache_hit_percentages(store::SimulationStore)
# - write_result!
# - read_result!
# - read_results
# - write_optimizer_stats!
# - read_optimizer_stats
# - get_dm_data
# - get_em_data
# - get_container_key_lookup

get_dm_data(store::SimulationStore) = store.dm_data
get_em_data(store::SimulationStore) = store.em_data


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct VariateTransformChain{
    VF <: VariateForm,
    ST <: VariateSpace,
    SF <: VariateSpace,
    N,
    NTT<:NTuple{N,VariateTransform{VF}}
} <: VariateTransform{VF,ST,SF}
    transforms::NTT
end

function VariateTransformChain(transforms::NTuple{N,VariateTransform{VF}}) where {N,VF}
    shape = varshape(first(values(transforms)))
    @argcheck all(isequal(shape), map(varshape, values(transforms)))
    ST = typeof(target_space(last(transforms)))
    SF = typeof(source_space(last(transforms)))
    NTT = typeof(transforms)
    VariateTransformChain{VF,ST,SF,N,NTT}(transforms)
end

import Base.∘

# ToDo: Utilize specialized merge operations at ends of chains:
∘(a::VariateTransform, b::VariateTransform) = VariateTransformChain((b, a))
∘(a::VariateTransform, b::VariateTransformChain) = VariateTransformChain((b, a.transforms...))
∘(a::VariateTransformChain, b::VariateTransform) = VariateTransformChain((b.transforms..., a))
∘(a::VariateTransformChain, b::VariateTransformChain) = VariateTransformChain((b.transforms..., a.transforms...))

∘(a::VariateTransform{VF}, b::IdentityVT{VF}) where VF = (@argcheck varshape(a) == varshape(b); a)
∘(a::IdentityVT{VF}, b::VariateTransform{VF}) where VF = (@argcheck varshape(a) == varshape(b); b)
∘(a::IdentityVT{VF}, b::IdentityVT{VF}) where VF = (@argcheck varshape(a) == varshape(b); a)

target_space(trafo::VariateTransformChain) = target_space(last(orig.transforms))
source_space(trafo::VariateTransformChain) = source_space(last(orig.transforms))

function _apply_trafo_chain(v::Any, prev_ladj::Real)
    (v = v, ladj = prev_ladj)
end

function _apply_trafo_chain(v::Any, prev_ladj::Real, t_first::VariateTransform, ts::VariateTransform...)
    r1 = apply_vartrafo_impl(t_first, v, prev_ladj)
    r2 = _apply_trafo_chain(r1.v, r1.ladj, ts...)
    (v = r2.v, ladj = r2.ladj)
end

Base.inv(trafo::VariateTransformChain) = VariateTransformChain(reverse(map(inv, trafo.transforms)))

ValueShapes.varshape(trafo::VariateTransformChain) = varshape(first(values(trafo.transforms)))

apply_vartrafo_impl(trafo::VariateTransformChain, v::Any, prev_ladj::Real) = _apply_trafo_chain(v, prev_ladj, trafo.transforms...)

apply_vartrafo_impl(trafo::InvVT{<:VariateTransformChain}, v::Any, prev_ladj::Real) = apply_vartrafo_impl(inv(trafo.orig), v, prev_ladj)

# ToDo: Better fix for ambiguity:
function apply_vartrafo_impl(trafo::VariateTransformChain{Univariate}, v::AbstractArray{<:Real,0}, prev_ladj::Real)
    r = apply_vartrafo_impl(trafo, v[], prev_ladj)
    (v = fill(r.v), ladj = r.ladj)
end


text:
using SparseArrays, Test
using CDCS

# Example page 4 of SeDuMi_Guide_105R5.pdf
@testset "Linear Programming example" begin
    c = [1.0, -1.0, 0.0, 0.0]
    A = [10.0 -7.0 -1.0 0.0
          1.0  0.5  0.0 1.0]
    b = [5.0, 3.0]
    sol, dual, z, status = cdcs(Matrix(A'), b, c, CDCS.Cone(0, 4), verbose=0)
    tol = 1e-2
    @test sol ≈ [47/24, 25//12, 0, 0] atol=tol rtol=tol
    @test dual ≈ [1/8, -1/4] atol=tol rtol=tol
    @test z ≈ [0, 0, 1/8, 1/4] atol=tol rtol=tol
    @test status["cost"] ≈ -1/8 atol=tol rtol=tol
end


text:
struct Scaler
    scale::Float64
    transx::Float64
    transy::Float64

    function Scaler(bbox::Array{Float64, 2})
        _scale = (abs(1 / (bbox[2] - bbox[1])))
        new(_scale, (0 - bbox[1] * _scale), (0 - bbox[1, 2] * _scale))
    end
end

function scaledpoint(p::Array{Float64,1}, scaler::Scaler)
    return [(p[1] * scaler.scale) + scaler.transx, (p[2] * scaler.scale) + scaler.transy]
end

function unscaledpoint(p::Array{Float64,1}, scaler::Scaler)
    return [(p[1] - scaler.transx) / scaler.scale, (p[2] - scaler.transy) / scaler.scale]
end

function scaledpoint2d(unscaled_point::Point2D, scaler::Scaler)
    return Point(
        ((getx(unscaled_point) * scaler.scale) + scaler.transx),
        ((gety(unscaled_point) * scaler.scale) + scaler.transy),
    )
end

function unscaledpoint2d(scaled_point::Point2D, scaler::Scaler)
    return Point(
        (getx(scaled_point) - scaler.transx) / scaler.scale,
        (gety(scaled_point) - scaler.transy) / scaler.scale,
    )
end

function unscaley(y::Float64, scaler::Scaler)
    return (y - scaler.transy) / scaler.scale
end

function unscalex(x::Float64, scaler::Scaler)
    return (x - scaler.transx) / scaler.scale
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
"""
This stream function applies periodic streaming, that is distribution
functions f_i that stream out of the domain will be placed on the opposite
side of the domain
"""
function stream(quadrature::Quadrature, f, f_new = copy(f))
    lx, ly, lq = size(f)

    @inbounds for x in 1:lx, y in 1:ly, f_idx in 1:lq
        next_x, next_y = stream_periodically_from(quadrature, x, y, lx, ly, f_idx)

        f_new[next_x, next_y, f_idx] = f[x, y, f_idx]
    end

    return f_new
end

stream!(q::Quadrature; f_new, f_old) = stream!(q, f_old, f_new)
function stream!(quadrature::Quadrature, f, f_new)
    lx, ly, lq = size(f)

    @inbounds for x in 1:lx, y in 1:ly, f_idx in 1:lq
        # Implement streaming using a gather approach
        from_x, from_y = stream_periodically_to(quadrature, x, y, lx, ly, f_idx)

        f_new[x, y, f_idx] = f[from_x, from_y, f_idx]
    end

    return
end

# function stream!(q::Quadrature, f, f_new)
#     for f_idx = 1 : length(q.weights)
#         f_new[:, :, f_idx] = circshift(f[:,:, f_idx], q.abscissae[:, f_idx]);
#     end
# end

"""
Choose the next indices which should be streamed to depending on the given
 x and y index and the direction index.
We use the global abscissae variable to determine the direction and make
 sure that the indices are bounded
"""
function stream_periodically_from(q::Quadrature, x, y, lx, ly, f_idx)
    # Note: to do circshift: we have to subtract
    next_x = x + q.abscissae[1, f_idx]
    if next_x > lx
        next_x -= lx
    elseif next_x < 1
        next_x += lx
    end

    next_y = y + q.abscissae[2, f_idx]
    if next_y > ly
        next_y -= ly
    elseif next_y < 1
        next_y += ly
    end

    return next_x, next_y
end

"""
Choose the next indices which should be streamed to depending on the given
 x and y index and the direction index.
We use the abscissae variable to determine the direction and make
 sure that the indices are bounded
"""
function stream_periodically_to(q::Quadrature, x, y, lx, ly, f_idx)
    next_x = x - q.abscissae[1, f_idx]
    next_y = y - q.abscissae[2, f_idx]

    return mod1(next_x, lx), mod1(next_y, ly)
end


text:
#=
XfoilWrapper.jl
    Written by Matthew White
    4/2/2019

    Input:  config file by name of Config.txt
    Output: Cn,Cm and polar file
=#
function xfoilWrapper(cnfg::String = "bin/config.txt")
    # import config file
    cnfg = DelimitedFiles.readdlm(cnfg, String)
    # airfoil
    uppercase(cnfg[1,2]) == "NACA" ? airfoil = cnfg[1,3] : airfoil = cnfg[1,2]
    re = cnfg[2,2] # Reynolds
    mach = cnfg[3,2] # Mach #
    changePanels = uppercase(cnfg[4,2]) .== "TRUE" # Manually change panels (otherwise runs PANE)
        panels = cnfg[5,2] # Panels to change to
    changeIter = uppercase(cnfg[6,2]) .== "TRUE" # change the number of iterations
        iter = cnfg[7,2] # iteration count to change to
    # AoA range
    aStart = parse(Float32, cnfg[8,2])
    aEnd = parse(Float32, cnfg[9,2])
    aStep = parse(Float32, cnfg[10,2])
    # XFOIL path
    xPath = cnfg[11,2]
    # Polar file path
    polarPath = cnfg[12,2]

    #   --- Open XFOIL Pipe ---
    p=open(pipeline(`$xPath`),"r+")
    run(`rm -rf $polarPath`) # remove any old polar file
    # Write to pipe
    length(airfoil) <= 5 ? write(p,"naca $airfoil\n") : write(p,"load $airfoil\n")
    changePanels ? write(p,"ppar $panels\n\n") : write(p,"pane\n")
    write(p,"oper\n")
    write(p,"visc $re\n")
    if changeIter write(p,"iter $iter\n") end
    write(p,"mach $mach\n")
    write(p,"pacc\n")
    write(p,"$polarPath\n\n")
    write(p,"aseq $aStart $aEnd $aStep\n")
    write(p,"\nquit")

    close(p)

    while !isfile(polarPath) end # wait unitl polar file has been created
    while filesize(polarPath) < 1000 end # wait until at least one value has been written

    #   --- Reading the Polar File ---
    prevSize = 0
    sec = 0
    n = 10
    while true # wait unit the polar file is done changing
        global polar = DelimitedFiles.readdlm(polarPath,Float64,skipstart = 12)
        if polar[end,1] > aEnd - aStep
            println("Done writing polar file")
            break #File is fully loaded
        else #check every n seconds for stable file size
            sleep(1)
            sec += 1
            if prevSize == filesize(polarPath) && sec == n
                println("Method could not converge for the final alphas")
                println("   Maximum angle of attack at $(polar[end,1]) degrees")
                break # after n seconds if filesize is constant then use current polar
            elseif sec >= n
                sec = 0 # reset sec every n seconds
            end
            prevSize = filesize(polarPath)
        end
    end

    # polar: alpha | Cl | CD | CDp | Cm | Top_Xtr | Bot_Xtr
    # Cn = Cl*cos(a) + Cd*sin(a)
    Cn = polar[:,2].*cosd.(polar[:,1]) + polar[:,3].*sind.(polar[:,1])
    Cm = polar[:,5]

    Cn,Cm

end


text:
using oneAPI.oneL0

@testset "driver" begin

drvs = drivers()
@assert !isempty(drvs)
drv = first(drvs)
show(devnull, drv)

api_version(drv)

properties(drv)
ipc_properties(drv)
#extension_properties(drv)

end

drv = first(drivers())

@testset "device" begin

devs = devices(drv)
@assert !isempty(devs)
dev = first(devs)
show(devnull, dev)

properties(dev)
compute_properties(dev)
module_properties(dev)
memory_properties(dev)
memory_access_properties(dev)
cache_properties(dev)
image_properties(dev)
p2p_properties(dev, dev)

end

dev = first(devices(drv))

@testset "context" begin

ctx = ZeContext(drv)
show(devnull, ctx)

#status(ctx)

end

ctx = ZeContext(drv)

@testset "command" begin

groups = command_queue_groups(dev)
@test !isempty(groups)

groups = compute_groups(dev)
group = first(groups)

queue = ZeCommandQueue(ctx, dev, group.ordinal)

list = ZeCommandList(ctx, dev, group.ordinal)
close(list)
execute!(queue, [list])
synchronize(queue)
reset(list)

list = ZeCommandList(ctx, dev, group.ordinal) do list
    @test list isa ZeCommandList
end

execute!(queue) do list
    @test list isa ZeCommandList
end

end

group = first(compute_groups(dev))
queue = ZeCommandQueue(ctx, dev, group.ordinal)

@testset "event" begin

ZeEventPool(ctx, 1)
ZeEventPool(ctx, 1, dev)

pool = ZeEventPool(ctx, 1)

event = pool[1]
@test !query(event)

signal(event)
ZeCommandList(ctx, dev, group.ordinal) do list
    append_signal!(list, event)
end
@test query(event)

wait(event, 1)
ZeCommandList(ctx, dev, group.ordinal) do list
    append_wait!(list, event)
end

reset(event)
ZeCommandList(ctx, dev, group.ordinal) do list
    append_reset!(list, event)
end

# timed_pool = ZeEventPool(ctx, 1; flags=oneL0.ZE_EVENT_POOL_FLAG_KERNEL_TIMESTAMP)
# timed_event = timed_pool[1]
# @test kernel_timestamp(timed_event).global.start == nothing
# @test kernel_timestamp(timed_event).context.start == nothing
# signal(timed_event) # FIXME: A kernel timestamp event can only be signaled from zeCommandListAppendLaunchKernel et al. functions
# @test kernel_timestamp(timed_event).global.start != nothing
# @test kernel_timestamp(timed_event).context.start != nothing

end

@testset "barrier" begin

pool = ZeEventPool(ctx, 1)
event = pool[1]

ZeCommandList(ctx, dev, group.ordinal) do list
    append_barrier!(list)
    append_barrier!(list, event)
    append_barrier!(list, event, event)
end

#device_barrier(dev)    # unsupported

end

@testset "module" begin

data = read(joinpath(@__DIR__, "dummy.spv"))
mod = ZeModule(ctx, dev, data)

@test length(kernels(mod)) == 2
@test haskey(kernels(mod), "foo")
@test !haskey(kernels(mod), "baz")
kernel = kernels(mod)["foo"]

suggest_groupsize(kernel, 1024)
groupsize!(kernel, 1)
groupsize!(kernel, (1,))
groupsize!(kernel, (1, 1))
groupsize!(kernel, (1, 1, 1))

kernel = kernels(mod)["bar"]
arguments(kernel)[1] = Int32(42)

@test indirect_access(kernel) == 0
indirect_access!(kernel, oneL0.ZE_KERNEL_INDIRECT_ACCESS_FLAG_DEVICE)
@test indirect_access(kernel) == oneL0.ZE_KERNEL_INDIRECT_ACCESS_FLAG_DEVICE

# oneapi-src/level-zero#55
if !parse(Bool, get(ENV, "ZE_ENABLE_PARAMETER_VALIDATION", "false"))
    attrs = source_attributes(kernel)
    @test isempty(attrs)
end

props = properties(kernel)
@test props.numKernelArgs == 1
@test props.requiredGroupSize isa oneL0.ZeDim3

@testset "kernel execution" begin

ZeCommandList(ctx, dev, group.ordinal) do list
    append_launch!(list, kernel, 1)
end

pool = ZeEventPool(ctx, 2)
signal_event = pool[1]
wait_event = pool[2]

execute!(queue) do list
    append_launch!(list, kernel, 1, signal_event, wait_event)
end
@test !query(signal_event)

signal(wait_event)
synchronize(queue)
@test query(signal_event)

end

end

@testset "memory" begin

buf = device_alloc(ctx, dev, 1024)
props = properties(buf)
@test props.device == dev
@test props.type == oneL0.ZE_MEMORY_TYPE_DEVICE
@test_throws ArgumentError convert(Ptr{Cvoid}, buf)
ptr = convert(ZePtr{Cvoid}, buf)
@test lookup_alloc(ctx, ptr) isa typeof(buf)
free(buf)

buf = host_alloc(ctx, 1024)
props = properties(buf)
@test props.type == oneL0.ZE_MEMORY_TYPE_HOST
ptr = convert(ZePtr{Cvoid}, buf)
@test lookup_alloc(ctx, ptr) isa typeof(buf)
ptr = convert(Ptr{Cvoid}, buf)
@test lookup_alloc(ctx, ptr) isa typeof(buf)
free(buf)

buf = shared_alloc(ctx, dev, 1024)
props = properties(buf)
@test props.type == oneL0.ZE_MEMORY_TYPE_SHARED
ptr = convert(ZePtr{Cvoid}, buf)
@test lookup_alloc(ctx, ptr) isa typeof(buf)
ptr = convert(Ptr{Cvoid}, buf)
@test lookup_alloc(ctx, ptr) isa typeof(buf)
free(buf)

end

@testset "copy" begin

let src = rand(Int, 1024)
    chk = ones(Int, length(src))

    dst = device_alloc(ctx, dev, sizeof(src))

    execute!(queue) do list
        append_copy!(list, pointer(dst), pointer(src), sizeof(src))
        append_barrier!(list)
        append_copy!(list, pointer(chk), pointer(dst), sizeof(src))
    end
    synchronize(queue)
    @test chk == src

    execute!(queue) do list
        pattern = [42]
        append_fill!(list, pointer(dst), pointer(pattern), sizeof(pattern), sizeof(src))
        append_barrier!(list)
        append_copy!(list, pointer(chk), pointer(dst), sizeof(src))
    end
    synchronize(queue)
    @test all(isequal(42), chk)

    free(dst)
end

for buf in [device_alloc(ctx, dev, 1024),
            host_alloc(ctx, 1024),
            shared_alloc(ctx, dev, 1024)]
    execute!(queue) do list
        append_prefetch!(list, pointer(buf), sizeof(buf))

        append_advise!(list, dev, pointer(buf), sizeof(buf),
                       oneL0.ZE_MEMORY_ADVICE_SET_READ_MOSTLY)
    end

    free(buf)
end

end

@testset "residency" begin

for buf in [device_alloc(ctx, dev, 1024),
            host_alloc(ctx, 1024),
            shared_alloc(ctx, dev, 1024)]
    make_resident(ctx, dev, buf)
    evict(ctx, dev, buf)
    make_resident(ctx, dev, buf, 1024)
    evict(ctx, dev, buf, 1024)
    free(buf)
end

end